From e704df2771e72ccdf9046967bb95545ddd735e67 Mon Sep 17 00:00:00 2001 From: rhaas Date: Tue, 15 Apr 2014 19:48:44 +0000 Subject: GRHydro: New C++ ePPM algorithm Tested with TOV star. Improved comments and pointed out a difference to F90 code which makes the C++ code more accurate. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@609 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_PPMReconstruct_drv_opt.cc | 107 ++++++---- src/GRHydro_PPM_opt.h | 41 ++++ src/GRHydro_ePPM.cc | 367 ++++++++++++++++++++++++++++++++++ src/make.code.defn | 1 + 4 files changed, 480 insertions(+), 36 deletions(-) create mode 100644 src/GRHydro_ePPM.cc diff --git a/src/GRHydro_PPMReconstruct_drv_opt.cc b/src/GRHydro_PPMReconstruct_drv_opt.cc index 9b70d2c..7ed18de 100644 --- a/src/GRHydro_PPMReconstruct_drv_opt.cc +++ b/src/GRHydro_PPMReconstruct_drv_opt.cc @@ -89,15 +89,13 @@ void CCTK_FNAME(GRHydro_PPMReconstruct_drv_opt)(cGH const * const restrict & res // trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","trivial"); not_trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","not_trivial"); - // ENHANCED PPM IS NOT IMPLEMENTED - // logical :: apply_enhanced_ppm - // ! if use_enhanced_ppm, allow old PPM on one level - // if (GRHydro_oppm_reflevel .eq. (-1) .or. & - // GRHydro_reflevel .ne. GRHydro_oppm_reflevel) then - // apply_enhanced_ppm = use_enhanced_ppm .ne. 0 - // else - // apply_enhanced_ppm = .false. - // end if + // if use_enhanced_ppm, allow old PPM on one level + bool apply_enhanced_ppm = false; + if (GRHydro_oppm_reflevel == (-1) || + *GRHydro_reflevel != GRHydro_oppm_reflevel) + apply_enhanced_ppm = use_enhanced_ppm != 0; + else + apply_enhanced_ppm = false; int i,j,k; @@ -178,54 +176,91 @@ void CCTK_FNAME(GRHydro_PPMReconstruct_drv_opt)(cGH const * const restrict & res double* restrict Bvcyplus, \ double* restrict Bvczplus, \ double* restrict psidcplus); - static ppm1d_cxx_ptr_t ppm1d_cxx_funcs[32]; + static ppm1d_cxx_ptr_t ppm1d_cxx_funcs[64]; ppm1d_cxx_ptr_t ppm1d_cxx_func; - #define MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect) ((do_temp<<4)|(do_ye<<3)|(do_mhd<<2)|(dc_flag<<1)|(do_ppm_detect)) - #define SETPOINTER(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect) \ - ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect)] = \ - GRHydro_ppm1d_cxx - + #define MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect,do_eppm) ((do_eppm<<5)|(do_temp<<4)|(do_ye<<3)|(do_mhd<<2)|(dc_flag<<1)|(do_ppm_detect)) + #define SETPOINTER(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect,do_eppm) \ + if (!do_eppm) \ + ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect, do_eppm)] = \ + GRHydro_ppm1d_cxx; \ + else \ + ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect, do_eppm)] = \ + GRHydro_eppm1d_cxx; + + // --- original PPM // first stuff without MHD - SETPOINTER(false,false,false,false,false); - SETPOINTER(false,false,false,false,true); + SETPOINTER(false,false,false,false,false,false); + SETPOINTER(false,false,false,false,true,false); // do_temp can only be true of do_ye is also true - SETPOINTER(true,true,false,false,false); - SETPOINTER(true,true,false,false,true); + SETPOINTER(true,true,false,false,false,false); + SETPOINTER(true,true,false,false,true,false); // but do_ye can be true and do_temp can be false - SETPOINTER(false,true,false,false,false); - SETPOINTER(false,true,false,false,true); + SETPOINTER(false,true,false,false,false,false); + SETPOINTER(false,true,false,false,true,false); // with MHD, but dc_flag false - SETPOINTER(false,false,true,false,false); - SETPOINTER(false,false,true,false,true); + SETPOINTER(false,false,true,false,false,false); + SETPOINTER(false,false,true,false,true,false); // do_temp can only be true of do_ye is also true - SETPOINTER(true,true,true,false,false); - SETPOINTER(true,true,true,false,true); + SETPOINTER(true,true,true,false,false,false); + SETPOINTER(true,true,true,false,true,false); // but do_ye can be true and do_temp can be false - SETPOINTER(false,true,true,false,false); - SETPOINTER(false,true,true,false,true); + SETPOINTER(false,true,true,false,false,false); + SETPOINTER(false,true,true,false,true,false); // with MHD, dc_flag true - SETPOINTER(false,false,true,true,false); - SETPOINTER(false,false,true,true,true); + SETPOINTER(false,false,true,true,false,false); + SETPOINTER(false,false,true,true,true,false); + // do_temp can only be true of do_ye is also true + SETPOINTER(true,true,true,true,false,false); + SETPOINTER(true,true,true,true,true,false); + // but do_ye can be true and do_temp can be false + SETPOINTER(false,true,true,true,false,false); + SETPOINTER(false,true,true,true,true,false); + + // ---- enhanced PPM + // first stuff without MHD + SETPOINTER(false,false,false,false,false,true); + SETPOINTER(false,false,false,false,true,true); + // do_temp can only be true of do_ye is also true + SETPOINTER(true,true,false,false,false,true); + SETPOINTER(true,true,false,false,true,true); + // but do_ye can be true and do_temp can be false + SETPOINTER(false,true,false,false,false,true); + SETPOINTER(false,true,false,false,true,true); + + // with MHD, but dc_flag false + SETPOINTER(false,false,true,false,false,true); + SETPOINTER(false,false,true,false,true,true); // do_temp can only be true of do_ye is also true - SETPOINTER(true,true,true,true,false); - SETPOINTER(true,true,true,true,true); + SETPOINTER(true,true,true,false,false,true); + SETPOINTER(true,true,true,false,true,true); // but do_ye can be true and do_temp can be false - SETPOINTER(false,true,true,true,false); - SETPOINTER(false,true,true,true,true); + SETPOINTER(false,true,true,false,false,true); + SETPOINTER(false,true,true,false,true,true); + // with MHD, dc_flag true + SETPOINTER(false,false,true,true,false,true); + SETPOINTER(false,false,true,true,true,true); + // do_temp can only be true of do_ye is also true + SETPOINTER(true,true,true,true,false,true); + SETPOINTER(true,true,true,true,true,true); + // but do_ye can be true and do_temp can be false + SETPOINTER(false,true,true,true,false,true); + SETPOINTER(false,true,true,true,true,true); + + const bool do_temp = *evolve_temper; const bool do_ye = *evolve_Y_e; const bool do_mhd = *evolve_MHD; const bool dc_flag = clean_divergence; const bool do_ppm_detect = ppm_detect; - ppm1d_cxx_func = ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect)]; + ppm1d_cxx_func = ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect,apply_enhanced_ppm)]; if(ppm1d_cxx_func == NULL) { CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, - "Combination do_temp = %d do_ye = %d do_mhd = %d dc_flag = %d do_ppm_detect = %d is not supported for PPM. Please instatate the respective template", - do_temp, do_ye, do_mhd, dc_flag, do_ppm_detect); + "Combination do_temp = %d do_ye = %d do_mhd = %d dc_flag = %d do_ppm_detect = %d do_eppm = %d is not supported for PPM. Please instantiate the respective template", + do_temp, do_ye, do_mhd, dc_flag, do_ppm_detect, apply_enhanced_ppm); } //index = i + j*nx + k*nx*ny; diff --git a/src/GRHydro_PPM_opt.h b/src/GRHydro_PPM_opt.h index 5107723..2a7f723 100644 --- a/src/GRHydro_PPM_opt.h +++ b/src/GRHydro_PPM_opt.h @@ -36,3 +36,44 @@ void GRHydro_ppm1d_cxx(const int nx, double* restrict Bvcyplus, double* restrict Bvczplus, double* restrict psidcplus); + +template +void GRHydro_eppm1d_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); + + diff --git a/src/GRHydro_ePPM.cc b/src/GRHydro_ePPM.cc new file mode 100644 index 0000000..8e5e7de --- /dev/null +++ b/src/GRHydro_ePPM.cc @@ -0,0 +1,367 @@ +#include +#include +#include +#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)) + +#define MIN3(a,b,c) (MIN(a, MIN(b, c))) +#define MIN4(a,b,c,d) (MIN(a, MIN(b, MIN(c, d)))) + +#define MAX5(a,b,c,d,e) (MAX(a, MAX(b, MAX(c, MAX(d, e))))) + +using namespace std; + + +/** + ePPM Reconstruction (see Reisswig et al. 2013) based on McCorquodale & Colella (2011). + Note: This routine differs from the Fortran routine (and the algorithm described in Reisswig et al. 2013) + in the following way: We do not treat the specific internal energy "epsilon" in any special way. + Rather, to avoid reconstruction to eps<0 for gamma-laws, we rely on + 'eps<0'-error checks after reconstruction (a similar treatment is also necessary for other reconstruction methods such as WENO5 and MP5). +*/ + + +/* + 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 double approx_at_cell_interface(const double* const restrict 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 double limit(const double* const restrict a, const double ah, const double C, const int i) +{ + if ((MIN(a[i],a[i+1]) <= ah) && (ah <= MAX(a[i],a[i+1]))) { + return ah; + } else { + const double D2a = 3.0 * ((a[i] + a[i+1]) - 2.0*ah ); + const double D2aL = ((a[i-1] + a[i+1]) - 2.0*a[i] ); + const double D2aR = ((a[i] + a[i+2]) - 2.0*a[i+1]); + const double D2aLim = copysign(1.0, D2a)*MIN3(C*fabs(D2aL), C*fabs(D2aR), fabs(D2a)) * 1.0/3.0; + if (D2a*D2aR >= 0 && D2a*D2aL >= 0) + return 0.5*(a[i]+a[i+1]) - D2aLim; + else + return 0.5*(a[i]+a[i+1]); + } + return 0; +} + + +/// Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34. +/// This does not use the check for deviations from a cubic. Thus, it gets away with only 3 stencil points! +static inline void monotonize(double* const restrict aminus, + const double* const restrict a, + double* const restrict aplus, + const double C, + const int i) +{ + double D2aLim = 0; + double rhi = 0; + const double daplus = aplus[i]-a[i]; + const double daminus = a[i]-aminus[i]; + if (daplus*daminus <= 0 || (a[i-2]-a[i])*(a[i]-a[i+2]) <= 0) { + const double D2a = - (12.0*a[i] - 6.0*(aminus[i]+aplus[i])); + const double D2aC = (a[i-1] + a[i+1]) - 2.0*a[i]; + const double D2aL = (a[i-2] + a[i] ) - 2.0*a[i-1]; + const double D2aR = (a[i] + a[i+2]) - 2.0*a[i+1]; + if (copysign(1.0, D2a) == copysign(1.0, D2aC) && copysign(1.0, D2a) == copysign(1.0, D2aL) && copysign(1.0, D2a) == copysign(1.0, D2aR)) + D2aLim = copysign(1.0, D2a) * MIN4(C*fabs(D2aL), C*fabs(D2aR), C*fabs(D2aC), fabs(D2a)); + if (!(fabs(D2a) <= 1e-12*MAX5(fabs(a[i-2]), fabs(a[i-1]), fabs(a[i]), fabs(a[i+1]), fabs(a[i+2])))) + rhi = D2aLim / D2a; + if (! (rhi >= 1.0 - 1e-12)) { + if (daplus*daminus < 0) { + aplus[i] = a[i] + daplus * rhi; + aminus[i] = a[i] - daminus * rhi; + } else if (fabs(daminus) >= 2.0*fabs(daplus)) { + aminus[i] = a[i] - (2.0*(1.0-rhi)*daplus + rhi*daminus); + } else if (fabs(daplus) >= 2.0*fabs(daminus)) { + aplus[i] = a[i] + (2.0*(1.0-rhi)*daminus + rhi*daplus); + } + } + //trivial_rp(i-1) = .false. + //trivial_rp(i) = .false. + } else { + if (fabs(daplus) >= 2.0*fabs(daminus)) + aplus[i] = a[i] + 2.0*daminus; + else if (fabs(daminus) >= 2.0*fabs(daplus)) + aminus[i] = a[i] - 2.0*daplus; + //trivial_rp(i-1) = .false. + //trivial_rp(i) = .false. + } + + return; +} + + +template +void GRHydro_eppm1d_cxx(const int nx, + const double dx, + const double* const restrict rho, + const double* const restrict velx, + const double* const restrict vely, + const double* const restrict velz, + const double* const restrict eps, + const double* const restrict press, + const double* const restrict temp, + const double* const restrict ye, + const double* const restrict Bvcx, + const double* const restrict Bvcy, + const double* const restrict Bvcz, + const double* const restrict psidc, + double* const restrict rhominus, + double* const restrict velxminus, + double* const restrict velyminus, + double* const restrict velzminus, + double* const restrict epsminus, + double* const restrict tempminus, + double* const restrict yeminus, + double* const restrict Bvcxminus, + double* const restrict Bvcyminus, + double* const restrict Bvczminus, + double* const restrict psidcminus, + double* const restrict rhoplus, + double* const restrict velxplus, + double* const restrict velyplus, + double* const restrict velzplus, + double* const restrict epsplus, + double* const restrict tempplus, + double* const restrict yeplus, + double* const restrict Bvcxplus, + double* const restrict Bvcyplus, + double* const restrict Bvczplus, + double* const restrict psidcplus) +{ + DECLARE_CCTK_PARAMETERS; + + + // We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, + // then checking for (13) of Colella & Sekora 2008 and applying + // (18) and (19) if (13) is not satisfied. This is done with LIMIT. + for (int i=1; i < nx-2; ++i) { + rhoplus[i] = approx_at_cell_interface(rho, i); + rhoplus[i] = limit(rho, rhoplus[i], enhanced_ppm_C2, i); + rhominus[i+1] = rhoplus[i]; + + epsplus[i] = approx_at_cell_interface(eps, i); + epsplus[i] = limit(eps, epsplus[i], enhanced_ppm_C2, i); + epsminus[i+1] = epsplus[i]; + + velxplus[i] = approx_at_cell_interface(velx, i); + velxplus[i] = limit(velx, velxplus[i], enhanced_ppm_C2, i); + velxminus[i+1] = velxplus[i]; + + velyplus[i] = approx_at_cell_interface(vely, i); + velyplus[i] = limit(vely, velyplus[i], enhanced_ppm_C2, i); + velyminus[i+1] = velyplus[i]; + + velzplus[i] = approx_at_cell_interface(velz, i); + velzplus[i] = limit(velz, velzplus[i], enhanced_ppm_C2, i); + velzminus[i+1] = velzplus[i]; + + if(do_ye) { + if(do_temp) { + tempplus[i] = approx_at_cell_interface(temp, i); + tempplus[i] = limit(temp, tempplus[i], enhanced_ppm_C2, i); + tempminus[i+1] = tempplus[i]; + } + yeplus[i] = approx_at_cell_interface(ye, i); + yeplus[i] = limit(ye, yeplus[i], enhanced_ppm_C2, i); + yeminus[i+1] = yeplus[i]; + } + if(do_mhd) { + Bvcxplus[i] = approx_at_cell_interface(Bvcx, i); + Bvcxplus[i] = limit(Bvcx, Bvcxplus[i], enhanced_ppm_C2, i); + Bvcxminus[i+1] = Bvcxplus[i]; + + Bvcyplus[i] = approx_at_cell_interface(Bvcy, i); + Bvcyplus[i] = limit(Bvcy, Bvcyplus[i], enhanced_ppm_C2, i); + Bvcyminus[i+1] = Bvcyplus[i]; + + Bvczplus[i] = approx_at_cell_interface(Bvcz, i); + Bvczplus[i] = limit(Bvcz, Bvczplus[i], enhanced_ppm_C2, i); + Bvczminus[i+1] = Bvczplus[i]; + + if(dc_flag) { + psidcplus[i] = approx_at_cell_interface(psidc, i); + psidcplus[i] = limit(psidc, psidcplus[i], enhanced_ppm_C2, i); + psidcminus[i+1] = psidcplus[i]; + } + } + + } + + // Monotonize + for(int i=GRHydro_stencil-1; i dpress(nx,0); + for(int i=1; i tilde_flatten(nx,0); + for(int i=2; 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 +#define INSTANTIATE(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect) \ +template void GRHydro_eppm1d_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); + +//template template + + +// 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 941afc1..ea823e4 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -84,6 +84,7 @@ SRCS = Utils.F90 \ Con2PrimM_fortran_interfaces.F90 \ GRHydro_Prim2Con.c \ GRHydro_PPM.cc \ + GRHydro_ePPM.cc \ GRHydro_PPMReconstruct_drv_opt.cc \ GRHydro_Wrappers.F90 -- cgit v1.2.3