add more workarounds in Con2Prim for hot EOS

* fix con2prim issues (pertaining to error messages)
* fix OMP CRITICAL issues in Prim2Con and Con2Prim

original commits by Christian Ott (cott)



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

2 files changed, 94 insertions, 51 deletions

diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index 2ba43cd..8279057 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -52,7 +52,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
   integer :: i, j, k, itracer, nx, ny, nz
   CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin
   logical :: epsnegative
-  character(len=100) warnline
+  character*256 :: warnline
   
   CCTK_INT :: type_bits, atmosphere
   CCTK_INT :: type2_bits
@@ -88,7 +88,8 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
 
   !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,&
-  !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp)
+  !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp,&
+  !$OMP warnline)
   do k = 1, nz 
     do j = 1, ny 
       do i = 1, nx
@@ -147,6 +148,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
                    press(i,j,k),keyerr,anyerr)
            else
               ! use polytropic EOS
+              keytemp = 0
               call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                    rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
 
@@ -189,7 +191,7 @@ 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), local_perc_ptol)
                if(GRHydro_C2P_failed(i,j,k).eq.2) then
-                  !$OMP CRITICAL(c2po1)
+                  !$OMP CRITICAL
                   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
@@ -201,7 +203,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
                      call CCTK_WARN(1,warnline)
                      call CCTK_WARN(0,"Aborting!!!")
                   endif
-                  !$OMP END CRITICAL(c2po1)
+                  !$OMP END CRITICAL
                endif
             endif
          endif
@@ -263,6 +265,38 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
 #endif          
 
         end if
+
+
+#if 1
+        ! this is needed in the current implementation of refluxing -- in this implementation,
+        ! the entire correction is applied to the coarse grid cell.
+        ! This leads to the cell getting too cold.
+        ! We work around this issue by enforcing that the temperature be
+        ! above 1.2 MeV at a density above 3.0e11. If this is the case, we inject heat.
+        if(evolve_temper.ne.0) then
+           if(GRHydro_eos_hot_eps_fix.ne.0.and.&
+                rho(i,j,k).gt.5.0d-7 .and. temperature(i,j,k).lt.1.2d0) then
+              !$OMP CRITICAL
+!              write(warnline,"(A40)") "Fixing temperature at:"
+!              call CCTK_WARN(1,warnline)
+              write(warnline,"(A10,4i6,1P3E15.6)") "fixtemp:", GRHydro_reflevel,i,j,k,&
+                   x(i,j,k),y(i,j,k),z(i,j,k)
+              call CCTK_WARN(1,warnline)
+              !$OMP END CRITICAL
+              keytemp = 1
+              temperature(i,j,k) = 1.5d0
+              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),press(i,j,k),keyerr,anyerr)
+              tau(i,j,k)  = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k) +press(i,j,k))&
+                   * w_lorentz(i,j,k)**2 - dens(i,j,k) - press(i,j,k)
+              keytemp = 0
+
+           endif
+        endif
+#endif
+
+
+
       end do
     end do
   end do
@@ -526,7 +560,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
   CCTK_REAL epsminl,pminl
   CCTK_REAL local_perc_ptol
 
-  character(len=200) warnline
+  character(len=256) warnline
   logical epsnegative
 
   integer :: failwarnmode 
@@ -538,7 +572,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
   integer :: nfudgemax,nf
   n=1;keytemp=0;anyerr=0;keyerr(1)=0
   temp0 = 0.0d0;xpress = 0.0d0
-  nf=0;nfudgemax=15
+  nf=0;nfudgemax=30
 ! end EOS Omni vars
 
   failinfomode = 1
@@ -553,7 +587,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
         y .lt. -1.0d-12 .or.&
         z .lt. -1.0d-12)) then
      failwarnmode = 2
-     failinfomode = 2
+     failinfomode = 2  
   endif
 
 !!$  Undensitize the variables 
@@ -573,17 +607,8 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
   pold = max(1.d-15,xpress)
   ! error handling
   if(anyerr.ne.0) then
-     !$OMP CRITICAL(c2p0)
+     !$OMP CRITICAL
      if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
-        call CCTK_WARN(failinfomode,"EOS error in c2p 0")
-        write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
-        call CCTK_WARN(failinfomode,warnline)
-        write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye
-        call CCTK_WARN(failinfomode,warnline)
-        write(warnline,"(A7,i8)") "code: ",keyerr(1)
-        call CCTK_WARN(failinfomode,warnline)
-        write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
-        call CCTK_WARN(failinfomode,warnline)
         if(keyerr(1).eq.667.and.GRHydro_eos_hot_eps_fix.ne.0) then
            ! Handling of the case when no new temperature can be
            ! found for a given epsilon. The amount of times
