Merge pull request #4351 from yoonso0-0/ffe/add_source_terms

Compute source terms in the ForceFree evolution system

src/Evolution/Systems/ForceFree/CMakeLists.txt

@@ -10,6 +10,7 @@ spectre_target_sources(
   PRIVATE
   Characteristics.cpp
   Fluxes.cpp
+  Sources.cpp
   )
 
 spectre_target_headers(
@@ -18,6 +19,7 @@ spectre_target_headers(
   HEADERS
   Characteristics.hpp
   Fluxes.hpp
+  Sources.hpp
   Tags.hpp
   )
 

src/Evolution/Systems/ForceFree/Sources.cpp

@@ -0,0 +1,145 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/ForceFree/Sources.hpp"
+
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Christoffel.hpp"
+#include "PointwiseFunctions/GeneralRelativity/IndexManipulation.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace ForceFree {
+
+namespace detail {
+void sources_impl(
+    const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+        source_tilde_e,
+    const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+        source_tilde_b,
+    const gsl::not_null<Scalar<DataVector>*> source_tilde_psi,
+    const gsl::not_null<Scalar<DataVector>*> source_tilde_phi,
+
+    // temp variables
+    const tnsr::I<DataVector, 3, Frame::Inertial>&
+        trace_spatial_christoffel_second,
+
+    // EM args
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_e,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
+    const Scalar<DataVector>& tilde_psi, const Scalar<DataVector>& tilde_phi,
+    const Scalar<DataVector>& tilde_q,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_current_density,
+    const double kappa_psi, const double kappa_phi,
+    // GR args
+    const Scalar<DataVector>& lapse,
+    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>& inv_spatial_metric,
+    const Scalar<DataVector>& sqrt_det_spatial_metric,
+    const tnsr::ii<DataVector, 3, Frame::Inertial>& extrinsic_curvature) {
+  // S(\tilde{E}^i)
+  raise_or_lower_index(source_tilde_e, d_lapse, inv_spatial_metric);
+  for (size_t i = 0; i < 3; ++i) {
+    source_tilde_e->get(i) -=
+        get(lapse) * trace_spatial_christoffel_second.get(i);
+    source_tilde_e->get(i) *= get(tilde_psi);
+    source_tilde_e->get(i) -= get(lapse) * get(sqrt_det_spatial_metric) *
+                              spatial_current_density.get(i);
+    for (size_t m = 0; m < 3; ++m) {
+      source_tilde_e->get(i) -= tilde_e.get(m) * d_shift.get(m, i);
+    }
+  }
+
+  // S(\tilde{B}^i)
+  raise_or_lower_index(source_tilde_b, d_lapse, inv_spatial_metric);
+  for (size_t i = 0; i < 3; ++i) {
+    source_tilde_b->get(i) -=
+        get(lapse) * trace_spatial_christoffel_second.get(i);
+    source_tilde_b->get(i) *= get(tilde_phi);
+    for (size_t m = 0; m < 3; ++m) {
+      source_tilde_b->get(i) -= tilde_b.get(m) * d_shift.get(m, i);
+    }
+  }
+
+  // S(\tilde{\psi})
+  trace(source_tilde_psi, extrinsic_curvature, inv_spatial_metric);
+  get(*source_tilde_psi) += kappa_psi;
+  get(*source_tilde_psi) *= -1.0 * get(tilde_psi);
+  get(*source_tilde_psi) += get(tilde_q);
+  get(*source_tilde_psi) *= get(lapse);
+  for (size_t m = 0; m < 3; ++m) {
+    get(*source_tilde_psi) += tilde_e.get(m) * d_lapse.get(m);
+  }
+
+  // S(\tilde{\phi})
+  trace(source_tilde_phi, extrinsic_curvature, inv_spatial_metric);
+  get(*source_tilde_phi) += kappa_phi;
+  get(*source_tilde_phi) *= -1.0 * get(lapse) * get(tilde_phi);
+  for (size_t m = 0; m < 3; ++m) {
+    get(*source_tilde_phi) += tilde_b.get(m) * d_lapse.get(m);
+  }
+}
+}  // namespace detail
+
+void Sources::apply(
+    const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+        source_tilde_e,
+    const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
+        source_tilde_b,
+    const gsl::not_null<Scalar<DataVector>*> source_tilde_psi,
+    const gsl::not_null<Scalar<DataVector>*> source_tilde_phi,
+    // EM variables
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_e,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
+    const Scalar<DataVector>& tilde_psi, const Scalar<DataVector>& tilde_phi,
+    const Scalar<DataVector>& tilde_q,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_current_density,
+    const double kappa_psi, const double kappa_phi,
+    // GR variables
+    const Scalar<DataVector>& lapse,
+    const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
+    const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
+    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) {
+  // temp variable to store metric derivative quantities
+  Variables<tmpl::list<
+      gr::Tags::SpatialChristoffelFirstKind<3, Frame::Inertial, DataVector>,
+      gr::Tags::SpatialChristoffelSecondKind<3, Frame::Inertial, DataVector>,
+      gr::Tags::TraceSpatialChristoffelSecondKind<3, Frame::Inertial,
+                                                  DataVector>>>
+      temp_tensors(get<0>(tilde_e).size());
+
+  // compute the product \gamma^jk \Gamma^i_{jk}
+  auto& spatial_christoffel_first_kind = get<
+      gr::Tags::SpatialChristoffelFirstKind<3, Frame::Inertial, DataVector>>(
+      temp_tensors);
+  gr::christoffel_first_kind(make_not_null(&spatial_christoffel_first_kind),
+                             d_spatial_metric);
+  auto& spatial_christoffel_second_kind = get<
+      gr::Tags::SpatialChristoffelSecondKind<3, Frame::Inertial, DataVector>>(
+      temp_tensors);
+  raise_or_lower_first_index(make_not_null(&spatial_christoffel_second_kind),
+                             spatial_christoffel_first_kind,
+                             inv_spatial_metric);
+  auto& trace_spatial_christoffel_second =
+      get<gr::Tags::TraceSpatialChristoffelSecondKind<3, Frame::Inertial,
+                                                      DataVector>>(
+          temp_tensors);
+  trace_last_indices(make_not_null(&trace_spatial_christoffel_second),
+                     spatial_christoffel_second_kind, inv_spatial_metric);
+
+  detail::sources_impl(source_tilde_e, source_tilde_b, source_tilde_psi,
+                       source_tilde_phi, trace_spatial_christoffel_second,
+                       tilde_e, tilde_b, tilde_psi, tilde_phi, tilde_q,
+                       spatial_current_density, kappa_psi, kappa_phi, lapse,
+                       d_lapse, d_shift, inv_spatial_metric,
+                       sqrt_det_spatial_metric, extrinsic_curvature);
+}
+
+}  // namespace ForceFree

