Skip to content

Commit 7512086

Browse files
nilsdeppeclaude
nilsdeppe
and
claude
committed
Add PPM reconstruction to ValenciaDivClean GRMHD
Adds PpmPrim, a piecewise parabolic method (PPM) reconstructor for the GRMHD ValenciaDivClean system, following the same structure as MonotonisedCentralPrim. PPM uses a 3-point stencil with the Colella & Woodward monotonicity limiter (ghost_zone_size = 2, no adaptive order). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
1 parent 2fe49c1 commit 7512086

7 files changed

Lines changed: 345 additions & 1 deletion

File tree

‎src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ spectre_target_sources(
77
MonotonicityPreserving5.cpp
88
MonotonisedCentral.cpp
99
PositivityPreservingAdaptiveOrder.cpp
10+
Ppm.cpp
1011
Reconstructor.cpp
1112
RegisterDerivedWithCharm.cpp
1213
Wcns5z.cpp
@@ -22,6 +23,7 @@ spectre_target_headers(
2223
MonotonicityPreserving5.hpp
2324
MonotonisedCentral.hpp
2425
PositivityPreservingAdaptiveOrder.hpp
26+
Ppm.hpp
2527
Reconstructor.hpp
2628
ReconstructWork.hpp
2729
ReconstructWork.tpp

‎src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Factory.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,5 +6,6 @@
66
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.hpp"
77
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp"
88
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp"
9+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Ppm.hpp"
910
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
1011
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.hpp"

‎src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Ppm.cpp‎

Lines changed: 160 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,160 @@
1+
// Distributed under the MIT License.
2+
// See LICENSE.txt for details.
3+
4+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Ppm.hpp"
5+
6+
#include <array>
7+
#include <cstddef>
8+
#include <memory>
9+
#include <pup.h>
10+
#include <utility>
11+
12+
#include "DataStructures/DataVector.hpp"
13+
#include "DataStructures/Index.hpp"
14+
#include "DataStructures/Tensor/Tensor.hpp"
15+
#include "DataStructures/Variables.hpp"
16+
#include "Domain/Structure/Direction.hpp"
17+
#include "Domain/Structure/DirectionalIdMap.hpp"
18+
#include "Domain/Structure/Element.hpp"
19+
#include "Domain/Structure/ElementId.hpp"
20+
#include "Evolution/DgSubcell/GhostData.hpp"
21+
#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
22+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp"
23+
#include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp"
24+
#include "NumericalAlgorithms/FiniteDifference/Ppm.hpp"
25+
#include "NumericalAlgorithms/Spectral/Mesh.hpp"
26+
#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
27+
#include "Utilities/GenerateInstantiations.hpp"
28+
#include "Utilities/Gsl.hpp"
29+
#include "Utilities/TMPL.hpp"
30+
31+
namespace grmhd::ValenciaDivClean::fd {
32+
PpmPrim::PpmPrim(const bool reconstruct_rho_times_temperature)
33+
: reconstruct_rho_times_temperature_(reconstruct_rho_times_temperature) {}
34+
35+
PpmPrim::PpmPrim(CkMigrateMessage* const msg) : Reconstructor(msg) {}
36+
37+
std::unique_ptr<Reconstructor> PpmPrim::get_clone() const {
38+
return std::make_unique<PpmPrim>(*this);
39+
}
40+
41+
void PpmPrim::pup(PUP::er& p) {
42+
Reconstructor::pup(p);
43+
p | reconstruct_rho_times_temperature_;
44+
}
45+
46+
// NOLINTNEXTLINE
47+
PUP::able::PUP_ID PpmPrim::my_PUP_ID = 0;
48+
49+
template <size_t ThermodynamicDim>
50+
void PpmPrim::reconstruct(
51+
const gsl::not_null<std::array<Variables<tags_list_for_reconstruct>, dim>*>
52+
vars_on_lower_face,
53+
const gsl::not_null<std::array<Variables<tags_list_for_reconstruct>, dim>*>
54+
vars_on_upper_face,
55+
const Variables<hydro::grmhd_tags<DataVector>>& volume_prims,
56+
const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
57+
const Element<dim>& element,
58+
const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>& ghost_data,
59+
const Mesh<dim>& subcell_mesh) const {
60+
DirectionalIdMap<dim, Variables<prims_to_reconstruct_tags>>
61+
neighbor_variables_data{};
62+
::fd::neighbor_data_as_variables<dim>(make_not_null(&neighbor_variables_data),
63+
ghost_data, ghost_zone_size(),
64+
subcell_mesh);
65+
reconstruct_prims_work<prims_to_reconstruct_tags>(
66+
vars_on_lower_face, vars_on_upper_face,
67+
[](auto upper_face_vars_ptr, auto lower_face_vars_ptr,
68+
const auto& volume_vars, const auto& ghost_cell_vars,
69+
const auto& subcell_extents, const size_t number_of_variables) {
70+
::fd::reconstruction::ppm(upper_face_vars_ptr, lower_face_vars_ptr,
71+
volume_vars, ghost_cell_vars, subcell_extents,
72+
number_of_variables);
73+
},
74+
volume_prims, eos, element, neighbor_variables_data, subcell_mesh,
75+
ghost_zone_size(), true, reconstruct_rho_times_temperature());
76+
}
77+
78+
template <size_t ThermodynamicDim>
79+
void PpmPrim::reconstruct_fd_neighbor(
80+
const gsl::not_null<Variables<tags_list_for_reconstruct>*> vars_on_face,
81+
const Variables<hydro::grmhd_tags<DataVector>>& subcell_volume_prims,
82+
const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
83+
const Element<dim>& element,
84+
const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>& ghost_data,
85+
const Mesh<dim>& subcell_mesh,
86+
const Direction<dim> direction_to_reconstruct) const {
87+
reconstruct_fd_neighbor_work<prims_to_reconstruct_tags,
88+
prims_to_reconstruct_tags>(
89+
vars_on_face,
90+
[](const auto tensor_component_on_face_ptr,
91+
const auto& tensor_component_volume,
92+
const auto& tensor_component_neighbor,
93+
const Index<dim>& subcell_extents,
94+
const Index<dim>& ghost_data_extents,
95+
const Direction<dim>& local_direction_to_reconstruct) {
96+
::fd::reconstruction::reconstruct_neighbor<
97+
Side::Lower, ::fd::reconstruction::detail::PpmReconstructor>(
98+
tensor_component_on_face_ptr, tensor_component_volume,
99+
tensor_component_neighbor, subcell_extents, ghost_data_extents,
100+
local_direction_to_reconstruct);
101+
},
102+
[](const auto tensor_component_on_face_ptr,
103+
const auto& tensor_component_volume,
104+
const auto& tensor_component_neighbor,
105+
const Index<dim>& subcell_extents,
106+
const Index<dim>& ghost_data_extents,
107+
const Direction<dim>& local_direction_to_reconstruct) {
108+
::fd::reconstruction::reconstruct_neighbor<
109+
Side::Upper, ::fd::reconstruction::detail::PpmReconstructor>(
110+
tensor_component_on_face_ptr, tensor_component_volume,
111+
tensor_component_neighbor, subcell_extents, ghost_data_extents,
112+
local_direction_to_reconstruct);
113+
},
114+
subcell_volume_prims, eos, element, ghost_data, subcell_mesh,
115+
direction_to_reconstruct, ghost_zone_size(), true,
116+
reconstruct_rho_times_temperature());
117+
}
118+
119+
bool PpmPrim::reconstruct_rho_times_temperature() const {
120+
return reconstruct_rho_times_temperature_;
121+
}
122+
123+
bool operator==(const PpmPrim& lhs, const PpmPrim& rhs) {
124+
return lhs.reconstruct_rho_times_temperature() ==
125+
rhs.reconstruct_rho_times_temperature();
126+
}
127+
128+
bool operator!=(const PpmPrim& lhs, const PpmPrim& rhs) {
129+
return not(lhs == rhs);
130+
}
131+
132+
#define THERMO_DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
133+
134+
#define INSTANTIATION(r, data) \
135+
template void PpmPrim::reconstruct( \
136+
gsl::not_null<std::array<Variables<tags_list_for_reconstruct>, 3>*> \
137+
vars_on_lower_face, \
138+
gsl::not_null<std::array<Variables<tags_list_for_reconstruct>, 3>*> \
139+
vars_on_upper_face, \
140+
const Variables<hydro::grmhd_tags<DataVector>>& volume_prims, \
141+
const EquationsOfState::EquationOfState<true, THERMO_DIM(data)>& eos, \
142+
const Element<3>& element, \
143+
const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \
144+
ghost_data, \
145+
const Mesh<3>& subcell_mesh) const; \
146+
template void PpmPrim::reconstruct_fd_neighbor( \
147+
gsl::not_null<Variables<tags_list_for_reconstruct>*> vars_on_face, \
148+
const Variables<hydro::grmhd_tags<DataVector>>& subcell_volume_prims, \
149+
const EquationsOfState::EquationOfState<true, THERMO_DIM(data)>& eos, \
150+
const Element<3>& element, \
151+
const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \
152+
ghost_data, \
153+
const Mesh<3>& subcell_mesh, \
154+
const Direction<3> direction_to_reconstruct) const;
155+
156+
GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))
157+
158+
#undef INSTANTIATION
159+
#undef THERMO_DIM
160+
} // namespace grmhd::ValenciaDivClean::fd