@@ -592,6 +617,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
            nf = 0
            do while(anyerr.ne.0.and.nf.le.nfudgemax)
               anyerr = 0
+              utau = utau + rho*abs(epsilon)*0.05d0*w_lorentz**2
               epsilon = epsilon + abs(epsilon)*0.05d0
               call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
                    rho,epsilon,temp,ye,xpress,keyerr,anyerr)
@@ -602,7 +628,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
            write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
            call CCTK_WARN(failinfomode,warnline)
            if(nf.gt.nfudgemax) then
-              call CCTK_WARN(failinfomode,"EOS error in c2p 1: injected heat too many times")
+              call CCTK_WARN(failinfomode,"EOS error in c2p 0a: injected heat too many times")
               write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
               call CCTK_WARN(failinfomode,warnline)
               write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye
@@ -614,12 +640,19 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
               call CCTK_WARN(failwarnmode,"Aborting!!!")
            endif
         else
+           call CCTK_WARN(failinfomode,"EOS error in c2p 0")
+           write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+           call CCTK_WARN(failinfomode,warnline)
+           write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye
+           call CCTK_WARN(failinfomode,warnline)
+           write(warnline,"(A7,i8)") "code: ",keyerr(1)
+           call CCTK_WARN(failinfomode,warnline)
+           write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
            call CCTK_WARN(failinfomode,warnline)
            call CCTK_WARN(failwarnmode,"Aborting!!!")
-
         endif
      endif
-     !$OMP END CRITICAL(c2p0)
+     !$OMP END CRITICAL
   endif
 
 !!$  Check that the variables have a chance of being physical
@@ -652,17 +685,8 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
                       rho,epsilon,temp,ye,xpress,keyerr,anyerr)
   ! error handling
   if(anyerr.ne.0) then
-     !$OMP CRITICAL(c2p1)
+     !$OMP CRITICAL
      if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
-        call CCTK_WARN(failinfomode,"EOS error in c2p 1")
-        write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
-        call CCTK_WARN(failinfomode,warnline)
-        write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye
-        call CCTK_WARN(failinfomode,warnline)
-        write(warnline,"(A7,i8)") "code: ",keyerr(1)
-        call CCTK_WARN(failinfomode,warnline)
-        write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
-        call CCTK_WARN(failinfomode,warnline)
         if(keyerr(1).eq.667.and.GRHydro_eos_hot_eps_fix.ne.0) then
            ! Handling of the case when no new temperature can be
            ! found for a given epsilon. The amount of times
@@ -671,6 +695,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
            nf = 0
            do while(anyerr.ne.0.and.nf.le.nfudgemax)
               anyerr = 0
+              utau = utau + rho*abs(epsilon)*0.05d0*w_lorentz**2
               epsilon = epsilon + abs(epsilon)*0.05d0
               call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
                    rho,epsilon,temp,ye,xpress,keyerr,anyerr)
@@ -693,11 +718,19 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
               call CCTK_WARN(failwarnmode,"Aborting!!!")
            endif
         else
+           call CCTK_WARN(failinfomode,"EOS error in c2p 1")
+           write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+           call CCTK_WARN(failinfomode,warnline)
+           write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye
+           call CCTK_WARN(failinfomode,warnline)
+           write(warnline,"(A7,i8)") "code: ",keyerr(1)
+           call CCTK_WARN(failinfomode,warnline)
+           write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
            call CCTK_WARN(failinfomode,warnline)
            call CCTK_WARN(failwarnmode,"Aborting!!!")
         endif
      endif
-     !$OMP END CRITICAL(c2p1)
+     !$OMP END CRITICAL
   endif
   f = pold - xpress
 
@@ -748,25 +781,19 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
          rho,epsilon,temp,ye,xpress,keyerr,anyerr)
     ! error handling
     if(anyerr.ne.0) then
-     !$OMP CRITICAL(c2p2)
+     !$OMP CRITICAL
        if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
-          call CCTK_WARN(failinfomode,"EOS error in c2p 2")
-          write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
-          call CCTK_WARN(failinfomode,warnline)
-          write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye
-          call CCTK_WARN(failinfomode,warnline)
-          write(warnline,"(A7,i8)") "code: ",keyerr(1)
-          call CCTK_WARN(failinfomode,warnline)
-          write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
-          call CCTK_WARN(failinfomode,warnline)
           if(keyerr(1).eq.667.and.GRHydro_eos_hot_eps_fix.ne.0) then
              ! Handling of the case when no new temperature can be
              ! found for a given epsilon. The amount of times
              ! we get to this place should be monitored, as this
              ! 'failsafe' mode creates artificial heat.
