* 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.



git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@213 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 43 insertions, 45 deletions

diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index 0158acf..fe57ca4 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -58,6 +58,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
   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 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
                  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 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
   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 @@ end subroutine Con2Prim_pt
 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 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
   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 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
   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 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
      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 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
   
   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 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
           call CCTK_WARN(failwarnmode,"Aborting!!!")
        endif
     endif
+
     f = pnew - xpress
 
   enddo
@@ -713,7 +712,6 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
     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