Fix some allocation error when running as an embeded code. Force heli…

…cal pertubation when run in memory and force the restart file to be written after the helpert is applied.

Sources/evolution.f90

@@ -70,6 +70,8 @@ MODULE evolution
       REAL (dp), DIMENSION(:), ALLOCATABLE :: xcdot
 !>  Change in force squared residule.
       REAL (dp)                            :: fsq_last
+!>  Pertubation was added.
+      LOGICAL                              :: pert_added
 
       CONTAINS
 
@@ -87,7 +89,7 @@ MODULE evolution
 !>  @param[in] ftol      Force residual tolerance.
 !-------------------------------------------------------------------------------
       SUBROUTINE converge_diagonal(wout_file, ftol)
-      USE siesta_namelist, ONLY: ladd_pert, l_output_alliter
+      USE siesta_namelist, ONLY: l_output_alliter
       USE perturbation, ONLY: add_perturb, niter_max
       USE dumping_mod, ONLY: write_output
 
@@ -138,17 +140,19 @@ SUBROUTINE converge_diagonal(wout_file, ftol)
             CALL write_output(wout_file, niter, .false.)
          END IF
 
-         l_iterate = fsq_ratio1 .le. 0.5_dp    .and.                           &
-                     fsq_total1 .gt. fsq_block .and.                           &
-                     niter      .lt. niter_max
+         l_iterate = (fsq_ratio1 .le. 0.5_dp    .and.                          &
+                      fsq_total1 .gt. fsq_block .and.                          &
+                      niter      .lt. niter_max) .or.                          &
+                     niter .eq. 1
 
-         IF (ladd_pert .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)
             CALL second0(t2)
             diag_add_pert_time=diag_add_pert_time+(t2-t1)
-            ladd_pert = .FALSE.
+            pert_added = .TRUE.
+            fsq_min = 1.0E20_dp
          END IF
          niter = niter + 1
       END DO
@@ -166,7 +170,7 @@ SUBROUTINE converge_diagonal(wout_file, ftol)
 !>  @param[in] ftol      Force residual tolerance.
 !-------------------------------------------------------------------------------
       SUBROUTINE converge_blocks(wout_file, ftol)
-      USE siesta_namelist, ONLY: l_output_alliter, ladd_pert
+      USE siesta_namelist, ONLY: l_output_alliter
       USE perturbation, ONLY: add_perturb, niter_max
       USE dumping_mod, ONLY: write_output
 
@@ -224,24 +228,25 @@ SUBROUTINE converge_blocks(wout_file, ftol)
          IF (nprecon .gt. 2 .and. ABS(1 - fsq_ratio1) .lt. 1.E-2_dp) THEN
             levm_scale = levm_scale/3
             nrow = nrow + 1
-            IF (nrow.EQ.3 .and. l_iterate) THEN
+!            IF (nrow.EQ.3 .and. l_iterate) THEN
 !               IF (iam .eq. 0 .and. lverbose) THEN
 !                  WRITE(*,1001), fsq_ratio1
 !               END IF
 !               l_iterate = .false.
-            END IF
-         ELSE 
+!            END IF
+         ELSE
             nrow = 0
          END IF
 
 !  In case we didn't add it in diag loop.
-         IF (ladd_pert .and. fsq_total1 .lt. 100*fsq_prec) THEN
+         IF (.not.pert_added .and. fsq_total1 .lt. 100*fsq_prec) THEN
            l_init_state = .true.
            CALL second0(t1)
            CALL add_perturb(xc, getwmhd)
            CALL second0(t2)
            block_add_pert_time = block_add_pert_time + (t2 - t1)
-           ladd_pert = .FALSE.
+           pert_added = .TRUE.
+           fsq_min = 1.0E20_dp
 !  To output right after pert is applied, set l_iterate = .false. here.
          END IF
 
@@ -316,15 +321,19 @@ SUBROUTINE init_evolution
       n1 = ns*mnmax
       neqs = ndims*n1
 
-      ALLOCATE(xc(neqs), col_scale(0:mpol,-ntor:ntor,ndims,ns), stat=istat)
-      CALL ASSERT(istat.eq.0,'Allocate xc failed!')
+      IF (.not.ALLOCATED(xc)) THEN
+         ALLOCATE(xc(neqs), col_scale(0:mpol,-ntor:ntor,ndims,ns), stat=istat)
+         CALL ASSERT(istat.eq.0,'Allocate xc failed!')
+      END IF
       xc = 0
       col_scale = 1
       CALL init_ptrs(xc, jvsupsmncf, jvsupumnsf, jvsupvmnsf,                   &
                          jvsupsmnsf, jvsupumncf, jvsupvmncf)
 
