diff --git a/docs/examples/DIIID_resistive_example/rdcon.in b/docs/examples/DIIID_resistive_example/rdcon.in index eef0ddbf1..36aebd085 100644 --- a/docs/examples/DIIID_resistive_example/rdcon.in +++ b/docs/examples/DIIID_resistive_example/rdcon.in @@ -75,8 +75,8 @@ bin_bal1=f ! Binary output for bal_flag poloidal functions out_bal2=f ! Ascii output for bal_flag functions bin_bal2=f ! Binary output for bal_flag functions - MRE_flag=t ! Netcdf output for various modified Rutherford equation terms (defined in mercier.f) - geom_flag=t ! Netcdf output for more flux-surface averaged terms (defined in mercier.f). To use, must set MRE_flag = t + MRE_flag=f ! Netcdf output for various modified Rutherford equation terms (defined in mercier.f) + geom_flag=f ! Netcdf output for more flux-surface averaged terms (defined in mercier.f). To use, must set MRE_flag = t / &UA_DIAGNOSE_LIST diff --git a/equil/direct.f b/equil/direct.f index b58e3de29..675b8f534 100644 --- a/equil/direct.f +++ b/equil/direct.f @@ -74,7 +74,7 @@ SUBROUTINE direct_run sq_in%fs(:,4)=SQRT(sq_in%xs) sq_in%name=" sq " sq_in%title=(/"psifac"," f ","mu0 p "," q "," rho "/) - CALL spline_fit(sq_in,"extrap") + CALL spline_fit_hermite(sq_in,"extrap",sq_in_slopes) psi_in%xs=rmin+(/(ir,ir=0,psi_in%mx)/)*(rmax-rmin)/psi_in%mx psi_in%ys=zmin+(/(iz,iz=0,psi_in%my)/)*(zmax-zmin)/psi_in%my CALL bicube_fit(psi_in,"extrap","extrap") @@ -179,6 +179,8 @@ SUBROUTINE direct_run c----------------------------------------------------------------------- sq%fs(ipsi,1)=bf%f*twopi sq%fs(ipsi,2)=bf%p + sq%fs1(ipsi,1)=bf%f1*twopi + sq%fs1(ipsi,2)=bf%p1 sq%fs(ipsi,3)=y_out(istep,1)*twopi*psio sq%fs(ipsi,4)=y_out(istep,3)*bf%f/twopi CALL spline_dealloc(ff) @@ -191,9 +193,8 @@ SUBROUTINE direct_run c----------------------------------------------------------------------- c fit surface quantities to cubic splines. c----------------------------------------------------------------------- - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) sq%name=" sq " - sq%title=(/" psi "," f "," p "," q "/) q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0) c----------------------------------------------------------------------- c revise q profile. @@ -205,10 +206,11 @@ SUBROUTINE direct_run DO ipsi=0,mpsi ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac + sq%fs1(ipsi,1)=sq%fs1(ipsi,1)/ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac ENDDO - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) ENDIF qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(1-sq%xs(mpsi)) c----------------------------------------------------------------------- diff --git a/equil/equil.f b/equil/equil.f index be7456d46..fbc12fe8f 100644 --- a/equil/equil.f +++ b/equil/equil.f @@ -39,7 +39,7 @@ SUBROUTINE equil_loadnamelists(op_psihigh, op_psilow) $ newq0,psihigh,psilow,input_only,jac_type,power_bp,power_r, $ power_b,jac_method,convert_type,power_flag,ns1, $ sp_pfac,sp_nx,sp_dx1,sp_dx2,use_galgrid,etol,enstep, - $ use_classic_splines + $ use_classic_splines,profile_source NAMELIST/equil_output/bin_2d,bin_eq_1d,bin_eq_2d,out_2d,out_eq_1d, $ out_eq_2d,bin_fl,out_fl,interp,gse_flag,dump_flag,verbose, $ out_ahg2msc @@ -150,6 +150,7 @@ SUBROUTINE equil_read(unit, op_psihigh, op_psilow) c----------------------------------------------------------------------- c read equilibrium data and diagnose. c----------------------------------------------------------------------- + sq_in_slopes=.FALSE. SELECT CASE(TRIM(eq_type)) CASE("fluxgrid") CALL read_eq_fluxgrid diff --git a/equil/equil_out.f b/equil/equil_out.f index 1069179be..bb4fb90d8 100644 --- a/equil/equil_out.f +++ b/equil/equil_out.f @@ -97,7 +97,9 @@ SUBROUTINE equil_out_diagnose(sq_flag,unit) sq_out%fs(:,2)=sq%fs(:,2) sq_out%fs(:,3)=sq%fs(:,4) sq_out%fs(:,4)=SQRT(sq_out%xs) - CALL spline_fit(sq_out,"extrap") + sq_out%fs1(:,1)=sq%fs1(:,1)/twopi + sq_out%fs1(:,2)=sq%fs1(:,2) + CALL spline_fit_hermite(sq_out,"extrap",sq_in_slopes) sq_out%name=" sq " sq_out%title=(/"psifac"," f ","mu0 p "," q "," rho "/) IF(bin_eq_1d)CALL bin_open(bin_2d_unit,"sq_out.bin","UNKNOWN", @@ -845,6 +847,11 @@ SUBROUTINE equil_out_dump WRITE(dump_unit)rzphi%name WRITE(dump_unit)rzphi%periodic c----------------------------------------------------------------------- +c write f, p knot derivatives and eqfun; optional for read_eq_dump. +c----------------------------------------------------------------------- + WRITE(dump_unit)sq_in_slopes(1:2),sq%fs1(:,1:2) + IF(eqfun%allocated)WRITE(dump_unit)eqfun%fs +c----------------------------------------------------------------------- c close dump file. c----------------------------------------------------------------------- CALL bin_close(dump_unit) diff --git a/equil/global.f b/equil/global.f index b57137f02..53594f553 100644 --- a/equil/global.f +++ b/equil/global.f @@ -16,6 +16,8 @@ MODULE global_mod CHARACTER(16) :: eq_type="fluxgrid" CHARACTER(128) :: eq_filename="fluxgrid.dat" + CHARACTER(16) :: profile_source="values" + LOGICAL, DIMENSION(4) :: sq_in_slopes=.FALSE. LOGICAL :: direct_flag=.FALSE. LOGICAL :: out_eq_1d=.FALSE.,bin_eq_1d=.FALSE. diff --git a/equil/inverse.f b/equil/inverse.f index 922b027e6..ac19224bb 100644 --- a/equil/inverse.f +++ b/equil/inverse.f @@ -54,7 +54,7 @@ SUBROUTINE inverse_run sq_in%fs(:,4)=SQRT(sq_in%xs) sq_in%name=" sq " sq_in%title=(/"psifac"," f ","mu0 p "," q "," rho "/) - CALL spline_fit(sq_in,"extrap") + CALL spline_fit_hermite(sq_in,"extrap",sq_in_slopes) rz_in%xs=sq_in%xs rz_in%ys=(/(itheta,itheta=0,rz_in%my)/)/REAL(rz_in%my,r8) rz_in%name=" rz " @@ -180,7 +180,7 @@ SUBROUTINE inverse_run CALL spline_alloc(spl,mtheta,5) DO ipsi=0,mpsi psifac=rzphi%xs(ipsi) - CALL spline_eval(sq_in,psifac,0) + CALL spline_eval(sq_in,psifac,1) spl%xs=rzphi%ys c----------------------------------------------------------------------- c interpolate to new surface. @@ -232,10 +232,12 @@ SUBROUTINE inverse_run c----------------------------------------------------------------------- sq%fs(ipsi,1)=sq_in%f(1)*twopi sq%fs(ipsi,2)=sq_in%f(2) + sq%fs1(ipsi,1)=sq_in%f1(1)*twopi + sq%fs1(ipsi,2)=sq_in%f1(2) sq%fs(ipsi,3)=spl%fsi(mtheta,3)*twopi*pi sq%fs(ipsi,4)=spl%fsi(mtheta,4)*sq%fs(ipsi,1)/(2*twopi*psio) ENDDO - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0) c----------------------------------------------------------------------- c revise q profile. @@ -247,10 +249,11 @@ SUBROUTINE inverse_run DO ipsi=0,mpsi ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac + sq%fs1(ipsi,1)=sq%fs1(ipsi,1)/ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac ENDDO - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) ENDIF qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(one-sq%xs(mpsi)) c----------------------------------------------------------------------- diff --git a/equil/local.f b/equil/local.f index b4e25380d..0cd9c436c 100644 --- a/equil/local.f +++ b/equil/local.f @@ -13,6 +13,7 @@ c 5. ascii_close. c 6. program_stop. c 7. floored_log. +c 8. read_var_len. c----------------------------------------------------------------------- c subprogram 0. local. c module declarations. @@ -259,4 +260,46 @@ FUNCTION floored_log(u) RESULT(ulog) c----------------------------------------------------------------------- RETURN END FUNCTION floored_log +c----------------------------------------------------------------------- +c subprogram 8. read_var_len. +c allocates variable length arrays from set-length read-in arrays. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE read_var_len(x_in, x_out) + + IMPLICIT NONE + + REAL(r8), DIMENSION(:), INTENT(IN) :: x_in ! Indexes from 1... + REAL(r8), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: x_out + INTEGER :: n, i + REAL(r8) :: default_value +c----------------------------------------------------------------------- +c assume last value of x_in is a sentinal +c----------------------------------------------------------------------- + default_value = x_in(SIZE(x_in)) +c----------------------------------------------------------------------- +c count leading non-sentinel values in x_in. +c----------------------------------------------------------------------- + n = 0 + DO i = 1, SIZE(x_in) - 1 + IF (x_in(i) /= default_value) THEN + n = n + 1 + ELSE + EXIT + END IF + END DO +c----------------------------------------------------------------------- +c allocate output array and copy if values exist. +c----------------------------------------------------------------------- + IF (n > 0) THEN + ALLOCATE(x_out(n)) + x_out(1:n) = x_in(1:n) + END IF +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE read_var_len END MODULE local_mod diff --git a/equil/read_eq.f b/equil/read_eq.f index 3ad187791..943973687 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -37,6 +37,8 @@ c 29. read_eq_chease4. c 30. read_eq_marklin_direct. c 31. read_eq_marklin_inverse. +c 32. read_eq_profiles_from_derivatives. +c 33. read_eq_profiles_warn. c----------------------------------------------------------------------- c subprogram 0. read_eq_mod. c module declarations. @@ -625,6 +627,7 @@ SUBROUTINE read_eq_efit REAL(r8) :: bcentr,cpasma,rgrid,rmaxis,rzero,ssibry1, $ ssibry2,ssimag1,ssimag2,xdim,xdum,zdim,zmaxis,zmid REAL(r8) :: tmp + REAL(r8), DIMENSION(:), ALLOCATABLE :: ffprim,pprime CHARACTER(LEN=7) :: shotstr, timestr c----------------------------------------------------------------------- c read equilibrium data. @@ -636,14 +639,15 @@ SUBROUTINE read_eq_efit READ(in_unit,'(5e16.9)')cpasma,ssimag2,xdum,rmaxis,xdum READ(in_unit,'(5e16.9)')zmaxis,xdum,ssibry2,xdum,xdum CALL spline_alloc(sq_in,nw-1,4) + ALLOCATE(ffprim(0:nw-1),pprime(0:nw-1)) READ(in_unit,'(5e16.9)')(sq_in%fs(i,1),i=0,nw-1) READ(in_unit,'(5e16.9)')(sq_in%fs(i,2),i=0,nw-1) - READ(in_unit,'(5e16.9)')(sq_in%fs(i,3),i=0,nw-1) - READ(in_unit,'(5e16.9)')(sq_in%fs(i,3),i=0,nw-1) + READ(in_unit,'(5e16.9)')(ffprim(i),i=0,nw-1) + READ(in_unit,'(5e16.9)')(pprime(i),i=0,nw-1) CALL bicube_alloc(psi_in,nw-1,nh-1,1) READ(in_unit,'(5e16.9)')((psi_in%fs(i,j,1),i=0,nw-1),j=0,nh-1) READ(in_unit,'(5e16.9)',iostat=ios)(sq_in%fs(i,3),i=0,nw-1) - IF(ios /= 0)sq_in%fs(i,3)=0 + IF(ios /= 0)sq_in%fs(:,3)=0 CALL ascii_close(in_unit) c----------------------------------------------------------------------- c translate to internal quantities. @@ -667,6 +671,9 @@ SUBROUTINE read_eq_efit sq_in%xs=(/(ia,ia=0,ma)/)/dfloat(ma) sq_in%fs(:,1)=ABS(sq_in%fs(:,1)) sq_in%fs(:,2)=MAX(sq_in%fs(:,2)*mu0,zero) + CALL read_eq_profiles_from_derivatives(ffprim*psio, + $ pprime*mu0*psio) + DEALLOCATE(ffprim,pprime) c----------------------------------------------------------------------- c copy and convert 2D quantities. c----------------------------------------------------------------------- @@ -808,8 +815,8 @@ END SUBROUTINE read_eq_ldp_d c----------------------------------------------------------------------- SUBROUTINE read_eq_ldp_i - INTEGER :: mx,my - REAL(r8), DIMENSION(:), POINTER :: psi,f,p,q + INTEGER :: mx,my,ios1,ios2 + REAL(r8), DIMENSION(:), POINTER :: psi,f,p,q,ffp,pp REAL(r8), DIMENSION(:,:), POINTER :: r,z c----------------------------------------------------------------------- c open data file, read scalar data. @@ -823,9 +830,9 @@ SUBROUTINE read_eq_ldp_i my=my-1 mx=mx-1 ALLOCATE(psi(0:mx),f(0:mx),p(0:mx),q(0:mx), - $ r(0:my,0:mx),z(0:my,0:mx)) + $ r(0:my,0:mx),z(0:my,0:mx),ffp(0:mx),pp(0:mx)) c----------------------------------------------------------------------- -c read binary data. +c read binary data; f*df/dpsi and dp/dpsi records are optional. c----------------------------------------------------------------------- READ(in_unit)psi READ(in_unit)f @@ -833,6 +840,8 @@ SUBROUTINE read_eq_ldp_i READ(in_unit)q READ(in_unit)r READ(in_unit)z + READ(in_unit,IOSTAT=ios1)ffp + READ(in_unit,IOSTAT=ios2)pp CALL bin_close(in_unit) c----------------------------------------------------------------------- c copy and revise 1D arrays. @@ -843,6 +852,11 @@ SUBROUTINE read_eq_ldp_i sq_in%fs(:,1)=f sq_in%fs(:,2)=p*mu0 sq_in%fs(:,3)=q + IF(ios1 == 0 .AND. ios2 == 0)THEN + CALL read_eq_profiles_from_derivatives(ffp*psio,pp*mu0*psio) + ELSEIF(TRIM(profile_source) /= "values")THEN + CALL read_eq_profiles_warn("no ffp, pp records in the file") + ENDIF if(psio.lt.0)psio=-psio c----------------------------------------------------------------------- c copy and revise 2D arrays. @@ -855,7 +869,7 @@ SUBROUTINE read_eq_ldp_i c----------------------------------------------------------------------- c process inverse equilibrium and deallocate local arrays. c----------------------------------------------------------------------- - DEALLOCATE(psi,f,p,q,r,z) + DEALLOCATE(psi,f,p,q,r,z,ffp,pp) CALL inverse_run c----------------------------------------------------------------------- c terminate. @@ -1425,7 +1439,8 @@ END SUBROUTINE read_eq_rtaylor c----------------------------------------------------------------------- SUBROUTINE read_eq_dump - INTEGER :: nqty + INTEGER :: nqty,ios + LOGICAL, DIMENSION(2) :: slopes LOGICAL, PARAMETER :: diagnose=.TRUE. INTEGER :: ix,iy,iqty,iside,ipsi @@ -1457,14 +1472,15 @@ SUBROUTINE read_eq_dump c write 1D splines. c----------------------------------------------------------------------- READ(dump_unit)mpsi,nqty + IF(sq%allocated)CALL spline_dealloc(sq) CALL spline_alloc(sq,mpsi,nqty) READ(dump_unit)sq%xs,sq%fs, $ sq%xpower,sq%x0,sq%title,sq%name,sq%periodic - CALL spline_fit(sq,"extrap") c----------------------------------------------------------------------- c read 2D splines. c----------------------------------------------------------------------- READ(dump_unit)mpsi,mtheta,nqty + IF(rzphi%allocated)CALL bicube_dealloc(rzphi) CALL bicube_alloc(rzphi,mpsi,mtheta,nqty) READ(dump_unit)rzphi%xs READ(dump_unit)rzphi%ys @@ -1479,6 +1495,23 @@ SUBROUTINE read_eq_dump READ(dump_unit)rzphi%name READ(dump_unit)rzphi%periodic c----------------------------------------------------------------------- +c read supplied f, p derivatives and eqfun if present; fit sq. +c----------------------------------------------------------------------- + READ(dump_unit,IOSTAT=ios)slopes,sq%fs1(:,1:2) + IF(ios == 0)sq_in_slopes(1:2)=slopes + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) + IF(eqfun%allocated)CALL bicube_dealloc(eqfun) + CALL bicube_alloc(eqfun,mpsi,mtheta,3) + READ(dump_unit,IOSTAT=ios)eqfun%fs + IF(ios == 0)THEN + eqfun%xs=sq%xs + eqfun%ys=rzphi%ys + eqfun%title=(/" b0 "," "," " /) + CALL bicube_fit(eqfun,"extrap","periodic") + ELSE + CALL bicube_dealloc(eqfun) + ENDIF +c----------------------------------------------------------------------- c restore x powers and fit. c----------------------------------------------------------------------- ALLOCATE(xfac(0:rzphi%mx)) @@ -1505,10 +1538,11 @@ SUBROUTINE read_eq_dump DO ipsi=0,mpsi ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac + sq%fs1(ipsi,1)=sq%fs1(ipsi,1)/ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac ENDDO - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) ENDIF qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(one-sq%xs(mpsi)) c----------------------------------------------------------------------- @@ -2378,4 +2412,116 @@ SUBROUTINE read_eq_marklin_inverse c----------------------------------------------------------------------- RETURN END SUBROUTINE read_eq_marklin_inverse +c----------------------------------------------------------------------- +c subprogram 32. read_eq_profiles_from_derivatives. +c profile_source = "integrate": replaces the f and mu0 p columns +c of sq_in by integrals of ffp = f*df/dpsifac and +c pp = d(mu0 p)/dpsifac inward from the boundary values, and sets +c ffp/f and pp as knot derivatives for spline_fit_hermite. Leaves +c sq_in unchanged for "values", or, with a warning, if the +c derivatives are unusable. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) + + REAL(r8), DIMENSION(0:), INTENT(IN) :: ffp,pp + + INTEGER :: mx + REAL(r8) :: sf,sp,sgn,dfmax,dpmax + REAL(r8), DIMENSION(0:SIZE(ffp)-1) :: f,p,f2,fint,pint + TYPE(spline_type) :: spl +c----------------------------------------------------------------------- +c select method. +c----------------------------------------------------------------------- + SELECT CASE(TRIM(profile_source)) + CASE("values") + RETURN + CASE("integrate") + CASE DEFAULT + CALL program_stop("Cannot recognize profile_source = " + $ //TRIM(profile_source)) + END SELECT + mx=sq_in%mx + f=sq_in%fs(:,1) + p=sq_in%fs(:,2) +c----------------------------------------------------------------------- +c check derivatives are finite and that f and p agree on their sign +c relative to the values (that common sign is applied below). +c----------------------------------------------------------------------- + IF(.NOT.(ALL(ABS(ffp) <= HUGE(one)) .AND. + $ ALL(ABS(pp) <= HUGE(one)) .AND. ALL(f /= 0)))THEN + CALL read_eq_profiles_warn("non-finite derivatives or f = 0") + RETURN + ENDIF + sf=SUM((f(1:mx)-f(0:mx-1)) + $ *(ffp(1:mx)/f(1:mx)+ffp(0:mx-1)/f(0:mx-1))) + sp=SUM((p(1:mx)-p(0:mx-1))*(pp(1:mx)+pp(0:mx-1))) + sgn=SIGN(one,sf+sp) + IF(sf*sp < 0)THEN + CALL read_eq_profiles_warn("derivative signs disagree with" + $ //" the values for only one of f and p") + RETURN + ENDIF +c----------------------------------------------------------------------- +c integrate derivatives inward from the boundary. +c----------------------------------------------------------------------- + CALL spline_alloc(spl,mx,2) + spl%xs=sq_in%xs + spl%fs(:,1)=sgn*ffp + spl%fs(:,2)=sgn*pp + CALL spline_fit(spl,"extrap") + CALL spline_int(spl) + f2=f(mx)**2-2*(spl%fsi(mx,1)-spl%fsi(:,1)) + pint=p(mx)-(spl%fsi(mx,2)-spl%fsi(:,2)) + CALL spline_dealloc(spl) + IF(ANY(f2 <= 0))THEN + CALL read_eq_profiles_warn("integrated f**2 <= 0") + RETURN + ENDIF + fint=SIGN(SQRT(f2),f(mx)) +c----------------------------------------------------------------------- +c report consistency of tabulated values and derivatives. +c----------------------------------------------------------------------- + dfmax=MAXVAL(ABS(fint-f))/MAXVAL(ABS(f)) + dpmax=MAXVAL(ABS(pint-p))/MAX(MAXVAL(ABS(p)),TINY(one)) + IF(verbose)WRITE(*,'(1x,2a,1p,2(a,es9.2))')"profile_source = ", + $ TRIM(profile_source),", max |df|/|f| =",dfmax, + $ ", max |dp|/|p| =",dpmax + IF(dfmax > 1e-4 .OR. dpmax > 5e-2)WRITE(*,'(1x,a)') + $ "WARNING: tabulated values and derivatives of f, p disagree" +c----------------------------------------------------------------------- +c store values and derivatives. +c----------------------------------------------------------------------- + sq_in%fs(:,1)=fint + sq_in%fs(:,2)=pint + sq_in%fs1(:,1)=sgn*ffp/fint + sq_in%fs1(:,2)=sgn*pp + sq_in_slopes(1:2)=.TRUE. +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE read_eq_profiles_from_derivatives +c----------------------------------------------------------------------- +c subprogram 33. read_eq_profiles_warn. +c warns that profile derivatives are not used. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE read_eq_profiles_warn(reason) + + CHARACTER(*), INTENT(IN) :: reason +c----------------------------------------------------------------------- +c write warning. +c----------------------------------------------------------------------- + WRITE(*,'(1x,3a)')"WARNING: profile_source = values used: ", + $ reason,"." +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE read_eq_profiles_warn END MODULE read_eq_mod diff --git a/equil/spline.f b/equil/spline.f index f7aca7896..4256cfc25 100644 --- a/equil/spline.f +++ b/equil/spline.f @@ -32,6 +32,7 @@ c 21. spline_thomas. c 22. spline_roots. c 23. spline_refine_root. +c 24. spline_fit_hermite. c----------------------------------------------------------------------- c subprogram 0. spline_type definition. c defines spline_type. @@ -1849,5 +1850,39 @@ SUBROUTINE spline_refine_root(spl, iqty, x, op_tol) c----------------------------------------------------------------------- RETURN END SUBROUTINE spline_refine_root +c----------------------------------------------------------------------- +c subprogram 24. spline_fit_hermite. +c fits cubic splines, keeping the caller's knot derivatives fs1 +c for the quantities flagged in keep (cubic Hermite interpolation). +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE spline_fit_hermite(spl,endmode,keep) + + TYPE(spline_type), INTENT(INOUT) :: spl + CHARACTER(*), INTENT(IN) :: endmode + LOGICAL, DIMENSION(spl%nqty), INTENT(IN) :: keep + + INTEGER :: iqty + REAL(r8), DIMENSION(0:spl%mx,spl%nqty) :: fs1 +c----------------------------------------------------------------------- +c fit all quantities, then restore the supplied derivatives. +c----------------------------------------------------------------------- + DO iqty=1,spl%nqty + IF(keep(iqty) .AND. ANY(spl%xpower(:,iqty) /= 0)) + $ CALL program_stop("spline_fit_hermite: "// + $ "supplied derivatives incompatible with xpower") + ENDDO + fs1=spl%fs1 + CALL spline_fit(spl,endmode) + DO iqty=1,spl%nqty + IF(keep(iqty))spl%fs1(:,iqty)=fs1(:,iqty) + ENDDO +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE spline_fit_hermite END MODULE spline_mod diff --git a/input/equil.in b/input/equil.in index 98cbb3a91..d922f0ed6 100644 --- a/input/equil.in +++ b/input/equil.in @@ -15,6 +15,7 @@ newq0=0 ! Grad-Shafranov solution invariant adjustment of the q profile to give the specified value of q at the axis. Default 0 uses input file value. etol = 1e-7 ! Tolerance for direct equilibrium reconstruction of flux surfaces/coordinates use_classic_splines = f ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions + profile_source="values" ! efit, and ldp_i files with FF', p' records. "values" (default): spline the tabulated f, p. "integrate": integrate FF', p' inward from the boundary and keep them as knot derivatives; for g-files with single-precision values; can shift pedestal surfaces input_only=f ! Generate information about the input and then quit with no further calculation use_galgrid=t ! Use the same grid as galerkin method when regrid=.True. diff --git a/rdcon/dcon.f b/rdcon/dcon.f index e1874bc50..9874298b0 100644 --- a/rdcon/dcon.f +++ b/rdcon/dcon.f @@ -77,6 +77,11 @@ SUBROUTINE dcon_dealloc(opt) DEALLOCATE(cellinfos%icell,cellinfos%iintvl, $ cellinfos%etypes,cellinfos%etypes_int,cellinfos%x1,cellinfos%x2) ENDIF + IF(MRE_flag)THEN + CALL spline_dealloc(mreterms) + DEALLOCATE(Zeff_) + DEALLOCATE(psi_N_Zeff_) + ENDIF c----------------------------------------------------------------------- c terminate. c----------------------------------------------------------------------- @@ -226,7 +231,8 @@ SUBROUTINE dcon_run $ tol_nr,tol_r,crossover,ucrit,singfac_min,singfac_max, $ cyl_flag,dmlim,lim_flag,sas_flag,sing_order,sort_type, $ gal_flag,regrid_flag,sing1_flag,qlow,qhigh, - $ sing_order_ceiling,degen_tol,coil,Zeff,reform_eq_with_psilim + $ sing_order_ceiling,degen_tol,coil,reform_eq_with_psilim, + $ Zeff, psi_N_Zeff NAMELIST/rdcon_output/interp,crit_break,out_bal1, $ bin_bal1,out_bal2,bin_bal2,out_metric,bin_metric,out_fmat, $ bin_fmat,out_gmat,bin_gmat,out_kmat,bin_kmat,out_sol, @@ -260,7 +266,6 @@ SUBROUTINE dcon_run CALL gal_read_input ENDIF CALL ascii_close(in_unit) - c----------------------------------------------------------------------- c open output files, read, process, and diagnose equilibrium. c----------------------------------------------------------------------- @@ -339,6 +344,8 @@ SUBROUTINE dcon_run c optionally compute modified Rutherford equation (MRE) terms. c----------------------------------------------------------------------- IF(MRE_flag)THEN + CALL read_var_len(Zeff,Zeff_) + CALL read_var_len(psi_N_Zeff,psi_N_Zeff_) CALL spline_alloc(mreterms,mpsi,30) mreterms%xs=sq%xs mreterms%fs=0 diff --git a/rdcon/dcon_mod.f b/rdcon/dcon_mod.f index 99bf7444c..74c161d1c 100644 --- a/rdcon/dcon_mod.f +++ b/rdcon/dcon_mod.f @@ -80,7 +80,10 @@ MODULE rdcon_mod $ sing_start=0,nzero=0 REAL(r8) :: thmax0=1,ucrit=1e4,tol_r=1e-5,tol_nr=1e-5, $ crossover=1e-2,mthsurf0=1,prefac=1.0, - $ plasma1=0.0,vacuum1=0.0,total1=0.0,Zeff=1.52 + $ plasma1=0.0,vacuum1=0.0,total1=0.0 + + REAL(r8), DIMENSION(:), ALLOCATABLE :: Zeff_, psi_N_Zeff_ + REAL(r8), DIMENSION(1:1000) :: Zeff=-1.0, psi_N_Zeff=-1.0 REAL(r8) :: dpsi_intvl=0.1,dpsi1_intvl=0.1 TYPE(spline_type) :: locstab, mreterms diff --git a/rdcon/mercier.f b/rdcon/mercier.f index f86faf262..b3e03ed23 100644 --- a/rdcon/mercier.f +++ b/rdcon/mercier.f @@ -37,12 +37,20 @@ SUBROUTINE mercier_scan $ eps_eff REAL(r8), DIMENSION(:), POINTER :: avg REAL(r8), DIMENSION(2) :: rmax_loc, rmin_loc, zmax_loc, zmin_loc - TYPE(spline_type), TARGET :: ff + TYPE(spline_type), TARGET :: ff, Zeff_spline c----------------------------------------------------------------------- c prepare spline types. c----------------------------------------------------------------------- IF(MRE_flag)THEN CALL spline_alloc(ff,mtheta,22) + IF(ALLOCATED(Zeff_) .AND. ALLOCATED(psi_N_Zeff_))THEN + CALL spline_alloc(Zeff_spline,SIZE(psi_N_Zeff_)-1,1) + Zeff_spline%xs=psi_N_Zeff_ + Zeff_spline%fs(:,1)=Zeff_ + CALL spline_fit(Zeff_spline,"extrap") + ELSE + CALL program_stop("Provide Zeff profile values") + ENDIF ELSE CALL spline_alloc(ff,mtheta,5) ENDIF @@ -207,8 +215,9 @@ SUBROUTINE mercier_scan c----------------------------------------------------------------------- ! mu_e_on_nu_e=(ftr/(1.d0-ftr))*(1.d0+0.533d0/Zeff) ! ^ Taking banana limit of eq. B17 (& B14) in Callen - mufrac=ftr*(1.d0+0.533d0/Zeff)/ - $ ((1.d0-ftr)+ftr*(1.d0+0.533d0/Zeff)) + CALL spline_eval(Zeff_spline,psifac,0) + mufrac=ftr*(1.d0+0.533d0/Zeff_spline%f(1))/ + $ ((1.d0-ftr)+ftr*(1.d0+0.533d0/Zeff_spline%f(1))) c----------------------------------------------------------------------- c evaluate geometric prefactors of MRE stability terms from c Hegna 1999 https://doi.org/10.1063/1.873661 @@ -279,6 +288,7 @@ SUBROUTINE mercier_scan 120 FORMAT(19(E30.15,1X)) ENDDO CALL spline_dealloc(ff) + IF(MRE_flag) CALL spline_dealloc(Zeff_spline) c----------------------------------------------------------------------- c terminate. c----------------------------------------------------------------------- diff --git a/rdcon/rdcon_netcdf.f b/rdcon/rdcon_netcdf.f index e605886ed..c0e06add9 100644 --- a/rdcon/rdcon_netcdf.f +++ b/rdcon/rdcon_netcdf.f @@ -131,7 +131,11 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et) CALL check( nf90_put_att(ncid,nf90_global,'betat', betat)) CALL check( nf90_put_att(ncid,nf90_global,'betan', betan)) CALL check( nf90_put_att(ncid,nf90_global,'volume', volume)) - CALL check( nf90_put_att(ncid,nf90_global,'Zeff', Zeff)) + IF(ALLOCATED(Zeff_) .AND. ALLOCATED(psi_N_Zeff_))THEN + CALL check( nf90_put_att(ncid,nf90_global,'Zeff', Zeff_)) + CALL check( nf90_put_att(ncid,nf90_global,'psi_N_Zeff', + $ psi_N_Zeff_)) + ENDIF CALL check( nf90_put_att(ncid,nf90_global,'bwall', bwall)) CALL check( nf90_put_att(ncid,nf90_global,'bt0', bt0)) CALL check( nf90_put_att(ncid,nf90_global,'q0', q0))