V3FIT
cyl_flux.f
1  MODULE cyl_flux
2  USE stel_kinds
3  USE stel_constants, ONLY: twopi, one, zero
4  IMPLICIT NONE
5 
6 !-------------------------------------------------------------------------------
7 ! Module for coordinate conversion, cylindrical <- -> VMEC flux
8 ! CONTAINS two subroutines:
9 ! flx2cyl - convert from VMEC flux to cylindrical - inverse Fourier Transform
10 ! cyl2flx - convert from cylindrical to VMEC flux - Root find using flx2cyl
11 !
12 ! cyl2flx CONTAINS two subroutines:
13 ! newt2d - implements a 2 dimensional Newton rootfind
14 ! get_flxcoord - interface with fewer arguments to flx2cyl
15 ! Note - since newt2d and get_flxcoord are CONTAINed in cyl2flx, they have
16 ! access to all the variables in cyl2flx.
17 !
18 ! 2011-09-06 JDH
19 ! Split this module from vmec_utils, to clarify interfaces and dependencies
20 !-------------------------------------------------------------------------------
21 
22  CONTAINS
23 
24 !*******************************************************************************
25 !-------------------------------------------------------------------------------
26 ! Convert from flux to cylindrical cooridnates
27 !-------------------------------------------------------------------------------
28  SUBROUTINE flx2cyl(rzl_array, c_flux, r_cyl, ns, ntor,
29  1 mpol, ntmax, lthreed, lasym, iflag,
30  2 mscale, nscale, Ru, Rv, Zu, Zv)
31  IMPLICIT NONE
32 C-----------------------------------------------
33 C D u m m y A r g u m e n t s
34 C-----------------------------------------------
35  INTEGER, INTENT(out) :: iflag
36  INTEGER, INTENT(in) :: ns, ntmax, mpol, ntor
37  LOGICAL :: lthreed, lasym
38  REAL(rprec), DIMENSION(ns,0:ntor,0:mpol-1,3*ntmax),
39  1 INTENT(in) :: rzl_array
40  REAL(rprec), INTENT(in) :: c_flux(3)
41  REAL(rprec), INTENT(out) :: r_cyl(3)
42  REAL(rprec), INTENT(in), OPTIONAL ::
43  1 mscale(0:mpol-1), nscale(0:ntor)
44  REAL(rprec), INTENT(out), OPTIONAL :: Ru, Rv, Zu, Zv
45 C-----------------------------------------------
46 C L o c a l P a r a m e t e r s
47 C-----------------------------------------------
48  INTEGER, PARAMETER :: rcc = 1
49 C-----------------------------------------------
50 C L o c a l V a r i a b l e s
51 C-----------------------------------------------
52  INTEGER :: rss, rsc, rcs, zsc, zcs, zcc, zss
53  INTEGER :: istat, jslo, jshi, mpol1, m, n
54  REAL(rprec), DIMENSION(0:ntor,0:mpol-1) ::
55  1 rmncc, rmnss, zmncs, zmnsc,
56  2 rmncs, rmnsc, zmncc, zmnss
57  REAL(rprec) :: wlo, whi, wlo_odd, whi_odd, hs1,
58  1 si, ui, vi, r11, z11
59  REAL(rprec) :: wmins(0:ntor,0:mpol-1),
60  1 wplus(0:ntor,0:mpol-1)
61  REAL(rprec) :: cosu, sinu, cosv, sinv,
62  1 cosmu(0:mpol-1), sinmu(0:mpol-1),
63  2 cosnv(0:ntor), sinnv(0:ntor),
64  3 cosnvn(0:ntor), sinnvn(0:ntor)
65  REAL(rprec) :: work1(0:mpol-1,12)
66  LOGICAL :: lru, lrv, lzu, lzv
67 !-------------------------------------------------------------------------------!
68 ! COMPUTES THE CYLINDRICAL COORDINATES R11 and Z11
69 ! AT A FIXED FLUX COORDINATE POINT si, ui(theta), vi(N*phi)
70 ! WHERE phi = geometric toroidal angle (0 to 2pi), N = no. field periods
71 !
72 ! INPUT:
73 ! c_flux: array of (si,ui,vi) values to convert to cylindrical coordinates
74 ! rzl_array: array of (r, z, lambda) Fourier coefficients for all radial, m, n values
75 ! ns, mpol,ntor, ntmax: radial, poloidal, toroidal, type (r,z,l) dimensions of rzl_array
76 ! mscale, nscale: option scale factors. Use ONLY if rzl_array comes from within VMEC.
77 ! If arising from WOUT file, mscale, nscale => 1 are not passed
78 !
79 ! OUTPUT:
80 ! r_cyl : R = r_cyl(1); N*PHI = r_cyl(2); Z = r_cyl(3)
81 !
82 ! OPTIONAL OUTPUT
83 ! Ru = dR/du; Rv = dR/dv = dR/dphi / N
84 ! Zu = dZ/du; Zv = dZ/dv = dZ/dphi / N
85 !
86 ! NOTE: User is responsible for multiplying Rv, Zv by N to get phi derivatives
87 !
88 !-------------------------------------------------------------------------------
89 
90  iflag = -1
91  si = c_flux(1); ui = c_flux(2); vi = c_flux(3)
92  r_cyl(2) = vi
93 
94 ! IF (si.lt.zero .or. si.gt.one) THEN
95  IF (si .lt. zero) THEN
96  WRITE(6, *)' In flx2cyl, s(flux) must be > 0'
97  RETURN
98  END IF
99 
100  lru = PRESENT(ru); lrv = PRESENT(rv)
101  lzu = PRESENT(zu); lzv = PRESENT(zv)
102  IF (lrv .and. .not. lthreed) rv = 0
103  IF (lzv .and. .not. lthreed) zv = 0
104 
105 !
106 ! FIND INTERPOLATED s VALUE AND COMPUTE INTERPOLATION WEIGHTS wlo, whi (for even m modes)
107 ! AND wlo_odd, whi_odd (for odd m modes).
108 ! FOR si <= 1, POINT IS INSIDE PLASMA;
109 ! FOR si > 1, TRY EXTRAPOLATION (WITH CONTINUOUS s, u DERIVATIVES INTO "vacuum" REGION
110 !
111  hs1 = one/(ns-1)
112  IF (si .le. one) THEN
113 
114  jslo = 1 + si/hs1
115  jshi = jslo+1
116  wlo = (hs1*(jshi-1) - si)/hs1
117  whi = 1 - wlo
118  IF (jslo .eq. ns) jshi = jslo
119 
120 !
121 ! USE Rmn, Zmn ~ SQRT(s) FOR ODD-m MODES, SO INTERPOLATE Xmn/SQRT(s)
122 !
123  whi_odd = whi*sqrt(si/(hs1*(jshi-1)))
124  IF (jslo .ne. 1) THEN
125  wlo_odd = wlo*sqrt(si/(hs1*(jslo-1)))
126  ELSE
127  wlo_odd = 0
128  whi_odd = sqrt(si/(hs1*(jshi-1)))
129  END IF
130 
131  ELSE
132 
133  jshi = ns
134  jslo = ns-1
135  wlo = -(si - 1)/hs1; wlo_odd = wlo
136  whi = 1 - wlo; whi_odd = whi
137 
138  ENDIF
139 
140 
141  mpol1 = mpol-1
142 
143  wmins(:,0:mpol1:2) = wlo
144  wplus(:,0:mpol1:2) = whi
145  wmins(:,1:mpol1:2) = wlo_odd
146  wplus(:,1:mpol1:2) = whi_odd
147 
148 
149 
150  IF (.not.lasym) THEN
151  IF (lthreed) THEN
152  IF (ntmax .ne. 2) stop 'ntmax != 2 in flx2cyl!'
153  rss = 2; zcs = 2
154  ELSE
155  IF (ntmax .ne. 1) stop 'ntmax != 1 in flx2cyl!'
156  END IF
157  ELSE
158  IF (lthreed) THEN
159  IF (ntmax .ne. 4) stop 'ntmax != 4 in flx2cyl!'
160  rss = 2; rsc = 3; rcs = 4
161  zcs = 2; zcc = 3; zss = 4
162  ELSE
163  IF (ntmax .ne. 2) stop 'ntmax != 2 in flx2cyl!'
164  rsc = 2; zcc = 2
165  END IF
166  END IF
167 
168  zsc = 1+ntmax; zcs = zcs+ntmax; zcc = zcc+ntmax; zss = zss+ntmax
169 
170  rmncc = wmins*rzl_array(jslo,:,:,rcc)
171  1 + wplus*rzl_array(jshi,:,:,rcc) !!COS(mu) COS(nv)
172  zmnsc = wmins*rzl_array(jslo,:,:,zsc)
173  1 + wplus*rzl_array(jshi,:,:,zsc) !!SIN(mu) COS(nv)
174 
175  IF (lthreed) THEN
176  rmnss = wmins*rzl_array(jslo,:,:,rss)
177  1 + wplus*rzl_array(jshi,:,:,rss) !!SIN(mu) SIN(nv)
178  zmncs = wmins*rzl_array(jslo,:,:,zcs)
179  1 + wplus*rzl_array(jshi,:,:,zcs) !!COS(mu) SIN(nv)
180  END IF
181 
182 !
183 ! SETUP TRIG ARRAYS
184 !
185  cosu = cos(ui); sinu = sin(ui)
186  cosv = cos(vi); sinv = sin(vi)
187 
188  cosmu(0) = 1; sinmu(0) = 0
189  cosnv(0) = 1; sinnv(0) = 0
190  DO m = 1, mpol1
191  cosmu(m) = cosmu(m-1)*cosu - sinmu(m-1)*sinu
192  sinmu(m) = sinmu(m-1)*cosu + cosmu(m-1)*sinu
193  END DO
194 
195  IF (PRESENT(mscale)) THEN
196  cosmu = cosmu*mscale; sinmu = sinmu*mscale
197  END IF
198 
199  DO n = 1, ntor
200  cosnv(n) = cosnv(n-1)*cosv - sinnv(n-1)*sinv
201  sinnv(n) = sinnv(n-1)*cosv + cosnv(n-1)*sinv
202  END DO
203 
204  IF (PRESENT(nscale)) THEN
205  cosnv = cosnv*nscale; sinnv = sinnv*nscale
206  END IF
207 
208  IF (lrv .or. lzv) THEN
209  DO n = 0, ntor
210  cosnvn(n) = n*cosnv(n)
211  sinnvn(n) =-n*sinnv(n)
212  END DO
213  END IF
214 
215  iflag = 0
216 
217 !
218 ! COMPUTE R11, Z11 IN REAL SPACE
219 !
220 ! FIRST, INVERSE TRANSFORM IN N-V SPACE, FOR FIXED M
221 !
222  DO m = 0, mpol1
223 
224  work1(m,1) = sum(rmncc(:,m)*cosnv(:))
225  work1(m,2) = sum(zmnsc(:,m)*cosnv(:))
226  IF (lru) work1(m,3) =-m*work1(m,1)
227  IF (lzu) work1(m,4) = m*work1(m,2)
228  IF (lthreed) THEN
229  IF (lrv) work1(m,5) = sum(rmncc(:,m)*sinnvn(:))
230  IF (lzv) work1(m,6) = sum(zmnsc(:,m)*sinnvn(:))
231  work1(m,7) = sum(rmnss(:,m)*sinnv(:))
232  work1(m,8) = sum(zmncs(:,m)*sinnv(:))
233  IF (lru) work1(m,9) = m*work1(m,7)
234  IF (lzu) work1(m,10) =-m*work1(m,8)
235  IF (lrv) work1(m,11) = sum(rmnss(:,m)*cosnvn(:))
236  IF (lzv) work1(m,12) = sum(zmncs(:,m)*cosnvn(:))
237  END IF
238 
239  END DO
240 
241 !
242 ! NEXT, INVERSE TRANSFORM IN M-U SPACE
243 !
244  IF (lthreed) THEN
245  r11 = sum(work1(:,1)*cosmu(:) + work1(:,7)*sinmu(:))
246  z11 = sum(work1(:,2)*sinmu(:) + work1(:,8)*cosmu(:))
247  IF (lru) ru = sum(work1(:,3)*sinmu(:) + work1(:,9)*cosmu(:))
248  IF (lzu) zu = sum(work1(:,4)*cosmu(:) + work1(:,10)*sinmu(:))
249  IF (lrv) rv = sum(work1(:,5)*cosmu(:) + work1(:,11)*sinmu(:))
250  IF (lzv) zv = sum(work1(:,6)*sinmu(:) + work1(:,12)*cosmu(:))
251  ELSE
252  r11 = sum(work1(:,1)*cosmu(:))
253  z11 = sum(work1(:,2)*sinmu(:))
254  IF (lru) ru = sum(work1(:,3)*sinmu(:))
255  IF (lzu) zu = sum(work1(:,4)*cosmu(:))
256  END IF
257 
258 
259  IF (.not.lasym) GOTO 1000
260 
261  rmnsc = wmins*rzl_array(jslo,:,:,rsc)
262  1 + wplus*rzl_array(jshi,:,:,rsc) !!SIN(mu) COS(nv)
263  zmncc = wmins*rzl_array(jslo,:,:,zcc)
264  1 + wplus*rzl_array(jshi,:,:,zcc) !!COS(mu) COS(nv)
265 
266  IF (lthreed) THEN
267  rmncs = wmins*rzl_array(jslo,:,:,rcs)
268  1 + wplus*rzl_array(jshi,:,:,rcs) !!COS(mu) SIN(nv)
269  zmnss = wmins*rzl_array(jslo,:,:,zss)
270  1 + wplus*rzl_array(jshi,:,:,zss) !!SIN(mu) SIN(nv)
271  END IF
272 
273 !
274 ! COMPUTE R11, Z11 IN REAL SPACE
275 !
276 ! FIRST, INVERSE TRANSFORM IN N-V SPACE, FOR FIXED M
277 !
278  DO m = 0, mpol1
279 
280  work1(m,1) = sum(rmnsc(:,m)*cosnv(:))
281  work1(m,2) = sum(zmncc(:,m)*cosnv(:))
282  IF (lru) work1(m,3) = m*work1(m,1)
283  IF (lzu) work1(m,4) =-m*work1(m,2)
284 
285  IF (lthreed) THEN
286  IF (lrv) work1(m,5) = sum(rmnsc(:,m)*sinnvn(:))
287  IF (lzv) work1(m,6) = sum(zmncc(:,m)*sinnvn(:))
288  work1(m,7) = sum(rmncs(:,m)*sinnv(:))
289  work1(m,8) = sum(zmnss(:,m)*sinnv(:))
290  IF (lru) work1(m,9) =-m*work1(m,7)
291  IF (lzu) work1(m,10) = m*work1(m,8)
292  IF (lrv) work1(m,11) = sum(rmncs(:,m)*cosnvn(:))
293  IF (lzv) work1(m,12) = sum(zmnss(:,m)*cosnvn(:))
294  END IF
295 
296  END DO
297 
298 !
299 ! NEXT, INVERSE TRANSFORM IN M-U SPACE
300 !
301  IF (lthreed) THEN
302  r11 = r11 + sum(work1(:,1)*sinmu(:) + work1(:,7)*cosmu(:))
303  z11 = z11 + sum(work1(:,2)*cosmu(:) + work1(:,8)*sinmu(:))
304  IF (lru) ru = ru +
305  1 sum(work1(:,3)*cosmu(:) + work1(:,9)*sinmu(:))
306  IF (lzu) zu = zu +
307  1 sum(work1(:,4)*sinmu(:) + work1(:,10)*cosmu(:))
308  IF (lrv) rv = rv +
309  1 sum(work1(:,5)*sinmu(:) + work1(:,11)*cosmu(:))
310  IF (lzv) zv = zv +
311  1 sum(work1(:,6)*cosmu(:) + work1(:,12)*sinmu(:))
312  ELSE
313  r11 = r11 + sum(work1(:,1)*sinmu(:))
314  z11 = z11 + sum(work1(:,2)*cosmu(:))
315  IF (lru) ru = ru + sum(work1(:,3)*cosmu(:))
316  IF (lzu) zu = zu + sum(work1(:,4)*sinmu(:))
317  END IF
318 
319  1000 CONTINUE
320 
321  r_cyl(1) = r11; r_cyl(3) = z11
322 
323  END SUBROUTINE flx2cyl
324 
325 !*******************************************************************************
326 !-------------------------------------------------------------------------------
327 ! convert from cylindrical to flux coordinates
328 !-------------------------------------------------------------------------------
329  SUBROUTINE cyl2flx(rzl_in, r_cyl, c_flx, ns_in, ntor_in, mpol_in,
330  1 ntmax_in, lthreed_in, lasym_in, info, nfe, fmin,
331  1 mscale, nscale, ru, zu, rv, zv)
332  IMPLICIT NONE
333 C-----------------------------------------------
334 C D u m m y A r g u m e n t s
335 C-----------------------------------------------
336  INTEGER, INTENT(out) :: info, nfe
337  INTEGER, INTENT(in) :: ns_in, ntor_in, mpol_in, ntmax_in
338  REAL(rprec), INTENT(in) :: r_cyl(3)
339  REAL(rprec), INTENT(inout) :: c_flx(3)
340  REAL(rprec), INTENT(in), TARGET ::
341  1 rzl_in(ns_in,0:ntor_in,0:mpol_in-1,3*ntmax_in)
342  REAL(rprec), TARGET, OPTIONAL ::
343  1 mscale(0:mpol_in-1), nscale(0:ntor_in)
344  REAL(rprec), INTENT(out), OPTIONAL :: ru, zu, rv, zv
345  REAL(rprec), INTENT(out) :: fmin
346  LOGICAL, INTENT(in) :: lthreed_in, lasym_in
347 C-----------------------------------------------
348 C L o c a l P a r a m e t e r s
349 C-----------------------------------------------
350  INTEGER, PARAMETER :: nvar = 2
351  real(rprec), PARAMETER :: ftol = 1.e-16_dp
352 C-----------------------------------------------
353 C L o c a l V a r i a b l e s
354 C-----------------------------------------------
355  REAL(rprec) :: xc_opt(nvar), r_cyl_out(3), fmin0
356  INTEGER :: iflag, itry, nfe_out
357 
358 C-----------------------------------------------
359 C Variables to communicate with internal subroutines newt2d and get_flxcoord
360 C-----------------------------------------------
361  INTEGER :: ns_loc, ntmax_loc, mpol_loc, ntor_loc
362  REAL(rprec) :: r_target, phi_target, z_target, fnorm
363  REAL(rprec), POINTER :: rzl_array(:,:,:,:)
364  REAL(rprec), POINTER :: mscale_loc(:), nscale_loc(:)
365  LOGICAL :: lthreed_loc, lasym_loc, lscale
366 
367 C-----------------------------------------------
368 ! LOCAL PARAMETERS:
369 ! ftol : nominally, set to 1.E-16. Gives a maximum (relative)
370 ! error in matching R and Z of sqrt(ftol), or 1.E-8.
371 ! To increase accuracy, ftol should be lowered, but this
372 ! may require more Newton iterations (slows code).
373 !
374 ! INPUT:
375 ! rzl_in : 4D array with r,z (lambda) Fourier coefficients vs. radius
376 !
377 ! r_cyl : vector specifying cylindrical point to match, R = r_cyl(1),
378 ! N*phi = r_cyl(2), Z = r_cyl(3)
379 ! NOTE: N*phi (N=no. field periods) is input, NOT phi!
380 ! ns_in : number of radial nodes in input array rzl_in
381 ! ntmax_in: number of different types of modes (cos-cos, sin-sin, etc)
382 ! mpol_in : number of poloidal modes (0:mpol_in-1)
383 ! ntor_in : number of toroidal modes = ntor_in+1 (0:ntor)
384 ! lthreed_in :true if this is a 3D plasma
385 ! lasym_in: true if this is an asymmetric plasma
386 ! mscale (nscale) :
387 ! optional scaling arrays for cos, sin arrays. Used
388 ! only if this routine is called from within VMEC.
389 !
390 ! OUTPUT:
391 ! nfe : number of function evaluations
392 ! info : = 0, no errors, -1, fmin > ftol (tolerance exceeded on output)
393 ! = -3, s > 1 outside plasma, probably
394 ! fmin : minimum value of f = (r - Rin)**2 + (z - Zin)**2 at c_flx
395 
396 ! INPUT/OUTPUT:
397 ! c_flx : array of flux coordinates (s = c_flx(1), u=theta= c_flx(2),
398 ! v = N*phi= c_flx(3))
399 ! on input, initial guess (recommend magnetic axis if "cold" start)
400 ! on output, s, u values corresponding to r_cyl
401 C-----------------------------------------------
402 ! Initialize global variables
403  rzl_array => rzl_in
404  lthreed_loc = lthreed_in; lasym_loc = lasym_in
405  mpol_loc = mpol_in; ntor_loc = ntor_in
406  ns_loc = ns_in; ntmax_loc = ntmax_in
407  lscale = PRESENT(mscale)
408  IF (lscale) THEN
409  mscale_loc => mscale; nscale_loc => nscale
410  END IF
411  r_target = r_cyl(1); phi_target = r_cyl(2); z_target = r_cyl(3)
412 
413 ! Initialize local variables
414  xc_opt(1) = c_flx(1); xc_opt(2) = c_flx(2)
415 
416 ! Avoid exact magnetic axis, which is singular point
417  IF (c_flx(1) .eq. zero) xc_opt(1) = one/(ns_loc-1)
418 
419  fnorm = r_target**2 + z_target**2
420  IF (fnorm .lt. epsilon(fnorm)) fnorm = 1
421  fnorm = one/fnorm
422 
423 
424  nfe = 0
425  fmin0 = 1
426 
427  DO itry = 1, 4
428 
429  CALL newt2d(xc_opt, fmin, ftol, nfe_out, nvar, info)
430  nfe = nfe + nfe_out
431 
432  IF (fmin.le.ftol .or. info.eq.-3) EXIT
433 !
434 ! JOG POINT (BY ROTATING ANGLE) TO IMPROVE CONVERGENCE
435 !
436  IF (fmin .gt. 1.e-3*fmin0) THEN
437  xc_opt(2) = xc_opt(2) + twopi/20
438  ELSE
439  xc_opt(2) = xc_opt(2) + twopi/40
440  END IF
441 
442  fmin0 = min(fmin, fmin0)
443 ! PRINT *,' ITRY = ', itry+1,' FMIN = ', fmin
444 
445  END DO
446 
447  c_flx(1) = xc_opt(1); c_flx(2) = xc_opt(2); c_flx(3) = phi_target
448 !SPH IF (info.eq.0 .and. c_flx(1).gt.one) c_flx(1) = one
449 
450  c_flx(2) = mod(c_flx(2), twopi)
451  DO WHILE (c_flx(2) .lt. zero)
452  c_flx(2) = c_flx(2) + twopi
453  END DO
454 
455 !
456 ! COMPUTE Ru, Zu, Rv, Zv IF REQUIRED
457 !
458 ! IF ((PRESENT(ru) .or. PRESENT(zu) .or.
459 ! 1 PRESENT(rv) .or. PRESENT(zv)) .and. info.eq.0)
460 ! 2 CALL flx2cyl(rzl_in, c_flx, r_cyl_out, ns_loc, ntor_loc,
461 ! 3 mpol_loc, ntmax_loc, lthreed_loc, lasym_loc,
462 ! 4 iflag, MSCALE=mscale, NSCALE=nscale,
463 ! 5 RU=ru, ZU=zu, RV=rv, ZV=zv)
464 
465  IF (PRESENT(ru) .or. PRESENT(zu) .or. PRESENT(rv) .or. &
466  & PRESENT(zv)) THEN
467  IF (info .eq. 0) THEN
468  CALL flx2cyl(rzl_in, c_flx, r_cyl_out, ns_loc, ntor_loc, &
469  & mpol_loc, ntmax_loc, lthreed_loc, lasym_loc, &
470  & iflag, mscale=mscale, nscale=nscale, &
471  & ru=ru, zu=zu, rv=rv, zv=zv)
472  ELSE ! insure that optional arguments are assigned. JDH 2014-03-11
473  IF (PRESENT(ru)) ru = 0
474  IF (PRESENT(zu)) zu = 0
475  IF (PRESENT(rv)) rv = 0
476  IF (PRESENT(zv)) zv = 0
477  ENDIF
478  ENDIF
479 
480  CONTAINS ! internal subprograms newt2d and get_flxcoord
481 
482 !-------------------------------------------------------------------------------
483 !-------------------------------------------------------------------------------
484  SUBROUTINE newt2d(xc_opt, fmin, ftol, nfe, nvar, iflag)
485  IMPLICIT NONE
486 C-----------------------------------------------
487 C D u m m y A r g u m e n t s
488 C-----------------------------------------------
489  INTEGER, INTENT(in) :: nvar
490  INTEGER, INTENT(out) :: nfe, iflag
491  REAL(rprec), INTENT(inout) :: xc_opt(nvar)
492  REAL(rprec), INTENT(in) :: ftol
493  REAL(rprec), INTENT(out) :: fmin
494 C-----------------------------------------------
495 C L o c a l V a r i a b l e s
496 C-----------------------------------------------
497  INTEGER, PARAMETER :: niter = 50
498  integer :: ieval, isgt1
499  REAL(rprec) :: c_flx(3), r_cyl_out(3), fvec(nvar), sflux,
500  1 uflux, eps0, eps, epu, xc_min(2), factor
501  REAL(rprec) :: x0(3), xs(3), xu(3), dels, delu, tau, fmin0,
502  1 ru1, zu1, edge_value, snew
503 !-------------------------------------------------------------------------------
504 ! INPUT/OUTPUT:
505 ! xc_opt: sflux = xc_opt(1), uflux = xc_opt(2) are the toroidal flux
506 ! coordinate and poloidal angle coordinate, respectively
507 ! iflag: = 0, successfully find s,u point
508 ! =-1, did not converge
509 ! =-3, sflux > 1, probably
510 !
511 ! LOCAL VARIABLES:
512 ! tau: d(R,Z)/d(s,u) (Jacobian)
513 ! isgt1: counter for number of times s>1
514 
515 ! FIND FLUX COORDINATES (s,u) WHICH CORRESPOND TO ZERO OF TARGET FUNCTION
516 !
517 ! F == (R - R_TARGET)**2 + (Z - Z_TARGET)**2
518 !
519 ! FOR A GIVEN CYLINDRICAL COORDINATE POINT (R_TARGET, N*PHI=PHI_TARGET, Z_TARGET)
520 !
521 ! Reference: S.E.Attenberger, W.A.Houlberg, S.P.Hirshman, J. Comp Phys 72 (1987) 435.
522 !
523 ! The algorithm used here modifies this slightly to improve "faltering" convergence
524 ! by choosing a steepest-descent path when the step size has been decreased sufficiently
525 ! without yielding a lower value of F.
526 !
527 !-------------------------------------------------------------------------------
528 
529  iflag = -1
530  eps0 = sqrt(epsilon(eps))
531  xc_min = xc_opt
532 
533  c_flx(3) = phi_target
534  fmin0 = 1.e10_dp
535  factor = 1
536  nfe = 0
537  edge_value = one + one/(ns_loc-1)
538  isgt1 = 0
539 
540  DO ieval = 1, niter
541  nfe = nfe + 1
542 
543  sflux = max(xc_opt(1), zero)
544 ! sflux = MIN(MAX(xc_opt(1), zero), one)
545  uflux = xc_opt(2)
546  c_flx(1) = sflux; c_flx(2) = uflux
547 
548 ! COMPUTE R,Z, Ru, Zu
549  CALL get_flxcoord(x0, c_flx, ru=ru1, zu=zu1)
550  xu(1) = ru1; xu(3) = zu1
551 
552 ! MAKE SURE sflux IS LARGE ENOUGH
553 ! TO COMPUTE d(sqrt(s))/ds ACCURATELY NEAR ORIGIN
554  IF (sflux .ge. 1000*eps0) THEN
555  eps = eps0
556  ELSE
557  eps = eps0*sflux
558  END IF
559 
560 ! COMPUTE Rs, Zs NUMERICALLY
561  eps = abs(eps)
562  IF (sflux .ge. 1-eps) eps = -eps
563  c_flx(1) = sflux + eps
564  CALL get_flxcoord(r_cyl_out, c_flx)
565  xs = (r_cyl_out - x0)/eps
566  c_flx(1) = sflux
567 
568  x0(1) = x0(1) - r_target
569  x0(3) = x0(3) - z_target
570  fmin = (x0(1)**2 + x0(3)**2)*fnorm
571 
572  IF (fmin .gt. fmin0) THEN
573  factor = (2*factor)/3
574  xc_opt = xc_min
575 ! REDIRECT ALONG STEEPEST-DESCENT PATH
576  IF (6*factor .lt. one) THEN
577  dels =-(x0(1)*xs(1) + x0(3)*xs(3))/(xs(1)**2 + xs(3)**2)
578  delu =-(x0(1)*xu(1) + x0(3)*xu(3))/(xu(1)**2 + xu(3)**2)
579  END IF
580  ELSE
581  fmin0 = fmin
582  factor = 1
583  xc_min = xc_opt
584 
585 ! NEWTON STEP
586  tau = xu(1)*xs(3) - xu(3)*xs(1)
587  IF (abs(tau) .le. abs(eps)*r_target**2) THEN
588  iflag = -2
589  EXIT
590  END IF
591  dels = ( x0(1)*xu(3) - x0(3)*xu(1))/tau
592  delu = (-x0(1)*xs(3) + x0(3)*xs(1))/tau
593  IF (fmin .gt. 1.e-3_dp) THEN
594  dels = dels/2; delu = delu/2
595  END IF
596 
597  END IF
598 
599  IF (fmin .le. ftol) EXIT
600 
601  IF (abs(dels) .gt. one) dels = sign(one, dels)
602 ! IF (ABS(delu) .gt. twopi/2) delu = SIGN(twopi/2, delu)
603 
604  snew = xc_opt(1) + dels*factor
605  IF (snew .lt. zero) THEN
606  xc_opt(1) = -snew/2 !Prevents oscillations around origin s=0
607  xc_opt(2) = xc_opt(2) + twopi/2
608  delu = -delu
609 ! factor = (-snew/2-xc_opt(1))/dels
610 ! xc_opt(1) = -snew/2
611  ELSE
612  xc_opt(1) = snew
613  END IF
614  xc_opt(2) = xc_opt(2) + delu*factor
615 
616  IF (xc_opt(1) .gt. edge_value) THEN
617  isgt1 = isgt1+1
618  IF (xc_opt(1) .gt. 2._dp) isgt1 = isgt1+1
619  IF (isgt1 .gt. 5) EXIT
620  END IF
621 
622  END DO
623 
624  IF (isgt1.gt.5) THEN
625  iflag = -3
626  xc_min = xc_opt
627  ELSE IF (xc_min(1) .gt. edge_value) THEN
628  iflag = -3
629  ELSE IF (fmin0 .le. ftol) THEN
630  iflag = 0
631  ELSE
632  iflag = -1
633  END IF
634 
635  fmin = fmin0
636  xc_opt = xc_min
637  xc_opt(2) = mod(xc_opt(2), twopi)
638 
639  END SUBROUTINE newt2d ! Internal subroutine to cyl2flx
640 
641 !-------------------------------------------------------------------------------
642 !-------------------------------------------------------------------------------
643  SUBROUTINE get_flxcoord(x1, c_flx, ru, zu) ! Internal subroutine to cyl2flx
644 C-----------------------------------------------
645 C D u m m y A r g u m e n t s
646 C-----------------------------------------------
647  REAL(rprec), INTENT(out) :: x1(3)
648  REAL(rprec), INTENT(in) :: c_flx(3)
649  REAL(rprec), INTENT(out), OPTIONAL :: ru, zu
650 C-----------------------------------------------
651 C L o c a l V a r i a b l e s
652 C-----------------------------------------------
653  INTEGER :: iflag
654 C-----------------------------------------------
655  IF (lscale) THEN
656  CALL flx2cyl(rzl_array, c_flx, x1, ns_loc, ntor_loc, mpol_loc,
657  1 ntmax_loc, lthreed_loc, lasym_loc, iflag,
658  2 mscale=mscale_loc, nscale=nscale_loc, ru=ru, zu=zu)
659  ELSE
660  CALL flx2cyl(rzl_array, c_flx, x1, ns_loc, ntor_loc, mpol_loc,
661  1 ntmax_loc, lthreed_loc, lasym_loc, iflag,
662  2 ru=ru, zu=zu)
663  END IF
664 
665  END SUBROUTINE get_flxcoord ! Internal subroutine to cyl2flx
666 
667  END SUBROUTINE cyl2flx
668 
669  END MODULE cyl_flux
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17