check for velocity-overshoots from reconstruction which otherwise would produce NaNs, e.g. in w_lorentz. We cannot really do something against this, so we simply limit w_lorentz to some really large value, or similar

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

3 files changed, 52 insertions, 20 deletions

diff --git a/src/Whisky_Eigenproblem_Marquina.F90 b/src/Whisky_Eigenproblem_Marquina.F90
index 34857e4..c618c2a 100644
--- a/src/Whisky_Eigenproblem_Marquina.F90
+++ b/src/Whisky_Eigenproblem_Marquina.F90
@@ -19,6 +19,8 @@
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 
+#define MARQUINA_SMALL 1.0d-10
+
  /*@@
    @routine    eigenproblem_marquina
    @date       Wed Feb 13 12:27:59 2002
@@ -72,7 +74,7 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   CCTK_REAL reivec1r(5),reivec2r(5),reivec3r(5),reivecpr(5),reivecmr(5)
   CCTK_REAL lam1l,lam2l,lam3l,lamml,lampl,lamm_nobetal,lamp_nobetal
   CCTK_REAL lam1r,lam2r,lam3r,lammr,lampr,lamm_nobetar,lamp_nobetar
-  CCTK_REAL lam1,lam2,lam3,lamm,lamp
+  CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamp_tmp
 
   CCTK_REAL cs2l,cs2r,one,two
   CCTK_REAL vlowxr,vlowyr,vlowzr,v2r,wr
@@ -118,6 +120,9 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   vlowzl = gxz*velxl + gyz*velyl + gzz*velzl
   v2l = vlowxl*velxl + vlowyl*velyl + vlowzl*velzl
 
+  if (v2l .ge. 1.d0-MARQUINA_SMALL) then
+    v2l = 1.d0-MARQUINA_SMALL
+  endif
 !!$  Assume consistent primitive data
 
   wl = one / sqrt(one - v2l)
@@ -135,11 +140,13 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   lam1l = velxl - betainvalp
   lam2l = lam1l
   lam3l = lam1l
+  lamp_tmp = cs2l*(one-v2l)*(u*(one-v2l*cs2l) - velxl**2*(one-cs2l))
+  if (lamp_tmp .le. MARQUINA_SMALL) then
+    lamp_tmp = MARQUINA_SMALL;
+  endif
   invfacl = one/(one-v2l*cs2l)
-  lamp_nobetal = (velxl*(one-cs2l) + sqrt(cs2l*(one-v2l)* &
-       (u*(one-v2l*cs2l) - velxl**2*(one-cs2l))))*invfacl
-  lamm_nobetal = (velxl*(one-cs2l) - sqrt(cs2l*(one-v2l)* &
-       (u*(one-v2l*cs2l) - velxl**2*(one-cs2l))))*invfacl
+  lamp_nobetal = (velxl*(one-cs2l) + sqrt(lamp_tmp))*invfacl
+  lamm_nobetal = (velxl*(one-cs2l) - sqrt(lamp_tmp))*invfacl
 
 !!$ BEGIN: Check for NaN value (2nd check)
   if (lamp_nobetal .ne. lamp_nobetal) then
@@ -169,8 +176,10 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   vlowzr = gxz*velxr + gyz*velyr + gzz*velzr
   v2r = vlowxr*velxr + vlowyr*velyr + vlowzr*velzr
 
+  if (v2r .ge. 1.d0-MARQUINA_SMALL) then
+    v2r = 1.d0-MARQUINA_SMALL
+  endif
 !!$  Assume consistent primitive data
-
   wr = one / sqrt(one - v2r)
 
 !!$ BEGIN: Check for NaN value (3rd check)
@@ -185,11 +194,13 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   lam1r = velxr - betainvalp
   lam2r = lam1r
   lam3r = lam1r
+  lamp_tmp = cs2r*(one-v2r)*(u*(one-v2r*cs2r) - velxr**2*(one-cs2r))
+  if (lamp_tmp .le. MARQUINA_SMALL) then
+    lamp_tmp = MARQUINA_SMALL;
+  endif
   invfacr = one/(one-v2r*cs2r)
-  lamp_nobetar = (velxr*(one-cs2r) + sqrt(cs2r*(one-v2r)* &
-       (u*(one-v2r*cs2r) - velxr**2*(one-cs2r))))*invfacr
-  lamm_nobetar = (velxr*(one-cs2r) - sqrt(cs2r*(one-v2r)* &
-       (u*(one-v2r*cs2r) - velxr**2*(one-cs2r))))*invfacr
+  lamp_nobetar = (velxr*(one-cs2r) + sqrt(lamp_tmp))*invfacr
+  lamm_nobetar = (velxr*(one-cs2r) - sqrt(lamp_tmp))*invfacr
 
 !!$ BEGIN: Check for NaN value (4th check)
   if (lamp_nobetar .ne. lamp_nobetar) then
diff --git a/src/Whisky_Marquina.F90 b/src/Whisky_Marquina.F90
index b1f09ca..35b240a 100644
--- a/src/Whisky_Marquina.F90
+++ b/src/Whisky_Marquina.F90
@@ -49,7 +49,7 @@ subroutine Whisky_Marquina(CCTK_ARGUMENTS)
     CCTK_REAL :: avg_alp,avg_beta,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
          avg_det,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,&
          pressp,pressm_i, &
-         w_lorentzp,w_lorentzm_i, usendh, psi4h
+         tmp_w_lorentzp, tmp_w_lorentzm_i, w_lorentzp,w_lorentzm_i, usendh, psi4h
     integer :: m
     integer :: i,j,k
     character(len=256) NaN_WarnLine
@@ -78,7 +78,8 @@ subroutine Whisky_Marquina(CCTK_ARGUMENTS)
     !$OMP PARALLEL DO PRIVATE(i,j,consp,consm_i,primp,primm_i,&
     !$OMP marquinaflux,avg_beta,avg_alp,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
     !$OMP psi4h,f_marquina,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,usendh,&
-    !$OMP w_lorentzp,w_lorentzm_i,fplus,fminus,m,avg_det)
+    !$OMP tmp_w_lorentzp, tmp_w_lorentzm_i,w_lorentzp,w_lorentzm_i,&
+    !$OMP fplus,fminus,m,avg_det)
     do k = whisky_stencil, cctk_lsh(3) - whisky_stencil
       do j = whisky_stencil, cctk_lsh(2) - whisky_stencil
         do i = whisky_stencil, cctk_lsh(1) - whisky_stencil
@@ -199,10 +200,16 @@ subroutine Whisky_Marquina(CCTK_ARGUMENTS)
 
 !!$left state
 
-            w_lorentzp = 1.d0 / sqrt(1.d0 - (gxxh*primp(2)*primp(2) + &
+            tmp_w_lorentzp = gxxh*primp(2)*primp(2) + &
                  gyyh*primp(3)*primp(3) + gzzh*primp(4)*primp(4) + &
                  2*gxyh*primp(2)*primp(3) + 2*gxzh*primp(2) *primp(4) + &
-                 2*gyzh*primp(3)*primp(4)))  
+                 2*gyzh*primp(3)*primp(4)
+            if (tmp_w_lorentzp .ge. 1.d0) then
+              w_lorentzp = Whisky_lorentz_overshoot_cutoff
+            else
+              w_lorentzp = 1.d0 / sqrt(1.d0 - tmp_w_lorentzp);
+            endif
+
 
 !!$ BEGIN: Check for NaN value (1st check)
   if (w_lorentzp .ne. w_lorentzp) then
@@ -217,11 +224,16 @@ subroutine Whisky_Marquina(CCTK_ARGUMENTS)
         
 !!$right state
 
-            w_lorentzm_i = 1.d0 / sqrt(1.d0 - (gxxh*primm_i(2)*primm_i(2) + &
+            tmp_w_lorentzm_i = gxxh*primm_i(2)*primm_i(2) + &
                  gyyh*primm_i(3)*primm_i(3) + gzzh*primm_i(4)*primm_i(4) + &
                  2*gxyh*primm_i(2)*primm_i(3) + &
                  2*gxzh*primm_i(2) *primm_i(4)+ &
-                 2*gyzh*primm_i(3)*primm_i(4)))  
+                 2*gyzh*primm_i(3)*primm_i(4)
+            if (tmp_w_lorentzm_i .ge. 1.d0) then
+              w_lorentzm_i = Whisky_lorentz_overshoot_cutoff
+            else
+              w_lorentzm_i = 1.d0 / sqrt(1.d0 - tmp_w_lorentzm_i);
+            endif
 
 !!$ BEGIN: Check for NaN value (2nd check)
   if (w_lorentzm_i .ne. w_lorentzm_i) then
diff --git a/src/Whisky_Prim2Con.F90 b/src/Whisky_Prim2Con.F90
index 3669516..eae7ec8 100644
--- a/src/Whisky_Prim2Con.F90
+++ b/src/Whisky_Prim2Con.F90
@@ -295,6 +295,11 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS)
              rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
              velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
              w_lorentzplus(i,j,k)) 
+  if (densminus(i,j,k) .ne. densminus(i,j,k)) then
+    !$OMP CRITICAL
+    call CCTK_WARN(0, "NaN in densminus(i,j,k) (Prim2Con)")
+    !$OMP END CRITICAL
+  endif
 
       end do
     end do
@@ -326,7 +331,7 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
   
   CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
   CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,&
-       deps, dpress, w, vlowx, vlowy, vlowz, sqrtdet   
+       deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet   
   CCTK_INT :: handle
 
 #ifdef _EOS_BASE_INC_
@@ -334,9 +339,13 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
 #endif
 #include "EOS_Base.inc"
   
-  w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz &
-       *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz&
-       *dvely*dvelz))  
+  w_tmp = gxx*dvelx*dvelx + gyy*dvely*dvely + gzz *dvelz*dvelz + &
+          2*gxy*dvelx*dvely + 2*gxz*dvelx*dvelz + 2*gyz*dvely*dvelz
+  if (w_tmp .ge. 1.d0) then
+    w = 1.d100
+  else
+    w = 1.d0 / sqrt(1.d0 - w_tmp)
+  endif
   
   dpress = EOS_Pressure(handle, drho, deps)
   deps = EOS_SpecificIntEnergy(handle, drho, dpress)