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@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) &&\