ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
remesh.f90
Go to the documentation of this file.
1 subroutine remesh(a, nrnew, npnew, meshno, cx, cy, xaxis, yaxis, &
2  nax, rax, sax)
3 !---------------------------------------------------------------------
4 ! form the system of new finite elements as flux coordinates using
5 ! the exact interpolation
6 ! xx,yy,psi : on exit contain the values on the new grid
7 ! nr,np : the number radial and poloidal points in the old grid
8 ! nrnew,npnew : ,, ,, in the new grid
9 ! xaxis : position of magnetic axis
10 ! nax1,nax2 : nodenumbers of element with magnetic axis
11 ! rax : r value of magnetic axis
12 !----------------------------------------------------------------------
13 
14  use itm_types
15  use mod_parameter
16  use mod_corners
17  use mod_dat
18  use mod_mesh
19  use mod_meshacc
20  use mod_nodes
21  use mod_helena_io, only : out_he, iu6
22  use mod_interpolation
23  use mod_output
24  use matrix_calculations
25  use math_functions
26  use mod_profiles
27  use helena_spline
28 
29  implicit none
30 
31  real(r8), intent(in) :: a, xaxis, yaxis, rax, sax
32  integer(itm_i4), intent(in) :: nrnew, npnew, nax
33  integer(itm_i4), intent(inout) :: meshno
34 
35  real(r8), intent(inout) :: cx, cy
36 
37  real(r8), parameter :: psitol = 1.e-8_r8
38  real(r8), parameter :: toltht = 1.e-10_r8
39  real(r8), parameter :: tolpsi = 1.e-10_r8
40 
41  real(r8) :: thtkn(npnew)
42  real(r8) :: factas
43  real(r8) :: rpsi, psid, dpsid, ddpsid
44  real(r8) :: psival, thtval
45  real(r8) :: rr, r, r2, r3
46  real(r8) :: psim, psimr, psip, psipr
47  real(r8) :: psimin, psimax
48  real(r8) :: ps, dpsir, dpsis
49  real(r8) :: a0, a1, a2, a3
50  real(r8) :: zx, dxr, dxs, zxr, zxs
51  real(r8) :: zy, dyr, dys, zyr, zys
52  real(r8) :: tht, tht0, tht1, tht2
53  real(r8) :: dd, ss
54  real(r8) :: x, xr, xs, xrs, xrr, xss
55  real(r8) :: y, yr, ys, yrs, yrr, yss
56  real(r8) :: zpsi, zpsir, zpsis, zpsirs, zpsirr, zpsiss
57  real(r8) :: rad
58  real(r8) :: thx, thy, thxx, thxy, thyy
59  real(r8) :: ths, thr, thrr, thrs, thss
60  real(r8) :: ptjac, ptjr, ptjs
61  real(r8) :: rt, st
62  real(r8) :: rpt, spt, xpt, ypt
63  real(r8) :: xyjac
64  real(r8) :: psix, psiy
65 ! real :: dum, ddum, dddum
66  real(r8) :: etap, ejac
67  real(r8) :: fvec1, fvec2
68  real(r8) :: fjac11, fjac12, fjac21, fjac22
69  real(r8) :: dis, dr, ds, dtot
70  real(r8) :: eqpsi
71  real(r8) :: s, psis, psir, psrs, psrr, psss
72  real(r8) :: rx, ry, sx, sy
73  real(r8) :: psixx, psiyy, psiyy0
74  real(r8) :: tn, tn2, cn
75  real(r8) :: tol
76  real(r8) :: trr, tss, ttht, tx, txr, txs, ty, tyr, tys, tpsi, tpsir, tpsis
77  real(r8) :: step, stepr, steps
78  real(r8) :: dis_old
79  real(r8) :: psi_bar, psi_min, psi_max, psi_cnt
80  real(r8) :: rr_exclusion_zone = +0.00_r8
81  real(r8) :: ds_old
82  integer(itm_i4) :: tjn, tjns, tnode, ttn1, ttn2, ttn3, ttn4
83  integer(itm_i4) :: nnsave1, nnsave2
84  integer(itm_i4) :: psi_prob, psi_prob_arr(1)
85  integer(itm_i4) :: nsym(2 * nr - 2), ssym(2 * nr - 2)
86  integer(itm_i4) :: i, j, k, n, ns, itn
87  integer(itm_i4) :: nn
88  integer(itm_i4) :: n1, n2, n3, n4
89  integer(itm_i4) :: node, in, jn, j0, n0
90  integer(itm_i4) :: ifail
91  integer(itm_i4) :: iprev, jprev
92  integer(itm_i4) :: itmax, ittest
93  integer(itm_i4) :: jindex
94  logical :: chagr(nrnew * npnew)
95  logical :: nobrack, found
96  type(spline_coefficients) :: f_spline
97  real(r8), dimension(3) :: abltg
98  real(r8) :: alfa_t, beta_t, tmp1
99  integer(itm_i4) :: typ_t
100  logical :: modified_psi
101  logical :: crashed
102 
103  crashed = .false.
104  factas = 1._r8
105  ds = 0.0_r8
106  if (ias == 0) factas = 2._r8
107 
108  do i = 1, nr - 1
109  nsym(i) = (i - 1) * (np - 1) + 1
110  ssym(i) = -1
111  nsym(nr - 1 + i) = (nr - 1) * (np - 1) - (i - 1) * (np - 1)
112  ssym(nr - 1 + i) = 1
113  end do
114 
115  if (imesh == 2) then
116  if (nr /= nrnew) then
117 !------------------- reallocation of radial grid for new mesh accumulation
118  deallocate(sg)
119  deallocate(dsg)
120  deallocate(ddsg)
121  allocate(sg(nrnew))
122  allocate(dsg(nrnew))
123  allocate(ddsg(nrnew))
124  call mesh_accumulation(nrnew, n_acc_points, s_acc(1 : n_acc_points), &
125  sig(1 : n_acc_points), weights(1 : n_acc_points), equidistant, &
126  sg, dsg, ddsg, .false.)
127  end if
128  do i = 1, nrnew
129  psikn(nrnew - i + 1) = sg(i)**2
130  dpsikn(nrnew - i + 1) = 2._r8 * sg(i) * dsg(i)
131  ddpsikn(nrnew - i + 1) = 2._r8 * sg(i) * ddsg(i) + 2._r8 &
132  * dsg(i)**2
133  radpsi(nrnew - i + 1) = dble(i - 1) / dble(nrnew - 1)
134  end do
135  else if (imesh == 3) then
136  if (nr /= nrnew) then
137 !------------------- reallocation of radial grid for new mesh accumulation
138  deallocate(sg)
139  deallocate(dsg)
140  deallocate(ddsg)
141  allocate(sg(nrnew))
142  allocate(dsg(nrnew))
143  allocate(ddsg(nrnew))
144  call mesh_accumulation(nrnew, n_acc_points, s_acc(1 : n_acc_points), &
145  sig(1 : n_acc_points), weights(1 : n_acc_points), equidistant, &
146  sg, dsg, ddsg, .true.)
147  end if
148  do i = 1, nrnew
149  psikn(nrnew - i + 1) = sg(i)**2
150  dpsikn(nrnew - i + 1) = 2._r8 * sg(i) * dsg(i)
151  ddpsikn(nrnew - i + 1) = 2._r8 * sg(i) * ddsg(i) + 2._r8 &
152  * dsg(i)**2
153  radpsi(nrnew - i + 1) = dble(i - 1) / dble(nrnew - 1)
154  end do
155  else
156  do i = 1, nrnew
157  rpsi = dble(i - 1) / dble(nrnew - 1)
158  call radial_mesh(rpsi, psid, dpsid, ddpsid)
159  psikn(nrnew - i + 1) = psid
160  dpsikn(nrnew - i + 1) = dpsid
161  ddpsikn(nrnew - i + 1) = ddpsid
162  radpsi(nrnew - i + 1) = rpsi
163  end do
164  end if
165  radpsi(nrnew) = 0._r8
166  psikn(nrnew) = 0._r8
167  do j = 1, npnew
168  thtkn(j) = (1._r8 + dble(ias)) * pi * dble(j - 1) / dble(npnew -1)
169  end do
170  meshno = meshno + 1
171 !------------------------- update old mesh for the first iteration ----
172  do i = 1, nrnew * npnew
173  chagr(i) = .false.
174  end do
175  do i = 1, nr * np
176  do k = 1, 4
177  xxold(k, i) = xx(k, i)
178  yyold(k, i) = yy(k, i)
179  psiold(4 * (i - 1) + k) = psi(4 * (i - 1) + k)
180  end do
181  end do
182 
183  psi_bar = 0._r8
184  psi_cnt = 0._r8
185  psi_min = 1.e30_r8
186  psi_max = -1.e30_r8
187  rr = rr_exclusion_zone
188  ss = 0._r8
189  do nn = (nr - 2) * (np - 1) + 1, (nr - 1) * (np - 1)
190  in = (nn - 1) / (np - 1) + 1
191  jn = mod(nn - 1, np - 1) + 1
192  call set_node_number(nn, n1, n2, n3, n4)
193 ! write(*,*) 'n1, n2, n3, n4 = ', n1, n2, n3, n4
194  call interpolation(0, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 - 1) + 1), &
195  psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), rr, ss, psim)
196 ! write(*,*) 'in, jn, psi = ', in, jn, psim
197  psi_bar = psi_bar + psim
198  psi_cnt = psi_cnt + 1._r8
199  psi_min = min(psi_min, psim)
200  psi_max = max(psi_max, psim)
201  end do
202  psi_bar = psi_bar / psi_cnt
203  if (verbosity > 1) write(iu6, *) 'psi_bar, psi_min, psi_max = ', psi_bar, &
204  psi_min, psi_max
205  modified_psi = .false.
206  if (psikn(nrnew) < psi_bar) then
207  i = 1
208  do while (psikn(i) > psi_bar)
209  i = i + 1
210  end do
211  i = i - 1
212  else
213  i=nrnew
214  end if
215  if ( i >= 1) then
216  if (psikn(i) < psi_max) then
217  modified_psi = .true.
218  if (verbosity > 3) write(iu6, *) "Before: ", psikn(1 : i)
219  psikn(1 : i) = 1 - (1 - psikn(1 : i)) * (1 - psi_max) / (1 - psikn(i))
220  if (verbosity > 3) write(iu6, *) "After: ", psikn(1 : i)
221  end if
222  end if
223  i = i + 1
224  if (i <= nrnew) then
225  if (psikn(i) > psi_min) then
226  modified_psi = .true.
227  if (verbosity > 3) write(iu6, *) "Before: ", psikn(i : nrnew)
228  psikn(i : nrnew) = psikn(i : nrnew) * psi_min / psikn(i)
229  if (verbosity > 3) write(iu6, *) "After: ", psikn(i : nrnew)
230  end if
231  end if
232 
233  if (verbosity > 0) write(iu6, *) 'modified_psi = ', modified_psi
234  if (modified_psi) then
235  call allocate_spline_coefficients(f_spline, nrnew)
236  alfa_t = 0._r8
237  beta_t = 0._r8
238  typ_t = 3
239  call spline(nrnew, radpsi, psikn, alfa_t, beta_t, typ_t, f_spline)
240  if (verbosity > 5) write(iu6, *) 'Compare psikn and its derivatives ' &
241  // 'with spline calculations of same'
242  do i = 1, nrnew
243  tmp1 = spwert(nrnew, radpsi(i), f_spline, radpsi, abltg, 2)
244  if (verbosity > 5) write(iu6,*) radpsi(i), psikn(i), tmp1, dpsikn(i), &
245  abltg(1), ddpsikn(i), abltg(2)
246  dpsikn(i)=abltg(1)
247  ddpsikn(i)=abltg(2)
248  end do
249  call deallocate_spline_coefficients(f_spline)
250  end if
251 
252 !-- deallocate old xx, yy, and psi and re-allocate with nrnew, npnew
253  deallocate(xx, yy, psi)
254  allocate(psi(4 * nrnew * npnew))
255  allocate(xx(4, nrnew * npnew))
256  allocate(yy(4, nrnew * npnew))
257 
258  xx = 0._r8
259  yy = 0._r8
260  psi = 0._r8
261 
262 !--------------------------- loop over all flux surfaces -----------
263  do i = 1, nrnew - 1
264  psival = psikn(i)
265  if (verbosity > 6) write(iu6, *) ' finding surface : ', psival, nax
266 !--------------------------- find starting point of flux contour
267  do n = nax, 1, 1 - np
268  call set_node_number(n, n1, n2, n3, n4)
269 !------------------------------------- quad. eq for r value at minimum -
270  rr = -1._r8
271  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 - 1) &
272  + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
273  rr, sax, psim, psimr, dpsis)
274  rr = +1._r8
275  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 - 1) &
276  + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
277  rr, sax, psip, psipr, dpsis)
278  tol = 1._r8 + 1.e-5_r8
279  call quad_equation(psim, psimr, psip, psipr, tol, r, ifail)
280  if (abs(r) > 1._r8) then
281  psimin = min(psim, psip)
282  psimax = max(psim, psip)
283  else
284  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 - 1) &
285  + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
286  r, sax, ps, dpsir, dpsis)
287  psimin = min(min(psim, ps), psip)
288  psimax = max(max(psim, ps), psip)
289  endif
290  if (psival >= psimin - psitol .and. psival <= psimax + psitol) then
291  a3 = (psim + psimr - psip + psipr) / 4._r8
292  a2 = (-psimr + psipr) / 4._r8
293  a1 = (-3._r8 * psim - psimr + 3._r8 * psip - psipr) / 4._r8
294  a0 = (2._r8 * psim + psimr + 2._r8 * psip - psipr) / 4._r8 &
295  - psival
296  call solve_cubic(a0, a1, a2, a3, rr, r2, r3, ifail)
297  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 - 1) &
298  + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
299  rr, sax, ps, dpsir, dpsis)
300  call interpolation(1, xxold(1, n1), xxold(1, n2), xxold(1, n3), &
301  xxold(1, n4), rr, sax, zx, dxr, dxs)
302  if (zx < xaxis) then
303  if (verbosity > 5) then
304  write(iu6, *) 'node on wrong side : ', zx, xaxis
305  write(iu6, *) 'Roots are ', rr, r2, r3
306  write(iu6, *) 'First solution for psi had value ',ps
307  endif
308  rr = r2
309  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 &
310  - 1) + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
311  rr, sax, ps, dpsir, dpsis)
312  call interpolation(1, xxold(1, n1), xxold(1, n2), xxold(1, n3), &
313  xxold(1, n4), rr, sax, zx, dxr, dxs)
314  if (verbosity > 5) then
315  write(iu6, *) 'Second solution for psi had value ', ps
316  write(iu6, *) 'zx, xaxis = ', zx, xaxis
317  endif
318  end if
319  if (abs(rr) <= 1._r8 + 1.e-5_r8 .and. zx >= xaxis) goto 45
320 !------------------ special case psi = 1._r8
321  if (i == 1 .and. abs(ps - psival) < 1.e-5_r8) then
322  rr = -1._r8
323  goto 45
324  end if
325  end if
326  end do
327  if (verbosity > 0) &
328  write(iu6, *) ' starting value not found : ', psival, rr, ps
329 ! If this happens program should be stopped!
330  45 continue
331 !--------------------------- starting position found ----------------
332  call interpolation(1, yyold(1, n1), yyold(1, n2), yyold(1, n3), &
333  yyold(1, n4), rr, sax, zy, dyr, dys)
334  tht0 = atan2(zy - yaxis, zx - xaxis)
335  if (tht0 < 0._r8) tht0 = tht0 + 2._r8 * pi
336 !----------------------------- trace flux surface to find theta values
337  dd = 0.25_r8 * (psival)**0.33_r8 / factas
338  iprev = 0
339  jprev = 0
340  ss = sax
341  nn = n
342  tht1 = tht0
343  tht2 = tht0
344  itmax = 250000 ! was 25000 ! was 2500
345  ittest = 0
346  nobrack = .true.
347  nnsave1 = nn
348  do j = 1, npnew
349 47 if (ias == 1) then
350  jindex = mod(int(tht0 / (2._r8 * pi) * npnew) + j, npnew) + 1
351  else
352  jindex = j
353  end if
354  thtval = thtkn(jindex)
355  if (crashed .and. verbosity > 6) &
356  write(iu6, *) "value of theta sought ", thtval
357  found = .false.
358 
359 !---------------------------------- treat theta = pi as special point
360  if (j == npnew .and. ias == 0) then
361 
362  do ns = 1, 2 * nr - 2
363  n = nsym(ns)
364  ss = ssym(ns)
365  call set_node_number(n, n1, n2, n3, n4)
366  rr = -1._r8
367  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 &
368  - 1) + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
369  rr, ss, psim, psimr, dpsis)
370  rr = 1._r8
371  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 &
372  - 1) + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
373  rr, ss, psip, psipr, dpsis)
374  a3 = (psim + psimr - psip + psipr) / 4._r8
375  a2 = (-psimr + psipr) / 4._r8
376  a1 = (-3._r8 * psim - psimr + 3._r8 * psip - psipr) / 4._r8
377  a0 = (2._r8 * psim + psimr + 2._r8 * psip - psipr) &
378  / 4._r8 - psival
379  call solve_cubic(a0, a1, a2, a3, rr, r2, r3, ifail)
380  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 &
381  - 1) + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
382  rr, ss, ps, dpsir, dpsis)
383  call interpolation(1, xxold(1, n1), xxold(1, n2), xxold(1, n3), &
384  xxold(1, n4), rr, ss, zx, dxr, dxs)
385  if (zx > xaxis) then
386  if (verbosity > 3) write(iu6, *) ' node on wrong side : ', zx, &
387  xaxis
388  rr = r2
389  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 &
390  - 1) + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 &
391  - 1) + 1), rr, ss, ps, dpsir, dpsis)
392  call interpolation(1, xxold(1, n1), xxold(1, n2), xxold(1, n3), &
393  xxold(1, n4), rr, ss, zx, dxr, dxs)
394  end if
395  if (abs(rr) <= 1._r8 + 1.e-5_r8 .and. zx < xaxis &
396  .and. zx >= -1._r8 - 1.e-8_r8) then
397  found = .true.
398  nn = n
399  nobrack= .false.
400  goto 145
401  end if
402  end do
403  end if
404  145 continue
405 
406  if (crashed .and. verbosity > 6) &
407  write(iu6, *) "found, nobrack = ", found, nobrack
408  if (crashed .and. verbosity > 6) &
409  write(iu6, *) "nn before itn loop ", nn
410  nnsave2=nn
411 
412  step = 0.5_r8
413 ! step=real(min(np,nr))/2.0_R8
414  dis_old = 0._r8
415  do itn = 1, itmax
416  in = (nn - 1) / (np - 1) + 1
417  jn = mod(nn - 1, np - 1) + 1
418 ! if (itn == 100) then
419 ! step = step / 10._r8
420 ! if (verbosity > 0) write(iu6, *) "step changed to ", step
421 ! end if
422  if (itn == 1000) then
423  step = step / 10._r8
424  if (verbosity > 0) write(iu6,*) "step changed to ", step
425  end if
426  if (itn == 10000) then
427  step = step / 10._r8
428  if (verbosity > 0) write(iu6,*) "step changed to ", step
429  end if
430  if (itn == 100000) then
431  step = step / 10._r8
432  if (verbosity > 0) write(iu6,*) "step changed to ", step
433  end if
434  stepr = step
435  steps = step
436  if (crashed .and. verbosity > 6) write(iu6,*) "itn = ", itn
437  call set_node_number(nn, n1, n2, n3, n4)
438  call interpolation(1, xxold(1, n1), xxold(1, n2), xxold(1, n3), &
439  xxold(1, n4), rr, ss, zx, zxr, zxs)
440  call interpolation(1, yyold(1, n1), yyold(1, n2), yyold(1, n3), &
441  yyold(1, n4), rr, ss, zy, zyr, zys)
442  call interpolation(1, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 - 1) &
443  + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
444  rr, ss, ps, dpsir, dpsis)
445  tht = atan2(zy - yaxis, zx - xaxis)
446  if (tht < 0._r8) tht = tht + 2._r8 * pi
447  if (crashed .and. verbosity > 6) &
448  write(iu6, *) "thtval, tht = ", thtval, tht
449  if (nobrack) then
450  tht1 = tht2
451  tht2 = tht
452  else if (found) then
453  tht1 = tht
454  end if
455  if (crashed .and. verbosity > 6) &
456  write(iu6, *) "theta error and tol = ", &
457  min(abs(thtval - tht), abs(thtval - tht + 2._r8 * pi), &
458  abs(thtval - tht - 2._r8 * pi)), toltht
459  if (crashed .and. verbosity > 6) &
460  write(iu6, *) "psi error and tol = ", abs(psival - ps), &
461  tolpsi
462  if ((abs(thtval - tht) < toltht .or. abs(thtval - tht &
463  + 2._r8 * pi) < toltht .or. abs(thtval - tht - 2._r8 * pi) &
464  < toltht) .and. abs(psival - ps) < tolpsi) then
465 !
466 !----------------------------- node located ---------------------
467 !
468 ! if (j == np) then
469 ! if (verbosity > 3) then
470 ! write(iu6, *) ' node located : '
471 ! write(iu6, 41) j, zx, zy, ps, tht, ps - psival, tht - thtval
472 ! end if
473 ! end if
474  if (crashed .and. verbosity > 6) write(iu6,*) "Found"
475  found = .true.
476  if (abs(rr) > 1._r8 + tolpsi .or. abs(ss) > 1._r8 &
477  + toltht) then
478  if (verbosity > 1) &
479  write(iu6, *) ' warning : ', rr, ss
480  end if
481  node = (i - 1) * npnew + jindex
482  call set_node_number(nn, n1, n2, n3, n4)
483  call interpolation(2, xxold(1, n1), xxold(1, n2), xxold(1, n3), &
484  xxold(1, n4), rr, ss, x, xr, xs, xrs, xrr, xss)
485  call interpolation(2, yyold(1, n1), yyold(1, n2), yyold(1, n3), &
486  yyold(1, n4), rr, ss, y, yr, ys, yrs, yrr, yss)
487  call interpolation(2, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 &
488  - 1) + 1), psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), &
489  rr, ss, zpsi, zpsir, zpsis, zpsirs, zpsirr, zpsiss)
490 !-------------------------------------------------------------------
491  rad = (x - xaxis)**2 + (y - yaxis)**2
492  thy = (x - xaxis) / rad
493  thx = -(y - yaxis) / rad
494  thxx = 2._r8 * (y - yaxis) * (x - xaxis) / rad**2
495  thyy = -thxx
496  thxy = ((y - yaxis)**2 - (x - xaxis)**2) / rad**2
497  ths = thx * xs + thy * ys
498  thr = thx * xr + thy * yr
499  thrr = thxx * xr * xr + 2._r8 * thxy * xr * yr + thx * xrr &
500  + thyy * yr * yr + thy * yrr
501  thrs = thxx * xr * xs + thxy * xr * ys + thx * xrs &
502  + thxy * yr * xs + thyy * yr * ys + thy * yrs
503  thss = thxx * xs * xs + 2._r8 * thxy * xs * ys + thx * xss &
504  + thyy * ys * ys + thy * yss
505  ptjac = zpsir * ths - zpsis * thr
506  rt = -zpsis / ptjac
507  st = zpsir / ptjac
508  ptjr = zpsirr * ths + zpsir * thrs - zpsirs * thr - zpsis &
509  * thrr
510  ptjs = zpsirs * ths + zpsir * thss - zpsiss * thr - zpsis &
511  * thrs
512  rpt = (-ptjr * ths / ptjac**2 + thrs / ptjac) * rt &
513  + (-ptjs * ths / ptjac**2 + thss / ptjac) * st
514  spt = (ptjr * thr / ptjac**2 - thrr / ptjac) * rt &
515  + (ptjs * thr / ptjac**2 - thrs / ptjac) * st
516  xpt = -xrr * ths * zpsis / ptjac**2 + xrs * (ths * zpsir &
517  + thr * zpsis) / ptjac**2 + xr * rpt + xs * spt - xss &
518  * thr * zpsir / ptjac**2
519  ypt = -yrr * ths * zpsis / ptjac**2 + yrs * (ths * zpsir &
520  + thr * zpsis) / ptjac**2 + yr * rpt + ys * spt - yss &
521  * thr * zpsir / ptjac**2
522  xyjac = xr * ys - xs * yr
523  psix = (ys * zpsir - yr * zpsis) / xyjac
524  psiy = (-xs * zpsir + xr * zpsis) / xyjac
525 ! call radial_mesh(radpsi(i), dum, ddum, dddum)
526  etap = 1._r8 / dpsikn(i)
527  ejac = etap * (psix * thy - psiy * thx)
528  xx(1, node) = x
529  yy(1, node) = y
530 !-------------------------------- watch minus sign from r orientation --
531  xx(2, node) = -(thy / ejac) / (2._r8 * (nrnew - 1))
532  yy(2, node) = -(-thx / ejac) / (2._r8 * (nrnew - 1))
533  xx(3, node) = (-etap * psiy / ejac) / (factas * (npnew &
534  - 1) / pi)
535  yy(3, node) = (etap * psix / ejac) / (factas * (npnew &
536  - 1) / pi)
537  xx(4, node) = -(xpt / etap) / (2._r8 * factas * (nrnew - 1) &
538  * (npnew - 1) / pi)
539  yy(4, node) = -(ypt / etap) / (2._r8 * factas * (nrnew - 1) &
540  * (npnew - 1) / pi)
541  psi(4 * (node - 1) + 1) = zpsi
542  psi(4 * (node - 1) + 2) = -dpsikn(i) / (2._r8 * (nrnew - 1.))
543  psi(4 * (node - 1) + 3) = 0._r8
544  psi(4 * (node - 1) + 4) = 0._r8
545  chagr(node) = .true.
546  if (itn > ittest) ittest = itn
547  if (crashed .and. verbosity > 6) &
548  write(iu6,*) "goto 50 after itn = ", itn
549  goto 50
550  else if ((tht1 <= thtval + 1.e-5_r8 .and. tht2 >= thtval - 1.e-5_r8) &
551  .or. (ias == 1 .and. tht1 > tht2 + 1.57_r8 &
552  .and. tht1 <= thtval + 1.e-5_r8 &
553  .and. tht2 + 2._r8 * pi >= thtval - 1.e-5_r8) &
554  .or. (ias == 1 .and. tht1 > tht2 + 1.57_r8 &
555  .and. tht1 - 2._r8 * pi <= thtval + 1.e-5_r8 &
556  .and. tht2 >= thtval - 1.e-5_r8)) then
557 !
558 !----------------------------- theta value bracketed ------------
559 !
560  if (crashed .and. verbosity > 6) write(iu6,*) "Bracketed"
561  fvec1 = -(zy - yaxis) + (zx - xaxis) * tan(thtval)
562  fvec2 = -ps + psival
563  if (crashed .and. verbosity > 6) &
564  write(iu6, *) "fvec1, fvec2 = ", fvec1, fvec2
565  if (crashed .and. verbosity > 6) &
566  write(iu6, *) "zyr, zxr, zys, zxs = ", zyr, zxr, zys, zxs
567  if (crashed .and. verbosity > 6) &
568  write(iu6, *) "thtval, tan(thtval) = ",thtval, tan(thtval)
569  fjac11 = zyr - zxr * tan(thtval)
570  fjac12 = zys - zxs * tan(thtval)
571  fjac21 = dpsir
572  fjac22 = dpsis
573  if (crashed .and. verbosity > 6) &
574  write(iu6, *) "fjac11, fjac12, fjac21, fjac22 = ", &
575  fjac11, fjac12, fjac21, fjac22
576  dis = fjac22 * fjac11 - fjac12 * fjac21
577  dr = (fjac22 * fvec1 - fjac12 * fvec2) / dis
578  ds = (fjac11 * fvec2 - fjac21 * fvec1) / dis
579  dtot = sqrt(dr * dr + ds * ds)
580  if (dis * dis_old < 0._r8) then
581  write(iu6, *) "Sign change in dis: old, new = ", dis_old, dis
582  write(iu6, *) nn
583  trr = rr
584  trr = rr_exclusion_zone
585  trr = 1 - (1 - rr_exclusion_zone) / 2._r8
586  trr = 0.75_r8
587  trr = 0.99_r8
588  trr = 0.999_r8
589  trr = 0.99999_r8
590  trr = 0.999999_r8
591  open(133, file='data.dump')
592  do tjn = 1, np - 1
593  tnode = (in - 1) * (np - 1) + tjn
594  call set_node_number(tnode, ttn1, ttn2, ttn3, ttn4)
595  do tjns = -100, 100
596  tss = real(tjns) / real(100)
597  call interpolation(1, xxold(1, ttn1), xxold(1, ttn2), &
598  xxold(1, ttn3), xxold(1, ttn4), trr, tss, tx, txr, txs)
599  call interpolation(1, yyold(1, ttn1), yyold(1, ttn2), &
600  yyold(1, ttn3), yyold(1, ttn4), trr, tss, ty, tyr, tys)
601  call interpolation(1, psiold(4 * (ttn1 - 1) + 1), &
602  psiold(4 * (ttn2 - 1) + 1), psiold(4 * (ttn3 - 1) + 1), &
603  psiold(4 * (ttn4 - 1) + 1), trr, tss, tpsi, tpsir, tpsis)
604  ttht = atan2(ty - yaxis, tx - xaxis)
605  write(133, *) tjn, tjns, tss, tx, ty, tpsi, ttht
606  end do
607  end do
608  close(133)
609  stop
610  end if
611  dis_old = dis
612  if (crashed .and. verbosity > 6) &
613  write(iu6,*) "dtot, dis = ", dtot, dis
614 ! trying out alternative logic
615  if (.false.) then
616  if (in == nr - 1) then
617  stepr = max(min(step, (1._r8 - rr) / 2._r8), 0._r8)
618  if (verbosity > 6) write(iu6,*) 'Bracketed: Step in r changed ' &
619  // 'to ',stepr
620  ds = ds * step / stepr
621  if (verbosity > 6) write(iu6, *) 'ds rescaled by ', &
622  step / stepr, ' to ', ds
623  if (abs(ds) > 10._r8) then
624  dr = dr / abs(ds)
625  if (verbosity > 6) write(iu6, *) 'dr rescaled by ', &
626  1._r8 / abs(ds),' to ', dr
627  end if
628  end if
629  dr = sign(min(abs(dr), stepr), dr)
630  ds = sign(min(abs(ds), steps), ds)
631  else
632  if (dtot > step) then
633  dr = dr / dtot * step
634  ds = ds / dtot * step
635  end if
636  if (in == nr - 1) then
637  dr = min(dr, (1._r8 - rr) / 2._r8)
638  if (verbosity > 6) write(iu6, *) 'Bracketed: dr changed to ',dr
639  end if
640  end if
641  rr = rr + dr
642  ss = ss + ds
643  if (crashed .and. verbosity > 0) &
644  write(iu6,*) "Incremented rr & ss to ", rr, ss
645  nobrack = .false.
646  else
647 !
648 !------------------------ theta not found yet, track flux surface!
649 
650  if (crashed .and. verbosity > 6) write(iu6, *) "Not bracketed"
651  nobrack = .true.
652  ds_old=ds
653  ds = - dd * dpsir / abs(dpsir)
654  if(itn.gt.1 .and. abs(ds + ds_old).lt.1e-10_r8) then
655  write(*,*) 'seem to have an oscillation in ds'
656  write(*,*) ds_old, ds, nn, rr, ss
657  step = step/2.0_r8
658  write(*,*) 'step changed to ',step
659  endif
660  if (crashed .and. verbosity > 6) write(iu6, *) "ds = ",ds
661  if (abs((psival - ps) / psival) > 0.01_r8) then
662  if (crashed .and. verbosity > 6) write(iu6, *) "moving in ps"
663  ds = 0._r8
664  dr = (psival - ps) / dpsir
665  if (in == nr - 1) then
666  stepr = max(min(step, (1._r8 - rr) / 2._r8), 0._r8)
667  if (verbosity > 6) write(iu6, *) 'Not bracketed (1): Step in ' &
668  // 'r changed to ', stepr
669  end if
670  dr = sign(min(abs(dr), stepr), dr)
671  if (dr == 0._r8) then
672  if (verbosity > 6) then
673  write(iu6, *) 'Only taking a step in r and dr == 0'
674  write(iu6, *) 'We will try to take a step in ds'
675  end if
676  if (psival > psi_bar) then
677  ds = steps / 10._r8
678  else
679  ds = -steps / 10._r8
680  end if
681  ds = sign(min(abs(ds), steps), ds)
682  if (verbosity > 6) write(iu6, *) 'ds = ', ds
683  end if
684  rr = rr + dr
685  ss = ss + ds
686  if (crashed .and. verbosity > 6) &
687  write(iu6, *) "Incremented rr to ", rr
688  else
689  dr = (psival - ps - dpsis * ds) / dpsir
690  dtot = sqrt(dr * dr + ds * ds)
691  dr = dr * dd / dtot
692  ds = ds * dd / dtot
693 ! trying out alternative logic
694  if (.false.) then
695  if (in == nr - 1) then
696  stepr = max(min(step, (1._r8 - rr) / 2._r8), 0._r8)
697  if (verbosity > 6) &
698  write(iu6, *) 'Not bracketed (2): Step in r changed to ', stepr
699  end if
700  dr = sign(min(abs(dr), stepr), dr)
701  ds = sign(min(abs(ds), steps), ds)
702  else
703  dtot = sqrt(dr * dr + ds * ds)
704  if (dtot > step) then
705  dr = dr / dtot * step
706  ds = ds / dtot * step
707  end if
708  if (in == nr - 1) then
709  dr = min(dr, (1._r8 - rr) / 2._r8)
710  if (verbosity > 6) &
711  write(iu6, *) 'Not bracketed (2): dr changed to ', dr
712  end if
713  end if
714  rr = rr + dr
715  ss = ss + ds
716  if (crashed .and. verbosity > 6) &
717  write(iu6, *) "Incremented rr & ss to ", rr, ss
718  end if
719  end if
720 
721  if (in == nr - 1 .and. verbosity > 6) then
722  write(iu6, *) 'Entered danger zone; in, rr, jn, ss, step = ', &
723  in, rr, jn, ss, step
724  write(iu6, *) 'Searching for psi, theta = ', psival, thtval
725  write(iu6, *) 'Current position is ', ps, tht
726  write(iu6, *) 'psi_min, psi_bar, psi_max = ', psi_min, psi_bar, &
727  psi_max
728  end if
729  if (crashed .and. verbosity > 6) &
730  write(iu6, *) 'in, jn = ', in, jn
731  if (ss < -1._r8) then
732  if (jprev == -1) then
733  nn = nn - 1
734  if (crashed .and. verbosity > 6) &
735  write(iu6, *) "Decrementing nn by 1 to ", nn
736  if (jn == 1) then
737  nn = nn + np - 1
738  if (crashed .and. verbosity > 6) &
739  write(iu6, *) "Skipping nn up to ", nn
740  end if
741  ss = ss + 2._r8
742  jprev = 0
743  else
744  ss = -1._r8
745  if (crashed .and. verbosity > 6) &
746  write(iu6, *) "Set ss = -1"
747  jprev = -1
748  end if
749  else if (ss > 1) then
750  if (jprev == 1) then
751  nn = nn + 1
752  if (crashed .and. verbosity > 6) &
753  write(iu6, *) "Incrementing nn by 1 to ", nn
754  if (jn == np - 1) then
755  nn = nn - np + 1
756  if (crashed .and. verbosity > 6) &
757  write(iu6, *) "Skipping nn down to ", nn
758  end if
759  ss = ss - 2._r8
760  jprev = 0
761  else
762  if (crashed .and. verbosity > 6) &
763  write(iu6, *) "Set ss = +1"
764  ss = 1._r8
765  jprev = 1
766  end if
767  end if
768  if (rr < -1._r8) then
769  if (in > 1 .and. iprev == -1) then
770  nn = nn - np + 1
771  if (crashed .and. verbosity > 6) &
772  write(iu6, *) "Decrementing nn by np to ", nn
773  rr = rr + 2._r8
774  iprev = 0
775  else
776  rr = -1._r8
777  if (crashed .and. verbosity > 6) &
778  write(iu6, *) "Set rr = -1"
779  iprev = -1
780  end if
781  else if (rr > 1._r8) then
782  if (in < nr - 1 .and. iprev == 1) then
783  nn = nn + np - 1
784  if (crashed .and. verbosity > 6) &
785  write(iu6, *) "Incrementing nn by np to ", nn
786  rr = rr - 2._r8
787  iprev = 0
788  else
789  rr = 1._r8
790  if (crashed .and. verbosity > 6) &
791  write(iu6, *) "Set rr = +1"
792  iprev = 1
793  end if
794  end if
795  end do
796  if (standard_output) &
797  write(out_he, *) ' fatal node not found : ',i, j, psival, thtval
798  write(iu6, *) ' fatal error: node not found in remesh !'
799  if (crashed) then
800  crashed = .false.
801  nn = nnsave2
802  stop !!!!
803  else
804  crashed = .true.
805  nn = nnsave2
806  goto 47
807  end if
808  50 end do
809  crashed = .false.
810  if (ittest > 2500) then
811  if (verbosity > 0) &
812  write(iu6, *) ' high max number of iterations : ', ittest
813  end if
814  end do
815 
816  if (ias == 1) then
817 !------------------------------------ copy tht = 0 to tht = 2pi
818  do i = 1, nrnew - 1
819  j = npnew
820  node = (i - 1) * npnew + j
821  j0 = 1
822  n0 = (i - 1) * npnew + j0
823  xx(1, node) = xx(1, n0)
824  xx(2, node) = xx(2, n0)
825  xx(3, node) = xx(3, n0)
826  xx(4, node) = xx(4, n0)
827  yy(1, node) = yy(1, n0)
828  yy(2, node) = yy(2, n0)
829  yy(3, node) = yy(3, n0)
830  yy(4, node) = yy(4, n0)
831  psi(4 * (node - 1) + 1) = psi(4 * (n0 - 1) + 1)
832  psi(4 * (node - 1) + 2) = psi(4 * (n0 - 1) + 2)
833  psi(4 * (node - 1) + 3) = psi(4 * (n0 - 1) + 3)
834  psi(4 * (node - 1) + 4) = psi(4 * (n0 - 1) + 4)
835  chagr(node) = .true.
836  end do
837  end if
838 
839  call set_node_number(nax, n1, n2, n3, n4)
840  s = sax
841  r = rax
842  call interpolation(2, xxold(1, n1), xxold(1, n2), xxold(1, n3), &
843  xxold(1, n4), r, s, x, xr, xs, xrs, xrr, xss)
844  call interpolation(2, yyold(1, n1), yyold(1, n2), yyold(1, n3), &
845  yyold(1, n4), r, s, y, yr, ys, yrs, yrr, yss)
846  call interpolation(2, psiold(4 * (n1 - 1) + 1), psiold(4 * (n2 - 1) + 1), &
847  psiold(4 * (n3 - 1) + 1), psiold(4 * (n4 - 1) + 1), r, s, eqpsi, &
848  psir, psis, psrs, psrr, psss)
849  ejac = xr * ys - xs * yr
850  ry = -xs / ejac
851  rx = ys / ejac
852  sy = xr / ejac
853  sx = -yr / ejac
854  psixx = psrr * rx * rx + 2._r8 * psrs * rx * sx + psss * sx * sx
855  psiyy = psrr * ry * ry + 2._r8 * psrs * ry * sy + psss * sy * sy
856  if (hbt) then
857  psiyy0 = a * (1._r8 + b * x * (1._r8 + eps * x / 2._r8)) - psixx
858  else
859  psiyy0 = a * (1._r8 + b * (1._r8 + eps * x)**2) - psixx
860  end if
861  cx = psixx / 2._r8
862  cy = psiyy0 / 2._r8
863  do j = 1, npnew
864  node = (nrnew - 1) * npnew + j
865  xx(1, node) = xaxis
866  yy(1, node) = yaxis
867  tn = tan(thtkn(j))
868  tn2 = tn**2
869  cn = cos(thtkn(j))
870  if (thtkn(j) == pi / 2._r8) then
871  xx(2, node) = 0._r8
872  yy(2, node) = -1._r8 / (sqrt(cy) * 2._r8 * (nrnew - 1))
873  xx(4, node) = 1._r8 / (sqrt(cy) * 2._r8 * factas * (nrnew - 1) &
874  * (npnew - 1) / pi)
875  yy(4, node) = 0._r8
876  elseif (thtkn(j) == (3._r8 * pi / 2._r8)) then
877  xx(2, node) = 0._r8
878  yy(2, node) = -1._r8 / (sqrt(cy) * 2._r8 * (nrnew - 1))
879  xx(4, node) = 1._r8 / (sqrt(cy) * 2._r8 * factas * (nrnew - 1) &
880  * (npnew - 1) / pi)
881  yy(4, node) = 0._r8
882  else
883  xx(2, node) = -sign(1._r8, cn) / (sqrt(cx + cy * tn2) * 2._r8 &
884  * (nrnew - 1))
885  yy(2, node) = -abs(tn) / (sqrt(cx + cy * tn2) * 2._r8 * (nrnew - 1))
886  xx(4, node) = (cx + cy * tn2)**(-1.5_r8) * cy * abs(tn) / (cn**2 &
887  * 2._r8* factas * (nrnew - 1) * (npnew - 1) / pi)
888  yy(4, node) = -cx * (cx + cy * tn2)**(-1.5_r8) / (cn * abs(cn) &
889  * 2._r8 * factas * (nrnew - 1) * (npnew - 1) / pi)
890  end if
891  if (thtkn(j) > pi) then
892  yy(2, node) = -yy(2, node)
893  xx(4, node) = -xx(4, node)
894  end if
895  xx(3, node) = 0._r8
896  yy(3, node) = 0._r8
897  psi(4 * (node - 1) + 1) = 0._r8
898  psi(4 * (node - 1) + 2) = 0._r8
899  psi(4 * (node - 1) + 3) = 0._r8
900  psi(4 * (node - 1) + 4) = 0._r8
901  chagr(node) = .true.
902  end do
903  do i = 1, nrnew * npnew
904  if (standard_output) then
905  if (.not. chagr(i)) write(out_he, *) 'node missed at i = ', i
906  end if
907  end do
908 
909 !-- re-calculate node numbers
910  if (nr /= nrnew .or. np /= npnew) then
911  deallocate(nodeno)
912  nr = nrnew
913  np = npnew
914  allocate(nodeno((nr - 1) * (np - 1), 4))
915  call node_numbers
916  end if
917 
918 ! 3 format('elongation on axis : cx, cy = ', 2e12.4)
919 ! 41 format(i3,2e12.4,4e14.6)
920 
921  return
922 
923 end subroutine remesh
mod_dat
Definition: mod_dat.f90:1
radial_mesh
subroutine radial_mesh(rpsi, zpsi, dzpsi, ddzpsi)
Definition: radial_mesh.f90:1
mod_mesh
Definition: mod_mesh.f90:1
mod_corners
Definition: mod_corners.f90:1
mod_interpolation::interpolation
subroutine interpolation(type_interp, xn1, xn2, xn3, xn4, r, s, x, xr, xs, xrs, xrr, xss, yn1, yn2, yn3, yn4, pn1, pn2, pn3, pn4, yr, ys, ps)
Definition: mod_interpolation.f90:10
matrix_calculations
Definition: matrix_calculations.f90:1
mod_profiles::allocate_spline_coefficients
subroutine allocate_spline_coefficients(spline, n)
Definition: mod_profiles.f90:27
mod_output
Definition: mod_output.f90:1
mod_interpolation
Definition: mod_interpolation.f90:1
remesh
subroutine remesh(a, nrnew, npnew, meshno, cx, cy, xaxis, yaxis, nax, rax, sax)
Definition: remesh.f90:1
spline
subroutine spline(N, X, Y, ALFA, BETA, TYP, A, B, C, D)
Definition: solution6.f90:655
spwert
REAL *8 function spwert(N, XWERT, A, B, C, D, X, ABLTG)
Definition: solution6.f90:155
mod_profiles
Definition: mod_profiles.f90:1
mod_nodes
Definition: mod_nodes.f90:1
math_functions
Definition: math_functions.f90:1
mod_parameter
Definition: mod_parameter.f90:1
math_functions::quad_equation
subroutine, public quad_equation(p_m, p_mr, p_p, p_pr, tol, r, ifail)
Definition: math_functions.f90:489
solve_cubic
subroutine solve_cubic(c0, c1, c2, c3, x1, x2, x3, ifail)
Definition: solve_cubic.f90:1
mod_profiles::deallocate_spline_coefficients
subroutine deallocate_spline_coefficients(spline)
Definition: mod_profiles.f90:43
mesh_accumulation
subroutine mesh_accumulation(nrtmp, nv, s_acc, sig, weights, equidistant, sg, dsg, ddsg, in_psi)
Definition: mesh_accumulation.f90:1
mod_meshacc
Definition: mod_meshacc.f90:1
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
node_numbers
subroutine node_numbers
Definition: node_numbers.f90:1
matrix_calculations::set_node_number
subroutine set_node_number(n_node, n1, n2, n3, n4)
Definition: matrix_calculations.f90:11
mod_helena_io
Definition: mod_helena_io.f90:1
helena_spline::spline_coefficients
Definition: helena_spline.f90:7
helena_spline
Definition: helena_spline.f90:1