From d6039ad0dd87d1668b9dccf0be2bc54e27bf6f6b Mon Sep 17 00:00:00 2001 From: Jooheon Yoo Date: Wed, 30 Oct 2024 11:00:57 -0400 Subject: [PATCH] Add StressEnergyTensor free function --- .../Hydro/ComovingMagneticField.cpp | 44 ++++ .../Hydro/ComovingMagneticField.hpp | 17 ++ src/PointwiseFunctions/Hydro/StressEnergy.cpp | 149 +++++++++++-- src/PointwiseFunctions/Hydro/StressEnergy.hpp | 24 +++ .../Hydro/Python/Test_StressEnergy.py | 199 ++++++++++++++++++ .../PointwiseFunctions/Hydro/TestFunctions.py | 59 ------ .../Hydro/Test_StressEnergy.cpp | 17 +- 7 files changed, 427 insertions(+), 82 deletions(-) diff --git a/src/PointwiseFunctions/Hydro/ComovingMagneticField.cpp b/src/PointwiseFunctions/Hydro/ComovingMagneticField.cpp index 98a08346d220..d6b6be3a3b62 100644 --- a/src/PointwiseFunctions/Hydro/ComovingMagneticField.cpp +++ b/src/PointwiseFunctions/Hydro/ComovingMagneticField.cpp @@ -10,6 +10,40 @@ namespace hydro { +template +void comoving_magnetic_field( + const gsl::not_null*> result, + const tnsr::I& spatial_velocity, + const tnsr::I& magnetic_field, + const Scalar& magnetic_field_dot_spatial_velocity, + const Scalar& lorentz_factor, const tnsr::I& shift, + const Scalar& 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 +tnsr::A comoving_magnetic_field( + const tnsr::I& spatial_velocity, + const tnsr::I& magnetic_field, + const Scalar& magnetic_field_dot_spatial_velocity, + const Scalar& lorentz_factor, const tnsr::I& shift, + const Scalar& lapse) { + tnsr::A result{}; + comoving_magnetic_field(make_not_null(&result), spatial_velocity, + magnetic_field, magnetic_field_dot_spatial_velocity, + lorentz_factor, shift, lapse); + return result; +} + template void comoving_magnetic_field_one_form( const gsl::not_null*> result, @@ -68,6 +102,16 @@ Scalar 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*> \ + comoving_magnetic_field_result, \ + const tnsr::I&, const tnsr::I&, \ + const Scalar&, const Scalar&, \ + const tnsr::I&, const Scalar&); \ + template tnsr::A comoving_magnetic_field( \ + const tnsr::I&, const tnsr::I&, \ + const Scalar&, const Scalar&, \ + const tnsr::I&, const Scalar&); \ template void comoving_magnetic_field_one_form( \ const gsl::not_null*> \ comoving_magnetic_field_one_form_result, \ diff --git a/src/PointwiseFunctions/Hydro/ComovingMagneticField.hpp b/src/PointwiseFunctions/Hydro/ComovingMagneticField.hpp index ab83cc152e3b..6252c1fdba50 100644 --- 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 +void comoving_magnetic_field( + const gsl::not_null*> result, + const tnsr::I& spatial_velocity, + const tnsr::I& magnetic_field, + const Scalar& magnetic_field_dot_spatial_velocity, + const Scalar& lorentz_factor, const tnsr::I& shift, + const Scalar& lapse); + +template +tnsr::A comoving_magnetic_field( + const tnsr::I& spatial_velocity, + const tnsr::I& magnetic_field, + const Scalar& magnetic_field_dot_spatial_velocity, + const Scalar& lorentz_factor, const tnsr::I& shift, + const Scalar& lapse); + +template void comoving_magnetic_field_one_form( const gsl::not_null*> result, const tnsr::i& spatial_velocity_one_form, diff --git a/src/PointwiseFunctions/Hydro/StressEnergy.cpp b/src/PointwiseFunctions/Hydro/StressEnergy.cpp index 6e8cdf8fdd08..7de1c0fad9da 100644 --- 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 +tnsr::A four_velocity(const tnsr::I& spatial_velocity, + const tnsr::I& shift, + const Scalar& lorentz_factor, + const Scalar& lapse) { + tnsr::A 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 four_velocity(const tnsr::I&, + const tnsr::I&, + const Scalar&, + const Scalar&); + +template tnsr::A four_velocity(const tnsr::I&, + const tnsr::I&, + const Scalar&, + const Scalar&); +} // namespace + namespace hydro { template @@ -70,23 +99,111 @@ void stress_trace(gsl::not_null*> result, (square(get(lorentz_factor)) * get(spatial_velocity_squared) + 1.); } +template +void stress_energy_tensor( + gsl::not_null*> result, + const Scalar& rest_mass_density, + const Scalar& specific_internal_energy, + const Scalar& pressure, const Scalar& lorentz_factor, + const Scalar& lapse, + const Scalar& comoving_magnetic_field_magnitude, + const tnsr::I& spatial_velocity, + const tnsr::I& shift, + const tnsr::I& magnetic_field, + const tnsr::ii& spatial_metric, + const tnsr::II& 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*>, const Scalar&, \ - const Scalar&, const Scalar&, \ - const Scalar&, const Scalar&, \ - const Scalar&); \ - template void momentum_density( \ - gsl::not_null*>, const Scalar&, \ - const Scalar&, const tnsr::I&, \ - const Scalar&, const tnsr::I&, \ - const Scalar&, const Scalar&); \ - template void stress_trace( \ - gsl::not_null*>, const Scalar&, \ - const Scalar&, const Scalar&, \ - const Scalar&, const Scalar&, \ - const Scalar&, const Scalar&); +#define INSTANTIATION(r, data) \ + template void energy_density( \ + gsl::not_null*>, const Scalar&, \ + const Scalar&, const Scalar&, \ + const Scalar&, const Scalar&, \ + const Scalar&); \ + template void momentum_density( \ + gsl::not_null*>, const Scalar&, \ + const Scalar&, const tnsr::I&, \ + const Scalar&, const tnsr::I&, \ + const Scalar&, const Scalar&); \ + template void stress_trace( \ + gsl::not_null*>, const Scalar&, \ + const Scalar&, const Scalar&, \ + const Scalar&, const Scalar&, \ + const Scalar&, const Scalar&); \ + template void stress_energy_tensor( \ + gsl::not_null*>, const Scalar&, \ + const Scalar&, const Scalar&, \ + const Scalar&, const Scalar&, \ + const Scalar&, const tnsr::I&, \ + const tnsr::I&, const tnsr::I&, \ + const tnsr::ii&, const tnsr::II&); GENERATE_INSTANTIATIONS(INSTANTIATION, (double, DataVector)) diff --git a/src/PointwiseFunctions/Hydro/StressEnergy.hpp b/src/PointwiseFunctions/Hydro/StressEnergy.hpp index 6d5f8be6d6f2..f9566bdaf5e1 100644 --- a/src/PointwiseFunctions/Hydro/StressEnergy.hpp +++ b/src/PointwiseFunctions/Hydro/StressEnergy.hpp @@ -132,4 +132,28 @@ void stress_trace(gsl::not_null*> result, const Scalar& magnetic_field_dot_spatial_velocity, const Scalar& 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 +void stress_energy_tensor( + gsl::not_null*> result, + const Scalar& rest_mass_density, + const Scalar& specific_internal_energy, + const Scalar& pressure, const Scalar& lorentz_factor, + const Scalar& lapse, + const Scalar& comoving_magnetic_field_magnitude, + const tnsr::I& spatial_velocity, + const tnsr::I& shift, + const tnsr::I& magnetic_field, + const tnsr::ii& spatial_metric, + const tnsr::II& inverse_spatial_metric); } // namespace hydro diff --git a/tests/Unit/PointwiseFunctions/Hydro/Python/Test_StressEnergy.py b/tests/Unit/PointwiseFunctions/Hydro/Python/Test_StressEnergy.py index 32e61c9134fa..5cec1ade5edb 100644 --- a/tests/Unit/PointwiseFunctions/Hydro/Python/Test_StressEnergy.py +++ b/tests/Unit/PointwiseFunctions/Hydro/Python/Test_StressEnergy.py @@ -1,2 +1,201 @@ # Distributed under the MIT License. # See LICENSE.txt for details. + +import numpy as np + + +def inverse_spacetime_metric(lapse, shift, inverse_spatial_metric): + result = np.zeros([4, 4]) + result[0, 0] = -1.0 / (lapse * lapse) + for i in range(3): + result[0, i + 1] = -shift[i] * result[0, 0] + result[i + 1, 0] = result[0, i + 1] + for i in range(3): + for j in range(3): + result[i + 1, j + 1] = inverse_spatial_metric[i, j] + ( + shift[i] * shift[j] * result[0, 0] + ) + return result + + +def four_velocity(spatial_velocity, shift, lorentz_factor, lapse): + result = np.zeros(4) + result[0] = lorentz_factor / lapse + result[1:] = lorentz_factor * (spatial_velocity - shift / lapse) + return result + + +def comoving_magnetic_field( + spatial_velocity, + magnetic_field, + magnetic_field_dot_spatial_velocity, + lorentz_factor, + shift, + lapse, +): + result = np.zeros(4) + result[0] = magnetic_field_dot_spatial_velocity / lapse + result[1:] = ( + magnetic_field + + lapse + * result[0] + * (lorentz_factor * (spatial_velocity - shift / lapse)) + ) / lorentz_factor + return result + + +def energy_density( + rest_mass_density, + specific_enthalpy, + pressure, + lorentz_factor, + magnetic_field_dot_spatial_velocity, + comoving_magnetic_field_squared, +): + return ( + rest_mass_density * specific_enthalpy * lorentz_factor**2 + - pressure + + comoving_magnetic_field_squared * (lorentz_factor**2 - 0.5) + - (lorentz_factor * magnetic_field_dot_spatial_velocity) ** 2 + ) + + +def momentum_density( + rest_mass_density, + specific_enthalpy, + spatial_velocity, + lorentz_factor, + magnetic_field, + magnetic_field_dot_spatial_velocity, + comoving_magnetic_field_squared, +): + return ( + rest_mass_density + * specific_enthalpy + * lorentz_factor**2 + * spatial_velocity + + comoving_magnetic_field_squared + * lorentz_factor**2 + * spatial_velocity + - magnetic_field_dot_spatial_velocity * magnetic_field + - magnetic_field_dot_spatial_velocity**2 + * lorentz_factor**2 + * spatial_velocity + ) + + +def stress_trace( + rest_mass_density, + specific_enthalpy, + pressure, + spatial_velocity_squared, + lorentz_factor, + magnetic_field_dot_spatial_velocity, + comoving_magnetic_field_squared, +): + return ( + 3.0 * pressure + + rest_mass_density * specific_enthalpy * (lorentz_factor**2 - 1.0) + + comoving_magnetic_field_squared + * (lorentz_factor**2 * spatial_velocity_squared + 0.5) + - magnetic_field_dot_spatial_velocity**2 + * (lorentz_factor**2 * spatial_velocity_squared + 1.0) + ) + + +def stress_energy_tensor( + rest_mass_density, + specific_internal_energy, + pressure, + lorentz_factor, + lapse, + comoving_magnetic_field_magnitude, + spatial_velocity, + shift, + magnetic_field, + spatial_metric, + inverse_spatial_metric, +): + inverse_spacetime_metric_v = inverse_spacetime_metric( + lapse, shift, inverse_spatial_metric + ) + + magnetic_field_dot_spatial_velocity = np.dot( + np.dot(spatial_metric, magnetic_field), spatial_velocity + ) + comoving_magnetic_field_v = comoving_magnetic_field( + spatial_velocity, + magnetic_field, + magnetic_field_dot_spatial_velocity, + lorentz_factor, + shift, + lapse, + ) + four_velocity_v = four_velocity( + spatial_velocity, shift, lorentz_factor, lapse + ) + rho_h_star = ( + (rest_mass_density + rest_mass_density * specific_internal_energy) + + pressure + + comoving_magnetic_field_magnitude**2 + ) + p_star = pressure + 0.5 * (comoving_magnetic_field_magnitude**2) + + result = np.zeros([4, 4]) + result[0, 0] = ( + (rho_h_star * (four_velocity_v[0] ** 2)) + + (p_star * inverse_spacetime_metric_v[0, 0]) + - (comoving_magnetic_field_v[0] ** 2) + ) + result[1, 0] = ( + (rho_h_star * four_velocity_v[0] * four_velocity_v[1]) + + (p_star * inverse_spacetime_metric_v[0, 1]) + - (comoving_magnetic_field_v[0] * comoving_magnetic_field_v[1]) + ) + result[1, 1] = ( + (rho_h_star * four_velocity_v[1] ** 2) + + (p_star * inverse_spacetime_metric_v[1, 1]) + - (comoving_magnetic_field_v[1] ** 2) + ) + result[0, 1] = result[1, 0] + result[2, 0] = ( + (rho_h_star * four_velocity_v[2] * four_velocity_v[0]) + + (p_star * inverse_spacetime_metric_v[2, 0]) + - (comoving_magnetic_field_v[2] * comoving_magnetic_field_v[0]) + ) + result[0, 2] = result[2, 0] + result[2, 1] = ( + (rho_h_star * four_velocity_v[2] * four_velocity_v[1]) + + (p_star * inverse_spacetime_metric_v[2, 1]) + - (comoving_magnetic_field_v[2] * comoving_magnetic_field_v[1]) + ) + result[1, 2] = result[2, 1] + result[2, 2] = ( + (rho_h_star * four_velocity_v[2] ** 2) + + (p_star * inverse_spacetime_metric_v[2, 2]) + - (comoving_magnetic_field_v[2] ** 2) + ) + result[3, 0] = ( + (rho_h_star * four_velocity_v[3] * four_velocity_v[0]) + + (p_star * inverse_spacetime_metric_v[3, 0]) + - (comoving_magnetic_field_v[3] * comoving_magnetic_field_v[0]) + ) + result[0, 3] = result[3, 0] + result[3, 1] = ( + (rho_h_star * four_velocity_v[3] * four_velocity_v[1]) + + (p_star * inverse_spacetime_metric_v[3, 1]) + - (comoving_magnetic_field_v[3] * comoving_magnetic_field_v[1]) + ) + result[1, 3] = result[3, 1] + result[3, 2] = ( + (rho_h_star * four_velocity_v[3] * four_velocity_v[2]) + + (p_star * inverse_spacetime_metric_v[3, 2]) + - (comoving_magnetic_field_v[3] * comoving_magnetic_field_v[2]) + ) + result[2, 3] = result[3, 2] + result[3, 3] = ( + (rho_h_star * four_velocity_v[3] ** 2) + + (p_star * inverse_spacetime_metric_v[3, 3]) + - (comoving_magnetic_field_v[3] ** 2) + ) + return result diff --git a/tests/Unit/PointwiseFunctions/Hydro/TestFunctions.py b/tests/Unit/PointwiseFunctions/Hydro/TestFunctions.py index 9256f17323fc..ecd2270756a7 100644 --- a/tests/Unit/PointwiseFunctions/Hydro/TestFunctions.py +++ b/tests/Unit/PointwiseFunctions/Hydro/TestFunctions.py @@ -112,62 +112,3 @@ def relativistic_specific_enthalpy( # End functions for testing SpecificEnthalpy.cpp - - -def energy_density( - rest_mass_density, - specific_enthalpy, - pressure, - lorentz_factor, - magnetic_field_dot_spatial_velocity, - comoving_magnetic_field_squared, -): - return ( - rest_mass_density * specific_enthalpy * lorentz_factor**2 - - pressure - + comoving_magnetic_field_squared * (lorentz_factor**2 - 0.5) - - (lorentz_factor * magnetic_field_dot_spatial_velocity) ** 2 - ) - - -def momentum_density( - rest_mass_density, - specific_enthalpy, - spatial_velocity, - lorentz_factor, - magnetic_field, - magnetic_field_dot_spatial_velocity, - comoving_magnetic_field_squared, -): - return ( - rest_mass_density - * specific_enthalpy - * lorentz_factor**2 - * spatial_velocity - + comoving_magnetic_field_squared - * lorentz_factor**2 - * spatial_velocity - - magnetic_field_dot_spatial_velocity * magnetic_field - - magnetic_field_dot_spatial_velocity**2 - * lorentz_factor**2 - * spatial_velocity - ) - - -def stress_trace( - rest_mass_density, - specific_enthalpy, - pressure, - spatial_velocity_squared, - lorentz_factor, - magnetic_field_dot_spatial_velocity, - comoving_magnetic_field_squared, -): - return ( - 3.0 * pressure - + rest_mass_density * specific_enthalpy * (lorentz_factor**2 - 1.0) - + comoving_magnetic_field_squared - * (lorentz_factor**2 * spatial_velocity_squared + 0.5) - - magnetic_field_dot_spatial_velocity**2 - * (lorentz_factor**2 * spatial_velocity_squared + 1.0) - ) diff --git a/tests/Unit/PointwiseFunctions/Hydro/Test_StressEnergy.cpp b/tests/Unit/PointwiseFunctions/Hydro/Test_StressEnergy.cpp index ae2769eff5da..b756627b3aeb 100644 --- a/tests/Unit/PointwiseFunctions/Hydro/Test_StressEnergy.cpp +++ b/tests/Unit/PointwiseFunctions/Hydro/Test_StressEnergy.cpp @@ -17,17 +17,20 @@ namespace hydro { SPECTRE_TEST_CASE("Unit.PointwiseFunctions.Hydro.StressEnergy", "[Unit][Hydro]") { pypp::SetupLocalPythonEnvironment local_python_env( - "PointwiseFunctions/Hydro/"); - const DataVector used_for_size(5); + "PointwiseFunctions/Hydro/Python/"); + const DataVector used_for_size(1); pypp::check_with_random_values<1>(&energy_density, - "TestFunctions", {"energy_density"}, + "Test_StressEnergy", {"energy_density"}, {{{0.0, 1.0}}}, used_for_size); pypp::check_with_random_values<1>(&momentum_density, - "TestFunctions", {"momentum_density"}, + "Test_StressEnergy", {"momentum_density"}, {{{0.0, 1.0}}}, used_for_size); - pypp::check_with_random_values<1>(&stress_trace, "TestFunctions", - {"stress_trace"}, {{{0.0, 1.0}}}, - used_for_size); + pypp::check_with_random_values<1>(&stress_trace, + "Test_StressEnergy", {"stress_trace"}, + {{{0.0, 1.0}}}, used_for_size); + pypp::check_with_random_values<1>( + &stress_energy_tensor, "Test_StressEnergy", + {"stress_energy_tensor"}, {{{0.0, 1.0}}}, used_for_size); } } // namespace hydro