diff options
| -rw-r--r-- | doc/documentation.tex | 2 | ||||
| -rw-r--r-- | param.ccl | 4 | ||||
| -rw-r--r-- | src/GRHydro_Con2Prim.F90 | 60 |
3 files changed, 55 insertions, 11 deletions
diff --git a/doc/documentation.tex b/doc/documentation.tex index b6f3ff9..1a8d6c2 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -1991,7 +1991,7 @@ GRAstro\_Hydro} codes gave an essential benchmark to aim for, and 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 @@ -314,6 +314,10 @@ real c2p_reset_pressure_to_value "The value to which the pressure is reset to wh 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" diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index f61eea7..0916852 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -319,8 +319,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) endif #endif - - end do end do end do @@ -368,7 +366,7 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & ! 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 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & 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,& |