diff --git a/ReconX/src/ppm.hxx b/ReconX/src/ppm.hxx index 63af84ce..088c69bf 100644 --- a/ReconX/src/ppm.hxx +++ b/ReconX/src/ppm.hxx @@ -18,98 +18,74 @@ using namespace std; using namespace Arith; using namespace Loop; -/* PPM reconstruction scheme. See Colella & Woodward (1984) (e.g. at - * https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf) */ -inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array -ppm(const GF3D2 &gf_var, - const array, 5> &cells, const CCTK_INT &dir, - const bool &gf_is_rho, const GF3D2 &gf_press, - const GF3D2 &gf_vel_dir, +/** + * @brief PPM reconstruction scheme. See Colella & Woodward (1984) (e.g. at + * https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf) + * + */ + +template +inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array +ppm(T gf_Imm, T gf_Im, T gf_I, T gf_Ip, T gf_Ipp, T press_Imm, T press_Im, + T press_Ip, T press_Ipp, T vel_dir_Im, T vel_dir_Ip, const bool &gf_is_rho, const reconstruct_params_t &reconstruct_params) { - // Unpack all cells in the stencil - const auto &Imm = cells[0]; - const auto &Im = cells[1]; - const auto &I = cells[2]; - const auto &Ip = cells[3]; - const auto &Ipp = cells[4]; // Unpack all PPM parameters const bool &ppm_shock_detection = reconstruct_params.ppm_shock_detection; const bool &ppm_zone_flattening = reconstruct_params.ppm_zone_flattening; - const CCTK_REAL &poly_k = reconstruct_params.poly_k; - const CCTK_REAL &poly_gamma = reconstruct_params.poly_gamma; - const CCTK_REAL &ppm_eta1 = reconstruct_params.ppm_eta1; - const CCTK_REAL &ppm_eta2 = reconstruct_params.ppm_eta2; - const CCTK_REAL &ppm_eps = reconstruct_params.ppm_eps; - const CCTK_REAL &ppm_eps_shock = reconstruct_params.ppm_eps_shock; - const CCTK_REAL &ppm_small = reconstruct_params.ppm_small; - const CCTK_REAL &ppm_omega1 = reconstruct_params.ppm_omega1; - const CCTK_REAL &ppm_omega2 = reconstruct_params.ppm_omega2; - - // Grid function at neighboring cells - const CCTK_REAL &gf_Imm = gf_var(Imm); - const CCTK_REAL &gf_Im = gf_var(Im); - const CCTK_REAL &gf_I = gf_var(I); - const CCTK_REAL &gf_Ip = gf_var(Ip); - const CCTK_REAL &gf_Ipp = gf_var(Ipp); + const T &poly_k = reconstruct_params.poly_k; + const T &poly_gamma = reconstruct_params.poly_gamma; + const T &ppm_eta1 = reconstruct_params.ppm_eta1; + const T &ppm_eta2 = reconstruct_params.ppm_eta2; + const T &ppm_eps = reconstruct_params.ppm_eps; + const T &ppm_eps_shock = reconstruct_params.ppm_eps_shock; + const T &ppm_small = reconstruct_params.ppm_small; + const T &ppm_omega1 = reconstruct_params.ppm_omega1; + const T &ppm_omega2 = reconstruct_params.ppm_omega2; // Helpers - const CCTK_REAL diff_Im = gf_I - gf_Imm; - const CCTK_REAL diff_I = gf_Ip - gf_Im; - const CCTK_REAL diff_Ip = gf_Ipp - gf_I; + const T diff_Im = gf_I - gf_Imm; + const T diff_I = gf_Ip - gf_Im; + const T diff_Ip = gf_Ipp - gf_I; - const CCTK_REAL delta_Im = 0.5 * diff_Im; - const CCTK_REAL delta_I = 0.5 * diff_I; - const CCTK_REAL delta_Ip = 0.5 * diff_Ip; + const T delta_Im = 0.5 * diff_Im; + const T delta_I = 0.5 * diff_I; + const T delta_Ip = 0.5 * diff_Ip; - const CCTK_REAL _2fabs_Im_Imm = 2 * fabs(gf_Im - gf_Imm); - const CCTK_REAL _2fabs_I_Im = 2 * fabs(gf_I - gf_Im); - const CCTK_REAL _2fabs_Ip_I = 2 * fabs(gf_Ip - gf_I); - const CCTK_REAL _2fabs_Ipp_Ip = 2 * fabs(gf_Ipp - gf_Ip); + const T _2fabs_Im_Imm = 2 * fabs(gf_Im - gf_Imm); + const T _2fabs_I_Im = 2 * fabs(gf_I - gf_Im); + const T _2fabs_Ip_I = 2 * fabs(gf_Ip - gf_I); + const T _2fabs_Ipp_Ip = 2 * fabs(gf_Ipp - gf_Ip); const bool same_sgn_Im = ((gf_I - gf_Im) * (gf_Im - gf_Imm) > 0); const bool same_sgn_I = ((gf_Ip - gf_I) * (gf_I - gf_Im) > 0); const bool same_sgn_Ip = ((gf_Ipp - gf_Ip) * (gf_Ip - gf_I) > 0); - const CCTK_REAL deltamod_Im = + const T deltamod_Im = same_sgn_Im ? sgn(delta_Im) * min(fabs(delta_Im), min(_2fabs_Im_Imm, _2fabs_I_Im)) : 0; - const CCTK_REAL deltamod_I = + const T deltamod_I = same_sgn_I ? sgn(delta_I) * min(fabs(delta_I), min(_2fabs_I_Im, _2fabs_Ip_I)) : 0; - const CCTK_REAL deltamod_Ip = + const T deltamod_Ip = same_sgn_Ip ? sgn(delta_Ip) * min(fabs(delta_Ip), min(_2fabs_Ip_I, _2fabs_Ipp_Ip)) : 0; /* Initial reconstructed states at the interfaces between cells Im/I ans I/Ip * NOTE: not const because they may change later */ - const CCTK_REAL gf_Im_I = + const T gf_Im_I = 0.5 * (gf_Im + gf_I) + (1. / 6.) * (deltamod_Im - deltamod_I); - const CCTK_REAL gf_I_Ip = + const T gf_I_Ip = 0.5 * (gf_I + gf_Ip) + (1. / 6.) * (deltamod_I - deltamod_Ip); - CCTK_REAL rc_low = gf_Im_I; - CCTK_REAL rc_up = gf_I_Ip; - - /* Pressure may or may not be needed, or may not be needed at all points - * (depending on whether shock detection and zone flattening are activated or - * not) */ - // const CCTK_REAL &press_Immm = gf_press(Immm); // Only used in the original - // zone flattening scheme - const CCTK_REAL &press_Imm = gf_press(Imm); - const CCTK_REAL &press_Im = gf_press(Im); - // const CCTK_REAL &press_I = gf_press(I); // Only used in the original - // zone flattening scheme - const CCTK_REAL &press_Ip = gf_press(Ip); - const CCTK_REAL &press_Ipp = gf_press(Ipp); - // const CCTK_REAL &press_Ippp = gf_press(Ippp); // Only used in the original - // zone flattening scheme - - const CCTK_REAL diff_press_I = press_Ip - press_Im; - const CCTK_REAL min_press_I = min(press_Im, press_Ip); + T rc_low = gf_Im_I; + T rc_up = gf_I_Ip; + + const T diff_press_I = press_Ip - press_Im; + const T min_press_I = min(press_Im, press_Ip); /* Shock detection (eqs. 1.15, 1.16, 1.17 with uniform grid spacing). * This is only applied to rho and only if the shock is marked as a contact @@ -117,20 +93,21 @@ ppm(const GF3D2 &gf_var, // FIXME: contact discontinuity check only valid for polytropic/ideal-fluid // EOS if (ppm_shock_detection and gf_is_rho) { - const CCTK_REAL k_gamma = poly_k * poly_gamma; + const T k_gamma = poly_k * poly_gamma; const bool contact_disc = (k_gamma * fabs(diff_I) * min_press_I >= fabs(diff_press_I) * min(gf_Ip, gf_Im)); if (contact_disc) { // This assumes gf_var is rho (gf_is_rho is true) - const CCTK_REAL d2rho_Im = gf_I - 2 * gf_Im + gf_Imm; - const CCTK_REAL d2rho_Ip = gf_Ipp - 2 * gf_Ip + gf_I; - const CCTK_REAL d2_prod = d2rho_Im * d2rho_Ip; - const bool cond2 = (fabs(diff_I) - ppm_eps_shock*min(fabs(gf_Ip), fabs(gf_Im))) > 0.; - - const CCTK_REAL eta_tilde_I = - ((d2_prod < 0) and cond2) ? ((-1. / 6.) * (d2rho_Ip - d2rho_Im) / diff_I) : 0; - const CCTK_REAL eta_I = - max(0., min(ppm_eta1 * (eta_tilde_I - ppm_eta2), 1.)); + const T d2rho_Im = gf_I - 2 * gf_Im + gf_Imm; + const T d2rho_Ip = gf_Ipp - 2 * gf_Ip + gf_I; + const T d2_prod = d2rho_Im * d2rho_Ip; + const bool cond2 = + (fabs(diff_I) - ppm_eps_shock * min(fabs(gf_Ip), fabs(gf_Im))) > 0.; + + const T eta_tilde_I = ((d2_prod < 0) and cond2) + ? ((-1. / 6.) * (d2rho_Ip - d2rho_Im) / diff_I) + : 0; + const T eta_I = max(0., min(ppm_eta1 * (eta_tilde_I - ppm_eta2), 1.)); rc_low = (1 - eta_I) * rc_low + eta_I * (gf_Im + 0.5 * deltamod_Im); rc_up = (1 - eta_I) * rc_up + eta_I * (gf_Ip - 0.5 * deltamod_Ip); @@ -144,52 +121,32 @@ ppm(const GF3D2 &gf_var, * be a major issue and is done in many GRMHD codes (e.g. WhiskyMHD, GRHydro, * Spritz, IllinoisGRMHD). */ if (ppm_zone_flattening) { - const CCTK_REAL w_I = ( (fabs(diff_press_I) > ppm_eps * min_press_I) and - (gf_vel_dir(Im) > gf_vel_dir(Ip))) - ? 1. - : 0; - const CCTK_REAL ftilde_I = (fabs(press_Ipp - press_Imm) < ppm_small) ? 1.0 : - max(0., 1. - w_I * max(0., ppm_omega2 * - ( (diff_press_I / (press_Ipp - press_Imm)) - ppm_omega1))); - - const CCTK_REAL one_minus_ftilde_I_gfI = (1 - ftilde_I)*gf_I; + const T w_I = ((fabs(diff_press_I) > ppm_eps * min_press_I) and + (vel_dir_Im > vel_dir_Ip)) + ? 1. + : 0; + const T ftilde_I = + (fabs(press_Ipp - press_Imm) < ppm_small) + ? 1.0 + : max(0., + 1. - w_I * max(0., ppm_omega2 * ((diff_press_I / + (press_Ipp - press_Imm)) - + ppm_omega1))); + + const T one_minus_ftilde_I_gfI = (1 - ftilde_I) * gf_I; rc_low = ftilde_I * rc_low + one_minus_ftilde_I_gfI; rc_up = ftilde_I * rc_up + one_minus_ftilde_I_gfI; - - // This would require one more ghost cell and it's not worth it - /*if (diff_press_I < 0) { - const CCTK_REAL diff_press_Ip = press_Ipp - press_I; - const CCTK_REAL w_Ip = (diff_press_Ip > ppm_eps*min(press_I, - press_Ipp) and gf_vel_dir(I) > gf_vel_dir(Ipp)) ? 1 : 0; const CCTK_REAL - ftilde_Ip = 1 - max(0., w_Ip*ppm_omega2*(diff_press_Ip/(press_Ippp - - press_Im) - ppm_omega1)); const CCTK_REAL f_I = max(ftilde_I, - ftilde_Ip); const CCTK_REAL one_minus_fI = 1 - f_I; const CCTK_REAL fI_gfI - = f_I*gf_I; rc_low = fI_gfI + one_minus_fI*rc_low; rc_up = fI_gfI + - one_minus_fI*rc_up; - } - - else { - const CCTK_REAL diff_press_Im = press_I - press_Imm; - const CCTK_REAL w_Im = (diff_press_Im > ppm_eps*min(press_Imm, - press_I) and gf_vel_dir(Imm) > gf_vel_dir(I)) ? 1 : 0; const CCTK_REAL - ftilde_Im = 1 - max(0., w_Im*ppm_omega2*(diff_press_Im/(press_Ip - - press_Immm) - ppm_omega1)); const CCTK_REAL f_I = max(ftilde_I, - ftilde_Im); const CCTK_REAL one_minus_fI = 1 - f_I; const CCTK_REAL fI_gfI - = f_I*gf_I; rc_low = fI_gfI + one_minus_fI*rc_low; rc_up = fI_gfI + - one_minus_fI*rc_up; - }*/ } // Monotonization (see eq. 1.10) if ((rc_up - gf_I) * (gf_I - rc_low) <= 0) { rc_low = gf_I; rc_up = gf_I; - } - else { - const CCTK_REAL diff_rc = rc_up - rc_low; - const CCTK_REAL diff_rc_sq = diff_rc * diff_rc; - const CCTK_REAL gf6_I = 6 * diff_rc * (gf_I - 0.5 * (rc_low + rc_up)); + } else { + const T diff_rc = rc_up - rc_low; + const T diff_rc_sq = diff_rc * diff_rc; + const T gf6_I = 6 * diff_rc * (gf_I - 0.5 * (rc_low + rc_up)); if (gf6_I > diff_rc_sq) { rc_low = 3 * gf_I - 2 * rc_up; @@ -201,9 +158,27 @@ ppm(const GF3D2 &gf_var, } // Return the lower and upper reconstructed states in cell I - const array rc = {rc_low, rc_up}; + const array rc = {rc_low, rc_up}; return rc; +} + +template +inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array +ppm_reconstruct(T gf_Immm, T gf_Imm, T gf_Im, T gf_Ip, T gf_Ipp, T gf_Ippp, + T press_Immm, T press_Imm, T press_Im, T press_Ip, T press_Ipp, + T press_Ippp, T vel_dir_Imm, T vel_dir_Im, T vel_dir_Ip, + T vel_dir_Ipp, const bool &gf_is_rho, + const reconstruct_params_t &reconstruct_params) { + + const auto rc_Im{ppm(gf_Immm, gf_Imm, gf_Im, gf_Ip, gf_Ipp, press_Immm, + press_Imm, press_Ip, press_Ipp, vel_dir_Imm, vel_dir_Ip, + gf_is_rho, reconstruct_params)}; + + const auto rc_Ip{ppm(gf_Imm, gf_Im, gf_Ip, gf_Ipp, gf_Ippp, press_Imm, + press_Im, press_Ipp, press_Ippp, vel_dir_Im, vel_dir_Ipp, + gf_is_rho, reconstruct_params)}; + return {rc_Im.at(1), rc_Ip.at(0)}; } } // namespace ReconX diff --git a/ReconX/src/reconstruct.hxx b/ReconX/src/reconstruct.hxx index 7b56fba0..56ee8502 100644 --- a/ReconX/src/reconstruct.hxx +++ b/ReconX/src/reconstruct.hxx @@ -63,17 +63,23 @@ reconstruct(const GF3D2 &gf_var, const PointDesc &p, } case reconstruction_t::ppm: { - const array, 5> cells_Im = {Immm, Imm, Im, Ip, Ipp}; - const array, 5> cells_Ip = {Imm, Im, Ip, Ipp, Ippp}; + return ppm_reconstruct( + gf_var(Immm), gf_var(Imm), gf_var(Im), gf_var(Ip), gf_var(Ipp), + gf_var(Ippp), gf_press(Immm), gf_press(Imm), gf_press(Im), gf_press(Ip), + gf_press(Ipp), gf_press(Ippp), gf_vel_dir(Imm), gf_vel_dir(Im), + gf_vel_dir(Ip), gf_vel_dir(Ipp), gf_is_rho, reconstruct_params); + } - const array rc_Im = - ppm(gf_var, cells_Im, dir, gf_is_rho, gf_press, gf_vel_dir, - reconstruct_params); - const array rc_Ip = - ppm(gf_var, cells_Ip, dir, gf_is_rho, gf_press, gf_vel_dir, - reconstruct_params); + case reconstruction_t::wenoz: { + return wenoz_reconstruct(gf_var(Immm), gf_var(Imm), gf_var(Im), gf_var(Ip), + gf_var(Ipp), gf_var(Ippp), + reconstruct_params.weno_eps); + } - return array{rc_Im.at(1), rc_Ip.at(0)}; + case reconstruction_t::mp5: { + return mp5_reconstruct(gf_var(Immm), gf_var(Imm), gf_var(Im), gf_var(Ip), + gf_var(Ipp), gf_var(Ippp), + reconstruct_params.mp5_alpha); } case reconstruction_t::eppm: { @@ -88,18 +94,6 @@ reconstruct(const GF3D2 &gf_var, const PointDesc &p, return array{rc_Im.at(1), rc_Ip.at(0)}; } - case reconstruction_t::wenoz: { - return wenoz_reconstruct(gf_var(Immm), gf_var(Imm), gf_var(Im), gf_var(Ip), - gf_var(Ipp), gf_var(Ippp), - reconstruct_params.weno_eps); - } - - case reconstruction_t::mp5: { - return mp5_reconstruct(gf_var(Immm), gf_var(Imm), gf_var(Im), gf_var(Ip), - gf_var(Ipp), gf_var(Ippp), - reconstruct_params.mp5_alpha); - } - default: assert(0); }