GRHydro: Option to reconstruct W*vel (W = Lorentz factor) instead of just the primitive velocity 'vel'.

This ensures that the Lorentz factor can never become unphysical during reconstruction.
This patch is necessary for NS collapse when ePPM+Refluxing is used.

The idea is due to Roland.

From: Christian Reisswig <reisswig@scriwalker.(none)>

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

3 files changed, 218 insertions, 130 deletions

diff --git a/param.ccl b/param.ccl
index a89d8cc..e66e529 100644
--- a/param.ccl
+++ b/param.ccl
@@ -231,6 +231,9 @@ real enhanced_ppm_C4 "Parameter for enhancecd ppm limiter. This is for the addit
    0:1 :: "must be greater than or equal to 0, and not larger than 1"
 } 0.0
 
+boolean reconstruct_Wv "Reconstruct the primitive velocity W_Lorentz*vel, rather than just vel." STEERABLE=ALWAYS
+{
+} no
 
 int eno_order "The order of accuracy of the ENO reconstruction"
 {
diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90
index 2ad0187..50da41c 100644
--- a/src/GRHydro_PPM.F90
+++ b/src/GRHydro_PPM.F90
@@ -106,6 +106,24 @@ subroutine SimplePPM_1d(handle,poly,&
 
   logical :: cond
 
+  CCTK_REAL, dimension(nx) :: vx,vy,vz
+  
+  
+  !!$ reconstruct w^i = w_lorentz*v^i to ensure slower than light
+  !!$ speeds?
+  if (reconstruct_Wv.ne.0) then
+    ! all velocity like quantities are now w_lorentz*velocity. We will convert
+    ! back to ordinary velocity at the end, using w_lorentz = sqrt(1 + g_{ij}
+    ! w^i w^j). 
+    vx = w_lorentz*velx
+    vy = w_lorentz*vely
+    vz = w_lorentz*velz
+  else
+    vx = velx
+    vy = vely
+    vz = velz
+  end if
+
 
 !!$  Initially, all the Riemann problems will be trivial
 trivial_rp = .true.
@@ -129,9 +147,9 @@ trivial_rp = .true.
 
       do i = 2, nx - 1
 	 drho(i) = 0.5d0 * (rho(i+1) - rho(i-1))
-	 dvelx(i) = 0.5d0 * (velx(i+1) - velx(i-1))
-	 dvely(i) = 0.5d0 * (vely(i+1) - vely(i-1))
-	 dvelz(i) = 0.5d0 * (velz(i+1) - velz(i-1))
+	 dvelx(i) = 0.5d0 * (vx(i+1) - vx(i-1))
+	 dvely(i) = 0.5d0 * (vy(i+1) - vy(i-1))
+	 dvelz(i) = 0.5d0 * (vz(i+1) - vz(i-1))
 	 dpress(i) = press(i+1) - press(i-1)
 	 d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1))! / 6.d0 / dx / dx
 	 ! since we use d2rho only for the condition d2rho(i+1)*d2rhoi(i-1)<0 
@@ -147,9 +165,9 @@ trivial_rp = .true.
 
       do i = 2, nx - 1
 	 STEEP(rho, drho, dmrho)
-	 STEEP(velx, dvelx, dmvelx)
-	 STEEP(vely, dvely, dmvely)
-	 STEEP(velz, dvelz, dmvelz)
+	 STEEP(vx, dvelx, dmvelx)
+	 STEEP(vy, dvely, dmvely)
+	 STEEP(vz, dvelz, dmvelz)
       end do
       if (poly .eq. 0) then
 	 do i = 2, nx - 1
@@ -163,13 +181,13 @@ trivial_rp = .true.
 	 rhoplus(i) = 0.5d0 * (rho(i) + rho(i+1)) + &
 	       (dmrho(i) - dmrho(i+1)) / 6.d0
 	 rhominus(i+1) = rhoplus(i)
-	 velxplus(i) = 0.5d0 * (velx(i) + velx(i+1)) + &
+	 velxplus(i) = 0.5d0 * (vx(i) + vx(i+1)) + &
 	       (dmvelx(i) - dmvelx(i+1)) / 6.d0
 	 velxminus(i+1) = velxplus(i)
-	 velyplus(i) = 0.5d0 * (vely(i) + vely(i+1)) + &
+	 velyplus(i) = 0.5d0 * (vy(i) + vy(i+1)) + &
 	       (dmvely(i) - dmvely(i+1)) / 6.d0
 	 velyminus(i+1) = velyplus(i)
