Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <optional>
8 : #include <utility>
9 :
10 : #include "DataStructures/ExtractPoint.hpp"
11 : #include "DataStructures/Tensor/Tensor.hpp"
12 : #include "Domain/Creators/Tags/Domain.hpp"
13 : #include "Domain/Domain.hpp"
14 : #include "Domain/ElementMap.hpp"
15 : #include "Domain/Structure/BlockId.hpp"
16 : #include "Domain/Structure/Direction.hpp"
17 : #include "Domain/Structure/DirectionMap.hpp"
18 : #include "Domain/Structure/DirectionalId.hpp"
19 : #include "Domain/Structure/DirectionalIdMap.hpp"
20 : #include "Domain/Structure/Element.hpp"
21 : #include "Domain/Structure/ElementId.hpp"
22 : #include "Domain/Tags.hpp"
23 : #include "Evolution/DgSubcell/GhostZoneLogicalCoordinates.hpp"
24 : #include "Evolution/DgSubcell/Mesh.hpp"
25 : #include "Evolution/DgSubcell/SliceTensor.hpp"
26 : #include "Evolution/DgSubcell/SubcellOptions.hpp"
27 : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
28 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
29 : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
30 : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
31 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
32 : #include "Utilities/ErrorHandling/Error.hpp"
33 : #include "Utilities/Gsl.hpp"
34 : #include "Utilities/TMPL.hpp"
35 :
36 : namespace evolution::dg::subcell {
37 : /*!
38 : * \brief Sets the `intrp::IrregularInterpolant`s for interpolating to ghost
39 : * zone data at block boundaries.
40 : *
41 : * The DG to FD interpolants are at full order of the DG grid. The FD to FD
42 : * interpolant is piecewise linear with no support for neighboring
43 : * elements. We will want to use high-order slope-limited FD interpolation in
44 : * the future, but that requires neighbor communication. A slightly simpler
45 : * approach would be to use high-order Lagrange interpolation, which still
46 : * requires neighbor communication but does not require any additional changes
47 : * to the reconstruction routines to work on non-uniform grids. This is what
48 : * the Multipatch-MHD code does, relying on the slope limiting from the ghost
49 : * zones to remove oscillations. I (Nils Deppe) am not sure I love that, but
50 : * it's worth a try since it should be pretty easy to do.
51 : *
52 : * \warning Currently assumes that neighboring DG/FD elements are on the same
53 : * refinement level and have the same DG mesh and subcell mesh.
54 : */
55 : template <size_t Dim, typename ReconstructorTag>
56 1 : struct SetInterpolators {
57 0 : using return_tags = tmpl::list<
58 : evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<Dim>,
59 : evolution::dg::subcell::Tags::InterpolatorsFromDgToNeighborFd<Dim>,
60 : evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>,
61 : evolution::dg::subcell::Tags::ExtensionDirections<Dim>>;
62 0 : using argument_tags =
63 : tmpl::list<::domain::Tags::Element<Dim>, ::domain::Tags::Domain<Dim>,
64 : domain::Tags::Mesh<Dim>, domain::Tags::Mesh<Dim>,
65 : evolution::dg::subcell::Tags::Mesh<Dim>,
66 : evolution::dg::subcell::Tags::Mesh<Dim>,
67 : ::domain::Tags::ElementMap<Dim, Frame::Grid>, ReconstructorTag,
68 : evolution::dg::subcell::Tags::SubcellOptions<Dim>>;
69 :
70 : template <typename ReconstructorType>
71 0 : static void apply(
72 : const gsl::not_null<
73 : DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>*>
74 : interpolators_fd_to_neighbor_fd_ptr,
75 : const gsl::not_null<
76 : DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>*>
77 : interpolators_dg_to_neighbor_fd_ptr,
78 : const gsl::not_null<
79 : DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>*>
80 : interpolators_neighbor_dg_to_fd_ptr,
81 : const gsl::not_null<
82 : DirectionMap<Dim, interpolators_detail::ExtensionDirection<Dim>>*>
83 : extension_direction_ptr,
84 : const Element<Dim>& element, const Domain<Dim>& domain,
85 : const Mesh<Dim>& my_dg_mesh,
86 : // Needs to be updated to support non-uniform h/p-refinement
87 : const Mesh<Dim>& neighbor_dg_mesh, const Mesh<Dim>& my_fd_mesh,
88 : // Needs to be updated to support non-uniform h/p-refinement
89 : const Mesh<Dim>& neighbor_fd_mesh,
90 : const ElementMap<Dim, Frame::Grid>& element_map,
91 : const ReconstructorType& reconstructor,
92 : const evolution::dg::subcell::SubcellOptions& subcell_options) {
93 : // Skip for elements that are DG-only: either in a DG-only block,
94 : // bordering a DG-only block, or on a non-subcell-compatible mesh.
95 : if (not fd::dg_mesh_supports_subcell(my_dg_mesh) or
96 : alg::found(subcell_options.only_dg_block_ids(),
97 : element.id().block_id()) or
98 : alg::any_of(
99 : element.neighbors(),
100 : [&subcell_options](const auto& direction_and_neighbors) {
101 : return alg::found(
102 : subcell_options.only_dg_block_ids(),
103 : direction_and_neighbors.second.ids().begin()->block_id());
104 : })) {
105 : return;
106 : }
107 : const bool enable_extension_directions =
108 : subcell_options.enable_extension_directions();
109 : if (enable_extension_directions) {
110 : *extension_direction_ptr = {};
111 : } // Initialize the extension directions to empty.
112 :
113 : const size_t number_of_ghost_zones = reconstructor.ghost_zone_size();
114 : const size_t my_block_id = element.id().block_id();
115 : for (const auto& direction_neighbors_in_direction : element.neighbors()) {
116 : const auto& direction = direction_neighbors_in_direction.first;
117 : const auto& neighbors_in_direction =
118 : direction_neighbors_in_direction.second;
119 : for (const ElementId<Dim>& neighbor_id : neighbors_in_direction) {
120 : const auto& orientation =
121 : neighbors_in_direction.orientation(neighbor_id);
122 : const auto direction_from_neighbor = orientation(direction.opposite());
123 : const size_t neighbor_block_id = neighbor_id.block_id();
124 : if (neighbor_block_id == my_block_id) {
125 : continue;
126 : }
127 : const auto& neighbor_block = domain.blocks()[neighbor_block_id];
128 : // InterpolatorsFromFdToNeighborFd &
129 : // InterpolatorsFromDgToNeighborFd
130 : // 1. Compute the grid coordinates of my neighbor's ghost zones.
131 : // 2. Compute the element logical coordinates of my neighbor's
132 : // ghost zones.
133 : // 3. Create interpolators
134 :
135 : if (not is_isotropic(neighbor_fd_mesh) and
136 : neighbor_fd_mesh.basis(Dim - 1) != Spectral::Basis::Cartoon) {
137 : ERROR("We assume an isotropic mesh but got "
138 : << neighbor_fd_mesh << " ElementID is " << element.id());
139 : }
140 : // Extra checks for cartoon meshes not checked above
141 : fd::verify_subcell_mesh(neighbor_fd_mesh, true);
142 :
143 : const auto get_logical_coords = [&element, &neighbor_id, &direction](
144 : const auto& map,
145 : const auto& grid_coords) {
146 : tnsr::I<DataVector, Dim, Frame::ElementLogical> logical_coords{
147 : get<0>(grid_coords).size()};
148 : for (size_t i = 0; i < get<0>(grid_coords).size(); ++i) {
149 : try {
150 : tnsr::I<double, Dim, Frame::ElementLogical> logical_coord =
151 : map.inverse(extract_point(grid_coords, i));
152 : for (size_t d = 0; d < Dim; ++d) {
153 : logical_coords.get(d)[i] = logical_coord.get(d);
154 : }
155 : } catch (const std::bad_optional_access& e) {
156 : ERROR(
157 : "Failed to get logical coordinates for neighbor's "
158 : "ghost zone grid coordinates. This could be because the "
159 : "ghost zones are not in the nearest neighbor but instead in "
160 : "the next-to-nearest neighbor. The code assumes all ghost "
161 : "zones, even on curved meshes, are in the nearest neighbors. "
162 : "The current element is "
163 : << element.id() << " and the neighbor id is " << neighbor_id
164 : << " in direction " << direction
165 : << " The neighbor grid coordinates are \n"
166 : << extract_point(grid_coords, i) << "\n");
167 : }
168 : }
169 : return logical_coords;
170 : };
171 :
172 : tnsr::I<DataVector, Dim, Frame::Grid> neighbor_grid_ghost_zone_coords{};
173 : // Get the neighbor's ghost zone coordinates in the grid
174 : // frame.
175 : if (const tnsr::I<DataVector, Dim, Frame::ElementLogical>
176 : neighbor_logical_ghost_zone_coords =
177 : evolution::dg::subcell::fd::ghost_zone_logical_coordinates(
178 : neighbor_fd_mesh, number_of_ghost_zones,
179 : direction_from_neighbor);
180 : neighbor_block.is_time_dependent()) {
181 : const ElementMap neighbor_element_map(
182 : neighbor_id,
183 : neighbor_block.moving_mesh_logical_to_grid_map().get_clone());
184 : neighbor_grid_ghost_zone_coords =
185 : neighbor_element_map(neighbor_logical_ghost_zone_coords);
186 : } else {
187 : const ElementMap neighbor_element_map(
188 : neighbor_id, neighbor_block.stationary_map().get_clone());
189 : const tnsr::I<DataVector, Dim, Frame::Inertial>
190 : neighbor_inertial_ghost_zone_coords =
191 : neighbor_element_map(neighbor_logical_ghost_zone_coords);
192 : for (size_t i = 0; i < Dim; ++i) {
193 : neighbor_grid_ghost_zone_coords[i] =
194 : neighbor_inertial_ghost_zone_coords[i];
195 : }
196 : }
197 : // Map the ghost zone grid coordinates back to our logical
198 : // coordinates.
199 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>
200 : neighbor_logical_ghost_zone_coords = get_logical_coords(
201 : element_map, neighbor_grid_ghost_zone_coords);
202 :
203 : // We need to check if the neighbor's ghost zone coordinates
204 : // are in the same element as the current element. If not, we
205 : // need to extend the mesh in the direction of the ghost zone
206 : // coordinates.
207 : // Note, we only need to do this check if we enable extending
208 : // the mesh to avoid extrapolation. If not, we simply
209 : // allow the extrapolation to happen.
210 :
211 : // 'needs_extension' is set to true if the neighbor's ghost zone
212 : // coordinates are outside the current element's mesh and if
213 : // we enable extending the mesh.
214 : bool needs_extension = false;
215 : std::optional<Direction<Dim>> direction_to_extend;
216 : if (enable_extension_directions) {
217 : for (size_t d = 0; d < Dim; ++d) {
218 : // small epsilon of 1e-10 used to ensure we are not unncessarily
219 : // flagging as problematic
220 : const double ext = 1. - (1. / my_fd_mesh.extents(d)) + 1.e-10;
221 : const auto& coords = neighbor_logical_ghost_zone_coords.get(d);
222 :
223 : for (size_t i = 0; i < coords.size(); ++i) {
224 : if (std::abs(coords[i]) > ext) {
225 : needs_extension = true;
226 : Direction<Dim> new_direction = Direction<Dim>{
227 : d, coords[i] > 0 ? Side::Upper : Side::Lower};
228 :
229 : if (!direction_to_extend.has_value()) {
230 : direction_to_extend = new_direction;
231 : } else if (direction_to_extend.value() != new_direction) {
232 : ERROR("Multiple directions to extend: existing = "
233 : << direction_to_extend.value()
234 : << ", new = " << new_direction);
235 : }
236 : break; // no reason to check remaining coords.
237 : }
238 : }
239 : }
240 : }
241 :
242 : if (needs_extension) {
243 : if (!direction_to_extend.has_value()) {
244 : ERROR(
245 : "Should have direction to extend if flagged "
246 : "as problematic!");
247 : }
248 : const auto& external_boundaries = element.external_boundaries();
249 : if (external_boundaries.find(direction_to_extend.value()) !=
250 : external_boundaries.end()) {
251 : ERROR(
252 : "Direction to extend is toward the "
253 : "faces of the Element that are external boundaries.");
254 : }
255 :
256 : auto new_basis = make_array<Dim>(my_fd_mesh.basis(0));
257 : auto new_extents = make_array<Dim>(my_fd_mesh.extents(0));
258 : auto new_quads = make_array<Dim>(my_fd_mesh.quadrature(0));
259 :
260 : const size_t problematic_dim =
261 : direction_to_extend.value().dimension();
262 :
263 : // note we are extending our current volume by including its own ghost
264 : // points in the problematic direction (direction to extend)
265 : // which means the logical coordinates of the ghost (to be sent)
266 : // must be transformed to accommodate the extended mesh in
267 : // direction to extend.
268 :
269 : const double rescale_factor =
270 : static_cast<double>(my_fd_mesh.extents(problematic_dim)) /
271 : (my_fd_mesh.extents(problematic_dim) + number_of_ghost_zones);
272 : double translation =
273 : static_cast<double>(number_of_ghost_zones) /
274 : (my_fd_mesh.extents(problematic_dim) + number_of_ghost_zones);
275 : // translation above is based on extending to Upper Side.
276 : if (direction_to_extend.value().side() == Side::Lower) {
277 : translation *= -1.;
278 : }
279 : auto new_neighbor_logical_ghost_zone_coords =
280 : neighbor_logical_ghost_zone_coords;
281 : for (size_t i = 0;
282 : i < new_neighbor_logical_ghost_zone_coords[0].size(); ++i) {
283 : new_neighbor_logical_ghost_zone_coords.get(problematic_dim)[i] *=
284 : rescale_factor;
285 : new_neighbor_logical_ghost_zone_coords.get(problematic_dim)[i] -=
286 : translation;
287 : }
288 :
289 : for (size_t d = 0; d < Dim; ++d) {
290 : gsl::at(new_basis, d) = my_fd_mesh.basis(d);
291 : gsl::at(new_quads, d) = my_fd_mesh.quadrature(d);
292 : if (d == problematic_dim) {
293 : gsl::at(new_extents, d) =
294 : my_fd_mesh.extents(d) + number_of_ghost_zones;
295 : } else {
296 : gsl::at(new_extents, d) = my_fd_mesh.extents(d);
297 : }
298 : }
299 : const Mesh<Dim> new_mesh{new_extents, new_basis, new_quads};
300 : (*interpolators_fd_to_neighbor_fd_ptr)[DirectionalId<Dim>{
301 : direction, neighbor_id}] = intrp::Irregular<Dim>{
302 : new_mesh, new_neighbor_logical_ghost_zone_coords,
303 : subcell_options.get_fd_to_fd_interp_order()};
304 : (*extension_direction_ptr)[direction] =
305 : interpolators_detail::ExtensionDirection<Dim>{
306 : direction_to_extend.value()};
307 : } else {
308 : (*interpolators_fd_to_neighbor_fd_ptr)[DirectionalId<Dim>{
309 : direction, neighbor_id}] = intrp::Irregular<Dim>{
310 : my_fd_mesh, neighbor_logical_ghost_zone_coords,
311 : subcell_options.get_fd_to_fd_interp_order()};
312 : }
313 : // Set up interpolators for our local element to our neighbor's
314 : // ghost zones.
315 : (*interpolators_dg_to_neighbor_fd_ptr)[DirectionalId<Dim>{
316 : direction, neighbor_id}] = intrp::Irregular<Dim>{
317 : my_dg_mesh, neighbor_logical_ghost_zone_coords};
318 :
319 : // InterpolatorsFromNeighborDgToFd: the interpolation from our
320 : // neighbor's DG grid to our FD ghost zones.
321 : //
322 : // 1. Compute the grid coordinates of my ghost zones.
323 : // 2. Compute neighbor's element logical coordinates of my ghost
324 : // zones
325 : // 3. Create interpolator for InterpolatorsFromNeighborDgToFd
326 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>
327 : my_logical_coords = logical_coordinates(my_fd_mesh);
328 : tnsr::I<DataVector, Dim, Frame::ElementLogical>
329 : my_logical_ghost_zone_coords =
330 : evolution::dg::subcell::slice_tensor_for_subcell(
331 : my_logical_coords, neighbor_fd_mesh.extents(),
332 : number_of_ghost_zones, direction,
333 : // We want to _set_ the interpolators, so just do a simple
334 : // slice.
335 : {});
336 : const double delta_xi =
337 : get<0>(my_logical_coords)[1] - get<0>(my_logical_coords)[0];
338 : // The sign accounts for whether we are shift along the
339 : // positive or negative axis.
340 : const double coordinate_shift =
341 : direction.sign() * delta_xi * number_of_ghost_zones;
342 : my_logical_ghost_zone_coords.get(direction.dimension()) +=
343 : coordinate_shift;
344 : const tnsr::I<DataVector, Dim, Frame::Grid> my_grid_ghost_zone_coords =
345 : element_map(my_logical_ghost_zone_coords);
346 : if (neighbor_block.is_time_dependent()) {
347 : const ElementMap neighbor_element_map(
348 : neighbor_id,
349 : neighbor_block.moving_mesh_logical_to_grid_map().get_clone());
350 : (*interpolators_neighbor_dg_to_fd_ptr)[DirectionalId<Dim>{
351 : direction, neighbor_id}] = intrp::Irregular<Dim>{
352 : neighbor_dg_mesh, get_logical_coords(neighbor_element_map,
353 : my_grid_ghost_zone_coords)};
354 : } else {
355 : const ElementMap neighbor_element_map(
356 : neighbor_id, neighbor_block.stationary_map().get_clone());
357 : const tnsr::I<DataVector, Dim, Frame::Inertial>
358 : view_my_grid_ghost_zone_coords{};
359 : for (size_t i = 0; i < Dim; ++i) {
360 : make_const_view(make_not_null(&view_my_grid_ghost_zone_coords[i]),
361 : my_grid_ghost_zone_coords[i], 0,
362 : my_grid_ghost_zone_coords[i].size());
363 : }
364 : (*interpolators_neighbor_dg_to_fd_ptr)[DirectionalId<Dim>{
365 : direction, neighbor_id}] = intrp::Irregular<Dim>{
366 : neighbor_dg_mesh,
367 : get_logical_coords(neighbor_element_map,
368 : view_my_grid_ghost_zone_coords)};
369 : }
370 : }
371 : }
372 : }
373 : };
374 : } // namespace evolution::dg::subcell
|
