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

[rhaas at tapir.caltech.edu]

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

Modified:
 /trunk/src/
  GRHydro_PPM.F90

Log:
 GRHydro: ePPM: Constrain reconstructed vel to below speed of light!
 
 From: Christian Reisswig <reisswig at scriwalker.(none)>

File Changes:

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

File [modified]: GRHydro_PPM.F90
Delta lines: +23 -0
===================================================================
--- trunk/src/GRHydro_PPM.F90	2012-07-17 17:08:22 UTC (rev 395)
+++ trunk/src/GRHydro_PPM.F90	2012-07-17 17:08:25 UTC (rev 396)
@@ -224,6 +224,7 @@
       endif
 
 
+
    if (PPM3) then
       
       !! We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, 
@@ -767,6 +768,24 @@
 
 
 
+!! 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) \
+       *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 &&\
+      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
@@ -775,6 +794,8 @@
 		  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)
 	    end do
 	    if (poly .eq. 0) then
 		  do i = 3, nx - 2
@@ -788,6 +809,8 @@
 		  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)
 	    end do
 	    if (poly .eq. 0) then
 		  do i = 4, nx - 3