-	 velzplus(i) = 0.5d0 * (velz(i) + velz(i+1)) + &
+	 velzplus(i) = 0.5d0 * (vz(i) + vz(i+1)) + &
 	       (dmvelz(i) - dmvelz(i+1)) / 6.d0
 	 velzminus(i+1) = velzplus(i)
       end do
@@ -235,16 +253,16 @@ trivial_rp = .true.
          LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i))
          rhominus(i+1) = rhoplus(i)
          
-         APPROX_AT_CELL_INTERFACE(velx, velxplus(i))
-         LIMIT(velx, velxplus(i), enhanced_ppm_C2, velxplus(i))
+         APPROX_AT_CELL_INTERFACE(vx, velxplus(i))
+         LIMIT(vx, velxplus(i), enhanced_ppm_C2, velxplus(i))
          velxminus(i+1) = velxplus(i)
          
-         APPROX_AT_CELL_INTERFACE(vely, velyplus(i))
-         LIMIT(vely, velyplus(i), enhanced_ppm_C2, velyplus(i))
+         APPROX_AT_CELL_INTERFACE(vy, velyplus(i))
+         LIMIT(vy, velyplus(i), enhanced_ppm_C2, velyplus(i))
          velyminus(i+1) = velyplus(i)
          
-         APPROX_AT_CELL_INTERFACE(velz, velzplus(i))
-         LIMIT(velz, velzplus(i), enhanced_ppm_C2, velzplus(i))
+         APPROX_AT_CELL_INTERFACE(vz, velzplus(i))
+         LIMIT(vz, velzplus(i), enhanced_ppm_C2, velzplus(i))
          velzminus(i+1) = velzplus(i)
       end do
       if (poly .eq. 0) then
@@ -264,16 +282,16 @@ trivial_rp = .true.
          LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i))
          rhominus(i+1) = rhoplus(i)
          
-         APPROX_AT_CELL_INTERFACE_STENCIL4(velx, velxplus(i))
-         LIMIT(velx, velxplus(i), enhanced_ppm_C2, velxplus(i))
+         APPROX_AT_CELL_INTERFACE_STENCIL4(vx, velxplus(i))
+         LIMIT(vx, velxplus(i), enhanced_ppm_C2, velxplus(i))
          velxminus(i+1) = velxplus(i)
          
-         APPROX_AT_CELL_INTERFACE_STENCIL4(vely, velyplus(i))
-         LIMIT(vely, velyplus(i), enhanced_ppm_C2, velyplus(i))
+         APPROX_AT_CELL_INTERFACE_STENCIL4(vy, velyplus(i))
+         LIMIT(vy, velyplus(i), enhanced_ppm_C2, velyplus(i))
          velyminus(i+1) = velyplus(i)
          
-         APPROX_AT_CELL_INTERFACE_STENCIL4(velz, velzplus(i))
-         LIMIT(velz, velzplus(i), enhanced_ppm_C2, velzplus(i))
+         APPROX_AT_CELL_INTERFACE_STENCIL4(vz, velzplus(i))
+         LIMIT(vz, velzplus(i), enhanced_ppm_C2, velzplus(i))
          velzminus(i+1) = velzplus(i)
       end do
       if (poly .eq. 0) then
@@ -444,7 +462,7 @@ trivial_rp = .true.
 !!$  Zone flattening. See appendix of C&W, p. 197-8.
   do i = 3, nx - 2
     dpress2 = press(i+2) - press(i-2)
-    dvel = velx(i+1) - velx(i-1)
+    dvel = vx(i+1) - vx(i-1)
     if ( (abs(dpress(i)) >  ppm_epsilon * min(press(i-1),press(i+1))) .and. &
          (dvel < 0.d0) ) then
       w = 1.d0
@@ -471,12 +489,12 @@ trivial_rp = .true.
 	    end if
 	    rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
 	    rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
-	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i)
-	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i)
-	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i)
-	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i)
-	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i)
-	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i)
+	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i)
+	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i)
+	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i)
+	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i)
+	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i)
+	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i)
 	    if (poly .eq. 0) then
 	       epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
 	       epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
@@ -492,12 +510,12 @@ trivial_rp = .true.
 	    end if
 	    rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
 	    rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
