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