V3FIT
evolution.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 evolution
13  USE v3_utilities, ONLY: assert
14  USE stel_kinds
15  USE stel_constants
16  USE shared_data
17  USE shared_functions
18  USE island_params, mpol=>mpol_i, ntor=>ntor_i, ns=>ns_i, &
19  nfp=>nfp_i, mnmax=>mnmax_i, hs=>hs_i
20  USE timer_mod
21  USE siesta_namelist, ONLY: eta_factor
22  USE descriptor_mod, ONLY: iam, nprocs
23  USE nscalingtools, ONLY: sksdbg, parsolver, parfunctisl, mpi_err, &
24  startglobrow, endglobrow, rcounts, disp
25  USE mpi_inc
26 
27 
28  IMPLICIT NONE
29 
30 !*******************************************************************************
31 ! Module Parameters
32 !*******************************************************************************
34  INTEGER, PARAMETER :: ndamp = 20
36  INTEGER, PARAMETER :: ncongrad_steps = 100
38  INTEGER, PARAMETER :: nprint_step = 1
39 
40 !*******************************************************************************
41 ! Module variables
42 !*******************************************************************************
44  LOGICAL :: l_EPAR
46  LOGICAL :: l_Gmres
48  LOGICAL :: l_conjmin
49 
51  INTEGER :: nfun_calls
52 
54  REAL (dp) :: delt_cg = 1
56  REAL (dp) :: fsqprev
58  REAL (dp) :: fsqprev1
60  REAL (dp) :: levm_loc
62  REAL (dp) :: fsq_min
64  REAL (dp) :: fsq_ratio
66  REAL (dp) :: fsq_ratio1
68  REAL (dp), DIMENSION(ndamp) :: otau
70  REAL (dp), DIMENSION(:), ALLOCATABLE :: xcdot
71 
72  CONTAINS
73 
74 !*******************************************************************************
75 ! UTILITY SUBROUTINES
76 !*******************************************************************************
77 !-------------------------------------------------------------------------------
86 !-------------------------------------------------------------------------------
87  SUBROUTINE converge_diagonal(wout_file, ftol)
88  USE siesta_namelist, ONLY: ladd_pert, l_output_alliter
89  USE perturbation, ONLY: add_perturb, niter_max
90  USE dumping_mod, ONLY: write_output
91 
92  IMPLICIT NONE
93 
94 ! Declare Arguments
95  CHARACTER (len=*), INTENT(in) :: wout_file
96  REAL(dp), INTENT(IN) :: ftol
97 
98 ! Local Variables
99  LOGICAL :: l_iterate
100  REAL (dp) :: fsq_block
101  REAL (dp) :: t1
102  REAL (dp) :: t2
103 
104 ! Local Parameters
105  REAL (dp),PARAMETER :: fsq_prec = 1.e-10_dp
106 
107 
108 ! Start of executable code.
109  fsq_block = max(ftol, fsq_prec)
110 
111 ! Call evolve until fsq_ratio stops decreasing.
112  l_iterate = .true.
113  l_gmres = .true.
114  nprecon_type = prediag
115 
116  IF (nprecon .ne. 0) THEN
117  nprecon_type = preblock
118  END IF
119  niter = 1
120 
121 ! Need to reset for reconstruction purposes.
122  levm_scale = 1
123  in_hess_nfunct = 0
124  out_hess_nfunct = 0
125 
126 ! Use column-scaled fsq_total here.
127  DO WHILE (l_iterate)
128  CALL second0(t1)
129  CALL evolve
130  CALL second0(t2)
131  diag_evolve_time=diag_evolve_time+(t2-t1)
132 
133  IF (l_output_alliter) THEN
134 ! Do not close wout_file. This session could be run from inside a
135 ! reconstruction instance.
136  CALL write_output(wout_file, niter, .false.)
137  END IF
138 
139  l_iterate = fsq_ratio1 .le. 0.5_dp .and. &
140  fsq_total1 .gt. fsq_block .and. &
141  niter .lt. niter_max
142 
143  IF (ladd_pert .and. fsq_total1.lt.100*fsq_block) THEN
144  l_init_state = .true.
145  CALL second0(t1)
146  CALL add_perturb(xc, getwmhd)
147  CALL second0(t2)
148  diag_add_pert_time=diag_add_pert_time+(t2-t1)
149  ladd_pert = .false.
150  END IF
151  niter = niter + 1
152  END DO
153 
154  END SUBROUTINE
155 
156 !-------------------------------------------------------------------------------
165 !-------------------------------------------------------------------------------
166  SUBROUTINE converge_blocks(wout_file, ftol)
167  USE siesta_namelist, ONLY: l_output_alliter, ladd_pert
168  USE perturbation, ONLY: add_perturb, niter_max
169  USE dumping_mod, ONLY: write_output
170 
171 ! Declare Arguments
172  CHARACTER (len=*), INTENT(in) :: wout_file
173  REAL (dp), INTENT(in) :: ftol
174 
175 ! Local Variables
176  INTEGER :: nrow
177  LOGICAL :: l_iterate
178  REAL (dp) :: t1
179  REAL (dp) :: t2
180 
181 ! Local Parameters
182  REAL (dp),PARAMETER :: fsq_prec = 1.e-10_dp
183 
184 ! Start of executable code.
185  l_gmres = .false.
186  nprecon_type = preblock
187  nprecon = 1
188 ! niter_max = niter + niter_max ! FIXME: MRC REMOVE LATER
189 ! Controls iteration sequence.
190  l_iterate = (niter.LT.niter_max) .AND. (fsq_total1.GT.ftol)
191  nrow = 0
192 
193  DO WHILE (l_iterate)
194  CALL second0(t1)
195  CALL evolve
196  CALL second0(t2)
197  block_evolve_time = block_evolve_time + (t2 - t1)
198 
199  IF (l_output_alliter) THEN
200 ! Do not close wout_file. This session could be run from inside a
201 ! reconstruction instance.
202  CALL write_output (wout_file, niter, .false.)
203  END IF
204 
205  l_iterate = niter .lt. niter_max .and. &
206  fsq_total1 .gt. ftol
207 
208  IF (.not.l_conjmin .and. .not.l_gmres) THEN
209  IF (iam .eq. 0 .and. lverbose) THEN
210  WRITE (*,1000)
211  END IF
212  l_iterate = .false.
213 ! ELSE IF (nprecon .gt. 3 .and. fsq_ratio .gt 1.E5_dp) THEN
214 ! FIXME: Disabled for now so we can test convergence.
215 ! IF (iam .eq. 0 .and. lverbose) THEN
216 ! WRITE(*,1001), fsq_ratio1
217 ! END IF
218 ! l_iterate=.FALSE.
219  END IF
220 
221 ! Stop is lack of progress.
222  IF (nprecon .gt. 2 .and. abs(1 - fsq_ratio1) .lt. 1.e-2_dp) THEN
223  levm_scale = levm_scale/3
224  nrow = nrow + 1
225  IF (nrow.EQ.3 .and. l_iterate) THEN
226 ! IF (iam .eq. 0 .and. lverbose) THEN
227 ! WRITE(*,1001), fsq_ratio1
228 ! END IF
229 ! l_iterate = .false.
230  END IF
231  ELSE
232  nrow = 0
233  END IF
234 
235 ! In case we didn't add it in diag loop.
236  IF (ladd_pert .and. fsq_total1 .lt. 100*fsq_prec) THEN
237  l_init_state = .true.
238  CALL second0(t1)
239  CALL add_perturb(xc, getwmhd)
240  CALL second0(t2)
241  block_add_pert_time = block_add_pert_time + (t2 - t1)
242  ladd_pert = .false.
243 ! To output right after pert is applied, set l_iterate = .false. here.
244  END IF
245 
246  nprecon = nprecon + 1
247  niter = niter + 1
248  END DO
249 
250 1000 FORMAT(' CONJ GRADIENT UNABLE TO MAKE PROGRESS')
251 1001 FORMAT(' FSQ RATIO: ',1p,e10.2, ' SIESTA CAN NOT CONVERGE FURTHER!')
252 
253  END SUBROUTINE
254 
255 !-------------------------------------------------------------------------------
257 !-------------------------------------------------------------------------------
258  SUBROUTINE init_evolution
259  USE hessian, ONLY: levmarq_param0, levmarq_param, l_compute_hessian
260  USE quantities, ONLY: jvsupsmncf,jvsupumnsf,jvsupvmnsf, &
261  jvsupsmnsf,jvsupumncf,jvsupvmncf, &
262  fsubsmncf, fsubumnsf, fsubvmnsf, &
263  fsubsmnsf, fsubumncf, fsubvmncf, &
264  fsupsmncf, fsupumnsf, fsupvmnsf, &
265  fsupsmnsf, fsupumncf, fsupvmncf
266  USE shared_functions, ONLY: init_ptrs
267 
268  IMPLICIT NONE
269 
270 ! Local Variables
271  INTEGER :: istat
272  INTEGER :: n1
273 
274 ! Start of executable code.
275 ! Natural boundary condition at origin. l_push_s and l_push_u can't both be
276 ! true. This causes a linearly dependent => 0 eigenvalue in the Hessian.
277 ! Control evolve displacements (m=1 for s,u; m=0 for v) at s=0.
278  l_push_s = .false.
279  l_push_u = .false.
280  l_push_v = .false.
281 
282 !avoid linear dependency
283  IF (l_push_u) THEN
284  l_push_s = .false.
285  END IF
286 ! l_push_edge = .true.
287  l_push_edge = .false.
288 
289  scale_s = hs
290  scale_u = sqrt(scale_s)
291 
292  ste = 0
293  l_epar = .false.
294 
295  nfun_calls = 0
296  fsqprev = -1; fsqprev1 = -1
297  fsq_min = 1.e20_dp
298  niter = 0
299  wtotal0 = -1
300 ! Initial Newton tolerance parameter.
301  etak_tol = 1.e-01_dp
302  delta_t = 1.0
303  l_linearize = .false.
304  l_getfsq = .true.
305 ! Set in getwmhd function.
306  l_getwmhd = .false.
307  levmarq_param = levmarq_param0
308  fsq_total = 1
309  fsq_total1 = 1
310 
311  CALL assert(ndims .eq. 3 .or. ndims .eq. 6,'WRONG ndims!')
312  CALL assert(mnmax .eq. (1 + mpol)*(2*ntor + 1),'WRONG mnmax!')
313 
314  n1 = ns*mnmax
315  neqs = ndims*n1
316 
317  ALLOCATE(xc(neqs), col_scale(0:mpol,-ntor:ntor,ndims,ns), stat=istat)
318  CALL assert(istat.eq.0,'Allocate xc failed!')
319  xc = 0
320  col_scale = 1
321  CALL init_ptrs(xc, jvsupsmncf, jvsupumnsf, jvsupvmnsf, &
322  jvsupsmnsf, jvsupumncf, jvsupvmncf)
323 
324  ALLOCATE(gc(neqs), gc_sup(neqs), stat=istat)
325  CALL assert(istat.eq.0,'Allocate gc failed!')
326  gc = 0
327  CALL init_ptrs(gc_sup, fsupsmncf, fsupumnsf, fsupvmnsf, &
328  fsupsmnsf, fsupumncf, fsupvmncf)
329  CALL init_ptrs(gc, fsubsmncf, fsubumnsf, fsubvmnsf, &
330  fsubsmnsf, fsubumncf, fsubvmncf)
331 
332  l_compute_hessian = .false.
333 
334  col_scale(:,:,1,:) = scale_s
335  col_scale(:,:,2,:) = scale_u
336  IF (lasym) THEN
337  col_scale(:,:,4,:) = scale_s
338  col_scale(:,:,5,:) = scale_u
339  END IF
340 
341  END SUBROUTINE
342 
343 !-------------------------------------------------------------------------------
349 !-------------------------------------------------------------------------------
350  SUBROUTINE evolve
351  USE quantities, ONLY: jvsupsmncf, jvsupumnsf, jvsupvmnsf, &
352  jvsupsmnsf, jvsupumncf, jvsupvmncf, &
353  fsubsmncf, fsubumnsf, fsubvmnsf, &
354  fsubsmnsf, fsubumncf, fsubvmncf
355  USE gmres, ONLY: gmres_fun
356  USE hessian
357  USE siesta_namelist, ONLY: eta_factor, lresistive, ftol, restart_ext, &
358  wout_file
359  USE perturbation, ONLY: niter_max
360  USE diagnostics_mod, ONLY: bgradp_rms, toroidal_flux0, bgradp, &
361  tflux
362  USE siesta_state, ONLY: update_state
363  USE blocktridiagonalsolver_s, ONLY: refactorhessian
364  USE timer_mod, ONLY: gmres_time
365  USE restart_mod, ONLY: restart_write
366  USE vmec_info, ONLY: vmec_curtor
367 
368  IMPLICIT NONE
369 
370 ! Local Variables.
371  REAL (dp) :: skston
372  REAL (dp) :: skstoff
373  INTEGER :: iprint
374  INTEGER :: m1
375  INTEGER :: i
376  INTEGER :: n1
377  INTEGER :: n2
378  REAL (dp) :: v2tot
379  REAL (dp) :: f1
380  REAL (dp) :: f2
381  REAL (dp) :: t1
382  REAL (dp) :: t2
383  REAL (dp) :: lm0
384  LOGICAL :: lprint
385  LOGICAL :: ldiagonal
386  LOGICAL :: lblock
387  INTEGER, DIMENSION(2) :: tnum_funct
388  CHARACTER(LEN=20) :: str_resistive
389 
390 ! Local Parameters
391  REAL (dp), DIMENSION(1), PARAMETER :: levscan = (/ 0.25_dp /)
392 
393 ! Start of executable code.
394  l_init_state = .true.
395  siesta_curtor = 0.0
396 
397 ! When running under a reconstruction context need to reset this parameter.
398  IF (niter .le. 1) THEN
399  levm_loc = 1
400  END IF
401 
402 ! Turns off conjugate gradient always.
403  l_gmres = .true.
404 
405 ! Determine levenberg and mu|| parameter guesses based on ||F||.
406  IF (fsq_total1 .lt. 1.e-6_dp) THEN
407 ! IF (fsq_total1 .lt. 1.E-8_dp) THEN
408 
409 ! Pseudo-transient continuation: Turn ogg for constant levmarq_param, muPar
410 ! option.
411  f2 = fsq_total1*1.e12_dp
412 
413 ! NOTE: Asym runs work better like this.
414  f2 = one/(one + f2)
415  IF (lasym) THEN
416  f2 = f2/4
417  END IF
418  f1 = fsq_total1**f2
419 
420  f2 = 1
421  IF (fsq_ratio1 .gt. 0.25_dp .and. &
422  fsq_ratio1 .lt. 1.5_dp) THEN
423  f2 = 0.5_dp
424  END IF
425  IF (fsq_ratio1 .gt. 1.5_dp) THEN
426  f2 = 1.5_dp
427  END IF
428  IF (fsq_ratio1 .gt. 10._dp) THEN
429  f2 = 3
430  END IF
431  levm_loc = max(0.1_dp, min(10._dp, f2*levm_loc))
432 
433  f1 = min(one, f1*levm_scale*levm_loc)
434 
435  levmarq_param = (levmarq_param0 - levm_ped)*f1 + levm_ped
436  IF (levmarq_param0 .eq. zero) THEN
437  levmarq_param = 0
438  END IF
439  mupar = (mupar0-mu_ped)*f1 + mu_ped
440  IF (mupar0 .eq. zero) THEN
441  mupar=0
442  END IF
443 
444  ELSE
445  levmarq_param = max(0.3_dp, levmarq_param0)
446  mupar = mupar0
447  END IF
448 
449 ! Except for the initial step where only initial unpreconditioned forces are
450 ! need, compute full or diagonal approximate hessian.
451  IF (niter .gt. 1) THEN
452  l_linearize = .true.
453  l_getfsq = .false.
454  l_applyprecon = .false.
455  l_getwmhd = .false.
456  ldiagonal = nprecon_type .eq. prediag .and. niter .eq. 2
457  lblock = nprecon_type .eq. preblock
458 
459  IF (ldiagonal .OR. lblock) THEN
460 
461  CALL second0(skston)
462  CALL compute_hessian_blocks(funct_island, ldiagonal)
463 
464  CALL second0(skstoff)
465  IF (ldiagonal) THEN
466  comp_diag_elements_time = comp_diag_elements_time &
467  + (skstoff - skston)
468  END IF
469  IF (lblock) THEN
470  compute_hessian_time = compute_hessian_time + (skstoff - skston)
471  END IF
472  END IF
473 
474  END IF
475 
476 ! Reset run-control logicals. These may have been reset to false in
477 ! Compute_Hessian => funct_island call.
478  l_init_state = .true.
479  l_applyprecon = .false.
480  l_linearize = .false.
481  l_getfsq = .true.
482  fsq_lin = -1
483 
484 ! Choose time-stepping algorithm. Update covariant Fourier components the
485 ! forces to evolve the contravariant components of the displacement.
486  IF (niter .eq. 1) THEN
487  l_init_state = .true.
488  l_getwmhd = .true.
489  CALL second0(skston)
490  CALL funct_island
491  CALL second0(skstoff)
492  evolve_funct_island_time = evolve_funct_island_time &
493  + (skstoff - skston)
494  l_getwmhd = .false.
495 
496  ELSE IF (l_gmres) THEN
497  CALL second0(skston)
498  IF (nprecon .NE. 1) THEN
499  etak_tol = min(0.1_dp, 1.e8_dp*fsq_total)
500  IF (fsq_total <= 0) THEN
501  etak_tol = 0.1_dp
502  END IF
503  END IF
504  CALL gmres_fun
505  CALL second0(skstoff)
506  gmres_time = gmres_time + (skstoff - skston)
507  END IF
508 
509  IF (.not.l_gmres .and. niter .gt. 1) THEN
510  CALL second0(skston)
511  lm0 = levmarq_param
512 ! FIXME: Turn off until RefactorHessian provided for LSCALAPACK = T
513  DO i = 1, 1 !1 + SIZE(levscan)
514  IF (i .GT. 1) THEN
515  levmarq_param = levscan(i - 1)*lm0
516  IF (levmarq_param .gt. levmarq_param0) THEN
517  levmarq_param = levmarq_param0*1.e-2_dp
518  END IF
519  IF (levmarq_param .le. levm_ped) THEN
520  levmarq_param = 100*levm_ped
521  END IF
522  IF (iam .EQ. 0 .and. lverbose) THEN
523  WRITE (*,1000) levmarq_param
524  END IF
525  CALL refactorhessian(levmarq_param)
526  END IF
527 
528  CALL conjugate_grad(ftol)
529  IF (l_conjmin) THEN
530  IF (i .gt. 1) THEN
531  levm_scale = levm_scale*levmarq_param/lm0
532  END IF
533  EXIT
534  END IF
535  END DO
536  CALL second0(skstoff)
537  conj_grad_time = conj_grad_time + (skstoff - skston)
538  END IF
539 
540  v2tot = sqrt(hs*sum(xc*xc))
541 
542  lprint = mod(niter, nprint_step) .eq. 0 .or. niter .eq. 1
543  l_update_state = .true.
544  CALL update_state(lprint, fsq_total1, zero)
545  l_update_state = .false.
546 
547 ! Compute force component residues |Fsubi|**2.
548  IF (fsq_total .lt. zero .and. iam .eq. 0) THEN
549  WRITE (*,1001) fsq_total
550  END IF
551  IF (fsqprev .gt. zero) THEN
552  fsq_ratio = fsq_total/fsqprev
553  fsq_ratio1 = fsq_total1/fsqprev1
554  ELSE
555  fsq_ratio = 0
556  fsq_ratio1 = 0
557  END IF
558  fsqprev = fsq_total
559  fsqprev1 = fsq_total1
560 
561 ! Convert output to REAL (VMEC) units, divide B-fields by b_factor, p by
562 ! p_factor, and WMHD by p_factor
563 
564 #if defined(MPI_OPT)
565  tnum_funct(1) = in_hess_nfunct
566  tnum_funct(2) = out_hess_nfunct
567  CALL mpi_allreduce(mpi_in_place, tnum_funct, 2, mpi_integer, &
568  mpi_max, siesta_comm, mpi_err)
569  nfun_calls = (tnum_funct(1) + tnum_funct(2))
570 #else
571  nfun_calls = (in_hess_nfunct + out_hess_nfunct)
572 #endif
573 
574  IF (niter .eq. 1) THEN
575  CALL bgradp(startglobrow, endglobrow)
576  CALL tflux
577  END IF
578 
579  IF (lprint) THEN
580 ! Needed to compute maximum forces for printout.
581  CALL gather_array(gc)
582 
583 #if defined(MPI_OPT)
584  CALL mpi_allreduce(mpi_in_place, siesta_curtor, 1, mpi_real8, &
585  mpi_sum, siesta_comm, mpi_err)
586 #endif
587 
588  IF (niter .gt. 1 .and. iam .eq. 0) THEN
589  DO iprint = 6, unit_out, unit_out - 6
590  IF (.not.lverbose .and. iprint .eq. 6) THEN
591  cycle
592  END IF
593  WRITE (iprint,1002)
594  END DO
595  CALL write_bounds(jvsupsmncf, fsubsmncf, 'X^S-max:', 'F_S-max:')
596  CALL write_bounds(jvsupumnsf, fsubumnsf, 'X^U-max:', 'F_U-max:')
597  CALL write_bounds(jvsupvmnsf, fsubvmnsf, 'X^V-max:', 'F_V-max:')
598  IF (lasym) THEN
599  DO iprint = 6, unit_out, unit_out - 6
600  IF (.not.lverbose .and. iprint .eq. 6) THEN
601  cycle
602  END IF
603  WRITE (iprint,1003)
604  END DO
605  CALL write_bounds(jvsupsmnsf, fsubsmnsf, 'X^S-max:', 'F_S-max:')
606  CALL write_bounds(jvsupumncf, fsubumncf, 'X^U-max:', 'F_U-max:')
607  CALL write_bounds(jvsupvmncf, fsubvmncf, 'X^V-max:', 'F_V-max:')
608  END IF
609 
610  IF (lverbose) THEN
611  WRITE (*,1004) bsbu_ratio_s, jsju_ratio_s, &
612  (bs0(i), i=2,6), (bu0(i), i=2,6)
613  IF (lasym) THEN
614  WRITE (*,1005) bsbu_ratio_a, jsju_ratio_a, &
615  (bs0(i), i=8,12), (bu0(i), i=8,12)
616  END IF
617  WRITE (*,*)
618  END IF
619 
620  DO iprint = 6, unit_out, unit_out - 6
621  IF (.not.lverbose .and. iprint .eq. 6) THEN
622  cycle
623  END IF
624  WRITE (iprint,1007) siesta_curtor, vmec_curtor
625  END DO
626  END IF
627 
628  DO iprint = 6, unit_out, unit_out - 6
629  IF ((.not.lverbose .and. iprint .eq. 6) .or. iam .ne. 0) THEN
630  cycle
631  END IF
632  IF (niter .EQ. 1) THEN
633  IF (lresistive) THEN
634  str_resistive = "RESISTIVE RUN "
635  ELSE
636  str_resistive = "NON-RESISTIVE RUN "
637  END IF
638  WRITE (iprint, 1006) str_resistive, eta_factor, lasym, &
639  l_push_s, l_push_u, l_push_v, l_push_edge
640  WRITE (iprint, 60) delta_t, etak_tol, l_hess_sym
641  WRITE (iprint, 65) levmarq_param0, mupar0
642  WRITE (iprint, 70) neqs, ns, mpol, ntor, nu_i, nv_i, &
643  ngmres_steps
644  END IF
645 
646  IF (fsq_lin .EQ. -1) THEN
647  WRITE (iprint, 100) &
648  niter, (wtotal-wtotal0)*1.e6_dp/wtotal0, &
649  fsq_total1, fsqvs, fsqvu, &
650  fsqvv, v2tot*delta_t, nfun_calls
651  ELSE
652  WRITE (iprint, 102) &
653  niter, (wtotal-wtotal0)*1.e6_dp/wtotal0, &
654  fsq_total1, fsq_lin, fsqvs, fsqvu, &
655  fsqvv, v2tot*delta_t, nfun_calls
656  END IF
657  END DO
658  END IF
659 
660  CALL second0(skston)
661  IF (fsq_total1 .LT. fsq_min) THEN
662  fsq_min = fsq_total1
663  CALL restart_write(restart_ext, wout_file)
664  END IF
665  CALL second0(skstoff)
666  evolve_restart_file_time=evolve_restart_file_time+(skstoff-skston)
667 
668 ! Add resistive part of B perturbation
669  CALL second0(skston)
670  IF (lresistive .AND. nprecon.GE.1 .AND. fsq_total1.GT.fsq_res) THEN
671  l_getfsq = .false.
672  CALL add_resistive_e
673  END IF
674  CALL second0(skstoff)
675  evolve_add_resistive_e_time=evolve_add_resistive_e_time+(skstoff-skston)
676 
677 1000 FORMAT(/,' Refactoring Hessian: LM = ',1pe12.3)
678 1001 FORMAT('Fsq_Total = ', 1pe12.3,' < 0!')
679 1002 FORMAT('SYMMETRIC DISPLACEMENTS AND FORCES')
680 1003 FORMAT('ASYMMETRIC DISPLACEMENTS AND FORCES')
681 1004 FORMAT(' |B^s-r12*B^u|/|B^s+r12*B^u| (m=1,r->0, sym) : ',1pe10.3,/, &
682  ' |J^s-r12*J^u|/|J^s+r12*J^u| (m=1,r->0, sym) : ',1pe10.3,/, &
683  ' JBSUPSH(JS=2-6,M=1)/r12 (sym): ',1p5e10.3,/, &
684  ' JBSUPUH(JS=2-6,M=1) : ',1p5e10.3)
685 1005 FORMAT(' |B^s+r12*B^u|/|B^s-r12*B^u| (m=1,r->0, asym) : ',1pe10.3,/, &
686  ' |J^s+r12*J^u|/|J^s-r12*J^u| (m=1,r->0, asym) : ',1pe10.3,/, &
687  '-JBSUPSH(JS=2-6,M=1)/r12 (asym): ',1p5e10.3,/, &
688  ' JBSUPUH(JS=2-6,M=1) : ',1p5e10.3)
689 1006 FORMAT(1x, a,' ETA_FACTOR: ',1pe10.2,' LASYM: ',l1, &
690  ' L_PUSH_S: ',l1,' L_PUSH_U: ',l1,' L_PUSH_V: ',l1, &
691  ' L_PUSH_EDGE: ',l1)
692 
693 50 FORMAT(' WMHD: ', 1pe13.6,' <BETA>: ',1pe11.4, &
694  ' TFLUX: ',1pe11.4,' B.GRAD-P (rms): ', 1pe11.4,/,21('-'))
695 60 FORMAT(' DELTA_T: ',1pe10.2,' ETA_K: ', &
696  1pe10.2,' HESSIAN SYM: ', l1)
697 65 FORMAT(' LEVMARQ_PARAM: ',1pe10.2,' MU_PAR: ',1pe10.2)
698 70 FORMAT(' NEQS: ', i6,' NS: ',i4,' MPOL: ',i4,' NTOR: ',i4, &
699  ' NTHETA: ', i4,' NZETA: ',i4, ' NGMRES-STEPS: ', i4,//, &
700  ' NITER (W-W0)/W0*1E6 F2(MHD) F2(LIN)', &
701  ' F2SUBS F2SUBU F2SUBV |V|rms NCALLS')
702 100 FORMAT(1x,i5,1x, 1pe12.4, 2x, 1x,1pe10.3, 2x,9 ('-'), &
703  4(1x,1pe10.3),i9)
704 102 FORMAT(1x,i5,1x, 1pe12.4, 2x, 6(1x,1pe10.3),i9)
705 112 FORMAT(' |B^s+r12*B^u|/|B^s-r12*B^u| (m=1,r->0, asym) : ',1pe10.3,/, &
706  ' |J^s+r12*J^u|/|J^s-r12*J^u| (m=1,r->0, asym) : ',1pe10.3,/, &
707  '-JBSUPSH(JS=2-6,M=1)/r12 (asym): ',1p5e10.3,/, &
708  ' JBSUPUH(JS=2-6,M=1) : ',1p5e10.3)
709 1007 FORMAT(' SIESTA Curtor : ',e12.4,' VMEC Curtor : 'e12.4,/)
710 
711  END SUBROUTINE
712 
713 
715  SUBROUTINE write_bounds(vsupXmn, fsubXmn, vlabel, flabel)
716 !-----------------------------------------------
717 ! D u m m y A r g u m e n t s
718 !-----------------------------------------------
719  REAL(dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(IN) &
720  :: vsupxmn, fsubxmn
721  CHARACTER(LEN=*), INTENT(IN) :: vlabel, flabel
722 !-----------------------------------------------
723 ! L o c a l V a r i a b l e s
724 !-----------------------------------------------
725  INTEGER, PARAMETER :: nm0=1, nn0=2, ns0=3
726  INTEGER :: imaxlocv(3), ilbound(3), imaxlocf(3), iprint
727  REAL(dp) :: vmax, fmax
728 !-----------------------------------------------
729 
730  ilbound = lbound(vsupxmn)-1
731  imaxlocv = maxloc(abs(vsupxmn)) + ilbound
732  vmax = vsupxmn(imaxlocv(1),imaxlocv(2),imaxlocv(3))
733  CALL assert(abs(vmax).EQ.maxval(abs(vsupxmn)),'vmax WRONG')
734  imaxlocv = imaxlocv - ilbound
735  imaxlocf = maxloc(abs(fsubxmn)) + ilbound
736  fmax = fsubxmn(imaxlocf(1),imaxlocf(2),imaxlocf(3))
737  CALL assert(abs(fmax).EQ.maxval(abs(fsubxmn)),'fmax WRONG')
738  imaxlocf = imaxlocf - ilbound
739 
740  DO iprint = 6, unit_out, unit_out-6
741  IF ((.NOT.lverbose .AND. iprint.EQ.6) .OR. iam.NE.0) cycle
742  WRITE(iprint, 45) &
743  trim(vlabel),vmax,imaxlocv(ns0),imaxlocv(nm0)-1, &
744  imaxlocv(nn0)-ntor-1, &
745  trim(flabel),fmax,imaxlocf(ns0),imaxlocf(nm0)-1, &
746  imaxlocf(nn0)-ntor-1
747  END DO
748 
749  45 FORMAT(2(1x,a,1x,1pe10.2,' AT JS: ',i4,' M: ',i4,' N: ',i4,2x))
750 
751  END SUBROUTINE write_bounds
752 
753 
755  SUBROUTINE conjugate_grad (ftol)
756  USE stel_kinds
757  USE hessian, ONLY: mupar
758  USE siesta_state, ONLY: update_state, update_state
759 !-----------------------------------------------
760 ! D u m m y A r g u m e n t s
761 !-----------------------------------------------
762  REAL(dp) :: ftol
763 !-----------------------------------------------
764 ! L o c a l V a r i a b l e s
765 !-----------------------------------------------
766  REAL(dp), PARAMETER :: p5 = 0.50_dp
767  INTEGER, PARAMETER :: ndamp = 10
768  INTEGER :: icount, iter
769  REAL(dp) :: b1, fac, fsq_start, fsq_min0, fsq_save, fsq_conv, &
770  wmhd_min, wmhd_start, mupars, dtau, delt_cg, otav
771  REAL(dp), ALLOCATABLE :: xcstart(:)
772 !-----------------------------------------------
773  IF (fsq_total1.LE.ftol .OR. nprecon.EQ.0) THEN
774  IF (nprecon .EQ. 0) nprecon = 1
775  RETURN
776  END IF
777 
778  l_getfsq = .true.
779  l_linearize = .false.
780  l_applyprecon = .true.
781  l_init_state = .true.
782  l_printoriginforces = .false.
783  l_getwmhd=.true.
784 
785 !
786 ! Linearize perturbations around B0 and p0
787 ! dB = curl(xc X B0)
788 ! dp = dp(p0,xc)
789 !
790 ! NON-LINEAR forces
791 ! J = J0+dJ, B = B0+dB, p = p0+dp
792 
793  ALLOCATE(xcdot(neqs), xc0(neqs), stat=icount)
794  CALL assert(icount.eq.0,'ALLOCATION FAILED IN CONJ_GRAD')
795 
796  xcdot = 0
797  xc = 0
798  xc0 = 0
799 
800  delt_cg = 1
801  dtau = 0
802  fsq_save = 0
803 
804  IF (iam .EQ. 0 .and. lverbose) print 90
805 90 FORMAT(1x,'CONJUGATE GRADIENT CONVERGENCE SUMMARY', &
806  /,' -------------',/,1x,'ITER',7x,'FSQ_NL',10x,'||X||')
807 
808 ! muParS = muPar; ! muPar=0
809 
810  DO icount = 1, 10*ncongrad_steps
811 !
812 ! LOOP OVER delt_cg TO FIND MINIMUM F(xc)
813 !
814  CALL funct_island
815 
816  CALL update_taudamp (fsq_save, delt_cg)
817 
818  IF (icount .EQ. 1) THEN
819  fsq_start = fsq_total1
820  fsq_min0 = fsq_total1
821  fsq_conv = fsq_min0
822  fsq_save = fsq_total
823  wmhd_start = wtotal
824  wmhd_min = wtotal
825  ALLOCATE (xcstart(SIZE(gc)))
826  xcstart = -gc
827  END IF
828 
829  IF (fsq_total1 .LT. fsq_min0) THEN
830  fsq_min0 = fsq_total1
831  xc0 = xc
832  IF (fsq_total1 .LT. ftol) EXIT
833  END IF
834 ! ELSE IF (wtotal .lt. wmhd_min) THEN
835 ! wmhd_min = wtotal
836 ! xc0 = xc
837 ! END IF
838  IF (iam.EQ.0 .AND. (icount.EQ.1 .OR. mod(icount,10).EQ.0) .and. lverbose) &
839  print 110, icount, fsq_total1, sqrt(sum(xc*xc))
840 
841  IF (fsq_total1 .GT. 1.5_dp*fsq_min0) THEN
842  delt_cg = 0.95_dp*delt_cg
843  IF (delt_cg .LT. 0.001_dp) EXIT
844  xc = xc0
845  xcdot = 0
846  IF (fsq_total1 .GT. 10*fsq_min0) THEN
847  delt_cg = delt_cg*p5
848  fsq_save = 0
849  cycle
850  END IF
851  END IF
852 
853  IF (mod(icount, 50) .EQ. 0) THEN
854  IF (fsq_min0 .GT. fsq_conv/2) EXIT
855  fsq_conv = fsq_min0
856  END IF
857 
858  CALL assert(.NOT.l_linearize,'LINEARIZATION TURNED ON!')
859  otav = sum(otau)/ndamp
860  dtau = delt_cg*otav/2
861 
862  b1 = one-dtau
863  fac = one/(one+dtau)
864 !
865 ! THIS IS THE TIME-STEPPING ALGORITHM. IT IS ESSENTIALLY A CONJUGATE
866 ! GRADIENT METHOD, WITHOUT THE LINE SEARCHES (FLETCHER-REEVES),
867 ! BASED ON A METHOD BY P. GARABEDIAN
868 !
869 ! CONJ GRAD
870  xcdot = fac*(b1*xcdot + delt_cg*gc)
871  xc = xc + delt_cg*xcdot
872 !
873 ! STEEPEST DESCENT (b1->0)
874 ! xc = b1*xc - delt_cg*gc
875 
876  END DO
877 
878 !
879 ! UPDATE FIELD, PRESSURE PERTS AND UPDATE NEW NON-LINEAR STATE
880 !
881  mupars = mupar; mupar = 0
882  l_printoriginforces = .false.
883 
884 ! IF (fsq_min0 .EQ. fsq_start) xc0 = xcstart
885 ! CALL LineSearch(xc0, fsq_min0)
886 
887  l_conjmin = (fsq_min0 .LT. fsq_start)
888  IF (l_conjmin) THEN
889  xc = xc0
890  ELSE
891  xc = 0
892  END IF
893 
894  IF (fsq_min0 .GT. 0.95_dp*fsq_start) l_conjmin=.false.
895  IF (.NOT.l_gmres .AND. l_conjmin) l_printoriginforces = .true.
896  l_applyprecon = .false.
897  CALL funct_island
898 
899  l_printoriginforces = .false.
900 
901  IF (mupars .NE. zero) fsq_min0 = fsq_total1
902 
903  IF (iam.EQ.0 .AND. l_conjmin) THEN
904  IF (lverbose) WRITE (6, 100) icount, fsq_start, fsq_min0
905  WRITE (unit_out,100) icount, fsq_start, fsq_min0
906  END IF
907 
908  mupar = mupars
909 
910  DEALLOCATE (xcdot, xc0, xcstart)
911  100 FORMAT(' ITER_CG: ', i3,' FSQ (START CG): ',1p,e12.3,' FSQ (END CG): ', 1pe12.3)
912  110 FORMAT(i5, 2(3x,1pe12.3))
913 
914  END SUBROUTINE conjugate_grad
915 
917 
918  SUBROUTINE update_taudamp (fsq_save, delt_cg)
919 !-----------------------------------------------
920 ! D u m m y A r g u m e n t s
921 !-----------------------------------------------
922  REAL(dp), INTENT(INOUT) :: fsq_save, delt_cg
923 !-----------------------------------------------
924 ! L o c a l V a r i a b l e s
925 !-----------------------------------------------
926  REAL(dp), PARAMETER :: mintau = 0.15_dp
927  REAL(dp) :: ratio
928 !-----------------------------------------------
929 
930  IF (fsq_save .GT. zero) THEN
931  ratio = abs(fsq_total/fsq_save)
932  ratio = min(mintau, ratio)
933  otau(1:ndamp-1) = otau(2:ndamp)
934  otau(ndamp) = ratio/delt_cg
935  fsq_save = fsq_total
936  ELSE
937  otau = mintau/delt_cg
938  END IF
939 
940  END SUBROUTINE update_taudamp
941 
942  END MODULE evolution
943 
shared_data::etak_tol
real(dp) etak_tol
FIXME: UNKNOWN.
Definition: shared_data.f90:117
shared_data::nprecon
integer nprecon
Preconditioner flag.
Definition: shared_data.f90:68
shared_data::out_hess_nfunct
integer out_hess_nfunct
FIXME UNKNOWN.
Definition: shared_data.f90:87
shared_data::l_push_edge
logical l_push_edge
Solve u,v components at s=1.
Definition: shared_data.f90:202
shared_data::lverbose
logical lverbose
Use verbose screen output.
Definition: shared_data.f90:234
shared_data::l_getfsq
logical l_getfsq
Compute |F|^2.
Definition: shared_data.f90:216
shared_functions
Module contained subroutines and functions updating MHD forces and Wmhd.
Definition: shared_functions.f90:4
shared_data::xc
real(dp), dimension(:), allocatable xc
1D array of Fourier mode displacement components.
Definition: shared_data.f90:97
shared_data::in_hess_nfunct
integer in_hess_nfunct
FIXME UNKNOWN.
Definition: shared_data.f90:85
shared_data::l_push_u
logical l_push_u
Solve for u component at origin.
Definition: shared_data.f90:206
shared_data::scale_u
real(dp) scale_u
FIXME: UNKNOWN.
Definition: shared_data.f90:149
hessian::gather_array
Definition: hessian.f90:31
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
shared_data::delta_t
real(dp) delta_t
Time step.
Definition: shared_data.f90:125
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
shared_data::l_push_s
logical l_push_s
Solve for s component at origin.
Definition: shared_data.f90:204
shared_data::bsbu_ratio_s
real(dp) bsbu_ratio_s
FIXME: UNKNOWN.
Definition: shared_data.f90:139
shared_data::l_applyprecon
logical l_applyprecon
Apply preconditioner.
Definition: shared_data.f90:218
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
shared_data::siesta_curtor
real(dp) siesta_curtor
Total toroidal current.
Definition: shared_data.f90:239
shared_data::preblock
integer, parameter preblock
Block preconditioning flag.
Definition: shared_data.f90:25
shared_data::l_init_state
logical l_init_state
Store initial field/pressure state.
Definition: shared_data.f90:222
shared_data::fsupumncf
real(dp), dimension(:,:,:), pointer fsupumncf
Contravariant force for stellarator symmetric u component on full grid.
Definition: shared_data.f90:192
shared_data::scale_s
real(dp) scale_s
FIXME: UNKNOWN.
Definition: shared_data.f90:147
shared_data::l_printoriginforces
logical l_printoriginforces
Print forces at the origin.
Definition: shared_data.f90:220
shared_data::fsupsmnsf
real(dp), dimension(:,:,:), pointer fsupsmnsf
Contravariant force for stellarator symmetric s component on full grid.
Definition: shared_data.f90:188
shared_data::fsqvv
real(dp) fsqvv
|F|^2 for v components.
Definition: shared_data.f90:131
v3_utilities::assert
Definition: v3_utilities.f:55
shared_data::levm_ped
real(dp), parameter levm_ped
FIXME: UNKNOWN.
Definition: shared_data.f90:29
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
shared_data::wtotal
real(dp) wtotal
MHD energy sum of magnetic and thermal.
Definition: shared_data.f90:121
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
shared_data::mu_ped
real(dp), parameter mu_ped
Pedestal value of levenberg/mu. Should be between 10^-5 and 10^-10.
Definition: shared_data.f90:31
siesta_namelist::eta_factor
real(dp) eta_factor
Resistivity value.
Definition: siesta_namelist.f90:148
siesta_namelist::wout_file
character(len=siesta_namelist_name_length) wout_file
Filename of the VMEC woutfile.
Definition: siesta_namelist.f90:189
shared_data::jsju_ratio_s
real(dp) jsju_ratio_s
FIXME: UNKNOWN.
Definition: shared_data.f90:141
shared_data::fsupumnsf
real(dp), dimension(:,:,:), pointer fsupumnsf
Contravariant force for stellarator non-symmetric u component on full grid.
Definition: shared_data.f90:194
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
shared_data::l_linearize
logical l_linearize
Use linearized forces.
Definition: shared_data.f90:210
restart_mod
Contains routines for writting the restart file.
Definition: restart_mod.f90:10
shared_data::fsupsmncf
real(dp), dimension(:,:,:), pointer fsupsmncf
Contravariant force for stellarator non-symmetric s component on full grid.
Definition: shared_data.f90:190
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
island_params::hs_i
real(dp) hs_i
Radial grid spacing. hs = s_i+1 - s_i.
Definition: island_params.f90:47
shared_data::nprecon_type
integer nprecon_type
Preconditioner type.
Definition: shared_data.f90:70
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
shared_data::gc
real(dp), dimension(:), allocatable, target gc
1D Array of Fourier mode MHD force components
Definition: shared_data.f90:99
shared_data::l_update_state
logical l_update_state
Update the ste array.
Definition: shared_data.f90:224
shared_data::fsq_res
real(dp), parameter fsq_res
Threshold force for turning off resistive perturbations.
Definition: shared_data.f90:27
shared_data::ste
real(dp), dimension(4) ste
Spectral Truncation RMS error,.
Definition: shared_data.f90:133
shared_data::fsqvu
real(dp) fsqvu
|F|^2 for u components.
Definition: shared_data.f90:129
shared_data::xc0
real(dp), dimension(:), allocatable xc0
Saved fouier displacements.
Definition: shared_data.f90:106
siesta_namelist::lresistive
logical lresistive
Use resistive perturbaton.
Definition: siesta_namelist.f90:126
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
shared_data::l_getwmhd
logical l_getwmhd
Compute MHD energy.
Definition: shared_data.f90:214
shared_data::fsq_total1
real(dp) fsq_total1
|F|^2 WITHOUT column scaling.
Definition: shared_data.f90:93
shared_data::ndims
integer ndims
Number of independent variables.
Definition: shared_data.f90:58
shared_data::fsupvmncf
real(dp), dimension(:,:,:), pointer fsupvmncf
Contravariant force for stellarator symmetric v component on full grid.
Definition: shared_data.f90:196
shared_data::prediag
integer, parameter prediag
Diagonal preconditioning flag.
Definition: shared_data.f90:23
shared_data::levm_scale
real(dp) levm_scale
FIXME: UNKNOWN.
Definition: shared_data.f90:119
restart_mod::restart_write
subroutine restart_write(restart_ext, wout_file)
Write the restart file.
Definition: restart_mod.f90:363
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::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
shared_data::fsqvs
real(dp) fsqvs
|F|^2 for s components.
Definition: shared_data.f90:127
siesta_namelist::ladd_pert
logical ladd_pert
Use helical perturbation.
Definition: siesta_namelist.f90:124
shared_data::fsupvmnsf
real(dp), dimension(:,:,:), pointer fsupvmnsf
Contravariant force for stellarator non-symmetric v component on full grid.
Definition: shared_data.f90:198
shared_data::wtotal0
real(dp) wtotal0
Saved MHD energy sum of magnetic and thermal.
Definition: shared_data.f90:123
shared_data::fsq_lin
real(dp) fsq_lin
Linear |F|^2.
Definition: shared_data.f90:115
shared_data::neqs
integer neqs
Number of elements in the xc array.
Definition: shared_data.f90:54
shared_data::ngmres_steps
integer, parameter ngmres_steps
Max number of gmres steps (10-100) should scale with ns.
Definition: shared_data.f90:21
shared_data::l_push_v
logical l_push_v
Solve for v component at origin.
Definition: shared_data.f90:208
siesta_namelist::l_output_alliter
logical l_output_alliter
Write output files on all iterations.
Definition: siesta_namelist.f90:136
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17