[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 354)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Wed Jun 20 11:14:28 CDT 2012
User: rhaas
Date: 2012/06/20 11:14 AM
Modified:
/trunk/
param.ccl
/trunk/doc/
documentation.tex
/trunk/src/
GRHydro_Con2Prim.F90
Log:
GRHydro: C2P: resort to bisection if Newton-Raphson fails. This is only
activated if user requested this.
Patch by Christian Reisswig.
File Changes:
Directory: /trunk/doc/
======================
File [modified]: documentation.tex
Delta lines: +1 -1
===================================================================
--- trunk/doc/documentation.tex 2012-06-06 14:37:40 UTC (rev 353)
+++ trunk/doc/documentation.tex 2012-06-20 16:14:28 UTC (rev 354)
@@ -1991,7 +1991,7 @@
encouragement that it was possible!
Thirdly, the support of the Cactus team, especially Tom Goodale,
-Gabrielle Allen and Thomas Radke made (and, especially Thomas, continue to make) life a
+Gabrielle Allen and Thomas Radke made (and continues to make) life a
lot easier.
Finally, for their work in coding, ideas and suggestions, or just
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2Prim.F90
Delta lines: +50 -10
===================================================================
--- trunk/src/GRHydro_Con2Prim.F90 2012-06-06 14:37:40 UTC (rev 353)
+++ trunk/src/GRHydro_Con2Prim.F90 2012-06-20 16:14:28 UTC (rev 354)
@@ -319,8 +319,6 @@
endif
#endif
-
-
end do
end do
end do
@@ -368,7 +366,7 @@
! begin EOS Omni vars
integer :: n,keytemp,anyerr,keyerr(1)
- real*8 :: xpress,xtemp,xye
+ real*8 :: xpress,xtemp,xye,tmp,plow
n=1;keytemp=0;anyerr=0;keyerr(1)=0
xpress=0.0d0;xtemp=0.0d0;xye=0.0d0
! end EOS Omni vars
@@ -387,6 +385,10 @@
call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
pold = max(1.d-10,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(pmin, sqrt(s2) - utau - udens)
!!$ Check that the variables have a chance of being physical
@@ -398,10 +400,10 @@
else
!$OMP CRITICAL
!!!! call CCTK_WARN(GRHydro_NaN_verbose, "c2p failed and being told not to reset the pressure")
+ !call CCTK_WARN(0, "c2p failed and being told not to reset the pressure")
GRHydro_C2P_failed = 1
!$OMP END CRITICAL
endif
-
endif
!!$ Calculate rho and epsilon
@@ -477,14 +479,52 @@
pold = pnew
- pnew = max(pold - f/df, pmin)
+ ! 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(pmin, sqrt(s2) - utau - udens)
+
+ if ((pnew .lt. plow .or. tmp .le. 0.0d0) .and. c2p_resort_to_bisection) 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(pmin, 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)
- w_lorentz = (utau + pnew + udens) / sqrt( (utau + pnew + udens)**2 - &
- s2)
- epsilon = (sqrt( (utau + pnew + udens)**2 - s2) - pnew * w_lorentz - &
+
+ rho = udens * sqrt( tmp)/(utau + pnew + udens)
+
+ w_lorentz = (utau + pnew + udens) / sqrt( tmp)
+ epsilon = (sqrt( tmp) - pnew * w_lorentz - &
udens)/udens
call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +4 -0
===================================================================
--- trunk/param.ccl 2012-06-06 14:37:40 UTC (rev 353)
+++ trunk/param.ccl 2012-06-20 16:14:28 UTC (rev 354)
@@ -314,6 +314,10 @@
0: :: "greater than zero"
} 1.e-20
+boolean c2p_resort_to_bisection "If the pressure guess is unphysical, should we try with bisection (slower, but more robust)"
+{
+} no
+
real GRHydro_eps_min "Minimum value of specific internal energy - this is now only used in GRHydro_InitData's GRHydro_Only_Atmo routine"
{
0: :: "Positive"
More information about the Commits
mailing list