[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 213)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Sat Feb 5 10:59:38 CST 2011
User: cott
Date: 2011/02/05 10:59 AM
Modified:
/trunk/src/
GRHydro_Con2Prim.F90
Log:
* update handling of the hot c2p business: add
fallback option when c2p does not converge -- in this case
the pointwise c2p is called again with a lower tolerance.
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2Prim.F90
Delta lines: +43 -45
===================================================================
--- trunk/src/GRHydro_Con2Prim.F90 2011-02-02 19:44:05 UTC (rev 212)
+++ trunk/src/GRHydro_Con2Prim.F90 2011-02-05 16:59:37 UTC (rev 213)
@@ -58,6 +58,7 @@
CCTK_INT :: type2_bits
CCTK_REAL :: local_min_tracer
+ CCTK_REAL :: local_perc_ptol
! begin EOS Omni vars
integer :: n,keytemp,anyerr,keyerr(1)
@@ -157,12 +158,40 @@
z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, &
GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
else
- call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2), &
- scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), &
- vel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
+ call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),&
+ scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),&
+ vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
+ press(i,j,k),w_lorentz(i,j,k), &
uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, &
- GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
+ GRHydro_reflevel, GRHydro_C2P_failed(i,j,k), GRHydro_perc_ptol)
+ if(GRHydro_C2P_failed(i,j,k).eq.2) then
+ ! this means c2p did not converge.
+ ! In this case, we attempt to call c2p with a reduced
+ ! accuracy requirement; if it fails again, we abort
+ GRHydro_C2P_failed(i,j,k) = 0
+ local_perc_ptol = GRHydro_perc_ptol*100.0d0
+ call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),&
+ scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),rho(i,j,k),&
+ vel(i,j,k,1),vel(i,j,k,2), vel(i,j,k,3),eps(i,j,k),&
+ temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
+ uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
+ z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, &
+ GRHydro_reflevel, GRHydro_C2P_failed(i,j,k), local_perc_ptol)
+ if(GRHydro_C2P_failed(i,j,k).eq.2) then
+ if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
+ call CCTK_WARN(1,"Convergence problem in c2p")
+ 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),&
+ temperature(i,j,k),y_e(i,j,k)
+ call CCTK_WARN(1,warnline)
+ call CCTK_WARN(0,"Aborting!!!")
+ endif
+ endif
+ endif
endif
if (epsnegative) then
@@ -272,7 +301,7 @@
n=1;keytemp=0;anyerr=0;keyerr(1)=0
xpress=0.0d0;xtemp=0.0d0;xye=0.0d0
! end EOS Omni vars
-
+
!!$ Undensitize the variables
udens = dens /sqrt(det)
@@ -466,7 +495,7 @@
subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, vely, &
velz, epsilon, temp, ye, press, w_lorentz, uxx, uxy, uxz, uyy, &
uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, &
- GRHydro_reflevel, GRHydro_C2P_failed)
+ GRHydro_reflevel, GRHydro_C2P_failed, local_perc_ptol)
implicit none
@@ -482,6 +511,7 @@
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 local_perc_ptol
character(len=200) warnline
logical epsnegative
@@ -495,11 +525,10 @@
n=1;keytemp=0;anyerr=0;keyerr(1)=0
temp0 = 0.0d0;xpress = 0.0d0
! end EOS Omni vars
-
+
failinfomode = 1
failwarnmode = 0
-
! set pmin and epsmin to something sensible:
pminl = 1.0d-28
epsminl = 1.0e-5
@@ -543,7 +572,6 @@
endif
endif
-
!!$ Check that the variables have a chance of being physical
if( (utau + pold + udens)**2 - s2 .le. 0.0d0) then
@@ -594,46 +622,16 @@
count = 0
pnew = pold
- do while ( ((abs(pnew - pold)/abs(pnew) .gt. GRHydro_perc_ptol) .and. &
+ do while ( ((abs(pnew - pold)/abs(pnew) .gt. local_perc_ptol) .and. &
(abs(pnew - pold) .gt. GRHydro_del_ptol)) .or. &
(count .lt. GRHydro_countmin))
count = count + 1
if (count > GRHydro_countmax) then
- GRHydro_C2P_failed = 1
-
- !$OMP CRITICAL
- if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
- call CCTK_WARN(failinfomode, 'count > GRHydro_countmax! ')
- write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
- call CCTK_WARN(failinfomode,warnline)
- write(warnline,'(a20,3g16.7)') 'xyz location: ',x,y,z
- call CCTK_WARN(failinfomode,warnline)
- write(warnline,'(a20,3i5)') 'ijk location: ',ii,jj,kk
- call CCTK_WARN(failinfomode,warnline)
- write(warnline,'(a20,g16.7)') 'radius: ',r
- call CCTK_WARN(failinfomode,warnline)
- call CCTK_WARN(failinfomode,"Setting the point to atmosphere")
- endif
- !$OMP END CRITICAL
-
- ! for safety, let's set the point to atmosphere
- rho = GRHydro_rho_min
- udens = rho
- dens = sqrt(det) * rho
- pnew = pminl
- epsilon = epsminl
- ! w_lorentz=1, so the expression for utau reduces to:
- utau = rho + rho*epsminl - udens
- sx = 0.d0
- sy = 0.d0
- sz = 0.d0
- s2 = 0.d0
- usx = 0.d0
- usy = 0.d0
- usz = 0.d0
- w_lorentz = 1.d0
- goto 51
+ ! non-convergence is now handled outside of this
+ ! routine
+ GRHydro_C2P_failed = 2
+ return
end if
call EOS_Omni_DPressByDRho(handle,keytemp,GRHydro_eos_rf_prec,n,&
@@ -679,6 +677,7 @@
call CCTK_WARN(failwarnmode,"Aborting!!!")
endif
endif
+
f = pnew - xpress
enddo
@@ -713,7 +712,6 @@
call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, &
rho,epsilon,temp,ye,xpress,keyerr,anyerr)
-
! error handling
if(anyerr.ne.0) then
if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
More information about the Commits
mailing list