V3FIT
runvmec.f
1  SUBROUTINE runvmec(ictrl_array, input_file0,
2  & lscreen, COMM_WORLD, reset_file_name)
3  USE vmec_main
4  USE vmec_params, ONLY: bad_jacobian_flag, more_iter_flag,
5  & norm_term_flag, successful_term_flag,
6  & restart_flag, readin_flag,
7  & timestep_flag, ns_error_flag,
8  & reset_jacdt_flag, lamscale
9  USE realspace
10  USE vmec_params, ONLY: ntmax
11  USE vacmod, ONLY: nuv, nuv3
12  USE timer_sub
13  USE parallel_include_module
14  USE parallel_vmec_module, ONLY: myenvvariables
15  USE parallel_vmec_module, ONLY: initrunvmec
16  USE parallel_vmec_module, ONLY: finalizerunvmec
17  USE parallel_vmec_module, ONLY: initsurfacecomm
18  USE parallel_vmec_module, ONLY: finalizesurfacecomm
19  USE parallel_vmec_module, ONLY: setvacuumcommunicator
20  USE blocktridiagonalsolver_bst, ONLY: initialize_bst
21  USE blocktridiagonalsolver_bst, ONLY: finalize_bst
22  USE xstuff
23  USE mpi_inc
24  IMPLICIT NONE
25 C-----------------------------------------------
26 C D u m m y A r g u m e n t s
27 C-----------------------------------------------
28  INTEGER, INTENT(inout), TARGET :: ictrl_array(5)
29  LOGICAL, INTENT(in) :: lscreen
30  CHARACTER(LEN=*), INTENT(in) :: input_file0
31  CHARACTER(LEN=*), OPTIONAL :: reset_file_name
32  INTEGER, INTENT(IN), OPTIONAL :: COMM_WORLD
33 C-----------------------------------------------
34 C L o c a l V a r i a b l e s
35 C-----------------------------------------------
36  INTEGER, POINTER :: ier_flag
37  INTEGER :: ictrl_flag, iseq_count
38  INTEGER :: ns_index, ns_min, nsval,
39  & ns_old=0, numsteps
40  INTEGER :: igrid, index_end, index_dat,
41  & jacob_off, niter_store
42  INTEGER, SAVE :: igrid0
43  INTEGER :: max_grid_size, flag
44  CHARACTER(LEN=120) :: input_file
45  LOGICAL :: lreset
46  REAL(dp) :: rvton, rvtoff, tiniton, tinitoff
47  REAL(dp) :: gridton, gridtoff
48  REAL(dp) :: bcastton, bcasttoff
49  REAL(dp), ALLOCATABLE, DIMENSION(:) :: bcastarr
50  INTEGER :: blklength, grid_id, i, js,
51  & nsmin, nsmax
52  CHARACTER(LEN=20) :: fname
53 
54 C-----------------------------------------------
55 !
56 ! ictrl_flag = ictrl_array(1)
57 ! flag that controls calling of various subroutines of vmec code
58 ! add together the values beow to utilize several subroutines with one call
59 !
60 ! value flag-name calls routines to...
61 ! ----- --------- ---------------------
62 ! 1 restart_flag reset internal run-control parameters (for example, if
63 ! jacobian was bad, to try a smaller time-step)
64 ! 2 readin_flag read in data from input_file and initialize parameters/arrays
65 ! which do not dependent on radial grid size
66 ! allocate internal grid-dependent arrays used by vmec;
67 ! initialize internal grid-dependent vmec profiles (xc, iota, etc);
68 ! setup loop for radial multi-grid meshes or, if ns_index = ictrl_array(4)
69 ! is > 0, use radial grid points specified by ns_array[ns_index]
70 ! 4 timestep_flag iterate vmec either by "niter" time steps or until ftol satisfied,
71 ! whichever comes first. If numsteps (see below) > 0, vmec will return
72 ! to caller after numsteps, rather than niter, steps.
73 ! 8 output_flag write out output files (wout, jxbout)
74 ! 16 cleanup_flag cleanup (deallocate arrays) - this terminates present run of the sequence
75 ! This flag will be ignored if the run might be continued. For example,
76 ! if ier_flag (see below) returns the value more_iter_flag, the cleanup
77 ! code will be skipped even if cleanup_flag is set, so that the run
78 ! could be continued on the next call to runvmec.
79 ! 32 reset_jacdt_flag Resets ijacobian flag and time step to delt0
80 !
81 ! thus, setting ictrl_flag = 1+2+4+8+16 will perform ALL the tasks thru cleanup_flag
82 ! in addition, if ns_index = 0 and numsteps = 0 (see below), vmec will control its own run history
83 !
84 ! ier_flag = ictrl_array(2)
85 ! specifies vmec error condition (if nonzero)
86 ! numsteps = ictrl_array(3)
87 ! number time steps to evolve the equilibrium. Iterations will stop EITHER if numsteps > 0 and
88 ! when the number of vmec iterations exceeds numsteps; OR if the ftol condition is satisfied,
89 ! whichever comes first. The timestep_flag must be set (in ictrl_flag) for this to be in effect.
90 ! If numsteps <= 0, then vmec will choose consecutive (and increasing) values from the ns_array
91 ! until ftol is satisfied on each successive multi-grid.
92 ! ns_index = ictrl_array(4)
93 ! if > 0 on entry, specifies index (in ns_array) of the radial grid to be used for the present iteration
94 ! phase. If ns_index <= 0, vmec will use the previous value of this index (if the ftol
95 ! condition was not satisfied during the last call to runvmec) or the next value of this index,
96 ! and it will iterate through each successive non-zero member of the ns_array until ftol-convergence
97 ! occurs on each multigrid.
98 ! on exit, contains last value of ns_array index used
99 ! iseq_count=ictrl_array(5)
100 ! specifies a unique sequence label for identifying output files in a sequential vmec run
101 C-----------------------------------------------
102  INTERFACE
103  SUBROUTINE initialize_radial(nsval, ns_old, delt0,
104  & lscreen, reset_file_name)
105  USE vmec_main
106  IMPLICIT NONE
107  INTEGER, INTENT(in) :: nsval
108  INTEGER, INTENT(inout) :: ns_old
109  CHARACTER(LEN=*), OPTIONAL :: reset_file_name
110  LOGICAL, INTENT(in) :: lscreen
111  REAL(rprec), INTENT(out) :: delt0
112  END SUBROUTINE initialize_radial
113  END INTERFACE
114 
115  runvmec_pass = runvmec_pass + 1
116  CALL second0(rvton)
117  CALL myenvvariables
118  CALL initrunvmec(comm_world,lfreeb)
119  lv3fitcall = l_v3fit
120  IF (lv3fitcall) THEN
121  IF (runvmec_pass.GT.1) THEN
122  CALL serial2parallel4x(xc,pxc)
123  CALL serial2parallel4x(xcdot,pxcdot)
124  CALL serial2parallel4x(xstore,pxstore)
125  CALL second0(bcastton)
126  CALL mpi_bcast(pxc,SIZE(pxc), mpi_real8, 0,
127  & runvmec_comm_world, mpi_err)
128  CALL mpi_bcast(pxcdot, SIZE(pxcdot), mpi_real8, 0,
129  & runvmec_comm_world, mpi_err)
130  CALL mpi_bcast(pxstore, SIZE(pxstore), mpi_real8, 0,
131  & runvmec_comm_world, mpi_err)
132  CALL mpi_bcast(iotas, SIZE(iotas), mpi_real8, 0,
133  & runvmec_comm_world, mpi_err)
134  CALL mpi_bcast(iotaf, SIZE(iotaf), mpi_real8, 0,
135  & runvmec_comm_world, mpi_err)
136  CALL mpi_bcast(phips, SIZE(phips), mpi_real8, 0,
137  & runvmec_comm_world, mpi_err)
138  CALL mpi_bcast(phipf, SIZE(phipf), mpi_real8, 0,
139  & runvmec_comm_world, mpi_err)
140  CALL mpi_bcast(chips, SIZE(chips), mpi_real8, 0,
141  & runvmec_comm_world, mpi_err)
142  CALL mpi_bcast(chipf, SIZE(chipf), mpi_real8, 0,
143  & runvmec_comm_world, mpi_err)
144  CALL mpi_bcast(mass, SIZE(mass), mpi_real8, 0,
145  & runvmec_comm_world, mpi_err)
146  CALL mpi_bcast(icurv, SIZE(icurv), mpi_real8, 0,
147  & runvmec_comm_world, mpi_err)
148  CALL mpi_bcast(lamscale, 1, mpi_real8, 0,
149  & runvmec_comm_world, mpi_err)
150 
151  nsmin = t1lglob; nsmax = t1rglob
152  DO js = nsmin, nsmax
153  pphip(:,js) = phips(js)
154  pchip(:,js) = chips(js)
155  END DO
156 
157  CALL second0(bcasttoff)
158  broadcast_time = broadcast_time + (bcasttoff - bcastton)
159  END IF
160  END IF
161 
162  ictrl_flag = ictrl_array(1)
163  numsteps = ictrl_array(3)
164  ier_flag => ictrl_array(2)
165  ns_index = ictrl_array(4)
166  iseq_count = ictrl_array(5)
167  CALL second0(timeon)
168 
169 !
170 ! PARSE input_file into path/input.ext
171 !
172  index_dat = index(input_file0, 'input.')
173  index_end = len_trim(input_file0)
174  IF (index_dat .gt. 0) THEN
175  input_file = trim(input_file0)
176  input_extension = input_file0(index_dat+6:index_end)
177  ELSE
178  input_extension = input_file0(1:index_end)
179  input_file = 'input.'//trim(input_extension)
180  END IF
181 
182 !
183 ! INITIALIZE PARAMETERS
184 !
185  lreset = (iand(ictrl_flag, restart_flag) .ne. 0)
186 
187  IF (lreset) THEN
188  CALL reset_params
189 ! res0 = -1 Done in reset_params
190  END IF
191 
192  IF (iand(ictrl_flag, reset_jacdt_flag) .NE. 0) THEN
193  ijacob = 0
194  delt0r = delt
195  END IF
196 
197  IF (iand(ictrl_flag, readin_flag) .NE. 0) THEN
198 !
199 ! READ INPUT FILE (INDATA NAMELIST), MGRID_FILE (VACUUM FIELD DATA)
200 !
201  CALL vsetup (iseq_count)
202 
203  CALL readin (input_file, iseq_count, ier_flag, lscreen)
204  max_grid_size = ns_array(multi_ns_grid)
205 
206  IF (ier_flag .NE. 0) GOTO 1000
207 !
208 ! COMPUTE NS-INVARIANT ARRAYS
209 !
210  CALL fixaray
211  END IF
212 
213 ! IF(lfreeb) CALL SetVacuumCommunicator(nuv, nuv3, max_grid_size)
214 
215  IF (lreset) THEN
216 !
217 ! COMPUTE INITIAL SOLUTION ON COARSE GRID
218 ! IF PREVIOUS SEQUENCE DID NOT CONVERGE WELL
219 !
220 ! IF (lreseta) THEN !NOTE: where externally, lreseta = T, set restart_flag bit
221 ! (ictrl_flag = IOR(ictrl_flag,restart_flag))
222  igrid0 = 1
223  ns_old = 0
224  IF (PRESENT(reset_file_name) .AND.
225  & len_trim(reset_file_name) .ne. 0) THEN
226  igrid0 = multi_ns_grid
227  END IF
228  IF (grank .EQ. 0) WRITE (nthreed, 30)
229  delt0r = delt
230  END IF
231 
232  30 FORMAT(' FSQR, FSQZ = Normalized Physical Force Residuals',/,
233  & ' fsqr, fsqz = Preconditioned Force Residuals',/,
234  & 1x,23('-'),/, ' BEGIN FORCE ITERATIONS',/,1x,23('-'),/)
235 
236  IF (all(ns_array .eq. 0) .and. ns_index .le. 0) THEN
237  ier_flag = ns_error_flag
238  GOTO 1000
239  END IF
240 
241  jacob_off = 0
242 
243  IF (iand(ictrl_flag, timestep_flag) .EQ. 0) GOTO 1000
244 
245  IF(lfreeb) CALL setvacuumcommunicator(nuv, nuv3, max_grid_size) !SAL 070719
246 
247  50 CONTINUE
248  iequi = 0
249  IF (lfreeb .and. jacob_off .eq. 1) ivac = 1 !!restart vacuum calculations
250 
251  ns_min = 3
252 
253  num_grids = multi_ns_grid
254  IF(.NOT.ALLOCATED(grid_procs)) THEN
255  ALLOCATE(grid_procs(num_grids))
256  ALLOCATE(grid_size(num_grids))
257  ALLOCATE(grid_time(num_grids))
258  ALLOCATE(f3d_time(num_grids))
259  ALLOCATE(f3d_num(num_grids))
260  IF (lfreeb) ALLOCATE(vgrid_time(num_grids))
261  END IF
262 
263  f3d_time = 0; f3d_num=0
264  blklength = (ntor + 1)*(mpol1 + 1)
265  !BEGIN - Main loop that will be parallelized - SKS
266  grid_id = 1
267  old_vacuum_time = 0
268 
269  DO igrid = igrid0 - jacob_off, multi_ns_grid
270  CALL second0(gridton)
271 
272  IF (igrid .lt. igrid0) THEN
273 ! TRY TO GET NON-SINGULAR JACOBIAN ON A 3 PT RADIAL MESH
274  nsval = 3; ivac = -1
275  ftolv = 1.e-4_dp
276  ELSE IF (ns_index .gt. 0) THEN
277  IF (ns_index .gt. SIZE(ns_array)) THEN
278  ier_flag = ns_error_flag
279  RETURN
280  END IF
281  nsval = ns_array(ns_index)
282  IF (nsval .le. 0) stop 'NSVAL <= 0: WRONG INDEX VALUE'
283  ftolv = ftol_array(ns_index)
284  niter = niter_array(ns_index)
285  ELSE
286  nsval = ns_array(igrid)
287  IF (nsval .lt. ns_min) cycle
288  ns_min = nsval
289  ictrl_array(4) = igrid
290  ftolv = ftol_array(igrid)
291  niter = niter_array(igrid)
292  END IF
293 
294  CALL second0(tiniton)
295  IF (parvmec .AND. ns_resltn .GE. 1) THEN
296  IF (lactive) THEN
297  CALL gather4xarray(pscalxc)
298  CALL gather4xarray(pxc)
299  END IF
300  CALL finalizesurfacecomm(ns_comm)
301 
302  CALL second0(bcastton)
303  CALL mpi_bcast(pscalxc, SIZE(pscalxc), mpi_real8, 0,
304  & runvmec_comm_world, mpi_err)
305  CALL mpi_bcast(pxc, SIZE(pxc), mpi_real8, 0,
306  & runvmec_comm_world, mpi_err)
307  CALL second0(bcasttoff)
308  broadcast_time = broadcast_time + (bcasttoff - bcastton)
309  END IF
310 
311  CALL initsurfacecomm(nsval, nzeta, ntheta3, ntmax, ntor, mpol1)
312  CALL second0(tinitoff)
313  init_parallel_time = init_parallel_time + (tinitoff-tiniton)
314 
315  grid_size(grid_id) = nsval
316  grid_procs(grid_id) = nranks
317 
318 ! JDH 2012-06-20. V3FIT fix, inserted with change from VMEC 8.48 -> 8.49
319 ! (Not sure just what in initialize_radial messes up convergence - happens slowly)
320 ! Logical l_v3fit is declared in vmec_input, available via vmec_main
321  IF (l_v3fit .AND. ns_old .ne. nsval) THEN
322  CALL initialize_radial(nsval, ns_old, delt0r, lscreen,
323  & reset_file_name)
324  ELSE IF (ns_old .le. nsval) THEN
325  CALL initialize_radial(nsval, ns_old, delt0r, lscreen,
326  & reset_file_name)
327  END IF
328 
329  CALL initialize_bst(.false., nsval, blklength)
330 
331 ! CONTROL NUMBER OF STEPS
332  IF (numsteps .GT. 0) THEN
333  niter_store = niter
334  niter = numsteps + iter2 - 1
335  END IF
336 
337  CALL eqsolve (ier_flag, lscreen)
338 
339  IF (numsteps .GT. 0) THEN
340  niter = niter_store
341  END IF
342 
343  IF (ier_flag .ne. norm_term_flag .and.
344  & ier_flag .ne. successful_term_flag .and.
345  & ier_flag .ne. more_iter_flag) EXIT
346  IF (numsteps .GT. 0 .or. ns_index .GT. 0) EXIT
347 
348 !
349 ! give up if it refuses to converge, M.Drevlak
350 ! it may help to end a vmec run in an optimization environment, if it
351 ! fails to converge in the first iterations of an ns_array sequence
352 ! within the set number of iterations specified by NITER.
353 ! The parameter fgiveup defaults to 30.
354 !
355 
356  IF (lgiveup .and. (fsqr .gt. ftolv*fgiveup .or.
357  & fsqz .gt. ftolv*fgiveup .or.
358  & fsql .gt. ftolv*fgiveup )) THEN
359  print *, "runvmec: giving up due to poor convergence"
360  EXIT
361  END IF
362 
363  CALL finalize_bst(.false.)
364  CALL second0(gridtoff)
365  grid_time(grid_id) = gridtoff - gridton
366  IF (lfreeb) THEN
367  IF (parvmec) THEN
368  vgrid_time(grid_id) = vacuum_time - old_vacuum_time
369  old_vacuum_time = vacuum_time
370  ELSE
371  vgrid_time(grid_id) = s_vacuum_time - old_vacuum_time
372  old_vacuum_time = s_vacuum_time
373  END IF
374  END IF
375  grid_id = grid_id + 1
376  END DO
377  !END - Main loop that will be parallelized - SKS
378 
379  100 CONTINUE
380 
381  IF (ier_flag .eq. bad_jacobian_flag .and. jacob_off .eq. 0) THEN
382  jacob_off = 1
383  GO TO 50
384  END IF
385 
386  CALL second0 (timeoff)
387  timer(tsum) = timer(tsum) + timeoff - timeon
388 !
389 ! WRITE OUTPUT TO THREED1, WOUT FILES; FREE MEMORY ALLOCATED GLOBALLY
390 !
391  1000 IF (lmoreiter .AND.
392  & ier_flag .EQ. more_iter_flag .AND.
393  & grank .EQ. 0) THEN ! J Geiger
394  print *, 'runvmec: Running some more iterations',
395  & ' -> Skipping call to fileout!'
396  ELSE IF (ier_flag .NE. more_iter_flag) THEN
397  IF (parvmec) THEN
398  CALL fileout_par(iseq_count, ictrl_flag, ier_flag, lscreen)
399  ELSE
400  CALL fileout(iseq_count, ictrl_flag, ier_flag, lscreen)
401  END IF
402  END IF
403 
404  IF(lv3fitcall) CALL finalizerunvmec(runvmec_comm_world)
405  CALL second0(rvtoff)
406  runvmec_time = runvmec_time + (rvtoff - rvton)
407 
408  END SUBROUTINE runvmec
409 !------------------------------------------------
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
Generated on Thu Mar 5 2020 15:49:27 for V3FIT by   doxygen 1.8.17