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)