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"