diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt b/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt index 19fa7221b400..35ea0c1e5727 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt @@ -12,6 +12,7 @@ spectre_target_sources( PunctureField.cpp PunctureFieldOrder0.cpp PunctureFieldOrder1.cpp + SelfForce.cpp ) spectre_target_headers( @@ -22,6 +23,7 @@ spectre_target_headers( SingletonChare.hpp Tags.hpp PunctureField.hpp + SelfForce.hpp Worldtube.hpp ) diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.cpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.cpp new file mode 100644 index 000000000000..a0930ff8c6ea --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.cpp @@ -0,0 +1,65 @@ + +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp" + +#include + +#include "DataStructures/Tensor/Tensor.hpp" +#include "Utilities/Gsl.hpp" + +namespace CurvedScalarWave::Worldtube { + +template +void self_force_acceleration( + gsl::not_null*> self_force_acc, + const Scalar& dt_psi_monopole, + const tnsr::i& psi_dipole, + const tnsr::I& particle_velocity, const double particle_charge, + const double particle_mass, const tnsr::AA& inverse_metric, + const Scalar& dilation_factor) { + const double factor = + particle_charge / particle_mass / square(get(dilation_factor)); + for (size_t i = 0; i < Dim; ++i) { + self_force_acc->get(i) = + (inverse_metric.get(i + 1, 0) - + particle_velocity.get(i) * inverse_metric.get(0, 0)) * + get(dt_psi_monopole) * factor; + for (size_t j = 0; j < Dim; ++j) { + self_force_acc->get(i) += + (inverse_metric.get(i + 1, j + 1) - + particle_velocity.get(i) * inverse_metric.get(0, j + 1)) * + psi_dipole.get(j) * factor; + } + } +} + +template +tnsr::I self_force_acceleration( + const Scalar& dt_psi_monopole, + const tnsr::i& psi_dipole, + const tnsr::I& particle_velocity, const double particle_charge, + const double particle_mass, const tnsr::AA& inverse_metric, + const Scalar& dilation_factor) { + tnsr::I self_force_acc{}; + self_force_acceleration(make_not_null(&self_force_acc), dt_psi_monopole, + psi_dipole, particle_velocity, particle_charge, + particle_mass, inverse_metric, dilation_factor); + return self_force_acc; +} + +// Instantiations +template void self_force_acceleration( + gsl::not_null*> self_force_acc, + const Scalar& dt_psi_monopole, const tnsr::i& psi_dipole, + const tnsr::I& particle_velocity, const double particle_charge, + const double particle_mass, const tnsr::AA& inverse_metric, + const Scalar& dilation_factor); + +template tnsr::I self_force_acceleration( + const Scalar& dt_psi_monopole, const tnsr::i& psi_dipole, + const tnsr::I& particle_velocity, const double particle_charge, + const double particle_mass, const tnsr::AA& inverse_metric, + const Scalar& dilation_factor); +} // namespace CurvedScalarWave::Worldtube diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp new file mode 100644 index 000000000000..e067cc0a1fd7 --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp @@ -0,0 +1,51 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/Tensor/Tensor.hpp" +#include "Utilities/Gsl.hpp" + +namespace CurvedScalarWave::Worldtube { + +/// @{ +/*! + * \brief Computes the coordinate acceleration due to the scalar self-force onto + * the charge. + * + * \details It is given by + * + * \begin{equation} + * (u^0)^2 \ddot{x}^i_p = \frac{q}{\mu}(g^{i + * \alpha} - \dot{x}^i_p g^{0 \alpha} ) \partial_\alpha \Psi^R + * \end{equation} + * + * where $\dot{x}^i_p$ is the position of the scalar charge, $\Psi^R$ is the + * regular field, $q$ is the particle's charge, $\mu$ is the particle's mass, + * $u^\alpha$ is the four-velocity and $g^{\alpha \beta}$ is the inverse + * spacetime metric in the inertial frame, evaluated at the position of the + * particle. An overdot denotes a derivative with respect to coordinate time. + * Greek indices are spacetime indices and Latin indices are purely spatial. + * Note that the coordinate geodesic acceleration is NOT included. + */ +template +void self_force_acceleration( + gsl::not_null*> self_force_acc, + const Scalar& dt_psi_monopole, + const tnsr::i& psi_dipole, + const tnsr::I& particle_velocity, double particle_charge, + double particle_mass, const tnsr::AA& inverse_metric, + const Scalar& dilation_factor); + +template +tnsr::I self_force_acceleration( + const Scalar& dt_psi_monopole, + const tnsr::i& psi_dipole, + const tnsr::I& particle_velocity, double particle_charge, + double particle_mass, const tnsr::AA& inverse_metric, + const Scalar& dilation_factor); + +/// @} +} // namespace CurvedScalarWave::Worldtube diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.cpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.cpp index b5fbc29f3c6e..294997bece8a 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.cpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.cpp @@ -7,6 +7,7 @@ #include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Variables.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp" #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp" #include "Utilities/Gsl.hpp" @@ -17,13 +18,28 @@ void UpdateAcceleration::apply( ::Tags::dt>>>*> dt_evolved_vars, const std::array, 2>& pos_vel, - const tnsr::I& geodesic_acc) { + const tuples::TaggedTuple, + gr::Tags::InverseSpacetimeMetric, + Tags::TimeDilationFactor>& background, + const tnsr::I& geodesic_acc, + const Scalar& dt_psi_monopole, + const tnsr::i& psi_dipole, + const double charge, const double mass, const bool apply_self_force) { + tnsr::I self_force_acc(0.); const auto& particle_velocity = pos_vel.at(1); + if (apply_self_force) { + const auto& inverse_metric = + get>(background); + const auto& dilation_factor = get(background); + self_force_acceleration(make_not_null(&self_force_acc), dt_psi_monopole, + psi_dipole, particle_velocity, charge, mass, + inverse_metric, dilation_factor); + } for (size_t i = 0; i < Dim; ++i) { get<::Tags::dt>>(*dt_evolved_vars).get(i)[0] = particle_velocity.get(i); get<::Tags::dt>>(*dt_evolved_vars).get(i)[0] = - geodesic_acc.get(i); + geodesic_acc.get(i) + self_force_acc.get(i); } } } // namespace CurvedScalarWave::Worldtube diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.hpp index 8cfe89099b2c..3861a4fd526b 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.hpp @@ -10,15 +10,21 @@ #include "DataStructures/Variables.hpp" #include "DataStructures/VariablesTag.hpp" #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp" +#include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp" #include "Parallel/AlgorithmExecution.hpp" #include "Parallel/GlobalCache.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "Utilities/Gsl.hpp" +#include "Utilities/TaggedTuple.hpp" namespace CurvedScalarWave::Worldtube { /*! - * \brief Computes the geodesic acceleration of the particle, see - * `Tags::GeodesicAccelerationCompute`. This mutator is run on the worldtube + * \brief Computes the final acceleration of the particle at this time step. + * \details If `apply_self_force` is `false`, the acceleration will simply be + * geodesic, see `gr::geodesic_acceleration`. Otherwise, the acceleration due + * to the scalar self-force is additionally applied to it, see + * `self_force_acceleration`. This mutator is run on the worldtube * singleton chare. */ struct UpdateAcceleration { @@ -27,14 +33,26 @@ struct UpdateAcceleration { tmpl::list, Tags::EvolvedVelocity>>; using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>; using return_tags = tmpl::list; - using argument_tags = tmpl::list, - Tags::GeodesicAcceleration>; + using argument_tags = tmpl::list< + Tags::ParticlePositionVelocity, Tags::BackgroundQuantities, + Tags::GeodesicAcceleration, + Stf::Tags::StfTensor<::Tags::dt, 0, Dim, + Frame::Inertial>, + Stf::Tags::StfTensor, + Tags::Charge, Tags::Mass, Tags::ApplySelfForce>; static void apply( gsl::not_null< Variables>, ::Tags::dt>>>*> dt_evolved_vars, const std::array, 2>& pos_vel, - const tnsr::I& geodesic_acc); + const tuples::TaggedTuple, + gr::Tags::InverseSpacetimeMetric, + Tags::TimeDilationFactor>& background, + const tnsr::I& geodesic_acc, + const Scalar& dt_psi_monopole, + const tnsr::i& psi_dipole, double charge, + double mass, bool apply_self_force); }; + } // namespace CurvedScalarWave::Worldtube diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt index 12a7fab87c08..64cecd9b5328 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt @@ -5,6 +5,7 @@ set(LIBRARY "Test_ScalarWaveWorldtube") set(LIBRARY_SOURCES Test_PunctureField.cpp + Test_SelfForce.cpp Test_Tags.cpp ElementActions/Test_SendToWorldtube.cpp ElementActions/Test_ReceiveWorldtubeData.cpp diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.py b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.py new file mode 100644 index 000000000000..93844d9c61b3 --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.py @@ -0,0 +1,23 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +import numpy as np + + +def self_force_acceleration( + dt_psi_monopole, + psi_dipole, + vel, + charge, + mass, + inverse_metric, + dilation, +): + # Prepend extra value so dimensions work out for einsum. + # The 0th component does not affect the final result + four_vel = np.concatenate((np.empty(1), vel), axis=0) + d_psi = np.concatenate(([dt_psi_monopole], psi_dipole), axis=0) + self_force_acc = np.einsum("ij,j", inverse_metric, d_psi) + self_force_acc -= np.einsum("i,j,j", four_vel, inverse_metric[0], d_psi) + self_force_acc *= charge / mass / dilation**2 + return self_force_acc[1:] diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateAcceleration.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateAcceleration.cpp index d7be09a5d8d4..382f615b38f0 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateAcceleration.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateAcceleration.cpp @@ -12,6 +12,7 @@ #include "DataStructures/DataBox/PrefixHelpers.hpp" #include "DataStructures/Variables.hpp" #include "DataStructures/VariablesTag.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp" #include "Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.hpp" #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp" #include "Framework/TestCreation.hpp" @@ -34,28 +35,61 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.Worldtube.UpdateAcceleration", MAKE_GENERATOR(gen); std::uniform_real_distribution<> dist(-1., 1.); const DataVector used_for_size(1); - const auto pos = make_with_random_values>( - make_not_null(&gen), dist, used_for_size); - const auto vel = make_with_random_values>( - make_not_null(&gen), dist, used_for_size); auto dt_evolved_vars = make_with_random_values( make_not_null(&gen), dist, used_for_size); - const auto acceleration = - make_with_random_values>(make_not_null(&gen), dist, 1); - auto box = db::create< - db::AddSimpleTags, dt_variables_tag, - Tags::GeodesicAcceleration>>( - std::array, 2>{{pos, vel}}, std::move(dt_evolved_vars), - acceleration); + const auto geodesic_acc = make_with_random_values>( + make_not_null(&gen), dist, 1); + const auto vel = make_with_random_values>( + make_not_null(&gen), dist, 1); + const auto pos = make_with_random_values>( + make_not_null(&gen), dist, 1); + const std::array, 2> pos_vel{{pos, vel}}; + const auto metric = make_with_random_values>( + make_not_null(&gen), dist, 1); + const auto inverse_metric = make_with_random_values>( + make_not_null(&gen), dist, 1); + const auto dilation = + make_with_random_values>(make_not_null(&gen), dist, 1); + const Tags::BackgroundQuantities::type background_quantities{ + metric, inverse_metric, dilation}; + const auto dt_psi_monopole = + make_with_random_values>(make_not_null(&gen), dist, 1); + const auto psi_dipole = make_with_random_values>( + make_not_null(&gen), dist, 1); + const bool apply_self_force = false; + auto box = db::create, + Tags::BackgroundQuantities, Tags::GeodesicAcceleration, + Stf::Tags::StfTensor<::Tags::dt, 0, Dim, + Frame::Inertial>, + Stf::Tags::StfTensor, + Tags::Charge, Tags::Mass, Tags::ApplySelfForce>>( + std::move(dt_evolved_vars), pos_vel, background_quantities, geodesic_acc, + dt_psi_monopole, psi_dipole, 1., 1., apply_self_force); db::mutate_apply(make_not_null(&box)); + const auto& dt_vars = db::get(box); + for (size_t i = 0; i < Dim; ++i) { + CHECK(get<::Tags::dt>>(dt_vars).get(i)[0] == + vel.get(i)); + CHECK(get<::Tags::dt>>(dt_vars).get(i)[0] == + geodesic_acc.get(i)); + } + + db::mutate( + [](const gsl::not_null apply_self_force_arg) { + *apply_self_force_arg = true; + }, + make_not_null(&box)); - const auto& dt_vars_after_mutate = db::get(box); + db::mutate_apply(make_not_null(&box)); + const auto self_force_acc = self_force_acceleration( + dt_psi_monopole, psi_dipole, vel, 1., 1., inverse_metric, dilation); for (size_t i = 0; i < Dim; ++i) { - CHECK(get<::Tags::dt>>(dt_vars_after_mutate) - .get(i)[0] == vel.get(i)); - CHECK(get<::Tags::dt>>(dt_vars_after_mutate) - .get(i)[0] == acceleration.get(i)); + CHECK(get<::Tags::dt>>(dt_vars).get(i)[0] == + vel.get(i)); + CHECK(get<::Tags::dt>>(dt_vars).get(i)[0] == + geodesic_acc.get(i) + self_force_acc.get(i)); } } } // namespace diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_SelfForce.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_SelfForce.cpp new file mode 100644 index 000000000000..b613e6e3cbeb --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_SelfForce.cpp @@ -0,0 +1,35 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include + +#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp" +#include "Framework/CheckWithRandomValues.hpp" +#include "Framework/SetupLocalPythonEnvironment.hpp" +#include "Framework/TestCreation.hpp" +#include "Framework/TestHelpers.hpp" +#include "Utilities/TMPL.hpp" + +namespace CurvedScalarWave::Worldtube { +namespace { + +void test_self_force_acceleration() { + pypp::SetupLocalPythonEnvironment local_python_env{ + "Evolution/Systems/CurvedScalarWave/Worldtube"}; + pypp::check_with_random_values<1>( + static_cast (*)( + const Scalar&, const tnsr::i&, + const tnsr::I&, const double, const double, + const tnsr::AA&, const Scalar&)>( + self_force_acceleration<3>), + "SelfForce", "self_force_acceleration", {{{-2.0, 2.0}}}, 1); +} + +SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.Worldtube.SelfForce", + "[Unit][Evolution]") { + test_self_force_acceleration(); +} +} // namespace +} // namespace CurvedScalarWave::Worldtube