From d93d83e385c78374dddd40fd99cd3dc10fd8b48b Mon Sep 17 00:00:00 2001 From: Francois Foucart Date: Thu, 31 Oct 2024 14:01:18 -0400 Subject: [PATCH] Take single MC step without communication --- .../MonteCarlo/EvolveMonteCarlo.hpp | 53 ++++--- .../MonteCarlo/Actions/CMakeLists.txt | 1 + .../Actions/InitializeMonteCarlo.hpp | 149 ++++++++++++++++++ .../MonteCarlo/Actions/TimeStepActions.hpp | 35 ++-- .../Particles/MonteCarlo/CellVolume.cpp | 4 +- .../Particles/MonteCarlo/CellVolume.hpp | 11 +- .../ImplicitMonteCarloCorrections.tpp | 6 +- .../MonteCarlo/NeutrinoInteractionTable.cpp | 11 +- src/Evolution/Particles/MonteCarlo/System.hpp | 11 +- src/Evolution/Particles/MonteCarlo/Tags.hpp | 2 +- .../Particles/MonteCarlo/TakeTimeStep.tpp | 6 +- .../MonteCarlo/TemplatedLocalFunctions.hpp | 2 +- .../MultiLinearSpanInterpolation.hpp | 29 +++- .../Particles/MonteCarlo/Test_CellVolume.cpp | 6 +- .../MonteCarlo/Test_TakeTimeStep.cpp | 4 +- .../MonteCarlo/Test_TimeStepAction.cpp | 4 +- 16 files changed, 267 insertions(+), 67 deletions(-) create mode 100644 src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp index e261164f660d..3511fbed7d8e 100644 --- a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp +++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp @@ -27,6 +27,7 @@ #include "Evolution/Initialization/Evolution.hpp" #include "Evolution/Initialization/Limiter.hpp" #include "Evolution/Initialization/SetVariables.hpp" +#include "Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp" #include "Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp" #include "Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp" #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp" @@ -69,7 +70,9 @@ #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp" #include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/ConstantM1.hpp" #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp" +#include "PointwiseFunctions/Hydro/LowerSpatialFourVelocity.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" +#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp" #include "Time/Actions/AdvanceTime.hpp" #include "Time/Actions/CleanHistory.hpp" #include "Time/Actions/RecordTimeStepperData.hpp" @@ -114,7 +117,6 @@ struct EvolutionMetavars { using TimeStepperBase = TimeStepper; static constexpr bool use_dg_subcell = true; - using initial_data_tag = evolution::initial_data::Tags::InitialData; using initial_data_list = grmhd::ValenciaDivClean::InitialData::initial_data_list; using equation_of_state_tag = hydro::Tags::GrmhdEquationOfState; @@ -136,6 +138,10 @@ struct EvolutionMetavars { ::Events::Tags::ObserverDetInvJacobianCompute< Frame::ElementLogical, Frame::Inertial>>; + using analytic_variables_tags = typename system::variables_tag::tags_list; + using analytic_compute = evolution::Tags::AnalyticSolutionsCompute< + volume_dim, analytic_variables_tags, false, initial_data_list>; + struct factory_creation : tt::ConformsTo { using factory_classes = tmpl::map< @@ -204,6 +210,12 @@ struct EvolutionMetavars { evolution::dg::subcell::BackgroundGrVars>, Actions::MutateApply, + Initialization::Actions::InitializeMCTags, + Initialization::Actions::AddComputeTags>>, // Initialization::TimeStepperHistory Parallel::Actions::TerminatePhase>; @@ -211,22 +223,25 @@ struct EvolutionMetavars { EvolutionMetavars, tmpl::list< Parallel::PhaseActions -//, -// Parallel::PhaseActions< -// Parallel::Phase::InitializeTimeStepperHistory, -// SelfStart::self_start_procedure>, + initialization_actions>, + //, + // Parallel::PhaseActions< + // Parallel::Phase::InitializeTimeStepperHistory, + // SelfStart::self_start_procedure>, -// Parallel::PhaseActions>, + // Parallel::PhaseActions>, -// Parallel::PhaseActions< -// Parallel::Phase::Evolve, -// tmpl::list>>>; + Parallel::PhaseActions< + Parallel::Phase::Evolve, + tmpl::list, + Parallel::Actions::TerminatePhase>> + // tmpl::list> >>; struct registration @@ -240,10 +255,10 @@ struct EvolutionMetavars { observers::ObserverWriter, dg_element_array>; - using const_global_cache_tags = - tmpl::list>; + using const_global_cache_tags = tmpl::list< + equation_of_state_tag, evolution::initial_data::Tags::InitialData, + Particles::MonteCarlo::Tags::InteractionRatesTable>; static constexpr Options::String help{ "Evolve Monte Carlo transport (without coupling to hydro).\n\n"}; diff --git a/src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt b/src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt index 2ab80257813e..d9e1b37ef3f3 100644 --- a/src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt +++ b/src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt @@ -5,5 +5,6 @@ spectre_target_headers( ${LIBRARY} INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src HEADERS + InitializeMonteCarlo.hpp TimeStepActions.hpp ) diff --git a/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp b/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp new file mode 100644 index 000000000000..21a1c59f67d3 --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp @@ -0,0 +1,149 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/Initialization/InitialData.hpp" +#include "Evolution/Particles/MonteCarlo/MortarData.hpp" +#include "Evolution/Particles/MonteCarlo/Packet.hpp" +#include "Evolution/Particles/MonteCarlo/Tags.hpp" +#include "Parallel/AlgorithmExecution.hpp" +#include "Parallel/GlobalCache.hpp" +#include "ParallelAlgorithms/Initialization/MutateAssign.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" +#include "Utilities/TMPL.hpp" + +/// \cond +namespace evolution::initial_data::Tags { +struct InitialData; +} // namespace evolution::initial_data::Tags + +namespace tuples { +template +class TaggedTuple; +} // namespace tuples + +namespace Parallel { +template +class GlobalCache; +} // namespace Parallel +/// \endcond + +namespace Initialization::Actions { + +/// \ingroup InitializationGroup +/// \brief Allocate variables needed for evolution of Monte Carlo transport +/// +/// Uses: +/// - evolution::dg::subcell::Tags::Mesh +/// +/// DataBox changes: +/// - Adds: +/// * Particles::MonteCarlo::Tags::PacketsOnElement +/// * Particles::MonteCarlo::Tags::RandomNumberGenerator +/// * Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission +/// +/// * Particles::MonteCarlo::Tags::CellLightCrossingTime +/// * Background hydro variables +/// * Particles::MonteCarlo::Tags::MortarDataTag +/// +/// - Removes: nothing +/// - Modifies: nothing +template +struct InitializeMCTags { + public: + using hydro_variables_tag = typename System::hydro_variables_tag; + + static constexpr size_t dim = System::volume_dim; + using simple_tags = + tmpl::list, + Particles::MonteCarlo::Tags::CellLightCrossingTime, + hydro_variables_tag, + Particles::MonteCarlo::Tags::MortarDataTag>; + + using compute_tags = tmpl::list<>; + + template + static Parallel::iterable_action_return_t apply( + db::DataBox& box, + const tuples::TaggedTuple& /*inboxes*/, + const Parallel::GlobalCache& /*cache*/, + const ArrayIndex& /*array_index*/, ActionList /*meta*/, + const ParallelComponent* const /*meta*/) { + const size_t num_grid_points = + db::get>(box) + .number_of_grid_points(); + using derived_classes = + tmpl::at; + using HydroVars = typename hydro_variables_tag::type; + call_with_dynamic_type( + &db::get(box), + [&box](const auto* const data_or_solution) { + static constexpr size_t dim = System::volume_dim; + const double initial_time = 0.0; // db::get<::Tags::Time>(box); + const size_t num_grid_points = + db::get>(box) + .number_of_grid_points(); + const auto& inertial_coords = db::get< + evolution::dg::subcell::Tags::Coordinates>( + box); + // Get hydro variables + HydroVars hydro_variables{num_grid_points}; + hydro_variables.assign_subset(evolution::Initialization::initial_data( + *data_or_solution, inertial_coords, initial_time, + typename hydro_variables_tag::tags_list{})); + Initialization::mutate_assign>( + make_not_null(&box), std::move(hydro_variables)); + }); + + typename Particles::MonteCarlo::Tags::PacketsOnElement::type all_packets; + Initialization::mutate_assign< + tmpl::list>( + make_not_null(&box), std::move(all_packets)); + + // Currently seeds with 0 for testing... + typename Particles::MonteCarlo::Tags::RandomNumberGenerator::type rng(0); + Initialization::mutate_assign< + tmpl::list>( + make_not_null(&box), std::move(rng)); + + // Currently hard-code energy at emission; should be set by option + typename Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission< + NeutrinoSpecies>::type packet_energy_at_emission = + make_with_value>( + DataVector{num_grid_points}, 1.e-12); + Initialization::mutate_assign< + tmpl::list>>(make_not_null(&box), + std::move(packet_energy_at_emission)); + + typename Particles::MonteCarlo::Tags::CellLightCrossingTime< + DataVector>::type cell_light_crossing_time(num_grid_points, 1.0); + Initialization::mutate_assign>>( + make_not_null(&box), std::move(cell_light_crossing_time)); + + // Initialize empty mortar data; do we need more at initialization stage? + using MortarData = + typename Particles::MonteCarlo::Tags::MortarDataTag::type; + MortarData mortar_data; + Initialization::mutate_assign< + tmpl::list>>( + make_not_null(&box), std::move(mortar_data)); + + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } +}; + +} // namespace Initialization::Actions diff --git a/src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp b/src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp index f31cca715fee..7d60ecb6dff7 100644 --- a/src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp +++ b/src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp @@ -41,7 +41,8 @@ struct TimeStepMutator { using return_tags = tmpl::list>; + Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission< + NeutrinoSpecies>>; // To do : check carefully DG vs Subcell quantities... everything should // be on the Subcell grid! using argument_tags = tmpl::list< @@ -60,11 +61,11 @@ struct TimeStepMutator { Frame::Inertial>, ::Tags::deriv, tmpl::size_t, Frame::Inertial>, - ::Tags::deriv, - tmpl::size_t, Frame::Inertial>, + ::Tags::deriv, tmpl::size_t, + Frame::Inertial>, gr::Tags::SpatialMetric, gr::Tags::InverseSpatialMetric, - gr::Tags::DetSpatialMetric, + gr::Tags::SqrtDetSpatialMetric, Particles::MonteCarlo::Tags::CellLightCrossingTime, evolution::dg::subcell::Tags::Mesh, evolution::dg::subcell::Tags::Coordinates, @@ -96,10 +97,10 @@ struct TimeStepMutator { const tnsr::i& d_lapse, const tnsr::iJ& d_shift, - const tnsr::iJJ& d_inv_spatial_metric, + const tnsr::ijj& d_spatial_metric, const tnsr::ii& spatial_metric, const tnsr::II& inv_spatial_metric, - const Scalar& determinant_spatial_metric, + const Scalar& sqrt_determinant_spatial_metric, const Scalar& cell_light_crossing_time, const Mesh& mesh, const tnsr::I& mesh_coordinates, const std::optional>& @@ -117,8 +118,8 @@ struct TimeStepMutator { const size_t num_ghost_zones = 1; // Get information stored in various databox containers in // the format expected by take_time_step_on_element - const double start_time = current_step_id.step_time().value(); - const double end_time = next_step_id.step_time().value(); + const double start_time = 0.0; // current_step_id.step_time().value(); + const double end_time = 0.1; // next_step_id.step_time().value(); Scalar det_jacobian_logical_to_inertial(lapse); get(det_jacobian_logical_to_inertial) = 1.0 / get(det_inverse_jacobian_logical_to_inertial); @@ -131,6 +132,22 @@ struct TimeStepMutator { const DirectionalIdMap>& cell_light_crossing_time_ghost = mortar_data.cell_light_crossing_time; + tnsr::iJJ d_inv_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++) { + d_inv_spatial_metric.get(k, i, j) -= + inv_spatial_metric.get(i, l) * inv_spatial_metric.get(j, m) * + d_spatial_metric.get(k, l, m); + } + } + } + } + } + TemplatedLocalFunctions templated_functions; templated_functions.take_time_step_on_element( packets, random_number_generator, single_packet_energy, start_time, @@ -138,7 +155,7 @@ struct TimeStepMutator { rest_mass_density, temperature, lorentz_factor, lower_spatial_four_velocity, lapse, shift, d_lapse, d_shift, d_inv_spatial_metric, spatial_metric, inv_spatial_metric, - determinant_spatial_metric, cell_light_crossing_time, mesh, + sqrt_determinant_spatial_metric, cell_light_crossing_time, mesh, mesh_coordinates, num_ghost_zones, mesh_velocity, inverse_jacobian_logical_to_inertial, det_jacobian_logical_to_inertial, inertial_to_fluid_jacobian, inertial_to_fluid_inverse_jacobian, diff --git a/src/Evolution/Particles/MonteCarlo/CellVolume.cpp b/src/Evolution/Particles/MonteCarlo/CellVolume.cpp index 89f2a71267ee..4f416722fff6 100644 --- a/src/Evolution/Particles/MonteCarlo/CellVolume.cpp +++ b/src/Evolution/Particles/MonteCarlo/CellVolume.cpp @@ -12,13 +12,13 @@ namespace Particles::MonteCarlo { void cell_proper_four_volume_finite_difference( const gsl::not_null*> cell_proper_four_volume, const Scalar& lapse, - const Scalar& determinant_spatial_metric, + const Scalar& sqrt_determinant_spatial_metric, const double time_step, const Mesh<3>& mesh, const Scalar& det_jacobian_logical_to_inertial) { const double cell_logical_volume = 8.0 / static_cast(mesh.number_of_grid_points()); cell_proper_four_volume->get() = - get(lapse) * get(determinant_spatial_metric) * time_step * + get(lapse) * get(sqrt_determinant_spatial_metric) * time_step * cell_logical_volume * get(det_jacobian_logical_to_inertial); } diff --git a/src/Evolution/Particles/MonteCarlo/CellVolume.hpp b/src/Evolution/Particles/MonteCarlo/CellVolume.hpp index f813836ef9a9..52db8d797b98 100644 --- a/src/Evolution/Particles/MonteCarlo/CellVolume.hpp +++ b/src/Evolution/Particles/MonteCarlo/CellVolume.hpp @@ -25,12 +25,11 @@ namespace Particles::MonteCarlo { /// in inertial coordinates, hence the need for /// det_jacobian_logical_to_inertial void cell_proper_four_volume_finite_difference( - gsl::not_null* > cell_proper_four_volume, - const Scalar& lapse, - const Scalar& determinant_spatial_metric, - double time_step, - const Mesh<3>& mesh, - const Scalar& det_jacobian_logical_to_inertial); + gsl::not_null*> cell_proper_four_volume, + const Scalar& lapse, + const Scalar& sqrt_determinant_spatial_metric, double time_step, + const Mesh<3>& mesh, + const Scalar& det_jacobian_logical_to_inertial); /// 3-volume of a cell in inertial coordinate. Note that this is /// the coordinate volume, not the proper volume. This quantity diff --git a/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp b/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp index c18cf3834f0f..699e76bfcf47 100644 --- a/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp +++ b/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp @@ -155,11 +155,13 @@ void TemplatedLocalFunctions:: // We have all we need to calculate beta here. We take the largest of the // two variations get(beta)[i] = fabs(get(neutrino_energy_1)[i] - get(neutrino_energy_0)[i]) / - fabs(get(fluid_energy_1)[i] - get(fluid_energy_0)[i]); + (fabs(get(fluid_energy_1)[i] - get(fluid_energy_0)[i]) + + 1.e-16); const double beta_lepton = fabs(get(neutrino_lepton_number_1)[i] - get(neutrino_lepton_number_0)[i]) / - fabs(get(fluid_lepton_number_1)[i] - get(fluid_lepton_number_0)[i]); + (fabs(get(fluid_lepton_number_1)[i] - get(fluid_lepton_number_0)[i]) + + 1.e-16); get(beta)[i] = std::max(get(beta)[i], beta_lepton); } diff --git a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp index aaccea76beee..cbbbf8b9e251 100644 --- a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp +++ b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp @@ -309,9 +309,12 @@ void NeutrinoInteractionTable:: double ye = 0.0; double log_rho = 0.0; double log_temp = 0.0; + const double minimum_density = exp(table_log_density[0]); for (size_t p = 0; p < get(electron_fraction).size(); p++) { ye = get(electron_fraction)[p]; - log_rho = log(get(rest_mass_density)[p]); + log_rho = get(rest_mass_density)[p] > minimum_density + ? log(get(rest_mass_density)[p]) + : log(minimum_density); log_temp = get(temperature)[p] > minimum_temperature ? log(get(temperature)[p]) : log(minimum_temperature); @@ -351,6 +354,12 @@ void NeutrinoInteractionTable:: square(square(temperature_correction_factor)); gsl::at(gsl::at(*scattering_opacity, ns), ng)[p] *= square(square(temperature_correction_factor)); + gsl::at(gsl::at(*absorption_opacity, ns), ng)[p] = + std::clamp(gsl::at(gsl::at(*absorption_opacity, ns), ng)[p], + min_kappa, max_kappa); + gsl::at(gsl::at(*scattering_opacity, ns), ng)[p] = + std::clamp(gsl::at(gsl::at(*scattering_opacity, ns), ng)[p], + min_kappa, max_kappa); } } } diff --git a/src/Evolution/Particles/MonteCarlo/System.hpp b/src/Evolution/Particles/MonteCarlo/System.hpp index 9bf599363b6e..c638ed2b45b8 100644 --- a/src/Evolution/Particles/MonteCarlo/System.hpp +++ b/src/Evolution/Particles/MonteCarlo/System.hpp @@ -9,6 +9,7 @@ #include "Evolution/Particles/MonteCarlo/Tags.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" +#include "PointwiseFunctions/Hydro/TagsDeclarations.hpp" #include "Utilities/TMPL.hpp" namespace Particles::MonteCarlo { @@ -31,14 +32,8 @@ struct System { ::Tags::Variables>; using flux_spacetime_variables_tag = ::Tags::Variables>; // Hydro tags needed for background metric - using hydro_variables_tag = ::Tags::Variables, - hydro::Tags::RestMassDensity, - hydro::Tags::ElectronFraction, - hydro::Tags::Temperature, - hydro::Tags::LowerSpatialFourVelocity>>; - using primitive_variables_tag = - hydro_variables_tag; //::Tags::Variables>; + using hydro_variables_tag = ::Tags::Variables>; + using primitive_variables_tag = hydro_variables_tag; using inverse_spatial_metric_tag = gr::Tags::InverseSpatialMetric; diff --git a/src/Evolution/Particles/MonteCarlo/Tags.hpp b/src/Evolution/Particles/MonteCarlo/Tags.hpp index dc76c4b7d7e5..e0274fc1ae89 100644 --- a/src/Evolution/Particles/MonteCarlo/Tags.hpp +++ b/src/Evolution/Particles/MonteCarlo/Tags.hpp @@ -53,7 +53,7 @@ struct InteractionRatesTable : db::SimpleTag { std::unique_ptr>; static constexpr bool pass_metavariables = false; using option_tags = - NeutrinoInteractionTable::options; + typename NeutrinoInteractionTable::options; static type create_from_options(const std::string filename) { std::unique_ptr> diff --git a/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp b/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp index 00838bf3f48b..06c715cbd48f 100644 --- a/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp +++ b/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp @@ -98,7 +98,7 @@ void TemplatedLocalFunctions:: const tnsr::iJJ& d_inv_spatial_metric, const tnsr::ii& spatial_metric, const tnsr::II& inv_spatial_metric, - const Scalar& determinant_spatial_metric, + const Scalar& sqrt_determinant_spatial_metric, const Scalar& cell_light_crossing_time, const Mesh<3>& mesh, @@ -137,8 +137,8 @@ void TemplatedLocalFunctions:: Scalar cell_inertial_three_volume = make_with_value>(lapse, 0.0); cell_proper_four_volume_finite_difference( - &cell_proper_four_volume, lapse, determinant_spatial_metric, time_step, - mesh, det_jacobian_logical_to_inertial); + &cell_proper_four_volume, lapse, sqrt_determinant_spatial_metric, + time_step, mesh, det_jacobian_logical_to_inertial); cell_inertial_coordinate_three_volume_finite_difference( &cell_inertial_three_volume, mesh, det_jacobian_logical_to_inertial); diff --git a/src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp b/src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp index b88054158b99..e62fbc5b67a8 100644 --- a/src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp +++ b/src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp @@ -85,7 +85,7 @@ struct TemplatedLocalFunctions { const tnsr::iJJ& d_inv_spatial_metric, const tnsr::ii& spatial_metric, const tnsr::II& inv_spatial_metric, - const Scalar& determinant_spatial_metric, + const Scalar& sqrt_determinant_spatial_metric, const Scalar& cell_light_crossing_time, const Mesh<3>& mesh, const tnsr::I& mesh_coordinates, size_t num_ghost_zones, diff --git a/src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp b/src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp index ad29e279b095..7de83e9076a3 100644 --- a/src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp +++ b/src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp @@ -234,19 +234,32 @@ size_t MultiLinearSpanInterpolation< << which_dimension << "\ncurrent_index: " << current_index << "\nnumber of points: " << number_of_points_[which_dimension] << "\ntarget point: " << target_points - << "\nrelative coordinate: " << relative_coordinate); + << "\nrelative coordinate: " << relative_coordinate + << "\nspacing: " << inverse_spacing_[which_dimension]); + + // Do not trigger errors for roundoff errors in index calculation + if (UNLIKELY(current_index + 1 == number_of_points_[which_dimension])) { + if (relative_coordinate * inverse_spacing_[which_dimension] - + static_cast(current_index) < + 1.e-13) { + current_index -= 1; + } + } // We are exceeding the table bounds: // Use linear extrapolation based of the highest // two points in the table - ASSERT(allow_extrapolation_abov_data_[which_dimension] or - UNLIKELY(current_index + 1 < number_of_points_[which_dimension]), - "Interpolation exceeds upper table bounds.\nwhich_dimension: " - << which_dimension << "\ncurrent_index: " << current_index - << "\nnumber of points: " << number_of_points_[which_dimension] - << "\ntarget point: " << target_points - << "\nrelative coordinate: " << relative_coordinate); + ASSERT( + allow_extrapolation_abov_data_[which_dimension] or + UNLIKELY(current_index + 1 < number_of_points_[which_dimension]), + "Interpolation exceeds upper table bounds.\nwhich_dimension: " + << which_dimension << "\ncurrent_index: " << current_index + << "\nnumber of points: " << number_of_points_[which_dimension] + << "\ntarget point: " << target_points + << "\nrelative coordinate: " << relative_coordinate << "\nposition: " + << relative_coordinate * inverse_spacing_[which_dimension] - + static_cast(number_of_points_[which_dimension] - 1)); // Enforce index ranges current_index = std::min(number_of_points_[which_dimension] - 2, diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp index 7e440d0a20b6..4c8cc0adfb8b 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp @@ -17,7 +17,7 @@ SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloCellVolume", const size_t dv_size = 27; DataVector zero_dv(dv_size, 0.0); Scalar lapse{DataVector(dv_size, 1.2)}; - Scalar determinant_spatial_metric{DataVector(dv_size, 0.9)}; + Scalar sqrt_determinant_spatial_metric{DataVector(dv_size, 0.9)}; Scalar det_jacobian_logical_to_inertial{DataVector(dv_size, 1.1)}; const double time_step = 0.6; @@ -25,8 +25,8 @@ SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloCellVolume", Scalar expected_cell_proper_four_volume{ DataVector(dv_size, 8.0 / 27.0 * 1.2 * 0.9 * 0.6 * 1.1)}; Particles::MonteCarlo::cell_proper_four_volume_finite_difference( - &cell_proper_four_volume, lapse, determinant_spatial_metric, time_step, - mesh, det_jacobian_logical_to_inertial); + &cell_proper_four_volume, lapse, sqrt_determinant_spatial_metric, + time_step, mesh, det_jacobian_logical_to_inertial); Scalar cell_inertial_three_volume{DataVector(dv_size, 0.0)}; Scalar expected_cell_inertial_three_volume{ diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp index 66dcdcc14166..df00fd2d40cc 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp @@ -58,7 +58,7 @@ void test_flat_space_time_step(const bool skip) { spatial_metric.get(0, 0) = 1.0; spatial_metric.get(1, 1) = 1.0; spatial_metric.get(2, 2) = 1.0; - Scalar determinant_spatial_metric(dv_size, 1.0); + Scalar sqrt_determinant_spatial_metric(dv_size, 1.0); tnsr::I shift = make_with_value>(lapse, 0.0); tnsr::i d_lapse = @@ -244,7 +244,7 @@ void test_flat_space_time_step(const bool skip) { electron_fraction, baryon_density, temperature, lorentz_factor, lower_spatial_four_velocity, lapse, shift, d_lapse, d_shift, d_inv_spatial_metric, spatial_metric, inv_spatial_metric, - determinant_spatial_metric, cell_light_crossing_time, mesh, + sqrt_determinant_spatial_metric, cell_light_crossing_time, mesh, mesh_coordinates, num_ghost_zones, mesh_velocity, inverse_jacobian_logical_to_inertial, det_jacobian_logical_to_inertial, jacobian_inertial_to_fluid, inverse_jacobian_inertial_to_fluid, diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp index 08b66b4fdb8f..a6d3063ca2be 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp @@ -190,7 +190,7 @@ void test_advance_packets(const bool skip) { spatial_metric.get(0, 0) = 1.0; spatial_metric.get(1, 1) = 1.0; spatial_metric.get(2, 2) = 1.0; - Scalar determinant_spatial_metric(n_pts, 1.0); + Scalar sqrt_determinant_spatial_metric(n_pts, 1.0); tnsr::I shift = make_with_value>(lapse, 0.0); tnsr::i d_lapse = @@ -320,7 +320,7 @@ void test_advance_packets(const bool skip) { d_inv_spatial_metric, spatial_metric, inv_spatial_metric, - determinant_spatial_metric, + sqrt_determinant_spatial_metric, mesh_coordinates, mesh_velocity, inverse_jacobian_logical_to_inertial,