Better domain errors, atmosphere density in ID, and velocity limiting in atmosphere.

src/Domain/DomainHelpers.cpp

@@ -349,8 +349,9 @@
     const gsl::not_null<
         std::vector<DirectionMap<VolumeDim, BlockNeighbor<VolumeDim>>>*>
         neighbors_of_all_blocks) {
-  ASSERT(orientations_of_all_blocks.size() == corners_of_all_blocks.size(),
-         "Each block must have an OrientationMap relative to an edifice.");
+  if (orientations_of_all_blocks.size() != corners_of_all_blocks.size()) {
+    ERROR("Each block must have an OrientationMap relative to an edifice.");
+  }
   size_t block_id = 0;
   for (auto& block : *neighbors_of_all_blocks) {
     const auto& corners_of_block = corners_of_all_blocks[block_id];
@@ -602,13 +603,17 @@
     const std::vector<domain::CoordinateMaps::Distribution>&
         radial_distribution,
     const ShellWedges which_wedges, const double opening_angle) {
-  ASSERT(not use_half_wedges or which_wedges == ShellWedges::All,
-         "If we are using half wedges we must also be using ShellWedges::All.");
-  ASSERT(radial_distribution.size() == 1 + radial_partitioning.size(),
-         "Specify a radial distribution for every spherical shell. You "
-         "specified "
-             << radial_distribution.size() << " items, but the domain has "
-             << 1 + radial_partitioning.size() << " shells.");
+  if (use_half_wedges and which_wedges != ShellWedges::All) {
+    ERROR(
+        "If we are using half wedges we must also be using ShellWedges::All.");
+  }
+  if (radial_distribution.size() != 1 + radial_partitioning.size()) {
+    ERROR(
+        "Specify a radial distribution for every spherical shell. You "
+        "specified "
+        << radial_distribution.size() << " items, but the domain has "
+        << 1 + radial_partitioning.size() << " shells.");
+  }
 
   const auto wedge_orientations = orientations_for_sphere_wrappings();
 
@@ -738,17 +743,22 @@
     const domain::CoordinateMaps::Distribution radial_distribution,
     const std::optional<double> distribution_value, const double sphericity,
     const double opening_angle) {
-  ASSERT(length_inner_cube < 0.5 * length_outer_cube,
-         "The outer cube is too small! The inner cubes will pierce the surface "
-         "of the outer cube.");
-  ASSERT(
-      abs(origin_preimage[0]) + length_inner_cube < 0.5 * length_outer_cube and
+  if (length_inner_cube >= 0.5 * length_outer_cube) {
+    ERROR("The outer cube (" << length_outer_cube
+                             << ") is too small! The inner cubes ("
+                             << length_inner_cube
+                             << ") will pierce the surface of the outer cube.");
+  }
+  if (not(abs(origin_preimage[0]) + length_inner_cube <
+              0.5 * length_outer_cube and
           abs(origin_preimage[1]) + length_inner_cube <
               0.5 * length_outer_cube and
           abs(origin_preimage[2]) + 0.5 * length_inner_cube <
-              0.5 * length_outer_cube,
-      "The current choice for `origin_preimage` results in the inner cubes "
-      "piercing the surface of the outer cube.");
+              0.5 * length_outer_cube)) {
+    ERROR(
+        "The current choice for `origin_preimage` results in the inner cubes "
+        "piercing the surface of the outer cube.");
+  }
   const auto frustum_orientations = orientations_for_sphere_wrappings();
   const double lower = 0.5 * length_inner_cube;
   const double top = 0.5 * length_outer_cube;

src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp

@@ -497,9 +497,8 @@ struct GhValenciaDivCleanTemplateBase<
                   volume_dim, Frame::ElementLogical, Frame::Inertial>,
               ::Events::Tags::ObserverMesh<volume_dim>>,
           gr::Tags::WeylElectricCompute<DataVector, 3, Frame::Inertial>,
-          gr::Tags::Psi4RealCompute<Frame::Inertial>
-
-          >,
+          gr::Tags::Psi4RealCompute<Frame::Inertial>,
+          ::Events::Tags::ObserverMeshVelocity<3>>,
       tmpl::conditional_t<
           use_dg_subcell,
           tmpl::list<evolution::dg::subcell::Tags::TciStatusCompute<volume_dim>,
@@ -569,7 +568,8 @@ struct GhValenciaDivCleanTemplateBase<
                      ::Events::Tags::ObserverJacobianCompute<
                          volume_dim, Frame::ElementLogical, Frame::Inertial>,
                      ::Events::Tags::ObserverDetInvJacobianCompute<
-                         Frame::ElementLogical, Frame::Inertial>>>,
+                         Frame::ElementLogical, Frame::Inertial>,
+                     ::Events::Tags::ObserverMeshVelocityCompute<3>>>,
       tmpl::list<analytic_compute, error_compute>,
       tmpl::list<::Tags::DerivCompute<
                      typename system::variables_tag,

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.tpp

@@ -39,9 +39,8 @@ class FunctionOfZ {
               const double electron_fraction, const EosType& equation_of_state,
               const double lorentz_max)
       : q_(tau / rest_mass_density_times_lorentz_factor),
-        r_squared_(momentum_density_squared /
-                   square(rest_mass_density_times_lorentz_factor)),
-        r_(std::sqrt(r_squared_)),
+        r_(std::sqrt(momentum_density_squared /
+                     square(rest_mass_density_times_lorentz_factor))),
         rest_mass_density_times_lorentz_factor_(
             rest_mass_density_times_lorentz_factor),
         electron_fraction_(electron_fraction),
@@ -75,17 +74,16 @@ class FunctionOfZ {
     }
   }
 
-  std::pair<double, double> root_bracket();
+  std::pair<double, double> root_bracket() const;
 
   Primitives primitives(double z) const;
 
   double operator()(double z) const;
 
-  bool has_no_root() { return state_is_unphysical_; };
+  bool has_no_root() const { return state_is_unphysical_; }
 
  private:
   double q_;
-  double r_squared_;
   double r_;
   bool state_is_unphysical_ = false;
   const double rest_mass_density_times_lorentz_factor_;
@@ -95,8 +93,8 @@ class FunctionOfZ {
 
 template <typename EosType, bool EnforcePhysicality>
 std::pair<double, double>
-FunctionOfZ<EosType, EnforcePhysicality>::root_bracket() {
-  auto k = r_ / (q_ + 1.);
+FunctionOfZ<EosType, EnforcePhysicality>::root_bracket() const {
+  const auto k = r_ / (q_ + 1.);
 
   const double rho_min = equation_of_state_.rest_mass_density_lower_bound();
 
@@ -108,17 +106,17 @@ FunctionOfZ<EosType, EnforcePhysicality>::root_bracket() {
   }
 
   // Compute bounds (C23)
-  double lower_bound = 0.5 * k / std::sqrt(std::abs(1. - 0.25 * k * k));
+  const double lower_bound = 0.5 * k / std::sqrt(std::abs(1. - 0.25 * k * k));
   // Ensure that upper_bound does not become degenerate when k ~ 0
   // Empirically, an offset of 1.e-8 has worked well.
-  double upper_bound = 1.e-8 + k / std::sqrt(std::abs(1. - k * k));
+  const double upper_bound = 1.e-8 + k / std::sqrt(std::abs(1. - k * k));
 
   return {lower_bound, upper_bound};
 }
 
 template <typename EosType, bool EnforcePhysicality>
 Primitives FunctionOfZ<EosType, EnforcePhysicality>::primitives(
-    double z) const {
+    const double z) const {
   // Compute Lorentz factor, note that z = lorentz * v
   const double w_hat = std::sqrt(1. + z * z);
 
@@ -158,20 +156,21 @@ Primitives FunctionOfZ<EosType, EnforcePhysicality>::primitives(
     // 1d-EOS. Instead, we recover the primitives and then reset the
     // specific internal energy and specific enthalpy using the EOS.
     p_hat =
-        get(equation_of_state_.pressure_from_density(Scalar<double>(rho_hat)));
+        get(equation_of_state_.pressure_from_density(Scalar<double>{rho_hat}));
   } else if constexpr (EosType::thermodynamic_dim == 2) {
     p_hat = get(equation_of_state_.pressure_from_density_and_energy(
-        Scalar<double>(rho_hat), Scalar<double>(epsilon_hat)));
+        Scalar<double>{rho_hat}, Scalar<double>{epsilon_hat}));
   } else if constexpr (EosType::thermodynamic_dim == 3) {
     p_hat = get(equation_of_state_.pressure_from_density_and_energy(
-        Scalar<double>(rho_hat), Scalar<double>(epsilon_hat),
-        Scalar<double>(electron_fraction_)));
+        Scalar<double>{rho_hat}, Scalar<double>{epsilon_hat},
+        Scalar<double>{electron_fraction_}));
   }
   return Primitives{rho_hat, w_hat, p_hat, epsilon_hat};
 }
 
 template <typename EosType, bool EnforcePhysicality>
