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

rhaas at tapir.caltech.edu rhaas at tapir.caltech.edu
Tue Apr 15 14:49:25 CDT 2014


User: rhaas
Date: 2014/04/15 02:49 PM

Added:
 /trunk/src/
  GRHydro_WENOReconstruct.cc, GRHydro_WENOReconstruct.hh

Modified:
 /trunk/
  param.ccl
 /trunk/src/
  GRHydro_HLLE.cc, make.code.defn

Log:
 GRHydro: add WENO reconstruct

File Changes:

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

File [modified]: GRHydro_HLLE.cc
Delta lines: +3 -6
===================================================================
--- trunk/src/GRHydro_HLLE.cc	2014-04-15 19:49:22 UTC (rev 623)
+++ trunk/src/GRHydro_HLLE.cc	2014-04-15 19:49:25 UTC (rev 624)
@@ -568,8 +568,7 @@
       	          	  velxminus[idxr],velyminus[idxr],velzminus[idxr],
       	          	  pressminus[idxr],savg_detr,
                           avg_alp,avg_betax,avg_betay,avg_betaz,blowxm,blowym,blowzm,b2m,ab0m,wm);
-      	    }
-      	    else if(fdir==2) {
+            } else if(fdir==2) {
       	      eigenvaluesM_cc(rhominus[idxr],
       	          	   velyminus[idxr],
       	          	   velzminus[idxr],
@@ -648,8 +647,7 @@
       	          	  pressminus[idxr],savg_detr,
                           avg_alp,avg_betaz,avg_betax,avg_betay,blowzm,blowxm,blowym,b2m,ab0m,wm);
       	    } 