src/Evolution/Systems/ForceFree/Sources.hpp

@@ -0,0 +1,119 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include "DataStructures/DataBox/Prefixes.hpp"
+#include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Evolution/Systems/ForceFree/Tags.hpp"
+#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+class DataVector;
+/// \endcond
+
+namespace ForceFree {
+
+namespace detail {
+void sources_impl(
+    gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_e,
+    gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_b,
+    gsl::not_null<Scalar<DataVector>*> source_tilde_psi,
+    gsl::not_null<Scalar<DataVector>*> source_tilde_phi,
+
+    // temp variables
+    const tnsr::I<DataVector, 3, Frame::Inertial>&
+        trace_spatial_christoffel_second,
+
+    // EM args
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_e,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
+    const Scalar<DataVector>& tilde_psi, const Scalar<DataVector>& tilde_phi,
+    const Scalar<DataVector>& tilde_q,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_current_density,
+    const double kappa_psi, const double kappa_phi,
+    // GR args
+    const Scalar<DataVector>& lapse,
+    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>& inv_spatial_metric,
+    const Scalar<DataVector>& sqrt_det_spatial_metric,
+    const tnsr::ii<DataVector, 3, Frame::Inertial>& extrinsic_curvature);
+}  // namespace detail
+
+/*!
+ * \brief Compute the source terms for the GRFFE system with divergence
+ * cleaning.
+ *
+ * \f{align*}
+ *  S(\tilde{E}^i) &= -\alpha \sqrt{\gamma} J^i - \tilde{E}^j \partial_j \beta^i
+ *    + \tilde{\psi} ( \gamma^{ij} \partial_j \alpha - \alpha \gamma^{jk}
+ *    \Gamma^i_{jk} ) \\
+ *  S(\tilde{B}^i) &= -\tilde{B}^j \partial_j \beta^i + \tilde{\phi} (
+ *    \gamma^{ij} \partial_j \alpha - \alpha \gamma^{jk} \Gamma^i_{jk} ) \\
+ *  S(\tilde{\psi}) &= \tilde{E}^k \partial_k \alpha + \alpha \tilde{q} - \alpha
+ *    \tilde{\phi} ( K + \kappa_\phi ) \\
+ *  S(\tilde{\phi}) &= \tilde{B}^k \partial_k \alpha - \alpha \tilde{\phi} (K +
+ *    \kappa_\phi ) \\
+ *  S(\tilde{q}) &= 0
+ * \f}
+ *
+ * where the conserved variables \f$\tilde{E}^i, \tilde{B}^i, \tilde{\psi},
+ * \tilde{\phi}, \tilde{q}\f$ are densitized electric field, magnetic field,
+ * magnetic divergence cleaning field, electric divergence cleaning field, and
+ * electric charge density.
+ *
+ * \f$J^i\f$ is the spatial electric current density, \f$\alpha\f$ is the lapse,
+ * \f$\beta^i\f$ is the shift, \f$\gamma^{ij}\f$ is the spatial metric,
+ * \f$\gamma\f$ is the determinant of spatial metric, \f$\Gamma^i_{jk}\f$ is the
+ * spatial Christoffel symbol, \f$K\f$ is the trace of extrinsic curvature.
+ * \f$\kappa_\phi\f$ and \f$\kappa_\psi\f$ are damping parameters associated
+ * with divergence cleaning of magnetic and electric fields, respectively.
+ *
+ */
+struct Sources {
+  using return_tags =
+      tmpl::list<::Tags::Source<Tags::TildeE>, ::Tags::Source<Tags::TildeB>,
+                 ::Tags::Source<Tags::TildePsi>,
+                 ::Tags::Source<Tags::TildePhi>>;
+
+  using argument_tags = tmpl::list<
+      // EM variables
+      Tags::TildeE, Tags::TildeB, Tags::TildePsi, Tags::TildePhi, Tags::TildeQ,
+      Tags::SpatialCurrentDensity, Tags::KappaPsi, Tags::KappaPhi,
+      // GR variables
+      gr::Tags::Lapse<>,
+      ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
+                    Frame::Inertial>,
+      ::Tags::deriv<gr::Tags::Shift<3, Frame::Inertial, DataVector>,
+                    tmpl::size_t<3>, Frame::Inertial>,
+      ::Tags::deriv<gr::Tags::SpatialMetric<3, Frame::Inertial, DataVector>,
+                    tmpl::size_t<3>, Frame::Inertial>,
+      gr::Tags::InverseSpatialMetric<3>, gr::Tags::SqrtDetSpatialMetric<>,
+      gr::Tags::ExtrinsicCurvature<3, Frame::Inertial, DataVector>>;
+
+  static void apply(
+      gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_e,
+      gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_b,
+      gsl::not_null<Scalar<DataVector>*> source_tilde_psi,
+      gsl::not_null<Scalar<DataVector>*> source_tilde_phi,
+      // EM variables
+      const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_e,
+      const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
+      const Scalar<DataVector>& tilde_psi, const Scalar<DataVector>& tilde_phi,
+      const Scalar<DataVector>& tilde_q,
+      const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_current_density,
+      const double kappa_psi, const double kappa_phi,
+      // GR variables
+      const Scalar<DataVector>& lapse,
+      const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
+      const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
+      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);
+};
+
+}  // namespace ForceFree

