Add StressEnergyTensor free function

jyoo1042 · Jan 17, 2025 · 018e625 · 018e625
1 parent 918181b
commit 018e625
Show file tree

Hide file tree

Showing 8 changed files with 358 additions and 85 deletions.
diff --git a/src/PointwiseFunctions/Hydro/CMakeLists.txt b/src/PointwiseFunctions/Hydro/CMakeLists.txt
@@ -51,10 +51,10 @@ target_link_libraries(
   DataStructures
   ErrorHandling
   H5
+  GeneralRelativity
   Options
   Serialization
   INTERFACE
-  GeneralRelativity
   )
 
 add_subdirectory(EquationsOfState)

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,8 +37,25 @@ namespace hydro {
  * $1/\sqrt{4\pi}$, following \cite Moesta2013dna .
  */
 template <typename DataType>
+void comoving_magnetic_field(
+    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,
+    gsl::not_null<tnsr::a<DataType, 3>*> result,
     const tnsr::i<DataType, 3>& spatial_velocity_one_form,
     const tnsr::i<DataType, 3>& magnetic_field_one_form,
     const Scalar<DataType>& magnetic_field_dot_spatial_velocity,
@@ -55,7 +72,7 @@ tnsr::a<DataType, 3> comoving_magnetic_field_one_form(
 
 template <typename DataType>
 void comoving_magnetic_field_squared(
-    const gsl::not_null<Scalar<DataType>*> result,
+    gsl::not_null<Scalar<DataType>*> result,
     const Scalar<DataType>& magnetic_field_squared,
     const Scalar<DataType>& magnetic_field_dot_spatial_velocity,
     const Scalar<DataType>& lorentz_factor);

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,78 @@ 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.;
+
+  for (size_t i = 0; i < 4; ++i) {
+    result->get(i, i) = rho_h_star * square(four_velocity_v.get(i)) +
+                        p_star * inverse_spacetime_metric.get(i, i) -
+                        square(comoving_magnetic_field_v.get(i));
+
+    for (size_t j = i + 1; j < 4; ++j) {
+      result->get(i, j) =
+          rho_h_star * four_velocity_v.get(i) * four_velocity_v.get(j) +
+          p_star * inverse_spacetime_metric.get(i, j) -
+          comoving_magnetic_field_v.get(i) * comoving_magnetic_field_v.get(j);
+    }
+  }
+}
+
 #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, $T^{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