+             write(warnline,"(4i5,1P10E15.6)") GRHydro_reflevel,ii,jj,kk,x,y,z
+             call CCTK_WARN(0,warnline)
              nf = 0
              do while(anyerr.ne.0.and.nf.le.nfudgemax)
                 anyerr = 0
+                utau = utau + rho*abs(epsilon)*0.05d0*w_lorentz**2
                 epsilon = epsilon + abs(epsilon)*0.05d0
                 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
                      rho,epsilon,temp,ye,xpress,keyerr,anyerr)
@@ -789,12 +816,19 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
                 call CCTK_WARN(failwarnmode,"Aborting!!!")
              endif
           else
+             call CCTK_WARN(failinfomode,"EOS error in c2p 2")
+             write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+             call CCTK_WARN(failinfomode,warnline)
+             write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye
+             call CCTK_WARN(failinfomode,warnline)
+             write(warnline,"(A7,i8)") "code: ",keyerr(1)
+             call CCTK_WARN(failinfomode,warnline)
+             write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
              call CCTK_WARN(failinfomode,warnline)
              call CCTK_WARN(failwarnmode,"Aborting!!!")
-
           endif
        endif
-       !$OMP END CRITICAL(c2p2)
+       !$OMP END CRITICAL
     endif
 
     f = pnew - xpress
@@ -833,7 +867,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
 
     ! error handling
     if(anyerr.ne.0) then
-       !$OMP CRITICAL(c2p3)
+       !$OMP CRITICAL
        if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
           call CCTK_WARN(failinfomode,"EOS error in c2p 3")
           write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
@@ -846,7 +880,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
           call CCTK_WARN(failinfomode,warnline)
           call CCTK_WARN(failwarnmode,"Aborting!!!")
        endif
-       !$OMP END CRITICAL(c2p3)
+       !$OMP END CRITICAL
     endif
     
     f = pold - xpress
@@ -1767,6 +1801,7 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS)
   
   implicit none
   
+  DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_ARGUMENTS
 
   integer :: i, j, k, nx, ny, nz
@@ -1785,6 +1820,14 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS)
 
         if (GRHydro_C2P_failed(i,j,k) == 1) then
 
+           if(con2prim_oct_hack.ne.0.and.&
+                (x(i,j,k) .lt. -1.0d-12 .or.&
+                y(i,j,k) .lt. -1.0d-12 .or.&
+                z(i,j,k) .lt. -1.0d-12)) then
+              cycle
+           endif
+
+
           !$OMP CRITICAL
           call CCTK_WARN(1,'Con2Prim failed; stopping the code.')
           call CCTK_WARN(1,'Even with mesh refinement, this point is not restricted from a finer level, so this is really an error')
diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90
index 8a428c9..fbe118b 100644
--- a/src/GRHydro_Prim2Con.F90
+++ b/src/GRHydro_Prim2Con.F90
@@ -254,7 +254,7 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, &
              drho,deps,temp,ye,dpress,keyerr,anyerr)  
         keytemp=0
         if(anyerr.ne.0) then
-           !$OMP CRITICAL(p2c1)
+           !$OMP CRITICAL
            call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2")
            write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
            call CCTK_WARN(1,warnline)
@@ -264,13 +264,13 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, &
            call CCTK_WARN(1,warnline)
            write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
            call CCTK_WARN(1,warnline)
-           !$OMP END CRITICAL(p2c1)
+           !$OMP END CRITICAL
         endif
      else
         ! This is a way of recovering even on finer refinement levels:
         ! Use the average temperature at the interface instead of the
         ! reconstructed specific internal energy.
-        !$OMP CRITICAL(p2c2)
+        !$OMP CRITICAL
         call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!")
         write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
         call CCTK_WARN(1,warnline)
@@ -280,14 +280,14 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, &
         call CCTK_WARN(1,warnline)
         write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
         call CCTK_WARN(1,warnline)
-        !$OMP END CRITICAL(p2c2)
+        !$OMP END CRITICAL
         keytemp=1
         temp = temp0
         call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
              drho,deps,temp,ye,dpress,keyerr,anyerr)  
         keytemp=0
         if(anyerr.ne.0) then
-           !$OMP CRITICAL(p2c3)
+           !$OMP CRITICAL
            call CCTK_WARN(1,"EOS error in prim2con_hot")
            write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
            call CCTK_WARN(1,warnline)
@@ -298,7 +298,7 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, &
            write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
            call CCTK_WARN(1,warnline)
            call CCTK_WARN(0,"Aborting!!!")
-           !$OMP END CRITICAL(p2c3)
+           !$OMP END CRITICAL
         endif
      endif
   endif