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

rhaas at tapir.caltech.edu rhaas at tapir.caltech.edu
Tue Jul 17 12:08:33 CDT 2012


User: rhaas
Date: 2012/07/17 12:08 PM

Modified:
 /trunk/src/
  GRHydro_Con2Prim.F90

Log:
 GRHydro: Con2Prim improvement to bisection algorithm. This fixes some NANs associated with NS collapse test case when full GR is used.
 
 From: Christian Reisswig <reisswig at scriwalker.(none)>

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_Con2Prim.F90
Delta lines: +175 -8
===================================================================
--- trunk/src/GRHydro_Con2Prim.F90	2012-07-17 17:08:31 UTC (rev 400)
+++ trunk/src/GRHydro_Con2Prim.F90	2012-07-17 17:08:33 UTC (rev 401)
@@ -118,6 +118,21 @@
     do j = 1, ny 
       do i = 1, nx
 
+#if 0
+        !$OMP CRITICAL
+         if (rho(i,j,k).le.0.0d0) then
+	    call CCTK_WARN(1,"rho less than zero at entry of Con2Prim")
+	    write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel
+	    call CCTK_WARN(1,warnline)
+	    write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k)
+	    call CCTK_WARN(1,warnline)
+	    write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k)
+	    call CCTK_WARN(1,warnline)
+	    call CCTK_WARN(0,"Aborting!!!")
+	 endif
+	 !$OMP END CRITICAL
+#endif
+
         !do not compute if in atmosphere or in excised region
         if ((atmosphere_mask(i,j,k) .ne. 0) .or. &
             (hydro_excision_mask(i,j,k) .ne. 0)) cycle
@@ -320,6 +335,46 @@
         endif
 #endif
 
+!! Some debugging convenience output
+# if 0
+	 !$OMP CRITICAL
+	 if (dens(i,j,k).ne.dens(i,j,k)) then
+	    call CCTK_WARN(1,"NAN at exit of Con2Prim")
+	    write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel
+	    call CCTK_WARN(1,warnline)
+	    write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k)
+	    call CCTK_WARN(1,warnline)
+	    write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k)
+	    call CCTK_WARN(1,warnline)
+	    call CCTK_WARN(0,"Aborting!!!")
+	 endif
+	 if (rho(i,j,k).ne.rho(i,j,k)) then
+	    call CCTK_WARN(1,"NAN in rho at exit of Con2Prim")
+	    write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel
+	    call CCTK_WARN(1,warnline)
+	    write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k)
+	    call CCTK_WARN(1,warnline)
+	    write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k)
+	    call CCTK_WARN(1,warnline)
+	    call CCTK_WARN(0,"Aborting!!!")
+	 endif
+	 !$OMP END CRITICAL
+	 
+	 !$OMP CRITICAL
+         if (rho(i,j,k).le.0.0d0) then
+	    call CCTK_WARN(1,"rho less than zero at exit of Con2Prim")
+	    write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel
+	    call CCTK_WARN(1,warnline)
+	    write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k)
+	    call CCTK_WARN(1,warnline)
+	    write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k)
+	    call CCTK_WARN(1,warnline)
+	    call CCTK_WARN(0,"Aborting!!!")
+	 endif
+	 !$OMP END CRITICAL
+#endif
+         
+
       end do
     end do
   end do
@@ -363,7 +418,7 @@
   CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, &
        w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin
   character(len=200) warnline
-  logical epsnegative
+  logical epsnegative, mustbisect
 
 ! begin EOS Omni vars
   integer :: n,keytemp,anyerr,keyerr(1)
@@ -385,12 +440,16 @@
 !!$  Set initial guess for pressure:
   call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
                       rho,epsilon,xtemp,xye,xpress,keyerr,anyerr)
-  pold = max(1.d-10,xpress)
+  !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)
+  
+  ! Start out with some reasonable initial guess
+  pold = max(plow+1.d-10,xpress)
 
+  mustbisect = .false.
 
 !!$  Check that the variables have a chance of being physical
 
@@ -423,6 +482,11 @@
 
   f = pold - xpress
 
+  if (f .ne. f) then 
+     ! Ok, this yielded nonsense, let's enforce bisection!
+     mustbisect = .true.
+  endif
+
 !!$Find the root
   
   count = 0
@@ -492,14 +556,15 @@
       tmp = (utau + pnew + udens)**2 - s2
       plow = max(pmin, sqrt(s2) - utau - udens)
       
-      if (pnew .lt. plow .or. tmp .le. 0.0d0) then
+      if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect .eq. .true.) 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
-	       
+	 
+	 mustbisect = .false.
       end if
     
     else
