User: cott Date: 2010/10/22 02:42 PM
Added: /EOS_Omni/ configuration.ccl /EOS_Omni/src/nuc_eos/ bisection.F90, eosmodule.F90, findtemp.F90, linterp_many.F90, make.code.defn, make.code.deps, nuc_eos.F90, readtable.F90
Modified: /EOS_Omni/src/ make.code.defn
Log: * add finite-T EOS routines (not called yet; but code compiles!)
File Changes:
Directory: /EOS_Omni/ =====================
File [added]: configuration.ccl Delta lines: +2 -0 =================================================================== --- EOS_Omni/configuration.ccl (rev 0) +++ EOS_Omni/configuration.ccl 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,2 @@ +REQUIRES HDF5 +
Directory: /EOS_Omni/src/nuc_eos/ =================================
File [added]: bisection.F90 Delta lines: +82 -0 =================================================================== --- EOS_Omni/src/nuc_eos/bisection.F90 (rev 0) +++ EOS_Omni/src/nuc_eos/bisection.F90 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,82 @@ +subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect) + + use eosmodule + + implicit none + + real*8 lr,lt0,y,eps0,lt + integer keyerrt + + integer keybisect + ! 1 -> operate on energy + ! 2 -> operate on entropy + + !temporary vars + real*8 lt1,lt2,ltmin,ltmax + real*8 f1,f2,fmid,dlt,ltmid + real*8 d1,d2,d3,tol + real*8 f1a,f2a + + integer bcount,i,itmax,maxbcount + + real*8 bivar + + bcount = 0 + maxbcount = 80 + + tol=1.d-9 ! need to find energy to less than 1 in 10^-9 + itmax=50 + + ltmax=logtemp(ntemp) + ltmin=logtemp(1) + + lt = lt0 + lt1 = dlog10(min(10.0d0**ltmax,1.10d0*(10.0d0**lt0))) + lt2 = dlog10(max(10.0d0**ltmin,0.90d0*(10.0d0**lt0))) + + call findthis(lr,lt1,y,f1a,bivar,d1,d2,d3) + call findthis(lr,lt2,y,f2a,bivar,d1,d2,d3) + + f1=f1a-eps0 + f2=f2a-eps0 + + keyerrt=0 + + + do while(f1*f2.ge.0.0d0) + bcount=bcount+1 + lt1=dlog10(min(10.0d0**ltmax,1.2d0*(10.0d0**lt1))) + lt2=dlog10(max(10.0d0**ltmin,0.8d0*(10.0d0**lt2))) + call findthis(lr,lt1,y,f1a,bivar,d1,d2,d3) + call findthis(lr,lt2,y,f2a,bivar,d1,d2,d3) + f1=f1a-eps0 + f2=f2a-eps0 + if(bcount.ge.maxbcount) then + keyerrt=667 + return + endif + enddo + + if(f1.lt.0.0d0) then + lt=lt1 + dlt=dlog10((10.0D0**lt2)-(10.0d0**lt1)) + else + lt=lt2 + dlt=dlog10((10.0D0**lt1)-(10.0d0**lt2)) + endif + + do i=1,itmax + dlt=dlog10((10.0d0**dlt)*0.5d0) + ltmid=dlog10(10.0d0**lt+10.0d0**dlt) + call findthis(lr,ltmid,y,f2a,bivar,d1,d2,d3) + fmid=f2a-eps0 + if(fmid.le.0.0d0) lt=ltmid + if(abs(1.0d0-f2a/eps0).lt.tol) then + lt=ltmid + return + endif + enddo + + + +end subroutine bisection
File [added]: eosmodule.F90 Delta lines: +60 -0 =================================================================== --- EOS_Omni/src/nuc_eos/eosmodule.F90 (rev 0) +++ EOS_Omni/src/nuc_eos/eosmodule.F90 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,60 @@ + module eosmodule + + implicit none + + integer,save :: nrho,ntemp,nye + + integer,save :: warn_from !warn from given reflevel + + real*8 :: energy_shift = 0.0d0 + + real*8 :: precision = 1.0d-9 + +! min-max values: + real*8 :: eos_rhomin,eos_rhomax + real*8 :: eos_yemin,eos_yemax + real*8 :: eos_tempmin,eos_tempmax + + real*8 :: t_max_hack = 240.0d0 + +! basics + integer, parameter :: nvars = 19 + real*8,allocatable :: alltables(:,:,:,:) + ! index variable mapping: + ! 1 -> logpress + ! 2 -> logenergy + ! 3 -> entropy + ! 4 -> munu + ! 5 -> cs2 + ! 6 -> dedT + ! 7 -> dpdrhoe + ! 8 -> dpderho + ! 9 -> muhat + ! 10 -> mu_e + ! 11 -> mu_p + ! 12 -> mu_n + ! 13 -> xa + ! 14 -> xh + ! 15 -> xn + ! 16 -> xp + ! 17 -> abar + ! 18 -> zbar + ! 19 -> gamma + + real*8,allocatable,save :: logrho(:) + real*8,allocatable,save :: logtemp(:) + real*8,allocatable,save :: ye(:) + +! constants + real*8,save :: mev_to_erg = 1.60217733d-6 + real*8,save :: amu_cgs = 1.66053873d-24 + real*8,save :: amu_mev = 931.49432d0 + real*8,save :: pi = 3.14159265358979d0 + real*8,save :: ggrav = 6.672d-8 + real*8,save :: temp_mev_to_kelvin = 1.1604447522806d10 + real*8,save :: clight = 2.99792458d10 + real*8,save :: kb_erg = 1.380658d-16 + real*8,save :: kb_mev = 8.61738568d-11 + + + end module eosmodule
File [added]: findtemp.F90 Delta lines: +207 -0 =================================================================== --- EOS_Omni/src/nuc_eos/findtemp.F90 (rev 0) +++ EOS_Omni/src/nuc_eos/findtemp.F90 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,207 @@ +subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps) + + use eosmodule + + implicit none + + real*8 lr,lt0,y,epsin + real*8 eps,lt,ldt + real*8 tol + real*8 d1,d2,d3 + real*8 eps0,eps1,lt1 + + real*8 ltn,ltmax,ltmin + 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=20 ! use at most 20 iterations, then bomb + + lt=lt0 + lt1=lt + + eps0=epsin + eps1=eps0 + + + ltmax=logtemp(ntemp) + ltmin=logtemp(1) + + ! Note: We are using Ewald's Lagrangian interpolator here! + + !preconditioning 1: do we already have the right temperature? + call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3) + if (abs(eps-eps0).lt.tol*abs(eps0)) then + return + endif + 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 + + ! 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 + 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 + + + if(i.ge.itmax) then + keyerrt=667 + call bisection(lr,lt0,y,eps0,lt,alltables(:,:,:,2),keyerrt,1) + if(keyerrt.eq.667) then + if(lt.ge.log10(t_max_hack)) then + ! handling for too high temperatures + lt = log10(t_max_hack) + keyerrt=0 + goto 12 + else if(abs(lt-log10(t_max_hack))/log10(t_max_hack).lt.0.025d0) then + lt0 = min(lt,log10(t_max_hack)) + keyerrt=0 + goto 12 + else + ! 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 + + lt0=min(lt,log10(t_max_hack)) + return + endif + +12 continue + + lt0=min(lt,log10(t_max_hack)) + + +end subroutine findtemp + +subroutine findtemp_entropy(lr,lt0,y,sin,keyerrt,rfeps) + +! This routine finds the new temperature based +! on rho, Y_e, entropy + + use eosmodule + + implicit none + + real*8 lr,lt0,y,sin + real*8 s,lt,ldt + real*8 tol + real*8 d1,d2,d3 + real*8 s0,s1,lt1 + + real*8 ltn,ltmax,ltmin + 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=20 ! use at most 20 iterations, then bomb + + lt=lt0 + lt1=lt + + s0=sin + s1=s0 + + ltmax=logtemp(ntemp) + ltmin=logtemp(1) + + !preconditioning 1: do we already have the right temperature? + call findthis(lr,lt,y,s,alltables(:,:,:,3),d1,d2,d3) + if (abs(s-s0).lt.tol*abs(s0)) then + return + endif + lt1=lt + s1=s + + + do i=1,itmax + !d2 is the derivative ds/dlogtemp; + ldt = -(s - s0)/d2 + ltn = lt+ldt + ltn = min(ltn,ltmax) + ltn = max(ltn,ltmin) + lt1=lt + lt=ltn + s1=s + call findthis(lr,lt,y,s,alltables(:,:,:,3),d1,d2,d3) + if (abs(s - s0).lt.tol*abs(s0)) then + exit + endif + !setup new d2 + + ! 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(s-s0).lt.1.0d-3*abs(s0)) then + d2 = (s-s1)/(lt-lt1) + endif + enddo + + + if(i.ge.itmax) then + keyerrt=667 + call bisection(lr,lt0,y,s0,lt,alltables(:,:,:,3),keyerrt,2) + 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" + write(*,"(i4,i4,1P10E19.10)") i,rl,lr,lt0,y,lt,s,s0,abs(s-s0)/s0 + write(*,*) "Tried calling bisection... didn't help... :-/" + write(*,*) "Bisection error: ",keyerrt + endif + + lt0=lt + return + endif + + + lt0=lt + + +end subroutine findtemp_entropy
File [added]: linterp_many.F90 Delta lines: +125 -0 =================================================================== --- EOS_Omni/src/nuc_eos/linterp_many.F90 (rev 0) +++ EOS_Omni/src/nuc_eos/linterp_many.F90 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,125 @@ + SUBROUTINE intp3d_many ( 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 = (yt(ny) - yt(1)) / FLOAT(ny-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 + INT( (y(n) - yt(1) - 1.e-10) * dyi ) + iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi ) +! + ix = MAX( 2, MIN( ix, nx ) ) + iy = MAX( 2, MIN( iy, ny ) ) + 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) = 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 +
File [added]: make.code.defn Delta lines: +7 -0 =================================================================== --- EOS_Omni/src/nuc_eos/make.code.defn (rev 0) +++ EOS_Omni/src/nuc_eos/make.code.defn 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,7 @@ +SRCS = eosmodule.F90 nuc_eos.F90 bisection.F90 \ + findtemp.F90 linterp_many.F90 readtable.F90 + +SUBDIRS = + +F90FLAGS+=-I$(HDF5_INC_DIRS) -I$(HDF5_LIB_DIRS) +
File [added]: make.code.deps Delta lines: +6 -0 =================================================================== --- EOS_Omni/src/nuc_eos/make.code.deps (rev 0) +++ EOS_Omni/src/nuc_eos/make.code.deps 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,6 @@ +bisection.F90.o: eosmodule.F90.o +findtemp.F90.o: eosmodule.F90.o +linterp_many.F90.o: eosmodule.F90.o +nuc_eos.F90.o: eosmodule.F90.o +readtable.F90.o: eosmodule.F90.o +
File [added]: nuc_eos.F90 Delta lines: +322 -0 =================================================================== --- EOS_Omni/src/nuc_eos/nuc_eos.F90 (rev 0) +++ EOS_Omni/src/nuc_eos/nuc_eos.F90 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,322 @@ +! ######################################################### +! +! Copyright C. D. Ott and Evan O'Connor, July 2009 +! +! UNITS: density g/cm^3 +! temperature MeV +! ye number fraction per baryon +! energy erg/g +! pressure dyn/cm^2 +! chemical potentials MeV +! entropy k_B / baryon +! cs2 cm^2/s^2 (not relativistic) +! +! keyerr --> error output; should be 0 +! rfeps --> root finding relative accuracy, set around 1.0d-10 +! keytemp: 0 -> coming in with eps +! 1 -> coming in with temperature +! 2 -> coming in with entropy +! +subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& + xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp,xabar,xzbar,xmu_e,xmu_n,xmu_p, & + xmuhat,keytemp,keyerr,rfeps) + + use eosmodule + implicit none + + real*8, intent(in) :: xrho,xye + real*8, intent(inout) :: xtemp,xenr,xent + real*8, intent(out) :: xprs,xcs2,xdedt + real*8, intent(out) :: xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp + real*8, intent(out) :: xabar,xzbar,xmu_e,xmu_n,xmu_p,xmuhat + real*8, intent(in) :: rfeps + integer, intent(in) :: keytemp + integer, intent(out) :: keyerr + + + ! local variables + real*8 :: lr,lt,y,xx,xeps,leps,xs + real*8 :: d1,d2,d3 + real*8 :: ff(nvars) + integer :: keyerrt = 0 + + if(xrho.gt.eos_rhomax) then + stop "nuc_eos: rho > rhomax" + endif + + if(xrho.lt.eos_rhomin) then + stop "nuc_eos: rho < rhomin" + endif + + if(xye.gt.eos_yemax) then + stop "nuc_eos: ye > yemax" + endif + + if(xye.lt.eos_yemin) then + 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 + stop "nuc_eos: temp < tempmin" + 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 + stop "Did not find temperature" + endif + xtemp = 10.0d0**lt + + elseif(keytemp.eq.2) then + !need to find temperature based on xent + xs = xent + call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps) + xtemp = 10.0d0**lt + endif + + ! have temperature, proceed: + call findall(lr,lt,y,ff) + + xprs = 10.0d0**ff(1) + + xeps = 10.0d0**ff(2) - energy_shift + if(keytemp.eq.1.or.keytemp.eq.2) then + xenr = xeps + endif + + if(keytemp.eq.1.or.keytemp.eq.0) then + xent = ff(3) + endif + + xcs2 = ff(5) + +! derivatives + xdedt = ff(6) + + xdpdrhoe = ff(7) + + xdpderho = ff(8) + +! chemical potentials + xmuhat = ff(9) + + xmu_e = ff(10) + + xmu_p = ff(11) + + xmu_n = ff(12) + +! compositions + xxa = ff(13) + + xxh = ff(14) + + xxn = ff(15) + + xxp = ff(16) + + xabar = ff(17) + + xzbar = ff(18) + + +end subroutine nuc_eos_full + +subroutine nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + + implicit none + real*8 xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe + real*8,parameter :: idealK1 = 1.2435d15 * (0.5d0**(4.d0/3.d0)) + real*8,parameter :: idealgamma = 1.41d0 + integer keytemp + + if(keytemp.eq.1) then +! energy wanted + xprs=idealK1*xrho**(idealgamma) + xenr=xprs/xrho/(idealgamma-1.d0) + xcs2=idealgamma*xprs/xrho + endif + + xprs = (idealgamma - 1.0d0) * xrho * xenr + + xdpderho = (idealgamma - 1.0d0 ) * xrho + xdpdrhoe = (idealgamma - 1.0d0 ) * xenr + + xcs2= xdpdrhoe+xdpderho*xprs/xrho**2 + +end subroutine nuc_low_eos + +subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,& + xdpderho,xdpdrhoe,xmunu,keytemp,keyerr,rfeps) + + use eosmodule + implicit none + + real*8, intent(in) :: xrho,xye + real*8, intent(inout) :: xtemp,xenr,xent + real*8, intent(in) :: rfeps + real*8, intent(out) :: xprs,xmunu,xcs2,xdedt + real*8, intent(out) :: xdpderho,xdpdrhoe + integer, intent(in) :: keytemp + integer, intent(out) :: keyerr + + ! local variables + 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) + xent = 4.0d0 + return + endif + + if(xye.gt.eos_yemax) then + stop "nuc_eos: ye > yemax" + endif + + if(xye.lt.eos_yemin) then + 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) + xent = 4.0d0 + 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 + !need to find temperature based on xent + xs = xent + call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps) + if(keyerrt.ne.0) then + keyerr = keyerrt + return + endif + xtemp = 10.0d0**lt + endif + + ! have temperature, proceed: + call findall_short(lr,lt,y,ff) + xprs = 10.0d0**ff(1) + + xeps = 10.0d0**ff(2) - energy_shift + if((keytemp.eq.1).or.(keytemp.eq.2)) then + xenr = xeps + endif + + if(keytemp.eq.1.or.keytemp.eq.0) then + xent = ff(3) + endif + + xmunu = ff(4) + + xcs2 = ff(5) + + xdedt = ff(6) + + xdpdrhoe = ff(7) + + xdpderho = ff(8) + +end subroutine nuc_eos_short + +subroutine findthis(lr,lt,y,value,array,d1,d2,d3) + + use eosmodule + + implicit none + + integer rip,rim + integer tip,tim + integer yip,yim + + real*8 lr,lt,y,value,d1,d2,d3 + real*8 array(*) + +! Ewald's interpolator + call intp3d(lr,lt,y,value,1,array,nrho,ntemp,nye,logrho,logtemp,ye,d1,d2,d3) + + +end subroutine findthis + + +subroutine findall(lr,lt,y,ff) + + use eosmodule + implicit none + + real*8 ff(nvars) + real*8 ffx(nvars,1) + real*8 lr,lt,y + integer i + +! Ewald's interpolator + call intp3d_many(lr,lt,y,ffx,1,alltables,& + nrho,ntemp,nye,nvars,logrho,logtemp,ye) + ff(:) = ffx(:,1) + + +end subroutine findall + + +subroutine findall_short(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 = 8 + + +! Ewald's interpolator + call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,1:8), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_short
File [added]: readtable.F90 Delta lines: +234 -0 =================================================================== --- EOS_Omni/src/nuc_eos/readtable.F90 (rev 0) +++ EOS_Omni/src/nuc_eos/readtable.F90 2010-10-22 19:42:05 UTC (rev 9) @@ -0,0 +1,234 @@ +subroutine readtable(eos_filename) +! This routine reads the table and initializes +! all variables in the module. + + use eosmodule + use hdf5 + + implicit none + + character(*) eos_filename + + character(len=100) message + +! HDF5 vars + integer(HID_T) file_id,dset_id,dspace_id + integer(HSIZE_T) dims1(1), dims3(3) + integer error,rank,accerr + integer i,j,k + + real*8 amu_cgs_andi + real*8 buffer1,buffer2,buffer3,buffer4 + accerr=0 + + write(*,*) "Reading Ott EOS Table" + + call h5open_f(error) + + call h5fopen_f (trim(adjustl(eos_filename)), H5F_ACC_RDONLY_F, file_id, error) + + write(6,*) trim(adjustl(eos_filename)) + +! read scalars + dims1(1)=1 + call h5dopen_f(file_id, "pointsrho", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_INTEGER, nrho, dims1, error) + call h5dclose_f(dset_id,error) + + if(error.ne.0) then + stop "Could not read EOS table file" + endif + + dims1(1)=1 + call h5dopen_f(file_id, "pointstemp", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_INTEGER, ntemp, dims1, error) + call h5dclose_f(dset_id,error) + + if(error.ne.0) then + stop "Could not read EOS table file" + endif + + dims1(1)=1 + call h5dopen_f(file_id, "pointsye", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_INTEGER, nye, dims1, error) + call h5dclose_f(dset_id,error) + + if(error.ne.0) then + stop "Could not read EOS table file" + endif + + write(message,"(a25,1P3i5)") "We have nrho ntemp nye: ", nrho,ntemp,nye + write(*,*) message + + allocate(alltables(nrho,ntemp,nye,nvars)) + + ! index variable mapping: + ! 1 -> logpress + ! 2 -> logenergy + ! 3 -> entropy + ! 4 -> munu + ! 5 -> cs2 + ! 6 -> dedT + ! 7 -> dpdrhoe + ! 8 -> dpderho + ! 9 -> muhat + ! 10 -> mu_e + ! 11 -> mu_p + ! 12 -> mu_n + ! 13 -> xa + ! 14 -> xh + ! 15 -> xn + ! 16 -> xp + ! 17 -> abar + ! 18 -> zbar + + + dims3(1)=nrho + dims3(2)=ntemp + dims3(3)=nye + call h5dopen_f(file_id, "logpress", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,1), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + call h5dopen_f(file_id, "logenergy", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,2), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + call h5dopen_f(file_id, "entropy", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,3), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + call h5dopen_f(file_id, "munu", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,4), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + call h5dopen_f(file_id, "cs2", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,5), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + call h5dopen_f(file_id, "dedt", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,6), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + call h5dopen_f(file_id, "dpdrhoe", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,7), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + call h5dopen_f(file_id, "dpderho", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,8), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + +! chemical potentials + call h5dopen_f(file_id, "muhat", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,9), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + call h5dopen_f(file_id, "mu_e", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,10), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + call h5dopen_f(file_id, "mu_p", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,11), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + call h5dopen_f(file_id, "mu_n", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,12), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + +! compositions + call h5dopen_f(file_id, "Xa", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,13), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + call h5dopen_f(file_id, "Xh", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,14), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + call h5dopen_f(file_id, "Xn", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,15), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + call h5dopen_f(file_id, "Xp", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,16), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + +! average nucleus + call h5dopen_f(file_id, "Abar", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,17), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + call h5dopen_f(file_id, "Zbar", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,18), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + +! Gamma + call h5dopen_f(file_id, "gamma", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,19), dims3, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + allocate(logrho(nrho)) + dims1(1)=nrho + call h5dopen_f(file_id, "logrho", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, logrho, dims1, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + allocate(logtemp(ntemp)) + dims1(1)=ntemp + call h5dopen_f(file_id, "logtemp", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, logtemp, dims1, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + allocate(ye(nye)) + dims1(1)=nye + call h5dopen_f(file_id, "ye", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, ye, dims1, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + + call h5dopen_f(file_id, "energy_shift", dset_id, error) + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, energy_shift, dims1, error) + call h5dclose_f(dset_id,error) + accerr=accerr+error + + if(accerr.ne.0) then + stop "Problem reading EOS table file" + endif + + + call h5fclose_f (file_id,error) + + call h5close_f (error) + + ! set min-max values: + + eos_rhomin = 10.0d0**logrho(1) + eos_rhomax = 10.0d0**logrho(nrho) + + eos_yemin = ye(1) + eos_yemax = ye(nye) + + eos_tempmin = 10.0d0**logtemp(1) + eos_tempmax = 10.0d0**logtemp(ntemp) + + write(6,*) "Done reading eos tables" + + +end subroutine readtable + +
Directory: /EOS_Omni/src/ =========================
File [modified]: make.code.defn Delta lines: +1 -1 =================================================================== --- EOS_Omni/src/make.code.defn 2010-08-26 19:26:54 UTC (rev 8) +++ EOS_Omni/src/make.code.defn 2010-10-22 19:42:05 UTC (rev 9) @@ -5,4 +5,4 @@ EOS_Omni_SingleVarCalls_harm.F90 EOS_Omni_Handles.c
# Subdirectories containing source files -SUBDIRS = +SUBDIRS = nuc_eos