‎src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Ppm.hpp‎

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
// Distributed under the MIT License.
2+
// See LICENSE.txt for details.
3+
4+
#pragma once
5+
6+
#include <array>
7+
#include <cstddef>
8+
#include <memory>
9+
#include <utility>
10+
11+
#include "DataStructures/DataBox/PrefixHelpers.hpp"
12+
#include "DataStructures/DataBox/Prefixes.hpp"
13+
#include "DataStructures/Tensor/TypeAliases.hpp"
14+
#include "Domain/Structure/DirectionalIdMap.hpp"
15+
#include "Domain/Tags.hpp"
16+
#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
17+
#include "Evolution/DgSubcell/Tags/Mesh.hpp"
18+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp"
19+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
20+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
21+
#include "Options/String.hpp"
22+
#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
23+
#include "PointwiseFunctions/Hydro/Tags.hpp"
24+
#include "Utilities/Serialization/CharmPupable.hpp"
25+
#include "Utilities/TMPL.hpp"
26+
27+
/// \cond
28+
class DataVector;
29+
template <size_t Dim>
30+
class Direction;
31+
template <size_t Dim>
32+
class Element;
33+
template <size_t Dim>
34+
class ElementId;
35+
namespace EquationsOfState {
36+
template <bool IsRelativistic, size_t ThermodynamicDim>
37+
class EquationOfState;
38+
} // namespace EquationsOfState
39+
template <size_t Dim>
40+
class Mesh;
41+
namespace gsl {
42+
template <typename T>
43+
class not_null;
44+
} // namespace gsl
45+
namespace PUP {
46+
class er;
47+
} // namespace PUP
48+
template <typename TagsList>
49+
class Variables;
50+
namespace evolution::dg::subcell {
51+
class GhostData;
52+
} // namespace evolution::dg::subcell
53+
/// \endcond
54+
55+
namespace grmhd::ValenciaDivClean::fd {
56+
/*!
57+
* \brief PPM (Piecewise Parabolic Method) reconstruction. See
58+
* ::fd::reconstruction::ppm() for details.
59+
*/
60+
class PpmPrim : public Reconstructor {
61+
private:
62+
// pressure -> temperature
63+
using prims_to_reconstruct_tags =
64+
tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
65+
hydro::Tags::ElectronFraction<DataVector>,
66+
hydro::Tags::Temperature<DataVector>,
67+
hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>,
68+
hydro::Tags::MagneticField<DataVector, 3>,
69+
hydro::Tags::DivergenceCleaningField<DataVector>>;
70+
71+
public:
72+
static constexpr size_t dim = 3;
73+
74+
struct ReconstructRhoTimesTemperature {
75+
using type = bool;
76+
static constexpr Options::String help = {
77+
"If 'true' then we reconstruct the rho*T, if 'false' we reconstruct "
78+
"T."};
79+
};
80+
81+
using options = tmpl::list<ReconstructRhoTimesTemperature>;
82+
static constexpr Options::String help{
83+
"PPM (Piecewise Parabolic Method) reconstruction scheme using primitive "
84+
"variables."};
85+
86+
PpmPrim() = default;
87+
PpmPrim(PpmPrim&&) = default;
88+
PpmPrim& operator=(PpmPrim&&) = default;
89+
PpmPrim(const PpmPrim&) = default;
90+
PpmPrim& operator=(const PpmPrim&) = default;
91+
~PpmPrim() override = default;
92+
93+
explicit PpmPrim(bool reconstruct_rho_times_temperature);
94+
95+
explicit PpmPrim(CkMigrateMessage* msg);
96+
97+
WRAPPED_PUPable_decl_base_template(Reconstructor, PpmPrim);
98+
99+
auto get_clone() const -> std::unique_ptr<Reconstructor> override;
100+
101+
static constexpr bool use_adaptive_order = false;
102+
103+
void pup(PUP::er& p) override;
104+
105+
size_t ghost_zone_size() const override { return 2; }
106+
107+
using reconstruction_argument_tags =
108+
tmpl::list<::Tags::Variables<hydro::grmhd_tags<DataVector>>,
109+
hydro::Tags::GrmhdEquationOfState, domain::Tags::Element<dim>,
110+
evolution::dg::subcell::Tags::GhostDataForReconstruction<dim>,
111+
evolution::dg::subcell::Tags::Mesh<dim>>;
112+
113+
template <size_t ThermodynamicDim>
114+
void reconstruct(
115+
gsl::not_null<std::array<Variables<tags_list_for_reconstruct>, dim>*>
116+
vars_on_lower_face,
117+
gsl::not_null<std::array<Variables<tags_list_for_reconstruct>, dim>*>
118+
vars_on_upper_face,
119+
const Variables<hydro::grmhd_tags<DataVector>>& volume_prims,
120+
const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
121+
const Element<dim>& element,
122+
const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>&
123+
ghost_data,
124+
const Mesh<dim>& subcell_mesh) const;
125+
126+
/// Called by an element doing DG when the neighbor is doing subcell.
127+
template <size_t ThermodynamicDim>
128+
void reconstruct_fd_neighbor(
129+
gsl::not_null<Variables<tags_list_for_reconstruct>*> vars_on_face,
130+
const Variables<hydro::grmhd_tags<DataVector>>& subcell_volume_prims,
131+
const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
132+
const Element<dim>& element,
133+
const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>&
134+
ghost_data,
135+
const Mesh<dim>& subcell_mesh,
136+
Direction<dim> direction_to_reconstruct) const;
137+
138+
bool reconstruct_rho_times_temperature() const override;
139+
140+
private:
141+
bool reconstruct_rho_times_temperature_{false};
142+
};
143+
144+
bool operator==(const PpmPrim& lhs, const PpmPrim& rhs);
145+
146+
bool operator!=(const PpmPrim& lhs, const PpmPrim& rhs);
147+
} // namespace grmhd::ValenciaDivClean::fd

