Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <random>
7 : #include <tuple>
8 :
9 : #include "DataStructures/DataBox/DataBox.hpp"
10 : #include "Domain/Structure/Element.hpp"
11 : #include "Domain/Tags.hpp"
12 : #include "Evolution/DgSubcell/ActiveGrid.hpp"
13 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
14 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
15 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
16 : #include "Evolution/Initialization/InitialData.hpp"
17 : #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
18 : #include "Evolution/Particles/MonteCarlo/MonteCarloOptions.hpp"
19 : #include "Evolution/Particles/MonteCarlo/MortarData.hpp"
20 : #include "Evolution/Particles/MonteCarlo/Packet.hpp"
21 : #include "Evolution/Particles/MonteCarlo/Tags.hpp"
22 : #include "Parallel/AlgorithmExecution.hpp"
23 : #include "Parallel/GlobalCache.hpp"
24 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
25 : #include "PointwiseFunctions/Hydro/Tags.hpp"
26 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
27 : #include "Utilities/TMPL.hpp"
28 :
29 : /// \cond
30 : namespace evolution::initial_data::Tags {
31 : struct InitialData;
32 : } // namespace evolution::initial_data::Tags
33 :
34 : namespace tuples {
35 : template <typename...>
36 : class TaggedTuple;
37 : } // namespace tuples
38 :
39 : namespace Parallel {
40 : template <typename Metavariables>
41 : class GlobalCache;
42 : } // namespace Parallel
43 : /// \endcond
44 :
45 : namespace Initialization::Actions {
46 :
47 : /// \ingroup InitializationGroup
48 : /// \brief Allocate variables needed for evolution of Monte Carlo transport
49 : ///
50 : /// Uses:
51 : /// - evolution::dg::subcell::Tags::Mesh<dim>
52 : /// - evolution::dg::subcell::Tags::Coordinates<dim, Frame::Inertial>
53 : /// - evolution::dg::subcell::Tags::ActiveGrid
54 : /// - ::Tags::Time
55 : /// - evolution::initial_data::Tags::InitialData
56 : /// - Particles::MonteCarlo::Tags::MonteCarloOptions<EnergyBins,
57 : /// NeutrinoSpecies>
58 : /// - domain::Tags::Element<dim>
59 : ///
60 : /// DataBox changes:
61 : /// - Adds:
62 : /// * Particles::MonteCarlo::Tags::PacketsOnElement
63 : /// * Particles::MonteCarlo::Tags::RandomNumberGenerator
64 : /// * Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
65 : /// NeutrinoSpecies>
66 : /// * Background hydro variables
67 : /// * Particles::MonteCarlo::Tags::CouplingTildeTau<DataVector>
68 : /// * Particles::MonteCarlo::Tags::CouplingTildeRhoYe<DataVector>
69 : /// * Particles::MonteCarlo::Tags::CouplingTildeS<DataVector,dim>
70 : /// * Particles::MonteCarlo::Tags::MortarDataTag<dim>
71 : /// * Particles::MonteCarlo::Tags::GhostZoneCouplingData<dim>
72 : /// * Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>
73 : ///
74 : /// - Removes: nothing
75 : /// - Modifies: nothing
76 : template <typename System, size_t EnergyBins, size_t NeutrinoSpecies,
77 : bool InitializeBackground>
78 1 : struct InitializeMCTags {
79 : public:
80 0 : using hydro_variables_tag = typename System::hydro_variables_tag;
81 :
82 0 : static constexpr size_t dim = System::volume_dim;
83 0 : using simple_tags =
84 : tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement,
85 : Particles::MonteCarlo::Tags::RandomNumberGenerator,
86 : Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
87 : NeutrinoSpecies>,
88 : hydro_variables_tag,
89 : Particles::MonteCarlo::Tags::CouplingTildeTau<DataVector>,
90 : Particles::MonteCarlo::Tags::CouplingTildeRhoYe<DataVector>,
91 : Particles::MonteCarlo::Tags::CouplingTildeS<DataVector, dim>,
92 : Particles::MonteCarlo::Tags::MortarDataTag<dim>,
93 : Particles::MonteCarlo::Tags::GhostZoneCouplingDataTag<dim>,
94 : Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>,
95 : evolution::dg::subcell::Tags::ActiveGrid>;
96 :
97 0 : using compute_tags = tmpl::list<>;
98 :
99 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
100 : typename ArrayIndex, typename ActionList,
101 : typename ParallelComponent>
102 0 : static Parallel::iterable_action_return_t apply(
103 : db::DataBox<DbTagsList>& box,
104 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
105 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
106 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
107 : const ParallelComponent* const /*meta*/) {
108 : if (db::get<evolution::dg::subcell::Tags::ActiveGrid>(box) !=
109 : evolution::dg::subcell::ActiveGrid::Subcell) {
110 : ERROR("MC requires all elements to use Subcell");
111 : }
112 : const Mesh<dim>& mesh =
113 : db::get<evolution::dg::subcell::Tags::Mesh<dim>>(box);
114 : const size_t num_grid_points = mesh.number_of_grid_points();
115 : // Number of ghost zones for MC is assumed to be 1 for now.
116 : const size_t num_ghost_zones = 1;
117 : size_t mesh_size_with_ghost_zones = 1;
118 : for (size_t d = 0; d < dim; d++) {
119 : mesh_size_with_ghost_zones *= (mesh.extents()[d] + 2 * num_ghost_zones);
120 : }
121 : const DataVector zero_dv_with_ghost_zones(mesh_size_with_ghost_zones, 0.0);
122 : const Scalar<DataVector> zero_scalar_with_ghost_zones =
123 : make_with_value<Scalar<DataVector>>(zero_dv_with_ghost_zones, 0.0);
124 : const tnsr::i<DataVector, dim, Frame::Inertial> zero_tnsr_with_ghost_zones =
125 : make_with_value<tnsr::i<DataVector, 3, Frame::Inertial>>(
126 : zero_dv_with_ghost_zones, 0.0);
127 : if constexpr (InitializeBackground) {
128 : using derived_classes =
129 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
130 : evolution::initial_data::InitialData>;
131 : using HydroVars = typename hydro_variables_tag::type;
132 : call_with_dynamic_type<void, derived_classes>(
133 : &db::get<evolution::initial_data::Tags::InitialData>(box),
134 : [&box, &num_grid_points](const auto* const data_or_solution) {
135 : static constexpr size_t dim = System::volume_dim;
136 : const double initial_time = db::get<::Tags::Time>(box);
137 : const auto& inertial_coords =
138 : db::get<evolution::dg::subcell::Tags::Coordinates<
139 : dim, Frame::Inertial>>(box);
140 : // Get hydro variables
141 : HydroVars hydro_variables{num_grid_points};
142 : hydro_variables.assign_subset(
143 : evolution::Initialization::initial_data(
144 : *data_or_solution, inertial_coords, initial_time,
145 : typename hydro_variables_tag::tags_list{}));
146 : Initialization::mutate_assign<tmpl::list<hydro_variables_tag>>(
147 : make_not_null(&box), std::move(hydro_variables));
148 : });
149 : }
150 :
151 : Initialization::mutate_assign<
152 : tmpl::list<Particles::MonteCarlo::Tags::CouplingTildeTau<DataVector>>>(
153 : make_not_null(&box), zero_scalar_with_ghost_zones);
154 : Initialization::mutate_assign<tmpl::list<
155 : Particles::MonteCarlo::Tags::CouplingTildeRhoYe<DataVector>>>(
156 : make_not_null(&box), zero_scalar_with_ghost_zones);
157 : Initialization::mutate_assign<tmpl::list<
158 : Particles::MonteCarlo::Tags::CouplingTildeS<DataVector, dim>>>(
159 : make_not_null(&box), zero_tnsr_with_ghost_zones);
160 :
161 : // Read global options for Monte-Carlo evolution
162 : const auto mc_options = db::get<
163 : Particles::MonteCarlo::Tags::MonteCarloOptions<NeutrinoSpecies>>(box);
164 : const auto& initial_packet_energy = mc_options.get_initial_packet_energy();
165 :
166 : typename Particles::MonteCarlo::Tags::PacketsOnElement::type all_packets;
167 : Initialization::mutate_assign<
168 : tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement>>(
169 : make_not_null(&box), std::move(all_packets));
170 :
171 : const unsigned long seed = std::random_device{}();
172 : typename Particles::MonteCarlo::Tags::RandomNumberGenerator::type rng(seed);
173 :
174 : Initialization::mutate_assign<
175 : tmpl::list<Particles::MonteCarlo::Tags::RandomNumberGenerator>>(
176 : make_not_null(&box), std::move(rng));
177 :
178 : // Initial energy of packets, read from MC options
179 : typename Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
180 : NeutrinoSpecies>::type packet_energy_at_emission =
181 : make_with_value<std::array<DataVector, NeutrinoSpecies>>(
182 : DataVector{num_grid_points}, 0.0);
183 : for (size_t s = 0; s < NeutrinoSpecies; s++) {
184 : packet_energy_at_emission[s] = initial_packet_energy[s];
185 : }
186 : Initialization::mutate_assign<
187 : tmpl::list<Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
188 : NeutrinoSpecies>>>(make_not_null(&box),
189 : std::move(packet_energy_at_emission));
190 :
191 : // Initialize mortar data and coupling data.
192 : // Currently assumes a single neighbor on each face (i.e. no h-refinement)
193 : using MortarData =
194 : typename Particles::MonteCarlo::Tags::MortarDataTag<dim>::type;
195 : MortarData mortar_data;
196 : using CouplingData =
197 : typename Particles::MonteCarlo::Tags::GhostZoneCouplingDataTag<
198 : dim>::type;
199 : CouplingData coupling_data;
200 : const Element<dim>& element = db::get<::domain::Tags::Element<dim>>(box);
201 : for (const auto& [direction, neighbors] : element.neighbors()) {
202 : const size_t sliced_mesh_size =
203 : mesh.slice_away(direction.dimension()).number_of_grid_points();
204 : const DataVector zero_dv_slice(sliced_mesh_size, 0.0);
205 : const Index<dim - 1> sliced_mesh_extents =
206 : mesh.slice_away(direction.dimension()).extents();
207 : size_t sliced_mesh_size_with_ghost_zone = 1;
208 : for (size_t d = 0; d < dim - 1; d++) {
209 : sliced_mesh_size_with_ghost_zone *= ( sliced_mesh_extents[d]
210 : + 2 * num_ghost_zones );
211 : }
212 : const DataVector zero_dv_ghost_zones(sliced_mesh_size_with_ghost_zone,
213 : 0.0);
214 : const tnsr::i<DataVector, dim, Frame::Inertial> zero_tnsr_ghost_zones =
215 : make_with_value<tnsr::i<DataVector, 3, Frame::Inertial>>(
216 : zero_dv_ghost_zones, 0.0);
217 :
218 : for (const auto& neighbor : neighbors) {
219 : const DirectionalId<dim> mortar_id{direction, neighbor};
220 : mortar_data.rest_mass_density.emplace(mortar_id, zero_dv_slice);
221 : mortar_data.electron_fraction.emplace(mortar_id, zero_dv_slice);
222 : mortar_data.temperature.emplace(mortar_id, zero_dv_slice);
223 : mortar_data.cell_light_crossing_time.emplace(mortar_id, zero_dv_slice);
224 : coupling_data.coupling_tilde_tau.emplace(mortar_id,
225 : zero_dv_ghost_zones);
226 : coupling_data.coupling_tilde_rho_ye.emplace(mortar_id,
227 : zero_dv_ghost_zones);
228 : coupling_data.coupling_tilde_s.emplace(mortar_id,
229 : zero_tnsr_ghost_zones);
230 : }
231 : }
232 : Initialization::mutate_assign<
233 : tmpl::list<Particles::MonteCarlo::Tags::MortarDataTag<dim>>>(
234 : make_not_null(&box), std::move(mortar_data));
235 : Initialization::mutate_assign<
236 : tmpl::list<Particles::MonteCarlo::Tags::GhostZoneCouplingDataTag<dim>>>(
237 : make_not_null(&box), std::move(coupling_data));
238 :
239 : using GhostZoneData =
240 : typename Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>::type;
241 : GhostZoneData ghost_zone_data{};
242 : Initialization::mutate_assign<
243 : tmpl::list<Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>>>(
244 : make_not_null(&box), std::move(ghost_zone_data));
245 :
246 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
247 : }
248 : };
249 :
250 : } // namespace Initialization::Actions
|
