From c6e137f1d93ca7d873bb74af727a117b77bcf65c Mon Sep 17 00:00:00 2001 From: Francois Foucart Date: Fri, 20 Dec 2024 14:58:12 -0500 Subject: [PATCH 1/3] Add functionality needed to have empty vars for MC system --- src/DataStructures/Variables.hpp | 5 +++ .../DgSubcell/Actions/Initialize.hpp | 3 +- src/Evolution/DgSubcell/BackgroundGrVars.hpp | 33 +++++++++++-------- .../Events/ObserveTimeStep.hpp | 8 ++++- tests/Unit/DataStructures/Test_Variables.cpp | 6 ++++ 5 files changed, 40 insertions(+), 15 deletions(-) diff --git a/src/DataStructures/Variables.hpp b/src/DataStructures/Variables.hpp index 17f355366ded..722126c2d096 100644 --- a/src/DataStructures/Variables.hpp +++ b/src/DataStructures/Variables.hpp @@ -594,6 +594,11 @@ class Variables> { template Variables(const T* /*pointer*/, const size_t /*size*/) {} static constexpr size_t size() { return 0; } + static constexpr size_t number_of_grid_points() { return 0; } + void assign_subset(const Variables>& /*unused*/) {} + void assign_subset(const tuples::TaggedTuple<>& /*unused*/) {} + // Initialization for empty variables should ignore any input. + void initialize(size_t /*number_of_grid_points*/) {} }; // 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..81bfb48edb72 100644 --- a/src/Evolution/DgSubcell/Actions/Initialize.hpp +++ b/src/Evolution/DgSubcell/Actions/Initialize.hpp @@ -103,7 +103,8 @@ struct SetSubcellGrid { subcell::Tags::ReconstructionOrder, evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd, evolution::dg::subcell::Tags::InterpolatorsFromDgToNeighborFd, - evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd>; + evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd, + typename System::variables_tag>; using compute_tags = tmpl::list, Tags::LogicalCoordinatesCompute, ::domain::Tags::MappedCoordinates< diff --git a/src/Evolution/DgSubcell/BackgroundGrVars.hpp b/src/Evolution/DgSubcell/BackgroundGrVars.hpp index 3863c5e6a6c4..abae1b80d90a 100644 --- a/src/Evolution/DgSubcell/BackgroundGrVars.hpp +++ b/src/Evolution/DgSubcell/BackgroundGrVars.hpp @@ -142,10 +142,13 @@ struct BackgroundGrVars : tt::ConformsTo { 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); + if constexpr (not std::is_same_v< + typename SubcellFaceGrVars::value_type::tags_list, + tmpl::list<>>){ + face_centered_impl(subcell_face_gr_vars, time, functions_of_time, + logical_to_grid_map, grid_to_inertial_map, + subcell_mesh, solution_or_data); + } } } @@ -158,19 +161,23 @@ 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/ParallelAlgorithms/Events/ObserveTimeStep.hpp b/src/ParallelAlgorithms/Events/ObserveTimeStep.hpp index 9b727f56ca56..ee7797de4c2a 100644 --- a/src/ParallelAlgorithms/Events/ObserveTimeStep.hpp +++ b/src/ParallelAlgorithms/Events/ObserveTimeStep.hpp @@ -175,7 +175,13 @@ class ObserveTimeStep : public Event { 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(); + // For empty vars, we use 1 grid point to avoid divisions by zero + // after reduction. + size_t number_of_grid_points = 1; + if constexpr (not std::is_same_v>) { + number_of_grid_points = 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/DataStructures/Test_Variables.cpp b/tests/Unit/DataStructures/Test_Variables.cpp index 397700c1bcd3..14a32cc3afb2 100644 --- a/tests/Unit/DataStructures/Test_Variables.cpp +++ b/tests/Unit/DataStructures/Test_Variables.cpp @@ -78,6 +78,12 @@ void test_empty_variables() { CHECK(get_output(tuple_with_empty_vars) == "(3,hello,{})"); CHECK(not contains_allocations(empty_vars)); + // Functions do nothing, but need to compile + empty_vars.initialize(0); + const Variables> input_vars; + empty_vars.assign_subset(input_vars); + CHECK(empty_vars.number_of_grid_points() == 0); + CHECK(empty_vars.size() == 0); } template From f6003921a5be576896ed103a59a3d38bde253a3f Mon Sep 17 00:00:00 2001 From: Francois Foucart Date: Tue, 29 Oct 2024 16:57:47 -0400 Subject: [PATCH 2/3] Base executable for MC --- .../RadiationTransport/CMakeLists.txt | 1 + .../MonteCarlo/CMakeLists.txt | 45 ++++ .../MonteCarlo/EvolveMonteCarlo.cpp | 33 +++ .../MonteCarlo/EvolveMonteCarlo.hpp | 254 ++++++++++++++++++ src/Evolution/Initialization/SetVariables.hpp | 21 +- .../MonteCarlo/Actions/CMakeLists.txt | 1 + .../Actions/InitializeMonteCarlo.hpp | 167 ++++++++++++ .../MonteCarlo/Actions/TimeStepActions.hpp | 31 ++- .../Particles/MonteCarlo/CMakeLists.txt | 1 + .../Particles/MonteCarlo/CellVolume.cpp | 4 +- .../Particles/MonteCarlo/CellVolume.hpp | 11 +- .../MonteCarlo/EvolvePacketsInElement.tpp | 16 ++ .../MonteCarlo/GhostZoneCommunication.hpp | 18 +- .../MonteCarlo/NeutrinoInteractionTable.cpp | 5 +- .../MonteCarlo/NeutrinoInteractionTable.hpp | 14 + src/Evolution/Particles/MonteCarlo/System.hpp | 41 +++ src/Evolution/Particles/MonteCarlo/Tags.hpp | 12 + .../Particles/MonteCarlo/TakeTimeStep.tpp | 6 +- .../MonteCarlo/TemplatedLocalFunctions.hpp | 2 +- .../MultiLinearSpanInterpolation.hpp | 29 +- .../Particles/MonteCarlo/Test_CellVolume.cpp | 14 +- .../MonteCarlo/Test_CommunicationTags.cpp | 4 + .../MonteCarlo/Test_TakeTimeStep.cpp | 24 +- .../MonteCarlo/Test_TimeStepAction.cpp | 14 +- 24 files changed, 701 insertions(+), 67 deletions(-) create mode 100644 src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt create mode 100644 src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp create mode 100644 src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp create mode 100644 src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp create mode 100644 src/Evolution/Particles/MonteCarlo/System.hpp 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..e1bdde9784fd --- /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..d36e4b2dda05 --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp @@ -0,0 +1,167 @@ +// 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/DgSubcell/Tags/Coordinates.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.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 +/// - evolution::dg::subcell::Tags::Coordinates +/// - evolution::initial_data::Tags::InitialData +/// +/// DataBox changes: +/// - Adds: +/// * Particles::MonteCarlo::Tags::PacketsOnElement +/// * Particles::MonteCarlo::Tags::RandomNumberGenerator +/// * Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission< +/// NeutrinoSpecies> +/// * 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 = 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)); + + const 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..0ced96571ea5 100644 --- a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp +++ b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp @@ -136,6 +136,16 @@ struct GhostDataMcPackets { const size_t d = max_distance_direction.value().dimension(); const Side side = max_distance_direction.value().side(); packet.coordinates.get(d) += (side == Side::Lower) ? (2.0) : (-2.0); + // Corner/edge treatment; move packet in a live point for now + // Definitely needs improvement... + for (size_t dd = 0; dd < Dim; dd++) { + if (dd != d && packet.coordinates.get(dd) < -1.0) { + packet.coordinates.get(dd) = -2.0 - packet.coordinates.get(dd); + } + if (dd != d && packet.coordinates.get(dd) > 1.0) { + packet.coordinates.get(dd) = 2.0 - packet.coordinates.get(dd); + } + } if (output.contains(max_distance_direction.value())) { output[max_distance_direction.value()].push_back(packet); } @@ -311,8 +321,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 +378,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/NeutrinoInteractionTable.cpp b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp index cbbbf8b9e251..e5894408c2bf 100644 --- a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp +++ b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp @@ -306,6 +306,7 @@ void NeutrinoInteractionTable:: const size_t n_ye_points = table_electron_fraction.size(); double temperature_correction_factor = 0.0; + const double dye = table_electron_fraction[1] - table_electron_fraction[0]; double ye = 0.0; double log_rho = 0.0; double log_temp = 0.0; @@ -323,8 +324,8 @@ void NeutrinoInteractionTable:: ? 1.0 : get(temperature)[p] / minimum_temperature; - ye = std::clamp(ye, table_electron_fraction[0], - table_electron_fraction[n_ye_points - 1]); + ye = std::clamp(ye, table_electron_fraction[0] + 1.e-8 * dye, + table_electron_fraction[n_ye_points - 1] - 1.e-8 * dye); log_rho = std::clamp(log_rho, table_log_density[0], table_log_density[n_rho_points - 1]); log_temp = std::clamp(log_temp, table_log_temperature[0], 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/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp index 7e440d0a20b6..61ac032d2d0b 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp @@ -15,18 +15,20 @@ SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloCellVolume", Spectral::Quadrature::CellCentered); 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 det_jacobian_logical_to_inertial{DataVector(dv_size, 1.1)}; + const DataVector zero_dv(dv_size, 0.0); + const Scalar lapse{DataVector(dv_size, 1.2)}; + const Scalar sqrt_determinant_spatial_metric + {DataVector(dv_size, 0.9)}; + const Scalar det_jacobian_logical_to_inertial + {DataVector(dv_size, 1.1)}; const double time_step = 0.6; Scalar cell_proper_four_volume{DataVector(dv_size, 0.0)}; 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_CommunicationTags.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp index a91b6f4c8f5e..46aaa87fdf8a 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp @@ -524,6 +524,10 @@ void test_send_receive_actions() { // topological coordinate of their new element. packet_east.coordinates[0] -= 2.0; packet_south.coordinates[1] += 2.0; + // Current hack for edges/corners + if constexpr (Dim > 2){ + packet_south.coordinates[2] = -2.0 - packet_south.coordinates[2]; + } { const auto& east_data = ActionTesting::get_inbox_tag< comp, Particles::MonteCarlo::McGhostZoneDataInboxTag< diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp index f15870ef41f6..15ebadbc052f 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp @@ -33,15 +33,15 @@ void test_flat_space_time_step() { const size_t num_ghost_zones = 1; const size_t dv_size = cube(size_1d); - DataVector zero_dv(dv_size, 0.0); + const DataVector zero_dv(dv_size, 0.0); const size_t dv_size_with_ghost = cube(size_1d + 2 * num_ghost_zones); - DataVector zero_dv_with_ghost(dv_size_with_ghost, 0.0); + const DataVector zero_dv_with_ghost(dv_size_with_ghost, 0.0); const size_t dv_size_in_ghost = square(size_1d) * num_ghost_zones; - DataVector zero_dv_in_ghost(dv_size_in_ghost, 0.0); - DataVector one_dv_in_ghost(dv_size_in_ghost, 1.0); + const DataVector zero_dv_in_ghost(dv_size_in_ghost, 0.0); + const DataVector one_dv_in_ghost(dv_size_in_ghost, 1.0); // Minkowski metric - Scalar lapse{DataVector(dv_size, 1.0)}; + const Scalar lapse{DataVector(dv_size, 1.0)}; tnsr::II inv_spatial_metric = make_with_value>(lapse, 0.0); inv_spatial_metric.get(0, 0) = 1.0; @@ -52,14 +52,14 @@ void test_flat_space_time_step() { 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); - tnsr::I shift = + const Scalar sqrt_determinant_spatial_metric(dv_size, 1.0); + const tnsr::I shift = make_with_value>(lapse, 0.0); - tnsr::i d_lapse = + const tnsr::i d_lapse = make_with_value>(lapse, 0.0); - tnsr::iJ d_shift = + const tnsr::iJ d_shift = make_with_value>(lapse, 0.0); - tnsr::iJJ d_inv_spatial_metric = + const tnsr::iJJ d_inv_spatial_metric = make_with_value>(lapse, 0.0); // Mesh velocity set to std::null for now @@ -90,7 +90,7 @@ void test_flat_space_time_step() { inverse_jacobian_logical_to_inertial.get(0, 0) = 1.0; inverse_jacobian_logical_to_inertial.get(1, 1) = 1.0; inverse_jacobian_logical_to_inertial.get(2, 2) = 1.0; - Scalar det_jacobian_logical_to_inertial(dv_size, 1.0); + const Scalar det_jacobian_logical_to_inertial(dv_size, 1.0); // Coordinates tnsr::I mesh_coordinates = @@ -238,7 +238,7 @@ void test_flat_space_time_step() { 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 2f2aec53b353..d7dca4d480ee 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp @@ -88,11 +88,11 @@ struct component { Frame::Inertial>, Tags::deriv, tmpl::size_t, Frame::Inertial>, - Tags::deriv, + Tags::deriv, tmpl::size_t, Frame::Inertial>, gr::Tags::SpatialMetric, gr::Tags::InverseSpatialMetric, - gr::Tags::DetSpatialMetric, + gr::Tags::SqrtDetSpatialMetric, evolution::dg::subcell::Tags::Coordinates, domain::Tags::MeshVelocity, evolution::dg::subcell::fd::Tags::InverseJacobianLogicalToInertial, @@ -185,15 +185,15 @@ void test_advance_packets() { 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 = make_with_value>(lapse, 0.0); tnsr::iJ d_shift = make_with_value>(lapse, 0.0); - tnsr::iJJ d_inv_spatial_metric = - make_with_value>(lapse, 0.0); + tnsr::ijj d_spatial_metric = + make_with_value>(lapse, 0.0); // Fluid variables Scalar rest_mass_density(zero_dv); @@ -313,10 +313,10 @@ void test_advance_packets() { shift, d_lapse, d_shift, - d_inv_spatial_metric, + d_spatial_metric, spatial_metric, inv_spatial_metric, - determinant_spatial_metric, + sqrt_determinant_spatial_metric, mesh_coordinates, mesh_velocity, inverse_jacobian_logical_to_inertial, From 9bbe256109dc1b125666f39f166fef2361990c1b Mon Sep 17 00:00:00 2001 From: Francois Foucart Date: Tue, 29 Oct 2024 16:57:47 -0400 Subject: [PATCH 3/3] Make MC executable work with input files (optically thick sphere) --- .../MonteCarlo/CMakeLists.txt | 1 + .../MonteCarlo/EvolveMonteCarlo.cpp | 8 +- .../MonteCarlo/EvolveMonteCarlo.hpp | 42 +++- .../Actions/InitializeMonteCarlo.hpp | 50 ++-- .../Particles/MonteCarlo/CMakeLists.txt | 4 + .../Particles/MonteCarlo/CellCrossingTime.cpp | 1 + .../Particles/MonteCarlo/CellCrossingTime.hpp | 30 +++ .../MonteCarlo/GhostZoneCommunication.hpp | 41 +++- .../ImplicitMonteCarloCorrections.tpp | 2 +- .../MonteCarlo/MonteCarloOptions.cpp | 27 +++ .../MonteCarlo/MonteCarloOptions.hpp | 48 ++++ .../MonteCarlo/NeutrinoMomentsFromMC.cpp | 32 +++ .../MonteCarlo/NeutrinoMomentsFromMC.hpp | 71 ++++++ src/Evolution/Particles/MonteCarlo/Tags.hpp | 19 ++ .../RadiationTransport/CMakeLists.txt | 1 + .../MonteCarlo/CMakeLists.txt | 33 +++ .../RadiationTransport/MonteCarlo/Factory.hpp | 12 + .../MonteCarlo/HomogeneousSphere.cpp | 216 ++++++++++++++++++ .../MonteCarlo/HomogeneousSphere.hpp | 191 ++++++++++++++++ .../MonteCarlo/Solutions.hpp | 11 + .../EquationsOfState/EquationOfState.hpp | 2 +- .../MonteCarlo/Test_CommunicationTags.cpp | 57 ++++- 22 files changed, 864 insertions(+), 35 deletions(-) create mode 100644 src/Evolution/Particles/MonteCarlo/MonteCarloOptions.cpp create mode 100644 src/Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp create mode 100644 src/Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.cpp create mode 100644 src/Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.hpp create mode 100644 src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/CMakeLists.txt create mode 100644 src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Factory.hpp create mode 100644 src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.cpp create mode 100644 src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.hpp create mode 100644 src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Solutions.hpp diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt b/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt index 8680009bc4e2..b9e8aa2bb25d 100644 --- a/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt +++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt @@ -33,6 +33,7 @@ target_link_libraries( LinearOperators MathFunctions MonteCarlo + MonteCarloSolutions Observer Options Parallel diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp index db4a08f04d8b..bb8683a91b2e 100644 --- a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp +++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp @@ -20,6 +20,12 @@ void register_neutrino_tables() { Particles::MonteCarlo::NeutrinoInteractionTable<16, 3>>{}); } +void register_mc_options() { + register_classes_with_charm( + tmpl::list, + Particles::MonteCarlo::MonteCarloOptions<3>>{}); +} + extern "C" void CkRegisterMainModule() { Parallel::charmxx::register_main_module>(); Parallel::charmxx::register_init_node_and_proc( @@ -28,6 +34,6 @@ extern "C" void CkRegisterMainModule() { &domain::FunctionsOfTime::register_derived_with_charm, &EquationsOfState::register_derived_with_charm, ®ister_factory_classes_with_charm>, - ®ister_neutrino_tables}, + ®ister_neutrino_tables, ®ister_mc_options}, {}); } diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp index e1bdde9784fd..4a9ff52c7a1f 100644 --- a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp +++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp @@ -13,6 +13,8 @@ #include "Evolution/ComputeTags.hpp" #include "Evolution/DgSubcell/Actions/Initialize.hpp" #include "Evolution/DgSubcell/BackgroundGrVars.hpp" +#include "Evolution/DgSubcell/Tags/ObserverCoordinates.hpp" +#include "Evolution/DgSubcell/Tags/ObserverMesh.hpp" #include "Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp" #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp" #include "Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp" @@ -29,8 +31,11 @@ #include "Evolution/Initialization/SetVariables.hpp" #include "Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp" #include "Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp" +#include "Evolution/Particles/MonteCarlo/CellCrossingTime.hpp" #include "Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp" #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp" +#include "Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp" +#include "Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.hpp" #include "Evolution/Particles/MonteCarlo/System.hpp" #include "Evolution/Particles/MonteCarlo/Tags.hpp" #include "Evolution/Systems/GrMhd/ValenciaDivClean/AllSolutions.hpp" @@ -57,6 +62,7 @@ #include "ParallelAlgorithms/Actions/MutateApply.hpp" #include "ParallelAlgorithms/Actions/TerminatePhase.hpp" #include "ParallelAlgorithms/Events/Factory.hpp" +#include "ParallelAlgorithms/Events/ObserveFields.hpp" #include "ParallelAlgorithms/EventsAndDenseTriggers/DenseTrigger.hpp" #include "ParallelAlgorithms/EventsAndDenseTriggers/DenseTriggers/Factory.hpp" #include "ParallelAlgorithms/EventsAndTriggers/Completion.hpp" @@ -68,7 +74,7 @@ #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/RadiationTransport/MonteCarlo/Factory.hpp" #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp" #include "PointwiseFunctions/Hydro/LowerSpatialFourVelocity.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -118,7 +124,7 @@ struct EvolutionMetavars { static constexpr bool use_dg_subcell = true; using initial_data_list = - grmhd::ValenciaDivClean::InitialData::initial_data_list; + RadiationTransport::MonteCarlo::Solutions::all_solutions; using equation_of_state_tag = hydro::Tags::GrmhdEquationOfState; struct SubcellOptions { @@ -129,12 +135,20 @@ struct EvolutionMetavars { }; using observe_fields = - tmpl::list, - domain::Tags::Coordinates>; + tmpl::list, + evolution::dg::subcell::Tags::ObserverCoordinatesCompute< + volume_dim, Frame::Grid>, + evolution::dg::subcell::Tags::ObserverCoordinatesCompute< + volume_dim, Frame::Inertial>, + Particles::MonteCarlo::Tags::InertialFrameEnergyDensity>; using non_tensor_compute_tags = - tmpl::list<::Events::Tags::ObserverMeshCompute, - ::Events::Tags::ObserverDetInvJacobianCompute< - Frame::ElementLogical, Frame::Inertial>>; + tmpl::list, + evolution::dg::subcell::Tags:: + ObserverJacobianAndDetInvJacobianCompute< + volume_dim, Frame::ElementLogical, Frame::Inertial>, + evolution::dg::subcell::Tags::ObserverInverseJacobianCompute< + volume_dim, Frame::ElementLogical, Frame::Inertial>>; using analytic_variables_tags = typename system::variables_tag::tags_list; using analytic_compute = evolution::Tags::AnalyticSolutionsCompute< @@ -145,7 +159,13 @@ struct EvolutionMetavars { using factory_classes = tmpl::map< tmpl::pair, tmpl::pair, domain_creators>, - tmpl::pair>>, + tmpl::pair, + dg::Events::ObserveFields>>>, tmpl::pair, tmpl::pair< grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField, @@ -187,6 +207,8 @@ struct EvolutionMetavars { NeutrinoSpecies>, Initialization::Actions::AddComputeTags>>, evolution::Actions::InitializeRunEventsAndDenseTriggers, @@ -197,6 +219,9 @@ struct EvolutionMetavars { tmpl::list< Parallel::PhaseActions, + Parallel::PhaseActions>, Parallel::PhaseActions< Parallel::Phase::Evolve, tmpl::list< @@ -237,6 +262,7 @@ struct EvolutionMetavars { dg_element_array>; using const_global_cache_tags = tmpl::list< + Particles::MonteCarlo::Tags::MonteCarloOptions, equation_of_state_tag, evolution::initial_data::Tags::InitialData, Particles::MonteCarlo::Tags::InteractionRatesTable>; diff --git a/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp b/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp index d36e4b2dda05..c5572953edd1 100644 --- a/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp +++ b/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp @@ -13,6 +13,7 @@ #include "Evolution/DgSubcell/Tags/Mesh.hpp" #include "Evolution/Initialization/InitialData.hpp" #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp" +#include "Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp" #include "Evolution/Particles/MonteCarlo/MortarData.hpp" #include "Evolution/Particles/MonteCarlo/Packet.hpp" #include "Evolution/Particles/MonteCarlo/Tags.hpp" @@ -48,6 +49,9 @@ namespace Initialization::Actions { /// - evolution::dg::subcell::Tags::Mesh /// - evolution::dg::subcell::Tags::Coordinates /// - evolution::initial_data::Tags::InitialData +/// - Particles::MonteCarlo::Tags::MonteCarloOptions +/// - domain::Tags::Element /// /// DataBox changes: /// - Adds: @@ -55,7 +59,6 @@ namespace Initialization::Actions { /// * Particles::MonteCarlo::Tags::RandomNumberGenerator /// * Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission< /// NeutrinoSpecies> -/// * Particles::MonteCarlo::Tags::CellLightCrossingTime /// * Background hydro variables /// * Particles::MonteCarlo::Tags::MortarDataTag /// * Particles::MonteCarlo::Tags::McGhostZoneDataTag @@ -73,7 +76,6 @@ struct InitializeMCTags { Particles::MonteCarlo::Tags::RandomNumberGenerator, Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission< NeutrinoSpecies>, - Particles::MonteCarlo::Tags::CellLightCrossingTime, hydro_variables_tag, Particles::MonteCarlo::Tags::MortarDataTag, Particles::MonteCarlo::Tags::McGhostZoneDataTag>; @@ -89,9 +91,9 @@ struct InitializeMCTags { 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(); + const Mesh& mesh = + db::get>(box); + const size_t num_grid_points = mesh.number_of_grid_points(); using derived_classes = tmpl::at; @@ -116,39 +118,55 @@ struct InitializeMCTags { make_not_null(&box), std::move(hydro_variables)); }); + // Read global options for Monte-Carlo evolution + const auto mc_options = db::get< + Particles::MonteCarlo::Tags::MonteCarloOptions>(box); + const auto& initial_packet_energy = mc_options.get_initial_packet_energy(); + typename Particles::MonteCarlo::Tags::PacketsOnElement::type all_packets; Initialization::mutate_assign< tmpl::list>( make_not_null(&box), std::move(all_packets)); - const unsigned long seed = - std::random_device{}(); // static_cast(time(NULL)); + const unsigned long seed = std::random_device{}(); 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 + // Initial energy of packets, read from MC options typename Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission< NeutrinoSpecies>::type packet_energy_at_emission = make_with_value>( - DataVector{num_grid_points}, 1.e-12); + DataVector{num_grid_points}, 0.0); + for (size_t s = 0; s < NeutrinoSpecies; s++) { + packet_energy_at_emission[s] = initial_packet_energy[s]; + } 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? + // Initialize mortar data. Currently assumes a single neighbor on each + // face (i.e. no h-refinement) using MortarData = typename Particles::MonteCarlo::Tags::MortarDataTag::type; MortarData mortar_data; + const Element& element = db::get<::domain::Tags::Element>(box); + for (const auto& [direction, neighbors] : element.neighbors()) { + const size_t sliced_mesh_size = + mesh.slice_away(direction.dimension()).number_of_grid_points(); + DataVector zero_dv_ghost_zones(sliced_mesh_size, 0.0); + for (const auto& neighbor : neighbors) { + const DirectionalId mortar_id{direction, neighbor}; + mortar_data.rest_mass_density.emplace(mortar_id, zero_dv_ghost_zones); + mortar_data.electron_fraction.emplace(mortar_id, zero_dv_ghost_zones); + mortar_data.temperature.emplace(mortar_id, zero_dv_ghost_zones); + mortar_data.cell_light_crossing_time.emplace(mortar_id, + zero_dv_ghost_zones); + } + } Initialization::mutate_assign< tmpl::list>>( make_not_null(&box), std::move(mortar_data)); diff --git a/src/Evolution/Particles/MonteCarlo/CMakeLists.txt b/src/Evolution/Particles/MonteCarlo/CMakeLists.txt index daca3ec8939d..6f6b60acd857 100644 --- a/src/Evolution/Particles/MonteCarlo/CMakeLists.txt +++ b/src/Evolution/Particles/MonteCarlo/CMakeLists.txt @@ -13,7 +13,9 @@ spectre_target_sources( CouplingTermsForPropagation.cpp EvolvePackets.cpp InverseJacobianInertialToFluidCompute.cpp + MonteCarloOptions.cpp NeutrinoInteractionTable.cpp + NeutrinoMomentsFromMC.cpp Packet.cpp Scattering.cpp TemplatedLocalFunctions.cpp @@ -34,8 +36,10 @@ spectre_target_headers( GhostZoneCommunicationTags.hpp ImplicitMonteCarloCorrections.tpp InverseJacobianInertialToFluidCompute.hpp + MonteCarloOptions.hpp MortarData.hpp NeutrinoInteractionTable.hpp + NeutrinoMomentsFromMC.hpp Packet.hpp Scattering.hpp System.hpp diff --git a/src/Evolution/Particles/MonteCarlo/CellCrossingTime.cpp b/src/Evolution/Particles/MonteCarlo/CellCrossingTime.cpp index 44929be86844..0a2548ae5bf4 100644 --- a/src/Evolution/Particles/MonteCarlo/CellCrossingTime.cpp +++ b/src/Evolution/Particles/MonteCarlo/CellCrossingTime.cpp @@ -25,6 +25,7 @@ void cell_light_crossing_time( inertial_coordinates.get(1)[step[1]] - inertial_coordinates.get(1)[0], inertial_coordinates.get(2)[step[2]] - inertial_coordinates.get(2)[0]}; + get(*cell_light_crossing_time) = DataVector(n_pts,0.0); // Estimate light-crossing time in the cell. for (size_t i = 0; i < n_pts; i++) { double& min_crossing_time = get(*cell_light_crossing_time)[i]; diff --git a/src/Evolution/Particles/MonteCarlo/CellCrossingTime.hpp b/src/Evolution/Particles/MonteCarlo/CellCrossingTime.hpp index 1a0b01291b32..05d9f187020b 100644 --- a/src/Evolution/Particles/MonteCarlo/CellCrossingTime.hpp +++ b/src/Evolution/Particles/MonteCarlo/CellCrossingTime.hpp @@ -4,6 +4,12 @@ #pragma once #include "DataStructures/Tensor/TypeAliases.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/DgSubcell/Tags/Coordinates.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/Particles/MonteCarlo/Tags.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "Utilities/Gsl.hpp" /// \cond namespace gsl { @@ -27,4 +33,28 @@ void cell_light_crossing_time( const tnsr::I& shift, const tnsr::II& inv_spatial_metric); +struct CellLightCrossingTimeCompute : Tags::CellLightCrossingTime, + db::ComputeTag { + using base = Tags::CellLightCrossingTime; + using return_type = typename base::type; + using argument_tags = tmpl::list< + evolution::dg::subcell::Tags::Mesh<3>, + evolution::dg::subcell::Tags::Coordinates<3, Frame::Inertial>, + gr::Tags::Lapse, + gr::Tags::Shift, + gr::Tags::InverseSpatialMetric>; + + static void function( + gsl::not_null cell_light_crossing_time_, + const Mesh<3>& mesh, + const tnsr::I& inertial_coordinates, + const Scalar& lapse, + const tnsr::I& shift, + const tnsr::II& inv_spatial_metric){ + cell_light_crossing_time(cell_light_crossing_time_, mesh, + inertial_coordinates, lapse, shift, inv_spatial_metric); + } +}; + + } // namespace Particles::MonteCarlo diff --git a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp index 0ced96571ea5..7f9128e768d8 100644 --- a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp +++ b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp @@ -23,6 +23,7 @@ #include "Evolution/DgSubcell/GhostData.hpp" #include "Evolution/DgSubcell/SliceData.hpp" #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp" +#include "Evolution/DgSubcell/Tags/Coordinates.hpp" #include "Evolution/DgSubcell/Tags/Interpolators.hpp" #include "Evolution/DgSubcell/Tags/Mesh.hpp" #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp" @@ -364,16 +365,31 @@ struct ReceiveDataForMcCommunication { }, make_not_null(&box)); } else { + const Mesh& subcell_mesh = + db::get>(box); + auto& mesh_coordinates = + db::get>(box); // TO DO: Deal with data coupling neutrinos back to fluid evolution db::mutate( - [&element, &received_data]( + [&element, &received_data, &subcell_mesh, &mesh_coordinates]( const gsl::not_null*> packet_list) { + const Index& extents = subcell_mesh.extents(); + std::array bottom_coord_mesh; + std::array step; + std::array dx_mesh; + for (size_t d = 0; d < Dim; d++) { + bottom_coord_mesh[d] = mesh_coordinates.get(d)[0]; + step[d] = (d == 0) ? 1 : step[d - 1] * extents[d - 1]; + dx_mesh[d] = + mesh_coordinates.get(d)[step[d]] - bottom_coord_mesh[d]; + } for (const auto& [direction, neighbors_in_direction] : element.neighbors()) { for (const auto& neighbor : neighbors_in_direction) { DirectionalId directional_element_id{direction, neighbor}; - const std::optional>& + std::optional>& received_data_packets = received_data[directional_element_id] .packets_entering_this_element; @@ -385,6 +401,26 @@ struct ReceiveDataForMcCommunication { } else { const size_t n_packets = received_data_packets.value().size(); for (size_t p = 0; p < n_packets; p++) { + // Find closest grid point to received packet at current + // time, using extents for live points only. + Packet& packet = received_data_packets.value()[p]; + std::array closest_point_index_3d; + for (size_t d = 0; d < Dim; 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]; + size_t shift = extents[0]; + for (size_t d = 1; d < Dim; d++) { + packet.index_of_closest_grid_point += + shift * closest_point_index_3d[d]; + shift *= extents[d]; + } + // Now add packets to list in current element packet_list->push_back(received_data_packets.value()[p]); } } @@ -393,7 +429,6 @@ struct ReceiveDataForMcCommunication { }, make_not_null(&box)); } - return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } }; diff --git a/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp b/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp index 699e76bfcf47..d76c7eba77e6 100644 --- a/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp +++ b/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp @@ -48,7 +48,7 @@ void TemplatedLocalFunctions:: interaction_table.get_neutrino_energies(); // Apply implicit MC corrections as needed - const size_t dv_size = rest_mass_density.size(); + const size_t dv_size = get(rest_mass_density).size(); // Calculate beta parameter (relative change of MC vs fluid variables) // For photon transport, diff --git a/src/Evolution/Particles/MonteCarlo/MonteCarloOptions.cpp b/src/Evolution/Particles/MonteCarlo/MonteCarloOptions.cpp new file mode 100644 index 000000000000..b0fd005c8198 --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/MonteCarloOptions.cpp @@ -0,0 +1,27 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp" + +#include +#include +#include +#include +#include +#include + +namespace Particles::MonteCarlo { + +template +PUP::able::PUP_ID MonteCarloOptions::my_PUP_ID = 0; // NOLINT + +template +void MonteCarloOptions::pup(PUP::er& p) { + PUP::able::pup(p); + p | initial_packet_energy_; +} + +} // namespace Particles::MonteCarlo + +template class Particles::MonteCarlo::MonteCarloOptions<2>; +template class Particles::MonteCarlo::MonteCarloOptions<3>; diff --git a/src/Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp b/src/Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp new file mode 100644 index 000000000000..c8b6ed8c994a --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp @@ -0,0 +1,48 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include + +#include "Options/String.hpp" +#include "Utilities/Serialization/CharmPupable.hpp" +#include "Utilities/TMPL.hpp" + +namespace Particles::MonteCarlo { + +template +class MonteCarloOptions : public PUP::able { + public: + explicit MonteCarloOptions( + const std::array initial_packet_energy) + : initial_packet_energy_(initial_packet_energy) {} + + static constexpr Options::String help = { + "Global options for Monte-Carlo evolution.\n" + "InitialPacketEnergy: [double, double, double] \n"}; + + struct InitialPacketEnergy { + using type = std::array; + static constexpr Options::String help{ + "Initial energy used to create packets"}; + }; + + using options = tmpl::list; + + explicit MonteCarloOptions(CkMigrateMessage* msg) : PUP::able(msg) {} + + using PUP::able::register_constructor; + void pup(PUP::er& p) override; + WRAPPED_PUPable_decl_template(MonteCarloOptions); + + const std::array& get_initial_packet_energy() const { + return initial_packet_energy_; + } + + private: + std::array initial_packet_energy_; +}; + +} // namespace Particles::MonteCarlo diff --git a/src/Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.cpp b/src/Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.cpp new file mode 100644 index 000000000000..08cbbf8d6852 --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.cpp @@ -0,0 +1,32 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.hpp" + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Evolution/Particles/MonteCarlo/CellVolume.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" + +namespace Particles::MonteCarlo { + +void inertial_frame_energy_density( + gsl::not_null*> fluid_frame_energy_density, + const std::vector& packets, const Scalar& lapse, + const Scalar& sqrt_determinant_spatial_metric, + const Mesh<3>& mesh, + const Scalar& det_jacobian_logical_to_inertial) { + Scalar cell_inertial_three_volume(lapse); + cell_inertial_coordinate_three_volume_finite_difference( + &cell_inertial_three_volume, mesh, det_jacobian_logical_to_inertial); + get(*fluid_frame_energy_density) = 0.0 * get(lapse); + for (const auto& packet : packets) { + const size_t& idx = packet.index_of_closest_grid_point; + get(*fluid_frame_energy_density)[idx] += + get(lapse)[idx] / get(cell_inertial_three_volume)[idx] / + get(sqrt_determinant_spatial_metric)[idx] * packet.momentum_upper_t * + packet.number_of_neutrinos; + } +} + +} // namespace Particles::MonteCarlo diff --git a/src/Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.hpp b/src/Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.hpp new file mode 100644 index 000000000000..fb8259e5afc2 --- /dev/null +++ b/src/Evolution/Particles/MonteCarlo/NeutrinoMomentsFromMC.hpp @@ -0,0 +1,71 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/DgSubcell/Tags/Jacobians.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/Particles/MonteCarlo/Packet.hpp" +#include "Evolution/Particles/MonteCarlo/Tags.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "Utilities/Gsl.hpp" + +/// \cond +namespace gsl { +template +class not_null; +} // namespace gsl + +class DataVector; + +template +class Mesh; +/// \endcond + +namespace Particles::MonteCarlo { + +void inertial_frame_energy_density( + gsl::not_null*> fluid_frame_energy_density, + const std::vector& packets, const Scalar& lapse, + const Scalar& sqrt_determinant_spatial_metric, + const Mesh<3>& mesh, + const Scalar& det_jacobian_logical_to_inertial); + +namespace Tags { +/// Simple tag containing the inertial frame energy +/// density on the grid for MC packets +struct InertialFrameEnergyDensity : db::SimpleTag { + using type = Scalar; +}; +} // namespace Tags + +struct InertialFrameEnergyDensityCompute : Tags::InertialFrameEnergyDensity, + db::ComputeTag { + using base = Tags::InertialFrameEnergyDensity; + using return_type = typename base::type; + using argument_tags = tmpl::list< + Particles::MonteCarlo::Tags::PacketsOnElement, + gr::Tags::Lapse, gr::Tags::SqrtDetSpatialMetric, + evolution::dg::subcell::Tags::Mesh<3>, + evolution::dg::subcell::fd::Tags::DetInverseJacobianLogicalToInertial>; + + static void function( + gsl::not_null inertial_frame_density, + const std::vector& packets, const Scalar& lapse, + const Scalar& sqrt_determinant_spatial_metric, + const Mesh<3>& mesh, + const Scalar& det_inv_jacobian_logical_to_inertial) { + Scalar det_jacobian_logical_to_inertial(lapse); + get(det_jacobian_logical_to_inertial) = + 1.0 / get(det_inv_jacobian_logical_to_inertial); + inertial_frame_energy_density(inertial_frame_density, packets, lapse, + sqrt_determinant_spatial_metric, mesh, + det_jacobian_logical_to_inertial); + } +}; + +} // namespace Particles::MonteCarlo diff --git a/src/Evolution/Particles/MonteCarlo/Tags.hpp b/src/Evolution/Particles/MonteCarlo/Tags.hpp index e0274fc1ae89..6954604cbc17 100644 --- a/src/Evolution/Particles/MonteCarlo/Tags.hpp +++ b/src/Evolution/Particles/MonteCarlo/Tags.hpp @@ -8,6 +8,7 @@ #include "DataStructures/DataBox/Tag.hpp" #include "DataStructures/Tensor/TypeAliases.hpp" +#include "Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp" #include "Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp" #include "Evolution/Particles/MonteCarlo/Packet.hpp" @@ -65,4 +66,22 @@ struct InteractionRatesTable : db::SimpleTag { } }; +template +struct MonteCarloOptions : db::SimpleTag { + using type = std::unique_ptr< + Particles::MonteCarlo::MonteCarloOptions>; + static constexpr bool pass_metavariables = false; + using option_tags = typename Particles::MonteCarlo::MonteCarloOptions< + NeutrinoSpecies>::options; + static type create_from_options( + const std::array initial_packet_energy) { + std::unique_ptr> + mc_options_ptr = std::make_unique< + Particles::MonteCarlo::MonteCarloOptions>( + initial_packet_energy); + return mc_options_ptr; + ; + } +}; + } // namespace Particles::MonteCarlo::Tags diff --git a/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/CMakeLists.txt b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/CMakeLists.txt index 9a5cb7938cb7..060886a54225 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/CMakeLists.txt +++ b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/CMakeLists.txt @@ -2,3 +2,4 @@ # See LICENSE.txt for details. add_subdirectory(M1Grey) +add_subdirectory(MonteCarlo) diff --git a/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/CMakeLists.txt b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/CMakeLists.txt new file mode 100644 index 000000000000..4bf586cde573 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/CMakeLists.txt @@ -0,0 +1,33 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +set(LIBRARY MonteCarloSolutions) + +add_spectre_library(${LIBRARY}) + +spectre_target_sources( + ${LIBRARY} + PRIVATE + HomogeneousSphere.cpp + ) + +spectre_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src + HEADERS + HomogeneousSphere.hpp + Factory.hpp + Solutions.hpp + ) + +target_link_libraries( + ${LIBRARY} + PUBLIC + Boost::boost + DataStructures + ErrorHandling + GeneralRelativitySolutions + Hydro + MonteCarlo + Options + ) diff --git a/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Factory.hpp b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Factory.hpp new file mode 100644 index 000000000000..bf6a5214e323 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Factory.hpp @@ -0,0 +1,12 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.hpp" +#include "Utilities/TMPL.hpp" + +namespace RadiationTransport::MonteCarlo::Solutions { +/// \brief List of all analytic solutions +using all_solutions = tmpl::list; +} // namespace RadiationTransport::MonteCarlo::Solutions diff --git a/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.cpp b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.cpp new file mode 100644 index 000000000000..8032f36e97d9 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.cpp @@ -0,0 +1,216 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.hpp" + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/MakeWithValue.hpp" + +namespace RadiationTransport::MonteCarlo::Solutions { + +namespace { +double get_mask(const tnsr::I& x, const double radius_bound) { + const double radius = sqrt((x.get(0) * x.get(0)) + (x.get(1) * x.get(1)) + + (x.get(2) * x.get(2))); + return radius > radius_bound ? 1.0 : 0.0; +} + +DataVector get_mask(const tnsr::I& x, + const double radius_bound) { + const DataVector radius = sqrt((x.get(0) * x.get(0)) + (x.get(1) * x.get(1)) + + (x.get(2) * x.get(2))); + DataVector mask = make_with_value(radius, 0.0); + for (size_t i = 0; i < mask.size(); i++) { + mask[i] = radius[i] > radius_bound ? 1.0 : 0.0; + } + return mask; +} +} // namespace + +HomogeneousSphere::HomogeneousSphere( + const double& radius, const std::array& densities, + const std::array& temperatures, + const std::array& electron_fractions, + std::unique_ptr> + local_eos) + : radius_(radius), + densities_(densities), + temperatures_(temperatures), + electron_fractions_(electron_fractions), + equation_of_state_(std::move(local_eos)) {} + +std::unique_ptr +HomogeneousSphere::get_clone() const { + return std::make_unique(*this); +} + +HomogeneousSphere::HomogeneousSphere(CkMigrateMessage* /*unused*/) {} + +HomogeneousSphere::HomogeneousSphere(const HomogeneousSphere& rhs) + : radius_(rhs.radius_), + densities_(rhs.densities_), + temperatures_(rhs.temperatures_), + electron_fractions_(rhs.electron_fractions_), + equation_of_state_(rhs.equation_of_state_->get_clone()) {} + +HomogeneousSphere& HomogeneousSphere::operator=(const HomogeneousSphere& rhs) { + radius_ = rhs.radius_; + densities_ = rhs.densities_; + temperatures_ = rhs.temperatures_; + electron_fractions_ = rhs.electron_fractions_; + equation_of_state_ = rhs.equation_of_state_->get_clone(); + return *this; +} + +void HomogeneousSphere::pup(PUP::er& p) { + p | radius_; + p | densities_; + p | temperatures_; + p | electron_fractions_; + p | equation_of_state_; +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const { + const DataType mask = get_mask(x, radius_); + return {Scalar{DataType{electron_fractions_[1] * mask + + electron_fractions_[0] * (1.0 - mask)}}}; +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const { + const DataType mask = get_mask(x, radius_); + return {Scalar{ + DataType{densities_[1] * mask + densities_[0] * (1.0 - mask)}}}; +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const { + const DataType mask = get_mask(x, radius_); + return {Scalar{ + DataType{temperatures_[1] * mask + temperatures_[0] * (1.0 - mask)}}}; +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const { + return {make_with_value>(x, 1.0)}; +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double t, + tmpl::list> /*meta*/) const { + auto primitives = this->variables( + x, t, + tmpl::list, + hydro::Tags::Temperature, + hydro::Tags::ElectronFraction>{}); + return equation_of_state_ + ->specific_internal_energy_from_density_and_temperature( + tuples::get>(primitives), + tuples::get>(primitives), + tuples::get>(primitives)); +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double t, + tmpl::list> /*meta*/) const { + auto primitives = this->variables( + x, t, + tmpl::list, + hydro::Tags::Temperature, + hydro::Tags::ElectronFraction>{}); + return equation_of_state_->pressure_from_density_and_temperature( + tuples::get>(primitives), + tuples::get>(primitives), + tuples::get>(primitives)); +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const { + return {make_with_value>(x, 0.0)}; +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const { + return {make_with_value>(x, 0.0)}; +} + +template +tuples::TaggedTuple> +HomogeneousSphere::variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const { + return {make_with_value>(x, 0.0)}; +} + +PUP::able::PUP_ID HomogeneousSphere::my_PUP_ID = 0; // NOLINT + +bool operator==(const HomogeneousSphere& lhs, const HomogeneousSphere& rhs) { + return (lhs.radius_ == rhs.radius_ && lhs.densities_ == rhs.densities_ && + lhs.temperatures_ == rhs.temperatures_ && + lhs.electron_fractions_ == rhs.electron_fractions_ && + *lhs.equation_of_state_ == *rhs.equation_of_state_); +} + +bool operator!=(const HomogeneousSphere& lhs, const HomogeneousSphere& rhs) { + return not(lhs == rhs); +} + +#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) +#define TAG(data) BOOST_PP_TUPLE_ELEM(1, data) + +#define INSTANTIATE_SCALARS(_, data) \ + template tuples::TaggedTuple > \ + HomogeneousSphere::variables( \ + const tnsr::I& x, double t, \ + tmpl::list > /*meta*/) const; + +GENERATE_INSTANTIATIONS(INSTANTIATE_SCALARS, (double, DataVector), + (hydro::Tags::DivergenceCleaningField, + hydro::Tags::RestMassDensity, + hydro::Tags::ElectronFraction, + hydro::Tags::Temperature, hydro::Tags::LorentzFactor, + hydro::Tags::SpecificInternalEnergy, + hydro::Tags::Pressure)) + +#define INSTANTIATE_VECTORS(_, data) \ + template tuples::TaggedTuple > \ + HomogeneousSphere::variables( \ + const tnsr::I& x, double t, \ + tmpl::list > /*meta*/) const; + +GENERATE_INSTANTIATIONS(INSTANTIATE_VECTORS, (double, DataVector), + (hydro::Tags::MagneticField, + hydro::Tags::SpatialVelocity)) + +#undef DTYPE +#undef TAG +#undef INSTANTIATE_SCALARS +#undef INSTANTIATE_VECTORS +} // namespace RadiationTransport::MonteCarlo::Solutions diff --git a/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.hpp b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.hpp new file mode 100644 index 000000000000..0420076b0de9 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/HomogeneousSphere.hpp @@ -0,0 +1,191 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include + +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Evolution/Particles/MonteCarlo/Tags.hpp" +#include "Options/String.hpp" +#include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "PointwiseFunctions/Hydro/TagsDeclarations.hpp" +#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" + +/// \cond +class DataVector; +/// \endcond + +namespace RadiationTransport::MonteCarlo::Solutions { + +/*! + * \brief Homogeneous sphere as fluid background to MC run + * + * Provides background fluid variables for a + * fluid with constant density, temperature, Ye + * in Minkowski spacetime. + * + */ +class HomogeneousSphere : public evolution::initial_data::InitialData, + public MarkAsAnalyticSolution { + public: + constexpr static bool IsRelativistic = true; + using equation_of_state_type = + EquationsOfState::EquationOfState; + + const EquationsOfState::EquationOfState& + equation_of_state() const { + return *equation_of_state_; + } + + /// The radius of the sphere + struct Radius { + using type = double; + static constexpr Options::String help = {"The radius of the sphere."}; + }; + /// The density inside and outside the sphere + struct Densities { + using type = std::array; + static constexpr Options::String help = {"Density inside and outside."}; + }; + /// The temperature inside and outside the sphere + struct Temperatures { + using type = std::array; + static constexpr Options::String help = {"Temperature inside and outside."}; + }; + /// The electron fraction inside and outside the sphere + struct ElectronFractions { + using type = std::array; + static constexpr Options::String help = {"Ye inside and outside."}; + }; + + using options = tmpl::list< + Radius, Densities, Temperatures, ElectronFractions, + hydro::OptionTags::InitialDataEquationOfState>; + static constexpr Options::String help = { + "Background for uniform sphere with constant rho, T, Ye"}; + + HomogeneousSphere() = default; + HomogeneousSphere(const HomogeneousSphere& /*rhs*/); + HomogeneousSphere& operator=(const HomogeneousSphere& /*rhs*/); + HomogeneousSphere(HomogeneousSphere&& /*rhs*/) = default; + HomogeneousSphere& operator=(HomogeneousSphere&& /*rhs*/) = default; + ~HomogeneousSphere() override = default; + + HomogeneousSphere( + const double& radius, const std::array& densities, + const std::array& temperatures, + const std::array& electron_fractions, + std::unique_ptr> + local_eos); + + auto get_clone() const + -> std::unique_ptr override; + + /// \cond + explicit HomogeneousSphere(CkMigrateMessage* msg); + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(HomogeneousSphere); + /// \endcond + + /// @{ + /// Retrieve fluid variables at `(x, t)` + template + auto variables(const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) + const -> tuples::TaggedTuple>; + + template + auto variables(const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) + const -> tuples::TaggedTuple>; + + template + auto variables(const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const + -> tuples::TaggedTuple>; + + template + auto variables(const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) + const -> tuples::TaggedTuple>; + + template + auto variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const + -> tuples::TaggedTuple>; + + template + auto variables(const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const + -> tuples::TaggedTuple>; + + template + auto variables(const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) + const -> tuples::TaggedTuple>; + + template + auto variables(const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) + const -> tuples::TaggedTuple>; + + template + auto variables( + const tnsr::I& x, double /*t*/, + tmpl::list> /*meta*/) const + -> tuples::TaggedTuple>; + /// @} + + /// Retrieve a collection of hydro variables at `(x, t)` + template + tuples::TaggedTuple variables(const tnsr::I& x, + double t, + tmpl::list /*meta*/) const { + static_assert(sizeof...(Tags) > 1, + "The generic template will recurse infinitely if only one " + "tag is being retrieved."); + return {get(variables(x, t, tmpl::list{}))...}; + } + + /// Retrieve the metric variables + template + tuples::TaggedTuple variables(const tnsr::I& x, double t, + tmpl::list /*meta*/) const { + return background_spacetime_.variables(x, t, tmpl::list{}); + } + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& /*p*/) override; + + private: + friend bool operator==(const HomogeneousSphere& lhs, + const HomogeneousSphere& rhs); + + double radius_ = std::numeric_limits::signaling_NaN(); + std::array densities_{ + {std::numeric_limits::signaling_NaN(), + std::numeric_limits::signaling_NaN()}}; + std::array temperatures_{ + {std::numeric_limits::signaling_NaN(), + std::numeric_limits::signaling_NaN()}}; + std::array electron_fractions_{ + {std::numeric_limits::signaling_NaN(), + std::numeric_limits::signaling_NaN()}}; + + std::unique_ptr equation_of_state_; + gr::Solutions::Minkowski<3> background_spacetime_{}; +}; + +bool operator!=(const HomogeneousSphere& lhs, const HomogeneousSphere& rhs); +} // namespace RadiationTransport::MonteCarlo::Solutions diff --git a/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Solutions.hpp b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Solutions.hpp new file mode 100644 index 000000000000..607e97eadeab --- /dev/null +++ b/src/PointwiseFunctions/AnalyticSolutions/RadiationTransport/MonteCarlo/Solutions.hpp @@ -0,0 +1,11 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +/*! + * \ingroup AnalyticSolutionsGroup + * \brief Holds classes providing background fluid/metric for MC + */ +namespace RadiationTransport::MonteCarlo::Solutions { +} // namespace RadiationTransport::MonteCarlo::Solutions diff --git a/src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp b/src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp index e9a0a886aaeb..7227a16f108c 100644 --- a/src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp +++ b/src/PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp @@ -568,7 +568,7 @@ class EquationOfState : public PUP::able { const EquationOfState& rhs) const = 0; virtual std::unique_ptr> - promote_to_3d_eos() { + promote_to_3d_eos() const { return this->get_clone(); } diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp index 46aaa87fdf8a..79d44b0c5bba 100644 --- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp +++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp @@ -66,6 +66,7 @@ struct component { using initial_tags = tmpl::list< ::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>, domain::Tags::Mesh, evolution::dg::subcell::Tags::Mesh, + evolution::dg::subcell::Tags::Coordinates, evolution::dg::subcell::Tags::ActiveGrid, domain::Tags::Element, Particles::MonteCarlo::Tags::McGhostZoneDataTag, Particles::MonteCarlo::Tags::PacketsOnElement, @@ -121,6 +122,7 @@ void test_send_receive_actions() { const Mesh subcell_mesh = evolution::dg::subcell::fd::mesh(dg_mesh); const evolution::dg::subcell::ActiveGrid active_grid = evolution::dg::subcell::ActiveGrid::Subcell; + const size_t mesh_size = subcell_mesh.extents()[0]; // ^ eta // +-+-+> xi @@ -201,6 +203,39 @@ void test_send_receive_actions() { get<0>(logical_coordinates(subcell_mesh)) * 2.0); std::vector packets_on_element{}; + // Logical coordinates + tnsr::I mesh_coordinates = + make_with_value>( + cell_light_crossing_time, 0.0); + for (size_t ix = 0; ix < mesh_size; ix++) { + const double x_coord = -1.0 + (0.5 + static_cast(ix)) / + static_cast(mesh_size) * 2.0; + if constexpr (Dim == 1) { + const size_t idx = ix; + mesh_coordinates.get(0)[idx] = x_coord; + } else { + for (size_t iy = 0; iy < mesh_size; iy++) { + const double y_coord = -1.0 + (0.5 + static_cast(iy)) / + static_cast(mesh_size) * 2.0; + if constexpr (Dim == 2) { + const size_t idx = ix + iy * mesh_size; + mesh_coordinates.get(0)[idx] = x_coord; + mesh_coordinates.get(1)[idx] = y_coord; + } else { + for (size_t iz = 0; iz < mesh_size; iz++) { + const double z_coord = -1.0 + (0.5 + static_cast(iz)) / + static_cast(mesh_size) * + 2.0; + const size_t idx = ix + iy * mesh_size + iz * mesh_size * mesh_size; + mesh_coordinates.get(0)[idx] = x_coord; + mesh_coordinates.get(1)[idx] = y_coord; + mesh_coordinates.get(2)[idx] = z_coord; + } + } + } + } + } + for (const auto& [direction, neighbor_ids] : neighbors) { (void)direction; for (const auto& neighbor_id : neighbor_ids) { @@ -211,7 +246,8 @@ void test_send_receive_actions() { &runner, ActionTesting::NodeId{0}, ActionTesting::LocalCoreId{0}, neighbor_id, {time_step_id, next_time_step_id, Mesh{}, Mesh{}, - active_grid, Element{}, NeighborDataMap{}, packets_on_element, + tnsr::I{}, active_grid, + Element{}, NeighborDataMap{}, packets_on_element, Variables{}, typename Particles::MonteCarlo::Tags::MortarDataTag::type{}, Scalar{}, Scalar{}, Scalar{}, @@ -252,11 +288,12 @@ void test_send_receive_actions() { const size_t species = 1; const double number_of_neutrinos = 2.0; - const size_t index_of_closest_grid_point = 0; const double t0 = 1.2; const double x0 = 0.3; const double y0 = 0.5; const double z0 = -0.7; + const size_t index_of_closest_grid_point = + (Dim == 1) ? 5 : ((Dim == 2) ? 5 + 6 * 9 : 5 + 6 * 9 + 1 * 81); const double p_upper_t0 = 1.1; const double p_x0 = 0.9; const double p_y0 = 0.7; @@ -290,9 +327,9 @@ void test_send_receive_actions() { Interps fd_to_neighbor_fd_interpolants{}; ActionTesting::emplace_array_component_and_initialize( &runner, ActionTesting::NodeId{0}, ActionTesting::LocalCoreId{0}, self_id, - {time_step_id, next_time_step_id, dg_mesh, subcell_mesh, active_grid, - element, neighbor_data, packets_on_element, evolved_vars, mortar_data, - rest_mass_density, electron_fraction, temperature, + {time_step_id, next_time_step_id, dg_mesh, subcell_mesh, mesh_coordinates, + active_grid, element, neighbor_data, packets_on_element, evolved_vars, + mortar_data, rest_mass_density, electron_fraction, temperature, cell_light_crossing_time, typename domain::Tags::NeighborMesh::type{}, fd_to_neighbor_fd_interpolants}); @@ -522,6 +559,8 @@ void test_send_receive_actions() { // Correct coordinates of packets sent east/south to get in the // topological coordinate of their new element. + // This is the unprocessed inbox, so the index of the closest + // point is still the same as on the old element. packet_east.coordinates[0] -= 2.0; packet_south.coordinates[1] += 2.0; // Current hack for edges/corners @@ -578,6 +617,14 @@ void test_send_receive_actions() { REQUIRE_FALSE(ActionTesting::next_action_if_ready( make_not_null(&runner), self_id)); + // Now calculates the correct index for the new points, as the + // processing of the Inbox resets the index of the closest grid point. + packet_east.index_of_closest_grid_point -= 4; + packet_south.index_of_closest_grid_point += 18; + if constexpr (Dim > 2) { + packet_south.index_of_closest_grid_point -= 81; + } + // Set up fake data coming from east neighbor { const std::optional>