diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/GRHydro_PPM.F90 | 25 |
1 files changed, 20 insertions, 5 deletions
diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90 index c434438..48b1d72 100644 --- a/src/GRHydro_PPM.F90 +++ b/src/GRHydro_PPM.F90 @@ -767,16 +767,31 @@ trivial_rp = .true. 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) &&\ |