From 4530e894a8002d16f7eb7c822a79c4038c6b7682 Mon Sep 17 00:00:00 2001 From: cott Date: Sat, 23 Mar 2013 20:10:20 +0000 Subject: * extrapolate to negative (!) temperatures, so that we can use this with the new Con2PrimHot treatment in GRHydro. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@78 8e189c6b-2ab8-4400-aa02-70a9cfce18b9 --- src/nuc_eos/findtemp.F90 | 107 +++++++++++++++--- src/nuc_eos/linterp_many.F90 | 254 ++++++++++++++++++++++++++++++++++++++++++- src/nuc_eos/nuc_eos.F90 | 184 ++++++++++++++++++++++++------- 3 files changed, 482 insertions(+), 63 deletions(-) diff --git a/src/nuc_eos/findtemp.F90 b/src/nuc_eos/findtemp.F90 index 646994d..125b0cd 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 ! precision to which we need to find temp itmax=20 ! use at most 20 iterations, then bomb lt=lt0 @@ -28,7 +28,6 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps) eps0=epsin eps1=eps0 - ltmax=logtemp(ntemp) ltmin=logtemp(1) @@ -42,25 +41,17 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps) 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 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps) 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 @@ subroutine findtemp_entropy(lr,lt0,y,sin,keyerrt,rfeps) 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 @@ subroutine findtemp_entropy(lr,lt0,y,sin,keyerrt,rfeps) write(*,*) "Tried calling bisection... didn't help... :-/" write(*,*) "Bisection error: ",keyerrt endif +#endif lt0=lt return @@ -207,3 +194,89 @@ subroutine findtemp_entropy(lr,lt0,y,sin,keyerrt,rfeps) 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 diff --git a/src/nuc_eos/linterp_many.F90 b/src/nuc_eos/linterp_many.F90 index d0ddb08..42b819f 100644 --- a/src/nuc_eos/linterp_many.F90 +++ b/src/nuc_eos/linterp_many.F90 @@ -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 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 -- cgit v1.2.3