-	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i)
-	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i)
-	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i)
-	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i)
-	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i)
-	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i)
+	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i)
+	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i)
+	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i)
+	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i)
+	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i)
+	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i)
 	    if (poly .eq. 0) then
 	       epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
 	       epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
@@ -767,50 +785,14 @@ trivial_rp = .true.
       endif   
 
 
-! Reconstruct quantity on plus face
-#define AVGP(a) \
-   0.5*(a(i)+a(i+1))
-
-! Reconstruct quantity on minus face
-#define AVGM(a) \
-   0.5*(a(i)+a(i-1))
-
-
-!! Compute scalar product v^i v_i on plus face
-#define VEL2P(dvelx,dvely,dvelz) \
-   (AVGP(gxx)*dvelx*dvelx + AVGP(gyy)*dvely*dvely + AVGP(gzz)  \
-       *dvelz*dvelz + 2*AVGP(gxy)*dvelx*dvely + 2*AVGP(gxz)*dvelx *dvelz + 2*AVGP(gyz) \
-       *dvely*dvelz)
-
-!! Compute scalar product v^i v_i on minus face
-#define VEL2M(dvelx,dvely,dvelz) \
-   (AVGM(gxx)*dvelx*dvelx + AVGM(gyy)*dvely*dvely + AVGM(gzz)  \
-       *dvelz*dvelz + 2*AVGM(gxy)*dvelx*dvely + 2*AVGM(gxz)*dvelx *dvelz + 2*AVGM(gyz) \
-       *dvely*dvelz)
-
-
-!! Check if velocity is below speed of light. If not, reduce reconstruction to first order!
-#define VELCHECK(vxminus,vx,vxplus,vyminus,vy,vyplus,vzminus,vz,vzplus) \
-   if (VEL2M(vxminus(i),vyminus(i),vzminus(i)) .ge. (1.0d0-enhanced_ppm_C4) .or. VEL2P(vxplus(i),vyplus(i),vzplus(i)) .ge. (1.0d0-enhanced_ppm_C4)) then &&\
-      vxminus(i) = vx(i)  &&\
-      vxplus(i) = vx(i)   &&\
-      vyminus(i) = vy(i)  &&\
-      vyplus(i) = vy(i)   &&\
-      vzminus(i) = vz(i)  &&\
-      vzplus(i) = vz(i)   &&\
-   endif
-
-
    if (use_enhanced_ppm .eq. 1) then
       !! Constrain parabolic profiles, PPM 2011/2008
       if (PPM3) then
 	    do i = 3, nx - 2
 		  MON_WITH_LOCAL_EXTREMUM(rhominus,rho,rhoplus, enhanced_ppm_C2)
-		  MON_WITH_LOCAL_EXTREMUM(velxminus,velx,velxplus, enhanced_ppm_C2)
-		  MON_WITH_LOCAL_EXTREMUM(velyminus,vely,velyplus, enhanced_ppm_C2)
-		  MON_WITH_LOCAL_EXTREMUM(velzminus,velz,velzplus, enhanced_ppm_C2)
-		  ! Ensure that velocity is below speed of light!
-		  VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus)
+		  MON_WITH_LOCAL_EXTREMUM(velxminus,vx,velxplus, enhanced_ppm_C2)
+		  MON_WITH_LOCAL_EXTREMUM(velyminus,vy,velyplus, enhanced_ppm_C2)
+		  MON_WITH_LOCAL_EXTREMUM(velzminus,vz,velzplus, enhanced_ppm_C2)
 	    end do
 	    if (poly .eq. 0) then
 		  do i = 3, nx - 2
@@ -821,11 +803,9 @@ trivial_rp = .true.
       else
 	    do i = 4, nx - 3
 		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(rhominus,rho,rhoplus, enhanced_ppm_C2, enhanced_ppm_C3)
-		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,velx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3)
-		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vely,velyplus, enhanced_ppm_C2, enhanced_ppm_C3)
-		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,velz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3)
-		  ! Ensure that velocity is below speed of light!
-		  VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus)
+		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,vx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3)
+		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vy,velyplus, enhanced_ppm_C2, enhanced_ppm_C3)
+		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,vz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3)
 	    end do
 	    if (poly .eq. 0) then
 		  do i = 4, nx - 3
@@ -844,12 +824,12 @@ trivial_rp = .true.
 	    end if
 	    rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
 	    rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
