V3FIT
siesta_currents.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 siesta_currents
13  USE v3_utilities, ONLY: assert
14  USE stel_kinds
15  USE stel_constants
16  USE descriptor_mod, ONLY: siesta_comm
17  USE shared_data, ONLY: l_update_state, lasym, ste, &
18  jsju_ratio_s, jsju_ratio_a
19  USE fourier
20  USE nscalingtools, ONLY: nranks, startglobrow, endglobrow, mpi_err
21  USE metrics, ONLY: tolowerf, tolowerh
22  USE quantities
23 
24  IMPLICIT NONE
25 
26  CONTAINS
27 
28 !*******************************************************************************
29 ! UTILITY SUBROUTINES
30 !*******************************************************************************
31 !-------------------------------------------------------------------------------
60 !-------------------------------------------------------------------------------
61  SUBROUTINE cv_currents(bsupsijh, bsupuijh, bsupvijh, &
62  ksupsijf, ksupuijf, ksupvijf, &
63  lmagen, lcurr)
64  USE island_params, ONLY: nu_i, rmajor_i, fourier_context
65  USE hessian, ONLY: mupar_norm
66  USE descriptor_mod, ONLY: iam
67  USE quantities, ONLY: bsq
68  USE timer_mod
69  USE mpi_inc
70  USE descriptor_mod, ONLY: siesta_comm
71  USE shared_data, ONLY: siesta_curtor, unit_out
72 
73  IMPLICIT NONE
74 
75 ! Declare Arguments
76  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: bsupsijh
77  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: bsupuijh
78  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: bsupvijh
79  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupsijf
80  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupuijf
81  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupvijf
82  LOGICAL, INTENT(in) :: lmagen
83  LOGICAL, INTENT(in) :: lcurr
84 
85 ! Local Variables
86  INTEGER :: istat
87  REAL (dp) :: beta
88  REAL (dp) :: ton
89  REAL (dp) :: toff
90  REAL (dp), DIMENSION(6) :: temp(6)
91  REAL (dp) :: tmp
92  INTEGER :: nsmin
93  INTEGER :: nsmax
94  INTEGER :: nmin
95  INTEGER :: nmax
96  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubsijh
97  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubuijh
98  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubvijh
99  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bs_filter
100  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bu_filter
101  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bv_filter
102  REAL (dp) :: edge_cur
103 
104 ! Start of executable code
105  CALL second0(ton)
106 
107  nmin = max(1, startglobrow)
108  nmax = min(ns, endglobrow)
109  nsmin = lbound(bsupsijh,3)
110  nsmax = ubound(bsupsijh,3)
111 
112 ! Calculate covariant BsubXh on half mesh
113  ALLOCATE(bsubsijh(ntheta,nzeta,nsmin:nsmax), &
114  bsubuijh(ntheta,nzeta,nsmin:nsmax), &
115  bsubvijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
116  CALL assert_eq(0, istat, 'Allocation1 failed in cv_currents')
117 
118  CALL tolowerh(bsupsijh, bsupuijh, bsupvijh, &
119  bsubsijh, bsubuijh, bsubvijh, nsmin, nsmax)
120 
121 ! Compute parallel damping scaling cofficient <bsupv^2/B.B> ~ k^2
122  IF (.not.ALLOCATED(mupar_norm)) THEN
123  ALLOCATE(mupar_norm(nsmin:nsmax), stat=istat)
124  CALL assert_eq(0, istat,'mupar_norm allocation failed!')
125  beta = max(1.e-5_dp, wp_i/wb_i)
126  mupar_norm = signjac*beta/(rmajor_i*rmajor_i) !/ vp_f(nsmin:nsmax) !add 1/jac correction factor in AddDamping
127  IF (iam .eq. 0 .and. lverbose) THEN
128  WRITE (*,1000) rmajor_i
129  END IF
130  END IF
131 
132  IF (lmagen) THEN
133  bsq(:,:,nmin:nmax) = bsupsijh(:,:,nmin:nmax)*bsubsijh(:,:,nmin:nmax) &
134  + bsupuijh(:,:,nmin:nmax)*bsubuijh(:,:,nmin:nmax) &
135  + bsupvijh(:,:,nmin:nmax)*bsubvijh(:,:,nmin:nmax)
136  IF (nmin .EQ. 1) THEN
137  bsq(:,:,1) = 0
138  END IF
139 
140  CALL assert(nmin .eq. 1 .or. &
141  all(bsq(:,:,nmin:nmax) .gt. zero), &
142  'BSQ <= 0 in cv_currents!')
143 
144  wb = sum(bsq(:,:,nmin:nmax)*jacobh(:,:,nmin:nmax)*wint(:,:,nmin:nmax))
145 
146  IF (parsolver) THEN
147 #if defined(MPI_OPT)
148 ! FIXME: All reduce is not deterministic. This causes a divergent run sequence.
149  CALL mpi_allreduce(mpi_in_place, wb, 1, mpi_real8, mpi_sum, &
150  siesta_comm, mpi_err)
151 #endif
152  END IF
153  wb = signjac*(twopi*pi*hs_i)*wb
154  END IF
155 
156 ! Calculate BsubXmnh and KsupXmnf harmonics. Deallocated in calling routine.
157  IF (l_update_state) THEN
158  ALLOCATE(bs_filter(ntheta,nzeta,nsmin:nsmax), &
159  bu_filter(ntheta,nzeta,nsmin:nsmax), &
160  bv_filter(ntheta,nzeta,nsmin:nsmax), stat=istat)
161  CALL assert_eq(0, istat, 'Allocation2 failed in CURLB')
162  END IF
163 
164  edge_cur = 0.0
165  CALL curlb(ksupsmnsf, ksupumncf, ksupvmncf, &
166  bsubsijh, bsubuijh, bsubvijh, &
167  bs_filter, bu_filter, bv_filter, edge_cur, &
168  jsju_ratio_s, f_sin, lcurr, nsmax, nsmin, siesta_curtor)
169  IF (lasym) THEN
170  CALL curlb(ksupsmncf, ksupumnsf, ksupvmnsf, &
171  bsubsijh, bsubuijh, bsubvijh, &
172  bs_filter, bu_filter, bv_filter, edge_cur, &
173  jsju_ratio_a, f_cos, lcurr, nsmax, nsmin, &
174  siesta_curtor)
175  END IF
176 
177  IF (nsmax .eq. ns) THEN
178  siesta_curtor = -fourier_context%orthonorm(m0,n0)*siesta_curtor &
179  / b_factor
180  END IF
181 
182  IF (lcurr) THEN
183  nsmax = min(ns, endglobrow + 1)
184  ELSE
185  nsmax = min(ns, endglobrow)
186  END IF
187 
188 ! Calculate full mesh, contravariant real-space current components
189 ! KsupXF = sqrt(g)*JsupXF.
190  CALL getcv(ksupsmnsf, ksupumncf, ksupvmncf, &
191  ksupsijf, ksupuijf, ksupvijf, f_sin, nsmin, nsmax)
192  IF (lasym) THEN
193  CALL getcv(ksupsmncf, ksupumnsf, ksupvmnsf, &
194  ksupsijf, ksupuijf, ksupvijf, f_cos, nsmin, nsmax)
195  END IF
196 
197 ! Covariant K components (with jacob factor still there) in real space needed
198 ! for current diffusion
199  IF (lcurr) THEN
200  CALL tolowerf(ksupsijf, ksupuijf, ksupvijf, &
201  ksubsijf, ksubuijf, ksubvijf, nsmin, nsmax)
202  IF (nsmin .eq. 1) THEN
203  ksubuijf(:,:,1) = 0
204  END IF
205  END IF
206 
207 ! Diagnostic output (only compute when state is updated)
208  diagno: IF (l_update_state) THEN
209 
210 ! Compute spectral truncation error measure (Use ksupsijf for temporary
211 ! filtered bsubXijh values)
212  IF (nsmax .eq. ns) THEN
213  tmp = sum(wint(:,:,ns)*(bsubvijh(:,:,ns)*bsubvijh(:,:,ns) + &
214  bsubuijh(:,:,ns)*bsubuijh(:,:,ns)))
215  IF (tmp .ne. zero) THEN
216  tmp = sqrt(edge_cur/tmp)
217  END IF
218  END IF
219 #if defined(MPI_OPT)
220  IF (parsolver) THEN
221  CALL mpi_bcast(tmp, 1, mpi_real8, nranks - 1, siesta_comm, mpi_err)
222  END IF
223 #endif
224  IF (iam .EQ. 0) THEN
225  IF (lverbose) THEN
226  WRITE (*, 1001) tmp
227  END IF
228  WRITE (unit_out, 1001) tmp
229  CALL flush(unit_out)
230  END IF
231 
232  temp(1) = sum(wint(:,:,nmin:nmax) * &
233  (bs_filter - bsubsijh(:,:,nmin:nmax))**2)
234  temp(2) = sum(wint(:,:,nmin:nmax) * &
235  (bs_filter + bsubsijh(:,:,nmin:nmax))**2)
236  temp(3) = sum(wint(:,:,nmin:nmax) * &
237  (bu_filter - bsubuijh(:,:,nmin:nmax))**2)
238  temp(4) = sum(wint(:,:,nmin:nmax) * &
239  (bu_filter + bsubuijh(:,:,nmin:nmax))**2)
240  temp(5) = sum(wint(:,:,nmin:nmax) * &
241  (bv_filter - bsubvijh(:,:,nmin:nmax))**2)
242  temp(6) = sum(wint(:,:,nmin:nmax) * &
243  (bv_filter + bsubvijh(:,:,nmin:nmax))**2)
244 #if defined(MPI_OPT)
245  IF (parsolver) THEN
246  CALL mpi_allreduce(mpi_in_place, temp, 6, mpi_real8, mpi_sum, &
247  siesta_comm, mpi_err)
248  END IF
249 #endif
250  IF (temp(2) .ne. zero) THEN
251  ste(2) = sqrt(temp(1)/temp(2))
252  END IF
253  IF (temp(4) .ne. zero) THEN
254  ste(3) = sqrt(temp(3)/temp(4))
255  END IF
256  IF (temp(6) .ne. zero) THEN
257  ste(4) = sqrt(temp(5)/temp(6))
258  END IF
259  DEALLOCATE (bs_filter, bu_filter, bv_filter)
260  END IF diagno
261 
262  DEALLOCATE (bsubsijh, bsubuijh, bsubvijh)
263 
264  CALL second0(toff)
265  cv_current_time = cv_current_time + (toff - ton)
266  time_current = time_current + (toff - ton)
267 
268 1000 FORMAT(' Rmajor_i = ',1p,e12.3)
269 1001 FORMAT(' RMS EDGE CURRENT ||KSUPS(NS)|| = ',1p,e10.3)
270 
271  END SUBROUTINE
272 
273 !-------------------------------------------------------------------------------
296 !-------------------------------------------------------------------------------
297  SUBROUTINE curlb(ksupsmnf, ksupumnf, ksupvmnf, &
298  bsubsijh, bsubuijh, bsubvijh, &
299  bs_filter, bu_filter, bv_filter, edge_cur, &
300  jsju_ratio, iparity, lcurr, nsmax, nsmin, curtor)
301  USE utilities, ONLY: curl_htof
302  USE island_params, ONLY: fourier_context
303 
304  IMPLICIT NONE
305 
306 ! Declare Arguments
307  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: ksupsmnf
308  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: ksupumnf
309  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: ksupvmnf
310  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(in) :: bsubsijh
311  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(in) :: bsubuijh
312  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: bsubvijh
313  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: bs_filter
314  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: bu_filter
315  REAL (dp), DIMENSION(:,:,:), INTENT(inout) :: bv_filter
316  REAL (dp), INTENT(inout) :: edge_cur
317  REAL (dp), INTENT(out) :: jsju_ratio
318  INTEGER, INTENT(in) :: iparity
319  LOGICAL, INTENT(in) :: lcurr
320  INTEGER, INTENT(in) :: nsmax
321  INTEGER, INTENT(in) :: nsmin
322  REAL (dp), INTENT(inout) :: curtor
323 
324 ! local variables
325  INTEGER :: fours
326  INTEGER :: fouruv
327  INTEGER :: fcomb
328  INTEGER :: nmin
329  INTEGER :: nmax
330  INTEGER :: istat
331  REAL (dp) :: r12
332  REAL (dp), DIMENSION(-ntor:ntor) :: temp
333  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubsmnh
334  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubumnh
335  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsubvmnh
336 
337 ! Start of executable code
338  nmin = max(1, startglobrow)
339  nmax = min(ns, endglobrow)
340 
341  IF (iparity .EQ. f_sin) THEN
342  fcomb = f_none
343  fours = f_sin
344  fouruv = f_cos
345  r12 = 0.5*hs_i
346  ELSE
347  fcomb = f_sum
348  fours = f_cos
349  fouruv = f_sin
350  r12 = -0.5*hs_i
351  END IF
352 
353  ALLOCATE(bsubsmnh(0:mpol,-ntor:ntor,nsmin:nsmax), &
354  bsubumnh(0:mpol,-ntor:ntor,nsmin:nsmax), &
355  bsubvmnh(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
356  CALL assert_eq(0, istat, 'Allocation1 failed in CURLB')
357 
358 ! Compute bsubxmnh
359  CALL fourier_context%tomnsp(bsubsijh, bsubsmnh, fours)
360  CALL fourier_context%tomnsp(bsubuijh, bsubumnh, fouruv)
361  CALL fourier_context%tomnsp(bsubvijh, bsubvmnh, fouruv)
362 
363  CALL curl_htof(bsubsmnh, bsubumnh, bsubvmnh, &
364  ksupsmnf, ksupumnf, ksupvmnf, &
365  iparity, nsmin, nsmax, lcurr, curtor)
366 
367 ! Diagnostics
368  diagno: IF (l_update_state) THEN
369  IF (nmax .eq. ns) THEN
370  edge_cur = edge_cur + sum(ksupsmnf(:,:,ns)*ksupsmnf(:,:,ns))
371  END IF
372 
373  CALL fourier_context%toijsp(bsubsmnh(:,:,nmin:nmax), bs_filter, &
374  fcomb, fours)
375  CALL fourier_context%toijsp(bsubumnh(:,:,nmin:nmax), bu_filter, &
376  fcomb, fouruv)
377  CALL fourier_context%toijsp(bsubvmnh(:,:,nmin:nmax), bv_filter, &
378  fcomb, fouruv)
379  IF (nmin .eq. 1) THEN
380  bv_filter(:,:,1) = bv_filter(:,:,2)
381  bsubvijh(:,:,1) = bsubvijh(:,:,2)
382  END IF
383 
384 ! K^s-r12*K^u at origin BEFORE bc applied
385  IF (nsmin .eq. 1) THEN
386  temp = ksupsmnf(m1,:,1) + r12*ksupumnf(m1,:,1)
387  jsju_ratio = sum(temp*temp)
388  IF (jsju_ratio .gt. zero) THEN
389  temp = ksupsmnf(m1,:,1) - r12*ksupumnf(m1,:,1)
390  jsju_ratio = sum(temp*temp)/jsju_ratio
391  jsju_ratio = sqrt(jsju_ratio)
392  END IF
393  END IF
394  END IF diagno
395 
396  DEALLOCATE (bsubsmnh, bsubumnh, bsubvmnh, stat=istat)
397  CALL assert_eq(0, istat, 'Deallocation failed in CURLB')
398 
399  END SUBROUTINE
400 
401 !-------------------------------------------------------------------------------
413 !-------------------------------------------------------------------------------
414  SUBROUTINE getcv(ksupsmnf, ksupumnf, ksupvmnf, &
415  ksupsijf, ksupuijf, ksupvijf, parity, nsmin, nsmax)
416  USE island_params, ONLY: fourier_context
417 
418  IMPLICIT NONE
419 
420 ! Declare Arguments
421  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: ksupsmnf
422  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: ksupumnf
423  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: ksupvmnf
424  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupsijf
425  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupuijf
426  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(inout) :: ksupvijf
427  INTEGER, INTENT(in) :: parity
428  INTEGER, INTENT(in) :: nsmin
429  INTEGER, INTENT(in) :: nsmax
430 
431 ! local variables
432  INTEGER :: fours
433  INTEGER :: fouruv
434  INTEGER :: fcomb
435 
436 ! Start of executable code
437  IF (parity .eq. f_sin) THEN
438  fcomb = f_none
439  fours = f_sin
440  fouruv = f_cos
441  ELSE
442  fcomb = f_sum
443  fours = f_cos
444  fouruv = f_sin
445  END IF
446 
447  CALL fourier_context%toijsp(ksupsmnf(:,:,nsmin:nsmax), ksupsijf, fcomb, &
448  fours)
449  CALL fourier_context%toijsp(ksupumnf(:,:,nsmin:nsmax), ksupuijf, fcomb, &
450  fouruv)
451  CALL fourier_context%toijsp(ksupvmnf(:,:,nsmin:nsmax), ksupvijf, fcomb, &
452  fouruv)
453 
454  END SUBROUTINE
455 
456  END MODULE
quantities
This file contains subroutines for allocating and initializing curvilinear magnetic covariant and pre...
Definition: quantities.f90:11
v3_utilities::assert_eq
Definition: v3_utilities.f:62
shared_data::siesta_curtor
real(dp) siesta_curtor
Total toroidal current.
Definition: shared_data.f90:239
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
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
island_params::rmajor_i
real(dp) rmajor_i
Major radius.
Definition: island_params.f90:64
shared_data::jsju_ratio_s
real(dp) jsju_ratio_s
FIXME: UNKNOWN.
Definition: shared_data.f90:141
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::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
shared_data::l_update_state
logical l_update_state
Update the ste array.
Definition: shared_data.f90:224
siesta_currents::curlb
subroutine curlb(ksupsmnf, ksupumnf, ksupvmnf,
Compute fourier components of contravariant currents on full radial grid.
Definition: siesta_currents.f90:298
metrics::tolowerf
subroutine tolowerf(xsupsij, xsupuij, xsupvij,
Converts to covariant component full grid.
Definition: metrics.f90:590
shared_data::ste
real(dp), dimension(4) ste
Spectral Truncation RMS error,.
Definition: shared_data.f90:133
siesta_currents::getcv
subroutine getcv(ksupsmnf, ksupumnf, ksupvmnf,
Transform full grid contravariant currents to real space.
Definition: siesta_currents.f90:415
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
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
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
shared_data::jsju_ratio_a
real(dp) jsju_ratio_a
FIXME: UNKNOWN.
Definition: shared_data.f90:145
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
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17