diff --git a/src/DataStructures/Variables.hpp b/src/DataStructures/Variables.hpp index 17f355366ded..f7e8d98a5a20 100644 --- a/src/DataStructures/Variables.hpp +++ b/src/DataStructures/Variables.hpp @@ -594,6 +594,8 @@ class Variables> { template Variables(const T* /*pointer*/, const size_t /*size*/) {} static constexpr size_t size() { return 0; } + void assign_subset(const Variables>& /*unused*/) {} + void assign_subset(const tuples::TaggedTuple<>& /*unused*/) {} }; // gcc8 screams when the empty Variables has pup as a member function, so we diff --git a/src/Evolution/DgSubcell/Actions/Initialize.hpp b/src/Evolution/DgSubcell/Actions/Initialize.hpp index 01cf8fb10c53..5ba41afe4ff9 100644 --- a/src/Evolution/DgSubcell/Actions/Initialize.hpp +++ b/src/Evolution/DgSubcell/Actions/Initialize.hpp @@ -46,6 +46,7 @@ #include "Utilities/CallWithDynamicType.hpp" #include "Utilities/ContainerHelpers.hpp" #include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/Overloader.hpp" #include "Utilities/TMPL.hpp" #include "Utilities/TaggedTuple.hpp" @@ -187,17 +188,10 @@ struct SetSubcellGrid { }, make_not_null(&box)); - db::mutate_apply< - tmpl::list, - tmpl::list<>>( - [&cell_is_not_on_external_boundary, &dg_mesh, - subcell_allowed_in_element, &subcell_mesh]( + const auto init_non_vars = + [&cell_is_not_on_external_boundary, subcell_allowed_in_element]( const gsl::not_null active_grid_ptr, const gsl::not_null did_rollback_ptr, - const auto active_vars_ptr, const gsl::not_null tci_decision_ptr, const gsl::not_null tci_calls_since_rollback_ptr, const gsl::not_null steps_since_tci_call_ptr) { @@ -210,17 +204,59 @@ struct SetSubcellGrid { subcell_enabled_at_external_boundary) and subcell_allowed_in_element) { *active_grid_ptr = ActiveGrid::Subcell; - active_vars_ptr->initialize(subcell_mesh.number_of_grid_points()); } else { *active_grid_ptr = ActiveGrid::Dg; - active_vars_ptr->initialize(dg_mesh.number_of_grid_points()); } *tci_decision_ptr = 0; *tci_calls_since_rollback_ptr = 0; *steps_since_tci_call_ptr = 0; + }; + + db::mutate_apply< + tmpl::flatten>, + tmpl::list<>, typename System::variables_tag>, + subcell::Tags::TciDecision, subcell::Tags::TciCallsSinceRollback, + subcell::Tags::StepsSinceTciCall>>, + tmpl::list<>>( + Overloader{ + [&init_non_vars]( + const gsl::not_null active_grid_ptr, + const gsl::not_null did_rollback_ptr, + const gsl::not_null tci_decision_ptr, + const gsl::not_null tci_calls_since_rollback_ptr, + const gsl::not_null steps_since_tci_call_ptr) { + init_non_vars(active_grid_ptr, did_rollback_ptr, tci_decision_ptr, + tci_calls_since_rollback_ptr, + steps_since_tci_call_ptr); + }, + [&init_non_vars, &dg_mesh, &subcell_mesh, + &cell_is_not_on_external_boundary, subcell_allowed_in_element]( + const gsl::not_null active_grid_ptr, + const gsl::not_null did_rollback_ptr, + const auto active_vars_ptr, + const gsl::not_null tci_decision_ptr, + const gsl::not_null tci_calls_since_rollback_ptr, + const gsl::not_null steps_since_tci_call_ptr) { + init_non_vars(active_grid_ptr, did_rollback_ptr, tci_decision_ptr, + tci_calls_since_rollback_ptr, + steps_since_tci_call_ptr); + if ((cell_is_not_on_external_boundary or + subcell_enabled_at_external_boundary) and + subcell_allowed_in_element) { + active_vars_ptr->initialize( + subcell_mesh.number_of_grid_points()); + } else { + active_vars_ptr->initialize(dg_mesh.number_of_grid_points()); + } + }, }, make_not_null(&box)); + if constexpr (System::has_primitive_and_conservative_vars) { db::mutate( [&dg_mesh, &subcell_mesh](const auto prim_vars_ptr, diff --git a/src/Evolution/DgSubcell/BackgroundGrVars.hpp b/src/Evolution/DgSubcell/BackgroundGrVars.hpp index 3863c5e6a6c4..6158c235fa1d 100644 --- a/src/Evolution/DgSubcell/BackgroundGrVars.hpp +++ b/src/Evolution/DgSubcell/BackgroundGrVars.hpp @@ -105,50 +105,62 @@ struct BackgroundGrVars : tt::ConformsTo { const bool did_rollback, const T& solution_or_data) { const size_t num_subcell_pts = subcell_mesh.number_of_grid_points(); - if (gsl::at(*subcell_face_gr_vars, 0).number_of_grid_points() != 0) { - // Evolution phase + bool evolution_phase = false; + if constexpr (not std::is_same_v< + typename SubcellFaceGrVars::value_type::tags_list, + tmpl::list<>>) { + evolution_phase = + (gsl::at(*subcell_face_gr_vars, 0).number_of_grid_points() != 0); + } - // Check if the mesh is actually moving i.e. block coordinate map is - // time-dependent. If not, we can skip the evaluation of GR variables - // since they may stay with their values assigned at the initialization - // phase. - const auto& element_id = element.id(); - const size_t block_id = element_id.block_id(); - const Block& block = domain.blocks()[block_id]; + if (evolution_phase) { + if constexpr (not std::is_same_v< + typename SubcellFaceGrVars::value_type::tags_list, + tmpl::list<>>) { + // Evolution phase - if (block.is_time_dependent()) { - if (did_rollback or not ComputeOnlyOnRollback) { - if (did_rollback) { - // Right after rollback, subcell GR vars are stored in the - // `inactive` one. - ASSERT(inactive_gr_vars->number_of_grid_points() == num_subcell_pts, - "The size of subcell GR variables (" - << inactive_gr_vars->number_of_grid_points() - << ") is not equal to the number of FD grid points (" - << subcell_mesh.number_of_grid_points() << ")."); + // Check if the mesh is actually moving i.e. block coordinate map is + // time-dependent. If not, we can skip the evaluation of GR variables + // since they may stay with their values assigned at the initialization + // phase. + const auto& element_id = element.id(); + const size_t block_id = element_id.block_id(); + const Block& block = domain.blocks()[block_id]; - cell_centered_impl(inactive_gr_vars, time, subcell_inertial_coords, - solution_or_data); + if (block.is_time_dependent()) { + if (did_rollback or not ComputeOnlyOnRollback) { + if (did_rollback) { + // Right after rollback, subcell GR vars are stored in the + // `inactive` one. + ASSERT( + inactive_gr_vars->number_of_grid_points() == num_subcell_pts, + "The size of subcell GR variables (" + << inactive_gr_vars->number_of_grid_points() + << ") is not equal to the number of FD grid points (" + << subcell_mesh.number_of_grid_points() << ")."); - } else { - // In this case the element didn't rollback but started from FD. - // Therefore subcell GR vars are in the `active` one. - ASSERT(active_gr_vars->number_of_grid_points() == num_subcell_pts, - "The size of subcell GR variables (" - << active_gr_vars->number_of_grid_points() - << ") is not equal to the number of FD grid points (" - << subcell_mesh.number_of_grid_points() << ")."); + cell_centered_impl(inactive_gr_vars, time, + subcell_inertial_coords, solution_or_data); - cell_centered_impl(active_gr_vars, time, subcell_inertial_coords, - solution_or_data); - } + } else { + // In this case the element didn't rollback but started from FD. + // Therefore subcell GR vars are in the `active` one. + ASSERT(active_gr_vars->number_of_grid_points() == num_subcell_pts, + "The size of subcell GR variables (" + << active_gr_vars->number_of_grid_points() + << ") is not equal to the number of FD grid points (" + << subcell_mesh.number_of_grid_points() << ")."); - face_centered_impl(subcell_face_gr_vars, time, functions_of_time, - logical_to_grid_map, grid_to_inertial_map, - subcell_mesh, solution_or_data); + cell_centered_impl(active_gr_vars, time, subcell_inertial_coords, + solution_or_data); + } + + face_centered_impl(subcell_face_gr_vars, time, functions_of_time, + logical_to_grid_map, grid_to_inertial_map, + subcell_mesh, solution_or_data); + } } } - } else { // Initialization phase (*inactive_gr_vars).initialize(num_subcell_pts); @@ -158,19 +170,22 @@ struct BackgroundGrVars : tt::ConformsTo { "The subcell mesh must have isotropic basis, quadrature. and " "extents but got " << subcell_mesh); - const size_t num_face_centered_mesh_grid_pts = - (subcell_mesh.extents(0) + 1) * subcell_mesh.extents(1) * - subcell_mesh.extents(2); - for (size_t d = 0; d < volume_dim; ++d) { - gsl::at(*subcell_face_gr_vars, d) - .initialize(num_face_centered_mesh_grid_pts); + if constexpr (not std::is_same_v< + typename SubcellFaceGrVars::value_type::tags_list, + tmpl::list<>>) { + const size_t num_face_centered_mesh_grid_pts = + (subcell_mesh.extents(0) + 1) * subcell_mesh.extents(1) * + subcell_mesh.extents(2); + for (size_t d = 0; d < volume_dim; ++d) { + gsl::at(*subcell_face_gr_vars, d) + .initialize(num_face_centered_mesh_grid_pts); + } + face_centered_impl(subcell_face_gr_vars, time, functions_of_time, + logical_to_grid_map, grid_to_inertial_map, + subcell_mesh, solution_or_data); } - cell_centered_impl(inactive_gr_vars, time, subcell_inertial_coords, solution_or_data); - face_centered_impl(subcell_face_gr_vars, time, functions_of_time, - logical_to_grid_map, grid_to_inertial_map, - subcell_mesh, solution_or_data); } } diff --git a/src/Evolution/Executables/RadiationTransport/CMakeLists.txt b/src/Evolution/Executables/RadiationTransport/CMakeLists.txt index 9a5cb7938cb7..060886a54225 100644 --- a/src/Evolution/Executables/RadiationTransport/CMakeLists.txt +++ b/src/Evolution/Executables/RadiationTransport/CMakeLists.txt @@ -2,3 +2,4 @@ # See LICENSE.txt for details. add_subdirectory(M1Grey) +add_subdirectory(MonteCarlo) diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt b/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt new file mode 100644 index 000000000000..8680009bc4e2 --- /dev/null +++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt @@ -0,0 +1,45 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +set(EXECUTABLE EvolveMonteCarlo) + +add_spectre_executable( + ${EXECUTABLE} + EXCLUDE_FROM_ALL + EvolveMonteCarlo.cpp + ) + +target_link_libraries( + ${EXECUTABLE} + PRIVATE + Actions + Charmxx::main + CoordinateMaps + DataStructures + DgSubcell + DiscontinuousGalerkin + DomainCreators + Events + EventsAndDenseTriggers + EventsAndTriggers + Evolution + GeneralRelativity + GeneralRelativitySolutions + GrMhdAnalyticData + GrMhdSolutions + H5 + Hydro + Informer + LinearOperators + MathFunctions + MonteCarlo + Observer + Options + Parallel + PhaseControl + Serialization + Time + Utilities + ValenciaDivClean + ) + diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp new file mode 100644 index 000000000000..db4a08f04d8b --- /dev/null +++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp @@ -0,0 +1,33 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp" + +#include + +#include "Domain/Creators/RegisterDerivedWithCharm.hpp" +#include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp" +#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" +#include "Parallel/CharmMain.tpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/RegisterDerivedWithCharm.hpp" +#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" + +void register_neutrino_tables() { + register_classes_with_charm( + tmpl::list, + Particles::MonteCarlo::NeutrinoInteractionTable<2, 3>, + Particles::MonteCarlo::NeutrinoInteractionTable<4, 3>, + Particles::MonteCarlo::NeutrinoInteractionTable<16, 3>>{}); +} + +extern "C" void CkRegisterMainModule() { + Parallel::charmxx::register_main_module>(); + Parallel::charmxx::register_init_node_and_proc( + {&domain::creators::register_derived_with_charm, + &domain::creators::time_dependence::register_derived_with_charm, + &domain::FunctionsOfTime::register_derived_with_charm, + &EquationsOfState::register_derived_with_charm, + ®ister_factory_classes_with_charm>, + ®ister_neutrino_tables}, + {}); +} diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp new file mode 100644 index 000000000000..8f90a3c79674 --- /dev/null +++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp @@ -0,0 +1,254 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include + +#include "Domain/Creators/Factory3D.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/Actions/RunEventsAndDenseTriggers.hpp" +#include "Evolution/Actions/RunEventsAndTriggers.hpp" +#include "Evolution/ComputeTags.hpp" +#include "Evolution/DgSubcell/Actions/Initialize.hpp" +#include "Evolution/DgSubcell/BackgroundGrVars.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp" +#include "Evolution/DiscontinuousGalerkin/BackgroundGrVars.hpp" +#include "Evolution/DiscontinuousGalerkin/DgElementArray.hpp" +#include "Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp" +#include "Evolution/DiscontinuousGalerkin/Initialization/QuadratureTag.hpp" +#include "Evolution/DiscontinuousGalerkin/Limiters/Minmod.hpp" +#include "Evolution/DiscontinuousGalerkin/Limiters/Tags.hpp" +#include "Evolution/Initialization/ConservativeSystem.hpp" +#include "Evolution/Initialization/DgDomain.hpp" +#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" +#include "Evolution/Particles/MonteCarlo/System.hpp" +#include "Evolution/Particles/MonteCarlo/Tags.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/AllSolutions.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/SwapGrTags.hpp" +#include "IO/Observer/Actions/RegisterEvents.hpp" +#include "IO/Observer/Helpers.hpp" +#include "IO/Observer/ObserverComponent.hpp" +#include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp" +#include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp" +#include "Options/Protocols/FactoryCreation.hpp" +#include "Options/String.hpp" +#include "Parallel/Local.hpp" +#include "Parallel/Phase.hpp" +#include "Parallel/PhaseControl/CheckpointAndExitAfterWallclock.hpp" +#include "Parallel/PhaseControl/ExecutePhaseChange.hpp" +#include "Parallel/PhaseControl/Factory.hpp" +#include "Parallel/PhaseControl/VisitAndReturn.hpp" +#include "Parallel/PhaseDependentActionList.hpp" +#include "Parallel/Protocols/RegistrationMetavariables.hpp" +#include "ParallelAlgorithms/Actions/AddComputeTags.hpp" +#include "ParallelAlgorithms/Actions/AddSimpleTags.hpp" +#include "ParallelAlgorithms/Actions/InitializeItems.hpp" +#include "ParallelAlgorithms/Actions/LimiterActions.hpp" +#include "ParallelAlgorithms/Actions/MutateApply.hpp" +#include "ParallelAlgorithms/Actions/TerminatePhase.hpp" +#include "ParallelAlgorithms/Events/Factory.hpp" +#include "ParallelAlgorithms/EventsAndDenseTriggers/DenseTrigger.hpp" +#include "ParallelAlgorithms/EventsAndDenseTriggers/DenseTriggers/Factory.hpp" +#include "ParallelAlgorithms/EventsAndTriggers/Completion.hpp" +#include "ParallelAlgorithms/EventsAndTriggers/Event.hpp" +#include "ParallelAlgorithms/EventsAndTriggers/EventsAndTriggers.hpp" +#include "ParallelAlgorithms/EventsAndTriggers/LogicalTriggers.hpp" +#include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp" +#include "PointwiseFunctions/AnalyticData/AnalyticData.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp" +#include "PointwiseFunctions/AnalyticData/Tags.hpp" +#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" +#include "Time/Actions/SelfStartActions.hpp" +#include "Time/Actions/UpdateU.hpp" +#include "Time/ChangeSlabSize/Action.hpp" +#include "Time/StepChoosers/Factory.hpp" +#include "Time/StepChoosers/StepChooser.hpp" +#include "Time/Tags/Time.hpp" +#include "Time/Tags/TimeStepId.hpp" +#include "Time/TimeSequence.hpp" +#include "Time/TimeSteppers/Factory.hpp" +#include "Time/TimeSteppers/LtsTimeStepper.hpp" +#include "Time/TimeSteppers/TimeStepper.hpp" +#include "Time/Triggers/TimeTriggers.hpp" +#include "Utilities/Functional.hpp" +#include "Utilities/ProtocolHelpers.hpp" +#include "Utilities/TMPL.hpp" + +/// \cond +namespace Frame { +struct Inertial; +} // namespace Frame +namespace PUP { +class er; +} // namespace PUP +namespace Parallel { +template +class CProxy_GlobalCache; +} // namespace Parallel +/// \endcond + +// NEED: +// Initial data + +template +struct EvolutionMetavars { + static constexpr size_t volume_dim = 3; + + using system = Particles::MonteCarlo::System; + using temporal_id = Tags::TimeStepId; + using TimeStepperBase = TimeStepper; + static constexpr bool use_dg_subcell = true; + + using initial_data_list = + grmhd::ValenciaDivClean::InitialData::initial_data_list; + using equation_of_state_tag = hydro::Tags::GrmhdEquationOfState; + + struct SubcellOptions { + using evolved_vars_tags = typename system::variables_tag::tags_list; + + static constexpr bool subcell_enabled = use_dg_subcell; + static constexpr bool subcell_enabled_at_external_boundary = true; + }; + + using observe_fields = + tmpl::list, + domain::Tags::Coordinates>; + using non_tensor_compute_tags = + tmpl::list<::Events::Tags::ObserverMeshCompute, + ::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< + tmpl::pair, + tmpl::pair, domain_creators>, + tmpl::pair>>, + tmpl::pair, + tmpl::pair< + grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField, + grmhd::AnalyticData::InitialMagneticFields:: + initial_magnetic_fields>, + tmpl::pair>, + tmpl::pair, + StepChoosers::standard_slab_choosers>, + tmpl::pair, + TimeSequences::all_time_sequences>, + tmpl::pair, + TimeSequences::all_time_sequences>, + tmpl::pair, + tmpl::pair>; + }; + + using observed_reduction_data_tags = + observers::collect_reduction_data_tags>>>; + + using dg_registration_list = + tmpl::list; + + using initialization_actions = tmpl::list< + Initialization::Actions::InitializeItems< + Initialization::TimeStepping, + evolution::dg::Initialization::Domain //, + >, + Initialization::Actions::AddSimpleTags< + evolution::dg::BackgroundGrVars>, + evolution::dg::subcell::Actions::SetSubcellGrid, + Initialization::Actions::AddSimpleTags< + evolution::dg::subcell::BackgroundGrVars>, + Actions::MutateApply, + Initialization::Actions::InitializeMCTags, + Initialization::Actions::AddComputeTags>>, + evolution::Actions::InitializeRunEventsAndDenseTriggers, + Parallel::Actions::TerminatePhase>; + + using dg_element_array = DgElementArray< + EvolutionMetavars, + tmpl::list< + Parallel::PhaseActions, + Parallel::PhaseActions< + Parallel::Phase::Evolve, + tmpl::list< + Actions::AdvanceTime, + evolution::Actions::RunEventsAndTriggers, + Actions::AdvanceTime, + evolution::Actions::RunEventsAndTriggers, + Actions::AdvanceTime, + evolution::Actions::RunEventsAndTriggers, + Particles::MonteCarlo::Actions::SendDataForMcCommunication< + volume_dim, + // No local time stepping + false, Particles::MonteCarlo::CommunicationStep::PreStep>, + Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication< + volume_dim, + Particles::MonteCarlo::CommunicationStep::PreStep>, + Particles::MonteCarlo::Actions::TakeTimeStep, + Particles::MonteCarlo::Actions::SendDataForMcCommunication< + volume_dim, + // No local time stepping + false, + Particles::MonteCarlo::CommunicationStep::PostStep>, + Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication< + volume_dim, + Particles::MonteCarlo::CommunicationStep::PostStep>, + PhaseControl::Actions::ExecutePhaseChange>>>>; + + struct registration + : tt::ConformsTo { + using element_registrars = + tmpl::map>; + }; + + using component_list = + tmpl::list, + observers::ObserverWriter, + dg_element_array>; + + 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"}; + + static constexpr std::array default_phase_order{ + {Parallel::Phase::Initialization, + Parallel::Phase::InitializeTimeStepperHistory, Parallel::Phase::Register, + Parallel::Phase::Evolve, Parallel::Phase::Exit}}; + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& /*p*/) {} +}; diff --git a/src/Evolution/Initialization/SetVariables.hpp b/src/Evolution/Initialization/SetVariables.hpp index f8a28f5930c1..53eef9d8ee31 100644 --- a/src/Evolution/Initialization/SetVariables.hpp +++ b/src/Evolution/Initialization/SetVariables.hpp @@ -155,15 +155,18 @@ struct SetVariables { using variables_tag = typename system::variables_tag; // Set initial data from analytic solution - using Vars = typename variables_tag::type; - db::mutate( - [&initial_time, &inertial_coords, - &solution_or_data](const gsl::not_null vars) { - vars->assign_subset(evolution::Initialization::initial_data( - solution_or_data, inertial_coords, initial_time, - typename Vars::tags_list{})); - }, - box); + if constexpr (not std::is_same_v>) { + using Vars = typename variables_tag::type; + db::mutate( + [&initial_time, &inertial_coords, + &solution_or_data](const gsl::not_null vars) { + vars->assign_subset(evolution::Initialization::initial_data( + solution_or_data, inertial_coords, initial_time, + typename Vars::tags_list{})); + }, + box); + } } } }; 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..f5ecaf298f3b --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp @@ -0,0 +1,163 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "Domain/Structure/Element.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/Initialization/InitialData.hpp" +#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.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 +/// * Particles::MonteCarlo::Tags::McGhostZoneDataTag +/// +/// - 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, + Particles::MonteCarlo::Tags::McGhostZoneDataTag>; + + 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)); + + unsigned long seed = + std::random_device{}(); // static_cast(time(NULL)); + typename Particles::MonteCarlo::Tags::RandomNumberGenerator::type rng(seed); + + 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)); + + using GhostZoneData = + typename Particles::MonteCarlo::Tags::McGhostZoneDataTag::type; + GhostZoneData ghost_zone_data{}; + Initialization::mutate_assign< + tmpl::list>>( + make_not_null(&box), std::move(ghost_zone_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..42bf52adaad2 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>& @@ -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/CMakeLists.txt b/src/Evolution/Particles/MonteCarlo/CMakeLists.txt index 0f8c3f106098..daca3ec8939d 100644 --- a/src/Evolution/Particles/MonteCarlo/CMakeLists.txt +++ b/src/Evolution/Particles/MonteCarlo/CMakeLists.txt @@ -38,6 +38,7 @@ spectre_target_headers( NeutrinoInteractionTable.hpp Packet.hpp Scattering.hpp + System.hpp Tags.hpp TakeTimeStep.tpp TemplatedLocalFunctions.hpp 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/EvolvePacketsInElement.tpp b/src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp index 3abfe9d8d072..effa473212c1 100644 --- a/src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp +++ b/src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp @@ -162,6 +162,22 @@ void TemplatedLocalFunctions::evolve_packets( initial_time = packet.time; dt_end_step = final_time - initial_time; + // Find closest grid point to packet at current time, using + // extents for live points only. + { + 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] - gsl::at(bottom_coord_mesh, d)) / + gsl::at(dx_mesh, d) + + 0.5); + } + 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]); + } + // Get quantities that we do NOT update if the packet // changes cell. // local_idx is the index on the mesh without ghost zones diff --git a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp index 2f9422a5848a..d11ba5cedcc3 100644 --- a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp +++ b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp @@ -311,8 +311,8 @@ struct ReceiveDataForMcCommunication { const DataVector& received_data_direction = received_data[directional_element_id] .ghost_zone_hydro_variables; - REQUIRE(received_data[directional_element_id] - .packets_entering_this_element == std::nullopt); + // REQUIRE(received_data[directional_element_id] + // .packets_entering_this_element == std::nullopt); if (mortar_data->rest_mass_density[directional_element_id] == std::nullopt) { continue; @@ -368,8 +368,8 @@ struct ReceiveDataForMcCommunication { received_data[directional_element_id] .packets_entering_this_element; // Temporary: currently no data for coupling to the fluid - REQUIRE(received_data[directional_element_id] - .ghost_zone_hydro_variables.size() == 0); + // REQUIRE(received_data[directional_element_id] + // .ghost_zone_hydro_variables.size() == 0); if (received_data_packets == std::nullopt) { continue; } else { 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/NeutrinoInteractionTable.hpp b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp index 69c057c841ae..db027ad2567f 100644 --- a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp +++ b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp @@ -11,6 +11,7 @@ #include "DataStructures/BoostMultiArray.hpp" #include "DataStructures/Tensor/TypeAliases.hpp" #include "NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp" +#include "Options/String.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/Serialization/CharmPupable.hpp" @@ -28,6 +29,19 @@ class NeutrinoInteractionTable : public PUP::able { /// Read table from disk and stores interaction rates. explicit NeutrinoInteractionTable(const std::string& filename); + static constexpr Options::String help = { + "Tabulated neutrino-matter interactions in NuLib format.\n" + "Emissivity, absorption and scattering opacities are \n" + "tabulated as a function of density, eletron fraction \n" + "and temperature."}; + + struct TableFilename { + using type = std::string; + static constexpr Options::String help{"File name of the NuLib table"}; + }; + + using options = tmpl::list; + /// Explicit instantiation from table values, for tests NeutrinoInteractionTable( std::vector table_data_, diff --git a/src/Evolution/Particles/MonteCarlo/System.hpp b/src/Evolution/Particles/MonteCarlo/System.hpp new file mode 100644 index 000000000000..c638ed2b45b8 --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/System.hpp @@ -0,0 +1,41 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/VariablesTag.hpp" +#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 { + +struct System { + static constexpr bool is_in_flux_conservative_form = false; + static constexpr bool has_primitive_and_conservative_vars = false; + static constexpr size_t volume_dim = 3; + // The EoS is used within the MC code itself, even if we provide + // the fluid variables as a background. + static constexpr size_t thermodynamic_dim = 3; + + using mc_variables_tag = ::Tags::Variables< + tmpl::list>; + using variables_tag = ::Tags::Variables>; // mc_variables_tag; + using flux_variables = tmpl::list<>; + using gradient_variables = tmpl::list<>; + // GR tags needed for background metric + using spacetime_variables_tag = + ::Tags::Variables>; + using flux_spacetime_variables_tag = ::Tags::Variables>; + // Hydro tags needed for background metric + using hydro_variables_tag = ::Tags::Variables>; + using primitive_variables_tag = hydro_variables_tag; + + using inverse_spatial_metric_tag = + gr::Tags::InverseSpatialMetric; +}; +} // namespace Particles::MonteCarlo diff --git a/src/Evolution/Particles/MonteCarlo/Tags.hpp b/src/Evolution/Particles/MonteCarlo/Tags.hpp index ed0fb3080df8..e0274fc1ae89 100644 --- a/src/Evolution/Particles/MonteCarlo/Tags.hpp +++ b/src/Evolution/Particles/MonteCarlo/Tags.hpp @@ -51,6 +51,18 @@ template struct InteractionRatesTable : db::SimpleTag { using type = std::unique_ptr>; + static constexpr bool pass_metavariables = false; + using option_tags = + typename NeutrinoInteractionTable::options; + static type create_from_options(const std::string filename) { + std::unique_ptr> + interaction_table_ptr = + std::make_unique>(filename); + return interaction_table_ptr; + ; + } }; } // namespace Particles::MonteCarlo::Tags 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/src/ParallelAlgorithms/Events/ObserveTimeStep.hpp b/src/ParallelAlgorithms/Events/ObserveTimeStep.hpp index 9b727f56ca56..63751cb189d0 100644 --- a/src/ParallelAlgorithms/Events/ObserveTimeStep.hpp +++ b/src/ParallelAlgorithms/Events/ObserveTimeStep.hpp @@ -165,17 +165,18 @@ class ObserveTimeStep : public Event { // We obtain the grid size from the variables, rather than the mesh, // so that this observer is not DG-specific. using argument_tags = - tmpl::list<::Tags::TimeStep, typename System::variables_tag>; + tmpl::list<::Tags::TimeStep>; //, typename System::variables_tag>; template void operator()(const TimeDelta& time_step, - const typename System::variables_tag::type& variables, + // const typename System::variables_tag::type& variables, Parallel::GlobalCache& cache, const ArrayIndex& array_index, const ParallelComponent* const /*meta*/, const ObservationValue& observation_value) const { - const size_t number_of_grid_points = variables.number_of_grid_points(); + const size_t number_of_grid_points = + 1; // variables.number_of_grid_points(); const double slab_size = time_step.slab().duration().value(); const double step_size = abs(time_step.value()); const double wall_time = sys::wall_time(); 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,