diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/IteratePunctureField.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/IteratePunctureField.hpp index e83e8aa39a993..cf63d9f8dd23e 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/IteratePunctureField.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/IteratePunctureField.hpp @@ -59,15 +59,31 @@ struct IteratePunctureField { const auto iterated_puncture_field) { tnsr::I iterated_acceleration{ {self_force_data[0], self_force_data[1], self_force_data[2]}}; + const size_t face_size = + get<0>(centered_face_coordinates.value()).size(); + if (not iterated_puncture_field->has_value()) { - iterated_puncture_field->emplace( - get<0>(centered_face_coordinates.value()).size()); + iterated_puncture_field->emplace(face_size); } puncture_field(make_not_null(&iterated_puncture_field->value()), centered_face_coordinates.value(), position_velocity[0], position_velocity[1], iterated_acceleration, 1., order); + Variables, + ::Tags::deriv, Frame::Inertial>>> + acceleration_terms(face_size); + acceleration_terms_1( + make_not_null(&acceleration_terms), + centered_face_coordinates.value(), position_velocity[0], + position_velocity[1], iterated_acceleration, self_force_data[3], + self_force_data[4], self_force_data[5], self_force_data[6], + self_force_data[7], self_force_data[8], self_force_data[9], + self_force_data[10], self_force_data[11], self_force_data[12], + self_force_data[13], self_force_data[14], 1.); + iterated_puncture_field->value() += acceleration_terms; iterated_puncture_field->value() *= charge; }, make_not_null(&box)); diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.cpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.cpp index 95e3930b282d6..5cb9f8cd6f5a4 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.cpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.cpp @@ -8,6 +8,7 @@ #include #include "Evolution/Systems/CurvedScalarWave/Worldtube/Inboxes.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp" #include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp" #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp" #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp" @@ -25,20 +26,183 @@ void IterateAccelerationTerms::apply( Tags::TimeDilationFactor>& background, const tnsr::I& geodesic_acc, const Scalar& psi_monopole, const Scalar& dt_psi_monopole, - const tnsr::i& psi_dipole, double charge, - std::optional mass, double time, std::optional turn_on_time, - std::optional turn_on_interval) { - const size_t data_size = 3; + const tnsr::i& psi_dipole, + const tnsr::i& dt_psi_dipole, + const double charge, const std::optional mass, const double time, + const std::optional turn_on_time, + const std::optional turn_on_interval, const size_t iteration) { + const size_t data_size = 15; acceleration_terms->get() = DataVector(data_size, 0.); - const auto& inverse_metric = - get>(background); - const auto& dilation_factor = get(background); - const auto& vel = pos_vel.at(1); - const auto self_force_acc = - self_force_acceleration(dt_psi_monopole, psi_dipole, vel, charge, - mass.value(), inverse_metric, dilation_factor); + auto acc = geodesic_acc; + double roll_on = 0.; + if (time > turn_on_time.value()) { + const auto& particle_position = pos_vel.at(0); + const auto& particle_velocity = pos_vel.at(1); + const auto& imetric = + get>(background); + const auto& metric = + get>(background); + const auto& christoffel = + get>(background); + const auto& trace_christoffel = + get>( + background); + const auto di_imetric = + spatial_derivative_inverse_ks_metric(particle_position); + const auto dij_imetric = + second_spatial_derivative_inverse_ks_metric(particle_position); + const auto di_trace_christoffel = + spatial_derivative_ks_contracted_christoffel(particle_position); + const auto di_metric = spatial_derivative_ks_metric(metric, di_imetric); + const auto dij_metric = second_spatial_derivative_metric( + metric, di_metric, di_imetric, dij_imetric); + const auto di_christoffel = spatial_derivative_christoffel( + di_metric, dij_metric, imetric, di_imetric); + const auto& u0 = get(get(background)); + + // The evolution of the mass is a second order effect so we only include it + // starting at the 2nd iteration. + const double evolved_mass = iteration == 1 + ? mass.value() + : mass.value() - charge * get(psi_monopole); + + const auto self_force_acc = self_force_acceleration( + dt_psi_monopole, psi_dipole, particle_velocity, charge, evolved_mass, + imetric, get(background)); + + const double t_minus_turnup = time - turn_on_time.value(); + roll_on = turn_on_function(t_minus_turnup, turn_on_interval.value()); + const double dt_roll_on = + dt_turn_on_function(t_minus_turnup, turn_on_interval.value()); + const double dt2_roll_on = + dt2_turn_on_function(t_minus_turnup, turn_on_interval.value()); + for (size_t i = 0; i < Dim; ++i) { + acc.get(i) += roll_on * self_force_acc.get(i); + } + const tnsr::A u{{u0, get<0>(particle_velocity) * u0, + get<1>(particle_velocity) * u0, + get<2>(particle_velocity) * u0}}; + const tnsr::a d_psiR{{get(dt_psi_monopole), get<0>(psi_dipole), + get<1>(psi_dipole), get<2>(psi_dipole)}}; + + // Eq.(55) + double dt2_psiR = get<0>(trace_christoffel) * get<0>(d_psiR); + for (size_t i = 0; i < Dim; ++i) { + dt2_psiR += -2. * imetric.get(0, i + 1) * dt_psi_dipole.get(i) + + trace_christoffel.get(i + 1) * psi_dipole.get(i); + } + dt2_psiR /= get<0, 0>(imetric); + + tnsr::a dt_d_psiR{{dt2_psiR, get<0>(dt_psi_dipole), + get<1>(dt_psi_dipole), + get<2>(dt_psi_dipole)}}; + for (size_t i = 0; i < Dim; ++i) { + get<0>(dt_d_psiR) += particle_velocity.get(i) * dt_psi_dipole.get(i); + } + + double dt_evolved_mass = get(dt_psi_monopole); + double dt2_evolved_mass = dt2_psiR; + for (size_t i = 0; i < Dim; ++i) { + dt_evolved_mass += psi_dipole.get(i) * particle_velocity.get(i); + dt2_evolved_mass += 2. * dt_psi_dipole.get(i) * particle_velocity.get(i) + + psi_dipole.get(i) * acc.get(i); + } + dt_evolved_mass *= iteration == 1 ? 0. : -charge; + dt2_evolved_mass *= iteration == 1 ? 0. : -charge; + + const double dt_mass_factor = -dt_evolved_mass / evolved_mass; + const double dt2_mass_factor = (2. * dt_evolved_mass * dt_evolved_mass - + evolved_mass * dt2_evolved_mass) / + square(evolved_mass); + const auto dt_christoffel = tenex::evaluate( + particle_velocity(ti::I) * di_christoffel(ti::i, ti::A, ti::b, ti::c)); + + const auto dt_imetric = tenex::evaluate( + particle_velocity(ti::I) * di_imetric(ti::i, ti::A, ti::B)); + const auto dt2_imetric = tenex::evaluate( + particle_velocity(ti::I) * particle_velocity(ti::J) * + dij_imetric(ti::i, ti::j, ti::A, ti::B) + + acc(ti::I) * di_imetric(ti::i, ti::A, ti::B)); + // Eq.(51a) + const auto dt_u = tenex::evaluate( + charge / evolved_mass / u0 * imetric(ti::A, ti::B) * d_psiR(ti::b) - + christoffel(ti::A, ti::b, ti::c) * u(ti::B) * u(ti::C) / u0); + // Eq.(51b) + const auto dt2_u = tenex::evaluate( + charge / evolved_mass / u0 * + (dt_imetric(ti::A, ti::B) * d_psiR(ti::b) + + imetric(ti::A, ti::B) * dt_d_psiR(ti::b) + + dt_mass_factor * imetric(ti::A, ti::B) * d_psiR(ti::b)) - + (dt_christoffel(ti::A, ti::b, ti::c) * u(ti::B) * u(ti::C) + + 2. * christoffel(ti::A, ti::b, ti::c) * dt_u(ti::B) * u(ti::C) + + get<0>(dt_u) * dt_u(ti::A)) / + u0); + const auto dt_u_rollon = tenex::evaluate(roll_on * dt_u(ti::A)); + const auto dt2_u_rollon = tenex::evaluate(dt_roll_on * dt_u(ti::A) + + roll_on * dt2_u(ti::A)); + + // Eq.(56) + tnsr::i d_dt2_psiR{0.}; + for (size_t i = 0; i < Dim; ++i) { + d_dt2_psiR.get(i) += + -di_imetric.get(i, 0, 0) * dt2_psiR + + di_trace_christoffel.get(i, 0) * get(dt_psi_monopole) + + trace_christoffel.get(0) * dt_psi_dipole.get(i); + for (size_t j = 0; j < Dim; ++j) { + d_dt2_psiR.get(i) += + -2. * di_imetric.get(i, 0, j + 1) * dt_psi_dipole.get(j) + + di_trace_christoffel.get(i, j + 1) * psi_dipole.get(j); + } + d_dt2_psiR.get(i) /= get<0, 0>(imetric); + } + + // Eq.(57) + double dt3_psiR = trace_christoffel.get(0) * dt2_psiR; + for (size_t i = 0; i < Dim; ++i) { + dt3_psiR += -2. * imetric.get(0, i + 1) * d_dt2_psiR.get(i) + + trace_christoffel.get(i + 1) * dt_psi_dipole.get(i); + } + dt3_psiR /= get<0, 0>(imetric); + tnsr::a dt2_d_psiR; + dt2_d_psiR.get(0) = dt3_psiR; + for (size_t i = 0; i < Dim; ++i) { + dt2_d_psiR.get(0) += 2. * d_dt2_psiR.get(i) * particle_velocity.get(i) + + dt_psi_dipole.get(i) * acc.get(i); + dt2_d_psiR.get(i + 1) = d_dt2_psiR.get(i); + } + auto f = self_force_per_mass(d_psiR, u, charge, evolved_mass, imetric); + auto dt_f = + dt_self_force_per_mass(d_psiR, dt_d_psiR, u, dt_u_rollon, charge, + evolved_mass, imetric, dt_imetric); + auto dt2_f = dt2_self_force_per_mass( + d_psiR, dt_d_psiR, dt2_d_psiR, u, dt_u_rollon, dt2_u_rollon, charge, + evolved_mass, imetric, dt_imetric, dt2_imetric); + for (size_t i = 0; i < 4; ++i) { + dt_f.get(i) += dt_mass_factor * f.get(i); + dt2_f.get(i) += + 2. * dt_mass_factor * dt_f.get(i) + dt2_mass_factor * f.get(i); + } + const auto f_roll_on = tenex::evaluate(roll_on * f(ti::A)); + const auto dt_f_roll_on = + tenex::evaluate(dt_roll_on * f(ti::A) + roll_on * dt_f(ti::A)); + const auto dt2_f_roll_on = tenex::evaluate( + dt2_roll_on * f(ti::A) + 2. * dt_roll_on * dt_f(ti::A) + + roll_on * dt2_f(ti::A)); + + const auto cov_f = + Du_self_force_per_mass(f_roll_on, dt_f_roll_on, u, christoffel); + const auto dt_cov_f = + dt_Du_self_force_per_mass(f_roll_on, dt_f_roll_on, dt2_f_roll_on, u, + dt_u_rollon, christoffel, dt_christoffel); + for (size_t i = 0; i < Dim; ++i) { + get(*acceleration_terms)[i + 3] = f_roll_on.get(i); + get(*acceleration_terms)[i + 6] = dt_f_roll_on.get(i); + get(*acceleration_terms)[i + 9] = cov_f.get(i); + get(*acceleration_terms)[i + 12] = dt_cov_f.get(i); + } + } for (size_t i = 0; i < Dim; ++i) { - get(*acceleration_terms)[i] = geodesic_acc.get(i) + self_force_acc.get(i); + get(*acceleration_terms)[i] = acc.get(i); } } } // namespace CurvedScalarWave::Worldtube diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.hpp index 7aca1e6ebd473..f84117104a362 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.hpp @@ -17,7 +17,11 @@ namespace CurvedScalarWave::Worldtube { /*! * \brief Computes the next iteration of the acceleration due to scalar self - * force from the current iteration of the regular field. + * force from the current iteration of the regular field, as well as the + * quantities required to compute the acceleration terms of the puncture field. + * + * \details Analytic expressions for the computed terms are given in Section V.B + * of \cite Wittek:2024gxn. */ struct IterateAccelerationTerms { static constexpr size_t Dim = 3; @@ -30,10 +34,12 @@ struct IterateAccelerationTerms { Stf::Tags::StfTensor<::Tags::dt, 0, Dim, Frame::Inertial>, Stf::Tags::StfTensor, + Stf::Tags::StfTensor<::Tags::dt, 1, Dim, + Frame::Inertial>, Tags::Charge, Tags::Mass, ::Tags::Time, Tags::SelfForceTurnOnTime, - Tags::SelfForceTurnOnInterval>; + Tags::SelfForceTurnOnInterval, Tags::CurrentIteration>; static void apply( - const gsl::not_null*> acceleration_terms, + gsl::not_null*> acceleration_terms, const std::array, 2>& pos_vel, const tuples::TaggedTuple< gr::Tags::SpacetimeMetric, @@ -43,10 +49,11 @@ struct IterateAccelerationTerms { Tags::TimeDilationFactor>& background, const tnsr::I& geodesic_acc, const Scalar& psi_monopole, const Scalar& dt_psi_monopole, - const tnsr::i& psi_dipole, double charge, + const tnsr::i& psi_dipole, + const tnsr::i& dt_psi_dipole, double charge, std::optional mass, double time, std::optional turn_on_time, - std::optional turn_on_interval); + std::optional turn_on_interval, size_t iteration); }; } // namespace CurvedScalarWave::Worldtube diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt index 3b87f4f0030ec..7e76d3a5ff480 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt @@ -18,6 +18,7 @@ set(LIBRARY_SOURCES SingletonActions/Test_ChangeSlabSize.cpp SingletonActions/Test_InitializeElementFacesGridCoordinates.cpp SingletonActions/Test_InitializeEvolvedVariables.cpp + SingletonActions/Test_IterateAccelerationTerms.cpp SingletonActions/Test_ObserveWorldtubeSolution.cpp SingletonActions/Test_UpdateAcceleration.cpp SingletonActions/Test_UpdateFunctionsOfTime.cpp diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.py b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.py index fceee439c6951..908a1936ee1de 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.py +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.py @@ -126,3 +126,151 @@ def dt_Du_self_force_per_mass( def turn_on_function(time, timescale): return 1.0 - np.exp(-((time / timescale) ** 4)) + + +def iterate_acceleration_terms( + vel, + psi_monopole, + dt_psi_monopole, + psi_dipole, + dt_psi_dipole, + roll_on, + dt_roll_on, + dt2_roll_on, + imetric, + christoffel, + trace_christoffel, + di_imetric, + dij_imetric, + di_trace_christoffel, + di_christoffel, + u0, + geodesic_acc, + charge, + mass, + current_iteration, +): + res = np.zeros(15) + + evolved_mass = ( + mass if current_iteration == 1 else mass - charge * psi_monopole + ) + dt_evolved_mass = dt_psi_monopole + np.einsum("i,i", vel, psi_dipole) + dt_evolved_mass = 0 if current_iteration == 1 else -charge * dt_evolved_mass + self_force_acc = self_force_acceleration( + dt_psi_monopole, psi_dipole, vel, charge, evolved_mass, imetric, u0 + ) + acc = geodesic_acc + roll_on * self_force_acc + res[:3] = acc + + d_psi = np.concatenate(([dt_psi_monopole], psi_dipole), axis=0) + u = np.concatenate(([u0], u0 * vel)) + f = self_force_per_mass(d_psi, u, charge, evolved_mass, imetric) + f_rollon = f * roll_on + res[3:6] = f_rollon[:3] + + dt_u = ( + charge / evolved_mass * np.einsum("ab,b", imetric, d_psi) + - np.einsum("abc,b,c", christoffel, u, u) + ) / u0 + dt_u_roll_on = dt_u * roll_on + dt_imetric = np.einsum("abc,a", di_imetric, vel) + dt2_psi = ( + -2.0 * np.einsum("i,i", imetric[0, 1:], dt_psi_dipole) + + trace_christoffel[0] * dt_psi_monopole + + np.einsum("i,i", trace_christoffel[1:], psi_dipole) + ) / imetric[0, 0] + dt_d_psi = np.concatenate(([dt2_psi], dt_psi_dipole)) + dt_d_psi[0] += np.einsum("i,i", vel, dt_psi_dipole) + dt_f = ( + charge + / evolved_mass + * ( + np.einsum("ab,b", dt_imetric, d_psi) + + np.einsum("a,b,b", dt_u_roll_on, u, d_psi) + + np.einsum("a,b,b", u, dt_u_roll_on, d_psi) + + np.einsum("ab,b", imetric, dt_d_psi) + + np.einsum("a,b,b", u, u, dt_d_psi) + ) + ) + dt_f -= dt_evolved_mass / evolved_mass * f + dt_f_rollon = dt_f * roll_on + f * dt_roll_on + res[6:9] = dt_f_rollon[:3] + + cov_f = Du_self_force_per_mass(f_rollon, dt_f_rollon, u, christoffel) + res[9:12] = cov_f[:3] + dt2_di_psi = ( + -2.0 * np.einsum("ij,j", di_imetric[:, 1:, 0], dt_psi_dipole) + + di_trace_christoffel[:, 0] * dt_psi_monopole + + trace_christoffel[0] * dt_psi_dipole + + np.einsum("ij,j", di_trace_christoffel[:, 1:], psi_dipole) + - di_imetric[:, 0, 0] * dt2_psi + ) / imetric[0, 0] + dt3_psi = ( + -2.0 * np.einsum("i,i", imetric[0, 1:], dt2_di_psi) + + trace_christoffel[0] * dt2_psi + + np.einsum("i,i", trace_christoffel[1:], dt_psi_dipole) + ) / imetric[0, 0] + dt2_d_psi = np.concatenate(([dt3_psi], dt2_di_psi)) + dt2_d_psi[0] += 2.0 * np.einsum("i,i", vel, dt2_di_psi) + np.einsum( + "i,i", acc, dt_psi_dipole + ) + dt_christoffel = np.einsum("iabc,i", di_christoffel, vel) + dt2_imetric = np.einsum("ijab,i,j", dij_imetric, vel, vel) + np.einsum( + "iab,i", di_imetric, acc + ) + dt2_u = ( + charge + / evolved_mass + * ( + np.einsum("ab,b", dt_imetric, d_psi) + + np.einsum("ab,b", imetric, dt_d_psi) + - dt_evolved_mass / evolved_mass * np.einsum("ab,b", imetric, d_psi) + ) + - np.einsum("abc,b,c", dt_christoffel, u, u) + - 2.0 * np.einsum("abc,b,c", christoffel, dt_u, u) + - dt_u * dt_u[0] + ) / u[0] + dt2_u_roll_on = dt_roll_on * dt_u + roll_on * dt2_u + dt2_f = dt2_self_force_per_mass( + d_psi, + dt_d_psi, + dt2_d_psi, + u, + dt_u_roll_on, + dt2_u_roll_on, + charge, + evolved_mass, + imetric, + dt_imetric, + dt2_imetric, + ) + dt2_evolved_mass = ( + dt2_psi + + 2.0 * np.einsum("i,i", vel, dt_psi_dipole) + + np.einsum("i,i", acc, psi_dipole) + ) + dt2_evolved_mass = ( + 0.0 if current_iteration == 1 else -charge * dt2_evolved_mass + ) + + dt2_f += ( + -2.0 * dt_evolved_mass / evolved_mass * dt_f + + f + * (2.0 * dt_evolved_mass**2 - dt2_evolved_mass * evolved_mass) + / evolved_mass**2 + ) + + dt2_f_roll_on = dt2_roll_on * f + 2.0 * dt_roll_on * dt_f + roll_on * dt2_f + dt_cov_f = dt_Du_self_force_per_mass( + f_rollon, + dt_f_rollon, + dt2_f_roll_on, + u, + dt_u_roll_on, + christoffel, + dt_christoffel, + ) + res[12:15] = dt_cov_f[:3] + + return res diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_IterateAccelerationTerms.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_IterateAccelerationTerms.cpp new file mode 100644 index 0000000000000..158b551bb3f30 --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_IterateAccelerationTerms.cpp @@ -0,0 +1,150 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include + +#include "Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.hpp" +#include "Framework/Pypp.hpp" +#include "Framework/SetupLocalPythonEnvironment.hpp" +#include "Framework/TestCreation.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/DataStructures/DataBox/TestHelpers.hpp" +#include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" +#include "Time/Tags/Time.hpp" +#include "Utilities/TMPL.hpp" + +namespace CurvedScalarWave::Worldtube { + +SPECTRE_TEST_CASE( + "Unit.Evolution.Systems.CSW.Worldtube.IterateAccelerationTerms", + "[Unit][Evolution]") { + const pypp::SetupLocalPythonEnvironment local_python_env{ + "Evolution/Systems/CurvedScalarWave/Worldtube"}; + MAKE_GENERATOR(gen); + static constexpr size_t Dim = 3; + const std::uniform_real_distribution dist(-1., 1.); + const std::uniform_real_distribution pos_dist(2., 10.); + const std::uniform_real_distribution vel_dist(-0.1, 0.1); + const gr::Solutions::KerrSchild kerr_schild(1.4, {{0.1, 0.2, 0.3}}, + {{0., 0., 0.}}); + const auto pos = make_with_random_values>( + make_not_null(&gen), make_not_null(&pos_dist), 1); + const auto vel = make_with_random_values>( + make_not_null(&gen), make_not_null(&vel_dist), 1); + + const auto psi_monopole = + make_with_random_values>(make_not_null(&gen), dist, 1); + 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 auto dt_psi_dipole = make_with_random_values>( + make_not_null(&gen), dist, 1); + const double charge = 0.1; + const double mass = 0.1; + const double time = 10.; + const double turn_on_time = 20.; + const double turn_on_interval = 1.; + const size_t current_iteration = 1; + auto box = db::create< + db::AddSimpleTags< + Tags::AccelerationTerms, Tags::ParticlePositionVelocity, + CurvedScalarWave::Tags::BackgroundSpacetime< + gr::Solutions::KerrSchild>, + Stf::Tags::StfTensor, + Stf::Tags::StfTensor<::Tags::dt, 0, Dim, + Frame::Inertial>, + Stf::Tags::StfTensor, + Stf::Tags::StfTensor<::Tags::dt, 1, Dim, + Frame::Inertial>, + Tags::Charge, Tags::Mass, ::Tags::Time, Tags::SelfForceTurnOnTime, + Tags::SelfForceTurnOnInterval, Tags::CurrentIteration>, + db::AddComputeTags, + Tags::GeodesicAccelerationCompute>>( + Scalar{}, std::array, 2>{pos, vel}, + kerr_schild, psi_monopole, dt_psi_monopole, psi_dipole, dt_psi_dipole, + charge, std::make_optional(mass), time, std::make_optional(turn_on_time), + std::make_optional(turn_on_interval), current_iteration); + db::mutate_apply(make_not_null(&box)); + + const auto& background = db::get>(box); + const auto& geodesic_acc = db::get>(box); + const auto& imetric = + get>(background); + const auto& metric = get>(background); + const auto& christoffel = + get>(background); + const auto& trace_christoffel = + get>( + background); + const auto di_imetric = spatial_derivative_inverse_ks_metric(pos); + const auto dij_imetric = second_spatial_derivative_inverse_ks_metric(pos); + const auto di_trace_christoffel = + spatial_derivative_ks_contracted_christoffel(pos); + const auto di_metric = spatial_derivative_ks_metric(metric, di_imetric); + const auto dij_metric = second_spatial_derivative_metric( + metric, di_metric, di_imetric, dij_imetric); + const auto di_christoffel = spatial_derivative_christoffel( + di_metric, dij_metric, imetric, di_imetric); + const auto& u0 = get(get(background)); + + double roll_on = 0.; + double dt_roll_on = 0.; + double dt2_roll_on = 0.; + + const auto expected_terms = pypp::call( + "SelfForce", "iterate_acceleration_terms", vel, psi_monopole, + dt_psi_monopole, psi_dipole, dt_psi_dipole, roll_on, dt_roll_on, + dt2_roll_on, imetric, christoffel, trace_christoffel, di_imetric, + dij_imetric, di_trace_christoffel, di_christoffel, u0, geodesic_acc, + charge, mass, current_iteration); + + CHECK_ITERABLE_APPROX(expected_terms, + get(db::get(box))); + + const double new_time = 21.123; + db::mutate<::Tags::Time>( + [&new_time](const gsl::not_null local_time) { + *local_time = new_time; + }, + make_not_null(&box)); + db::mutate_apply(make_not_null(&box)); + + const double t_minus_turnup = new_time - turn_on_time; + roll_on = turn_on_function(t_minus_turnup, turn_on_interval); + dt_roll_on = dt_turn_on_function(t_minus_turnup, turn_on_interval); + dt2_roll_on = dt2_turn_on_function(t_minus_turnup, turn_on_interval); + + const auto expected_terms_2 = pypp::call( + "SelfForce", "iterate_acceleration_terms", vel, psi_monopole, + dt_psi_monopole, psi_dipole, dt_psi_dipole, roll_on, dt_roll_on, + dt2_roll_on, imetric, christoffel, trace_christoffel, di_imetric, + dij_imetric, di_trace_christoffel, di_christoffel, u0, geodesic_acc, + charge, mass, current_iteration); + + CHECK_ITERABLE_APPROX(expected_terms_2, + get(db::get(box))); + + static constexpr size_t new_iteration = 7; + db::mutate( + [](const gsl::not_null local_current_iteration) { + *local_current_iteration = new_iteration; + }, + make_not_null(&box)); + db::mutate_apply(make_not_null(&box)); + + const auto expected_terms_3 = pypp::call( + "SelfForce", "iterate_acceleration_terms", vel, psi_monopole, + dt_psi_monopole, psi_dipole, dt_psi_dipole, roll_on, dt_roll_on, + dt2_roll_on, imetric, christoffel, trace_christoffel, di_imetric, + dij_imetric, di_trace_christoffel, di_christoffel, u0, geodesic_acc, + charge, mass, new_iteration); + CHECK_ITERABLE_APPROX(expected_terms_3, + get(db::get(box))); +} +} // namespace CurvedScalarWave::Worldtube