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