From 468964738e8cf4871c812ab75d46ca1d6c367e28 Mon Sep 17 00:00:00 2001 From: cianciosa Date: Wed, 21 Feb 2024 13:41:44 -0500 Subject: [PATCH] Fix issues using the SIESTA run context in V3FIT. --- Sources/evolution.f90 | 2 +- Sources/gmres.f90 | 54 +++++++++++++++++++++---------------- Sources/grid_extension.f90 | 5 ++-- Sources/hessian.f90 | 9 ++++--- Sources/metrics.f90 | 33 +++++++++++------------ Sources/nscalingtools.f90 | 2 +- Sources/pchelms.f90 | 4 +-- Sources/perturbation.f90 | 2 +- Sources/quantities.f90 | 49 ++++++++++++++++++++------------- Sources/siesta_namelist.f90 | 2 +- Sources/siesta_run.f90 | 43 +++++++++++++++++++++++++++++ Sources/vmec_info.f90 | 4 +-- 12 files changed, 136 insertions(+), 73 deletions(-) diff --git a/Sources/evolution.f90 b/Sources/evolution.f90 index b74069c..2b498a8 100644 --- a/Sources/evolution.f90 +++ b/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) diff --git a/Sources/gmres.f90 b/Sources/gmres.f90 index 2d99287..0d0a56f 100644 --- a/Sources/gmres.f90 +++ b/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) diff --git a/Sources/grid_extension.f90 b/Sources/grid_extension.f90 index dc425c6..0566bc4 100644 --- a/Sources/grid_extension.f90 +++ b/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 diff --git a/Sources/hessian.f90 b/Sources/hessian.f90 index 3b51f9f..f2b5fd0 100644 --- a/Sources/hessian.f90 +++ b/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 diff --git a/Sources/metrics.f90 b/Sources/metrics.f90 index 7faac84..f051955 100644 --- a/Sources/metrics.f90 +++ b/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 diff --git a/Sources/nscalingtools.f90 b/Sources/nscalingtools.f90 index ead0d77..987de88 100644 --- a/Sources/nscalingtools.f90 +++ b/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. !-------------------------------------------------------------------------- diff --git a/Sources/pchelms.f90 b/Sources/pchelms.f90 index dbd4597..e59ac53 100644 --- a/Sources/pchelms.f90 +++ b/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 (*,*) ' = ', & twopi*asubvmncf(1,ntor + 1,ns)*fourier_context%orthonorm(0,0) WRITE (*,*) ' = ', & diff --git a/Sources/perturbation.f90 b/Sources/perturbation.f90 index e4470dc..9ba6151 100644 --- a/Sources/perturbation.f90 +++ b/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)) diff --git a/Sources/quantities.f90 b/Sources/quantities.f90 index 5647d0f..ba168ef 100644 --- a/Sources/quantities.f90 +++ b/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 diff --git a/Sources/siesta_namelist.f90 b/Sources/siesta_namelist.f90 index a4d9c31..8893222 100644 --- a/Sources/siesta_namelist.f90 +++ b/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 diff --git a/Sources/siesta_run.f90 b/Sources/siesta_run.f90 index 289ed5c..60ca28a 100644 --- a/Sources/siesta_run.f90 +++ b/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 diff --git a/Sources/vmec_info.f90 b/Sources/vmec_info.f90 index 237582e..f8591f9 100644 --- a/Sources/vmec_info.f90 +++ b/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)