From f9e1d1e1874385e1e1a53eef03838649059f0fdd Mon Sep 17 00:00:00 2001 From: Francois Foucart Date: Mon, 18 Mar 2024 13:44:25 -0700 Subject: [PATCH] Add test to evolve packet on Kerr background Also add ghost zones to MC evolution in Kerr test --- .../Particles/MonteCarlo/CMakeLists.txt | 1 + .../MonteCarlo/Test_EvolvePackets.cpp | 220 +++++++++++++++++- 2 files changed, 218 insertions(+), 3 deletions(-) diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt b/tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt index 8117d6b1c1510..3fb408bb45dbd 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt +++ b/tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt @@ -31,6 +31,7 @@ target_link_libraries( DgSubcell GeneralRelativity GeneralRelativityHelpers + GeneralRelativitySolutions H5 Hydro HydroHelpers diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_EvolvePackets.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_EvolvePackets.cpp index 6163046f46616..de2f587eeb5c6 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_EvolvePackets.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_EvolvePackets.cpp @@ -10,12 +10,15 @@ #include "Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp" #include "Framework/TestHelpers.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Units.hpp" using hydro::units::nuclear::proton_mass; -SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloEvolution", - "[Unit][Evolution]") { +namespace { + +void test_evolve_minkowski() { const size_t mesh_1d_size = 2; const Mesh<3> mesh(mesh_1d_size, Spectral::Basis::FiniteDifference, Spectral::Quadrature::CellCentered); @@ -103,7 +106,6 @@ SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloEvolution", Particles::MonteCarlo::Packet packet(1, 1.0, 0, 0.0, -0.75, -1.0, -1.0, 1.0, 1.0, 0.0, 0.0); - packet.renormalize_momentum(inv_spatial_metric, lapse); CHECK(packet.momentum_upper_t == 1.0); @@ -170,3 +172,215 @@ SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloEvolution", CHECK(fabs(get(coupling_rho_ye)[ext_idx_0 + 1] + proton_mass * 5.0e-61) < 1.e-72); } +tnsr::I spatial_coords_logical( + const Mesh<3>& mesh) { + const DataVector used_for_size(mesh.number_of_grid_points()); + auto x = make_with_value>( + used_for_size, 0.0); + const Index<3>& extents = mesh.extents(); + Index<3> cur_extents{0, 0, 0}; + for (cur_extents[0] = 0; cur_extents[0] < extents[0]; cur_extents[0]++) { + for (cur_extents[1] = 0; cur_extents[1] < extents[1]; cur_extents[1]++) { + for (cur_extents[2] = 0; cur_extents[2] < extents[2]; cur_extents[2]++) { + const size_t storage_index = mesh.storage_index(cur_extents); + x.get(0)[storage_index] = + -1.0 + 2.0 * (static_cast(cur_extents[0]) + 0.5) / + static_cast(extents[0]); + x.get(1)[storage_index] = + -1.0 + 2.0 * (static_cast(cur_extents[1]) + 0.5) / + static_cast(extents[1]); + x.get(2)[storage_index] = + -1.0 + 2.0 * (static_cast(cur_extents[2]) + 0.5) / + static_cast(extents[2]); + } + } + } + return x; +} + // inertial coordinates are set such that the element used for the + // evolution is centered on (6,0,0), and the jacobian matrix + // is the identity matrix. +tnsr::I spatial_coords_inertial( + tnsr::I logical_coords) { + auto x = make_with_value>( + logical_coords.get(0), 0.0); + x.get(0) = logical_coords.get(0) + 6.0; + x.get(1) = logical_coords.get(1); + x.get(2) = logical_coords.get(2); + return x; +} + // Evolve a single Monte-Carlo packet along Kerr geodesic + // with check that final position matches analytical expectations +void test_evolve_kerr(const size_t mesh_size_1d) { + MAKE_GENERATOR(generator); + + // Parameters for KerrSchild solution + const double mass = 1.01; + const std::array spin{{0.0, 0.0, 0.0}}; + const std::array center{{0.0, 0.0, 0.0}}; + const double t = 1.3; + // Evaluate solution + gr::Solutions::KerrSchild solution(mass, spin, center); + + // Set domain and coordintes + const Mesh<3> mesh(mesh_size_1d, Spectral::Basis::FiniteDifference, + Spectral::Quadrature::CellCentered); + const size_t num_ghost_zones = 1; + const size_t extended_mesh_1d_size = mesh_size_1d + 2 * num_ghost_zones; + const Mesh<3> extended_mesh(extended_mesh_1d_size, + Spectral::Basis::FiniteDifference, + Spectral::Quadrature::CellCentered); + const size_t dv_size = mesh.number_of_grid_points(); + const size_t extended_dv_size = extended_mesh.number_of_grid_points(); + DataVector zero_dv(dv_size, 0.0); + DataVector extended_zero_dv(extended_dv_size, 0.0); + const auto mesh_coordinates = spatial_coords_logical(mesh); + const auto inertial_coordinates = spatial_coords_inertial(mesh_coordinates); + + // Mesh velocity set to std::null for now + const std::optional> mesh_velocity = + std::nullopt; + + // Compute metric quantities + const auto vars = solution.variables( + inertial_coordinates, t, + typename gr::Solutions::KerrSchild::tags{}); + const auto& lapse = get>(vars); + const auto& deriv_lapse = get>(vars); + const auto& shift = + get>(vars); + const auto& deriv_shift = get>(vars); + const auto& spatial_metric = + get>(vars); + const auto& inverse_spatial_metric = + get>(vars); + Scalar cell_light_crossing_time = + make_with_value>(zero_dv, 1.0); + const auto& deriv_spatial_metric = + get>(vars); + + tnsr::iJJ deriv_inverse_spatial_metric = + make_with_value>(lapse, 0.0); + for (size_t i = 0; i < 3; i++) { + for (size_t j = i; j < 3; j++) { + for (size_t k = 0; k < 3; k++) { + for (size_t l = 0; l < 3; l++) { + for (size_t m = 0; m < 3; m++) { + deriv_inverse_spatial_metric.get(k, i, j) -= + inverse_spatial_metric.get(i, l) * + inverse_spatial_metric.get(j, m) * + deriv_spatial_metric.get(k, l, m); + } + } + } + } + } + + // Fluid quantities + Scalar lorentz_factor = + make_with_value>(lapse, 1.0); + tnsr::i lower_spatial_four_velocity = + make_with_value>(lapse, 0.0); + + // Initialize packet on the x-axis + const double dx = 2.0 / mesh_size_1d; + Particles::MonteCarlo::Packet packet(0, 1.0, 5, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, + 1.0, 0.0); + // Self-consistency: update index of closest point and p^t + std::array closest_point_index_3d{0, 0, 0}; + for (size_t d = 0; d < 3; d++) { + gsl::at(closest_point_index_3d, d) = + std::floor((packet.coordinates[d] - mesh_coordinates.get(d)[0]) / + (2.0 / mesh_size_1d) + + 0.5); + } + const Index<3>& extents = mesh.extents(); + packet.index_of_closest_grid_point = + closest_point_index_3d[0] + + extents[0] * + (closest_point_index_3d[1] + extents[1] * closest_point_index_3d[2]); + packet.renormalize_momentum(inverse_spatial_metric, lapse); + + // Jacobians set to identity for now + InverseJacobian + inverse_jacobian_inertial_to_fluid = make_with_value< + InverseJacobian>(lapse, + 0.0); + inverse_jacobian_inertial_to_fluid.get(0, 0) = 1.0; + inverse_jacobian_inertial_to_fluid.get(1, 1) = 1.0; + inverse_jacobian_inertial_to_fluid.get(2, 2) = 1.0; + inverse_jacobian_inertial_to_fluid.get(3, 3) = 1.0; + Jacobian + + jacobian_inertial_to_fluid = make_with_value< + Jacobian>(lapse, 0.0); + jacobian_inertial_to_fluid.get(0, 0) = 1.0; + jacobian_inertial_to_fluid.get(1, 1) = 1.0; + jacobian_inertial_to_fluid.get(2, 2) = 1.0; + jacobian_inertial_to_fluid.get(3, 3) = 1.0; + // Logical to inertial inverse jacobian, also identity for now + InverseJacobian + inverse_jacobian = + make_with_value>(lapse, 0.0); + inverse_jacobian.get(0, 0) = 1.0; + inverse_jacobian.get(1, 1) = 1.0; + inverse_jacobian.get(2, 2) = 1.0; + + // Interaction rates + const std::array energy_at_bin_center = {2.0, 5.0}; + std::array, 2> absorption_opacity = { + std::array{{extended_zero_dv, extended_zero_dv}}, + std::array{{extended_zero_dv, extended_zero_dv}}}; + std::array, 2> scattering_opacity = { + std::array{{extended_zero_dv, extended_zero_dv}}, + std::array{{extended_zero_dv, extended_zero_dv}}}; + + Scalar coupling_tilde_tau = + make_with_value>(extended_zero_dv, 0.0); + Scalar coupling_rho_ye = + make_with_value>(extended_zero_dv, 0.0); + tnsr::i coupling_tilde_s = + make_with_value>(extended_zero_dv, + 0.0); + + std::vector packets{packet}; + Particles::MonteCarlo::TemplatedLocalFunctions<2, 2> MonteCarloStruct; + // Getting CFL constant + const double cfl_constant = 0.25; + const double final_time = 1.0; + const double dt_step = cfl_constant * dx; + double current_time = 0.0; + packets[0].renormalize_momentum(inverse_spatial_metric, lapse); + + //time evolution with step size dt + while (current_time < final_time) { + double dt = std::min(dt_step, final_time - current_time); + MonteCarloStruct.evolve_packets( + &packets, &generator, &coupling_tilde_tau, &coupling_tilde_s, + &coupling_rho_ye, current_time + dt, mesh, mesh_coordinates, + num_ghost_zones, absorption_opacity, scattering_opacity, + energy_at_bin_center, lorentz_factor, lower_spatial_four_velocity, + lapse, shift, deriv_lapse, deriv_shift, deriv_inverse_spatial_metric, + spatial_metric, inverse_spatial_metric, cell_light_crossing_time, + mesh_velocity, inverse_jacobian, jacobian_inertial_to_fluid, + inverse_jacobian_inertial_to_fluid); + current_time += dt; + } + const double final_r = sqrt(pow(packets[0].coordinates.get(0) + 6.0, 2) + + pow(packets[0].coordinates.get(1), 2)+pow(packets[0].coordinates.get(2),2)); + + //Check r against geodesic evolution value obtained from python code + CHECK(std::abs(final_r - 6.298490) < 1e-2); +} + +} // namespace + +SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloEvolution", + "[Unit][Evolution]") { + test_evolve_minkowski(); + test_evolve_kerr(15); +}