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