V3FIT
angle_constraints.f90
1  MODULE angle_constraints
2  USE vmec_main, ONLY: ns, mpol, ntor, dp, mpol1, lthreed, lasym
3  USE vmec_params, ONLY: signgs, ntmax, rcc, rss, zsc, zcs, rsc, rcs, zss, zcc
4  USE precon2d, ONLY: ictrl_prec2d
5  IMPLICIT NONE
6  INTEGER, PARAMETER :: pexp=4, m0=0, m1=1, m2=2, m3=3
7  LOGICAL, PARAMETER :: lorigin=.false.
8  INTEGER :: mrho, m, istat
9  REAL(dp), PARAMETER :: p5=0.5_dp, zero=0
10  REAL(dp), ALLOCATABLE :: t1m(:), t2m(:), cos_HB(:), sin_HB(:)
11  REAL(dp), ALLOCATABLE :: rz_array0(:,:,:), xtempa(:)
12  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: arhod, arhom, brhod, brhom, &
13  ard2, arm2, azd2, azm2, arhod2, arhom2, &
14  brd2, brm2, bzd2, bzm2, brhod2, brhom2
15  REAL(dp), DIMENSION(:), ALLOCATABLE :: crhod, sin2u, cos2u, sfact
16  REAL(dp) :: sqp5
17 #ifdef _HBANGLE
18  CONTAINS
19 
20 ! CALLED FROM FIXARAY
21  SUBROUTINE init_multipliers
22  IMPLICIT NONE
23  REAL(dp) :: dnorm, t0
24 ! NOTE: rho = SUM(m) rhomn*max(m,1)**pexp cos(mu-nv), etc
25 ! the extra m**pexp factor is to give good scaling
26 
27  IF (ALLOCATED(t1m)) RETURN
28 
29  sqp5 = sqrt(p5)
30  mrho = mpol1
31 !HB Angle constraint
32 ! mrho = mpol1-1
33 
34  ALLOCATE(t1m(0:mpol), t2m(0:mpol), cos_hb(0:mpol1), sin_hb(0:mpol1), stat=istat)
35  IF (istat .NE. 0) stop 'Allocation error in init_multipliers!'
36 
37 !
38 ! REF: Hirshman/Breslau paper
39 !
40  t1m(m0) = 0; t2m(m0) = 0
41  DO m = 1, mpol1+1
42  t0 = max(1,m-1)
43  t1m(m) = t0/m
44  t0 = m+1
45  t2m(m) = t0/m
46  t1m(m) = t1m(m)**pexp
47  t2m(m) = t2m(m)**pexp
48  END DO
49 
50 ! t1m(3) = 0
51 !SPH-TEST (IMPROVES CONDITION #?)
52 ! t1m(m2) = t2m(m2)
53 
54 !HB Constraint
55 ! t2m(mpol1) = 0; t2m(mpol) = 0
56 
57 ! dnorm = MAXVAL(t1m); dnorm = MAX(dnorm,MAXVAL(t2m))
58 ! t1m = t1m/dnorm; t2m = t2m/dnorm
59 
60  cos_hb(m0) = 1; sin_hb(m0) = 0
61 
62  DO m = 1, mrho
63  dnorm = sqrt(t1m(m+1)**2 + t2m(m-1)**2)
64  IF (dnorm .EQ. zero) cycle
65  cos_hb(m) = t1m(m+1)/dnorm
66  sin_hb(m) = t2m(m-1)/dnorm
67  END DO
68 
69  IF (mrho .NE. mpol1) THEN
70  cos_hb(mpol1) = 0; sin_hb(mpol1) = 1
71  END IF
72 
73 !DOESN'T CONVERGE IF THIS GOES BEFORE cos,sin_HB calculations (F_0 too larger???)
74 ! t2m(m0) = t1m(m2)
75 
76  END SUBROUTINE init_multipliers
77 
78  SUBROUTINE free_multipliers
79  IF (ALLOCATED(t1m)) DEALLOCATE (t1m, t2m, cos_hb, sin_hb)
80  IF (ALLOCATED(rz_array0)) DEALLOCATE (rz_array0)
81  IF (ALLOCATED(xtempa)) DEALLOCATE (xtempa)
82  IF (ALLOCATED(arhod)) DEALLOCATE(arhod, arhom, brhod, brhom, crhod, &
83  ard2, arm2, azd2, azm2, arhod2, arhom2, &
84  brd2, brm2, bzd2, bzm2, brhod2, brhom2)
85  IF (ALLOCATED(cos2u)) DEALLOCATE(cos2u, sin2u)
86  IF (ALLOCATED(sfact)) DEALLOCATE(sfact)
87 
88  END SUBROUTINE free_multipliers
89 
90  SUBROUTINE store_init_array(rzl_array)
91  USE vmec_main, ONLY: neqs2, nznt, nzeta, cosmui, sinmui, cosmu, sinmu, &
92  nzeta, ntheta2, ntheta3, hs
93  REAL(dp), DIMENSION(ns*(1+ntor),0:mpol1,3*ntmax), INTENT(inout) :: rzl_array
94  INTEGER :: nsp1, l, js
95 
96  IF (ALLOCATED(rz_array0)) DEALLOCATE (rz_array0)
97  ALLOCATE(rz_array0(ns*(1+ntor),0:mpol1,2*ntmax))
98  rz_array0 = rzl_array(:,:,1:2*ntmax)
99  rzl_array(:,:,1:2*ntmax) = 0
100 
101  CALL init_multipliers
102  CALL get_rep_mismatch(rz_array0, rzl_array)
103 
104  IF (ALLOCATED(xtempa)) DEALLOCATE (xtempa)
105  ALLOCATE (xtempa(neqs2)) !Used as temp storage of xc in funct3d
106 
107  nsp1=ns+1
108  IF (ALLOCATED(arhod)) DEALLOCATE(arhod, arhom, brhod, brhom, crhod, &
109  ard2, arm2, azd2, azm2, arhod2, arhom2, &
110  brd2, brm2, bzd2, bzm2, brhod2, brhom2)
111  ALLOCATE(arhod(nsp1,2),arhom(nsp1,2),brhod(nsp1,2),brhom(nsp1,2), &
112  ard2(nsp1,2), arm2(nsp1,2), arhod2(nsp1,2), arhom2(nsp1,2),&
113  azd2(nsp1,2), azm2(nsp1,2), &
114  brd2(nsp1,2), brm2(nsp1,2), brhod2(nsp1,2), brhom2(nsp1,2),&
115  bzd2(nsp1,2), bzm2(nsp1,2), crhod(nsp1))
116  IF (ALLOCATED(cos2u)) DEALLOCATE(cos2u, sin2u)
117  ALLOCATE (cos2u(nznt), sin2u(nznt))
118  DO l=1,nzeta
119  cos2u(l:nznt:nzeta) = cosmui(:,m2)
120  sin2u(l:nznt:nzeta) = sinmui(:,m2)
121  END DO
122 
123  IF (ALLOCATED(sfact)) DEALLOCATE(sfact)
124  ALLOCATE(sfact(ns))
125  DO js=1,ns
126  sfact(js) = hs*(js-1)
127  END DO
128 
129  END SUBROUTINE store_init_array
130 
131  SUBROUTINE getrz (rz_array)
132 ! USE vmec_main, ONLY: iter2
133 ! USE xstuff, ONLY: xc, xstore
134  IMPLICIT NONE
135  REAL(dp), DIMENSION(ns*(1+ntor),0:mpol1,2*ntmax), INTENT(inout) :: rz_array
136  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: rhocc, rhoss, rhocs, rhosc
137  REAL(dp), DIMENSION(:), ALLOCATABLE :: r0c, r0s, z0c, z0s
138  INTEGER :: nsnt, mrho1, istat
139 
140 ! INPUT: rho quasi-polar components (Fourier modes) stored in R
141 ! R,Z(m=0) centroid components stored in Z(m=0)
142 ! OUTPUT: Cylindrical component Fourier modes of R, Z
143 ! REF: HB Paper, Eq 32
144 !
145 ! FOR STELLARATOR SYMMETRY, rho = cos(mu-nv), R' = sum(m^pexp Rmn cos(mu-nv)), Z' = sum(m^pexp Zmn sin(mu-nv)
146 ! R' = (rhocc + rhoss) cosu
147 ! Z' = (rhocc + rhoss) sinu
148 !
149 ! FOR ASYMMETRY, ADD rho = sin(mu-nv) TERMS
150 ! R' = (rhosc + rhocs) cosu
151 ! Z' = (rhosc + rhocs) sinu
152 !
153 
154  nsnt = ns*(1+ntor)
155  mrho1 = mrho+1
156 ! IF (mrho .NE. mpol1) STOP 'mrho != mpol1'
157 
158 ! Enforce asymptotic behavior near axis (HB paper after Eq (34))
159  IF (lorigin) THEN
160  rz_array(2:nsnt:ns,m2,1:ntmax) = sqp5*rz_array(3:nsnt:ns,m2,1:ntmax)
161  IF (mpol1 .GE. m3) THEN
162  rz_array(2:nsnt:ns,m1:m3:2,1:ntmax) = p5*rz_array(3:nsnt:ns,m1:m3:2,1:ntmax)
163  rz_array(2:nsnt:ns,m3+1:,1:ntmax) = 0
164  END IF
165  END IF
166 
167  ALLOCATE (rhocc(nsnt,0:mpol1+1), rhoss(nsnt,0:mpol1+1), &
168  r0c(nsnt), z0s(nsnt),stat=istat)
169  IF (istat .NE. 0) stop 'Allocation Error #1 in GETRZ'
170  rhocc(:,0:mrho) = rz_array(:,0:mrho,rcc)
171  r0c = rz_array(:,m0,zsc+ntmax) !Used for temp storage
172  rhocc(:,mrho1) = 0; rhocc(:,mpol1+1) = 0
173 
174  IF (lthreed) THEN
175  rhoss(:,1:mrho) = rz_array(:,1:mrho,rss)
176  rhoss(:,m0) = 0
177  rhoss(:,mrho1) = 0; rhoss(:,0) = 0
178  z0s = rz_array(:,m0,zcs+ntmax)
179  END IF
180 
181  modes: DO m = 0, mpol1
182  IF (m .EQ. m0) THEN
183  rz_array(:,m0,rcc) = r0c + t2m(m0)*rhocc(:,m1)
184  rz_array(:,m0,zsc+ntmax) = 0
185  IF (lthreed) THEN
186  rz_array(:,m0,zcs+ntmax) = z0s - t2m(m0)*rhoss(:,m1)*signgs
187  rz_array(:,m0,rss) = 0
188  END IF
189  ELSE
190  rz_array(:,m,rcc) = (t1m(m)*rhocc(:,m-1) + t2m(m)*rhocc(:,m+1))
191  rz_array(:,m,zsc+ntmax) =-(t1m(m)*rhocc(:,m-1) - t2m(m)*rhocc(:,m+1))*signgs
192  IF (lthreed) THEN
193  rz_array(:,m,rss) = (t1m(m)*rhoss(:,m-1) + t2m(m)*rhoss(:,m+1))
194  rz_array(:,m,zcs+ntmax) = (t1m(m)*rhoss(:,m-1) - t2m(m)*rhoss(:,m+1))*signgs
195  END IF
196  ENDIF
197  END DO modes
198 
199  DEALLOCATE (rhocc, rhoss, r0c, z0s)
200 
201  IF (.NOT.lasym) GOTO 1002
202 
203  ALLOCATE (rhosc(nsnt,0:mrho1), rhocs(nsnt,0:mrho1), &
204  r0s(nsnt), z0c(nsnt), stat=istat)
205  IF (istat .NE. 0) stop 'Allocation Error #2 in GETRZ'
206 
207  rhosc(:,0:mrho) = rz_array(:,0:mrho,rsc)
208  rhosc(:,mrho1) = 0
209  z0c = rz_array(:,m0,zcc+ntmax)
210  IF (lthreed) THEN
211  rhocs(:,0:mrho) = rz_array(:,0:mrho,rcs)
212  rhocs(:,mrho1) = 0; rhocs(:,0) = 0
213  r0s = rz_array(:,m0,zss+ntmax)
214  END IF
215 
216  modea: DO m = 0, mpol1
217  IF (m .EQ. 0) THEN
218  rz_array(:,m0,zcc+ntmax) = z0c - t2m(m0)*rhosc(:,m1)*signgs
219  rz_array(:,m0,rsc) = 0
220  IF (lthreed) THEN
221  rz_array(:,m0,rcs) = r0s + t2m(m0)*rhocs(:,m1)
222  rz_array(:,m0,zss+ntmax) = 0
223  END IF
224  ELSE
225  rz_array(:,m,rsc) = (t1m(m)*rhosc(:,m-1) + t2m(m)*rhosc(:,m+1))
226  rz_array(:,m,zcc+ntmax) = (t1m(m)*rhosc(:,m-1) - t2m(m)*rhosc(:,m+1))*signgs
227  IF (lthreed) THEN
228  rz_array(:,m,rcs) = (t1m(m)*rhocs(:,m-1) + t2m(m)*rhocs(:,m+1))
229  rz_array(:,m,zss+ntmax) =-(t1m(m)*rhocs(:,m-1) - t2m(m)*rhocs(:,m+1))*signgs
230  END IF
231  ENDIF
232  END DO modea
233 
234  DEALLOCATE (rhosc, rhocs, r0s, z0c)
235 
236  1002 CONTINUE
237 
238  !Add initial (scaled) boundary values if needed
239  IF (ictrl_prec2d .NE. 3) rz_array = rz_array + rz_array0
240 
241  END SUBROUTINE getrz
242 
243  SUBROUTINE getfrho (frho_array)
244  USE vmec_main, ONLY: iter2, fnorm, hs
245  REAL(dp), DIMENSION(ns*(1+ntor),0:mpol1,2*ntmax), INTENT(inout) :: &
246  frho_array
247  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: frcc, frss, fzcs, fzsc, &
248  frsc, frcs, fzcc, fzss
249  INTEGER :: nsnt, mrho1, istat
250 
251 ! INPUT: Frho contains cylindrical components of FR, FZ
252 ! OUTPUT: Frho(1:ntmax) contains Quasi-polar components the MHD forces
253 ! Frho(ntmax+1:2*ntmax) contains centroid (m=0) forces
254  nsnt = ns*(1+ntor)
255  mrho1 = mrho+1
256 
257  ALLOCATE (frcc(nsnt,0:mrho1), frss(nsnt,0:mrho1), &
258  fzsc(nsnt,0:mrho1), fzcs(nsnt,0:mrho1), stat=istat)
259  IF (istat .NE. 0) stop 'Allocation error #1 in getfrho!'
260 
261  frcc(:,0:mrho) = frho_array(:,0:mrho,rcc); frcc(:,mrho1)=0
262  frho_array(:,:,rcc)=0
263 
264  fzsc(:,0:mrho) = frho_array(:,0:mrho,zsc+ntmax); fzsc(:,mrho1) = 0
265  frho_array(:,:,zsc+ntmax)=0
266 
267  IF (lthreed) THEN
268  frss(:,0:mrho) = frho_array(:,0:mrho,rss); frss(:,mrho1) = 0
269  frho_array(:,:,rss)=0
270  fzcs(:,0:mrho) = frho_array(:,0:mrho,zcs+ntmax); fzcs(:,mrho1) = 0
271  frho_array(:,:,zcs+ntmax)=0
272  END IF
273 
274  modes: DO m = 0, mrho
275  IF (m .EQ. m0) THEN
276  frho_array(:,m0,zsc+ntmax) = frcc(:,m0) !storage
277  frho_array(:,m0,rcc) = cos_hb(m0)*(frcc(:,m1) - signgs*fzsc(:,m1))
278  IF (.NOT.lthreed) cycle
279  frho_array(:,m0,zcs+ntmax) = fzcs(:,m0)
280  ELSE
281  frho_array(:,m,rcc) = cos_hb(m)*(frcc(:,m+1) - signgs*fzsc(:,m+1)) &
282  + sin_hb(m)*(frcc(:,m-1) + signgs*fzsc(:,m-1))
283  IF (.NOT.lthreed) cycle
284  frho_array(:,m,rss) = cos_hb(m)*(frss(:,m+1) + signgs*fzcs(:,m+1)) &
285  + sin_hb(m)*(frss(:,m-1) - signgs*fzcs(:,m-1))
286  ENDIF
287  END DO modes
288 
289  DEALLOCATE (frcc, frss, fzsc, fzcs)
290  IF (lthreed) THEN
291  IF (any(frho_array(:,m0,rss) .NE. zero)) stop 'FRHO(m0,rss) != 0'
292  END IF
293 
294  IF (.NOT. lasym) GOTO 1000
295 
296  ALLOCATE (frsc(nsnt,0:mrho1), frcs(nsnt,0:mrho1), &
297  fzcc(nsnt,0:mrho1), fzss(nsnt,0:mrho1), stat=istat)
298  IF (istat .NE. 0) stop 'Allocation error #2 in getfrho!'
299 
300  frsc(:,0:mrho) = frho_array(:,0:mrho,rsc); frsc(:,mrho1) = 0
301  frho_array(:,:,rsc) = 0
302  fzcc(:,0:mrho) = frho_array(:,0:mrho,zcc+ntmax); fzcc(:,mrho1) = 0
303  frho_array(:,:,zcc+ntmax) = 0
304  IF (lthreed) THEN
305  frcs(:,0:mrho) = frho_array(:,0:mrho,rcs); frcs(:,mrho1) = 0
306  frho_array(:,:,rcs)=0
307  fzss(:,0:mrho) = frho_array(:,0:mrho,zss+ntmax); fzss(:,mrho1) = 0
308  frho_array(:,:,zss+ntmax)=0
309  END IF
310 
311  modea: DO m = 0, mrho
312  IF (m .EQ. m0) THEN
313  frho_array(:,m0,zcc+ntmax) = fzcc(:,m0)
314  IF (.NOT.lthreed) cycle
315  frho_array(:,m0,zss+ntmax) = frcs(:,m0)
316  frho_array(:,m0,rcs) =cos_hb(m0)*(frcs(:,m1) - signgs*fzss(:,m1))
317  ELSE
318  frho_array(:,m,rsc) = cos_hb(m)*(frsc(:,m+1) + signgs*fzcc(:,m+1)) &
319  + sin_hb(m)*(frsc(:,m-1) - signgs*fzcc(:,m-1))
320  IF (.NOT.lthreed) cycle
321  frho_array(:,m,rcs) = cos_hb(m)*(frcs(:,m+1) - signgs*fzss(:,m+1)) &
322  + sin_hb(m)*(frcs(:,m-1) + signgs*fzss(:,m-1))
323  ENDIF
324  END DO modea
325 
326  DEALLOCATE (frsc, frcs, fzcc, fzss)
327 
328  IF (lasym) THEN
329  IF (any(frho_array(:,m0,rsc) .NE. zero)) stop 'FRHO(m=0,rsc) != 0'
330  END IF
331 
332  1000 CONTINUE
333 
334 ! Use asymptotic behavior near axis, rather than evolution (SPH010214)
335  IF (lorigin) THEN
336  frho_array(2:nsnt:ns,m1:,1:ntmax) = 0
337  END IF
338 
339 !ADD UNIQUE POLAR AXIS "CONSTRAINT"
340  frho_array(:,m1,1:ntmax) = 0 !Unique angle
341 
342 
343  END SUBROUTINE getfrho
344 
345 
346  SUBROUTINE scalfor_rho(gcr, gcz)
347  USE vmec_main
348  USE vmec_params
349  IMPLICIT NONE
350 !-----------------------------------------------
351 ! D u m m y A r g u m e n t s
352 !-----------------------------------------------
353  REAL(dp), DIMENSION(ns,0:ntor,0:mpol1,ntmax), INTENT(inout) :: &
354  gcr, gcz
355 !-----------------------------------------------
356 ! L o c a l V a r i a b l e s
357 !-----------------------------------------------
358  REAL(dp), PARAMETER :: edge_pedestal=0.05_dp
359  INTEGER :: m , mp, n, js, jmax, jmin4(0:mnsize-1)
360  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: arho, brho, drho, &
361  arho2,brho2,drho2
362  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: acen, bcen, dcen
363  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: gcen
364  REAL(dp) :: tar, taz, tc, tc1
365 !-----------------------------------------------
366  jmax = ns
367  IF (ivac .lt. 1) jmax = ns1
368 
369  ALLOCATE (acen(ns,0:ntor), bcen(ns,0:ntor), dcen(ns,0:ntor), &
370  gcen(ns,0:ntor,4))
371 !
372 ! FIRST, SCALE m=0 R,Z COMPONENTS
373 !
374  gcen=0
375  gcen(:,:,1) = gcz(:,:,m0,zsc) !used for r0n storage
376  IF (lasym) gcen(:,:,3) = gcz(:,:,m0,zcc)
377  IF (lthreed) THEN
378  gcen(:,:,2) = gcz(:,:,m0,zcs)
379  IF (lasym) gcen(:,:,4) = gcz(:,:,m0,zss)
380  END IF
381 
382  IF (sum(gcz*gcz) .NE. sum(gcen*gcen)) stop 'ERROR #1 IN SCALFOR_RHO'
383 
384  CALL scalaxis(arm, ard, crd, acen, bcen, dcen, jmax, gcen(:,:,1), 0)
385  IF (lasym) CALL scalaxis(azm, azd, crd, acen, bcen, dcen, jmax, gcen(:,:,3), 0)
386  IF (lthreed) THEN
387  CALL scalaxis(azm, azd, crd, acen, bcen, dcen, jmax, gcen(:,:,2), 1)
388  IF (lasym) CALL scalaxis(arm, ard, crd, acen, bcen, dcen, jmax, gcen(:,:,4), 1)
389  END IF
390 
391  gcz = 0
392 !
393 ! RESTORE SCALED CENTROID VALUES
394 !
395  gcz(:,:,m0,zsc) = gcen(:,:,1) !R0n-cc
396  IF (lasym) gcz(:,:,m0,zcc) = gcen(:,:,3) !Z0n-cc
397  IF (lthreed) THEN
398  gcz(:,:,m0,zcs) = gcen(:,:,2) !Z0n-cs
399  IF (lasym) gcz(:,:,m0,zss) = gcen(:,:,4) !R0n-cs
400  END IF
401 
402 !
403 ! NEXT, SCALE QUASI-POLAR COMPONENTS
404 !
405  ALLOCATE (arho(ns,0:ntor,0:mpol1), brho(ns,0:ntor,0:mpol1), &
406  drho(ns,0:ntor,0:mpol1), arho2(ns,0:ntor,0:mpol1), &
407  brho2(ns,0:ntor,0:mpol1),drho2(ns,0:ntor,0:mpol1))
408 
409  DO m = 0, mpol1
410  mp = mod(m,2)+1
411  DO n = 0, ntor
412  DO js = 1, ns
413  arho(js,n,m) = -(arhom(js+1,mp) + brhom(js+1,mp)*m**2)
414  brho(js,n,m) = -(arhom(js,mp) + brhom(js,mp) *m**2)
415  drho(js,n,m) = -(arhod(js,mp) + brhod(js,mp) *m**2 &
416  + crhod(js)*(n*nfp)**2)
417  arho2(js,n,m)= -p5*(arhom2(js+1,mp)+ brhom2(js+1,mp)*m**2)
418  brho2(js,n,m)= -p5*(arhom2(js,mp) + brhom2(js,mp) *m**2)
419  drho2(js,n,m)= -p5*(arhod2(js,mp) + brhod2(js,mp) *m**2)
420  END DO
421  END DO
422  END DO
423 
424  CALL avg_rho(arho, arho2)
425  CALL avg_rho(brho, brho2)
426  CALL avg_rho(drho, drho2)
427 
428  IF (jmax .ge. ns) THEN
429 !
430 ! SMALL EDGE PEDESTAL NEEDED TO IMPROVE CONVERGENCE
431 ! IN PARTICULAR, NEEDED TO ACCOUNT FOR POTENTIAL ZERO
432 ! EIGENVALUE DUE TO NEUMANN (GRADIENT) CONDITION AT EDGE
433 !
434  drho(ns,:,:) = (1+edge_pedestal)*drho(ns,:,:)
435 ! drho(ns,:,3:) = (1+2*edge_pedestal)*drho(ns,:,3:)
436  END IF
437 
438  jmin4 = 2
439  CALL tridslv (arho,drho,brho,gcr(:,:,:,rcc),jmin4,jmax, &
440  mnsize-1,ns,1)
441  IF (lthreed) CALL tridslv (arho,drho,brho,gcr(:,:,:,rss), &
442  jmin4,jmax,mnsize-1,ns,1)
443  IF (lasym) THEN
444  CALL tridslv (arho,drho,brho,gcr(:,:,:,rsc),jmin4,jmax, &
445  mnsize-1,ns,1)
446  IF (lthreed) CALL tridslv (arho,drho,brho,gcr(:,:,:,rcs), &
447  jmin4,jmax,mnsize-1,ns,1)
448  END IF
449 
450  DEALLOCATE (arho, brho, drho, arho2, brho2, drho2, gcen)
451 
452  END SUBROUTINE scalfor_rho
453 
454  SUBROUTINE scalaxis(axm, axd, cxd, acen, bcen, dcen, jmax, gcen, iflag)
455  USE vmec_main
456  USE vmec_params
457  REAL(dp), PARAMETER :: edge_pedestal=0.15_dp, fac=0.25_dp
458 ! REAL(dp), PARAMETER :: edge_pedestal=0.05_dp, fac=0.25_dp
459  REAL(dp), INTENT(IN), DIMENSION(ns+1,2) :: axm, axd
460  REAL(dp), INTENT(IN), DIMENSION(ns+1) :: cxd
461  REAL(dp), DIMENSION(ns,0:ntor) :: acen, bcen, dcen
462  REAL(dp), INTENT(INOUT) :: gcen(ns,0:ntor)
463  REAL(dp) :: mult_fac
464  INTEGER, INTENT(IN) :: jmax, iflag
465  INTEGER :: mp, n, js, mnsize0, jmin4(0:mnsize-1)
466 
467  mp=1
468  DO n = 0, ntor
469  DO js = 1, jmax
470  acen(js,n) = -axm(js+1,mp)
471  bcen(js,n) = -axm(js,mp)
472  dcen(js,n) = -(axd(js,mp)+cxd(js)*(n*nfp)**2)
473  END DO
474  END DO
475 
476  IF (jmax .ge. ns) THEN
477 !
478 ! SMALL EDGE PEDESTAL NEEDED TO IMPROVE CONVERGENCE
479 ! IN PARTICULAR, NEEDED TO ACCOUNT FOR POTENTIAL ZERO
480 ! EIGENVALUE DUE TO NEUMANN (GRADIENT) CONDITION AT EDGE
481 !
482  dcen(ns,:) = (1+edge_pedestal)*dcen(ns,:)
483 !
484 ! STABILIZATION ALGORITHM FOR ZC_00(NS)
485 ! FOR UNSTABLE CASE, HAVE TO FLIP SIGN OF -FAC -> +FAC FOR CONVERGENCE
486 ! COEFFICIENT OF < Ru (R Pvac)> ~ -fac*(z-zeq) WHERE fac (EIGENVALUE, OR
487 ! FIELD INDEX) DEPENDS ON THE EQUILIBRIUM MAGNETIC FIELD AND CURRENT,
488 ! AND zeq IS THE EQUILIBRIUM EDGE VALUE OF Z00
489  IF (iflag .eq. 1) THEN
490 !
491 ! METHOD 1: SUBTRACT (INSTABILITY) Pedge ~ fac*z/hs FROM PRECONDITIONER AT EDGE
492 !
493  mult_fac = min(fac, fac*hs*15)
494  dcen(ns,0) = dcen(ns,0)*(1-mult_fac)/(1+edge_pedestal)
495  END IF
496 
497  ENDIF
498 
499 
500  mnsize0 = 1+ntor
501  jmin4 = jmin3
502  IF (iresidue.GE.0 .AND. iresidue.LT.3) jmin4(0)=2
503 
504  CALL tridslv(acen,dcen,bcen,gcen,jmin4,jmax,mnsize0-1,ns,1)
505 
506  END SUBROUTINE scalaxis
507 
508  SUBROUTINE avg_rho(ax, ax2)
509  USE vmec_main
510  IMPLICIT NONE
511  REAL(dp), INTENT(inout), DIMENSION(ns*(1+ntor), 0:mpol1) :: ax, ax2
512  REAL(dp), ALLOCATABLE :: ax1(:,:)
513  INTEGER :: m
514 
515  ALLOCATE (ax1(ns*(1+ntor),0:mpol1))
516 
517  ax1(:,m0) = cos_hb(m0)*t1m(m1)*(ax(:,m1)+ax2(:,m1))
518  DO m=1,mpol1-1
519  ax1(:,m) = cos_hb(m)*(t1m(m+1)*ax(:,m+1)+t2m(m-1)*ax2(:,m-1)) &
520  + sin_hb(m)*(t2m(m-1)*ax(:,m-1)+t1m(m+1)*ax2(:,m+1))
521  END DO
522 ! ax1(:,m2) = ax1(:,m2)+sin_HB(m2)*t2m(m1)*ax2(:,m2-1)
523  ax1(:,mpol1) = sin_hb(mpol1)*t2m(mpol1-1)*ax(:,mpol1-1) &
524  + cos_hb(mpol1)*t2m(mpol1-1)*ax2(:,mpol1-1)
525  ax = ax1
526 
527  DEALLOCATE (ax1)
528 
529  END SUBROUTINE avg_rho
530 
531  SUBROUTINE precondn_rho
532  USE vmec_main, ONLY: ard, arm, brd, brm, azd, azm, bzd, bzm, crd
533  USE fbal, ONLY: rzu_fac, rru_fac, frcc_fac, fzsc_fac
534  REAL(dp), PARAMETER :: one=1
535 
536  arhod = ard + azd
537  arhom = arm + azm
538  brhod = brd + bzd
539  brhom = brm + bzm
540  crhod = crd + crd
541 
542  arhod2 = ard2-azd2; brhod2 = (brd2-bzd2)
543  arhom2 = arm2-azm2; brhom2 = (brm2-bzm2)
544 
545  rzu_fac = (rzu_fac+rru_fac); rru_fac=rzu_fac
546  frcc_fac(2:ns-1) = one/rzu_fac(2:ns-1)
547  fzsc_fac(2:ns-1) = -frcc_fac(2:ns-1)
548 
549  END SUBROUTINE precondn_rho
550 
551  SUBROUTINE get_rep_mismatch(rz0_array, rho_array)
552  USE vmec_main, ONLY: irst, hs
553  IMPLICIT NONE
554  REAL(dp), PARAMETER :: p5=0.5_dp, p25=p5*p5
555  REAL(dp), DIMENSION(ns*(1+ntor),0:mpol1,2*ntmax), INTENT(in) :: rz0_array
556  REAL(dp), DIMENSION(ns*(1+ntor),0:mpol1,2*ntmax), INTENT(out) :: rho_array
557  INTEGER :: m, nsnt, mrho1, js, n, ntc
558  REAL(dp) :: match, delta, t1(ns*(1+ntor)), t2(ns*(1+ntor)), &
559  temp1(ns*(1+ntor)), temp2(ns*(1+ntor)), es
560 !
561 ! COMPUTES mismatch BETWEEN INITIAL R-Z REPRESENTATION AND QPOLAR FORM
562 !
563  nsnt = SIZE(rz0_array,1)
564  mrho1 = max(mrho+1,mpol1)
565  temp1 = 0; temp2 = 0
566 
567 !
568 ! STORE m=0 AXIS data
569 !
570  rho_array(:,m0,zsc+ntmax) = rz0_array(:,m0,rcc)
571  IF (lthreed) THEN
572  rho_array(:,m0,zcs+ntmax) = rz0_array(:,m0,zcs+ntmax)
573  END IF
574 
575  DO m = 0, mrho
576  IF (m .LT. mpol1-1) THEN
577  t1 = (rz0_array(:,m+1,rcc) - signgs*rz0_array(:,m+1,zsc+ntmax))/t1m(m+1)
578  IF (m .LE. 1) THEN
579  t2 = t1
580  ELSE
581  t2 = (rz0_array(:,m-1,rcc) + signgs*rz0_array(:,m-1,zsc+ntmax))/t2m(m-1)
582  END IF
583  ELSE IF (m .GT. 1) THEN
584  t2 = (rz0_array(:,m-1,rcc) + signgs*rz0_array(:,m-1,zsc+ntmax))/t2m(m-1)
585  t1 = t2
586  END IF
587  rho_array(:,m,rcc) = p25*(t1 + t2)
588  temp1 = temp1 + (t1 - t2)**2
589  temp2 = temp2 + (t1 + t2)**2
590 
591  IF (.NOT.lthreed) cycle
592 
593  IF (m .LT. mpol1-1) THEN
594  t1 = (rz0_array(:,m+1,rss) + signgs*rz0_array(:,m+1,zcs+ntmax))/t1m(m+1)
595  IF (m .LE. 1) THEN
596  t2 = t1
597  ELSE
598  t2 = (rz0_array(:,m-1,rss) - signgs*rz0_array(:,m-1,zcs+ntmax))/t2m(m-1)
599  END IF
600  ELSE
601  t2 = (rz0_array(:,m-1,rss) - rz0_array(:,m-1,rcs+ntmax))/t2m(m-1)
602  t1 = t2
603  END IF
604  rho_array(:,m,rss) = p25*(t1 + t2)
605  temp1 = temp1 + (t1 - t2)**2
606  temp2 = temp2 + (t1 + t2)**2
607  END DO
608 
609  rho_array(:,m0,zsc+ntmax) = rho_array(:,m0,zsc+ntmax) - t2m(m0)*rho_array(:,m1,rcc)
610  IF (lthreed) THEN
611  rho_array(:,m0,zcs+ntmax) = rho_array(:,m0,zcs+ntmax) + t2m(m0)*rho_array(:,m1,rss)*signgs
612  END IF
613 
614 !
615 ! NON-SYMMETRIC CONTRIBUTIONS
616 !
617  IF (.NOT.lasym) GOTO 1000
618 
619  rho_array(:,m0,zcc+ntmax) = rz0_array(:,m0,zcc+ntmax)
620  IF (lthreed) THEN
621  rho_array(:,m0,zss+ntmax) = rz0_array(:,m0,rcs)
622  END IF
623 
624  DO m = 0, mrho
625  IF (m .LT. mpol1-1) THEN
626  t1 = (rz0_array(:,m+1,rsc) + signgs*rz0_array(:,m+1,zcc+ntmax))/t1m(m+1)
627  IF (m .LE. 1) THEN
628  t2 = t1
629  ELSE
630  t2 = (rz0_array(:,m-1,rsc) - signgs*rz0_array(:,m-1,zcc+ntmax))/t2m(m-1)
631  END IF
632  ELSE
633  t2 = (rz0_array(:,m-1,rsc) - signgs*rz0_array(:,m-1,zcc+ntmax))/t2m(m-1)
634  t1 = t2
635  END IF
636  rho_array(:,m,rsc) = p25*(t1 + t2)
637  temp1 = temp1 + (t1 - t2)**2
638  temp2 = temp2 + (t1 + t2)**2
639 
640  IF (.not.lthreed) cycle
641 
642  IF (m .LT. mpol1-1) THEN
643  t1 = (rz0_array(:,m+1,rcs) - signgs*rz0_array(:,m+1,zss+ntmax))/t1m(m+1)
644  IF (m .LE. 1) THEN
645  t2 = t1
646  ELSE
647  t2 = (rz0_array(:,m-1,rcs) + signgs*rz0_array(:,m-1,zss+ntmax))/t2m(m-1)
648  END IF
649  ELSE
650  t2 = (rz0_array(:,m-1,rcs) + signgs*rz0_array(:,m-1,zss+ntmax))/t2m(m-1)
651  t1 = t2
652  END IF
653  rho_array(:,m,rcs) = p25*(t1 + t2)
654  temp1 = temp1 + (t1 - t2)**2
655  temp2 = temp2 + (t1 + t2)**2
656  END DO
657 
658  rho_array(:,m0,zcc+ntmax) = rho_array(:,m0,zcc+ntmax) + t2m(m0)*rho_array(:,m1,rsc)
659  IF (lthreed) THEN
660  rho_array(:,m0,zss+ntmax) = rho_array(:,m0,zss+ntmax) - t2m(m0)*rho_array(:,m1,rcs)*signgs
661  END IF
662 
663  1000 CONTINUE
664 
665  IF (irst .EQ. 1) GOTO 2000
666 
667  DO js = 2, ns-1
668  es = (js-1)*hs
669  DO n = 0, ntor
670  ntc = js+ns*ntor
671  DO m = 0, mrho
672  IF (mod(m,2) .EQ. 0) THEN
673  rho_array(ntc, m, 1:ntmax) = sqrt(es)*rho_array(ns*(1+ntor), m, 1:ntmax)
674  ELSE
675  rho_array(ntc, m, 1:ntmax) = es*rho_array(ns*(1+ntor), m, 1:ntmax)
676  END IF
677  END DO
678  END DO
679  END DO
680 
681  2000 CONTINUE
682 
683  WRITE (*, '(a)') ' Quasi-polar representation mismatch vs radius'
684  DO js=2,ns,ns-2
685  delta = sum(temp1(js:nsnt:ns))
686  match = sum(temp2(js:nsnt:ns))
687  IF (match .NE. 0._dp) WRITE (*, '(1x,i4,1p,e10.2)')js,delta/match
688  END DO
689 
690 ! rho_array=0 !Use this if delta at the edge is not good
691  rz_array0=0 !Initial jacobian can change sign: check
692 
693 ! IF (ALLOCATED(rho1)) DEALLOCATE(rho1)
694 ! ALLOCATE (rho1(nsnt,ntmax))
695 ! rho1 = rho_array(:,m1,1:ntmax)
696 
697  END SUBROUTINE get_rep_mismatch
698 #endif
699  END MODULE angle_constraints
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17