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