diff options
Diffstat (limited to 'src/nuc_eos/findtemp.F90')
| -rw-r--r-- | src/nuc_eos/findtemp.F90 | 107 |
1 files changed, 90 insertions, 17 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 |