Merge pull request #38 from ORNL-Fusion/sparse

Fix issues using the SIESTA run context in V3FIT.

Sources/evolution.f90

@@ -145,7 +145,7 @@ SUBROUTINE converge_diagonal(wout_file, ftol)
                       niter      .lt. niter_max) .or.                          &
                      (niter .eq. 1 .and. niter_max .ne. 1)
 
-         IF (.not.pert_added .and. fsq_total1.lt.100*fsq_block) THEN
+         IF (.not.pert_added .and. fsq_total1 .lt. 100*fsq_block) THEN
             l_init_state = .true.
             CALL second0(t1)
             CALL add_perturb(xc, getwmhd)

Sources/gmres.f90

@@ -151,32 +151,40 @@ SUBROUTINE gmres_fun
 !
 !     THIS IS STANDARD GMRES CALL
 !
-      CALL second0(skston)
-      xc = xcmin
-      CALL gmres_wrap (xc, brhs)
-      CALL second0(skstoff)
-      gmres_wrap_time=gmres_wrap_time+(skstoff-skston)
-
-      IF (fsq_total1 .LT. fsq_min .OR. ALL(xcmin .EQ. zero)) THEN
-         xcmin = xc
-         fsq_min = fsq_total1
-         IF (lm0 .GT. zero) scale_fac = levmarq_param/lm0
-      END IF
+         CALL second0(skston)
+         xc = xcmin
+         CALL gmres_wrap (xc, brhs)
+         CALL second0(skstoff)
+         gmres_wrap_time=gmres_wrap_time+(skstoff-skston)
+
+         IF (fsq_total1 .LT. fsq_min .OR. ALL(xcmin .EQ. zero)) THEN
+            xcmin = xc
+            fsq_min = fsq_total1
+            IF (lm0 .GT. zero) THEN
+               scale_fac = levmarq_param/lm0
+            END IF
+         END IF
 !!!END GMRES CONVERGENCE
 !FOR NOW, THIS IS ONLY IMPLEMENTED FOR PARSOLVER=TRUE: MUST IMPLEMENT
 !RefactorHessian FOR LSCALAPACK=T
-      IF (.FALSE.) THEN
-!      IF (.NOT.l_Diagonal_Only .AND. PARSOLVER) THEN
-         IF (fsq_min.GE.0.95_dp*fmhd .AND. iter.LT.SIZE(levscan)) THEN
-            iter = iter+1
-            levmarq_param = lm0*levscan(iter)
-            IF (levmarq_param .LT. levm_ped) levmarq_param = 100*levm_ped
-            IF (levmarq_param .GE. levmarq_param0) levmarq_param = levmarq_param0/10._dp**iter
-            IF (iam.EQ.0 .AND. lverbose) PRINT 930,' Refactoring Hessian: LM = ',levmarq_param
-            CALL RefactorHessian(levmarq_param)
-            GOTO 920
+         IF (.FALSE.) THEN
+!        IF (.NOT.l_Diagonal_Only .AND. PARSOLVER) THEN
+            IF (fsq_min.GE.0.95_dp*fmhd .AND. iter.LT.SIZE(levscan)) THEN
+               iter = iter+1
+               levmarq_param = lm0*levscan(iter)
+               IF (levmarq_param .LT. levm_ped) THEN
+                  levmarq_param = 100*levm_ped
+               END IF
+               IF (levmarq_param .GE. levmarq_param0) THEN
+                  levmarq_param = levmarq_param0/10._dp**iter
+               END IF
+               IF (iam.EQ.0 .AND. lverbose) THEN
+                  PRINT 930,' Refactoring Hessian: LM = ',levmarq_param
+               END IF
+               CALL RefactorHessian(levmarq_param)
+               GOTO 920
+            END IF
          END IF
-      END IF
       END IF LGMRES
 
  912  FORMAT(i5,3(3x,1pe12.3))
@@ -231,7 +239,7 @@ SUBROUTINE gmres_fun
       muPar = muParS
 
       DEALLOCATE (xcmin, brhs)
-      
+
       END SUBROUTINE gmres_fun
 
       SUBROUTINE init_gmres0 (icntl, cntl, etak, m, n)

Sources/grid_extension.f90

@@ -11,7 +11,8 @@ MODULE grid_extension
       USE vmec_info, ONLY: jcurrumnc, jcurrvmnc
       USE stel_kinds
       USE stel_constants
