unify code: make primitive2conservativegeneral() use prim2conpolytype()

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

1 files changed, 22 insertions, 73 deletions

diff --git a/src/Whisky_Prim2Con.F90 b/src/Whisky_Prim2Con.F90
index 2c823f5..4ab24f2 100644
--- a/src/Whisky_Prim2Con.F90
+++ b/src/Whisky_Prim2Con.F90
@@ -360,8 +360,10 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
     w = 1.d0 / sqrt(1.d0 - w_tmp)
   endif
   
-  dpress = EOS_Pressure(handle, drho, deps)
-  deps = EOS_SpecificIntEnergy(handle, drho, dpress)
+  if (handle .ge. 0) then
+    dpress = EOS_Pressure(handle, drho, deps)
+    deps = EOS_SpecificIntEnergy(handle, drho, dpress)
+  end if
 
   vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
   vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
@@ -580,78 +582,25 @@ subroutine primitive2conservativegeneral(CCTK_ARGUMENTS)
         call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl)
         call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr)
 
-        w_lorentzminus(i,j,k) = 1.d0 / sqrt(1.d0 - &
-             (gxxl * velxminus(i,j,k)**2 + &
-              gyyl * velyminus(i,j,k)**2 + &
-              gzzl * velzminus(i,j,k)**2 + &
-              2.d0 * gxyl * velxminus(i,j,k) * velyminus(i,j,k) + &
-              2.d0 * gxzl * velxminus(i,j,k) * velzminus(i,j,k) + &
-              2.d0 * gyzl * velyminus(i,j,k) * velzminus(i,j,k) ))
-        vlowx = gxxl * velxminus(i,j,k) + &
-                gxyl * velyminus(i,j,k) + &
-                gxzl * velzminus(i,j,k)
-        vlowy = gxyl * velxminus(i,j,k) + &
-                gyyl * velyminus(i,j,k) + &
-                gyzl * velzminus(i,j,k)
-        vlowz = gxzl * velxminus(i,j,k) + &
-                gyzl * velyminus(i,j,k) + &
-                gzzl * velzminus(i,j,k)
-        densminus(i,j,k) = sqrt(avg_detl) * rhominus(i,j,k) * &
-             w_lorentzminus(i,j,k)
-        sxminus(i,j,k) =  sqrt(avg_detl) * &
-             ( rhominus(i,j,k) * &
-               (1.d0 + epsminus(i,j,k)) + pressminus(i,j,k) ) * &
-               w_lorentzminus(i,j,k)**2 * vlowx
-        syminus(i,j,k) =  sqrt(avg_detl) * &
-             ( rhominus(i,j,k) * &
-               (1.d0 + epsminus(i,j,k)) + pressminus(i,j,k) ) * &
-               w_lorentzminus(i,j,k)**2 * vlowy
-        szminus(i,j,k) =  sqrt(avg_detl) * &
-             ( rhominus(i,j,k) * &
-               (1.d0 + epsminus(i,j,k)) + pressminus(i,j,k) ) * &
-               w_lorentzminus(i,j,k)**2 * vlowz
-        tauminus(i,j,k) = sqrt(avg_detl) * &
-             ( ( rhominus(i,j,k) * &
-                 (1.d0 + epsminus(i,j,k)) + pressminus(i,j,k) ) * &
-                 w_lorentzminus(i,j,k)**2 - pressminus(i,j,k) ) - &
-             densminus(i,j,k)
-
-        w_lorentzplus(i,j,k) = 1.d0 / sqrt(1.d0 - &
-             (gxxr * velxplus(i,j,k)**2 + &
-              gyyr * velyplus(i,j,k)**2 + &
-              gzzr * velzplus(i,j,k)**2 + &
-              2.d0 * gxyr * velxplus(i,j,k) * velyplus(i,j,k) + &
-              2.d0 * gxzr * velxplus(i,j,k) * velzplus(i,j,k) + &
-              2.d0 * gyzr * velyplus(i,j,k) * velzplus(i,j,k) ))
-        vlowx = gxxr * velxplus(i,j,k) + &
-                gxyr * velyplus(i,j,k) + &
-                gxzr * velzplus(i,j,k)
-        vlowy = gxyr * velxplus(i,j,k) + &
-                gyyr * velyplus(i,j,k) + &
-                gyzr * velzplus(i,j,k)
-        vlowz = gxzr * velxplus(i,j,k) + &
-                gyzr * velyplus(i,j,k) + &
-                gzzr * velzplus(i,j,k)
-        densplus(i,j,k) = sqrt(avg_detr) * rhoplus(i,j,k) * &
-             w_lorentzplus(i,j,k)
-        sxplus(i,j,k) =  sqrt(avg_detr) * &
-             ( rhoplus(i,j,k) * &
-               (1.d0 + epsplus(i,j,k)) + pressplus(i,j,k) ) * &
-               w_lorentzplus(i,j,k)**2 * vlowx
-        syplus(i,j,k) =  sqrt(avg_detr) * &
-             ( rhoplus(i,j,k) * &
-               (1.d0 + epsplus(i,j,k)) + pressplus(i,j,k) ) * &
-               w_lorentzplus(i,j,k)**2 * vlowy
-        szplus(i,j,k) =  sqrt(avg_detr) * &
-             ( rhoplus(i,j,k) * &
-               (1.d0 + epsplus(i,j,k)) + pressplus(i,j,k) ) * &
-               w_lorentzplus(i,j,k)**2 * vlowz
-        tauplus(i,j,k) = sqrt(avg_detr) * &
-             ( ( rhoplus(i,j,k) * &
-                 (1.d0 + epsplus(i,j,k)) + pressplus(i,j,k) ) * &
-                 w_lorentzplus(i,j,k)**2 - pressplus(i,j,k) ) - &
-             densplus(i,j,k)
+        call prim2conpolytype(-1, gxxl,gxyl,gxzl,&
+             gyyl,gyzl,gzzl, &
+             avg_detl,densminus(i,j,k),sxminus(i,j,k),&
+             syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), &
+             velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
+             epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))
 
+        call prim2conpolytype(-1, gxxr,gxyr,gxzr,&
+             gyyr,gyzr,gzzr, &
+             avg_detr, densplus(i,j,k),sxplus(i,j,k),&
+             syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
+             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
   end do