From d9c0fe4ff6adfddc226c426d0398d95d03ece204 Mon Sep 17 00:00:00 2001 From: Nikolas Date: Sat, 23 Nov 2024 11:03:59 +0100 Subject: [PATCH] Add turn on function to self force acceleration --- .../SingletonActions/UpdateAcceleration.cpp | 16 ++++++---- .../SingletonActions/UpdateAcceleration.hpp | 10 +++++-- .../ElementActions/Test_Iterations.cpp | 13 +++++++-- .../ElementActions/Test_SendToWorldtube.cpp | 12 ++++++-- .../Test_UpdateAcceleration.cpp | 29 +++++++++++++++---- 5 files changed, 62 insertions(+), 18 deletions(-) diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.cpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.cpp index dac0515a86e4b..04fcf1a0a74b0 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.cpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.cpp @@ -25,25 +25,31 @@ void UpdateAcceleration::apply( gr::Tags::TraceSpacetimeChristoffelSecondKind, Tags::TimeDilationFactor>& background, const tnsr::I& geodesic_acc, - const Scalar& dt_psi_monopole, + const Scalar& psi_monopole, const Scalar& dt_psi_monopole, const tnsr::i& psi_dipole, const double charge, const std::optional mass, - const size_t max_iterations) { + const size_t max_iterations, const double time, + const std::optional turn_on_time, + const std::optional turn_on_interval) { tnsr::I self_force_acc(0.); const auto& particle_velocity = pos_vel.at(1); - if (max_iterations > 0) { + double roll_on = 0.; + if (max_iterations > 0 and time > turn_on_time.value()) { + roll_on = + turn_on_function(time - turn_on_time.value(), turn_on_interval.value()); const auto& inverse_metric = get>(background); const auto& dilation_factor = get(background); + const double evolved_mass = mass.value() - charge * get(psi_monopole); self_force_acceleration(make_not_null(&self_force_acc), dt_psi_monopole, - psi_dipole, particle_velocity, charge, mass.value(), + psi_dipole, particle_velocity, charge, evolved_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) + self_force_acc.get(i); + geodesic_acc.get(i) + roll_on * 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 614f309368c29..42a7ef1f6683c 100644 --- a/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.hpp @@ -36,10 +36,12 @@ struct UpdateAcceleration { using argument_tags = tmpl::list< Tags::ParticlePositionVelocity, Tags::BackgroundQuantities, Tags::GeodesicAcceleration, + Stf::Tags::StfTensor, Stf::Tags::StfTensor<::Tags::dt, 0, Dim, Frame::Inertial>, Stf::Tags::StfTensor, - Tags::Charge, Tags::Mass, Tags::MaxIterations>; + Tags::Charge, Tags::Mass, Tags::MaxIterations, ::Tags::Time, + Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>; static void apply( gsl::not_null< Variables>, @@ -53,9 +55,11 @@ struct UpdateAcceleration { gr::Tags::TraceSpacetimeChristoffelSecondKind, Tags::TimeDilationFactor>& background, const tnsr::I& geodesic_acc, - const Scalar& dt_psi_monopole, + const Scalar& psi_monopole, const Scalar& dt_psi_monopole, const tnsr::i& psi_dipole, double charge, - std::optional mass, size_t max_iterations); + std::optional mass, size_t max_iterations, double time, + std::optional turn_on_time, + std::optional turn_on_interval); }; } // namespace CurvedScalarWave::Worldtube diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp index d462f4e03486a..96bdc20455650 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp @@ -121,7 +121,8 @@ struct MockMetavariables { domain::Tags::Domain, CurvedScalarWave::Tags::BackgroundSpacetime, Tags::ExcisionSphere, Tags::WorldtubeRadius, Tags::ExpansionOrder, - Tags::MaxIterations, Tags::Charge, Tags::Mass>; + Tags::MaxIterations, Tags::Charge, Tags::Mass, ::Tags::Time, + Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>; }; void test_iterations(const size_t max_iterations) { @@ -163,11 +164,14 @@ void test_iterations(const size_t max_iterations) { const auto& initial_refinements = shell.initial_refinement_levels(); const auto& initial_extents = shell.initial_extents(); + // unused but the tag is needed to compile + const double time = std::numeric_limits::signaling_NaN(); tuples::TaggedTuple< domain::Tags::Domain, CurvedScalarWave::Tags::BackgroundSpacetime, Tags::ExcisionSphere, Tags::WorldtubeRadius, Tags::ExpansionOrder, - Tags::MaxIterations, Tags::Charge, Tags::Mass> + Tags::MaxIterations, Tags::Charge, Tags::Mass, ::Tags::Time, + Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval> tuple_of_opts{shell.create_domain(), kerr_schild, excision_sphere, @@ -175,7 +179,10 @@ void test_iterations(const size_t max_iterations) { expansion_order, max_iterations, charge, - std::make_optional(mass)}; + std::make_optional(mass), + time, + std::nullopt, + std::nullopt}; ActionTesting::MockRuntimeSystem runner{std::move(tuple_of_opts)}; const auto element_ids = initial_element_ids(initial_refinements); const auto& blocks = shell_domain.blocks(); diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp index 28dc78bb6a6ff..a25bf2112a2d3 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp @@ -123,7 +123,8 @@ struct MockMetavariables { using const_global_cache_tags = tmpl::list, Tags::ExcisionSphere, Tags::WorldtubeRadius, Tags::ExpansionOrder, - Tags::MaxIterations, Tags::Charge, Tags::Mass>; + Tags::MaxIterations, Tags::Charge, Tags::Mass, ::Tags::Time, + Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>; }; // This test checks that `SendToWorldtube` integrates the regular field on the @@ -142,6 +143,8 @@ SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.SendToWorldtube", "[Unit]") { const size_t initial_extent = 10; const size_t face_size = initial_extent * initial_extent; const auto quadrature = Spectral::Quadrature::GaussLobatto; + // unused but the tag is needed to compile + const double time = std::numeric_limits::signaling_NaN(); // we create several differently refined shells so a different number of // elements sends data for (const auto& [expansion_order, initial_refinement, worldtube_radius] : @@ -166,13 +169,18 @@ SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.SendToWorldtube", "[Unit]") { // self force and therefore iterative scheme is turned off tuples::TaggedTuple, Tags::ExcisionSphere, Tags::WorldtubeRadius, Tags::ExpansionOrder, - Tags::MaxIterations, Tags::Charge, Tags::Mass> + Tags::MaxIterations, Tags::Charge, Tags::Mass, + ::Tags::Time, Tags::SelfForceTurnOnTime, + Tags::SelfForceTurnOnInterval> tuple_of_opts{shell.create_domain(), excision_sphere, excision_sphere.radius(), expansion_order, static_cast(0), 0.1, + std::nullopt, + time, + std::nullopt, std::nullopt}; ActionTesting::MockRuntimeSystem runner{std::move(tuple_of_opts)}; const auto element_ids = initial_element_ids(initial_refinements); 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 3b1adb7624a4e..e9db3d94fae2c 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateAcceleration.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateAcceleration.cpp @@ -21,6 +21,7 @@ #include "Helpers/DataStructures/MakeWithRandomValues.hpp" #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" #include "ParallelAlgorithms/Actions/MutateApply.hpp" +#include "Time/Tags/Time.hpp" #include "Utilities/TMPL.hpp" namespace CurvedScalarWave::Worldtube { @@ -56,20 +57,33 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.Worldtube.UpdateAcceleration", make_with_random_values>(make_not_null(&gen), dist, 1); const Tags::BackgroundQuantities::type background_quantities{ metric, inverse_metric, christoffel, trace_christoffel, dilation}; + 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 double mass = 0.1; + const double charge = 0.2; + + const double time = 105.; + const double turn_on_time = 100.; + const double turn_on_interval = 10.; const size_t max_iterations = 0; auto box = db::create, Tags::BackgroundQuantities, Tags::GeodesicAcceleration, + Stf::Tags::StfTensor, Stf::Tags::StfTensor<::Tags::dt, 0, Dim, Frame::Inertial>, Stf::Tags::StfTensor, - Tags::Charge, Tags::Mass, Tags::MaxIterations>>( + Tags::Charge, Tags::Mass, Tags::MaxIterations, ::Tags::Time, + Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>>( std::move(dt_evolved_vars), pos_vel, background_quantities, geodesic_acc, - dt_psi_monopole, psi_dipole, 1., std::make_optional(1.), max_iterations); + psi_monopole, dt_psi_monopole, psi_dipole, charge, + std::make_optional(mass), max_iterations, time, + std::make_optional(turn_on_time), std::make_optional(turn_on_interval)); db::mutate_apply(make_not_null(&box)); const auto& dt_vars = db::get(box); @@ -87,13 +101,18 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.Worldtube.UpdateAcceleration", make_not_null(&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); + + const double roll_on = + turn_on_function(time - turn_on_time, turn_on_interval); + const double evolved_mass = mass - charge * get(psi_monopole); + const auto self_force_acc = + self_force_acceleration(dt_psi_monopole, psi_dipole, vel, charge, + evolved_mass, inverse_metric, dilation); 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) + self_force_acc.get(i)); + geodesic_acc.get(i) + roll_on * self_force_acc.get(i)); } } } // namespace