diff --git a/AsterX/src/unit_tests/minmod.cxx b/AsterX/src/unit_tests/minmod.cxx index 57531c33..9ee3aa57 100644 --- a/AsterX/src/unit_tests/minmod.cxx +++ b/AsterX/src/unit_tests/minmod.cxx @@ -19,70 +19,117 @@ void AsterXTests::test_minmod(std::mt19937_64 &engine, int repetitions) { for (int i = 0; i < repetitions; i++) { - CCTK_VINFO("Testing minmod repetition %i", i); - { - bool success{true}; + const auto p1{positive_distrib(engine)}; + const auto p2{positive_distrib(engine)}; + const auto p3{positive_distrib(engine)}; + const auto p4{positive_distrib(engine)}; - const auto p1{positive_distrib(engine)}; - const auto p2{positive_distrib(engine)}; + const auto n1{negative_distrib(engine)}; + const auto n2{negative_distrib(engine)}; + const auto n3{negative_distrib(engine)}; + const auto n4{negative_distrib(engine)}; - const auto n1{negative_distrib(engine)}; - const auto n2{negative_distrib(engine)}; + const auto r1{real_distrib(engine)}; + const auto r2{real_distrib(engine)}; - const auto r1{real_distrib(engine)}; - const auto r2{real_distrib(engine)}; + const auto positive_minmod{minmod(p1, p2)}; + const auto negative_minmod{minmod(n1, n2)}; - const auto positive_minmod{minmod(p1, p2)}; - const auto negative_minmod{minmod(n1, n2)}; + const auto mixed_minmod_1{minmod(n1, p1)}; + const auto mixed_minmod_2{minmod(n2, p2)}; - const auto mixed_minmod_1{minmod(n1, p1)}; - const auto mixed_minmod_2{minmod(n2, p2)}; + const auto real_minmod_1{minmod(r1, r2)}; - const auto real_minmod(minmod(r1, r2)); + const auto min_p1p2{min(p1, p2)}; + const auto max_n1n2{max(n1, n2)}; - const auto min_p1p2{min(p1, p2)}; - const auto max_n1n2{max(n1, n2)}; + const auto max_abs{max(abs(r1), abs(r2))}; - const auto max_abs{max(abs(r1), abs(r2))}; + CCTK_VINFO("Testing minmod of positive numbers, rep. %i", i); + { if (!isapprox(positive_minmod, min_p1p2)) { CCTK_VINFO(" FAILED. Reason: minmod with positive " - "numbers expected %.f16 and got %.f16", - min_p1p2, positive_minmod); - success &= false; + "numbers %.f16 and %.f16 expected to be %.f16 but resulted " + "in %.f16", + p1, p2, min_p1p2, positive_minmod); + } else { + CCTK_VINFO(" PASSED."); } + } + CCTK_VINFO("Testing minmod of negative numbers, rep. %i", i); + { if (!isapprox(negative_minmod, max_n1n2)) { CCTK_VINFO(" FAILED. Reason: minmod with negative " - "numbers expected %.f16 and got %.f16", - min_p1p2, negative_minmod); - success &= false; + "numbers %.f16 and %.f16 expected to be %.f16 but resulted " + "in %.f16", + n1, n2, max_n1n2, negative_minmod); + } else { + CCTK_VINFO(" PASSED."); } + } + CCTK_VINFO("Testing minmod of mixed sign numbers, rep. %i", i); + { if (!isapprox(mixed_minmod_1, CCTK_REAL(0))) { - CCTK_VINFO(" FAILED. Reason: minmod with mixed sign " - "numbers numbers expected 0.0 and got %.f16", - mixed_minmod_1); - success &= false; + CCTK_VINFO( + " FAILED. Reason: minmod with mixed sign " + "numbers %.f16 and %.f16 expected tp be 0.0 but resulted in %.f16", + n1, p1, mixed_minmod_1); + } else { + CCTK_VINFO(" PASSED."); } - if (!isapprox(mixed_minmod_2, CCTK_REAL(0))) { - CCTK_VINFO(" FAILED. Reason: minmod with mixed sign " - "numbers numbers expected 0.0 and got %.f16", - mixed_minmod_2); - success &= false; + if (!isapprox(mixed_minmod_1, CCTK_REAL(0))) { + CCTK_VINFO( + " FAILED. Reason: minmod with mixed sign " + "numbers %.f16 and %.f16 expected tp be 0.0 but resulted in %.f16", + n2, p2, mixed_minmod_2); + } else { + CCTK_VINFO(" PASSED."); } + } - if (!(real_minmod < max_abs || isapprox(real_minmod, max_abs))) { + CCTK_VINFO("Testing minmod(a, b) <= max(|a|,|b|), rep. %i", i); + { + if (!(real_minmod_1 < max_abs || isapprox(real_minmod_1, max_abs))) { CCTK_VINFO(" FAILED. Reason: minmod with reals %.f16 and %.f16 " "returned %.f, which is not bounded by %.f16", - r1, r2, real_minmod, max_abs); - success &= false; + r1, r2, real_minmod_1, max_abs); + } else { + CCTK_VINFO(" PASSED."); } + } + + CCTK_VINFO("Testing minmod(a, b, c, d) = minmod(minmod(a, b), minmod(c, " + "d)) for positive arguments, rep. %i", + i); + { + const auto expected_min{min(p1, min(p2, min(p3, p4)))}; + const auto computed_min(minmod(minmod(p1, p2), minmod(p3, p4))); + + if (!isapprox(expected_min, computed_min)) { + CCTK_VINFO(" FAILED. Reason: expected %.f16 and got %.f16", + expected_min, computed_min); + } else { + CCTK_VINFO(" PASSED."); + } + } + + CCTK_VINFO("Testing minmod(a, b, c, d) = minmod(minmod(a, b), minmod(c, " + "d)) for negative arguments, rep. %i", + i); + { + const auto expected_max{max(n1, max(n2, max(n3, n4)))}; + const auto computed_max(minmod(minmod(n1, n2), minmod(n3, n4))); - if (success) { + if (!isapprox(expected_max, computed_max)) { + CCTK_VINFO(" FAILED. Reason: expected %.f16 and got %.f16", + expected_max, computed_max); + } else { CCTK_VINFO(" PASSED."); } } } -} +} \ No newline at end of file diff --git a/AsterX/test/Balsara1_shocktube_godunov.par b/AsterX/test/Balsara1_shocktube_godunov.par new file mode 100644 index 00000000..f595315d --- /dev/null +++ b/AsterX/test/Balsara1_shocktube_godunov.par @@ -0,0 +1,103 @@ +ActiveThorns = " + CarpetX + IOUtil + ODESolvers + ADMBase + HydroBase + TmunuBase + AsterSeeds + AsterX + EOSX +" + +$nlevels = 1 +$ncells = 100 + +Cactus::cctk_show_schedule = yes + +Cactus::presync_mode = "mixed-error" + +Cactus::terminate = "time" +Cactus::cctk_final_time = 0.40 + +ADMBase::initial_data = "Cartesian Minkowski" +ADMBase::initial_lapse = "one" +ADMBase::initial_shift = "zero" +ADMBase::initial_dtlapse = "none" +ADMBase::initial_dtshift = "none" + +CarpetX::verbose = no + +CarpetX::xmin = -0.5 +CarpetX::ymin = -0.5 +CarpetX::zmin = -0.5 + +CarpetX::xmax = +0.5 +CarpetX::ymax = +0.5 +CarpetX::zmax = +0.5 + +CarpetX::boundary_x = "linear extrapolation" +CarpetX::boundary_y = "linear extrapolation" +CarpetX::boundary_z = "linear extrapolation" + +CarpetX::boundary_upper_x = "linear extrapolation" +CarpetX::boundary_upper_y = "linear extrapolation" +CarpetX::boundary_upper_z = "linear extrapolation" + +CarpetX::ncells_x = $ncells +CarpetX::ncells_y = 2 +CarpetX::ncells_z = 2 + +CarpetX::max_num_levels = $nlevels +CarpetX::regrid_every = 100000 +CarpetX::blocking_factor_x = 1 +CarpetX::blocking_factor_y = 1 +CarpetX::blocking_factor_z = 1 +CarpetX::regrid_error_threshold = 0.01 + +CarpetX::prolongation_type = "ddf" +CarpetX::ghost_size = 3 +CarpetX::dtfac = 0.25 + +AsterSeeds::test_type = "1DTest" +AsterSeeds::test_case = "Balsara1" + +ReconX::reconstruction_method = "godunov" +AsterX::flux_type = "HLLE" +AsterX::c2p_tol = 1e-8 +AsterX::max_iter = 100 +AsterX::debug_mode = "yes" + +Con2PrimFactory::unit_test = "yes" + +EOSX::evol_eos_name = "IdealGas" +EOSX::gl_gamma = 2.0 +EOSX::poly_gamma = 2.0 +EOSX::rho_max = 10000 +EOSX::eps_max = 10000 + +ODESolvers::method = "RK4" + +IO::out_dir = $parfile +IO::out_every = 10 +IO::parfile_write = no + +CarpetX::out_metadata = no +CarpetX::out_norm_vars = " " +CarpetX::out_norm_omit_unstable = yes + +CarpetX::out_silo_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" + +CarpetX::out_tsv_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" \ No newline at end of file diff --git a/AsterX/test/Balsara1_shocktube_minmod.par b/AsterX/test/Balsara1_shocktube_minmod.par new file mode 100644 index 00000000..40c56ae7 --- /dev/null +++ b/AsterX/test/Balsara1_shocktube_minmod.par @@ -0,0 +1,103 @@ +ActiveThorns = " + CarpetX + IOUtil + ODESolvers + ADMBase + HydroBase + TmunuBase + AsterSeeds + AsterX + EOSX +" + +$nlevels = 1 +$ncells = 50 + +Cactus::cctk_show_schedule = yes + +Cactus::presync_mode = "mixed-error" + +Cactus::terminate = "time" +Cactus::cctk_final_time = 0.40 + +ADMBase::initial_data = "Cartesian Minkowski" +ADMBase::initial_lapse = "one" +ADMBase::initial_shift = "zero" +ADMBase::initial_dtlapse = "none" +ADMBase::initial_dtshift = "none" + +CarpetX::verbose = no + +CarpetX::xmin = -0.5 +CarpetX::ymin = -0.5 +CarpetX::zmin = -0.5 + +CarpetX::xmax = +0.5 +CarpetX::ymax = +0.5 +CarpetX::zmax = +0.5 + +CarpetX::boundary_x = "linear extrapolation" +CarpetX::boundary_y = "linear extrapolation" +CarpetX::boundary_z = "linear extrapolation" + +CarpetX::boundary_upper_x = "linear extrapolation" +CarpetX::boundary_upper_y = "linear extrapolation" +CarpetX::boundary_upper_z = "linear extrapolation" + +CarpetX::ncells_x = $ncells +CarpetX::ncells_y = 2 +CarpetX::ncells_z = 2 + +CarpetX::max_num_levels = $nlevels +CarpetX::regrid_every = 100000 +CarpetX::blocking_factor_x = 1 +CarpetX::blocking_factor_y = 1 +CarpetX::blocking_factor_z = 1 +CarpetX::regrid_error_threshold = 0.01 + +CarpetX::prolongation_type = "ddf" +CarpetX::ghost_size = 3 +CarpetX::dtfac = 0.25 + +AsterSeeds::test_type = "1DTest" +AsterSeeds::test_case = "Balsara1" + +ReconX::reconstruction_method = "minmod" +AsterX::flux_type = "HLLE" +AsterX::c2p_tol = 1e-8 +AsterX::max_iter = 100 +AsterX::debug_mode = "yes" + +Con2PrimFactory::unit_test = "yes" + +EOSX::evol_eos_name = "IdealGas" +EOSX::gl_gamma = 2.0 +EOSX::poly_gamma = 2.0 +EOSX::rho_max = 10000 +EOSX::eps_max = 10000 + +ODESolvers::method = "RK4" + +IO::out_dir = $parfile +IO::out_every = 10 +IO::parfile_write = no + +CarpetX::out_metadata = no +CarpetX::out_norm_vars = " " +CarpetX::out_norm_omit_unstable = yes + +CarpetX::out_silo_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" + +CarpetX::out_tsv_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" \ No newline at end of file diff --git a/AsterX/test/Balsara1_shocktube_monocentral.par b/AsterX/test/Balsara1_shocktube_monocentral.par new file mode 100644 index 00000000..259f7a24 --- /dev/null +++ b/AsterX/test/Balsara1_shocktube_monocentral.par @@ -0,0 +1,103 @@ +ActiveThorns = " + CarpetX + IOUtil + ODESolvers + ADMBase + HydroBase + TmunuBase + AsterSeeds + AsterX + EOSX +" + +$nlevels = 1 +$ncells = 50 + +Cactus::cctk_show_schedule = yes + +Cactus::presync_mode = "mixed-error" + +Cactus::terminate = "time" +Cactus::cctk_final_time = 0.40 + +ADMBase::initial_data = "Cartesian Minkowski" +ADMBase::initial_lapse = "one" +ADMBase::initial_shift = "zero" +ADMBase::initial_dtlapse = "none" +ADMBase::initial_dtshift = "none" + +CarpetX::verbose = no + +CarpetX::xmin = -0.5 +CarpetX::ymin = -0.5 +CarpetX::zmin = -0.5 + +CarpetX::xmax = +0.5 +CarpetX::ymax = +0.5 +CarpetX::zmax = +0.5 + +CarpetX::boundary_x = "linear extrapolation" +CarpetX::boundary_y = "linear extrapolation" +CarpetX::boundary_z = "linear extrapolation" + +CarpetX::boundary_upper_x = "linear extrapolation" +CarpetX::boundary_upper_y = "linear extrapolation" +CarpetX::boundary_upper_z = "linear extrapolation" + +CarpetX::ncells_x = $ncells +CarpetX::ncells_y = 2 +CarpetX::ncells_z = 2 + +CarpetX::max_num_levels = $nlevels +CarpetX::regrid_every = 100000 +CarpetX::blocking_factor_x = 1 +CarpetX::blocking_factor_y = 1 +CarpetX::blocking_factor_z = 1 +CarpetX::regrid_error_threshold = 0.01 + +CarpetX::prolongation_type = "ddf" +CarpetX::ghost_size = 3 +CarpetX::dtfac = 0.25 + +AsterSeeds::test_type = "1DTest" +AsterSeeds::test_case = "Balsara1" + +ReconX::reconstruction_method = "monocentral" +AsterX::flux_type = "HLLE" +AsterX::c2p_tol = 1e-8 +AsterX::max_iter = 100 +AsterX::debug_mode = "yes" + +Con2PrimFactory::unit_test = "yes" + +EOSX::evol_eos_name = "IdealGas" +EOSX::gl_gamma = 2.0 +EOSX::poly_gamma = 2.0 +EOSX::rho_max = 10000 +EOSX::eps_max = 10000 + +ODESolvers::method = "RK4" + +IO::out_dir = $parfile +IO::out_every = 10 +IO::parfile_write = no + +CarpetX::out_metadata = no +CarpetX::out_norm_vars = " " +CarpetX::out_norm_omit_unstable = yes + +CarpetX::out_silo_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" + +CarpetX::out_tsv_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" \ No newline at end of file diff --git a/AsterX/test/Balsara1_shocktube_mp5.par b/AsterX/test/Balsara1_shocktube_mp5.par new file mode 100644 index 00000000..7411eae4 --- /dev/null +++ b/AsterX/test/Balsara1_shocktube_mp5.par @@ -0,0 +1,103 @@ +ActiveThorns = " + CarpetX + IOUtil + ODESolvers + ADMBase + HydroBase + TmunuBase + AsterSeeds + AsterX + EOSX +" + +$nlevels = 1 +$ncells = 50 + +Cactus::cctk_show_schedule = yes + +Cactus::presync_mode = "mixed-error" + +Cactus::terminate = "time" +Cactus::cctk_final_time = 0.40 + +ADMBase::initial_data = "Cartesian Minkowski" +ADMBase::initial_lapse = "one" +ADMBase::initial_shift = "zero" +ADMBase::initial_dtlapse = "none" +ADMBase::initial_dtshift = "none" + +CarpetX::verbose = no + +CarpetX::xmin = -0.5 +CarpetX::ymin = -0.5 +CarpetX::zmin = -0.5 + +CarpetX::xmax = +0.5 +CarpetX::ymax = +0.5 +CarpetX::zmax = +0.5 + +CarpetX::boundary_x = "linear extrapolation" +CarpetX::boundary_y = "linear extrapolation" +CarpetX::boundary_z = "linear extrapolation" + +CarpetX::boundary_upper_x = "linear extrapolation" +CarpetX::boundary_upper_y = "linear extrapolation" +CarpetX::boundary_upper_z = "linear extrapolation" + +CarpetX::ncells_x = $ncells +CarpetX::ncells_y = 2 +CarpetX::ncells_z = 2 + +CarpetX::max_num_levels = $nlevels +CarpetX::regrid_every = 100000 +CarpetX::blocking_factor_x = 1 +CarpetX::blocking_factor_y = 1 +CarpetX::blocking_factor_z = 1 +CarpetX::regrid_error_threshold = 0.01 + +CarpetX::prolongation_type = "ddf" +CarpetX::ghost_size = 3 +CarpetX::dtfac = 0.25 + +AsterSeeds::test_type = "1DTest" +AsterSeeds::test_case = "Balsara1" + +ReconX::reconstruction_method = "mp5" +AsterX::flux_type = "HLLE" +AsterX::c2p_tol = 1e-8 +AsterX::max_iter = 100 +AsterX::debug_mode = "yes" + +Con2PrimFactory::unit_test = "yes" + +EOSX::evol_eos_name = "IdealGas" +EOSX::gl_gamma = 2.0 +EOSX::poly_gamma = 2.0 +EOSX::rho_max = 10000 +EOSX::eps_max = 10000 + +ODESolvers::method = "RK4" + +IO::out_dir = $parfile +IO::out_every = 10 +IO::parfile_write = no + +CarpetX::out_metadata = no +CarpetX::out_norm_vars = " " +CarpetX::out_norm_omit_unstable = yes + +CarpetX::out_silo_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" + +CarpetX::out_tsv_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" \ No newline at end of file diff --git a/AsterX/test/Balsara1_shocktube_wenoz.par b/AsterX/test/Balsara1_shocktube_wenoz.par new file mode 100644 index 00000000..d4e10932 --- /dev/null +++ b/AsterX/test/Balsara1_shocktube_wenoz.par @@ -0,0 +1,103 @@ +ActiveThorns = " + CarpetX + IOUtil + ODESolvers + ADMBase + HydroBase + TmunuBase + AsterSeeds + AsterX + EOSX +" + +$nlevels = 1 +$ncells = 50 + +Cactus::cctk_show_schedule = yes + +Cactus::presync_mode = "mixed-error" + +Cactus::terminate = "time" +Cactus::cctk_final_time = 0.40 + +ADMBase::initial_data = "Cartesian Minkowski" +ADMBase::initial_lapse = "one" +ADMBase::initial_shift = "zero" +ADMBase::initial_dtlapse = "none" +ADMBase::initial_dtshift = "none" + +CarpetX::verbose = no + +CarpetX::xmin = -0.5 +CarpetX::ymin = -0.5 +CarpetX::zmin = -0.5 + +CarpetX::xmax = +0.5 +CarpetX::ymax = +0.5 +CarpetX::zmax = +0.5 + +CarpetX::boundary_x = "linear extrapolation" +CarpetX::boundary_y = "linear extrapolation" +CarpetX::boundary_z = "linear extrapolation" + +CarpetX::boundary_upper_x = "linear extrapolation" +CarpetX::boundary_upper_y = "linear extrapolation" +CarpetX::boundary_upper_z = "linear extrapolation" + +CarpetX::ncells_x = $ncells +CarpetX::ncells_y = 2 +CarpetX::ncells_z = 2 + +CarpetX::max_num_levels = $nlevels +CarpetX::regrid_every = 100000 +CarpetX::blocking_factor_x = 1 +CarpetX::blocking_factor_y = 1 +CarpetX::blocking_factor_z = 1 +CarpetX::regrid_error_threshold = 0.01 + +CarpetX::prolongation_type = "ddf" +CarpetX::ghost_size = 3 +CarpetX::dtfac = 0.25 + +AsterSeeds::test_type = "1DTest" +AsterSeeds::test_case = "Balsara1" + +ReconX::reconstruction_method = "wenoz" +AsterX::flux_type = "HLLE" +AsterX::c2p_tol = 1e-8 +AsterX::max_iter = 100 +AsterX::debug_mode = "yes" + +Con2PrimFactory::unit_test = "yes" + +EOSX::evol_eos_name = "IdealGas" +EOSX::gl_gamma = 2.0 +EOSX::poly_gamma = 2.0 +EOSX::rho_max = 10000 +EOSX::eps_max = 10000 + +ODESolvers::method = "RK4" + +IO::out_dir = $parfile +IO::out_every = 10 +IO::parfile_write = no + +CarpetX::out_metadata = no +CarpetX::out_norm_vars = " " +CarpetX::out_norm_omit_unstable = yes + +CarpetX::out_silo_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" + +CarpetX::out_tsv_vars = " + HydroBase::rho + HydroBase::vel + HydroBase::eps + HydroBase::press + HydroBase::Bvec +" \ No newline at end of file diff --git a/ReconX/src/minmod.hxx b/ReconX/src/minmod.hxx index 8847f355..ce1e9f79 100644 --- a/ReconX/src/minmod.hxx +++ b/ReconX/src/minmod.hxx @@ -1,25 +1,21 @@ #ifndef RECONX_MINMOD_HXX #define RECONX_MINMOD_HXX -#include - #include -#include -#include - -#include -#include -#include "reconx_utils.hxx" +#include +#include namespace ReconX { -using namespace std; -using namespace Arith; +using std::array; -template +template inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST T minmod(const T &x, const T &y) { + using std::fabs; + using std::signbit; + if (signbit(x) != signbit(y)) return T(0); if (fabs(x) < fabs(y)) @@ -28,6 +24,25 @@ inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST T minmod(const T &x, return y; } +template +inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array +minmod_reconstruct(T q_Imm, T q_Im, T q_Ip, T q_Ipp) { + // reconstructs values of Im and Ip at the common face between these + // two cells + + const T var_slope_p{q_Ipp - q_Ip}; + const T var_slope_c{q_Ip - q_Im}; + const T var_slope_m{q_Im - q_Imm}; + + // reconstructed Im on its "plus/right" side + const T var_m{q_Im + minmod(var_slope_c, var_slope_m) / 2}; + + // reconstructed Ip on its "minus/left" side + const T var_p{q_Ip - minmod(var_slope_p, var_slope_c) / 2}; + + return {var_m, var_p}; +} + } // namespace ReconX #endif // RECONX_MINMOD_HXX diff --git a/ReconX/src/monocentral.hxx b/ReconX/src/monocentral.hxx index 5b360489..65a67ffa 100644 --- a/ReconX/src/monocentral.hxx +++ b/ReconX/src/monocentral.hxx @@ -1,31 +1,49 @@ #ifndef RECONX_MONOCENTRAL_HXX #define RECONX_MONOCENTRAL_HXX -#include - #include -#include -#include - -#include -#include #include "reconx_utils.hxx" +#include +#include + namespace ReconX { -using namespace std; -using namespace Arith; +using std::array; -template +template inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST T monocentral(const T &x, const T &y) { + + using std::fabs; + using std::min; + if (sgn(x) != sgn(y)) return 0; else return sgn(x) * min(2 * fabs(x), min(2 * fabs(y), fabs(x + y) / 2)); } +template +inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array +monocentral_reconstruct(T q_Imm, T q_Im, T q_Ip, T q_Ipp) { + // reconstructs values of Im and Ip at the common face between these + // two cells + + // reconstructed Im on its "plus/right" side + T var_slope_p{q_Ip - q_Im}; + T var_slope_m{q_Im - q_Imm}; + const T var_m{q_Im + monocentral(var_slope_p, var_slope_m) / 2}; + + // reconstructed Ip on its "minus/left" side + var_slope_p = q_Ipp - q_Ip; + var_slope_m = q_Ip - q_Im; + const T var_p{q_Ip - monocentral(var_slope_p, var_slope_m) / 2}; + + return {var_m, var_p}; +} + } // namespace ReconX #endif // RECONX_MINMOD_HXX diff --git a/ReconX/src/mp5.hxx b/ReconX/src/mp5.hxx index 06d2ad7f..1ad5cb97 100644 --- a/ReconX/src/mp5.hxx +++ b/ReconX/src/mp5.hxx @@ -1,95 +1,96 @@ #ifndef RECONX_MP5_HXX #define RECONX_MP5_HXX -#include - #include -#include -#include - -#include -#include -#include "reconx_utils.hxx" #include "minmod.hxx" -namespace ReconX { - -using namespace std; -using namespace Arith; -using namespace Loop; +#include +namespace ReconX { // Compute the median of three numbers -template inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST -T median(T &x, T &y, T &z) { - return x + minmod(y-x, z-x); +template +inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST T median(T &x, T &y, + T &z) { + return x + minmod(y - x, z - x); } -/* Fifth-order monotonicity preserving (MP5) scheme of Suresh and Huynh 1997. -Paper: "Accurate monotonicity-preserving schemes with Runge–Kutta time stepping, J. Comput. Phys. 136 (1997) 83–99." -*/ - +/** + * Fifth-order monotonicity preserving (MP5) scheme of Suresh and Huynh + * 1997. Paper: "Accurate monotonicity-preserving schemes with Runge–Kutta time + * stepping, J. Comput. Phys. 136 (1997) 83–99." + */ +template inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST CCTK_REAL -mp5(const GF3D2 &gf_var, - const array, 5> &cells, - const reconstruct_params_t &reconstruct_params) { - - // Unpack all MP5 parameters - const CCTK_REAL mp5_alpha = reconstruct_params.mp5_alpha; - - // 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]; - - // 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 CCTK_REAL ul = (2*gf_Imm - 13*gf_Im + 47*gf_I + - 27*gf_Ip - 3*gf_Imm) / 60.0; - - const CCTK_REAL deltam = gf_I-gf_Im; - const CCTK_REAL deltap = gf_Ip-gf_I; - - const CCTK_REAL ump = gf_I + minmod(deltap, mp5_alpha*deltam ); - - if((ul - gf_I)*(ul - ump) < 0) { - return ul; - } - else { - const CCTK_REAL dm = gf_Imm + gf_I - 2*gf_Im; - const CCTK_REAL d = gf_Im + gf_Ip - 2*gf_I; - const CCTK_REAL dp = gf_I + gf_Ipp - 2*gf_Ip; +mp5(T gf_Imm, T gf_Im, T gf_I, T gf_Ip, T gf_Ipp, T mp5_alpha) { + + using std::max; + using std::min; - const CCTK_REAL dmp = minmod(minmod(4*d-dp, 4*dp-d), - minmod(d, dp)); - const CCTK_REAL dmm = minmod(minmod(4*dm - d, 4*d - dm), - minmod(dm, d)); + // Eq. (2.1) + const T ul{(2 * gf_Imm - 13 * gf_Im + 47 * gf_I + 27 * gf_Ip - 3 * gf_Ipp) / + 60.0}; - const CCTK_REAL ulc = gf_I + 0.5*deltam + 4.0/3.0*dmm; - const CCTK_REAL umd = 0.5*(gf_I + gf_Ip) - 0.5*dmp; + const T deltam{gf_I - gf_Im}; + const T deltap{gf_Ip - gf_I}; - const CCTK_REAL uul = gf_I + mp5_alpha*deltam; + // Eq. (2.12) + const T ump{gf_I + minmod(deltap, mp5_alpha * deltam)}; - const CCTK_REAL umin = std::max( - std::min(gf_I, std::min(gf_Ip, umd)), - std::min(gf_I, std::min(uul, ulc))); - const CCTK_REAL umax = std::min( - std::max(gf_I, std::max(gf_Ip, umd)), - std::max(gf_I, std::max(uul, ulc))); + // Condition Eq. (2.30) + if ((ul - gf_I) * (ul - ump) <= 1.0e-10) { + return ul; + } else { - return median(umin, ul, umax); + // Eq. (2.19) + const T dm{gf_Imm + gf_I - 2 * gf_Im}; + const T d{gf_Im + gf_Ip - 2 * gf_I}; + const T dp{gf_I + gf_Ipp - 2 * gf_Ip}; + // Eq. (2.27) + const T dmp{minmod(minmod(4 * d - dp, 4 * dp - d), minmod(d, dp))}; + const T dmm{minmod(minmod(4 * dm - d, 4 * d - dm), minmod(dm, d))}; + + // Eq. (2.8) + const T uul{gf_I + mp5_alpha * deltam}; + + // Eq. (2.16) + const T uav{(gf_I + gf_Ip) / 2.0}; + + // Eq. (2.28) + const T umd{uav - dmp / 2.0}; + + // Eq. (2.29) + const T ulc{gf_I + deltam / 2.0 + 4.0 / 3.0 * dmm}; + + // Eq. (2.24 a) + const T umin{max(min(gf_I, min(gf_Ip, umd)), min(gf_I, min(uul, ulc)))}; + + // Eq. (2.24 b) + const T umax{min(max(gf_I, max(gf_Ip, umd)), max(gf_I, max(uul, ulc)))}; + + // Eq. (2.26) + return median(ul, umin, umax); } - return 0; + return T{0}; +} +template +inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array +mp5_reconstruct(T gf_Immm, T gf_Imm, T gf_Im, T gf_Ip, T gf_Ipp, T gf_Ippp, + T mp5_alpha) { + // for the left cell, the plus side has sequence: Immm, Imm, Im, Ip, Ipp + // for the left cell, the minus side has sequence: Ipp, Ip, Im, Im, Immm + // here, we need the plus side + const auto rc_Im{mp5(gf_Immm, gf_Imm, gf_Im, gf_Ip, gf_Ipp, mp5_alpha)}; + + // for the right cell, the plus side has sequence: Imm, Im, Ip, Ipp, Ippp + // for the right cell, the minus side has sequence: Ippp, Ipp, Ip, Im, Imm + // here, we need the minus side + const auto rc_Ip{mp5(gf_Ippp, gf_Ipp, gf_Ip, gf_Im, gf_Imm, mp5_alpha)}; + + return {rc_Im, rc_Ip}; } } // namespace ReconX diff --git a/ReconX/src/ppm.hxx b/ReconX/src/ppm.hxx index 63af84ce..4879e7f2 100644 --- a/ReconX/src/ppm.hxx +++ b/ReconX/src/ppm.hxx @@ -1,115 +1,89 @@ #ifndef RECONX_PPM_HXX #define RECONX_PPM_HXX -#include - #include -#include -#include - -#include -#include #include "reconx_utils.hxx" +#include +#include + namespace ReconX { -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, +using std::array; + +/** + * @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]; + + using std::fabs; + using std::max; + using std::min; // 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 +91,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 +119,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 +156,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 937140c2..759c0bc5 100644 --- a/ReconX/src/reconstruct.hxx +++ b/ReconX/src/reconstruct.hxx @@ -4,8 +4,8 @@ #include #include -#include -#include +// #include +// #include #include "monocentral.hxx" #include "minmod.hxx" @@ -14,9 +14,11 @@ #include "wenoz.hxx" #include "mp5.hxx" +#include + namespace ReconX { -using namespace std; +using std::array; using namespace Arith; // enum class for different reconstruction routines @@ -49,50 +51,36 @@ reconstruct(const GF3D2 &gf_var, const PointDesc &p, switch (reconstruction) { case reconstruction_t::Godunov: { - CCTK_REAL var_m = gf_var(Im); - CCTK_REAL var_p = gf_var(Ip); - return array{var_m, var_p}; + return {gf_var(Im), gf_var(Ip)}; } case reconstruction_t::minmod: { - // reconstructs values of Im and Ip at the common face between these - // two cells - CCTK_REAL var_slope_p = gf_var(Ipp) - gf_var(Ip); - CCTK_REAL var_slope_c = gf_var(Ip) - gf_var(Im); - CCTK_REAL var_slope_m = gf_var(Im) - gf_var(Imm); - // reconstructed Im on its "plus/right" side - CCTK_REAL var_m = gf_var(Im) + minmod(var_slope_c, var_slope_m) / 2; - // reconstructed Ip on its "minus/left" side - CCTK_REAL var_p = gf_var(Ip) - minmod(var_slope_p, var_slope_c) / 2; - return array{var_m, var_p}; + return minmod_reconstruct(gf_var(Imm), gf_var(Im), gf_var(Ip), gf_var(Ipp)); } case reconstruction_t::monocentral: { - // reconstructs values of Im and Ip at the common face between these - // two cells - // reconstructed Im on its "plus/right" side - CCTK_REAL var_slope_p = gf_var(Ip) - gf_var(Im); - CCTK_REAL var_slope_m = gf_var(Im) - gf_var(Imm); - CCTK_REAL var_m = gf_var(Im) + monocentral(var_slope_p, var_slope_m) / 2; - // reconstructed Ip on its "minus/left" side - var_slope_p = gf_var(Ipp) - gf_var(Ip); - var_slope_m = gf_var(Ip) - gf_var(Im); - CCTK_REAL var_p = gf_var(Ip) - monocentral(var_slope_p, var_slope_m) / 2; - return array{var_m, var_p}; + return monocentral_reconstruct(gf_var(Imm), gf_var(Im), gf_var(Ip), + gf_var(Ipp)); } 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[1], rc_Ip[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: { @@ -106,37 +94,7 @@ reconstruct(const GF3D2 &gf_var, const PointDesc &p, eppm(gf_var, cells_Ip, gf_is_press, gf_press, gf_vel_dir, reconstruct_params); - return array{rc_Im[1], rc_Ip[0]}; - } - - case reconstruction_t::wenoz: { - const array, 5> cells_Im = {Immm, Imm, Im, Ip, Ipp}; - const array, 5> cells_Ip = {Imm, Im, Ip, Ipp, Ippp}; - - const array rc_Im = - wenoz(gf_var, cells_Im, reconstruct_params); - const array rc_Ip = - wenoz(gf_var, cells_Ip, reconstruct_params); - - return array{rc_Im[1], rc_Ip[0]}; - } - - case reconstruction_t::mp5: { - // for the left cell, the plus side has sequence: Immm, Imm, Im, Ip, Ipp - // for the left cell, the minus side has sequence: Ipp, Ip, Im, Im, Immm - // here, we need the plus side - const array, 5> cells_Im = {Immm, Imm, Im, Ip, Ipp}; - - // for the right cell, the plus side has sequence: Imm, Im, Ip, Ipp, Ippp - // for the right cell, the minus side has sequence: Ippp, Ipp, Ip, Im, Imm - // here, we need the minus side - - const array, 5> cells_Ip = {Ippp, Ipp, Ip, Im, Imm}; - - const CCTK_REAL rc_Im = mp5(gf_var, cells_Im, reconstruct_params); - const CCTK_REAL rc_Ip = mp5(gf_var, cells_Ip, reconstruct_params); - - return array{rc_Im, rc_Ip}; + return array{rc_Im.at(1), rc_Ip.at(0)}; } default: diff --git a/ReconX/src/wenoz.hxx b/ReconX/src/wenoz.hxx index a015170d..8a3eeed9 100644 --- a/ReconX/src/wenoz.hxx +++ b/ReconX/src/wenoz.hxx @@ -1,52 +1,32 @@ #ifndef RECONX_WENOZ_HXX #define RECONX_WENOZ_HXX -#include - #include -#include -#include -#include -#include +#include -#include "reconx_utils.hxx" +#include +#include namespace ReconX { -using namespace std; -using namespace Arith; -using namespace Loop; - -/* WENO-Z: Performs the reconstruction of a given variable using 5th order - * WENO-Z method based on Borges et al., ''An improved weighted essentially - * non-oscillatory scheme for hyperbolic conservation laws'', 2008) */ -// Also, see the Spritz code for details - -inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array -wenoz(const GF3D2 &gf_var, - const array, 5> &cells, - const reconstruct_params_t &reconstruct_params) { - - // Unpack all WENO-Z parameters - const CCTK_REAL weno_eps = reconstruct_params.weno_eps; - - // 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]; - - // 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); +using std::array; + +/** + * WENO-Z: Performs the reconstruction of a given variable using 5th order + * WENO-Z method based on Borges et al., "An improved weighted essentially + * non-oscillatory scheme for hyperbolic conservation laws", 2008). Also, see + * the Spritz code for details + */ +template +inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array +wenoz(T gf_Imm, T gf_Im, T gf_I, T gf_Ip, T gf_Ipp, T weno_eps) { + + using Arith::vec; + using std::abs; // Computing the smoothness indicators (Borges et al. 2008) - const vec betaZ{ + const vec betaZ{ (13.0 / 12.0) * pow2(gf_Imm - 2.0 * gf_Im + gf_I) + (1.0 / 4.0) * pow2(gf_Imm - 4.0 * gf_Im + 3.0 * gf_I), @@ -57,20 +37,20 @@ wenoz(const GF3D2 &gf_var, (1.0 / 4.0) * pow2(3.0 * gf_I - 4.0 * gf_Ip + gf_Ipp)}; // Defining tau5 based on eq. 25 of (Borges et al. 2008) - const CCTK_REAL tau5 = abs(betaZ(0) - betaZ(2)); + const T tau5{abs(betaZ(0) - betaZ(2))}; // Unnormalized weights // Optimal weights are chosen according to Del Zanna et al. 2007 - vec aux_alphaZ; + vec aux_alphaZ; aux_alphaZ(0) = 1.0 + tau5 / (betaZ(0) + weno_eps); aux_alphaZ(1) = 1.0 + tau5 / (betaZ(1) + weno_eps); aux_alphaZ(2) = 1.0 + tau5 / (betaZ(2) + weno_eps); - // const vec wt = {5.0 / 16.0, 10.0 / 16.0, 1.0 / 16.0}; + // const vec wt = {5.0 / 16.0, 10.0 / 16.0, 1.0 / 16.0}; // Original weights as suggested in (Borges et al. 2008) - const vec wt = {3.0 / 10.0, 3.0 / 5.0, 1.0 / 10.0}; + const vec wt{3.0 / 10.0, 3.0 / 5.0, 1.0 / 10.0}; - vec, 3> alphaZ; + vec, 3> alphaZ; // for minus side alphaZ(0)(0) = wt(0) * aux_alphaZ(0); @@ -83,12 +63,12 @@ wenoz(const GF3D2 &gf_var, alphaZ(2)(1) = wt(0) * aux_alphaZ(2); // Normalized weights for the reconstruction (Del Zanna et al. 2007) - const vec omega_denom([&](int j) ARITH_INLINE { + const vec omega_denom([&](int j) ARITH_INLINE { return alphaZ(0)(j) + alphaZ(1)(j) + alphaZ(2)(j); }); - const vec, 3> omegaZ([&](int j) ARITH_INLINE { - return vec( + const vec, 3> omegaZ([&](int j) ARITH_INLINE { + return vec( [&](int f) ARITH_INLINE { return alphaZ(j)(f) / omega_denom(f); }); }); @@ -96,30 +76,40 @@ wenoz(const GF3D2 &gf_var, // interfaces /* Spritz Weights: - const CCTK_REAL var_m = + const T var_m = (omegaZ(2)(0) / 8.0) * (3.0 * gf_Ipp - 10.0 * gf_Ip + 15.0 * gf_I) + (omegaZ(1)(0) / 8.0) * (-1.0 * gf_Ip + 6.0 * gf_I + 3.0 * gf_Im) + (omegaZ(0)(0) / 8.0) * (3.0 * gf_I + 6.0 * gf_Im - 1.0 * gf_Imm); - const CCTK_REAL var_p = + const T var_p = (omegaZ(0)(1) / 8.0) * (3.0 * gf_Imm - 10.0 * gf_Im + 15.0 * gf_I) + (omegaZ(1)(1) / 8.0) * (-1.0 * gf_Im + 6.0 * gf_I + 3.0 * gf_Ip) + (omegaZ(2)(1) / 8.0) * (3.0 * gf_I + 6.0 * gf_Ip - 1.0 * gf_Ipp); */ // GRHydro Weights: - const CCTK_REAL var_m = + const T var_m{ (omegaZ(2)(0) / 6.0) * (2.0 * gf_Ipp - 7.0 * gf_Ip + 11.0 * gf_I) + (omegaZ(1)(0) / 6.0) * (-1.0 * gf_Ip + 5.0 * gf_I + 2.0 * gf_Im) + - (omegaZ(0)(0) / 6.0) * (2.0 * gf_I + 5.0 * gf_Im - 1.0 * gf_Imm); + (omegaZ(0)(0) / 6.0) * (2.0 * gf_I + 5.0 * gf_Im - 1.0 * gf_Imm)}; - const CCTK_REAL var_p = + const T var_p{ (omegaZ(0)(1) / 6.0) * (2.0 * gf_Imm - 7.0 * gf_Im + 11.0 * gf_I) + (omegaZ(1)(1) / 6.0) * (-1.0 * gf_Im + 5.0 * gf_I + 2.0 * gf_Ip) + - (omegaZ(2)(1) / 6.0) * (2.0 * gf_I + 5.0 * gf_Ip - 1.0 * gf_Ipp); + (omegaZ(2)(1) / 6.0) * (2.0 * gf_I + 5.0 * gf_Ip - 1.0 * gf_Ipp)}; + + return {var_m, var_p}; +} + +template +inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST array +wenoz_reconstruct(T gf_Immm, T gf_Imm, T gf_Im, T gf_Ip, T gf_Ipp, T gf_Ippp, + T weno_eps) { + + const auto rc_Im{wenoz(gf_Immm, gf_Imm, gf_Im, gf_Ip, gf_Ipp, weno_eps)}; + const auto rc_Ip{wenoz(gf_Imm, gf_Im, gf_Ip, gf_Ipp, gf_Ippp, weno_eps)}; - const array var_rc = {var_m, var_p}; - return var_rc; + return {rc_Im.at(1), rc_Ip.at(0)}; } } // namespace ReconX