V3FIT
vacuum.f
1  SUBROUTINE vacuum_par (rmnc, rmns, zmns, zmnc, xm, xn,
2  & plascur, rbtor, wint, ns, ivac_skip,
3  & ivac, mnmax, ier_flag, lscreen)
4  USE vacmod
5  USE vparams, ONLY: nthreed, zero, one, mu0
6  USE vmec_params, ONLY: norm_term_flag, phiedge_error_flag
7  USE vmec_input, ONLY: lrfp ! JDH Added 2013-11-25, to test for RFP
8  USE vmec_main, ONLY: nznt, irst
9  USE parallel_include_module
10  USE timer_sub
11  IMPLICIT NONE
12 C-----------------------------------------------
13 C D u m m y A r g u m e n t s
14 C-----------------------------------------------
15  INTEGER :: ns, ivac_skip, ivac, mnmax, ier_flag
16  REAL(dp) :: plascur, rbtor
17  REAL(dp), DIMENSION(mnmax), INTENT(in) ::
18  & rmnc, rmns, zmns, zmnc, xm, xn
19  REAL(dp), DIMENSION(nuv3), INTENT(in) :: wint
20  LOGICAL :: lscreen
21 C-----------------------------------------------
22 C L o c a l V a r i a b l e s
23 C-----------------------------------------------
24  INTEGER :: mn, n, n1, m, i, info=0, j
25  REAL(dp), DIMENSION(:), POINTER :: potcos, potsin
26  REAL(dp), ALLOCATABLE :: potu(:), potv(:)
27  REAL(dp), ALLOCATABLE :: amatrix(:)
28  REAL(dp):: dn2, dm2, cosmn, sinmn, huv, hvv,
29  & det, bsubuvac, fac, ton, toff
30  REAL(dp) :: tmp1(2), tmp2(2)
31 C-----------------------------------------------
32 !
33 ! THIS ROUTINE COMPUTES .5 * B**2 ON THE VACUUM / PLASMA SURFACE
34 ! BASED ON THE PROGRAM BY P. MERKEL [J. Comp. Phys. 66, 83 (1986)]
35 ! AND MODIFIED BY W. I. VAN RIJ AND S. P. HIRSHMAN (1987)
36 
37 ! THE USER MUST SUPPLY THE FILE << MGRID >> WHICH INCLUDES THE MAGNETIC
38 ! FIELD DATA TO BE READ BY THE SUBROUTINE BECOIL
39 ! THE "VACUUM.INC" FILE IS DEFINED IN VMEC.UNIX
40 !
41 !
42  CALL second0(tvacon)
43  IF (.NOT.ALLOCATED(potvac)) stop 'POTVAC not ALLOCATED in VACCUM'
44 
45  ALLOCATE (amatrix(mnpd2*mnpd2), potu(nuv3), potv(nuv3), stat=i)
46  IF (i .NE. 0) stop 'Allocation error in vacuum'
47 
48  potsin => potvac(1:mnpd)
49  potcos => potvac(1+mnpd:)
50 
51  ALLOCATE (bexu(nuv3), bexv(nuv3), bexn(nuv3),
52  & bexni(nuv3), r1b(nuv), rub(nuv3), rvb(nuv3),
53  & z1b(nuv), zub(nuv3), zvb(nuv3), auu(nuv3), auv(nuv3),
54  & avv(nuv3), snr(nuv3), snv(nuv3), snz(nuv3), drv(nuv3),
55  & guu_b(nuv3), guv_b(nuv3), gvv_b(nuv3), rzb2(nuv),
56  & rcosuv(nuv), rsinuv(nuv), stat=i)
57  IF (i .NE. 0) stop 'Allocation error in vacuum'
58 
59 !
60 ! INDEX OF LOCAL VARIABLES
61 !
62 ! rmnc,rmns,zmns,zmnc: Surface Fourier coefficients (m,n) of R,Z
63 ! xm,xn: m, n values corresponding to rc,zs array
64 ! bsqvac: B**2/2 at the vacuum INTERFACE
65 ! plascur: net toroidal current
66 ! rbtor : net (effective) poloidal current (loop integrated R*Btor)
67 ! mnmax: number of R, Z modes in Fourier series of R,Z
68 ! ivac_skip: regulates whether full (=0) or incremental (>0)
69 ! update of matrix elements is necessary
70 !
71 !
72 ! compute and store mean magnetic fields (due to
73 ! toroidal plasma current and EXTERNAL tf-coils)
74 ! note: these are fixed for a constant current iteration
75 !
76 ! bfield = rbtor*grad(zeta) + plascur*grad("theta") - grad(potential)
77 !
78 ! where "theta" is computed using Biot-Savart law for filaments
79 ! Here, the potential term is needed to satisfy B dot dS = 0 and has the form:
80 !
81 ! potential = SUM potsin*SIN(mu - nv) + potcos*COS(mu - nv)
82 !
83  CALL second0(ton)
84  IF (.NOT. ALLOCATED(tanu)) CALL precal (wint)
85  CALL surface (rmnc, rmns, zmns, zmnc, xm, xn, mnmax)
86  CALL second0(toff)
87  timer_vac(tsurf) = timer_vac(tsurf) + (toff-ton)
88 
89  ton = toff
90  CALL bextern (plascur, wint, lscreen)
91  CALL second0(toff)
92  timer_vac(tbext) = timer_vac(tbext)+(toff-ton)
93 !
94 ! Determine scalar magnetic potential POTVAC
95 !
96  CALL second0(ton)
97  CALL scalpot (potvac, amatrix, ivac_skip)
98  CALL second0(toff)
99  timer_vac(tscal) = timer_vac(tscal) + (toff-ton)
100 
101  ton = toff
102  CALL solver (amatrix, potvac, mnpd2, 1, info)
103  CALL second0(toff)
104  timer_vac(tsolver) = timer_vac(tsolver) + (toff-ton)
105  solver_time = solver_time + (toff - ton)
106 
107  IF (info .NE. 0) stop 'Error in solver in VACUUM'
108 !
109 ! compute tangential covariant (sub u,v) and contravariant
110 ! (super u,v) magnetic field components on the plasma surface
111 !
112  potu(nuv3min:nuv3max) = 0; potv(nuv3min:nuv3max) = 0
113 
114  mn = 0
115  DO n = -nf, nf
116  dn2 = -(n*nfper)
117  n1 = abs(n)
118  DO m = 0, mf
119  mn = mn + 1
120  dm2 = m
121  j = 0
122  DO i = nuv3min, nuv3max
123  j = j + 1
124  cosmn = potsin(mn)*cosmni(j,mn)/(pi2*pi2*wint(i))
125  potu(i) = potu(i) + dm2*cosmn
126  potv(i) = potv(i) + dn2*cosmn
127  END DO
128  IF (.NOT.lasym) cycle
129  j = 0
130  DO i = nuv3min, nuv3max
131  j = j + 1
132  sinmn = potcos(mn)*sinmni(j,mn)/(pi2*pi2*wint(i))
133  potu(i) = potu(i) - dm2*sinmn
134  potv(i) = potv(i) - dn2*sinmn
135  END DO
136  END DO
137  END DO
138 
139  DO i = nuv3min, nuv3max
140  bsubu_sur(i) = potu(i) + bexu(i) !Covariant components
141  bsubv_sur(i) = potv(i) + bexv(i)
142  huv = p5*guv_b(i)*(nfper)
143  hvv = gvv_b(i)*(nfper*nfper)
144  det = one/(guu_b(i)*hvv-huv*huv)
145  bsupu_sur(i) = (hvv*bsubu_sur(i)-huv*bsubv_sur(i))*det !Contravariant components
146  bsupv_sur(i) = ((-huv*bsubu_sur(i))+guu_b(i)*bsubv_sur(i))*det
147  bsqvac(i) = p5*(bsubu_sur(i)*bsupu_sur(i)
148  & + bsubv_sur(i)*bsupv_sur(i)) !.5*|Bvac|**2
149  brv(i) = rub(i)*bsupu_sur(i) + rvb(i)*bsupv_sur(i)
150  bphiv(i) = r1b(i)*bsupv_sur(i)
151  bzv(i) = zub(i)*bsupu_sur(i) + zvb(i)*bsupv_sur(i)
152  END DO
153 !
154 ! PRINT OUT VACUUM PARAMETERS
155 !
156  IF (ivac .EQ. 0) THEN
157  ivac = ivac + 1
158  IF (vrank .EQ. 0) THEN
159  IF (lscreen) WRITE (*, 200) nfper, mf, nf, nu, nv
160  WRITE (nthreed, 200) nfper, mf, nf, nu, nv
161  END IF
162  200 FORMAT(/,2x,'In VACUUM, np =',i3,2x,'mf =',i3,2x,'nf =',i3,
163  & ' nu =',i3,2x,'nv = ',i4)
164 
165  bsubuvac = 0
166  bsubvvac = 0
167  DO i = nuv3min, nuv3max
168  bsubuvac = bsubuvac + bsubu_sur(i)*wint(i)
169  bsubvvac = bsubvvac + bsubv_sur(i)*wint(i)
170  END DO
171  tmp1(1) = bsubuvac
172  tmp1(2)=bsubvvac
173  CALL second0(ton)
174 
175  IF (vlactive) THEN
176  CALL mpi_allreduce(tmp1, tmp2, 2, mpi_real8, mpi_sum,
177  & vac_comm, mpi_err)
178  END IF
179 
180  CALL second0(toff)
181  allreduce_time = allreduce_time + (toff - ton)
182  bsubuvac = tmp2(1); bsubvvac = tmp2(2)
183  bsubuvac = bsubuvac*signgs*pi2
184 
185  fac = 1.e-6_dp/mu0
186  IF (vrank .EQ. 0) THEN
187  IF (lscreen ) THEN
188  WRITE (*,1000) bsubuvac*fac, plascur*fac, bsubvvac, rbtor
189  END IF
190  WRITE (nthreed, 1000) bsubuvac*fac, plascur*fac, bsubvvac,
191  & rbtor
192  END IF
193  1000 FORMAT(2x,'2*pi * a * -BPOL(vac) = ',1p,e10.2,
194  & ' TOROIDAL CURRENT = ',e10.2,/,2x,'R * BTOR(vac) = ',
195  & e10.2,' R * BTOR(plasma) = ',e10.2)
196 ! JDH Add test for RFP. 2013-11-25
197 ! IF (rbtor*bsubvvac .lt. zero) ier_flag = phiedge_error_flag
198 ! IF (ABS((plascur - bsubuvac)/rbtor) .gt. 1.e-2_dp)
199 ! 1 ier_flag = 10
200  IF (rbtor*bsubvvac .LT. zero) THEN
201  IF (lrfp) THEN
202  IF (vrank .EQ. 0) THEN
203  IF (lscreen) WRITE(*,1100)
204  WRITE(nthreed,1100)
205  END IF
206  ELSE
207  ier_flag = phiedge_error_flag
208  ENDIF
209  ENDIF
210  IF (abs((plascur - bsubuvac)/rbtor) .GT. 5.e-2_dp) THEN
211  IF (lrfp) THEN
212  IF (vrank .EQ. 0) THEN
213  IF (lscreen) WRITE(*,1200)
214  WRITE(nthreed,1200)
215  END IF
216  ELSE
217  ier_flag = 10
218  ENDIF
219  ENDIF
220 ! END JDH Add test for RFP. 2013-11-25
221  ENDIF
222 1100 FORMAT('lrfp is TRUE. Ignore phiedge sign problem')
223 1200 FORMAT('lrfp is TRUE. Proceed with convergence')
224 
225  IF (ALLOCATED(bexu))
226  & DEALLOCATE (bexu, bexv, bexn, bexni, r1b, rub, rvb, z1b, zub,
227  & zvb, auu, auv, avv, snr, snv, snz, drv, guu_b, guv_b, gvv_b,
228  & rzb2, rcosuv, rsinuv, stat=i)
229  IF (i .NE. 0) stop 'Deallocation error in vacuum'
230 
231  DEALLOCATE (amatrix, potu, potv, stat=i)
232  IF (i .NE. 0) stop 'Deallocation error in vacuum'
233 
234  CALL second0(ton)
235 
236  IF (vlactive) THEN
237  CALL mpi_allgatherv(mpi_in_place, numjs_vac, mpi_real8, bsqvac,
238  & counts_vac, disps_vac, mpi_real8, vac_comm,
239  & mpi_err)
240  END IF
241 
242  CALL second0(toff)
243  timer_vac(tallgv) = timer_vac(tallgv) + (toff-ton)
244 
245  tvacoff = toff
246  vacuum_time = vacuum_time + (tvacoff - tvacon)
247 
248  END SUBROUTINE vacuum_par
249 
250 
251  SUBROUTINE vacuum(rmnc, rmns, zmns, zmnc, xm, xn,
252  & plascur, rbtor, wint, ns, ivac_skip, ivac,
253  & mnmax, ier_flag, lscreen)
254  USE vacmod
255  USE vparams, ONLY: nthreed, zero, one, mu0
256  USE vmec_params, ONLY: norm_term_flag, phiedge_error_flag
257  USE vmec_input, ONLY: lrfp ! JDH Added 2013-11-25, to test for RFP
258 
259  IMPLICIT NONE
260 C-----------------------------------------------
261 C D u m m y A r g u m e n t s
262 C-----------------------------------------------
263  INTEGER :: ns, ivac_skip, ivac, mnmax, ier_flag
264  REAL(dp) :: plascur, rbtor
265  REAL(dp), DIMENSION(mnmax), INTENT(in) ::
266  & rmnc, rmns, zmns, zmnc, xm, xn
267  REAL(dp), DIMENSION(ns, nuv3), INTENT(in) :: wint
268  LOGICAL :: lscreen
269 C-----------------------------------------------
270 C L o c a l V a r i a b l e s
271 C-----------------------------------------------
272  INTEGER :: i
273  REAL(dp), ALLOCATABLE :: tmpwint(:)
274 C-----------------------------------------------
275  ALLOCATE (tmpwint(nuv3), stat=i)
276  IF (i .NE. 0) stop 'Allocation error in vacuum'
277 
278  tmpwint(:) = wint(ns, :)
279 
280  CALL vacuum_par (rmnc, rmns, zmns, zmnc, xm, xn,
281  & plascur, rbtor, tmpwint, ns, ivac_skip,
282  & ivac, mnmax, ier_flag, lscreen)
283 
284  DEALLOCATE (tmpwint, stat = i)
285 
286  END SUBROUTINE vacuum
Generated on Thu Mar 5 2020 15:49:27 for V3FIT by   doxygen 1.8.17