diff options
Diffstat (limited to 'src/nuc_eos/nuc_eos.F90')
| -rw-r--r-- | src/nuc_eos/nuc_eos.F90 | 184 |
1 files changed, 142 insertions, 42 deletions
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90 index 3420ee0..1dfb9e8 100644 --- a/src/nuc_eos/nuc_eos.F90 +++ b/src/nuc_eos/nuc_eos.F90 @@ -69,16 +69,18 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& if(keytemp.eq.1) then if(xtemp.gt.eos_tempmax) then - call CCTK_WARN (0, "nuc_eos: temp > tempmax") + keyerr = 105 + return endif if(xtemp.lt.eos_tempmin) then - call CCTK_WARN (0, "nuc_eos: temp < tempmin") + keyerr = 106 + return endif endif lr = log10(xrho) - lt = log10(xtemp) + lt = log10(max(xtemp,eos_tempmin)) y = xye xeps = xenr + energy_shift leps = log10(max(xeps,1.0d0)) @@ -221,7 +223,7 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& integer, intent(out) :: keyerr ! local variables - real*8 :: lr,lt,y,xx,xeps,leps,xs,xpressure + real*8 :: lr,lt,y,xt,xx,xeps,leps,xs,xpressure real*8 :: d1,d2,d3,ff(8) integer :: keyerrt integer :: keyerrr @@ -256,14 +258,13 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& endif if(xtemp.lt.eos_tempmin) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) - xent = 0.0d0 + keyerr = 106 return endif endif lr = log10(xrho) - lt = log10(xtemp) + lt = log10(max(xtemp,eos_tempmin)) y = xye keyerr = 0 @@ -273,12 +274,22 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& xeps = xenr + energy_shift leps = log10(max(xeps,1.0d0)) call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.ne.0) then + if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + ! turns out, we can still converge, but + ! too a bogus, possibly negative temperature + xt = xtemp + call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + if(keyerrt.eq.0) then + keyerrt = 668 + xtemp = xt + else + keyerrt = 667 + return + endif keyerr = keyerrt - return + else + xtemp = 10.0d0**lt endif - xtemp = 10.0d0**lt - elseif(keytemp.eq.2) then !need to find temperature based on xent xs = xent @@ -302,7 +313,11 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& endif ! have rho,temp,ye; proceed: - call findall_short(lr,lt,y,ff) + if(keyerr.ne.668) then + call findall_short(lr,lt,y,ff) + else + call findall_short_lowT(lr,xt,y,ff) + endif !unless we want xprs to be constant (keytemp==3), reset xprs if(.not.keytemp.eq.3) then @@ -347,7 +362,7 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& ! local variables real*8 :: xcs2,xdpderho,xdpdrhoe - real*8 :: lr,lt,y,xx,xeps,leps,xs + real*8 :: lr,lt,y,xt,xx,xeps,leps,xs real*8 :: d1,d2,d3,ff(8) integer :: keyerrt @@ -379,7 +394,7 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& endif if(xtemp.lt.eos_tempmin) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + keyerr = 106 return endif endif @@ -390,29 +405,44 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& call CCTK_WARN (0, "eos_nuc_press does not support keytemp other than 0 and 1") endif - lr = log10(xrho) - lt = log10(xtemp) + lr = log10(xrho) + lt = log10(max(xtemp,eos_tempmin)) y = xye if(keytemp.eq.0) then !need to find temperature based on xeps xeps = xenr + energy_shift leps = log10(max(xeps,1.0d0)) call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.ne.0) then + if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + ! turns out, we can still converge, but + ! too a bogus, possibly negative temperature + xt = xtemp + call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + if(keyerrt.eq.0) then + keyerrt = 668 + xtemp = xt + else + keyerrt = 667 + return + endif keyerr = keyerrt - return + else + xtemp = 10.0d0**lt endif - xtemp = 10.0d0**lt endif ! have temperature, proceed: - call findall_press_eps(lr,lt,y,ff) + if(keyerr.ne.668) then + call findall_press_eps(lr,lt,y,ff) + else + call findall_press_eps_lowT(lr,xt,y,ff) + endif xprs = 10.0d0**ff(1) xenr = 10.0d0**ff(2) - energy_shift end subroutine nuc_eos_press_eps -subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& - xdpderho,& +subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpderho,& + xdpdrhoe,& keytemp,keyerr,rfeps) use eosmodule @@ -427,8 +457,8 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& ! local variables real*8 :: xcs2,xprs - real*8 :: lr,lt,y,xx,xeps,leps,xs - real*8 :: d1,d2,d3,ff(8) + real*8 :: lr,lt,y,xt,xx,xeps,leps,xs + real*8 :: d1,d2,d3,ff(2) integer :: keyerrt keyerrt = 0 @@ -459,38 +489,52 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& endif if(xtemp.lt.eos_tempmin) then - call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + keyerr = 106 return endif endif lr = log10(xrho) - lt = log10(xtemp) + lt = log10(max(xtemp,eos_tempmin)) 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 + xeps = xenr + energy_shift + leps = log10(max(xeps,1.0d0)) call findtemp(lr,lt,y,leps,keyerrt,rfeps) - if(keyerrt.ne.0) then + if(keyerrt.eq.667 .and. lt .lt. 10.0d0) then + ! turns out, we can still converge, but + ! too a bogus, possibly negative temperature + xt = xtemp + call findtemp_low(lr,xt,y,leps,keyerrt,rfeps) + if(keyerrt.eq.0) then + keyerrt = 668 + xtemp = xt + else + keyerrt = 667 + return + endif keyerr = keyerrt - return + else + xtemp = 10.0d0**lt endif - xtemp = 10.0d0**lt - elseif(keytemp.gt.1) then call CCTK_WARN (0, "eos_nuc_press does not support keytemp > 1") endif ! have temperature, proceed: - call findall_dpdr_dpde(lr,lt,y,ff) - xdpdrhoe = 10.0d0**ff(1) - xdpderho = 10.0d0**ff(2) + if(keyerr.ne.668) then + call findall_dpdr_dpde(lr,lt,y,ff) + else + call findall_dpdr_dpde_lowT(lr,xtemp,y,ff) + endif - if(keytemp.eq.1) then + xdpdrhoe = ff(1) + xdpderho = ff(2) + + if(keytemp.eq.1.and.keyerr.ne.668) then call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3) xeps = 10.0d0**xeps - energy_shift xenr = xeps @@ -515,7 +559,6 @@ subroutine findthis(lr,lt,y,value,array,d1,d2,d3) ! Ewald's interpolator call intp3d(lr,lt,y,value,1,array,nrho,ntemp,nye,logrho,logtemp,ye,d1,d2,d3) - end subroutine findthis @@ -557,14 +600,33 @@ subroutine findall_short(lr,lt,y,ff) end subroutine findall_short - -subroutine findall_press_eps(lr,lt,y,ff) +subroutine findall_short_lowT(lr,t,y,ff) use eosmodule implicit none real*8 ffx(8,1) real*8 ff(8) + real*8 lr,t,y + integer i + integer :: nvarsx = 8 + + +! Ewald's interpolator + call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,1:8), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_short_lowT + + +subroutine findall_press_eps(lr,lt,y,ff) + + use eosmodule + implicit none + + real*8 ffx(2,1) + real*8 ff(2) real*8 lr,lt,y integer i integer :: nvarsx = 2 @@ -577,13 +639,32 @@ subroutine findall_press_eps(lr,lt,y,ff) end subroutine findall_press_eps +subroutine findall_press_eps_lowT(lr,t,y,ff) + + use eosmodule + implicit none + + real*8 ffx(2,1) + real*8 ff(2) + real*8 lr,t,y + integer i + integer :: nvarsx = 2 + +! Ewald's interpolator + call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,1:2), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_press_eps_lowT + + subroutine findall_dpdr_dpde(lr,lt,y,ff) use eosmodule implicit none - real*8 ffx(8,1) - real*8 ff(8) + real*8 ffx(2,1) + real*8 ff(2) real*8 lr,lt,y integer i integer :: nvarsx = 2 @@ -595,3 +676,22 @@ subroutine findall_dpdr_dpde(lr,lt,y,ff) ff(:) = ffx(:,1) end subroutine findall_dpdr_dpde + +subroutine findall_dpdr_dpde_lowT(lr,t,y,ff) + + use eosmodule + implicit none + + real*8 ffx(2,1) + real*8 ff(2) + real*8 lr,t,y + integer i + integer :: nvarsx = 2 + + +! Ewald's interpolator + call intp3d_many_linearTlow(lr,t,y,ffx,1,alltables(:,:,:,7:8), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_dpdr_dpde_lowT |