V3FIT
siesta_state.f90
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_state
13 #define DUMP_STATE
14  USE stel_kinds
15  USE quantities
16  USE descriptor_mod, ONLY: iam, siesta_comm
17  USE diagnostics_mod
18  USE shared_data, ONLY: ste, bs0, bu0, bsbu_ratio_s, bsbu_ratio_a, &
19  unit_out
20  USE timer_mod, ONLY: time_update_state
21  USE siesta_init, ONLY: init_state
22  USE nscalingtools, ONLY: startglobrow, endglobrow, mpi_err, &
23  parsolver
24  USE mpi_inc
25  USE fourier, ONLY: f_sin, f_cos
26 
27  IMPLICIT NONE
28 
29 !*******************************************************************************
30 ! Module variables
31 !*******************************************************************************
33  LOGICAL, PUBLIC :: lfirst = .true.
34 
35  CONTAINS
36 
37 !*******************************************************************************
38 ! UTILITY SUBROUTINES
39 !*******************************************************************************
40 !-------------------------------------------------------------------------------
49 !-------------------------------------------------------------------------------
50  SUBROUTINE update_state(lprint, fsq_total, ftol)
51 
52  IMPLICIT NONE
53 
54 ! Declare Arguments
55  LOGICAL, INTENT(in) :: lprint
56  REAL (dp), INTENT(in) :: fsq_total
57  REAL (dp), INTENT(in) :: ftol
58 
59 ! local variables
60  REAL (dp) :: ton
61  REAL (dp) :: toff
62  INTEGER :: nsmin
63  INTEGER :: nsmax
64 
65 ! Start of executable code
66  CALL second0(ton)
67 
68  CALL updatefields(jbsupsmnsh, jbsupumnch, jbsupvmnch, jpmnch, &
69  djbsupsmnsh, djbsupumnch, djbsupvmnch, djpmnch)
70  IF (lasym) THEN
71  CALL updatefields(jbsupsmnch, jbsupumnsh, jbsupvmnsh, jpmnsh, &
72  djbsupsmnch, djbsupumnsh, djbsupvmnsh, djpmnsh)
73  END IF
74 
75 ! Gather new total fields onto all processors
76  IF (parsolver) THEN
77  CALL gatherfields
78  END IF
79 
80  CALL second0(toff)
81  time_update_state = time_update_state + (toff - ton)
82 
83 ! Reset perturbations
84  CALL clear_field_perts
85 
86  CALL second0(toff)
87  time_update_state = time_update_state + (toff - ton)
88 
89  IF (.not.lprint) THEN
90  RETURN
91  END IF
92 
93  CALL second0(ton)
94 
95  CALL update_diagnostics(jbsupsmnsh, jbsupumnch, jbsupumnch, jpmnch, &
96  bs0(1:6), bu0(1:6), bsbu_ratio_s, pwr_spec_s, &
97  f_sin)
98  IF (lasym) THEN
99  CALL update_diagnostics(jbsupsmnch, jbsupumnsh, jbsupvmnsh, jpmnsh, &
100  bs0(7:12), bu0(7:12), bsbu_ratio_a, &
101  pwr_spec_a, f_cos)
102  END IF
103  lfirst = .false.
104 
105  nsmin = max(1, startglobrow)
106  nsmax = min(endglobrow, ns)
107 
108  CALL init_state(.true.)
109 
110 ! Compute co-variant and contravariant components of current (times jacobian)
111 ! need for divJ, BdotJ
112  lcurr_init = .true.
113  CALL divb(nsmin, nsmax)
114  CALL divj(nsmin, nsmax)
115  CALL bgradp(nsmin, nsmax)
116  CALL tflux
117  CALL bdotj_par
118  lcurr_init = .false.
119 
120  toroidal_flux = toroidal_flux - toroidal_flux0
121  IF (iam .EQ. 0) THEN
122  IF (lverbose) THEN
123  WRITE (*,1000) ste(1), ste(2), ste(3), ste(4), divb_rms, &
124  toroidal_flux, wp/wb, bgradp_rms, max_bgradp, &
125  min_bgradp, bdotj_rms, bdotj2_rms, divj_rms
126  END IF
127  WRITE (unit_out,1000) ste(1), ste(2), ste(3), ste(4), divb_rms, &
128  toroidal_flux, wp/wb, bgradp_rms, max_bgradp, &
129  min_bgradp, bdotj_rms, bdotj2_rms, divj_rms
130  END IF
131 
132  CALL second0(toff)
133  time_update_state = time_update_state + (toff - ton)
134 
135 1000 FORMAT(' SPECTRAL TRUNC ERROR - p: ',1pe11.3,' B_s: ',1pe11.3, &
136  ' B_u: ',1pe11.3,' B_v: ',1pe11.3,/, &
137  ' DIV-B (rms): ',1pe11.3, ' DEL_TFLUX: ',1pe11.3,/, &
138  ' <BETA>: ', 1pe11.3,' B.GRAD-P (rms): ',1pe11.3, &
139  ' B.GRAD-P (max): ',1pe11.3,' B.GRAD-P (min): ',1pe11.3,/, &
140  ' (J*B)/|JxB| (rms): ', 1pe11.3,' (J_par)/|J_tot| (rms): ', &
141  1pe11.3,' DIV-J (rms): ',1pe11.3)
142 
143  END SUBROUTINE
144 
145 !-------------------------------------------------------------------------------
159 !-------------------------------------------------------------------------------
160  SUBROUTINE updatefields(jbsupsmnh, jbsupumnh, jbsupvmnh, jpmnh, &
161  djbsupsmnh, djbsupumnh, djbsupvmnh, djpmnh)
162 
163  IMPLICIT NONE
164 
165 ! Declare Arguments
166  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: jbsupsmnh
167  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: jbsupumnh
168  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: jbsupvmnh
169  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: jpmnh
170  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: djbsupsmnh
171  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: djbsupumnh
172  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: djbsupvmnh
173  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: djpmnh
174 
175 ! Local Variables
176  INTEGER :: n1
177  INTEGER :: n2
178 
179 ! Start of executable code.
180  n1 = startglobrow
181  n2 = endglobrow
182 
183  jbsupsmnh(:,:,n1:n2) = jbsupsmnh(:,:,n1:n2) + djbsupsmnh(:,:,n1:n2)
184  jbsupumnh(:,:,n1:n2) = jbsupumnh(:,:,n1:n2) + djbsupumnh(:,:,n1:n2)
185  jbsupvmnh(:,:,n1:n2) = jbsupvmnh(:,:,n1:n2) + djbsupvmnh(:,:,n1:n2)
186  jpmnh(:,:,n1:n2) = jpmnh(:,:,n1:n2) + djpmnh(:,:,n1:n2)
187 
188  END SUBROUTINE
189 
190 !-------------------------------------------------------------------------------
202 !-------------------------------------------------------------------------------
203  SUBROUTINE update_diagnostics(jbsupsmnh, jbsupumnh, jbsupvmnh, &
204  jpmnh, bs, bu, bsbu_ratio, &
205  pwr_spec, iparity)
206  USE island_params, ns=>ns_i, hs=>hs_i
207  USE fourier, ONLY: m1
208 
209  IMPLICIT NONE
210 
211 ! Declare Arguments
212  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(IN) :: jbsupsmnh
213  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(IN) :: jbsupumnh
214  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(IN) :: jbsupvmnh
215  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(IN) :: jpmnh
216  REAL (dp), DIMENSION(6), INTENT(out) :: bs
217  REAL (dp), DIMENSION(6), INTENT(out) :: bu
218  REAL (dp), INTENT(out) :: bsbu_ratio
219  REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns,4), INTENT(inout) :: pwr_spec
220  INTEGER, INTENT(in) :: iparity
221 
222 ! Local Variables
223  REAL (dp) :: r12
224  REAL (dp) :: diag_b
225  REAL (dp) :: diag_p
226  REAL (dp), DIMENSION(12) :: d1
227  INTEGER :: js
228  INTEGER :: itype
229  INTEGER :: n1
230  INTEGER :: n2
231 
232 ! Local Parameters
233  CHARACTER (len=4), DIMENSION(2), PARAMETER :: &
234  def = (/'SYM ', 'ASYM'/)
235 
236 ! Start of executable code.
237  n1 = startglobrow
238  n2 = endglobrow
239 
240  IF (iparity .eq. f_sin) THEN
241  itype = 1
242  r12 = hs/2
243  ELSE
244  itype = 2
245  r12 = -hs/2
246  END IF
247 
248 #undef _TEST_STATE
249 !#define _TEST_STATE
250 #if defined(_TEST_STATE)
251  IF (iam .eq. 0) THEN
252  WRITE (*,1000) def(itype)
253  DO js = -ntor,ntor
254  WRITE (*,1001) js, jbsupsmnh(m1,js,2), r12*jbsupumnh(m1,js,2)
255  END DO
256  END IF
257 #endif
258 
259  bs(1) = sqrt(sum((jbsupsmnh(m1,:,2) - r12*jbsupumnh(m1,:,2))**2))/r12
260  bu(1) = sqrt(sum((jbsupsmnh(m1,:,2) + r12*jbsupumnh(m1,:,2))**2))/r12
261 
262 #if defined(MPI_OPT)
263  IF (parsolver) THEN
264  d1(1) = bs(1)
265  d1(2) = bu(1)
266  CALL mpi_bcast(d1, 2, mpi_real8, 0, siesta_comm, mpi_err)
267  bs(1) = d1(1)
268  bu(1) = d1(2)
269  END IF
270 #endif
271  IF (abs(bu(1)) .GT. 1.e-10_dp) THEN
272  bsbu_ratio = bs(1)/bu(1)
273  d1 = 0
274  DO js=n1, min(6, n2)
275  d1(js) = sqrt(sum(jbsupsmnh(m1,:,js)**2))/abs(r12)
276  d1(js+6) = sqrt(sum(jbsupumnh(m1,:,js)**2))
277  END DO
278 #if defined(MPI_OPT)
279  IF (parsolver) THEN
280  CALL mpi_allreduce(mpi_in_place, d1, 12, mpi_real8, mpi_sum, &
281  siesta_comm, mpi_err)
282  END IF
283 #endif
284  bs(1:6) = d1(1:6)
285  bu(1:6) = d1(7:12)
286  ELSE
287  bs(1:6) = 0
288  bu(1:6) = 0
289  bsbu_ratio_s = 1
290  END IF
291 
292 #if defined(DUMP_STATE)
293  IF (lfirst) THEN
294  pwr_spec(:,:,n1:n2,1) = jbsupsmnh(:,:,n1:n2)
295  pwr_spec(:,:,n1:n2,2) = jbsupumnh(:,:,n1:n2)
296  pwr_spec(:,:,n1:n2,3) = jbsupvmnh(:,:,n1:n2)
297  pwr_spec(:,:,n1:n2,4) = jpmnh(:,:,n1:n2)
298  ELSE
299  d1(1) = sum((jbsupsmnh(:,:,n1:n2) - pwr_spec(:,:,n1:n2,1))**2 + &
300  (jbsupumnh(:,:,n1:n2) - pwr_spec(:,:,n1:n2,2))**2 + &
301  (jbsupvmnh(:,:,n1:n2) - pwr_spec(:,:,n1:n2,3))**2)
302  d1(2) = sum((jpmnh(:,:,n1:n2) - pwr_spec(:,:,n1:n2,4))**2)
303 #if defined(MPI_OPT)
304  IF (parsolver) THEN
305  CALL mpi_allreduce(mpi_in_place, d1, 2, mpi_real8, mpi_sum, &
306  siesta_comm, mpi_err)
307  END IF
308 #endif
309  diag_b = d1(1)
310  diag_p = d1(2)
311  IF (iam .EQ. 0) THEN
312  DO js = 6, unit_out, unit_out - 6
313  IF (.NOT.lverbose .AND. js.EQ.6) cycle
314  WRITE(js, 1002) def(itype), diag_b, diag_p
315  END DO
316  END IF
317  END IF
318 
319 1000 FORMAT('ipar: ',a,/,' n B^s r12*B^u')
320 1001 FORMAT(i4,1p,2e14.4)
321 1002 FORMAT(1x,'POWER SPECTRA(',a,') -- dB: ',1p,e10.3,' dP: ',1p,e10.3)
322 #endif
323  END SUBROUTINE
324 
325 !-------------------------------------------------------------------------------
327 !-------------------------------------------------------------------------------
328  SUBROUTINE clear_field_perts
329 
330  IMPLICIT NONE
331 
332 ! Start of executable code.
333  djbsupsmnsh = 0
334  djbsupumnch = 0
335  djbsupvmnch = 0
336  djpmnch = 0
337 
338  IF (lasym) THEN
339  djbsupsmnch = 0
340  djbsupumnsh = 0
341  djbsupvmnsh = 0
342  djpmnsh = 0
343  END IF
344 
345  END SUBROUTINE
346 
347  END MODULE
quantities
This file contains subroutines for allocating and initializing curvilinear magnetic covariant and pre...
Definition: quantities.f90:11
shared_data::bsbu_ratio_s
real(dp) bsbu_ratio_s
FIXME: UNKNOWN.
Definition: shared_data.f90:139
siesta_state
This file contains subroutines for aupdating from t to t + delta_t the magnetic field and pressure as...
Definition: siesta_state.f90:12
shared_data::bu0
real(dp), dimension(12) bu0
FIXME: UNKNOWN.
Definition: shared_data.f90:137
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::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
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
shared_data::unit_out
integer, parameter unit_out
File output io unit.
Definition: shared_data.f90:39
island_params::hs_i
real(dp) hs_i
Radial grid spacing. hs = s_i+1 - s_i.
Definition: island_params.f90:47
shared_data::ste
real(dp), dimension(4) ste
Spectral Truncation RMS error,.
Definition: shared_data.f90:133
shared_data::bs0
real(dp), dimension(12) bs0
FIXME: UNKNOWN.
Definition: shared_data.f90:135
shared_data::bsbu_ratio_a
real(dp) bsbu_ratio_a
FIXME: UNKNOWN.
Definition: shared_data.f90:143
fourier::m1
integer, parameter m1
m = 1 mode.
Definition: fourier.f90:60
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
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17