V3FIT
pchelms.f90
Go to the documentation of this file.
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
11 !*******************************************************************************
12  MODULE pchelms
13  USE stel_kinds
14  USE island_params, mpol=>mpol_i, ntor=>ntor_i, ns=>ns_i
15  USE island_params, ntheta=>nu_i, nzeta=>nv_i, hs=>hs_i
16  USE island_params, ohs => ohs_i, mnmax=>mnmax_i, nfp=>nfp_i
17  USE hessian, ONLY: inithess, compute_hessian_blocks, asym_index, levmarq_param
18  USE blocktridiagonalsolver_s, ONLY: initialize, getranks
19  USE nscalingtools, ONLY: myenvvariables, parsolver, disp, startglobrow, endglobrow, rcounts, mpi_err
20  USE shared_data, ONLY: gc, xc, neqs, gc0, xc0, fsq_total, mblk_size, ndims
21  USE shared_data, ONLY: asubsmnsf, asubumncf, asubvmncf, &
22  fsupsmnsf, fsupumncf, fsupvmncf, &
23  asubsmncf, asubumnsf, asubvmnsf, &
24  fsupsmncf, fsupumnsf, fsupvmnsf
25  USE shared_data, ONLY: torflux, polflux, nsp, r1_i, z1_i, ru_i, &
26  zu_i, rv_i, zv_i, lcolscale
27  USE descriptor_mod, ONLY: iam, nprocs, lscalapack, inhessian, siesta_comm
28  USE v3_utilities, ONLY: assert
29 
30  USE vmec_info, ONLY: jcurrumnc, jcurrvmnc, jcurrumns, jcurrvmns
31  USE fourier
32  USE metrics, ONLY: tolowerh
33  USE quantities, ONLY: jacobh, jbsupsmnsh, jbsupumnch, jbsupvmnch, &
34  jbsupsmnch, jbsupumnsh, jbsupvmnsh
35 
36  USE mpi_inc
37 
38  IMPLICIT NONE
39 
40 !*******************************************************************************
41 ! pchelms module parameters
42 !*******************************************************************************
46  LOGICAL, PARAMETER, PRIVATE :: l_asedge = .true.
47 
48 ! FIXME: Remove these. They should either be passed in as arguments.
49  INTEGER, PRIVATE :: nsmin
50  INTEGER, PRIVATE :: nsmax
51  INTEGER, PRIVATE :: ns_match
52 
53  !Added SKS
54  LOGICAL, PRIVATE :: linit
55  LOGICAL, PRIVATE :: lhessian
56  REAL (dp), PRIVATE :: bnorm = -1
57  REAL (dp), PRIVATE :: line_bu
58 
59  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: bsupsijh
60  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: bsupuijh
61  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: bsupvijh
62  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: bsubsijh
63  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: bsubuijh
64  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: bsubvijh
65  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: jacobmnch
66  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: jacobmnsh
67 
68  CONTAINS
69 
70 !*******************************************************************************
71 ! UTILITY SUBROUTINES
72 !*******************************************************************************
73 !-------------------------------------------------------------------------------
75 !-------------------------------------------------------------------------------
76  SUBROUTINE run_pchelms
77  USE shared_data, ONLY: lverbose, siesta_curtor
78  USE diagnostics_mod, ONLY: divb, divb_rms
79  USE vmec_info, ONLY: vmec_curtor
80 
81  IMPLICIT NONE
82 
83 ! local variables
84  INTEGER :: istat
85  INTEGER :: js
86  REAL (dp) :: ton
87  REAL (dp) :: toff
88  LOGICAL :: flag
89  LOGICAL :: lcolscale_save
90 
91 ! local parameters
92  REAL (dp), PARAMETER :: levmarq_param_inital = 1.e-6_dp
93 
94 ! Start of executable code
95  siesta_curtor = 0.0
96 
97 ! Turn of column scaling. Save the value of the column scaling flag so it can
98 ! be reset afterwards.
99  lcolscale_save = lcolscale
100  lcolscale = .false.
101 
102  IF (lverbose .and. iam .EQ. 0) THEN
103  WRITE (*,1000)
104  END IF
105 
106  CALL init_pchelms
107 
108  CALL second0(ton)
109 
110 ! Add initial diffusive term to preconditioner instead
111  levmarq_param = levmarq_param_inital
112 
113 ! Note initHess must have already been called prior to this point.
114  CALL compute_hessian_blocks(compute_forces_lin, .false.)
115  CALL second0(toff)
116  IF (lverbose .and. iam .EQ. 0) THEN
117  WRITE (*,1001) asym_index, toff - ton
118  END IF
119 
120  CALL second0(ton)
121 ! Call conjugate gradient pchelms
122  CALL gmres_pchelms
123  CALL second0(toff)
124 
125  IF (lverbose .and. iam .EQ. 0) THEN
126  WRITE (*,1002) toff - ton
127  END IF
128 
129 ! Update solution (initial guess + perturbation) and check final forces
130  CALL backslv
131 
132 ! Dump output: MUST call backslv to get updated B's before dumping
133  DO js = 1, ns
134  CALL dump_a(js, 333)
135  CALL dump_b(js, 666)
136  END DO
137 
138  nsmin = startglobrow
139  nsmax = min(ns, endglobrow + 1)
140  CALL divb(nsmin, nsmax)
141 
142 #if defined(MPI_OPT)
143  CALL mpi_allreduce(mpi_in_place, siesta_curtor, 1, mpi_real8, mpi_sum, &
144  siesta_comm, mpi_err)
145 #endif
146 
147  IF (lverbose .and. iam .EQ. 0) THEN
148  WRITE (*,1003) divb_rms, siesta_curtor, vmec_curtor
149  END IF
150 
151  CALL check_current
152 
153  CALL cleanup
154 
155  lcolscale = lcolscale_save
156 
157 1000 FORMAT(/,'-----------------------------------------------', &
158  /,'STARTING EXTERNAL B-FIELD CALCULATION (PCHELMS)', &
159  /,'-----------------------------------------------')
160 1001 FORMAT(' ASYMMETRY INDEX: ',1p,e12.4,/, &
161  ' PCHELMS HESSIAN COMPUTATION TIME: ',f10.1,'s')
162 1002 FORMAT(' PCHELMS HESSIAN SOLVER TIME: ',f10.1,'s')
163 1003 FORMAT(/,' |DivB| (normed) = ',1p,e12.3,/, &
164  /,' SIESTA Curtor : ',e12.4,' VMEC Curtor : 'e12.4, &
165  /,'-----------------------------------------------', &
166  /,'ENDING PCHELMS', &
167  /,'-----------------------------------------------')
168 
169  END SUBROUTINE
170 
171 !-------------------------------------------------------------------------------
185 !-------------------------------------------------------------------------------
186  SUBROUTINE curla_pchelms(Asubsmnf, Asubumnf, Asubvmnf, &
187  jbsupsmnh, jbsupumnh, jbsupvmnh, iparity)
188  USE utilities, ONLY: curl_ftoh, set_bndy_fouier_m0, set_bndy_full_origin
189  USE island_params, ONLY: fourier_context
190 
191  IMPLICIT NONE
192 
193 ! Declare Arguments
194  REAL (dp), DIMENSION(:,:,:), POINTER :: Asubsmnf
195  REAL (dp), DIMENSION(:,:,:), POINTER :: Asubumnf
196  REAL (dp), DIMENSION(:,:,:), POINTER :: Asubvmnf
197  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: jbsupsmnh
198  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: jbsupumnh
199  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: jbsupvmnh
200  INTEGER, INTENT(in) :: iparity
201 
202 ! local variables
203  INTEGER :: fcomb
204  INTEGER :: fours
205  INTEGER :: fouruv
206  INTEGER :: moff
207  INTEGER :: noff
208  INTEGER :: nmin
209  INTEGER :: nmatch
210 
211  INTEGER :: m,n
212 
213 ! Start of executable code
214  moff = lbound(asubsmnf, 1)
215  noff = ntor + lbound(asubsmnf, 2)
216  nmin = nsmin
217 
218  IF (iparity .EQ. f_sin) THEN
219  fcomb = f_none
220  fours = f_sin
221  fouruv = f_cos
222  ELSE
223  fcomb = f_sum
224  fours = f_cos
225  fouruv = f_sin
226  END IF
227 
228  CALL set_bndy_fouier_m0(asubsmnf, asubumnf, asubvmnf, iparity)
229  IF (nsmin .EQ. 1) THEN
230  CALL set_bndy_full_origin(asubsmnf, asubumnf, asubvmnf)
231  END IF
232 
233 ! Dirichelet bc (fixed) for A_u, A_v (A_s floats) consistent with Poynting
234 ! flux = 0.
235 ! IF (nsmax .eq. ns .and. &
236 ! lHessian .and. &
237 ! .not.lInit) THEN
238 ! IF (.not.l_AsEdge) THEN
239 ! asubsmnf(:,:,ns) = 0
240 ! END IF
241 ! asubumnf(:,:,ns) = 0
242 ! asubvmnf(:,:,ns) = 0
243 ! END IF
244 
245 ! Need this to keep Hessian symmetric
246  IF (nmin .le. ns_match .and. lhessian .and. .not.linit) THEN
247  nmatch = min(ns_match, nsmax)
248  asubumnf(m1 + moff,:,nmin:nmatch) = 0
249  asubvmnf(m1 + moff,:,nmin:nmatch) = 0
250  END IF
251 
252  CALL curl_ftoh(asubsmnf, asubumnf, asubvmnf, &
253  jbsupsmnh, jbsupumnh, jbsupvmnh, iparity, nmin, nsmax)
254 
255  nmin = startglobrow
256 
257 ! Convert to real space.
258  CALL fourier_context%toijsp(jbsupsmnh(:,:,nmin:nsmax), bsupsijh, fcomb, &
259  fours)
260  CALL fourier_context%toijsp(jbsupumnh(:,:,nmin:nsmax), bsupuijh, fcomb, &
261  fouruv)
262  CALL fourier_context%toijsp(jbsupvmnh(:,:,nmin:nsmax), bsupvijh, fcomb, &
263  fouruv)
264 
265  END SUBROUTINE
266 
267 !-------------------------------------------------------------------------------
283 !-------------------------------------------------------------------------------
284  SUBROUTINE curlb_pchelms(ksupsmnf, ksupumnf, ksupvmnf, asubsmnf, &
285  bsupsmnh, iparity, curtor)
286  USE utilities, ONLY: curl_htof, gradienthalf, set_bndy_fouier_m0
287  USE v3_utilities, ONLY: assert_eq
288  USE island_params, ONLY: fourier_context
289 
290  IMPLICIT NONE
291 
292 ! Declare Arguments
293  REAL(dp), POINTER, DIMENSION(:,:,:) :: ksupsmnf
294  REAL(dp), POINTER, DIMENSION(:,:,:) :: ksupumnf
295  REAL(dp), POINTER, DIMENSION(:,:,:) :: ksupvmnf
296  REAL(dp), POINTER, DIMENSION(:,:,:) :: asubsmnf
297  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(IN) :: bsupsmnh
298  INTEGER, INTENT(IN) :: iparity
299  REAL (dp), INTENT(inout) :: curtor
300 
301 ! Local variables
302  INTEGER :: fours
303  INTEGER :: fouruv
304  INTEGER :: moff
305  INTEGER :: noff
306  INTEGER :: nl
307  INTEGER :: nh
308  INTEGER :: istat
309  INTEGER :: nmin
310  INTEGER :: nmax
311  INTEGER :: mp
312  INTEGER :: np
313  INTEGER :: m
314  INTEGER :: n
315  INTEGER :: sparity
316  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: dA_sds
317  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: bsubsmnh
318  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: bsubumnh
319  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: bsubvmnh
320  REAL (dp) :: Diff
321  REAL (dp), DIMENSION(nsp) :: t1(nsp)
322 
323 ! Start of executable code
324  IF (iparity .eq. f_sin) THEN
325  fours = f_sin
326  fouruv = f_cos
327  sparity = 1
328  ELSE
329  fours = f_cos
330  fouruv = f_sin
331  sparity = -1
332  END IF
333 
334  nmin = startglobrow
335  nmax = min(ns, endglobrow+1)
336  moff = lbound(ksupsmnf,1)
337  noff = lbound(ksupsmnf,2) + ntor
338 
339  ALLOCATE(bsubsmnh(0:mpol,-ntor:ntor,nmin:nmax), &
340  bsubumnh(0:mpol,-ntor:ntor,nmin:nmax), &
341  bsubvmnh(0:mpol,-ntor:ntor,nmin:nmax), stat=istat)
342  CALL assert(istat .eq. 0,'Allocation1 failed in CURLB_PCHELMS')
343 
344  CALL fourier_context%tomnsp(bsubsijh, bsubsmnh(:,:,nmin:nmax), fours) ! Compute BSUBXMNH
345  CALL fourier_context%tomnsp(bsubuijh, bsubumnh(:,:,nmin:nmax), fouruv)
346  CALL fourier_context%tomnsp(bsubvijh, bsubvmnh(:,:,nmin:nmax), fouruv)
347 
348  CALL set_bndy_fouier_m0(bsubsmnh, bsubumnh, bsubvmnh, fours)
349 
350 ! Sub s term has sin parity so the subu term has cos parity.
351  IF (nmax .eq. ns .and. iparity .eq. f_sin) THEN
352  line_bu = bsubumnh(0,0,ns)
353  END IF
354 
355  CALL curl_htof(bsubsmnh, bsubumnh, bsubvmnh, &
356  ksupsmnf, ksupumnf, ksupvmnf, &
357  iparity, nmin, nmax, .false., curtor)
358 
359  nmax = endglobrow
360  DEALLOCATE (bsubsmnh, bsubumnh, bsubvmnh, stat=istat)
361  CALL assert_eq(0, istat, 'Deallocation failed in CURLB_PCHELMS')
362 #if 0
363 ! Add smoothing diffusion for A_s in Hessian:
364  IF (lhessian) THEN
365  np = max(1, nmin - 1)
366  mp = min(ns, nmax + 1)
367  ALLOCATE (da_sds(0:mpol,-ntor:ntor,np:mp))
368  CALL gradienthalf(da_sds, asubsmnf)
369  diff = 1.e-4_dp*ohs
370  ksupsmnf(:,:,nmin:mp-1) = ksupsmnf(:,:,nmin:mp-1) &
371  + diff*(da_sds(:,:,nmin+1:mp) - &
372  da_sds(:,:,nmin:mp-1))
373  IF (nmax .eq. ns) THEN
374  ksupsmnf(:,:,ns) = ksupsmnf(:,:,ns) - 2*diff*da_sds(:,:,ns)
375  END IF
376  IF (iparity .eq. isym) THEN
377  ksupsmnf(m0+moff,n0+noff,nmin:nmax) = 0
378  END IF
379  DEALLOCATE(da_sds)
380  END IF
381 
382  nh = min(nmax, nsp)
383  IF (nh .GE. nmin) THEN
384 !Force b^s -> 0 for js <= nsp (add (B^s)**2 to energy)
385  nl = nmin
386  diff = 0.1_dp
387  DO m = 0,mpol
388  mp = m*sparity
389  DO n = -ntor,ntor
390  np = n*sparity*nfp
391  t1(nl:nh) = diff*(bsupsmnh(m,n,nl:nh) &
392  + bsupsmnh(m,n,nl+1:nh+1))
393  ksupumnf(m+moff,n+noff,nl:nh) = ksupumnf(m+moff,n+noff,nl:nh) &
394  - np*t1(nl:nh)
395  ksupvmnf(m+moff,n+noff,nl:nh) = ksupvmnf(m+moff,n+noff,nl:nh) &
396  + mp*t1(nl:nh)
397  END DO
398  END DO
399  END IF
400 #endif
401 
402 ! Origin boundary conditions: If asubs(0,:,1) = 0 in curla, must set
403 ! ksupsmnf(0,:,1) = 0 here for hessian symmetry.
404 ! Need factor of 1/2 for hessian symmetry at js=1
405 ! IF (nmin .eq. 1) THEN
406 ! ksupsmnf(:,:,1) = 0
407 ! ksupumnf(:,:,1) = 0
408 ! ksupvmnf(m0+moff,:,1) = 0.5*ksupvmnf(m0+moff,:,1)
409 ! ksupvmnf(m1+moff,:,1) = 0
410 ! END IF
411 ! IF (nmax .eq. ns) THEN
412 ! ksupsmnf(:,:,ns) = 0.5*ksupsmnf(:,:,ns)
413 ! ksupumnf(:,:,ns) = 0.5*ksupumnf(:,:,ns)
414 ! ksupvmnf(:,:,ns) = 0.5*ksupvmnf(:,:,ns)
415 ! END IF
416 
417 ! Flip sign to get correct diagonal sign (>0).
418  ksupsmnf(:,:,nmin:nmax) = -ksupsmnf(:,:,nmin:nmax)
419  ksupumnf(:,:,nmin:nmax) = -ksupumnf(:,:,nmin:nmax)
420  ksupvmnf(:,:,nmin:nmax) = -ksupvmnf(:,:,nmin:nmax)
421 
422  END SUBROUTINE
423 
424 !-------------------------------------------------------------------------------
439 !-------------------------------------------------------------------------------
440  SUBROUTINE cyl2vmec_a(A_r, A_p, A_z, cA_s, cA_u, cA_v)
441 
442  IMPLICIT NONE
443 
444 ! Declare Arguments
445  REAL (dp), DIMENSION(:,:), INTENT(IN) :: A_r
446  REAL (dp), DIMENSION(:,:), INTENT(IN) :: A_p
447  REAL (dp), DIMENSION(:,:), INTENT(IN) :: A_z
448  REAL (dp), DIMENSION(:,:), INTENT(OUT) :: cA_s
449  REAL (dp), DIMENSION(:,:), INTENT(OUT) :: cA_u
450  REAL (dp), DIMENSION(:,:), INTENT(OUT) :: cA_v
451 
452 ! Local variables
453  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: rs
454  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: zs
455  INTEGER :: u
456  INTEGER :: v
457  INTEGER :: js
458 
459 ! Start of executable code
460  ALLOCATE(rs(ntheta, nzeta))
461  ALLOCATE(zs(ntheta, nzeta))
462  rs = (r1_i(:,:,ns) - r1_i(:,:,ns - 1))*ohs
463  zs = (z1_i(:,:,ns) - z1_i(:,:,ns - 1))*ohs
464 
465  ca_s(:,:) = a_r(:,:)*rs(:,:) + a_z(:,:)*zs(:,:)
466  ca_u(:,:) = a_r(:,:)*ru_i(:,:,ns) + a_z(:,:)*zu_i(:,:,ns)
467  ca_v(:,:) = a_r(:,:)*rv_i(:,:,ns) + a_z(:,:)*zv_i(:,:,ns) &
468  + a_p(:,:)*r1_i(:,:,ns)
469 
470  DEALLOCATE(rs,zs)
471 
472  END SUBROUTINE
473 
474 !-------------------------------------------------------------------------------
493 !-------------------------------------------------------------------------------
494  SUBROUTINE init_a(cA_s, cA_u, cA_v, A_s, A_u, A_v, parity)
495  USE metrics, ONLY: phipf_i, chipf_i
496  USE utilities, ONLY: set_bndy_fouier_m0, set_bndy_full_origin
497  USE island_params, ONLY: fourier_context
498 
499  IMPLICIT NONE
500 
501 ! Declare Arguments
502  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: cA_s
503  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: cA_u
504  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: cA_v
505  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(out) :: A_s
506  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(out) :: A_u
507  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(out) :: A_v
508  INTEGER, INTENT(IN) :: parity
509 
510 ! Local variables
511  INTEGER :: js
512  INTEGER :: m
513  INTEGER :: n
514  INTEGER :: nmatch
515  INTEGER :: fours
516  INTEGER :: fouruv
517  REAL (dp) :: t1
518  REAL (dp) :: p1
519  REAL (dp) :: q1
520 
521 ! Start of executable code
522 ! Only do this once, so no need to parallelize
523  a_s = 0
524  a_u = 0
525  a_v = 0
526 
527  IF (parity .EQ. f_sin) THEN
528 ! e_s (sin), e_u (cos), e_v (cos)
529  fours = f_sin
530  fouruv = f_cos
531  ELSE
532 ! e_s (cos), e_u (sin), e_v (sin)
533  fours = f_cos
534  fouruv = f_sin
535  END IF
536 
537  CALL fourier_context%tomnsp(ca_s, a_s(:,:,ns:ns), fours)
538  CALL fourier_context%tomnsp(ca_u, a_u(:,:,ns:ns), fouruv)
539  CALL fourier_context%tomnsp(ca_v, a_v(:,:,ns:ns), fouruv)
540 
541 ! Extapolate the vector potential from the edge to the center.
542  nmatch = 0
543  t1 = 1
544  DO js = 1, ns
545  p1 = real(js-1,dp)/(ns-1)
546  p1 = p1*p1
547  q1 = 1 !p1 ~ flux
548  DO m = 0, mpol
549  q1 = q1*p1 !p1**m
550  DO n = -ntor, ntor
551  a_s(m,n,js) = a_s(m,n,ns)*q1
552  a_u(m,n,js) = a_u(m,n,ns)*q1
553  a_v(m,n,js) = a_v(m,n,ns)*q1
554  IF (nmatch .gt. 0) THEN
555  a_s(m,n,js) = t1*a_s(m,n,js) + (1 - t1)*a_s(m,n,nmatch)
556  a_u(m,n,js) = t1*a_u(m,n,js) + (1 - t1)*a_u(m,n,nmatch)
557  a_v(m,n,js) = t1*a_v(m,n,js) + (1 - t1)*a_v(m,n,nmatch)
558  END IF
559  END DO
560  END DO
561  END DO
562 
563  CALL set_bndy_fouier_m0(a_s, a_u, a_v, fours)
564  CALL set_bndy_full_origin(a_s, a_u, a_v)
565 
566  END SUBROUTINE
567 
568 !-------------------------------------------------------------------------------
584 !-------------------------------------------------------------------------------
585  SUBROUTINE init_f(Fsupsmn, Fsupumn, Fsupvmn, jcurrumn, jcurrvmn)
586 
587  IMPLICIT NONE
588 
589 ! Declare Arguments
590  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(out) :: Fsupsmn
591  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(out) :: Fsupumn
592  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(out) :: Fsupvmn
593  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(in) :: jcurrumn
594  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(in) :: jcurrvmn
595 
596 ! Local variables
597  INTEGER :: js
598  REAL (dp) :: t1
599 
600 ! Start of executable code
601 !NOTE: mu0 factored into j's in metric; (note: F -> -F in CURLB_PCHELMS,
602 !to get correct sign of hessian diagonal elements)
603  nsmin = startglobrow
604  nsmax = endglobrow
605 
606  CALL assert(all(jcurrumn(:,:,nsp+1:).EQ.0._dp), 'jcurumn != 0 in vacuum!')
607  CALL assert(all(jcurrvmn(:,:,nsp+1:).EQ.0._dp), 'jcurvmn != 0 in vacuum!')
608 
609  fsupsmn(:,:,nsmin:nsmax) = 0
610  fsupumn(:,:,nsmin:nsmax) = jcurrumn(:,:,nsmin:nsmax)
611  fsupvmn(:,:,nsmin:nsmax) = jcurrvmn(:,:,nsmin:nsmax)
612 
613 ! Smooth edge current gradient.
614  DO js = nsp - 10, nsp
615  t1 = 1 - real(js-1,dp)/(nsp-1)
616  fsupumn(:,:,js) = t1*fsupumn(:,:,js)
617  fsupvmn(:,:,js) = t1*fsupvmn(:,:,js)
618  END DO
619 
620 ! Fsupumn(:,:,1) = 0
621 ! IF (iparity .eq. 0) THEN
622 ! Fsupvmn(m1:,:,1) = 0
623 ! ELSE
624 ! Fsupvmn(m1:,:,1) = 0
625 ! END IF
626 
627 !IMPROVES CONVERGENCE BY SUPPRESSING jac*B^u(m=1) NEAR AXIS (AT LEAST FOR W7X CASE)
628 ! IF (ns_match .GT. 0) RETURN
629 ! Fsupumn(m1,:,1:nsp/4) = 0
630 ! Fsupvmn(m1,:,1:nsp/4) = 0
631 
632  END SUBROUTINE
633 
634 !-------------------------------------------------------------------------------
648 !-------------------------------------------------------------------------------
649  SUBROUTINE compare_current(Fsupumn, Fsupvmn, jcurrumn, jcurrvmn)
650  USE island_params, ONLY: fourier_context
651  USE v3_utilities, ONLY: assert_eq
652 
653  IMPLICIT NONE
654 
655 ! Declare Arguments
656  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: Fsupumn
657  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: Fsupvmn
658  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: jcurrumn
659  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: jcurrvmn
660 
661 ! Local variables
662  INTEGER :: ns_mid
663  INTEGER :: m
664  INTEGER :: n
665  INTEGER :: moff
666  INTEGER :: noff
667  REAL (dp) :: t1
668 
669 ! Start of executable code
670 ! Note: mu0 factored into j's in metric; minus sign because r = b - Ax
671 ! Only do this once, so no need to parallelize (note: F -> -F to get
672 ! correct sign of hessian diagonal elements)
673  nsmin = startglobrow
674  nsmax = endglobrow
675  ns_mid = (nsmin + nsmax)/2
676 
677  CALL assert_eq(lbound(fsupumn,1), lbound(jcurrumn,1), ' LBOUND1 FAILED')
678  CALL assert_eq(lbound(fsupumn,2), lbound(jcurrumn,2), ' LBOUND2 FAILED')
679 
680  WRITE (1000 + iam, 1000) ns_mid
681  DO n = -ntor, ntor
682  noff = n + ntor + lbound(fsupumn,2)
683  DO m = 0, mpol
684  moff = m + lbound(fsupumn,1)
685  t1 = fourier_context%orthonorm(m, n)
686  WRITE (1000 + iam, 1001) m, n, &
687  jcurrumn(moff,noff,ns_mid)*t1, &
688  -fsupumn(moff,noff,ns_mid)*t1, &
689  jcurrvmn(moff,noff,ns_mid)*t1, &
690  -fsupvmn(moff,noff,ns_mid)*t1
691  END DO
692  END DO
693 
694 1000 FORMAT(' JS = ',i4,/ &
695  ' M N J^u F^u J^v F^v')
696 1001 FORMAT(2i6, 1p,4e12.4)
697 
698  END SUBROUTINE
699 
700 !-------------------------------------------------------------------------------
705 !-------------------------------------------------------------------------------
706  SUBROUTINE dump_a(js, iunit)
707  USE shared_data, ONLY: lasym
708 
709  IMPLICIT NONE
710 
711 ! Declare Arguments
712  INTEGER, INTENT(in) :: js
713  INTEGER, INTENT(in) :: iunit
714 
715 ! Local variables
716  INTEGER :: m
717  INTEGER :: n
718  INTEGER :: moff
719  INTEGER :: noff
720  INTEGER :: mp
721  INTEGER :: np
722  INTEGER :: nb
723  INTEGER :: is
724 
725 ! Start of executable code
726 ! xc = xc+xc0 ON ENTRY
727  nsmin = lbound(asubsmnsf,3)
728  nsmax = ubound(asubsmnsf,3)
729  moff = lbound(asubsmnsf,1)
730  noff = lbound(asubsmnsf,2)
731 
732  IF (js .lt. nsmin .or. js .gt. nsmax) THEN
733  RETURN
734  END IF
735  nb = min(2, ntor)
736 
737  IF (js .EQ. 1) THEN
738  WRITE(iunit + 10, 1000)
739  DO n = -nb, nb
740  np = n + ntor + noff
741 
742  DO m = 0, 1
743  mp = m + moff
744  IF (n .lt. 0 .and. m .eq. 0) cycle
745  WRITE(iunit + 10, 1001) n, m
746  DO is = 1, ns
747  WRITE(iunit + 10,1002) is, &
748  asubsmnsf(mp,np,is) * &
749  fourier_context%orthonorm(m,n), &
750  asubumncf(mp,np,is) * &
751  fourier_context%orthonorm(m,n), &
752  asubvmncf(mp,np,is) * &
753  fourier_context%orthonorm(m,n)
754  END DO
755  END DO
756  WRITE(iunit + 10,*)
757  END DO
758  END IF
759 
760  WRITE(iunit, 1003) js, neqs
761  DO n = -ntor, ntor
762  np = n + ntor + noff
763  DO m = 0, mpol
764  mp = m + moff
765  WRITE(iunit, 1004) m, n, asubsmnsf(mp,np,js) * &
766  fourier_context%orthonorm(m,n), &
767  asubumncf(mp,np,js) * &
768  fourier_context%orthonorm(m,n), &
769  asubvmncf(mp,np,js) * &
770  fourier_context%orthonorm(m,n), &
771  jcurrumnc(m,n,js) * &
772  fourier_context%orthonorm(m,n), &
773  jcurrvmnc(m,n,js) * &
774  fourier_context%orthonorm(m,n)
775  IF (.not.lasym) cycle
776  WRITE(iunit + 100, 1004) m, n, asubsmncf(mp,np,js) * &
777  fourier_context%orthonorm(m,n), &
778  asubumnsf(mp,np,js) * &
779  fourier_context%orthonorm(m,n), &
780  asubvmnsf(mp,np,js) * &
781  fourier_context%orthonorm(m,n), &
782  jcurrumns(m,n,js) * &
783  fourier_context%orthonorm(m,n), &
784  jcurrvmns(m,n,js) * &
785  fourier_context%orthonorm(m,n)
786  END DO
787  END DO
788 
789  WRITE(iunit,*)
790 
791 1000 FORMAT(' RAD A_s A_u A_v')
792 1001 FORMAT(' N=',i2,' M=',i2)
793 1002 FORMAT(i4,3(1p,e12.4))
794 1003 FORMAT(' RADIUS(JS): ',i4,' NEQS: ',i6, &
795  /,' M N A_s A_u A_v ', &
796  'jcurru jcurrv')
797 1004 FORMAT(2i8, 5(1p,e12.4))
798 
799  END SUBROUTINE
800 
801 !-------------------------------------------------------------------------------
806 !-------------------------------------------------------------------------------
807  SUBROUTINE dump_b(js, iunit)
808  USE metrics, ONLY: r1_i, chipf_i, phipf_i
809  USE hessian, ONLY: gather_array
810  USE siesta_namelist, ONLY: nsin
811  USE vmec_info
812  USE shared_data, ONLY: lasym
813  USE island_params, ONLY: fourier_context
814 
815  IMPLICIT NONE
816 
817 ! Declare Arguments
818  INTEGER, INTENT(IN) :: js
819  INTEGER, INTENT(IN) :: iunit
820 
821 ! Local variables
822  INTEGER :: m
823  INTEGER :: n
824  INTEGER :: is
825  INTEGER :: i
826  INTEGER :: j
827  INTEGER :: n1
828  INTEGER :: n2
829  REAL (dp) :: rjac
830  REAL (dp) :: rbsupv
831  REAL (dp) :: rho
832  REAL (dp) :: t1
833 
834 ! Start of executable code
835 ! xc = xc+xc0 ON ENTRY
836 ! Here, bsupX = jac*bsupX
837  n1 = startglobrow
838  n2 = endglobrow
839 
840  IF (js .EQ. 1) THEN
841  ALLOCATE(jacobmnch(0:mpol,-ntor:ntor,ns))
842  CALL fourier_context%tomnsp(jacobh, jacobmnch, f_cos)
843  IF (lasym) THEN
844  ALLOCATE(jacobmnsh(0:mpol,-ntor:ntor,ns))
845  CALL fourier_context%tomnsp(jacobh, jacobmnsh, f_sin)
846  END IF
847 
848  CALL gather_array(jbsupsmnsh)
849  CALL gather_array(jbsupumnch)
850  CALL gather_array(jbsupvmnch)
851  IF (lasym) THEN
852  CALL gather_array(jbsupsmnch)
853  CALL gather_array(jbsupumnsh)
854  CALL gather_array(jbsupvmnsh)
855  END IF
856 
857  WRITE(iunit + 10,1000)
858 
859  IF (endglobrow .EQ. ns) THEN
860  is = ns
861  DO m = 0, mpol
862  DO n = -ntor, ntor
863  IF (n.lt.0 .AND. m.eq.0) cycle
864  t1 = -fourier_context%orthonorm(m,n)
865  WRITE(iunit + 10, 1001) m, n, &
866  (1.5_dp*jbsupsmnsh(m,n,is) - &
867  0.5_dp*jbsupsmnsh(m,n,is - 1))*t1, &
868  (1.5_dp*jbsupumnch(m,n,is) - &
869  0.5_dp*jbsupumnch(m,n,is - 1))*t1, &
870  (1.5_dp*jbsupvmnch(m,n,is) - &
871  0.5_dp*jbsupvmnch(m,n,is - 1))*t1
872  END DO
873  END DO
874 
875  WRITE(iunit + 10,1002)
876  DO n = -1,1
877  DO m = 0, 1
878  IF (n.LT.0 .AND. m.EQ.0) cycle
879  t1 = -fourier_context%orthonorm(m,n)
880  WRITE(iunit + 10,1003) n, m
881  DO is = 2, ns
882  WRITE(iunit + 10, 1004) is, jbsupsmnsh(m,n,is)*t1, &
883  jbsupumnch(m,n,is)*t1, &
884  jbsupvmnch(m,n,is)*t1
885  END DO
886  WRITE(iunit+10,*)
887  END DO
888  END DO
889 
890  CALL flush(iunit + 10)
891  END IF !endglobrow=ns
892 
893  RETURN
894  END IF !js=1
895 
896  IF (iam .EQ. 0) THEN
897  WRITE(iunit, 1005) js, neqs
898 
899  WRITE(iunit, 1006)
900  IF (lasym) THEN
901  WRITE(iunit + 100, 1007)
902  END IF
903 
904  DO n = -ntor, ntor
905  DO m = 0, mpol
906  t1 = -fourier_context%orthonorm(m,n)
907  WRITE(iunit, 1008) m, n, jbsupsmnsh(m,n,js)*t1, &
908  jbsupumnch(m,n,js)*t1, &
909  jbsupvmnch(m,n,js)*t1, &
910  jacobmnch(m,n,js)*t1
911  IF (lasym) THEN
912  WRITE(iunit+100, 1008) m, n, jbsupsmnch(m,n,js)*t1, &
913  jbsupumnsh(m,n,js)*t1, &
914  jbsupvmnsh(m,n,js)*t1, &
915  jacobmnsh(m,n,js)*t1
916  END IF
917  END DO
918  END DO
919 
920  WRITE(iunit,*)
921  CALL flush(iunit)
922  END IF !iam=0
923 
924  IF (js .NE. ns) THEN
925  RETURN
926  END IF
927 
928 ! Estimate iota profile on extended mesh. Note bsupx is jac*bupx here and minus
929 ! sign for sign of jacobian. FIXME: Sign of the jacobian is hard coded here.
930 ! This should be read from the wout file.
931 ! chipf_i(1) = 0.0
932 ! phipf_i(1) = 0.0
933 ! DO is = 2, ns - 1
934  DO is = nsin + 1, ns - 1, 1
935  chipf_i(is) = -(jbsupumnch(0,0,is + 1) + jbsupumnch(0,0,is))*0.5
936  phipf_i(is) = -(jbsupvmnch(0,0,is + 1) + jbsupvmnch(0,0,is))*0.5
937  END DO
938  IF (ns .gt. nsin) THEN
939  chipf_i(ns) = -2.5*jbsupumnch(0,0,ns) + 1.5*jbsupumnch(0,0,ns - 1)
940  phipf_i(ns) = -2.5*jbsupvmnch(0,0,ns) + 1.5*jbsupvmnch(0,0,ns - 1)
941  END IF
942 
943 ! Check Edge value
944  IF (iam .EQ. 0) THEN
945 
946  rjac = 0
947  rbsupv = 0
948  DO n = -ntor, ntor
949  DO m = 0, mpol
950  rjac = rjac + jacobmnch(m,n,ns)*fourier_context%orthonorm(m,n)
951  rbsupv = rbsupv &
952  + jbsupvmnch(m,n,ns)*fourier_context%orthonorm(m,n)
953  END DO
954  END DO
955 
956  rbsupv = rbsupv*r1_i(1,1,ns)/rjac
957 
958  WRITE (*,1009) rjac, jacobh(1,1,ns), rbsupv
959 
960 ! Estimate iota profile. Note bsupx is jac*bsupx here and minus sign for sign
961 ! of jacobian. FIXME: Sign of the jacobian is hard coded here. This should be
962 ! read from the wout file.
963  WRITE (*, 1010)
964  DO is = 2, ns!, 10 !Skip Every 10
965 ! DO is = 2, ns
966  rho = (is-1.5_dp)*hs
967  WRITE (*,1011) rho, 0.5_dp*(chipf_i(is) + chipf_i(is - 1)), &
968  -jbsupumnch(0,0,is), &
969  -(asubvmncf(1,ntor + 1,is - 1) - &
970  asubvmncf(1,ntor + 1,is))*ohs, &
971  0.5_dp*(phipf_i(is) + phipf_i(is - 1)), &
972  -jbsupvmnch(0,0,is), &
973  -(asubumncf(1,ntor + 1,is) - &
974  asubumncf(1,ntor + 1,is-1))*ohs
975  END DO
976  END IF !iam = 0
977 
978  chipf_i(1) = 0.0
979  phipf_i(1) = 0.0
980  DO is = 2, nsin - 1
981  chipf_i(is) = -(jbsupumnch(0,0,is + 1) + jbsupumnch(0,0,is))*0.5
982  phipf_i(is) = -(jbsupvmnch(0,0,is + 1) + jbsupvmnch(0,0,is))*0.5
983  END DO
984 
985 1000 FORMAT(' EDGE VALUES OF jB^X FOR ALL M,N',/, &
986  ' M N jB^s jB^u jB^v')
987 1001 FORMAT(2i4,3(1p,e10.2))
988 1002 FORMAT(/,' RAD jB^s jB^u jB^v')
989 1003 FORMAT(' N=',i2,' M=',i2)
990 1004 FORMAT(i4,3(1p,e11.3))
991 1005 FORMAT('HALF RADIUS (JS): ',i4, ' NEQS: ',i6)
992 1006 FORMAT(' M N jB^s(sin) jB^u(cos) jB^v(cos) jacob(cos)')
993 1007 FORMAT(' M N jB^s(cos) jB^u(sin) jB^v(sin) jacob(sin)')
994 1008 FORMAT(2i8, 5(1p,e12.4))
995 1009 FORMAT(' jac(u=0,v=0,ns): ',1p,e12.3,' jacobh: ',1p,e12.3,/ &
996  ' R*B^v(u=0,v=0,ns) : ',1p,e12.3)
997 1010 FORMAT(/,' RAD CHIP(WOUT) CHIP(PCH) CHIP(VP)' &
998  ' PHIP(WOUT) PHIP(PCH) PHIP(VP)')
999 1011 FORMAT(f7.3,1p,6e14.3)
1000 
1001  END SUBROUTINE
1002 
1003 !-------------------------------------------------------------------------------
1005 !-------------------------------------------------------------------------------
1006  SUBROUTINE getfsq_pchelms
1007  USE shared_data, ONLY: fsupsmnsf, fsupumncf, fsupvmncf, &
1008  fsupsmncf, fsupumnsf, fsupvmnsf, lasym
1009  USE mpi_inc
1010 
1011  IMPLICIT NONE
1012 
1013 ! Start of executable code
1014  nsmin = startglobrow
1015  nsmax = endglobrow
1016 
1017  fsq_total = sum(fsupsmnsf(:,:,nsmin:nsmax)**2 + &
1018  fsupumncf(:,:,nsmin:nsmax)**2 + &
1019  fsupvmncf(:,:,nsmin:nsmax)**2)
1020  IF (lasym) THEN
1021  fsq_total = fsq_total &
1022  + sum(fsupsmncf(:,:,nsmin:nsmax)**2 + &
1023  fsupumnsf(:,:,nsmin:nsmax)**2 + &
1024  fsupvmnsf(:,:,nsmin:nsmax)**2)
1025  END IF
1026 #if defined(MPI_OPT)
1027  CALL mpi_allreduce(mpi_in_place, fsq_total, 1, mpi_real8, mpi_sum, &
1028  siesta_comm, mpi_err)
1029 #endif
1030  IF (bnorm .GT. zero) fsq_total = fsq_total/bnorm
1031 
1032  END SUBROUTINE
1033 
1034 !-------------------------------------------------------------------------------
1036 !-------------------------------------------------------------------------------
1037  SUBROUTINE cleanup
1038  USE shared_data, ONLY: lasym
1039 
1040  IMPLICIT NONE
1041 
1042 ! Start of executable code
1043  DEALLOCATE(xc0, gc0)
1044 
1045  IF (ALLOCATED(jacobmnch)) THEN
1046  DEALLOCATE(jacobmnch)
1047  END IF
1048  IF (lasym .and. ALLOCATED(jacobmnsh)) THEN
1049  DEALLOCATE(jacobmnsh)
1050  END IF
1051 
1052  DEALLOCATE(bsupsijh, bsupuijh, bsupvijh)
1053  DEALLOCATE(bsubsijh, bsubuijh, bsubvijh)
1054 
1055  xc = 0
1056  gc = 0
1057 
1058  END SUBROUTINE
1059 
1060 !-------------------------------------------------------------------------------
1063 !-------------------------------------------------------------------------------
1064  SUBROUTINE check_current
1065  USE metrics, ONLY: jsupvdota
1066 
1067  IMPLICIT NONE
1068 
1069 ! Local variables
1070  INTEGER :: i
1071  REAL (dp) :: currv_int
1072 
1073 ! Start of executable code
1074 ! Average over toroidal angle too (should not dependent on toroidal angle)
1075  IF (endglobrow .eq. ns) THEN
1076  currv_int = hs*(sum(jcurrvmnc(0,0,1:ns - 1)) + jcurrvmnc(0,0,ns)/2)
1077  WRITE (*,1000) line_bu, currv_int, jsupvdota
1078  END IF
1079 
1080 1000 FORMAT(' Line Integral B_u: ',1p,e12.4,4x, &
1081  ' Surface-integrated current: ',1p,e12.4,4x, &
1082  ' VMEC Surface-integrated current: ',1p,e12.4)
1083 
1084  END SUBROUTINE
1085 
1086 !-------------------------------------------------------------------------------
1088 !-------------------------------------------------------------------------------
1089  SUBROUTINE compute_forces_lin
1090 
1091  IMPLICIT NONE
1092 
1093 ! Start of executable code
1094  lhessian = .true.
1095  CALL compute_forces(.false.)
1096  lhessian = .false.
1097 
1098  END SUBROUTINE
1099 
1100 !-------------------------------------------------------------------------------
1102 !-------------------------------------------------------------------------------
1103  SUBROUTINE backslv
1104  USE metrics, ONLY: r1_i
1105  USE hessian, ONLY: gather_array
1106  USE shared_data, ONLY: lasym
1107 
1108  IMPLICIT NONE
1109 
1110 ! Start of executable code
1111 ! update solution (already gathered to all processors)
1112  xc = xc + xc0
1113 
1114  CALL init_f(fsupsmnsf, fsupumncf, fsupvmncf, &
1115  jcurrumnc, jcurrvmnc)
1116  IF (lasym) THEN
1117  CALL init_f(fsupsmncf, fsupumnsf, fsupvmnsf, &
1118  jcurrumns, jcurrvmns)
1119  END IF
1120  gc0 = gc
1121 
1122  CALL compute_forces(.true.)
1123  CALL gather_array(gc)
1124 
1125  IF (iam .EQ. 0) THEN
1126  WRITE(*,1000) fsq_total, maxval(abs(gc)), sqrt(sum(xc*xc))
1127  END IF
1128 
1129 1000 FORMAT(/,' Back Solve Check', &
1130  /,' |F|^2: ',1p,e10.3,' MAX |F|: ', 1pe10.3, ' |X|: ',1pe10.3)
1131 
1132  END SUBROUTINE
1133 
1134 !-------------------------------------------------------------------------------
1136 !-------------------------------------------------------------------------------
1137  SUBROUTINE init_pchelms
1138  USE blocktridiagonalsolver_s, ONLY: initialize
1139  USE descriptor_mod, ONLY: iam, nprocs
1140  USE nscalingtools, ONLY: initremap
1141  USE shared_functions, ONLY: init_ptrs
1142  USE shared_data, ONLY: lasym
1143 
1144  IMPLICIT NONE
1145 
1146 ! Local variables
1147  INTEGER :: n1
1148  INTEGER :: istat
1149 
1150 ! Start of executable code
1151  IF (lasym) THEN
1152  ndims = 6
1153  ELSE
1154  ndims = 3
1155  END IF
1156 
1157  mblk_size = ndims*mnmax
1158  CALL initialize(ns, mblk_size)
1159 
1160  n1 = ns*mnmax
1161  neqs = ndims*n1
1162 
1163  CALL init_ptrs(xc, asubsmnsf, asubumncf, asubvmncf, &
1164  asubsmncf, asubumnsf, asubvmnsf)
1165 
1166  CALL init_ptrs(gc, fsupsmnsf, fsupumncf, fsupvmncf, &
1167  fsupsmncf, fsupumnsf, fsupvmnsf)
1168 
1169 ! Initialize forces
1170  CALL initremap(mpol, ntor, ns, nprocs, iam)
1171 
1172  CALL compute_forces(.true.)
1173  xc = 0
1174 
1175  END SUBROUTINE
1176 
1177 !-------------------------------------------------------------------------------
1182 !-------------------------------------------------------------------------------
1183  SUBROUTINE compute_forces(lNLinear)
1184  USE shared_data, ONLY: asubsmnsf, asubumncf, asubvmncf, &
1185  fsupsmnsf, fsupumncf, fsupvmncf, &
1186  asubsmncf, asubumnsf, asubvmnsf, &
1187  fsupsmncf, fsupumnsf, fsupvmnsf, neqs, &
1188  siesta_curtor, lasym
1189  USE bmw_run, ONLY: bmw_exec
1190  USE siesta_namelist, ONLY: wout_file, mgrid_file
1191  USE island_params, ONLY: fourier_context
1192  USE v3_utilities, ONLY: assert_eq
1193 
1194  IMPLICIT NONE
1195 
1196 ! Declare Arguments
1197  LOGICAL, INTENT(in) :: lNLinear
1198 
1199 ! Local variables
1200  INTEGER :: istat
1201  REAL (dp) :: dphi
1202  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: A_r
1203  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: A_p
1204  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: A_z
1205  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: cA_s
1206  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: cA_u
1207  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: cA_v
1208 
1209 ! Start of executable code
1210  nsmin = startglobrow
1211  nsmax = min(ns,endglobrow + 1)
1212  ns_match = 0 !nsp /4
1213 
1214 ! FIXME: All this initalizion should be somewhere else.
1215  IF (bnorm .LT. 0) THEN
1216  linit = .true.
1217  IF (.NOT. ALLOCATED(xc0)) THEN
1218  ALLOCATE(xc0(neqs), gc0(neqs))
1219  CALL assert(ALLOCATED(jbsupsmnsh), 'B^S not allocated!')
1220  IF (lasym) THEN
1221  CALL assert(ALLOCATED(jbsupsmnch), 'B^S not allocated!')
1222  END IF
1223  ALLOCATE(bsupsijh(ntheta,nzeta,nsmin:nsmax), &
1224  bsupuijh(ntheta,nzeta,nsmin:nsmax), &
1225  bsupvijh(ntheta,nzeta,nsmin:nsmax), &
1226  bsubsijh(ntheta,nzeta,nsmin:nsmax), &
1227  bsubuijh(ntheta,nzeta,nsmin:nsmax), &
1228  bsubvijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
1229  CALL assert_eq(0, istat, 'Allocation failed in COMPUTE_FORCES')
1230  END IF
1231 
1232  CALL init_f(fsupsmnsf, fsupumncf, fsupvmncf, &
1233  jcurrumnc, jcurrvmnc)
1234  IF (lasym) THEN
1235  CALL init_f(fsupsmncf, fsupumnsf, fsupvmnsf, &
1236  jcurrumns, jcurrvmns)
1237  END IF
1238  gc0 = gc
1239 
1240  CALL getfsq_pchelms
1241  bnorm = fsq_total
1242 
1243  dphi = (2*pi)/(nfp*nzeta)
1244  ALLOCATE(a_r(ntheta,nzeta,1), &
1245  a_p(ntheta,nzeta,1), &
1246  a_z(ntheta,nzeta,1))
1247 
1248  CALL bmw_exec(mgrid_file, wout_file, r1_i, z1_i, dphi, ns, &
1249  a_r, a_p, a_z)
1250 
1251  ALLOCATE(ca_s(ntheta,nzeta,1), &
1252  ca_u(ntheta,nzeta,1), &
1253  ca_v(ntheta,nzeta,1))
1254  CALL cyl2vmec_a(a_r(:,:,1), a_p(:,:,1), a_z(:,:,1), &
1255  ca_s(:,:,1), ca_u(:,:,1), ca_v(:,:,1))
1256  DEALLOCATE(a_r, a_p, a_z)
1257 
1258  CALL init_a(ca_s, ca_u, ca_v, asubsmnsf, asubumncf, asubvmncf, f_sin)
1259  IF (lasym) THEN
1260  CALL init_a(ca_s, ca_u, ca_v, asubsmncf, asubumnsf, asubvmnsf, &
1261  f_cos)
1262  END IF
1263  DEALLOCATE(ca_s, ca_u, ca_v)
1264 
1265  xc0 = xc
1266  END IF
1267 
1268  gc = 0
1269 
1270 ! Get sqrt(g)*CURLA.grad(s,u,v) [sqrt(g) B^i] on half grid in real space, for
1271 ! input AsubXmn on full grid.
1272  nsmin = max(1, startglobrow - 1)
1273  nsmax = min(ns, endglobrow + 1)
1274  CALL curla_pchelms(asubsmnsf, asubumncf, asubvmncf, &
1275  jbsupsmnsh, jbsupumnch, jbsupvmnch, f_sin)
1276  IF (lasym) THEN
1277  CALL curla_pchelms(asubsmncf, asubumnsf, asubvmnsf, &
1278  jbsupsmnch, jbsupumnsh, jbsupvmnsh, f_cos)
1279  END IF
1280 
1281 ! Add 1/jac factor
1282  nsmin = startglobrow
1283 
1284  bsupsijh = bsupsijh/jacobh(:,:,nsmin:nsmax)
1285  bsupuijh = bsupuijh/jacobh(:,:,nsmin:nsmax)
1286  bsupvijh = bsupvijh/jacobh(:,:,nsmin:nsmax)
1287 
1288 ! Convert to covariant component of B_i on half grid.
1289  CALL tolowerh(bsupsijh, bsupuijh, bsupvijh, &
1290  bsubsijh, bsubuijh, bsubvijh, nsmin, nsmax)
1291 
1292 ! Get -sqrt(g)*CURLB dot grad(s,u,v) Fourier coefficients [sqrt(g)*J^i] on full
1293 ! grid and store in FsupXmnf (sign is flipped to make positive contributions to
1294 ! Hessian)
1295  nsmin = startglobrow
1296  nsmax = min(ns, endglobrow + 1)
1297  CALL curlb_pchelms(fsupsmnsf, fsupumncf, fsupvmncf, asubsmnsf, &
1298  jbsupsmnsh, f_sin, siesta_curtor)
1299  IF (lasym) THEN
1300  CALL curlb_pchelms(fsupsmncf, fsupumnsf, fsupvmnsf, asubsmncf, &
1301  jbsupsmnch, f_cos, siesta_curtor)
1302  END IF
1303 
1304  IF (nsmax .eq. ns) THEN
1305  siesta_curtor = -fourier_context%orthonorm(m0,n0)*siesta_curtor
1306  END IF
1307 
1308  nsmax = endglobrow
1309 
1310 ! Dump debugging info.
1311 ! FIXME: This should only get written out if compiled in debug mode. Check
1312 ! linit to see if this is only written once.
1313  IF (linit) THEN
1314  CALL compare_current(fsupumncf, fsupvmncf, jcurrumnc, jcurrvmnc)
1315  IF (lasym) THEN
1316  CALL compare_current(fsupumnsf, fsupvmnsf, jcurrumns, jcurrvmns)
1317  END IF
1318  END IF
1319 
1320  IF (linit .OR. lnlinear) THEN
1321  gc = gc + gc0
1322  END IF
1323 
1324 ! Compute ONLY -Ax, NOT residue= -Ax + b (for preconditioner and GMRES)
1325 ! Subtract [sqrt(g)*J]mn (input from VMEC) to get net force (residual=b-Ax)
1326  CALL boundaryconditions(fsupsmnsf, fsupumncf, fsupvmncf, f_sin)
1327  IF (lasym) THEN
1328  CALL boundaryconditions(fsupsmncf, fsupumnsf, fsupvmnsf, f_cos)
1329  END IF
1330 
1331 ! Subtract A*x0 since A_u, A_v are fixed at the edge b' = b + A*x0 (new
1332 ! right-hand side)
1333  IF (linit) THEN
1334  gc0 = gc
1335  END IF
1336 
1337  IF (lnlinear .OR. linit) THEN
1338  linit = .false.
1339  CALL getfsq_pchelms
1340  END IF
1341 
1342  END SUBROUTINE
1343 
1344 !-------------------------------------------------------------------------------
1351 !-------------------------------------------------------------------------------
1352  SUBROUTINE boundaryconditions(Fsupsmn, Fsupumn, Fsupvmn, iparity)
1353  USE hessian, ONLY: l_compute_hessian
1354  USE utilities, ONLY: set_bndy_fouier_m0, set_bndy_full_origin
1355 
1356  IMPLICIT NONE
1357 
1358 ! Declare Arguments
1359  REAL(dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(inout) :: Fsupsmn
1360  REAL(dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(inout) :: Fsupumn
1361  REAL(dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(inout) :: Fsupvmn
1362  INTEGER, INTENT(in) :: iparity
1363 
1364 ! Local variables
1365  INTEGER :: nmatch
1366 
1367 ! Start of executable code
1368  IF (endglobrow .EQ. ns) THEN
1369  IF (.not.l_asedge) THEN
1370  fsupsmn(:,:,ns) = 0
1371  END IF
1372  fsupumn(:,:,ns) = 0
1373  fsupvmn(:,:,ns) = 0
1374  END IF
1375 
1376  CALL set_bndy_fouier_m0(fsupsmn, fsupumn, fsupvmn, iparity)
1377  IF (startglobrow .eq. 1) THEN
1378  CALL set_bndy_full_origin(fsupsmn, fsupumn, fsupvmn)
1379  END IF
1380 
1381 ! Preserve surfaces for js < ns_match
1382  IF (ns_match .GE. nsmin) THEN
1383  nmatch = min(ns_match, nsmax)
1384  fsupumn(m1,:,nsmin:nmatch) = 0
1385  fsupvmn(m1,:,nsmin:nmatch) = 0
1386  END IF
1387 
1388  END SUBROUTINE
1389 
1390 !-------------------------------------------------------------------------------
1394 !-------------------------------------------------------------------------------
1395  SUBROUTINE gmres_pchelms
1396  USE descriptor_mod, ONLY: siesta_comm, iam, nprocs
1397  USE gmres_lib, ONLY: gmres_info, gmres_par, gmres_ser
1398  USE hessian, ONLY: apply_precond
1399  USE shared_data, ONLY: lverbose, unit_out, gmres_no_peak, iortho
1400 
1401  IMPLICIT NONE
1402 
1403 ! Local variables
1404  TYPE (gmres_info) :: gi
1405  INTEGER :: n
1406  INTEGER :: m
1407  INTEGER :: js
1408  REAL (dp), DIMENSION(:), ALLOCATABLE :: b
1409  REAL (dp), DIMENSION(:), ALLOCATABLE :: x0
1410 
1411 ! Local Parameters
1412  INTEGER, PARAMETER :: noPrec = 0
1413  integer, PARAMETER :: leftprec = 1
1414  INTEGER, PARAMETER :: rightPrec = 2
1415  INTEGER, PARAMETER :: dblePrec = 3
1416  REAL (dp), PARAMETER :: one = 1
1417  REAL (dp), PARAMETER :: ftol = 1.e-18_dp
1418 
1419 ! Start of executable code
1420  ALLOCATE(b(neqs), x0(neqs))
1421 
1422  xc = 0
1423 
1424  CALL compute_forces(.true.)
1425  IF (iam .EQ. 0) THEN
1426  IF (lverbose) THEN
1427  WRITE (*,*)
1428  WRITE (*,1001) 1, fsq_total, sqrt(sum(xc0**2))
1429  END IF
1430  WRITE (unit_out, 1000) neqs
1431  WRITE (unit_out, 1001) 1, fsq_total, sqrt(sum(xc0**2))
1432  END IF
1433 
1434  b = gc0
1435  x0 = xc
1436  x0 = 0
1437 
1438  CALL init_dgmres(gi%icntl, gi%cntl)
1439 
1440 ! Tune GMRES parameters
1441 ! Tolerance
1442  gi%cntl(1) = sqrt(ftol)
1443 ! Write warnings to fort.21
1444  gi%icntl(2) = 21
1445 ! Save the convergence history in file fort.20
1446  gi%icntl(3) = 20
1447 ! Preconditioning input flags (note: different than revcom flags in driver)
1448  gi%icntl(4) = rightprec ! leftPrec, dblePrec, noPrec
1449 ! Orthogonalization type.
1450  gi%icntl(5) = iortho
1451 ! Initial guess (use it if = 1, otherwise ignore it)
1452  gi%icntl(6) = 1
1453 ! Maximum number of iterations at each step (~ns/5)
1454  gi%icntl(7) = 20
1455 ! Default
1456  gi%icntl(8) = 1
1457 ! Steps for peek at progress during rev com loop.
1458  gi%icntl(9) = 1
1459 ! linear solver (output on screen).
1460  gi%l_nonlinear = .false.
1461 
1462 
1463  n = neqs
1464  gi%m = 200
1465 
1466 ! SIESTA process mpi rank.
1467  gi%iam = iam
1468 ! Number of siesta processes.
1469  gi%nprocs=nprocs
1470  gi%ngmres_type = gmres_no_peak
1471  gi%ftol = ftol
1472  IF (parsolver) THEN
1473  gi%endglobrow = endglobrow
1474  gi%startglobrow = startglobrow
1475  gi%mblk_size = mblk_size
1476  gi%rcounts => rcounts
1477  gi%disp => disp
1478 
1479  gi%my_comm = siesta_comm
1480  gi%my_comm_world = siesta_comm
1481 
1482  CALL gmres_par(n, gi, matvec_par_pchelms, apply_precond, &
1483  getnlforce_pchelms, x0, b)
1484  ELSE
1485  CALL gmres_ser(n, gi, matvec_pchelms, apply_precond, &
1486  getnlforce_pchelms, x0, b)
1487  END IF
1488 
1489  xc = x0
1490 
1491  DEALLOCATE(x0, b)
1492 
1493 1000 FORMAT(' NEQS : ',i6)
1494 1001 FORMAT(' ITER: ',i6,3x,'|F|^2: ',1pe12.3,3x,' |X| = ',1pe12.3)
1495 
1496  END SUBROUTINE
1497 
1498 !-------------------------------------------------------------------------------
1504 !-------------------------------------------------------------------------------
1505  SUBROUTINE matvec_pchelms (p, Ap, ndim)
1506  USE v3_utilities, ONLY: assert_eq
1507 
1508  IMPLICIT NONE
1509 
1510 ! Declare Arguments
1511  REAL(dp), INTENT(IN) :: p(ndim)
1512  REAL(dp), INTENT(OUT) :: Ap(ndim)
1513  INTEGER, INTENT(IN) :: ndim
1514 
1515 ! Start of executable code
1516  CALL assert_eq(neqs, ndim, ' neqs != ndim')
1517 
1518  xc = p
1519  CALL compute_forces(.false.)
1520  ap = -gc
1521 
1522  END SUBROUTINE
1523 
1524 !-------------------------------------------------------------------------------
1530 !-------------------------------------------------------------------------------
1531  SUBROUTINE matvec_par_pchelms (ploc, Ap, nloc)
1532  USE hessian, ONLY: eps_factor
1533  USE blocktridiagonalsolver_s, ONLY: parmatvec, padsides
1534  USE v3_utilities, ONLY: assert_eq
1535 
1536  IMPLICIT NONE
1537 
1538 ! Declare Arguments
1539  REAL(dp), DIMENSION(nloc), INTENT(IN) :: ploc
1540  REAL(dp), DIMENSION(nloc), INTENT(OUT) :: Ap
1541  INTEGER, INTENT(IN) :: nloc
1542 
1543 ! local variables
1544  INTEGER :: myrowstart
1545  INTEGER :: myrowend
1546  INTEGER :: istat
1547  REAL (dp), DIMENSION(:), ALLOCATABLE :: p
1548 
1549 ! Start of executable code
1550  istat = (endglobrow - startglobrow + 1)*mblk_size
1551  CALL assert_eq(istat, nloc, 'nloc wrong in matvec_par')
1552 
1553  myrowstart = (startglobrow - 1)*mblk_size + 1
1554  myrowend = myrowstart + nloc - 1
1555 
1556  ALLOCATE (p(ns*mblk_size), stat=istat)
1557  CALL assert_eq(0, istat, 'Allocation error in matvec_par')
1558 
1559  p(myrowstart:myrowend) = ploc
1560  CALL padsides(p, mblk_size, 1, 1)
1561 
1562  xc = p
1563  CALL compute_forces(.false.)
1564  ap = -gc(myrowstart:myrowend)
1565 
1566  DEALLOCATE (p)
1567 
1568  END SUBROUTINE
1569 
1570 !-------------------------------------------------------------------------------
1576 !-------------------------------------------------------------------------------
1577  SUBROUTINE getnlforce_pchelms(xcstate, fsq_nl, bnorm)
1578 
1579  IMPLICIT NONE
1580 
1581 ! Declare Arguments
1582  REAL (dp), DIMENSION(neqs), INTENT(IN) :: xcstate
1583  REAL (dp), INTENT(OUT) :: fsq_nl
1584  REAL (dp), INTENT(IN) :: bnorm
1585 
1586 ! Start of executable code
1587  WRITE (*,1000)
1588  fsq_nl = 1
1589  RETURN
1590 
1591 1000 FORMAT('PCHELMS should not get here.')
1592 
1593  END SUBROUTINE
1594 
1595  END MODULE pchelms
shared_data::r1_i
real(dp), dimension(:,:,:), allocatable r1_i
R coordinates of the computational grid.
Definition: shared_data.f90:153
shared_data::lverbose
logical lverbose
Use verbose screen output.
Definition: shared_data.f90:234
shared_data::zv_i
real(dp), dimension(:,:,:), allocatable zv_i
dZ/dv coordinates of the computational grid.
Definition: shared_data.f90:163
shared_functions
Module contained subroutines and functions updating MHD forces and Wmhd.
Definition: shared_functions.f90:4
shared_data::xc
real(dp), dimension(:), allocatable xc
1D array of Fourier mode displacement components.
Definition: shared_data.f90:97
shared_data::iortho
integer iortho
Orthogonalization in GMRES.
Definition: shared_data.f90:81
hessian::gather_array
Definition: hessian.f90:31
shared_data::fsq_total
real(dp) fsq_total
|F|^2 WITH column scaling.
Definition: shared_data.f90:91
quantities
This file contains subroutines for allocating and initializing curvilinear magnetic covariant and pre...
Definition: quantities.f90:11
bmw_run
BMW is a code for extending fields belond the VMEC domain in a manner that ensures divergence free fi...
Definition: bmw_run.f:27
v3_utilities::assert_eq
Definition: v3_utilities.f:62
pchelms::curlb_pchelms
subroutine curlb_pchelms(ksupsmnf, ksupumnf, ksupvmnf, asubsmnf,
Compute real-space components of contravariant jac*J^i on full radial grid from curl(B).
Definition: pchelms.f90:285
shared_data::siesta_curtor
real(dp) siesta_curtor
Total toroidal current.
Definition: shared_data.f90:239
shared_data::torflux
real(dp), dimension(:), allocatable torflux
Toroidal flux profile.
Definition: shared_data.f90:169
shared_data::asubvmnsf
real(dp), dimension(:,:,:), pointer asubvmnsf
Covariant vector potential for non-stellator symmetric v component on full grid.
Definition: shared_data.f90:186
shared_data::asubumncf
real(dp), dimension(:,:,:), pointer asubumncf
Covariant vector potential for stellator symmetric u component on full grid.
Definition: shared_data.f90:178
island_params::phipf_i
real(dp), dimension(:), allocatable phipf_i
Radial toroidal flux derivative.
Definition: island_params.f90:68
shared_data::fsupumncf
real(dp), dimension(:,:,:), pointer fsupumncf
Contravariant force for stellarator symmetric u component on full grid.
Definition: shared_data.f90:192
island_params::fourier_context
type(fourier_class), pointer fourier_context
Fourier transform object.
Definition: island_params.f90:76
shared_data::fsupsmnsf
real(dp), dimension(:,:,:), pointer fsupsmnsf
Contravariant force for stellarator symmetric s component on full grid.
Definition: shared_data.f90:188
fourier
Module contains subroutines for computing FFT on parallel radial intervals. Converts quantities betwe...
Definition: fourier.f90:13
pchelms::gmres_pchelms
subroutine gmres_pchelms
Setup and run GMRES solver for the helmholtz problem.
Definition: pchelms.f90:1396
fourier::m0
integer, parameter m0
m = 0 mode.
Definition: fourier.f90:58
siesta_namelist::mgrid_file
character(len=siesta_namelist_name_length) mgrid_file
Filename of the VMEC woutfile.
Definition: siesta_namelist.f90:193
shared_data::asubsmncf
real(dp), dimension(:,:,:), pointer asubsmncf
Covariant vector potential for non-stellator symmetric s component on full grid.
Definition: shared_data.f90:176
v3_utilities::assert
Definition: v3_utilities.f:55
fourier::f_sin
integer, parameter f_sin
Sine parity.
Definition: fourier.f90:27
shared_data::asubsmnsf
real(dp), dimension(:,:,:), pointer asubsmnsf
Covariant vector potential for stellator symmetric s component on full grid.
Definition: shared_data.f90:173
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
pchelms::l_asedge
logical, parameter, private l_asedge
Controls if the edge Aubs values are evolved or fixed.
Definition: pchelms.f90:46
siesta_namelist::wout_file
character(len=siesta_namelist_name_length) wout_file
Filename of the VMEC woutfile.
Definition: siesta_namelist.f90:189
pchelms::dump_a
subroutine dump_a(js, iunit)
Write out vector potential to file.
Definition: pchelms.f90:707
shared_data::lcolscale
logical lcolscale
Apply column scaling to hessian.
Definition: shared_data.f90:228
shared_data::fsupumnsf
real(dp), dimension(:,:,:), pointer fsupumnsf
Contravariant force for stellarator non-symmetric u component on full grid.
Definition: shared_data.f90:194
shared_data::gmres_no_peak
integer, parameter gmres_no_peak
GMRES no_peak improvement.
Definition: shared_data.f90:35
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
shared_data::zu_i
real(dp), dimension(:,:,:), allocatable zu_i
dZ/du coordinates of the computational grid.
Definition: shared_data.f90:159
pchelms::compare_current
subroutine compare_current(Fsupumn, Fsupvmn, jcurrumn, jcurrvmn)
Compare Curl(Curl(A)) with the expected VMEC and vacuum currents.
Definition: pchelms.f90:650
shared_data::polflux
real(dp), dimension(:), allocatable polflux
Poloidal flux profile.
Definition: shared_data.f90:171
pchelms::cyl2vmec_a
subroutine cyl2vmec_a(A_r, A_p, A_z, cA_s, cA_u, cA_v)
Convert cylindical vector potential to contravariant components.
Definition: pchelms.f90:441
fourier::n0
integer, parameter n0
n = 0 mode.
Definition: fourier.f90:65
gmres_lib::gmres_info
Definition: gmres_lib.f:7
pchelms::init_a
subroutine init_a(cA_s, cA_u, cA_v, A_s, A_u, A_v, parity)
Initialize vector potential.
Definition: pchelms.f90:495
pchelms::getnlforce_pchelms
subroutine getnlforce_pchelms(xcstate, fsq_nl, bnorm)
Non linear force callback.
Definition: pchelms.f90:1578
pchelms::run_pchelms
subroutine run_pchelms
Run the pchelms code to solve for inital jbsup values.
Definition: pchelms.f90:77
shared_data::fsupsmncf
real(dp), dimension(:,:,:), pointer fsupsmncf
Contravariant force for stellarator non-symmetric s component on full grid.
Definition: shared_data.f90:190
shared_data::unit_out
integer, parameter unit_out
File output io unit.
Definition: shared_data.f90:39
shared_data::lasym
logical lasym
Use non-stellarator symmetry.
Definition: shared_data.f90:230
island_params::hs_i
real(dp) hs_i
Radial grid spacing. hs = s_i+1 - s_i.
Definition: island_params.f90:47
blocktridiagonalsolver_s
Module contains subroutines for solver for block tri-diagonal matrices.
Definition: blocktridiagonalsolver_s.f90:7
siesta_namelist::nsin
integer nsin
Radial size of the plasma grid.
Definition: siesta_namelist.f90:162
shared_data::gc
real(dp), dimension(:), allocatable, target gc
1D Array of Fourier mode MHD force components
Definition: shared_data.f90:99
island_params::chipf_i
real(dp), dimension(:), allocatable chipf_i
Radial poloidal flux derivative.
Definition: island_params.f90:70
shared_data::asubvmncf
real(dp), dimension(:,:,:), pointer asubvmncf
Covariant vector potential for stellator symmetric v component on full grid.
Definition: shared_data.f90:183
shared_data::rv_i
real(dp), dimension(:,:,:), allocatable rv_i
dR/dv coordinates of the computational grid.
Definition: shared_data.f90:161
pchelms
This file solves the helmholtz equation to set inital fields that match vmec and vacuum currents from...
Definition: pchelms.f90:12
shared_data::nsp
integer nsp
Total radial grid size in the VMEC region.
Definition: shared_data.f90:64
shared_data::xc0
real(dp), dimension(:), allocatable xc0
Saved fouier displacements.
Definition: shared_data.f90:106
pchelms::boundaryconditions
subroutine boundaryconditions(Fsupsmn, Fsupumn, Fsupvmn, iparity)
Apply boundary conditions.
Definition: pchelms.f90:1353
metrics::tolowerh
subroutine tolowerh(xsupsij, xsupuij, xsupvij,
Converts to covariant component half grid.
Definition: metrics.f90:539
fourier::f_sum
integer, parameter f_sum
Sum fouier real space transformed quantity. This is used when a non-stellarator symmetric parity is b...
Definition: fourier.f90:43
fourier::m1
integer, parameter m1
m = 1 mode.
Definition: fourier.f90:60
pchelms::matvec_pchelms
subroutine matvec_pchelms(p, Ap, ndim)
Serial callback function for MatVec GMRES routine.
Definition: pchelms.f90:1506
shared_data::ndims
integer ndims
Number of independent variables.
Definition: shared_data.f90:58
shared_data::fsupvmncf
real(dp), dimension(:,:,:), pointer fsupvmncf
Contravariant force for stellarator symmetric v component on full grid.
Definition: shared_data.f90:196
island_params::ohs_i
real(dp) ohs_i
Radial grid derivative factor. d X/ ds = (x_i+1 - x_i)/(s_i+1 - s_i) (1) Where ohs = 1/(s_i+1 - s_i) ...
Definition: island_params.f90:51
shared_data::mblk_size
integer mblk_size
Block size. (mpol + 1)*(2*ntor + 1)*ndims.
Definition: shared_data.f90:62
shared_data::z1_i
real(dp), dimension(:,:,:), allocatable z1_i
Z coordinates of the computational grid.
Definition: shared_data.f90:155
metrics
Module for reading in VMEC input and computing metric elements on the SIESTA (sqrt-flux) radial grid....
Definition: metrics.f90:13
fourier::f_cos
integer, parameter f_cos
Cosine parity.
Definition: fourier.f90:25
siesta_namelist
This file contains all the variables and maximum sizes of the inputs for a SIESTA namelist input file...
Definition: siesta_namelist.f90:103
shared_data::ru_i
real(dp), dimension(:,:,:), allocatable ru_i
dR/du coordinates of the computational grid.
Definition: shared_data.f90:157
shared_data
This file contains variables and parameters used by many modules in SIESTA.
Definition: shared_data.f90:10
shared_data::asubumnsf
real(dp), dimension(:,:,:), pointer asubumnsf
Covariant vector potential for non-stellator symmetric u component on full grid.
Definition: shared_data.f90:181
shared_data::fsupvmnsf
real(dp), dimension(:,:,:), pointer fsupvmnsf
Contravariant force for stellarator non-symmetric v component on full grid.
Definition: shared_data.f90:198
pchelms::init_f
subroutine init_f(Fsupsmn, Fsupumn, Fsupvmn, jcurrumn, jcurrvmn)
Initialize expected currents.
Definition: pchelms.f90:586
shared_data::neqs
integer neqs
Number of elements in the xc array.
Definition: shared_data.f90:54
shared_data::gc0
real(dp), dimension(:), allocatable gc0
Saved fouier MHD forces.
Definition: shared_data.f90:108
pchelms::matvec_par_pchelms
subroutine matvec_par_pchelms(ploc, Ap, nloc)
Parallel callback function for MatVec GMRES routine.
Definition: pchelms.f90:1532
pchelms::curla_pchelms
subroutine curla_pchelms(Asubsmnf, Asubumnf, Asubvmnf,
Compute real-space components of contravariant jac*B^i on half radial grid from curl(A).
Definition: pchelms.f90:187
fourier::f_none
integer, parameter f_none
Do not sum fouier real space transformed quantity. This is used when a stellarator symmetric parity i...
Definition: fourier.f90:32
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17