From a68b6409a8daf2a0d3e56315ed09aae07d19cd5a Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 10:40:49 -0400 Subject: [PATCH 01/13] switch from double precision to dp type in xnet_types --- tools/starkiller-helmholtz/actual_eos.F90 | 569 +++++++++++----------- tools/starkiller-helmholtz/eos_type.F90 | 122 +++-- 2 files changed, 345 insertions(+), 346 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 3a9159ae..2bb806b2 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -1,5 +1,6 @@ module actual_eos_module + use xnet_types, only: dp use eos_type_module character (len=64), public :: eos_name = "helmholtz" @@ -8,92 +9,92 @@ module actual_eos_module logical, allocatable :: do_coulomb logical, allocatable :: input_is_constant - double precision, parameter, private :: ZERO = 0.0d0, HALF = 0.5d0, TWO = 2.0d0 + real(dp), parameter, private :: ZERO = 0.0_dp, HALF = 0.5_dp, TWO = 2.0_dp !..for the tables, in general integer, parameter, private :: imax = 541, jmax = 201 integer, allocatable :: itmax, jtmax - double precision, allocatable :: d(:), t(:) + real(dp), allocatable :: d(:), t(:) - double precision, allocatable :: tlo, thi, tstp, tstpi - double precision, allocatable :: dlo, dhi, dstp, dstpi + real(dp), allocatable :: tlo, thi, tstp, tstpi + real(dp), allocatable :: dlo, dhi, dstp, dstpi - double precision, allocatable :: ttol, dtol + real(dp), allocatable :: ttol, dtol !..for the helmholtz free energy tables - double precision, allocatable :: f(:,:), fd(:,:), & - ft(:,:), fdd(:,:), ftt(:,:), & - fdt(:,:), fddt(:,:), fdtt(:,:), & - fddtt(:,:) + real(dp), allocatable :: f(:,:), fd(:,:), & + ft(:,:), fdd(:,:), ftt(:,:), & + fdt(:,:), fddt(:,:), fdtt(:,:), & + fddtt(:,:) !..for the pressure derivative with density ables - double precision, allocatable :: dpdf(:,:), dpdfd(:,:), & - dpdft(:,:), dpdfdt(:,:) + real(dp), allocatable :: dpdf(:,:), dpdfd(:,:), & + dpdft(:,:), dpdfdt(:,:) !..for chemical potential tables - double precision, allocatable :: ef(:,:), efd(:,:), & - eft(:,:), efdt(:,:) + real(dp), allocatable :: ef(:,:), efd(:,:), & + eft(:,:), efdt(:,:) !..for the number density tables - double precision, allocatable :: xf(:,:), xfd(:,:), & - xft(:,:), xfdt(:,:) + real(dp), allocatable :: xf(:,:), xfd(:,:), & + xft(:,:), xfdt(:,:) !..for storing the differences - double precision, allocatable :: dt_sav(:), dt2_sav(:), & - dti_sav(:), dt2i_sav(:), & - dd_sav(:), dd2_sav(:), & - ddi_sav(:), dd2i_sav(:) + real(dp), allocatable :: dt_sav(:), dt2_sav(:), & + dti_sav(:), dt2i_sav(:), & + dd_sav(:), dd2_sav(:), & + ddi_sav(:), dd2i_sav(:) integer, parameter :: max_newton = 100 ! 2006 CODATA physical constants ! Math constants - double precision, parameter :: pi = 3.1415926535897932384d0 + real(dp), parameter :: pi = 3.1415926535897932384_dp ! Physical constants - double precision, parameter :: h = 6.6260689633d-27 - double precision, parameter :: hbar = 0.5d0 * h/pi - double precision, parameter :: qe = 4.8032042712d-10 - double precision, parameter :: avo_eos = 6.0221417930d23 - double precision, parameter :: clight = 2.99792458d10 - double precision, parameter :: kerg = 1.380650424d-16 - double precision, parameter :: ev2erg_eos = 1.60217648740d-12 - double precision, parameter :: kev = kerg/ev2erg_eos - double precision, parameter :: amu = 1.66053878283d-24 - double precision, parameter :: me_eos = 9.1093821545d-28 - double precision, parameter :: rbohr = hbar*hbar/(me_eos * qe * qe) - double precision, parameter :: fine = qe*qe/(hbar*clight) - - double precision, parameter :: ssol = 5.67051d-5 - double precision, parameter :: asol = 4.0d0 * ssol / clight - double precision, parameter :: weinlam = h*clight/(kerg * 4.965114232d0) - double precision, parameter :: weinfre = 2.821439372d0*kerg/h + real(dp), parameter :: h = 6.6260689633e-27_dp + real(dp), parameter :: hbar = 0.5_dp * h/pi + real(dp), parameter :: qe = 4.8032042712e-10_dp + real(dp), parameter :: avo_eos = 6.0221417930e23_dp + real(dp), parameter :: clight = 2.99792458e10_dp + real(dp), parameter :: kerg = 1.380650424e-16_dp + real(dp), parameter :: ev2erg_eos = 1.60217648740e-12_dp + real(dp), parameter :: kev = kerg/ev2erg_eos + real(dp), parameter :: amu = 1.66053878283e-24_dp + real(dp), parameter :: me_eos = 9.1093821545e-28_dp + real(dp), parameter :: rbohr = hbar*hbar/(me_eos * qe * qe) + real(dp), parameter :: fine = qe*qe/(hbar*clight) + + real(dp), parameter :: ssol = 5.67051e-5_dp + real(dp), parameter :: asol = 4.0_dp * ssol / clight + real(dp), parameter :: weinlam = h*clight/(kerg * 4.965114232_dp) + real(dp), parameter :: weinfre = 2.821439372_dp*kerg/h ! Astronomical constants - double precision, parameter :: ly = 9.460528d17 - double precision, parameter :: pc = 3.261633d0 * ly + real(dp), parameter :: ly = 9.460528e17_dp + real(dp), parameter :: pc = 3.261633_dp * ly ! Some other useful combinations of the constants - double precision, parameter :: sioncon = (2.0d0 * pi * amu * kerg)/(h*h) - double precision, parameter :: forth = 4.0d0/3.0d0 - double precision, parameter :: forpi = 4.0d0 * pi - double precision, parameter :: kergavo = kerg * avo_eos - double precision, parameter :: ikavo = 1.0d0/kergavo - double precision, parameter :: asoli3 = asol/3.0d0 - double precision, parameter :: light2 = clight * clight + real(dp), parameter :: sioncon = (2.0_dp * pi * amu * kerg)/(h*h) + real(dp), parameter :: forth = 4.0_dp/3.0_dp + real(dp), parameter :: forpi = 4.0_dp * pi + real(dp), parameter :: kergavo = kerg * avo_eos + real(dp), parameter :: ikavo = 1.0_dp/kergavo + real(dp), parameter :: asoli3 = asol/3.0_dp + real(dp), parameter :: light2 = clight * clight ! Constants used for the Coulomb corrections - double precision, parameter :: a1 = -0.898004d0 - double precision, parameter :: b1 = 0.96786d0 - double precision, parameter :: c1 = 0.220703d0 - double precision, parameter :: d1 = -0.86097d0 - double precision, parameter :: e1 = 2.5269d0 - double precision, parameter :: a2 = 0.29561d0 - double precision, parameter :: b2 = 1.9885d0 - double precision, parameter :: c2 = 0.288675d0 - double precision, parameter :: onethird = 1.0d0/3.0d0 - double precision, parameter :: esqu = qe * qe + real(dp), parameter :: a1 = -0.898004_dp + real(dp), parameter :: b1 = 0.96786_dp + real(dp), parameter :: c1 = 0.220703_dp + real(dp), parameter :: d1 = -0.86097_dp + real(dp), parameter :: e1 = 2.5269_dp + real(dp), parameter :: a2 = 0.29561_dp + real(dp), parameter :: b2 = 1.9885_dp + real(dp), parameter :: c2 = 0.288675_dp + real(dp), parameter :: onethird = 1.0_dp/3.0_dp + real(dp), parameter :: esqu = qe * qe !$acc declare & !$acc create(tlo, thi, dlo, dhi) & @@ -174,82 +175,82 @@ subroutine actual_eos(input, state) type (eos_t), intent(inout) :: state !..rows to store EOS data - double precision :: temp_row, & - den_row, & - abar_row, & - zbar_row, & - ye_row, & - etot_row, & - ptot_row, & - cv_row, & - cp_row, & - xne_row, & - xnp_row, & - etaele_row, & - detadt_row, & - pele_row, & - ppos_row, & - dpd_row, & - dpt_row, & - dpa_row, & - dpz_row, & - ded_row, & - det_row, & - dea_row, & - dez_row, & - stot_row, & - dsd_row, & - dst_row, & - htot_row, & - dhd_row, & - dht_row, & - dpe_row, & - dpdr_e_row, & - gam1_row, & - cs_row + real(dp) :: temp_row, & + den_row, & + abar_row, & + zbar_row, & + ye_row, & + etot_row, & + ptot_row, & + cv_row, & + cp_row, & + xne_row, & + xnp_row, & + etaele_row, & + detadt_row, & + pele_row, & + ppos_row, & + dpd_row, & + dpt_row, & + dpa_row, & + dpz_row, & + ded_row, & + det_row, & + dea_row, & + dez_row, & + stot_row, & + dsd_row, & + dst_row, & + htot_row, & + dhd_row, & + dht_row, & + dpe_row, & + dpdr_e_row, & + gam1_row, & + cs_row !..declare local variables logical :: single_iter, double_iter, converged integer :: var, dvar, var1, var2, iter - double precision :: v_want - double precision :: v1_want, v2_want - double precision :: xnew, xtol, dvdx, smallx, error, v - double precision :: v1, v2, dv1dt, dv1dr, dv2dt,dv2dr, delr, error1, error2, told, rold, tnew, rnew, v1i, v2i - - double precision :: x,y,zz,zzi,deni,tempi,xni,dxnidd,dxnida, & - dpepdt,dpepdd,deepdt,deepdd,dsepdd,dsepdt, & - dpraddd,dpraddt,deraddd,deraddt,dpiondd,dpiondt, & - deiondd,deiondt,dsraddd,dsraddt,dsiondd,dsiondt, & - dse,dpe,dsp,kt,ktinv,prad,erad,srad,pion,eion, & - sion,xnem,pele,eele,sele,pres,ener,entr,dpresdd, & - dpresdt,denerdd,denerdt,dentrdd,dentrdt,cv,cp, & - gam1,gam2,gam3,chit,chid,nabad,sound,etaele, & - detadt,detadd,xnefer,dxnedt,dxnedd,s, & - temp,den,abar,zbar,ytot1,ye + real(dp) :: v_want + real(dp) :: v1_want, v2_want + real(dp) :: xnew, xtol, dvdx, smallx, error, v + real(dp) :: v1, v2, dv1dt, dv1dr, dv2dt,dv2dr, delr, error1, error2, told, rold, tnew, rnew, v1i, v2i + + real(dp) :: x,y,zz,zzi,deni,tempi,xni,dxnidd,dxnida, & + dpepdt,dpepdd,deepdt,deepdd,dsepdd,dsepdt, & + dpraddd,dpraddt,deraddd,deraddt,dpiondd,dpiondt, & + deiondd,deiondt,dsraddd,dsraddt,dsiondd,dsiondt, & + dse,dpe,dsp,kt,ktinv,prad,erad,srad,pion,eion, & + sion,xnem,pele,eele,sele,pres,ener,entr,dpresdd, & + dpresdt,denerdd,denerdt,dentrdd,dentrdt,cv,cp, & + gam1,gam2,gam3,chit,chid,nabad,sound,etaele, & + detadt,detadd,xnefer,dxnedt,dxnedd,s, & + temp,den,abar,zbar,ytot1,ye !..for the interpolations - integer :: iat,jat - double precision :: free,df_d,df_t,df_tt,df_dt - double precision :: xt,xd,mxt,mxd, & - si0t,si1t,si2t,si0mt,si1mt,si2mt, & - si0d,si1d,si2d,si0md,si1md,si2md, & - dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt, & - dsi0d,dsi1d,dsi2d,dsi0md,dsi1md,dsi2md, & - ddsi0t,ddsi1t,ddsi2t,ddsi0mt,ddsi1mt,ddsi2mt, & - z,din,fi(36) + integer :: iat,jat + real(dp) :: free,df_d,df_t,df_tt,df_dt + real(dp) :: xt,xd,mxt,mxd, & + si0t,si1t,si2t,si0mt,si1mt,si2mt, & + si0d,si1d,si2d,si0md,si1md,si2md, & + dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt, & + dsi0d,dsi1d,dsi2d,dsi0md,dsi1md,dsi2md, & + ddsi0t,ddsi1t,ddsi2t,ddsi0mt,ddsi1mt,ddsi2mt, & + z,din,fi(36) !..for the coulomb corrections - double precision :: dsdd,dsda,lami,inv_lami,lamida,lamidd, & - plasg,plasgdd,plasgdt,plasgda,plasgdz, & - ecoul,decouldd,decouldt,decoulda,decouldz, & - pcoul,dpcouldd,dpcouldt,dpcoulda,dpcouldz, & - scoul,dscouldd,dscouldt,dscoulda,dscouldz + real(dp) :: dsdd,dsda,lami,inv_lami,lamida,lamidd, & + plasg,plasgdd,plasgdt,plasgda,plasgdz, & + ecoul,decouldd,decouldt,decoulda,decouldz, & + pcoul,dpcouldd,dpcouldt,dpcoulda,dpcouldz, & + scoul,dscouldd,dscouldt,dscoulda,dscouldz - double precision :: p_temp, e_temp + real(dp) :: p_temp, e_temp - double precision :: smallt, smalld + real(dp) :: smallt, smalld !$gpu @@ -324,41 +325,41 @@ subroutine actual_eos(input, state) endif - ptot_row = 0.0d0 - dpt_row = 0.0d0 - dpd_row = 0.0d0 - dpa_row = 0.0d0 - dpz_row = 0.0d0 - dpe_row = 0.0d0 - dpdr_e_row = 0.0d0 + ptot_row = 0.0_dp + dpt_row = 0.0_dp + dpd_row = 0.0_dp + dpa_row = 0.0_dp + dpz_row = 0.0_dp + dpe_row = 0.0_dp + dpdr_e_row = 0.0_dp - etot_row = 0.0d0 - det_row = 0.0d0 - ded_row = 0.0d0 - dea_row = 0.0d0 - dez_row = 0.0d0 + etot_row = 0.0_dp + det_row = 0.0_dp + ded_row = 0.0_dp + dea_row = 0.0_dp + dez_row = 0.0_dp - stot_row = 0.0d0 - dst_row = 0.0d0 - dsd_row = 0.0d0 + stot_row = 0.0_dp + dst_row = 0.0_dp + dsd_row = 0.0_dp - htot_row = 0.0d0 - dhd_row = 0.0d0 - dht_row = 0.0d0 + htot_row = 0.0_dp + dhd_row = 0.0_dp + dht_row = 0.0_dp - pele_row = 0.0d0 - ppos_row = 0.0d0 + pele_row = 0.0_dp + ppos_row = 0.0_dp - xne_row = 0.0d0 - xnp_row = 0.0d0 + xne_row = 0.0_dp + xnp_row = 0.0_dp - etaele_row = 0.0d0 - detadt_row = 0.0d0 + etaele_row = 0.0_dp + detadt_row = 0.0_dp - cv_row = 0.0d0 - cp_row = 0.0d0 - cs_row = 0.0d0 - gam1_row = 0.0d0 + cv_row = 0.0_dp + cp_row = 0.0_dp + cs_row = 0.0_dp + gam1_row = 0.0_dp converged = .false. @@ -371,24 +372,24 @@ subroutine actual_eos(input, state) abar = abar_row zbar = zbar_row - ytot1 = 1.0d0 / abar + ytot1 = 1.0_dp / abar ye = ye_row din = ye * den !..initialize - deni = 1.0d0/den - tempi = 1.0d0/temp + deni = 1.0_dp/den + tempi = 1.0_dp/temp kt = kerg * temp - ktinv = 1.0d0/kt + ktinv = 1.0_dp/kt !..radiation section: prad = asoli3 * temp * temp * temp * temp - dpraddd = 0.0d0 - dpraddt = 4.0d0 * prad*tempi + dpraddd = 0.0_dp + dpraddt = 4.0_dp * prad*tempi - erad = 3.0d0 * prad*deni + erad = 3.0_dp * prad*deni deraddd = -erad*deni - deraddt = 3.0d0 * dpraddt*deni + deraddt = 3.0_dp * dpraddt*deni srad = (prad*deni + erad)*tempi dsraddd = (dpraddd*deni - prad*deni*deni + deraddd)*tempi @@ -403,9 +404,9 @@ subroutine actual_eos(input, state) dpiondd = dxnidd * kt dpiondt = xni * kerg - eion = 1.5d0 * pion*deni - deiondd = (1.5d0 * dpiondd - eion)*deni - deiondt = 1.5d0 * dpiondt*deni + eion = 1.5_dp * pion*deni + deiondd = (1.5_dp * dpiondd - eion)*deni + deiondt = 1.5_dp * dpiondt*deni x = abar*abar*sqrt(abar) * deni/avo_eos s = sioncon * temp @@ -416,7 +417,7 @@ subroutine actual_eos(input, state) - kergavo * deni * ytot1 dsiondt = (dpiondt*deni + deiondt)*tempi - & (pion*deni + eion) * tempi*tempi & - + 1.5d0 * kergavo * tempi*ytot1 + + 1.5_dp * kergavo * tempi*ytot1 x = avo_eos*kerg/abar !..electron-positron section: @@ -471,10 +472,10 @@ subroutine actual_eos(input, state) fi(36) = fddtt(iat+1,jat+1) !..various differences - xt = max( (temp - t(jat))*dti_sav(jat), 0.0d0) - xd = max( (din - d(iat))*ddi_sav(iat), 0.0d0) - mxt = 1.0d0 - xt - mxd = 1.0d0 - xd + xt = max( (temp - t(jat))*dti_sav(jat), 0.0_dp) + xd = max( (din - d(iat))*ddi_sav(iat), 0.0_dp) + mxt = 1.0_dp - xt + mxd = 1.0_dp - xd !..the six density and six temperature basis functions si0t = psi0(xt) @@ -608,7 +609,7 @@ subroutine actual_eos(input, state) dpepdd = h3( fi, & si0t, si1t, si0mt, si1mt, & si0d, si1d, si0md, si1md) - dpepdd = max(ye * dpepdd,0.0d0) + dpepdd = max(ye * dpepdd,0.0_dp) !..look in the electron chemical potential table only once fi(1) = ef(iat,jat) @@ -672,7 +673,7 @@ subroutine actual_eos(input, state) x = h3( fi, & si0t, si1t, si0mt, si1mt, & dsi0d, dsi1d, dsi0md, dsi1md) - x = max(x,0.0d0) + x = max(x,0.0_dp) dxnedd = ye * x !..derivative with respect to temperature @@ -690,8 +691,8 @@ subroutine actual_eos(input, state) x = din * din pele = x * df_d dpepdt = x * df_dt - ! dpepdd = ye * (x * df_dd + 2.0d0 * din * df_d) - s = dpepdd/ye - 2.0d0 * din * df_d + ! dpepdd = ye * (x * df_dd + 2.0_dp * din * df_d) + s = dpepdd/ye - 2.0_dp * din * df_d x = ye * ye sele = -df_t * ye @@ -704,21 +705,21 @@ subroutine actual_eos(input, state) !..coulomb section: !..initialize - pcoul = 0.0d0 - dpcouldd = 0.0d0 - dpcouldt = 0.0d0 - dpcoulda = 0.0d0 - dpcouldz = 0.0d0 - ecoul = 0.0d0 - decouldd = 0.0d0 - decouldt = 0.0d0 - decoulda = 0.0d0 - decouldz = 0.0d0 - scoul = 0.0d0 - dscouldd = 0.0d0 - dscouldt = 0.0d0 - dscoulda = 0.0d0 - dscouldz = 0.0d0 + pcoul = 0.0_dp + dpcouldd = 0.0_dp + dpcouldt = 0.0_dp + dpcoulda = 0.0_dp + dpcouldz = 0.0_dp + ecoul = 0.0_dp + decouldd = 0.0_dp + decouldt = 0.0_dp + decoulda = 0.0_dp + decouldz = 0.0_dp + scoul = 0.0_dp + dscouldd = 0.0_dp + dscouldt = 0.0_dp + dscoulda = 0.0_dp + dscouldz = 0.0_dp !..uniform background corrections only !..from yakovlev & shalybkov 1989 @@ -729,8 +730,8 @@ subroutine actual_eos(input, state) dsdd = z * dxnidd dsda = z * dxnida - lami = 1.0d0/s**onethird - inv_lami = 1.0d0/lami + lami = 1.0_dp/s**onethird + inv_lami = 1.0_dp/lami z = -onethird * lami lamidd = z * dsdd/s lamida = z * dsda/s @@ -740,20 +741,20 @@ subroutine actual_eos(input, state) plasgdd = z * lamidd plasgda = z * lamida plasgdt = -plasg*ktinv * kerg - plasgdz = 2.0d0 * plasg/zbar + plasgdz = 2.0_dp * plasg/zbar ! TURN ON/OFF COULOMB if ( do_coulomb ) then !...yakovlev & shalybkov 1989 equations 82, 85, 86, 87 - if (plasg .ge. 1.0D0) then - x = plasg**(0.25d0) + if (plasg .ge. 1.0_dp) then + x = plasg**(0.25_dp) y = avo_eos * ytot1 * kerg ecoul = y * temp * (a1*plasg + b1*x + c1/x + d1) pcoul = onethird * den * ecoul - scoul = -y * (3.0d0*b1*x - 5.0d0*c1/x & - + d1 * (log(plasg) - 1.0d0) - e1) + scoul = -y * (3.0_dp*b1*x - 5.0_dp*c1/x & + + d1 * (log(plasg) - 1.0_dp) - e1) - y = avo_eos*ytot1*kt*(a1 + 0.25d0/plasg*(b1*x - c1/x)) + y = avo_eos*ytot1*kt*(a1 + 0.25_dp/plasg*(b1*x - c1/x)) decouldd = y * plasgdd decouldt = y * plasgdt + ecoul/temp decoulda = y * plasgda - ecoul/abar @@ -766,33 +767,33 @@ subroutine actual_eos(input, state) dpcouldz = y * decouldz y = -avo_eos*kerg/(abar*plasg)* & - (0.75d0*b1*x+1.25d0*c1/x+d1) + (0.75_dp*b1*x+1.25_dp*c1/x+d1) dscouldd = y * plasgdd dscouldt = y * plasgdt dscoulda = y * plasgda - scoul/abar dscouldz = y * plasgdz !...yakovlev & shalybkov 1989 equations 102, 103, 104 - else if (plasg .lt. 1.0D0) then + else if (plasg .lt. 1.0_dp) then x = plasg*sqrt(plasg) y = plasg**b2 z = c2 * x - onethird * a2 * y pcoul = -pion * z - ecoul = 3.0d0 * pcoul/den - scoul = -avo_eos/abar*kerg*(c2*x -a2*(b2-1.0d0)/b2*y) + ecoul = 3.0_dp * pcoul/den + scoul = -avo_eos/abar*kerg*(c2*x -a2*(b2-1.0_dp)/b2*y) - s = 1.5d0*c2*x/plasg - onethird*a2*b2*y/plasg + s = 1.5_dp*c2*x/plasg - onethird*a2*b2*y/plasg dpcouldd = -dpiondd*z - pion*s*plasgdd dpcouldt = -dpiondt*z - pion*s*plasgdt - s = 3.0d0/den + s = 3.0_dp/den decouldd = s * dpcouldd - ecoul/den decouldt = s * dpcouldt decoulda = s * dpcoulda decouldz = s * dpcouldz s = -avo_eos*kerg/(abar*plasg)* & - (1.5d0*c2*x-a2*(b2-1.0d0)*y) + (1.5_dp*c2*x-a2*(b2-1.0_dp)*y) dscouldd = s * plasgdd dscouldt = s * plasgdt dscoulda = s * plasgda - scoul/abar @@ -807,21 +808,21 @@ subroutine actual_eos(input, state) if (p_temp .le. ZERO .or. e_temp .le. ZERO) then - pcoul = 0.0d0 - dpcouldd = 0.0d0 - dpcouldt = 0.0d0 - dpcoulda = 0.0d0 - dpcouldz = 0.0d0 - ecoul = 0.0d0 - decouldd = 0.0d0 - decouldt = 0.0d0 - decoulda = 0.0d0 - decouldz = 0.0d0 - scoul = 0.0d0 - dscouldd = 0.0d0 - dscouldt = 0.0d0 - dscoulda = 0.0d0 - dscouldz = 0.0d0 + pcoul = 0.0_dp + dpcouldd = 0.0_dp + dpcouldt = 0.0_dp + dpcoulda = 0.0_dp + dpcouldz = 0.0_dp + ecoul = 0.0_dp + decouldd = 0.0_dp + decouldt = 0.0_dp + decoulda = 0.0_dp + decouldz = 0.0_dp + scoul = 0.0_dp + dscouldd = 0.0_dp + dscouldt = 0.0_dp + dscoulda = 0.0_dp + dscouldz = 0.0_dp end if end if @@ -852,19 +853,19 @@ subroutine actual_eos(input, state) chid = dpresdd*zzi cv = denerdt x = zz * chit/(temp * cv) - gam3 = x + 1.0d0 + gam3 = x + 1.0_dp gam1 = chit*x + chid nabad = x/gam1 - gam2 = 1.0d0/(1.0d0 - nabad) + gam2 = 1.0_dp/(1.0_dp - nabad) cp = cv * gam1/chid - z = 1.0d0 + (ener + light2)*zzi + z = 1.0_dp + (ener + light2)*zzi sound = clight * sqrt(gam1/z) !..maxwell relations; each is zero if the consistency is perfect x = den * den - dse = temp*dentrdt/denerdt - 1.0d0 - dpe = (denerdd*x + temp*dpresdt)/pres - 1.0d0 - dsp = -dentrdd*x/dpresdt - 1.0d0 + dse = temp*dentrdt/denerdt - 1.0_dp + dpe = (denerdd*x + temp*dpresdt)/pres - 1.0_dp + dsp = -dentrdd*x/dpresdt - 1.0_dp ptot_row = pres dpt_row = dpresdt @@ -885,10 +886,10 @@ subroutine actual_eos(input, state) dht_row = denerdt + dpresdt / den pele_row = pele - ppos_row = 0.0d0 + ppos_row = 0.0_dp xne_row = xnefer - xnp_row = 0.0d0 + xnp_row = 0.0_dp etaele_row = etaele detadt_row = detadt @@ -1156,9 +1157,9 @@ subroutine actual_eos_init implicit none - double precision :: dth, dt2, dti, dt2i - double precision :: dd, dd2, ddi, dd2i - double precision :: tsav, dsav + real(dp) :: dth, dt2, dti, dt2i + real(dp) :: dd, dd2, ddi, dd2i + real(dp) :: tsav, dsav integer :: i, j integer :: status @@ -1214,27 +1215,27 @@ subroutine actual_eos_init input_is_constant = .true. do_coulomb = .true. - ttol = 1.0d-8 - dtol = 1.0d-8 + ttol = 1.0e-8_dp + dtol = 1.0e-8_dp !.. read the helmholtz free energy table itmax = imax jtmax = jmax - tlo = 3.0d0 - thi = 13.0d0 + tlo = 3.0_dp + thi = 13.0_dp tstp = (thi - tlo)/float(jmax-1) - tstpi = 1.0d0/tstp - dlo = -12.0d0 - dhi = 15.0d0 + tstpi = 1.0_dp/tstp + dlo = -12.0_dp + dhi = 15.0_dp dstp = (dhi - dlo)/float(imax-1) - dstpi = 1.0d0/dstp + dstpi = 1.0_dp/dstp do j=1,jmax tsav = tlo + (j-1)*tstp - t(j) = 10.0d0**(tsav) + t(j) = 10.0_dp**(tsav) do i=1,imax dsav = dlo + (i-1)*dstp - d(i) = 10.0d0**(dsav) + d(i) = 10.0_dp**(dsav) end do end do @@ -1304,8 +1305,8 @@ subroutine actual_eos_init do j = 1, jmax-1 dth = t(j+1) - t(j) dt2 = dth * dth - dti = 1.0d0/dth - dt2i = 1.0d0/dt2 + dti = 1.0_dp/dth + dt2i = 1.0_dp/dt2 dt_sav(j) = dth dt2_sav(j) = dt2 dti_sav(j) = dti @@ -1314,8 +1315,8 @@ subroutine actual_eos_init do i = 1, imax-1 dd = d(i+1) - d(i) dd2 = dd * dd - ddi = 1.0d0/dd - dd2i = 1.0d0/dd2 + ddi = 1.0_dp/dd + dd2i = 1.0_dp/dd2 dd_sav(i) = dd dd2_sav(i) = dd2 ddi_sav(i) = ddi @@ -1354,85 +1355,85 @@ end subroutine actual_eos_init ! psi0 and its derivatives pure function psi0(z) result(psi0r) !$acc routine seq - double precision, intent(in) :: z - double precision :: psi0r + real(dp), intent(in) :: z + real(dp) :: psi0r !$gpu - psi0r = z**3 * ( z * (-6.0d0*z + 15.0d0) -10.0d0) + 1.0d0 + psi0r = z**3 * ( z * (-6.0_dp*z + 15.0_dp) -10.0_dp) + 1.0_dp end function psi0 pure function dpsi0(z) result(dpsi0r) !$acc routine seq - double precision, intent(in) :: z - double precision :: dpsi0r + real(dp), intent(in) :: z + real(dp) :: dpsi0r !$gpu - dpsi0r = z**2 * ( z * (-30.0d0*z + 60.0d0) - 30.0d0) + dpsi0r = z**2 * ( z * (-30.0_dp*z + 60.0_dp) - 30.0_dp) end function dpsi0 pure function ddpsi0(z) result(ddpsi0r) !$acc routine seq - double precision, intent(in) :: z - double precision :: ddpsi0r + real(dp), intent(in) :: z + real(dp) :: ddpsi0r !$gpu - ddpsi0r = z* ( z*( -120.0d0*z + 180.0d0) -60.0d0) + ddpsi0r = z* ( z*( -120.0_dp*z + 180.0_dp) -60.0_dp) end function ddpsi0 ! psi1 and its derivatives pure function psi1(z) result(psi1r) !$acc routine seq - double precision, intent(in) :: z - double precision :: psi1r + real(dp), intent(in) :: z + real(dp) :: psi1r !$gpu - psi1r = z* ( z**2 * ( z * (-3.0d0*z + 8.0d0) - 6.0d0) + 1.0d0) + psi1r = z* ( z**2 * ( z * (-3.0_dp*z + 8.0_dp) - 6.0_dp) + 1.0_dp) end function psi1 pure function dpsi1(z) result(dpsi1r) !$acc routine seq - double precision, intent(in) :: z - double precision :: dpsi1r + real(dp), intent(in) :: z + real(dp) :: dpsi1r !$gpu - dpsi1r = z*z * ( z * (-15.0d0*z + 32.0d0) - 18.0d0) +1.0d0 + dpsi1r = z*z * ( z * (-15.0_dp*z + 32.0_dp) - 18.0_dp) +1.0_dp end function dpsi1 pure function ddpsi1(z) result(ddpsi1r) !$acc routine seq - double precision, intent(in) :: z - double precision :: ddpsi1r + real(dp), intent(in) :: z + real(dp) :: ddpsi1r !$gpu - ddpsi1r = z * (z * (-60.0d0*z + 96.0d0) -36.0d0) + ddpsi1r = z * (z * (-60.0_dp*z + 96.0_dp) -36.0_dp) end function ddpsi1 ! psi2 and its derivatives pure function psi2(z) result(psi2r) !$acc routine seq - double precision, intent(in) :: z - double precision :: psi2r + real(dp), intent(in) :: z + real(dp) :: psi2r !$gpu - psi2r = 0.5d0*z*z*( z* ( z * (-z + 3.0d0) - 3.0d0) + 1.0d0) + psi2r = 0.5_dp*z*z*( z* ( z * (-z + 3.0_dp) - 3.0_dp) + 1.0_dp) end function psi2 pure function dpsi2(z) result(dpsi2r) !$acc routine seq - double precision, intent(in) :: z - double precision :: dpsi2r + real(dp), intent(in) :: z + real(dp) :: dpsi2r !$gpu - dpsi2r = 0.5d0*z*( z*(z*(-5.0d0*z + 12.0d0) - 9.0d0) + 2.0d0) + dpsi2r = 0.5_dp*z*( z*(z*(-5.0_dp*z + 12.0_dp) - 9.0_dp) + 2.0_dp) end function dpsi2 pure function ddpsi2(z) result(ddpsi2r) !$acc routine seq - double precision, intent(in) :: z - double precision :: ddpsi2r + real(dp), intent(in) :: z + real(dp) :: ddpsi2r !$gpu - ddpsi2r = 0.5d0*(z*( z * (-20.0d0*z + 36.0d0) - 18.0d0) + 2.0d0) + ddpsi2r = 0.5_dp*(z*( z * (-20.0_dp*z + 36.0_dp) - 18.0_dp) + 2.0_dp) end function ddpsi2 ! biquintic hermite polynomial function pure function h5(fi,w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md) result(h5r) !$acc routine seq - double precision, intent(in) :: fi(36) - double precision, intent(in) :: w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md - double precision :: h5r + real(dp), intent(in) :: fi(36) + real(dp), intent(in) :: w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md + real(dp) :: h5r !$gpu @@ -1461,46 +1462,46 @@ end function h5 ! psi0 & derivatives pure function xpsi0(z) result(xpsi0r) !$acc routine seq - double precision, intent(in) :: z - double precision :: xpsi0r + real(dp), intent(in) :: z + real(dp) :: xpsi0r !$gpu - xpsi0r = z * z * (2.0d0*z - 3.0d0) + 1.0 + xpsi0r = z * z * (2.0_dp*z - 3.0_dp) + 1.0 end function xpsi0 pure function xdpsi0(z) result(xdpsi0r) !$acc routine seq - double precision, intent(in) :: z - double precision :: xdpsi0r + real(dp), intent(in) :: z + real(dp) :: xdpsi0r !$gpu - xdpsi0r = z * (6.0d0*z - 6.0d0) + xdpsi0r = z * (6.0_dp*z - 6.0_dp) end function xdpsi0 ! psi1 & derivatives pure function xpsi1(z) result(xpsi1r) !$acc routine seq - double precision, intent(in) :: z - double precision :: xpsi1r + real(dp), intent(in) :: z + real(dp) :: xpsi1r !$gpu - xpsi1r = z * ( z * (z - 2.0d0) + 1.0d0) + xpsi1r = z * ( z * (z - 2.0_dp) + 1.0_dp) end function xpsi1 pure function xdpsi1(z) result(xdpsi1r) !$acc routine seq - double precision, intent(in) :: z - double precision :: xdpsi1r + real(dp), intent(in) :: z + real(dp) :: xdpsi1r !$gpu - xdpsi1r = z * (3.0d0*z - 4.0d0) + 1.0d0 + xdpsi1r = z * (3.0_dp*z - 4.0_dp) + 1.0_dp end function xdpsi1 ! bicubic hermite polynomial function pure function h3(fi,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) result(h3r) !$acc routine seq - double precision, intent(in) :: fi(36) - double precision, intent(in) :: w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md - double precision :: h3r + real(dp), intent(in) :: fi(36) + real(dp), intent(in) :: w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md + real(dp) :: h3r !$gpu - h3r = fi(1) *w0d*w0t + fi(2) *w0md*w0t & + h3r = fi(1) *w0d*w0t + fi(2) *w0md*w0t & + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt & + fi(5) *w0d*w1t + fi(6) *w0md*w1t & + fi(7) *w0d*w1mt + fi(8) *w0md*w1mt & diff --git a/tools/starkiller-helmholtz/eos_type.F90 b/tools/starkiller-helmholtz/eos_type.F90 index 69b47c7c..a8d04249 100644 --- a/tools/starkiller-helmholtz/eos_type.F90 +++ b/tools/starkiller-helmholtz/eos_type.F90 @@ -1,15 +1,13 @@ module eos_type_module - use xnet_types, only : rt => dp + use xnet_types, only: dp use xnet_util, only: xnet_terminate implicit none - private :: rt - - double precision, parameter, private :: ONE = 1.0d0 - double precision, parameter, private :: ZERO = 0.0d0 - double precision, parameter, private :: small_x = 0.0d0 + real(dp), parameter, private :: ONE = 1.0_dp + real(dp), parameter, private :: ZERO = 0.0_dp + real(dp), parameter, private :: small_x = 0.0_dp integer, parameter :: eos_input_rt = 1 ! rho, T are inputs integer, parameter :: eos_input_rh = 2 ! rho, h are inputs @@ -45,22 +43,22 @@ module eos_type_module ! Minimum and maximum thermodynamic quantities permitted by the EOS. - real(rt), allocatable :: mintemp - real(rt), allocatable :: maxtemp - real(rt), allocatable :: mindens - real(rt), allocatable :: maxdens - real(rt), allocatable :: minx - real(rt), allocatable :: maxx - real(rt), allocatable :: minye - real(rt), allocatable :: maxye - real(rt), allocatable :: mine - real(rt), allocatable :: maxe - real(rt), allocatable :: minp - real(rt), allocatable :: maxp - real(rt), allocatable :: mins - real(rt), allocatable :: maxs - real(rt), allocatable :: minh - real(rt), allocatable :: maxh + real(dp), allocatable :: mintemp + real(dp), allocatable :: maxtemp + real(dp), allocatable :: mindens + real(dp), allocatable :: maxdens + real(dp), allocatable :: minx + real(dp), allocatable :: maxx + real(dp), allocatable :: minye + real(dp), allocatable :: maxye + real(dp), allocatable :: mine + real(dp), allocatable :: maxe + real(dp), allocatable :: minp + real(dp), allocatable :: maxp + real(dp), allocatable :: mins + real(dp), allocatable :: maxs + real(dp), allocatable :: minh + real(dp), allocatable :: maxh !$acc declare & !$acc create(mintemp, maxtemp, mindens, maxdens, minx, maxx, minye, maxye) & @@ -109,42 +107,42 @@ module eos_type_module type :: eos_t - real(rt) :: rho - real(rt) :: T - real(rt) :: p - real(rt) :: e - real(rt) :: h - real(rt) :: s - - real(rt) :: dpdT - real(rt) :: dpdr - real(rt) :: dedT - real(rt) :: dedr - real(rt) :: dhdT - real(rt) :: dhdr - real(rt) :: dsdT - real(rt) :: dsdr - real(rt) :: dpde - real(rt) :: dpdr_e - - real(rt) :: cv - real(rt) :: cp - real(rt) :: xne - real(rt) :: xnp - real(rt) :: eta - real(rt) :: detadt - real(rt) :: pele - real(rt) :: ppos - real(rt) :: mu - real(rt) :: mu_e - real(rt) :: y_e - real(rt) :: gam1 - real(rt) :: cs - - real(rt) :: abar - real(rt) :: zbar - - real(rt) :: conductivity + real(dp) :: rho + real(dp) :: T + real(dp) :: p + real(dp) :: e + real(dp) :: h + real(dp) :: s + + real(dp) :: dpdT + real(dp) :: dpdr + real(dp) :: dedT + real(dp) :: dedr + real(dp) :: dhdT + real(dp) :: dhdr + real(dp) :: dsdT + real(dp) :: dsdr + real(dp) :: dpde + real(dp) :: dpdr_e + + real(dp) :: cv + real(dp) :: cp + real(dp) :: xne + real(dp) :: xnp + real(dp) :: eta + real(dp) :: detadt + real(dp) :: pele + real(dp) :: ppos + real(dp) :: mu + real(dp) :: mu_e + real(dp) :: y_e + real(dp) :: gam1 + real(dp) :: cs + + real(dp) :: abar + real(dp) :: zbar + + real(dp) :: conductivity end type eos_t @@ -240,7 +238,7 @@ subroutine eos_get_small_temp(small_temp_out) implicit none - real(rt), intent(out) :: small_temp_out + real(dp), intent(out) :: small_temp_out !$gpu @@ -256,7 +254,7 @@ subroutine eos_get_small_dens(small_dens_out) implicit none - real(rt), intent(out) :: small_dens_out + real(dp), intent(out) :: small_dens_out !$gpu @@ -272,7 +270,7 @@ subroutine eos_get_max_temp(max_temp_out) implicit none - real(rt), intent(out) :: max_temp_out + real(dp), intent(out) :: max_temp_out !$gpu @@ -288,7 +286,7 @@ subroutine eos_get_max_dens(max_dens_out) implicit none - real(rt), intent(out) :: max_dens_out + real(dp), intent(out) :: max_dens_out !$gpu From 0d2d035d8d5a28d5ee257187e691f5855f1bdd91 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 10:42:12 -0400 Subject: [PATCH 02/13] remove !$gpu directive --- tools/starkiller-helmholtz/actual_eos.F90 | 18 ------------------ tools/starkiller-helmholtz/eos_type.F90 | 14 -------------- 2 files changed, 32 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 2bb806b2..2478bced 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -252,8 +252,6 @@ subroutine actual_eos(input, state) real(dp) :: smallt, smalld - !$gpu - call eos_get_small_temp(smallt) call eos_get_small_dens(smalld) @@ -1357,7 +1355,6 @@ pure function psi0(z) result(psi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: psi0r - !$gpu psi0r = z**3 * ( z * (-6.0_dp*z + 15.0_dp) -10.0_dp) + 1.0_dp end function psi0 @@ -1365,7 +1362,6 @@ pure function dpsi0(z) result(dpsi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: dpsi0r - !$gpu dpsi0r = z**2 * ( z * (-30.0_dp*z + 60.0_dp) - 30.0_dp) end function dpsi0 @@ -1373,7 +1369,6 @@ pure function ddpsi0(z) result(ddpsi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: ddpsi0r - !$gpu ddpsi0r = z* ( z*( -120.0_dp*z + 180.0_dp) -60.0_dp) end function ddpsi0 @@ -1382,7 +1377,6 @@ pure function psi1(z) result(psi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: psi1r - !$gpu psi1r = z* ( z**2 * ( z * (-3.0_dp*z + 8.0_dp) - 6.0_dp) + 1.0_dp) end function psi1 @@ -1390,7 +1384,6 @@ pure function dpsi1(z) result(dpsi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: dpsi1r - !$gpu dpsi1r = z*z * ( z * (-15.0_dp*z + 32.0_dp) - 18.0_dp) +1.0_dp end function dpsi1 @@ -1398,7 +1391,6 @@ pure function ddpsi1(z) result(ddpsi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: ddpsi1r - !$gpu ddpsi1r = z * (z * (-60.0_dp*z + 96.0_dp) -36.0_dp) end function ddpsi1 @@ -1407,7 +1399,6 @@ pure function psi2(z) result(psi2r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: psi2r - !$gpu psi2r = 0.5_dp*z*z*( z* ( z * (-z + 3.0_dp) - 3.0_dp) + 1.0_dp) end function psi2 @@ -1415,7 +1406,6 @@ pure function dpsi2(z) result(dpsi2r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: dpsi2r - !$gpu dpsi2r = 0.5_dp*z*( z*(z*(-5.0_dp*z + 12.0_dp) - 9.0_dp) + 2.0_dp) end function dpsi2 @@ -1423,7 +1413,6 @@ pure function ddpsi2(z) result(ddpsi2r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: ddpsi2r - !$gpu ddpsi2r = 0.5_dp*(z*( z * (-20.0_dp*z + 36.0_dp) - 18.0_dp) + 2.0_dp) end function ddpsi2 @@ -1435,8 +1424,6 @@ pure function h5(fi,w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md) resul real(dp), intent(in) :: w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md real(dp) :: h5r - !$gpu - h5r = fi(1) *w0d*w0t + fi(2) *w0md*w0t & + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt & + fi(5) *w0d*w1t + fi(6) *w0md*w1t & @@ -1464,7 +1451,6 @@ pure function xpsi0(z) result(xpsi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: xpsi0r - !$gpu xpsi0r = z * z * (2.0_dp*z - 3.0_dp) + 1.0 end function xpsi0 @@ -1472,7 +1458,6 @@ pure function xdpsi0(z) result(xdpsi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: xdpsi0r - !$gpu xdpsi0r = z * (6.0_dp*z - 6.0_dp) end function xdpsi0 @@ -1482,7 +1467,6 @@ pure function xpsi1(z) result(xpsi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: xpsi1r - !$gpu xpsi1r = z * ( z * (z - 2.0_dp) + 1.0_dp) end function xpsi1 @@ -1490,7 +1474,6 @@ pure function xdpsi1(z) result(xdpsi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: xdpsi1r - !$gpu xdpsi1r = z * (3.0_dp*z - 4.0_dp) + 1.0_dp end function xdpsi1 @@ -1500,7 +1483,6 @@ pure function h3(fi,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) result(h3r) real(dp), intent(in) :: fi(36) real(dp), intent(in) :: w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md real(dp) :: h3r - !$gpu h3r = fi(1) *w0d*w0t + fi(2) *w0md*w0t & + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt & + fi(5) *w0d*w1t + fi(6) *w0md*w1t & diff --git a/tools/starkiller-helmholtz/eos_type.F90 b/tools/starkiller-helmholtz/eos_type.F90 index a8d04249..ceeee8ee 100644 --- a/tools/starkiller-helmholtz/eos_type.F90 +++ b/tools/starkiller-helmholtz/eos_type.F90 @@ -156,8 +156,6 @@ subroutine copy_eos_t(to_eos, from_eos) type(eos_t) :: to_eos, from_eos - !$gpu - to_eos % rho = from_eos % rho to_eos % T = from_eos % T to_eos % p = from_eos % p @@ -206,8 +204,6 @@ subroutine clean_state(state) type (eos_t), intent(inout) :: state - !$gpu - state % T = min(maxtemp, max(mintemp, state % T)) state % rho = min(maxdens, max(mindens, state % rho)) @@ -240,8 +236,6 @@ subroutine eos_get_small_temp(small_temp_out) real(dp), intent(out) :: small_temp_out - !$gpu - small_temp_out = mintemp end subroutine eos_get_small_temp @@ -256,8 +250,6 @@ subroutine eos_get_small_dens(small_dens_out) real(dp), intent(out) :: small_dens_out - !$gpu - small_dens_out = mindens end subroutine eos_get_small_dens @@ -272,8 +264,6 @@ subroutine eos_get_max_temp(max_temp_out) real(dp), intent(out) :: max_temp_out - !$gpu - max_temp_out = maxtemp end subroutine eos_get_max_temp @@ -288,8 +278,6 @@ subroutine eos_get_max_dens(max_dens_out) real(dp), intent(out) :: max_dens_out - !$gpu - max_dens_out = maxdens end subroutine eos_get_max_dens @@ -304,8 +292,6 @@ function eos_input_has_var(input, ivar) result(has) integer, intent(in) :: input, ivar logical :: has - !$gpu - has = .false. select case (ivar) From 058d3c78786d38ee6cfbcc4e150536a3949eebf4 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 10:43:59 -0400 Subject: [PATCH 03/13] add implicit none to actual_eos --- tools/starkiller-helmholtz/actual_eos.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 2478bced..6cd7f632 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -3,6 +3,8 @@ module actual_eos_module use xnet_types, only: dp use eos_type_module + implicit none + character (len=64), public :: eos_name = "helmholtz" ! Runtime parameters From 84d90c09580be813f2abcea46d1f041e3d4f12a7 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 10:44:50 -0400 Subject: [PATCH 04/13] remove pure function (for now) --- tools/starkiller-helmholtz/actual_eos.F90 | 30 +++++++++++------------ 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 6cd7f632..788e07ef 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -1353,21 +1353,21 @@ end subroutine actual_eos_init ! quintic hermite polynomial functions ! psi0 and its derivatives - pure function psi0(z) result(psi0r) + function psi0(z) result(psi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: psi0r psi0r = z**3 * ( z * (-6.0_dp*z + 15.0_dp) -10.0_dp) + 1.0_dp end function psi0 - pure function dpsi0(z) result(dpsi0r) + function dpsi0(z) result(dpsi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: dpsi0r dpsi0r = z**2 * ( z * (-30.0_dp*z + 60.0_dp) - 30.0_dp) end function dpsi0 - pure function ddpsi0(z) result(ddpsi0r) + function ddpsi0(z) result(ddpsi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: ddpsi0r @@ -1375,21 +1375,21 @@ pure function ddpsi0(z) result(ddpsi0r) end function ddpsi0 ! psi1 and its derivatives - pure function psi1(z) result(psi1r) + function psi1(z) result(psi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: psi1r psi1r = z* ( z**2 * ( z * (-3.0_dp*z + 8.0_dp) - 6.0_dp) + 1.0_dp) end function psi1 - pure function dpsi1(z) result(dpsi1r) + function dpsi1(z) result(dpsi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: dpsi1r dpsi1r = z*z * ( z * (-15.0_dp*z + 32.0_dp) - 18.0_dp) +1.0_dp end function dpsi1 - pure function ddpsi1(z) result(ddpsi1r) + function ddpsi1(z) result(ddpsi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: ddpsi1r @@ -1397,21 +1397,21 @@ pure function ddpsi1(z) result(ddpsi1r) end function ddpsi1 ! psi2 and its derivatives - pure function psi2(z) result(psi2r) + function psi2(z) result(psi2r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: psi2r psi2r = 0.5_dp*z*z*( z* ( z * (-z + 3.0_dp) - 3.0_dp) + 1.0_dp) end function psi2 - pure function dpsi2(z) result(dpsi2r) + function dpsi2(z) result(dpsi2r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: dpsi2r dpsi2r = 0.5_dp*z*( z*(z*(-5.0_dp*z + 12.0_dp) - 9.0_dp) + 2.0_dp) end function dpsi2 - pure function ddpsi2(z) result(ddpsi2r) + function ddpsi2(z) result(ddpsi2r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: ddpsi2r @@ -1420,7 +1420,7 @@ end function ddpsi2 ! biquintic hermite polynomial function - pure function h5(fi,w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md) result(h5r) + function h5(fi,w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md) result(h5r) !$acc routine seq real(dp), intent(in) :: fi(36) real(dp), intent(in) :: w0t,w1t,w2t,w0mt,w1mt,w2mt,w0d,w1d,w2d,w0md,w1md,w2md @@ -1449,14 +1449,14 @@ end function h5 ! cubic hermite polynomial functions ! psi0 & derivatives - pure function xpsi0(z) result(xpsi0r) + function xpsi0(z) result(xpsi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: xpsi0r xpsi0r = z * z * (2.0_dp*z - 3.0_dp) + 1.0 end function xpsi0 - pure function xdpsi0(z) result(xdpsi0r) + function xdpsi0(z) result(xdpsi0r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: xdpsi0r @@ -1465,14 +1465,14 @@ end function xdpsi0 ! psi1 & derivatives - pure function xpsi1(z) result(xpsi1r) + function xpsi1(z) result(xpsi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: xpsi1r xpsi1r = z * ( z * (z - 2.0_dp) + 1.0_dp) end function xpsi1 - pure function xdpsi1(z) result(xdpsi1r) + function xdpsi1(z) result(xdpsi1r) !$acc routine seq real(dp), intent(in) :: z real(dp) :: xdpsi1r @@ -1480,7 +1480,7 @@ pure function xdpsi1(z) result(xdpsi1r) end function xdpsi1 ! bicubic hermite polynomial function - pure function h3(fi,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) result(h3r) + function h3(fi,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) result(h3r) !$acc routine seq real(dp), intent(in) :: fi(36) real(dp), intent(in) :: w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md From 8b85679ffa2c339e6880e87a1b42684e83f3fdcb Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 10:57:30 -0400 Subject: [PATCH 05/13] make starkiller modules private by default --- tools/starkiller-helmholtz/actual_eos.F90 | 15 ++-- tools/starkiller-helmholtz/eos_type.F90 | 94 +++++++++++------------ 2 files changed, 53 insertions(+), 56 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 788e07ef..1fd0588c 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -4,6 +4,7 @@ module actual_eos_module use eos_type_module implicit none + private character (len=64), public :: eos_name = "helmholtz" @@ -11,8 +12,6 @@ module actual_eos_module logical, allocatable :: do_coulomb logical, allocatable :: input_is_constant - real(dp), parameter, private :: ZERO = 0.0_dp, HALF = 0.5_dp, TWO = 2.0_dp - !..for the tables, in general integer, parameter, private :: imax = 541, jmax = 201 integer, allocatable :: itmax, jtmax @@ -47,7 +46,7 @@ module actual_eos_module dd_sav(:), dd2_sav(:), & ddi_sav(:), dd2i_sav(:) - integer, parameter :: max_newton = 100 + integer, parameter :: max_newton = 100 ! 2006 CODATA physical constants @@ -110,15 +109,13 @@ module actual_eos_module !$acc create(dd_sav, dd2_sav, ddi_sav, dd2i_sav) & !$acc create(do_coulomb, input_is_constant) -public actual_eos, actual_eos_init, actual_eos_finalize, eos_supports_input_type + public :: actual_eos, actual_eos_init, actual_eos_finalize, eos_supports_input_type contains function eos_supports_input_type(input) result(supported) - use eos_type_module - implicit none integer, intent(in) :: input @@ -806,7 +803,7 @@ subroutine actual_eos(input, state) p_temp = prad + pion + pele + pcoul e_temp = erad + eion + eele + ecoul - if (p_temp .le. ZERO .or. e_temp .le. ZERO) then + if (p_temp .le. 0.0_dp .or. e_temp .le. 0.0_dp) then pcoul = 0.0_dp dpcouldd = 0.0_dp @@ -1040,8 +1037,8 @@ subroutine actual_eos(input, state) ! Don't let the temperature or density change by more ! than a factor of two - tnew = max(HALF * told, min(tnew, TWO * told)) - rnew = max(HALF * rold, min(rnew, TWO * rold)) + tnew = max(0.5_dp * told, min(tnew, 2.0_dp * told)) + rnew = max(0.5_dp * rold, min(rnew, 2.0_dp * rold)) ! Don't let us freeze or evacuate tnew = max(smallt, tnew) diff --git a/tools/starkiller-helmholtz/eos_type.F90 b/tools/starkiller-helmholtz/eos_type.F90 index ceeee8ee..32ce3e9c 100644 --- a/tools/starkiller-helmholtz/eos_type.F90 +++ b/tools/starkiller-helmholtz/eos_type.F90 @@ -4,66 +4,66 @@ module eos_type_module use xnet_util, only: xnet_terminate implicit none + private - real(dp), parameter, private :: ONE = 1.0_dp - real(dp), parameter, private :: ZERO = 0.0_dp - real(dp), parameter, private :: small_x = 0.0_dp - - integer, parameter :: eos_input_rt = 1 ! rho, T are inputs - integer, parameter :: eos_input_rh = 2 ! rho, h are inputs - integer, parameter :: eos_input_tp = 3 ! T, p are inputs - integer, parameter :: eos_input_rp = 4 ! rho, p are inputs - integer, parameter :: eos_input_re = 5 ! rho, e are inputs - integer, parameter :: eos_input_ps = 6 ! p, s are inputs - integer, parameter :: eos_input_ph = 7 ! p, h are inputs - integer, parameter :: eos_input_th = 8 ! T, h are inputs + integer, parameter, public :: eos_input_rt = 1 ! rho, T are inputs + integer, parameter, public :: eos_input_rh = 2 ! rho, h are inputs + integer, parameter, public :: eos_input_tp = 3 ! T, p are inputs + integer, parameter, public :: eos_input_rp = 4 ! rho, p are inputs + integer, parameter, public :: eos_input_re = 5 ! rho, e are inputs + integer, parameter, public :: eos_input_ps = 6 ! p, s are inputs + integer, parameter, public :: eos_input_ph = 7 ! p, h are inputs + integer, parameter, public :: eos_input_th = 8 ! T, h are inputs ! these are used to allow for a generic interface to the ! root finding - integer, parameter :: itemp = 1 - integer, parameter :: idens = 2 - integer, parameter :: iener = 3 - integer, parameter :: ienth = 4 - integer, parameter :: ientr = 5 - integer, parameter :: ipres = 6 + integer, parameter, public :: itemp = 1 + integer, parameter, public :: idens = 2 + integer, parameter, public :: iener = 3 + integer, parameter, public :: ienth = 4 + integer, parameter, public :: ientr = 5 + integer, parameter, public :: ipres = 6 ! error codes - integer, parameter :: ierr_general = 1 - integer, parameter :: ierr_input = 2 - integer, parameter :: ierr_iter_conv = 3 - integer, parameter :: ierr_neg_e = 4 - integer, parameter :: ierr_neg_p = 5 - integer, parameter :: ierr_neg_h = 6 - integer, parameter :: ierr_neg_s = 7 - integer, parameter :: ierr_iter_var = 8 - integer, parameter :: ierr_init = 9 - integer, parameter :: ierr_init_xn = 10 - integer, parameter :: ierr_out_of_bounds = 11 - integer, parameter :: ierr_not_implemented = 12 + integer, parameter, public :: ierr_general = 1 + integer, parameter, public :: ierr_input = 2 + integer, parameter, public :: ierr_iter_conv = 3 + integer, parameter, public :: ierr_neg_e = 4 + integer, parameter, public :: ierr_neg_p = 5 + integer, parameter, public :: ierr_neg_h = 6 + integer, parameter, public :: ierr_neg_s = 7 + integer, parameter, public :: ierr_iter_var = 8 + integer, parameter, public :: ierr_init = 9 + integer, parameter, public :: ierr_init_xn = 10 + integer, parameter, public :: ierr_out_of_bounds = 11 + integer, parameter, public :: ierr_not_implemented = 12 ! Minimum and maximum thermodynamic quantities permitted by the EOS. - real(dp), allocatable :: mintemp - real(dp), allocatable :: maxtemp - real(dp), allocatable :: mindens - real(dp), allocatable :: maxdens - real(dp), allocatable :: minx - real(dp), allocatable :: maxx - real(dp), allocatable :: minye - real(dp), allocatable :: maxye - real(dp), allocatable :: mine - real(dp), allocatable :: maxe - real(dp), allocatable :: minp - real(dp), allocatable :: maxp - real(dp), allocatable :: mins - real(dp), allocatable :: maxs - real(dp), allocatable :: minh - real(dp), allocatable :: maxh + real(dp), allocatable, public :: mintemp + real(dp), allocatable, public :: maxtemp + real(dp), allocatable, public :: mindens + real(dp), allocatable, public :: maxdens + real(dp), allocatable, public :: minx + real(dp), allocatable, public :: maxx + real(dp), allocatable, public :: minye + real(dp), allocatable, public :: maxye + real(dp), allocatable, public :: mine + real(dp), allocatable, public :: maxe + real(dp), allocatable, public :: minp + real(dp), allocatable, public :: maxp + real(dp), allocatable, public :: mins + real(dp), allocatable, public :: maxs + real(dp), allocatable, public :: minh + real(dp), allocatable, public :: maxh !$acc declare & !$acc create(mintemp, maxtemp, mindens, maxdens, minx, maxx, minye, maxye) & !$acc create(mine, maxe, minp, maxp, mins, maxs, minh, maxh) + public :: clean_state, print_state, eos_get_small_temp, eos_get_small_dens + public :: eos_get_max_temp, eos_get_max_dens, eos_input_has_var + ! A generic structure holding thermodynamic quantities and their derivatives, ! plus some other quantities of interest. @@ -105,7 +105,7 @@ module eos_type_module ! dpdr_e -- d pressure / d rho |_energy ! conductivity -- thermal conductivity (in erg/cm/K/sec) - type :: eos_t + type, public :: eos_t real(dp) :: rho real(dp) :: T From f24033ca6fc0bdb582af76e81dabb68145bdcbb3 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 13:04:58 -0400 Subject: [PATCH 06/13] comment portions of helmholtz dedicated to iterations (e.g. temp from energy) -- all of our calls are eos_input_rt (dens,temp) --- tools/starkiller-helmholtz/actual_eos.F90 | 452 +++++++++++----------- 1 file changed, 226 insertions(+), 226 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 1fd0588c..3116aa45 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -249,10 +249,10 @@ subroutine actual_eos(input, state) real(dp) :: p_temp, e_temp - real(dp) :: smallt, smalld + !real(dp) :: smallt, smalld - call eos_get_small_temp(smallt) - call eos_get_small_dens(smalld) + !call eos_get_small_temp(smallt) + !call eos_get_small_dens(smalld) temp_row = state % T den_row = state % rho @@ -262,65 +262,65 @@ subroutine actual_eos(input, state) ! Initial setup for iterations - single_iter = .false. - double_iter = .false. + !single_iter = .false. + !double_iter = .false. - if (input .eq. eos_input_rt) then + !if (input .eq. eos_input_rt) then - ! Nothing to do here. + ! ! Nothing to do here. - elseif (input .eq. eos_input_rh) then + !elseif (input .eq. eos_input_rh) then - single_iter = .true. - v_want = state % h - var = ienth - dvar = itemp + ! single_iter = .true. + ! v_want = state % h + ! var = ienth + ! dvar = itemp - elseif (input .eq. eos_input_tp) then + !elseif (input .eq. eos_input_tp) then - single_iter = .true. - v_want = state % p - var = ipres - dvar = idens + ! single_iter = .true. + ! v_want = state % p + ! var = ipres + ! dvar = idens - elseif (input .eq. eos_input_rp) then + !elseif (input .eq. eos_input_rp) then - single_iter = .true. - v_want = state % p - var = ipres - dvar = itemp + ! single_iter = .true. + ! v_want = state % p + ! var = ipres + ! dvar = itemp - elseif (input .eq. eos_input_re) then + !elseif (input .eq. eos_input_re) then - single_iter = .true. - v_want = state % e - var = iener - dvar = itemp + ! single_iter = .true. + ! v_want = state % e + ! var = iener + ! dvar = itemp - elseif (input .eq. eos_input_ps) then + !elseif (input .eq. eos_input_ps) then - double_iter = .true. - v1_want = state % p - v2_want = state % s - var1 = ipres - var2 = ientr + ! double_iter = .true. + ! v1_want = state % p + ! v2_want = state % s + ! var1 = ipres + ! var2 = ientr - elseif (input .eq. eos_input_ph) then + !elseif (input .eq. eos_input_ph) then - double_iter = .true. - v1_want = state % p - v2_want = state % h - var1 = ipres - var2 = ienth + ! double_iter = .true. + ! v1_want = state % p + ! v2_want = state % h + ! var1 = ipres + ! var2 = ienth - elseif (input .eq. eos_input_th) then + !elseif (input .eq. eos_input_th) then - single_iter = .true. - v_want = state % h - var = ienth - dvar = idens + ! single_iter = .true. + ! v_want = state % h + ! var = ienth + ! dvar = idens - endif + !endif ptot_row = 0.0_dp dpt_row = 0.0_dp @@ -358,11 +358,11 @@ subroutine actual_eos(input, state) cs_row = 0.0_dp gam1_row = 0.0_dp - converged = .false. + !converged = .false. - if (input .eq. eos_input_rt) converged = .true. + !if (input .eq. eos_input_rt) converged = .true. - do iter = 1, max_newton + !do iter = 1, max_newton temp = temp_row den = den_row @@ -896,167 +896,167 @@ subroutine actual_eos(input, state) cs_row = sound gam1_row = gam1 - if (converged) then - - exit - - elseif (single_iter) then - - if (dvar .eq. itemp) then - - x = temp_row - smallx = smallt - xtol = ttol - - if (var .eq. ipres) then - v = ptot_row - dvdx = dpt_row - elseif (var .eq. iener) then - v = etot_row - dvdx = det_row - elseif (var .eq. ientr) then - v = stot_row - dvdx = dst_row - elseif (var .eq. ienth) then - v = htot_row - dvdx = dht_row - else - exit - endif - - else ! dvar == density - - x = den_row - smallx = smalld - xtol = dtol - - if (var .eq. ipres) then - v = ptot_row - dvdx = dpd_row - elseif (var .eq. iener) then - v = etot_row - dvdx = ded_row - elseif (var .eq. ientr) then - v = stot_row - dvdx = dsd_row - elseif (var .eq. ienth) then - v = htot_row - dvdx = dhd_row - else - exit - endif - - endif - - ! Now do the calculation for the next guess for T/rho - - xnew = x - (v - v_want) / dvdx - - ! Don't let the temperature/density change by more than a factor of two - xnew = max(0.5 * x, min(xnew, 2.0 * x)) - - ! Don't let us freeze/evacuate - xnew = max(smallx, xnew) - - ! Store the new temperature/density - - if (dvar .eq. itemp) then - temp_row = xnew - else - den_row = xnew - endif - - ! Compute the error from the last iteration - - error = abs( (xnew - x) / x ) - - if (error .lt. xtol) converged = .true. - - elseif (double_iter) then - - ! Figure out which variables we're using - - told = temp_row - rold = den_row - - if (var1 .eq. ipres) then - v1 = ptot_row - dv1dt = dpt_row - dv1dr = dpd_row - elseif (var1 .eq. iener) then - v1 = etot_row - dv1dt = det_row - dv1dr = ded_row - elseif (var1 .eq. ientr) then - v1 = stot_row - dv1dt = dst_row - dv1dr = dsd_row - elseif (var1 .eq. ienth) then - v1 = htot_row - dv1dt = dht_row - dv1dr = dhd_row - else - exit - endif - - if (var2 .eq. ipres) then - v2 = ptot_row - dv2dt = dpt_row - dv2dr = dpd_row - elseif (var2 .eq. iener) then - v2 = etot_row - dv2dt = det_row - dv2dr = ded_row - elseif (var2 .eq. ientr) then - v2 = stot_row - dv2dt = dst_row - dv2dr = dsd_row - elseif (var2 .eq. ienth) then - v2 = htot_row - dv2dt = dht_row - dv2dr = dhd_row - else - exit - endif - - ! Two functions, f and g, to iterate over - v1i = v1_want - v1 - v2i = v2_want - v2 - - ! - ! 0 = f + dfdr * delr + dfdt * delt - ! 0 = g + dgdr * delr + dgdt * delt - ! - - ! note that dfi/dT = - df/dT - delr = (-v1i*dv2dt + v2i*dv1dt) / (dv2dr*dv1dt - dv2dt*dv1dr) - - rnew = rold + delr - - tnew = told + (v1i - dv1dr*delr) / dv1dt - - ! Don't let the temperature or density change by more - ! than a factor of two - tnew = max(0.5_dp * told, min(tnew, 2.0_dp * told)) - rnew = max(0.5_dp * rold, min(rnew, 2.0_dp * rold)) - - ! Don't let us freeze or evacuate - tnew = max(smallt, tnew) - rnew = max(smalld, rnew) - - ! Store the new temperature and density - den_row = rnew - temp_row = tnew - - ! Compute the errors - error1 = abs( (rnew - rold) / rold ) - error2 = abs( (tnew - told) / told ) - - if (error1 .LT. dtol .and. error2 .LT. ttol) converged = .true. + !if (converged) then + + ! exit + + !elseif (single_iter) then + + ! if (dvar .eq. itemp) then + + ! x = temp_row + ! smallx = smallt + ! xtol = ttol + + ! if (var .eq. ipres) then + ! v = ptot_row + ! dvdx = dpt_row + ! elseif (var .eq. iener) then + ! v = etot_row + ! dvdx = det_row + ! elseif (var .eq. ientr) then + ! v = stot_row + ! dvdx = dst_row + ! elseif (var .eq. ienth) then + ! v = htot_row + ! dvdx = dht_row + ! else + ! exit + ! endif + + ! else ! dvar == density + + ! x = den_row + ! smallx = smalld + ! xtol = dtol + + ! if (var .eq. ipres) then + ! v = ptot_row + ! dvdx = dpd_row + ! elseif (var .eq. iener) then + ! v = etot_row + ! dvdx = ded_row + ! elseif (var .eq. ientr) then + ! v = stot_row + ! dvdx = dsd_row + ! elseif (var .eq. ienth) then + ! v = htot_row + ! dvdx = dhd_row + ! else + ! exit + ! endif + + ! endif + + ! ! Now do the calculation for the next guess for T/rho + + ! xnew = x - (v - v_want) / dvdx + + ! ! Don't let the temperature/density change by more than a factor of two + ! xnew = max(0.5 * x, min(xnew, 2.0 * x)) + + ! ! Don't let us freeze/evacuate + ! xnew = max(smallx, xnew) + + ! ! Store the new temperature/density + + ! if (dvar .eq. itemp) then + ! temp_row = xnew + ! else + ! den_row = xnew + ! endif + + ! ! Compute the error from the last iteration + + ! error = abs( (xnew - x) / x ) + + ! if (error .lt. xtol) converged = .true. + + !elseif (double_iter) then + + ! ! Figure out which variables we're using + + ! told = temp_row + ! rold = den_row + + ! if (var1 .eq. ipres) then + ! v1 = ptot_row + ! dv1dt = dpt_row + ! dv1dr = dpd_row + ! elseif (var1 .eq. iener) then + ! v1 = etot_row + ! dv1dt = det_row + ! dv1dr = ded_row + ! elseif (var1 .eq. ientr) then + ! v1 = stot_row + ! dv1dt = dst_row + ! dv1dr = dsd_row + ! elseif (var1 .eq. ienth) then + ! v1 = htot_row + ! dv1dt = dht_row + ! dv1dr = dhd_row + ! else + ! exit + ! endif + + ! if (var2 .eq. ipres) then + ! v2 = ptot_row + ! dv2dt = dpt_row + ! dv2dr = dpd_row + ! elseif (var2 .eq. iener) then + ! v2 = etot_row + ! dv2dt = det_row + ! dv2dr = ded_row + ! elseif (var2 .eq. ientr) then + ! v2 = stot_row + ! dv2dt = dst_row + ! dv2dr = dsd_row + ! elseif (var2 .eq. ienth) then + ! v2 = htot_row + ! dv2dt = dht_row + ! dv2dr = dhd_row + ! else + ! exit + ! endif + + ! ! Two functions, f and g, to iterate over + ! v1i = v1_want - v1 + ! v2i = v2_want - v2 + + ! ! + ! ! 0 = f + dfdr * delr + dfdt * delt + ! ! 0 = g + dgdr * delr + dgdt * delt + ! ! + + ! ! note that dfi/dT = - df/dT + ! delr = (-v1i*dv2dt + v2i*dv1dt) / (dv2dr*dv1dt - dv2dt*dv1dr) + + ! rnew = rold + delr + + ! tnew = told + (v1i - dv1dr*delr) / dv1dt + + ! ! Don't let the temperature or density change by more + ! ! than a factor of two + ! tnew = max(0.5_dp * told, min(tnew, 2.0_dp * told)) + ! rnew = max(0.5_dp * rold, min(rnew, 2.0_dp * rold)) + + ! ! Don't let us freeze or evacuate + ! tnew = max(smallt, tnew) + ! rnew = max(smalld, rnew) - endif + ! ! Store the new temperature and density + ! den_row = rnew + ! temp_row = tnew - enddo + ! ! Compute the errors + ! error1 = abs( (rnew - rold) / rold ) + ! error2 = abs( (tnew - told) / told ) + + ! if (error1 .LT. dtol .and. error2 .LT. ttol) converged = .true. + + !endif + + !enddo state % T = temp_row state % rho = den_row @@ -1107,41 +1107,41 @@ subroutine actual_eos(input, state) state % cs = sqrt(state % gam1 * state % p / state % rho) - if (input_is_constant) then + !if (input_is_constant) then - if (input .eq. eos_input_rh) then + ! if (input .eq. eos_input_rh) then - state % h = v_want + ! state % h = v_want - elseif (input .eq. eos_input_tp) then + ! elseif (input .eq. eos_input_tp) then - state % p = v_want + ! state % p = v_want - elseif (input .eq. eos_input_rp) then + ! elseif (input .eq. eos_input_rp) then - state % p = v_want + ! state % p = v_want - elseif (input .eq. eos_input_re) then + ! elseif (input .eq. eos_input_re) then - state % e = v_want + ! state % e = v_want - elseif (input .eq. eos_input_ps) then + ! elseif (input .eq. eos_input_ps) then - state % p = v1_want - state % s = v2_want + ! state % p = v1_want + ! state % s = v2_want - elseif (input .eq. eos_input_ph) then + ! elseif (input .eq. eos_input_ph) then - state % p = v1_want - state % h = v2_want + ! state % p = v1_want + ! state % h = v2_want - elseif (input .eq. eos_input_th) then + ! elseif (input .eq. eos_input_th) then - state % h = v_want + ! state % h = v_want - endif + ! endif - endif + !endif end subroutine actual_eos From 8192e47bc24cbc8a61226241680ae5a21d568b8e Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 13:13:44 -0400 Subject: [PATCH 07/13] comment branching for coulumb corrections --- assume they are always on --- tools/starkiller-helmholtz/actual_eos.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 3116aa45..2c9baf28 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -741,7 +741,7 @@ subroutine actual_eos(input, state) plasgdz = 2.0_dp * plasg/zbar ! TURN ON/OFF COULOMB - if ( do_coulomb ) then + !if ( do_coulomb ) then !...yakovlev & shalybkov 1989 equations 82, 85, 86, 87 if (plasg .ge. 1.0_dp) then x = plasg**(0.25_dp) @@ -822,7 +822,7 @@ subroutine actual_eos(input, state) dscouldz = 0.0_dp end if - end if + !end if !..sum all the components pres = prad + pion + pele + pcoul From a280009210163f18b5b02a34b5c6a69bce3e85c8 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 13:17:24 -0400 Subject: [PATCH 08/13] Revert "comment branching for coulumb corrections --- assume they are always on" This reverts commit 8192e47bc24cbc8a61226241680ae5a21d568b8e. --- tools/starkiller-helmholtz/actual_eos.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 2c9baf28..3116aa45 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -741,7 +741,7 @@ subroutine actual_eos(input, state) plasgdz = 2.0_dp * plasg/zbar ! TURN ON/OFF COULOMB - !if ( do_coulomb ) then + if ( do_coulomb ) then !...yakovlev & shalybkov 1989 equations 82, 85, 86, 87 if (plasg .ge. 1.0_dp) then x = plasg**(0.25_dp) @@ -822,7 +822,7 @@ subroutine actual_eos(input, state) dscouldz = 0.0_dp end if - !end if + end if !..sum all the components pres = prad + pion + pele + pcoul From 98820747adfe0824bdfcbf2183088ead760533ef Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 13:17:33 -0400 Subject: [PATCH 09/13] Revert "comment portions of helmholtz dedicated to iterations (e.g. temp from energy) -- all of our calls are eos_input_rt (dens,temp)" This reverts commit f24033ca6fc0bdb582af76e81dabb68145bdcbb3. --- tools/starkiller-helmholtz/actual_eos.F90 | 452 +++++++++++----------- 1 file changed, 226 insertions(+), 226 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 3116aa45..1fd0588c 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -249,10 +249,10 @@ subroutine actual_eos(input, state) real(dp) :: p_temp, e_temp - !real(dp) :: smallt, smalld + real(dp) :: smallt, smalld - !call eos_get_small_temp(smallt) - !call eos_get_small_dens(smalld) + call eos_get_small_temp(smallt) + call eos_get_small_dens(smalld) temp_row = state % T den_row = state % rho @@ -262,65 +262,65 @@ subroutine actual_eos(input, state) ! Initial setup for iterations - !single_iter = .false. - !double_iter = .false. + single_iter = .false. + double_iter = .false. - !if (input .eq. eos_input_rt) then + if (input .eq. eos_input_rt) then - ! ! Nothing to do here. + ! Nothing to do here. - !elseif (input .eq. eos_input_rh) then + elseif (input .eq. eos_input_rh) then - ! single_iter = .true. - ! v_want = state % h - ! var = ienth - ! dvar = itemp + single_iter = .true. + v_want = state % h + var = ienth + dvar = itemp - !elseif (input .eq. eos_input_tp) then + elseif (input .eq. eos_input_tp) then - ! single_iter = .true. - ! v_want = state % p - ! var = ipres - ! dvar = idens + single_iter = .true. + v_want = state % p + var = ipres + dvar = idens - !elseif (input .eq. eos_input_rp) then + elseif (input .eq. eos_input_rp) then - ! single_iter = .true. - ! v_want = state % p - ! var = ipres - ! dvar = itemp + single_iter = .true. + v_want = state % p + var = ipres + dvar = itemp - !elseif (input .eq. eos_input_re) then + elseif (input .eq. eos_input_re) then - ! single_iter = .true. - ! v_want = state % e - ! var = iener - ! dvar = itemp + single_iter = .true. + v_want = state % e + var = iener + dvar = itemp - !elseif (input .eq. eos_input_ps) then + elseif (input .eq. eos_input_ps) then - ! double_iter = .true. - ! v1_want = state % p - ! v2_want = state % s - ! var1 = ipres - ! var2 = ientr + double_iter = .true. + v1_want = state % p + v2_want = state % s + var1 = ipres + var2 = ientr - !elseif (input .eq. eos_input_ph) then + elseif (input .eq. eos_input_ph) then - ! double_iter = .true. - ! v1_want = state % p - ! v2_want = state % h - ! var1 = ipres - ! var2 = ienth + double_iter = .true. + v1_want = state % p + v2_want = state % h + var1 = ipres + var2 = ienth - !elseif (input .eq. eos_input_th) then + elseif (input .eq. eos_input_th) then - ! single_iter = .true. - ! v_want = state % h - ! var = ienth - ! dvar = idens + single_iter = .true. + v_want = state % h + var = ienth + dvar = idens - !endif + endif ptot_row = 0.0_dp dpt_row = 0.0_dp @@ -358,11 +358,11 @@ subroutine actual_eos(input, state) cs_row = 0.0_dp gam1_row = 0.0_dp - !converged = .false. + converged = .false. - !if (input .eq. eos_input_rt) converged = .true. + if (input .eq. eos_input_rt) converged = .true. - !do iter = 1, max_newton + do iter = 1, max_newton temp = temp_row den = den_row @@ -896,167 +896,167 @@ subroutine actual_eos(input, state) cs_row = sound gam1_row = gam1 - !if (converged) then - - ! exit - - !elseif (single_iter) then - - ! if (dvar .eq. itemp) then - - ! x = temp_row - ! smallx = smallt - ! xtol = ttol - - ! if (var .eq. ipres) then - ! v = ptot_row - ! dvdx = dpt_row - ! elseif (var .eq. iener) then - ! v = etot_row - ! dvdx = det_row - ! elseif (var .eq. ientr) then - ! v = stot_row - ! dvdx = dst_row - ! elseif (var .eq. ienth) then - ! v = htot_row - ! dvdx = dht_row - ! else - ! exit - ! endif - - ! else ! dvar == density - - ! x = den_row - ! smallx = smalld - ! xtol = dtol - - ! if (var .eq. ipres) then - ! v = ptot_row - ! dvdx = dpd_row - ! elseif (var .eq. iener) then - ! v = etot_row - ! dvdx = ded_row - ! elseif (var .eq. ientr) then - ! v = stot_row - ! dvdx = dsd_row - ! elseif (var .eq. ienth) then - ! v = htot_row - ! dvdx = dhd_row - ! else - ! exit - ! endif - - ! endif - - ! ! Now do the calculation for the next guess for T/rho - - ! xnew = x - (v - v_want) / dvdx - - ! ! Don't let the temperature/density change by more than a factor of two - ! xnew = max(0.5 * x, min(xnew, 2.0 * x)) - - ! ! Don't let us freeze/evacuate - ! xnew = max(smallx, xnew) - - ! ! Store the new temperature/density - - ! if (dvar .eq. itemp) then - ! temp_row = xnew - ! else - ! den_row = xnew - ! endif - - ! ! Compute the error from the last iteration - - ! error = abs( (xnew - x) / x ) - - ! if (error .lt. xtol) converged = .true. - - !elseif (double_iter) then - - ! ! Figure out which variables we're using - - ! told = temp_row - ! rold = den_row - - ! if (var1 .eq. ipres) then - ! v1 = ptot_row - ! dv1dt = dpt_row - ! dv1dr = dpd_row - ! elseif (var1 .eq. iener) then - ! v1 = etot_row - ! dv1dt = det_row - ! dv1dr = ded_row - ! elseif (var1 .eq. ientr) then - ! v1 = stot_row - ! dv1dt = dst_row - ! dv1dr = dsd_row - ! elseif (var1 .eq. ienth) then - ! v1 = htot_row - ! dv1dt = dht_row - ! dv1dr = dhd_row - ! else - ! exit - ! endif - - ! if (var2 .eq. ipres) then - ! v2 = ptot_row - ! dv2dt = dpt_row - ! dv2dr = dpd_row - ! elseif (var2 .eq. iener) then - ! v2 = etot_row - ! dv2dt = det_row - ! dv2dr = ded_row - ! elseif (var2 .eq. ientr) then - ! v2 = stot_row - ! dv2dt = dst_row - ! dv2dr = dsd_row - ! elseif (var2 .eq. ienth) then - ! v2 = htot_row - ! dv2dt = dht_row - ! dv2dr = dhd_row - ! else - ! exit - ! endif - - ! ! Two functions, f and g, to iterate over - ! v1i = v1_want - v1 - ! v2i = v2_want - v2 - - ! ! - ! ! 0 = f + dfdr * delr + dfdt * delt - ! ! 0 = g + dgdr * delr + dgdt * delt - ! ! - - ! ! note that dfi/dT = - df/dT - ! delr = (-v1i*dv2dt + v2i*dv1dt) / (dv2dr*dv1dt - dv2dt*dv1dr) - - ! rnew = rold + delr - - ! tnew = told + (v1i - dv1dr*delr) / dv1dt - - ! ! Don't let the temperature or density change by more - ! ! than a factor of two - ! tnew = max(0.5_dp * told, min(tnew, 2.0_dp * told)) - ! rnew = max(0.5_dp * rold, min(rnew, 2.0_dp * rold)) - - ! ! Don't let us freeze or evacuate - ! tnew = max(smallt, tnew) - ! rnew = max(smalld, rnew) + if (converged) then + + exit + + elseif (single_iter) then + + if (dvar .eq. itemp) then + + x = temp_row + smallx = smallt + xtol = ttol + + if (var .eq. ipres) then + v = ptot_row + dvdx = dpt_row + elseif (var .eq. iener) then + v = etot_row + dvdx = det_row + elseif (var .eq. ientr) then + v = stot_row + dvdx = dst_row + elseif (var .eq. ienth) then + v = htot_row + dvdx = dht_row + else + exit + endif + + else ! dvar == density + + x = den_row + smallx = smalld + xtol = dtol + + if (var .eq. ipres) then + v = ptot_row + dvdx = dpd_row + elseif (var .eq. iener) then + v = etot_row + dvdx = ded_row + elseif (var .eq. ientr) then + v = stot_row + dvdx = dsd_row + elseif (var .eq. ienth) then + v = htot_row + dvdx = dhd_row + else + exit + endif + + endif + + ! Now do the calculation for the next guess for T/rho + + xnew = x - (v - v_want) / dvdx + + ! Don't let the temperature/density change by more than a factor of two + xnew = max(0.5 * x, min(xnew, 2.0 * x)) + + ! Don't let us freeze/evacuate + xnew = max(smallx, xnew) + + ! Store the new temperature/density + + if (dvar .eq. itemp) then + temp_row = xnew + else + den_row = xnew + endif + + ! Compute the error from the last iteration + + error = abs( (xnew - x) / x ) + + if (error .lt. xtol) converged = .true. + + elseif (double_iter) then + + ! Figure out which variables we're using + + told = temp_row + rold = den_row + + if (var1 .eq. ipres) then + v1 = ptot_row + dv1dt = dpt_row + dv1dr = dpd_row + elseif (var1 .eq. iener) then + v1 = etot_row + dv1dt = det_row + dv1dr = ded_row + elseif (var1 .eq. ientr) then + v1 = stot_row + dv1dt = dst_row + dv1dr = dsd_row + elseif (var1 .eq. ienth) then + v1 = htot_row + dv1dt = dht_row + dv1dr = dhd_row + else + exit + endif + + if (var2 .eq. ipres) then + v2 = ptot_row + dv2dt = dpt_row + dv2dr = dpd_row + elseif (var2 .eq. iener) then + v2 = etot_row + dv2dt = det_row + dv2dr = ded_row + elseif (var2 .eq. ientr) then + v2 = stot_row + dv2dt = dst_row + dv2dr = dsd_row + elseif (var2 .eq. ienth) then + v2 = htot_row + dv2dt = dht_row + dv2dr = dhd_row + else + exit + endif + + ! Two functions, f and g, to iterate over + v1i = v1_want - v1 + v2i = v2_want - v2 + + ! + ! 0 = f + dfdr * delr + dfdt * delt + ! 0 = g + dgdr * delr + dgdt * delt + ! + + ! note that dfi/dT = - df/dT + delr = (-v1i*dv2dt + v2i*dv1dt) / (dv2dr*dv1dt - dv2dt*dv1dr) + + rnew = rold + delr + + tnew = told + (v1i - dv1dr*delr) / dv1dt + + ! Don't let the temperature or density change by more + ! than a factor of two + tnew = max(0.5_dp * told, min(tnew, 2.0_dp * told)) + rnew = max(0.5_dp * rold, min(rnew, 2.0_dp * rold)) + + ! Don't let us freeze or evacuate + tnew = max(smallt, tnew) + rnew = max(smalld, rnew) + + ! Store the new temperature and density + den_row = rnew + temp_row = tnew + + ! Compute the errors + error1 = abs( (rnew - rold) / rold ) + error2 = abs( (tnew - told) / told ) + + if (error1 .LT. dtol .and. error2 .LT. ttol) converged = .true. - ! ! Store the new temperature and density - ! den_row = rnew - ! temp_row = tnew + endif - ! ! Compute the errors - ! error1 = abs( (rnew - rold) / rold ) - ! error2 = abs( (tnew - told) / told ) - - ! if (error1 .LT. dtol .and. error2 .LT. ttol) converged = .true. - - !endif - - !enddo + enddo state % T = temp_row state % rho = den_row @@ -1107,41 +1107,41 @@ subroutine actual_eos(input, state) state % cs = sqrt(state % gam1 * state % p / state % rho) - !if (input_is_constant) then + if (input_is_constant) then - ! if (input .eq. eos_input_rh) then + if (input .eq. eos_input_rh) then - ! state % h = v_want + state % h = v_want - ! elseif (input .eq. eos_input_tp) then + elseif (input .eq. eos_input_tp) then - ! state % p = v_want + state % p = v_want - ! elseif (input .eq. eos_input_rp) then + elseif (input .eq. eos_input_rp) then - ! state % p = v_want + state % p = v_want - ! elseif (input .eq. eos_input_re) then + elseif (input .eq. eos_input_re) then - ! state % e = v_want + state % e = v_want - ! elseif (input .eq. eos_input_ps) then + elseif (input .eq. eos_input_ps) then - ! state % p = v1_want - ! state % s = v2_want + state % p = v1_want + state % s = v2_want - ! elseif (input .eq. eos_input_ph) then + elseif (input .eq. eos_input_ph) then - ! state % p = v1_want - ! state % h = v2_want + state % p = v1_want + state % h = v2_want - ! elseif (input .eq. eos_input_th) then + elseif (input .eq. eos_input_th) then - ! state % h = v_want + state % h = v_want - ! endif + endif - !endif + endif end subroutine actual_eos From 95a60019053285ba1c606675103fa8b8bc8e7829 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 13:24:44 -0400 Subject: [PATCH 10/13] Create pared down version of EoS to only return vars needed by XNet --- tools/starkiller-helmholtz/actual_eos.F90 | 384 +++++++++++++++++++++- 1 file changed, 378 insertions(+), 6 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 1fd0588c..d0e590ec 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -80,6 +80,7 @@ module actual_eos_module real(dp), parameter :: sioncon = (2.0_dp * pi * amu * kerg)/(h*h) real(dp), parameter :: forth = 4.0_dp/3.0_dp real(dp), parameter :: forpi = 4.0_dp * pi + real(rt), parameter :: forthpi = forth * pi real(dp), parameter :: kergavo = kerg * avo_eos real(dp), parameter :: ikavo = 1.0_dp/kergavo real(dp), parameter :: asoli3 = asol/3.0_dp @@ -110,6 +111,7 @@ module actual_eos_module !$acc create(do_coulomb, input_is_constant) public :: actual_eos, actual_eos_init, actual_eos_finalize, eos_supports_input_type + public :: xnet_actual_eos, actual_eos_eta, actual_eos_cv contains @@ -217,7 +219,7 @@ subroutine actual_eos(input, state) real(dp) :: xnew, xtol, dvdx, smallx, error, v real(dp) :: v1, v2, dv1dt, dv1dr, dv2dt,dv2dr, delr, error1, error2, told, rold, tnew, rnew, v1i, v2i - real(dp) :: x,y,zz,zzi,deni,tempi,xni,dxnidd,dxnida, & + real(dp) :: x,y,z,zz,zzi,deni,tempi,xni,dxnidd,dxnida, & dpepdt,dpepdd,deepdt,deepdd,dsepdd,dsepdt, & dpraddd,dpraddt,deraddd,deraddt,dpiondd,dpiondt, & deiondd,deiondt,dsraddd,dsraddt,dsiondd,dsiondt, & @@ -226,19 +228,18 @@ subroutine actual_eos(input, state) dpresdt,denerdd,denerdt,dentrdd,dentrdt,cv,cp, & gam1,gam2,gam3,chit,chid,nabad,sound,etaele, & detadt,detadd,xnefer,dxnedt,dxnedd,s, & - temp,den,abar,zbar,ytot1,ye + temp,den,abar,zbar,ytot1,ye,din !..for the interpolations integer :: iat,jat real(dp) :: free,df_d,df_t,df_tt,df_dt - real(dp) :: xt,xd,mxt,mxd, & + real(dp) :: xt,xd,mxt,mxd,fi(36), & si0t,si1t,si2t,si0mt,si1mt,si2mt, & si0d,si1d,si2d,si0md,si1md,si2md, & dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt, & dsi0d,dsi1d,dsi2d,dsi0md,dsi1md,dsi2md, & - ddsi0t,ddsi1t,ddsi2t,ddsi0mt,ddsi1mt,ddsi2mt, & - z,din,fi(36) + ddsi0t,ddsi1t,ddsi2t,ddsi0mt,ddsi1mt,ddsi2mt !..for the coulomb corrections real(dp) :: dsdd,dsda,lami,inv_lami,lamida,lamidd, & @@ -248,7 +249,6 @@ subroutine actual_eos(input, state) scoul,dscouldd,dscouldt,dscoulda,dscouldz real(dp) :: p_temp, e_temp - real(dp) :: smallt, smalld call eos_get_small_temp(smallt) @@ -1146,6 +1146,378 @@ subroutine actual_eos(input, state) end subroutine actual_eos + subroutine xnet_actual_eos(input, state) + + ! This is a pruned version of actual_eos that only calculates those + ! quantities needed by XNet: electron chemical potential, its derivative + ! w.r.t. temperature, and specific heat. + + !$acc routine seq + + implicit none + + !..input arguments + integer, intent(in ) :: input + type (eos_t), intent(inout) :: state + + !..declare local variables + real(dp) :: cv,etaele,detadt, & + temp,den,abar,zbar,ye + + temp = state % T + den = state % rho + abar = state % abar + zbar = state % zbar + + ye = state % y_e + + call actual_eos_cv(temp,den,abar,zbar,ye,cv) + call actual_eos_eta(temp,den,ye,etaele,detadt) + + state % cv = cv + state % eta = etaele + state % detadt = detadt + + end subroutine xnet_actual_eos + + + subroutine actual_eos_eta(temp,den,ye,etaele,detadt) + + ! This is a pruned version of actual_eos that only calculates those + ! quantities needed by XNet: electron chemical potential, its derivative + ! w.r.t. temperature, and specific heat. + + !$acc routine seq + + implicit none + + !..input arguments + real(rt), intent(in ) :: temp,den,ye + real(rt), intent(out) :: etaele,detadt + + !..declare local variables + real(rt) :: din + + !..for the interpolations + integer :: iat,jat + real(rt) :: xt,xd,mxt,mxd,dfi(16), & + si0t,si1t,si2t,si0mt,si1mt,si2mt, & + si0d,si1d,si2d,si0md,si1md,si2md, & + dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt + + din = ye * den + + !..electron-positron section: + !..hash locate this temperature and density + jat = int((log10(temp) - tlo)*tstpi) + 1 + jat = max(1,min(jat,jtmax-1)) + iat = int((log10(din) - dlo)*dstpi) + 1 + iat = max(1,min(iat,itmax-1)) + + !..various differences + xt = max( (temp - t(jat))*dti_sav(jat), 0.0_rt) + xd = max( (din - d(iat))*ddi_sav(iat), 0.0_rt) + mxt = 1.0_dp - xt + mxd = 1.0_dp - xd + + !..now get the pressure derivative with density, chemical potential, and + !..electron positron number densities + !..get the interpolation weight functions + si0t = xpsi0(xt) + si1t = xpsi1(xt)*dt_sav(jat) + + si0mt = xpsi0(mxt) + si1mt = -xpsi1(mxt)*dt_sav(jat) + + si0d = xpsi0(xd) + si1d = xpsi1(xd)*dd_sav(iat) + + si0md = xpsi0(mxd) + si1md = -xpsi1(mxd)*dd_sav(iat) + + !..derivatives of weight functions + dsi0t = xdpsi0(xt)*dti_sav(jat) + dsi1t = xdpsi1(xt) + + dsi0mt = -xdpsi0(mxt)*dti_sav(jat) + dsi1mt = xdpsi1(mxt) + + !..look in the electron chemical potential table only once + dfi(1) = ef(iat,jat) + dfi(2) = ef(iat+1,jat) + dfi(3) = ef(iat,jat+1) + dfi(4) = ef(iat+1,jat+1) + dfi(5) = eft(iat,jat) + dfi(6) = eft(iat+1,jat) + dfi(7) = eft(iat,jat+1) + dfi(8) = eft(iat+1,jat+1) + dfi(9) = efd(iat,jat) + dfi(10) = efd(iat+1,jat) + dfi(11) = efd(iat,jat+1) + dfi(12) = efd(iat+1,jat+1) + dfi(13) = efdt(iat,jat) + dfi(14) = efdt(iat+1,jat) + dfi(15) = efdt(iat,jat+1) + dfi(16) = efdt(iat+1,jat+1) + + !..electron chemical potential etaele + etaele = h3( dfi, & + si0t, si1t, si0mt, si1mt, & + si0d, si1d, si0md, si1md) + + !..derivative with respect to temperature + detadt = h3( dfi, & + dsi0t, dsi1t, dsi0mt, dsi1mt, & + si0d, si1d, si0md, si1md) + + end subroutine actual_eos_eta + + + subroutine actual_eos_cv(temp,den,abar,zbar,ye,cv) + + ! This is a pruned version of actual_eos that only calculates those + ! quantities needed by XNet: electron chemical potential, its derivative + ! w.r.t. temperature, and specific heat. + + implicit none + + !$acc routine seq + + !..input arguments + real(rt), intent(in ) :: temp,den,abar,zbar,ye + real(rt), intent(out) :: cv + + !..declare local variables + real(rt) :: x,y,z,deni,tempi,xni, & + deepdt,dsepdt, & + dpraddt,deraddt,dpiondt, & + deiondt, & + kt,ktinv,prad,erad,pion,eion, & + pele,eele,sele, & + s, & + ytot1,din + + + !..for the interpolations + integer :: iat,jat + real(rt) :: free,df_d,df_t,df_tt + real(rt) :: xt,xd,mxt,mxd,fi(36),dfi(16), & + si0t,si1t,si2t,si0mt,si1mt,si2mt, & + si0d,si1d,si2d,si0md,si1md,si2md, & + dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt, & + dsi0d,dsi1d,dsi2d,dsi0md,dsi1md,dsi2md, & + ddsi0t,ddsi1t,ddsi2t,ddsi0mt,ddsi1mt,ddsi2mt + + !..for the coulomb corrections + real(rt) :: lami,inv_lami, & + plasg,plasgdt, & + ecoul,decouldt, & + pcoul,dpcouldt + + real(rt) :: p_temp, e_temp + + !$gpu + + ytot1 = 1.0_dp/abar + din = ye * den + + !..initialize + deni = 1.0_dp/den + tempi = 1.0_dp/temp + kt = kerg * temp + ktinv = 1.0_dp/kt + + !..radiation section: + prad = asoli3 * temp * temp * temp * temp + dpraddt = 4.0_dp * prad*tempi + deraddt = 3.0_dp * dpraddt*deni + + !..ion section: + xni = avo_eos * ytot1 * den + pion = xni * kt + dpiondt = xni * kerg + deiondt = 1.5_dp * dpiondt*deni + + !..electron-positron section: + !..hash locate this temperature and density + jat = int((log10(temp) - tlo)*tstpi) + 1 + jat = max(1,min(jat,jtmax-1)) + iat = int((log10(din) - dlo)*dstpi) + 1 + iat = max(1,min(iat,itmax-1)) + + !..various differences + xt = max( (temp - t(jat))*dti_sav(jat), 0.0_rt) + xd = max( (din - d(iat))*ddi_sav(iat), 0.0_rt) + mxt = 1.0_dp - xt + mxd = 1.0_dp - xd + + !..the six density and six temperature basis functions + si0t = psi0(xt) + si1t = psi1(xt)*dt_sav(jat) + si2t = psi2(xt)*dt2_sav(jat) + + si0mt = psi0(mxt) + si1mt = -psi1(mxt)*dt_sav(jat) + si2mt = psi2(mxt)*dt2_sav(jat) + + si0d = psi0(xd) + si1d = psi1(xd)*dd_sav(iat) + si2d = psi2(xd)*dd2_sav(iat) + + si0md = psi0(mxd) + si1md = -psi1(mxd)*dd_sav(iat) + si2md = psi2(mxd)*dd2_sav(iat) + + !..derivatives of the weight functions + !dsi0t = dpsi0(xt)*dti_sav(jat) + !dsi1t = dpsi1(xt) + !dsi2t = dpsi2(xt)*dt_sav(jat) + + !dsi0mt = -dpsi0(mxt)*dti_sav(jat) + !dsi1mt = dpsi1(mxt) + !dsi2mt = -dpsi2(mxt)*dt_sav(jat) + + dsi0d = dpsi0(xd)*ddi_sav(iat) + dsi1d = dpsi1(xd) + dsi2d = dpsi2(xd)*dd_sav(iat) + + dsi0md = -dpsi0(mxd)*ddi_sav(iat) + dsi1md = dpsi1(mxd) + dsi2md = -dpsi2(mxd)*dd_sav(iat) + + !..second derivatives of the weight functions + ddsi0t = ddpsi0(xt)*dt2i_sav(jat) + ddsi1t = ddpsi1(xt)*dti_sav(jat) + ddsi2t = ddpsi2(xt) + + ddsi0mt = ddpsi0(mxt)*dt2i_sav(jat) + ddsi1mt = -ddpsi1(mxt)*dti_sav(jat) + ddsi2mt = ddpsi2(mxt) + + !..access the table locations only once + fi(1) = f(iat,jat) + fi(2) = f(iat+1,jat) + fi(3) = f(iat,jat+1) + fi(4) = f(iat+1,jat+1) + fi(5) = ft(iat,jat) + fi(6) = ft(iat+1,jat) + fi(7) = ft(iat,jat+1) + fi(8) = ft(iat+1,jat+1) + fi(9) = ftt(iat,jat) + fi(10) = ftt(iat+1,jat) + fi(11) = ftt(iat,jat+1) + fi(12) = ftt(iat+1,jat+1) + fi(13) = fd(iat,jat) + fi(14) = fd(iat+1,jat) + fi(15) = fd(iat,jat+1) + fi(16) = fd(iat+1,jat+1) + fi(17) = fdd(iat,jat) + fi(18) = fdd(iat+1,jat) + fi(19) = fdd(iat,jat+1) + fi(20) = fdd(iat+1,jat+1) + fi(21) = fdt(iat,jat) + fi(22) = fdt(iat+1,jat) + fi(23) = fdt(iat,jat+1) + fi(24) = fdt(iat+1,jat+1) + fi(25) = fddt(iat,jat) + fi(26) = fddt(iat+1,jat) + fi(27) = fddt(iat,jat+1) + fi(28) = fddt(iat+1,jat+1) + fi(29) = fdtt(iat,jat) + fi(30) = fdtt(iat+1,jat) + fi(31) = fdtt(iat,jat+1) + fi(32) = fdtt(iat+1,jat+1) + fi(33) = fddtt(iat,jat) + fi(34) = fddtt(iat+1,jat) + fi(35) = fddtt(iat,jat+1) + fi(36) = fddtt(iat+1,jat+1) + + !..the free energy + !free = h5( fi, & + ! si0t, si1t, si2t, si0mt, si1mt, si2mt, & + ! si0d, si1d, si2d, si0md, si1md, si2md) + + !..derivative with respect to density + df_d = h5( fi, & + si0t, si1t, si2t, si0mt, si1mt, si2mt, & + dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md) + + !..derivative with respect to temperature + !df_t = h5( fi, & + ! dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, & + ! si0d, si1d, si2d, si0md, si1md, si2md) + + !..derivative with respect to temperature**2 + df_tt = h5( fi, & + ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt, & + si0d, si1d, si2d, si0md, si1md, si2md) + + !..the desired electron-positron thermodynamic quantities + x = din * din + pele = x * df_d + + !sele = -df_t * ye + dsepdt = -df_tt * ye + + !eele = ye*free + temp * sele + deepdt = temp * dsepdt + + !..coulomb section: + !..uniform background corrections only + !..from yakovlev & shalybkov 1989 + !..plasg is the plasma coupling parameter + z = forth * pi + s = z * xni + + lami = 1.0_dp/s**onethird + inv_lami = 1.0_dp/lami + + plasg = zbar*zbar*esqu*ktinv*inv_lami + plasgdt = -plasg*ktinv * kerg + + ! TURN ON/OFF COULOMB + !...yakovlev & shalybkov 1989 equations 82, 85, 86, 87 + if (plasg .ge. 1.0_dp) then + x = plasg**(0.25_dp) + y = avo_eos * ytot1 * kerg + ecoul = y * temp * (a1*plasg + b1*x + c1/x + d1) + pcoul = onethird * den * ecoul + + y = avo_eos*ytot1*kt*(a1 + 0.25_dp/plasg*(b1*x - c1/x)) + decouldt = y * plasgdt + ecoul/temp + + !...yakovlev & shalybkov 1989 equations 102, 103, 104 + else if (plasg .lt. 1.0_dp) then + x = plasg*sqrt(plasg) + y = plasg**b2 + z = c2 * x - onethird * a2 * y + pcoul = -pion * z + ecoul = 3.0_dp * pcoul/den + + s = 1.5_dp*c2*x/plasg - onethird*a2*b2*y/plasg + dpcouldt = -dpiondt*z - pion*s*plasgdt + + s = 3.0_dp/den + decouldt = s * dpcouldt + end if + + ! Disable Coulomb corrections if they cause + ! the energy or pressure to go negative. + + p_temp = prad + pion + pele + pcoul + e_temp = 0.0_dp + !e_temp = erad + eion + eele + ecoul + + if (p_temp .le. ZERO ) then + !if (p_temp .le. ZERO .or. e_temp .le. ZERO) then + decouldt = 0.0_dp + end if + + !..the specific heat at constant volume (c&g 9.92) + cv = deraddt + deiondt + deepdt + decouldt + + end subroutine actual_eos_cv + subroutine actual_eos_init From b22ce4a16117f7bcac76522a18d556c9bc037ccf Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 13:28:16 -0400 Subject: [PATCH 11/13] switch XNet EOS interface to starkiller to call the XNet-specific version --- source/xnet_eos_starkiller.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/xnet_eos_starkiller.f90 b/source/xnet_eos_starkiller.f90 index 9e1f0b2a..c6b08efe 100644 --- a/source/xnet_eos_starkiller.f90 +++ b/source/xnet_eos_starkiller.f90 @@ -31,7 +31,7 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) Use xnet_controls, Only: idiag, iheat, iscrn, lun_diag Use xnet_types, Only: dp - Use actual_eos_module, Only: actual_eos + Use actual_eos_module, Only: xnet_actual_eos Use eos_type_module, Only: eos_input_rt, eos_t Implicit None @@ -63,7 +63,7 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) eos_state%zbar = ye*abar ! Call the eos - Call actual_eos(eos_input_rt,eos_state) + Call xnet_actual_eos(eos_input_rt,eos_state) ! Convert units from ergs/g to MeV/nucleon and K to GK etae = eos_state%eta From 1fb40299f5fd43868117bd0ad55146dce0051924 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 13:32:14 -0400 Subject: [PATCH 12/13] fix some oversights --- tools/starkiller-helmholtz/actual_eos.F90 | 46 +++++++++++------------ 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index d0e590ec..614a4007 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -80,7 +80,7 @@ module actual_eos_module real(dp), parameter :: sioncon = (2.0_dp * pi * amu * kerg)/(h*h) real(dp), parameter :: forth = 4.0_dp/3.0_dp real(dp), parameter :: forpi = 4.0_dp * pi - real(rt), parameter :: forthpi = forth * pi + real(dp), parameter :: forthpi = forth * pi real(dp), parameter :: kergavo = kerg * avo_eos real(dp), parameter :: ikavo = 1.0_dp/kergavo real(dp), parameter :: asoli3 = asol/3.0_dp @@ -1192,15 +1192,15 @@ subroutine actual_eos_eta(temp,den,ye,etaele,detadt) implicit none !..input arguments - real(rt), intent(in ) :: temp,den,ye - real(rt), intent(out) :: etaele,detadt + real(dp), intent(in ) :: temp,den,ye + real(dp), intent(out) :: etaele,detadt !..declare local variables - real(rt) :: din + real(dp) :: din !..for the interpolations integer :: iat,jat - real(rt) :: xt,xd,mxt,mxd,dfi(16), & + real(dp) :: xt,xd,mxt,mxd,dfi(16), & si0t,si1t,si2t,si0mt,si1mt,si2mt, & si0d,si1d,si2d,si0md,si1md,si2md, & dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt @@ -1284,11 +1284,11 @@ subroutine actual_eos_cv(temp,den,abar,zbar,ye,cv) !$acc routine seq !..input arguments - real(rt), intent(in ) :: temp,den,abar,zbar,ye - real(rt), intent(out) :: cv + real(dp), intent(in ) :: temp,den,abar,zbar,ye + real(dp), intent(out) :: cv !..declare local variables - real(rt) :: x,y,z,deni,tempi,xni, & + real(dp) :: x,y,z,deni,tempi,xni, & deepdt,dsepdt, & dpraddt,deraddt,dpiondt, & deiondt, & @@ -1300,8 +1300,8 @@ subroutine actual_eos_cv(temp,den,abar,zbar,ye,cv) !..for the interpolations integer :: iat,jat - real(rt) :: free,df_d,df_t,df_tt - real(rt) :: xt,xd,mxt,mxd,fi(36),dfi(16), & + real(dp) :: free,df_d,df_t,df_tt + real(dp) :: xt,xd,mxt,mxd,fi(36),dfi(16), & si0t,si1t,si2t,si0mt,si1mt,si2mt, & si0d,si1d,si2d,si0md,si1md,si2md, & dsi0t,dsi1t,dsi2t,dsi0mt,dsi1mt,dsi2mt, & @@ -1309,14 +1309,12 @@ subroutine actual_eos_cv(temp,den,abar,zbar,ye,cv) ddsi0t,ddsi1t,ddsi2t,ddsi0mt,ddsi1mt,ddsi2mt !..for the coulomb corrections - real(rt) :: lami,inv_lami, & + real(dp) :: lami,inv_lami, & plasg,plasgdt, & ecoul,decouldt, & pcoul,dpcouldt - real(rt) :: p_temp, e_temp - - !$gpu + real(dp) :: p_temp, e_temp ytot1 = 1.0_dp/abar din = ye * den @@ -1849,19 +1847,19 @@ function xdpsi1(z) result(xdpsi1r) end function xdpsi1 ! bicubic hermite polynomial function - function h3(fi,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) result(h3r) + function h3(dfi,w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md) result(h3r) !$acc routine seq - real(dp), intent(in) :: fi(36) + real(dp), intent(in) :: dfi(16) real(dp), intent(in) :: w0t,w1t,w0mt,w1mt,w0d,w1d,w0md,w1md real(dp) :: h3r - h3r = fi(1) *w0d*w0t + fi(2) *w0md*w0t & - + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt & - + fi(5) *w0d*w1t + fi(6) *w0md*w1t & - + fi(7) *w0d*w1mt + fi(8) *w0md*w1mt & - + fi(9) *w1d*w0t + fi(10) *w1md*w0t & - + fi(11) *w1d*w0mt + fi(12) *w1md*w0mt & - + fi(13) *w1d*w1t + fi(14) *w1md*w1t & - + fi(15) *w1d*w1mt + fi(16) *w1md*w1mt + h3r = dfi(1) *w0d*w0t + dfi(2) *w0md*w0t & + + dfi(3) *w0d*w0mt + dfi(4) *w0md*w0mt & + + dfi(5) *w0d*w1t + dfi(6) *w0md*w1t & + + dfi(7) *w0d*w1mt + dfi(8) *w0md*w1mt & + + dfi(9) *w1d*w0t + dfi(10) *w1md*w0t & + + dfi(11) *w1d*w0mt + dfi(12) *w1md*w0mt & + + dfi(13) *w1d*w1t + dfi(14) *w1md*w1t & + + dfi(15) *w1d*w1mt + dfi(16) *w1md*w1mt end function h3 subroutine actual_eos_finalize From 9175d97bc88aa2829308e9760e5c49077f1001b8 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Thu, 27 Jun 2024 13:33:00 -0400 Subject: [PATCH 13/13] typos --- tools/starkiller-helmholtz/actual_eos.F90 | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/tools/starkiller-helmholtz/actual_eos.F90 b/tools/starkiller-helmholtz/actual_eos.F90 index 614a4007..5af166e9 100644 --- a/tools/starkiller-helmholtz/actual_eos.F90 +++ b/tools/starkiller-helmholtz/actual_eos.F90 @@ -1215,8 +1215,8 @@ subroutine actual_eos_eta(temp,den,ye,etaele,detadt) iat = max(1,min(iat,itmax-1)) !..various differences - xt = max( (temp - t(jat))*dti_sav(jat), 0.0_rt) - xd = max( (din - d(iat))*ddi_sav(iat), 0.0_rt) + xt = max( (temp - t(jat))*dti_sav(jat), 0.0_dp) + xd = max( (din - d(iat))*ddi_sav(iat), 0.0_dp) mxt = 1.0_dp - xt mxd = 1.0_dp - xd @@ -1344,8 +1344,8 @@ subroutine actual_eos_cv(temp,den,abar,zbar,ye,cv) iat = max(1,min(iat,itmax-1)) !..various differences - xt = max( (temp - t(jat))*dti_sav(jat), 0.0_rt) - xd = max( (din - d(iat))*ddi_sav(iat), 0.0_rt) + xt = max( (temp - t(jat))*dti_sav(jat), 0.0_dp) + xd = max( (din - d(iat))*ddi_sav(iat), 0.0_dp) mxt = 1.0_dp - xt mxd = 1.0_dp - xd @@ -1506,8 +1506,7 @@ subroutine actual_eos_cv(temp,den,abar,zbar,ye,cv) e_temp = 0.0_dp !e_temp = erad + eion + eele + ecoul - if (p_temp .le. ZERO ) then - !if (p_temp .le. ZERO .or. e_temp .le. ZERO) then + if (p_temp .le. 0.0_dp ) then decouldt = 0.0_dp end if