Add source terms for ValenciaDivClean with Cartoon

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: michaeldmurphy1

src/Evolution/Systems/GrMhd/ValenciaDivClean/Sources.cpp

@@ -13,9 +13,13 @@
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
 #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "NumericalAlgorithms/Spectral/Quadrature.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Christoffel.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "PointwiseFunctions/Hydro/ComovingMagneticField.hpp"
 #include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/Hydro/TransportVelocity.hpp"
 #include "Utilities/ConstantExpressions.hpp"
 #include "Utilities/Gsl.hpp"
 
@@ -65,6 +69,9 @@ struct PressureStar : db::SimpleTag {
 struct EnthalpyTimesDensityWSquaredPlusBSquared : db::SimpleTag {
   using type = Scalar<DataVector>;
 };
+struct ComovingMagneticFieldOneForm : db::SimpleTag {
+  using type = tnsr::a<DataVector, 3, Frame::Inertial>;
+};
 }  // namespace
 
 namespace grmhd::ValenciaDivClean {
@@ -165,6 +172,97 @@ void sources_impl(
     get(*source_tilde_phi) += tilde_b.get(m) * d_lapse.get(m);
   }
 }
+
+void cartoon_sources_impl(
+    const gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
+        source_tilde_s,
+    const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+        source_tilde_b,
+
+    const Scalar<DataVector>& pressure_star,
+    const tnsr::i<DataVector, 3, Frame::Inertial>& magnetic_field_one_form,
+    const Scalar<DataVector>& magnetic_field_dot_spatial_velocity,
+
+    const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
+    const Scalar<DataVector>& tilde_phi,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
+    const Scalar<DataVector>& lorentz_factor, const Scalar<DataVector>& lapse,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
+    const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
+    const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
+    const Scalar<DataVector>& sqrt_det_spatial_metric,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
+    const Spectral::Quadrature cartoon_quadrature) {
+#ifdef SPECTRE_DEBUG
+  for (size_t i = 0; i < get<0>(inertial_coords).size(); ++i) {
+    ASSERT(not equal_within_roundoff(
+               0.0, get<0>(inertial_coords)[i],
+               std::numeric_limits<double>::epsilon() * 100.0,
+               max(get<0>(inertial_coords))),
+           "Computing cartoon source terms with x=0 in domain, got "
+           "inertial_coords = "
+               << inertial_coords);
+  }
+#endif  // SPECTRE_DEBUG
+  if (cartoon_quadrature == Spectral::Quadrature::SphericalSymmetry) {
+    get<0>(*source_tilde_s) += 2.0 * get(lapse) * get(sqrt_det_spatial_metric) *
+                               get(pressure_star) / get<0>(inertial_coords);
+    get<0>(*source_tilde_b) +=
+        get(lapse) *
+        (get<1, 1>(inv_spatial_metric) + get<2, 2>(inv_spatial_metric)) *
+        get(tilde_phi) / get<0>(inertial_coords);
+  } else {
+    ASSERT(
+        cartoon_quadrature == Spectral::Quadrature::AxialSymmetry,
+        "Got unexpected quadrature for Cartoon basis: " << cartoon_quadrature);
+
+    Variables<tmpl::list<hydro::Tags::SpatialVelocityOneForm<DataVector, 3>,
+                         ComovingMagneticFieldOneForm,
+                         hydro::Tags::TransportVelocity<DataVector, 3>>>
+        temp_tensors(get<0>(inertial_coords).size());
+
+    auto& spatial_velocity_one_form =
+        get<hydro::Tags::SpatialVelocityOneForm<DataVector, 3>>(temp_tensors);
+    raise_or_lower_index(make_not_null(&spatial_velocity_one_form),
+                         spatial_velocity, spatial_metric);
+
+    auto& comoving_magnetic_field_one_form =
+        get<ComovingMagneticFieldOneForm>(temp_tensors);
+    hydro::comoving_magnetic_field_one_form(
+        make_not_null(&comoving_magnetic_field_one_form),
+        spatial_velocity_one_form, magnetic_field_one_form,
+        magnetic_field_dot_spatial_velocity, lorentz_factor, shift, lapse);
+
+    auto& transport_velocity =
+        get<hydro::Tags::TransportVelocity<DataVector, 3>>(temp_tensors);
+    hydro::transport_velocity(make_not_null(&transport_velocity),
+                              spatial_velocity, lapse, shift);
+
+    get<0>(*source_tilde_s) +=
+        (get(lapse) * get(sqrt_det_spatial_metric) * get(pressure_star) +
+         get(lapse) * get<2>(tilde_b) *
+             (-get<3>(comoving_magnetic_field_one_form)) / get(lorentz_factor) +
+         (get<2>(tilde_s)) * get<2>(transport_velocity)) /
+        get<0>(inertial_coords);
+    get<2>(*source_tilde_s) +=
+        (get(lapse) * get<2>(tilde_b) *
+             (get<1>(comoving_magnetic_field_one_form)) / get(lorentz_factor) +
+         (-get<0>(tilde_s)) * get<2>(transport_velocity)) /
+        get<0>(inertial_coords);
+
+    get<0>(*source_tilde_b) +=
+        (get(lapse) * (get<2, 2>(inv_spatial_metric)) * get(tilde_phi) +
+         get(lapse) * (-get<2>(spatial_velocity)) * get<2>(tilde_b) +
+         (get<2>(tilde_b)) * get<2>(transport_velocity)) /
+        get<0>(inertial_coords);
+    get<2>(*source_tilde_b) +=
+        (get(lapse) * (-get<2, 0>(inv_spatial_metric)) * get(tilde_phi) +
+         get(lapse) * (get<0>(spatial_velocity)) * get<2>(tilde_b) +
+         (-get<0>(tilde_b)) * get<2>(transport_velocity)) /
+        get<0>(inertial_coords);
+  }
+}
 }  // namespace detail
 
 void ComputeSources::apply(
@@ -186,14 +284,17 @@ void ComputeSources::apply(
     const Scalar<DataVector>& specific_internal_energy,
     const Scalar<DataVector>& lorentz_factor,
     const Scalar<DataVector>& pressure, const Scalar<DataVector>& lapse,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
     const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
     const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
     const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
     const tnsr::ijj<DataVector, 3, Frame::Inertial>& d_spatial_metric,
     const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
     const Scalar<DataVector>& sqrt_det_spatial_metric,
     const tnsr::ii<DataVector, 3, Frame::Inertial>& extrinsic_curvature,
-    const double constraint_damping_parameter) {
+    const double constraint_damping_parameter,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
+    const Mesh<3>& dg_mesh) {
   Variables<
       tmpl::list<TildeSUp, DensitizedStress, MagneticFieldOneForm,
                  hydro::Tags::MagneticFieldDotSpatialVelocity<DataVector>,
@@ -260,5 +361,14 @@ void ComputeSources::apply(
 
       rest_mass_density, electron_fraction, pressure, specific_internal_energy,
       extrinsic_curvature, constraint_damping_parameter);
+
+  if (dg_mesh.basis(2) == Spectral::Basis::Cartoon) {
+    detail::cartoon_sources_impl(
+        source_tilde_s, source_tilde_b, pressure_star, magnetic_field_oneform,
+        magnetic_field_dot_spatial_velocity, tilde_s, tilde_b, tilde_phi,
+        spatial_velocity, lorentz_factor, lapse, shift, spatial_metric,
+        inv_spatial_metric, sqrt_det_spatial_metric, inertial_coords,
+        dg_mesh.quadrature(2));
+  }
 }
 }  // namespace grmhd::ValenciaDivClean

