Skip to content

Commit

Permalink
Take single MC step without communication
Browse files Browse the repository at this point in the history
  • Loading branch information
Francois Foucart committed Nov 14, 2024
1 parent 1bd140a commit d93d83e
Show file tree
Hide file tree
Showing 16 changed files with 267 additions and 67 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include "Evolution/Initialization/Evolution.hpp"
#include "Evolution/Initialization/Limiter.hpp"
#include "Evolution/Initialization/SetVariables.hpp"
#include "Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp"
#include "Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp"
#include "Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp"
#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
Expand Down Expand Up @@ -69,7 +70,9 @@
#include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/ConstantM1.hpp"
#include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
#include "PointwiseFunctions/Hydro/LowerSpatialFourVelocity.hpp"
#include "PointwiseFunctions/Hydro/Tags.hpp"
#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
#include "Time/Actions/AdvanceTime.hpp"
#include "Time/Actions/CleanHistory.hpp"
#include "Time/Actions/RecordTimeStepperData.hpp"
Expand Down Expand Up @@ -114,7 +117,6 @@ struct EvolutionMetavars {
using TimeStepperBase = TimeStepper;
static constexpr bool use_dg_subcell = true;

using initial_data_tag = evolution::initial_data::Tags::InitialData;
using initial_data_list =
grmhd::ValenciaDivClean::InitialData::initial_data_list;
using equation_of_state_tag = hydro::Tags::GrmhdEquationOfState;
Expand All @@ -136,6 +138,10 @@ struct EvolutionMetavars {
::Events::Tags::ObserverDetInvJacobianCompute<
Frame::ElementLogical, Frame::Inertial>>;

using analytic_variables_tags = typename system::variables_tag::tags_list;
using analytic_compute = evolution::Tags::AnalyticSolutionsCompute<
volume_dim, analytic_variables_tags, false, initial_data_list>;

struct factory_creation
: tt::ConformsTo<Options::protocols::FactoryCreation> {
using factory_classes = tmpl::map<
Expand Down Expand Up @@ -204,29 +210,38 @@ struct EvolutionMetavars {
evolution::dg::subcell::BackgroundGrVars<system, EvolutionMetavars,
true, false>>,
Actions::MutateApply<grmhd::ValenciaDivClean::subcell::SwapGrTags>,
Initialization::Actions::InitializeMCTags<system, EnergyBins,
NeutrinoSpecies>,
Initialization::Actions::AddComputeTags<tmpl::list<
hydro::Tags::LowerSpatialFourVelocityCompute,
Particles::MonteCarlo::InverseJacobianInertialToFluidCompute,
domain::Tags::JacobianCompute<4, Frame::Inertial, Frame::Fluid>>>,
// Initialization::TimeStepperHistory<EvolutionMetavars>
Parallel::Actions::TerminatePhase>;

using dg_element_array = DgElementArray<
EvolutionMetavars,
tmpl::list<
Parallel::PhaseActions<Parallel::Phase::Initialization,
initialization_actions>
//,
// Parallel::PhaseActions<
// Parallel::Phase::InitializeTimeStepperHistory,
// SelfStart::self_start_procedure<step_actions, system>>,
initialization_actions>,
//,
// Parallel::PhaseActions<
// Parallel::Phase::InitializeTimeStepperHistory,
// SelfStart::self_start_procedure<step_actions, system>>,

// Parallel::PhaseActions<Parallel::Phase::Register,
// tmpl::list<dg_registration_list,
// Parallel::Actions::TerminatePhase>>,
// Parallel::PhaseActions<Parallel::Phase::Register,
// tmpl::list<dg_registration_list,
// Parallel::Actions::TerminatePhase>>,

// Parallel::PhaseActions<
// Parallel::Phase::Evolve,
// tmpl::list<evolution::Actions::RunEventsAndTriggers,
// Actions::ChangeSlabSize, //step_actions,
// Actions::AdvanceTime,
// PhaseControl::Actions::ExecutePhaseChange>>>>;
Parallel::PhaseActions<
Parallel::Phase::Evolve,
tmpl::list<Particles::MonteCarlo::Actions::TakeTimeStep<
EnergyBins, NeutrinoSpecies>,
Parallel::Actions::TerminatePhase>>
// tmpl::list<evolution::Actions::RunEventsAndTriggers,
// Actions::ChangeSlabSize, //step_actions,
// Actions::AdvanceTime,
// PhaseControl::Actions::ExecutePhaseChange>>
>>;

struct registration
Expand All @@ -240,10 +255,10 @@ struct EvolutionMetavars {
observers::ObserverWriter<EvolutionMetavars>,
dg_element_array>;

using const_global_cache_tags =
tmpl::list<equation_of_state_tag, initial_data_tag,
Particles::MonteCarlo::Tags::InteractionRatesTable<
EnergyBins, NeutrinoSpecies>>;
using const_global_cache_tags = tmpl::list<
equation_of_state_tag, evolution::initial_data::Tags::InitialData,
Particles::MonteCarlo::Tags::InteractionRatesTable<EnergyBins,
NeutrinoSpecies>>;

static constexpr Options::String help{
"Evolve Monte Carlo transport (without coupling to hydro).\n\n"};
Expand Down
1 change: 1 addition & 0 deletions src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -5,5 +5,6 @@ spectre_target_headers(
${LIBRARY}
INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
HEADERS
InitializeMonteCarlo.hpp
TimeStepActions.hpp
)
149 changes: 149 additions & 0 deletions src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,149 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

#pragma once

#include <tuple>

#include "DataStructures/DataBox/DataBox.hpp"
#include "Domain/Tags.hpp"
#include "Evolution/Initialization/InitialData.hpp"
#include "Evolution/Particles/MonteCarlo/MortarData.hpp"
#include "Evolution/Particles/MonteCarlo/Packet.hpp"
#include "Evolution/Particles/MonteCarlo/Tags.hpp"
#include "Parallel/AlgorithmExecution.hpp"
#include "Parallel/GlobalCache.hpp"
#include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
#include "PointwiseFunctions/Hydro/Tags.hpp"
#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
#include "Utilities/TMPL.hpp"

/// \cond
namespace evolution::initial_data::Tags {
struct InitialData;
} // namespace evolution::initial_data::Tags

namespace tuples {
template <typename...>
class TaggedTuple;
} // namespace tuples

namespace Parallel {
template <typename Metavariables>
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<dim>
///
/// DataBox changes:
/// - Adds:
/// * Particles::MonteCarlo::Tags::PacketsOnElement
/// * Particles::MonteCarlo::Tags::RandomNumberGenerator
/// * Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission
/// <NeutrinoSpecies>
/// * Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>
/// * Background hydro variables
/// * Particles::MonteCarlo::Tags::MortarDataTag<dim>
///
/// - Removes: nothing
/// - Modifies: nothing
template <typename System, size_t EnergyBins, size_t NeutrinoSpecies>
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::PacketsOnElement,
Particles::MonteCarlo::Tags::RandomNumberGenerator,
Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
NeutrinoSpecies>,
Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>,
hydro_variables_tag,
Particles::MonteCarlo::Tags::MortarDataTag<dim>>;

using compute_tags = tmpl::list<>;

template <typename DbTagsList, typename... InboxTags, typename Metavariables,
typename ArrayIndex, typename ActionList,
typename ParallelComponent>
static Parallel::iterable_action_return_t apply(
db::DataBox<DbTagsList>& box,
const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
const Parallel::GlobalCache<Metavariables>& /*cache*/,
const ArrayIndex& /*array_index*/, ActionList /*meta*/,
const ParallelComponent* const /*meta*/) {
const size_t num_grid_points =

Check failure on line 83 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Debug)

variable 'num_grid_points' is not initialized

Check failure on line 83 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Release)

variable 'num_grid_points' is not initialized
db::get<evolution::dg::subcell::Tags::Mesh<dim>>(box)

Check failure on line 84 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Debug)

no member named 'dg' in namespace 'evolution'

Check failure on line 84 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Release)

no member named 'dg' in namespace 'evolution'
.number_of_grid_points();
using derived_classes =
tmpl::at<typename Metavariables::factory_creation::factory_classes,
evolution::initial_data::InitialData>;
using HydroVars = typename hydro_variables_tag::type;
call_with_dynamic_type<void, derived_classes>(
&db::get<evolution::initial_data::Tags::InitialData>(box),
[&box](const auto* const data_or_solution) {
static constexpr size_t dim = System::volume_dim;
const double initial_time = 0.0; // db::get<::Tags::Time>(box);
const size_t num_grid_points =

Check failure on line 95 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Debug)

variable 'num_grid_points' is not initialized

Check failure on line 95 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Release)

variable 'num_grid_points' is not initialized
db::get<evolution::dg::subcell::Tags::Mesh<dim>>(box)

Check failure on line 96 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Debug)

no member named 'dg' in namespace 'evolution'

Check failure on line 96 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Release)

no member named 'dg' in namespace 'evolution'
.number_of_grid_points();
const auto& inertial_coords = db::get<
evolution::dg::subcell::Tags::Coordinates<dim, Frame::Inertial>>(

Check failure on line 99 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Debug)

no member named 'dg' in namespace 'evolution'

Check failure on line 99 in src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Release)