-double FunctionOfZ<EosType, EnforcePhysicality>::operator()(double z) const {
+double FunctionOfZ<EosType, EnforcePhysicality>::operator()(
+    const double z) const {
   const auto [rho_hat, w_hat, p_hat, epsilon_hat] = primitives(z);
   // Equation (C5)
   const double a_hat = p_hat / (rho_hat * (1.0 + epsilon_hat));
@@ -193,7 +192,7 @@ std::optional<PrimitiveRecoveryData> KastaunEtAlHydro::apply(
     const grmhd::ValenciaDivClean::PrimitiveFromConservativeOptions&
         primitive_from_conservative_options) {
   // Master function see Equation (44)
-  auto f_of_z =
+  const auto f_of_z =
       KastaunEtAlHydro_detail::FunctionOfZ<EosType, EnforcePhysicality>{
           tau,
           momentum_density_squared,
@@ -233,7 +232,7 @@ std::optional<PrimitiveRecoveryData> KastaunEtAlHydro::apply(
       pressure,
       specific_internal_energy,
       (rest_mass_density * (1. + specific_internal_energy) + pressure) *
-          (1. + z * z),
+          square(lorentz_factor),
       electron_fraction};
 }
 }  // namespace grmhd::ValenciaDivClean::PrimitiveRecoverySchemes