User: cott
Date: 2013/03/23 03:10 PM
Modified:
/trunk/src/nuc_eos/
findtemp.F90, linterp_many.F90, nuc_eos.F90
Log:
* extrapolate to negative (!) temperatures, so that we can use this
with the new Con2PrimHot treatment in GRHydro.
File Changes:
Directory: /trunk/src/nuc_eos/
==============================
File [modified]: findtemp.F90
Delta lines: +90 -17
===================================================================
--- trunk/src/nuc_eos/findtemp.F90 2013-03-23 20:08:12 UTC (rev 77)
+++ trunk/src/nuc_eos/findtemp.F90 2013-03-23 20:10:20 UTC (rev 78)
@@ -16,10 +16,10 @@
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 ! precision to which we need to find temp
itmax=20 ! use at most 20 iterations, then bomb
lt=lt0
@@ -28,7 +28,6 @@
eps0=epsin
eps1=eps0
-
ltmax=logtemp(ntemp)
ltmin=logtemp(1)
@@ -42,25 +41,17 @@
lt1=lt
eps1=eps
-! write(*,"(i4,1P12E19.10)") 0,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0,d2
-! write(*,"(i4,1P12E19.10)") 0,lr,lt0,y,ltmin,ltmax
-
do i=1,itmax
!d2 is the derivative deps/dlogtemp;
ldt = -(eps - eps0)/d2
-! if(ldt.gt.0.d0) ltmin=lt
-! if(ldt.le.0.d0) ltmax=lt
ltn = lt+ldt
ltn = min(ltn,ltmax)
ltn = max(ltn,ltmin)
lt1=lt
lt=ltn
eps1=eps
-! write(*,"(i4,1P12E19.10)") i,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0,d2,ldt
call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3)
-! call findthis(lr,lt,y,eps,epst,d1,d2,d3)
if (abs(eps - eps0).lt.tol*abs(eps0)) then
-! write(*,"(1P12E19.10)") tol,abs(eps-eps0)/eps0
exit
endif
!setup new d2
@@ -72,12 +63,6 @@
if(abs(eps-eps0).lt.1.0d-3*abs(eps0)) then
d2 = (eps-eps1)/(lt-lt1)
endif
-! if(i.ge.10) then
-! write(*,*) "EOS: Did not converge in findtemp!"
-! write(*,*) "rl,logrho,logtemp0,ye,lt,eps,eps0,abs(eps-eps0)/eps0"
-! if(i.gt.5) then
-! write(*,"(i4,1P10E22.14)") i,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0
-! endif
enddo
@@ -190,6 +175,7 @@
if(i.ge.itmax) then
keyerrt=667
call bisection(lr,lt0,y,s0,lt,alltables(:,:,:,3),keyerrt,2)
+#if 0
if(keyerrt.eq.667) then
write(*,*) "EOS: Did not converge in findtemp_entropy!"
write(*,*) "rl,logrho,logtemp0,ye,lt,s,s0,abs(s-s0)/s0"
@@ -197,6 +183,7 @@
write(*,*) "Tried calling bisection... didn't help... :-/"
write(*,*) "Bisection error: ",keyerrt
endif
+#endif
lt0=lt
return
@@ -207,3 +194,89 @@
end subroutine findtemp_entropy
+
+
+subroutine findtemp_low(lr,t0,y,epsin,keyerrt,rfeps)
+
+ ! this routine is for finding fake temperatures
+ ! outside of the table range; we use linear
+ ! extrapolation in temperature
+
+ use eosmodule
+
+ implicit none
+
+ real*8 lr,t0,y,epsin
+ real*8 eps,t,dt
+ real*8 tol
+ real*8 dlepsdt
+ real*8 eps0,eps1,t1
+
+ real*8 tn,tmin,tmax
+ real*8 tinput,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
+ itmax=100 ! use at most 20 iterations, then bomb
+
+ t=t0
+ t1=t
+
+ eps0=epsin
+ eps1=eps0
+
+ tmax=10.0d0**logtemp(2)
+ tmin=-20.0d0
+
+ !preconditioning 1: do we already have the right temperature?
+ call intp3d_linearTlow(lr,t,y,eps,1,alltables(:,:,:,2),nrho,&
+ ntemp,nye,logrho,logtemp,ye,dlepsdt)
+
+ if (abs(eps-eps0).lt.tol*abs(eps0)) then
+ return
+ endif
+ t1=t
+ eps1=eps
+
+ do i=1,itmax
+ dt = -(eps - eps0)/dlepsdt
+ tn = t+dt
+ if(tn >= tmax) then
+ tn = tmax*0.9d0
+ endif
+ tn = min(tn,tmax)
+ tn = max(tn,tmin)
+ t1=t
+ t=tn
+ eps1=eps
+
+ call intp3d_linearTlow(lr,t,y,eps,1,alltables(:,:,:,2),nrho,&
+ ntemp,nye,logrho,logtemp,ye,dlepsdt)
+
+ if (abs(eps - eps0).lt.tol*abs(eps0)) then
+ exit
+ endif
+
+ ! if we are closer than 10^-2 to the
+ ! root (eps-eps0)=0, we are switching to
+ ! the secant method, since the table is rather coarse and the
+ ! derivatives may be garbage.
+ if(abs(eps-eps0).lt.1.0d-3*abs(eps0)) then
+ dlepsdt = (eps-eps1)/(t-t1)
+ endif
+ enddo
+
+ t0 = t
+
+ if(i.ge.itmax) then
+ keyerrt=667
+ return
+ endif
+
+
+end subroutine findtemp_low
File [modified]: linterp_many.F90
Delta lines: +250 -4
===================================================================
--- trunk/src/nuc_eos/linterp_many.F90 2013-03-23 20:08:12 UTC (rev 77)
+++ trunk/src/nuc_eos/linterp_many.F90 2013-03-23 20:10:20 UTC (rev 78)
@@ -1,5 +1,3 @@
-#include "cctk.h"
-
SUBROUTINE intp3d_many ( x, y, z, f, kt, ft, nx, ny, nz, nvars, xt, yt, zt)
!
implicit none
@@ -45,7 +43,7 @@
real*8 dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi
integer n,ix,iy,iz
- IF (kt .GT. ktx) call CCTK_WARN(0, '***KTX**')
+ IF (kt .GT. ktx) STOP'***KTX**'
!
!
!------ determine spacing parameters of (equidistant!!!) table
@@ -67,7 +65,7 @@
!
!------- loop over all points to be interpolated
!
- do n = 1, kt
+ dO n = 1, kt
!
!------- determine location in (equidistant!!!) table
!
@@ -125,3 +123,251 @@
end SUBROUTINE intp3d_many
+
+ SUBROUTINE intp3d_linearTlow ( x, y, z, f, kt, ft, nx, ny, nz, xt, yt, zt, &
+ d2)
+!
+ implicit none
+!
+!---------------------------------------------------------------------
+!
+! purpose: interpolation of a function of three variables in an
+! equidistant(!!!) table.
+!
+! method: 8-point Lagrange linear interpolation formula
+!
+! x input vector of first variable
+! y temperature
+! z input vector of third variable
+!
+! f output vector of interpolated function values
+!
+! kt vector length of input and output vectors
+!
+! ft 3d array of tabulated function values
+! nx x-dimension of table
+! ny y-dimension of table
+! nz z-dimension of table
+! xt vector of x-coordinates of table
+! yt vector of y-coordinates of table
+! zt vector of z-coordinates of table
+!
+!---------------------------------------------------------------------
+
+ integer kt,nx,ny,nz
+ real*8 :: ft(nx,ny,nz)
+
+ real*8 x(kt),y(kt),z(kt),f(kt)
+ real*8 xt(nx),yt(ny),zt(nz)
+ real*8 d1,d2,d3
+!
+!
+ integer,parameter :: ktx = 1
+ real*8 fh(ktx,8), delx(ktx), dely(ktx), delz(ktx), &
+ a1(ktx), a2(ktx), a3(ktx), a4(ktx), &
+ a5(ktx), a6(ktx), a7(ktx), a8(ktx)
+
+ real*8 dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi
+ integer n,ix,iy,iz
+
+ IF (kt .GT. ktx) STOP'***KTX**'
+!
+!
+!------ determine spacing parameters of (equidistant!!!) table
+!
+ dx = (xt(nx) - xt(1)) / FLOAT(nx-1)
+ dy = 10.0d0**yt(2) - 10.0d0**yt(1)
+ dz = (zt(nz) - zt(1)) / FLOAT(nz-1)
+!
+ dxi = 1. / dx
+ dyi = 1. / dy
+ dzi = 1. / dz
+!
+ dxyi = dxi * dyi
+ dxzi = dxi * dzi
+ dyzi = dyi * dzi
+!
+ dxyzi = dxi * dyi * dzi
+!
+!
+!------- loop over all points to be interpolated
+!
+ dO n = 1, kt
+!
+!------- determine location in (equidistant!!!) table
+!
+ ix = 2 + INT( (x(n) - xt(1) - 1.e-10) * dxi )
+ iy = 2
+ iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi )
+!
+ ix = MAX( 2, MIN( ix, nx ) )
+ iz = MAX( 2, MIN( iz, nz ) )
+!
+! write(*,*) iy-1,iy,iy+1
+!
+!------- set-up auxiliary arrays for Lagrange interpolation
+!
+ delx(n) = xt(ix) - x(n)
+ dely(n) = 10.0d0**yt(iy) - y(n)
+ delz(n) = zt(iz) - z(n)
+!
+ fh(n,1) = ft(ix , iy , iz)
+ fh(n,2) = ft(ix-1, iy , iz)
+ fh(n,3) = ft(ix , iy-1, iz)
+ fh(n,4) = ft(ix , iy , iz-1)
+ fh(n,5) = ft(ix-1, iy-1, iz)
+ fh(n,6) = ft(ix-1, iy , iz-1)
+ fh(n,7) = ft(ix , iy-1, iz-1)
+ fh(n,8) = ft(ix-1, iy-1, iz-1)
+!
+!------ set up coefficients of the interpolation polynomial and
+! evaluate function values
+ !
+ a1(n) = fh(n,1)
+ a2(n) = dxi * ( fh(n,2) - fh(n,1) )
+ a3(n) = dyi * ( fh(n,3) - fh(n,1) )
+ a4(n) = dzi * ( fh(n,4) - fh(n,1) )
+ a5(n) = dxyi * ( fh(n,5) - fh(n,2) - fh(n,3) + fh(n,1) )
+ a6(n) = dxzi * ( fh(n,6) - fh(n,2) - fh(n,4) + fh(n,1) )
+ a7(n) = dyzi * ( fh(n,7) - fh(n,3) - fh(n,4) + fh(n,1) )
+ a8(n) = dxyzi * ( fh(n,8) - fh(n,1) + fh(n,2) + fh(n,3) + &
+ fh(n,4) - fh(n,5) - fh(n,6) - fh(n,7) )
+!
+ d2 = -a3(n)
+ f(n) = a1(n) + a2(n) * delx(n) &
+ + a3(n) * dely(n) &
+ + a4(n) * delz(n) &
+ + a5(n) * delx(n) * dely(n) &
+ + a6(n) * delx(n) * delz(n) &
+ + a7(n) * dely(n) * delz(n) &
+ + a8(n) * delx(n) * dely(n) * delz(n)
+!
+
+ enddo
+!
+
+ end SUBROUTINE intp3d_linearTlow
+
+ SUBROUTINE intp3d_many_linearTLow ( x, y, z, f, kt, ft, nx, ny, nz, nvars, xt, yt, zt)
+!
+ implicit none
+!
+!---------------------------------------------------------------------
+!
+! purpose: interpolation of a function of three variables in an
+! equidistant(!!!) table.
+!
+! method: 8-point Lagrange linear interpolation formula
+!
+! x input vector of first variable
+! y input vector of second variable
+! z input vector of third variable
+!
+! f output vector of interpolated function values
+!
+! kt vector length of input and output vectors
+!
+! ft 3d array of tabulated function values
+! nx x-dimension of table
+! ny y-dimension of table
+! nz z-dimension of table
+! xt vector of x-coordinates of table
+! yt vector of y-coordinates of table
+! zt vector of z-coordinates of table
+!
+!---------------------------------------------------------------------
+
+ integer kt,nx,ny,nz,iv,nvars
+ real*8 :: ft(nx,ny,nz,nvars)
+
+ real*8 x(kt),y(kt),z(kt),f(kt,nvars)
+ real*8 xt(nx),yt(ny),zt(nz)
+ real*8 d1,d2,d3
+!
+!
+ integer,parameter :: ktx = 1
+ real*8 fh(ktx,8,nvars), delx(ktx), dely(ktx), delz(ktx), &
+ a1(ktx,nvars), a2(ktx,nvars), a3(ktx,nvars), a4(ktx,nvars), &
+ a5(ktx,nvars), a6(ktx,nvars), a7(ktx,nvars), a8(ktx,nvars)
+
+ real*8 dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi
+ integer n,ix,iy,iz
+
+ IF (kt .GT. ktx) STOP'***KTX**'
+!
+!
+!------ determine spacing parameters of (equidistant!!!) table
+!
+ dx = (xt(nx) - xt(1)) / FLOAT(nx-1)
+ dy = 10.0d0**yt(2) - 10.0d0**yt(1)
+ dz = (zt(nz) - zt(1)) / FLOAT(nz-1)
+!
+ dxi = 1. / dx
+ dyi = 1. / dy
+ dzi = 1. / dz
+!
+ dxyi = dxi * dyi
+ dxzi = dxi * dzi
+ dyzi = dyi * dzi
+!
+ dxyzi = dxi * dyi * dzi
+!
+!
+!------- loop over all points to be interpolated
+!
+ dO n = 1, kt
+!
+!------- determine location in (equidistant!!!) table
+!
+ ix = 2 + INT( (x(n) - xt(1) - 1.e-10) * dxi )
+ iy = 2
+ iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi )
+!
+ ix = MAX( 2, MIN( ix, nx ) )
+ iz = MAX( 2, MIN( iz, nz ) )
+!
+! write(*,*) iy-1,iy,iy+1
+!
+!------- set-up auxiliary arrays for Lagrange interpolation
+!
+ delx(n) = xt(ix) - x(n)
+ dely(n) = 10.0d0**yt(iy) - y(n)
+ delz(n) = zt(iz) - z(n)
+!
+ do iv = 1, nvars
+ fh(n,1,iv) = ft(ix , iy , iz, iv )
+ fh(n,2,iv) = ft(ix-1, iy , iz, iv )
+ fh(n,3,iv) = ft(ix , iy-1, iz, iv )
+ fh(n,4,iv) = ft(ix , iy , iz-1, iv)
+ fh(n,5,iv) = ft(ix-1, iy-1, iz, iv )
+ fh(n,6,iv) = ft(ix-1, iy , iz-1, iv)
+ fh(n,7,iv) = ft(ix , iy-1, iz-1, iv)
+ fh(n,8,iv) = ft(ix-1, iy-1, iz-1, iv)
+!
+!------ set up coefficients of the interpolation polynomial and
+! evaluate function values
+ !
+ a1(n,iv) = fh(n,1,iv)
+ a2(n,iv) = dxi * ( fh(n,2,iv) - fh(n,1,iv) )
+ a3(n,iv) = dyi * ( fh(n,3,iv) - fh(n,1,iv) )
+ a4(n,iv) = dzi * ( fh(n,4,iv) - fh(n,1,iv) )
+ a5(n,iv) = dxyi * ( fh(n,5,iv) - fh(n,2,iv) - fh(n,3,iv) + fh(n,1,iv) )
+ a6(n,iv) = dxzi * ( fh(n,6,iv) - fh(n,2,iv) - fh(n,4,iv) + fh(n,1,iv) )
+ a7(n,iv) = dyzi * ( fh(n,7,iv) - fh(n,3,iv) - fh(n,4,iv) + fh(n,1,iv) )
+ a8(n,iv) = dxyzi * ( fh(n,8,iv) - fh(n,1,iv) + fh(n,2,iv) + fh(n,3,iv) + &
+ fh(n,4,iv) - fh(n,5,iv) - fh(n,6,iv) - fh(n,7,iv) )
+!
+ f(n,iv) = a1(n,iv) + a2(n,iv) * delx(n) &
+ + a3(n,iv) * dely(n) &
+ + a4(n,iv) * delz(n) &
+ + a5(n,iv) * delx(n) * dely(n) &
+ + a6(n,iv) * delx(n) * delz(n) &
+ + a7(n,iv) * dely(n) * delz(n) &
+ + a8(n,iv) * delx(n) * dely(n) * delz(n)
+!
+ enddo
+
+ enddo
+!
+
+ end SUBROUTINE intp3d_many_linearTlow
File [modified]: nuc_eos.F90
Delta lines: +142 -42
===================================================================
--- trunk/src/nuc_eos/nuc_eos.F90 2013-03-23 20:08:12 UTC (rev 77)
+++ trunk/src/nuc_eos/nuc_eos.F90 2013-03-23 20:10:20 UTC (rev 78)
@@ -69,16 +69,18 @@
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 @@
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 @@
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 @@
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 @@
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 @@
! 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 @@
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 @@
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 @@
! 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 @@
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 @@
! 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 @@
end subroutine findall_short
+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(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
@@ -577,13 +639,32 @@
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 @@
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