-      USE shared_data, ONLY: lasym, r1_i, z1_i, ru_i, zu_i, rv_i, zv_i, nsp
+      USE shared_data, ONLY: lasym, r1_i, z1_i, ru_i, zu_i, rv_i, zv_i, nsp,   &
+                             lverbose
       USE island_params, ns=>ns_i, mpol=>mpol_i, ntor=>ntor_i, &
                          ntheta=>nu_i, nzeta=>nv_i, mnmax=>mnmax_i, nsv=>ns_v, &
                          ohs=>ohs_i, nfp=>nfp_i
@@ -720,7 +721,7 @@ FUNCTION read_vessel_file()
   100 CONTINUE
 
       REWIND(UNIT=iou, iostat=read_vessel_file)
-      IF (read_vessel_file .NE. 0) THEN
+      IF (read_vessel_file .NE. 0 .and. lverbose) THEN
          WRITE (*,*) vessel_file
          RETURN
       END IF

Sources/hessian.f90

@@ -1341,9 +1341,12 @@ SUBROUTINE CheckEigenvalues_Serial(nblock, bsize)
       IF (JOBZ .NE. 'N') DEALLOCATE(Z)
 
       CALL second0(toff)
-      PRINT 100,' Eigenvalue TIME: ', (toff-ton), ' INFO: ', info,      &
-                ' Eigenvalues written to FORT.',4000+ncount
- 100  FORMAT(a,1p,e10.2,2(a,i4))      
+      IF (lverbose) THEN
+         PRINT 100,' Eigenvalue TIME: ', (toff-ton), ' INFO: ', info,          &
+                   ' Eigenvalues written to FORT.',4000+ncount
+      END IF
+
+ 100  FORMAT(a,1p,e10.2,2(a,i4))
 
       END SUBROUTINE CheckEigenvalues_Serial
 

Sources/metrics.f90

@@ -16,14 +16,7 @@ MODULE metrics
       USE island_params
       USE shared_data, ONLY: lasym, r1_i, z1_i, ru_i, zu_i, rv_i, zv_i,  &
                              lverbose, jsupvdotA, nsp
-      USE read_wout_mod, ns_vmec=>ns, mpol_vmec=>mpol, ntor_vmec=>ntor,  &
-          rmnc_vmec=>rmnc, zmns_vmec=>zmns, lmns_vmec=>lmns,             &
-          xm_vmec=>xm, xn_vmec=>xn, chipf_vmec=>chipf,                   &  ! MRC 4/1/2016
-          rmns_vmec=>rmns, zmnc_vmec=>zmnc, lmnc_vmec=>lmnc,             &  ! MRC 12/1/2016
-          phipf_vmec=>phipf, presf_vmec=>presf, nfp_vmec=>nfp,           &
-          wb_vmec=>wb, wp_vmec=>wp, gamma_vmec=>gamma,                   &
-          volume_vmec=>volume, raxis_vmec=>raxis, lasym_vmec=>lasym,     &
-          iasym_vmec=>iasym
+      USE v3_utilities
       USE descriptor_mod, ONLY: iam
       USE timer_mod
       USE utilities, ONLY: to_full_mesh
@@ -226,10 +219,14 @@ SUBROUTINE LoadGrid(istat)
       INTEGER, INTENT(out) :: istat
 
 !  Start executable code.
-      ALLOCATE (r1_i(nu_i, nv_i, ns_i), z1_i(nu_i, nv_i, ns_i),         &
-                ru_i(nu_i, nv_i, ns_i), zu_i(nu_i, nv_i, ns_i),         &
-                rv_i(nu_i, nv_i, ns_i), zv_i(nu_i, nv_i, ns_i),         &
-                stat = istat)
+      IF (.not.ALLOCATED(r1_i)) THEN
+         ALLOCATE (r1_i(nu_i, nv_i, ns_i), z1_i(nu_i, nv_i, ns_i),             &
+                   ru_i(nu_i, nv_i, ns_i), zu_i(nu_i, nv_i, ns_i),             &
+                   rv_i(nu_i, nv_i, ns_i), zv_i(nu_i, nv_i, ns_i),             &
+                   stat = istat)
+      ELSE
+         istat = 0
+      END IF
 
       CALL rz_to_ijsp(rmnc_i, zmns_i, .false.)
       IF (lasym) THEN
