V3FIT
siesta_init.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 !
9 !*******************************************************************************
10  MODULE siesta_init
11  USE v3_utilities, ONLY: assert, assert_eq
12  USE stel_kinds
13  USE descriptor_mod, ONLY: siesta_comm, iam
14  USE fourier, ONLY: f_none, f_sin, f_cos, f_du, f_dv, f_sum, m0, m1
15  USE island_params, ONLY: fourier_context
16  USE nscalingtools, ONLY: startglobrow, endglobrow, mpi_err
17  USE quantities
18  USE shared_data, ONLY: l_update_state, l_pedge
19  USE mpi_inc
20 
21  IMPLICIT NONE
22 
23  PRIVATE
24 !*******************************************************************************
25 ! utilities module variables
26 !*******************************************************************************
27 ! FIXME: Remove these if possible.
28  INTEGER :: nsmin
29  INTEGER :: nsmax
30  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: pfilter
31  REAL (dp) :: pedge
32 
33  PUBLIC :: init_state
34 
35  CONTAINS
36 
37 !-------------------------------------------------------------------------------
43 !-------------------------------------------------------------------------------
44  SUBROUTINE init_state(lcurrent_only, lpar_in)
45  USE bhtobf
46  USE utilities, ONLY: gradientfull
47  USE metrics, ONLY: tolowerf
48  USE shared_data, ONLY: fsq_total, l_getfsq, ste, buv_res, &
49  l_update_state
50  USE hessian, ONLY: mupar
51  USE siesta_currents, ONLY: cv_currents
52  USE island_params, ONLY: ohs=>ohs_i
53 
54  IMPLICIT NONE
55 
56 ! Declare Arguments
57  LOGICAL, INTENT(IN) :: lcurrent_only
58  LOGICAL, INTENT(IN), OPTIONAL :: lpar_in
59 
60 ! local variables
61  INTEGER :: istat
62  INTEGER :: n1
63  INTEGER :: n2
64  INTEGER :: nmax
65  INTEGER :: endr
66  INTEGER :: startr
67  REAL (dp) :: ton
68  REAL (dp) :: toff
69  REAL (dp) :: temp(2)
70  LOGICAL :: ladd_pert0
71  LOGICAL :: lpar
72 
73 ! Start of executable code
74 ! jbupXmnPh and jpmnPh are initially computed on full ns mesh, and are gathered
75 ! to that full ns mesh in update_state
76  istat = 0
77  IF (PRESENT(lpar_in)) THEN
78  lpar = lpar_in
79  ELSE
80  lpar = .true.
81  END IF
82 
83  IF (.not.lpar) THEN
84  startr = startglobrow
85  endr = endglobrow
86  startglobrow = 1
87  endglobrow = ns
88  END IF
89  nsmin = max(1, startglobrow - 1)
90  nsmax = min(ns, endglobrow + 2)
91 
92 !Unnecessary to recompute magnetic fields or currents: already stored at this point!
93  ladd_pert0 = ALLOCATED(buv_res)
94  IF (.not.ladd_pert0) THEN
95  CALL init_allocate_arrays(lpar)
96 
97 ! Calculate unperturbed half mesh jac*bsupX and jac*pressure from nsmin to
98 ! nsmax. Sum symmetric and (for lasym) asymmetric parts. Returning from
99 ! update_state, jbsupXmn, jpmn have been gathered onto all processors.
100  CALL getfields(jbsupsmnsh, jbsupumnch, jbsupvmnch, jpmnch, f_sin)
101  IF (lasym) THEN
102  CALL getfields(jbsupsmnch, jbsupumnsh, jbsupvmnsh, jpmnsh, f_cos)
103  END IF
104 
105 ! Convert jac*bsubXh to bsupXh and jac*pres to presh in real space. On exit,
106 ! bsupXf and pf0 ae valid on [nsmin, nsmax - 1]
107  CALL bhalftobfull(bsupsijh0, bsupuijh0, bsupvijh0, &
108  bsupsijf0, bsupuijf0, bsupvijf0, &
109  pijh0, pijf0)
110 
111 ! Compute and store unperturbed current components on full mesh,
112 ! KsupX = sqrt(g)*JsupX
113  CALL cv_currents(bsupsijh0, bsupuijh0, bsupvijh0, &
114  ksupsijf0, ksupuijf0, ksupvijf0, &
115  l_getfsq, .true.)
116 
117  IF (lcurrent_only) THEN
118  RETURN
119  END IF
120  END IF
121 
122 ! Need this for resistive diffusion (|| resistivity) and mupar damping
123  IF (ladd_pert0 .or. mupar .ne. zero) THEN
124  nmax = min(endglobrow + 1, ns)
125 ! Calculate K || B
126  CALL tolowerf(bsupsijf0,bsupuijf0,bsupvijf0, &
127  bsubsijf, bsubuijf, bsubvijf, &
128  nsmin, nmax)
129 ! IN UPDATE-BFIELD, KSUB = jacob*JSUB: PARALLEL CURRENT
130  IF (ladd_pert0) THEN
131  ksubsijf(:,:,nsmin:nmax) = bsubsijf(:,:,nsmin:nmax)
132  ksubuijf(:,:,nsmin:nmax) = bsubuijf(:,:,nsmin:nmax)
133  ksubvijf(:,:,nsmin:nmax) = bsubvijf(:,:,nsmin:nmax)
134  RETURN
135  END IF
136  END IF
137 
138 ! Compute spectrally filtered dp/du and dp/dv on half radial grid [nsmin:nsmax]
139  IF (l_update_state) THEN
140  ALLOCATE(pfilter(ntheta, nzeta, nsmin:nsmax))
141  END IF
142  CALL filterpressure(f_cos)
143  IF (lasym) THEN
144  CALL filterpressure(f_sin)
145  END IF
146 
147 ! Compute radial pressure derivative pijf0_ds at full mesh points
148 ! [nsmin:nsmax-1]
149  CALL gradientfull(pijf0_ds, pijh0)
150 
151 ! SPH: one-sided (m=1) derivative at origin yields factor of 2 (1/(hs/2)).
152 ! pijh(:,:,1) should contain the value of jpmnch(m0,:,2) from Init_Fields.
153  IF (nsmin .EQ. 1) THEN
154  pijf0_ds(:,:,1) = 2.0*(pijh0(:,:,2) - pijh0(:,:,1))*ohs
155  END IF
156 
157 ! SPH11-3-16 preserve s=1 as iso-pressure contour
158  IF (nsmax .eq. ns .and. l_pedge) THEN
159  pijf0(:,:,ns) = pedge
160  pijf0_ds(:,:,ns) = 0
161  END IF
162 
163 ! Fourier filter (Nyquist) check
164  IF (.not. l_update_state) THEN
165  RETURN
166  END IF
167 
168 ! pijh is unfiltered
169  n1 = max(1, startglobrow)
170  n2 = min(ns, endglobrow)
171  IF (n1 .EQ. 1) THEN
172  pfilter(:,:,1) = 0
173  END IF
174  temp(1) = sum((pijh0(:,:,n1:n2) - pfilter(:,:,n1:n2))**2)
175  temp(2) = sum((pijh0(:,:,n1:n2) + pfilter(:,:,n1:n2))**2)
176 #if defined(MPI_OPT)
177  CALL mpi_allreduce(mpi_in_place, temp, 2, mpi_real8, mpi_sum, &
178  siesta_comm, mpi_err)
179 #endif
180  IF (temp(2) .ne. zero) THEN
181  ste(1) = sqrt(temp(1)/temp(2))
182  END IF
183  DEALLOCATE(pfilter)
184 
185  CALL checkfft(jbsupsmnsh, jbsupumnch, jbsupvmnch, jpmnch, f_sin)
186  IF (lasym) THEN
187  CALL checkfft(jbsupsmnch, jbsupumnsh, jbsupvmnsh, jpmnsh, f_cos)
188  END IF
189 
190  IF (.not.lpar) THEN
191  startglobrow = startr
192  endglobrow = endr
193  END IF
194 
195  END SUBROUTINE
196 
198  SUBROUTINE init_allocate_arrays(lpar)
199  USE shared_data, ONLY: l_par_state
200 !-----------------------------------------------
201 ! D u m m y A r g u m e n t s
202 !-----------------------------------------------
203  LOGICAL, INTENT(in) :: lpar
204 !-----------------------------------------------
205 ! L o c a l V a r i a b l e s
206 !-----------------------------------------------
207  LOGICAL :: lrealloc, lalloc
208  INTEGER :: istat, n1, n2, n3, nsmin, nsmax, &
209  nsmin1, nsmax1
210 !-----------------------------------------------
211  lalloc = ALLOCATED(jvsupsijf)
212  lrealloc = .true.
213  istat = 0
214  nsmin =max(1,startglobrow); nsmax =min(endglobrow+1,ns)
215  nsmin1=max(1,startglobrow-1); nsmax1=min(endglobrow+2,ns)
216  n1 = ntheta; n2 = nzeta
217 
218  IF (lpar) THEN
219  n3 = nsmax1-nsmin1+1
220  l_par_state=.true.
221  ELSE
222  n3 = ns
223  l_par_state=.false.
224  END IF
225 
226  IF (lalloc) THEN
227  IF (SIZE(bsupsijh0,3) .EQ. n3) lrealloc=.false.
228  IF (lrealloc) THEN
229  DEALLOCATE (jvsupsijf, jvsupuijf, jvsupvijf, &
230  bsupsijf0, bsupuijf0, bsupvijf0, &
231  bsupsijf, bsupuijf, bsupvijf, &
232  bsupsijh0, bsupuijh0, bsupvijh0, &
233  bsubsijf, bsubuijf, bsubvijf, &
234  ksupsijf0, ksupuijf0, ksupvijf0, &
235  ksupsijf, ksupuijf, ksupvijf, &
236  ksubsijf, ksubuijf, ksubvijf, &
237  pijf0, pijh0, pijf0_ds, &
238  pijh0_du, pijh0_dv, stat=istat)
239  CALL assert_eq(0, istat, &
240  'Deallocation error #1 in INIT_ALLOCATE_ARRAYS')
241  END IF
242  END IF
243 
244  IF (lrealloc) THEN
245  n3 = ns
246  ALLOCATE (jvsupsijf(n1,n2,n3), jvsupuijf(n1,n2,n3), & ! Full mesh quantities (real space)
247  jvsupvijf(n1,n2,n3)) ! jacobian*(V_s , V_u , V_v)
248  IF (lpar) THEN
249  ALLOCATE (bsupsijh0(n1,n2,nsmin1:nsmax1), & ! B^s, B^u, B^v (half)
250  bsupuijh0(n1,n2,nsmin1:nsmax1), &
251  bsupvijh0(n1,n2,nsmin1:nsmax1), &
252  bsupsijf0(n1,n2,nsmin1:nsmax1), & ! B^s, B^u, B^v (full)
253  bsupuijf0(n1,n2,nsmin1:nsmax1), &
254  bsupvijf0(n1,n2,nsmin1:nsmax1), &
255  ksupsijf0(n1,n2,nsmin1:nsmax), & ! K^s, K^u, K^v
256  ksupuijf0(n1,n2,nsmin1:nsmax), &
257  ksupvijf0(n1,n2,nsmin1:nsmax), &
258  ksubsijf(n1,n2,nsmin1:nsmax), & ! Full mesh quantities (real space)
259  ksubuijf(n1,n2,nsmin1:nsmax), & ! K_s, K_u, K_v
260  ksubvijf(n1,n2,nsmin1:nsmax), &
261  bsubsijf(n1,n2,nsmin1:nsmax), & ! Full mesh quantities (real space)
262  bsubuijf(n1,n2,nsmin1:nsmax), & ! B_s, B_u, B_v
263  bsubvijf(n1,n2,nsmin1:nsmax), &
264  ksupsijf(n1,n2,nsmin:nsmax), &
265  ksupuijf(n1,n2,nsmin:nsmax), &
266  ksupvijf(n1,n2,nsmin:nsmax), &
267  bsupsijf(n1,n2,nsmin:nsmax), &
268  bsupuijf(n1,n2,nsmin:nsmax), &
269  bsupvijf(n1,n2,nsmin:nsmax), &
270  pijh0(n1,n2,nsmin1:nsmax1), &
271  pijh0_du(n1,n2,nsmin1:nsmax), &
272  pijh0_dv(n1,n2,nsmin1:nsmax), &
273  pijf0(n1,n2,nsmin1:nsmax1), &
274  pijf0_ds(n1,n2,nsmin1:nsmax), stat=istat)
275  ELSE
276  ALLOCATE (ksupsijf0(n1,n2,n3), ksupuijf0(n1,n2,n3), &
277  ksupvijf0(n1,n2,n3), ksupsijf(n1,n2,n3), &
278  ksupuijf(n1,n2,n3), ksupvijf(n1,n2,n3), &
279  ksubsijf(n1,n2,n3), &
280  ksubuijf(n1,n2,n3), &
281  ksubvijf(n1,n2,n3), &
282  bsubsijf(n1,n2,n3), &
283  bsubuijf(n1,n2,n3), &
284  bsubvijf(n1,n2,n3), &
285  pijf0(n1,n2,n3), pijh0(n1,n2,n3), &
286  pijh0_du(n1,n2,n3), pijh0_dv(n1,n2,n3), &
287  pijf0_ds(n1,n2,n3), bsupsijf0(n1,n2,n3), &
288  bsupuijf0(n1,n2,n3),bsupvijf0(n1,n2,n3), &
289  bsupsijf(n1,n2,n3), bsupuijf(n1,n2,n3), &
290  bsupvijf(n1,n2,n3), bsupsijh0(n1,n2,n3), &
291  bsupuijh0(n1,n2,n3),bsupvijh0(n1,n2,n3), &
292  stat=istat)
293  END IF
294  CALL assert_eq(0, istat, &
295  'Allocation error #1 in Init_Allocate_Arrays')
296  jvsupsijf = 0; jvsupuijf = 0; jvsupvijf = 0 !Need in add_resistivity loop
297  END IF
298 
299 
300  lalloc = ALLOCATED(ksupsmnsf)
301  lrealloc = .true.
302 
303  IF (lpar) THEN
304  n1 = mpol+1; n2 = 2*ntor+1; n3 = nsmax-nsmin1+1
305  ELSE
306  n1 = mpol+1; n2 = 2*ntor+1; n3 = ns
307  END IF
308 
309  IF (lalloc) THEN
310  IF (SIZE(ksupsmnsf,3) .EQ. n3) THEN
311  lrealloc = .false.
312  END IF
313  IF (lrealloc) THEN
314  DEALLOCATE(ksupsmnsf, ksupumncf, ksupvmncf, &
315  djpmnch, djbsupsmnsh, djbsupumnch, djbsupvmnch, &
316  stat=istat)
317  CALL assert_eq(0, istat, &
318  'Deallocation error #2 in Init_Allocate_Arrays')
319  IF (lasym) THEN
320  DEALLOCATE(ksupsmncf, ksupumnsf, ksupvmnsf, &
321  djpmnsh, djbsupsmnch, djbsupumnsh, djbsupvmnsh, &
322  stat=istat)
323  CALL assert_eq(0, istat, &
324  'Deallocation error #3 in Init_Allocate_Arrays')
325  END IF
326  END IF
327  END IF
328 
329  IF (lrealloc) THEN
330  IF (lpar) THEN
331  ALLOCATE(ksupsmnsf(0:mpol,-ntor:ntor,nsmin1:nsmax), &
332  ksupumncf(0:mpol,-ntor:ntor,nsmin1:nsmax), &
333  ksupvmncf(0:mpol,-ntor:ntor,nsmin1:nsmax), &
334  djpmnch(0:mpol,-ntor:ntor,nsmin:nsmax), &
335  djbsupsmnsh(0:mpol,-ntor:ntor,nsmin:nsmax), &
336  djbsupumnch(0:mpol,-ntor:ntor,nsmin:nsmax), &
337  djbsupvmnch(0:mpol,-ntor:ntor,nsmin:nsmax), &
338  stat=istat)
339  IF (lasym) THEN
340  ALLOCATE(ksupsmncf(0:mpol,-ntor:ntor,nsmin1:nsmax), &
341  ksupumnsf(0:mpol,-ntor:ntor,nsmin1:nsmax), &
342  ksupvmnsf(0:mpol,-ntor:ntor,nsmin1:nsmax), &
343  djpmnsh(0:mpol,-ntor:ntor,nsmin:nsmax), &
344  djbsupsmnch(0:mpol,-ntor:ntor,nsmin:nsmax), &
345  djbsupumnsh(0:mpol,-ntor:ntor,nsmin:nsmax), &
346  djbsupvmnsh(0:mpol,-ntor:ntor,nsmin:nsmax), &
347  stat=istat)
348  END IF
349  ELSE
350  ALLOCATE(ksupsmnsf(0:mpol,-ntor:ntor,n3), &
351  ksupumncf(0:mpol,-ntor:ntor,n3), &
352  ksupvmncf(0:mpol,-ntor:ntor,n3), &
353  djpmnch(0:mpol,-ntor:ntor,n3), &
354  djbsupsmnsh(0:mpol,-ntor:ntor,n3), &
355  djbsupumnch(0:mpol,-ntor:ntor,n3), &
356  djbsupvmnch(0:mpol,-ntor:ntor,n3), &
357  stat=istat)
358  IF (lasym) THEN
359  ALLOCATE(ksupsmncf(0:mpol,-ntor:ntor,n3), &
360  ksupumnsf(0:mpol,-ntor:ntor,n3), &
361  ksupvmnsf(0:mpol,-ntor:ntor,n3), &
362  djpmnsh(0:mpol,-ntor:ntor,n3), &
363  djbsupsmnch(0:mpol,-ntor:ntor,n3), &
364  djbsupumnsh(0:mpol,-ntor:ntor,n3), &
365  djbsupvmnsh(0:mpol,-ntor:ntor,n3), &
366  stat=istat)
367  END IF
368  END IF
369  CALL assert_eq(0, istat, &
370  'Allocation error #2 in Init_Allocate_Arrays')
371  djpmnch = 0
372  djbsupsmnsh = 0
373  djbsupumnch = 0
374  djbsupvmnch = 0
375  IF (lasym) THEN
376  djpmnsh = 0
377  djbsupsmnch = 0
378  djbsupumnsh = 0
379  djbsupvmnsh = 0
380  END IF
381  END IF
382 
383  END SUBROUTINE
384 
385 !-------------------------------------------------------------------------------
395 !-------------------------------------------------------------------------------
396  SUBROUTINE getfields(jbsupsmnh, jbsupumnh, jbsupvmnh, jpmnh, iparity)
397 
398  IMPLICIT NONE
399 
400 ! Declare Arguments
401  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupsmnh
402  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupumnh
403  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupvmnh
404  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jpmnh
405  INTEGER, INTENT(in) :: iparity
406 
407 ! local variables
408  INTEGER :: fours
409  INTEGER :: fouruvp
410  INTEGER :: fcomb
411 
412 ! Start of executable code
413 ! Check origin values
414  IF (nsmin .EQ. 1) THEN
415  CALL assert(all(jbsupsmnh(m1,:,1) .eq. jbsupsmnh(m1,:,2)), &
416  'bmnsh(1) != bmnsh(2)')
417  CALL assert(all(jbsupsmnh(m0,:,1) .eq. zero), 'bmnsh(m0:,1) != 0')
418  CALL assert(all(jbsupsmnh(m1+1:,:,1) .eq. zero), 'bmnsh(m2:,1) != 0')
419  CALL assert(all(jbsupumnh(m0:m1,:,1) .eq. jbsupumnh(m0:m1,:,2)), &
420  'bmnuh(1) != bmnuh(2)')
421  CALL assert(all(jbsupumnh(m1+1:,:,1) .eq. zero), 'bmnuh(m2:,1) != 0')
422  CALL assert(all(jbsupvmnh(m0,:,1) .eq. jbsupvmnh(m0,:,2)), &
423  'bmnvh(m0,1) != bmnvh(m0,2)')
424  CALL assert(all(jbsupvmnh(m1:,:,1) .eq. zero), 'bmnvh(m1:,1) != 0')
425  CALL assert(all(jpmnh(m0,:,1) .eq. jpmnh(m0,:,2)), &
426  'pmnh(m0,1) != pmnh(m0,2)')
427  CALL assert(all(jpmnh(m1:,:,1) .eq. zero), 'pmnh(m1:,1) != 0')
428  END IF
429 
430  IF (iparity .EQ. f_sin) THEN
431  fcomb = f_none
432  fours = f_sin
433  fouruvp = f_cos
434  ELSE
435  fcomb = f_sum
436  fours = f_cos
437  fouruvp = f_sin
438  END IF
439 
440  CALL fourier_context%toijsp(jbsupsmnh(:,:,nsmin:nsmax), bsupsijh0, &
441  fcomb, fours)
442  CALL fourier_context%toijsp(jbsupumnh(:,:,nsmin:nsmax), bsupuijh0, &
443  fcomb, fouruvp)
444  CALL fourier_context%toijsp(jbsupvmnh(:,:,nsmin:nsmax), bsupvijh0, &
445  fcomb, fouruvp)
446  CALL fourier_context%toijsp(jpmnh(:,:,nsmin:nsmax), pijh0, fcomb, fouruvp)
447 
448  END SUBROUTINE
449 
450 !-------------------------------------------------------------------------------
457 !-------------------------------------------------------------------------------
458  SUBROUTINE filterpressure(parity)
459  USE island_params, ONLY: fourier_context
460 
461  IMPLICIT NONE
462 
463 ! Declare Arguments
464  INTEGER, INTENT(in) :: parity
465 
466 ! local variables
467  INTEGER :: fcomb
468  INTEGER :: istat
469  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: pmnh
470 
471 ! Start of executable code
472  IF (parity .eq. f_cos) THEN
473  fcomb = f_none
474  ELSE
475  fcomb = f_sum
476  END IF
477 
478 ! pmnh contains Nyquist-filtered (to MPOL, NTOR) Fourier harmonics of p on
479 ! half mesh
480  ALLOCATE(pmnh(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
481  CALL assert_eq(0, istat, 'Allocation error in FilterPressure')
482 
483  CALL fourier_context%tomnsp(pijh0, pmnh, parity)
484 
485  IF (parity .eq. f_cos .and. nsmax .eq. ns) THEN
486  pedge = pmnh(m0,0,ns)*fourier_context%orthonorm(0,0)
487  END IF
488 
489 ! Spectrally filtered angle derivatives of p in real space (half mesh).
490  CALL fourier_context%toijsp(pmnh, pijh0_du, ior(f_du, fcomb), parity)
491  CALL fourier_context%toijsp(pmnh, pijh0_dv, ior(f_dv, fcomb), parity)
492 
493  IF (l_update_state) THEN
494  CALL fourier_context%toijsp(pmnh, pfilter, fcomb, parity)
495  END IF
496 
497  DEALLOCATE(pmnh)
498 
499  END SUBROUTINE
500 
501 !-------------------------------------------------------------------------------
509 !-------------------------------------------------------------------------------
510  SUBROUTINE checkfft(jbsupsmnh, jbsupumnh, jbsupvmnh, jpmnh, iparity)
511 
512  IMPLICIT NONE
513 
514 ! Declare Arguments
515  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupsmnh
516  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupumnh
517  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupvmnh
518  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jpmnh
519  INTEGER, INTENT(in) :: iparity
520 
521 ! Local Variables
522  INTEGER :: fours
523  INTEGER :: fouruv
524  INTEGER :: fcomb
525  INTEGER :: sparity
526  INTEGER :: m
527  INTEGER :: mp
528  INTEGER :: n
529  INTEGER :: np
530  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: FFT_TEST
531  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: pmnh
532  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: temp
533  REAL (dp) :: rtest
534 
535 ! Local Parameters
536  REAL (dp), PARAMETER :: rtol = 1.e-14_dp
537 
538 ! Start of executable code.
539 #undef _TEST_INIT
540 !#define _TEST_INIT
541 #if defined(_TEST_INIT)
542  IF (iparity .EQ. f_sin) THEN
543  fours = f_sin
544  fouruv = f_cos
545  sparity = 1
546  ELSE
547  fours = f_cos
548  fouruv = f_sin
549  sparity = -1
550  END IF
551 
552 !check independent variables
553  ALLOCATE (fft_test(0:mpol,-ntor:ntor,nsmin:nsmax), &
554  temp(ntheta,nzeta,nsmin:nsmax), &
555  pmnh(0:mpol,-ntor:ntor,nsmin:nsmax), stat=m)
556  CALL assert(m.EQ.0,' Allocation failed in CheckFFT')
557  temp = bsupsijh0*jacobh(:,:,nsmin:nsmax)
558  CALL tomnsp(temp, fft_test, fours)
559  rtest = maxval(abs(fft_test - jbsupsmnh(:,:,nsmin:nsmax)))
560  CALL assert(rtest.LT.rtol, ' jbsupsmnh FFT error in siesta_init')
561  temp = bsupuijh0*jacobh(:,:,nsmin:nsmax)
562  CALL tomnsp(temp, fft_test, fouruv)
563  rtest = maxval(abs(fft_test - jbsupumnh(:,:,nsmin:nsmax)))
564  CALL assert(rtest.LT.rtol,' jbsupumnh FFT error in siesta_init')
565  temp = bsupvijh0*jacobh(:,:,nsmin:nsmax)
566  CALL tomnsp(temp, fft_test, fouruv)
567  rtest = maxval(abs(fft_test - jbsupvmnh(:,:,nsmin:nsmax)))
568  CALL assert(rtest.LT.rtol,' jbsupvmnh FFT error in siesta_init')
569  temp = pijh0*jacobh(:,:,nsmin:nsmax)
570  CALL tomnsp(temp, fft_test, fouruv)
571  rtest = maxval(abs(fft_test - jpmnh(:,:,nsmin:nsmax)))
572  CALL assert(rtest.LT.rtol,' jpmnh FFT error in siesta_init')
573 
574 ! check angle derivatives of pressure
575  CALL tomnsp(pijh0, pmnh, fouruv)
576  CALL tomnsp(pijh0_du, fft_test, fouruv)
577  DO m = 0, mpol
578  mp = -sparity*m
579  CALL assert(maxval(abs(fft_test(m,:,:) - mp*pmnh(m,:,:))) .LE. rtol, &
580  ' pijh0_du FFT error in siesta_init')
581  END DO
582 
583  CALL tomnsp(pijh0_dv, fft_test, fouruv)
584  DO n = -ntor, ntor
585  np = -sparity*n*nfp
586  CALL assert(maxval(abs(fft_test(:,n,:) - np*pmnh(:,n,:))) .LE. rtol, &
587  ' pijh0_dv FFT error in siesta_init')
588  END DO
589 
590  DEALLOCATE(fft_test, pmnh, temp)
591 #endif
592 
593  END SUBROUTINE
594 
595 !-------------------------------------------------------------------------------
597 !-------------------------------------------------------------------------------
598  SUBROUTINE checkpressure
599 
600  IMPLICIT NONE
601 
602 ! Start of executable code
603  IF (lverbose .and. minval(pijh0) .lt. zero) THEN
604  WRITE (*,1000)
605  END IF
606 
607 ! This may stabilize convergence (lower f2 exponent or raise lev_marq)
608  WHERE (pijh0 .LT. zero)
609  pijh0 = zero
610  END WHERE
611 
612 1000 FORMAT(' pijh < 0 in init_state. ===> Recommend lowering f2 exponent', &
613  ' in evolution')
614 
615  END SUBROUTINE
616 
617  END MODULE siesta_init
shared_data::l_getfsq
logical l_getfsq
Compute |F|^2.
Definition: shared_data.f90:216
shared_data::l_par_state
logical l_par_state
Parallel allocated quantities? FIXME: check this.
Definition: shared_data.f90:226
shared_data::l_pedge
logical, parameter l_pedge
Preserve s=1 as iso-pressure surface.
Definition: shared_data.f90:43
siesta_init::init_allocate_arrays
subroutine init_allocate_arrays(lpar)
subroutine for allocating unperturbed array
Definition: siesta_init.f90:199
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
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
siesta_init::init_state
subroutine, public init_state(lcurrent_only, lpar_in)
Initialize equilibrium state.
Definition: siesta_init.f90:45
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::buv_res
real(dp), dimension(:,:,:), allocatable buv_res
Resonant magnetic field perturbation.
Definition: shared_data.f90:104
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
shared_data::l_update_state
logical l_update_state
Update the ste array.
Definition: shared_data.f90:224
metrics::tolowerf
subroutine tolowerf(xsupsij, xsupuij, xsupvij,
Converts to covariant component full grid.
Definition: metrics.f90:590
fourier::f_du
integer, parameter f_du
Poloidal derivative flag.
Definition: fourier.f90:36
shared_data::ste
real(dp), dimension(4) ste
Spectral Truncation RMS error,.
Definition: shared_data.f90:133
siesta_init::getfields
subroutine getfields(jbsupsmnh, jbsupumnh, jbsupvmnh, jpmnh, ipari
Get magnetic field and pressure on the real space mesh.
Definition: siesta_init.f90:397
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
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_init
Initializes unperturbed siesta fields and pressure in real space.
Definition: siesta_init.f90:10
shared_data
This file contains variables and parameters used by many modules in SIESTA.
Definition: shared_data.f90:10
fourier::f_dv
integer, parameter f_dv
Toroidal derivative flag.
Definition: fourier.f90:38
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