VMEC 8.52
3D Equilibrium solver with nested flux surfaces.
Loading...
Searching...
No Matches
totzsp.f90
Go to the documentation of this file.
1
3
32SUBROUTINE totzsps(rzl_array, r11, ru1, rv1, z11, zu1, zv1, lu1, lv1, rcn1, zcn1)
33
34 USE vmec_main
35 USE stel_kinds, only: rprec
36 USE vmec_params, ONLY: jmin1, jlam, ntmax, rcc, rss, zsc, zcs, m0, m1, n0
37
38 IMPLICIT NONE
39
40 REAL(rprec), DIMENSION(ns,0:ntor,0:mpol1,3*ntmax), TARGET, INTENT(inout) :: rzl_array
41 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: r11
42 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: ru1
43 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: rv1
44 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: z11
45 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: zu1
46 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: zv1
47 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: lu1
48 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: lv1
49 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: rcn1
50 REAL(rprec), DIMENSION(ns*nzeta*ntheta3,0:1), INTENT(out) :: zcn1
51
52 INTEGER :: n, m, mparity, k, i, j1, l, j1l, nsl
53 INTEGER :: ioff, joff, mj, ni, nsz
54 REAL(rprec), DIMENSION(:,:,:), POINTER :: rmncc, rmnss, zmncs, zmnsc, lmncs, lmnsc
55
56 ! WORK1 Array of inverse transforms in toroidal angle (zeta), for all radial positions
57 ! NOTE: ORDERING OF LAST INDEX IS DIFFERENT HERE THAN IN PREVIOUS VMEC2000 VERSIONS
58 REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: work1
59 REAL(rprec) :: cosmux, sinmux
60
61 ! This is the inverse Fourier transform from the (odd-m-scaled-by-1/sqrt(s))-xc
62 ! given in rzl_array to real space. Tangential derivatives are also computed.
63
64 ! WHEN COMPUTING PRECONDITIONER, USE FASTER HESSIAN VERSION (totzsps_hess) INSTEAD.
65 rmncc => rzl_array(:,:,:,rcc) ! COS(mu) COS(nv)
66 zmnsc => rzl_array(:,:,:,zsc+ntmax) ! SIN(mu) COS(nv)
67 lmnsc => rzl_array(:,:,:,zsc+2*ntmax) ! SIN(mu) COS(nv)
68 IF (lthreed) THEN
69 rmnss => rzl_array(:,:,:,rss) ! SIN(mu) SIN(nv)
70 zmncs => rzl_array(:,:,:,zcs+ntmax) ! COS(mu) SIN(nv)
71 lmncs => rzl_array(:,:,:,zcs+2*ntmax) ! COS(mu) SIN(nv)
72
73! #ifndef _HBANGLE
74 ! CONVERT FROM INTERNAL XC REPRESENTATION FOR m=1 MODES, R+(stored at rss) = .5(rss + zcs),
75 ! R-(stored at zcs) = .5(rss - zcs), TO EXTERNAL ("PHYSICAL") rss, zcs FORMS. NEED THIS EVEN
76 ! WHEN COMPUTING HESSIAN FOR FREE BOUNDARY (rmnss, zmncs at JS=NS needed in vacuum call)
77
78 CALL convert_sym (rmnss, zmncs)
79! #end /* ndef _HBANGLE */
80
81 END IF
82
83 ! v8.50: Norm for preconditioned R,Z forces: scale to boundary value only
84 ! v8.51 Restore hs dependence (1:ns, not just ns)
85 ! fnorm1 = one/SUM(rzl_array(1:ns,:,m1:,1:2*ntmax)**2)
86
87 ! ORIGIN EXTRAPOLATION (JS=1) FOR M=1 MODES
88 ! --> TODO: why? quantities on axis should not have a poloidal dependence ???
89 ! NOTE: PREVIOUS VERSIONS OF VMEC USED TWO-POINT EXTRAPOLATION FOR R,Z.
90 ! HOWEVER,THIS CAN NOT BE USED TO COMPUTE THE TRI-DIAG 2D PRECONDITIONER.
91 rzl_array(1,:,m1,:) = rzl_array(2,:,m1,:) ! now this is a constant extrapolation to the axis (????)
92
93 ioff = lbound(rmncc,2) ! starting index along n
94 joff = lbound(rmncc,3) ! starting index along m
95
96 ! ORIGIN EXTRAPOLATION OF M=0 MODES FOR LAMBDA
97 IF (lthreed .and. jlam(m0).gt.1) then
98 lmncs(1,:,m0+joff) = lmncs(2,:,m0+joff)
99 end if
100
101 nsz = ns*nzeta
102 ALLOCATE (work1(nsz,12), stat=i)
103 IF (i .ne. 0) stop 'Allocation error in VMEC2000 totzsps'
104
105 r11 = 0; ru1 = 0; rv1 = 0; rcn1 = 0
106 z11 = 0; zu1 = 0; zv1 = 0; zcn1 = 0
107 lu1 = 0; lv1 = 0
108
109 ! COMPUTE R, Z, AND LAMBDA IN REAL SPACE
110 ! NOTE: LU = d(Lam)/du, LV = -d(Lam)/dv
111 DO m = 0, mpol1
112 mparity = mod(m,2)
113 mj = m+joff
114 work1 = 0
115
116 ! m>1 contributions start at js=2, not 1 (axis)
117 ! TODO: why has the axis even an m=1 contribution?
118 j1 = jmin1(m)
119
120 ! INVERSE TRANSFORM IN N-ZETA, FOR FIXED M
121 DO n = 0, ntor
122 ni = n+ioff
123 DO k = 1, nzeta
124 l = ns*(k-1)
125
126 j1l = j1+l
127 nsl = ns+l
128
129 work1(j1l:nsl, 1) = work1(j1l:nsl, 1) + rmncc(j1:ns,ni,mj)*cosnv(k,n)
130 work1(j1l:nsl, 6) = work1(j1l:nsl, 6) + zmnsc(j1:ns,ni,mj)*cosnv(k,n)
131 work1(j1l:nsl,10) = work1(j1l:nsl,10) + lmnsc(j1:ns,ni,mj)*cosnv(k,n)
132
133 IF (.not.lthreed) cycle
134
135 work1(j1l:nsl, 4) = work1(j1l:nsl, 4) + rmnss(j1:ns,ni,mj)*cosnvn(k,n)
136 work1(j1l:nsl, 7) = work1(j1l:nsl, 7) + zmncs(j1:ns,ni,mj)*cosnvn(k,n)
137 work1(j1l:nsl,11) = work1(j1l:nsl,11) + lmncs(j1:ns,ni,mj)*cosnvn(k,n)
138
139 work1(j1l:nsl, 2) = work1(j1l:nsl, 2) + rmnss(j1:ns,ni,mj)*sinnv(k,n)
140 work1(j1l:nsl, 5) = work1(j1l:nsl, 5) + zmncs(j1:ns,ni,mj)*sinnv(k,n)
141 work1(j1l:nsl, 9) = work1(j1l:nsl, 9) + lmncs(j1:ns,ni,mj)*sinnv(k,n)
142
143 work1(j1l:nsl, 3) = work1(j1l:nsl, 3) + rmncc(j1:ns,ni,mj)*sinnvn(k,n)
144 work1(j1l:nsl, 8) = work1(j1l:nsl, 8) + zmnsc(j1:ns,ni,mj)*sinnvn(k,n)
145 work1(j1l:nsl,12) = work1(j1l:nsl,12) + lmnsc(j1:ns,ni,mj)*sinnvn(k,n)
146 END DO
147 END DO
148
149 ! INVERSE TRANSFORM IN M-THETA, FOR ALL RADIAL, ZETA VALUES
150 l = 0
151 DO i = 1, ntheta2
152
153 ! start and end indices of poloidal slice of current flux surface
154 j1l = l+1
155 nsl = l+nsz
156
157 l = l + nsz
158
159 ! basis functions for spectral constraint
160 cosmux = xmpq(m,1)*cosmu(i,m)
161 sinmux = xmpq(m,1)*sinmu(i,m)
162
163 r11(j1l:nsl,mparity) = r11(j1l:nsl,mparity) + work1(1:nsz,1)*cosmu(i,m)
164 ru1(j1l:nsl,mparity) = ru1(j1l:nsl,mparity) + work1(1:nsz,1)*sinmum(i,m)
165
166! #ifndef _HBANGLE
167 rcn1(j1l:nsl,mparity) = rcn1(j1l:nsl,mparity) + work1(1:nsz,1)*cosmux ! rmncc*cosnv*cosmux
168! #end /* ndef _HBANGLE */
169
170 z11(j1l:nsl,mparity) = z11(j1l:nsl,mparity) + work1(1:nsz,6)*sinmu(i,m)
171 zu1(j1l:nsl,mparity) = zu1(j1l:nsl,mparity) + work1(1:nsz,6)*cosmum(i,m)
172
173! #ifndef _HBANGLE
174 zcn1(j1l:nsl,mparity) = zcn1(j1l:nsl,mparity) + work1(1:nsz,6)*sinmux ! zmnsc*cosnv*sinmux
175! #end /* ndef _HBANGLE */
176
177 lu1(j1l:nsl,mparity) = lu1(j1l:nsl,mparity) + work1(1:nsz,10)*cosmum(i,m)
178
179 IF (.not.lthreed) cycle
180
181 r11(j1l:nsl,mparity) = r11(j1l:nsl,mparity) + work1(1:nsz,2)*sinmu(i,m)
182 ru1(j1l:nsl,mparity) = ru1(j1l:nsl,mparity) + work1(1:nsz,2)*cosmum(i,m)
183
184! #ifndef _HBANGLE
185 rcn1(j1l:nsl,mparity) = rcn1(j1l:nsl,mparity) + work1(1:nsz,2)*sinmux ! rmnss*sinnv*sinmux
186! #end /* ndef _HBANGLE */
187
188 rv1(j1l:nsl,mparity) = rv1(j1l:nsl,mparity) + work1(1:nsz,3)*cosmu(i,m) &
189 + work1(1:nsz,4)*sinmu(i,m)
190 z11(j1l:nsl,mparity) = z11(j1l:nsl,mparity) + work1(1:nsz,5)*cosmu(i,m)
191
192 zu1(j1l:nsl,mparity) = zu1(j1l:nsl,mparity) + work1(1:nsz,5)*sinmum(i,m)
193
194! #ifndef _HBANGLE
195 zcn1(j1l:nsl,mparity) = zcn1(j1l:nsl,mparity) + work1(1:nsz,5)*cosmux ! zmncs*sinnv*cosmux
196! #end /* ndef _HBANGLE */
197
198 zv1(j1l:nsl,mparity) = zv1(j1l:nsl,mparity) + work1(1:nsz,7)*cosmu(i,m) &
199 + work1(1:nsz,8)*sinmu(i,m)
200
201 lu1(j1l:nsl,mparity) = lu1(j1l:nsl,mparity) + work1(1:nsz,9)*sinmum(i,m)
202 lv1(j1l:nsl,mparity) = lv1(j1l:nsl,mparity) - ( work1(1:nsz,11)*cosmu(i,m) &
203 + work1(1:nsz,12)*sinmu(i,m) )
204 END DO
205
206 END DO
207
208 DEALLOCATE (work1)
209
210END SUBROUTINE totzsps
211
216SUBROUTINE convert_sym(rmnss, zmncs)
217
218 USE vmec_main
219 USE vmec_params, only: m1
220 USE stel_kinds, only: rprec
221
222 IMPLICIT NONE
223
224 REAL(rprec), DIMENSION(ns,0:ntor,0:mpol1), INTENT(inout) :: rmnss, zmncs
225
226 REAL(rprec), DIMENSION(ns,0:ntor) :: temp
227
228 ! CONVERT FROM INTERNAL REPRESENTATION TO "PHYSICAL" RMNSS, ZMNCS FOURIER FORM
229 ! rmnss = (RMNSS+ZMNCS)
230 ! zmncs = (RMNSS-ZMNCS)
231 IF (lconm1) then
232
233 ! This essentially reverts the operation performed at the end of readin().
234
235 temp(:,:) = rmnss(:,:,m1) !This is internal
236 rmnss(:,:,m1) = temp(:,:) + zmncs(:,:,m1) !Now these are physical
237 zmncs(:,:,m1) = temp(:,:) - zmncs(:,:,m1)
238
239 end if
240
241END SUBROUTINE convert_sym
242
256SUBROUTINE totzspa(rzl_array, r11, ru1, rv1, z11, zu1, zv1, lu1, lv1, rcn1, zcn1)
257
258 USE vmec_main
259 USE stel_kinds, only: rprec
260 USE vmec_params, ONLY: jmin1, jlam, ntmax, rcs, rsc, zcc, zss, m0, m1, n0
261
262 IMPLICIT NONE
263
264 REAL(rprec), DIMENSION(ns,0:ntor,0:mpol1,3*ntmax), TARGET, INTENT(inout) :: rzl_array
265 REAL(rprec), DIMENSION(ns*nzeta,ntheta3,0:1), INTENT(out) :: &
266 r11, ru1, rv1, z11, zu1, zv1, lu1, lv1, rcn1, zcn1
267
268 INTEGER :: m, n, mparity, k, i, l, j1, j1l, nsl
269 INTEGER :: ioff, joff, mj, ni
270 REAL(rprec), DIMENSION(:,:,:), POINTER :: rmncs, rmnsc, zmncc, zmnss, lmncc, lmnss
271 REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: work1
272 REAL(rprec) :: cosmux, sinmux
273
274 ! WHEN COMPUTING PRECONDITIONER, USE FASTER HESSIAN VERSION (totzsps_hess) INSTEAD.
275
276 rmnsc => rzl_array(:,:,:,rsc) !!SIN(mu) COS(nv)
277 zmncc => rzl_array(:,:,:,zcc+ntmax) !!COS(mu) COS(nv)
278 lmncc => rzl_array(:,:,:,zcc+2*ntmax) !!COS(mu) COS(nv)
279 IF (lthreed) THEN
280 rmncs => rzl_array(:,:,:,rcs) !!COS(mu) SIN(nv)
281 zmnss => rzl_array(:,:,:,zss+ntmax) !!SIN(mu) SIN(nv)
282 lmnss => rzl_array(:,:,:,zss+2*ntmax) !!SIN(mu) SIN(nv)
283 END IF
284
285 ! CONVERT FROM INTERNAL XC REPRESENTATION FOR m=1 MODES, R+(at rsc) = .5(rsc + zcc),
286 ! R-(at zcc) = .5(rsc - zcc), TO REQUIRED rsc, zcc FORMS
287
288! #ifndef _HBANGLE
289 CALL convert_asym (rmnsc, zmncc)
290! #end /* ndef _HBANGLE */
291
292 ioff = lbound(rmnsc,2)
293 joff = lbound(rmnsc,3)
294
295 ALLOCATE (work1(ns*nzeta,12), stat=i)
296 IF (i .ne. 0) stop 'Allocation error in VMEC totzspa'
297
298 ! INITIALIZATION BLOCK
299 r11 = 0; ru1 = 0; rv1 = 0
300 z11 = 0; zu1 = 0; zv1 = 0
301 lu1 = 0; lv1 = 0
302 rcn1 = 0; zcn1 = 0
303
304 IF (jlam(m0) .gt. 1) lmncc(1,:,m0+joff) = lmncc(2,:,m0+joff)
305
306 DO m = 0, mpol1
307 mparity = mod(m,2)
308 mj = m+joff
309 work1 = 0
310 j1 = jmin1(m)
311 DO n = 0, ntor
312 ni = n+ioff
313 DO k = 1, nzeta
314 l = ns*(k-1)
315 j1l = j1+l
316 nsl = ns+l
317 work1(j1l:nsl,1) = work1(j1l:nsl,1) + rmnsc(j1:ns,ni,mj)*cosnv(k,n)
318 work1(j1l:nsl,6) = work1(j1l:nsl,6) + zmncc(j1:ns,ni,mj)*cosnv(k,n)
319 work1(j1l:nsl,10) = work1(j1l:nsl,10) + lmncc(j1:ns,ni,mj)*cosnv(k,n)
320
321 IF (.not.lthreed) cycle
322
323 work1(j1l:nsl,2) = work1(j1l:nsl,2) + rmncs(j1:ns,ni,mj)*sinnv(k,n)
324 work1(j1l:nsl,3) = work1(j1l:nsl,3) + rmnsc(j1:ns,ni,mj)*sinnvn(k,n)
325 work1(j1l:nsl,4) = work1(j1l:nsl,4) + rmncs(j1:ns,ni,mj)*cosnvn(k,n)
326 work1(j1l:nsl,5) = work1(j1l:nsl,5) + zmnss(j1:ns,ni,mj)*sinnv(k,n)
327 work1(j1l:nsl,7) = work1(j1l:nsl,7) + zmnss(j1:ns,ni,mj)*cosnvn(k,n)
328 work1(j1l:nsl,8) = work1(j1l:nsl,8) + zmncc(j1:ns,ni,mj)*sinnvn(k,n)
329 work1(j1l:nsl,9) = work1(j1l:nsl,9) + lmnss(j1:ns,ni,mj)*sinnv(k,n)
330 work1(j1l:nsl,11) = work1(j1l:nsl,11) + lmnss(j1:ns,ni,mj)*cosnvn(k,n)
331 work1(j1l:nsl,12) = work1(j1l:nsl,12) + lmncc(j1:ns,ni,mj)*sinnvn(k,n)
332 END DO
333 END DO
334
335 ! INVERSE TRANSFORM IN M-THETA
336 DO i = 1, ntheta2
337
338 cosmux = xmpq(m,1)*cosmu(i,m)
339 sinmux = xmpq(m,1)*sinmu(i,m)
340
341 r11(:,i,mparity) = r11(:,i,mparity) + work1(:,1)*sinmu(i,m)
342 ru1(:,i,mparity) = ru1(:,i,mparity) + work1(:,1)*cosmum(i,m)
343 z11(:,i,mparity) = z11(:,i,mparity) + work1(:,6)*cosmu(i,m)
344 zu1(:,i,mparity) = zu1(:,i,mparity) + work1(:,6)*sinmum(i,m)
345 lu1(:,i,mparity) = lu1(:,i,mparity) + work1(:,10)*sinmum(i,m)
346
347! #ifndef _HBANGLE
348 rcn1(:,i,mparity) = rcn1(:,i,mparity) + work1(:,1)*sinmux
349 zcn1(:,i,mparity) = zcn1(:,i,mparity) + work1(:,6)*cosmux
350! #end /* ndef _HBANGLE */
351
352 IF (.not.lthreed) cycle
353
354 r11(:,i,mparity) = r11(:,i,mparity) + work1(:,2)*cosmu(i,m)
355 ru1(:,i,mparity) = ru1(:,i,mparity) + work1(:,2)*sinmum(i,m)
356 z11(:,i,mparity) = z11(:,i,mparity) + work1(:,5)*sinmu(i,m)
357 zu1(:,i,mparity) = zu1(:,i,mparity) + work1(:,5)*cosmum(i,m)
358 lu1(:,i,mparity) = lu1(:,i,mparity) + work1(:,9)*cosmum(i,m)
359
360! #ifndef _HBANGLE
361 rcn1(:,i,mparity) = rcn1(:,i,mparity) + work1(:,2)*cosmux
362 zcn1(:,i,mparity) = zcn1(:,i,mparity) + work1(:,5)*sinmux
363! #end /* ndef _HBANGLE */
364
365 rv1(:,i,mparity) = rv1(:,i,mparity) + work1(:,3)*sinmu(i,m) + work1(:,4)*cosmu(i,m)
366 zv1(:,i,mparity) = zv1(:,i,mparity) + work1(:,7)*sinmu(i,m) + work1(:,8)*cosmu(i,m)
367 lv1(:,i,mparity) = lv1(:,i,mparity) - (work1(:,11)*sinmu(i,m)+work1(:,12)*cosmu(i,m))
368 END DO
369 END DO
370
371 DEALLOCATE (work1)
372
373END SUBROUTINE totzspa
374
379SUBROUTINE convert_asym(rmnsc, zmncc)
380
381 USE vmec_main
382 USE stel_kinds, only: rprec
383
384 IMPLICIT NONE
385
386 REAL(rprec), DIMENSION(ns,0:ntor,0:mpol1), INTENT(inout) :: rmnsc, zmncc
387
388 REAL(rprec), DIMENSION(ns,0:ntor) :: temp
389
390 ! CONVERT FROM INTERNAL REPRESENTATION TO "PHYSICAL" RMNSC, ZMNCC FOURIER FORM
391 ! rmnsc(in) = .5(RMNSC+ZMNCC), zmncc(in) = .5(RMNSC-ZMNCC) -> 0
392 IF (lconm1) then
393
394 temp(:,:) = rmnsc(:,:,1) !THIS IS INTERNAL
395 rmnsc(:,:,1) = temp(:,:) + zmncc(:,:,1) !NOW THE rmnsc,zmncc ARE PHYSICAL
396 zmncc(:,:,1) = temp(:,:) - zmncc(:,:,1)
397
398 end if
399
400END SUBROUTINE convert_asym
stel_kinds
Definition stel_kinds.f90:2
stel_kinds::rprec
integer, parameter rprec
Definition stel_kinds.f90:10
vmec_main
Definition vmec_main.f90:2
vmec_main::lconm1
logical lconm1
Definition vmec_main.f90:143
vmec_main::xmpq
real(rprec), dimension(0:mpol1d, 3) xmpq
spectral condensation weighting factors
Definition vmec_main.f90:81
vmec_main::lthreed
logical lthreed
Definition vmec_main.f90:142
vmec_params
Definition vmec_params.f90:2
vmec_params::rcs
integer rcs
Definition vmec_params.f90:48
vmec_params::n0
integer, parameter n0
from totzsp
Definition vmec_params.f90:68
vmec_params::zcc
integer zcc
Definition vmec_params.f90:51
vmec_params::jlam
integer, dimension(0:mpold), parameter jlam
starting js(m) values for which Lambda is evolved
Definition vmec_params.f90:17
vmec_params::jmin1
integer, dimension(0:mpold), parameter jmin1
starting js(m) values where R,Z are non-zero
Definition vmec_params.f90:15
vmec_params::zss
integer zss
Definition vmec_params.f90:52
vmec_params::m1
integer, parameter m1
from totzsp
Definition vmec_params.f90:67
vmec_params::zsc
integer zsc
Definition vmec_params.f90:49
vmec_params::rss
integer rss
Definition vmec_params.f90:46
vmec_params::zcs
integer zcs
Definition vmec_params.f90:50
vmec_params::m0
integer, parameter m0
from totzsp
Definition vmec_params.f90:66
vmec_params::ntmax
integer ntmax
number of contributing Fourier basis function (can be 1, 2 or 4); assigned in read_indata()
Definition vmec_params.f90:43
vmec_params::rcc
integer rcc
Definition vmec_params.f90:45
vmec_params::rsc
integer rsc
Definition vmec_params.f90:47
convert_sym
subroutine convert_sym(rmnss, zmncs)
Convert from internal representation to "physical" rmnss, zmncs Fourier form.
Definition totzsp.f90:217
totzspa
subroutine totzspa(rzl_array, r11, ru1, rv1, z11, zu1, zv1, lu1, lv1, rcn1, zcn1)
Inverse-Fourier-transform anti-symmetric part of geometry from Fourier space to real space.
Definition totzsp.f90:257
totzsps
subroutine totzsps(rzl_array, r11, ru1, rv1, z11, zu1, zv1, lu1, lv1, rcn1, zcn1)
Inverse-Fourier-transform symmetric part of geometry from Fourier space to real space.
Definition totzsp.f90:33
convert_asym
subroutine convert_asym(rmnsc, zmncc)
Convert from internal representation to "physical" rmnsc, zmncc Fourier form.
Definition totzsp.f90:380