[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 355)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Wed Jun 20 11:15:16 CDT 2012
User: rhaas
Date: 2012/06/20 11:15 AM
Modified:
/trunk/src/
GRHydro_Con2Prim.F90
Log:
C2P_hot: Resort to bisection as well (if user requested so).
Patch by Christian Reisswig.
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2Prim.F90
Delta lines: +44 -3
===================================================================
--- trunk/src/GRHydro_Con2Prim.F90 2012-06-20 16:14:28 UTC (rev 354)
+++ trunk/src/GRHydro_Con2Prim.F90 2012-06-20 16:15:15 UTC (rev 355)
@@ -491,7 +491,7 @@
tmp = (utau + pnew + udens)**2 - s2
plow = max(pmin, sqrt(s2) - utau - udens)
- if ((pnew .lt. plow .or. tmp .le. 0.0d0) .and. c2p_resort_to_bisection) then
+ if (pnew .lt. plow .or. tmp .le. 0.0d0) then
! Ok, Newton-Raphson ended up finding something unphysical.
! Let's try to find our root via bisection (which converges slower but is more robust)
@@ -624,7 +624,7 @@
CCTK_REAL GRHydro_C2P_failed
CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, &
w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin
- CCTK_REAL epsminl,pminl
+ CCTK_REAL epsminl,pminl,plow,tmp
CCTK_REAL local_perc_ptol
character(len=256) warnline
@@ -672,6 +672,10 @@
call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
rho,epsilon,temp,ye,xpress,keyerr,anyerr)
pold = max(1.d-15,xpress)
+ ! This is the lowest admissible pressure, otherwise we are off the physical regime
+ ! Sometimes, we may end up with bad initial guesses. In that case we need to start off with something
+ ! reasonable at least
+ plow = max(pminl, sqrt(s2) - utau - udens)
! error handling
if(anyerr.ne.0) then
!$OMP CRITICAL
@@ -861,8 +865,45 @@
dpressbydeps*depsbydpress
pold = pnew
- pnew = max(pold - f/df, pminl)
+ ! Try to obtain new pressure via Newton-Raphson.
+
+ pnew = pold - f/df
+
+ ! Check if Newton-Raphson resulted in something reasonable!
+
+ if (c2p_resort_to_bisection) then
+
+ tmp = (utau + pnew + udens)**2 - s2
+ plow = max(pminl, sqrt(s2) - utau - udens)
+
+ if (pnew .lt. plow .or. tmp .le. 0.0d0) then
+
+ ! Ok, Newton-Raphson ended up finding something unphysical.
+ ! Let's try to find our root via bisection (which converges slower but is more robust)
+
+ pnew = (plow + pold) / 2
+ tmp = (utau + pnew + udens)**2 - s2
+
+ end if
+
+ else
+
+ ! This is the standard algorithm without resorting to bisection.
+
+ pnew = max(pminl, pnew)
+ tmp = (utau + pnew + udens)**2 - s2
+
+ endif
+
+ ! Check if we are still in the physical range now.
+ ! If not, we set the C2P failed mask, and set pnew = pold (which will exit the loop).
+
+ if ((tmp .le. 0.0d0) .and. GRHydro_C2P_failed .eq. 0) then
+ GRHydro_C2P_failed = 1
+ pnew = pold
+ endif
+
!!$ Recalculate primitive variables and function
rho = udens * sqrt( (utau + pnew + udens)**2 - s2)/(utau + pnew + udens)
More information about the Commits
mailing list