-      ALLOCATE(gc(neqs), gc_sup(neqs), stat=istat)
-      CALL ASSERT(istat.eq.0,'Allocate gc failed!')
+      IF (.not.ALLOCATED(gc)) THEN
+         ALLOCATE(gc(neqs), gc_sup(neqs), stat=istat)
+         CALL ASSERT(istat.eq.0,'Allocate gc failed!')
+      END IF
       gc = 0
       CALL init_ptrs(gc_sup, fsupsmncf, fsupumnsf, fsupvmnsf,                  &
                              fsupsmnsf, fsupumncf, fsupvmncf)

Sources/hessian.f90

@@ -991,7 +991,7 @@ SUBROUTINE Compute_Hessian_Blocks_Thomas (xc, gc, func)
          CALL CheckEigenvalues_Serial(ns, mblk_size)
          CALL CheckConditionNumber_Serial(ns,mblk_size,anorm,rcond,info)
          CALL second0(toff)                  
-         IF (INFO .EQ. 0) THEN
+         IF (INFO .EQ. 0 .and. lverbose) THEN
             PRINT '(1x,3(a,1p,e12.3))','RCOND = ', rcond,               &
             ' ||A|| = ', ANORM,' TIME: ', toff-ton
          END IF

Sources/quantities.f90

@@ -253,8 +253,10 @@ MODULE quantities
 !*******************************************************************************
 !-------------------------------------------------------------------------------
 !>  @brief Intialize quantities.
+!>
+!>  @param[in] restarted Flag to signal that we are in restart mode.
 !-------------------------------------------------------------------------------
-      SUBROUTINE init_quantities
+      SUBROUTINE init_quantities(restarted)
       USE descriptor_mod, ONLY: LSCALAPACK
       USE island_params, ONLY: fourier_context
       USE fourier, ONLY: f_none, f_cos, f_sin, f_sum, m0
@@ -264,10 +266,12 @@ SUBROUTINE init_quantities
       USE prof_mod, ONLY: SetUpScalingAllGather
 #endif
       USE blocktridiagonalsolver_s, ONLY: Initialize
-      USE siesta_namelist, ONLY: lrestart
 
       IMPLICIT NONE
 
+!  Declare Arguments
+      LOGICAL                                  :: restarted
+
 !  Local variables.
       INTEGER                                  :: istat
       INTEGER                                  :: l
@@ -287,7 +291,7 @@ SUBROUTINE init_quantities
          endglobrow = ns
       END IF
 
-      IF (.not.lrestart) THEN
+      IF (.not.restarted) THEN
 !  Allocate and initialize dependent variables
          CALL alloc_quantities
 

Sources/restart_mod.f90

