[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 355)

[rhaas at tapir.caltech.edu]

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)