V3FIT
siesta_force.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 !
36 !*******************************************************************************
37  MODULE siesta_force
38  USE stel_kinds
39  USE stel_constants
40  USE v3_utilities, ONLY: assert
41  USE descriptor_mod, ONLY: siesta_comm, iam, nprocs
42  USE fourier
43  USE island_params, ONLY: nu_i, hs_i
44  USE timer_mod
45  USE shared_data, ONLY: l_linearize, l_getwmhd, lverbose, l_pedge
46  USE nscalingtools, ONLY: startglobrow, endglobrow, mpi_err
47  USE utilities, ONLY: gradientfull, to_full_mesh
48  USE quantities
49  USE mpi_inc
50 
51  IMPLICIT NONE
52 
53  CONTAINS
54 
55 !*******************************************************************************
56 ! UTILITY SUBROUTINES
57 !*******************************************************************************
58 !-------------------------------------------------------------------------------
68 !-------------------------------------------------------------------------------
69  SUBROUTINE update_force
70  USE shared_data, ONLY: gc, col_scale
71  USE siesta_currents, ONLY: cv_currents
72  USE hessian, ONLY: apply_colscale
73  USE bhtobf
74 
75  IMPLICIT NONE
76 
77 ! local variables
78  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: bsupsijh
79  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: bsupuijh
80  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: bsupvijh
81  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: pijh
82  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: pijf1
83  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: KxBsij
84  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: KxBuij
85  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: KxBvij
86  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: pardamp
87  INTEGER :: istat
88  INTEGER :: n1
89  INTEGER :: n2
90  REAL (dp) :: ton
91  REAL (dp) :: toff
92  REAL (dp) :: skston
93  REAL (dp) :: skstoff
94  INTEGER :: nsmin
95  INTEGER :: nsmax
96 
97 ! Start of executable code
98  CALL second0(skston)
99  nsmin = max(1, startglobrow)
100  nsmax = min(ns, endglobrow + 1)
101 
102  ALLOCATE(bsupsijh(ntheta,nzeta,nsmin:nsmax), &
103  bsupuijh(ntheta,nzeta,nsmin:nsmax), &
104  bsupvijh(ntheta,nzeta,nsmin:nsmax), &
105  pijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
106  CALL assert_eq(0, istat, 'Allocation1 failed in UPDATE_FORCE')
107 
108 ! Compute real-space contravariant components of B and p. Fields are updated
109 ! in update state, not here. Unperturbed fields are unchanged until
110 ! update_state call. Perturbed fields were calculated in update_bfield and
111 ! update_press.
112  CALL incfields(jbsupsmnsh, jbsupumnch, jbsupvmnch, jpmnch, &
113  djbsupsmnsh, djbsupumnch, djbsupvmnch, djpmnch, &
114  bsupsijh, bsupuijh, bsupvijh, pijh, f_sin)
115  IF (lasym) THEN
116  CALL incfields(jbsupsmnch, jbsupumnsh, jbsupvmnsh, jpmnsh, &
117  djbsupsmnch, djbsupumnsh, djbsupvmnsh, djpmnsh, &
118  bsupsijh, bsupuijh, bsupvijh, pijh, f_cos)
119  END IF
120 
121 ! Update thermal energy (pressure based). pijh contains the jacobian term at
122 ! this point.
123  wpres: IF (l_getwmhd) THEN
124  CALL assert(.not.inhessian, &
125  'l_getwmhd must be set to FALSE in Hessian')
126  n2 = min(ns, endglobrow)
127  n1 = max(2, nsmin)
128  wp = signjac*(twopi*twopi*hs_i)*sum(pijh(:,:,n1:n2)*wint(:,:,n1:n2))
129  n1 = nsmin
130  n2 = nsmax
131 #if defined(MPI_OPT)
132  IF (parsolver) THEN
133 ! FIXME: All reduce is not deterministic. This causes a divergent run sequence.
134  CALL mpi_allreduce(mpi_in_place, wp, 1, mpi_real8, mpi_sum, &
135  siesta_comm, mpi_err)
136  END IF
137 #endif
138  END IF wpres
139 
140  ALLOCATE(pijf1(ntheta,nzeta,nsmin:nsmax), stat=istat)
141  CALL assert_eq(0, istat, 'Allocation2 failed in UPDATE_FORCE')
142 
143 ! Update full mesh components of bsup* and pressure in real-space by averaging.
144 ! This removes the jacobian term.
145  CALL bhalftobfull(bsupsijh, bsupuijh, bsupvijh, &
146  bsupsijf, bsupuijf, bsupvijf, pijh, pijf1)
147 
148 ! Calculate contravariant components of J = curl B (Ksup = gsqrt*Jsup)
149 ! and update the magnetic energy if nonlinear force is being computed.
150  CALL cv_currents(bsupsijh, bsupuijh, bsupvijh, &
151  ksupsijf, ksupuijf, ksupvijf, &
152  .not.l_linearize .OR. l_getwmhd, .false.)
153 
154  DEALLOCATE(bsupsijh, bsupuijh, bsupvijh, stat=istat)
155 
156  nsmax = min(ns, endglobrow)
157 
158  ALLOCATE(kxbsij(ntheta,nzeta,nsmin:nsmax), &
159  kxbuij(ntheta,nzeta,nsmin:nsmax), &
160  kxbvij(ntheta,nzeta,nsmin:nsmax), stat=istat)
161  CALL assert_eq(0, istat, 'Allocation3 failed in UPDATE_FORCE')
162 
163 ! KxBsij used as a scratch array.
164  CALL initpardamping(kxbsij, pardamp)
165 
166 ! Compute curvilinear Lorentz force components
167  IF (l_linearize) THEN
168 ! dJ X B0
169  CALL lorentz(bsupsijf0, bsupuijf0, bsupvijf0, &
170  ksupsijf, ksupuijf, ksupvijf, &
171  kxbsij, kxbuij, kxbvij, f_none)
172 ! J0 X dB
173  CALL lorentz(bsupsijf, bsupuijf, bsupvijf, &
174  ksupsijf0, ksupuijf0, ksupvijf0, &
175  kxbsij, kxbuij, kxbvij, f_sum)
176  ELSE
177 ! (J0+dJ) X (B0+dB)
178  CALL lorentz(bsupsijf, bsupuijf, bsupvijf, &
179  ksupsijf, ksupuijf, ksupvijf, &
180  kxbsij, kxbuij, kxbvij, f_none)
181  END IF
182 
183  IF (nsmax .eq. ns .and. .not.l_vessel) THEN
184  CALL assert(all(bsupsijf(:,:,ns) .eq. zero), &
185  'bsupsijf(ns) != 0 in UPDATE_FORCE')
186  END IF
187 
188  CALL getmhdforce(fsubsmncf, fsubumnsf, fsubvmnsf, &
189  pijh, kxbsij, kxbuij, kxbvij, pardamp, f_cos)
190  IF (lasym) THEN
191  CALL getmhdforce(fsubsmnsf, fsubumncf, fsubvmncf, &
192  pijh, kxbsij, kxbuij, kxbvij, pardamp, f_sin)
193  END IF
194 
195  CALL apply_colscale(gc, col_scale, nsmin, nsmax)
196 
197  DEALLOCATE(pijf1, kxbsij, kxbuij, kxbvij, pijh, stat=istat)
198  IF (ALLOCATED(pardamp)) THEN
199  DEALLOCATE(pardamp)
200  END IF
201 
202  CALL second0(skstoff)
203  time_update_force = time_update_force + (skstoff - skston)
204 
205  END SUBROUTINE
206 
207 !-------------------------------------------------------------------------------
229 !-------------------------------------------------------------------------------
230  SUBROUTINE incfields(jbsupsmnh, jbsupumnh, jbsupvmnh, jpmnh, &
231  djbsupsmnh, djbsupumnh, djbsupvmnh, djpmnh, &
232  jbsupsijh, jbsupuijh, jbsupvijh, jpijh, &
233  parity)
234  USE island_params, ONLY: fourier_context
235 
236  IMPLICIT NONE
237 
238 ! Declare Arguments
239  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupsmnh
240  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupumnh
241  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupvmnh
242  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jpmnh
243  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: djbsupsmnh
244  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: djbsupumnh
245  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: djbsupvmnh
246  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: djpmnh
247  REAL (dp), DIMENSION(:,:,:), INTENT(out) :: jbsupsijh
248  REAL (dp), DIMENSION(:,:,:), INTENT(out) :: jbsupuijh
249  REAL (dp), DIMENSION(:,:,:), INTENT(out) :: jbsupvijh
250  REAL (dp), DIMENSION(:,:,:), INTENT(out) :: jpijh
251  INTEGER, INTENT(in) :: parity
252 
253 ! local variables
254  INTEGER :: fours
255  INTEGER :: fouruv
256  INTEGER :: fcomb
257  INTEGER :: istat
258  INTEGER :: nsmin
259  INTEGER :: nsmax
260  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: jbsupsmnh_i
261  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: jbsupumnh_i
262  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: jbsupvmnh_i
263  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: jpmnh_i
264 
265 ! Start of executable code
266  nsmin = max(1, startglobrow)
267  nsmax = min(ns, endglobrow + 1)
268 
269  IF (parity .EQ. f_sin) THEN
270  fcomb = f_none
271  fours = f_sin
272  fouruv = f_cos
273  ELSE
274  fcomb = f_sum
275  fours = f_cos
276  fouruv = f_sin
277  END IF
278 
279 ! Allocate space for total values of bsup*mnh_i and pmnh_i.
280  ALLOCATE(jbsupsmnh_i(0:mpol,-ntor:ntor,nsmin:nsmax), &
281  jbsupumnh_i(0:mpol,-ntor:ntor,nsmin:nsmax), &
282  jbsupvmnh_i(0:mpol,-ntor:ntor,nsmin:nsmax), &
283  jpmnh_i(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
284  CALL assert_eq(0, istat, 'Allocation failed in IncFields')
285 
286 ! For nonlinear force, add perturbation to unperturbed value. For linear force,
287 ! compute JUST perturbed fields here.
288  jbsupsmnh_i = djbsupsmnh
289  jbsupumnh_i = djbsupumnh
290  jbsupvmnh_i = djbsupvmnh
291  jpmnh_i = djpmnh
292  IF (.not.l_linearize) THEN
293  jbsupsmnh_i = jbsupsmnh_i + jbsupsmnh(:,:,nsmin:nsmax)
294  jbsupumnh_i = jbsupumnh_i + jbsupumnh(:,:,nsmin:nsmax)
295  jbsupvmnh_i = jbsupvmnh_i + jbsupvmnh(:,:,nsmin:nsmax)
296  jpmnh_i = jpmnh_i + jpmnh(:,:,nsmin:nsmax)
297  END IF
298 
299 ! js=1 origin values. Note: This assumed that the magnetic axis and the grid
300 ! axis align. Maybe the cause of poor convergence.
301  IF (nsmin .EQ. 1) THEN
302  jbsupsmnh_i(:,:,1) = 0
303  jbsupsmnh_i(m1,:,1) = jbsupsmnh_i(m1,:,2)
304  jbsupumnh_i(:,:,1) = 0
305  jbsupumnh_i(m0:m1,:,1) = jbsupumnh_i(m0:m1,:,2)
306  jbsupvmnh_i(:,:,1) = 0
307  jbsupvmnh_i(m0,:,1) = jbsupvmnh_i(m0,:,2)
308  jpmnh_i(:,:,1) = 0
309  jpmnh_i(m0,:,1) = jpmnh_i(m0,:,2)
310  END IF
311 
312 ! Calculate real-space jbsup*ihh and jpijh half mesh.
313  CALL fourier_context%toijsp(jbsupsmnh_i, jbsupsijh, fcomb, fours)
314  CALL fourier_context%toijsp(jbsupumnh_i, jbsupuijh, fcomb, fouruv)
315  CALL fourier_context%toijsp(jbsupvmnh_i, jbsupvijh, fcomb, fouruv)
316  CALL fourier_context%toijsp(jpmnh_i, jpijh, fcomb, fouruv)
317 
318  DEALLOCATE(jbsupsmnh_i, jbsupumnh_i, jbsupvmnh_i, jpmnh_i)
319 
320  END SUBROUTINE
321 
322 !-------------------------------------------------------------------------------
346 !-------------------------------------------------------------------------------
347  SUBROUTINE lorentz(bsupsijf, bsupuijf, bsupvijf, &
348  ksupsijf, ksupuijf, ksupvijf, &
349  KxBsij, KxBuij, KxBvij, fcomb)
350 
351  IMPLICIT NONE
352 
353 ! Declare Arguments
354  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: bsupsijf
355  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: bsupuijf
356  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: bsupvijf
357  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: ksupsijf
358  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: ksupuijf
359  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: ksupvijf
360  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: KxBsij
361  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: KxBuij
362  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: KxBvij
363  INTEGER, INTENT(in) :: fcomb
364 
365 ! local variables
366  INTEGER :: nsmin
367  INTEGER :: nsmax
368 
369 ! Start of executable code
370  nsmin = max(1, startglobrow)
371  nsmax = min(ns, endglobrow)
372 
373 ! Initialize to avoid NaN errors.
374  IF (fcomb .eq. f_none) THEN
375  kxbsij = 0
376  kxbuij = 0
377  kxbvij = 0
378  END IF
379 
380 
381  kxbsij = kxbsij &
382  + ksupuijf(:,:,nsmin:nsmax)*bsupvijf(:,:,nsmin:nsmax) &
383  - ksupvijf(:,:,nsmin:nsmax)*bsupuijf(:,:,nsmin:nsmax)
384  kxbuij = kxbuij &
385  + ksupvijf(:,:,nsmin:nsmax)*bsupsijf(:,:,nsmin:nsmax) &
386  - ksupsijf(:,:,nsmin:nsmax)*bsupvijf(:,:,nsmin:nsmax)
387  kxbvij = kxbvij &
388  + ksupsijf(:,:,nsmin:nsmax)*bsupuijf(:,:,nsmin:nsmax) &
389  - ksupuijf(:,:,nsmin:nsmax)*bsupsijf(:,:,nsmin:nsmax)
390 
391  END SUBROUTINE
392 
393 !-------------------------------------------------------------------------------
427 !-------------------------------------------------------------------------------
428  SUBROUTINE getmhdforce(fsubsmnf, fsubumnf, fsubvmnf, pijh, &
429  KxBsij, KxBuij, KxBvij, pardamp, parity)
430  USE island_params, ONLY: ohs=>ohs_i, fourier_context
431  USE utilities, ONLY: set_bndy_half_to_full, set_bndy_half_to_full_ds, &
432  set_bndy_full_origin
433 
434  IMPLICIT NONE
435 
436 ! Declare Arguments
437  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: fsubsmnf
438  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: fsubumnf
439  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: fsubvmnf
440  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: pijh
441  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: KxBsij
442  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: KxBuij
443  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: KxBvij
444  REAL (dp), DIMENSION(:,:), ALLOCATABLE, INTENT(inout) :: pardamp
445  INTEGER, INTENT(in) :: parity
446 
447 ! local variables
448  INTEGER :: fours
449  INTEGER :: fouruv
450  INTEGER :: sparity
451  INTEGER :: m
452  INTEGER :: n
453  INTEGER :: moff
454  INTEGER :: noff
455  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: pmnh
456  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: pmnf
457  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: pmnf_ds
458  INTEGER :: nsmin
459  INTEGER :: nsmax
460  INTEGER :: istat
461 
462 ! Start of executable code
463  nsmax = min(ns, endglobrow + 1)
464  nsmin = max(1, startglobrow)
465 
466  ALLOCATE(pmnf_ds(0:mpol,-ntor:ntor,nsmin:nsmax), &
467  pmnf(0:mpol,-ntor:ntor,nsmin:nsmax), &
468  pmnh(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
469  CALL assert_eq(0, istat, 'Allocation failed before GetMHDForce')
470 
471  nsmax = min(ns, endglobrow)
472 
473  IF (parity .EQ. f_cos) THEN
474  fours = f_cos
475  fouruv = f_sin
476  sparity = -1
477  ELSE
478  fours = f_sin
479  fouruv = f_cos
480  sparity = 1
481  END IF
482 
483 ! Harmonics of covariant Lorentz force components on full grid.
484  CALL fourier_context%tomnsp(kxbsij, fsubsmnf(:,:,nsmin:nsmax), fours)
485  CALL fourier_context%tomnsp(kxbuij, fsubumnf(:,:,nsmin:nsmax), fouruv)
486  CALL fourier_context%tomnsp(kxbvij, fsubvmnf(:,:,nsmin:nsmax), fouruv)
487 
488  IF (nsmin .eq. 1) THEN
489  CALL set_bndy_full_origin(fsubsmnf, fsubumnf, fsubvmnf)
490  END IF
491 
492 ! Harmonics of true pressure. The jacobian factor has been divided out of pijh.
493  CALL fourier_context%tomnsp(pijh, pmnh, fours)
494  CALL to_full_mesh(pmnh, pmnf)
495  CALL set_bndy_half_to_full(pmnh, pmnf, nsmin, fours)
496 
497 ! Radial derivative of pressure.
498  CALL gradientfull(pmnf_ds, pmnh)
499  CALL set_bndy_half_to_full_ds(pmnh, pmnf_ds, nsmin, fours)
500 
501  IF (nsmax.EQ.ns .AND. l_pedge) THEN
502  pmnf(:,:,ns) = 0
503  pmnf(m0,0,ns) = pmnh(m0,0,ns)
504  END IF
505 
506 ! Add pressure gradient to the lorentz force.
507  fsubsmnf(:,:,nsmin:nsmax) = fsubsmnf(:,:,nsmin:nsmax) - pmnf_ds
508  DO m = 0, mpol
509  moff = m + lbound(fsubumnf, 1)
510  fsubumnf(moff,:,nsmin:nsmax) = fsubumnf(moff,:,nsmin:nsmax) &
511  - m*sparity*pmnf(m,:,nsmin:nsmax)
512  END DO
513  DO n = -ntor, ntor
514  noff = n + ntor + lbound(fsubvmnf, 2)
515  fsubvmnf(:,noff,nsmin:nsmax) = fsubvmnf(:,noff,nsmin:nsmax) &
516  - n*nfp*sparity*pmnf(:,n,nsmin:nsmax)
517  END DO
518 
519  DEALLOCATE(pmnf, pmnf_ds, pmnh)
520 
521  IF (nsmin .eq. 1) THEN
522  CALL set_bndy_full_origin(fsubsmnf, fsubumnf, fsubvmnf)
523  END IF
524 
525  CALL get_force_harmonics(pardamp, fsubsmnf, fsubumnf, fsubvmnf, parity)
526 
527  END SUBROUTINE
528 
529 !-------------------------------------------------------------------------------
534 !-------------------------------------------------------------------------------
535  SUBROUTINE initpardamping(parscale, pardamp)
536  USE hessian, ONLY: l_compute_hessian, mupar, mupar_norm
537 
538 ! Declare Arguments
539  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: parscale
540  REAL (dp), DIMENSION(:,:), ALLOCATABLE, INTENT(inout) :: pardamp
541 
542 ! local variables
543  INTEGER :: nsmin
544  INTEGER :: nsmax
545 
546 ! Start of executable code
547  IF (mupar .eq. zero .or. .not.l_compute_hessian) THEN
548  RETURN
549  ELSE
550  nsmin = max(1, startglobrow)
551  nsmax = min(ns, endglobrow)
552 
553  CALL assert_eq(nsmin, lbound(parscale, 3), 'PARSCALE LBOUND WRONG!')
554  CALL assert_eq(nsmax, ubound(parscale, 3), 'PARSCALE UBOUND WRONG!')
555  ALLOCATE(pardamp(nsmin:nsmax,3))
556 
557 ! Update mupar_norm in cv_currents and put in 1/jacobh factor. Simpler and
558 ! faster approximation, keep only diagonal displacement terms from v dot D.
559  mupar = abs(mupar)
560 
561  parscale = bsubsijf(:,:,nsmin:nsmax)*bsubsijf(:,:,nsmin:nsmax) &
562  / jacobf(:,:,nsmin:nsmax)
563  CALL surfaverage(pardamp(nsmin:nsmax,1), parscale, nsmin, nsmax)
564  pardamp(:,1) = pardamp(:,1)*mupar*mupar_norm(nsmin:nsmax)
565 
566  parscale = bsubuijf(:,:,nsmin:nsmax)*bsubuijf(:,:,nsmin:nsmax) &
567  / jacobf(:,:,nsmin:nsmax)
568  CALL surfaverage(pardamp(nsmin:nsmax,2), parscale, nsmin, nsmax)
569  pardamp(:,2) = pardamp(:,2)*mupar*mupar_norm(nsmin:nsmax)
570 
571  parscale = bsubvijf(:,:,nsmin:nsmax)*bsubvijf(:,:,nsmin:nsmax) &
572  / jacobf(:,:,nsmin:nsmax)
573  CALL surfaverage(pardamp(nsmin:nsmax,3), parscale, nsmin, nsmax)
574  pardamp(:,3) = pardamp(:,3)*mupar*mupar_norm(nsmin:nsmax)
575  END IF
576 
577  END SUBROUTINE
578 
579 !-------------------------------------------------------------------------------
593 !-------------------------------------------------------------------------------
594  SUBROUTINE get_force_harmonics(pardamp, f_smnf, f_umnf, f_vmnf, parity)
595  USE hessian, ONLY: mupar, l_compute_hessian
596  USE shared_data, ONLY: nprecon, l_printoriginforces, &
597  l_push_s, l_push_u, l_push_v, &
598  l_push_edge, col_scale
599  USE nscalingtools, ONLY: firstpartition_ge_2, lastpartition_ge_2
600 ! These were computed in siesta_displacment and are used for || damping
601  USE quantities, ONLY: gvsupsmncf => fsupsmncf, gvsupumnsf => fsupumnsf, &
602  gvsupvmnsf => fsupvmnsf, gvsupsmnsf => fsupsmnsf, &
603  gvsupumncf => fsupumncf, gvsupvmncf => fsupvmncf
604  USE utilities, ONLY: set_bndy_full_origin, set_bndy_fouier_m0
605 
606  IMPLICIT NONE
607 
608 ! Declare Arguments
609  REAL (dp), DIMENSION(:,:), ALLOCATABLE, INTENT(inout) :: pardamp
610  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: f_smnf
611  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: f_umnf
612  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: f_vmnf
613  INTEGER, INTENT(in) :: parity
614 
615 ! local variables
616  REAL (dp) :: skston
617  REAL (dp) :: skstoff
618  REAL (dp), DIMENSION(10) :: tmps
619  INTEGER :: nsmin
620  INTEGER :: nsmax
621  INTEGER :: moff
622  INTEGER :: noff
623 
624 ! local parameters
625  REAL (dp), PARAMETER :: p5 = 0.5_dp
626 
627 ! Start of executable code
628  CALL second0(skston)
629 
630  nsmin = max(1, startglobrow)
631  nsmax = min(ns,endglobrow)
632 
633  moff = 1
634  noff = 1 + ntor
635 
636  CALL assert(all(f_smnf(m0+moff,-ntor+noff:noff-1,nsmin:nsmax) .eq. zero),&
637  'f_smnf != 0, m=0, n<0')
638  CALL assert(all(f_umnf(m0+moff,-ntor+noff:noff-1,nsmin:nsmax) .eq. zero),&
639  'f_umnf != 0, m=0, n<0')
640  CALL assert(all(f_vmnf(m0+moff,-ntor+noff:noff-1,nsmin:nsmax) .eq. zero),&
641  'f_vmnf != 0, m=0, n<0')
642  IF (ALLOCATED(pardamp)) THEN
643  IF (parity .eq. f_cos) THEN
644  CALL addpardamping(pardamp, &
645  f_smnf(:,:,nsmin:nsmax), &
646  f_umnf(:,:,nsmin:nsmax), &
647  f_vmnf(:,:,nsmin:nsmax), &
648  gvsupsmncf(:,:,nsmin:nsmax), &
649  gvsupumnsf(:,:,nsmin:nsmax), &
650  gvsupvmnsf(:,:,nsmin:nsmax))
651  CALL assert(all(f_umnf(m0+moff,noff,:) .eq. zero), &
652  'f_umnf(m0,0) != 0')
653  CALL assert(all(f_vmnf(m0+moff,noff,:) .eq. zero), &
654  'f_vmnf(m0,0) != 0')
655  ELSE
656  CALL addpardamping(pardamp, &
657  f_smnf(:,:,nsmin:nsmax), &
658  f_umnf(:,:,nsmin:nsmax), &
659  f_vmnf(:,:,nsmin:nsmax), &
660  gvsupsmnsf(:,:,nsmin:nsmax), &
661  gvsupumncf(:,:,nsmin:nsmax), &
662  gvsupvmncf(:,:,nsmin:nsmax))
663  CALL assert(all(f_smnf(m0+moff,noff,:) .eq. zero), &
664  'f_smnf(m0,0) != 0')
665  END IF
666  END IF
667 
668 ! GATHER BDY FORCES AND PRINT OUT ON PROC=0
669  print_o: IF (l_printoriginforces) THEN
670  tmps = 0
671 ! Assume s=0 is on the 0th processor.
672  IF (nsmin .le. 2 .and. nsmax .ge. 2) THEN
673  tmps(1) = sqrt(sum(f_smnf(m1+moff,:,2)**2))
674  tmps(2) = sqrt(sum(f_smnf(m0+moff,:,2)**2) + &
675  sum(f_smnf(m2+moff:,:,2)**2))
676  tmps(3) = sqrt(sum(f_umnf(m1+moff,:,2)**2))
677  tmps(4) = sqrt(sum(f_umnf(m0+moff,:,2)**2) + &
678  sum(f_umnf(m2+moff:,:,2)**2))
679  tmps(5) = sqrt(sum(f_vmnf(m0+moff,:,2)**2))
680  tmps(6) = sqrt(sum(f_vmnf(m1+moff:,:,2)**2))
681  END IF
682 
683  IF (nsmin .le. ns - 1 .and. nsmax .ge. ns - 1) THEN
684  tmps(7) = sqrt(sum(f_umnf(:,:,ns-1)**2))
685  tmps(8) = sqrt(sum(f_vmnf(:,:,ns-1)**2))
686  END IF
687 
688  IF (nsmax .eq. ns) THEN
689  tmps(9) = sqrt(sum(f_umnf(:,:,ns)**2))
690  tmps(10) = sqrt(sum(f_vmnf(:,:,ns)**2))
691  END IF
692 
693 #if defined(MPI_OPT)
694  CALL mpi_allreduce(mpi_in_place, tmps, 10, mpi_real8, mpi_sum, &
695  siesta_comm, mpi_err)
696 #endif
697 
698  IF (iam .eq. 0 .and. lverbose) THEN
699  IF (parity .eq. f_cos) THEN
700  WRITE (*,1000)
701  ELSE
702  WRITE (*,1001)
703  END IF
704  WRITE (*,1002) sqrt(sum(f_smnf(m1+moff,:,1)**2)), tmps(1), tmps(2)
705  WRITE (*,1003) sqrt(sum(f_umnf(m1+moff,:,1)**2)), tmps(3), tmps(4)
706  WRITE (*,1004) sqrt(sum(f_vmnf(m0+moff,:,1)**2)), tmps(5), tmps(6)
707  WRITE (*,1005) tmps(9), tmps(7)
708  WRITE (*,1006) tmps(10), tmps(8)
709  END IF
710  END IF print_o
711 
712 ! Origin boundary conditions. 1/2 factors below are needed to preserve hessian
713 ! symmetry.
714  IF (nsmin .eq. 1) THEN
715 
716  IF (.not.l_push_s) THEN
717  f_smnf(m1 + moff,:,1) = 0
718  ELSE
719  f_smnf(m1 + moff,:,1) = p5*f_smnf(m1 + moff,:,1)
720  END IF
721 
722 ! IF (.not.l_push_u) THEN
723 ! f_umnf(m1 + moff,:,1) = 0
724 ! ELSE
725 ! f_umnf(m1 + moff,:,1) = p5*f_umnf(m1 + moff,:,1)
726 ! END IF
727 
728  IF (.not.l_push_v) THEN
729  f_vmnf(:m1 + moff,:,1) = 0
730  ELSE
731  f_vmnf(:m1 + moff,:,1) = p5*f_vmnf(:m1 + moff,:,1)
732  END IF
733  END IF
734 
735 ! Edge boundary conditions.
736  edge: IF (nsmax .EQ. ns) THEN
737  f_smnf(:,:,ns) = 0
738  IF (l_pedge) THEN
739  f_umnf(:,:,ns) = 0
740  END IF
741  IF (.not.l_push_edge) THEN
742  f_umnf(:,:,ns) = 0
743  f_vmnf(:,:,ns) = 0
744  ELSE
745  f_umnf(:,:,ns) = p5*f_umnf(:,:,ns)
746  f_vmnf(:,:,ns) = p5*f_vmnf(:,:,ns)
747  END IF
748  END IF edge
749 
750  CALL second0(skstoff)
751  get_force_harmonics_time = get_force_harmonics_time + (skstoff - skston)
752 
753 1000 FORMAT(/,1x,'SYMMETRIC FORCES: ')
754 1001 FORMAT(/,1x,'ASYMMETRIC FORCES: ')
755 1002 FORMAT(' fs(1,m=1): ',1p,e12.3, &
756  ' fs(2,m=1): ',1pe12.3,' fs(2,m!=1):',1pe12.3)
757 1003 FORMAT(' fu(1,m=1): ',1p,e12.3, &
758  ' fu(2,m=1): ',1pe12.3,' fu(2,m!=1):',1pe12.3)
759 1004 FORMAT(' fv(1,m=0): ',1p,e12.3, &
760  ' fv(2,m=0): ',1pe12.3,' fv(2,m>0): ',1pe12.3)
761 1005 FORMAT(' fu(ns-1) : ',1p,e12.3,' fu(ns) : ',1pe12.3)
762 1006 FORMAT(' fv(ns-1) : ',1p,e12.3,' fv(ns) : ',1pe12.3,/)
763 
764  END SUBROUTINE
765 
766 !-------------------------------------------------------------------------------
782 !-------------------------------------------------------------------------------
783  SUBROUTINE addpardamping(pardamp, f_smnf, f_umnf, f_vmnf, &
784  jvsupsmnf, jvsupumnf, jvsupvmnf)
785  USE shared_data, ONLY: col_scale
786 
787  IMPLICIT NONE
788 
789 ! Declare Arguments
790  REAL (dp), DIMENSION(:,:), INTENT(in) :: pardamp
791  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: f_smnf
792  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: f_umnf
793  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: f_vmnf
794  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: jvsupsmnf
795  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: jvsupumnf
796  REAL (dp), DIMENSION(:,:,:), INTENT(in) :: jvsupvmnf
797 
798 ! local variables
799  INTEGER :: js
800 
801 ! Start of executable code
802 ! Update mupar_norm in cv_currents and put in 1/jacobh factor. Simpler, faster
803 ! approcimately diagonal form.
804  DO js = 1, SIZE(pardamp,1)
805  f_smnf(:,:,js) = f_smnf(:,:,js) + jvsupsmnf(:,:,js)*pardamp(js,1)
806  f_umnf(:,:,js) = f_umnf(:,:,js) + jvsupumnf(:,:,js)*pardamp(js,2)
807  f_vmnf(:,:,js) = f_vmnf(:,:,js) + jvsupvmnf(:,:,js)*pardamp(js,3)
808  END DO
809 
810  END SUBROUTINE
811 
812  END MODULE
shared_data::nprecon
integer nprecon
Preconditioner flag.
Definition: shared_data.f90:68
shared_data::l_push_edge
logical l_push_edge
Solve u,v components at s=1.
Definition: shared_data.f90:202
shared_data::lverbose
logical lverbose
Use verbose screen output.
Definition: shared_data.f90:234
shared_data::l_push_u
logical l_push_u
Solve for u component at origin.
Definition: shared_data.f90:206
siesta_force::get_force_harmonics
subroutine get_force_harmonics(pardamp, f_smnf, f_umnf, f_vmnf, pa
Final computation of the MHD covariant Fourier force components.
Definition: siesta_force.f90:595
shared_data::l_pedge
logical, parameter l_pedge
Preserve s=1 as iso-pressure surface.
Definition: shared_data.f90:43
quantities
This file contains subroutines for allocating and initializing curvilinear magnetic covariant and pre...
Definition: quantities.f90:11
shared_data::l_push_s
logical l_push_s
Solve for s component at origin.
Definition: shared_data.f90:204
siesta_force::lorentz
subroutine lorentz(bsupsijf, bsupuijf, bsupvijf,
Compute covariant (sub) real-space components of the Lorentz K X B force.
Definition: siesta_force.f90:348
v3_utilities::assert_eq
Definition: v3_utilities.f:62
siesta_force::incfields
subroutine incfields(jbsupsmnh, jbsupumnh, jbsupvmnh, jpmnh,
Compute nonlinear or linearized contravariant magnetic field components and pressure in real space.
Definition: siesta_force.f90:231
island_params::fourier_context
type(fourier_class), pointer fourier_context
Fourier transform object.
Definition: island_params.f90:76
shared_data::l_printoriginforces
logical l_printoriginforces
Print forces at the origin.
Definition: shared_data.f90:220
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
shared_data::col_scale
real(dp), dimension(:,:,:,:), allocatable col_scale
Column scaling factors.
Definition: shared_data.f90:110
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
siesta_force::update_force
subroutine update_force
Update MHD forces on full radial mesh.
Definition: siesta_force.f90:70
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
shared_data::l_linearize
logical l_linearize
Use linearized forces.
Definition: shared_data.f90:210
siesta_force
Compute the JxB - Grad(p) covariant force components. The plasma is in equilibrium when the force of ...
Definition: siesta_force.f90:37
island_params::hs_i
real(dp) hs_i
Radial grid spacing. hs = s_i+1 - s_i.
Definition: island_params.f90:47
fourier::m2
integer, parameter m2
m = 2 mode.
Definition: fourier.f90:62
shared_data::gc
real(dp), dimension(:), allocatable, target gc
1D Array of Fourier mode MHD force components
Definition: shared_data.f90:99
siesta_force::initpardamping
subroutine initpardamping(parscale, pardamp)
Compute scaling factors for parallel flow damping.
Definition: siesta_force.f90:536
siesta_force::addpardamping
subroutine addpardamping(pardamp, f_smnf, f_umnf, f_vmnf,
Add parallel flow damping terms to forces (for Hessian calculation) to suppress null space at resonan...
Definition: siesta_force.f90:784
shared_data::l_getwmhd
logical l_getwmhd
Compute MHD energy.
Definition: shared_data.f90:214
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
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
siesta_currents::cv_currents
subroutine cv_currents(bsupsijh, bsupuijh, bsupvijh,
Compute current density from Curl(B).
Definition: siesta_currents.f90:62
siesta_currents
Takes the Curl(B) to return the contravariant current density. Contravariant components of the magnet...
Definition: siesta_currents.f90:12
fourier::f_cos
integer, parameter f_cos
Cosine parity.
Definition: fourier.f90:25
siesta_force::getmhdforce
subroutine getmhdforce(fsubsmnf, fsubumnf, fsubvmnf, pijh,
Compute covariant (sub) components of the MHD force J X B - grad(p).
Definition: siesta_force.f90:429
shared_data
This file contains variables and parameters used by many modules in SIESTA.
Definition: shared_data.f90:10
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
shared_data::l_push_v
logical l_push_v
Solve for v component at origin.
Definition: shared_data.f90:208
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17