tests/Unit/Evolution/Systems/ForceFree/CMakeLists.txt

@@ -7,6 +7,7 @@ set(LIBRARY_SOURCES
   BoundaryConditions/Test_Periodic.cpp
   Test_Characteristics.cpp
   Test_Fluxes.cpp
+  Test_Sources.cpp
   Test_Tags.cpp
   )
 

tests/Unit/Evolution/Systems/ForceFree/Sources.py

@@ -0,0 +1,52 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+import numpy as np
+
+
+def trace_spatial_Gamma_second_kind(inv_spatial_metric, d_spatial_metric):
+    # returns \gamma^{jk} \Gamma^i_{jk}
+    term_one = np.einsum("ab, iab", inv_spatial_metric, d_spatial_metric)
+    term_two = np.einsum("ab, aib", inv_spatial_metric, d_spatial_metric)
+    return -0.5 * np.einsum("ia, a", inv_spatial_metric, term_one) + np.einsum(
+        "ia, a", inv_spatial_metric, term_two)
+
+
+def source_tilde_e(tilde_e, tilde_b, tilde_psi, tilde_phi, tilde_q,
+                   current_density, kappa_psi, kappa_phi, lapse, d_lapse,
+                   d_shift, d_spatial_metric, inv_spatial_metric,
+                   sqrt_det_spatial_metric, extrinsic_curvature):
+    return -lapse * sqrt_det_spatial_metric * current_density - np.einsum(
+        "a, ai", tilde_e,
+        d_shift) + tilde_psi * (np.einsum("a, ia", d_lapse, inv_spatial_metric)
+                                - lapse * trace_spatial_Gamma_second_kind(
+                                    inv_spatial_metric, d_spatial_metric))
+
+
+def source_tilde_b(tilde_e, tilde_b, tilde_psi, tilde_phi, tilde_q,
+                   current_density, kappa_psi, kappa_phi, lapse, d_lapse,
+                   d_shift, d_spatial_metric, inv_spatial_metric,
+                   sqrt_det_spatial_metric, extrinsic_curvature):
+    return -np.einsum("a, ai", tilde_b, d_shift) + tilde_phi * (
+        np.einsum("a, ia", d_lapse, inv_spatial_metric) - lapse *
+        trace_spatial_Gamma_second_kind(inv_spatial_metric, d_spatial_metric))
+
+
+def source_tilde_phi(tilde_e, tilde_b, tilde_psi, tilde_phi, tilde_q,
+                     current_density, kappa_psi, kappa_phi, lapse, d_lapse,
+                     d_shift, d_spatial_metric, inv_spatial_metric,
+                     sqrt_det_spatial_metric, extrinsic_curvature):
+
+    return np.einsum(
+        "a, a", tilde_b, d_lapse) - lapse * tilde_phi * (np.einsum(
+            "ab, ab", inv_spatial_metric, extrinsic_curvature) + kappa_phi)
+
+
+def source_tilde_psi(tilde_e, tilde_b, tilde_psi, tilde_phi, tilde_q,
+                     current_density, kappa_psi, kappa_phi, lapse, d_lapse,
+                     d_shift, d_spatial_metric, inv_spatial_metric,
+                     sqrt_det_spatial_metric, extrinsic_curvature):
+    return np.einsum(
+        "a, a", tilde_e,
+        d_lapse) + lapse * tilde_q - lapse * tilde_psi * (np.einsum(
+            "ab, ab", inv_spatial_metric, extrinsic_curvature) + kappa_psi)

tests/Unit/Evolution/Systems/ForceFree/Test_Sources.cpp

@@ -0,0 +1,20 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include "DataStructures/DataVector.hpp"
+#include "Evolution/Systems/ForceFree/Sources.hpp"
+#include "Framework/CheckWithRandomValues.hpp"
+#include "Framework/SetupLocalPythonEnvironment.hpp"
+
+SPECTRE_TEST_CASE("Unit.Evolution.Systems.ForceFree.Sources",
+                  "[Unit][Evolution]") {
+  pypp::SetupLocalPythonEnvironment local_python_env{
+      "Evolution/Systems/ForceFree"};
+
+  pypp::check_with_random_values<1>(&ForceFree::Sources::apply, "Sources",
+                                    {"source_tilde_e", "source_tilde_b",
+                                     "source_tilde_psi", "source_tilde_phi"},
+                                    {{{0.0, 1.0}}}, DataVector{5});
+}