Merge pull request #6455 from nilsdeppe/bns_improvements_11

Allow skipping GRMHD RHS computation in atmosphere/wavezone

src/Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp

@@ -4,6 +4,7 @@
 #pragma once
 
 #include <cstddef>
+#include <cstdlib>
 #include <optional>
 #include <ostream>
 #include <type_traits>
@@ -14,6 +15,7 @@
 #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
+#include "Evolution/DiscontinuousGalerkin/TimeDerivativeDecisions.hpp"
 #include "Evolution/PassVariables.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/MetricIdentityJacobian.hpp"
@@ -112,8 +114,7 @@ void volume_terms(
       "are defined, and that at least one of them is a non-empty list of "
       "tags.");
 
-  using flux_variables =
-      tmpl::list<FluxVariablesTags...>;
+  using flux_variables = tmpl::list<FluxVariablesTags...>;
 
   // Compute d_i u_\alpha for nonconservative products
   if constexpr (has_partial_derivs) {
@@ -124,25 +125,26 @@ void volume_terms(
   // For now just zero dt_vars. If this is a performance bottle neck we
   // can re-evaluate in the future.
   dt_vars_ptr->initialize(mesh.number_of_grid_points(), 0.0);
+  evolution::dg::TimeDerivativeDecisions<Dim> time_derivative_decisions{};
 
   // Compute volume du/dt and fluxes
   if constexpr (std::is_base_of_v<evolution::PassVariables,
                                   ComputeVolumeTimeDerivativeTerms>) {
     if constexpr (sizeof...(FluxVariablesTags) != 0) {
-      ComputeVolumeTimeDerivativeTerms::apply(
+      time_derivative_decisions = ComputeVolumeTimeDerivativeTerms::apply(
           dt_vars_ptr, volume_fluxes, temporaries,
           get<::Tags::deriv<PartialDerivTags, tmpl::size_t<Dim>,
                             Frame::Inertial>>(*partial_derivs)...,
           time_derivative_args...);
     } else {
-      ComputeVolumeTimeDerivativeTerms::apply(
+      time_derivative_decisions = ComputeVolumeTimeDerivativeTerms::apply(
           dt_vars_ptr, temporaries,
           get<::Tags::deriv<PartialDerivTags, tmpl::size_t<Dim>,
                             Frame::Inertial>>(*partial_derivs)...,
           time_derivative_args...);
     }
   } else {
-    ComputeVolumeTimeDerivativeTerms::apply(
+    time_derivative_decisions = ComputeVolumeTimeDerivativeTerms::apply(
         make_not_null(&get<::Tags::dt<VariablesTags>>(*dt_vars_ptr))...,
         make_not_null(&get<::Tags::Flux<FluxVariablesTags, tmpl::size_t<Dim>,
                                         Frame::Inertial>>(*volume_fluxes))...,
@@ -154,6 +156,9 @@ void volume_terms(
 
   // Add volume terms for moving meshes
   if (mesh_velocity.has_value()) {
+    if (not time_derivative_decisions.compute_flux_divergence) {
+      goto end_of_flux_mesh_velocity;  // NOLINT(cppcoreguidelines-avoid-goto)
+    }
     tmpl::for_each<flux_variables>([&div_mesh_velocity, &dt_vars_ptr,
                                     &evolved_vars, &mesh_velocity,
                                     &volume_fluxes](auto tag_v) {
@@ -193,6 +198,7 @@ void volume_terms(
             get(*div_mesh_velocity);
       }
     });
+  end_of_flux_mesh_velocity:
 
     // We add the mesh velocity to all equations that don't have flux terms.
     // This doesn't need to be equal to the equations that have partial
@@ -237,6 +243,9 @@ void volume_terms(
   // Add the flux divergence term to du_\alpha/dt, which must be done
   // after the corrections for the moving mesh are made.
   if constexpr (has_fluxes) {
+    if (not time_derivative_decisions.compute_flux_divergence) {
+      goto end_of_flux_div_calculation;  // NOLINT(cppcoreguidelines-avoid-goto)
+    }
     if (dg_formulation == ::dg::Formulation::StrongInertial) {
       divergence(div_fluxes, *volume_fluxes, mesh,
                  logical_to_inertial_inverse_jacobian);
@@ -289,6 +298,8 @@ void volume_terms(
             }
           }
         });
+  end_of_flux_div_calculation:
+    (void)5;  // null statement required after label at end of block until C++23
   } else {
     (void)dg_formulation;
   }

src/Evolution/Systems/Burgers/TimeDerivativeTerms.cpp

@@ -8,10 +8,11 @@
 #include "Utilities/Gsl.hpp"
 
 namespace Burgers {
-void TimeDerivativeTerms::apply(
+evolution::dg::TimeDerivativeDecisions<1> TimeDerivativeTerms::apply(
     const gsl::not_null<Scalar<DataVector>*> /*non_flux_terms_dt_vars*/,
     const gsl::not_null<tnsr::I<DataVector, 1, Frame::Inertial>*> flux_u,
     const Scalar<DataVector>& u) {
   get<0>(*flux_u) = 0.5 * square(get(u));
+  return {true};
 }
 }  // namespace Burgers

src/Evolution/Systems/Burgers/TimeDerivativeTerms.hpp

@@ -5,6 +5,7 @@
 
 #include "DataStructures/DataBox/Prefixes.hpp"
 #include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Evolution/DiscontinuousGalerkin/TimeDerivativeDecisions.hpp"
 #include "Utilities/TMPL.hpp"
 
 /// \cond
@@ -30,7 +31,7 @@ struct TimeDerivativeTerms {
   using temporary_tags = tmpl::list<>;
   using argument_tags = tmpl::list<Tags::U>;
 
-  static void apply(
+  static evolution::dg::TimeDerivativeDecisions<1> apply(
       // Time derivatives returned by reference. No source terms or
       // nonconservative products, so not used. All the tags in the
       // variables_tag in the system struct.

src/Evolution/Systems/CurvedScalarWave/TimeDerivative.cpp

@@ -16,7 +16,7 @@
 
 namespace CurvedScalarWave {
 template <size_t Dim>
-void TimeDerivative<Dim>::apply(
+evolution::dg::TimeDerivativeDecisions<Dim> TimeDerivative<Dim>::apply(
     const gsl::not_null<Scalar<DataVector>*> dt_psi,
     const gsl::not_null<Scalar<DataVector>*> dt_pi,
     const gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> dt_phi,
@@ -62,6 +62,7 @@ void TimeDerivative<Dim>::apply(
                   gamma2() * lapse() * (d_psi(ti::i) - phi(ti::i)) -
                   pi() * deriv_lapse(ti::i) +
                   phi(ti::j) * deriv_shift(ti::i, ti::J));
+  return {true};
 }
 }  // namespace CurvedScalarWave
 // Generate explicit instantiations of partial_derivatives function as well as

src/Evolution/Systems/CurvedScalarWave/TimeDerivative.hpp

@@ -6,6 +6,7 @@
 #include <cstddef>
 
 #include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Evolution/DiscontinuousGalerkin/TimeDerivativeDecisions.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Characteristics.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
@@ -70,7 +71,7 @@ struct TimeDerivative {
                  gr::Tags::TraceExtrinsicCurvature<DataVector>,
                  Tags::ConstraintGamma1, Tags::ConstraintGamma2>;
 
-  static void apply(
+  static evolution::dg::TimeDerivativeDecisions<Dim> apply(
       gsl::not_null<Scalar<DataVector>*> dt_psi,
       gsl::not_null<Scalar<DataVector>*> dt_pi,
       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> dt_phi,

src/Evolution/Systems/ForceFree/TimeDerivativeTerms.cpp

@@ -16,7 +16,7 @@
 
 namespace ForceFree {
 
-void TimeDerivativeTerms::apply(
+evolution::dg::TimeDerivativeDecisions<3> TimeDerivativeTerms::apply(
     const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
         non_flux_terms_dt_tilde_e,
     const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
@@ -107,6 +107,7 @@ void TimeDerivativeTerms::apply(
                        tilde_psi, tilde_phi, tilde_q, *tilde_j_drift, kappa_psi,
                        kappa_phi, lapse, d_lapse, d_shift, inv_spatial_metric,
                        extrinsic_curvature);
+  return {true};
 }
 
 }  // namespace ForceFree

src/Evolution/Systems/ForceFree/TimeDerivativeTerms.hpp

@@ -5,6 +5,7 @@
 
 #include "DataStructures/DataBox/DataBoxTag.hpp"
 #include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Evolution/DiscontinuousGalerkin/TimeDerivativeDecisions.hpp"
 #include "Evolution/Systems/ForceFree/Tags.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
 #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
@@ -67,7 +68,7 @@ struct TimeDerivativeTerms {
       ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
                     Frame::Inertial>>;
 
-  static void apply(
+  static evolution::dg::TimeDerivativeDecisions<3> apply(
       // Time derivatives returned by reference.
       gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
           non_flux_terms_dt_tilde_e,

src/Evolution/Systems/GeneralizedHarmonic/TimeDerivative.hpp

@@ -8,6 +8,7 @@
 #include "DataStructures/Tensor/TypeAliases.hpp"
 #include "Domain/Tags.hpp"
 #include "Domain/TagsTimeDependent.hpp"
+#include "Evolution/DiscontinuousGalerkin/TimeDerivativeDecisions.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/ConstraintDamping/Tags.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/DuDtTempTags.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Gauges.hpp"
@@ -140,7 +141,7 @@ struct TimeDerivative {
                                                Frame::Inertial>,
                  domain::Tags::MeshVelocity<Dim, Frame::Inertial>>;
 
-  static void apply(
+  static evolution::dg::TimeDerivativeDecisions<Dim> apply(
       gsl::not_null<tnsr::aa<DataVector, Dim>*> dt_spacetime_metric,
       gsl::not_null<tnsr::aa<DataVector, Dim>*> dt_pi,
       gsl::not_null<tnsr::iaa<DataVector, Dim>*> dt_phi,

src/Evolution/Systems/GeneralizedHarmonic/TimeDerivative.tpp

@@ -11,6 +11,7 @@
 #include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"
 #include "DataStructures/Tensor/EagerMath/Trace.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
+#include "Evolution/DiscontinuousGalerkin/TimeDerivativeDecisions.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/ConstraintDamping/Tags.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/DuDtTempTags.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Dispatch.hpp"
@@ -30,7 +31,8 @@
 
 namespace gh {
 template <class AllSolutionsForChristoffelAnalytic, size_t Dim>
-void TimeDerivative<AllSolutionsForChristoffelAnalytic, Dim>::apply(
+evolution::dg::TimeDerivativeDecisions<Dim>
+TimeDerivative<AllSolutionsForChristoffelAnalytic, Dim>::apply(
     const gsl::not_null<tnsr::aa<DataVector, Dim>*> dt_spacetime_metric,
     const gsl::not_null<tnsr::aa<DataVector, Dim>*> dt_pi,
     const gsl::not_null<tnsr::iaa<DataVector, Dim>*> dt_phi,
@@ -406,5 +408,6 @@ void TimeDerivative<AllSolutionsForChristoffelAnalytic, Dim>::apply(
       }
     }
   }
+  return {true};
 }
 }  // namespace gh

src/Evolution/Systems/GrMhd/GhValenciaDivClean/TimeDerivativeTerms.hpp

@@ -4,13 +4,15 @@
 #pragma once
 
 #include <cstddef>
+#include <limits>
 #include <utility>
 
 #include "DataStructures/DataBox/PrefixHelpers.hpp"
 #include "DataStructures/DataBox/Prefixes.hpp"
 #include "DataStructures/DataVector.hpp"
 #include "DataStructures/TaggedContainers.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
+#include "Evolution/DiscontinuousGalerkin/TimeDerivativeDecisions.hpp"
 #include "Evolution/PassVariables.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Harmonic.hpp"
 #include "Evolution/Systems/GeneralizedHarmonic/System.hpp"
@@ -21,6 +23,8 @@
 #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp"
 #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
 #include "Evolution/Systems/GrMhd/ValenciaDivClean/TimeDerivativeTerms.hpp"
+#include "Evolution/VariableFixing/FixToAtmosphere.hpp"
+#include "Evolution/VariableFixing/Tags.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
 // Tag obtained from gh::TimeDerivative needs to be complete here to
 // be used in TemporaryReference.
@@ -68,7 +72,7 @@ struct TimeDerivativeTermsImpl<
     tmpl::list<TraceReversedStressResultTags...>,
     tmpl::list<TraceReversedStressArgumentTags...>> {
   template <typename TemporaryTagsList, typename... ExtraTags>
-  static void apply(
+  static evolution::dg::TimeDerivativeDecisions<3> apply(
       const gsl::not_null<
           Variables<tmpl::list<GhDtTags..., ValenciaDtTags...>>*>
           dt_vars_ptr,
@@ -92,6 +96,39 @@ struct TimeDerivativeTermsImpl<
                      get<Tags::detail::TemporaryReference<GhArgTags>>(
                          arguments)...);
 
+    // If we are in the atmosphere, then we can skip the evolution
+    // of the GRMHD system completely. This inevitably depends on parameters
+    // of the FixToAtmosphere code, so we grab those directly.
+    if (const auto& fix_to_atmosphere = get<Tags::detail::TemporaryReference<
+            ::Tags::VariableFixer<::VariableFixing::FixToAtmosphere<3>>>>(
+            arguments);
+        max(get(get<Tags::detail::TemporaryReference<
+                    hydro::Tags::RestMassDensity<DataVector>>>(arguments))) <=
+        std::min({fix_to_atmosphere.density_of_atmosphere(),
+                  (fix_to_atmosphere.velocity_limiting().has_value()
+                       ? fix_to_atmosphere.velocity_limiting()
+                             ->atmosphere_density_cutoff
+                       : std::numeric_limits<double>::infinity()),
+                  (fix_to_atmosphere.kappa_limiting().has_value()
+                       ? fix_to_atmosphere.kappa_limiting()->density_lower_bound
+                       : std::numeric_limits<double>::infinity())}) *
+            (1.0 + 10.0 * std::numeric_limits<double>::epsilon())) {
+      // Point into the right memory, then set it to zero.
+      ASSERT(
+          max(get(get<tmpl::front<tmpl::list<ValenciaDtTags...>>>(
+              *dt_vars_ptr))) == 0.0,
+          "GH+GRMHD assumes the time derivatives are set to zero in general."
+          " If this is no longer the case, please set them to zero in "
+          "atmosphere by changing the code where this ASSERT was triggered.");
+      // Code that we could use to set the sources to zero if needed.
+      // Variables<tmpl::list<ValenciaDtTags...>> dt_div_clean(
+      //     get<tmpl::front<tmpl::list<ValenciaDtTags...>>>(*dt_vars_ptr)[0]
+      //         .data(),
+      //     0.0);
+      fluxes_ptr->initialize(fluxes_ptr->number_of_grid_points(), 0.0);
+      return evolution::dg::TimeDerivativeDecisions<3>{false};
+    }
+
     if (get<Tags::detail::TemporaryReference<gh::gauges::Tags::GaugeCondition>>(
             arguments)
             .is_harmonic()) {
@@ -174,6 +211,7 @@ struct TimeDerivativeTermsImpl<
         get<grmhd::GhValenciaDivClean::Tags::TraceReversedStressEnergy>(
             *temps_ptr),
         get<gr::Tags::Lapse<DataVector>>(*temps_ptr));
+    return evolution::dg::TimeDerivativeDecisions<3>{true};
   }
 };
 }  // namespace detail
@@ -263,15 +301,18 @@ struct TimeDerivativeTerms : evolution::PassVariables {
           gh_temp_tags, valencia_temp_tags, valencia_extra_temp_tags,
           trace_reversed_stress_result_tags, extra_temp_tags>>,
       gr::Tags::SpatialMetric<DataVector, 3>>;
-  using argument_tags =
-      tmpl::remove<tmpl::remove<tmpl::append<gh_arg_tags,
+  using argument_tags = tmpl::remove<
+      tmpl::remove<tmpl::append<gh_arg_tags,
+
+                                valencia_arg_tags,
 
-                                             valencia_arg_tags>,
-                                gr::Tags::SpatialMetric<DataVector, 3>>,
-                   d_spatial_metric>;
+                                tmpl::list<::Tags::VariableFixer<
+                                    ::VariableFixing::FixToAtmosphere<3>>>>,
+                   gr::Tags::SpatialMetric<DataVector, 3>>,
+      d_spatial_metric>;
 
   template <typename... Args>
-  static void apply(
+  static evolution::dg::TimeDerivativeDecisions<3> apply(
       const gsl::not_null<Variables<dt_tags>*> dt_vars_ptr,
       const gsl::not_null<Variables<db::wrap_tags_in<
           ::Tags::Flux, typename ValenciaDivClean::System::flux_variables,
@@ -317,7 +358,7 @@ struct TimeDerivativeTerms : evolution::PassVariables {
       }
     }
 
-    detail::TimeDerivativeTermsImpl<
+    return detail::TimeDerivativeTermsImpl<
         gh_dt_tags, valencia_dt_tags, valencia_flux_tags, gh_temp_tags,
         valencia_temp_tags, gh_gradient_tags, gh_arg_tags, valencia_arg_tags,
         typename grmhd::ValenciaDivClean::TimeDerivativeTerms::argument_tags,

src/Evolution/Systems/GrMhd/GhValenciaDivClean/VolumeTermsInstantiation.cpp

@@ -8,6 +8,7 @@
 #include "Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp"
 #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp"
 #include "Evolution/Systems/GrMhd/GhValenciaDivClean/TimeDerivativeTerms.hpp"
+#include "Evolution/VariableFixing/FixToAtmosphere.hpp"
 
 namespace evolution::dg::Actions::detail {
 template void volume_terms<::grmhd::GhValenciaDivClean::TimeDerivativeTerms>(
@@ -72,5 +73,6 @@ template void volume_terms<::grmhd::GhValenciaDivClean::TimeDerivativeTerms>(
     const Scalar<DataVector>& rest_mass_density,
     const Scalar<DataVector>& electron_fraction,
     const Scalar<DataVector>& specific_internal_energy,
-    const double& constraint_damping_parameter);
+    const double& constraint_damping_parameter,
+    const ::VariableFixing::FixToAtmosphere<3>& fix_to_atmosphere);
 }  // namespace evolution::dg::Actions::detail