V3FIT
funct3d.f
1  SUBROUTINE funct3d_par (lscreen, ier_flag)
2  USE vmec_main
3  USE vacmod, ONLY: bsqvac, bsqvac0, raxis_nestor, zaxis_nestor,
4  & nuv, nuv3
5  USE vmec_params, ONLY: ntmax, norm_term_flag
6  USE realspace
7  USE vforces
8  USE xstuff
9  USE timer_sub
10  USE precon2d, ONLY: ictrl_prec2d, lhess_exact, l_edge
11  USE vparams, ONLY: twopi
12  USE totzsp_mod
13  USE tomnsp_mod
14  USE timer_sub
15  USE parallel_include_module
16  USE parallel_vmec_module, ONLY: saxlastntype, zerolastntype,
17  & saxpbylastntype
18  USE blocktridiagonalsolver, ONLY: l_colscale
19 
20  IMPLICIT NONE
21 C-----------------------------------------------
22 C D u m m y A r g u m e n t s
23 C-----------------------------------------------
24  INTEGER, INTENT(inout) :: ier_flag
25  LOGICAL, INTENT(in) :: lscreen
26 C-----------------------------------------------
27 C L o c a l V a r i a b l e s
28 C-----------------------------------------------
29  INTEGER :: l0pi, l, ivacskip
30  INTEGER :: nvskip0 = 0
31  REAL(dp), DIMENSION(mnmax) ::
32  & rmnc, zmns, lmns, rmns, zmnc, lmnc
33  REAL(dp), DIMENSION(:,:,:), POINTER :: lu, lv
34  REAL(dp) :: presf_ns, delr_mse, delt0
35  REAL(dp) :: tbroadon, tbroadoff
36  REAL(dp), EXTERNAL :: pmass
37  INTEGER :: i, j, k, nsmin, nsmax, m
38  REAL(dp), ALLOCATABLE, DIMENSION(:) :: bcastbuf
39  INTEGER, DIMENSION(4) :: bbuf
40 C-----------------------------------------------
41  CALL second0 (tfunon)
42 !
43 ! POINTER ALIASES
44 !
45 
46  nfunct3d = nfunct3d + 1
47  lu => pczmn; lv => pcrmn
48 
49 
50 ! CONVERT ODD M TO 1/SQRT(S) INTERNAL REPRESENTATION
51  active1: IF (lactive) THEN
52  IF (ictrl_prec2d .EQ. 3) THEN
53  CALL saxpbylastntype(one, pxc, one, pxcdot, pgc)
54  CALL saxlastntype(pgc, pscalxc, pgc)
55  ELSE IF (ictrl_prec2d.EQ.1 .AND. l_colscale) THEN
56  pgc = (pxc-pxsave)*pcol_scale + pxsave
57  CALL saxlastntype(pgc, pscalxc, pgc)
58  ELSE
59  CALL saxlastntype(pxc, pscalxc, pgc)
60  END IF
61 
62 ! RIGID BODY SHIFT OF RMNCC(JS.GT.1,0,0) BY DELR_MSE= R00-RAXMSE
63 !
64 
65 ! INVERSE FOURIER TRANSFORM TO S,THETA,ZETA SPACE
66 ! R, Z, AND LAMBDA ARRAYS IN FOURIER SPACE
67 ! FIRST, DO SYMMETRIC [ F(u,v) = F(-u,-v) ] PIECES
68 ! ON THE RANGE u = 0,pi and v = 0,2*pi
69 !
70 
71  CALL totzsps_par (pgc, pr1, pru, prv, pz1, pzu, pzv, lu, lv,
72  & prcon, pzcon, ier_flag)
73 
74 !
75 ! ANTI-SYMMETRIC CONTRIBUTIONS TO INVERSE TRANSFORMS
76 !
77  IF (lasym) THEN
78  CALL totzspa_par (pgc, parmn, pbrmn, pextra3, pazmn, pbzmn,
79  & pextra4, pblmn, pclmn, pextra1, pextra2)
80 
81 ! SUM SYMMETRIC, ANTISYMMETRIC PIECES APPROPRIATELY
82 ! TO GET R, Z, L, (AND RCON, ZCON) ON FULL RANGE OF u (0 to 2*pi)
83  CALL symrzl_par (pr1, pru, prv, pz1, pzu, pzv, lu, lv,
84  & prcon, pzcon, parmn, pbrmn, pextra3, pazmn,
85  & pbzmn, pextra4, pblmn, pclmn, pextra1,
86  & pextra2)
87  END IF
88 
89 ! l0pi = ns*(1 + nzeta*(ntheta2 - 1)) !u = pi, v = 0, js = ns
90 ! router = r1(ns,0) + r1(ns,1)
91 ! rinner = r1(l0pi,0) + r1(l0pi,1)
92  r00 = pr1(1,1,0)
93  z00 = pz1(1,1,0)
94 
95 
96 !
97 ! COMPUTE CONSTRAINT RCON, ZCON
98 !
99  nsmin=tlglob; nsmax=trglob
100  DO l = nsmin, nsmax
101  prcon(:,l,0) = prcon(:,l,0) + prcon(:,l,1)*psqrts(:,l)
102  pzcon(:,l,0) = pzcon(:,l,0) + pzcon(:,l,1)*psqrts(:,l)
103  pru0(:,l) = pru(:,l,0) + pru(:,l,1)*psqrts(:,l)
104  pzu0(:,l) = pzu(:,l,0) + pzu(:,l,1)*psqrts(:,l)
105  END DO
106 
107 ! COMPUTE RCON0, ZCON0 FOR FIXED BOUNDARY BY SCALING EDGE VALUES
108 ! SCALE BY POWER OF SQRTS, RATHER THAN USE rcon0 = rcon, etc. THIS
109 ! PREVENTS A DISCONTINUITY WHEN RESTARTING FIXED BOUNDARY WITH NEW RCON0....
110 !
111 ! NOTE: IN ORDER TO MAKE INITIAL CONSTRAINT FORCES SAME FOR FREE/FIXED
112 ! BOUNDARY, WE SET RCON0,ZCON0 THE SAME INITIALLY, BUT TURN THEM OFF
113 ! SLOWLY IN FREE-BOUNDARY VACUUM LOOP (BELOW)
114 !
115  IF (ictrl_prec2d .EQ. 2) THEN
116  DO l = nsmin, nsmax
117  prcon0(:,l) = prcon(:,l,0)
118  pzcon0(:,l) = pzcon(:,l,0)
119  END DO
120  ELSE IF (iter2 .EQ. iter1 .AND.
121  & ivac .LE. 0 .AND.
122  & ictrl_prec2d .EQ. 0) THEN
123 #if defined(MPI_OPT)
124  ALLOCATE(bcastbuf(2*nznt))
125  bcastbuf(1:nznt)=prcon(:,ns,0)
126  bcastbuf(nznt+1:2*nznt)=pzcon(:,ns,0)
127  CALL second0(tbroadon)
128  CALL mpi_bcast(bcastbuf,SIZE(bcastbuf),mpi_real8,nranks-1,
129  & ns_comm,mpi_err)
130  CALL second0(tbroadoff)
131  broadcast_time = broadcast_time + (tbroadoff-tbroadon)
132  prcon(:,ns,0)=bcastbuf(1:nznt)
133  pzcon(:,ns,0)=bcastbuf(nznt+1:2*nznt)
134  DEALLOCATE(bcastbuf)
135 #endif
136  DO l = nsmin, nsmax
137  prcon0(:,l) = prcon(:,ns,0)*psqrts(:,l)**2
138  pzcon0(:,l) = pzcon(:,ns,0)*psqrts(:,l)**2
139  END DO
140  END IF
141 
142 !
143 ! COMPUTE S AND THETA DERIVATIVE OF R AND Z AND JACOBIAN ON HALF-GRID
144 !
145  CALL jacobian_par
146 
147 ! COMPUTE COVARIANT COMPONENTS OF B, MAGNETIC AND KINETIC
148 ! PRESSURE, AND METRIC ELEMENTS ON HALF-GRID
149 
150  CALL second0(tbcovon)
151  CALL bcovar_par(lu, lv, pxc, ier_flag)
152  CALL second0(tbcovoff)
153  bcovar_time=bcovar_time+(tbcovoff - tbcovon)
154 
155  END IF active1
156 
157  CALL mpi_bcast( ier_flag, 1, mpi_integer, 0,
158  & runvmec_comm_world, mpi_err) !SAL 070719
159 
160 #if defined(MPI_OPT)
161  bbuf(1)=irst; bbuf(2)=iequi; bbuf(3)=ivac; bbuf(4)=iter2
162  CALL mpi_bcast(bbuf,4,mpi_integer,0,runvmec_comm_world,mpi_err)
163  irst=bbuf(1); iequi=bbuf(2); ivac=bbuf(3); iter2=bbuf(4)
164  CALL mpi_bcast(lfreeb,1,mpi_logical,0,runvmec_comm_world,mpi_err)
165  IF (ier_flag .ne. norm_term_flag) RETURN !SAL 070719
166 #endif
167 
168  IF (irst.EQ.2 .AND. iequi.EQ.0) THEN
169  CALL zerolastntype(pgc)
170  GOTO 100
171  END IF
172 
173  timer(tbcov) = timer(tbcov) + (tbcovoff - tbcovon)
174 
175 ! COMPUTE VACUUM MAGNETIC PRESSURE AT PLASMA EDGE
176 ! NOTE: FOR FREE BOUNDARY RUNS, THE VALUE OF RBTOR=R*BTOR
177 ! AT THE PLASMA EDGE SHOULD BE ADJUSTED TO APPROXIMATELY
178 ! EQUAL THE VACUUM VALUE. THIS CAN BE DONE BY CHANGING
179 ! EITHER PHIEDGE OR THE INITIAL CROSS SECTION ACCORDING
180 ! TO THE SCALING LAW R*BTOR .EQ. PHIEDGE/(R1 * Z1).
181 
182  IF (lfreeb .AND.
183  & iter2 .GT. 1 .AND.
184  & iequi .EQ. 0) THEN
185 
186  IF (ictrl_prec2d.LE.1 .AND. (fsqr + fsqz).LE.1.e-3_dp)
187  & ivac = ivac+1 !decreased from e-1 to e-3 - sph12/04
188 
189  IF (nvskip0 .EQ. 0) nvskip0 = max(1, nvacskip)
190 
191  ivac0: IF (ivac .GE. 0) THEN
192 !SPH OFF: 6.20.17
193 ! IF INITIALLY ON, TURN OFF rcon0, zcon0 SLOWLY
194  IF (lactive) THEN
195  IF (ictrl_prec2d .EQ. 2) THEN
196  prcon0(:,nsmin:nsmax) = 0; pzcon0(:,nsmin:nsmax) = 0
197  ELSE IF (ictrl_prec2d .EQ. 0) THEN
198  prcon0(:,nsmin:nsmax) = 0.9_dp*prcon0(:,nsmin:nsmax)
199  pzcon0(:,nsmin:nsmax) = 0.9_dp*pzcon0(:,nsmin:nsmax)
200  END IF
201  ENDIF
202  CALL second0 (tvacon)
203  ivacskip = mod(iter2-iter1,nvacskip)
204  IF (ivac .LE. 2) ivacskip = 0
205 
206 ! EXTEND NVACSKIP AS EQUILIBRIUM CONVERGES
207  IF (ivacskip .EQ. 0) THEN
208  nvacskip = one/max(1.e-1_dp, 1.e11_dp*(fsqr+fsqz))
209  nvacskip = max(nvacskip, nvskip0)
210  END IF
211 
212 !
213 ! NORMALLY, WHEN COMPUTING THE HESSIAN, IT IS SUFFICIENT TO
214 ! COMPUTE THE VARIATIONS IN THE "EXACT" SOLUTION, NOT THE ENTIRE
215 ! FIELD PERIOD SUM. THUS, FOR ictrl_prec2d >= 2, SET ivacskip = 1
216 ! FOR ictrl_prec2d = 1 (RUN WITH PRECONDITIONER APPLIED), MUST
217 ! COMPUTE EXACT VACUUM RESPONSE NOW.
218 !
219 ! THE EXCEPTION TO THIS IS IF WE ARE TESTING THE HESSIAN (lHess_exact=T),
220 ! THEN MUST USE FULL VACUUM CALCULATION TO COMPUTE IT (ivacskip=0)
221 !
222 ! lHess_exact = .FALSE.
223 
224  IF (ictrl_prec2d .NE. 0) THEN
225  IF (lhess_exact .OR. ictrl_prec2d.EQ.2) THEN !Accurate Hessian
226  ivacskip = 0
227  ELSE
228  ivacskip = 1 !Fast vacuum calculation used to compute Hessian
229  ENDIF
230  ENDIF
231 
232 ! NOTE: pgc contains correct edge values of r,z,l arrays
233 ! convert_sym, convert_asym have been applied to m=1 modes
234 
235  CALL convert_par(rmnc,zmns,lmns,rmns,zmnc,lmnc,pgc)
236 
237 ! DO NOT UPDATE THIS WHEN USING PRECONDITIONER: BREAKS TRI-DIAGONAL STRUCTURE
238  IF (ictrl_prec2d.EQ.0 .OR. ictrl_prec2d.EQ.2) THEN
239  raxis_nestor(1:nzeta) = pr1(1:nzeta,1,0)
240  zaxis_nestor(1:nzeta) = pz1(1:nzeta,1,0)
241 
242 #if defined(MPI_OPT)
243  ALLOCATE (bcastbuf(2*nzeta))
244  bcastbuf(1:nzeta) = raxis_nestor(1:nzeta)
245  bcastbuf(nzeta+1:2*nzeta) = zaxis_nestor(1:nzeta)
246  CALL second0(tbroadon)
247  CALL mpi_bcast(bcastbuf,SIZE(bcastbuf),mpi_real8,0,
248  & runvmec_comm_world,mpi_err)
249  CALL second0(tbroadoff)
250  broadcast_time = broadcast_time + (tbroadoff - tbroadon)
251  raxis_nestor(1:nzeta) = bcastbuf(1:nzeta)
252  zaxis_nestor(1:nzeta) = bcastbuf(nzeta+1:2*nzeta)
253  DEALLOCATE (bcastbuf)
254 #endif
255  END IF
256 
257 #if defined(MPI_OPT)
258  ALLOCATE (bcastbuf(2))
259  bcastbuf(1)=rbtor
260  bcastbuf(2)=ctor
261  IF (lactive) THEN
262  CALL second0(tbroadon)
263  CALL mpi_bcast(bcastbuf,SIZE(bcastbuf),mpi_real8,
264  & nranks-1,ns_comm,mpi_err)
265  CALL second0(tbroadoff)
266  broadcast_time = broadcast_time + (tbroadoff -tbroadon)
267  END IF
268 
269  CALL second0(tbroadon)
270  IF (vlactive) THEN
271  CALL mpi_bcast(bcastbuf,SIZE(bcastbuf),mpi_real8,0,
272  & vac_comm,mpi_err)
273  END IF
274  CALL second0(tbroadoff)
275  broadcast_time = broadcast_time + (tbroadoff -tbroadon)
276  rbtor=bcastbuf(1)
277  ctor=bcastbuf(2)
278  DEALLOCATE (bcastbuf)
279 #endif
280 
281  IF (vlactive) THEN
282  IF (ictrl_prec2d .NE. 3 .OR.
283  & l_edge) THEN
284  CALL vacuum_par (rmnc, rmns, zmns, zmnc, xm, xn,
285  & ctor, rbtor, pwint_ns, ns, ivacskip,
286  & ivac, mnmax, ier_flag, lscreen)
287  IF (ictrl_prec2d .EQ. 2) bsqvac0 = bsqvac
288  ELSE
289  bsqvac = bsqvac0; ier_flag = 0
290  END IF
291  END IF
292 
293  IF (vnranks .LT. nranks) THEN
294 #if defined(MPI_OPT)
295  CALL mpi_bcast(bsqvac,SIZE(bsqvac),mpi_real8,0,
296  & ns_comm,mpi_err)
297 #endif
298  END IF
299 
300  IF (ier_flag .NE. 0) THEN
301  RETURN
302  END IF
303 !
304 ! RESET FIRST TIME FOR SOFT START
305 !
306  IF (ivac .EQ. 1) THEN
307  irst = 2; delt0 = delt
308  CALL restart_iter(delt0)
309  irst = 1
310  END IF
311 
312 !
313 ! IN CASE PRESSURE IS NOT ZERO AT EXTRAPOLATED EDGE...
314 ! UNCOMMENT ALL "RPRES" COMMENTS HERE AND IN BCOVAR, FORCES ROUTINES
315 ! IF NON-VARIATIONAL FORCES ARE DESIRED
316 !
317 ! presf_ns = 1.5_dp*pres(ns) - 0.5_dp*pres(ns1)
318 ! MUST NOT BREAK TRI-DIAGONAL RADIAL COUPLING: OFFENDS PRECONDITIONER!
319  presf_ns = pmass(hs*(ns-1.5_dp))
320  IF (presf_ns .NE. zero) THEN
321  presf_ns = (pmass(1._dp)/presf_ns) * pres(ns)
322  END IF
323 
324  DO l = 1, nznt
325  bsqsav(l,3) = 1.5_dp*pbzmn_o(l,ns)
326  & - 0.5_dp*pbzmn_o(l,ns-1)
327  pgcon(l,ns) = bsqvac(l) + presf_ns
328  rbsq(l) = pgcon(l,ns)*(pr1(l,ns,0) + pr1(l,ns,1))*ohs
329  dbsq(l) = abs(pgcon(l,ns)-bsqsav(l,3))
330  END DO
331 
332  IF (ivac .EQ. 1) THEN
333  IF (vlactive) THEN
334  bsqsav(:nznt,1) = pbzmn_o(:,ns)
335  bsqsav(:nznt,2) = bsqvac(:nznt)
336 #if defined(MPI_OPT)
337  CALL mpi_bcast(bsqsav(:,1),nznt,mpi_real8,
338  & nranks-1,ns_comm,mpi_err)
339 #endif
340  END IF
341  ELSE IF (ictrl_prec2d .NE. 3) THEN
342 #if defined(MPI_OPT)
343  CALL mpi_bcast(bsqsav(:,1),nznt,mpi_real8,
344  & 0,ns_comm,mpi_err)
345 #endif
346  END IF
347 
348  CALL second0 (tvacoff)
349  timer(tvac) = timer(tvac) + (tvacoff - tvacon)
350  IF (ictrl_prec2d .GE. 2) THEN
351  timer(tvac_2d) = timer(tvac_2d)+ (tvacoff - tvacon)
352  END IF
353  END IF ivac0
354  END IF
355 !
356 ! COMPUTE CONSTRAINT FORCE
357 !
358  active2: IF (lactive) THEN
359  IF (iequi .NE. 1) THEN
360  DO l = nsmin, nsmax
361  pextra1(:,l,0) = (prcon(:,l,0) - prcon0(:,l))*pru0(:,l)
362  & + (pzcon(:,l,0) - pzcon0(:,l))*pzu0(:,l)
363  END DO
364  CALL alias_par (pgcon, pextra1(:,:,0), pgc, pgc(1+mns),
365  & pgc(1+2*mns), pextra1(:,:,1))
366  ELSE
367  IF (lrecon) THEN
368  pxc(:ns) = pxc(:ns) + delr_mse
369  END IF
370  GOTO 100
371  END IF
372 
373 !
374 ! COMPUTE MHD FORCES ON INTEGER-MESH
375 !
376  CALL forces_par
377 
378 ! SYMMETRIZE FORCES (in u-v space)
379 !
380 
381  IF (lasym) THEN
382  CALL symforce_par (parmn, pbrmn, pcrmn, pazmn, pbzmn,
383  & pczmn, pblmn, pclmn, prcon, pzcon, pr1,
384  & pru, prv, pz1, pzu, pzv, pextra3,
385  & pextra4, pextra1, pextra2)
386  END IF
387 
388 !
389 ! FOURIER-TRANSFORM MHD FORCES TO (M,N)-SPACE
390 !
391 
392  CALL tomnsps_par (pgc, parmn, pbrmn, pcrmn, pazmn, pbzmn,
393  & pczmn, pblmn, pclmn, prcon, pzcon)
394 
395  IF (lasym) THEN
396  CALL tomnspa_par (pgc, pr1, pru, prv, pz1, pzu, pzv,
397  & pextra3, pextra4, pextra1, pextra2)
398  END IF
399 
400 !================================================================
401 !
402 ! COMPUTE FORCE RESIDUALS (RAW AND PRECONDITIONED)
403 !
404 !================================================================
405  CALL second0 (treson)
406 
407  CALL saxlastntype(pgc, pscalxc, pgc)
408 
409  CALL residue_par(pgc, pgc(1+irzloff), pgc(1+2*irzloff))
410 
411  END IF active2
412 
413 !NEED THIS ON ALL PROCESSORS IN GROUP (NOT JUST ACTIVE ONES) FOR STOPPING CRITERION IN EVOLVE
414 #if defined(MPI_OPT)
415  IF (gnranks .GT. nranks) THEN
416  ALLOCATE(bcastbuf(6))
417  bcastbuf(1) = fsqr; bcastbuf(2) = fsqr1
418  bcastbuf(3) = fsqz; bcastbuf(4) = fsqz1
419  bcastbuf(5) = fsql; bcastbuf(6) = fsql1
420  CALL second0(tbroadon)
421  CALL mpi_bcast(bcastbuf,SIZE(bcastbuf),mpi_real8,0,
422  & runvmec_comm_world,mpi_err)
423  CALL second0(tbroadoff)
424  broadcast_time = broadcast_time + (tbroadoff -tbroadon)
425  fsqr = bcastbuf(1); fsqr1 = bcastbuf(2)
426  fsqz = bcastbuf(3); fsqz1 = bcastbuf(4)
427  fsql = bcastbuf(5); fsql1 = bcastbuf(6)
428  DEALLOCATE(bcastbuf)
429  END IF
430 #endif
431 
432 ! Force new initial axis guess IF ALLOWED (l_moveaxis=T)
433  IF (lmove_axis .and.
434  & iter2 .eq .1 .and.
435  & (fsqr + fsqz +fsql) .gt. 1.e2_dp) THEN
436  irst = 4
437  END IF
438 
439  CALL second0 (tresoff)
440  timer(tres) = timer(tres) + (tresoff - treson)
441 
442  100 CONTINUE
443 
444  CALL second0 (tfunoff)
445  timer(tfun) = timer(tfun) + (tfunoff - tfunon)
446  IF (ictrl_prec2d .GE. 2) THEN
447  timer(tfun_2d) = timer(tfun_2d) + (tfunoff - tfunon)
448  END IF
449 
450  END SUBROUTINE funct3d_par
451 
452  SUBROUTINE funct3d (lscreen, ier_flag)
453  USE vmec_main
454 #ifdef _VACUUM2
455  USE vac2_vacmod, ONLY: mf, nf, bsqvac
456 #else
457  USE vacmod, ONLY: bsqvac, bsqvac0, raxis_nestor, zaxis_nestor
458 #endif
459  USE vmec_params, ONLY: ntmax
460  USE realspace
461  USE vforces
462  USE xstuff
463  USE timer_sub
464  USE precon2d, ONLY: ictrl_prec2d, lhess_exact, l_edge
465  USE vparams, ONLY: twopi
466  USE totzsp_mod
467  USE tomnsp_mod
468  USE timer_sub
469 #ifdef _HBANGLE
470  USE angle_constraints, ONLY: getrz, getfrho, xtempa
471 #endif
472  IMPLICIT NONE
473 !-----------------------------------------------
474 ! D u m m y A r g u m e n t s
475 !-----------------------------------------------
476  INTEGER, INTENT(inout) :: ier_flag
477  LOGICAL, INTENT(in) :: lscreen
478 !-----------------------------------------------
479 ! L o c a l V a r i a b l e s
480 !-----------------------------------------------
481  INTEGER :: l0pi, l, lk, ivacskip
482  INTEGER :: nvskip0 = 0
483  REAL(dp), DIMENSION(mnmax) ::
484  1 rmnc, zmns, lmns, rmns, zmnc, lmnc
485  REAL(dp), DIMENSION(:), POINTER :: lu, lv
486  REAL(dp) :: presf_ns, delr_mse, delt0
487  REAL(dp), EXTERNAL :: pmass
488 !-----------------------------------------------
489 !
490 ! POINTER ALIASES
491 !
492  lu => czmn; lv => crmn
493 
494  CALL second0 (tfunon)
495 
496 ! CONVERT ODD M TO 1/SQRT(S) INTERNAL REPRESENTATION
497 !
498 #ifdef _HBANGLE
499 !Overwrites rzl_array, but that is OK since gc = rzl_array in CALL, xc preserved
500  xtempa = xc
501  CALL getrz(xc)
502 #endif
503  IF (ictrl_prec2d .EQ. 3) THEN
504  gc(:neqs) = scalxc(:neqs)*(xc(:neqs)+xcdot(:neqs))
505  ELSE
506  gc(:neqs) = scalxc(:neqs)*xc(:neqs)
507  END IF
508 
509 !
510 ! INVERSE FOURIER TRANSFORM TO S,THETA,ZETA SPACE
511 ! R, Z, AND LAMBDA ARRAYS IN FOURIER SPACE
512 ! FIRST, DO SYMMETRIC [ F(u,v) = F(-u,-v) ] PIECES
513 ! ON THE RANGE u = 0,pi and v = 0,2*pi
514 !
515 
516  CALL totzsps (gc, r1, ru, rv, z1, zu, zv, lu, lv, rcon, zcon)
517 !
518 ! ANTI-SYMMETRIC CONTRIBUTIONS TO INVERSE TRANSFORMS
519 !
520  IF (lasym) THEN
521  CALL totzspa (gc, armn, brmn, extra3, azmn, bzmn, extra4,
522  & blmn, clmn, extra1, extra2)
523 
524 ! SUM SYMMETRIC, ANTISYMMETRIC PIECES APPROPRIATELY
525 ! TO GET R, Z, L, (AND RCON, ZCON) ON FULL RANGE OF u (0 to 2*pi)
526  CALL symrzl (r1, ru, rv, z1, zu, zv, lu, lv, rcon, zcon,
527  & armn, brmn, extra3, azmn, bzmn, extra4, blmn,
528  & clmn, extra1, extra2)
529  END IF
530 
531  l0pi = ns*(1 + nzeta*(ntheta2 - 1)) !u = pi, v = 0, js = ns
532  router = r1(ns,0) + r1(ns,1)
533  rinner = r1(l0pi,0) + r1(l0pi,1)
534  r00 = r1(1,0)
535  z00 = z1(1,0)
536 
537 !
538 ! COMPUTE CONSTRAINT RCON, ZCON
539 !
540 #ifndef _HBANGLE
541  rcon(:nrzt,0) = rcon(:nrzt,0) + rcon(:nrzt,1)*sqrts(:nrzt)
542  zcon(:nrzt,0) = zcon(:nrzt,0) + zcon(:nrzt,1)*sqrts(:nrzt)
543 #endif
544  ru0(:nrzt) = ru(:nrzt,0) + ru(:nrzt,1)*sqrts(:nrzt)
545  zu0(:nrzt) = zu(:nrzt,0) + zu(:nrzt,1)*sqrts(:nrzt)
546 
547 !
548 ! COMPUTE RCON0, ZCON0 FOR FIXED BOUNDARY BY SCALING EDGE VALUES
549 ! SCALE BY POWER OF SQRTS, RATHER THAN USE rcon0 = rcon, etc. THIS
550 ! PREVENTS A DISCONTINUITY WHEN RESTARTING FIXED BOUNDARY WITH NEW RCON0....
551 !
552 ! NOTE: IN ORDER TO MAKE INITIAL CONSTRAINT FORCES SAME FOR FREE/FIXED
553 ! BOUNDARY, WE SET RCON0,ZCON0 THE SAME INITIALLY, BUT TURN THEM OFF
554 ! SLOWLY IN FREE-BOUNDARY VACUUM LOOP (BELOW)
555 !
556 #ifndef _HBANGLE
557  IF (ictrl_prec2d .EQ. 2) THEN
558  rcon0(1:nrzt) = rcon(1:nrzt,0)
559  zcon0(1:nrzt) = zcon(1:nrzt,0)
560  ALLOCATE(bsqvac0(nznt))
561  ELSE IF (iter2 .EQ. iter1 .AND.
562  & ivac .le. 0 .AND.
563  & ictrl_prec2d .EQ. 0) THEN
564  DO l = 1, ns
565  rcon0(l:nrzt:ns) = rcon(ns:nrzt:ns,0)*sqrts(l:nrzt:ns)**2
566  zcon0(l:nrzt:ns) = zcon(ns:nrzt:ns,0)*sqrts(l:nrzt:ns)**2
567  END DO
568  END IF
569 #endif
570 !
571 ! COMPUTE S AND THETA DERIVATIVE OF R AND Z AND JACOBIAN ON HALF-GRID
572 !
573  CALL jacobian
574  IF (irst .EQ. 2 .AND.
575  & iequi .EQ. 0) THEN
576  gc=0
577  GOTO 100
578  END IF
579 
580 !
581 ! COMPUTE COVARIANT COMPONENTS OF B, MAGNETIC AND KINETIC
582 ! PRESSURE, AND METRIC ELEMENTS ON HALF-GRID
583 !
584  CALL second0 (tbcovon)
585  CALL bcovar (lu, lv)
586 
587  CALL second0 (tbcovoff)
588  timer(tbcov) = timer(tbcov) + (tbcovoff - tbcovon)
589 
590 ! COMPUTE VACUUM MAGNETIC PRESSURE AT PLASMA EDGE
591 ! NOTE: FOR FREE BOUNDARY RUNS, THE VALUE OF RBTOR=R*BTOR
592 ! AT THE PLASMA EDGE SHOULD BE ADJUSTED TO APPROXIMATELY
593 ! EQUAL THE VACUUM VALUE. THIS CAN BE DONE BY CHANGING
594 ! EITHER PHIEDGE OR THE INITIAL CROSS SECTION ACCORDING
595 ! TO THE SCALING LAW R*BTOR .EQ. PHIEDGE/(R1 * Z1).
596 
597  IF (lfreeb .and. iter2.gt.1 .and. iequi.eq.0) THEN
598  IF (ictrl_prec2d .le. 1 .and.
599  & (fsqr + fsqz) .le. 1.e-3_dp) THEN
600  ivac = ivac+1 !decreased from e-1 to e-3 - sph12/04
601  END IF
602  IF (nvskip0 .eq. 0) THEN
603  nvskip0 = max(1, nvacskip)
604  END IF
605  ivac0: IF (ivac .ge. 0) THEN
606 !SPH OFF: 6.20.17
607 ! IF INITIALLY ON, TURN OFF rcon0, zcon0 SLOWLY
608  IF (ictrl_prec2d .eq. 2) THEN
609  rcon0 = 0; zcon0 = 0
610  ELSE IF (ictrl_prec2d .eq. 0) THEN
611  rcon0 = 0.9_dp*rcon0; zcon0 = 0.9_dp*zcon0
612  END IF
613  CALL second0 (tvacon)
614  ivacskip = mod(iter2-iter1,nvacskip)
615  IF (ivac .LE. 2) THEN
616  ivacskip = 0
617  END IF
618 
619 ! EXTEND NVACSKIP AS EQUILIBRIUM CONVERGES
620  IF (ivacskip .eq. 0) THEN
621  nvacskip = one/max(1.e-1_dp, 1.e11_dp*(fsqr+fsqz))
622  nvacskip = max(nvacskip, nvskip0)
623  END IF
624 
625 !
626 ! NORMALLY, WHEN COMPUTING THE HESSIAN, IT IS SUFFICIENT TO
627 ! COMPUTE THE VARIATIONS IN THE "EXACT" SOLUTION, NOT THE ENTIRE
628 ! FIELD PERIOD SUM. THUS, FOR ictrl_prec2d >= 2, SET ivacskip = 1
629 ! FOR ictrl_prec2d = 1 (RUN WITH PRECONDITIONER APPLIED), MUST
630 ! COMPUTE EXACT VACUUM RESPONSE NOW.
631 !
632 ! THE EXCEPTION TO THIS IS IF WE ARE TESTING THE HESSIAN (lHess_exact=T),
633 ! THEN MUST USE FULL VACUUM CALCULATION TO COMPUTE IT (ivacskip=0)
634 !
635 ! lHess_exact = .FALSE.
636  IF (ictrl_prec2d .NE. 0) THEN
637  IF (lhess_exact .OR. ictrl_prec2d.EQ.2) THEN !Accurate Hessian
638  ivacskip = 0
639  ELSE
640  ivacskip = 1 !Fast vacuum calculation used to compute Hessian
641  END IF
642  ENDIF
643 
644 ! NOTE: gc contains correct edge values of r,z,l arrays
645 ! convert_sym, convert_asym have been applied to m=1 modes
646  CALL convert (rmnc, zmns, lmns, rmns, zmnc, lmnc, gc, ns)
647  IF (ictrl_prec2d.NE.3 .OR. l_edge) THEN
648 #ifdef _VACUUM2
649  CALL vac2_vacuum(rmnc, rmns, zmns, zmnc, xm, xn, ctor,
650  & ivacskip, mnmax)
651 
652 #else
653 ! DO NOT UPDATE THIS WHEN USING PRECONDITIONER: BREAKS TRI-DIAGONAL STRUCTURE
654  IF (ictrl_prec2d .EQ. 0 .OR.
655  & ictrl_prec2d .EQ. 2) THEN
656  raxis_nestor(1:nzeta) = r1(1:ns*nzeta:ns,0)
657  zaxis_nestor(1:nzeta) = z1(1:ns*nzeta:ns,0)
658  END IF
659 
660  CALL vacuum (rmnc, rmns, zmns, zmnc, xm, xn,
661  & ctor, rbtor, wint, ns, ivacskip, ivac,
662  & mnmax, ier_flag, lscreen)
663  IF (ictrl_prec2d .EQ. 2) THEN
664  bsqvac0 = bsqvac
665  END IF
666 #endif
667  ELSE
668  bsqvac = bsqvac0
669  ier_flag = 0
670  END IF
671 
672  IF (ier_flag .NE. 0) THEN
673  GOTO 100
674  END IF
675 !
676 ! RESET FIRST TIME FOR SOFT START
677 !
678  IF (ivac .eq. 1) THEN
679  irst = 2; delt0 = delt
680 #ifdef _HBANGLE
681  gc = xc; xc = xtempa
682 #endif
683  CALL restart_iter(delt0)
684 #ifdef _HBANGLE
685  xc = gc
686 #endif
687  irst = 1
688  END IF
689 
690 !
691 ! IN CASE PRESSURE IS NOT ZERO AT EXTRAPOLATED EDGE...
692 ! UNCOMMENT ALL "RPRES" COMMENTS HERE AND IN BCOVAR, FORCES ROUTINES
693 ! IF NON-VARIATIONAL FORCES ARE DESIRED
694 !
695 ! presf_ns = 1.5_dp*pres(ns) - 0.5_dp*pres(ns1)
696 ! MUST NOT BREAK TRI-DIAGONAL RADIAL COUPLING: OFFENDS PRECONDITIONER!
697  presf_ns = pmass(hs*(ns-1.5_dp))
698  IF (presf_ns .ne. zero) THEN
699  presf_ns = (pmass(one)/presf_ns) * pres(ns)
700  END IF
701 
702  lk = 0
703 ! gcon(:nrzt) = r1(:nrzt,0)+sqrts(:nrzt)*r1(:nrzt,1)
704 ! gcon(1+nrzt) = 0
705  DO l = ns, nrzt, ns
706  lk = lk + 1
707  bsqsav(lk,3) = 1.5_dp*bzmn_o(l) - 0.5_dp*bzmn_o(l-1)
708 #ifdef _ANIMEC
709  gcon(l) = bsqvac(lk) + pperp_ns(lk)
710 #else
711  gcon(l) = bsqvac(lk) + presf_ns
712 #endif
713 
714  rbsq(lk) = gcon(l)*(r1(l,0) + r1(l,1))*ohs
715  dbsq(lk) = abs(gcon(l)-bsqsav(lk,3))
716  END DO
717 !
718 ! COMPUTE m=0,n=0 EDGE "pedestals"
719 !
720 !! alphaR = hs*hs*ard(ns,1)
721 !! IF (alphaR .ne. zero) alphaR =
722 !! & hs*SUM(wint(ns:nrzt:ns)*zu0(ns:nrzt:ns)*rbsq)/alphaR
723 
724 !! PRINT *,' alphaR/r1(ns) = ', alphaR/gcon(ns)
725 
726  IF (ivac .eq. 1) THEN
727  bsqsav(:nznt,1) = bzmn_o(ns:nrzt:ns)
728  bsqsav(:nznt,2) = bsqvac(:nznt)
729  END IF
730  CALL second0 (tvacoff)
731  timer(tvac) = timer(tvac) + (tvacoff - tvacon)
732  IF (ictrl_prec2d .GE. 2) THEN
733  timer(tvac_2d) = timer(tvac_2d)+ (tvacoff - tvacon)
734  END IF
735 
736  END IF ivac0
737  END IF
738 
739 !
740 ! COMPUTE CONSTRAINT FORCE
741 !
742 #ifdef _HBANGLE
743  gcon = 0
744  IF (iequi .EQ. 1) THEN
745  GOTO 100
746  END IF
747 #else
748  IF (iequi .NE. 1) THEN
749  extra1(:nrzt,0) = (rcon(:nrzt,0) - rcon0(:nrzt))*ru0(:nrzt)
750  & + (zcon(:nrzt,0) - zcon0(:nrzt))*zu0(:nrzt)
751  CALL alias (gcon, extra1(:,0), gc, gc(1+mns), gc(1+2*mns),
752  & extra1(:,1))
753  ELSE
754  IF (lrecon) xc(:ns) = xc(:ns) + delr_mse
755  GOTO 100
756  END IF
757 #endif
758 !
759 ! COMPUTE MHD FORCES ON INTEGER-MESH
760 !
761  CALL forces
762 
763 !
764 ! FFT TEST
765 !
766 #ifdef _TEST_FOURIER
767 ! xc = xc*scalxc
768 ! CALL tomnsps_t (gc, xc, r1, ru, rv, z1, zu, zv)
769 ! IF (lasym) CALL tomnspa_t (gc, xc, r1, ru, rv, z1, zu, zv)
770 ! STOP
771 #endif
772 !
773 ! SYMMETRIZE FORCES (in u-v space): NOTE - gc IS SMALL BY FACTOR 2
774 ! IF lasym=T
775 !
776  IF (lasym) THEN
777  CALL symforce (armn, brmn, crmn, azmn, bzmn,
778  & czmn, blmn, clmn, rcon, zcon, r1, ru, rv,
779  & z1, zu, zv, extra3, extra4, extra1, extra2)
780 !NOT NECESSARY (EVEN THOUGH CORRECT) gc = 2*gc
781  END IF
782 
783 !
784 ! FOURIER-TRANSFORM MHD FORCES TO (M,N)-SPACE
785 !
786  CALL tomnsps (gc, armn, brmn, crmn, azmn, bzmn, czmn,
787  & blmn, clmn, rcon, zcon)
788 
789  IF (lasym) THEN
790  CALL tomnspa (gc, r1, ru, rv, z1, zu, zv,
791  & extra3, extra4, extra1, extra2)
792  END IF
793 
794 !================================================================
795 !
796 ! COMPUTE FORCE RESIDUALS (RAW AND PRECONDITIONED)
797 !
798 !================================================================
799  CALL second0 (treson)
800 
801  gc = gc * scalxc !!IS THIS CORRECT: SPH010214?
802 #ifdef _HBANGLE
803  CALL getfrho(gc)
804 #endif
805  CALL residue (gc, gc(1+irzloff), gc(1+2*irzloff))
806 ! Force new initial axis guess IF ALLOWED (l_moveaxis=T)
807  IF (lmove_axis .and.
808  & iter2 .eq. 1 .and.
809  & (fsqr + fsqz + fsql) .gt. 1.e2_dp) THEN
810  irst = 4
811  END IF
812 
813  CALL second0 (tresoff)
814  timer(tres) = timer(tres) + (tresoff - treson)
815 
816  100 CONTINUE
817 
818 #ifdef _HBANGLE
819  xc = xtempa
820 #endif
821  CALL second0 (tfunoff)
822  timer(tfun) = timer(tfun) + (tfunoff - tfunon)
823  IF (ictrl_prec2d .GE. 2) THEN
824  timer(tfun_2d) = timer(tfun_2d) + (tfunoff - tfunon)
825  END IF
826 
827  END SUBROUTINE funct3d
blocktridiagonalsolver
Solver for block tri-diagonal matrices. [Kalyan S. Perumalla, ORNL, 2009-2011].
Definition: blocktridiagonalsolver.f90:8
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17