-          }
-          else {
+          } else {
   	    
       	    if (fdir==1) {
       	      eigenvalues_cc(rhominus[idxr],
@@ -682,8 +680,7 @@
       	          	  fminus[0],fminus[1],fminus[2],fminus[3],fminus[4],
       	          	  velxminus[idxr],pressminus[idxr],savg_detr,
       	          	  avg_alp,avg_beta);	    
-      	    }
-      	    else if(fdir==2) {
+            } else if(fdir==2) {
       	      eigenvalues_cc(rhominus[idxr],
       	          	   velyminus[idxr],
       	          	   velzminus[idxr],

File [added]: GRHydro_WENOReconstruct.cc
Delta lines: +195 -0
===================================================================
--- trunk/src/GRHydro_WENOReconstruct.cc	                        (rev 0)
+++ trunk/src/GRHydro_WENOReconstruct.cc	2014-04-15 19:49:25 UTC (rev 624)
@@ -0,0 +1,195 @@
+#include <cmath>
+#include <iostream>
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+#include "GRHydro_Reconstruct_drv_impl.hh"
+#include "GRHydro_WENOReconstruct.hh"
+
+using namespace std;
+
+template <typename T> static inline T SQR (T const & x) { return x*x; }
+
+/**
+   WENO5 reconstruction operator.
+   Supports standard WENO5 (with and without adaptive epsilon), and WENO-z.
+*/
+template <bool do_wenoz, bool do_adaptive_epsilon>
+template <int dir>
+void GRHydro_WENOReconstruct1d_cxx<do_wenoz,do_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)
+{
+   DECLARE_CCTK_PARAMETERS;
+   
+#define A(i_) (a[ijk[i_]])
+#define Aplus(i_) (aplus[ijk[i_]])
+#define Aminus(i_) (aminus[ijk[i_]])
+   
+   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)
+                         };
+                       
+   
+         
+   
+      assert(! (do_wenoz && do_adaptive_epsilon) && "Adaptive_epsilon not supported for WENO-Z");
+
+      if (do_wenoz)
+      {
+         static const CCTK_REAL 
+                         weno_coeffs[3][5] = { { 2.0/6.0, -7.0/6.0, 11.0/6.0, 0,        0 }, 
+                                               { 0,       -1.0/6.0, 5.0/6.0,  2.0/6.0,  0 },
+                                               { 0,        0,       2.0/6.0,  5.0/6.0, -1.0/6.0 } };
+      
+         const CCTK_REAL beta1 = 13.0/12.0*SQR(A(0)-2.0*A(1)+A(2)) + 1.0/4.0*SQR(A(0)-4.0*A(1)+3.0*A(2));
+         const CCTK_REAL beta2 = 13.0/12.0*SQR(A(1)-2.0*A(2)+A(3)) + 1.0/4.0*SQR(A(1)-A(3));
+         const CCTK_REAL beta3 = 13.0/12.0*SQR(A(2)-2.0*A(3)+A(4)) + 1.0/4.0*SQR(3.0*A(2)-4.0*A(3)+A(4));
+            
+            
+         //    Compute weights according to weno-z alorithm
+         const CCTK_REAL wbarplus1 = 1.0/10.0 * (1.0 + abs(beta1-beta3) / (weno_eps + beta1));
+         const CCTK_REAL wbarplus2 = 3.0/5.0 * (1.0 + abs(beta1-beta3) / (weno_eps + beta2));
+         const CCTK_REAL wbarplus3 = 3.0/10.0 * (1.0 + abs(beta1-beta3) / (weno_eps + beta3));
+
+         const CCTK_REAL wbarminus1 = 3.0/10.0 * (1.0 + abs(beta1-beta3) / (weno_eps + beta1));
+         const CCTK_REAL wbarminus2 = 3.0/5.0 * (1.0 + abs(beta1-beta3) / (weno_eps + beta2));
+         const CCTK_REAL wbarminus3 = 1.0/10.0 * (1.0 + abs(beta1-beta3) / (weno_eps + beta3));
+         
+         const CCTK_REAL iwbarplussum = 1.0 / (wbarplus1 + wbarplus2 + wbarplus3);
+         
+         const CCTK_REAL wplus1 = wbarplus1 * iwbarplussum;
+         const CCTK_REAL wplus2 = wbarplus2 * iwbarplussum;
+         const CCTK_REAL wplus3 = wbarplus3 * iwbarplussum;
+         
+         const CCTK_REAL iwbarminussum = 1.0 / (wbarminus1 + wbarminus2 + wbarminus3);
+         
+         const CCTK_REAL wminus1 = wbarminus1 * iwbarminussum;
+         const CCTK_REAL wminus2 = wbarminus2 * iwbarminussum;
+         const CCTK_REAL wminus3 = wbarminus3 * iwbarminussum;
+         
+         //    Calculate the reconstruction
+         Aplus(2) = 0;
+         Aminus(2) = 0;
+         for (int j=0; j < 5; ++j) {
+               Aplus(2) += (wplus1 * weno_coeffs[0][j]
+                          + wplus2 * weno_coeffs[1][j]
+                          + wplus3 * weno_coeffs[2][j]) * A(j);
+               Aminus(2) += (wminus1 * weno_coeffs[2][4-j]
+                           + wminus2 * weno_coeffs[1][4-j]
+                           + wminus3 * weno_coeffs[0][4-j]) * A(j);
+         }
+      } else {
+         
+         static const CCTK_REAL beta_shu[3][6] = { { 4.0/3.0,  -19.0/3.0, 25.0/3.0, 11.0/3.0, -31.0/3.0, 10.0/3.0 },
+                                      { 4.0/3.0,  -13.0/3.0, 13.0/3.0, 5.0/3.0,  -13.0/3.0, 4.0/3.0 },
+                                      { 10.0/3.0, -31.0/3.0, 25.0/3.0, 11.0/3.0, -19.0/3.0, 4.0/3.0 } };
+         static const CCTK_REAL weno_coeffs[3][5] = { { 3.0/8.0, -5.0/4.0, 15.0/8.0, 0,      0 },
+                                         { 0,       -1.0/8.0, 3.0/4.0,  3.0/8.0, 0 },
+                                         { 0,       0,        3.0/8.0,  3.0/4.0, -1.0/8.0 } };
+                                      
+         //    Compute smoothness indicators
+         //    This is from Tchekhovskoy et al 2007 (WHAM code paper).
+         CCTK_REAL beta1  = beta_shu[0][0]*SQR(A(0))
+                  + beta_shu[0][1]*A(0)*A(1)
+                  + beta_shu[0][2]*SQR(A(1))
+                  + beta_shu[0][3]*A(0)*A(2)
+                  + beta_shu[0][4]*A(1)*A(2)
+                  + beta_shu[0][5]*SQR(A(2));
+         
+         CCTK_REAL beta2  = beta_shu[1][0]*SQR(A(1))
+                  + beta_shu[1][1]*A(1)*A(2)
+                  + beta_shu[1][2]*SQR(A(2))
+                  + beta_shu[1][3]*A(1)*A(3)
+                  + beta_shu[1][4]*A(2)*A(3)
+                  + beta_shu[1][5]*SQR(A(3));
+         
+         CCTK_REAL beta3  = beta_shu[2][0]*SQR(A(2))
+                  + beta_shu[2][1]*A(2)*A(3)
+                  + beta_shu[2][2]*SQR(A(3))
+                  + beta_shu[2][3]*A(2)*A(4)
+                  + beta_shu[2][4]*A(3)*A(4)
+                  + beta_shu[2][5]*SQR(A(4));
+         
+         
+         if (do_adaptive_epsilon) {
+            const CCTK_REAL vnorm = (SQR(A(0)) + SQR(A(1)) + SQR(A(2)) + SQR(A(3)) + SQR(A(4)));
+               
+            beta1 += 100.0*weno_eps*(vnorm + 1.0);
+            beta2 += 100.0*weno_eps*(vnorm + 1.0);
+            beta3 += 100.0*weno_eps*(vnorm + 1.0);
+               
+            const CCTK_REAL ibetanorm = 1.0 / (beta1 + beta2 + beta3);
+               
+            beta1 *= ibetanorm;
+            beta2 *= ibetanorm;
+            beta3 *= ibetanorm;
+         }
+         
+         const CCTK_REAL wbarplus1 = 1.0/16.0 / SQR(weno_eps + beta1);
+         const CCTK_REAL wbarplus2 = 5.0/8.0 / SQR(weno_eps + beta2);
+         const CCTK_REAL wbarplus3 = 5.0/16.0 / SQR(weno_eps + beta3);
+         
+         const CCTK_REAL iwbarplussum = 1.0 / (wbarplus1 + wbarplus2 + wbarplus3);
+         
+         const CCTK_REAL wplus1 = wbarplus1 * iwbarplussum;
+         const CCTK_REAL wplus2 = wbarplus2 * iwbarplussum;
+         const CCTK_REAL wplus3 = wbarplus3 * iwbarplussum;
+
+         const CCTK_REAL wbarminus1 = 5.0/16.0 / SQR(weno_eps + beta1);
+         const CCTK_REAL wbarminus2 = 5.0/8.0 / SQR(weno_eps + beta2);
+         const CCTK_REAL wbarminus3 = 1.0/16.0 / SQR(weno_eps + beta3);
+         
+         const CCTK_REAL iwbarminussum = 1.0 / (wbarminus1 + wbarminus2 + wbarminus3);
+         
+         const CCTK_REAL wminus1 = wbarminus1 * iwbarminussum;
+         const CCTK_REAL wminus2 = wbarminus2 * iwbarminussum;
+         const CCTK_REAL wminus3 = wbarminus3 * iwbarminussum;
+                                         
+         //    Calculate the reconstruction
+         Aplus(2) = 0;
+         Aminus(2) = 0;
+         for (int j=0; j < 5; ++j) {
+               Aplus(2) += (wplus1 * weno_coeffs[0][j]
+                          + wplus2 * weno_coeffs[1][j]
+                          + wplus3 * weno_coeffs[2][j]) * A(j);
+               Aminus(2) += (wminus1 * weno_coeffs[2][4-j]
+                           + wminus2 * weno_coeffs[1][4-j]
+                           + wminus3 * weno_coeffs[0][4-j]) * A(j);
+         }
+      }
+   }
+}
+
+// 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_WENOReconstruct1d_cxx<false,false>;
+template class GRHydro_WENOReconstruct1d_cxx<false,true >;
+template class GRHydro_WENOReconstruct1d_cxx<true ,false>;
+
+// the typedefs are required since the preprocessor chokes on the comma in the
+// template specialization. Alternative:
+// #define COMMA ,
+// or #define / #undef a la Carpet's typecase.hh
+typedef GRHydro_WENOReconstruct1d_cxx<false,false>
+  GRHydro_WENOReconstruct1d_cxx00;
+typedef GRHydro_WENOReconstruct1d_cxx<false,true >
+  GRHydro_WENOReconstruct1d_cxx01;
+typedef GRHydro_WENOReconstruct1d_cxx<true ,false>
+  GRHydro_WENOReconstruct1d_cxx10;
+typedef GRHydro_WENOReconstruct1d_cxx<true ,true >
+  GRHydro_WENOReconstruct1d_cxx11;
+
+INSTANTIATE_RECONSTRUCTION_OPERATOR(GRHydro_WENOReconstruct1d_cxx00)
+INSTANTIATE_RECONSTRUCTION_OPERATOR(GRHydro_WENOReconstruct1d_cxx01)
+INSTANTIATE_RECONSTRUCTION_OPERATOR(GRHydro_WENOReconstruct1d_cxx10)
+INSTANTIATE_RECONSTRUCTION_OPERATOR(GRHydro_WENOReconstruct1d_cxx11)

