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(a)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
