GRHydro: reimplement PPM in C

* add prototype of a re-implementation of PPM in C. So far
  we can do old PPM without and with temperature/Y_e reconstruction.
* also add a version of the reconstruction driver that is optimized
  for the use of the new PPM routine.
* no change to standard code version; one needs to manually replace
  source files to get the new stuff.

* new version is factor 3 faster than old version (based on preliminary
  measurements).

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

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

1 files changed, 56 insertions, 0 deletions

diff --git a/src/GRHydro_PPM.h b/src/GRHydro_PPM.h
new file mode 100644
index 0000000..dcc63ad
--- /dev/null
+++ b/src/GRHydro_PPM.h
@@ -0,0 +1,56 @@
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#if 1
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#endif
+
+#define MIN(a,b) (((a)<(b))?(a):(b))
+#define MAX(a,b) (((a)>(b))?(a):(b))
+
+static inline void steep(double *x, double *dx, double* 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,
+			      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;
+}	     	    
+
+
+