GRHydro: add templated (yay\!) PPM routine; can do old ppm, but for MHD, hot, cold whatever

From: Christian Ott <cott@tapir.caltech.edu>

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@573 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

2 files changed, 415 insertions, 0 deletions

diff --git a/src/GRHydro_PPM.cc b/src/GRHydro_PPM.cc
new file mode 100644
index 0000000..e1f426f
--- /dev/null
+++ b/src/GRHydro_PPM.cc
@@ -0,0 +1,414 @@
+#include <cmath>
+#include <algorithm>
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#define MIN(a,b) (((a)<(b))?(a):(b))
+#define MAX(a,b) (((a)>(b))?(a):(b))
+
+using namespace std;
+
+/*
+  Cases that must be considered:
+  * basic hydro
+  * hydro + temperature + ye
+  * hydro + ye
+  * basic mhd
+  * mhd + temperature + ye 
+  * mhd + ye 
+  * mppm (not supported right now)
+  * not supporting trivial_rp
+  * with or without divergence cleaning
+ */
+
+
+static inline void steep(const double *restrict x, 
+			 const double *restrict dx, 
+			 double *restrict dmx, 
+			 const int i) {
+  if ( (x[i+1] - x[i]) * (x[i]-x[i-1]) > 0.0 ) {
+    dmx[i] = copysign(1.0,dx[i]) * MIN(fabs(dx[i]),
+					MIN(2.0*fabs(x[i]-x[i-1]),
+					     2.0*fabs(x[i+1]-x[i])));
+  } else {
+    dmx[i] = 0.0;
+  }
+}
+
+static inline double approx_at_cell_interface(double* a, const int i) {
+  return 7.0/12.0*(a[i]+a[i+1]) - 1.0/12.0*(a[i-1]+a[i+2]);
+}
+
+static inline void monotonize(double* restrict xminus,
+			      const double* restrict x,
+			      double* restrict xplus,
+			      const int i) {
+
+  if (  !(xplus[i]==x[i] && x[i]==xminus[i]) 
+	&& ( (xplus[i]-x[i])*(x[i]-xminus[i]) <= 0.0 ) ) 
+    {
+      xminus[i] = x[i];
+      xplus[i] = x[i];
+    }  else if( 6.0 * (xplus[i]-xminus[i]) * 
+		(x[i]-0.5*(xplus[i]+xminus[i])) >
+		(xplus[i]-xminus[i])*(xplus[i]-xminus[i]) )
+    {
+      xminus[i] = 3.0*x[i]-2.0*xplus[i]; 
+    } else if( 6.0 * (xplus[i]-xminus[i]) * 
+	       (x[i]-0.5*(xplus[i]+xminus[i])) <
+	       -(xplus[i]-xminus[i])*(xplus[i]-xminus[i]) ) 
+    {
+      xplus[i] = 3.0*x[i]-2.0*xminus[i]; 
+    }
+  
+  return;
+}	     	    
+
+
+
+template<bool do_temp, bool do_ye, bool do_mhd, 
+	 bool dc_flag, bool do_ppm_detect>
+static void ppm1d_cxx(const int nx, 
+		      const double dx, 
+		      const double* restrict rho, 
+		      const double* restrict velx,
+		      const double* restrict vely, 
+		      const double* restrict velz, 
+		      const double* restrict eps, 
+		      const double* restrict press,
+		      const double* restrict temp,
+		      const double* restrict ye,
+		      const double* restrict Bvcx,
+		      const double* restrict Bvcy,
+		      const double* restrict Bvcz,
+		      const double* restrict psidc,
+		      double* restrict rhominus, 
+		      double* restrict velxminus, 
+		      double* restrict velyminus,
+		      double* restrict velzminus, 
+		      double* restrict epsminus, 
+		      double* restrict tempminus, 
+		      double* restrict yeminus, 
+		      double* restrict Bvcxminus,
+		      double* restrict Bvcyminus,
+		      double* restrict Bvczminus,
+		      double* restrict psidcminus,
+		      double* restrict rhoplus, 
+		      double* restrict velxplus, 
+		      double* restrict velyplus,
+		      double* restrict velzplus, 
+		      double* restrict epsplus,
+		      double* restrict tempplus, 
+		      double* restrict yeplus,
+		      double* restrict Bvcxplus,
+		      double* restrict Bvcyplus,
+		      double* restrict Bvczplus,
+		      double* restrict psidcplus)
+{
+  DECLARE_CCTK_PARAMETERS;
+
+  double drho[nx], dpress[nx], deps[nx], d2rho[nx]; 
+  double dmrho[nx], dmeps[nx];
+  double dmtemp[nx], dmye[nx];
+  double dtemp[nx], dye[nx];
+  double dvelx[nx], dvely[nx], dvelz[nx];
+  double dmvelx[nx], dmvely[nx], dmvelz[nx];
+  double dBvcx[nx], dBvcy[nx], dBvcz[nx], dpsidc[nx];
+  double dmBvcx[nx], dmBvcy[nx], dmBvcz[nx], dmpsidc[nx];
+  double tilde_flatten[nx];
+  
+  const double onesixth = 1.0/6.0;
+
+  //  Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178
+  //  This is the expression for an even grid.
+  for(int i=1;i < nx-1;i++) {
+    drho[i]   = 0.5 * (rho[i+1]-rho[i-1]);
+    d2rho[i]  = rho[i+1] - 2.0 * rho[i] + rho[i-1];
+    dpress[i] = (press[i+1]-press[i-1]);
+    deps[i]   = 0.5 * (eps[i+1]-eps[i-1]);
+    dvelx[i]  = 0.5 * (velx[i+1]-velx[i-1]);
+    dvely[i]  = 0.5 * (vely[i+1]-vely[i-1]);
+    dvelz[i]  = 0.5 * (velz[i+1]-velz[i-1]);
+    if(do_ye) {
+      if(do_temp) {
+	dtemp[i]  = 0.5 * (temp[i+1]-temp[i-1]);
+      }
+      dye[i]    = 0.5 * (ye[i+1]-ye[i-1]);
+    }
+    if(do_mhd) {
+      dBvcx[i]  = 0.5 * (Bvcx[i+1]-Bvcx[i-1]);
+      dBvcy[i]  = 0.5 * (Bvcy[i+1]-Bvcy[i-1]);
+      dBvcz[i]  = 0.5 * (Bvcz[i+1]-Bvcz[i-1]);
+      if(dc_flag) {
+	dpsidc[i] = 0.5 * (psidc[i+1]-psidc[i-1]);
+      }
+    }
+  }
+
+  //  Steepened slope. See (1.8) of Colella and Woodward, p.178
+  for(int i=1;i<nx-1;i++) {
+    steep(rho,drho,dmrho,i);
+    steep(eps,deps,dmeps,i);
+    steep(velx,dvelx,dmvelx,i);
+    steep(vely,dvely,dmvely,i);
+    steep(velz,dvelz,dmvelz,i);
+    if(do_ye) {
+      if(do_temp) {
+	steep(temp,dtemp,dmtemp,i);
+      }
+      steep(ye,dye,dmye,i);
+    }
+    if(do_mhd) {
+      steep(Bvcx,dBvcx,dmBvcx,i);
+      steep(Bvcy,dBvcy,dmBvcy,i);
+      steep(Bvcz,dBvcz,dmBvcz,i);
+      if(dc_flag) {
+	steep(psidc,dpsidc,dmpsidc,i);
+      }
+    }
+  }
+
+  // Initial boundary states. See (1.9) of Colella and Woodward, p.178
+  for(int i=1;i<nx-2;i++) {
+    rhoplus[i] = 0.5 * (rho[i] + rho[i+1]) + 
+      (dmrho[i] - dmrho[i+1]) * onesixth;
+    epsplus[i] = 0.5 * (eps[i] + eps[i+1]) + 
+      (dmeps[i] - dmeps[i+1]) * onesixth;
+    velxplus[i] = 0.5 * (velx[i] + velx[i+1]) + 
+      (dmvelx[i] - dmvelx[i+1]) * onesixth;
+    velyplus[i] = 0.5 * (vely[i] + vely[i+1]) + 
+      (dmvely[i] - dmvely[i+1]) * onesixth;
+    velzplus[i] = 0.5 * (velz[i] + velz[i+1]) + 
+      (dmvelz[i] - dmvelz[i+1]) * onesixth;
+    if(do_ye) {
+      if(do_temp) {
+	tempplus[i] = 0.5 * (temp[i] + temp[i+1]) + 
+	  (dmtemp[i] - dmtemp[i+1]) * onesixth;
+      }
+      yeplus[i] = 0.5 * (ye[i] + ye[i+1]) + 
+	(dmye[i] - dmye[i+1]) * onesixth;
+    }
+    if(do_mhd) {
+      Bvcxplus[i] = 0.5 * (Bvcx[i] + Bvcx[i+1]) + 
+	(dmBvcx[i] - dmBvcx[i+1]) * onesixth;
+      Bvcyplus[i] = 0.5 * (Bvcy[i] + Bvcy[i+1]) + 
+	(dmBvcy[i] - dmBvcy[i+1]) * onesixth;
+      Bvczplus[i] = 0.5 * (Bvcz[i] + Bvcz[i+1]) + 
+	(dmBvcz[i] - dmBvcz[i+1]) * onesixth;
+      if(dc_flag) {
+	psidcplus[i] = 0.5 * (psidc[i] + psidc[i+1]) + 
+	  (dmpsidc[i] - dmpsidc[i+1]) * onesixth;
+      }
+    }
+  }
+
+  // fill minus states
+  for(int i=1;i<nx-2;i++) {
+    rhominus[i+1] = rhoplus[i];
+    epsminus[i+1] = epsplus[i];
+    velxminus[i+1] = velxplus[i];
+    velyminus[i+1] = velyplus[i];
+    velzminus[i+1] = velzplus[i];
+    if(do_ye) {
+      if(do_temp) {
+	tempminus[i+1] = tempplus[i];
+      }
+      yeminus[i+1] = yeplus[i];
+    }
+    if(do_mhd) {
+      Bvcxminus[i+1] = Bvcxplus[i];
+      Bvcyminus[i+1] = Bvcyplus[i];
+      Bvczminus[i+1] = Bvczplus[i];
+      if(dc_flag) {
+	psidcminus[i+1] = psidcplus[i];
+      }
+    }
+  }
+
+  /*
+    Discontinuity steepening. See (1.14-17) of C&W.
+    This is the detect routine which mat be 
+    activated with the ppm_detect parameter
+    Note that this part really also depends on the grid being even. 
+    Note also that we don''t have access to the gas constant gamma.
+    So this is just dropped from eq. (3.2) of C&W.
+    We can get around this by just rescaling the constant k0 (ppm_k0 here).
+  */
+
+  if(do_ppm_detect) {
+    for(int i=2;i<nx-2;i++) {
+      double etatilde = 0.0;
+      if ( d2rho[i+1]*d2rho[i-1] < 0.0 
+	   && ( fabs(rho[i+1]-rho[i-1]) - ppm_epsilon_shock 
+		* MIN(fabs(rho[i+1]), 
+		      fabs(rho[i-1])) > 0.0) )
+	{
+	  etatilde = (rho[i-2] - rho[i+2] + 4.0 * drho[i]) / (drho[i] * 12.0);
+	}
+      double eta = MAX(0.0, MIN(1.0, ppm_eta1 * (etatilde - ppm_eta2)));
+      if (ppm_k0 * fabs(drho[i]) * MIN(press[i-1],press[i+1]) 
+	  < fabs(dpress[i]) * MIN(rho[i-1], rho[i+1])) 
+	{
+	  eta = 0.0;
+	}
+      rhominus[i] = rhominus[i] * (1.0 - eta) + 
+	(rho[i-1] + 0.5 * dmrho[i-1]) * eta;
+      rhoplus[i] = rhoplus[i] * (1.0 - eta) + 
+	(rho[i+1] - 0.5 * dmrho[i+1]) * eta;
+    }
+  } 
+  
+  // flattening
+  for(int i=2;i<nx-2;i++) {
+    double dpress2 = press[i+2] - press[i-2];
+    double dvel = velx[i+1] - velx[i-1];
+    double w=0.0;
+    if ( (fabs(dpress[i]) >  ppm_epsilon * MIN(press[i-1],press[i+1])) 
+	 && (dvel < 0.0) ) 
+      {
+	w = 1.0;
+      } 
+    if (fabs(dpress2) < ppm_small) 
+      {
+	tilde_flatten[i] = 1.0;
+      } 
+    else
+      {
+	tilde_flatten[i] = MAX(0.0, 1.0 - w * MAX(0.0, 
+			   ppm_omega2 * (dpress[i] / dpress2 - ppm_omega1)));
+      }
+  } 
+  
+  for(int i=2;i<nx-2;i++) {
+    double flatten = tilde_flatten[i];
+    rhoplus[i] = flatten * rhoplus[i] + (1.0 - flatten) * rho[i];
+    rhominus[i] = flatten * rhominus[i] + (1.0 - flatten) * rho[i];
+    epsplus[i] = flatten * epsplus[i] + (1.0 - flatten) * eps[i];
+    epsminus[i] = flatten * epsminus[i] + (1.0 - flatten) * eps[i];
+    velxplus[i] = flatten * velxplus[i] + (1.0 - flatten) * velx[i];
+    velxminus[i] = flatten * velxminus[i] + (1.0 - flatten) * velx[i];
+    velyplus[i] = flatten * velyplus[i] + (1.0 - flatten) * vely[i];
+    velyminus[i] = flatten * velyminus[i] + (1.0 - flatten) * vely[i];
+    velzplus[i] = flatten * velzplus[i] + (1.0 - flatten) * velz[i];
+    velzminus[i] = flatten * velzminus[i] + (1.0 - flatten) * velz[i];
+
+    if(do_ye) {
+      if(do_temp) {
+	tempplus[i] = flatten * tempplus[i] + (1.0 - flatten) * temp[i];
+	tempminus[i] = flatten * tempminus[i] + (1.0 - flatten) * temp[i];
+      }
+      yeplus[i] = flatten * yeplus[i] + (1.0 - flatten) * ye[i];
+      yeminus[i] = flatten * yeminus[i] + (1.0 - flatten) * ye[i];
+    }
+
+    if(do_mhd) {
+      Bvcxplus[i] = flatten * Bvcxplus[i] + (1.0 - flatten) * Bvcx[i];
+      Bvcxminus[i] = flatten * Bvcxminus[i] + (1.0 - flatten) * Bvcx[i];
+      Bvcyplus[i] = flatten * Bvcyplus[i] + (1.0 - flatten) * Bvcy[i];
+      Bvcyminus[i] = flatten * Bvcyminus[i] + (1.0 - flatten) * Bvcy[i];
+      Bvczplus[i] = flatten * Bvczplus[i] + (1.0 - flatten) * Bvcz[i];
+      Bvczminus[i] = flatten * Bvczminus[i] + (1.0 - flatten) * Bvcz[i];
+      if(dc_flag) {
+	psidcplus[i] = flatten * psidcplus[i] + (1.0 - flatten) * psidc[i];
+	psidcminus[i] = flatten * psidcminus[i] + (1.0 - flatten) * psidc[i];
+      }
+    }
+  } // flattening
+
+  for(int i=2;i<nx-2;i++) {
+    monotonize(rhominus,rho,rhoplus,i);
+    monotonize(epsminus,eps,epsplus,i);
+    monotonize(velxminus,velx,velxplus,i);
+    monotonize(velyminus,vely,velyplus,i);
+    monotonize(velzminus,velz,velzplus,i);
+    if(do_ye) {
+      if(do_temp) {
+	monotonize(tempminus,temp,tempplus,i);
+      }
+      monotonize(yeminus,ye,yeplus,i);
+    }
+    if(do_mhd) {
+      monotonize(Bvcxminus,Bvcx,Bvcxplus,i);
+      monotonize(Bvcyminus,Bvcy,Bvcyplus,i);
+      monotonize(Bvczminus,Bvcz,Bvczplus,i);
+      if(dc_flag) {
+	monotonize(psidcminus,psidc,psidcplus,i);
+      }
+    }
+  }
+
+  return;
+}
+
+//template template<bool do_temp, bool do_ye, bool do_mhd, bool dc_flag, bool do_ppm_detect>
+#define INSTANTIATE(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect)			\
+template static void ppm1d_cxx<do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect>(const int nx, \
+		      const double dx,					\
+		      const double* restrict rho,			\
+		      const double* restrict velx,			\
+		      const double* restrict vely,			\
+		      const double* restrict velz,			\
+		      const double* restrict eps,			\
+		      const double* restrict press,			\
+		      const double* restrict temp,			\
+		      const double* restrict ye,			\
+		      const double* restrict Bvcx,			\
+		      const double* restrict Bvcy,			\
+		      const double* restrict Bvcz,			\
+		      const double* restrict psidc,			\
+		      double* restrict rhominus,			\
+		      double* restrict velxminus,			\
+		      double* restrict velyminus,			\
+		      double* restrict velzminus,			\
+		      double* restrict epsminus,			\
+		      double* restrict tempminus,			\
+		      double* restrict yeminus,				\
+		      double* restrict Bvcxminus,			\
+		      double* restrict Bvcyminus,			\
+		      double* restrict Bvczminus,			\
+		      double* restrict psidcminus,			\
+		      double* restrict rhoplus,				\
+		      double* restrict velxplus,			\
+		      double* restrict velyplus,			\
+		      double* restrict velzplus,			\
+		      double* restrict epsplus,				\
+		      double* restrict tempplus,			\
+		      double* restrict yeplus,				\
+		      double* restrict Bvcxplus,			\
+		      double* restrict Bvcyplus,			\
+		      double* restrict Bvczplus,			\
+		      double* restrict psidcplus);
+
+//template template<bool do_temp, bool do_ye, bool do_mhd, bool dc_flag, bool do_ppm_detect>
+
+
+// first stuff without MHD
+INSTANTIATE(false,false,false,false,false)
+INSTANTIATE(false,false,false,false,true)
+// do_temp can only be true of do_ye is also true
+INSTANTIATE(true,true,false,false,false)
+INSTANTIATE(true,true,false,false,true)
+// but do_ye can be true and do_temp can be false
+INSTANTIATE(false,true,false,false,false)
+INSTANTIATE(false,true,false,false,true)
+
+// with MHD, but dc_flag false
+INSTANTIATE(false,false,true,false,false)
+INSTANTIATE(false,false,true,false,true)
+// do_temp can only be true of do_ye is also true
+INSTANTIATE(true,true,true,false,false)
+INSTANTIATE(true,true,true,false,true)
+// but do_ye can be true and do_temp can be false
+INSTANTIATE(false,true,true,false,false)
+INSTANTIATE(false,true,true,false,true)
+
+// with MHD, dc_flag true
+INSTANTIATE(false,false,true,true,false)
+INSTANTIATE(false,false,true,true,true)
+// do_temp can only be true of do_ye is also true
+INSTANTIATE(true,true,true,true,false)
+INSTANTIATE(true,true,true,true,true)
+// but do_ye can be true and do_temp can be false
+INSTANTIATE(false,true,true,true,false)
+INSTANTIATE(false,true,true,true,true)
diff --git a/src/make.code.defn b/src/make.code.defn
index 527713e..21181df 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -82,6 +82,7 @@ SRCS = 	Utils.F90 \
         Con2Prim_fortran_interfaces.F90 \
         Con2PrimM_fortran_interfaces.F90 \
         GRHydro_Prim2Con.c \
+	GRHydro_PPM.cc \
         GRHydro_Wrappers.F90
 	
 # Subdirectories containing source files