-	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i)
-	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i)
-	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i)
-	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i)
-	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i)
-	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i)
+	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i)
+	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i)
+	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i)
+	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i)
+	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i)
+	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i)
 	    if (poly .eq. 0) then
 	       epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
 	       epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
@@ -865,12 +845,12 @@ trivial_rp = .true.
 	    end if
 	    rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
 	    rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
-	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i)
-	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i)
-	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i)
-	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i)
-	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i)
-	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i)
+	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i)
+	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i)
+	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i)
+	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i)
+	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i)
+	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i)
 	    if (poly .eq. 0) then
 	       epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
 	       epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
@@ -885,9 +865,9 @@ trivial_rp = .true.
      do i = GRHydro_stencil, nx - GRHydro_stencil + 1
 	 ! original Colella&Woodward monotonicity preservation
 	 MON(rhominus,rho,rhoplus)
-	 MON(velxminus,velx,velxplus)
-	 MON(velyminus,vely,velyplus)
-	 MON(velzminus,velz,velzplus)
+	 MON(velxminus,vx,velxplus)
+	 MON(velyminus,vy,velyplus)
+	 MON(velzminus,vz,velzplus)
      end do
      if (poly .eq. 0) then
 	 do i = GRHydro_stencil, nx - GRHydro_stencil + 1
@@ -919,20 +899,20 @@ trivial_rp = .true.
       if (cond) then
         rhominus(i)=rho(i)
         rhoplus(i)=rho(i)
-        velxminus(i)=velx(i)
-        velxplus(i)=velx(i)
-        velyminus(i)=vely(i)
-        velyplus(i)=vely(i)
-        velzminus(i)=velz(i)
-        velzplus(i)=velz(i)
+        velxminus(i)=vx(i)
+        velxplus(i)=vx(i)
+        velyminus(i)=vy(i)
+        velyplus(i)=vy(i)
+        velzminus(i)=vz(i)
+        velzplus(i)=vz(i)
         rhominus(i-1)=rho(i)
         rhoplus(i-1)=rho(i)
-        velxminus(i-1)=velx(i)
-        velxplus(i-1)=velx(i)
-        velyminus(i-1)=vely(i)
-        velyplus(i-1)=vely(i)
-        velzminus(i-1)=velz(i)
-        velzplus(i-1)=velz(i)
+        velxminus(i-1)=vx(i)
+        velxplus(i-1)=vx(i)
+        velyminus(i-1)=vy(i)
+        velyplus(i-1)=vy(i)
+        velzminus(i-1)=vz(i)
+        velzplus(i-1)=vz(i)
         if (poly .eq. 0) then
           epsminus(i)=eps(i)
           epsplus(i)=eps(i)
@@ -946,9 +926,9 @@ trivial_rp = .true.
         end if
         if (cond) then
           call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i))
-          call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i))
-          call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i))
-          call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i))
+          call PPM_TVD(vx(i-1), vx(i), vx(i+1), velxminus(i), velxplus(i))
+          call PPM_TVD(vy(i-1), vy(i), vy(i+1), velyminus(i), velyplus(i))
+          call PPM_TVD(vz(i-1), vz(i), vz(i+1), velzminus(i), velzplus(i))
           if (poly .eq. 0) then
             call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i))
           end if
@@ -961,20 +941,20 @@ trivial_rp = .true.
       if (cond) then
         rhominus(i)=rho(i)
         rhoplus(i)=rho(i)
-        velxminus(i)=velx(i)
-        velxplus(i)=velx(i)
-        velyminus(i)=vely(i)
-        velyplus(i)=vely(i)
-        velzminus(i)=velz(i)
-        velzplus(i)=velz(i)
+        velxminus(i)=vx(i)
+        velxplus(i)=vx(i)
+        velyminus(i)=vy(i)
+        velyplus(i)=vy(i)
+        velzminus(i)=vz(i)
+        velzplus(i)=vz(i)
         rhominus(i+1)=rho(i)
         rhoplus(i+1)=rho(i)
-        velxminus(i+1)=velx(i)
-        velxplus(i+1)=velx(i)
-        velyminus(i+1)=vely(i)
-        velyplus(i+1)=vely(i)
-        velzminus(i+1)=velz(i)
-        velzplus(i+1)=velz(i)
+        velxminus(i+1)=vx(i)
+        velxplus(i+1)=vx(i)
+        velyminus(i+1)=vy(i)
+        velyplus(i+1)=vy(i)
+        velzminus(i+1)=vz(i)
+        velzplus(i+1)=vz(i)
         if (poly .eq. 0) then
           epsminus(i)=eps(i)
           epsplus(i)=eps(i)
