From 46028a9d88c4f90529ce408fd791599041e9859e Mon Sep 17 00:00:00 2001 From: Nikolas Date: Fri, 17 Feb 2023 17:56:10 +0100 Subject: [PATCH] Add worldtube boundary conditions to CSW --- .../BoundaryConditions/CMakeLists.txt | 2 + .../BoundaryConditions/Worldtube.cpp | 156 ++++++++++++ .../BoundaryConditions/Worldtube.hpp | 146 +++++++++++ .../BoundaryConditions/Test_Worldtube.cpp | 227 ++++++++++++++++++ .../BoundaryConditions/Worldtube.py | 141 +++++++++++ .../Systems/CurvedScalarWave/CMakeLists.txt | 1 + 6 files changed, 673 insertions(+) create mode 100644 src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.cpp create mode 100644 src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.hpp create mode 100644 tests/Unit/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Test_Worldtube.cpp create mode 100644 tests/Unit/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.py diff --git a/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/CMakeLists.txt b/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/CMakeLists.txt index 54847d84b3d9..5ffe82f25f82 100644 --- a/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/CMakeLists.txt +++ b/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/CMakeLists.txt @@ -7,6 +7,7 @@ spectre_target_sources( ConstraintPreservingSphericalRadiation.cpp BoundaryCondition.cpp DemandOutgoingCharSpeeds.cpp + Worldtube.cpp ) spectre_target_headers( @@ -17,4 +18,5 @@ spectre_target_headers( BoundaryCondition.hpp DemandOutgoingCharSpeeds.hpp Factory.hpp + Worldtube.hpp ) diff --git a/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.cpp b/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.cpp new file mode 100644 index 000000000000..49965029454a --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.cpp @@ -0,0 +1,156 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.hpp" + +#include +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Evolution/Systems/CurvedScalarWave/Characteristics.hpp" +#include "Evolution/Systems/CurvedScalarWave/Tags.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp" +#include "Utilities/Gsl.hpp" +namespace CurvedScalarWave::BoundaryConditions { + +template +Worldtube::Worldtube(CkMigrateMessage* const msg) + : BoundaryCondition(msg) {} + +template +std::unique_ptr +Worldtube::get_clone() const { + return std::make_unique(*this); +} + +template +void Worldtube::pup(PUP::er& p) { + BoundaryCondition::pup(p); +} + +template +// NOLINTNEXTLINE +PUP::able::PUP_ID Worldtube::my_PUP_ID = 0; + +template +std::optional Worldtube::dg_time_derivative( + gsl::not_null*> dt_psi_correction, + gsl::not_null*> dt_pi_correction, + gsl::not_null*> dt_phi_correction, + const std::optional>& face_mesh_velocity, + const tnsr::i& normal_covector, + const tnsr::I& normal_vector, + const Scalar& /*psi*/, const Scalar& /*pi*/, + const tnsr::i& phi, const Scalar& lapse, + const tnsr::I& shift, + const tnsr::II& /*inverse_spatial_metric*/, + const Scalar& gamma1, const Scalar& gamma2, + const Scalar& /*dt_psi*/, const tnsr::i& d_psi, + const tnsr::ij& d_phi, + const Variables>>& /*worldtube_vars*/) const { + auto char_speeds = + characteristic_speeds(gamma1, lapse, shift, normal_covector); + if (face_mesh_velocity.has_value()) { + const auto face_speed = dot_product(normal_covector, *face_mesh_velocity); + for (size_t i = 0; i < 4; ++i) { + char_speeds.at(i) -= get(face_speed); + } + } + *dt_psi_correction = + tenex::evaluate(normal_vector(ti::I) * (d_psi(ti::i) - phi(ti::i))); + *dt_phi_correction = tenex::evaluate( + normal_vector(ti::J) * (d_phi(ti::j, ti::i) - d_phi(ti::i, ti::j))); + + for (size_t i = 0; i < get(*dt_psi_correction).size(); ++i) { + get(*dt_psi_correction)[i] *= std::min(0., gsl::at(char_speeds[0], i)); + for (size_t j = 0; j < Dim; ++j) { + dt_phi_correction->get(j)[i] *= std::min(0., gsl::at(char_speeds[1], i)); + } + } + get(*dt_pi_correction) = get(gamma2) * get(*dt_psi_correction); + return {}; +} + +template +std::optional Worldtube::dg_ghost( + const gsl::not_null*> psi, + const gsl::not_null*> pi, + const gsl::not_null*> phi, + const gsl::not_null*> lapse, + const gsl::not_null*> shift, + const gsl::not_null*> + inverse_spatial_metric, + const gsl::not_null*> gamma1, + const gsl::not_null*> gamma2, + const gsl::not_null*> + inverse_spatial_metric_2, + + const std::optional>& + /*face_mesh_velocity*/, + const tnsr::i& normal_covector, + const tnsr::I& normal_vector, + const Scalar& psi_interior, + const Scalar& pi_interior, + const tnsr::i& phi_interior, + const Scalar& lapse_interior, + const tnsr::I& shift_interior, + const tnsr::II& + inverse_spatial_metric_interior, + const Scalar& gamma1_interior, + const Scalar& gamma2_interior, + const Scalar& /*dt_psi*/, + const tnsr::i& /*d_psi*/, + const tnsr::ij& /*d_phi*/, + const Variables>>& + worldtube_vars) const { + ASSERT(worldtube_vars.number_of_grid_points() == get(lapse_interior).size(), + "The worldtube solution has the wrong number of grid points: " + << worldtube_vars.number_of_grid_points()); + *lapse = lapse_interior; + *shift = shift_interior; + *inverse_spatial_metric = inverse_spatial_metric_interior; + *inverse_spatial_metric_2 = inverse_spatial_metric_interior; + *gamma1 = gamma1_interior; + *gamma2 = gamma2_interior; + + const auto& psi_worldtube = get(worldtube_vars); + const auto& pi_worldtube = get(worldtube_vars); + const auto& phi_worldtube = get>(worldtube_vars); + + // `VMinus` is calculated from the worldtube solution, the other + // characteristic fields are set to the corresponding values of the evolved + // variables. This is equivalent to giving no boundary conditions on those + // fields. + Variables, Tags::VPlus, Tags::VMinus>> + char_fields_mixed(get(gamma1_interior).size()); + + auto& v_psi = get(char_fields_mixed); + auto& v_zero = get>(char_fields_mixed); + auto& v_plus = get(char_fields_mixed); + auto& v_minus = get(char_fields_mixed); + + // use allocation for temporary + v_psi = dot_product(normal_vector, phi_interior); + v_zero = tenex::evaluate(phi_interior(ti::i) - + normal_covector(ti::i) * v_psi()); + v_plus = tenex::evaluate(pi_interior() + v_psi() - + gamma2_interior() * psi_interior()); + v_psi = psi_interior; + v_minus = tenex::evaluate(pi_worldtube() - + normal_vector(ti::I) * phi_worldtube(ti::i) - + gamma2_interior() * psi_worldtube()); + + evolved_fields_from_characteristic_fields(psi, pi, phi, gamma2_interior, + v_psi, v_zero, v_plus, v_minus, + normal_covector); + return {}; +} + +template class Worldtube<3>; +} // namespace CurvedScalarWave::BoundaryConditions diff --git a/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.hpp b/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.hpp new file mode 100644 index 000000000000..112a1699ea38 --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.hpp @@ -0,0 +1,146 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Evolution/BoundaryConditions/Type.hpp" +#include "Evolution/Systems/CurvedScalarWave/BoundaryConditions/BoundaryCondition.hpp" +#include "Evolution/Systems/CurvedScalarWave/Tags.hpp" +#include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp" +#include "Parallel/CharmPupable.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace CurvedScalarWave::BoundaryConditions { + +/*! + * \brief Sets boundary conditions for the elements abutting the worldtube using + * a combination of constraint-preserving boundary conditions and the local + * solution evolved inside the worldtube. + * + * \details After extensive experimentation we found this set of boundary + * conditions to be optimal. They are formulated in terms of characteristic + * fields. + * + * Boundary conditions for \f$w^0_\psi\f$ \f$w^0_i\f$ are formulated + * by demanding that there are no constraint violations flowing into the + * numerical domain and are formulated as corrections to the time derivative of + * the evolved variables directly. The derivation is described in + * \ref ConstraintPreservingSphericalRadiation . + * + * Boundary conditions for \f$\w^-\f$ are formulated by evaluating the + * analytical solution of the worldtube at the grid points of each element + * abutting the worldtube. The data is updated and saved to + * \ref CurvedScalarWave::Worldtube::Tags::WorldtubeSolution each time step. + * It is then treated like a ghost field and applied with the chosen numerical + * flux. So far only the upwind flux has been tried. + * + * We tried several other combinations such as setting all fields from the + * worldtube solution but found that this caused major constraint violations to + * flow out of the worldtube. + * + * \note We found that, depending on the worldtube size, a fairly high value for + * \f$\gamma_2\f$ of ca. 10 is required near the worldtube to ensure a stable + * evolution when using these boundary conditions. + */ +template +class Worldtube final : public BoundaryConditions::BoundaryCondition { + public: + using options = tmpl::list<>; + static constexpr Options::String help{ + "Boundary conditions set by the worldtube. w^- will be set by the " + "internal worldtube solution, w^psi and w^0_i are fixed by constraint " + "preserving boundary conditions on the time derivative."}; + + Worldtube() = default; + explicit Worldtube(CkMigrateMessage* msg); + + WRAPPED_PUPable_decl_base_template( + domain::BoundaryConditions::BoundaryCondition, Worldtube); + + auto get_clone() const -> std::unique_ptr< + domain::BoundaryConditions::BoundaryCondition> override; + + static constexpr evolution::BoundaryConditions::Type bc_type = + evolution::BoundaryConditions::Type::GhostAndTimeDerivative; + + void pup(PUP::er& p) override; + + using dg_interior_temporary_tags = tmpl::list< + gr::Tags::Lapse, + gr::Tags::Shift, + gr::Tags::InverseSpatialMetric, + Tags::ConstraintGamma1, Tags::ConstraintGamma2>; + + using dg_gridless_tags = + tmpl::list>; + + using dg_interior_evolved_variables_tags = + tmpl::list>; + using dg_interior_dt_vars_tags = tmpl::list<::Tags::dt>; + using dg_interior_deriv_vars_tags = tmpl::list< + ::Tags::deriv, Frame::Inertial>, + ::Tags::deriv, tmpl::size_t, Frame::Inertial>>; + + std::optional dg_time_derivative( + gsl::not_null*> dt_psi_correction, + gsl::not_null*> dt_pi_correction, + gsl::not_null*> + dt_phi_correction, + const std::optional>& face_mesh_velocity, + const tnsr::i& normal_covector, + const tnsr::I& normal_vector, + const Scalar& /*psi*/, const Scalar& /*pi*/, + const tnsr::i& phi, const Scalar& lapse, + const tnsr::I& shift, + const tnsr::II& /*inverse_spatial_metric*/, + const Scalar& gamma1, const Scalar& gamma2, + const Scalar& /*dt_psi*/, + const tnsr::i& d_psi, + const tnsr::ij& d_phi, + const Variables>>& /*worldtube_vars*/) const; + + std::optional dg_ghost( + const gsl::not_null*> psi, + const gsl::not_null*> pi, + const gsl::not_null*> phi, + const gsl::not_null*> lapse, + const gsl::not_null*> shift, + const gsl::not_null*> + inverse_spatial_metric, + const gsl::not_null*> gamma1, + const gsl::not_null*> gamma2, + const gsl::not_null*> + inverse_spatial_metric_2, + + const std::optional>& + /*face_mesh_velocity*/, + const tnsr::i& normal_covector, + const tnsr::I& normal_vector, + const Scalar& psi_interior, + const Scalar& pi_interior, + const tnsr::i& phi_interior, + const Scalar& lapse_interior, + const tnsr::I& shift_interior, + const tnsr::II& + inverse_spatial_metric_interior, + const Scalar& gamma1_interior, + const Scalar& gamma2_interior, + const Scalar& /*dt_psi*/, + const tnsr::i& d_psi, + const tnsr::ij& /*d_phi*/, + const Variables>>& + worldtube_vars) const; +}; +} // namespace CurvedScalarWave::BoundaryConditions diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Test_Worldtube.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Test_Worldtube.cpp new file mode 100644 index 000000000000..66b868eb5ce1 --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Test_Worldtube.cpp @@ -0,0 +1,227 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/EagerMath/Magnitude.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.hpp" +#include "Evolution/Systems/CurvedScalarWave/BoundaryCorrections/UpwindPenalty.hpp" +#include "Evolution/Systems/CurvedScalarWave/Characteristics.hpp" +#include "Evolution/Systems/CurvedScalarWave/System.hpp" +#include "Evolution/Systems/CurvedScalarWave/Tags.hpp" +#include "Framework/SetupLocalPythonEnvironment.hpp" +#include "Framework/TestCreation.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "Helpers/Evolution/DiscontinuousGalerkin/BoundaryConditions.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace { + +// here we check that the dg_ghost field returns evolved variables which +// correspond to `VMinus` being calculated by the worldtube solution and all +// other fields calculated by the interior solution. +void test_dg_ghost(const gsl::not_null gen) { + static constexpr size_t Dim = 3; + std::uniform_real_distribution<> dist(-10., 10.); + gr::Solutions::KerrSchild kerr_schild{2., {0.1, 0.2, 0.3}, {0.4, 0.6, 0.9}}; + const size_t num_points = 1000.; + const auto coords = + make_with_random_values>(gen, dist, num_points); + const auto spacetime_variables = kerr_schild.variables( + coords, 0., + tmpl::list< + gr::Tags::Lapse, + gr::Tags::Shift, + gr::Tags::InverseSpatialMetric>{}); + const auto& lapse = get>(spacetime_variables); + const auto& shift = get>( + spacetime_variables); + const auto& inverse_spatial_metric = + get>( + spacetime_variables); + + auto normal_covector = + make_with_random_values>(gen, dist, num_points); + const auto covector_mag = magnitude(normal_covector, inverse_spatial_metric); + for (size_t i = 0; i < Dim; ++i) { + normal_covector.get(i) /= covector_mag.get(); + } + const auto normal_vector = tenex::evaluate( + inverse_spatial_metric(ti::I, ti::J) * normal_covector(ti::j)); + CHECK_ITERABLE_APPROX(dot_product(normal_covector, normal_vector).get(), + DataVector(num_points, 1.)); + + const auto psi_interior = + make_with_random_values>(gen, dist, num_points); + const auto pi_interior = + make_with_random_values>(gen, dist, num_points); + const auto phi_interior = + make_with_random_values>(gen, dist, num_points); + const auto gamma1 = + make_with_random_values>(gen, dist, num_points); + const auto gamma2 = + make_with_random_values>(gen, dist, num_points); + const auto dt_psi = + make_with_random_values>(gen, dist, num_points); + const auto d_psi = + make_with_random_values>(gen, dist, num_points); + const auto d_phi = + make_with_random_values>(gen, dist, num_points); + const auto worldtube_vars = make_with_random_values>>>(gen, dist, num_points); + + Scalar psi_res(num_points); + Scalar pi_res(num_points); + tnsr::i phi_res(num_points); + Scalar lapse_res(num_points); + tnsr::I shift_res(num_points); + Scalar gamma1_res(num_points); + Scalar gamma2_res(num_points); + tnsr::II inverse_spatial_metric_res_1( + num_points); + tnsr::II inverse_spatial_metric_res_2( + num_points); + + CurvedScalarWave::BoundaryConditions::Worldtube worldtube_bcs{}; + worldtube_bcs.dg_ghost( + make_not_null(&psi_res), make_not_null(&pi_res), make_not_null(&phi_res), + make_not_null(&lapse_res), make_not_null(&shift_res), + make_not_null(&inverse_spatial_metric_res_1), make_not_null(&gamma1_res), + make_not_null(&gamma2_res), make_not_null(&inverse_spatial_metric_res_2), + std::nullopt, normal_covector, normal_vector, psi_interior, pi_interior, + phi_interior, lapse, shift, inverse_spatial_metric, gamma1, gamma2, + dt_psi, d_psi, d_phi, worldtube_vars); + + CHECK(lapse_res == lapse); + CHECK(shift_res == shift); + CHECK(inverse_spatial_metric_res_1 == inverse_spatial_metric); + CHECK(inverse_spatial_metric_res_2 == inverse_spatial_metric); + CHECK(gamma1_res == gamma1); + CHECK(gamma2_res == gamma2); + + const auto char_fields_interior = CurvedScalarWave::characteristic_fields( + gamma2, inverse_spatial_metric, psi_interior, pi_interior, phi_interior, + normal_covector); + + const auto char_fields_worldtube = CurvedScalarWave::characteristic_fields( + gamma2, inverse_spatial_metric, + get(worldtube_vars), + get(worldtube_vars), + get>(worldtube_vars), normal_covector); + + const auto char_fields_res = CurvedScalarWave::characteristic_fields( + gamma2, inverse_spatial_metric, psi_res, pi_res, phi_res, + normal_covector); + Approx approx = Approx::custom().epsilon(1.e-12).scale(1.); + CHECK_ITERABLE_CUSTOM_APPROX( + get(char_fields_interior), + get(char_fields_res), approx); + CHECK_ITERABLE_CUSTOM_APPROX( + get>(char_fields_interior), + get>(char_fields_res), approx); + CHECK_ITERABLE_CUSTOM_APPROX( + get(char_fields_interior), + get(char_fields_res), approx); + CHECK_ITERABLE_CUSTOM_APPROX( + get(char_fields_worldtube), + get(char_fields_res), approx); +} + +// this takes a variables on the C++ side filled with 1 and returns only +// 1 on the python side, point by point. There is little point in trying to +// transfer the variables to the python side because it is unclear how the +// different components should be packaged and ordered. +template +struct ConvertWorldtubeSolution { + using unpacked_container = double; + using packed_container = Variables< + tmpl::list>>; + using packed_type = packed_container; + static constexpr size_t face_size = NumGridPoints * NumGridPoints; + + static packed_container create_container() { + return packed_container{face_size, 1.}; + } + + static inline unpacked_container unpack(const packed_container /*packed*/, + const size_t /*grid_point_index*/) { + return 1.; + } + + static inline void pack(const gsl::not_null packed, + const unpacked_container /*unpacked*/, + const size_t /*grid_point_index*/) { + *packed = packed_container{face_size, 1.}; + } + + static inline size_t get_size(const packed_container& packed) { + return packed.number_of_grid_points(); + } +}; + +void test_python(const gsl::not_null gen) { + static constexpr size_t Dim = 3; + std::uniform_real_distribution<> dist(-10., 10.); + static constexpr size_t num_grid_points = 5; + auto box = db::create>>( + ConvertWorldtubeSolution::create_container()); + + namespace helpers = TestHelpers::evolution::dg; + helpers::test_boundary_condition_with_python< + CurvedScalarWave::BoundaryConditions::Worldtube, + CurvedScalarWave::BoundaryConditions::BoundaryCondition, + CurvedScalarWave::System, + tmpl::list>, + tmpl::list>, + tmpl::list>>( + gen, + "Evolution.Systems.CurvedScalarWave.BoundaryConditions." + "Worldtube", + tuples::TaggedTuple< + helpers::Tags::PythonFunctionForErrorMessage<>, + helpers::Tags::PythonFunctionName< + ::Tags::dt>, + helpers::Tags::PythonFunctionName< + ::Tags::dt>, + helpers::Tags::PythonFunctionName< + ::Tags::dt>>, + helpers::Tags::PythonFunctionName, + helpers::Tags::PythonFunctionName, + helpers::Tags::PythonFunctionName>, + helpers::Tags::PythonFunctionName>, + helpers::Tags::PythonFunctionName>, + helpers::Tags::PythonFunctionName< + gr::Tags::InverseSpatialMetric>, + helpers::Tags::PythonFunctionName< + CurvedScalarWave::Tags::ConstraintGamma1>, + helpers::Tags::PythonFunctionName< + CurvedScalarWave::Tags::ConstraintGamma2>>{ + "error", "dt_psi_worldtube", "dt_pi_worldtube", "dt_phi_worldtube", + "psi", "pi", "phi", "lapse", "shift", "inverse_spatial_metric", + "gamma1", "gamma2"}, + "Worldtube:\n", Index{num_grid_points}, box, + tuples::TaggedTuple<>{}); +} + +} // namespace +SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.BoundaryConditions.Worldtube", + "[Unit][Evolution]") { + pypp::SetupLocalPythonEnvironment local_python_env{""}; + MAKE_GENERATOR(gen); + test_dg_ghost(make_not_null(&gen)); + test_python(make_not_null(&gen)); +} diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.py b/tests/Unit/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.py new file mode 100644 index 000000000000..449d1315b57e --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/BoundaryConditions/Worldtube.py @@ -0,0 +1,141 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +import numpy as np +from Evolution.Systems.CurvedScalarWave import Characteristics + + +def error(face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, + lapse, shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, + d_phi, worldtube_vars, dim): + return None + + +def dt_psi_worldtube(face_mesh_velocity, normal_covector, normal_vector, psi, + pi, phi, lapse, shift, inverse_spatial_metric, gamma1, + gamma2, dt_psi, d_psi, d_phi, worldtube_vars, dim): + char_speed_psi = Characteristics.char_speed_vpsi(gamma1, lapse, shift, + normal_covector) + + if face_mesh_velocity is not None: + char_speed_psi -= np.dot(normal_covector, face_mesh_velocity) + + return np.dot(normal_vector, d_psi - phi) * min(0., char_speed_psi) + + +def dt_phi_worldtube(face_mesh_velocity, normal_covector, normal_vector, psi, + pi, phi, lapse, shift, inverse_spatial_metric, gamma1, + gamma2, dt_psi, d_psi, d_phi, worldtube_vars, dim): + char_speed_zero = Characteristics.char_speed_vzero(gamma1, lapse, shift, + normal_covector) + if face_mesh_velocity is not None: + char_speed_zero -= np.dot(normal_covector, face_mesh_velocity) + return np.einsum("ij,j", d_phi.T - d_phi, normal_vector) * min( + 0, char_speed_zero) + + +def dt_pi_worldtube(face_mesh_velocity, normal_covector, normal_vector, psi, + pi, phi, lapse, shift, inverse_spatial_metric, gamma1, + gamma2, dt_psi, d_psi, d_phi, worldtube_vars, dim): + return gamma2 * dt_psi_worldtube( + face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, + lapse, shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, + d_phi, worldtube_vars, dim) + + +def psi(face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, + lapse, shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, + d_phi, worldtube_vars, dim): + return psi + + +def pi(face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, lapse, + shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, d_phi, + worldtube_vars, dim): + # The worldtube data on the C++ side is a Variables filled with ones. + # Since we do not want a Variables on the python side, we hard-code it here + wt_psi = 1. + wt_pi = 1. + wt_phi = np.ones(3) + + vpsi_interior = Characteristics.char_field_vpsi(gamma2, + inverse_spatial_metric, + psi, pi, phi, + normal_covector) + vzero_interior = Characteristics.char_field_vzero(gamma2, + inverse_spatial_metric, + psi, pi, phi, + normal_covector) + vplus_interior = Characteristics.char_field_vplus(gamma2, + inverse_spatial_metric, + psi, pi, phi, + normal_covector) + vminus_wt = Characteristics.char_field_vminus(gamma2, + inverse_spatial_metric, + wt_psi, wt_pi, wt_phi, + normal_covector) + + return Characteristics.evol_field_pi(gamma2, vpsi_interior, vzero_interior, + vplus_interior, vminus_wt, + normal_covector) + + +def phi(face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, + lapse, shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, + d_phi, worldtube_vars, dim): + # The worldtube data on the C++ side is a Variables filled with ones. + # Since we do not want a Variables on the python side, we hard-code it here + wt_psi = 1. + wt_pi = 1. + wt_phi = np.ones(3) + + vpsi_interior = Characteristics.char_field_vpsi(gamma2, + inverse_spatial_metric, + psi, pi, phi, + normal_covector) + vzero_interior = Characteristics.char_field_vzero(gamma2, + inverse_spatial_metric, + psi, pi, phi, + normal_covector) + vplus_interior = Characteristics.char_field_vplus(gamma2, + inverse_spatial_metric, + psi, pi, phi, + normal_covector) + vminus_wt = Characteristics.char_field_vminus(gamma2, + inverse_spatial_metric, + wt_psi, wt_pi, wt_phi, + normal_covector) + return Characteristics.evol_field_phi(gamma2, vpsi_interior, + vzero_interior, vplus_interior, + vminus_wt, normal_covector) + + +def lapse(face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, + lapse, shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, + d_phi, worldtube_vars, dim): + return lapse + + +def shift(face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, + lapse, shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, + d_phi, worldtube_vars, dim): + return shift + + +def gamma1(face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, + lapse, shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, + d_phi, worldtube_vars, dim): + return gamma1 + + +def gamma2(face_mesh_velocity, normal_covector, normal_vector, psi, pi, phi, + lapse, shift, inverse_spatial_metric, gamma1, gamma2, dt_psi, d_psi, + d_phi, worldtube_vars, dim): + return gamma2 + + +def inverse_spatial_metric(face_mesh_velocity, normal_covector, normal_vector, + psi, pi, phi, lapse, shift, inverse_spatial_metric, + gamma1, gamma2, dt_psi, d_psi, d_phi, + worldtube_vars, dim): + return inverse_spatial_metric diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt index 9e6be3f60d7f..b052d615fb7e 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt @@ -7,6 +7,7 @@ set(LIBRARY_SOURCES BoundaryCorrections/Test_UpwindPenalty.cpp BoundaryConditions/Test_ConstraintPreservingSphericalRadiation.cpp BoundaryConditions/Test_DemandOutgoingCharSpeeds.cpp + BoundaryConditions/Test_Worldtube.cpp Test_BackgroundSpacetime.cpp Test_CalculateGrVars.cpp Test_Constraints.cpp