[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 397)

[rhaas at tapir.caltech.edu]

User: rhaas
Date: 2012/07/17 12:08 PM

Modified:
 /trunk/src/
  GRHydro_PPM.F90

Log:
 GRHydro: ePPM: Corrected a logic bug in vel2 constraint. Also, use averaged metric.
 
 From: Christian Reisswig <reisswig at scriwalker.(none)>

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_PPM.F90
Delta lines: +20 -5
===================================================================
--- trunk/src/GRHydro_PPM.F90	2012-07-17 17:08:25 UTC (rev 396)
+++ trunk/src/GRHydro_PPM.F90	2012-07-17 17:08:26 UTC (rev 397)
@@ -767,16 +767,31 @@
       endif   
 
 
+! Reconstruct quantity on plus face
+#define AVGP(a) \
+   0.5*(a(i)+a(i+1))
 
-!! Compute scalar product v^i v_i
-#define VEL2(dvelx,dvely,dvelz) \
-   (gxx(i)*dvelx*dvelx + gyy(i)*dvely*dvely + gzz(i)  \
-       *dvelz*dvelz + 2*gxy(i)*dvelx*dvely + 2*gxz(i)*dvelx *dvelz + 2*gyz(i) \
+! 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 (1.0d0 .ge. VEL2(vxminus(i),vyminus(i),vzminus(i)) .or. 1.0d0 .ge. VEL2(vxplus(i),vyplus(i),vzplus(i))) then &&\
+   if (VEL2M(vxminus(i),vyminus(i),vzminus(i)) .ge. 1.0d0 .or. VEL2P(vxplus(i),vyplus(i),vzplus(i)) .ge. 1.0d0) then &&\
       vxminus(i) = vx(i)  &&\
       vxplus(i) = vx(i)   &&\
       vyminus(i) = vy(i)  &&\