diff --git a/LIBSTELL/Sources/STELLOPT/stellopt_input_mod.f90 b/LIBSTELL/Sources/STELLOPT/stellopt_input_mod.f90 index 05f52766a..727b84573 100644 --- a/LIBSTELL/Sources/STELLOPT/stellopt_input_mod.f90 +++ b/LIBSTELL/Sources/STELLOPT/stellopt_input_mod.f90 @@ -355,6 +355,7 @@ MODULE stellopt_input_mod target_Jstar, sigma_Jstar, NumJstar,& target_helicity, sigma_helicity, helicity,& target_helicity_old, sigma_helicity_old, & + target_quasiiso, sigma_quasiiso, & target_resjac, sigma_resjac, xm_resjac, xn_resjac,& target_separatrix, sigma_separatrix, & r_separatrix, z_separatrix, phi_separatrix, & @@ -917,6 +918,8 @@ SUBROUTINE init_stellopt_input helicity = CMPLX(0.0,0.0) target_helicity_old(:) = 0.0 sigma_helicity_old(:) = bigno + target_quasiiso(:) = 0.0 + sigma_quasiiso(:) = bigno target_resjac(:) = 0.0 sigma_resjac(:) = bigno xn_resjac(:) = 0 @@ -1095,6 +1098,7 @@ SUBROUTINE read_stellopt_input(filename, istat) target_dkes(1) = 0.0; sigma_dkes(1) = bigno target_dkes(2) = 0.0; sigma_dkes(2) = bigno target_helicity(1) = 0.0; sigma_helicity(1) = bigno + target_quasiiso(1) = 0.0; sigma_quasiiso(1) = bigno target_Jstar(1) = 0.0; sigma_Jstar(1) = bigno target_dkes_Erdiff(1) = 0.0; sigma_dkes_Erdiff(1) = bigno target_dkes_alpha(1) = 0.0; sigma_dkes_alpha(1) = bigno @@ -1689,6 +1693,20 @@ SUBROUTINE write_optimum_namelist(iunit,istat) 'SIGMA_HELICITY_OLD(',ik,') = ',sigma_helicity_old(ik) END DO END IF + IF (ANY(sigma_quasiiso < bigno)) THEN + WRITE(iunit,'(A)') '!----------------------------------------------------------------------' + WRITE(iunit,'(A)') '! BOOZER QUASI-ISODYNAMIC METRIC' + WRITE(iunit,'(A)') '!----------------------------------------------------------------------' + n=0 + DO ik = 1,UBOUND(sigma_quasiiso,DIM=1) + IF(sigma_quasiiso(ik) < bigno) n=ik + END DO + DO ik = 1, n + IF (sigma_quasiiso(ik) < bigno) WRITE(iunit,"(2(2X,A,I3.3,A,ES22.12E3))") & + 'TARGET_QUASIISO(',ik,') = ',target_quasiiso(ik), & + 'SIGMA_QUASIISO(',ik,') = ',sigma_quasiiso(ik) + END DO + END IF IF (ANY(sigma_resjac < bigno)) THEN WRITE(iunit,'(A)') '!----------------------------------------------------------------------' WRITE(iunit,'(A)') '! BOOZER Resonant Modes' diff --git a/LIBSTELL/Sources/STELLOPT/stellopt_targets.f90 b/LIBSTELL/Sources/STELLOPT/stellopt_targets.f90 index 547883c0e..2f894ce04 100644 --- a/LIBSTELL/Sources/STELLOPT/stellopt_targets.f90 +++ b/LIBSTELL/Sources/STELLOPT/stellopt_targets.f90 @@ -157,6 +157,7 @@ MODULE stellopt_targets REAL(rprec), DIMENSION(nsd) :: target_magwell, sigma_magwell REAL(rprec), DIMENSION(nsd) :: target_helicity, sigma_helicity REAL(rprec), DIMENSION(nsd) :: target_helicity_old, sigma_helicity_old + REAL(rprec), DIMENSION(nsd) :: target_quasiiso, sigma_quasiiso COMPLEX :: helicity REAL(rprec), DIMENSION(nsd) :: target_resjac, sigma_resjac, & xm_resjac, xn_resjac @@ -273,6 +274,7 @@ MODULE stellopt_targets INTEGER, PARAMETER :: jtarget_neo = 603 INTEGER, PARAMETER :: jtarget_Jstar = 604 INTEGER, PARAMETER :: jtarget_helicity = 605 + INTEGER, PARAMETER :: jtarget_quasiiso = 6051 INTEGER, PARAMETER :: jtarget_resjac = 606 INTEGER, PARAMETER :: jtarget_txport = 607 INTEGER, PARAMETER :: jtarget_dkes = 608 @@ -417,6 +419,8 @@ SUBROUTINE write_targets(iunit,var_num) WRITE(iunit, out_format) 'Trapped Particle J*' CASE(jtarget_helicity) WRITE(iunit, out_format) 'Boozer Spectrum Helicity' + CASE(jtarget_quasiiso) + WRITE(iunit, out_format) 'Boozer Quasi-isodynamic metric' CASE(jtarget_txport) WRITE(iunit, out_format) 'Turbulent Transport' CASE(jtarget_orbit) diff --git a/SHARE/make_macports.inc b/SHARE/make_macports.inc index d2c226953..db58d3db1 100644 --- a/SHARE/make_macports.inc +++ b/SHARE/make_macports.inc @@ -104,7 +104,8 @@ endif ####################################################################### # GENE Options ####################################################################### - LGENE = F +ifneq ("$(wildcard $(GENE_PATH))","") + LGENE = T GENE_INC = -I$(GENE_PATH) GENE_DIR = $(GENE_PATH) LIB_GENE = libgene.a @@ -113,30 +114,42 @@ endif -L$(FFTWHOME)/lib -lfftw3 \ -L$(SLEPC_DIR)/$(PETSC_ARCH)/lib -lslepc \ -L$(PETSC_DIR)/$(PETSC_ARCH)/lib -lpetsc -lX11 +else + LGENE = F +endif ####################################################################### # COILOPT++ Options ####################################################################### - LCOILOPT = F +ifneq ("$(wildcard $(COILOPT_PATH))","") + LCOILOPT = T COILOPT_INC = -I$(COILOPT_PATH) COILOPTPP_DIR = $(COILOPT_PATH) LIB_COILOPTPP = libcoilopt++.a COILOPT_LIB = $(COILOPT_PATH)/$(LIB_COILOPTPP) \ -L$(GSLHOME)/lib -lgsl -lgslcblas -lstdc++ -lmpi_cxx +else + LCOILOPT = F +endif ####################################################################### # TERPSICHORE Options ####################################################################### - LTERPSICHORE= F +ifneq ("$(wildcard $(TERPSICHORE_PATH))","") + LTERPSICHORE= T TERPSICHORE_INC = -I$(TERPSICHORE_PATH) TERPSICHORE_DIR = $(TERPSICHORE_PATH) LIB_TERPSICHORE = libterpsichore.a TERPSICHORE_LIB = $(TERPSICHORE_DIR)/$(LIB_TERPSICHORE) +else + LTERPSICHORE = F +endif ####################################################################### # TRAVIS Options ####################################################################### - LTRAVIS= F +ifneq ("$(wildcard $(TRAVIS_PATH))","") + LTRAVIS = T TRAVIS_DIR = $(TRAVIS_PATH) LIB_TRAVIS = libtravis64_sopt.a LIB_MCONF = libmconf64.a @@ -144,35 +157,49 @@ endif TRAVIS_LIB = $(TRAVIS_DIR)lib/$(LIB_TRAVIS) \ $(TRAVIS_DIR)faddeeva_package/lib/$(LIB_FADDEEVA) \ $(TRAVIS_DIR)magconf/lib/$(LIB_MCONF) -lc++ +else + LTRAVIS = F +endif ####################################################################### # REGCOIL Options ####################################################################### - LREGCOIL= F +ifneq ("$(wildcard $(REGCOIL_PATH))","") + LREGCOIL= T REGCOIL_DIR = $(REGCOIL_PATH) REGCOIL_INC = -I$(REGCOIL_DIR) LIB_REGCOIL = libregcoil.a REGCOIL_LIB = $(REGCOIL_DIR)/$(LIB_REGCOIL) -fopenmp +else + LREGCOIL = F +endif ####################################################################### # SFINCS Options ####################################################################### - - LSFINCS = F +ifneq ("$(wildcard $(SFINCS_PATH))","") + LSFINCS = T SFINCS_DIR = $(SFINCS_PATH) SFINCS_INC = -I$(SFINCS_DIR) LIB_SFINCS = libsfincs.a SFINCS_LIB = $(SFINCS_DIR)/$(LIB_SFINCS) \ -L$(PETSC_DIR)/$(PETSC_ARCH)/lib -lpetsc -lX11 +else + LSFINCS = F +endif ####################################################################### # Available Energy Options ####################################################################### - LAEOPT= F +ifneq ("$(wildcard $(AEOPT_PATH))","") + LAEOPT = T AEOPT_DIR = $(AEOPT_PATH) AEOPT_INC = -I$(AEOPT_DIR) LIB_AEOPT = libtrapAE.a AEOPT_LIB = $(AEOPT_PATH)/$(LIB_AEOPT) +else + LAEOPT = F +endif ####################################################################### # LIBSTELL Shared Options diff --git a/STELLOPTV2/ObjectList b/STELLOPTV2/ObjectList index b2b460980..2f027307f 100644 --- a/STELLOPTV2/ObjectList +++ b/STELLOPTV2/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +chisq_quasiiso.o \ stellopt_read_cws.o \ stellopt_write_header.o \ chisq_dkes_alpha.o \ diff --git a/STELLOPTV2/STELLOPTV2.dep b/STELLOPTV2/STELLOPTV2.dep index 0cf028067..a51389ca0 100644 --- a/STELLOPTV2/STELLOPTV2.dep +++ b/STELLOPTV2/STELLOPTV2.dep @@ -279,6 +279,14 @@ chisq_helicity_ornl.o : \ equil_vals.o \ $(LIB_DIR)/$(LOCTYPE)/read_boozer_mod.o + +chisq_quasiiso.o : \ + stellopt_runtime.o \ + $(LIB_DIR)/$(LOCTYPE)/stellopt_targets.o \ + equil_utils.o \ + $(LIB_DIR)/$(LOCTYPE)/read_boozer_mod.o \ + $(LIB_DIR)/$(LOCTYPE)/vmec_input.o + chisq_resjac.o : \ stellopt_runtime.o \ $(LIB_DIR)/$(LOCTYPE)/stellopt_targets.o \ diff --git a/STELLOPTV2/Sources/Chisq/chisq_quasiiso.f90 b/STELLOPTV2/Sources/Chisq/chisq_quasiiso.f90 new file mode 100644 index 000000000..702cda87c --- /dev/null +++ b/STELLOPTV2/Sources/Chisq/chisq_quasiiso.f90 @@ -0,0 +1,193 @@ +!----------------------------------------------------------------------- +! Subroutine: chisq_quasiiso +! Authors: S. Lazerson (samuel.lazerson@gauss-fusion.com) +! Date: 10/08/2024 +! Description: This chisquared functional replicates the QI +! functional described in: +! https://doi.org/10.1017/S002237782300065X +! and +! https://doi.org/10.1103/PRXEnergy.3.023010 +!----------------------------------------------------------------------- + SUBROUTINE chisq_quasiiso(target,sigma,niter,iflag) +!----------------------------------------------------------------------- +! Libraries +!----------------------------------------------------------------------- + USE stellopt_runtime + USE stellopt_targets + USE equil_utils, ONLY: get_equil_iota + USE read_boozer_mod, ONLY: bmnc_b, ixn_b, ixm_b, mnboz_b, ns_b, nfp_b + USE vmec_input, ONLY: mpol, ntor +! USE stel_kinds, ONLY: rprec + +!----------------------------------------------------------------------- +! Input/Output Variables +! +!----------------------------------------------------------------------- + IMPLICIT NONE + REAL(rprec), INTENT(in) :: target(nsd) + REAL(rprec), INTENT(in) :: sigma(nsd) + INTEGER, INTENT(in) :: niter + INTEGER, INTENT(inout) :: iflag + +!----------------------------------------------------------------------- +! Local Variables +! +!----------------------------------------------------------------------- + INTEGER :: ik, mn, ier, i, l, m, lmin, lmax, i1, i2, lh, lr + INTEGER, PARAMETER :: nalpha = 65 ! must be odd + INTEGER, PARAMETER :: ntheta0 = 5 + INTEGER, PARAMETER :: nlambda = 4 + REAL(rprec) :: s_temp, iota0, phi, theta, deltaphi, ftemp, norm, & + Bmax, Bmin, Bmir, lambda, theta0, phimin,phimax + REAL(rprec), DIMENSION(:), ALLOCATABLE :: modb, modbs, dl, & + integral_A, integral_B + REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: J_I, J_C + + LOGICAL, PARAMETER :: lwrite_out = .False. + + +!----------------------------------------------------------------------- +! BEGIN SUBROUTINE +!----------------------------------------------------------------------- + IF (iflag < 0) RETURN + ik = COUNT(sigma < bigno) + IF (iflag == 1) WRITE(iunit_out,'(A,2(2X,I3.3))') 'QUASIISO ',ik,4 + IF (iflag == 1) WRITE(iunit_out,'(A)') 'TARGET SIGMA VAL K' + IF (niter >= 0) THEN + ! Allocate helpers along field line + ALLOCATE(modb(nalpha),modbs(nalpha),dl(nalpha),integral_A(nalpha),& + integral_B(nalpha)) + ! Allocate J arrays + ALLOCATE(J_C(nlambda,ntheta0),J_I(nlambda,ntheta0)) + ! Loop over surfaces + DO ik = 1, nsd + IF (ABS(sigma(ik)) >= bigno) CYCLE + ! Get iota for the surface + s_temp = DBLE(ik)/DBLE(ns_b) + ier = 0 + CALL get_equil_iota(s_temp,iota0,ier) + ! Loop over inital values of theta0 + J_C = 0.0; J_I = 0.0 + DO i = 1, ntheta0 + ! Calculate the initial distance to follow a field line + deltaphi = pi2/iota0 + modb = 0.0 + modbs = 0.0 + theta0 = pi2*(i-1)/ntheta0 + DO l = 1, nalpha + phi = deltaphi*(l-1)/(nalpha-1) + theta = theta0 + iota0*phi + DO mn = 1, mnboz_b + modb(l) = modb(l)+bmnc_b(mn,ik)*COS(DBLE(ixm_b(mn))*theta+DBLE(ixn_b(mn))*phi/nfp_b) + END DO + END DO + IF (lwrite_out) WRITE(320,*) modb + ! Find min/max modB values + lmin = MINLOC(modb,DIM=1) + lmax = MAXLOC(modb,DIM=1) + ! Recalculate length of fieldline + phimin = deltaphi*(lmin-1)/(nalpha-1) + phimax = deltaphi*(lmax-1)/(nalpha-1) + deltaphi = ABS(phimax-phimin) + ! Now recalculate modb centered around lmin + modb = 0 + DO l = 1, nalpha + phi = phimin-deltaphi + 2*deltaphi*(l-1)/(nalpha-1) + !phi = deltaphi*(l-lmin-1)/(nalpha-1) - deltaphi/2 + theta = theta0 + iota0*phi + DO mn = 1, mnboz_b + modb(l) = modb(l)+bmnc_b(mn,ik)*COS(DBLE(ixm_b(mn))*theta+DBLE(ixn_b(mn))*phi/nfp_b) + END DO + END DO + IF (lwrite_out) WRITE(321,*) modb + ! Find the Bmin and half point + lh = FLOOR(nalpha*0.5)+1 + lmin = MINLOC(modb,DIM=1) + lmin = MIN(MAX(lmin,2),nalpha-1) + ! Now shift lmin to the middle + modb = CSHIFT(modb,lmin-lh) !BOOSH + ! Now recalc lh as lmin + lmin = MINLOC(modb,DIM=1) + lmin = MIN(MAX(lmin,2),nalpha-1) + ! Squash the array + modbs = modb + DO l = lmin, 1, -1 + IF (modbs(l)Bmir) + 1 + i2 = COUNT(modbs(lmin:nalpha)0 Well)') + elif (plot_name == 'QUASIISO_evolution'): + x = self.stel_data.QUASIISO_K + y = self.stel_data.QUASIISO_VAL + t = self.stel_data.QUASIISO_TARGET + d = self.stel_data.QUASIISO_SIGMA + self.ax2.errorbar(x[0,:],t[0,:],yerr=d[0,:],fmt='ok',fillstyle='none',label='Target') + self.ax2.plot(x[0,:],y[0,:],'o',fillstyle='none',label='Initial',color='red') + for i in range(1,niter-1,1): + self.ax2.plot(x[i,:],y[i,:],'.k',fillstyle='none') + self.ax2.plot(x[niter-1,:],y[niter-1,:],'o',fillstyle='none',label='Final',color='green') + self.ax2.set_xlabel('Radial Grid') + self.ax2.set_ylabel('QI Metric') + self.ax2.set_title('Quasi-isodynamic Metric') elif (plot_name == 'CURVATURE_P2'): #x = self.stel_data.TXPORT_S #y = self.stel_data.CURVATURE_P2_P2 diff --git a/pySTEL/libstell/stellopt.py b/pySTEL/libstell/stellopt.py index c3648d68a..2762347e0 100644 --- a/pySTEL/libstell/stellopt.py +++ b/pySTEL/libstell/stellopt.py @@ -33,7 +33,7 @@ def __init__(self): 'NEO', 'DKES', 'DKES_ERDIFF', 'DKES_ALPHA', 'TXPORT', \ 'ORBIT', 'HELICITY', 'HELICITY_FULL', 'JSTAR', 'RESJAC', \ 'COIL_BNORM', 'REGCOIL_CHI2_B', 'CURVATURE_P2', 'GAMMA_C', \ - 'KINK'] + 'KINK', 'QUASIISO'] def read_stellopt_map(self,filename='map.dat'): """Reads a STELLOPT MAP output file