V3FIT
siesta_run.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 !
10 !*******************************************************************************
11 
12  MODULE siesta_run
13  USE stel_kinds
14  USE shared_data, ONLY: unit_out
15  USE timer_mod
16 
17  IMPLICIT NONE
18 
19 !*******************************************************************************
20 ! DERIVED-TYPE DECLARATIONSf
21 ! 1) siesta_run base class
22 !
23 !*******************************************************************************
24  TYPE siesta_run_class
26  REAL (dp) :: time_on
27  END TYPE
28 
29  CONTAINS
30 
31 !*******************************************************************************
32 ! CONSTRUCTION SUBROUTINES
33 !*******************************************************************************
34 !-------------------------------------------------------------------------------
43 !-------------------------------------------------------------------------------
44  FUNCTION siesta_run_construct(run_comm, verbose)
45  USE nscalingtools, ONLY: myenvvariables, parsolver, siesta_comm
46  USE descriptor_mod, ONLY: lscalapack, inhessian, iam, nprocs, &
47  nprow, npcol, icontxt, icontxt_1xp, &
48  icontxt_px1, icontxt_global, isroot, &
49  myrow, mycol
50 #if defined(MPI_OPT)
51  USE prof_mod, ONLY: profinit
52 #endif
53  USE shared_data, ONLY: nprecon, lverbose
54  USE perturbation, ONLY: init_data
55  USE hessian, ONLY: hesspass
56  USE siesta_state, ONLY: lfirst
57  USE siesta_error
58  USE diagnostics_mod, ONLY: toroidal_flux0
59  USE blocktridiagonalsolver_s, ONLY: initialize, getranks
60  USE siesta_namelist, ONLY: mpolin, ntorin
61  USE metrics, ONLY: set_grid_sizes
62 
63  IMPLICIT NONE
64 
65 ! Declare Arguments
66  TYPE (siesta_run_class), POINTER :: siesta_run_construct
67  INTEGER, INTENT(in) :: run_comm
68  LOGICAL, INTENT(in) :: verbose
69 
70 ! local variables.
71  REAL (dp) :: ton
72  REAL (dp) :: toff
73 
74 ! Start of executable code
75  ALLOCATE(siesta_run_construct)
76 
77  CALL siesta_error_clear_all
78 
79  lverbose = verbose
80  lfirst = .true.
81  toroidal_flux0=0
82 
83  siesta_comm = run_comm
84  hesspass = 0
85  parsolver=.false.
86 
87 #if defined(MPI_OPT)
88  lscalapack = .true.
89  inhessian = .false.
90 
91  CALL myenvvariables
92  IF (parsolver) THEN
93  lscalapack=.false.
94  CALL initialize(0, 0) !Must call this to initialize mpi timer
95  CALL getranks(iam, nprocs)
96  CALL sks_timers
97  END IF
98 
99  IF (lscalapack) THEN
100  CALL profinit()
101 
102 ! setup blacs
103  CALL blacs_pinfo(iam, nprocs)
104  CALL blacs_setup(iam, nprocs)
105  DO nprow = int(sqrt(dble(nprocs))) + 1, 1, -1
106  npcol = int( nprocs/nprow )
107  IF (nprow*npcol .eq. nprocs) EXIT
108  END DO
109 
110  CALL blacs_get(0, 0, icontxt)
111  CALL blacs_gridinit(icontxt, 'C', nprow, npcol)
112 
113  CALL blacs_get(0, 0, icontxt_1xp)
114  CALL blacs_gridinit(icontxt_1xp, 'R', 1, nprow*npcol)
115 
116  CALL blacs_get(0, 0, icontxt_px1)
117  CALL blacs_gridinit(icontxt_px1, 'C', nprow*npcol, 1)
118 
119  icontxt_global = icontxt_1xp
120 
121  CALL blacs_gridinfo(icontxt, nprow, npcol, myrow, mycol)
122  isroot = (myrow .eq. 0) .and. (mycol .eq. 0)
123  IF (isroot .and. lverbose) THEN
124  WRITE (*,1000) nprow, npcol
125  WRITE (*,1001) icontxt_global, icontxt, icontxt_1xp, icontxt_px1
126  END IF
127  END IF
128 #else
129  lscalapack = .false.
130  iam = 0
131  nprocs = 1
132 #endif
133 
134  CALL second0(ton)
135  siesta_run_construct%time_on = ton
136 
137  nprecon = init_data()
138  CALL set_grid_sizes(mpolin, ntorin)
139 
140  CALL second0(toff)
141 
142  init_data_time = toff - ton
143 
144 1000 FORMAT('nprow,npcol ',2(x,i5))
145 1001 FORMAT('icontxt_global,icontxt,icontxt_1xp,icontxt_px1 ',4(x,i5))
146 
147  END FUNCTION siesta_run_construct
148 
149 !*******************************************************************************
150 ! DESTRUCTION SUBROUTINES
151 !*******************************************************************************
152 !-------------------------------------------------------------------------------
161 !-------------------------------------------------------------------------------
162  SUBROUTINE siesta_run_destruct(this, finalize_mpi, close_wout)
163  USE nscalingtools, ONLY: output_timings, gettimes, parsolver, &
164  parfunctisl, pargmres, sksdbg, tofu, &
165  finalizeremap
166  USE blocktridiagonalsolver_s, ONLY: finalize
167  USE descriptor_mod, ONLY: iam, icontxt, lscalapack, myrow, mycol, nprocs
168  USE shared_data, ONLY: nprecon, xc, gc, fsq_total1, mblk_size, &
169  niter, col_scale, lverbose, ngmres_type, &
170  lcolscale, gc_sup
171  USE hessian, ONLY: levmarq_param, mupar, asym_index, dealloc_hessian
172  USE quantities, ONLY: fbdy, dealloc_quantities
173  USE dumping_mod, ONLY: write_output
174  USE diagnostics_mod, ONLY: write_profiles, dealloc_diagnostics
175  USE siesta_namelist, ONLY: wout_file
176  USE island_params, ONLY: ns_i
177  USE metrics, ONLY: cleanup_metric_elements, dealloc_full_lower_metrics
178 #if defined(MPI_OPT)
179  USE prof_mod, ONLY: profstat
180 #endif
181  IMPLICIT NONE
182 
183 ! Declare Arguments
184  TYPE (siesta_run_class), POINTER :: this
185  LOGICAL, INTENT(in) :: finalize_mpi
186  LOGICAL, INTENT(in) :: close_wout
187 
188 ! local variables.
189  REAL (dp) :: ton
190  REAL (dp) :: toff
191  INTEGER :: i
192 
193 ! Start of executable code
194  CALL second0(toff)
195  time_total = toff - this%time_on
196 
197  IF (output_timings) THEN
198  CALL gettimes
199  END IF
200 
201  IF (iam .EQ. 0) THEN
202  DO i = 6, unit_out, unit_out-6
203  IF (.NOT.lverbose .AND. i.EQ.6) cycle
204  WRITE (i, 1000) nprecon, levmarq_param, mupar, asym_index
205  WRITE (i,'(a,i5)') ' Number processors: ', nprocs
206  WRITE (i, 1001)
207  WRITE (i, 1002) time_total, fbdy(1), &
208  time_init, fbdy(2), &
209  time_diag_prec, fbdy(3), &
210  time_block_prec, fbdy(4), &
211  time_factor_blocks, fbdy(5), &
212  time_toijsp, fbdy(6), &
213  time_tomnsp, &
214  gmres_time, &
215  conj_grad_time, &
216  time_update_pres, fbdy(7), &
217  time_update_bfield, fbdy(8), &
218  time_current, fbdy(9), &
219  get_force_harmonics_time, fbdy(10), &
220  time_update_force, fbdy(11), &
221  time_update_upperv, fbdy(12), &
222  time_update_state, fbdy(13), &
223  time_funci, &
224  time_apply_precon, &
225  (diag_add_pert_time + block_add_pert_time), &
226  time_init_state
227  WRITE (i,*)
228  IF (parsolver) THEN
229  WRITE (i,*) 'PARSOLVER=T : NSCALED'
230  ELSE IF (lscalapack) THEN
231  WRITE (i,*) 'PARSOLVER=F : SCALAPACK'
232  ELSE
233  WRITE (i,*) 'PARSOLVER=F : SERIAL'
234  END IF
235  WRITE (i,*) 'PARFUNCTISL :', parfunctisl
236  WRITE (i,*) 'COLUMN SCALING :', lcolscale
237  WRITE (i,'(1x,a,L2,a,i2)') 'PARGMRES :', pargmres, &
238  ' GMRES_TYPE: ', ngmres_type
239 
240  WRITE (i,*) 'OUTPUT_TIMINGS :', output_timings
241 
242  WRITE (i,*)
243  WRITE (i,101)' TIME DIVB : ', time_divb
244  WRITE (i,101)' TIME DIVJ : ', time_divj
245  WRITE (i,101)' TIME BGRADP : ', time_bgradp
246  WRITE (i,101)' TIME BDOTJ : ', time_bdotj
247  WRITE (i,*)
248  WRITE (i,102)' M (block size) :', mblk_size
249  WRITE (i,102)' N (block rows) :', ns_i
250  WRITE (i,102)' P (processors) :', nprocs
251  END DO
252  ENDIF
253  101 FORMAT(a,1p,e10.3)
254  102 FORMAT(a,i5)
255 
256  CALL write_output(wout_file, niter, close_wout)
257  CALL write_profiles(fsq_total1) ! SPH: write pmn, bsupXmn, ksubXmn, jvsupXmn profiles
258  IF (iam .EQ. 0) THEN
259  IF (lverbose) print *,' Writing output to "siesta_profiles.txt" is finished!'
260  CLOSE (unit=unit_out)
261  ENDIF
262 
263 !CLEAN UP ARRAYS
264  CALL cleanup_metric_elements
265 
266  CALL dealloc_quantities
267  CALL dealloc_hessian
268  CALL dealloc_diagnostics
269  CALL dealloc_full_lower_metrics
270  IF (ALLOCATED(xc)) DEALLOCATE(xc, gc, gc_sup)
271  IF (ALLOCATED(col_scale)) DEALLOCATE(col_scale)
272 
273 #if defined(MPI_OPT)
274  IF (lscalapack) THEN
275  CALL blacs_barrier(icontxt, 'All')
276  CALL blacs_gridexit(icontxt)
277 
278  CALL blacs_exit(0)
279  IF ((myrow.EQ.0) .AND. (mycol.EQ.0)) THEN
280  CALL profstat
281  END IF
282 
283  ELSE
284  IF (sksdbg) THEN
285  WRITE(tofu,*) 'Called finalizeRemap and Finalize'
286  FLUSH(tofu)
287  END IF
288  CALL finalizeremap
289  CALL finalize(finalize_mpi)
290 
291  END IF
292 #endif
293 
294  DEALLOCATE(this)
295 
296 1000 FORMAT(/,' nprecon: ',i3,' LM parameter: ',1pe9.2,' mu||: ',1pe9.2, &
297  ' Symmetry Index: ',1pe9.2)
298 1001 FORMAT(/,'==============================',14x,'=======================', &
299  /, &
300  /,' TIMING INFORMATION ',14x,' RMS BOUNDARY FORCES', &
301  /, &
302  /,'==============================',14x,'=======================')
303 1002 FORMAT(' Total runtime : ', f12.3,15x,'fs(1,m=1) :',1pe10.2/, &
304  ' Initialization : ',0pf12.3,15x,'fs(2,m=1) :',1pe10.2/, &
305  ' Diagonal prec : ',0pf12.3,15x,'fs(2,m!=1) :',1pe10.2/, &
306  ' Compute blocks : ',0pf12.3,15x,'fu(1,m=1) :',1pe10.2/, &
307  ' Factor blocks : ',0pf12.3,15x,'fu(2,m=1) :',1pe10.2/, &
308  ' Toijsp : ',0pf12.3,15x,'fu(2,m!=1) :',1pe10.2/, &
309  ' Tomnsp : ',0pf12.3,/, &
310  ' GMRES : ',0pf12.3,/, &
311  ' Conj Gradient : ',0pf12.3,//, &
312  ' SUBROUTINES ',/, &
313  ' Update Pressure: ',0pf12.3,15x,'fv(1,m=0) :',1pe10.2/, &
314  ' Update Bfield : ',0pf12.3,15x,'fv(2,m=0) :',1pe10.2/, &
315  ' CV Currents : ',0pf12.3,15x,'fv(2,m!=0) :',1pe10.2/, &
316  ' Force Harmonics: ',0pf12.3,15x,'fu(ns) :',1pe10.2/, &
317  ' Update Force : ',0pf12.3,15x,'fu(ns-1) :',1pe10.2/, &
318  ' Update UpperV : ',0pf12.3,15x,'fv(ns) :',1pe10.2/, &
319  ' Update State : ',0pf12.3,15x,'fv(ns-1) :',1pe10.2/, &
320  ' Funct Island : ',0pf12.3/, &
321  ' Apply Precon : ',0pf12.3/, &
322  ' Add Perturb : ',0pf12.3/, &
323  ' Init State : ',0pf12.3, &
324  /,'==============================',14x,'=======================')
325  END SUBROUTINE siesta_run_destruct
326 
327 !*******************************************************************************
328 ! SETTER SUBROUTINES
329 !*******************************************************************************
330 !-------------------------------------------------------------------------------
337 !-------------------------------------------------------------------------------
338  SUBROUTINE siesta_run_set_vmec(this)
339  USE siesta_namelist, ONLY: nsin, mpolin, ntorin, nfpin, wout_file, &
340  l_vessel
341  USE metrics, ONLY: init_metric_elements, loadgrid, sqrtg
342  USE quantities, ONLY: init_quantities, init_fields
343  USE island_params, ns=>ns_i, mpol=>mpol_i, ntor=>ntor_i
344  USE shared_data, ONLY: nsp
345  USE vmec_info
346  USE grid_extension, ONLY: grid_extender, read_vessel_file
347  USE pchelms
348  USE evolution, ONLY: init_evolution
349  USE nscalingtools, ONLY: initremap
350 
351  IMPLICIT NONE
352 
353 ! Declare Arguments
354  TYPE (siesta_run_class), INTENT(inout) :: this
355 
356 ! local variables
357  INTEGER :: istat
358 
359 ! Start of executable code
360  CALL vmec_info_set_wout(wout_file, nsin, mpolin, ntorin, nfpin)
361 
362 ! CONSTRUCT R, Z, L REAL-SPACE ARRAYS ON SQRT(FLUX) - "POLAR" - MESH AND
363 ! COMPUTE METRIC ELEMENTS AND JACOBIAN
364  CALL loadgrid(istat)
365  CALL assert_eq(0, istat, 'LoadRZL error in siesta_run_set_vmec')
366 
367 ! If the l_vessel is true try to read a vessel file and extend the grid.
368 ! Otherwise proceed with fixed boundary equilibrium.
369  IF (l_vessel .and. read_vessel_file() .eq. 0) THEN
370  CALL grid_extender(wout_file, 'quad', istat)
371  ELSE
372  l_vessel = .false.
373  END IF
374 
375  CALL init_metric_elements()
376  CALL init_quantities !Initializes BCYCLIC
377  CALL init_evolution !neqs is set here
378  CALL initremap(mpol, ntor, ns, nprocs, iam)
379  CALL inithess
380 
381  IF (l_vessel) THEN
382  CALL run_pchelms
383  END IF
384  CALL init_fields
385 
386  DEALLOCATE (sqrtg)
387 
388  END SUBROUTINE
389 
390 !-------------------------------------------------------------------------------
398 !-------------------------------------------------------------------------------
399  SUBROUTINE siesta_run_set_restart(this)
400  USE restart_mod, ONLY: restart_read
401  USE siesta_namelist, ONLY: lrestart, mpolin, ntorin, nsin, restart_ext, &
402  wout_file, nfpin, nsin_ext, l_vessel
403  USE shared_data, ONLY: nprecon
404  USE metrics, ONLY: init_metric_elements
405  USE quantities, ONLY: init_quantities
406  USE descriptor_mod, ONLY: iam, nprocs
407  USE evolution, ONLY: init_evolution
408  USE nscalingtools, ONLY: initremap
409  USE island_params, ONLY: ns=>ns_i, hs=>hs_i, ohs=>ohs_i, nsh
410  USE hessian, ONLY: inithess
411 
412  IMPLICIT NONE
413 
414 ! Declare Arguments
415  TYPE (siesta_run_class), INTENT(inout) :: this
416 
417 ! Start of executable code
418  IF (l_vessel) THEN
419  ns = nsin + nsin_ext
420  ELSE
421  ns = nsin
422  END IF
423  nsh = ns - 1
424  ohs = nsin - 1
425  hs = 1.0/ohs
426 
427  nprecon = restart_read(restart_ext, wout_file, mpolin, ntorin, nfpin, &
428  ns)
429 
430  CALL init_metric_elements()
431  CALL init_quantities !Initializes BCYCLIC
432  CALL init_evolution !neqs is set here
433  CALL initremap(mpolin, ntorin, ns, nprocs, iam)
434  CALL inithess
435 
436  END SUBROUTINE
437 
438 !*******************************************************************************
439 ! UTILITY SUBROUTINES
440 !*******************************************************************************
441 !-------------------------------------------------------------------------------
451 !-------------------------------------------------------------------------------
452  SUBROUTINE siesta_run_converge(this)
453  USE evolution, ONLY: converge_diagonal, converge_blocks
454  USE descriptor_mod, ONLY: diagonaldone, iam
455  USE siesta_namelist, ONLY: ftol, wout_file
456 
457  IMPLICIT NONE
458 
459 ! Declare Arguments
460  TYPE (siesta_run_class), INTENT(inout) :: this
461 
462 ! Start of executable code
463 ! Converge initial residues with diagonal preconditioner
464  diagonaldone = .false.
465  CALL converge_diagonal(wout_file, ftol)
466  diagonaldone = .true.
467 
468 ! Converge using block preconditioner
469  CALL converge_blocks(wout_file, ftol)
470 
471  END SUBROUTINE
472 
473  END MODULE
siesta_namelist::lrestart
logical lrestart
Restart SIESTA from pervious run.
Definition: siesta_namelist.f90:128
shared_data::nprecon
integer nprecon
Preconditioner flag.
Definition: shared_data.f90:68
shared_data::lverbose
logical lverbose
Use verbose screen output.
Definition: shared_data.f90:234
siesta_namelist::mpolin
integer mpolin
Number of poloidal modes.
Definition: siesta_namelist.f90:166
shared_data::xc
real(dp), dimension(:), allocatable xc
1D array of Fourier mode displacement components.
Definition: shared_data.f90:97
metrics::cleanup_metric_elements
subroutine cleanup_metric_elements
Deallocate memory containing metric elements on the half mesh.
Definition: metrics.f90:642
siesta_namelist::l_vessel
logical l_vessel
If extended grid is to be used using an available vessel file.
Definition: siesta_namelist.f90:140
shared_data::gc_sup
real(dp), dimension(:), allocatable, target gc_sup
1D Array of Fourier mode MHD force components, FIXME Check if this is really needed.
Definition: shared_data.f90:102
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
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
metrics::dealloc_full_lower_metrics
subroutine dealloc_full_lower_metrics
Deallocate memory containing metric elements on the full mesh.
Definition: metrics.f90:666
shared_data::ngmres_type
integer ngmres_type
GMRES control flag.
Definition: shared_data.f90:74
siesta_namelist::nfpin
integer nfpin
Number of field periods to use. -1 means set this to the value in the wout file.
Definition: siesta_namelist.f90:171
siesta_namelist::ntorin
integer ntorin
Number of toroidal modes.
Definition: siesta_namelist.f90:168
siesta_namelist::restart_ext
character(len=siesta_namelist_name_length) restart_ext
Name of the restart file extension.
Definition: siesta_namelist.f90:191
shared_data::col_scale
real(dp), dimension(:,:,:,:), allocatable col_scale
Column scaling factors.
Definition: shared_data.f90:110
shared_data::niter
integer niter
Total number of iteration to run.
Definition: shared_data.f90:60
siesta_namelist::wout_file
character(len=siesta_namelist_name_length) wout_file
Filename of the VMEC woutfile.
Definition: siesta_namelist.f90:189
shared_data::lcolscale
logical lcolscale
Apply column scaling to hessian.
Definition: shared_data.f90:228
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
restart_mod
Contains routines for writting the restart file.
Definition: restart_mod.f90:10
pchelms::run_pchelms
subroutine run_pchelms
Run the pchelms code to solve for inital jbsup values.
Definition: pchelms.f90:77
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
blocktridiagonalsolver_s
Module contains subroutines for solver for block tri-diagonal matrices.
Definition: blocktridiagonalsolver_s.f90:7
siesta_namelist::ftol
real(dp) ftol
Force tolarance.
Definition: siesta_namelist.f90:146
metrics::loadgrid
subroutine loadgrid(istat)
Load the R, Z and lambda arrays from VMEC.
Definition: metrics.f90:218
siesta_namelist::nsin
integer nsin
Radial size of the plasma grid.
Definition: siesta_namelist.f90:162
shared_data::gc
real(dp), dimension(:), allocatable, target gc
1D Array of Fourier mode MHD force components
Definition: shared_data.f90:99
metrics::init_metric_elements
subroutine init_metric_elements()
Initialize the metric elements.
Definition: metrics.f90:103
siesta_error
This module contains all the code needed to define error.
Definition: siesta_error.f90:10
pchelms
This file solves the helmholtz equation to set inital fields that match vmec and vacuum currents from...
Definition: pchelms.f90:12
shared_data::nsp
integer nsp
Total radial grid size in the VMEC region.
Definition: shared_data.f90:64
grid_extension::read_vessel_file
integer function read_vessel_file()
Read the vessel file.
Definition: grid_extension.f90:596
restart_mod::restart_read
integer function restart_read(restart_ext, wout_file, mpolin, ntorin, nfpi
Reads the restart file.
Definition: restart_mod.f90:164
shared_data::fsq_total1
real(dp) fsq_total1
|F|^2 WITHOUT column scaling.
Definition: shared_data.f90:93
grid_extension
This file contains utilities for extending the vmec grid.
Definition: grid_extension.f90:10
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
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
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
siesta_namelist::nsin_ext
integer nsin_ext
Radial size of the extended grid.
Definition: siesta_namelist.f90:164
metrics::set_grid_sizes
subroutine set_grid_sizes(mpol_in, ntor_in)
Set grid sizes.
Definition: metrics.f90:127
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17