[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