[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 623)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Tue Apr 15 14:49:23 CDT 2014
User: rhaas
Date: 2014/04/15 02:49 PM
Added:
/trunk/src/
GRHydro_MP5Reconstruct.cc, GRHydro_MP5Reconstruct.hh
Modified:
/trunk/
param.ccl
/trunk/src/
make.code.defn
Log:
GRHydro: add MP5 reconstruct
File Changes:
Directory: /trunk/src/
======================
File [added]: GRHydro_MP5Reconstruct.cc
Delta lines: +114 -0
===================================================================
--- trunk/src/GRHydro_MP5Reconstruct.cc (rev 0)
+++ trunk/src/GRHydro_MP5Reconstruct.cc 2014-04-15 19:49:22 UTC (rev 623)
@@ -0,0 +1,114 @@
+#include <cmath>
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+#include "GRHydro_Reconstruct_drv_impl.hh"
+#include "GRHydro_MP5Reconstruct.hh"
+
+using namespace std;
+
+template <typename T> static inline T SQR (T const & x) { return x*x; }
+
+#define MIN(a,b) (((a)<(b))?(a):(b))
+#define MAX(a,b) (((a)>(b))?(a):(b))
+
+#define MIN3(a,b,c) (MIN(a, MIN(b,c)))
+#define MIN4(a,b,c,d) (MIN(a, MIN(b,MIN(c,d))))
+#define MAX3(a,b,c) (MAX(a, MAX(b,c)))
+
+#define MINMOD(x,y) \
+ (0.5*(copysign(1.0,x) + copysign(1.0,y)) * MIN(fabs(x), fabs(y)))
+
+#define MINMOD4(w,x,y,z) \
+ (0.125*( copysign(1.0,w)+copysign(1.0,x) )*fabs( (copysign(1.0,w)+copysign(1.0,y)) * (copysign(1.0,w)+copysign(1.0,z)) )*MIN4(fabs(w), fabs(x), fabs(y), fabs(z)))
+
+
+static inline CCTK_REAL MP5(const CCTK_REAL am2,
+ const CCTK_REAL am1,
+ const CCTK_REAL a,
+ const CCTK_REAL ap1,
+ const CCTK_REAL ap2,
+ const CCTK_REAL anorm,
+ const CCTK_REAL mp5_eps,
+ const CCTK_REAL mp5_alpha
+ )
+{
+ const CCTK_REAL vl = (2.0*am2 - 13.0*am1 + 47.0*a + 27.0*ap1 - 3.0*ap2) * 1.0/60.0;
+ const CCTK_REAL vmp = a + MINMOD( ap1-a, mp5_alpha*(a-am1) );
+ if ((vl-a)*(vl-vmp) <= mp5_eps*anorm)
+ return vl;
+ else {
+ const CCTK_REAL djm1 = am2 -2.0*am1 + a;
+ const CCTK_REAL dj = am1 -2.0*a + ap1;
+ const CCTK_REAL djp1 = a -2.0*ap1 + ap2;
+ const CCTK_REAL dm4jph = MINMOD4(4.0*dj-djp1, 4.0*djp1-dj, dj, djp1);
+ const CCTK_REAL dm4jmh = MINMOD4(4.0*dj-djm1, 4.0*djm1-dj, dj, djm1);
+ const CCTK_REAL vul = a + mp5_alpha*(a-am1);
+ const CCTK_REAL vav = 0.5*(a+ap1);
+ const CCTK_REAL vmd = vav - 0.5*dm4jph;
+ const CCTK_REAL vlc = a + 0.5*(a-am1) + 4.0/3.0*dm4jmh;
+ const CCTK_REAL vmin = MAX(MIN3(a,ap1,vmd), MIN3(a,vul,vlc));
+ const CCTK_REAL vmax = MIN(MAX3(a,ap1,vmd), MAX3(a,vul,vlc));
+ return vl + MINMOD(vmin-vl, vmax-vl);
+ }
+ return 0;
+}
+
+
+/**
+ MP5 reconstruction operator.
+*/
+template <bool do_MP5_adaptive_epsilon>
+template <int dir>
+inline void
+GRHydro_MP5Reconstruct1d_cxx<do_MP5_adaptive_epsilon>::
+ apply(const int nx,
+ const CCTK_REAL* const restrict a,
+ CCTK_REAL* const restrict aminus,
+ CCTK_REAL* const restrict aplus,
+ const cGH* const cctkGH,
+ const int j, const int k
+ )
+{
+#define A(i_) (a[ijk[i_]])
+#define Aplus(i_) (aplus[ijk[i_]])
+#define Aminus(i_) (aminus[ijk[i_]])
+
+ DECLARE_CCTK_PARAMETERS;
+
+ for (int i=GRHydro_stencil-1; i < nx-GRHydro_stencil+1; ++i)
+ {
+ const int ijk[5] = {
+ dir ==0 ? CCTK_GFINDEX3D(cctkGH, i-2, j, k) : dir ==1 ? CCTK_GFINDEX3D(cctkGH, j, i-2, k) : CCTK_GFINDEX3D(cctkGH, j, k, i-2),
+ dir ==0 ? CCTK_GFINDEX3D(cctkGH, i-1, j, k) : dir ==1 ? CCTK_GFINDEX3D(cctkGH, j, i-1, k) : CCTK_GFINDEX3D(cctkGH, j, k, i-1),
+ dir ==0 ? CCTK_GFINDEX3D(cctkGH, i , j, k) : dir ==1 ? CCTK_GFINDEX3D(cctkGH, j, i , k) : CCTK_GFINDEX3D(cctkGH, j, k, i ),
+ dir ==0 ? CCTK_GFINDEX3D(cctkGH, i+1, j, k) : dir ==1 ? CCTK_GFINDEX3D(cctkGH, j, i+1, k) : CCTK_GFINDEX3D(cctkGH, j, k, i+1),
+ dir ==0 ? CCTK_GFINDEX3D(cctkGH, i+2, j, k) : dir ==1 ? CCTK_GFINDEX3D(cctkGH, j, i+2, k) : CCTK_GFINDEX3D(cctkGH, j, k, i+2)
+ };
+
+ if (!do_MP5_adaptive_epsilon) {
+
+ Aplus(2) = MP5(A(0), A(1), A(2), A(3), A(4), 1.0, mp5_eps, mp5_alpha);
+ Aminus(2) = MP5(A(4), A(3), A(2), A(1), A(0), 1.0, mp5_eps, mp5_alpha);
+
+ } else {
+
+ const CCTK_REAL anorm = sqrt(SQR(A(0)) + SQR(A(1)) + SQR(A(2)) + SQR(A(3)) + SQR(A(4)));
+
+ Aplus(2) = MP5(A(0), A(1), A(2), A(3), A(4), anorm, mp5_eps, mp5_alpha);
+ Aminus(2) = MP5(A(4), A(3), A(2), A(1), A(0), anorm, mp5_eps, mp5_alpha);
+
+ }
+
+ }
+}
+
+// instantiate all copies we need, this way different operators can be compiled
+// in parallel. This must match the select routine in GRHydro_Reconstruct.cc
+template class GRHydro_MP5Reconstruct1d_cxx<false>;
+template class GRHydro_MP5Reconstruct1d_cxx<true>;
+
+INSTANTIATE_RECONSTRUCTION_OPERATOR(GRHydro_MP5Reconstruct1d_cxx<false>)
+INSTANTIATE_RECONSTRUCTION_OPERATOR(GRHydro_MP5Reconstruct1d_cxx<true >)
File [added]: GRHydro_MP5Reconstruct.hh
Delta lines: +25 -0
===================================================================
--- trunk/src/GRHydro_MP5Reconstruct.hh (rev 0)
+++ trunk/src/GRHydro_MP5Reconstruct.hh 2014-04-15 19:49:22 UTC (rev 623)
@@ -0,0 +1,25 @@
+#ifndef _GRHYDRO_MP5RECONSTRUCT_HH
+#define _GRHYDRO_MP5RECONSTRUCT_HH
+
+#include "cctk.h"
+
+using namespace std;
+
+
+/**
+ MP5 reconstruction operator.
+*/
+template <bool do_MP5_adaptive_epsilon>
+struct GRHydro_MP5Reconstruct1d_cxx {
+
+template <int dir>
+static inline void apply(const int nx,
+ const CCTK_REAL* const restrict a,
+ CCTK_REAL* const restrict aminus,
+ CCTK_REAL* const restrict aplus,
+ const cGH* const cctkGH,
+ const int j, const int k
+ );
+}; // end struct
+
+#endif // _GRHYDRO_MP5RECONSTRUCT_HH
File [modified]: make.code.defn
Delta lines: +1 -0
===================================================================
--- trunk/src/make.code.defn 2014-04-15 19:49:20 UTC (rev 622)
+++ trunk/src/make.code.defn 2014-04-15 19:49:22 UTC (rev 623)
@@ -87,6 +87,7 @@
GRHydro_PPM.cc \
GRHydro_ePPM.cc \
GRHydro_PPMReconstruct_drv_opt.cc \
+ GRHydro_MP5Reconstruct.cc \
GRHydro_Wrappers.F90 \
GRHydro_HLLE.cc
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +3 -0
===================================================================
--- trunk/param.ccl 2014-04-15 19:49:20 UTC (rev 622)
+++ trunk/param.ccl 2014-04-15 19:49:22 UTC (rev 623)
@@ -252,6 +252,9 @@
0:* :: "0.0 or very small and positive. 1e-10 is suggested by Suresh&Huynh. TOV star requires 0.0"
} 0.0
+boolean mp5_adaptive_eps "Same as WENO adaptive epsilon: adaptively reduce mp5_eps according to norm of stencil. Original algorithm does not use this."
+{
+} no
boolean use_enhanced_ppm "Use the enhanced ppm reconstruction method proposed by Colella & Sekora 2008 and McCorquodale & Colella 2011" STEERABLE=RECOVER
{
More information about the Commits
mailing list