@@ -988,9 +968,9 @@ trivial_rp = .true.
         end if
         if (cond) then
           call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i))
-          call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i))
-          call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i))
-          call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i))
+          call PPM_TVD(vx(i-1), vx(i), vx(i+1), velxminus(i), velxplus(i))
+          call PPM_TVD(vy(i-1), vy(i), vy(i+1), velyminus(i), velyplus(i))
+          call PPM_TVD(vz(i-1), vz(i), vz(i+1), velzminus(i), velzplus(i))
           if (poly .eq. 0) then
             call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i))
           end if
@@ -998,6 +978,50 @@ trivial_rp = .true.
       end if
     end if
   end do
+  
+  !!$ transform back to vel if Wv was used!
+  if (reconstruct_Wv.ne.0) then
+    do i = grhydro_stencil, nx - grhydro_stencil + 1
+      ! divide out the Loretnz factor obtained from w_lorentz =
+      ! sqrt(1+g_{ij} w^i w^j) for both the
+      ! plus and minus quantities this should by construction ensure
+      ! that any Lorentz factor calculated from them later on is
+      ! physical (ie. > 1.d0)
+      agxx = 0.5d0*( gxx(i) + gxx(i-1) )
+      agxy = 0.5d0*( gxy(i) + gxy(i-1) )
+      agxz = 0.5d0*( gxz(i) + gxz(i-1) )
+      agyy = 0.5d0*( gyy(i) + gyy(i-1) )
+      agyz = 0.5d0*( gyz(i) + gyz(i-1) )
+      agzz = 0.5d0*( gzz(i) + gzz(i-1) )
+      w = sqrt( 1.d0 + agxx*velxminus(i)*velxminus(i) + &
+                       agyy*velyminus(i)*velyminus(i) &
+              + agzz*velzminus(i)*velzminus(i) + &
+                2.d0*agxy*velxminus(i)*velyminus(i) &
+              + 2.d0*agxz*velxminus(i)*velzminus(i) + &
+                2.d0*agyz*velyminus(i)*velzminus(i) )
+      velxminus(i) = velxminus(i)/w
+      velyminus(i) = velyminus(i)/w
+      velzminus(i) = velzminus(i)/w
+      
+      agxx = 0.5d0*( gxx(i) + gxx(i+1) )
+      agxy = 0.5d0*( gxy(i) + gxy(i+1) )
+      agxz = 0.5d0*( gxz(i) + gxz(i+1) )
+      agyy = 0.5d0*( gyy(i) + gyy(i+1) )
+      agyz = 0.5d0*( gyz(i) + gyz(i+1) )
+      agzz = 0.5d0*( gzz(i) + gzz(i+1) )
+      w = sqrt( 1.d0 + agxx*velxplus(i)*velxplus(i) + &
+                agyy*velyplus(i)*velyplus(i) &
+              + agzz*velzplus(i)*velzplus(i) + &
+               2.d0*agxy*velxplus(i)*velyplus(i) &
+              + 2.d0*agxz*velxplus(i)*velzplus(i) + &
+               2.d0*agyz*velyplus(i)*velzplus(i) )
+      velxplus(i) = velxplus(i)/w
+      velyplus(i) = velyplus(i)/w
+      velzplus(i) = velzplus(i)/w
+    end do
+  end if
+  
+  
 return
 
 end subroutine SimplePPM_1d
diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index 238808a..b0db4b5 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -36,6 +36,8 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
 
   ! save memory when MP is not used
   ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
+  TARGET gaa, gab, gac, gbb, gbc, gcc
+  TARGET gxx, gxy, gxz, gyy, gyz, gzz
   TARGET lvel, vel
   
   DECLARE_CCTK_ARGUMENTS
@@ -45,14 +47,28 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
   CCTK_INT  :: i,j,k
   CCTK_REAL :: local_min_tracer
   CCTK_INT :: type_bits, not_trivial
+  CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz, w
 
   ! save memory when MP is not used
   CCTK_INT :: GRHydro_UseGeneralCoordinates
+  CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
   CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup
 
   if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+    g11 => gaa
