Make FixToAtmosphere fix the electron fraction in atmosphere

Also add fix_to_atmosphere calls in reconstruction

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<true, ThermodynamicDim>& eos,
     const Element<dim>& element,
     const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>& ghost_data,
-    const Mesh<dim>& subcell_mesh) const {
+    const Mesh<dim>& subcell_mesh,
+    const VariableFixing::FixToAtmosphere<dim>& 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 <size_t ThermodynamicDim, typename TagsList>
@@ -135,6 +138,7 @@ void MonotonisedCentralPrim::reconstruct_fd_neighbor(
     const Element<dim>& element,
     const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>& ghost_data,
     const Mesh<dim>& subcell_mesh,
+    const VariableFixing::FixToAtmosphere<dim>& fix_to_atmosphere,
     const Direction<dim> 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<dim>& fix_to_atmosphere) const; \
   template void MonotonisedCentralPrim::reconstruct_fd_neighbor(            \
       gsl::not_null<Variables<tags_list_for_reconstruct_fd_neighbor>*>      \
           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<dim>& fix_to_atmosphere,        \
       const Direction<3> direction_to_reconstruct) const;
 
 GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))

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<dim>,
                  evolution::dg::subcell::Tags::GhostDataForReconstruction<dim>,
-                 evolution::dg::subcell::Tags::Mesh<dim>>;
+                 evolution::dg::subcell::Tags::Mesh<dim>,
+                 ::Tags::VariableFixer<VariableFixing::FixToAtmosphere<dim>>>;
 
   template <size_t ThermodynamicDim, typename TagsList>
   void reconstruct(
@@ -114,7 +117,8 @@ class MonotonisedCentralPrim : public Reconstructor {
       const Element<dim>& element,
       const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>&
           ghost_data,
-      const Mesh<dim>& subcell_mesh) const;
+      const Mesh<dim>& subcell_mesh,
+      const VariableFixing::FixToAtmosphere<dim>& fix_to_atmosphere) const;
 
   /// Called by an element doing DG when the neighbor is doing subcell.
   template <size_t ThermodynamicDim, typename TagsList>
@@ -129,6 +133,7 @@ class MonotonisedCentralPrim : public Reconstructor {
       const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>&
           ghost_data,
       const Mesh<dim>& subcell_mesh,
+      const VariableFixing::FixToAtmosphere<dim>& fix_to_atmosphere,
       const Direction<dim> direction_to_reconstruct) const;
 };
 

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<true, ThermodynamicDim>& eos,
     const Element<dim>& element,
     const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>& ghost_data,
-    const Mesh<dim>& subcell_mesh) const {
+    const Mesh<dim>& subcell_mesh,
+    const VariableFixing::FixToAtmosphere<dim>& 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<tmpl::list<gr::Tags::SpacetimeMetric<DataVector, 3>>,
                          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<dim>& element,
     const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>& ghost_data,
     const Mesh<dim>& subcell_mesh,
-    const Direction<dim> direction_to_reconstruct) const {
+    const VariableFixing::FixToAtmosphere<dim>& fix_to_atmosphere,
+    const Direction<dim>& 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<dim>& fix_to_atmosphere) const; \
   template void                                                             \
   PositivityPreservingAdaptiveOrderPrim::reconstruct_fd_neighbor(           \
       gsl::not_null<Variables<tags_list_for_reconstruct_fd_neighbor>*>      \
@@ -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<dim>& fix_to_atmosphere,        \
+      const Direction<3>& direction_to_reconstruct) const;
 
 GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))
 

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<dim>,
                  evolution::dg::subcell::Tags::GhostDataForReconstruction<dim>,
-                 evolution::dg::subcell::Tags::Mesh<dim>>;
+                 evolution::dg::subcell::Tags::Mesh<dim>,
+                 ::Tags::VariableFixer<VariableFixing::FixToAtmosphere<dim>>>;
 
   template <size_t ThermodynamicDim, typename TagsList>
   void reconstruct(
@@ -174,7 +177,8 @@ class PositivityPreservingAdaptiveOrderPrim : public Reconstructor {
       const Element<dim>& element,
       const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>&
           ghost_data,
-      const Mesh<dim>& subcell_mesh) const;
+      const Mesh<dim>& subcell_mesh,
+      const VariableFixing::FixToAtmosphere<dim>& fix_to_atmosphere) const;
 
   /// Called by an element doing DG when the neighbor is doing subcell.
   template <size_t ThermodynamicDim, typename TagsList>
@@ -189,7 +193,8 @@ class PositivityPreservingAdaptiveOrderPrim : public Reconstructor {
       const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>&
           ghost_data,
       const Mesh<dim>& subcell_mesh,
-      const Direction<dim> direction_to_reconstruct) const;
+      const VariableFixing::FixToAtmosphere<dim>& fix_to_atmosphere,
+      const Direction<dim>& direction_to_reconstruct) const;
 
  private:
   // NOLINTNEXTLINE(readability-redundant-declaration)

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<PrimsTagsSentByNeighbor>>&
         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

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<PrimsTagsSentByNeighbor>>&
         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

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<true, ThermodynamicDim>& eos,
     const Element<dim>& element,
     const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>& ghost_data,
-    const Mesh<dim>& subcell_mesh) const {
+    const Mesh<dim>& subcell_mesh,
+    const VariableFixing::FixToAtmosphere<dim>& 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 <size_t ThermodynamicDim, typename TagsList>
@@ -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<dim>& 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<dim>& fix_to_atmosphere) const; \
   template void Wcns5zPrim::reconstruct_fd_neighbor(                        \
       gsl::not_null<Variables<tags_list_for_reconstruct_fd_neighbor>*>      \
           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<dim>& fix_to_atmosphere,        \
       const Direction<3>& direction_to_reconstruct) const;
 
 GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))