@@ -480,17 +477,17 @@ SUBROUTINE check_metrics
 
 !  Start executable code.
       CALL ASSERT(ALL(ABS(hss*gssf + hsu*gsuf + hsv*gsvf - 1) .lt. tolarance), &
-                      's Diagonal metric element failed.')
+                  's Diagonal metric element failed.')
       CALL ASSERT(ALL(ABS(hsu*gsuf + huu*guuf + huv*guvf - 1) .lt. tolarance), &
-                      'u Diagonal metric element failed.')
+                  'u Diagonal metric element failed.')
       CALL ASSERT(ALL(ABS(hsv*gsvf + huv*guvf + hvv*gvvf - 1) .lt. tolarance), &
-                      'v Diagonal metric element failed.')
+                  'v Diagonal metric element failed.')
       CALL ASSERT(ALL(ABS(hss*gsuf + hsu*guuf + hsv*guvf) .lt. tolarance),     &
-                      'su Off diagonal metric element failed.')
+                  'su Off diagonal metric element failed.')
       CALL ASSERT(ALL(ABS(hss*gsvf + hsu*guvf + hsv*gvvf) .lt. tolarance),     &
-                      'sv Off diagonal metric element failed.')
+                  'sv Off diagonal metric element failed.')
       CALL ASSERT(ALL(ABS(hsu*gsvf + huu*guvf + huv*gvvf) .lt. tolarance),     &
-                      'uv Off diagonal metric element failed.')
+                  'uv Off diagonal metric element failed.')
 
       END SUBROUTINE
 

Sources/nscalingtools.f90

@@ -74,7 +74,7 @@ SUBROUTINE  MyEnvVariables
 !--------------------------------------------------------------------------
 ! Default values of environment values that users are allowed to change
 !--------------------------------------------------------------------------
-!PARSOLVER=.TRUE.
+PARSOLVER=.FALSE.
 PARFUNCTISL=.TRUE.
 !--------------------------------------------------------------------------
 

Sources/pchelms.f90

@@ -21,7 +21,7 @@ MODULE pchelms
                                endglobrow, rcounts, MPI_ERR
       USE shared_data, ONLY: gc, xc, neqs, gc0, xc0, fsq_total, mblk_size, ndims
       USE shared_data, ONLY: asubsmnsf, asubumncf, asubvmncf,                  &
-                             asubsmncf, asubumnsf, asubvmnsf
+                             asubsmncf, asubumnsf, asubvmnsf, lverbose
       USE quantities, ONLY: fsupsmnsf, fsupumncf, fsupvmncf,                   &
                             fsupsmncf, fsupumnsf, fsupvmnsf
       USE shared_data, ONLY: nsp, r1_i, z1_i, ru_i,                            &
@@ -1332,7 +1332,7 @@ SUBROUTINE COMPUTE_FORCES(lNLinear)
                         f_cos)
          END IF
 
-         IF (iam .eq. 0) THEN
+         IF (iam .eq. 0 .and. lverbose) THEN
             WRITE (*,*) '<chiedge> = ',                                        &
                twopi*asubvmncf(1,ntor + 1,ns)*fourier_context%orthonorm(0,0)
             WRITE (*,*) '<phiedge> = ',                                        &

Sources/perturbation.f90

@@ -381,7 +381,7 @@ FUNCTION getwmhd(p)
       lresistive = lresist_save
       eta_factor = eta_factor_save
 
