[Commits] [svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos/ (Rev. 89)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Wed Mar 12 22:01:00 CDT 2014
User: rhaas
Date: 2014/03/12 10:01 PM
Modified:
/trunk/src/nuc_eos/
eosmodule.F90, findtemp.F90, linterp.F, linterp_many.F90, nuc_eos.F90
Log:
* Major bug fix in EOS_Omni/nuc_eos. In the code (before the fix),
offending specific internal energies where limited and the EOS
returned seemingly fine, but thermodynamically inconsistent data.
This has been fixed. This fix also makes isolated microphysical NS
simulations not work anymore. Why this is the case, we are still
investigating
From: Christian Ott <cott at tapir.caltech.edu>
File Changes:
Directory: /trunk/src/nuc_eos/
==============================
File [modified]: eosmodule.F90
Delta lines: +4 -0
===================================================================
--- trunk/src/nuc_eos/eosmodule.F90 2014-03-13 03:00:46 UTC (rev 88)
+++ trunk/src/nuc_eos/eosmodule.F90 2014-03-13 03:01:00 UTC (rev 89)
@@ -25,6 +25,7 @@
! basics
integer, parameter :: nvars = 19
real*8,pointer :: alltables(:,:,:,:)
+ real*8,allocatable :: epstable(:,:,:)
! index variable mapping:
! 1 -> logpress
! 2 -> logenergy
@@ -87,6 +88,8 @@
logtemp => logtemp_
ye => ye_
+ allocate(epstable(nrho,ntemp,nye))
+ epstable(1:nrho,1:ntemp,1:nye) = 10.0d0**alltables(1:nrho,1:ntemp,1:nye,2)
energy_shift = energy_shift_
if(do_energy_shift.ne.1) then
@@ -103,5 +106,6 @@
eos_tempmin = 10.0d0**logtemp(1)
eos_tempmax = 10.0d0**logtemp(ntemp)
+
end subroutine setup_eosmodule
File [modified]: findtemp.F90
Delta lines: +126 -65
===================================================================
--- trunk/src/nuc_eos/findtemp.F90 2014-03-13 03:00:46 UTC (rev 88)
+++ trunk/src/nuc_eos/findtemp.F90 2014-03-13 03:01:00 UTC (rev 89)
@@ -1,3 +1,107 @@
+subroutine findtemp_deb(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 ! precision to which we need to find temp
+ 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
+
+ do i=1,itmax
+ !d2 is the derivative deps/dlogtemp;
+ ldt = -(eps - eps0)/d2
+ ltn = lt+ldt
+ ltn = min(ltn,ltmax)
+ ltn = max(ltn,ltmin)
+ lt1=lt
+ lt=ltn
+ eps1=eps
+ call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3)
+ if (abs(eps - eps0).lt.tol*abs(eps0)) 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(eps-eps0).lt.1.0d-3*abs(eps0)) then
+ d2 = (eps-eps1)/(lt-lt1)
+ 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
+#if 0
+ ! 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
+ endif
+
+ lt0=min(lt,log10(t_max_hack))
+ return
+ endif
+
+12 continue
+
+ lt0=min(lt,log10(t_max_hack))
+
+
+end subroutine findtemp_deb
+
subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps)
use eosmodule
@@ -206,77 +310,34 @@
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 lr,t0,y,epsin,rfeps,dlepsdt
+ integer :: keyerrt
- real*8 tn,tmin,tmax
- real*8 tinput,rfeps
- integer :: rl = 0
- integer itmax,i,keyerrt
- integer ii,jj,kk
-
- keyerrt=0
+ real*8 :: eps2, eps1
+ real*8 :: t2, t1
+ real*8 :: m
- tol=rfeps ! need to find energy to less than 1 in 10^-10
- itmax=100 ! use at most 20 iterations, then bomb
+ keyerrt = 0
+ t2 = 10.0d0**logtemp(2)
+ t1 = 10.0d0**logtemp(1)
- t=t0
- t1=t
+ call intp3d_linearTlow(lr,t2,y,eps2,1,epstable,nrho,&
+ ntemp,nye,logrho,logtemp,ye,dlepsdt)
- eps0=epsin
- eps1=eps0
-
- tmax=10.0d0**logtemp(2)
- tmin=-20.0d0
+ call intp3d_linearTlow(lr,t1,y,eps1,1,epstable,nrho,&
+ ntemp,nye,logrho,logtemp,ye,dlepsdt)
- !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)
+ m = (t2-t1)/(eps2-eps1)
+ t0 = (epsin-eps1) * m + t1
- 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
+#if 0
+ !$OMP CRITICAL
+ write(6,*) t0,epsin,eps2,eps1
+ call intp3d_linearTlow(lr,t0,y,eps2,1,epstable,nrho,&
+ ntemp,nye,logrho,logtemp,ye,dlepsdt)
+ write(6,"(1P10E15.6)") abs(epsin-eps2)/epsin
+ !$OMP END CRITICAL
+#endif
- 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.F
Delta lines: +3 -3
===================================================================
--- trunk/src/nuc_eos/linterp.F 2014-03-13 03:00:46 UTC (rev 88)
+++ trunk/src/nuc_eos/linterp.F 2014-03-13 03:01:00 UTC (rev 89)
@@ -61,9 +61,9 @@
dy = (yt(ny) - yt(1)) / FLOAT(ny-1)
dz = (zt(nz) - zt(1)) / FLOAT(nz-1)
c
- dxi = 1. / dx
- dyi = 1. / dy
- dzi = 1. / dz
+ dxi = 1.0d0 / dx
+ dyi = 1.0d0 / dy
+ dzi = 1.0d0 / dz
c
dxyi = dxi * dyi
dxzi = dxi * dzi
File [modified]: linterp_many.F90
Delta lines: +9 -9
===================================================================
--- trunk/src/nuc_eos/linterp_many.F90 2014-03-13 03:00:46 UTC (rev 88)
+++ trunk/src/nuc_eos/linterp_many.F90 2014-03-13 03:01:00 UTC (rev 89)
@@ -54,9 +54,9 @@
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
+ dxi = 1.0d0 / dx
+ dyi = 1.0d0 / dy
+ dzi = 1.0d0 / dz
!
dxyi = dxi * dyi
dxzi = dxi * dzi
@@ -181,9 +181,9 @@
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
+ dxi = 1.0d0 / dx
+ dyi = 1.0d0 / dy
+ dzi = 1.0d0 / dz
!
dxyi = dxi * dyi
dxzi = dxi * dzi
@@ -304,9 +304,9 @@
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
+ dxi = 1.0d0 / dx
+ dyi = 1.0d0 / dy
+ dzi = 1.0d0 / dz
!
dxyi = dxi * dyi
dxzi = dxi * dzi
File [modified]: nuc_eos.F90
Delta lines: +48 -15
===================================================================
--- trunk/src/nuc_eos/nuc_eos.F90 2014-03-13 03:00:46 UTC (rev 88)
+++ trunk/src/nuc_eos/nuc_eos.F90 2014-03-13 03:01:00 UTC (rev 89)
@@ -83,13 +83,17 @@
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
- call findtemp(lr,lt,y,leps,keyerrt,rfeps)
+ if(xeps .gt. 0.0d0) then
+ leps = log10(xeps)
+ call findtemp(lr,lt,y,leps,keyerrt,rfeps)
+ else
+ keyerrt = 667
+ endif
if(keyerrt.ne.0) then
call CCTK_WARN (0, "Did not find temperature")
endif
@@ -272,18 +276,23 @@
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.eq.667 .and. lt .lt. 10.0d0) then
+ if(xeps .gt. 0.0d0) then
+ leps = log10(xeps)
+ call findtemp(lr,lt,y,leps,keyerrt,rfeps)
+ else
+ keyerrt = 667
+ endif
+ if(keyerrt.eq.667) 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)
+ call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps)
if(keyerrt.eq.0) then
keyerrt = 668
xtemp = xt
else
keyerrt = 667
+ keyerr = keyerrt
return
endif
keyerr = keyerrt
@@ -365,6 +374,7 @@
real*8 :: lr,lt,y,xt,xx,xeps,leps,xs
real*8 :: d1,d2,d3,ff(8)
integer :: keyerrt
+ character(len=256) :: warnline
keyerrt = 0
@@ -411,20 +421,38 @@
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.eq.667 .and. lt .lt. 10.0d0) then
+ if(xeps .gt. 0.0d0) then
+ leps = log10(xeps)
+ call findtemp(lr,lt,y,leps,keyerrt,rfeps)
+ else
+ keyerrt = 667
+ endif
+! if(xenr .lt. -9.00725347549139d20) then
+! if(xenr.lt.-2.28659899d20) then
+! write(warnline,"(A4,i5,1P10E15.6)") "UGGx",keyerrt,lr,lt,xt,xtemp,y,xeps,xenr
+! call CCTK_WARN(1,warnline)
+! endif
+ if(keyerrt.eq.667) 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)
+ call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps)
if(keyerrt.eq.0) then
keyerrt = 668
+ if(xeps.lt.0.0d0) then
+ keyerrt = 668 ! impossible to do anything
+ endif
xtemp = xt
else
keyerrt = 667
+ keyerr = keyerrt
return
endif
+! if(xenr .lt. -9.00725347549139d20) then
+! if(xenr.lt.-2.28659899d20) then
+! write(warnline,"(A4,i5,1P10E15.6)") "UGG",keyerrt,lr,lt,xt,xtemp,y,xeps,xenr
+! call CCTK_WARN(0,warnline)
+! endif
keyerr = keyerrt
else
xtemp = 10.0d0**lt
@@ -502,20 +530,25 @@
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.eq.667 .and. lt .lt. 10.0d0) then
+ if(xeps .gt. 0.0d0) then
+ leps = log10(xeps)
+ call findtemp(lr,lt,y,leps,keyerrt,rfeps)
+ else
+ keyerrt = 667
+ endif
+ if(keyerrt.eq.667) 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)
+ call findtemp_low(lr,xt,y,xeps,keyerrt,rfeps)
if(keyerrt.eq.0) then
keyerrt = 668
xtemp = xt
else
keyerrt = 667
+ keyerr = keyerrt
return
- endif
+ endif
keyerr = keyerrt
else
xtemp = 10.0d0**lt
More information about the Commits
mailing list