diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.cpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.cpp index ca08d9a94113..dbdb0096e8a7 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.cpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.cpp @@ -27,6 +27,7 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/FiniteDifference/MonotonisedCentral.hpp" #include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp" #include "NumericalAlgorithms/FiniteDifference/Unlimited.hpp" @@ -67,7 +68,8 @@ void MonotonisedCentralPrim::reconstruct( const EquationsOfState::EquationOfState& eos, const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const { + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const { using prim_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags::primitive_grmhd_reconstruction_tags; using all_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags:: @@ -121,7 +123,8 @@ void MonotonisedCentralPrim::reconstruct( shift, spacetime_metric); }, volume_prims, volume_spacetime_and_cons_vars, eos, element, - neighbor_variables_data, subcell_mesh, ghost_zone_size(), true); + neighbor_variables_data, subcell_mesh, ghost_zone_size(), true, + fix_to_atmosphere); } template @@ -135,6 +138,7 @@ void MonotonisedCentralPrim::reconstruct_fd_neighbor( const Element& element, const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction direction_to_reconstruct) const { using prim_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags::primitive_grmhd_reconstruction_tags; @@ -222,7 +226,7 @@ void MonotonisedCentralPrim::reconstruct_fd_neighbor( }, subcell_volume_prims, subcell_volume_spacetime_metric, eos, element, ghost_data, subcell_mesh, direction_to_reconstruct, ghost_zone_size(), - true); + true, fix_to_atmosphere); } bool operator==(const MonotonisedCentralPrim& /*lhs*/, @@ -250,7 +254,8 @@ bool operator!=(const MonotonisedCentralPrim& lhs, const Element<3>& element, \ const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ - const Mesh<3>& subcell_mesh) const; \ + const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; \ template void MonotonisedCentralPrim::reconstruct_fd_neighbor( \ gsl::not_null*> \ vars_on_face, \ @@ -263,6 +268,7 @@ bool operator!=(const MonotonisedCentralPrim& lhs, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, \ const Direction<3> direction_to_reconstruct) const; GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp index c9eafffbe27a..033bd3d0a0a6 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp @@ -23,6 +23,8 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "Options/String.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -101,7 +103,8 @@ class MonotonisedCentralPrim : public Reconstructor { typename System::variables_tag, hydro::Tags::GrmhdEquationOfState, domain::Tags::Element, evolution::dg::subcell::Tags::GhostDataForReconstruction, - evolution::dg::subcell::Tags::Mesh>; + evolution::dg::subcell::Tags::Mesh, + ::Tags::VariableFixer>>; template void reconstruct( @@ -114,7 +117,8 @@ class MonotonisedCentralPrim : public Reconstructor { const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const; + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; /// Called by an element doing DG when the neighbor is doing subcell. template @@ -129,6 +133,7 @@ class MonotonisedCentralPrim : public Reconstructor { const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction direction_to_reconstruct) const; }; diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp index 6977f4737fb5..c779ff979c3c 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp @@ -29,6 +29,7 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "NumericalAlgorithms/FiniteDifference/MonotonicityPreserving5.hpp" #include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp" @@ -120,7 +121,8 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct( const EquationsOfState::EquationOfState& eos, const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const { + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const { using all_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags:: primitive_grmhd_and_spacetime_reconstruction_tags; @@ -178,15 +180,15 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct( shift, spacetime_metric); }, volume_prims, volume_spacetime_and_cons_vars, eos, element, - neighbor_variables_data, subcell_mesh, ghost_zone_size(), false); + neighbor_variables_data, subcell_mesh, ghost_zone_size(), false, + fix_to_atmosphere); reconstruct_prims_work>, non_positive_tags>( vars_on_lower_face, vars_on_upper_face, - [this]( - auto upper_face_vars_ptr, auto lower_face_vars_ptr, - const auto& volume_vars, const auto& ghost_cell_vars, - const auto& subcell_extents, const size_t number_of_variables) { + [this](auto upper_face_vars_ptr, auto lower_face_vars_ptr, + const auto& volume_vars, const auto& ghost_cell_vars, + const auto& subcell_extents, const size_t number_of_variables) { reconstruct_(upper_face_vars_ptr, lower_face_vars_ptr, volume_vars, ghost_cell_vars, subcell_extents, number_of_variables, four_to_the_alpha_5_, @@ -227,7 +229,8 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct( shift, spacetime_metric); }, volume_prims, volume_spacetime_and_cons_vars, eos, element, - neighbor_variables_data, subcell_mesh, ghost_zone_size(), true); + neighbor_variables_data, subcell_mesh, ghost_zone_size(), true, + fix_to_atmosphere); } // The current implementation does not use positivity-preserving @@ -245,7 +248,8 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct_fd_neighbor( const Element& element, const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, - const Direction direction_to_reconstruct) const { + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, + const Direction& direction_to_reconstruct) const { using prim_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags::primitive_grmhd_reconstruction_tags; using all_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags:: @@ -336,7 +340,7 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct_fd_neighbor( }, subcell_volume_prims, subcell_volume_spacetime_metric, eos, element, ghost_data, subcell_mesh, direction_to_reconstruct, ghost_zone_size(), - true); + true, fix_to_atmosphere); } bool operator==(const PositivityPreservingAdaptiveOrderPrim& lhs, @@ -372,7 +376,8 @@ bool operator!=(const PositivityPreservingAdaptiveOrderPrim& lhs, const Element<3>& element, \ const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ - const Mesh<3>& subcell_mesh) const; \ + const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; \ template void \ PositivityPreservingAdaptiveOrderPrim::reconstruct_fd_neighbor( \ gsl::not_null*> \ @@ -386,7 +391,8 @@ bool operator!=(const PositivityPreservingAdaptiveOrderPrim& lhs, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ const Mesh<3>& subcell_mesh, \ - const Direction<3> direction_to_reconstruct) const; + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, \ + const Direction<3>& direction_to_reconstruct) const; GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp index 633fca930ac9..af09912185f3 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp @@ -20,6 +20,8 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "Options/Auto.hpp" #include "Options/Context.hpp" @@ -159,7 +161,8 @@ class PositivityPreservingAdaptiveOrderPrim : public Reconstructor { typename System::variables_tag, hydro::Tags::GrmhdEquationOfState, domain::Tags::Element, evolution::dg::subcell::Tags::GhostDataForReconstruction, - evolution::dg::subcell::Tags::Mesh>; + evolution::dg::subcell::Tags::Mesh, + ::Tags::VariableFixer>>; template void reconstruct( @@ -174,7 +177,8 @@ class PositivityPreservingAdaptiveOrderPrim : public Reconstructor { const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const; + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; /// Called by an element doing DG when the neighbor is doing subcell. template @@ -189,7 +193,8 @@ class PositivityPreservingAdaptiveOrderPrim : public Reconstructor { const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, - const Direction direction_to_reconstruct) const; + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, + const Direction& direction_to_reconstruct) const; private: // NOLINTNEXTLINE(readability-redundant-declaration) diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.hpp index 4f4c2cc2c1c5..553ad0ec9aaa 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.hpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.hpp @@ -14,6 +14,7 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/TagsDeclarations.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp" #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp" #include "Utilities/TMPL.hpp" @@ -97,7 +98,8 @@ void reconstruct_prims_work( const DirectionalIdMap<3, Variables>& neighbor_data, const Mesh<3>& subcell_mesh, size_t ghost_zone_size, - bool compute_conservatives); + bool compute_conservatives, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere); /*! * \brief Reconstructs \f$\rho, p, Wv^i, B^i\f$, \f$\Phi\f$, the spacetime @@ -132,5 +134,6 @@ void reconstruct_fd_neighbor_work( const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, const Mesh<3>& subcell_mesh, const Direction<3>& direction_to_reconstruct, - size_t ghost_zone_size, bool compute_conservatives); + size_t ghost_zone_size, bool compute_conservative, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere); } // namespace grmhd::GhValenciaDivClean::fd diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.tpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.tpp index 624a39373a6d..a29e880dbaa2 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.tpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/ReconstructWork.tpp @@ -27,6 +27,7 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/ConservativeFromPrimitive.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/FiniteDifference/Reconstruct.tpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" @@ -59,7 +60,8 @@ void reconstruct_prims_work( const DirectionalIdMap<3, Variables>& neighbor_data, const Mesh<3>& subcell_mesh, const size_t ghost_zone_size, - const bool compute_conservatives) { + const bool compute_conservatives, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere) { ASSERT(Mesh<3>(subcell_mesh.extents(0), subcell_mesh.basis(0), subcell_mesh.quadrature(0)) == subcell_mesh, "The subcell mesh should be isotropic but got " << subcell_mesh); @@ -224,9 +226,11 @@ void reconstruct_prims_work( if (compute_conservatives) { ValenciaDivClean::fd::compute_conservatives_for_reconstruction( - make_not_null(&gsl::at(*vars_on_lower_face, i)), eos); + make_not_null(&gsl::at(*vars_on_lower_face, i)), eos, + fix_to_atmosphere); ValenciaDivClean::fd::compute_conservatives_for_reconstruction( - make_not_null(&gsl::at(*vars_on_upper_face, i)), eos); + make_not_null(&gsl::at(*vars_on_upper_face, i)), eos, + fix_to_atmosphere); } } } @@ -255,7 +259,8 @@ void reconstruct_fd_neighbor_work( const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, const Mesh<3>& subcell_mesh, const Direction<3>& direction_to_reconstruct, - const size_t ghost_zone_size, const bool compute_conservatives) { + const size_t ghost_zone_size, const bool compute_conservatives, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere) { const DirectionalId<3> mortar_id{ direction_to_reconstruct, *element.neighbors().at(direction_to_reconstruct).begin()}; @@ -363,8 +368,8 @@ void reconstruct_fd_neighbor_work( spacetime_vars_for_grmhd(vars_on_face); } if (compute_conservatives) { - ValenciaDivClean::fd::compute_conservatives_for_reconstruction(vars_on_face, - eos); + ValenciaDivClean::fd::compute_conservatives_for_reconstruction( + vars_on_face, eos, fix_to_atmosphere); } } } // namespace grmhd::GhValenciaDivClean::fd diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.cpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.cpp index 96c742c1be30..9e0d869c431a 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.cpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.cpp @@ -25,6 +25,7 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp" #include "NumericalAlgorithms/FiniteDifference/Reconstruct.tpp" @@ -92,7 +93,8 @@ void Wcns5zPrim::reconstruct( const EquationsOfState::EquationOfState& eos, const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const { + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const { using all_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags:: primitive_grmhd_and_spacetime_reconstruction_tags; @@ -144,7 +146,8 @@ void Wcns5zPrim::reconstruct( shift, spacetime_metric); }, volume_prims, volume_spacetime_and_cons_vars, eos, element, - neighbor_variables_data, subcell_mesh, ghost_zone_size(), true); + neighbor_variables_data, subcell_mesh, ghost_zone_size(), true, + fix_to_atmosphere); } template @@ -158,6 +161,7 @@ void Wcns5zPrim::reconstruct_fd_neighbor( const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, const Mesh<3>& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction<3>& direction_to_reconstruct) const { using prim_tags_for_reconstruction = grmhd::GhValenciaDivClean::Tags::primitive_grmhd_reconstruction_tags; @@ -242,7 +246,7 @@ void Wcns5zPrim::reconstruct_fd_neighbor( }, subcell_volume_prims, subcell_volume_spacetime_metric, eos, element, ghost_data, subcell_mesh, direction_to_reconstruct, ghost_zone_size(), - true); + true, fix_to_atmosphere); } bool operator==(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs) { @@ -273,7 +277,8 @@ bool operator!=(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs) { const Element<3>& element, \ const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ - const Mesh<3>& subcell_mesh) const; \ + const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; \ template void Wcns5zPrim::reconstruct_fd_neighbor( \ gsl::not_null*> \ vars_on_face, \ @@ -286,6 +291,7 @@ bool operator!=(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs) { const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, \ const Direction<3>& direction_to_reconstruct) const; GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp index 4210cfa46f2c..efebea90b5db 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Wcns5z.hpp @@ -19,6 +19,8 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "Options/Auto.hpp" #include "Options/Context.hpp" @@ -137,7 +139,8 @@ class Wcns5zPrim : public Reconstructor { typename System::variables_tag, hydro::Tags::GrmhdEquationOfState, domain::Tags::Element, evolution::dg::subcell::Tags::GhostDataForReconstruction, - evolution::dg::subcell::Tags::Mesh>; + evolution::dg::subcell::Tags::Mesh, + ::Tags::VariableFixer>>; template void reconstruct( @@ -150,7 +153,8 @@ class Wcns5zPrim : public Reconstructor { const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const; + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; template void reconstruct_fd_neighbor( @@ -164,6 +168,7 @@ class Wcns5zPrim : public Reconstructor { const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction& direction_to_reconstruct) const; private: diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.cpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.cpp index 226276f941e0..cef96c3458fe 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.cpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.cpp @@ -44,6 +44,7 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" @@ -79,6 +80,8 @@ DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( const Mesh<3>& subcell_mesh = db::get>(box); const Mesh<3>& dg_mesh = db::get>(box); + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere = + db::get<::Tags::VariableFixer>>(box); const auto& subcell_options = db::get>(box); const auto& evolved_vars = db::get(box); @@ -112,8 +115,9 @@ DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( call_with_dynamic_type( &base_boundary_correction, [&box, &dg_mesh, &mortars_to_reconstruct_to, &neighbor_package_data, - &ghost_subcell_data, &recons, &subcell_mesh, &subcell_options, - &volume_prims, &volume_spacetime_vars](const auto* gh_grmhd_correction) { + &ghost_subcell_data, &recons, &subcell_mesh, &fix_to_atmosphere, + &subcell_options, &volume_prims, + &volume_spacetime_vars](const auto* gh_grmhd_correction) { using DerivedCorrection = std::decay_t; const auto& boundary_correction = dynamic_cast(*gh_grmhd_correction); @@ -156,13 +160,14 @@ DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( call_with_dynamic_type( - &recons, [&element, &eos, &mortar_id, &ghost_subcell_data, - &subcell_mesh, &vars_on_face, &volume_prims, - &volume_spacetime_vars](const auto& reconstructor) { + &recons, + [&element, &eos, &mortar_id, &ghost_subcell_data, &subcell_mesh, + &fix_to_atmosphere, &vars_on_face, &volume_prims, + &volume_spacetime_vars](const auto& reconstructor) { reconstructor->reconstruct_fd_neighbor( make_not_null(&vars_on_face), volume_prims, volume_spacetime_vars, eos, element, ghost_subcell_data, - subcell_mesh, mortar_id.direction()); + subcell_mesh, fix_to_atmosphere, mortar_id.direction()); }); // Get the mesh velocity if needed diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.cpp index 1d899a39d72b..2c87f1bb3690 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.cpp @@ -20,6 +20,7 @@ #include "Evolution/DgSubcell/GhostData.hpp" #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/FiniteDifference/MonotonicityPreserving5.hpp" #include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" @@ -61,7 +62,8 @@ void MonotonicityPreserving5Prim::reconstruct( const EquationsOfState::EquationOfState& eos, const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const { + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const { DirectionalIdMap> neighbor_variables_data{}; ::fd::neighbor_data_as_variables(make_not_null(&neighbor_variables_data), @@ -79,7 +81,7 @@ void MonotonicityPreserving5Prim::reconstruct( epsilon_); }, volume_prims, eos, element, neighbor_variables_data, subcell_mesh, - ghost_zone_size(), true); + ghost_zone_size(), true, fix_to_atmosphere); } template @@ -90,6 +92,7 @@ void MonotonicityPreserving5Prim::reconstruct_fd_neighbor( const Element& element, const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction direction_to_reconstruct) const { reconstruct_fd_neighbor_work( @@ -121,7 +124,7 @@ void MonotonicityPreserving5Prim::reconstruct_fd_neighbor( local_direction_to_reconstruct, alpha_, epsilon_); }, subcell_volume_prims, eos, element, ghost_data, subcell_mesh, - direction_to_reconstruct, ghost_zone_size(), true); + direction_to_reconstruct, ghost_zone_size(), true, fix_to_atmosphere); } bool operator==(const MonotonicityPreserving5Prim& lhs, @@ -147,7 +150,8 @@ bool operator!=(const MonotonicityPreserving5Prim& lhs, const Element<3>& element, \ const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ - const Mesh<3>& subcell_mesh) const; \ + const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; \ template void MonotonicityPreserving5Prim::reconstruct_fd_neighbor( \ gsl::not_null*> vars_on_face, \ const Variables>& subcell_volume_prims, \ @@ -156,6 +160,7 @@ bool operator!=(const MonotonicityPreserving5Prim& lhs, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, \ const Direction<3> direction_to_reconstruct) const; GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.hpp index 5c80f63c060a..ba983b7f087f 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonicityPreserving5.hpp @@ -20,6 +20,8 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "Options/String.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -121,7 +123,8 @@ class MonotonicityPreserving5Prim : public Reconstructor { tmpl::list<::Tags::Variables>, hydro::Tags::GrmhdEquationOfState, domain::Tags::Element, evolution::dg::subcell::Tags::GhostDataForReconstruction, - evolution::dg::subcell::Tags::Mesh>; + evolution::dg::subcell::Tags::Mesh, + ::Tags::VariableFixer>>; template void reconstruct( @@ -134,7 +137,8 @@ class MonotonicityPreserving5Prim : public Reconstructor { const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const; + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; template void reconstruct_fd_neighbor( @@ -145,6 +149,7 @@ class MonotonicityPreserving5Prim : public Reconstructor { const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction direction_to_reconstruct) const; private: diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.cpp index f81f46e18271..d5cd167e15fc 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.cpp @@ -20,6 +20,7 @@ #include "Evolution/DgSubcell/GhostData.hpp" #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/FiniteDifference/MonotonisedCentral.hpp" #include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" @@ -51,7 +52,8 @@ void MonotonisedCentralPrim::reconstruct( const EquationsOfState::EquationOfState& eos, const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const { + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const { DirectionalIdMap> neighbor_variables_data{}; ::fd::neighbor_data_as_variables(make_not_null(&neighbor_variables_data), @@ -67,7 +69,7 @@ void MonotonisedCentralPrim::reconstruct( ghost_cell_vars, subcell_extents, number_of_variables); }, volume_prims, eos, element, neighbor_variables_data, subcell_mesh, - ghost_zone_size(), true); + ghost_zone_size(), true, fix_to_atmosphere); } template @@ -78,6 +80,7 @@ void MonotonisedCentralPrim::reconstruct_fd_neighbor( const Element& element, const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction direction_to_reconstruct) const { reconstruct_fd_neighbor_work( @@ -109,7 +112,7 @@ void MonotonisedCentralPrim::reconstruct_fd_neighbor( local_direction_to_reconstruct); }, subcell_volume_prims, eos, element, ghost_data, subcell_mesh, - direction_to_reconstruct, ghost_zone_size(), true); + direction_to_reconstruct, ghost_zone_size(), true, fix_to_atmosphere); } bool operator==(const MonotonisedCentralPrim& /*lhs*/, @@ -135,7 +138,8 @@ bool operator!=(const MonotonisedCentralPrim& lhs, const Element<3>& element, \ const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ - const Mesh<3>& subcell_mesh) const; \ + const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; \ template void MonotonisedCentralPrim::reconstruct_fd_neighbor( \ gsl::not_null*> vars_on_face, \ const Variables>& subcell_volume_prims, \ @@ -144,6 +148,7 @@ bool operator!=(const MonotonisedCentralPrim& lhs, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, \ const Direction<3> direction_to_reconstruct) const; GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp index dfa95aecf38f..75599aed2752 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/MonotonisedCentral.hpp @@ -18,6 +18,8 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "Options/String.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -98,7 +100,8 @@ class MonotonisedCentralPrim : public Reconstructor { tmpl::list<::Tags::Variables>, hydro::Tags::GrmhdEquationOfState, domain::Tags::Element, evolution::dg::subcell::Tags::GhostDataForReconstruction, - evolution::dg::subcell::Tags::Mesh>; + evolution::dg::subcell::Tags::Mesh, + ::Tags::VariableFixer>>; template void reconstruct( @@ -111,7 +114,8 @@ class MonotonisedCentralPrim : public Reconstructor { const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const; + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; /// Called by an element doing DG when the neighbor is doing subcell. template @@ -123,6 +127,7 @@ class MonotonisedCentralPrim : public Reconstructor { const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction direction_to_reconstruct) const; }; diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp index 6a9ecd977d33..29e05001e5c6 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.cpp @@ -24,6 +24,7 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp" #include "NumericalAlgorithms/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp" @@ -106,7 +107,8 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct( const EquationsOfState::EquationOfState& eos, const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, - const Mesh<3>& subcell_mesh) const { + const Mesh<3>& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const { DirectionalIdMap> neighbor_variables_data{}; ::fd::neighbor_data_as_variables(make_not_null(&neighbor_variables_data), @@ -129,7 +131,7 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct( std::numeric_limits::signaling_NaN())); }, volume_prims, eos, element, neighbor_variables_data, subcell_mesh, - ghost_zone_size(), false); + ghost_zone_size(), false, fix_to_atmosphere); reconstruct_prims_work( vars_on_lower_face, vars_on_upper_face, [this](auto upper_face_vars_ptr, auto lower_face_vars_ptr, @@ -144,7 +146,7 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct( std::numeric_limits::signaling_NaN())); }, volume_prims, eos, element, neighbor_variables_data, subcell_mesh, - ghost_zone_size(), true); + ghost_zone_size(), true, fix_to_atmosphere); } template @@ -155,6 +157,7 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct_fd_neighbor( const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, const Mesh<3>& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction<3> direction_to_reconstruct) const { reconstruct_fd_neighbor_work( @@ -190,7 +193,7 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct_fd_neighbor( std::numeric_limits::signaling_NaN())); }, subcell_volume_prims, eos, element, ghost_data, subcell_mesh, - direction_to_reconstruct, ghost_zone_size(), false); + direction_to_reconstruct, ghost_zone_size(), false, fix_to_atmosphere); reconstruct_fd_neighbor_work( vars_on_face, [this](const auto tensor_component_on_face_ptr, @@ -224,7 +227,7 @@ void PositivityPreservingAdaptiveOrderPrim::reconstruct_fd_neighbor( std::numeric_limits::signaling_NaN())); }, subcell_volume_prims, eos, element, ghost_data, subcell_mesh, - direction_to_reconstruct, ghost_zone_size(), true); + direction_to_reconstruct, ghost_zone_size(), true, fix_to_atmosphere); } bool operator==(const PositivityPreservingAdaptiveOrderPrim& lhs, @@ -258,7 +261,8 @@ bool operator!=(const PositivityPreservingAdaptiveOrderPrim& lhs, const Element<3>& element, \ const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ - const Mesh<3>& subcell_mesh) const; \ + const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; \ template void \ PositivityPreservingAdaptiveOrderPrim::reconstruct_fd_neighbor( \ gsl::not_null*> vars_on_face, \ @@ -268,6 +272,7 @@ bool operator!=(const PositivityPreservingAdaptiveOrderPrim& lhs, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, \ const Direction<3> direction_to_reconstruct) const; GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp index 4a3fb2abc771..c0a3d08f2697 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PositivityPreservingAdaptiveOrder.hpp @@ -16,6 +16,8 @@ #include "Evolution/DgSubcell/Tags/Mesh.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "Options/Auto.hpp" #include "Options/Context.hpp" @@ -157,7 +159,8 @@ class PositivityPreservingAdaptiveOrderPrim : public Reconstructor { tmpl::list<::Tags::Variables>, hydro::Tags::GrmhdEquationOfState, domain::Tags::Element, evolution::dg::subcell::Tags::GhostDataForReconstruction, - evolution::dg::subcell::Tags::Mesh>; + evolution::dg::subcell::Tags::Mesh, + ::Tags::VariableFixer>>; template void reconstruct( @@ -172,7 +175,8 @@ class PositivityPreservingAdaptiveOrderPrim : public Reconstructor { const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const; + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; template void reconstruct_fd_neighbor( @@ -183,6 +187,7 @@ class PositivityPreservingAdaptiveOrderPrim : public Reconstructor { const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction direction_to_reconstruct) const; private: diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp index ba9b6f7999bb..79be81294ab9 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp @@ -11,6 +11,7 @@ #include "DataStructures/Tensor/TypeAliases.hpp" #include "Domain/Structure/DirectionalIdMap.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/TagsDeclarations.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp" #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp" #include "Utilities/TMPL.hpp" @@ -106,10 +107,11 @@ void reconstruct_prims_work( const DirectionalIdMap<3, Variables>& neighbor_data, const Mesh<3>& subcell_mesh, size_t ghost_zone_size, - bool compute_conservatives); + bool compute_conservatives, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere); /*! - * \brief Reconstructs the `PrimTagsForReconstruction` and if + * \brief Reconstructs the `PrimTagsForRecons truction` and if * `compute_conservatives` is `true` computes the Lorentz factor, upper spatial * velocity, specific internal energy, and the conserved variables. * @@ -133,5 +135,6 @@ void reconstruct_fd_neighbor_work( const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, const Mesh<3>& subcell_mesh, const Direction<3>& direction_to_reconstruct, - const size_t ghost_zone_size, bool compute_conservatives); + size_t ghost_zone_size, bool compute_conservatives, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere); } // namespace grmhd::ValenciaDivClean::fd diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp index 738179d5188a..74b421cd74f7 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp @@ -17,14 +17,15 @@ #include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Variables.hpp" #include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionMap.hpp" #include "Domain/Structure/DirectionalId.hpp" #include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/DirectionMap.hpp" #include "Domain/Structure/Element.hpp" #include "Domain/Structure/ElementId.hpp" #include "Evolution/DgSubcell/GhostData.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/ConservativeFromPrimitive.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" @@ -36,7 +37,8 @@ namespace grmhd::ValenciaDivClean::fd { template void compute_conservatives_for_reconstruction( const gsl::not_null*> vars_on_face, - const EquationsOfState::EquationOfState& eos) { + const EquationsOfState::EquationOfState& eos, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere) { // Computes: // 1. W v^i // 2. Lorentz factor as sqrt(1 + Wv^i Wv^j\gamma_{ij}) @@ -71,21 +73,33 @@ void compute_conservatives_for_reconstruction( } // pointers to primitive variables - const auto& rest_mass_density = + auto& rest_mass_density = get>(*vars_on_face); - const auto& electron_fraction = + auto& electron_fraction = get>(*vars_on_face); auto& pressure = get>(*vars_on_face); auto& specific_internal_energy = get>(*vars_on_face); - const auto& temperature = - get>(*vars_on_face); + auto& temperature = get>(*vars_on_face); + // The variables have to be fixed to atmosphere + // before the EoS calls, because the EoSs will generically + // fail for bad values of the rest_mass_density, temperature, and + // electron_fraction + fix_to_atmosphere(make_not_null(&rest_mass_density), + make_not_null(&specific_internal_energy), + make_not_null(&spatial_velocity), + make_not_null(&lorentz_factor), make_not_null(&pressure), + make_not_null(&temperature), + make_not_null(&electron_fraction), spatial_metric, eos); + ASSERT(min(rest_mass_density.get()) >= 0.0, + "rest mass density should be positive"); + ASSERT(min(electron_fraction.get()) >= 0.0, + "electron fraction should be positive"); // EoS calls based on reconstructed primitives if constexpr (ThermodynamicDim == 1) { specific_internal_energy = - eos.specific_internal_energy_from_density( - rest_mass_density); + eos.specific_internal_energy_from_density(rest_mass_density); pressure = eos.pressure_from_density(rest_mass_density); } else if constexpr (ThermodynamicDim == 2) { specific_internal_energy = @@ -134,7 +148,8 @@ void reconstruct_prims_work( const DirectionalIdMap<3, Variables>& neighbor_data, const Mesh<3>& subcell_mesh, const size_t ghost_zone_size, - const bool compute_conservatives) { + const bool compute_conservatives, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere) { ASSERT(Mesh<3>(subcell_mesh.extents(0), subcell_mesh.basis(0), subcell_mesh.quadrature(0)) == subcell_mesh, "The subcell mesh should be isotropic but got " << subcell_mesh); @@ -229,9 +244,11 @@ void reconstruct_prims_work( for (size_t i = 0; compute_conservatives and i < 3; ++i) { compute_conservatives_for_reconstruction( - make_not_null(&gsl::at(*vars_on_lower_face, i)), eos); + make_not_null(&gsl::at(*vars_on_lower_face, i)), eos, + fix_to_atmosphere); compute_conservatives_for_reconstruction( - make_not_null(&gsl::at(*vars_on_upper_face, i)), eos); + make_not_null(&gsl::at(*vars_on_upper_face, i)), eos, + fix_to_atmosphere); } } @@ -245,7 +262,8 @@ void reconstruct_fd_neighbor_work( const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, const Mesh<3>& subcell_mesh, const Direction<3>& direction_to_reconstruct, - const size_t ghost_zone_size, const bool compute_conservatives) { + const size_t ghost_zone_size, const bool compute_conservatives, + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere) { const DirectionalId<3> mortar_id{ direction_to_reconstruct, *element.neighbors().at(direction_to_reconstruct).begin()}; @@ -315,7 +333,8 @@ void reconstruct_fd_neighbor_work( }); if (compute_conservatives) { - compute_conservatives_for_reconstruction(vars_on_face, eos); + compute_conservatives_for_reconstruction(vars_on_face, eos, + fix_to_atmosphere); } } } // namespace grmhd::ValenciaDivClean::fd diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.cpp index 672b54dcc053..1f5e14e95f69 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.cpp @@ -24,6 +24,7 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.tpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "NumericalAlgorithms/FiniteDifference/NeighborDataAsVariables.hpp" #include "NumericalAlgorithms/FiniteDifference/Reconstruct.tpp" @@ -84,7 +85,8 @@ void Wcns5zPrim::reconstruct( const EquationsOfState::EquationOfState& eos, const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, - const Mesh<3>& subcell_mesh) const { + const Mesh<3>& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const { DirectionalIdMap> neighbor_variables_data{}; ::fd::neighbor_data_as_variables(make_not_null(&neighbor_variables_data), @@ -101,7 +103,7 @@ void Wcns5zPrim::reconstruct( epsilon_, max_number_of_extrema_); }, volume_prims, eos, element, neighbor_variables_data, subcell_mesh, - ghost_zone_size(), true); + ghost_zone_size(), true, fix_to_atmosphere); } template @@ -112,6 +114,7 @@ void Wcns5zPrim::reconstruct_fd_neighbor( const Element<3>& element, const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& ghost_data, const Mesh<3>& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction<3> direction_to_reconstruct) const { reconstruct_fd_neighbor_work( @@ -139,7 +142,7 @@ void Wcns5zPrim::reconstruct_fd_neighbor( local_direction_to_reconstruct, epsilon_, max_number_of_extrema_); }, subcell_volume_prims, eos, element, ghost_data, subcell_mesh, - direction_to_reconstruct, ghost_zone_size(), true); + direction_to_reconstruct, ghost_zone_size(), true, fix_to_atmosphere); } bool operator==(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs) { @@ -168,7 +171,8 @@ bool operator!=(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs) { const Element<3>& element, \ const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ - const Mesh<3>& subcell_mesh) const; \ + const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; \ template void Wcns5zPrim::reconstruct_fd_neighbor( \ gsl::not_null*> vars_on_face, \ const Variables>& subcell_volume_prims, \ @@ -177,6 +181,7 @@ bool operator!=(const Wcns5zPrim& lhs, const Wcns5zPrim& rhs) { const DirectionalIdMap<3, evolution::dg::subcell::GhostData>& \ ghost_data, \ const Mesh<3>& subcell_mesh, \ + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, \ const Direction<3> direction_to_reconstruct) const; GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.hpp index 9ee355610ac4..d3de3afee66f 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Wcns5z.hpp @@ -16,6 +16,8 @@ #include "Evolution/DgSubcell/Tags/Mesh.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "Options/String.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -131,7 +133,8 @@ class Wcns5zPrim : public Reconstructor { tmpl::list<::Tags::Variables>, hydro::Tags::GrmhdEquationOfState, domain::Tags::Element, evolution::dg::subcell::Tags::GhostDataForReconstruction, - evolution::dg::subcell::Tags::Mesh>; + evolution::dg::subcell::Tags::Mesh, + ::Tags::VariableFixer>>; template void reconstruct( @@ -144,7 +147,8 @@ class Wcns5zPrim : public Reconstructor { const Element& element, const DirectionalIdMap& ghost_data, - const Mesh& subcell_mesh) const; + const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere) const; template void reconstruct_fd_neighbor( @@ -155,6 +159,7 @@ class Wcns5zPrim : public Reconstructor { const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, + const VariableFixing::FixToAtmosphere& fix_to_atmosphere, const Direction direction_to_reconstruct) const; private: diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.cpp index 79140f48ac1d..80325692aa4b 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.cpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.cpp @@ -77,6 +77,8 @@ DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( const Mesh<3>& subcell_mesh = db::get>(box); const Mesh<3>& dg_mesh = db::get>(box); + const VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere = + db::get<::Tags::VariableFixer>>(box); const auto& subcell_options = db::get>(box); @@ -99,7 +101,8 @@ DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( &neighbor_package_data, &ghost_subcell_data, &recons, &subcell_mesh, &subcell_options, - &volume_prims]( + &volume_prims, + &fix_to_atmosphere]( auto derived_correction_v) { using DerivedCorrection = tmpl::type_from; if (typeid(boundary_correction) == typeid(DerivedCorrection)) { @@ -154,12 +157,13 @@ DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( call_with_dynamic_type( - &recons, - [&element, &eos, &mortar_id, &ghost_subcell_data, &subcell_mesh, - &vars_on_face, &volume_prims](const auto& reconstructor) { + &recons, [&element, &eos, &mortar_id, &ghost_subcell_data, + &subcell_mesh, &vars_on_face, &volume_prims, + &fix_to_atmosphere](const auto& reconstructor) { reconstructor->reconstruct_fd_neighbor( make_not_null(&vars_on_face), volume_prims, eos, element, - ghost_subcell_data, subcell_mesh, mortar_id.direction()); + ghost_subcell_data, subcell_mesh, fix_to_atmosphere, + mortar_id.direction()); }); // Get the mesh velocity if needed diff --git a/src/Evolution/VariableFixing/FixToAtmosphere.cpp b/src/Evolution/VariableFixing/FixToAtmosphere.cpp index 5a50506c6b2f..6481994d1927 100644 --- a/src/Evolution/VariableFixing/FixToAtmosphere.cpp +++ b/src/Evolution/VariableFixing/FixToAtmosphere.cpp @@ -11,6 +11,7 @@ #include "Options/ParseError.hpp" #include "Utilities/GenerateInstantiations.hpp" +#include "Parallel/Printf/Printf.hpp" namespace VariableFixing { template @@ -63,7 +64,7 @@ void FixToAtmosphere::operator()( const gsl::not_null*> lorentz_factor, const gsl::not_null*> pressure, const gsl::not_null*> temperature, - const Scalar& electron_fraction, + const gsl::not_null*> electron_fraction, const tnsr::ii& spatial_metric, const EquationsOfState::EquationOfState& equation_of_state) const { @@ -84,12 +85,12 @@ void FixToAtmosphere::operator()( // For 2D & 3D EoS, we also need to limit the temperature / energy if constexpr (ThermodynamicDim > 1) { - bool changed_temperature = false; + bool recompute_eos_dependent_quantities=false; if (const double min_temperature = equation_of_state.temperature_lower_bound(); get(*temperature)[i] < min_temperature) { get(*temperature)[i] = min_temperature; - changed_temperature = true; + recompute_eos_dependent_quantities = true; } // We probably need a better maximum temperature as well, but this is not @@ -98,10 +99,24 @@ void FixToAtmosphere::operator()( equation_of_state.temperature_upper_bound(); get(*temperature)[i] > max_temperature) { get(*temperature)[i] = max_temperature; - changed_temperature = true; + recompute_eos_dependent_quantities = true; + } + // For 3D we also need to fix the electron fraction + if constexpr (ThermodynamicDim == 3) { + if (get(*electron_fraction)[i] > + equation_of_state.electron_fraction_upper_bound() or + get(*electron_fraction)[i] < + equation_of_state.electron_fraction_lower_bound()) { + get(*electron_fraction)[i] = + get(equation_of_state + .equilibrium_electron_fraction_from_density_temperature( + Scalar{get(*rest_mass_density)[i]}, + Scalar{get(*temperature)[i]})); + recompute_eos_dependent_quantities = true; + } } - if (changed_temperature) { + if (recompute_eos_dependent_quantities) { if constexpr (ThermodynamicDim == 2) { specific_internal_energy->get()[i] = get(equation_of_state @@ -118,12 +133,12 @@ void FixToAtmosphere::operator()( .specific_internal_energy_from_density_and_temperature( Scalar{rest_mass_density->get()[i]}, Scalar{get(*temperature)[i]}, - Scalar{get(electron_fraction)[i]})); + Scalar{get(*electron_fraction)[i]})); pressure->get()[i] = get(equation_of_state.pressure_from_density_and_temperature( Scalar{rest_mass_density->get()[i]}, Scalar{temperature->get()[i]}, - Scalar{get(electron_fraction)[i]})); + Scalar{get(*electron_fraction)[i]})); } } } @@ -137,14 +152,16 @@ void FixToAtmosphere::set_density_to_atmosphere( const gsl::not_null*> specific_internal_energy, const gsl::not_null*> temperature, const gsl::not_null*> pressure, - const Scalar& electron_fraction, + const gsl::not_null*> electron_fraction, const EquationsOfState::EquationOfState& equation_of_state, const size_t grid_index) const { const Scalar atmosphere_density{density_of_atmosphere_}; rest_mass_density->get()[grid_index] = get(atmosphere_density); get(*temperature)[grid_index] = equation_of_state.temperature_lower_bound(); - + get(*electron_fraction)[grid_index] = get( + equation_of_state.equilibrium_electron_fraction_from_density_temperature( + atmosphere_density, Scalar{get(*temperature)[grid_index]})); if constexpr (ThermodynamicDim == 1) { pressure->get()[grid_index] = get(equation_of_state.pressure_from_density(atmosphere_density)); @@ -166,14 +183,14 @@ void FixToAtmosphere::set_density_to_atmosphere( specific_internal_energy->get()[grid_index] = get(equation_of_state .specific_internal_energy_from_density_and_temperature( - Scalar{get(*rest_mass_density)[grid_index]}, - Scalar{get(*temperature)[grid_index]}, - Scalar{get(electron_fraction)[grid_index]})); + Scalar{atmosphere_density}, + Scalar{atmosphere_temperature}, + Scalar{get(*electron_fraction)[grid_index]})); pressure->get()[grid_index] = get(equation_of_state.pressure_from_density_and_temperature( - Scalar{get(*rest_mass_density)[grid_index]}, - Scalar{get(*temperature)[grid_index]}, - Scalar{get(electron_fraction)[grid_index]})); + Scalar{atmosphere_density}, + Scalar{atmosphere_temperature}, + Scalar{get(*electron_fraction)[grid_index]})); } } } @@ -241,6 +258,8 @@ GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) #undef INSTANTIATION +#undef INSTANTIATION + #define INSTANTIATION(r, data) \ template void FixToAtmosphere::operator()( \ const gsl::not_null*> rest_mass_density, \ @@ -250,7 +269,7 @@ GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) const gsl::not_null*> lorentz_factor, \ const gsl::not_null*> pressure, \ const gsl::not_null*> temperature, \ - const Scalar& electron_fraction, \ + const gsl::not_null*> electron_fraction, \ const tnsr::ii& spatial_metric, \ const EquationsOfState::EquationOfState& \ equation_of_state) const; diff --git a/src/Evolution/VariableFixing/FixToAtmosphere.hpp b/src/Evolution/VariableFixing/FixToAtmosphere.hpp index a6878dcadcf0..ef9c4b313716 100644 --- a/src/Evolution/VariableFixing/FixToAtmosphere.hpp +++ b/src/Evolution/VariableFixing/FixToAtmosphere.hpp @@ -128,9 +128,9 @@ class FixToAtmosphere { hydro::Tags::SpatialVelocity, hydro::Tags::LorentzFactor, hydro::Tags::Pressure, - hydro::Tags::Temperature>; - using argument_tags = tmpl::list, - gr::Tags::SpatialMetric, + hydro::Tags::Temperature, + hydro::Tags::ElectronFraction>; + using argument_tags = tmpl::list, hydro::Tags::EquationOfStateBase>; // for use in `db::mutate_apply` @@ -143,7 +143,7 @@ class FixToAtmosphere { gsl::not_null*> lorentz_factor, gsl::not_null*> pressure, gsl::not_null*> temperature, - const Scalar& electron_fraction, + gsl::not_null*> electron_fraction, const tnsr::ii& spatial_metric, const EquationsOfState::EquationOfState& equation_of_state) const; @@ -155,7 +155,7 @@ class FixToAtmosphere { gsl::not_null*> specific_internal_energy, gsl::not_null*> temperature, gsl::not_null*> pressure, - const Scalar& electron_fraction, + gsl::not_null*> electron_fraction, const EquationsOfState::EquationOfState& equation_of_state, size_t grid_index) const; diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp index 451bba9d66b9..a1d713a1f651 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp @@ -61,6 +61,8 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/ConservativeFromPrimitive.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "Parallel/Tags/Metavariables.hpp" @@ -263,6 +265,7 @@ double test(const size_t num_dg_pts) { evolution::Tags::BoundaryCorrection< grmhd::GhValenciaDivClean::System>, hydro::Tags::GrmhdEquationOfState, + ::Tags::VariableFixer>, typename System::primitive_variables_tag, dt_variables_tag, variables_tag, evolution::dg::subcell::Tags::GhostDataForReconstruction<3>, @@ -305,7 +308,9 @@ double test(const size_t num_dg_pts) { std::unique_ptr< grmhd::GhValenciaDivClean::BoundaryCorrections::BoundaryCorrection>{ std::make_unique()}, - soln.equation_of_state().promote_to_3d_eos(), dg_prim_vars, + soln.equation_of_state().promote_to_3d_eos(), + VariableFixing::FixToAtmosphere<3>{1.0e-16, 1.1e-16, 1.0e-15, 1.0}, + dg_prim_vars, // Set incorrect size for dt variables because they should get resized. Variables{}, initial_variables, neighbor_data, 1.0, std::move(element_map), moving_mesh_map.get_clone(), diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp index 50a3ab6c3bc0..ba3c462f321d 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_TimeDerivative.cpp @@ -61,6 +61,8 @@ #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/ConservativeFromPrimitive.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" @@ -613,6 +615,7 @@ double test(const size_t num_dg_pts, std::optional expansion_velocity, evolution::Tags::BoundaryCorrection< grmhd::GhValenciaDivClean::System>, hydro::Tags::GrmhdEquationOfState, + ::Tags::VariableFixer>, typename System::primitive_variables_tag, dt_variables_tag, variables_tag, evolution::dg::subcell::Tags::GhostDataForReconstruction<3>, @@ -675,7 +678,9 @@ double test(const size_t num_dg_pts, std::optional expansion_velocity, std::make_unique, ValenciaDivClean::BoundaryCorrections::Hll>>()}, - soln.equation_of_state().promote_to_3d_eos(), cell_centered_prim_vars, + soln.equation_of_state().promote_to_3d_eos(), + VariableFixing::FixToAtmosphere<3>{1.0e-16, 1.1e-16, 1.0e-15, 1.0}, + cell_centered_prim_vars, // Set incorrect size for dt variables because they should get resized. Variables{}, initial_variables, neighbor_data, dummy_reconstruction_order, 1.0, mortar_data, diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp index 1a6a2650f8e1..203b20675057 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp @@ -49,6 +49,8 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "PointwiseFunctions/AnalyticSolutions/GrMhd/BondiMichel.hpp" @@ -249,6 +251,7 @@ double test(const size_t num_dg_pts) { evolution::dg::subcell::Tags::Mesh<3>, fd::Tags::Reconstructor, evolution::Tags::BoundaryCorrection, hydro::Tags::GrmhdEquationOfState, + ::Tags::VariableFixer>, typename System::spacetime_variables_tag, typename System::primitive_variables_tag, variables_tag, evolution::dg::subcell::Tags::OnSubcellFaces< @@ -268,8 +271,9 @@ double test(const size_t num_dg_pts) { grmhd::ValenciaDivClean::BoundaryCorrections::BoundaryCorrection>{ std::make_unique< grmhd::ValenciaDivClean::BoundaryCorrections::Hll>()}, - soln.equation_of_state().promote_to_3d_eos(), dg_spacetime_vars, - dg_prim_vars, + soln.equation_of_state().promote_to_3d_eos(), + VariableFixing::FixToAtmosphere<3>{1.0e-16, 1.1e-16, 1.0e-15, 1.0}, + dg_spacetime_vars, dg_prim_vars, typename variables_tag::type{dg_mesh.number_of_grid_points()}, face_centered_gr_tags(subcell_mesh, time, element_map, moving_mesh_map, functions_of_time, soln), diff --git a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_TimeDerivative.cpp b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_TimeDerivative.cpp index 8b62143265ae..7480910310a7 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_TimeDerivative.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/Test_TimeDerivative.cpp @@ -57,6 +57,8 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" #include "Evolution/TagsDomain.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" +#include "Evolution/VariableFixing/Tags.hpp" #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp" #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" @@ -382,6 +384,7 @@ std::array test(const size_t num_dg_pts, domain::Tags::Mesh<3>, fd::Tags::Reconstructor, evolution::Tags::BoundaryCorrection, hydro::Tags::GrmhdEquationOfState, + ::Tags::VariableFixer>, typename System::spacetime_variables_tag, typename System::primitive_variables_tag, dt_variables_tag, variables_tag, @@ -450,6 +453,7 @@ std::array test(const size_t num_dg_pts, std::make_unique< grmhd::ValenciaDivClean::BoundaryCorrections::Hll>()}, soln.equation_of_state().promote_to_3d_eos(), + VariableFixing::FixToAtmosphere<3>{1.0e-16, 1.1e-16, 1.0e-15, 1.0}, cell_centered_spacetime_vars, cell_centered_prim_vars, Variables{ subcell_mesh.number_of_grid_points()}, diff --git a/tests/Unit/Evolution/VariableFixing/Test_FixToAtmosphere.cpp b/tests/Unit/Evolution/VariableFixing/Test_FixToAtmosphere.cpp index 8efea9ae439f..dc7f653cbcad 100644 --- a/tests/Unit/Evolution/VariableFixing/Test_FixToAtmosphere.cpp +++ b/tests/Unit/Evolution/VariableFixing/Test_FixToAtmosphere.cpp @@ -28,8 +28,7 @@ void test_variable_fixer( auto specific_internal_energy = equation_of_state.specific_internal_energy_from_density(density); auto temperature = equation_of_state.temperature_from_density(density); - const Scalar electron_fraction{ - DataVector{get(density).size(), 0.5}}; + Scalar electron_fraction{DataVector{get(density).size(), 0.5}}; Scalar lorentz_factor{ DataVector{5.0 / 3.0, 7.0710678118654752, 1.8898223650461359}}; @@ -42,7 +41,7 @@ void test_variable_fixer( spatial_metric.get(i, i) = 2.0; } variable_fixer(&density, &specific_internal_energy, &spatial_velocity, - &lorentz_factor, &pressure, &temperature, electron_fraction, + &lorentz_factor, &pressure, &temperature, &electron_fraction, spatial_metric, equation_of_state); Scalar expected_density{DataVector{1.e-12, 2.e-11, 4.e-12}}; @@ -84,8 +83,7 @@ void test_variable_fixer( density, specific_internal_energy); auto temperature = equation_of_state.temperature_from_density_and_energy( density, specific_internal_energy); - const Scalar electron_fraction{ - DataVector{get(density).size(), 0.5}}; + Scalar electron_fraction{DataVector{get(density).size(), 0.5}}; Scalar lorentz_factor{DataVector{ 5. / 3., 7.0710678118654752, 1.8898223650461359, 7.0710678118654752}}; @@ -100,7 +98,7 @@ void test_variable_fixer( spatial_metric.get(i, i) = 2.; } variable_fixer(&density, &specific_internal_energy, &spatial_velocity, - &lorentz_factor, &pressure, &temperature, electron_fraction, + &lorentz_factor, &pressure, &temperature, &electron_fraction, spatial_metric, equation_of_state); Scalar expected_density{ @@ -140,6 +138,79 @@ void test_variable_fixer( CHECK_ITERABLE_APPROX(spatial_velocity, expected_spatial_velocity); } +template +void test_variable_fixer( + const VariableFixing::FixToAtmosphere& variable_fixer, + const EquationsOfState::EquationOfState& equation_of_state) { + Scalar density{DataVector{2.e-12, 2.e-11, 4.e-12, 2.e-11}}; + Scalar temperature{DataVector{2.0, 1.0, 2.0, 1.0}}; + + const Scalar bounded_electron_fraction{ + DataVector{.3, .2, 1.0, 0.0}}; + auto pressure = equation_of_state.pressure_from_density_and_temperature( + density, temperature, bounded_electron_fraction); + auto specific_internal_energy = + equation_of_state.specific_internal_energy_from_density_and_temperature( + density, temperature, bounded_electron_fraction); + Scalar electron_fraction{DataVector{.3, .2, 1.6, -.4}}; + + Scalar lorentz_factor{DataVector{ + 5. / 3., 7.0710678118654752, 1.8898223650461359, 7.0710678118654752}}; + CHECK(get(lorentz_factor).size() == get(density).size()); + auto spatial_velocity = + make_with_value>(density, 0.); + spatial_velocity.get(0) = DataVector{0.8, 0.7, 0.6, 0.7}; + CHECK(spatial_velocity.get(0).size() == get(density).size()); + auto spatial_metric = + make_with_value>(density, 0.); + for (size_t i = 0; i < Dim; ++i) { + spatial_metric.get(i, i) = 2.; + } + variable_fixer(&density, &specific_internal_energy, &spatial_velocity, + &lorentz_factor, &pressure, &temperature, &electron_fraction, + spatial_metric, equation_of_state); + + Scalar expected_density{ + DataVector{1.e-12, 2.e-11, 4.e-12, 2.e-11}}; + // Temperature reset to zero in atmosphere + const Scalar expected_temperature{DataVector{0.0, 1.0, 2.0, 1.0}}; + // The [0] component is atmosphere, value reset to equilibrium + // The [2] and [3] components are above and below the allowed values + // so are reset to "equalibrium" which in this case is just .1 + const Scalar expected_electron_fraction{ + DataVector{.1, .2, .1, .1}}; + auto expected_pressure = equation_of_state.pressure_from_density_and_energy( + expected_density, expected_temperature, expected_electron_fraction); + auto expected_specific_internal_energy = + equation_of_state.specific_internal_energy_from_density_and_temperature( + expected_density, expected_temperature, expected_electron_fraction); + + Scalar expected_lorentz_factor{DataVector{ + 1., 7.0710678118654752, 1.0000000001020408, 7.0710678118654752}}; + auto expected_spatial_velocity = + make_with_value>(density, 0.); + expected_spatial_velocity.get(0)[1] = 0.7; + + // The [2] component of the expected velocity is: + // velocity *= max_velocity_magnitude_ * (rho - rho_cut) / + // (trans_rho - rho_cut) / |v^i| + expected_spatial_velocity.get(0)[2] = 0.6 * (4.e-12 - 3.e-12) / + (1.e-11 - 3.e-12) * 1.e-4 / + sqrt(0.6 * 0.6 * 2.); + // The [3] component is the same as the [1] component + expected_spatial_velocity.get(0)[3] = expected_spatial_velocity.get(0)[1]; + + CHECK(get(expected_lorentz_factor).size() == get(expected_density).size()); + CHECK_ITERABLE_APPROX(density, expected_density); + CHECK_ITERABLE_APPROX(pressure, expected_pressure); + CHECK_ITERABLE_APPROX(temperature, expected_temperature); + CHECK_ITERABLE_APPROX(specific_internal_energy, + expected_specific_internal_energy); + CHECK_ITERABLE_APPROX(lorentz_factor, expected_lorentz_factor); + CHECK_ITERABLE_APPROX(spatial_velocity, expected_spatial_velocity); + CHECK_ITERABLE_APPROX(electron_fraction, expected_electron_fraction); +} + template void test_variable_fixer() { // Test for representative 1-d equation of state @@ -163,6 +234,16 @@ void test_variable_fixer() { test_serialization(variable_fixer); test_variable_fixer(fixer_from_options, ideal_fluid); + + // Test for not-so representative 3-d equation of state + // Replace with a realisitic 3-d EoS when that's easy + // enough to do + const EquationsOfState::Barotropic3D> + barotropic_3d_eos{EquationsOfState::PolytropicFluid{100.0, 2.0}}; + test_variable_fixer(variable_fixer, barotropic_3d_eos); + test_serialization(variable_fixer); + + test_variable_fixer(fixer_from_options, ideal_fluid); } } // namespace diff --git a/tests/Unit/Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp b/tests/Unit/Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp index 6486975c3530..122d28f78b98 100644 --- a/tests/Unit/Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp +++ b/tests/Unit/Helpers/Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/PrimReconstructor.hpp @@ -203,7 +203,8 @@ void test_prim_reconstructor_impl( using cons_tags = typename ghmhd::System::variables_tag::tags_list; using spacetime_tags = ::grmhd::GhValenciaDivClean::Tags::spacetime_reconstruction_tags; - + const VariableFixing::FixToAtmosphere<3> fix_to_atmosphere{1.0e-16, 1.1e-16, + 1.0e-15, 1.0}; const Mesh<3> subcell_mesh{points_per_dimension, Spectral::Basis::FiniteDifference, Spectral::Quadrature::CellCentered}; @@ -293,7 +294,8 @@ void test_prim_reconstructor_impl( .reconstruct(make_not_null(&vars_on_lower_face), make_not_null(&vars_on_upper_face), make_not_null(&reconstruction_order), volume_prims, - volume_cons_vars, eos, element, ghost_data, subcell_mesh); + volume_cons_vars, eos, element, ghost_data, subcell_mesh, + fix_to_atmosphere); for (size_t d = 0; d < 3; ++d) { CAPTURE(d); for (size_t i = 0; i < gsl::at(reconstruction_order_storage, d).size(); @@ -307,7 +309,8 @@ void test_prim_reconstructor_impl( dynamic_cast(reconstructor) .reconstruct(make_not_null(&vars_on_lower_face), make_not_null(&vars_on_upper_face), volume_prims, - volume_cons_vars, eos, element, ghost_data, subcell_mesh); + volume_cons_vars, eos, element, ghost_data, subcell_mesh, + fix_to_atmosphere); } for (size_t dim = 0; dim < 3; ++dim) { @@ -457,19 +460,19 @@ void test_prim_reconstructor_impl( num_pts_on_mortar}; if (dim != 2) { dynamic_cast(reconstructor) - .reconstruct_fd_neighbor(make_not_null(&upper_side_vars_on_mortar), - volume_prims, volume_spacetime_vars, eos, - element, ghost_data, subcell_mesh, - Direction<3>{dim, Side::Upper}); + .reconstruct_fd_neighbor( + make_not_null(&upper_side_vars_on_mortar), volume_prims, + volume_spacetime_vars, eos, element, ghost_data, subcell_mesh, + fix_to_atmosphere, Direction<3>{dim, Side::Upper}); } Variables lower_side_vars_on_mortar{ num_pts_on_mortar}; dynamic_cast(reconstructor) - .reconstruct_fd_neighbor(make_not_null(&lower_side_vars_on_mortar), - volume_prims, volume_spacetime_vars, eos, - element, ghost_data, subcell_mesh, - Direction<3>{dim, Side::Lower}); + .reconstruct_fd_neighbor( + make_not_null(&lower_side_vars_on_mortar), volume_prims, + volume_spacetime_vars, eos, element, ghost_data, subcell_mesh, + fix_to_atmosphere, Direction<3>{dim, Side::Lower}); tmpl::for_each>( [dim, &expected_lower_face_values, &expected_upper_face_values, diff --git a/tests/Unit/Helpers/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PrimReconstructor.hpp b/tests/Unit/Helpers/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PrimReconstructor.hpp index a1ddcbfbcb5d..7d11266c8d4d 100644 --- a/tests/Unit/Helpers/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PrimReconstructor.hpp +++ b/tests/Unit/Helpers/Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/PrimReconstructor.hpp @@ -32,6 +32,7 @@ #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" +#include "Evolution/VariableFixing/FixToAtmosphere.hpp" #include "Framework/TestHelpers.hpp" #include "NumericalAlgorithms/Spectral/Basis.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" @@ -114,6 +115,8 @@ void test_prim_reconstructor_impl( const Mesh<3> subcell_mesh{points_per_dimension, Spectral::Basis::FiniteDifference, Spectral::Quadrature::CellCentered}; + const VariableFixing::FixToAtmosphere<3> fix_to_atmosphere{1.0e-16, 1.1e-16, + 1.0e-15, 1.0}; auto logical_coords = logical_coordinates(subcell_mesh); // Make the logical coordinates different in each direction for (size_t i = 1; i < 3; ++i) { @@ -234,7 +237,7 @@ void test_prim_reconstructor_impl( .reconstruct(make_not_null(&vars_on_lower_face), make_not_null(&vars_on_upper_face), make_not_null(&reconstruction_order), volume_prims, eos, - element, ghost_data, subcell_mesh); + element, ghost_data, subcell_mesh, fix_to_atmosphere); for (size_t d = 0; d < 3; ++d) { CAPTURE(d); for (size_t i = 0; i < gsl::at(reconstruction_order_storage, d).size(); @@ -248,7 +251,7 @@ void test_prim_reconstructor_impl( dynamic_cast(reconstructor) .reconstruct(make_not_null(&vars_on_lower_face), make_not_null(&vars_on_upper_face), volume_prims, eos, - element, ghost_data, subcell_mesh); + element, ghost_data, subcell_mesh, fix_to_atmosphere); } for (size_t dim = 0; dim < 3; ++dim) { @@ -384,7 +387,8 @@ void test_prim_reconstructor_impl( dynamic_cast(reconstructor) .reconstruct_fd_neighbor(make_not_null(&upper_side_vars_on_mortar), volume_prims, eos, element, ghost_data, - subcell_mesh, Direction<3>{dim, Side::Upper}); + subcell_mesh, fix_to_atmosphere, + Direction<3>{dim, Side::Upper}); Variables lower_side_vars_on_mortar{ num_pts_on_mortar}; @@ -404,7 +408,8 @@ void test_prim_reconstructor_impl( dynamic_cast(reconstructor) .reconstruct_fd_neighbor(make_not_null(&lower_side_vars_on_mortar), volume_prims, eos, element, ghost_data, - subcell_mesh, Direction<3>{dim, Side::Lower}); + subcell_mesh, fix_to_atmosphere, + Direction<3>{dim, Side::Lower}); tmpl::for_each>( [dim, &expected_lower_face_values, &expected_upper_face_values,