@@ -196,17 +196,20 @@ FUNCTION restart_read(restart_ext, wout_file, mpolin, ntorin, nfpin,     &
                             jbsupumnsh, jbsupumnch,                            &
                             jbsupvmnsh, jbsupvmnch,                            &
                             jpmnsh,     jpmnch,                                &
-                            b_factor, p_factor, alloc_quantities
+                            b_factor, p_factor, alloc_quantities,              &
+                            dealloc_quantities
       USE island_params, ONLY: chipf => chipf_i, phipf => phipf_i,             &
      &                         wb => wb_i, wp => wp_i, nfp_i, gamma,           &
      &                         gnorm => gnorm_i, rmajor => rmajor_i,           &
      &                         fourier_context, nu_i, nv_i
       USE vmec_info, ONLY: rmnc => rmnc_i, zmns => zmns_i,                     &
      &                     rmns => rmns_i, zmnc => zmnc_i,                     &
-     &                     vmec_info_construct_island
+     &                     vmec_info_construct_island,                         &
+     &                     vmec_info_destruct_island
       USE metrics, ONLY: LoadGrid
       USE fourier, ONLY: m0, m1, fourier_class
       USE stel_constants, ONLY: one
+      USE siesta_namelist, ONLY: lrestart
 
       IMPLICIT NONE
 
@@ -237,6 +240,7 @@ FUNCTION restart_read(restart_ext, wout_file, mpolin, ntorin, nfpin,     &
 !  Start of executable code
       restart_read = 1
 
+      CALL dealloc_quantities
       CALL alloc_quantities
 
       filename = 'siesta_' // TRIM(restart_ext) // '.nc'
@@ -272,13 +276,18 @@ FUNCTION restart_read(restart_ext, wout_file, mpolin, ntorin, nfpin,     &
 !  The namelist input file may turn the asymmetric terms on and off.
       ALLOCATE(tempmn_r(0:mpol,-ntor:ntor,ns))
       ALLOCATE(temp_r(ns))
+      CALL vmec_info_destruct_island
       CALL vmec_info_construct_island(mpolin, ntorin, nsin, lasym)
 
-      ALLOCATE(chipf(nsin))
+      IF (.not.ALLOCATED(chipf)) THEN
+         ALLOCATE(chipf(nsin))
+      END IF
       CALL cdf_read(ncid, vn_chipf, temp_r)
       CALL interpit_1d(temp_r, chipf, ns, nsin, .false., 1)
 
-      ALLOCATE(phipf(nsin))
+      IF (.not.ALLOCATED(phipf)) THEN
+         ALLOCATE(phipf(nsin))
+      END IF
       CALL cdf_read(ncid, vn_phipf, temp_r)
       CALL interpit_1d(temp_r, phipf, ns, nsin, .false., 1)
 

Sources/siesta_run.f90

@@ -435,7 +435,7 @@ SUBROUTINE siesta_run_set_vmec_temp(this, load_wout)
       END IF
 
       CALL init_metric_elements()
-      CALL init_quantities !Initializes BCYCLIC
+      CALL init_quantities(.false.) !Initializes BCYCLIC
       CALL init_evolution !neqs is set here
       CALL initRemap(mpol, ntor, ns, nprocs, iam)
       CALL InitHess
@@ -491,7 +491,7 @@ SUBROUTINE siesta_run_set_restart(this)
                              ns, ntor_modes(-ntorin:ntorin))
 
       CALL init_metric_elements()
-      CALL init_quantities !Initializes BCYCLIC
+      CALL init_quantities(.true.) !Initializes BCYCLIC
       CALL init_evolution !neqs is set here
       CALL initRemap(mpolin, ntorin, ns, nprocs, iam)
       CALL InitHess
@@ -582,10 +582,11 @@ FUNCTION siesta_run_get_1d(this, param_name, index)
 !>  @param[inout] this A @ref siesta_run_class instance.
 !-------------------------------------------------------------------------------
       SUBROUTINE siesta_run_converge(this)
-      USE evolution, ONLY: converge_diagonal, converge_blocks
+      USE evolution, ONLY: converge_diagonal, converge_blocks, pert_added
       USE descriptor_mod, ONLY: DIAGONALDONE
-      USE siesta_namelist, ONLY: ftol, wout_file
+      USE siesta_namelist, ONLY: ftol, wout_file, ladd_pert
       USE utilities, ONLY: test_utilities
+      USE shared_data, ONLY: xc
 
       IMPLICIT NONE
 
@@ -596,12 +597,15 @@ SUBROUTINE siesta_run_converge(this)
       LOGICAL                                 :: passed
 
 !  Start of executable code
+      pert_added = .not.ladd_pert
+
       IF (test_utilities()) THEN
          WRITE (*,*) 'Failed Diverence Test.'
 !         STOP
       END IF
 
 !  Converge initial residues with diagonal preconditioner
+      xc = 0
       DIAGONALDONE = .false.
       CALL converge_diagonal(wout_file, ftol)
       DIAGONALDONE = .true.

Sources/vmec_info.f90

@@ -219,11 +219,13 @@ SUBROUTINE vmec_info_destruct_island
 !  Start of executable code
 
 !  Stellarator symmetric quantities.
-      DEALLOCATE(rmnc_i)
-      DEALLOCATE(zmns_i)
-      DEALLOCATE(lmns_i)
-      DEALLOCATE(jcurrumnc)
-      DEALLOCATE(jcurrvmnc)
+      IF (ALLOCATED(rmnc_i)) THEN
+         DEALLOCATE(rmnc_i)
+         DEALLOCATE(zmns_i)
+         DEALLOCATE(lmns_i)
+         DEALLOCATE(jcurrumnc)
+         DEALLOCATE(jcurrvmnc)
+      END IF
 
 !  Asymmetric quantities.
       IF (ALLOCATED(rmns_i)) THEN