@@ -524,6 +589,48 @@
         
     rho = udens * sqrt( tmp)/(utau + pnew + udens)
 
+!! Some debugging convenience output
+#if 0
+    if (rho .ne. rho) then
+    !$OMP CRITICAL
+      write(warnline,'(a38, i16)') 'NAN in rho! Root finding iteration: ', count
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'radius: ',r
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'udens: ', udens
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'dens: ', dens
+      call CCTK_WARN(1,warnline)
+      !write(warnline,'(a20,g16.7)') 'rho-old: ', rhoold
+      !call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'pnew: ', pnew
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'pold: ', pold
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'xpress: ', xpress
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'f: ', f
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'df: ', df
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'f/df: ', f/df
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a30,g16.7)') '(utau + pnew + udens)**2 - s2: ', (utau + pnew + udens)**2 - s2
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a30,g16.7)') '(utau + pnew + udens): ', (utau + pnew + udens)
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'drhobydpress: ', drhobydpress
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'depsbydpress: ', depsbydpress
+      call CCTK_WARN(0,warnline)
+   !$OMP END CRITICAL
+  endif
+#endif
+
     w_lorentz = (utau + pnew + udens) / sqrt( tmp)
     epsilon = (sqrt( tmp) - pnew * w_lorentz - &
          udens)/udens
@@ -533,6 +640,11 @@
     
     f = pnew - xpress
 
+    if (f .ne. f) then 
+       ! Ok, this yielded nonsense, let's enforce bisection!
+       mustbisect = .true.
+    endif
+
   enddo
   
 !!$  Polish the root
@@ -604,6 +716,47 @@
     epsnegative = .true.
   endif
 
+#if 0
+  if (rho .le. 0.0d0) then
+    !$OMP CRITICAL
+      write(warnline,'(a38, i16)') 'Epsilon negative!'
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'radius: ',r
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'udens: ', udens
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'dens: ', dens
+      call CCTK_WARN(1,warnline)
+      !write(warnline,'(a20,g16.7)') 'rho-old: ', rhoold
+      !call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'pnew: ', pnew
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'pold: ', pold
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'xpress: ', xpress
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'f: ', f
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'df: ', df
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'f/df: ', f/df
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a30,g16.7)') '(utau + pnew + udens)**2 - s2: ', (utau + pnew + udens)**2 - s2
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a30,g16.7)') '(utau + pnew + udens): ', (utau + pnew + udens)
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'drhobydpress: ', drhobydpress
+      call CCTK_WARN(1,warnline)
+      write(warnline,'(a20,g16.7)') 'depsbydpress: ', depsbydpress
+      call CCTK_WARN(0,warnline)
+   !$OMP END CRITICAL
+  endif
+#endif
+
 end subroutine Con2Prim_pt
 
 subroutine Con2Prim_pt_hot(cctk_iteration, ii,jj,kk,handle, dens, &
@@ -629,7 +782,7 @@
   CCTK_REAL local_perc_ptol
 
   character(len=256) warnline
-  logical epsnegative
+  logical epsnegative, mustbisect
 
   integer :: failwarnmode 
   integer :: failinfomode 
@@ -643,6 +796,8 @@
   nf=0;nfudgemax=30
 ! end EOS Omni vars
 
+  mustbisect = .false.
+
   failinfomode = 1
   failwarnmode = 0
 
@@ -672,11 +827,13 @@
   temp0 = temp
   call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
                       rho,epsilon,temp,ye,xpress,keyerr,anyerr)
-  pold = max(1.d-15,xpress)
+  !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)
+  ! Start out with some reasonable initial guess
+  pold = max(plow+1.d-10,xpress)
   ! error handling
   if(anyerr.ne.0) then
      !$OMP CRITICAL
@@ -836,6 +993,10 @@
   endif
   f = pold - xpress
 
+  if (f .ne. f) then 
+     ! Ok, this yielded nonsense, let's enforce bisection!
+     mustbisect = .true.
+  endif
 
 !!$Find the root
   
@@ -878,14 +1039,15 @@
       tmp = (utau + pnew + udens)**2 - s2
       plow = max(pminl, sqrt(s2) - utau - udens)
       
-      if (pnew .lt. plow .or. tmp .le. 0.0d0) then
+      if (pnew .lt. plow .or. tmp .le. 0.0d0 .or. mustbisect .eq. .true.) 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
-	       
+	 
+	 mustbisect = .false.
       end if
     
     else
@@ -985,6 +1147,11 @@
 
     f = pnew - xpress
 
+    if (f .ne. f) then 
+      ! Ok, this yielded nonsense, let's enforce bisection!
+      mustbisect = .true.
+    endif
+
   enddo
 
   



More information about the Commits mailing list