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