‎src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@ namespace grmhd::ValenciaDivClean::fd {
1414
class MonotonicityPreserving5Prim;
1515
class MonotonisedCentralPrim;
1616
class PositivityPreservingAdaptiveOrderPrim;
17+
class PpmPrim;
1718
class Wcns5zPrim;
1819
/// \endcond
1920

@@ -36,7 +37,7 @@ class Reconstructor : public PUP::able {
3637

3738
using creatable_classes =
3839
tmpl::list<MonotonicityPreserving5Prim, MonotonisedCentralPrim,
39-
PositivityPreservingAdaptiveOrderPrim, Wcns5zPrim>;
40+
PositivityPreservingAdaptiveOrderPrim, PpmPrim, Wcns5zPrim>;
4041

4142
virtual std::unique_ptr<Reconstructor> get_clone() const = 0;
4243

‎tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ set(LIBRARY_SOURCES
1616
FiniteDifference/Test_MonotonicityPreserving5.cpp
1717
FiniteDifference/Test_MonotonisedCentral.cpp
1818
FiniteDifference/Test_PositivityPreservingAdaptiveOrder.cpp
19+
FiniteDifference/Test_Ppm.cpp
1920
FiniteDifference/Test_Tag.cpp
2021
FiniteDifference/Test_Wcns5z.cpp
2122
Subcell/Test_ComputeFluxes.cpp

‎tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Test_Ppm.cpp‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
// Distributed under the MIT License.
2+
// See LICENSE.txt for details.
3+
4+
#include "Framework/TestingFramework.hpp"
5+
6+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Factory.hpp"
7+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Ppm.hpp"
8+
#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Tag.hpp"
9+
#include "Framework/TestCreation.hpp"
10+
#include "Helpers/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PrimReconstructor.hpp"
11+
12+
SPECTRE_TEST_CASE("Unit.Evolution.Systems.GrMhd.ValenciaDivClean.Fd.PpmPrim",
13+
"[Unit][Evolution]") {
14+
namespace helpers = TestHelpers::grmhd::ValenciaDivClean::fd;
15+
// We only test reconstructing T because for low-order reconstructors this
16+
// fails because the nonlinear terms (quadratics) aren't captured
17+
// accurately enough.
18+
const grmhd::ValenciaDivClean::fd::PpmPrim ppm_recons{false};
19+
helpers::test_prim_reconstructor(5, ppm_recons);
20+
const auto ppm_from_options_base = TestHelpers::test_factory_creation<
21+
grmhd::ValenciaDivClean::fd::Reconstructor,
22+
grmhd::ValenciaDivClean::fd::OptionTags::Reconstructor>(
23+
"PpmPrim:\n"
24+
" ReconstructRhoTimesTemperature: false\n");
25+
auto* const ppm_from_options =
26+
dynamic_cast<const grmhd::ValenciaDivClean::fd::PpmPrim*>(
27+
ppm_from_options_base.get());
28+
REQUIRE(ppm_from_options != nullptr);
29+
CHECK(*ppm_from_options == ppm_recons);
30+
const grmhd::ValenciaDivClean::fd::PpmPrim ppm_recons_recon_T{true};
31+
CHECK_FALSE(ppm_recons_recon_T == ppm_recons);
32+
}

0 commit comments

Comments
 (0)