-1000  FORMAT(' Adding helical magnetic field perturbations'/' 10^6 X Del-W',      &
+1000  FORMAT(' Adding helical magnetic field perturbations'/' 10^6 X Del-W',   &
              '    mres    nres     HelPert     rad  |m*chip+n*phip|',          &
              '     iota   radial width')
 1001  FORMAT(1p,e12.3,2(i8),3x,1pe10.2,0p,f8.2,f11.2,4x,2(f11.2))

Sources/quantities.f90

@@ -342,8 +342,10 @@ SUBROUTINE init_quantities(restarted)
          wint(l,:,:) = fourier_context%cosmui(l,m0)
       END DO
 
-      ALLOCATE (vp_f(ns), stat=istat)
-      CALL assert_eq(0, istat, 'Allocation error in init_quantities')
+      IF (.not.ALLOCATED(vp_f)) THEN
+         ALLOCATE (vp_f(ns), stat=istat)
+         CALL assert_eq(0, istat, 'Allocation error in init_quantities')
+      END IF
       CALL SurfAverage(vp_f, jacobf, 1, ns)
 
 !  Compar volumes for a sanity check that the siesta grid and orginal VMEC
@@ -796,37 +798,46 @@ SUBROUTINE alloc_quantities
 
 !  Start of executable code
 !  Half mesh quantities (harmonics), B^s (sine), B^u (cosine), B^v (cosine)
-      ALLOCATE(jbsupsmnsh(0:mpol,-ntor:ntor,ns),                               &
-               jbsupumnch(0:mpol,-ntor:ntor,ns),                               &
-               jbsupvmnch(0:mpol,-ntor:ntor,ns),                               &
-               jpmnch    (0:mpol,-ntor:ntor,ns), stat=istat)
-      CALL assert_eq(0, istat, 'Allocation #1 failed in alloc_quantities')
+      IF (.not.ALLOCATED(jbsupsmnsh)) THEN
+         ALLOCATE(jbsupsmnsh(0:mpol,-ntor:ntor,ns),                            &
+                  jbsupumnch(0:mpol,-ntor:ntor,ns),                            &
+                  jbsupvmnch(0:mpol,-ntor:ntor,ns),                            &
+                  jpmnch    (0:mpol,-ntor:ntor,ns), stat=istat)
+         CALL assert_eq(0, istat, 'Allocation #1 failed in alloc_quantities')
+      END IF
       jbsupsmnsh = zero
       jbsupumnch = zero
       jbsupvmnch = zero
       jpmnch = zero
 
-      ALLOCATE(pwr_spec_s(0:mpol,-ntor:ntor,ns,4),                             &
-               pwr_spec_a(0:mpol,-ntor:ntor,ns,4), stat=istat)
-      CALL assert_eq(0, istat, 'Allocation #2 failed in alloc_quantities')
+      IF (.not.ALLOCATED(pwr_spec_s)) THEN
+         ALLOCATE(pwr_spec_s(0:mpol,-ntor:ntor,ns,4),                          &
+                  pwr_spec_a(0:mpol,-ntor:ntor,ns,4), stat=istat)
+         CALL assert_eq(0, istat, 'Allocation #2 failed in alloc_quantities')
+      END IF
 
       IF (lasym) THEN
 !  Half mesh quantities (harmonics), B^s (sine), B^u (cosine), B^v (cosine)
-         ALLOCATE(jbsupsmnch(0:mpol,-ntor:ntor,ns),                            &
-                  jbsupumnsh(0:mpol,-ntor:ntor,ns),                            &
-                  jbsupvmnsh(0:mpol,-ntor:ntor,ns),                            &
-                  jpmnsh    (0:mpol,-ntor:ntor,ns), stat=istat)
-         CALL assert_eq(0, istat, 'Allocation #3 failed in alloc_quantities')
+         IF (.not.ALLOCATED(jbsupsmnsh)) THEN
+            ALLOCATE(jbsupsmnch(0:mpol,-ntor:ntor,ns),                         &
+                     jbsupumnsh(0:mpol,-ntor:ntor,ns),                         &
+                     jbsupvmnsh(0:mpol,-ntor:ntor,ns),                         &
+                     jpmnsh    (0:mpol,-ntor:ntor,ns), stat=istat)
+            CALL assert_eq(0, istat,                                           &
+     &                     'Allocation #3 failed in alloc_quantities')
+         END IF
          jbsupsmnch = zero
          jbsupumnsh = zero
          jbsupvmnsh = zero
          jpmnsh = zero
       END IF
-
-      ALLOCATE(bsq(ntheta,nzeta,ns), jacobh(ntheta,nzeta,ns),                  &
-               jacobf(ntheta,nzeta,ns), wint(ntheta,nzeta,ns), stat=istat)
+
+      IF (.not.ALLOCATED(bsq)) THEN
+         ALLOCATE(bsq(ntheta,nzeta,ns), jacobh(ntheta,nzeta,ns),               &
+                  jacobf(ntheta,nzeta,ns), wint(ntheta,nzeta,ns), stat=istat)
+         CALL assert_eq(0, istat, 'Allocation #4 failed in alloc_quantities')
+      END IF
       bsq = 0.0
-      CALL assert_eq(0, istat, 'Allocation #4 failed in alloc_quantities')
 
       END SUBROUTINE
 

Sources/siesta_namelist.f90

@@ -240,7 +240,7 @@ MODULE siesta_namelist
       SUBROUTINE siesta_namelist_read(namelist_file)
       USE safe_open_mod
       USE v3_utilities
-      USE Hessian, ONLY: levmarq_param, mupar, levmarq_param0, mupar0
+      USE Hessian, ONLY: levmarq_param0, mupar0
 
       IMPLICIT NONE
 

Sources/siesta_run.f90

@@ -26,6 +26,8 @@ MODULE siesta_run
       INTEGER, PARAMETER :: siesta_run_control_mpi = 1
 !>  Bit position to close the wout file.
       INTEGER, PARAMETER :: siesta_run_control_wout = 2
+!>  Bit position to sycn the namelist inputs.
+      INTEGER, PARAMETER :: siesta_run_sync_namelist = 3
 
 !*******************************************************************************
 !  DERIVED-TYPE DECLARATIONSf
@@ -46,6 +48,7 @@ MODULE siesta_run
          PROCEDURE :: set_1d => siesta_run_set_1d
          GENERIC   :: set => set_1d
          PROCEDURE :: converge => siesta_run_converge
+         PROCEDURE :: sync => siesta_run_sync
          FINAL     :: siesta_run_destruct
       END TYPE
 
@@ -416,6 +419,8 @@ SUBROUTINE siesta_run_set_vmec_temp(this, load_wout)
       INTEGER                                 :: istat
 
 !  Start of executable code
+      CALL this%sync
+
       CALL vmec_info_set_wout(wout_file, nsin, mpolin, ntorin, nfpin,          &
      &                        ntor_modes(-ntorin:ntorin), load_wout)
 
@@ -478,6 +483,8 @@ SUBROUTINE siesta_run_set_restart(this)
       CLASS (siesta_run_class), INTENT(inout) :: this
 
 !  Start of executable code
+      CALL this%sync
+
       IF (l_vessel) THEN
          ns = nsin + nsin_ext
       ELSE
@@ -523,6 +530,8 @@ SUBROUTINE siesta_run_set_1d(this, param_name, value, index)
 
          CASE ('helpert')
             helpert(index) = value
+            this%control_state = IBSET(this%control_state,                     &
+                                       siesta_run_sync_namelist)
 
          CASE DEFAULT
             CALL siesta_error_set_error(siesta_error_param,                    &
@@ -615,4 +624,38 @@ SUBROUTINE siesta_run_converge(this)
 
       END SUBROUTINE
 
+!*******************************************************************************
+!  MPI SUBROUTINES
+!*******************************************************************************
+!-------------------------------------------------------------------------------
+!>  @brief Sync child equilibria.
+!>
+!>  @param[inout] this A @ref siesta_run_class instance.
+!-------------------------------------------------------------------------------
+      SUBROUTINE siesta_run_sync(this)
+      USE siesta_namelist, ONLY: helpert
+      USE nscalingtools, ONLY: SIESTA_COMM, MPI_ERR
+#if defined(MPI_OPT)
+      USE mpi_inc
+#endif
+
+      IMPLICIT NONE
+
+!  Declate Arguments
+      CLASS (siesta_run_class), INTENT(inout) :: this
+
+!  Start of executable code
+#if defined(MPI_OPT)
+      CALL MPI_BCAST(this%control_state, 1, MPI_INTEGER, 0, SIESTA_COMM,       &
+                     MPI_ERR)
+      IF (BTEST(this%control_state, siesta_run_sync_namelist)) THEN
+         CALL MPI_BCAST(helpert, SIZE(helpert), MPI_REAL8, 0, SIESTA_COMM,     &
+                        MPI_ERR)
+         this%control_state = IBCLR(this%control_state,                        &
+                                    siesta_run_sync_namelist)
+      END IF
+#endif
+
+      END SUBROUTINE
+
       END MODULE

Sources/vmec_info.f90

@@ -146,6 +146,7 @@ SUBROUTINE vmec_info_construct_island(mpol, ntor, ns, is_asym)
 !  Start of executable code
 
 !  Stellarator symmetric quantities.
+      CALL vmec_info_destruct_island
       ALLOCATE(rmnc_i(0:mpol,-ntor:ntor,ns))
       ALLOCATE(zmns_i(0:mpol,-ntor:ntor,ns))
       ALLOCATE(lmns_i(0:mpol,-ntor:ntor,ns + 1))
@@ -324,8 +325,7 @@ SUBROUTINE vmec_info_set_wout(wout_file, ns_in, mpol_in, ntor_in,        &
       ALLOCATE(phipf_i(ns), chipf_i(ns), presf_i(ns), stat=istat)
       CALL vmec_info_spline_oned_array(chipf_vmec, chipf_i, istat)
       CALL vmec_info_spline_oned_array(phipf_vmec, phipf_i, istat)
-      presf_vmec = mu0*presf_vmec
-      CALL vmec_info_spline_oned_array(presf_vmec, presf_i, istat)
+      CALL vmec_info_spline_oned_array(mu0*presf_vmec, presf_i, istat)
 
 !  Pessure should never be negative.
       WHERE (presf_i .lt. 0)