From ede8035802c6794b835af5c12982b4d3359ca56a Mon Sep 17 00:00:00 2001 From: cott Date: Sun, 21 Nov 2010 04:02:46 +0000 Subject: * add a new routine that computes dpde|rho and dpdrho|e at once. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/EOS_Omni@23 8e189c6b-2ab8-4400-aa02-70a9cfce18b9 --- interface.ccl | 16 ++++++ src/EOS_Omni_MultiVarCalls.F90 | 120 +++++++++++++++++++++++++++++++++++++++++ src/nuc_eos/findtemp.F90 | 6 ++- src/nuc_eos/nuc_eos.F90 | 104 +++++++++++++++++++++++++++++++++++ 4 files changed, 244 insertions(+), 2 deletions(-) diff --git a/interface.ccl b/interface.ccl index 90e6151..1903bcd 100644 --- a/interface.ccl +++ b/interface.ccl @@ -52,6 +52,22 @@ void FUNCTION EOS_Omni_DPressByDRho(CCTK_INT IN eoskey, \ PROVIDES FUNCTION EOS_Omni_DPressByDRho WITH EOS_Omni_EOS_DPressByDRho LANGUAGE Fortran +void FUNCTION EOS_Omni_dpderho_dpdrhoe(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_REAL IN rf_precision, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY dpderho, \ + CCTK_REAL OUT ARRAY dpdrhoe, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +PROVIDES FUNCTION EOS_Omni_dpderho_dpdrhoe WITH EOS_Omni_EOS_dpderho_dpdrhoe LANGUAGE Fortran + + void FUNCTION EOS_Omni_cs2(CCTK_INT IN eoskey, \ CCTK_INT IN havetemp, \ diff --git a/src/EOS_Omni_MultiVarCalls.F90 b/src/EOS_Omni_MultiVarCalls.F90 index ec35cfe..25be0a2 100644 --- a/src/EOS_Omni_MultiVarCalls.F90 +++ b/src/EOS_Omni_MultiVarCalls.F90 @@ -82,3 +82,123 @@ subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,& enddo end subroutine EOS_Omni_EOS_short + + +subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,& + rho,eps,temp,ye,dpderho,dpdrhoe,keyerr,anyerr) + + use EOS_Omni_Module + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(out) :: keyerr(npoints) + CCTK_INT, intent(out) :: anyerr + CCTK_REAL, intent(in) :: rf_precision + CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) + CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) + CCTK_REAL, intent(out) :: dpderho(npoints) + CCTK_REAL, intent(out) :: dpdrhoe(npoints) + + ! local vars + integer :: i + character(256) :: warnstring + real*8 :: hybrid_local_gamma + real*8 :: hybrid_local_k_cgs + real*8 :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2 + ! temporary vars for nuc_eos + real*8 :: xrho,xye,xtemp,xenr,xent + real*8 :: xprs,xmunu,xcs2 + real*8 :: xdedt,xdpderho,xdpdrhoe + + + anyerr = 0 + keyerr(:) = 0 + + select case (eoskey) + case (1) + ! polytropic EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) / & + (poly_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + dpdrhoe(i) = press_gf * poly_k_cgs * & + poly_gamma * inv_rho_gf * & + (rho(i)*inv_rho_gf) ** (poly_gamma - 1.d0) + dpderho(i) = 0.0d0 + enddo + case (2) + ! gamma-law EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & + (gl_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + dpdrhoe(i) = (gl_gamma-1.0d0) * & + eps(i) + dpderho(i) = (gl_gamma - 1.0d0) * & + rho(i) + enddo + case (3) + ! hybrid EOS + do i=1,npoints + if(rho(i).gt.hybrid_rho_nuc) then + hybrid_local_gamma = hybrid_gamma2 + hybrid_local_k_cgs = hybrid_k2_cgs + else + hybrid_local_gamma = hybrid_gamma1 + hybrid_local_k_cgs = hybrid_k1_cgs + endif + hybrid_dp_poly = hybrid_local_gamma * press_gf * & + hybrid_local_k_cgs * rho(i)**(hybrid_local_gamma - 1.0d0) * & + inv_rho_gf**hybrid_local_gamma + + hybrid_dp_th1 = - hybrid_local_gamma * press_gf * hybrid_local_k_cgs * & + (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) * & + rho(i)**(hybrid_local_gamma - 1.d0) * inv_rho_gf**hybrid_local_gamma + + hybrid_dp_th2 = (hybrid_gamma_th - 1.d0) * eps(i) & + - (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & + (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & + hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) + + dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2 + dpderho(i) = (hybrid_gamma_th - 1.0d0) * rho(i) + enddo + case (4) + do i=1,npoints + xrho = rho(i) * inv_rho_gf + xtemp = temp(i) + xye = ye(i) + xenr = eps(i) * inv_eps_gf + call nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& + xdpderho,& + keytemp,keyerr(i),rf_precision) + + if(keyerr(i).ne.0) then + anyerr = 1 + endif + + if(keytemp.eq.1) then + eps(i) = xenr * eps_gf + else + temp(i) = xtemp + endif + + dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf + dpderho(i) = xdpderho * press_gf * inv_eps_gf + enddo + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_WARN(0,warnstring) + end select + + end subroutine EOS_Omni_EOS_dpderho_dpdrhoe diff --git a/src/nuc_eos/findtemp.F90 b/src/nuc_eos/findtemp.F90 index d208a5f..646994d 100644 --- a/src/nuc_eos/findtemp.F90 +++ b/src/nuc_eos/findtemp.F90 @@ -16,10 +16,10 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps) integer :: rl = 0 integer itmax,i,keyerrt integer ii,jj,kk - + keyerrt=0 - tol=rfeps ! need to find energy to less than 1 in 10^-10 + tol=rfeps ! need to find energy to less than 1 in 10^10 itmax=20 ! use at most 20 iterations, then bomb lt=lt0 @@ -95,12 +95,14 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps) keyerrt=0 goto 12 else +#if 0 ! total failure write(*,*) "EOS: Did not converge in findtemp!" write(*,*) "rl,logrho,logtemp0,ye,lt,eps,eps0,abs(eps-eps0)/eps0" write(*,"(i4,i4,1P10E19.10)") i,rl,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0 write(*,*) "Tried calling bisection... didn't help... :-/" write(*,*) "Bisection error: ",keyerrt +#endif endif endif diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90 index e84ac82..399d89b 100644 --- a/src/nuc_eos/nuc_eos.F90 +++ b/src/nuc_eos/nuc_eos.F90 @@ -348,6 +348,91 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& end subroutine nuc_eos_press_eps +subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& + xdpderho,& + keytemp,keyerr,rfeps) + + use eosmodule + implicit none + + real*8, intent(in) :: xrho,xye + real*8, intent(inout) :: xtemp,xenr + real*8, intent(in) :: rfeps + real*8, intent(out) :: xdpdrhoe,xdpderho + integer, intent(in) :: keytemp + integer, intent(out) :: keyerr + + ! local variables + real*8 :: xcs2,xprs + real*8 :: lr,lt,y,xx,xeps,leps,xs + real*8 :: d1,d2,d3,ff(8) + integer :: keyerrt = 0 + + if(xrho.gt.eos_rhomax) then + stop "nuc_eos: rho > rhomax" + endif + + if(xrho.lt.eos_rhomin*1.2d0) then + call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + return + endif + + if(xye.gt.eos_yemax) then + stop "nuc_eos: ye > yemax" + endif + + if(xye.lt.eos_yemin) then + keyerr = 102 + write(6,*) "ye: ", xye + stop "nuc_eos ye < yemin" + endif + + if(keytemp.eq.1) then + if(xtemp.gt.eos_tempmax) then + stop "nuc_eos: temp > tempmax" + endif + + if(xtemp.lt.eos_tempmin) then + call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + return + endif + endif + + lr = log10(xrho) + lt = log10(xtemp) + y = xye + xeps = xenr + energy_shift + leps = log10(max(xeps,1.0d0)) + + keyerr = 0 + + if(keytemp.eq.0) then + !need to find temperature based on xeps + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + if(keyerrt.ne.0) then + keyerr = keyerrt + return + endif + xtemp = 10.0d0**lt + + elseif(keytemp.eq.2) then + stop "eos_nuc_press does not support keytemp.eq.2" + endif + + ! have temperature, proceed: + call findall_dpdr_dpde(lr,lt,y,ff) + xdpdrhoe = 10.0d0**ff(1) + xdpderho = 10.0d0**ff(2) + + + if((keytemp.eq.1).or.(keytemp.eq.2)) then + call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3) + xeps = 10.0d0**xeps - energy_shift + xenr = xeps + endif + +end subroutine nuc_eos_dpdr_dpde + subroutine findthis(lr,lt,y,value,array,d1,d2,d3) @@ -426,3 +511,22 @@ subroutine findall_press_eps(lr,lt,y,ff) ff(:) = ffx(:,1) end subroutine findall_press_eps + +subroutine findall_dpdr_dpde(lr,lt,y,ff) + + use eosmodule + implicit none + + real*8 ffx(8,1) + real*8 ff(8) + real*8 lr,lt,y + integer i + integer :: nvarsx = 2 + + +! Ewald's interpolator + call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,7:8), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_dpdr_dpde -- cgit v1.2.3