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

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@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