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/Coordinates.hpp"
27 : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
28 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
29 : #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationStep.hpp"
30 : #include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
31 : #include "Evolution/Particles/MonteCarlo/MortarData.hpp"
32 : #include "Evolution/Particles/MonteCarlo/Tags.hpp"
33 : #include "Parallel/AlgorithmExecution.hpp"
34 : #include "Parallel/GlobalCache.hpp"
35 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
36 : #include "Time/TimeStepId.hpp"
37 : #include "Utilities/ErrorHandling/Assert.hpp"
38 : #include "Utilities/Gsl.hpp"
39 :
40 : /// \cond
41 : namespace Tags {
42 : struct TimeStepId;
43 : } // namespace Tags
44 : /// \endcond
45 :
46 : namespace Particles::MonteCarlo::Actions {
47 :
48 : /// Mutator to get required volume data for communication
49 : /// before a MC step; i.e. data sent from live points
50 : /// to ghost points in neighbors.
51 1 : struct GhostDataMutatorPreStep {
52 0 : using return_tags = tmpl::list<>;
53 0 : using argument_tags = tmpl::list<
54 : hydro::Tags::RestMassDensity<DataVector>,
55 : hydro::Tags::ElectronFraction<DataVector>,
56 : hydro::Tags::Temperature<DataVector>,
57 : Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>>;
58 0 : static const size_t number_of_vars = 4;
59 :
60 0 : static DataVector apply(const Scalar<DataVector>& rest_mass_density,
61 : const Scalar<DataVector>& electron_fraction,
62 : const Scalar<DataVector>& temperature,
63 : const Scalar<DataVector>& cell_light_crossing_time) {
64 : const size_t dv_size = get(rest_mass_density).size();
65 : DataVector buffer{dv_size * number_of_vars};
66 : std::copy(
67 : get(rest_mass_density).data(),
68 : std::next(get(rest_mass_density).data(), static_cast<int>(dv_size)),
69 : buffer.data());
70 : std::copy(
71 : get(electron_fraction).data(),
72 : std::next(get(electron_fraction).data(), static_cast<int>(dv_size)),
73 : std::next(buffer.data(), static_cast<int>(dv_size)));
74 : std::copy(get(temperature).data(),
75 : std::next(get(temperature).data(), static_cast<int>(dv_size)),
76 : std::next(buffer.data(), static_cast<int>(dv_size * 2)));
77 : std::copy(get(cell_light_crossing_time).data(),
78 : std::next(get(cell_light_crossing_time).data(),
79 : static_cast<int>(dv_size)),
80 : std::next(buffer.data(), static_cast<int>(dv_size * 3)));
81 : return buffer;
82 : }
83 : };
84 :
85 : /// Mutator to get required volume data for communication
86 : /// after a MC step; i.e. data sent from ghost points
87 : /// in neighbors to live points evolving the fluid. The
88 : /// data contains information about the back reaction of
89 : /// neutrinos on the fluid
90 : template <size_t Dim>
91 1 : struct GhostDataMutatorPostStep {
92 0 : using return_tags = tmpl::list<>;
93 0 : using argument_tags =
94 : tmpl::list<Particles::MonteCarlo::Tags::CouplingTildeTau<DataVector>,
95 : Particles::MonteCarlo::Tags::CouplingTildeRhoYe<DataVector>,
96 : Particles::MonteCarlo::Tags::CouplingTildeS<DataVector, Dim>>;
97 0 : static const size_t number_of_components = 2 + Dim;
98 :
99 0 : static DataVector apply(const Scalar<DataVector>& coupling_tilde_tau,
100 : const Scalar<DataVector>& coupling_tilde_rho_ye,
101 : const tnsr::i<DataVector, Dim>& coupling_tilde_s) {
102 : const size_t dv_size = get(coupling_tilde_tau).size();
103 : DataVector buffer{dv_size * number_of_components};
104 : std::copy(
105 : get(coupling_tilde_tau).data(),
106 : std::next(get(coupling_tilde_tau).data(), static_cast<int>(dv_size)),
107 : buffer.data());
108 : std::copy(
109 : get(coupling_tilde_rho_ye).data(),
110 : std::next(get(coupling_tilde_rho_ye).data(), static_cast<int>(dv_size)),
111 : std::next(buffer.data(), static_cast<int>(dv_size)));
112 : for (size_t d = 0; d < Dim; d++) {
113 : std::copy(
114 : coupling_tilde_s.get(d).data(),
115 : std::next(coupling_tilde_s.get(d).data(), static_cast<int>(dv_size)),
116 : std::next(buffer.data(), static_cast<int>(dv_size * (2 + d))));
117 : }
118 : return buffer;
119 : }
120 : };
121 :
122 : /// Mutator that returns packets currently in ghost zones in a
123 : /// DirectionMap<Dim,std::vector<Particles::MonteCarlo::Packet>>
124 : /// and remove them from the list of packets of the current
125 : /// element
126 : template <size_t Dim>
127 1 : struct GhostDataMcPackets {
128 0 : using return_tags = tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement>;
129 0 : using argument_tags = tmpl::list<::domain::Tags::Element<Dim>>;
130 :
131 0 : static DirectionMap<Dim, std::vector<Particles::MonteCarlo::Packet>> apply(
132 : const gsl::not_null<std::vector<Particles::MonteCarlo::Packet>*> packets,
133 : const Element<Dim>& element) {
134 : DirectionMap<Dim, std::vector<Particles::MonteCarlo::Packet>> output{};
135 : for (const auto& [direction, neighbors_in_direction] :
136 : element.neighbors()) {
137 : output[direction] = std::vector<Particles::MonteCarlo::Packet>{};
138 : }
139 : // Loop over packets. If a packet is out of the current element,
140 : // remove it from the element and copy it to the appropriate
141 : // neighbor list if it exists.
142 : // Note: At this point, we only handle simple boundaries
143 : // (one neighbor per direction, with the logical coordinates
144 : // being transformed by adding/removing 2). We assume that
145 : // packets that go out of the element but do not match any
146 : // neighbor have reached the domain boundary and can be removed.
147 : size_t n_packets = packets->size();
148 : for (size_t p = 0; p < n_packets; p++) {
149 : Packet& packet = (*packets)[p];
150 : std::optional<Direction<Dim>> max_distance_direction = std::nullopt;
151 : double max_distance = 0.0;
152 : // TO DO: Deal with dimensions higher than the grid dimension
153 : // e.g. for axisymmetric evolution
154 : for (size_t d = 0; d < Dim; d++) {
155 : if (packet.coordinates.get(d) < -1.0) {
156 : const double distance = -1.0 - packet.coordinates.get(d);
157 : if (max_distance < distance) {
158 : max_distance = distance;
159 : max_distance_direction = Direction<Dim>(d, Side::Lower);
160 : }
161 : }
162 : if (packet.coordinates.get(d) > 1.0) {
163 : const double distance = packet.coordinates.get(d) - 1.0;
164 : if (max_distance < distance) {
165 : max_distance = distance;
166 : max_distance_direction = Direction<Dim>(d, Side::Upper);
167 : }
168 : }
169 : }
170 : // If a packet should be moved along an edge / corner, we move
171 : // it into the direct neighbor it is effectively closest to
172 : // (in topological coordinates). This may need improvements.
173 : if (max_distance_direction != std::nullopt) {
174 : const size_t d = max_distance_direction.value().dimension();
175 : const Side side = max_distance_direction.value().side();
176 : packet.coordinates.get(d) += (side == Side::Lower) ? (2.0) : (-2.0);
177 : // Corner/edge treatment; move packet in a live point for now
178 : // Definitely needs improvement...
179 : for (size_t dd = 0; dd < Dim; dd++) {
180 : if (dd != d && packet.coordinates.get(dd) < -1.0) {
181 : packet.coordinates.get(dd) = -2.0 - packet.coordinates.get(dd);
182 : }
183 : if (dd != d && packet.coordinates.get(dd) > 1.0) {
184 : packet.coordinates.get(dd) = 2.0 - packet.coordinates.get(dd);
185 : }
186 : }
187 : if (output.contains(max_distance_direction.value())) {
188 : output[max_distance_direction.value()].push_back(packet);
189 : }
190 : std::swap((*packets)[p], (*packets)[n_packets - 1]);
191 : packets->pop_back();
192 : p--;
193 : n_packets--;
194 : }
195 : }
196 : return output;
197 : }
198 : };
199 :
200 : // Take the data for coupling MC to hydro gathered from neighboring
201 : // ghost zones (stored after communication in GhostZoneCouplingDataTag)
202 : // and add them to the coupling terms in live zones. After calling this
203 : // action, the coupling terms CouplingTildeTau, CouplingTildeRhoYe,
204 : // CouplingTildeS should contain the changes to the energy, RhoYe, and
205 : // momentum of the fluid due to MC packets evolved locally and those
206 : // that evolved within ghost zones on neighboring elements.
207 : // THIS FUNCTION IS CURRENTLY ONLY IMPLENTED IN 3D!
208 : // As other parts of the MC communication, it is also incompatible with
209 : // mesh refinement.
210 : template <size_t Dim>
211 0 : struct CombineCouplingDataPostStep {
212 0 : using return_tags =
213 : tmpl::list<Particles::MonteCarlo::Tags::CouplingTildeTau<DataVector>,
214 : Particles::MonteCarlo::Tags::CouplingTildeRhoYe<DataVector>,
215 : Particles::MonteCarlo::Tags::CouplingTildeS<DataVector, Dim>>;
216 0 : using argument_tags =
217 : tmpl::list<Particles::MonteCarlo::Tags::GhostZoneCouplingDataTag<Dim>,
218 : ::domain::Tags::Element<Dim>,
219 : evolution::dg::subcell::Tags::Mesh<Dim>>;
220 0 : static void apply(
221 : gsl::not_null<Scalar<DataVector>*> coupling_tilde_tau,
222 : gsl::not_null<Scalar<DataVector>*> coupling_tilde_rho_ye,
223 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
224 : coupling_tilde_s,
225 : const Particles::MonteCarlo::GhostZoneCouplingData<Dim>& coupling_data,
226 : const Element<Dim>& element, const Mesh<Dim>& subcell_mesh) {
227 : // Check that we do not use this function for 1D/2D runs
228 : ASSERT(Dim == 3, "CombineCouplingDataPostStep only coded in 3D so far");
229 :
230 : const size_t num_ghost_zones = 1;
231 : const Index<Dim> local_extents = subcell_mesh.extents();
232 : Index<Dim> ghost_extents = local_extents;
233 : for (size_t d = 0; d < 3; d++) {
234 : ghost_extents[d] += 2 * num_ghost_zones;
235 : }
236 : // Loop over neighboring elements
237 : for (const auto& [direction, neighbors_in_direction] :
238 : element.neighbors()) {
239 : for (const auto& neighbor : neighbors_in_direction) {
240 : DirectionalId<Dim> directional_element_id{direction, neighbor};
241 : if (coupling_data.coupling_tilde_tau.at(directional_element_id) ==
242 : std::nullopt) {
243 : continue;
244 : }
245 : const size_t dimension = directional_element_id.direction().dimension();
246 : const Side side = directional_element_id.direction().side();
247 : // Extents of coupling data: mesh with ghost points, sliced
248 : // in direction 'dimension'
249 : Index<Dim> ghost_zone_extents = ghost_extents;
250 : ghost_zone_extents[dimension] = num_ghost_zones;
251 : // These loops only work in 3D; needs fixing for other dimensions
252 : for (size_t i = 0; i < ghost_zone_extents[0]; ++i) {
253 : for (size_t j = 0; j < ghost_zone_extents[1]; ++j) {
254 : for (size_t k = 0; k < ghost_zone_extents[2]; ++k) {
255 : Index<Dim> index_3d{i, j, k};
256 : // Index of point in ghost zone data
257 : const size_t ghost_index =
258 : collapsed_index(index_3d, ghost_zone_extents);
259 : // Index of corresponding live point within the full
260 : // mesh, with ghost zones.
261 : // For two ghost zones e.g., the 1d mesh looks like
262 : // x x o o o o o o x x
263 : // with neighbors
264 : // o o o o o x x x x o o o o o
265 : // with x = GZ, o = live points. On the lower face,
266 : // we thus add data from the first ghost point to
267 : // index num_ghost_zones. On the upper face, the first
268 : // ghost point goes to index local_extents (the size of
269 : // the mesh without ghost zone).
270 : index_3d[dimension] =
271 : (side == Side::Lower)
272 : ? index_3d[dimension] + num_ghost_zones
273 : : local_extents[dimension] + index_3d[dimension];
274 : const size_t extended_index =
275 : collapsed_index(index_3d, ghost_extents);
276 : // Add gathered data to existing coupling terms
277 : coupling_tilde_tau->get()[extended_index] +=
278 : coupling_data.coupling_tilde_tau.at(directional_element_id)
279 : .value()[ghost_index];
280 : coupling_tilde_rho_ye->get()[extended_index] +=
281 : coupling_data.coupling_tilde_rho_ye.at(directional_element_id)
282 : .value()[ghost_index];
283 : for (size_t d = 0; d < Dim; d++) {
284 : coupling_tilde_s->get(d)[extended_index] +=
285 : coupling_data.coupling_tilde_s.at(directional_element_id)
286 : .value()
287 : .get(d)[ghost_index];
288 : }
289 : }
290 : }
291 : }
292 : }
293 : }
294 : }
295 : };
296 :
297 : /// Action responsible for the Send operation of ghost
298 : /// zone communication for Monte-Carlo transport.
299 : /// If CommStep == PreStep, this sends the fluid and
300 : /// metric variables needed for evolution of MC packets.
301 : /// If CommStep == PostStep, this sends packets that
302 : /// have moved from one element to another as well
303 : /// as coupling information from the ghost zone to the
304 : /// live points.
305 : template <size_t Dim, bool LocalTimeStepping,
306 : Particles::MonteCarlo::CommunicationStep CommStep>
307 1 : struct SendDataForMcCommunication {
308 0 : using inbox_tags =
309 : tmpl::list<Particles::MonteCarlo::McGhostZoneDataInboxTag<Dim, CommStep>>;
310 :
311 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
312 : typename ActionList, typename ParallelComponent,
313 : typename Metavariables>
314 0 : static Parallel::iterable_action_return_t apply(
315 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
316 : Parallel::GlobalCache<Metavariables>& cache,
317 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
318 : const ParallelComponent* const /*meta*/) {
319 : static_assert(not LocalTimeStepping,
320 : "Monte Carlo is not tested with local time stepping.");
321 :
322 : ASSERT(db::get<evolution::dg::subcell::Tags::ActiveGrid>(box) ==
323 : evolution::dg::subcell::ActiveGrid::Subcell,
324 : "The SendDataForReconstructionPreStep action in MC "
325 : "can only be called when Subcell is the active scheme.");
326 : db::mutate<Particles::MonteCarlo::Tags::McGhostZoneDataTag<Dim>>(
327 : [](const auto ghost_data_ptr) {
328 : // Clear the previous neighbor data and add current local data
329 : ghost_data_ptr->clear();
330 : },
331 : make_not_null(&box));
332 :
333 : auto& receiver_proxy =
334 : Parallel::get_parallel_component<ParallelComponent>(cache);
335 : const size_t ghost_zone_size = 1;
336 :
337 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
338 : const Mesh<Dim>& subcell_mesh =
339 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(box);
340 : const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box);
341 : Index<Dim> extents_with_ghost_zone = subcell_mesh.extents();
342 : for (size_t d = 0; d < Dim; d++) {
343 : extents_with_ghost_zone[d] += 2 * ghost_zone_size;
344 : }
345 :
346 : // Fill volume data that should be fed to neighbor ghost zones
347 : // PreStep, we send fluid data from the live points to the ghost points.
348 : // PostStep, we send coupling data from the ghost points to the live
349 : // points (note the different DataVector sizes; the PostStep data has
350 : // different dimensions because it is taken from DataVectors including
351 : // ghost points.
352 : DataVector volume_data_to_slice =
353 : CommStep == Particles::MonteCarlo::CommunicationStep::PreStep
354 : ? db::mutate_apply(GhostDataMutatorPreStep{}, make_not_null(&box))
355 : : db::mutate_apply(GhostDataMutatorPostStep<Dim>{},
356 : make_not_null(&box));
357 : const DirectionMap<Dim, DataVector> all_sliced_data =
358 : CommStep == Particles::MonteCarlo::CommunicationStep::PreStep
359 : ? evolution::dg::subcell::slice_data(
360 : volume_data_to_slice, subcell_mesh.extents(), ghost_zone_size,
361 : element.internal_boundaries(), 0,
362 : db::get<evolution::dg::subcell::Tags::
363 : InterpolatorsFromFdToNeighborFd<Dim>>(box))
364 : : evolution::dg::subcell::slice_data(
365 : volume_data_to_slice, extents_with_ghost_zone,
366 : ghost_zone_size, element.internal_boundaries(), 0,
367 : db::get<evolution::dg::subcell::Tags::
368 : InterpolatorsFromFdToNeighborFd<Dim>>(box));
369 : const DirectionMap<Dim, std::vector<Particles::MonteCarlo::Packet>>
370 : all_packets_ghost_zone =
371 : CommStep == Particles::MonteCarlo::CommunicationStep::PostStep
372 : ? db::mutate_apply(GhostDataMcPackets<Dim>{},
373 : make_not_null(&box))
374 : : DirectionMap<Dim,
375 : std::vector<Particles::MonteCarlo::Packet>>{};
376 :
377 : for (const auto& [direction, neighbors_in_direction] :
378 : element.neighbors()) {
379 : for (const ElementId<Dim>& neighbor : neighbors_in_direction) {
380 : const auto& orientation = neighbors_in_direction.orientation(neighbor);
381 : const auto direction_from_neighbor = orientation(direction.opposite());
382 : std::optional<std::vector<Particles::MonteCarlo::Packet>>
383 : packets_to_send = std::nullopt;
384 : DataVector subcell_data_to_send{};
385 :
386 : // PreStep and PostStep both send slice data to neighboring domains
387 : const auto& sliced_data_in_direction = all_sliced_data.at(direction);
388 : subcell_data_to_send = DataVector{sliced_data_in_direction.size()};
389 : std::copy(sliced_data_in_direction.begin(),
390 : sliced_data_in_direction.end(), subcell_data_to_send.begin());
391 : if (CommStep == Particles::MonteCarlo::CommunicationStep::PostStep) {
392 : packets_to_send = all_packets_ghost_zone.at(direction);
393 : if (packets_to_send.value().empty()) {
394 : packets_to_send = std::nullopt;
395 : }
396 : }
397 :
398 : McGhostZoneData<Dim> data{subcell_data_to_send, packets_to_send};
399 :
400 : Parallel::receive_data<
401 : Particles::MonteCarlo::McGhostZoneDataInboxTag<Dim, CommStep>>(
402 : receiver_proxy[neighbor], time_step_id,
403 : std::pair{DirectionalId<Dim>{direction_from_neighbor, element.id()},
404 : std::move(data)});
405 : }
406 : }
407 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
408 : }
409 : };
410 :
411 : /// Action responsible for the Receive operation of ghost
412 : /// zone communication for Monte-Carlo transport.
413 : /// If CommStep == PreStep, this gets the fluid and
414 : /// metric variables needed for evolution of MC packets.
415 : /// If CommStep == PostStep, this gets packets that
416 : /// have moved into the current element.
417 : /// The current code does not yet deal with data required
418 : /// to calculate the backreaction on the fluid.
419 : template <size_t Dim, CommunicationStep CommStep>
420 1 : struct ReceiveDataForMcCommunication {
421 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
422 : typename ActionList, typename ParallelComponent,
423 : typename Metavariables>
424 0 : static Parallel::iterable_action_return_t apply(
425 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
426 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
427 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
428 : const ParallelComponent* const /*meta*/) {
429 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
430 : const auto number_of_expected_messages = element.neighbors().size();
431 : if (UNLIKELY(number_of_expected_messages == 0)) {
432 : // We have no neighbors, so just continue without doing any work
433 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
434 : }
435 :
436 : const auto& current_time_step_id = db::get<::Tags::TimeStepId>(box);
437 : std::map<TimeStepId, DirectionalIdMap<
438 : Dim, Particles::MonteCarlo::McGhostZoneData<Dim>>>&
439 : inbox = tuples::get<Particles::MonteCarlo::McGhostZoneDataInboxTag<
440 : Metavariables::volume_dim, CommStep>>(inboxes);
441 : const auto& received = inbox.find(current_time_step_id);
442 : // Check we have at least some data from correct time, and then check that
443 : // we have received all data
444 : if (received == inbox.end() or
445 : received->second.size() != number_of_expected_messages) {
446 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
447 : }
448 :
449 : // Now that we have received all the data, copy it over as needed.
450 : DirectionalIdMap<Dim, Particles::MonteCarlo::McGhostZoneData<Dim>>
451 : received_data = std::move(inbox[current_time_step_id]);
452 : inbox.erase(current_time_step_id);
453 :
454 : // Copy the received PreStep data in the MortarData container
455 : if (CommStep == Particles::MonteCarlo::CommunicationStep::PreStep) {
456 : db::mutate<Particles::MonteCarlo::Tags::MortarDataTag<Dim>>(
457 : [&element,
458 : &received_data](const gsl::not_null<MortarData<Dim>*> mortar_data) {
459 : for (const auto& [direction, neighbors_in_direction] :
460 : element.neighbors()) {
461 : for (const auto& neighbor : neighbors_in_direction) {
462 : DirectionalId<Dim> directional_element_id{direction, neighbor};
463 : const DataVector& received_data_direction =
464 : received_data[directional_element_id]
465 : .ghost_zone_hydro_variables;
466 : if (mortar_data->rest_mass_density[directional_element_id] ==
467 : std::nullopt) {
468 : continue;
469 : }
470 : const size_t mortar_data_size =
471 : mortar_data->rest_mass_density[directional_element_id]
472 : .value()
473 : .size();
474 : ASSERT(received_data_direction.size() == mortar_data_size * 4,
475 : "Inconsistent sizes between inbox and mortar data");
476 : std::copy(received_data_direction.data(),
477 : std::next(received_data_direction.data(),
478 : static_cast<int>(mortar_data_size)),
479 : mortar_data->rest_mass_density[directional_element_id]
480 : .value()
481 : .data());
482 : std::copy(std::next(received_data_direction.data(),
483 : static_cast<int>(mortar_data_size)),
484 : std::next(received_data_direction.data(),
485 : static_cast<int>(mortar_data_size * 2)),
486 : mortar_data->electron_fraction[directional_element_id]
487 : .value()
488 : .data());
489 : std::copy(std::next(received_data_direction.data(),
490 : static_cast<int>(mortar_data_size * 2)),
491 : std::next(received_data_direction.data(),
492 : static_cast<int>(mortar_data_size * 3)),
493 : mortar_data->temperature[directional_element_id]
494 : .value()
495 : .data());
496 : std::copy(std::next(received_data_direction.data(),
497 : static_cast<int>(mortar_data_size * 3)),
498 : std::next(received_data_direction.data(),
499 : static_cast<int>(mortar_data_size * 4)),
500 : mortar_data
501 : ->cell_light_crossing_time[directional_element_id]
502 : .value()
503 : .data());
504 : }
505 : }
506 : },
507 : make_not_null(&box));
508 : } else {
509 : const Mesh<Dim>& subcell_mesh =
510 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(box);
511 : db::mutate<Particles::MonteCarlo::Tags::GhostZoneCouplingDataTag<Dim>,
512 : Particles::MonteCarlo::Tags::PacketsOnElement>(
513 : [&element, &received_data, &subcell_mesh](
514 : const gsl::not_null<GhostZoneCouplingData<Dim>*> coupling_data,
515 : const gsl::not_null<std::vector<Particles::MonteCarlo::Packet>*>
516 : packet_list) {
517 : // Two parts here: first deal with coupling data brought back
518 : // from ghost zones to live points; then handle packet
519 : // communication.
520 : const Index<Dim>& extents = subcell_mesh.extents();
521 : for (const auto& [direction, neighbors_in_direction] :
522 : element.neighbors()) {
523 : for (const auto& neighbor : neighbors_in_direction) {
524 : DirectionalId<Dim> directional_element_id{direction, neighbor};
525 : // Move boundary data to the coupling_data structure for later
526 : // use
527 : const DataVector& received_data_direction =
528 : received_data[directional_element_id]
529 : .ghost_zone_hydro_variables;
530 : if (coupling_data->coupling_tilde_tau[directional_element_id] !=
531 : std::nullopt) {
532 : const size_t coupling_data_size =
533 : coupling_data->coupling_tilde_tau[directional_element_id]
534 : .value()
535 : .size();
536 : ASSERT(received_data_direction.size() ==
537 : coupling_data_size * (2 + Dim),
538 : "Inconsistent sizes between inbox and coupling data");
539 : std::copy(
540 : received_data_direction.data(),
541 : std::next(received_data_direction.data(),
542 : static_cast<int>(coupling_data_size)),
543 : coupling_data->coupling_tilde_tau[directional_element_id]
544 : .value()
545 : .data());
546 : std::copy(std::next(received_data_direction.data(),
547 : static_cast<int>(coupling_data_size)),
548 : std::next(received_data_direction.data(),
549 : static_cast<int>(coupling_data_size * 2)),
550 : coupling_data
551 : ->coupling_tilde_rho_ye[directional_element_id]
552 : .value()
553 : .data());
554 : for (size_t d = 0; d < Dim; d++) {
555 : std::copy(
556 : std::next(
557 : received_data_direction.data(),
558 : static_cast<int>(coupling_data_size * (2 + d))),
559 : std::next(
560 : received_data_direction.data(),
561 : static_cast<int>(coupling_data_size * (3 + d))),
562 : coupling_data->coupling_tilde_s[directional_element_id]
563 : .value()
564 : .get(d)
565 : .data());
566 : }
567 : }
568 : // Get packet data
569 : std::optional<std::vector<Particles::MonteCarlo::Packet>>&
570 : received_data_packets =
571 : received_data[directional_element_id]
572 : .packets_entering_this_element;
573 : if (received_data_packets == std::nullopt) {
574 : continue;
575 : } else {
576 : const size_t n_packets = received_data_packets.value().size();
577 : for (size_t p = 0; p < n_packets; p++) {
578 : // Find closest grid point to received packet at current
579 : // time, using extents for live points only.
580 : Packet& packet = received_data_packets.value()[p];
581 : Index<Dim> closest_point_index_3d;
582 : for (size_t d = 0; d < Dim; d++) {
583 : closest_point_index_3d[d] =
584 : std::floor((packet.coordinates[d] + 1.0) *
585 : static_cast<double>(extents[d]) / 2.0);
586 : if (UNLIKELY(closest_point_index_3d[d] == extents[d])) {
587 : closest_point_index_3d[d]--;
588 : }
589 : }
590 : packet.index_of_closest_grid_point =
591 : collapsed_index(closest_point_index_3d, extents);
592 : // Now add packets to list in current element
593 : packet_list->push_back(received_data_packets.value()[p]);
594 : }
595 : }
596 : }
597 : }
598 : },
599 : make_not_null(&box));
600 : // Combine ghost zone data with live points data. Currently only done in
601 : // 3D
602 : if constexpr (Dim == 3) {
603 : db::mutate_apply(CombineCouplingDataPostStep<Dim>{},
604 : make_not_null(&box));
605 : }
606 : }
607 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
608 : }
609 : };
610 :
611 : } // namespace Particles::MonteCarlo::Actions
|
