V3FIT
quantities.f90
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
10 !*******************************************************************************
11  MODULE quantities
12  USE stel_kinds
13  USE stel_constants
14  USE v3_utilities, ONLY: assert, assert_eq
15  USE island_params, ns=>ns_i, ntheta=>nu_i, nzeta=>nv_i, &
16  mpol=>mpol_i, ntor=>ntor_i, nuv=>nuv_i, mnmax=>mnmax_i, &
17  nfp=>nfp_i
18  USE descriptor_mod, ONLY: iam, nprocs, siesta_comm
19  USE shared_data, ONLY: lasym, lverbose
20  USE nscalingtools, ONLY: startglobrow, endglobrow, parsolver, mpi_err
21  USE utilities, ONLY: to_half_mesh, to_full_mesh
22  USE siesta_namelist, ONLY: l_vessel
23 
24  IMPLICIT NONE
25 
26 !*******************************************************************************
27 ! Module variables
28 !*******************************************************************************
30  REAL (dp) :: b_factor
32  REAL (dp) :: p_factor
34  REAL (dp) :: signjac
36  REAL (dp) :: wp
38  REAL (dp) :: wb
40  REAL (dp) :: wp0
41 
44  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jbsupsmnsh
47  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jbsupumnch
50  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jbsupvmnch
53  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jpmnch
56  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupsmnsf
59  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupumncf
62  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupvmncf
65  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: djpmnch
68  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: djbsupsmnsh
71  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: djbsupumnch
74  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: djbsupvmnch
75 
78  REAL (dp), POINTER, DIMENSION(:,:,:) :: jvsupsmncf
81  REAL (dp), POINTER, DIMENSION(:,:,:) :: jvsupumnsf
84  REAL (dp), POINTER, DIMENSION(:,:,:) :: jvsupvmnsf
87  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsubsmncf
90  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsubumnsf
93  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsubvmnsf
96  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsupsmncf
99  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsupumnsf
102  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsupvmnsf
103 
106  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jbsupsmnch
109  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jbsupumnsh
112  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jbsupvmnsh
115  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jpmnsh
118  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupsmncf
121  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupumnsf
124  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupvmnsf
127  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: djpmnsh
130  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: djbsupsmnch
133  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: djbsupumnsh
136  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: djbsupvmnsh
137 
140  REAL (dp), POINTER, DIMENSION(:,:,:) :: jvsupsmnsf
143  REAL (dp), POINTER, DIMENSION(:,:,:) :: jvsupumncf
146  REAL (dp), POINTER, DIMENSION(:,:,:) :: jvsupvmncf
149  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsubsmnsf
152  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsubumncf
155  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsubvmncf
158  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsupsmnsf
161  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsupumncf
164  REAL (dp), POINTER, DIMENSION(:,:,:) :: fsupvmncf
165 
167  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: pwr_spec_s
169  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: pwr_spec_a
170 
172  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jvsupsijf
174  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jvsupuijf
176  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jvsupvijf
178  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jacobh
180  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: jacobf
182  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: wint
183 
185  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupsijf0
187  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupuijf0
189  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupvijf0
192  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupsijh0
195  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupuijh0
198  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupvijh0
201  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupsijf
204  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupuijf
207  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupvijf
209  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubsijf
211  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubuijf
213  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubvijf
215  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsq
217  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksubsijf
219  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksubuijf
221  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksubvijf
223  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupsijf0
225  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupuijf0
227  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupvijf0
229  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupsijf
231  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupuijf
233  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: ksupvijf
234 
236  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: pijh0
238  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: pijh0_du
240  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: pijh0_dv
242  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: pijf0
244  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: pijf0_ds
245 
247  REAL (dp) :: fbdy(13)
248 
249  CONTAINS
250 
251 !*******************************************************************************
252 ! UTILITY SUBROUTINES
253 !*******************************************************************************
254 !-------------------------------------------------------------------------------
256 !-------------------------------------------------------------------------------
257  SUBROUTINE init_quantities
258  USE descriptor_mod, ONLY: lscalapack
259  USE island_params, ONLY: fourier_context
260  USE fourier, ONLY: f_none, f_cos, f_sin, f_sum, m0
261  USE metrics, ONLY: sqrtg
262  USE shared_data, ONLY: mblk_size, ndims, lasym
263 #if defined(MPI_OPT)
264  USE prof_mod, ONLY: setupscalingallgather
265 #endif
266  USE blocktridiagonalsolver_s, ONLY: initialize
267  USE siesta_namelist, ONLY: lrestart
268 
269  IMPLICIT NONE
270 
271 ! Local variables.
272  INTEGER :: istat
273  INTEGER :: l
274  REAL (dp) :: sum1
275  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: tempmn
276 
277 ! Start of executable code.
278  IF (parsolver) THEN
279  CALL initialize(ns, mblk_size)
280  ELSE IF (lscalapack) THEN
281 #if defined(MPI_OPT)
282  CALL setupscalingallgather(mblk_size)
283 #endif
284  ELSE
285  startglobrow = 1
286  endglobrow = ns
287  END IF
288 
289  IF (.not.lrestart) THEN
290 ! Allocate and initialize dependent variables
291  CALL alloc_quantities
292 
293  p_factor = 1.0_dp/abs(wb_i)
294  b_factor = sqrt(p_factor)
295  gnorm_i = abs(wb_i)/(wb_i + wp_i/(gamma - 1.0))
296  END IF
297 
298 ! Reshapes sqrtg in local form
299  jacobh(:,:,:) = reshape(sqrtg, shape(jacobh))
300 
301 !AVOID DIVIDE BY ZERO AND STORE AVERAGED VALUE (OVER THETA)
302 !AT FIRST 1/2 GRID PT AT ORIGIN IN JACOBH(1)
303  CALL fourier_context%get_origin(jacobh, m0, lasym)
304  signjac = jacobh(1,1,2)/abs(jacobh(1,1,2))
305 
306 ! Filter jacobian ro preserve the unperturbed (constant in s) part of the
307 ! pressure.
308  ALLOCATE(tempmn(0:mpol,-ntor:ntor,SIZE(jacobh,3)))
309  CALL fourier_context%tomnsp(jacobh, tempmn, f_cos)
310  CALL fourier_context%toijsp(tempmn, jacobh, f_none, f_cos)
311  IF (lasym) THEN
312  CALL fourier_context%tomnsp(jacobh, tempmn, f_sin)
313  CALL fourier_context%toijsp(tempmn, jacobh, f_sum, f_sin)
314  END IF
315  DEALLOCATE(tempmn)
316 
317  CALL assert(all(jacobh*signjac .gt. 0), 'FILTERED JACOBIAN CHANGED SIGN!')
318 
319  IF (iam .eq. 0 .and. lverbose) THEN
320  sum1 = sum((jacobh(:,:,2:) - jacobf(:,:,2:))**2 / &
321  (jacobh(:,:,2:) + jacobf(:,:,2:))**2)
322  WRITE (*,1000) sqrt(sum1/SIZE(jacobh))
323  END IF
324 
325  CALL to_full_mesh(jacobh, jacobf)
326 ! Must be consistent with variational form of pressure evolution equation. Do
327 ! not use 2pt extrapoltaion.
328  jacobf(:,:,1) = jacobh(:,:,1)
329  jacobf(:,:,ns)= jacobh(:,:,ns)
330 
331  DO l = 1, ntheta
332  wint(l,:,:) = fourier_context%cosmui(l,m0)
333  END DO
334 
335  ALLOCATE (vp_f(ns), stat=istat)
336  CALL assert_eq(0, istat, 'Allocation error in init_quantities')
337  CALL surfaverage(vp_f, jacobf, 1, ns)
338 
339 ! Compar volumes for a sanity check that the siesta grid and orginal VMEC
340 ! volume match.
341  sum1 = 0.5*hs_i*(sum(vp_f(2:ns)) + sum(vp_f(1:ns - 1)))*signjac
342  sum1 = twopi*twopi*sum1
343  IF (abs((sum1 - volume_i)/volume_i) .GT. 1.e-3_dp) THEN
344  IF (iam .EQ. 0 .and. lverbose) THEN
345  WRITE (*,1001) volume_i, sum1
346  END IF
347  END IF
348 
349 1000 FORMAT(' JACOBIAN SPECTRAL TRUNCATION ERROR: ',1pe8.2)
350 1001 FORMAT(' VMEC VOLUME: ', 1pe8.2, ' SIESTA VOLUME: ',1pe8.2)
351 
352  END SUBROUTINE
353 
354 !-------------------------------------------------------------------------------
359 !-------------------------------------------------------------------------------
360  SUBROUTINE init_fields
361  USE island_params, ONLY: fourier_context
362  USE fourier, ONLY: f_sin, f_cos, m0, m1
363  USE vmec_info, ONLY: lmns_i, lmnc_i, iflipj
364  USE metrics, ONLY: sqrtg
365 
366  IMPLICIT NONE
367 
368 ! Local variables.
369  INTEGER :: istat
370  INTEGER :: nsmin
371  INTEGER :: nsmax
372  INTEGER :: nloc
373  INTEGER :: js
374  REAL(dp), DIMENSION(:), ALLOCATABLE :: phiph
375  REAL(dp), DIMENSION(:), ALLOCATABLE :: chiph
376  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: presif
377  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: presih
378 
379 ! Start of executable code.
380  nsmin = max(1, startglobrow)
381  nsmax = min(endglobrow, ns)
382 
383  ALLOCATE(phiph(ns), chiph(ns))
384 
385  phipf_i = signjac*iflipj*phipf_i
386  chipf_i = signjac*iflipj*chipf_i
387 
388  phiph(2:ns) = (phipf_i(2:ns) + phipf_i(1:nsh))/2
389  chiph(2:ns) = (chipf_i(2:ns) + chipf_i(1:nsh))/2
390  phiph(1) = 0
391  chiph(1) = 0
392 
393  IF (.not.l_vessel) THEN
394  CALL recompute_lambda(lmns_i, lmnc_i, jacobh, &
395  fourier_context%orthonorm, phiph, chiph, &
396  nsmin, nsmax)
397  END IF
398 
399  CALL init_bfield(jbsupsmnsh, jbsupumnch, jbsupvmnch, &
400  lmns_i, phiph, chiph, nsmin, nsmax, f_sin)
401  IF (lasym) THEN
402  CALL init_bfield(jbsupsmnch, jbsupumnsh, jbsupvmnsh, &
403  lmnc_i, phiph, chiph, nsmin, nsmax, f_cos)
404  END IF
405 
406  DEALLOCATE (phiph, chiph, stat=istat)
407  CALL assert(istat.EQ.0,'Deallocate error #1 in init_fields')
408 
409 ! Initialize half mesh pressure
410  nloc = max(1,startglobrow - 1)
411  ALLOCATE(presih(ntheta,nzeta,nloc:nsmax), &
412  presif(ntheta,nzeta,nloc:nsmax), stat=istat)
413  CALL assert(istat.EQ.0, 'Allocate error #1 in init_fields')
414 
415  DO js = nloc, nsmax
416  presif(:,:,js) = p_factor*presf_i(js)
417  END DO
418 
419  CALL to_half_mesh(presif, presih)
420 
421 ! Initialize (internally) jacob*(real pressure), which is evolved
422  presih(:,:,nsmin:nsmax) = &
423  & presih(:,:,nsmin:nsmax)*jacobh(:,:,nsmin:nsmax)
424 
425 ! Initialize internal pressure harmonics (half mesh)
426  CALL fourier_context%tomnsp(presih(:,:,nsmin:nsmax), &
427  jpmnch(:,:,nsmin:nsmax), f_cos)
428  jpmnch(:,:,1) = 0
429  jpmnch(m0,:,1) = jpmnch(m0,:,2)
430  IF (lasym) THEN
431  CALL fourier_context%tomnsp(presih(:,:,nsmin:nsmax), &
432  jpmnsh(:,:,nsmin:nsmax), f_sin)
433  jpmnsh(:,:,1) = 0
434  jpmnsh(m0,:,1) = jpmnsh(m0,:,2)
435  END IF
436 
437  DEALLOCATE(presif, presih, stat=istat)
438  CALL assert_eq(istat, 0, 'Deallocate error #2 in init_fields')
439 
440 ! Gather jbsupmnX's, jpmn onto all processors. This same as in update_state but
441 ! for initial values.
442  CALL gatherfields
443 
444  END SUBROUTINE
445 
446 !-------------------------------------------------------------------------------
458 !-------------------------------------------------------------------------------
459  SUBROUTINE init_bfield(jbsupsmnh, jbsupumnh, jbsupvmnh, &
460  lmn, phiph, chiph, nsmin, nsmax, parity)
461  USE fourier, ONLY: m0, m1, n0, f_sin
462  USE siesta_namelist, ONLY: l_vessel
463  USE utilities, ONLY: set_bndy_fouier_m0
464  USE island_params, ONLY: fourier_context
465 
466  IMPLICIT NONE
467 
468 ! Declare Arguments
469  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: jbsupsmnh
470  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: jbsupumnh
471  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: jbsupvmnh
472  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(in) :: lmn
473  REAL (dp), DIMENSION(:), INTENT(in) :: phiph
474  REAL (dp), DIMENSION(:), INTENT(in) :: chiph
475  INTEGER, INTENT(in) :: nsmin
476  INTEGER, INTENT(in) :: nsmax
477  INTEGER, INTENT(in) :: parity
478 
479 ! local variables
480  INTEGER :: m
481  INTEGER :: n
482  INTEGER :: i
483  INTEGER :: sparity
484  INTEGER :: mp
485  INTEGER :: np
486 
487 ! Start of executable code.
488  IF (parity .eq. f_sin) THEN
489  sparity = 1
490  ELSE
491  sparity = -1
492  END IF
493 
494  IF (.not.l_vessel) THEN
495 
496  DO i = nsmin, nsmax
497  DO n = -ntor, ntor
498  np = n*nfp*sparity
499  DO m = 0, mpol
500  mp = m*sparity
501  jbsupumnh(m,n,i) = phiph(i)*(-np*lmn(m,n,i))
502  jbsupvmnh(m,n,i) = phiph(i)*( mp*lmn(m,n,i))
503  END DO
504  END DO
505 
506  IF (parity .eq. f_sin) THEN
507  jbsupumnh(m0,n0,i) = chiph(i) &
508  / fourier_context%orthonorm(m0,n0)
509  jbsupvmnh(m0,n0,i) = phiph(i) &
510  / fourier_context%orthonorm(m0,n0)
511  END IF
512  END DO
513 
514  jbsupsmnh = 0
515  END IF
516 
517  jbsupsmnh(:,:,nsmin:nsmax) = b_factor*jbsupsmnh(:,:,nsmin:nsmax)
518  jbsupumnh(:,:,nsmin:nsmax) = b_factor*jbsupumnh(:,:,nsmin:nsmax)
519  jbsupvmnh(:,:,nsmin:nsmax) = b_factor*jbsupvmnh(:,:,nsmin:nsmax)
520 
521  CALL set_bndy_fouier_m0(jbsupsmnh, jbsupumnh, jbsupvmnh, parity)
522 
523 ! Store m=0 component of jbsupvmnh (and m=0,1 of jbsupumnh) at origin of
524 ! half-grid js=1.
525  IF (nsmin .eq. 1) THEN
526  jbsupsmnh(:,:,1) = 0
527  jbsupsmnh(m1,:,1) = jbsupsmnh(m1,:,2)
528  jbsupumnh(:,:,1) = 0
529  jbsupumnh(m0:m1,:,1) = jbsupumnh(m0:m1,:,2)
530  jbsupvmnh(m0,:,1) = jbsupvmnh(m0,:,2)
531  jbsupvmnh(m1:,:,1) = 0
532  END IF
533 
534  END SUBROUTINE
535 
536 !-------------------------------------------------------------------------------
578 !-------------------------------------------------------------------------------
579  SUBROUTINE recompute_lambda(lmns, lmnc, jacobh, orthonorm, &
580  phiph, chiph, nsmin, nsmax)
581  USE island_params, ns=>ns_i, nuv=>nuv_i, mnmax=>mnmax_i, &
582  nfp => nfp_i, ntheta => nu_i, nzeta => nv_i
583  USE metrics, ONLY: guu, guv, gvv
584 
585  IMPLICIT NONE
586 
587 ! Declare Arguments
588  REAL (dp), DIMENSION(mnmax,ns), INTENT(out) :: lmns
589  REAL (dp), DIMENSION(mnmax,ns), INTENT(out) :: lmnc
590  REAL (dp), DIMENSION(nuv,ns), INTENT(in) :: jacobh
591  REAL (dp), DIMENSION(mnmax), INTENT(in) :: orthonorm
592  REAL (dp), DIMENSION(ns), INTENT(in) :: phiph
593  REAL (dp), DIMENSION(ns), INTENT(in) :: chiph
594  INTEGER, INTENT(in) :: nsmin
595  INTEGER, INTENT(in) :: nsmax
596 
597 ! Local Variables
598  INTEGER :: js
599  INTEGER :: m
600  INTEGER :: n
601  INTEGER :: mp
602  INTEGER :: np
603  INTEGER :: mn
604  INTEGER :: mnp
605  INTEGER :: lk
606  INTEGER :: lu
607  INTEGER :: lv
608  INTEGER :: info
609  INTEGER :: isign
610  INTEGER :: n2
611  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: amat
612  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: brhs
613  REAL (dp), DIMENSION(:), ALLOCATABLE :: atemp
614  REAL (dp), DIMENSION(:), ALLOCATABLE :: btemp
615  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: cosmni(:,:)
616  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: cosmnp(:,:)
617  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: sinmni(:,:)
618  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: sinmnp(:,:)
619  REAL (dp) :: ton
620  REAL (dp) :: toff
621 
622 ! Start of executable code
623  CALL second0(ton)
624 
625  lmns(:,1) = 0
626  n2 = 1
627  IF (lasym) THEN
628  lmnc(:,1) = 0
629  n2 = 2
630  END IF
631 
632  ALLOCATE(atemp(nuv), btemp(nuv), cosmni(nuv,mnmax), cosmnp(nuv,mnmax), &
633  brhs(mnmax,n2), amat(mnmax*n2, mnmax*n2), stat=info)
634  CALL assert(info.EQ.0,'Allocation error in Recompute_Lambda')
635 
636  IF (lasym) THEN
637  ALLOCATE (sinmni(nuv,mnmax), sinmnp(nuv,mnmax), stat=info)
638  CALL assert(info.EQ.0,'Allocation error in Recompute_Lambda')
639  END IF
640 
641 ! Compute matrix elements surface by surface.
642  mn = 0
643  DO n = -ntor, ntor
644  isign = 1
645  IF (n .lt. 0) THEN
646  isign = -1
647  END IF
648 
649  np = abs(n)
650  DO m = 0, mpol
651  mn = mn + 1
652  lk = 0
653 
654  DO lv = 1, nzeta
655  DO lu = 1, ntheta
656  lk = lk + 1
657 
658  cosmni(lk,mn) = fourier_context%cosmui(lu,m) * &
659  fourier_context%cosnv(lv,np) &
660  - isign*fourier_context%sinmui(lu,m) * &
661  fourier_context%sinnv(lv,np)
662  cosmnp(lk,mn) = fourier_context%cosmu(lu,m) * &
663  fourier_context%cosnv(lv,np) &
664  - isign*fourier_context%sinmu(lu,m) * &
665  fourier_context%sinnv(lv,np)
666 
667  IF (lasym) THEN
668  sinmni(lk,mn) = fourier_context%sinmui(lu,m) * &
669  fourier_context%cosnv(lv,np) &
670  + isign*fourier_context%cosmui(lu,m) * &
671  fourier_context%sinnv(lv,np)
672  sinmnp(lk,mn) = fourier_context%sinmu(lu,m) * &
673  fourier_context%cosnv(lv,np) &
674  + isign*fourier_context%cosmu(lu,m) * &
675  fourier_context%sinnv(lv,np)
676  END IF
677  END DO
678  END DO
679 
680  IF (m .eq. 0 .and. n .lt. 0) THEN
681  cosmni(:,mn) = 0
682  cosmnp(:,mn) = 0
683  IF (lasym) THEN
684  sinmni(:,mn) = 0
685  sinmnp(:,mn) = 0
686  END IF
687  END IF
688  END DO
689  END DO
690 
691 ! Compute lambda from sqrt(g)*J^s == mB_v - nB_u = 0
692  DO js = max(nsmin,2), nsmax
693  mn = 0
694 
695  DO n = -ntor, ntor
696  DO m = 0, mpol
697 
698 ! Right hand side terms sinmi and -cosmi parts.
699  mn = mn + 1
700  btemp = (-m* (chiph(js)*guv(:,js) + phiph(js)*gvv(:,js)) + &
701  n*nfp*(chiph(js)*guu(:,js) + phiph(js)*guv(:,js))) &
702  / jacobh(:,js)
703 
704  brhs(mn,1) = sum(cosmni(:,mn)*btemp)
705  IF (lasym) THEN
706  brhs(mn,2) = -sum(sinmni(:,mn)*btemp)
707  END IF
708 
709 ! Matrix elements. Coefficients of lambda(ms, js)
710  mnp = 0
711 
712  DO np = -ntor, ntor
713  DO mp = 0, mpol
714  mnp = mnp + 1
715  atemp = (m*mp*gvv(:,js) + n*nfp*np*nfp*guu(:,js) - &
716  nfp*(m*np + mp*n)*guv(:,js))/jacobh(:,js)
717 
718  amat(mn,mnp) = sum(cosmni(:,mn)*cosmnp(:,mnp)*atemp)
719  IF (mnp .eq. mn .and. amat(mn,mnp) .eq. zero) THEN
720  amat(mn,mnp) = signjac
721  END IF
722 
723  IF (lasym) THEN
724  amat(mn,mnp+mnmax) = &
725  -sum(cosmni(:,mn)*sinmnp(:,mnp)*atemp)
726  amat(mn+mnmax,mnp) = &
727  -sum(sinmni(:,mn)*cosmnp(:,mnp)*atemp)
728  amat(mn+mnmax,mnp+mnmax) = &
729  sum(sinmni(:,mn)*sinmnp(:,mnp)*atemp)
730  IF (mnp .eq. mn .and. &
731  amat(mn+mnmax,mnp+mnmax).eq. zero) THEN
732  amat(mn+mnmax,mnp+mnmax) = signjac
733  END IF
734  END IF
735  END DO
736  END DO
737  END DO
738  END DO
739 
740  CALL solver(amat, brhs, mnmax*n2, 1, info)
741  CALL assert(info.EQ.0,'INFO != 0 IN RECOMPUTE_LAMBDA')
742 
743 ! Divide by phiph -> old style
744  lmns(:,js) = brhs(:,1)/(orthonorm*phiph(js))
745 !! lmns(:,js) = brhs(:,1)/orthonorm/phiph(js)
746  IF (lasym) THEN
747  lmnc(:,js) = brhs(:,2)/(orthonorm*phiph(js))
748 !! lmnc(:,js) = brhs(:,2)/orthonorm/phiph(js)
749  END IF
750  END DO
751 
752  CALL assert(mn .eq. mnmax .and. mnp .eq. mnmax, &
753  'mn or mnp != mnmax in RECOMPUTE_LAMBDA')
754 
755  DEALLOCATE(atemp, btemp, cosmni, cosmnp, brhs, amat)
756  IF (lasym) THEN
757  DEALLOCATE(sinmni, sinmnp)
758  END IF
759 
760  CALL second0(toff)
761  IF (iam .EQ. 0 .and. lverbose) THEN
762  WRITE (*,1000) toff - ton
763  END IF
764 
765 1000 FORMAT(' RECOMPUTE_LAMBDA - TIME: ',f10.2,' S')
766  END SUBROUTINE
767 
768 !-------------------------------------------------------------------------------
773 !-------------------------------------------------------------------------------
774  SUBROUTINE alloc_quantities
775 
776  IMPLICIT NONE
777 
778 ! local variables
779  INTEGER :: istat
780 
781 ! Start of executable code
782 ! Half mesh quantities (harmonics), B^s (sine), B^u (cosine), B^v (cosine)
783  ALLOCATE(jbsupsmnsh(0:mpol,-ntor:ntor,ns), &
784  jbsupumnch(0:mpol,-ntor:ntor,ns), &
785  jbsupvmnch(0:mpol,-ntor:ntor,ns), &
786  jpmnch(0:mpol,-ntor:ntor,ns), stat=istat)
787  CALL assert_eq(0, istat, 'Allocation #1 failed in alloc_quantities')
788  jbsupsmnsh = zero
789  jbsupumnch = zero
790  jbsupvmnch = zero
791  jpmnch = zero
792 
793  ALLOCATE(pwr_spec_s(0:mpol,-ntor:ntor,ns,4), &
794  pwr_spec_a(0:mpol,-ntor:ntor,ns,4), stat=istat)
795  CALL assert_eq(0, istat, 'Allocation #2 failed in alloc_quantities')
796 
797  IF (lasym) THEN
798 ! Half mesh quantities (harmonics), B^s (sine), B^u (cosine), B^v (cosine)
799  ALLOCATE(jbsupsmnch(0:mpol,-ntor:ntor,ns), &
800  jbsupumnsh(0:mpol,-ntor:ntor,ns), &
801  jbsupvmnsh(0:mpol,-ntor:ntor,ns), &
802  jpmnsh(0:mpol,-ntor:ntor,ns), stat=istat)
803  CALL assert_eq(0, istat, 'Allocation #3 failed in alloc_quantities')
804  jbsupsmnch = zero
805  jbsupumnsh = zero
806  jbsupvmnsh = zero
807  jpmnsh = zero
808  END IF
809 
810  ALLOCATE(bsq(ntheta,nzeta,ns), jacobh(ntheta,nzeta,ns), &
811  jacobf(ntheta,nzeta,ns), wint(ntheta,nzeta,ns), stat=istat)
812  bsq = 0.0
813  CALL assert_eq(0, istat, 'Allocation #4 failed in alloc_quantities')
814 
815  END SUBROUTINE
816 
817 !-------------------------------------------------------------------------------
822 !-------------------------------------------------------------------------------
823  SUBROUTINE dealloc_quantities
824 
825  IMPLICIT NONE
826 
827 ! Start of executable code
828 ! Quantities allocated in alloc_quantities
829  DEALLOCATE(jpmnch, jbsupsmnsh, jbsupumnch, jbsupvmnch)
830  DEALLOCATE(pwr_spec_s, pwr_spec_a)
831 
832 ! Other quantities
833  DEALLOCATE(djpmnch, djbsupsmnsh, djbsupumnch, djbsupvmnch, &
834  ksupsmnsf, ksupumncf, ksupvmncf)
835  DEALLOCATE(jacobh, jacobf, jvsupsijf, jvsupuijf, jvsupvijf, &
836  ksubsijf, ksubuijf, ksubvijf, pijh0, pijf0, &
837  pijf0_ds, bsupuijf0, vp_f, bsq, wint, &
838  bsupvijf0, bsupsijf0, pijh0_du, pijh0_dv)
839  IF (ALLOCATED(bsupsijf0)) THEN
840  DEALLOCATE(bsupsijf0, bsupuijf0, bsupvijf0)
841  END IF
842  IF (ALLOCATED(bsupsijh0)) THEN
843  DEALLOCATE(bsupsijh0, bsupuijh0, bsupvijh0)
844  END IF
845  IF (ALLOCATED(ksupsijf0)) THEN
846  DEALLOCATE(ksupsijf0, ksupuijf0, ksupvijf0)
847  END IF
848  IF (ALLOCATED(bsubsijf)) THEN
849  DEALLOCATE(bsubsijf, bsubuijf, bsubvijf)
850  END IF
851 
852 ! Asymmetric quantities.
853  IF (lasym) THEN
854 ! Quantities allocated in alloc_quantities
855  DEALLOCATE(jpmnsh, jbsupsmnch, jbsupumnsh, jbsupvmnsh)
856 ! Other quantities
857  DEALLOCATE(djpmnsh, djbsupsmnch, djbsupumnsh, djbsupvmnsh, &
858  ksupsmncf, ksupumnsf, ksupvmnsf)
859  END IF
860 
861  END SUBROUTINE
862 
863 !-------------------------------------------------------------------------------
865 !-------------------------------------------------------------------------------
866  SUBROUTINE toupper_forces
867  USE island_params, ONLY: ntheta=>nu_i, nzeta=>nv_i, &
868  mpol=>mpol_i, ntor=>ntor_i, &
869  fourier_context
870  USE fourier, ONLY: f_none, f_cos, f_sin, f_sum
871  USE metrics, ONLY: toupper
872  USE utilities, ONLY: set_bndy_full_origin
873 
874 ! Local Variables
875  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: fsubsijf
876  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: fsubuijf
877  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: fsubvijf
878  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: fsupsijf
879  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: fsupuijf
880  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: fsupvijf
881  INTEGER :: nsmin
882  INTEGER :: nsmax
883  INTEGER :: istat
884 
885 ! Start of executable code.
886  nsmin = max(1, startglobrow)
887  nsmax = min(ns, endglobrow)
888 
889  ALLOCATE(fsubsijf(ntheta,nzeta,nsmin:nsmax), &
890  fsubuijf(ntheta,nzeta,nsmin:nsmax), &
891  fsubvijf(ntheta,nzeta,nsmin:nsmax), &
892  fsupsijf(ntheta,nzeta,nsmin:nsmax), &
893  fsupuijf(ntheta,nzeta,nsmin:nsmax), &
894  fsupvijf(ntheta,nzeta,nsmin:nsmax), stat=istat)
895  CALL assert(istat.eq.0,' Allocation error in toupper_forces')
896 
897  CALL fourier_context%toijsp(fsubsmncf(:,:,nsmin:nsmax), fsubsijf, &
898  f_none, f_cos)
899  CALL fourier_context%toijsp(fsubumnsf(:,:,nsmin:nsmax), fsubuijf, &
900  f_none, f_sin)
901  CALL fourier_context%toijsp(fsubvmnsf(:,:,nsmin:nsmax), fsubvijf, &
902  f_none, f_sin)
903 
904  IF (lasym) THEN
905  CALL fourier_context%toijsp(fsubsmnsf(:,:,nsmin:nsmax), fsubsijf, &
906  f_sum, f_sin)
907  CALL fourier_context%toijsp(fsubumncf(:,:,nsmin:nsmax), fsubuijf, &
908  f_sum, f_cos)
909  CALL fourier_context%toijsp(fsubvmncf(:,:,nsmin:nsmax), fsubvijf, &
910  f_sum, f_cos)
911  END IF
912 
913  IF (nsmin .eq. 1) THEN
914  fsubuijf(:,:,1) = 0
915  END IF
916 
917 ! Compute contravariant forces needed to compute <||F||^2>.
918  CALL toupper(fsubsijf, fsubuijf, fsubvijf, &
919  fsupsijf, fsupuijf, fsupvijf, nsmin, nsmax)
920  DEALLOCATE(fsubsijf, fsubuijf, fsubvijf)
921 
922  CALL fourier_context%tomnsp(fsupsijf, fsupsmncf(:,:,nsmin:nsmax), f_cos)
923  CALL fourier_context%tomnsp(fsupuijf, fsupumnsf(:,:,nsmin:nsmax), f_sin)
924  CALL fourier_context%tomnsp(fsupvijf, fsupvmnsf(:,:,nsmin:nsmax), f_sin)
925 
926  IF (nsmin .eq. 1) THEN
927  CALL set_bndy_full_origin(fsupsmncf, fsupumnsf, fsupvmnsf)
928  END IF
929 
930  IF (lasym) THEN
931  CALL fourier_context%tomnsp(fsupsijf, fsupsmnsf(:,:,nsmin:nsmax), &
932  f_sin)
933  CALL fourier_context%tomnsp(fsupuijf, fsupumncf(:,:,nsmin:nsmax), &
934  f_cos)
935  CALL fourier_context%tomnsp(fsupvijf, fsupvmncf(:,:,nsmin:nsmax), &
936  f_cos)
937 
938  IF (nsmin .eq. 1) THEN
939  CALL set_bndy_full_origin(fsupsmnsf, fsupumncf, fsupvmncf)
940  END IF
941  END IF
942 
943  DEALLOCATE(fsupsijf, fsupuijf, fsupvijf, stat=istat)
944 
945  END SUBROUTINE
946 
947 !-------------------------------------------------------------------------------
949 !-------------------------------------------------------------------------------
950  SUBROUTINE gatherfields
951 
952  IMPLICIT NONE
953 
954 ! Start of executable code
955  IF (parsolver) THEN
956  CALL gather_fields(jbsupsmnsh, jbsupumnch, jbsupvmnch, jpmnch)
957  IF (lasym) THEN
958  CALL gather_fields(jbsupsmnch, jbsupumnsh, jbsupvmnsh, jpmnsh)
959  END IF
960  END IF
961 
962  END SUBROUTINE
963 
964 !-------------------------------------------------------------------------------
974 !-------------------------------------------------------------------------------
975  SUBROUTINE gather_fields(jbsupsmnh, jbsupumnh, jbsupvmnh, jpmnh)
976  USE hessian, ONLY: gather_array
977 
978  IMPLICIT NONE
979 
980 ! Start of executable code
981  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: jbsupsmnh
982  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: jbsupumnh
983  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: jbsupvmnh
984  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: jpmnh
985 
986 ! Start of executable code.
987  CALL assert_eq(ns, SIZE(jbsupsmnh,3), &
988  'Radial array size needs to be ns in GATHER_FIELDS')
989 
990  CALL gather_array(jbsupsmnh)
991  CALL gather_array(jbsupumnh)
992  CALL gather_array(jbsupvmnh)
993  CALL gather_array(jpmnh)
994 
995  END SUBROUTINE
996 
997 !-------------------------------------------------------------------------------
1007 !-------------------------------------------------------------------------------
1008  SUBROUTINE surfaverage(average, q3d, nsmin, nsmax)
1009  USE stel_kinds, ONLY: dp
1010  USE island_params, ONLY: ns=>ns_i, ntheta=>nu_i, nzeta=>nv_i
1011 
1012  IMPLICIT NONE
1013 
1014 ! Declare Arguments.
1015  REAL(dp), DIMENSION(nsmin:nsmax), INTENT(out) :: average
1016  REAL(dp), DIMENSION(ntheta,nzeta,nsmin:nsmax), INTENT(in) :: q3d
1017  INTEGER, INTENT(in) :: nsmin
1018  INTEGER, INTENT(in) :: nsmax
1019 
1020 ! Local variables.
1021  INTEGER :: js
1022 
1023 ! Start of executable code.
1024  DO js = nsmin, nsmax
1025  average(js) = sum(q3d(:,:,js)*wint(:,:,js))
1026  END DO
1027 
1028  END SUBROUTINE
1029 
1030  END MODULE
siesta_namelist::lrestart
logical lrestart
Restart SIESTA from pervious run.
Definition: siesta_namelist.f90:128
shared_data::lverbose
logical lverbose
Use verbose screen output.
Definition: shared_data.f90:234
siesta_namelist::l_vessel
logical l_vessel
If extended grid is to be used using an available vessel file.
Definition: siesta_namelist.f90:140
hessian::gather_array
Definition: hessian.f90:31
quantities
This file contains subroutines for allocating and initializing curvilinear magnetic covariant and pre...
Definition: quantities.f90:11
metrics::sqrtg
real(dp), dimension(:,:), allocatable sqrtg
Jacobian on half grid.
Definition: metrics.f90:38
v3_utilities::assert_eq
Definition: v3_utilities.f:62
island_params::fourier_context
type(fourier_class), pointer fourier_context
Fourier transform object.
Definition: island_params.f90:76
fourier
Module contains subroutines for computing FFT on parallel radial intervals. Converts quantities betwe...
Definition: fourier.f90:13
fourier::m0
integer, parameter m0
m = 0 mode.
Definition: fourier.f90:58
v3_utilities::assert
Definition: v3_utilities.f:55
fourier::f_sin
integer, parameter f_sin
Sine parity.
Definition: fourier.f90:27
metrics::guu
real(dp), dimension(:,:), allocatable guu
Lower metric tensor half grid. e_u . e_u.
Definition: metrics.f90:46
metrics::gvv
real(dp), dimension(:,:), allocatable gvv
Lower metric tensor half grid. e_v . e_v.
Definition: metrics.f90:50
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
fourier::n0
integer, parameter n0
n = 0 mode.
Definition: fourier.f90:65
shared_data::lasym
logical lasym
Use non-stellarator symmetry.
Definition: shared_data.f90:230
metrics::toupper
subroutine toupper(xsubsij, xsubuij, xsubvij,
Converts to contravariant component full grid.
Definition: metrics.f90:485
blocktridiagonalsolver_s
Module contains subroutines for solver for block tri-diagonal matrices.
Definition: blocktridiagonalsolver_s.f90:7
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
shared_data::ndims
integer ndims
Number of independent variables.
Definition: shared_data.f90:58
shared_data::mblk_size
integer mblk_size
Block size. (mpol + 1)*(2*ntor + 1)*ndims.
Definition: shared_data.f90:62
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
This file contains variables and parameters used by many modules in SIESTA.
Definition: shared_data.f90:10
island_params::vp_f
real(dp), dimension(:), allocatable vp_f
Volume of a radial grid surface.
Definition: island_params.f90:74
metrics::guv
real(dp), dimension(:,:), allocatable guv
Lower metric tensor half grid. e_u . e_v.
Definition: metrics.f90:48
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