Merge pull request #33 from ORNL-Fusion/sparse

Sparse

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 .and. niter_max .ne. 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