From 0bbcb0b8b3efaf2dccd237aa5b11469cf4100d69 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Mon, 1 Jul 2024 12:48:24 -0400 Subject: [PATCH 1/4] scalar/vector interfaces for y_moment; and changed argument list to require {x,a,z}ext instead of optional izb --- source/model_input_ascii.f90 | 24 +++---- source/xnet_abundances.f90 | 120 ++++++++++++++++++++++++++------- source/xnet_eos_bahcall.f90 | 22 ++---- source/xnet_eos_helm.f90 | 22 ++---- source/xnet_eos_starkiller.f90 | 26 +++---- source/xnet_integrate.f90 | 5 +- source/xnet_nse.f90 | 9 ++- source/xnet_screening.f90 | 4 +- 8 files changed, 140 insertions(+), 92 deletions(-) diff --git a/source/model_input_ascii.f90 b/source/model_input_ascii.f90 index 72e81db2..5f26b33d 100644 --- a/source/model_input_ascii.f90 +++ b/source/model_input_ascii.f90 @@ -161,13 +161,17 @@ Subroutine load_initial_abundances( inab_file, abund_desc, ierr, mask_in ) ! Initialize ierr = 0 - ! Allow auxiliary nucleus - ! Set to 0 to force normalization of abundances - ! and exclude aux nuclear species - iaux = 1 Do izb = zb_lo, zb_hi yestart(izb) = 0.0 ystart(:,izb) = 0.0 + + ! Allow auxiliary nucleus + ! Set to 0 to force normalization of abundances + ! and exclude aux nuclear species + iaux(izb) = 1 + xext(izb) = 0.0 + aext(izb) = 1.0 + zext(izb) = 0.0 EndDo Do izb = zb_lo, zb_hi @@ -190,12 +194,9 @@ Subroutine load_initial_abundances( inab_file, abund_desc, ierr, mask_in ) If ( t9start(izb) <= t9nse .or. yestart(izb) <= 0.0 .or. yestart(izb) >= 1.0 ) Then Call read_inab_file(inab_file(izone),abund_desc(izb),yein,yin,xext_loc,aext_loc,zext_loc,xnorm,ierr) - If ( iaux(izb)==0 ) then ! Turn off aux nucleus - yin = yin/xnorm + If ( iaux(izb) == 0 ) then ! Turn off aux nucleus + yin = yin / xnorm Write(lun_diag,"(a)") 'Normalizing initial abundances.' - xext(izb) = 0.0 - aext(izb) = 1.0 - zext(izb) = 0.0 Else xext(izb) = xext_loc aext(izb) = aext_loc @@ -203,7 +204,7 @@ Subroutine load_initial_abundances( inab_file, abund_desc, ierr, mask_in ) Endif ! If Ye is not provided in the initial abundance file explicitly, calculate it here If ( yein <= 0.0 .or. yein >= 1.0 ) & - & Call y_moment(yin,yein,ytot,abar,zbar,z2bar,zibar,izb) + & Call y_moment(yin,yein,ytot,abar,zbar,z2bar,zibar,xext(izb),aext(izb),zext(izb)) yestart(izb) = yein @@ -223,9 +224,6 @@ Subroutine load_initial_abundances( inab_file, abund_desc, ierr, mask_in ) Call nse_solve(rhostart(izb),t9start(izb),yestart(izb)) ystart(:,izb) = ynse(:) iaux(izb) = 0 - xext(izb) = 0.d0 - aext(izb) = 1.d0 - zext(izb) = 1.d0 Else ystart(:,izb) = yin(:) EndIf diff --git a/source/xnet_abundances.f90 b/source/xnet_abundances.f90 index 142b570c..0167edbe 100644 --- a/source/xnet_abundances.f90 +++ b/source/xnet_abundances.f90 @@ -20,51 +20,125 @@ Module xnet_abundances Real(dp), Allocatable :: aext(:) ! Mass number of auxiliary/external species Real(dp), Allocatable :: zext(:) ! Charge number of auxiliary/external species + Interface y_moment + Module Procedure y_moment_scalar + Module Procedure y_moment_vector + End Interface y_moment + + Private :: y_moment_internal + Contains - Subroutine y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,izb) + Subroutine y_moment_internal(y,ye,ytot,abar,zbar,z2bar,zibar,xext,yext,zext) !----------------------------------------------------------------------------------------------- ! This routine calculates moments of the abundance distribution for the EOS. !----------------------------------------------------------------------------------------------- Use nuclear_data, Only: ny, aa, zz, zz2, zzi - Use xnet_controls, Only: idiag, lun_diag Use xnet_constants, Only: thbim1 Use xnet_types, Only: dp Implicit None ! Input variables Real(dp), Intent(in) :: y(ny) - Integer, Intent(in),Optional :: izb + Real(dp), Intent(in) :: xext, yext, zext ! Output variables Real(dp), Intent(out) :: ye, ytot, abar, zbar, z2bar, zibar ! Local variables - Real(dp) :: atot, ztot, yext, xext_loc,zext_loc, aext_loc - - If (present(izb)) Then - yext = xext(izb)/aext(izb) - xext_loc=xext(izb) - aext_loc=aext(izb) - zext_loc=zext(izb) - Else - yext=0.0 - xext_loc=0.0 - aext_loc=1.0 - zext_loc=0.0 - Endif + Real(dp) :: atot, ztot + ! Calculate abundance moments ytot = sum(y) + yext - atot = sum(y(:) * aa(:)) + xext_loc - ztot = sum(y * zz) + yext*zext_loc + atot = sum(y * aa) + xext + ztot = sum(y * zz) + yext * zext abar = atot / ytot zbar = ztot / ytot - z2bar = ( sum(y * zz2) + yext * zext_loc * zext_loc ) / ytot - zibar = ( sum(y * zzi) + yext * zext_loc**thbim1 ) / ytot - ye = ztot / atot - If ( idiag >= 3 ) Write(lun_diag,"(a4,7es23.15)") 'YMom',ytot,abar,zbar,z2bar,zibar,ye, atot + z2bar = ( sum(y * zz2) + yext * zext * zext ) / ytot + zibar = ( sum(y * zzi) + yext * zext**thbim1 ) / ytot + ye = ztot / atot + + Return + End Subroutine y_moment_internal + + Subroutine y_moment_scalar(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) + !----------------------------------------------------------------------------------------------- + ! This is the interface for the scalar version. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny + Use xnet_controls, Only: idiag, lun_diag + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: y(ny) + Real(dp), Intent(in) :: xext, aext, zext + + ! Output variables + Real(dp), Intent(out) :: ye, ytot, abar, zbar, z2bar, zibar + + ! Local variables + Real(dp) :: yext + + ! Calculate abundance moments + yext = xext / aext + Call y_moment_internal(y,ye,ytot,abar,zbar,z2bar,zibar,xext,yext,zext) + If ( idiag >= 3 ) Write(lun_diag,"(a4,6es23.15)") 'YMom',ytot,abar,zbar,z2bar,zibar,ye + + Return + End Subroutine y_moment_scalar + + Subroutine y_moment_vector(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext,mask_in) + !----------------------------------------------------------------------------------------------- + ! This is the interface for the vector version. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny + Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: y(ny,zb_lo:zb_hi) + Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi) + + ! Output variables + Real(dp), Intent(out) :: ye(zb_lo:zb_hi), ytot(zb_lo:zb_hi), abar(zb_lo:zb_hi) + Real(dp), Intent(out) :: zbar(zb_lo:zb_hi), z2bar(zb_lo:zb_hi), zibar(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Real(dp) :: yext + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + ! Calculate abundance moments + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + yext = xext(izb) / aext(izb) + Call y_moment_internal(y(:,izb),ye(izb),ytot(izb), & + & abar(izb),zbar(izb),z2bar(izb),zibar(izb),xext(izb),yext,zext(izb)) + EndIf + EndDo + + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Write(lun_diag,"(a4,6es23.15)") 'YMom', & + ytot(izb),abar(izb),zbar(izb),z2bar(izb),zibar(izb),ye(izb) + EndIf + EndDo + EndIf Return - End Subroutine y_moment + End Subroutine y_moment_vector End Module xnet_abundances diff --git a/source/xnet_eos_bahcall.f90 b/source/xnet_eos_bahcall.f90 index 70bb85c3..88a30e7a 100644 --- a/source/xnet_eos_bahcall.f90 +++ b/source/xnet_eos_bahcall.f90 @@ -18,7 +18,7 @@ Subroutine eos_initialize Return End Subroutine eos_initialize - Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,izb) + Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine updates the Equation of State for changes in temperature and density. ! @@ -41,8 +41,7 @@ Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,izb) Implicit None ! Input variables - Real(dp), Intent(in) :: t9, rho, y(ny) - Integer, Intent(in),Optional :: izb + Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext ! Output variables Real(dp), Intent(out) :: ye, cv, efermkt, defermktdt9 @@ -56,11 +55,7 @@ Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,izb) Real(dp) :: deiondt9, deraddt9, decouldt9 ! Calculate Ye and other needed moments of the abundance distribution - If (present(izb)) Then - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,izb) - Else - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar) - Endif + Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) If ( iscrn > 0 .or. iheat > 0 ) Then @@ -102,7 +97,7 @@ Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,izb) Return End Subroutine eos_interface - Subroutine eos_screen(t9,rho,y,efermkt,defermktdt9,ztilde,zinter,lambda0,gammae,dztildedt9,izb) + Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine Bahcall's approach to calculate the factors needed for screening from the input ! temperature, density and composition. @@ -119,8 +114,7 @@ Subroutine eos_screen(t9,rho,y,efermkt,defermktdt9,ztilde,zinter,lambda0,gammae, Implicit None ! Input variables - Real(dp), Intent(in) :: t9, rho, y(ny), efermkt, defermktdt9 - Integer, Intent(in),Optional :: izb + Real(dp), Intent(in) :: t9, rho, y(ny), efermkt, defermktdt9, xext, aext, zext ! Output variables Real(dp), Intent(out) :: ztilde, zinter, lambda0, gammae, dztildedt9 @@ -130,11 +124,7 @@ Subroutine eos_screen(t9,rho,y,efermkt,defermktdt9,ztilde,zinter,lambda0,gammae, Real(dp) :: sratio, ae, dsratiodefermkt ! Calculate Ye and other needed moments of the abundance distribution - If (present(izb)) Then - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,izb) - Else - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar) - Endif + Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24) Call salpeter_ratio(efermkt,sratio,dsratiodefermkt) diff --git a/source/xnet_eos_helm.f90 b/source/xnet_eos_helm.f90 index 04b0a126..ebbbf408 100644 --- a/source/xnet_eos_helm.f90 +++ b/source/xnet_eos_helm.f90 @@ -21,7 +21,7 @@ Subroutine eos_initialize Return End Subroutine eos_initialize - Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) + Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext) Use nuclear_data, Only: ny Use xnet_abundances, Only: y_moment Use xnet_constants, Only: avn, epmev @@ -30,8 +30,7 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) Implicit None ! Input variables - Real(dp), Intent(in) :: t9, rho, y(ny) - Integer, Intent(in),Optional :: izb + Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext ! Ouput variables Real(dp), Intent(out) :: ye, cv, etae, detaedt9 @@ -40,11 +39,7 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) Real(dp) :: ytot, abar, zbar, z2bar, zibar ! Calculate Ye - If (present(izb)) Then - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,izb) - Else - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar) - Endif + Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) If ( iscrn > 0 .or. iheat > 0 ) Then @@ -73,7 +68,7 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) Return End Subroutine eos_interface - Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,izb) + Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine uses the current composition and prior updates to the Equation of State to ! calculate the factors needed for screening. @@ -87,8 +82,7 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild Implicit None ! Input variables - Real(dp), Intent(in) :: t9, rho, y(ny), etae, detaedt9 - Integer, Intent(in),Optional :: izb + Real(dp), Intent(in) :: t9, rho, y(ny), etae, detaedt9, xext, aext, zext ! Output variables Real(dp), Intent(out) :: ztilde, zinter, lambda0, gammae, dztildedt9 @@ -98,11 +92,7 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild Real(dp) :: sratio, ae, dsratiodeta ! Calculate Ye and other needed moments of the abundance distribution - If (present(izb)) Then - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,izb) - Else - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar) - Endif + Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24) Call salpeter_ratio(etae,sratio,dsratiodeta) diff --git a/source/xnet_eos_starkiller.f90 b/source/xnet_eos_starkiller.f90 index af3fec4e..fe2e1d45 100644 --- a/source/xnet_eos_starkiller.f90 +++ b/source/xnet_eos_starkiller.f90 @@ -65,7 +65,7 @@ Subroutine eosx(t9,rho,ye,abar,zbar,cv,etae,detaedt9) End Subroutine eosx - Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) + Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine updates the equation of state for changes in temperature and density. !----------------------------------------------------------------------------------------------- @@ -76,21 +76,16 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) Implicit None ! Input variables - Real(dp), Intent(in) :: t9, rho, y(ny) - Integer, Intent(in), Optional :: izb + Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext + ! Ouput variables Real(dp), Intent(out) :: ye, cv, etae, detaedt9 ! Local variables Real(dp) :: ytot, abar, zbar, z2bar, zibar - Real(dp) :: xext, aext, zext, xnet, xnorm ! Calculate Ye - If (present(izb)) Then - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,izb) - Else - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar) - Endif + Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) ! Call the eos Call eosx(t9,rho,ye,abar,zbar,cv,etae,detaedt9) @@ -99,7 +94,7 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,izb) Return End Subroutine eos_interface - Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,izb) + Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine uses the current composition and prior updates to the Equation of State to ! calculate the factors needed for screening. @@ -112,8 +107,8 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild Implicit None ! Input variables - Real(dp), Intent(in) :: t9, rho, y(ny), etae, detaedt9 - Integer, Intent(in), Optional :: izb + Real(dp), Intent(in) :: t9, rho, y(ny), etae, detaedt9, xext, aext, zext + ! Output variables Real(dp), Intent(out) :: ztilde, zinter, lambda0, gammae, dztildedt9 @@ -122,12 +117,7 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild Real(dp) :: sratio, ae, dsratiodeta ! Calculate Ye and other needed moments of the abundance distribution - If (present(izb)) Then - ! Takes auxiliary species into account, if it exists - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,izb) - Else - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar) - EndIf + Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5) Call salpeter_ratio(etae,sratio,dsratiodeta) diff --git a/source/xnet_integrate.f90 b/source/xnet_integrate.f90 index 247f3460..f93bd72e 100644 --- a/source/xnet_integrate.f90 +++ b/source/xnet_integrate.f90 @@ -278,7 +278,7 @@ Subroutine update_eos(mask_in) !----------------------------------------------------------------------------------------------- ! This routine updates the dependent thermodynamic state variables by interfacing with the EoS. !----------------------------------------------------------------------------------------------- - Use xnet_abundances, Only: yt + Use xnet_abundances, Only: yt, xext, aext, zext Use xnet_conditions, Only: t9t, rhot, yet, cv, etae, detaedt9 Use xnet_controls, Only: zb_lo, zb_hi, lzactive Use xnet_eos, Only: eos_interface @@ -305,7 +305,8 @@ Subroutine update_eos(mask_in) Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - call eos_interface(t9t(izb),rhot(izb),yt(:,izb),yet(izb),cv(izb),etae(izb),detaedt9(izb),izb) + call eos_interface(t9t(izb),rhot(izb),yt(:,izb),yet(izb),cv(izb),etae(izb),detaedt9(izb), & + & xext(izb),aext(izb),zext(izb)) EndIf EndDo diff --git a/source/xnet_nse.f90 b/source/xnet_nse.f90 index aaaf6c4a..fce2e9c4 100644 --- a/source/xnet_nse.f90 +++ b/source/xnet_nse.f90 @@ -1242,6 +1242,7 @@ Subroutine nse_screen Real(dp) :: theta12iw(izmax-1), theta12is(izmax-1), theta12si(izmax-1) Real(dp) :: fhs(0:izmax+1), fhi(0:izmax+1), gammaz(0:izmax+1) Real(dp) :: cv, etae, detaedt9, ztot, ztilde, zinter, lambda0, gammae, dztildedt9, s0 + Real(dp) :: xext,aext,zext start_timer = xnet_wtime() timer_nsescrn = timer_nsescrn - start_timer @@ -1251,8 +1252,12 @@ Subroutine nse_screen ElseIf ( .not. use_CP98 ) Then ! Call EOS to get plasma quantities - call eos_interface(t9nse,rhonse,ynse,ztot,cv,etae,detaedt9) - Call eos_screen(t9nse,rhonse,ynse,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9) + xext = 0.0 + aext = 1.0 + zext = 0.0 + call eos_interface(t9nse,rhonse,ynse,ztot,cv,etae,detaedt9,xext,aext,zext) + Call eos_screen(t9nse,rhonse,ynse,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9, & + xext,aext,zext) !--------------------------------------------------------------------------------------------- ! Calculate screening energies as a function of Z, for prescriptions that follow this approach diff --git a/source/xnet_screening.f90 b/source/xnet_screening.f90 index 820b4cc2..45631817 100644 --- a/source/xnet_screening.f90 +++ b/source/xnet_screening.f90 @@ -141,7 +141,7 @@ Subroutine screening(mask_in) !----------------------------------------------------------------------------------------------- Use nuclear_data, Only: izmax, nname, zseq, zseq53, zseqi Use reaction_data, Only: n2i, n3i, n4i, nreac - Use xnet_abundances, Only: yt + Use xnet_abundances, Only: yt, xext, aext, zext Use xnet_constants, Only: bi, bip1, cds, kbi, thbim2 Use xnet_conditions, Only: rhot, t9t, etae, detaedt9 Use xnet_controls, Only: idiag, iheat, iscrn, lun_diag, szbatch, zb_lo, zb_hi, lzactive @@ -183,7 +183,7 @@ Subroutine screening(mask_in) ! Call EOS to get plasma quantities call eos_screen(t9t(izb),rhot(izb),yt(:,izb),etae(izb),detaedt9(izb), & - & ztilde,zinter,lambda0,gammae,dztildedt9,izb) + & ztilde,zinter,lambda0,gammae,dztildedt9,xext(izb),aext(izb),zext(izb)) ! Calculate screening energies as a function of Z, for prescriptions that follow this approach gammaz(0) = 0.0d0 From e5d66f4e77f8c2c9e596b66b1124b5e3cbf5d1e7 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Mon, 1 Jul 2024 13:24:46 -0400 Subject: [PATCH 2/4] scalar/vector interface for t9rhofind --- source/xnet_conditions.f90 | 72 ++++++++++++++++++++++++++++++++++---- 1 file changed, 65 insertions(+), 7 deletions(-) diff --git a/source/xnet_conditions.f90 b/source/xnet_conditions.f90 index e1b81527..91f2686d 100644 --- a/source/xnet_conditions.f90 +++ b/source/xnet_conditions.f90 @@ -38,24 +38,82 @@ Module xnet_conditions Real(dp), Allocatable :: tstart(:),tstop(:),tdelstart(:),t9start(:),rhostart(:),yestart(:) Real(dp), Allocatable :: th(:,:),t9h(:,:),rhoh(:,:),yeh(:,:) + Interface t9rhofind + Module Procedure t9rhofind_scalar + Module Procedure t9rhofind_vector + End Interface t9rhofind + Contains - Subroutine t9rhofind(kstep,tf,nf,t9f,rhof,mask_in) + Subroutine t9rhofind_scalar(kstep,tf,nf,t9f,rhof,ns,ts,t9s,rhos) + !----------------------------------------------------------------------------------------------- + ! This routine calculates t9 and rho as a function of time, either via interpolation or from an + ! analytic expression. (scalar interface) + !----------------------------------------------------------------------------------------------- + Use, Intrinsic :: iso_fortran_env, Only: lun_stdout=>output_unit + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Integer, Intent(in) :: kstep, ns + Real(dp), Intent(in) :: tf, ts(:), t9s(:), rhos(:) + + ! Output variables + Integer, Intent(out) :: nf + Real(dp), Intent(out) :: t9f, rhof + + ! Local variables + Real(dp) :: dt, rdt, dt9, drho + Integer :: n + + ! For constant conditions (ns = 1), set temperature and density + If ( ns == 1 ) Then + t9f = t9s(1) + rhof = rhos(1) + nf = 1 + + ! Otherwise, calculate T9 and rho by interpolation + Else + Do n = 1, ns + If ( tf <= ts(n) ) Exit + EndDo + nf = n + If ( n > 1 .and. n <= ns ) Then + rdt = 1.0 / (ts(n)-ts(n-1)) + dt = tf - ts(n-1) + dt9 = t9s(n) - t9s(n-1) + drho = rhos(n) - rhos(n-1) + t9f = dt*rdt*dt9 + t9s(n-1) + rhof = dt*rdt*drho + rhos(n-1) + ElseIf ( n == 1 ) Then + t9f = t9s(1) + rhof = rhos(1) + Else + t9f = t9s(ns) + rhof = rhos(ns) + Write(lun_stdout,*) 'Time beyond thermodynamic range',tf,' >',ts(ns) + EndIf + EndIf + + Return + End Subroutine t9rhofind_scalar + + Subroutine t9rhofind_vector(kstep,tf,nf,t9f,rhof,mask_in) !----------------------------------------------------------------------------------------------- ! This routine calculates t9 and rho as a function of time, either via interpolation or from an - ! analytic expression. + ! analytic expression. (vector interface) !----------------------------------------------------------------------------------------------- - Use xnet_controls, Only: lun_diag, lun_stdout, nzevolve, zb_lo, zb_hi, lzactive + Use xnet_controls, Only: lun_diag, lun_stdout, zb_lo, zb_hi, lzactive Use xnet_types, Only: dp Implicit None ! Input variables Integer, Intent(in) :: kstep - Real(dp), Intent(in) :: tf(nzevolve) + Real(dp), Intent(in) :: tf(zb_lo:zb_hi) ! Input/Output variables - Integer, Intent(inout) :: nf(nzevolve) - Real(dp), Intent(inout) :: t9f(nzevolve), rhof(nzevolve) + Integer, Intent(inout) :: nf(zb_lo:zb_hi) + Real(dp), Intent(inout) :: t9f(zb_lo:zb_hi), rhof(zb_lo:zb_hi) ! Optional variables Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) @@ -107,6 +165,6 @@ Subroutine t9rhofind(kstep,tf,nf,t9f,rhof,mask_in) EndDo Return - End Subroutine t9rhofind + End Subroutine t9rhofind_vector End Module xnet_conditions From e6631ab5572746f121d0cc13cb19fa13fbc0fd96 Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Mon, 1 Jul 2024 13:51:23 -0400 Subject: [PATCH 3/4] scalar/vector versions of eos_interface, eos_screen, and salpeter_ratio --- source/xnet_eos_starkiller.f90 | 220 +++++++++++++++++++++++++++++++-- 1 file changed, 211 insertions(+), 9 deletions(-) diff --git a/source/xnet_eos_starkiller.f90 b/source/xnet_eos_starkiller.f90 index fe2e1d45..0c223d91 100644 --- a/source/xnet_eos_starkiller.f90 +++ b/source/xnet_eos_starkiller.f90 @@ -6,7 +6,24 @@ !*************************************************************************************************** Module xnet_eos + Use xnet_types, Only: dp Implicit None + Real(dp), Allocatable :: ye(:), ytot(:), abar(:), zbar(:), z2bar(:), zibar(:), sratio(:) + + Interface eos_interface + Module Procedure eos_interface_scalar + Module Procedure eos_interface_vector + End Interface + + Interface eos_screen + Module Procedure eos_screen_scalar + Module Procedure eos_screen_vector + End Interface + + Interface salpeter_ratio + Module Procedure salpeter_ratio_scalar + Module Procedure salpeter_ratio_vector + End Interface Contains @@ -14,10 +31,28 @@ Subroutine eos_initialize !----------------------------------------------------------------------------------------------- ! This routine initializes starkiller !----------------------------------------------------------------------------------------------- + Use xnet_controls, Only: nzevolve + Use actual_eos_module, Only: actual_eos_init Implicit None + Call actual_eos_init() + Allocate (ye(nzevolve)) + Allocate (ytot(nzevolve)) + Allocate (abar(nzevolve)) + Allocate (zbar(nzevolve)) + Allocate (z2bar(nzevolve)) + Allocate (zibar(nzevolve)) + Allocate (sratio(nzevolve)) + ye = 0.0 + ytot = 0.0 + abar = 0.0 + zbar = 0.0 + z2bar = 0.0 + zibar = 0.0 + sratio = 0.0 + Return End Subroutine eos_initialize @@ -65,7 +100,7 @@ Subroutine eosx(t9,rho,ye,abar,zbar,cv,etae,detaedt9) End Subroutine eosx - Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext) + Subroutine eos_interface_scalar(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine updates the equation of state for changes in temperature and density. !----------------------------------------------------------------------------------------------- @@ -89,12 +124,68 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext) ! Call the eos Call eosx(t9,rho,ye,abar,zbar,cv,etae,detaedt9) + If ( idiag >= 3 ) Write(lun_diag,"(a,6es23.15)") 'EOS',t9,rho,ye,cv,etae,detaedt9 Return - End Subroutine eos_interface + End Subroutine eos_interface_scalar + + Subroutine eos_interface_vector(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine updates the equation of state for changes in temperature and density. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny + Use xnet_abundances, Only: y_moment + Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi) + Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi) + + ! Ouput variables + Real(dp), Intent(out) :: ye(zb_lo:zb_hi), cv(zb_lo:zb_hi) + Real(dp), Intent(out) :: etae(zb_lo:zb_hi), detaedt9(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + ! Calculate Ye + Call y_moment(y,ye,ytot(zb_lo:zb_hi), & + & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), & + & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask) + + ! Call the eos + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Call eosx(t9(izb),rho(izb),ye(izb),abar(izb),zbar(izb),cv(izb),etae(izb),detaedt9(izb)) + EndIf + EndDo + + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Write(lun_diag,"(a,6es23.15)") 'EOS',t9(izb),rho(izb),ye(izb),cv(izb),etae(izb),detaedt9(izb) + EndIf + EndDo + EndIf + + Return + End Subroutine eos_interface_vector - Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) + Subroutine eos_screen_scalar(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine uses the current composition and prior updates to the Equation of State to ! calculate the factors needed for screening. @@ -114,16 +205,16 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild ! Local variables Real(dp) :: ye, ytot, abar, zbar, z2bar, zibar - Real(dp) :: sratio, ae, dsratiodeta + Real(dp) :: sratio ! Calculate Ye and other needed moments of the abundance distribution Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5) - Call salpeter_ratio(etae,sratio,dsratiodeta) + Call salpeter_ratio(etae,sratio,dztildedt9) ztilde = sqrt(z2bar + zbar*sratio) ! DGC, Eq. 4 If ( iheat > 0 ) Then - dztildedt9 = 0.5*zbar/ztilde * dsratiodeta*detaedt9 + dztildedt9 = 0.5*zbar/ztilde * dztildedt9*detaedt9 EndIf ! Calculate plasma quantities @@ -132,9 +223,78 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild & t9,rho,ye,z2bar,zbar,sratio,ztilde,ztilde*lambda0,gammae Return - End Subroutine eos_screen + End Subroutine eos_screen_scalar - Subroutine salpeter_ratio(eta,ratio,dratiodeta) + Subroutine eos_screen_vector(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine uses the current composition and prior updates to the Equation of State to + ! calculate the factors needed for screening. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny + Use xnet_abundances, Only: y_moment + Use xnet_controls, Only: idiag, iheat, lun_diag, zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Use xnet_util, Only: plasma + Implicit None + + ! Input variables + Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi) + Real(dp), Intent(in) :: etae(zb_lo:zb_hi), detaedt9(zb_lo:zb_hi) + Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi) + + ! Output variables + Real(dp), Intent(out) :: ztilde(zb_lo:zb_hi), zinter(zb_lo:zb_hi) + Real(dp), Intent(out) :: lambda0(zb_lo:zb_hi), gammae(zb_lo:zb_hi) + Real(dp), Intent(out) :: dztildedt9(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + ! Calculate Ye + Call y_moment(y,ye(zb_lo:zb_hi),ytot(zb_lo:zb_hi), & + & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), & + & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask) + + ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5) + Call salpeter_ratio(etae,sratio(zb_lo:zb_hi),dztildedt9,mask_in = mask) + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + ztilde(izb) = sqrt(z2bar(izb) + sratio(izb)*zbar(izb)) ! DGC, Eq. 4 + If ( iheat > 0 ) Then + dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*detaedt9(izb) + EndIf + + ! Calculate plasma quantities + Call plasma(t9(izb),rho(izb),ytot(izb),ye(izb),zbar(izb), & + & zibar(izb),ztilde(izb),zinter(izb),lambda0(izb),gammae(izb)) + EndIf + EndDo + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', & + & t9(izb),rho(izb),ye(izb),z2bar(izb),zbar(izb),sratio(izb), & + & ztilde(izb),ztilde(izb)*lambda0(izb),gammae(izb) + EndIf + EndDo + EndIf + + Return + End Subroutine eos_screen_vector + + Subroutine salpeter_ratio_scalar(eta,ratio,dratiodeta) !----------------------------------------------------------------------------------------------- ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening. ! eta is the ratio of electron chemical potential to kT. @@ -173,6 +333,48 @@ Subroutine salpeter_ratio(eta,ratio,dratiodeta) EndIf Return - End Subroutine salpeter_ratio + End Subroutine salpeter_ratio_scalar + + Subroutine salpeter_ratio_vector(eta,ratio,dratiodeta,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening. + ! eta is the ratio of electron chemical potential to kT. + ! f'/f is also defined as theta_e in DeWitt+ (1973). + ! + ! Calculation uses Fermi function relation d/dx f_(k+1) = (k+1) f_k and the rational function + ! expansions of Fukushima (2015; AMC 259 708) for the F-D integrals of order 1/2, -1/2, and -3/2. + !----------------------------------------------------------------------------------------------- + Use xnet_controls, Only: zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: eta(zb_lo:zb_hi) + + ! Output variables + Real(dp), Intent(out) :: ratio(zb_lo:zb_hi), dratiodeta(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Call salpeter_ratio_scalar(eta(izb),ratio(izb),dratiodeta(izb)) + EndIf + EndDo + + Return + End Subroutine salpeter_ratio_vector End Module xnet_eos From e07142cccd9bc43a37c41861d418c58fb38506ea Mon Sep 17 00:00:00 2001 From: Austin Harris Date: Mon, 1 Jul 2024 14:57:23 -0400 Subject: [PATCH 4/4] scalar/vector interfaces for other EoS interface modules --- source/xnet_eos_bahcall.f90 | 319 ++++++++++++++++++++++++++++----- source/xnet_eos_helm.f90 | 297 +++++++++++++++++++++++++----- source/xnet_eos_starkiller.f90 | 10 +- 3 files changed, 528 insertions(+), 98 deletions(-) diff --git a/source/xnet_eos_bahcall.f90 b/source/xnet_eos_bahcall.f90 index 88a30e7a..b6331f96 100644 --- a/source/xnet_eos_bahcall.f90 +++ b/source/xnet_eos_bahcall.f90 @@ -6,7 +6,24 @@ !*************************************************************************************************** Module xnet_eos + Use xnet_types, Only: dp Implicit None + Real(dp), Allocatable :: ye(:), ytot(:), abar(:), zbar(:), z2bar(:), zibar(:), sratio(:) + + Interface eos_interface + Module Procedure eos_interface_scalar + Module Procedure eos_interface_vector + End Interface + + Interface eos_screen + Module Procedure eos_screen_scalar + Module Procedure eos_screen_vector + End Interface eos_screen + + Interface salpeter_ratio + Module Procedure salpeter_ratio_scalar + Module Procedure salpeter_ratio_vector + End Interface Contains @@ -14,49 +31,53 @@ Subroutine eos_initialize !----------------------------------------------------------------------------------------------- ! This routine initializes the EoS. !----------------------------------------------------------------------------------------------- + Use xnet_controls, Only: nzevolve Implicit None + + Allocate (ye(nzevolve)) + Allocate (ytot(nzevolve)) + Allocate (abar(nzevolve)) + Allocate (zbar(nzevolve)) + Allocate (z2bar(nzevolve)) + Allocate (zibar(nzevolve)) + Allocate (sratio(nzevolve)) + ye = 0.0 + ytot = 0.0 + abar = 0.0 + zbar = 0.0 + z2bar = 0.0 + zibar = 0.0 + sratio = 0.0 + Return End Subroutine eos_initialize - Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext) + Subroutine eosx(t9,rho,ye,ytot,abar,zbar,z53bar,cv,efermkt,defermktdt9) !----------------------------------------------------------------------------------------------- - ! This routine updates the Equation of State for changes in temperature and density. - ! - ! Bahcall introduced a fit for the plasma parameters needed to calculate the screening factors - ! for nuclear reaction based on the easier to calculate efermkt (the ratio of the electron - ! Fermi Energy to the themal energy) instead of the more general expresions based on etae (the - ! ratio of the electron chemical potential to the thermal energy), which generally require a - ! complex Equation of State to calculate. For basic cases, XNet uses Bahcall's approach - ! to simplify the - ! - ! Note: the factors efermkt and defermktdt9 (the derivative of efermkt wrt to temperature) in the - ! subroutine interface are stored in XNet's conditions module as etae and detaedt9. For more - ! general EoS, etae and detaedt9 are their proper selves. + ! This routine interfaces with and calls the underlying EoS. !----------------------------------------------------------------------------------------------- - Use nuclear_data, Only: ny, zz53 - Use xnet_abundances, Only: y_moment Use xnet_constants, Only: asig, avn, bok, e2, pi, pi2, third, emass, clt, ele_en, hbar - Use xnet_controls, Only: idiag, iheat, iscrn, lun_diag + Use xnet_controls, Only: iheat, iscrn Use xnet_types, Only: dp Implicit None ! Input variables - Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext + Real(dp), Intent(in) :: t9, rho, ye, ytot, abar, zbar, z53bar - ! Output variables - Real(dp), Intent(out) :: ye, cv, efermkt, defermktdt9 + ! Ouput variables + Real(dp), Intent(out) :: cv, efermkt, defermktdt9 ! Local variables Real(dp), Parameter :: a1 = 0.898004, b1 = 0.96786, c1 = 0.220703, d1 = -0.86097 Real(dp), Parameter :: a2 = -0.86602540378, b2 = 0.29561, c2 = 1.9885 - Real(dp) :: ytot, bkt, abar, zbar, z2bar, zibar, z53bar + Real(dp) :: bkt Real(dp) :: etae_mb, ae, gam, gam14, gam32, gamc2, rel_ef Real(dp) :: eion, erad, ecoul Real(dp) :: deiondt9, deraddt9, decouldt9 - ! Calculate Ye and other needed moments of the abundance distribution - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) - + cv = 0.0 + efermkt = 0.0 + defermktdt9 = 0.0 If ( iscrn > 0 .or. iheat > 0 ) Then bkt = bok*t9 @@ -73,7 +94,6 @@ Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext) deraddt9 = 4.0*erad / t9 ae = (3.0 / (4.0*pi*avn*rho*ye))**third ! electron-sphere radius - z53bar = sum(zz53 * y) / ytot gam = z53bar*e2 / (ae*bkt) ! ionic Coulomb coupling parameter If ( gam >= 1.0 ) Then gam14 = gam**0.25 @@ -86,18 +106,122 @@ Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext) decouldt9 = ecoul/t9 - bok * ytot * (1.5*a2*gam32 + b2*c2*gamc2) EndIf cv = deiondt9 + deraddt9 + decouldt9 - Else - efermkt = 0.0 - defermktdt9 = 0.0 - cv = 0.0 EndIf + End Subroutine eosx + + Subroutine eos_interface_scalar(rho,t9,y,ye,cv,efermkt,defermktdt9,xext,aext,zext) + !----------------------------------------------------------------------------------------------- + ! This routine updates the Equation of State for changes in temperature and density. + ! + ! Bahcall introduced a fit for the plasma parameters needed to calculate the screening factors + ! for nuclear reaction based on the easier to calculate efermkt (the ratio of the electron + ! Fermi Energy to the themal energy) instead of the more general expresions based on etae (the + ! ratio of the electron chemical potential to the thermal energy), which generally require a + ! complex Equation of State to calculate. For basic cases, XNet uses Bahcall's approach + ! to simplify the + ! + ! Note: the factors efermkt and defermktdt9 (the derivative of efermkt wrt to temperature) in the + ! subroutine interface are stored in XNet's conditions module as etae and detaedt9. For more + ! general EoS, etae and detaedt9 are their proper selves. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny, zz53 + Use xnet_abundances, Only: y_moment + Use xnet_controls, Only: idiag, lun_diag + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext + + ! Output variables + Real(dp), Intent(out) :: ye, cv, efermkt, defermktdt9 + + ! Local variables + Real(dp) :: ytot, abar, zbar, z2bar, zibar, z53bar + + ! Calculate Ye + Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) + z53bar = sum(zz53 * y) / ytot + + ! Call the eos + Call eosx(t9,rho,ye,ytot,abar,zbar,z53bar,cv,efermkt,defermktdt9) + If ( idiag >= 3 ) Write(lun_diag,"(a,6es23.15)") 'EOS',t9,rho,ye,cv,efermkt,defermktdt9 Return - End Subroutine eos_interface + End Subroutine eos_interface_scalar - Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) + Subroutine eos_interface_vector(t9,rho,y,ye,cv,efermkt,defermktdt9,xext,aext,zext,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine updates the Equation of State for changes in temperature and density. + ! + ! Bahcall introduced a fit for the plasma parameters needed to calculate the screening factors + ! for nuclear reaction based on the easier to calculate efermkt (the ratio of the electron + ! Fermi Energy to the themal energy) instead of the more general expresions based on etae (the + ! ratio of the electron chemical potential to the thermal energy), which generally require a + ! complex Equation of State to calculate. For basic cases, XNet uses Bahcall's approach + ! to simplify the + ! + ! Note: the factors efermkt and defermktdt9 (the derivative of efermkt wrt to temperature) in the + ! subroutine interface are stored in XNet's conditions module as etae and detaedt9. For more + ! general EoS, etae and detaedt9 are their proper selves. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny, zz53 + Use xnet_abundances, Only: y_moment + Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi) + Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi) + + ! Ouput variables + Real(dp), Intent(out) :: ye(zb_lo:zb_hi), cv(zb_lo:zb_hi) + Real(dp), Intent(out) :: efermkt(zb_lo:zb_hi), defermktdt9(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Real(dp) :: z53bar + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + ! Calculate Ye + Call y_moment(y,ye,ytot(zb_lo:zb_hi), & + & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), & + & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask) + + ! Call the eos + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + z53bar = sum(zz53 * y(:,izb)) / ytot(izb) + Call eosx(t9(izb),rho(izb),ye(izb),ytot(izb),abar(izb),zbar(izb), & + & z53bar,cv(izb),efermkt(izb),defermktdt9(izb)) + EndIf + EndDo + + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Write(lun_diag,"(a,6es23.15)") 'EOS',t9(izb),rho(izb),ye(izb),cv(izb),efermkt(izb),defermktdt9(izb) + EndIf + EndDo + EndIf + + Return + End Subroutine eos_interface_vector + + Subroutine eos_screen_scalar(t9,rho,y,efermkt,defermktdt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine Bahcall's approach to calculate the factors needed for screening from the input ! temperature, density and composition. @@ -107,10 +231,9 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild !----------------------------------------------------------------------------------------------- Use nuclear_data, Only: ny Use xnet_abundances, Only: y_moment - Use xnet_constants, Only: avn, bok, clt, e2, ele_en, emass, hbar, pi, pi2, third, two3rd, & - & thbim2, twm2bi - Use xnet_controls, Only: idiag, iheat, lun_diag + Use xnet_controls, Only: idiag, lun_diag Use xnet_types, Only: dp + Use xnet_util, Only: plasma Implicit None ! Input variables @@ -127,27 +250,89 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24) - Call salpeter_ratio(efermkt,sratio,dsratiodefermkt) + Call salpeter_ratio(efermkt,sratio,dztildedt9) ztilde = sqrt(z2bar + zbar*sratio) - If ( iheat > 0 ) Then - dztildedt9 = 0.5*zbar/ztilde * dsratiodefermkt*defermktdt9 - Else - dztildedt9 = 0.0 - EndIf + dztildedt9 = 0.5*zbar/ztilde * dztildedt9*defermktdt9 ! Calculate plasma quantities - bkt = bok*t9 - lambda0 = sqrt(4.0*pi*rho*avn*ytot) * (e2/bkt)**1.5 ! DGC, Eq. 3 - ae = (3.0 / (4.0*pi*avn*rho*ye))**third ! electron-sphere radius - gammae = e2 / (ae*bkt) ! electron Coulomb coupling parameter - zinter = zibar / (ztilde**thbim2 * zbar**twm2bi) - If ( idiag >= 2 ) Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', & + Call plasma(t9,rho,ytot,ye,zbar,zibar,ztilde,zinter,lambda0,gammae) + If ( idiag >= 3 ) Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', & & t9,rho,ye,z2bar,zbar,sratio,ztilde,ztilde*lambda0,gammae Return - End Subroutine eos_screen + End Subroutine eos_screen_scalar + + Subroutine eos_screen_vector(t9,rho,y,efermkt,defermktdt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine Bahcall's approach to calculate the factors needed for screening from the input + ! temperature, density and composition. + ! + ! Note: the efermkt and defermktdt9 in the subroutine interface are stored in XNet's conditions + ! module as etae and detaedt9. For more general EoS, etae and detaedt9 are their proper selves. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny + Use xnet_abundances, Only: y_moment + Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Use xnet_util, Only: plasma + Implicit None + + ! Input variables + Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi) + Real(dp), Intent(in) :: efermkt(zb_lo:zb_hi), defermktdt9(zb_lo:zb_hi) + Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi) + + ! Output variables + Real(dp), Intent(out) :: ztilde(zb_lo:zb_hi), zinter(zb_lo:zb_hi) + Real(dp), Intent(out) :: lambda0(zb_lo:zb_hi), gammae(zb_lo:zb_hi) + Real(dp), Intent(out) :: dztildedt9(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + ! Calculate Ye + Call y_moment(y,ye(zb_lo:zb_hi),ytot(zb_lo:zb_hi), & + & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), & + & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask) + + ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5) + Call salpeter_ratio(efermkt,sratio(zb_lo:zb_hi),dztildedt9,mask_in = mask) + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + ztilde(izb) = sqrt(z2bar(izb) + sratio(izb)*zbar(izb)) ! DGC, Eq. 4 + dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*defermktdt9(izb) + + ! Calculate plasma quantities + Call plasma(t9(izb),rho(izb),ytot(izb),ye(izb),zbar(izb), & + & zibar(izb),ztilde(izb),zinter(izb),lambda0(izb),gammae(izb)) + EndIf + EndDo + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', & + & t9(izb),rho(izb),ye(izb),z2bar(izb),zbar(izb),sratio(izb), & + & ztilde(izb),ztilde(izb)*lambda0(izb),gammae(izb) + EndIf + EndDo + EndIf + + Return + End Subroutine eos_screen_vector - Subroutine salpeter_ratio(efmkt,ratio,dratiodefmkt) + Subroutine salpeter_ratio_scalar(efmkt,ratio,dratiodefmkt) !----------------------------------------------------------------------------------------------- ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening, using ! a fit to Figure 24 in that paper. efmkt is the ratio of electron chemical potential to kT. @@ -186,6 +371,44 @@ Subroutine salpeter_ratio(efmkt,ratio,dratiodefmkt) EndIf Return - End Subroutine salpeter_ratio + End Subroutine salpeter_ratio_scalar + + Subroutine salpeter_ratio_vector(efmkt,ratio,dratiodefmkt,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening, using + ! a fit to Figure 24 in that paper. efmkt is the ratio of electron chemical potential to kT. + !----------------------------------------------------------------------------------------------- + Use xnet_controls, Only: zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: efmkt(zb_lo:zb_hi) + + ! Output variables + Real(dp), Intent(out) :: ratio(zb_lo:zb_hi), dratiodefmkt(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Call salpeter_ratio_scalar(efmkt(izb),ratio(izb),dratiodefmkt(izb)) + EndIf + EndDo + + Return + End Subroutine salpeter_ratio_vector End Module xnet_eos diff --git a/source/xnet_eos_helm.f90 b/source/xnet_eos_helm.f90 index ebbbf408..aff85f84 100644 --- a/source/xnet_eos_helm.f90 +++ b/source/xnet_eos_helm.f90 @@ -6,8 +6,25 @@ !*************************************************************************************************** Module xnet_eos + Use xnet_types, Only: dp Implicit None include 'vector_eos.dek' + Real(dp), Allocatable :: ye(:), ytot(:), abar(:), zbar(:), z2bar(:), zibar(:), sratio(:) + + Interface eos_interface + Module Procedure eos_interface_scalar + Module Procedure eos_interface_vector + End Interface + + Interface eos_screen + Module Procedure eos_screen_scalar + Module Procedure eos_screen_vector + End Interface + + Interface salpeter_ratio + Module Procedure salpeter_ratio_scalar + Module Procedure salpeter_ratio_vector + End Interface Contains @@ -15,32 +32,47 @@ Subroutine eos_initialize !----------------------------------------------------------------------------------------------- ! This routine updates the Equation of State for changes in temperature and density. !----------------------------------------------------------------------------------------------- + Use xnet_controls, Only: nzevolve Implicit None + Call read_helm_table + Allocate (ye(nzevolve)) + Allocate (ytot(nzevolve)) + Allocate (abar(nzevolve)) + Allocate (zbar(nzevolve)) + Allocate (z2bar(nzevolve)) + Allocate (zibar(nzevolve)) + Allocate (sratio(nzevolve)) + ye = 0.0 + ytot = 0.0 + abar = 0.0 + zbar = 0.0 + z2bar = 0.0 + zibar = 0.0 + sratio = 0.0 + Return End Subroutine eos_initialize - Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext) - Use nuclear_data, Only: ny - Use xnet_abundances, Only: y_moment - Use xnet_constants, Only: avn, epmev - Use xnet_controls, Only: idiag, iheat, iscrn, lun_diag + Subroutine eosx(t9,rho,ye,abar,zbar,cv,etae,detaedt9) + !----------------------------------------------------------------------------------------------- + ! This routine interfaces with and calls the underlying EoS. + !----------------------------------------------------------------------------------------------- + Use xnet_constants, Only: amu + Use xnet_controls, Only: iheat, iscrn Use xnet_types, Only: dp Implicit None ! Input variables - Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext + Real(dp), Intent(in) :: t9, rho, ye, abar, zbar ! Ouput variables - Real(dp), Intent(out) :: ye, cv, etae, detaedt9 - - ! Local variables - Real(dp) :: ytot, abar, zbar, z2bar, zibar - - ! Calculate Ye - Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) + Real(dp), Intent(out) :: cv, etae, detaedt9 + cv = 0.0 + etae = 0.0 + detaedt9 = 0.0 If ( iscrn > 0 .or. iheat > 0 ) Then ! Load input variables for the eos @@ -55,30 +87,108 @@ Subroutine eos_interface(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext) Call helmeos ! Convert units from ergs/g to MeV/nucleon and K to GK + cv = cv_row(1) * amu * 1e9 etae = etaele_row(1) - detaedt9 = detat_row(1)*1e9 - cv = cv_row(1)*1.0d9/epmev/avn - Else - etae = 0.0 - detaedt9 = 0.0 - cv = 0.0 + detaedt9 = detat_row(1) * 1e9 EndIf + + End Subroutine eosx + + Subroutine eos_interface_scalar(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext) + !----------------------------------------------------------------------------------------------- + ! This routine updates the equation of state for changes in temperature and density. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny + Use xnet_abundances, Only: y_moment + Use xnet_controls, Only: idiag, lun_diag + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: t9, rho, y(ny), xext, aext, zext + + ! Ouput variables + Real(dp), Intent(out) :: ye, cv, etae, detaedt9 + + ! Local variables + Real(dp) :: ytot, abar, zbar, z2bar, zibar + + ! Calculate Ye + Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) + + ! Call the eos + Call eosx(t9,rho,ye,abar,zbar,cv,etae,detaedt9) + If ( idiag >= 3 ) Write(lun_diag,"(a,6es23.15)") 'EOS',t9,rho,ye,cv,etae,detaedt9 Return - End Subroutine eos_interface + End Subroutine eos_interface_scalar + + Subroutine eos_interface_vector(t9,rho,y,ye,cv,etae,detaedt9,xext,aext,zext,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine updates the equation of state for changes in temperature and density. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny + Use xnet_abundances, Only: y_moment + Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi) + Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi) + + ! Ouput variables + Real(dp), Intent(out) :: ye(zb_lo:zb_hi), cv(zb_lo:zb_hi) + Real(dp), Intent(out) :: etae(zb_lo:zb_hi), detaedt9(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + ! Calculate Ye + Call y_moment(y,ye,ytot(zb_lo:zb_hi), & + & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), & + & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask) + + ! Call the eos + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Call eosx(t9(izb),rho(izb),ye(izb),abar(izb),zbar(izb),cv(izb),etae(izb),detaedt9(izb)) + EndIf + EndDo + + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Write(lun_diag,"(a,6es23.15)") 'EOS',t9(izb),rho(izb),ye(izb),cv(izb),etae(izb),detaedt9(izb) + EndIf + EndDo + EndIf + + Return + End Subroutine eos_interface_vector - Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) + Subroutine eos_screen_scalar(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext) !----------------------------------------------------------------------------------------------- ! This routine uses the current composition and prior updates to the Equation of State to ! calculate the factors needed for screening. !----------------------------------------------------------------------------------------------- Use nuclear_data, Only: ny Use xnet_abundances, Only: y_moment - Use xnet_constants, Only: avn, bok, clt, e2, ele_en, emass, hbar, pi, pi2, third, two3rd, & - & thbim2, twm2bi Use xnet_controls, Only: idiag, iheat, lun_diag Use xnet_types, Only: dp + Use xnet_util, Only: plasma Implicit None ! Input variables @@ -88,37 +198,97 @@ Subroutine eos_screen(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztild Real(dp), Intent(out) :: ztilde, zinter, lambda0, gammae, dztildedt9 ! Local variables - Real(dp) :: ye, ytot, bkt, abar, zbar, z2bar, zibar - Real(dp) :: sratio, ae, dsratiodeta + Real(dp) :: ye, ytot, abar, zbar, z2bar, zibar + Real(dp) :: sratio ! Calculate Ye and other needed moments of the abundance distribution Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar,xext,aext,zext) - ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24) - Call salpeter_ratio(etae,sratio,dsratiodeta) - ztilde = sqrt(z2bar + zbar*sratio) - If ( iheat > 0 ) Then - dztildedt9 = 0.5*zbar/ztilde * dsratiodeta*detaedt9 - Else - dztildedt9 = 0.0 - EndIf + ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5) + Call salpeter_ratio(etae,sratio,dztildedt9) + ztilde = sqrt(z2bar + zbar*sratio) ! DGC, Eq. 4 + dztildedt9 = 0.5*zbar/ztilde * dztildedt9*detaedt9 ! Calculate plasma quantities - bkt = bok*t9 - lambda0 = sqrt(4.0*pi*rho*avn*ytot) * (e2/bkt)**1.5 ! DGC, Eq. 3 - ae = (3.0 / (4.0*pi*avn*rho*ye))**third ! electron-sphere radius - gammae = e2 / (ae*bkt) ! electron Coulomb coupling parameter - zinter = zibar / (ztilde**thbim2 * zbar**twm2bi) + Call plasma(t9,rho,ytot,ye,zbar,zibar,ztilde,zinter,lambda0,gammae) If ( idiag >= 3 ) Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', & & t9,rho,ye,z2bar,zbar,sratio,ztilde,ztilde*lambda0,gammae Return - End Subroutine eos_screen + End Subroutine eos_screen_scalar + + Subroutine eos_screen_vector(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae,dztildedt9,xext,aext,zext,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine uses the current composition and prior updates to the Equation of State to + ! calculate the factors needed for screening. + !----------------------------------------------------------------------------------------------- + Use nuclear_data, Only: ny + Use xnet_abundances, Only: y_moment + Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Use xnet_util, Only: plasma + Implicit None + + ! Input variables + Real(dp), Intent(in) :: t9(zb_lo:zb_hi), rho(zb_lo:zb_hi), y(ny,zb_lo:zb_hi) + Real(dp), Intent(in) :: etae(zb_lo:zb_hi), detaedt9(zb_lo:zb_hi) + Real(dp), Intent(in) :: xext(zb_lo:zb_hi), aext(zb_lo:zb_hi), zext(zb_lo:zb_hi) + + ! Output variables + Real(dp), Intent(out) :: ztilde(zb_lo:zb_hi), zinter(zb_lo:zb_hi) + Real(dp), Intent(out) :: lambda0(zb_lo:zb_hi), gammae(zb_lo:zb_hi) + Real(dp), Intent(out) :: dztildedt9(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + ! Calculate Ye + Call y_moment(y,ye(zb_lo:zb_hi),ytot(zb_lo:zb_hi), & + & abar(zb_lo:zb_hi),zbar(zb_lo:zb_hi),z2bar(zb_lo:zb_hi),zibar(zb_lo:zb_hi), & + & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask) + + ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5) + Call salpeter_ratio(etae,sratio(zb_lo:zb_hi),dztildedt9,mask_in = mask) + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + ztilde(izb) = sqrt(z2bar(izb) + sratio(izb)*zbar(izb)) ! DGC, Eq. 4 + dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*detaedt9(izb) + + ! Calculate plasma quantities + Call plasma(t9(izb),rho(izb),ytot(izb),ye(izb),zbar(izb), & + & zibar(izb),ztilde(izb),zinter(izb),lambda0(izb),gammae(izb)) + EndIf + EndDo + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', & + & t9(izb),rho(izb),ye(izb),z2bar(izb),zbar(izb),sratio(izb), & + & ztilde(izb),ztilde(izb)*lambda0(izb),gammae(izb) + EndIf + EndDo + EndIf + + Return + End Subroutine eos_screen_vector - Subroutine salpeter_ratio(eta,ratio,dratiodeta) + Subroutine salpeter_ratio_scalar(eta,ratio,dratiodeta) !----------------------------------------------------------------------------------------------- ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening. ! eta is the ratio of electron chemical potential to kT. + ! f'/f is also defined as theta_e in DeWitt+ (1973). ! ! Calculation uses Fermi function relation d/dx f_(k+1) = (k+1) f_k and the rational function ! expansions of Fukushima (2015; AMC 259 708) for the F-D integrals of order 1/2, -1/2, and -3/2. @@ -142,7 +312,7 @@ Subroutine salpeter_ratio(eta,ratio,dratiodeta) fermim = fdm1h(eta) fermip = fd1h(eta) - ! Evalutate the Salpeter ratio + ! Evalutate the Salpeter ratio (extra factor of 1/2 from FD integral definition) ratio = 0.5 * fermim/fermip If ( iheat > 0 ) Then dfmdeta = -0.5 * fdm3h(eta) @@ -153,7 +323,48 @@ Subroutine salpeter_ratio(eta,ratio,dratiodeta) EndIf Return - End Subroutine salpeter_ratio + End Subroutine salpeter_ratio_scalar -End Module xnet_eos + Subroutine salpeter_ratio_vector(eta,ratio,dratiodeta,mask_in) + !----------------------------------------------------------------------------------------------- + ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening. + ! eta is the ratio of electron chemical potential to kT. + ! f'/f is also defined as theta_e in DeWitt+ (1973). + ! + ! Calculation uses Fermi function relation d/dx f_(k+1) = (k+1) f_k and the rational function + ! expansions of Fukushima (2015; AMC 259 708) for the F-D integrals of order 1/2, -1/2, and -3/2. + !----------------------------------------------------------------------------------------------- + Use xnet_controls, Only: zb_lo, zb_hi, lzactive + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: eta(zb_lo:zb_hi) + + ! Output variables + Real(dp), Intent(out) :: ratio(zb_lo:zb_hi), dratiodeta(zb_lo:zb_hi) + + ! Optional variables + Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) + + ! Local variables + Integer :: izb + Logical, Pointer :: mask(:) + + If ( present(mask_in) ) Then + mask(zb_lo:) => mask_in + Else + mask(zb_lo:) => lzactive(zb_lo:zb_hi) + EndIf + If ( .not. any(mask) ) Return + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Call salpeter_ratio_scalar(eta(izb),ratio(izb),dratiodeta(izb)) + EndIf + EndDo + + Return + End Subroutine salpeter_ratio_vector +End Module xnet_eos diff --git a/source/xnet_eos_starkiller.f90 b/source/xnet_eos_starkiller.f90 index 0c223d91..321825b9 100644 --- a/source/xnet_eos_starkiller.f90 +++ b/source/xnet_eos_starkiller.f90 @@ -213,9 +213,7 @@ Subroutine eos_screen_scalar(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24; DGC, Eq. 5) Call salpeter_ratio(etae,sratio,dztildedt9) ztilde = sqrt(z2bar + zbar*sratio) ! DGC, Eq. 4 - If ( iheat > 0 ) Then - dztildedt9 = 0.5*zbar/ztilde * dztildedt9*detaedt9 - EndIf + dztildedt9 = 0.5*zbar/ztilde * dztildedt9*detaedt9 ! Calculate plasma quantities Call plasma(t9,rho,ytot,ye,zbar,zibar,ztilde,zinter,lambda0,gammae) @@ -232,7 +230,7 @@ Subroutine eos_screen_vector(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae !----------------------------------------------------------------------------------------------- Use nuclear_data, Only: ny Use xnet_abundances, Only: y_moment - Use xnet_controls, Only: idiag, iheat, lun_diag, zb_lo, zb_hi, lzactive + Use xnet_controls, Only: idiag, lun_diag, zb_lo, zb_hi, lzactive Use xnet_types, Only: dp Use xnet_util, Only: plasma Implicit None @@ -272,9 +270,7 @@ Subroutine eos_screen_vector(t9,rho,y,etae,detaedt9,ztilde,zinter,lambda0,gammae Do izb = zb_lo, zb_hi If ( mask(izb) ) Then ztilde(izb) = sqrt(z2bar(izb) + sratio(izb)*zbar(izb)) ! DGC, Eq. 4 - If ( iheat > 0 ) Then - dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*detaedt9(izb) - EndIf + dztildedt9(izb) = 0.5*zbar(izb)/ztilde(izb) * dztildedt9(izb)*detaedt9(izb) ! Calculate plasma quantities Call plasma(t9(izb),rho(izb),ytot(izb),ye(izb),zbar(izb), &