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 : evolution::dg::borders_only_dg_block(element, only_dg_block_ids)) {
101 : return;
102 : }
103 : const bool enable_extension_directions =
104 : subcell_options.enable_extension_directions();
105 : if (enable_extension_directions) {
106 : *extension_direction_ptr = {};
107 : } // Initialize the extension directions to empty.
108 :
109 : const size_t number_of_ghost_zones = reconstructor.ghost_zone_size();
110 : const size_t my_block_id = element.id().block_id();
111 : for (const auto& direction_neighbors_in_direction : element.neighbors()) {
112 : const auto& direction = direction_neighbors_in_direction.first;
113 : const auto& neighbors_in_direction =
114 : direction_neighbors_in_direction.second;
115 : for (const ElementId<Dim>& neighbor_id : neighbors_in_direction) {
116 : const auto& orientation =
117 : neighbors_in_direction.orientation(neighbor_id);
118 : const auto direction_from_neighbor = orientation(direction.opposite());
119 : const size_t neighbor_block_id = neighbor_id.block_id();
120 : if (neighbor_block_id == my_block_id) {
121 : continue;
122 : }
123 : const auto& neighbor_block = domain.blocks()[neighbor_block_id];
124 : // InterpolatorsFromFdToNeighborFd &
125 : // InterpolatorsFromDgToNeighborFd
126 : // 1. Compute the grid coordinates of my neighbor's ghost zones.
127 : // 2. Compute the element logical coordinates of my neighbor's
128 : // ghost zones.
129 : // 3. Create interpolators
130 :
131 : if (not is_isotropic(neighbor_fd_mesh) and
132 : neighbor_fd_mesh.basis(Dim - 1) != Spectral::Basis::Cartoon) {
133 : ERROR("We assume an isotropic mesh but got "
134 : << neighbor_fd_mesh << " ElementID is " << element.id());
135 : }
136 : // Extra checks for cartoon meshes not checked above
137 : fd::verify_subcell_mesh(neighbor_fd_mesh, true);
138 :
139 : const auto get_logical_coords = [&element, &neighbor_id, &direction](
140 : const auto& map,
141 : const auto& grid_coords) {
142 : tnsr::I<DataVector, Dim, Frame::ElementLogical> logical_coords{
143 : get<0>(grid_coords).size()};
144 : for (size_t i = 0; i < get<0>(grid_coords).size(); ++i) {
145 : try {
146 : tnsr::I<double, Dim, Frame::ElementLogical> logical_coord =
147 : map.inverse(extract_point(grid_coords, i));
148 : for (size_t d = 0; d < Dim; ++d) {
149 : logical_coords.get(d)[i] = logical_coord.get(d);
150 : }
151 : } catch (const std::bad_optional_access& e) {
152 : ERROR(
153 : "Failed to get logical coordinates for neighbor's "
154 : "ghost zone grid coordinates. This could be because the "
155 : "ghost zones are not in the nearest neighbor but instead in "
156 : "the next-to-nearest neighbor. The code assumes all ghost "
157 : "zones, even on curved meshes, are in the nearest neighbors. "
158 : "The current element is "
159 : << element.id() << " and the neighbor id is " << neighbor_id
160 : << " in direction " << direction
161 : << " The neighbor grid coordinates are \n"
162 : << extract_point(grid_coords, i) << "\n");
163 : }
164 : }
165 : return logical_coords;
166 : };
167 :
168 : tnsr::I<DataVector, Dim, Frame::Grid> neighbor_grid_ghost_zone_coords{};
169 : // Get the neighbor's ghost zone coordinates in the grid
170 : // frame.
171 : if (const tnsr::I<DataVector, Dim, Frame::ElementLogical>
172 : neighbor_logical_ghost_zone_coords =
173 : evolution::dg::subcell::fd::ghost_zone_logical_coordinates(
174 : neighbor_fd_mesh, number_of_ghost_zones,
175 : direction_from_neighbor);
176 : neighbor_block.is_time_dependent()) {
177 : const ElementMap neighbor_element_map(
178 : neighbor_id,
179 : neighbor_block.moving_mesh_logical_to_grid_map().get_clone());
180 : neighbor_grid_ghost_zone_coords =
181 : neighbor_element_map(neighbor_logical_ghost_zone_coords);
182 : } else {
183 : const ElementMap neighbor_element_map(
184 : neighbor_id, neighbor_block.stationary_map().get_clone());
185 : const tnsr::I<DataVector, Dim, Frame::Inertial>
186 : neighbor_inertial_ghost_zone_coords =
187 : neighbor_element_map(neighbor_logical_ghost_zone_coords);
188 : for (size_t i = 0; i < Dim; ++i) {
189 : neighbor_grid_ghost_zone_coords[i] =
190 : neighbor_inertial_ghost_zone_coords[i];
191 : }
192 : }
193 : // Map the ghost zone grid coordinates back to our logical
194 : // coordinates.
195 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>
196 : neighbor_logical_ghost_zone_coords = get_logical_coords(
197 : element_map, neighbor_grid_ghost_zone_coords);
198 :
199 : // We need to check if the neighbor's ghost zone coordinates
200 : // are in the same element as the current element. If not, we
201 : // need to extend the mesh in the direction of the ghost zone
202 : // coordinates.
203 : // Note, we only need to do this check if we enable extending
204 : // the mesh to avoid extrapolation. If not, we simply
205 : // allow the extrapolation to happen.
206 :
207 : // 'needs_extension' is set to true if the neighbor's ghost zone
208 : // coordinates are outside the current element's mesh and if
209 : // we enable extending the mesh.
210 : bool needs_extension = false;
211 : std::optional<Direction<Dim>> direction_to_extend;
212 : if (enable_extension_directions) {
213 : for (size_t d = 0; d < Dim; ++d) {
214 : // small epsilon of 1e-10 used to ensure we are not unncessarily
215 : // flagging as problematic
216 : const double ext = 1. - (1. / my_fd_mesh.extents(d)) + 1.e-10;
217 : const auto& coords = neighbor_logical_ghost_zone_coords.get(d);
218 :
219 : for (size_t i = 0; i < coords.size(); ++i) {
220 : if (std::abs(coords[i]) > ext) {
221 : needs_extension = true;
222 : Direction<Dim> new_direction = Direction<Dim>{
223 : d, coords[i] > 0 ? Side::Upper : Side::Lower};
224 :
225 : if (!direction_to_extend.has_value()) {
226 : direction_to_extend = new_direction;
227 : } else if (direction_to_extend.value() != new_direction) {
228 : ERROR("Multiple directions to extend: existing = "
229 : << direction_to_extend.value()
230 : << ", new = " << new_direction);
231 : }
232 : break; // no reason to check remaining coords.
233 : }
234 : }
235 : }
236 : }
237 :
238 : if (needs_extension) {
239 : if (!direction_to_extend.has_value()) {
240 : ERROR(
241 : "Should have direction to extend if flagged "
242 : "as problematic!");
243 : }
244 : const auto& external_boundaries = element.external_boundaries();
245 : if (external_boundaries.find(direction_to_extend.value()) !=
246 : external_boundaries.end()) {
247 : ERROR(
248 : "Direction to extend is toward the "
249 : "faces of the Element that are external boundaries.");
250 : }
251 :
252 : auto new_basis = make_array<Dim>(my_fd_mesh.basis(0));
253 : auto new_extents = make_array<Dim>(my_fd_mesh.extents(0));
254 : auto new_quads = make_array<Dim>(my_fd_mesh.quadrature(0));
255 :
256 : const size_t problematic_dim =
257 : direction_to_extend.value().dimension();
258 :
259 : // note we are extending our current volume by including its own ghost
260 : // points in the problematic direction (direction to extend)
261 : // which means the logical coordinates of the ghost (to be sent)
262 : // must be transformed to accommodate the extended mesh in
263 : // direction to extend.
264 :
265 : const double rescale_factor =
266 : static_cast<double>(my_fd_mesh.extents(problematic_dim)) /
267 : (my_fd_mesh.extents(problematic_dim) + number_of_ghost_zones);
268 : double translation =
269 : static_cast<double>(number_of_ghost_zones) /
270 : (my_fd_mesh.extents(problematic_dim) + number_of_ghost_zones);
271 : // translation above is based on extending to Upper Side.
272 : if (direction_to_extend.value().side() == Side::Lower) {
273 : translation *= -1.;
274 : }
275 : auto new_neighbor_logical_ghost_zone_coords =
276 : neighbor_logical_ghost_zone_coords;
277 : for (size_t i = 0;
278 : i < new_neighbor_logical_ghost_zone_coords[0].size(); ++i) {
279 : new_neighbor_logical_ghost_zone_coords.get(problematic_dim)[i] *=
280 : rescale_factor;
281 : new_neighbor_logical_ghost_zone_coords.get(problematic_dim)[i] -=
282 : translation;
283 : }
284 :
285 : for (size_t d = 0; d < Dim; ++d) {
286 : gsl::at(new_basis, d) = my_fd_mesh.basis(d);
287 : gsl::at(new_quads, d) = my_fd_mesh.quadrature(d);
288 : if (d == problematic_dim) {
289 : gsl::at(new_extents, d) =
290 : my_fd_mesh.extents(d) + number_of_ghost_zones;
291 : } else {
292 : gsl::at(new_extents, d) = my_fd_mesh.extents(d);
293 : }
294 : }
295 : const Mesh<Dim> new_mesh{new_extents, new_basis, new_quads};
296 : (*interpolators_fd_to_neighbor_fd_ptr)[DirectionalId<Dim>{
297 : direction, neighbor_id}] = intrp::Irregular<Dim>{
298 : new_mesh, new_neighbor_logical_ghost_zone_coords,
299 : subcell_options.get_fd_to_fd_interp_order()};
300 : (*extension_direction_ptr)[direction] =
301 : interpolators_detail::ExtensionDirection<Dim>{
302 : direction_to_extend.value()};
303 : } else {
304 : (*interpolators_fd_to_neighbor_fd_ptr)[DirectionalId<Dim>{
305 : direction, neighbor_id}] = intrp::Irregular<Dim>{
306 : my_fd_mesh, neighbor_logical_ghost_zone_coords,
307 : subcell_options.get_fd_to_fd_interp_order()};
308 : }
309 : // Set up interpolators for our local element to our neighbor's
310 : // ghost zones.
311 : (*interpolators_dg_to_neighbor_fd_ptr)[DirectionalId<Dim>{
312 : direction, neighbor_id}] = intrp::Irregular<Dim>{
313 : my_dg_mesh, neighbor_logical_ghost_zone_coords};
314 :
315 : // InterpolatorsFromNeighborDgToFd: the interpolation from our
316 : // neighbor's DG grid to our FD ghost zones.
317 : //
318 : // 1. Compute the grid coordinates of my ghost zones.
319 : // 2. Compute neighbor's element logical coordinates of my ghost
320 : // zones
321 : // 3. Create interpolator for InterpolatorsFromNeighborDgToFd
322 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>
323 : my_logical_coords = logical_coordinates(my_fd_mesh);
324 : tnsr::I<DataVector, Dim, Frame::ElementLogical>
325 : my_logical_ghost_zone_coords =
326 : evolution::dg::subcell::slice_tensor_for_subcell(
327 : my_logical_coords, neighbor_fd_mesh.extents(),
328 : number_of_ghost_zones, direction,
329 : // We want to _set_ the interpolators, so just do a simple
330 : // slice.
331 : {});
332 : const double delta_xi =
333 : get<0>(my_logical_coords)[1] - get<0>(my_logical_coords)[0];
334 : // The sign accounts for whether we are shift along the
335 : // positive or negative axis.
336 : const double coordinate_shift =
337 : direction.sign() * delta_xi * number_of_ghost_zones;
338 : my_logical_ghost_zone_coords.get(direction.dimension()) +=
339 : coordinate_shift;
340 : const tnsr::I<DataVector, Dim, Frame::Grid> my_grid_ghost_zone_coords =
341 : element_map(my_logical_ghost_zone_coords);
342 : if (neighbor_block.is_time_dependent()) {
343 : const ElementMap neighbor_element_map(
344 : neighbor_id,
345 : neighbor_block.moving_mesh_logical_to_grid_map().get_clone());
346 : (*interpolators_neighbor_dg_to_fd_ptr)[DirectionalId<Dim>{
347 : direction, neighbor_id}] = intrp::Irregular<Dim>{
348 : neighbor_dg_mesh, get_logical_coords(neighbor_element_map,
349 : my_grid_ghost_zone_coords)};
350 : } else {
351 : const ElementMap neighbor_element_map(
352 : neighbor_id, neighbor_block.stationary_map().get_clone());
353 : const tnsr::I<DataVector, Dim, Frame::Inertial>
354 : view_my_grid_ghost_zone_coords{};
355 : for (size_t i = 0; i < Dim; ++i) {
356 : make_const_view(make_not_null(&view_my_grid_ghost_zone_coords[i]),
357 : my_grid_ghost_zone_coords[i], 0,
358 : my_grid_ghost_zone_coords[i].size());
359 : }
360 : (*interpolators_neighbor_dg_to_fd_ptr)[DirectionalId<Dim>{
361 : direction, neighbor_id}] = intrp::Irregular<Dim>{
362 : neighbor_dg_mesh,
363 : get_logical_coords(neighbor_element_map,
364 : view_my_grid_ghost_zone_coords)};
365 : }
366 : }
367 : }
368 : }
369 : };
370 : } // namespace evolution::dg::subcell
|
