V3FIT
guess_axis.f
1  SUBROUTINE guess_axis_par(r1, z1, ru0, zu0, lscreen)
2  USE vmec_main
3  USE vmec_params, ONLY: nscale, signgs
4  USE realspace, ONLY: psqrts
5  USE parallel_include_module
6  IMPLICIT NONE
7 C-----------------------------------------------
8 C D u m m y A r g u m e n t s
9 C-----------------------------------------------
10  REAL(dp), DIMENSION(nzeta,ntheta3,ns,0:1), INTENT(inout) :: r1, z1
11  REAL(dp), DIMENSION(nzeta,ntheta3,ns),INTENT(inout) :: ru0, zu0
12  LOGICAL, INTENT(in) :: lscreen
13 C-----------------------------------------------
14 C L o c a l P a r a m e t e r s
15 C-----------------------------------------------
16  INTEGER, PARAMETER :: limpts = 61
17  REAL(dp), PARAMETER :: p5 = 0.5_dp, two = 2
18 C-----------------------------------------------
19 C L o c a l V a r i a b l e s
20 C-----------------------------------------------
21  INTEGER :: i, j, k
22  INTEGER :: iv, iu, iu_r, ivminus, nlim, ns12, klim, n
23  REAL(dp), DIMENSION(nzeta) :: rcom, zcom
24  REAL(dp), DIMENSION(ntheta1) :: r1b, z1b, rub, zub
25  REAL(dp), DIMENSION(ntheta1) :: r12, z12
26  REAL(dp), DIMENSION(ntheta1) :: rs, zs, tau, ru12, zu12, tau0
27  REAL(dp) :: rlim, zlim
28  REAL(dp) :: rmax, rmin, zmax, zmin, dzeta
29  REAL(dp) :: ds, mintau, mintemp
30  INTEGER :: blksize, numjs, left, right, bcastrank
31  INTEGER, ALLOCATABLE, DIMENSION(:) :: counts, disps
32  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: send_buf
33  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: send_buf2
34  REAL(dp), ALLOCATABLE, DIMENSION(:) :: recv_buf
35  REAL(dp) :: tbroadon, tbroadoff, tguesson, tguessoff
36  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: bcastbuf
37  REAL(dp), ALLOCATABLE :: tmp(:)
38  REAL(dp) :: tmp2(nzeta,2)
39 C-----------------------------------------------
40 !
41 ! COMPUTES GUESS FOR MAGNETIC AXIS IF USER GUESS
42 ! LEADS TO INITIAL SIGN CHANGE OF JACOBIAN. DOES A GRID
43 ! SEARCH (irgrid, izgrid) IN EACH PHI-PLANE FOR POINTS WHICH
44 ! YIELD A VALUE FOR THE JACOBIAN WITH THE CORRECT SIGN (SIGNGS)
45 ! CHOOSES THE AXIS POSITION SO THE MIN VALUE OF THE JACOBIAN IS MAXIMIZED
46 !
47  CALL second0(tguesson)
48  ns12 = (ns + 1)/2
49 
50 #if defined(MPI_OPT)
51  IF (nranks.GT.1) THEN
52  bcastrank = -1
53  IF (ns.EQ.nranks) THEN
54  bcastrank=ns12
55  ELSE IF (ns.GT.nranks) THEN
56  !!! Can make this log P (for later) : SKS : NOV 13, 2014
57  DO i=1, nranks
58  IF (tlglob_arr(i) .LE. ns12 .AND.
59  & ns12 .LE. trglob_arr(i)) THEN
60  bcastrank = i - 1
61  EXIT
62  END IF
63  END DO
64  ELSE
65  IF (rank.EQ.0) THEN
66  WRITE(*,*) 'Something wrong in guess_axis'
67  END IF
68  CALL stopmpi(666)
69  END IF
70 
71  ALLOCATE(bcastbuf(nzeta,ntheta3,6))
72  bcastbuf(:,:,1) = r1(:,:,ns12,0)
73  bcastbuf(:,:,2) = r1(:,:,ns12,1)
74  bcastbuf(:,:,3) = z1(:,:,ns12,0)
75  bcastbuf(:,:,4) = z1(:,:,ns12,1)
76  bcastbuf(:,:,5) = ru0(:,:,ns12)
77  bcastbuf(:,:,6) = zu0(:,:,ns12)
78 
79  CALL second0(tbroadon)
80  CALL mpi_bcast(bcastbuf, 6*nznt, mpi_real8, bcastrank,
81  & ns_comm, mpi_err)
82  CALL second0(tbroadoff)
83  broadcast_time = broadcast_time + (tbroadoff-tbroadon)
84 
85  r1(:,:,ns12,0) = bcastbuf(:,:,1)
86  r1(:,:,ns12,1) = bcastbuf(:,:,2)
87  z1(:,:,ns12,0) = bcastbuf(:,:,3)
88  z1(:,:,ns12,1) = bcastbuf(:,:,4)
89  ru0(:,:,ns12) = bcastbuf(:,:,5)
90  zu0(:,:,ns12) = bcastbuf(:,:,6)
91 
92  bcastbuf(:,:,1) = r1(:,:,ns,0)
93  bcastbuf(:,:,2) = r1(:,:,ns,1)
94  bcastbuf(:,:,3) = z1(:,:,ns,0)
95  bcastbuf(:,:,4) = z1(:,:,ns,1)
96  bcastbuf(:,:,5) = ru0(:,:,ns)
97  bcastbuf(:,:,6) = zu0(:,:,ns)
98 
99  CALL second0(tbroadon)
100  CALL mpi_bcast(bcastbuf, 6*nznt, mpi_real8, nranks - 1,
101  & ns_comm, mpi_err)
102  CALL second0(tbroadoff)
103  broadcast_time = broadcast_time + (tbroadoff - tbroadon)
104 
105  r1(:,:,ns,0) = bcastbuf(:,:,1)
106  r1(:,:,ns,1) = bcastbuf(:,:,2)
107  z1(:,:,ns,0) = bcastbuf(:,:,3)
108  z1(:,:,ns,1) = bcastbuf(:,:,4)
109  ru0(:,:,ns) = bcastbuf(:,:,5)
110  zu0(:,:,ns) = bcastbuf(:,:,6)
111  DEALLOCATE(bcastbuf)
112 
113  CALL second0(tbroadon)
114  CALL mpi_bcast(psqrts(1,ns12), 1, mpi_real8, bcastrank,
115  & ns_comm, mpi_err)
116  CALL second0(tbroadoff)
117  broadcast_time = broadcast_time + (tbroadoff - tbroadon)
118 
119  ALLOCATE(tmp(2*nzeta))
120  tmp(1:nzeta) = r1(:,1,1,0)
121  tmp(nzeta + 1:2*nzeta) = z1(:,1,1,0)
122  CALL second0(tbroadon)
123  CALL mpi_bcast(tmp, 2*nzeta, mpi_real8, 0, ns_comm, mpi_err)
124  CALL second0(tbroadoff)
125  broadcast_time = broadcast_time + (tbroadoff - tbroadon)
126  r1(:,1,1,0) = tmp(1:nzeta)
127  z1(:,1,1,0) = tmp(nzeta + 1:2*nzeta)
128  DEALLOCATE(tmp)
129  END IF
130 #endif
131 
132  planes: DO iv = 1, nzeta
133  IF (.not.lasym .and. iv .gt. nzeta/2 + 1) THEN
134  rcom(iv) = rcom(nzeta + 2 - iv)
135  zcom(iv) =-zcom(nzeta + 2 - iv)
136  cycle
137  END IF
138  r1b(:ntheta3) = r1(iv,:,ns,0) + r1(iv,:,ns,1)
139  z1b(:ntheta3) = z1(iv,:,ns,0) + z1(iv,:,ns,1)
140  r12(:ntheta3) = r1(iv,:,ns12,0)+r1(iv,:,ns12,1)*psqrts(1,ns12)
141  z12(:ntheta3) = z1(iv,:,ns12,0)+z1(iv,:,ns12,1)*psqrts(1,ns12)
142  rub(:ntheta3) = ru0(iv,:,ns)
143  zub(:ntheta3) = zu0(iv,:,ns)
144  ru12(:ntheta3) = p5*(ru0(iv,:,ns) + ru0(iv,:,ns12))
145  zu12(:ntheta3) = p5*(zu0(iv,:,ns) + zu0(iv,:,ns12))
146 
147  IF (.not.lasym) THEN
148 !
149 ! USE Z(v,-u) = -Z(twopi-v,u), R(v,-u) = R(twopi-v,u)
150 ! TO DO EXTEND R,Z, etc. OVER ALL THETA (NOT JUST 0,PI)
151 !
152  ivminus = mod(nzeta + 1 - iv,nzeta) + 1 !!(twopi-v)
153  DO iu = 1+ntheta2, ntheta1
154  iu_r = ntheta1 + 2 - iu
155  r1b(iu) = r1(ivminus,iu_r,ns,0) + r1(ivminus,iu_r,ns,1)
156  z1b(iu) =-(z1(ivminus,iu_r,ns,0) +
157  & z1(ivminus,iu_r,ns,1))
158  r12(iu) = r1(ivminus,iu_r,ns12,0)
159  & + r1(ivminus,iu_r,ns12,1)*psqrts(1,ns12)
160  z12(iu) =-(z1(ivminus,iu_r,ns12,0) +
161  & z1(ivminus,iu_r,ns12,1)*psqrts(1,ns12))
162  rub(iu) =-ru0(ivminus,iu_r,ns)
163  zub(iu) = zu0(ivminus,iu_r,ns)
164  ru12(iu) =-p5*(ru0(ivminus,iu_r,ns) +
165  & ru0(ivminus,iu_r,ns12))
166  zu12(iu) = p5*(zu0(ivminus,iu_r,ns) +
167  & zu0(ivminus,iu_r,ns12))
168  END DO
169  END IF
170 !
171 ! Scan over r-z grid for interior point
172 !
173  rmin = minval(r1b)
174  rmax = maxval(r1b)
175  zmin = minval(z1b)
176  zmax = maxval(z1b)
177  rcom(iv) = (rmax + rmin)/2
178  zcom(iv) = (zmax + zmin)/2
179 
180 !
181 ! Estimate jacobian based on boundary and 1/2 surface
182 !
183  ds = (ns - ns12)*hs
184  DO iu = 1, ntheta1
185  rs(iu) = (r1b(iu) - r12(iu))/ds + r1(iv,1,1,0)
186  zs(iu) = (z1b(iu) - z12(iu))/ds + z1(iv,1,1,0)
187  tau0(iu) = ru12(iu)*zs(iu) - zu12(iu)*rs(iu)
188  END DO
189 
190  mintau = 0
191 
192  DO nlim = 1, limpts
193  zlim = zmin + ((zmax - zmin)*(nlim-1))/(limpts - 1)
194  IF (.not.lasym .and.
195  & (iv .eq. 1 .or. iv .eq. nzeta/2 + 1)) THEN
196  zlim = 0
197  IF (nlim .gt. 1) EXIT
198  END IF
199 !
200 ! Find value of magnetic axis that maximizes the minimum jacobian value
201 !
202  DO klim = 1, limpts
203  rlim = rmin + ((rmax - rmin)*(klim - 1))/(limpts - 1)
204  tau = signgs*(tau0 - ru12(:)*zlim + zu12(:)*rlim)
205  mintemp = minval(tau)
206  IF (mintemp .gt. mintau) THEN
207  mintau = mintemp
208  rcom(iv) = rlim
209  zcom(iv) = zlim
210 ! If up-down symmetric and lasym=T, need this to pick z = 0
211  ELSE IF (mintemp .eq. mintau) THEN
212  IF (abs(zcom(iv)).gt.abs(zlim)) THEN
213  zcom(iv) = zlim
214  END IF
215  END IF
216  END DO
217  END DO
218 
219  END DO planes
220 
221 !Distribute to all processors, not just NS_COMM
222 #if defined(MPI_OPT)
223  tmp2(:,1) = rcom
224  tmp2(:,2) = zcom
225  CALL mpi_bcast(tmp2, 2*nzeta, mpi_real8, 0,
226  1 runvmec_comm_world, mpi_err)
227  rcom = tmp2(:,1)
228  zcom = tmp2(:,2)
229 #endif
230 !
231 ! FOURIER TRANSFORM RCOM, ZCOM
232 !
233  dzeta = two/nzeta
234  DO n = 0, ntor
235  raxis_cc(n) = dzeta*sum(cosnv(:,n)*rcom(:))/nscale(n)
236  zaxis_cs(n) =-dzeta*sum(sinnv(:,n)*zcom(:))/nscale(n)
237  raxis_cs(n) =-dzeta*sum(sinnv(:,n)*rcom(:))/nscale(n)
238  zaxis_cc(n) = dzeta*sum(cosnv(:,n)*zcom(:))/nscale(n)
239  IF (n .eq. 0 .or. n .eq. nzeta/2) THEN
240  raxis_cc(n) = p5*raxis_cc(n)
241  zaxis_cc(n) = p5*zaxis_cc(n)
242  END IF
243  END DO
244 
245  IF (grank == 0 .and. lscreen) THEN
246  print *,' ---- Improved AXIS Guess ----'
247  print *,' RAXIS_CC = ',raxis_cc(0:ntor)
248  print *,' ZAXIS_CS = ',zaxis_cs(0:ntor)
249  print *,' -----------------------------'
250  END IF
251 
252  CALL second0(tguessoff)
253  guess_axis_time = guess_axis_time + (tguessoff - tguesson)
254 
255  END SUBROUTINE guess_axis_par
256 
257  SUBROUTINE guess_axis(r1, z1, ru0, zu0)
258  USE vmec_main
259  USE vmec_params, ONLY: nscale, signgs
260  USE realspace, ONLY: sqrts
261  USE parallel_include_module
262  IMPLICIT NONE
263 C-----------------------------------------------
264 C D u m m y A r g u m e n t s
265 C-----------------------------------------------
266  REAL(dp), DIMENSION(ns,nzeta,ntheta3,0:1),
267  1 INTENT(in) :: r1, z1
268  REAL(dp), DIMENSION(ns,nzeta,ntheta3), INTENT(in) :: ru0, zu0
269 C-----------------------------------------------
270 C L o c a l P a r a m e t e r s
271 C-----------------------------------------------
272  INTEGER, PARAMETER :: limpts = 61
273  REAL(dp), PARAMETER :: p5 = 0.5_dp, two = 2
274 C-----------------------------------------------
275 C L o c a l V a r i a b l e s
276 C-----------------------------------------------
277  INTEGER :: i, j, k
278  INTEGER :: iv, iu, iu_r, ivminus, nlim, ns12, klim, n
279  REAL(dp), DIMENSION(nzeta) :: rcom, zcom
280  REAL(dp), DIMENSION(ntheta1) :: r1b, z1b, rub, zub
281  REAL(dp), DIMENSION(ntheta1) :: r12, z12
282  REAL(dp), DIMENSION(ntheta1) :: rs, zs, tau, ru12, zu12, tau0
283  REAL(dp) :: rlim, zlim
284  REAL(dp) :: rmax, rmin, zmax, zmin, dzeta
285  REAL(dp) :: ds, mintau, mintemp
286  REAL(dp) :: tguesson, tguessoff
287 C-----------------------------------------------
288 !
289 ! COMPUTES GUESS FOR MAGNETIC AXIS IF USER GUESS
290 ! LEADS TO INITIAL SIGN CHANGE OF JACOBIAN. DOES A GRID
291 ! SEARCH (irgrid, izgrid) IN EACH PHI-PLANE FOR POINTS WHICH
292 ! YIELD A VALUE FOR THE JACOBIAN WITH THE CORRECT SIGN (SIGNGS)
293 ! CHOOSES THE AXIS POSITION SO THE MIN VALUE OF THE JACOBIAN IS MAXIMIZED
294 !
295  CALL second0(tguesson)
296  ns12 = (ns+1)/2
297 
298  planes: DO iv = 1, nzeta
299  IF (.not.lasym .and. iv.gt.nzeta/2+1) THEN
300  rcom(iv) = rcom(nzeta+2-iv)
301  zcom(iv) =-zcom(nzeta+2-iv)
302  cycle
303  END IF
304  r1b(:ntheta3) = r1(ns,iv,:,0) + r1(ns,iv,:,1)
305  z1b(:ntheta3) = z1(ns,iv,:,0) + z1(ns,iv,:,1)
306  r12(:ntheta3) = r1(ns12,iv,:,0) + r1(ns12,iv,:,1)*sqrts(ns12)
307  z12(:ntheta3) = z1(ns12,iv,:,0) + z1(ns12,iv,:,1)*sqrts(ns12)
308  rub(:ntheta3) = ru0(ns,iv,:)
309  zub(:ntheta3) = zu0(ns,iv,:)
310  ru12(:ntheta3) = p5*(ru0(ns,iv,:) + ru0(ns12,iv,:))
311  zu12(:ntheta3) = p5*(zu0(ns,iv,:) + zu0(ns12,iv,:))
312 
313  IF (.not.lasym) THEN
314 !
315 ! USE Z(v,-u) = -Z(twopi-v,u), R(v,-u) = R(twopi-v,u)
316 ! TO DO EXTEND R,Z, etc. OVER ALL THETA (NOT JUST 0,PI)
317 !
318  ivminus = mod(nzeta + 1 - iv,nzeta) + 1 !!(twopi-v)
319  DO iu = 1+ntheta2, ntheta1
320  iu_r = ntheta1 + 2 - iu
321  r1b(iu) = r1(ns,ivminus,iu_r,0) + r1(ns,ivminus,iu_r,1)
322  z1b(iu) =-(z1(ns,ivminus,iu_r,0) + z1(ns,ivminus,iu_r,1))
323  r12(iu) = r1(ns12,ivminus,iu_r,0)
324  & + r1(ns12,ivminus,iu_r,1)*sqrts(ns12)
325  z12(iu) =-(z1(ns12,ivminus,iu_r,0)
326  & + z1(ns12,ivminus,iu_r,1)*sqrts(ns12))
327  rub(iu) =-ru0(ns,ivminus,iu_r)
328  zub(iu) = zu0(ns,ivminus,iu_r)
329  ru12(iu) = -p5*(ru0(ns,ivminus,iu_r)
330  & + ru0(ns12,ivminus,iu_r))
331  zu12(iu) = p5*(zu0(ns,ivminus,iu_r)
332  & + zu0(ns12,ivminus,iu_r))
333  END DO
334 
335  END IF
336 
337 !
338 ! Scan over r-z grid for interior point
339 !
340  rmin = minval(r1b)
341  rmax = maxval(r1b)
342  zmin = minval(z1b)
343  zmax = maxval(z1b)
344  rcom(iv) = (rmax + rmin)/2
345  zcom(iv) = (zmax + zmin)/2
346 
347 !
348 ! Estimate jacobian based on boundary and 1/2 surface
349 !
350  ds = (ns - ns12)*hs
351  DO iu = 1, ntheta1
352  rs(iu) = (r1b(iu) - r12(iu))/ds + r1(1,iv,1,0)
353  zs(iu) = (z1b(iu) - z12(iu))/ds + z1(1,iv,1,0)
354  tau0(iu) = ru12(iu)*zs(iu) - zu12(iu)*rs(iu)
355  END DO
356 
357  mintau = 0
358 
359  DO nlim = 1, limpts
360  zlim = zmin + ((zmax - zmin)*(nlim - 1))/(limpts - 1)
361  IF (.not.lasym .and.
362  & (iv .eq. 1 .or. iv .eq. nzeta/2 + 1)) THEN
363  zlim = 0
364  IF (nlim .gt. 1) EXIT
365  END IF
366 !
367 ! Find value of magnetic axis that maximizes the minimum jacobian value
368 !
369  DO klim = 1, limpts
370  rlim = rmin + ((rmax - rmin)*(klim - 1))/(limpts - 1)
371  tau = signgs*(tau0 - ru12(:)*zlim + zu12(:)*rlim)
372  mintemp = minval(tau)
373  IF (mintemp .gt. mintau) THEN
374  mintau = mintemp
375  rcom(iv) = rlim
376  zcom(iv) = zlim
377 ! If up-down symmetric and lasym=T, need this to pick z = 0
378  ELSE IF (mintemp .eq. mintau) THEN
379  IF (abs(zcom(iv)).gt.abs(zlim)) zcom(iv) = zlim
380  END IF
381  END DO
382  END DO
383 
384  END DO planes
385 
386 !
387 ! FOURIER TRANSFORM RCOM, ZCOM
388 !
389  dzeta = two/nzeta
390  DO n = 0, ntor
391  raxis_cc(n) = dzeta*sum(cosnv(:,n)*rcom(:))/nscale(n)
392  zaxis_cs(n) =-dzeta*sum(sinnv(:,n)*zcom(:))/nscale(n)
393  raxis_cs(n) =-dzeta*sum(sinnv(:,n)*rcom(:))/nscale(n)
394  zaxis_cc(n) = dzeta*sum(cosnv(:,n)*zcom(:))/nscale(n)
395  IF (n .eq. 0 .or. n .eq. nzeta/2) THEN
396  raxis_cc(n) = p5*raxis_cc(n)
397  zaxis_cc(n) = p5*zaxis_cc(n)
398  END IF
399  END DO
400 
401 ! 100 FORMAT(' n = ',i4,' raxis = ',1pe10.3,' zaxis = ',1pe10.3)
402 
403  CALL second0(tguessoff)
404  s_guess_axis_time = s_guess_axis_time + (tguessoff - tguesson)
405 
406  END SUBROUTINE guess_axis
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17