diff options
| author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2012-07-17 17:08:26 +0000 |
|---|---|---|
| committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2012-07-17 17:08:26 +0000 |
| commit | 375c87d95be79ae5b04c2fa464fff89c41c98dd4 (patch) | |
| tree | 2080b61a8993608551f8be3c7592b74603a5420c /src | |
| parent | f0618cd697af346f5fcd0ee0e53875826723915a (diff) | |
GRHydro: ePPM: Corrected a logic bug in vel2 constraint. Also, use averaged metric.
From: Christian Reisswig <reisswig@scriwalker.(none)>
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@397 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
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) &&\ |