Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <optional>
7 : #include <tuple>
8 : #include <utility>
9 :
10 : #include "DataStructures/DataBox/DataBox.hpp"
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/Index.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Domain/Structure/Direction.hpp"
15 : #include "Domain/Structure/DirectionalId.hpp"
16 : #include "Domain/Structure/DirectionalIdMap.hpp"
17 : #include "Domain/Structure/Element.hpp"
18 : #include "Domain/Structure/ElementId.hpp"
19 : #include "Domain/Structure/OrientationMapHelpers.hpp"
20 : #include "Domain/Tags.hpp"
21 : #include "Domain/Tags/NeighborMesh.hpp"
22 : #include "Evolution/DgSubcell/ActiveGrid.hpp"
23 : #include "Evolution/DgSubcell/GhostData.hpp"
24 : #include "Evolution/DgSubcell/SliceData.hpp"
25 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
26 : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
27 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
28 : #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp"
29 : #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
30 : #include "Evolution/Particles/MonteCarlo/MortarData.hpp"
31 : #include "Evolution/Particles/MonteCarlo/Tags.hpp"
32 : #include "Parallel/AlgorithmExecution.hpp"
33 : #include "Parallel/GlobalCache.hpp"
34 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
35 : #include "Time/TimeStepId.hpp"
36 : #include "Utilities/ErrorHandling/Assert.hpp"
37 : #include "Utilities/Gsl.hpp"
38 :
39 : /// \cond
40 : namespace Tags {
41 : struct TimeStepId;
42 : } // namespace Tags
43 : /// \endcond
44 :
45 : namespace Particles::MonteCarlo::Actions {
46 :
47 : /// Mutator to get required volume data for communication
48 : /// before a MC step; i.e. data sent from live points
49 : /// to ghost points in neighbors.
50 1 : struct GhostDataMutatorPreStep {
51 0 : using return_tags = tmpl::list<>;
52 0 : using argument_tags = tmpl::list<
53 : hydro::Tags::RestMassDensity<DataVector>,
54 : hydro::Tags::ElectronFraction<DataVector>,
55 : hydro::Tags::Temperature<DataVector>,
56 : Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>>;
57 0 : static const size_t number_of_vars = 4;
58 :
59 0 : static DataVector apply(const Scalar<DataVector>& rest_mass_density,
60 : const Scalar<DataVector>& electron_fraction,
61 : const Scalar<DataVector>& temperature,
62 : const Scalar<DataVector>& cell_light_crossing_time) {
63 : const size_t dv_size = get(rest_mass_density).size();
64 : DataVector buffer{dv_size * number_of_vars};
65 : std::copy(
66 : get(rest_mass_density).data(),
67 : std::next(get(rest_mass_density).data(), static_cast<int>(dv_size)),
68 : buffer.data());
69 : std::copy(
70 : get(electron_fraction).data(),
71 : std::next(get(electron_fraction).data(), static_cast<int>(dv_size)),
72 : std::next(buffer.data(), static_cast<int>(dv_size)));
73 : std::copy(get(temperature).data(),
74 : std::next(get(temperature).data(), static_cast<int>(dv_size)),
75 : std::next(buffer.data(), static_cast<int>(dv_size * 2)));
76 : std::copy(get(cell_light_crossing_time).data(),
77 : std::next(get(cell_light_crossing_time).data(),
78 : static_cast<int>(dv_size)),
79 : std::next(buffer.data(), static_cast<int>(dv_size * 3)));
80 : return buffer;
81 : }
82 : };
83 :
84 : /// Mutator that returns packets currently in ghost zones in a
85 : /// DirectionMap<Dim,std::vector<Particles::MonteCarlo::Packet>>
86 : /// and remove them from the list of packets of the current
87 : /// element
88 : template <size_t Dim>
89 1 : struct GhostDataMcPackets {
90 0 : using return_tags = tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement>;
91 0 : using argument_tags = tmpl::list<::domain::Tags::Element<Dim>>;
92 :
93 0 : static DirectionMap<Dim, std::vector<Particles::MonteCarlo::Packet>> apply(
94 : const gsl::not_null<std::vector<Particles::MonteCarlo::Packet>*> packets,
95 : const Element<Dim>& element) {
96 : DirectionMap<Dim, std::vector<Particles::MonteCarlo::Packet>> output{};
97 : for (const auto& [direction, neighbors_in_direction] :
98 : element.neighbors()) {
99 : output[direction] = std::vector<Particles::MonteCarlo::Packet>{};
100 : }
101 : // Loop over packets. If a packet is out of the current element,
102 : // remove it from the element and copy it to the appropriate
103 : // neighbor list if it exists.
104 : // Note: At this point, we only handle simple boundaries
105 : // (one neighbor per direction, with the logical coordinates
106 : // being transformed by adding/removing 2). We assume that
107 : // packets that go out of the element but do not match any
108 : // neighbor have reached the domain boundary and can be removed.
109 : size_t n_packets = packets->size();
110 : for (size_t p = 0; p < n_packets; p++) {
111 : Packet& packet = (*packets)[p];
112 : std::optional<Direction<Dim>> max_distance_direction = std::nullopt;
113 : double max_distance = 0.0;
114 : // TO DO: Deal with dimensions higher than the grid dimension
115 : // e.g. for axisymmetric evolution
116 : for (size_t d = 0; d < Dim; d++) {
117 : if (packet.coordinates.get(d) < -1.0) {
118 : const double distance = -1.0 - packet.coordinates.get(d);
119 : if (max_distance < distance) {
120 : max_distance = distance;
121 : max_distance_direction = Direction<Dim>(d, Side::Lower);
122 : }
123 : }
124 : if (packet.coordinates.get(d) > 1.0) {
125 : const double distance = packet.coordinates.get(d) - 1.0;
126 : if (max_distance < distance) {
127 : max_distance = distance;
128 : max_distance_direction = Direction<Dim>(d, Side::Upper);
129 : }
130 : }
131 : }
132 : // If a packet should be moved along an edge / corner, we move
133 : // it into the direct neighbor it is effectively closest to
134 : // (in topological coordinates). This may need improvements.
135 : if (max_distance_direction != std::nullopt) {
136 : const size_t d = max_distance_direction.value().dimension();
137 : const Side side = max_distance_direction.value().side();
138 : packet.coordinates.get(d) += (side == Side::Lower) ? (2.0) : (-2.0);
139 : if (output.contains(max_distance_direction.value())) {
140 : output[max_distance_direction.value()].push_back(packet);
141 : }
142 : std::swap((*packets)[p], (*packets)[n_packets - 1]);
143 : packets->pop_back();
144 : p--;
145 : n_packets--;
146 : }
147 : }
148 : return output;
149 : }
150 : };
151 :
152 : /// Action responsible for the Send operation of ghost
153 : /// zone communication for Monte-Carlo transport.
154 : /// If CommStep == PreStep, this sends the fluid and
155 : /// metric variables needed for evolution of MC packets.
156 : /// If CommStep == PostStep, this sends packets that
157 : /// have moved from one element to another.
158 : /// The current code does not yet deal with data required
159 : /// to calculate the backreaction on the fluid.
160 : template <size_t Dim, bool LocalTimeStepping,
161 : Particles::MonteCarlo::CommunicationStep CommStep>
162 1 : struct SendDataForMcCommunication {
163 0 : using inbox_tags =
164 : tmpl::list<Particles::MonteCarlo::McGhostZoneDataInboxTag<Dim, CommStep>>;
165 :
166 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
167 : typename ActionList, typename ParallelComponent,
168 : typename Metavariables>
169 0 : static Parallel::iterable_action_return_t apply(
170 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
171 : Parallel::GlobalCache<Metavariables>& cache,
172 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
173 : const ParallelComponent* const /*meta*/) {
174 : static_assert(not LocalTimeStepping,
175 : "Monte Carlo is not tested with local time stepping.");
176 :
177 : ASSERT(db::get<evolution::dg::subcell::Tags::ActiveGrid>(box) ==
178 : evolution::dg::subcell::ActiveGrid::Subcell,
179 : "The SendDataForReconstructionPreStep action in MC "
180 : "can only be called when Subcell is the active scheme.");
181 : db::mutate<Particles::MonteCarlo::Tags::McGhostZoneDataTag<Dim>>(
182 : [](const auto ghost_data_ptr) {
183 : // Clear the previous neighbor data and add current local data
184 : ghost_data_ptr->clear();
185 : },
186 : make_not_null(&box));
187 :
188 : auto& receiver_proxy =
189 : Parallel::get_parallel_component<ParallelComponent>(cache);
190 : const size_t ghost_zone_size = 1;
191 :
192 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
193 : const Mesh<Dim>& subcell_mesh =
194 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(box);
195 : const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box);
196 :
197 : // Fill volume data that should be fed to neighbor ghost zones
198 : DataVector volume_data_to_slice =
199 : CommStep == Particles::MonteCarlo::CommunicationStep::PreStep
200 : ? db::mutate_apply(GhostDataMutatorPreStep{}, make_not_null(&box))
201 : : DataVector{};
202 : const DirectionMap<Dim, DataVector> all_sliced_data =
203 : CommStep == Particles::MonteCarlo::CommunicationStep::PreStep
204 : ? evolution::dg::subcell::slice_data(
205 : volume_data_to_slice, subcell_mesh.extents(), ghost_zone_size,
206 : element.internal_boundaries(), 0,
207 : db::get<evolution::dg::subcell::Tags::
208 : InterpolatorsFromFdToNeighborFd<Dim>>(box))
209 : : DirectionMap<Dim, DataVector>{};
210 :
211 : const DirectionMap<Dim, std::vector<Particles::MonteCarlo::Packet>>
212 : all_packets_ghost_zone =
213 : CommStep == Particles::MonteCarlo::CommunicationStep::PostStep
214 : ? db::mutate_apply(GhostDataMcPackets<Dim>{},
215 : make_not_null(&box))
216 : : DirectionMap<Dim,
217 : std::vector<Particles::MonteCarlo::Packet>>{};
218 :
219 : // TO DO:
220 : // Need to deal with coupling data, which is brought back from
221 : // neighbor's ghost zones to the processor with the live cell!
222 :
223 : for (const auto& [direction, neighbors_in_direction] :
224 : element.neighbors()) {
225 : const auto& orientation = neighbors_in_direction.orientation();
226 : const auto direction_from_neighbor = orientation(direction.opposite());
227 : for (const ElementId<Dim>& neighbor : neighbors_in_direction) {
228 : std::optional<std::vector<Particles::MonteCarlo::Packet>>
229 : packets_to_send = std::nullopt;
230 : DataVector subcell_data_to_send{};
231 :
232 : if (CommStep == Particles::MonteCarlo::CommunicationStep::PreStep) {
233 : const auto& sliced_data_in_direction = all_sliced_data.at(direction);
234 : subcell_data_to_send = DataVector{sliced_data_in_direction.size()};
235 : std::copy(sliced_data_in_direction.begin(),
236 : sliced_data_in_direction.end(),
237 : subcell_data_to_send.begin());
238 : }
239 : if (CommStep == Particles::MonteCarlo::CommunicationStep::PostStep) {
240 : packets_to_send = all_packets_ghost_zone.at(direction);
241 : if (packets_to_send.value().empty()) {
242 : packets_to_send = std::nullopt;
243 : }
244 : }
245 :
246 : McGhostZoneData<Dim> data{subcell_data_to_send, packets_to_send};
247 :
248 : Parallel::receive_data<
249 : Particles::MonteCarlo::McGhostZoneDataInboxTag<Dim, CommStep>>(
250 : receiver_proxy[neighbor], time_step_id,
251 : std::pair{DirectionalId<Dim>{direction_from_neighbor, element.id()},
252 : std::move(data)});
253 : }
254 : }
255 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
256 : }
257 : };
258 :
259 : /// Action responsible for the Receive operation of ghost
260 : /// zone communication for Monte-Carlo transport.
261 : /// If CommStep == PreStep, this gets the fluid and
262 : /// metric variables needed for evolution of MC packets.
263 : /// If CommStep == PostStep, this gets packets that
264 : /// have moved into the current element.
265 : /// The current code does not yet deal with data required
266 : /// to calculate the backreaction on the fluid.
267 : template <size_t Dim, CommunicationStep CommStep>
268 1 : struct ReceiveDataForMcCommunication {
269 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
270 : typename ActionList, typename ParallelComponent,
271 : typename Metavariables>
272 0 : static Parallel::iterable_action_return_t apply(
273 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
274 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
275 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
276 : const ParallelComponent* const /*meta*/) {
277 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
278 : const auto number_of_expected_messages = element.neighbors().size();
279 : if (UNLIKELY(number_of_expected_messages == 0)) {
280 : // We have no neighbors, so just continue without doing any work
281 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
282 : }
283 :
284 : const auto& current_time_step_id = db::get<::Tags::TimeStepId>(box);
285 : std::map<TimeStepId, DirectionalIdMap<
286 : Dim, Particles::MonteCarlo::McGhostZoneData<Dim>>>&
287 : inbox = tuples::get<Particles::MonteCarlo::McGhostZoneDataInboxTag<
288 : Metavariables::volume_dim, CommStep>>(inboxes);
289 : const auto& received = inbox.find(current_time_step_id);
290 : // Check we have at least some data from correct time, and then check that
291 : // we have received all data
292 : if (received == inbox.end() or
293 : received->second.size() != number_of_expected_messages) {
294 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
295 : }
296 :
297 : // Now that we have received all the data, copy it over as needed.
298 : DirectionalIdMap<Dim, Particles::MonteCarlo::McGhostZoneData<Dim>>
299 : received_data = std::move(inbox[current_time_step_id]);
300 : inbox.erase(current_time_step_id);
301 :
302 : // Copy the received PreStep data in the MortarData container
303 : if (CommStep == Particles::MonteCarlo::CommunicationStep::PreStep) {
304 : db::mutate<Particles::MonteCarlo::Tags::MortarDataTag<Dim>>(
305 : [&element,
306 : &received_data](const gsl::not_null<MortarData<Dim>*> mortar_data) {
307 : for (const auto& [direction, neighbors_in_direction] :
308 : element.neighbors()) {
309 : for (const auto& neighbor : neighbors_in_direction) {
310 : DirectionalId<Dim> directional_element_id{direction, neighbor};
311 : const DataVector& received_data_direction =
312 : received_data[directional_element_id]
313 : .ghost_zone_hydro_variables;
314 : REQUIRE(received_data[directional_element_id]
315 : .packets_entering_this_element == std::nullopt);
316 : if (mortar_data->rest_mass_density[directional_element_id] ==
317 : std::nullopt) {
318 : continue;
319 : }
320 : const size_t mortar_data_size =
321 : mortar_data->rest_mass_density[directional_element_id]
322 : .value()
323 : .size();
324 : std::copy(received_data_direction.data(),
325 : std::next(received_data_direction.data(),
326 : static_cast<int>(mortar_data_size)),
327 : mortar_data->rest_mass_density[directional_element_id]
328 : .value()
329 : .data());
330 : std::copy(std::next(received_data_direction.data(),
331 : static_cast<int>(mortar_data_size)),
332 : std::next(received_data_direction.data(),
333 : static_cast<int>(mortar_data_size * 2)),
334 : mortar_data->electron_fraction[directional_element_id]
335 : .value()
336 : .data());
337 : std::copy(std::next(received_data_direction.data(),
338 : static_cast<int>(mortar_data_size * 2)),
339 : std::next(received_data_direction.data(),
340 : static_cast<int>(mortar_data_size * 3)),
341 : mortar_data->temperature[directional_element_id]
342 : .value()
343 : .data());
344 : std::copy(std::next(received_data_direction.data(),
345 : static_cast<int>(mortar_data_size * 3)),
346 : std::next(received_data_direction.data(),
347 : static_cast<int>(mortar_data_size * 4)),
348 : mortar_data
349 : ->cell_light_crossing_time[directional_element_id]
350 : .value()
351 : .data());
352 : }
353 : }
354 : },
355 : make_not_null(&box));
356 : } else {
357 : // TO DO: Deal with data coupling neutrinos back to fluid evolution
358 : db::mutate<Particles::MonteCarlo::Tags::PacketsOnElement>(
359 : [&element, &received_data](
360 : const gsl::not_null<std::vector<Particles::MonteCarlo::Packet>*>
361 : packet_list) {
362 : for (const auto& [direction, neighbors_in_direction] :
363 : element.neighbors()) {
364 : for (const auto& neighbor : neighbors_in_direction) {
365 : DirectionalId<Dim> directional_element_id{direction, neighbor};
366 : const std::optional<std::vector<Particles::MonteCarlo::Packet>>&
367 : received_data_packets =
368 : received_data[directional_element_id]
369 : .packets_entering_this_element;
370 : // Temporary: currently no data for coupling to the fluid
371 : REQUIRE(received_data[directional_element_id]
372 : .ghost_zone_hydro_variables.size() == 0);
373 : if (received_data_packets == std::nullopt) {
374 : continue;
375 : } else {
376 : const size_t n_packets = received_data_packets.value().size();
377 : for (size_t p = 0; p < n_packets; p++) {
378 : packet_list->push_back(received_data_packets.value()[p]);
379 : }
380 : }
381 : }
382 : }
383 : },
384 : make_not_null(&box));
385 : }
386 :
387 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
388 : }
389 : };
390 :
391 : } // namespace Particles::MonteCarlo::Actions
|