+    g12 => gab
+    g13 => gac
+    g22 => gbb
+    g23 => gbc
+    g33 => gcc
     vup => lvel
   else
+    g11 => gxx
+    g12 => gxy
+    g13 => gxz
+    g22 => gyy
+    g23 => gyz
+    g33 => gzz
     vup => vel
   end if
 
@@ -101,7 +117,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
     call CCTK_WARN(0, "Flux direction not x,y,z")
   end if
 
-  !$OMP PARALLEL DO PRIVATE(i,j,k)
+  !$OMP PARALLEL DO PRIVATE(i,j,k, agxx, agxy, agxz, agyy, agyz, agzz, w)
   do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1
     do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1
       do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1
@@ -112,12 +128,57 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
                  rhominus(i,j,k) = rho(i,j,k)
                  epsplus(i,j,k)  = eps(i,j,k)
                  epsminus(i,j,k) = eps(i,j,k)
-                 velxplus(i,j,k) = vup(i,j,k,1)
-                 velyplus(i,j,k) = vup(i,j,k,2)
-                 velzplus(i,j,k) = vup(i,j,k,3)
-                 velxminus(i,j,k) = vup(i,j,k,1)
-                 velyminus(i,j,k) = vup(i,j,k,2)
-                 velzminus(i,j,k) = vup(i,j,k,3)
+                 if (reconstruct_Wv.ne.0) then
+                    ! divide out the Loretnz factor for both the
+		    ! plus and minus quantities this should by construction ensure
+		    ! that any Lorentz factor calculated from them later on is
+		    ! physical (ie. > 1.d0)
+		    velxplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1)
+                    velyplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,2)
+                    velzplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,3)
+                    agxx = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset))
+		    agxy = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset))
+		    agxz = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset))
+		    agyy = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset))
+		    agyz = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset))
+		    agzz = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset))
+		    w = sqrt( 1.d0 + agxx*velxplus(i,j,k)*velxplus(i,j,k) + &
+			      agyy*velyplus(i,j,k)*velyplus(i,j,k) &
+			   + agzz*velzplus(i,j,k)*velzplus(i,j,k) + &
+			      2.d0*agxy*velxplus(i,j,k)*velyplus(i,j,k) &
+			   + 2.d0*agxz*velxplus(i,j,k)*velzplus(i,j,k) + &
+			      2.d0*agyz*velyplus(i,j,k)*velzplus(i,j,k) )
+		    velxplus(i,j,k) = velxplus(i,j,k)/w
+		    velyplus(i,j,k) = velyplus(i,j,k)/w
+		    velzplus(i,j,k) = velzplus(i,j,k)/w
+                    
+                    velxminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1)
+                    velyminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,2)
+                    velzminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,3)
+		    agxx = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset))
+		    agxy = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset))
+		    agxz = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset))
+		    agyy = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset))
+		    agyz = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset))
+		    agzz = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset))
+		    w = sqrt( 1.d0 + agxx*velxminus(i,j,k)*velxminus(i,j,k) + &
+				    agyy*velyminus(i,j,k)*velyminus(i,j,k) &
+			   + agzz*velzminus(i,j,k)*velzminus(i,j,k) + &
+			      2.d0*agxy*velxminus(i,j,k)*velyminus(i,j,k) &
+			   + 2.d0*agxz*velxminus(i,j,k)*velzminus(i,j,k) + &
+			      2.d0*agyz*velyminus(i,j,k)*velzminus(i,j,k) )
+                    velxminus(i,j,k) = velxminus(i,j,k)/w
+                    velyminus(i,j,k) = velyminus(i,j,k)/w
+                    velzminus(i,j,k) = velzminus(i,j,k)/w
+                 else
+                    ! This is the standard way of doing it
+                    velxplus(i,j,k) = vup(i,j,k,1)
+                    velyplus(i,j,k) = vup(i,j,k,2)
+                    velzplus(i,j,k) = vup(i,j,k,3)
+                    velxminus(i,j,k) = vup(i,j,k,1)
+                    velyminus(i,j,k) = vup(i,j,k,2)
+                    velzminus(i,j,k) = vup(i,j,k,3)
+                 endif
                  if(evolve_y_e.ne.0) then
                     y_e_plus(i,j,k)  = y_e(i,j,k)
                     y_e_minus(i,j,k) = y_e(i,j,k)