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