src/Evolution/Systems/GrMhd/ValenciaDivClean/Sources.hpp

@@ -4,15 +4,20 @@
 #pragma once
 
 #include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Domain/Tags.hpp"
+#include "Evolution/DgSubcell/Tags/Coordinates.hpp"
 #include "Evolution/Systems/GrMhd/ValenciaDivClean/TagsDeclarations.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "NumericalAlgorithms/Spectral/Quadrature.hpp"
 #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
 #include "PointwiseFunctions/Hydro/Tags.hpp"
 #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
 #include "Utilities/TMPL.hpp"
 
 /// \cond
 class DataVector;
+template <size_t Dim>
+class Mesh;
 
 namespace Tags {
 template <typename>
@@ -23,6 +28,26 @@ struct Source;
 namespace grmhd {
 namespace ValenciaDivClean {
 namespace detail {
+void cartoon_sources_impl(
+    gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> source_tilde_s,
+    gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_b,
+
+    const Scalar<DataVector>& pressure_star,
+    const tnsr::i<DataVector, 3, Frame::Inertial>& magnetic_field_one_form,
+    const Scalar<DataVector>& magnetic_field_dot_spatial_velocity,
+
+    const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
+    const Scalar<DataVector>& tilde_phi,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
+    const Scalar<DataVector>& lorentz_factor, const Scalar<DataVector>& lapse,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
+    const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
+    const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
+    const Scalar<DataVector>& sqrt_det_spatial_metric,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
+    Spectral::Quadrature cartoon_quadrature);
+
 void sources_impl(
     gsl::not_null<Scalar<DataVector>*> source_tilde_tau,
     gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> source_tilde_s,
@@ -113,6 +138,31 @@ void sources_impl(
  * of the divergence-free (no-monopole) condition \f$\Phi = \partial_i {\tilde
  * B}^i = 0\f$ .
  *
+ * When using the cartoon method, the flux of vector-quantity conserved
+ * variables result in terms that are not flux-like but can be interpreted as
+ * sources. For spherical symmetry these extra source terms are
+ * \f{align*}
+ * S(\tilde{S}_i) &= \frac{2}{x} \alpha \sqrt{\gamma} (p + p_m) \\
+ * S(\tilde{B}^i) &= \frac{1}{x} \alpha (\gamma^{yy} + \gamma^{zz}) \tilde{\Phi}
+ * \f}
+ *
+ * and for axial symmetry they are
+ * \f{align*}
+ * S(\tilde{S}_i) &=
+ *   \frac{1}{x} \left( \alpha \sqrt{\gamma} (p + p_m) +
+ *   (\tilde{S}_z - \tilde{S}_x ) v^z_\mathrm{tr} +
+ *   \frac{\alpha}{W} (b_x - b_z) \tilde{B}^z  \right) \\
+ * S(\tilde{B}^i) &=
+ *   \frac{1}{x} \left( \alpha (\gamma^{zz}- \gamma^{zx}) \tilde{\Phi}
+ *   + \alpha (v^x - v^z)  \tilde{B}^z +
+ *   (\tilde{B}^z - \tilde{B}^x) v^z_\mathrm{tr} \right),
+ * \f}
+ *
+ * where \f$v^i_\mathrm{tr} = \alpha v^i - \beta^i\f$ is the transport
+ * velocity. These terms are associated with the $x$ component of
+ * $\tilde{S}_j$ and $\tilde{B}^j$ because that is the divergence direction we
+ * have cartoonified.
+ *
  * \note For the electron fraction side, the source term is currently set to
  * \f$S(\tilde{Y}_e) = 0\f$ where the conserved variable \f$\tilde{Y}_e\f$ is a
  * generalized electron fraction. Implementing the source term using neutrino
@@ -140,6 +190,7 @@ struct ComputeSources {
                  hydro::Tags::SpecificInternalEnergy<DataVector>,
                  hydro::Tags::LorentzFactor<DataVector>,
                  hydro::Tags::Pressure<DataVector>, gr::Tags::Lapse<DataVector>,
+                 gr::Tags::Shift<DataVector, 3>,
                  ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
                                Frame::Inertial>,
                  ::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
@@ -150,7 +201,9 @@ struct ComputeSources {
                  gr::Tags::InverseSpatialMetric<DataVector, 3>,
                  gr::Tags::SqrtDetSpatialMetric<DataVector>,
                  gr::Tags::ExtrinsicCurvature<DataVector, 3>,
-                 grmhd::ValenciaDivClean::Tags::ConstraintDampingParameter>;
+                 grmhd::ValenciaDivClean::Tags::ConstraintDampingParameter,
+                 evolution::dg::subcell::Tags::Coordinates<3, Frame::Inertial>,
+                 domain::Tags::Mesh<3>>;
 
   static void apply(
       gsl::not_null<Scalar<DataVector>*> source_tilde_tau,
@@ -169,14 +222,17 @@ struct ComputeSources {
       const Scalar<DataVector>& specific_internal_energy,
       const Scalar<DataVector>& lorentz_factor,
       const Scalar<DataVector>& pressure, const Scalar<DataVector>& lapse,
+      const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
       const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
       const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
       const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
       const tnsr::ijj<DataVector, 3, Frame::Inertial>& d_spatial_metric,
       const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
       const Scalar<DataVector>& sqrt_det_spatial_metric,
       const tnsr::ii<DataVector, 3, Frame::Inertial>& extrinsic_curvature,
-      double constraint_damping_parameter);
+      double constraint_damping_parameter,
+      const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
+      const Mesh<3>& dg_mesh);
 };
 }  // namespace ValenciaDivClean
 }  // namespace grmhd

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

@@ -62,6 +62,7 @@ target_link_libraries(
   Hydro
   HydroHelpers
   Importers
+  LinearOperators
   Utilities
   ValenciaDivClean
   )

tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Sources.py

@@ -66,6 +66,7 @@ def source_tilde_tau(
     lorentz_factor,
     pressure,
     lapse,
+    shift,
     d_lapse,
     d_shift,
     spatial_metric,
@@ -106,6 +107,7 @@ def source_tilde_s(
     lorentz_factor,
     pressure,
     lapse,
+    shift,
     d_lapse,
     d_shift,
     spatial_metric,
@@ -148,6 +150,7 @@ def source_tilde_b(
     lorentz_factor,
     pressure,
     lapse,
+    shift,
     d_lapse,
     d_shift,
     spatial_metric,
@@ -186,6 +189,7 @@ def source_tilde_phi(
     lorentz_factor,
     pressure,
     lapse,
+    shift,
     d_lapse,
     d_shift,
     spatial_metric,