File [added]: GRHydro_WENOReconstruct.hh
Delta lines: +25 -0
===================================================================
--- trunk/src/GRHydro_WENOReconstruct.hh	                        (rev 0)
+++ trunk/src/GRHydro_WENOReconstruct.hh	2014-04-15 19:49:25 UTC (rev 624)
@@ -0,0 +1,25 @@
+#ifndef _GRHYDRO_WENORECONSTRUCT_HH
+#define _GRHYDRO_WENORECONSTRUCT_HH
+
+#include "cctk.h"
+
+using namespace std;
+
+/**
+   WENO5 reconstruction operator.
+   Supports standard WENO5 (with and without adaptive epsilon), and WENO-z.
+*/
+template <bool do_wenoz, bool do_adaptive_epsilon>
+struct GRHydro_WENOReconstruct1d_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_WENORECONSTRUCT_HH

File [modified]: make.code.defn
Delta lines: +1 -0
===================================================================
--- trunk/src/make.code.defn	2014-04-15 19:49:22 UTC (rev 623)
+++ trunk/src/make.code.defn	2014-04-15 19:49:25 UTC (rev 624)
@@ -87,6 +87,7 @@
 	GRHydro_PPM.cc \
 	GRHydro_ePPM.cc \
         GRHydro_PPMReconstruct_drv_opt.cc \
+        GRHydro_WENOReconstruct.cc \
         GRHydro_MP5Reconstruct.cc \
         GRHydro_Wrappers.F90 \
 	GRHydro_HLLE.cc

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

File [modified]: param.ccl
Delta lines: +1 -1
===================================================================
--- trunk/param.ccl	2014-04-15 19:49:22 UTC (rev 623)
+++ trunk/param.ccl	2014-04-15 19:49:25 UTC (rev 624)
@@ -289,7 +289,7 @@
   5  :: "Fifth-order"
 } 5
 
-real weno_eps "WENO epsilon parameter"
+real weno_eps "WENO epsilon parameter. For WENO-z, 1e-40 is recommended" STEERABLE=ALWAYS
 {
   0:*  :: "small and positive"
 } 1e-26



More information about the Commits mailing list