Add StressEnergyTensor free function

jyoo1042 · Nov 8, 2024 · d49618a · d49618a
1 parent b0e38ad
commit d49618a
Show file tree

Hide file tree

Showing 7 changed files with 427 additions and 82 deletions.
diff --git a/src/PointwiseFunctions/Hydro/ComovingMagneticField.cpp b/src/PointwiseFunctions/Hydro/ComovingMagneticField.cpp
@@ -10,6 +10,40 @@
 
 namespace hydro {
 
+template <typename DataType>
+void comoving_magnetic_field(
+    const gsl::not_null<tnsr::A<DataType, 3>*> result,
+    const tnsr::I<DataType, 3>& spatial_velocity,
+    const tnsr::I<DataType, 3>& magnetic_field,
+    const Scalar<DataType>& magnetic_field_dot_spatial_velocity,
+    const Scalar<DataType>& lorentz_factor, const tnsr::I<DataType, 3>& shift,
+    const Scalar<DataType>& lapse) {
+  get<0>(*result) = get(magnetic_field_dot_spatial_velocity) / get(lapse);
+
+  for (size_t i = 0; i < 3; ++i) {
+    result->get(i + 1) =
+        (magnetic_field.get(i) +
+         get(lapse) * result->get(0) *
+             (get(lorentz_factor) *
+              (spatial_velocity.get(i) - shift.get(i) / get(lapse)))) /
+        get(lorentz_factor);
+  }
+}
+
+template <typename DataType>
+tnsr::A<DataType, 3> comoving_magnetic_field(
+    const tnsr::I<DataType, 3>& spatial_velocity,
+    const tnsr::I<DataType, 3>& magnetic_field,
+    const Scalar<DataType>& magnetic_field_dot_spatial_velocity,
+    const Scalar<DataType>& lorentz_factor, const tnsr::I<DataType, 3>& shift,
+    const Scalar<DataType>& lapse) {
+  tnsr::A<DataType, 3> result{};
+  comoving_magnetic_field(make_not_null(&result), spatial_velocity,
+                          magnetic_field, magnetic_field_dot_spatial_velocity,
+                          lorentz_factor, shift, lapse);
+  return result;
+}
+
 template <typename DataType>
 void comoving_magnetic_field_one_form(
     const gsl::not_null<tnsr::a<DataType, 3>*> result,
@@ -68,6 +102,16 @@ Scalar<DataType> comoving_magnetic_field_squared(
 
 #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
 #define INSTANTIATION(r, data)                                               \
+  template void comoving_magnetic_field(                                     \
+      const gsl::not_null<tnsr::A<DTYPE(data), 3>*>                          \
+          comoving_magnetic_field_result,                                    \
+      const tnsr::I<DTYPE(data), 3>&, const tnsr::I<DTYPE(data), 3>&,        \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,                \
+      const tnsr::I<DTYPE(data), 3>&, const Scalar<DTYPE(data)>&);           \
+  template tnsr::A<DTYPE(data), 3> comoving_magnetic_field(                  \
+      const tnsr::I<DTYPE(data), 3>&, const tnsr::I<DTYPE(data), 3>&,        \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,                \
+      const tnsr::I<DTYPE(data), 3>&, const Scalar<DTYPE(data)>&);           \
   template void comoving_magnetic_field_one_form(                            \
       const gsl::not_null<tnsr::a<DTYPE(data), 3>*>                          \
           comoving_magnetic_field_one_form_result,                           \

diff --git a/src/PointwiseFunctions/Hydro/ComovingMagneticField.hpp b/src/PointwiseFunctions/Hydro/ComovingMagneticField.hpp
@@ -37,6 +37,23 @@ namespace hydro {
  * $1/\sqrt{4\pi}$, following \cite Moesta2013dna .
  */
 template <typename DataType>
+void comoving_magnetic_field(
+    const gsl::not_null<tnsr::A<DataType, 3>*> result,
+    const tnsr::I<DataType, 3>& spatial_velocity,
+    const tnsr::I<DataType, 3>& magnetic_field,
+    const Scalar<DataType>& magnetic_field_dot_spatial_velocity,
+    const Scalar<DataType>& lorentz_factor, const tnsr::I<DataType, 3>& shift,
+    const Scalar<DataType>& lapse);
+
+template <typename DataType>
+tnsr::A<DataType, 3> comoving_magnetic_field(
+    const tnsr::I<DataType, 3>& spatial_velocity,
+    const tnsr::I<DataType, 3>& magnetic_field,
+    const Scalar<DataType>& magnetic_field_dot_spatial_velocity,
+    const Scalar<DataType>& lorentz_factor, const tnsr::I<DataType, 3>& shift,
+    const Scalar<DataType>& lapse);
+
+template <typename DataType>
 void comoving_magnetic_field_one_form(
     const gsl::not_null<tnsr::a<DataType, 3>*> result,
     const tnsr::i<DataType, 3>& spatial_velocity_one_form,

diff --git a/src/PointwiseFunctions/Hydro/StressEnergy.cpp b/src/PointwiseFunctions/Hydro/StressEnergy.cpp
@@ -4,11 +4,40 @@
 #include "PointwiseFunctions/Hydro/StressEnergy.hpp"
 
 #include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
+#include "PointwiseFunctions/GeneralRelativity/InverseSpacetimeMetric.hpp"
+#include "PointwiseFunctions/Hydro/ComovingMagneticField.hpp"
 #include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Gsl.hpp"
 
+namespace {
+template <typename DataType>
+tnsr::A<DataType, 3> four_velocity(const tnsr::I<DataType, 3>& spatial_velocity,
+                                   const tnsr::I<DataType, 3>& shift,
+                                   const Scalar<DataType>& lorentz_factor,
+                                   const Scalar<DataType>& lapse) {
+  tnsr::A<DataType, 3> result{};
+  get<0>(result) = get(lorentz_factor) / get(lapse);
+  for (size_t i = 0; i < 3; ++i) {
+    result.get(i + 1) = get(lorentz_factor) *
+                        (spatial_velocity.get(i) - shift.get(i) / get(lapse));
+  }
+  return result;
+}
+
+template tnsr::A<DataVector, 3> four_velocity(const tnsr::I<DataVector, 3>&,
+                                              const tnsr::I<DataVector, 3>&,
+                                              const Scalar<DataVector>&,
+                                              const Scalar<DataVector>&);
+
+template tnsr::A<double, 3> four_velocity(const tnsr::I<double, 3>&,
+                                          const tnsr::I<double, 3>&,
+                                          const Scalar<double>&,
+                                          const Scalar<double>&);
+}  // namespace
+
 namespace hydro {
 
 template <typename DataType>
@@ -70,23 +99,111 @@ void stress_trace(gsl::not_null<Scalar<DataType>*> result,
           (square(get(lorentz_factor)) * get(spatial_velocity_squared) + 1.);
 }
 
+template <typename DataType>
+void stress_energy_tensor(
+    gsl::not_null<tnsr::AA<DataType, 3>*> result,
+    const Scalar<DataType>& rest_mass_density,
+    const Scalar<DataType>& specific_internal_energy,
+    const Scalar<DataType>& pressure, const Scalar<DataType>& lorentz_factor,
+    const Scalar<DataType>& lapse,
+    const Scalar<DataType>& comoving_magnetic_field_magnitude,
+    const tnsr::I<DataType, 3>& spatial_velocity,
+    const tnsr::I<DataType, 3>& shift,
+    const tnsr::I<DataType, 3>& magnetic_field,
+    const tnsr::ii<DataType, 3>& spatial_metric,
+    const tnsr::II<DataType, 3>& inverse_spatial_metric) {
+  const auto inverse_spacetime_metric =
+      gr::inverse_spacetime_metric(lapse, shift, inverse_spatial_metric);
+
+  const auto magnetic_field_dot_spatial_velocity =
+      dot_product(magnetic_field, spatial_velocity, spatial_metric);
+
+  const auto comoving_magnetic_field_v = comoving_magnetic_field(
+      spatial_velocity, magnetic_field, magnetic_field_dot_spatial_velocity,
+      lorentz_factor, shift, lapse);
+
+  const auto four_velocity_v =
+      four_velocity(spatial_velocity, shift, lorentz_factor, lapse);
+
+  const auto rho_h_star =
+      (get(rest_mass_density) +
+       get(rest_mass_density) * get(specific_internal_energy)) +
+      get(pressure) + square(get(comoving_magnetic_field_magnitude));
+
+  const auto p_star =
+      get(pressure) + square(get(comoving_magnetic_field_magnitude)) / 2.;
+
+  get<0, 0>(*result) = (rho_h_star * square(get<0>(four_velocity_v))) +
+                       (p_star * get<0, 0>(inverse_spacetime_metric)) -
+                       (square(get<0>(comoving_magnetic_field_v)));
+
+  get<1, 0>(*result) =
+      (rho_h_star * get<1>(four_velocity_v) * get<0>(four_velocity_v)) +
+      (p_star * get<1, 0>(inverse_spacetime_metric)) -
+      (get<1>(comoving_magnetic_field_v) * get<0>(comoving_magnetic_field_v));
+
+  get<1, 1>(*result) = (rho_h_star * square(get<1>(four_velocity_v))) +
+                       (p_star * get<1, 1>(inverse_spacetime_metric)) -
+                       (square(get<1>(comoving_magnetic_field_v)));
+
+  get<2, 0>(*result) =
+      (rho_h_star * get<2>(four_velocity_v) * get<0>(four_velocity_v)) +
+      (p_star * get<2, 0>(inverse_spacetime_metric)) -
+      (get<2>(comoving_magnetic_field_v) * get<0>(comoving_magnetic_field_v));
+
+  get<2, 1>(*result) =
+      (rho_h_star * get<2>(four_velocity_v) * get<1>(four_velocity_v)) +
+      (p_star * get<2, 1>(inverse_spacetime_metric)) -
+      (get<2>(comoving_magnetic_field_v) * get<1>(comoving_magnetic_field_v));
+
+  get<2, 2>(*result) = (rho_h_star * square(get<2>(four_velocity_v))) +
+                       (p_star * get<2, 2>(inverse_spacetime_metric)) -
+                       (square(get<2>(comoving_magnetic_field_v)));
+
+  get<3, 0>(*result) =
+      (rho_h_star * get<3>(four_velocity_v) * get<0>(four_velocity_v)) +
+      (p_star * get<3, 0>(inverse_spacetime_metric)) -
+      (get<3>(comoving_magnetic_field_v) * get<0>(comoving_magnetic_field_v));
+
+  get<3, 1>(*result) =
+      (rho_h_star * get<3>(four_velocity_v) * get<1>(four_velocity_v)) +
+      (p_star * get<3, 1>(inverse_spacetime_metric)) -
+      (get<3>(comoving_magnetic_field_v) * get<1>(comoving_magnetic_field_v));
+
+  get<3, 2>(*result) =
+      (rho_h_star * get<3>(four_velocity_v) * get<2>(four_velocity_v)) +
+      (p_star * get<3, 2>(inverse_spacetime_metric)) -
+      (get<3>(comoving_magnetic_field_v) * get<2>(comoving_magnetic_field_v));
+
+  get<3, 3>(*result) = (rho_h_star * square(get<3>(four_velocity_v))) +
+                       (p_star * get<3, 3>(inverse_spacetime_metric)) -
+                       (square(get<3>(comoving_magnetic_field_v)));
+}
+
 #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
-#define INSTANTIATION(r, data)                                             \
-  template void energy_density(                                            \
-      gsl::not_null<Scalar<DTYPE(data)>*>, const Scalar<DTYPE(data)>&,     \
-      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,              \
-      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,              \
-      const Scalar<DTYPE(data)>&);                                         \
-  template void momentum_density(                                          \
-      gsl::not_null<tnsr::I<DTYPE(data), 3>*>, const Scalar<DTYPE(data)>&, \
-      const Scalar<DTYPE(data)>&, const tnsr::I<DTYPE(data), 3>&,          \
-      const Scalar<DTYPE(data)>&, const tnsr::I<DTYPE(data), 3>&,          \
-      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&);             \
-  template void stress_trace(                                              \
-      gsl::not_null<Scalar<DTYPE(data)>*>, const Scalar<DTYPE(data)>&,     \
-      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,              \
-      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,              \
-      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&);
+#define INSTANTIATION(r, data)                                              \
+  template void energy_density(                                             \
+      gsl::not_null<Scalar<DTYPE(data)>*>, const Scalar<DTYPE(data)>&,      \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,               \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,               \
+      const Scalar<DTYPE(data)>&);                                          \
+  template void momentum_density(                                           \
+      gsl::not_null<tnsr::I<DTYPE(data), 3>*>, const Scalar<DTYPE(data)>&,  \
+      const Scalar<DTYPE(data)>&, const tnsr::I<DTYPE(data), 3>&,           \
+      const Scalar<DTYPE(data)>&, const tnsr::I<DTYPE(data), 3>&,           \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&);              \
+  template void stress_trace(                                               \
+      gsl::not_null<Scalar<DTYPE(data)>*>, const Scalar<DTYPE(data)>&,      \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,               \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,               \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&);              \
+  template void stress_energy_tensor(                                       \
+      gsl::not_null<tnsr::AA<DTYPE(data), 3>*>, const Scalar<DTYPE(data)>&, \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,               \
+      const Scalar<DTYPE(data)>&, const Scalar<DTYPE(data)>&,               \
+      const Scalar<DTYPE(data)>&, const tnsr::I<DTYPE(data), 3>&,           \
+      const tnsr::I<DTYPE(data), 3>&, const tnsr::I<DTYPE(data), 3>&,       \
+      const tnsr::ii<DTYPE(data), 3>&, const tnsr::II<DTYPE(data), 3>&);
 
 GENERATE_INSTANTIATIONS(INSTANTIATION, (double, DataVector))
 

diff --git a/src/PointwiseFunctions/Hydro/StressEnergy.hpp b/src/PointwiseFunctions/Hydro/StressEnergy.hpp
@@ -132,4 +132,28 @@ void stress_trace(gsl::not_null<Scalar<DataType>*> result,
                   const Scalar<DataType>& magnetic_field_dot_spatial_velocity,
                   const Scalar<DataType>& comoving_magnetic_field_squared);
 
+/*!
+ * \brief Stress Energy Tesnor, $S^{ab}=
+ * (\rho h)^{*} u^a u ^b + p^{*} g^{ab} - b^{a} b^{b}$,
+ *
+ * where $(\rho h)^{*} = \rho h + b^{2}$ and $p^{*} = p + b^{2}/2$
+ * are the enthalpy density and fluid pressure augmented by contributions of
+ * magnetic pressure $p_{mag}$ = b^{2}/2, respectively.
+ *
+ * $b$ refers to magnetic field measured in the comoving frame of the fluid
+ * $b^{a} = ^{*}F^{ab} u_{b}$.
+ */
+template <typename DataType>
+void stress_energy_tensor(
+    gsl::not_null<tnsr::AA<DataType, 3>*> result,
+    const Scalar<DataType>& rest_mass_density,
+    const Scalar<DataType>& specific_internal_energy,
+    const Scalar<DataType>& pressure, const Scalar<DataType>& lorentz_factor,
+    const Scalar<DataType>& lapse,
+    const Scalar<DataType>& comoving_magnetic_field_magnitude,
+    const tnsr::I<DataType, 3>& spatial_velocity,
+    const tnsr::I<DataType, 3>& shift,
+    const tnsr::I<DataType, 3>& magnetic_field,
+    const tnsr::ii<DataType, 3>& spatial_metric,
+    const tnsr::II<DataType, 3>& inverse_spatial_metric);
 }  // namespace hydro