no member named 'dg' in namespace 'evolution'
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<tmpl::list<hydro_variables_tag>>(
make_not_null(&box), std::move(hydro_variables));
});

typename Particles::MonteCarlo::Tags::PacketsOnElement::type all_packets;
Initialization::mutate_assign<
tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement>>(
make_not_null(&box), std::move(all_packets));

// Currently seeds with 0 for testing...
typename Particles::MonteCarlo::Tags::RandomNumberGenerator::type rng(0);
Initialization::mutate_assign<
tmpl::list<Particles::MonteCarlo::Tags::RandomNumberGenerator>>(
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<std::array<DataVector, NeutrinoSpecies>>(
DataVector{num_grid_points}, 1.e-12);
Initialization::mutate_assign<
tmpl::list<Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
NeutrinoSpecies>>>(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<tmpl::list<
Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>>>(
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<dim>::type;
MortarData mortar_data;
Initialization::mutate_assign<
tmpl::list<Particles::MonteCarlo::Tags::MortarDataTag<dim>>>(
make_not_null(&box), std::move(mortar_data));

return {Parallel::AlgorithmExecution::Continue, std::nullopt};
}
};

} // namespace Initialization::Actions
35 changes: 26 additions & 9 deletions src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,8 @@ struct TimeStepMutator {
using return_tags =
tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement,
Particles::MonteCarlo::Tags::RandomNumberGenerator,
Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<3>>;
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<
Expand All @@ -60,11 +61,11 @@ struct TimeStepMutator {
Frame::Inertial>,
::Tags::deriv<gr::Tags::Shift<DataVector, Dim>, tmpl::size_t<Dim>,
Frame::Inertial>,
::Tags::deriv<gr::Tags::InverseSpatialMetric<DataVector, Dim>,
tmpl::size_t<Dim>, Frame::Inertial>,
::Tags::deriv<gr::Tags::SpatialMetric<DataVector, Dim>, tmpl::size_t<Dim>,
Frame::Inertial>,
gr::Tags::SpatialMetric<DataVector, Dim, Frame::Inertial>,
gr::Tags::InverseSpatialMetric<DataVector, Dim, Frame::Inertial>,
gr::Tags::DetSpatialMetric<DataVector>,
gr::Tags::SqrtDetSpatialMetric<DataVector>,
Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>,
evolution::dg::subcell::Tags::Mesh<Dim>,
evolution::dg::subcell::Tags::Coordinates<Dim, Frame::ElementLogical>,
Expand Down Expand Up @@ -96,10 +97,10 @@ struct TimeStepMutator {

const tnsr::i<DataVector, Dim, Frame::Inertial>& d_lapse,
const tnsr::iJ<DataVector, Dim, Frame::Inertial>& d_shift,
const tnsr::iJJ<DataVector, Dim, Frame::Inertial>& d_inv_spatial_metric,
const tnsr::ijj<DataVector, Dim, Frame::Inertial>& d_spatial_metric,
const tnsr::ii<DataVector, Dim, Frame::Inertial>& spatial_metric,
const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric,
const Scalar<DataVector>& determinant_spatial_metric,
const Scalar<DataVector>& sqrt_determinant_spatial_metric,
const Scalar<DataVector>& cell_light_crossing_time, const Mesh<Dim>& mesh,
const tnsr::I<DataVector, Dim, Frame::ElementLogical>& mesh_coordinates,
const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
Expand All @@ -117,8 +118,8 @@ struct TimeStepMutator {
const size_t num_ghost_zones = 1;
// Get information stored in various databox containers in
// the format expected by take_time_step_on_element
const double start_time = current_step_id.step_time().value();
const double end_time = next_step_id.step_time().value();
const double start_time = 0.0; // current_step_id.step_time().value();
const double end_time = 0.1; // next_step_id.step_time().value();
Scalar<DataVector> det_jacobian_logical_to_inertial(lapse);
get(det_jacobian_logical_to_inertial) =
1.0 / get(det_inverse_jacobian_logical_to_inertial);
Expand All @@ -131,14 +132,30 @@ struct TimeStepMutator {
const DirectionalIdMap<Dim, std::optional<DataVector>>&
cell_light_crossing_time_ghost = mortar_data.cell_light_crossing_time;

tnsr::iJJ<DataVector, 3, Frame::Inertial> d_inv_spatial_metric =
make_with_value<tnsr::iJJ<DataVector, 3, Frame::Inertial>>(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<EnergyBins, NeutrinoSpecies> templated_functions;
templated_functions.take_time_step_on_element(
packets, random_number_generator, single_packet_energy, start_time,
end_time, equation_of_state, interaction_table, electron_fraction,
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,
Expand Down
4 changes: 2 additions & 2 deletions src/Evolution/Particles/MonteCarlo/CellVolume.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,13 +12,13 @@ namespace Particles::MonteCarlo {
void cell_proper_four_volume_finite_difference(
const gsl::not_null<Scalar<DataVector>*> cell_proper_four_volume,
const Scalar<DataVector>& lapse,
const Scalar<DataVector>& determinant_spatial_metric,
const Scalar<DataVector>& sqrt_determinant_spatial_metric,
const double time_step, const Mesh<3>& mesh,
const Scalar<DataVector>& det_jacobian_logical_to_inertial) {
const double cell_logical_volume =
8.0 / static_cast<double>(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);
}

Expand Down
11 changes: 5 additions & 6 deletions src/Evolution/Particles/MonteCarlo/CellVolume.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<Scalar<DataVector>* > cell_proper_four_volume,
const Scalar<DataVector>& lapse,
const Scalar<DataVector>& determinant_spatial_metric,
double time_step,
const Mesh<3>& mesh,
const Scalar<DataVector>& det_jacobian_logical_to_inertial);
gsl::not_null<Scalar<DataVector>*> cell_proper_four_volume,
const Scalar<DataVector>& lapse,
const Scalar<DataVector>& sqrt_determinant_spatial_metric, double time_step,
const Mesh<3>& mesh,
const Scalar<DataVector>& 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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -155,11 +155,13 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::
// We have all we need to calculate beta here. We take the largest of the
// two variations
get(beta)[i] = fabs(get(neutrino_energy_1)[i] - get(neutrino_energy_0)[i]) /
fabs(get(fluid_energy_1)[i] - get(fluid_energy_0)[i]);
(fabs(get(fluid_energy_1)[i] - get(fluid_energy_0)[i]) +
1.e-16);
const double beta_lepton =
fabs(get(neutrino_lepton_number_1)[i] -
get(neutrino_lepton_number_0)[i]) /
fabs(get(fluid_lepton_number_1)[i] - get(fluid_lepton_number_0)[i]);
(fabs(get(fluid_lepton_number_1)[i] - get(fluid_lepton_number_0)[i]) +
1.e-16);
get(beta)[i] =
std::max(get(beta)[i], beta_lepton);
}
Expand Down
Loading

0 comments on commit d93d83e

Please sign in to comment.