Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <type_traits>
9 : #include <unordered_map>
10 : #include <utility>
11 : #include <vector>
12 :
13 : #include "DataStructures/DataBox/MetavariablesTag.hpp"
14 : #include "DataStructures/SliceVariables.hpp"
15 : #include "DataStructures/TaggedTuple.hpp"
16 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
17 : #include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"
18 : #include "DataStructures/Tensor/Tensor.hpp"
19 : #include "DataStructures/Variables.hpp"
20 : #include "DataStructures/VariablesTag.hpp"
21 : #include "Domain/Creators/Tags/Domain.hpp"
22 : #include "Domain/Creators/Tags/InitialExtents.hpp"
23 : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
24 : #include "Domain/Domain.hpp"
25 : #include "Domain/ElementMap.hpp"
26 : #include "Domain/FaceNormal.hpp"
27 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
28 : #include "Domain/FunctionsOfTime/Tags.hpp"
29 : #include "Domain/InterfaceLogicalCoordinates.hpp"
30 : #include "Domain/Structure/CreateInitialMesh.hpp"
31 : #include "Domain/Structure/Direction.hpp"
32 : #include "Domain/Structure/DirectionMap.hpp"
33 : #include "Domain/Structure/DirectionalId.hpp"
34 : #include "Domain/Structure/Element.hpp"
35 : #include "Domain/Structure/IndexToSliceAt.hpp"
36 : #include "Domain/Structure/OrientationMapHelpers.hpp"
37 : #include "Domain/Tags.hpp"
38 : #include "Domain/Tags/FaceNormal.hpp"
39 : #include "Domain/Tags/Faces.hpp"
40 : #include "Domain/Tags/NeighborMesh.hpp"
41 : #include "Domain/Tags/SurfaceJacobian.hpp"
42 : #include "Elliptic/DiscontinuousGalerkin/Penalty.hpp"
43 : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
44 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
45 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
46 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
47 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
48 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
49 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
50 : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
51 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
52 : #include "Utilities/Algorithm.hpp"
53 : #include "Utilities/CallWithDynamicType.hpp"
54 : #include "Utilities/ErrorHandling/Assert.hpp"
55 : #include "Utilities/Gsl.hpp"
56 : #include "Utilities/ProtocolHelpers.hpp"
57 : #include "Utilities/TMPL.hpp"
58 :
59 : namespace elliptic::dg {
60 :
61 : namespace detail {
62 : template <size_t Dim>
63 : void initialize_coords_and_jacobians(
64 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
65 : logical_coords,
66 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> inertial_coords,
67 : gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
68 : Frame::Inertial>*>
69 : inv_jacobian,
70 : gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
71 : gsl::not_null<Scalar<DataVector>*> det_jacobian,
72 : gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
73 : Frame::Inertial>*>
74 : det_times_inv_jacobian,
75 : const Mesh<Dim>& mesh, const ElementMap<Dim, Frame::Inertial>& element_map,
76 : const domain::FunctionsOfTimeMap& functions_of_time);
77 : } // namespace detail
78 :
79 : /*!
80 : * \brief Initialize the background-independent geometry for the elliptic DG
81 : * operator
82 : *
83 : * ## Geometric aliasing
84 : *
85 : * The geometric quantities such as Jacobians are evaluated on the DG grid.
86 : * Since we know them analytically, we could also evaluate them on a
87 : * higher-order grid or with a stronger quadrature (Gauss instead of
88 : * Gauss-Lobatto) to combat geometric aliasing. See discussions in
89 : * \cite Vincent2019qpd and \cite Fischer2021voj .
90 : */
91 : template <size_t Dim>
92 1 : struct InitializeGeometry {
93 0 : using return_tags = tmpl::list<
94 : domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
95 : domain::Tags::NeighborMesh<Dim>, domain::Tags::ElementMap<Dim>,
96 : domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
97 : domain::Tags::Coordinates<Dim, Frame::Inertial>,
98 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
99 : Frame::Inertial>,
100 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
101 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>,
102 : domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
103 : Frame::Inertial>>;
104 0 : using argument_tags =
105 : tmpl::list<domain::Tags::InitialExtents<Dim>,
106 : domain::Tags::InitialRefinementLevels<Dim>,
107 : domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime,
108 : elliptic::dg::Tags::Quadrature>;
109 0 : using volume_tags = argument_tags;
110 0 : static void apply(
111 : gsl::not_null<Mesh<Dim>*> mesh, gsl::not_null<Element<Dim>*> element,
112 : gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*> neighbor_meshes,
113 : gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
114 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
115 : logical_coords,
116 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> inertial_coords,
117 : gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
118 : Frame::Inertial>*>
119 : inv_jacobian,
120 : gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
121 : gsl::not_null<Scalar<DataVector>*> det_jacobian,
122 : gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
123 : Frame::Inertial>*>
124 : det_times_inv_jacobian,
125 : const std::vector<std::array<size_t, Dim>>& initial_extents,
126 : const std::vector<std::array<size_t, Dim>>& initial_refinement,
127 : const Domain<Dim>& domain,
128 : const domain::FunctionsOfTimeMap& functions_of_time,
129 : Spectral::Quadrature i1_quadrature, const ElementId<Dim>& element_id);
130 : };
131 :
132 : template <size_t Dim>
133 0 : struct ProjectGeometry : tt::ConformsTo<::amr::protocols::Projector> {
134 0 : using return_tags = tmpl::list<
135 : domain::Tags::ElementMap<Dim>,
136 : domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
137 : domain::Tags::Coordinates<Dim, Frame::Inertial>,
138 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
139 : Frame::Inertial>,
140 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
141 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>,
142 : domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
143 : Frame::Inertial>>;
144 0 : using argument_tags =
145 : tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
146 : domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime>;
147 0 : using volume_tags =
148 : tmpl::list<domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime>;
149 :
150 : // p-refinement
151 0 : static void apply(
152 : const gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
153 : const gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
154 : logical_coords,
155 : const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
156 : inertial_coords,
157 : const gsl::not_null<InverseJacobian<
158 : DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
159 : inv_jacobian,
160 : const gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
161 : const gsl::not_null<Scalar<DataVector>*> det_jacobian,
162 : const gsl::not_null<InverseJacobian<
163 : DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
164 : det_times_inv_jacobian,
165 : const Mesh<Dim>& mesh, const Element<Dim>& /*element*/,
166 : const Domain<Dim>& /*domain*/,
167 : const domain::FunctionsOfTimeMap& functions_of_time,
168 : const std::pair<Mesh<Dim>, Element<Dim>>& old_mesh_and_element) {
169 : if (mesh == old_mesh_and_element.first) {
170 : // Neighbors may have changed, but this element hasn't. Nothing to do.
171 : return;
172 : }
173 : // The element map doesn't change with p-refinement
174 : detail::initialize_coords_and_jacobians(
175 : logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian,
176 : det_jacobian, det_times_inv_jacobian, mesh, *element_map,
177 : functions_of_time);
178 : }
179 :
180 : // h-refinement
181 : template <typename... ParentOrChildrenItemsType>
182 0 : static void apply(
183 : const gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
184 : const gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
185 : logical_coords,
186 : const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
187 : inertial_coords,
188 : const gsl::not_null<InverseJacobian<
189 : DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
190 : inv_jacobian,
191 : const gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
192 : const gsl::not_null<Scalar<DataVector>*> det_jacobian,
193 : const gsl::not_null<InverseJacobian<
194 : DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
195 : det_times_inv_jacobian,
196 : const Mesh<Dim>& mesh, const Element<Dim>& element,
197 : const Domain<Dim>& domain,
198 : const domain::FunctionsOfTimeMap& functions_of_time,
199 : const ParentOrChildrenItemsType&... /*parent_or_children_items*/) {
200 : const auto& element_id = element.id();
201 : const auto& block = domain.blocks()[element_id.block_id()];
202 : *element_map = ElementMap<Dim, Frame::Inertial>{element_id, block};
203 : detail::initialize_coords_and_jacobians(
204 : logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian,
205 : det_jacobian, det_times_inv_jacobian, mesh, *element_map,
206 : functions_of_time);
207 : }
208 : };
209 :
210 : namespace detail {
211 : // Compute derivative of the Jacobian numerically
212 : template <size_t Dim>
213 : void deriv_unnormalized_face_normals_impl(
214 : gsl::not_null<DirectionMap<Dim, tnsr::ij<DataVector, Dim>>*>
215 : deriv_unnormalized_face_normals,
216 : const Mesh<Dim>& mesh, const Element<Dim>& element,
217 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
218 : Frame::Inertial>& inv_jacobian);
219 :
220 : // Get element-logical coordinates of the mortar collocation points
221 : template <size_t Dim>
222 : tnsr::I<DataVector, Dim, Frame::ElementLogical> mortar_logical_coordinates(
223 : const Mesh<Dim - 1>& mortar_mesh,
224 : const ::dg::MortarSize<Dim - 1>& mortar_size,
225 : const Direction<Dim>& direction);
226 :
227 : template <size_t Dim>
228 : void mortar_jacobian(
229 : gsl::not_null<Scalar<DataVector>*> mortar_jacobian,
230 : gsl::not_null<Scalar<DataVector>*> perpendicular_element_size,
231 : const Mesh<Dim - 1>& mortar_mesh,
232 : const ::dg::MortarSize<Dim - 1>& mortar_size,
233 : const Direction<Dim>& direction,
234 : const tnsr::I<DataVector, Dim, Frame::ElementLogical>&
235 : mortar_logical_coords,
236 : const std::optional<tnsr::II<DataVector, Dim>>& inv_metric_on_mortar,
237 : const ElementMap<Dim, Frame::Inertial>& element_map,
238 : const domain::FunctionsOfTimeMap& functions_of_time);
239 : } // namespace detail
240 :
241 : /// Initialize the geometry on faces and mortars for the elliptic DG operator
242 : ///
243 : /// To normalize face normals this function needs the inverse background metric.
244 : /// Pass the tag representing the inverse background metric to the
245 : /// `InvMetricTag` template parameter, and the tag representing the analytic
246 : /// background from which it can be retrieved to the `BackgroundTag` template
247 : /// parameter. Set `InvMetricTag` and `BackgroundTag` to `void` to normalize
248 : /// face normals with the Euclidean magnitude.
249 : ///
250 : /// Mortar Jacobians are added only on nonconforming internal element
251 : /// boundaries, i.e., when `Spectral::needs_projection()` is true.
252 : ///
253 : /// The `::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<Dim>>` is only added
254 : /// on external boundaries, for use by boundary conditions.
255 : template <size_t Dim, typename InvMetricTag, typename BackgroundTag>
256 1 : struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> {
257 0 : using return_tags = tmpl::append<
258 : domain::make_faces_tags<
259 : Dim,
260 : tmpl::list<domain::Tags::Direction<Dim>,
261 : domain::Tags::Coordinates<Dim, Frame::Inertial>,
262 : domain::Tags::FaceNormal<Dim>,
263 : domain::Tags::FaceNormalVector<Dim>,
264 : domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>,
265 : domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
266 : Frame::Inertial>,
267 : // This is the volume inverse Jacobian on the face grid
268 : // points, multiplied by the determinant of the _face_
269 : // Jacobian (the tag above)
270 : domain::Tags::DetTimesInvJacobian<
271 : Dim, Frame::ElementLogical, Frame::Inertial>,
272 : // Possible optimization: The derivative of the face normal
273 : // could be omitted for some systems, but its memory usage
274 : // is probably insignificant since it's only added on
275 : // external boundaries.
276 : ::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<Dim>,
277 : tmpl::size_t<Dim>, Frame::Inertial>>>,
278 : tmpl::list<
279 : ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
280 : ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
281 : ::Tags::Mortars<domain::Tags::Coordinates<Dim, Frame::Inertial>, Dim>,
282 : ::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
283 : Frame::ElementLogical, Frame::Inertial>,
284 : Dim>,
285 : ::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>>;
286 0 : using argument_tags = tmpl::append<
287 : tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
288 : domain::Tags::NeighborMesh<Dim>, domain::Tags::ElementMap<Dim>,
289 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
290 : Frame::Inertial>,
291 : domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime,
292 : elliptic::dg::Tags::PenaltyParameter>,
293 : tmpl::conditional_t<
294 : std::is_same_v<BackgroundTag, void>, tmpl::list<>,
295 : tmpl::list<BackgroundTag, Parallel::Tags::Metavariables>>>;
296 0 : using volume_tags = tmpl::append<
297 : tmpl::list<domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime,
298 : elliptic::dg::Tags::PenaltyParameter>,
299 : tmpl::conditional_t<
300 : std::is_same_v<BackgroundTag, void>, tmpl::list<>,
301 : tmpl::list<BackgroundTag, Parallel::Tags::Metavariables>>>;
302 : template <typename... AmrData>
303 0 : static void apply(
304 : const gsl::not_null<DirectionMap<Dim, Direction<Dim>>*> face_directions,
305 : const gsl::not_null<DirectionMap<Dim, tnsr::I<DataVector, Dim>>*>
306 : faces_inertial_coords,
307 : const gsl::not_null<DirectionMap<Dim, tnsr::i<DataVector, Dim>>*>
308 : face_normals,
309 : const gsl::not_null<DirectionMap<Dim, tnsr::I<DataVector, Dim>>*>
310 : face_normal_vectors,
311 : const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
312 : face_normal_magnitudes,
313 : const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
314 : face_jacobians,
315 : const gsl::not_null<DirectionMap<
316 : Dim, InverseJacobian<DataVector, Dim, Frame::ElementLogical,
317 : Frame::Inertial>>*>
318 : face_jacobian_times_inv_jacobian,
319 : const gsl::not_null<DirectionMap<Dim, tnsr::ij<DataVector, Dim>>*>
320 : deriv_unnormalized_face_normals,
321 : const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_meshes,
322 : const gsl::not_null<::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>*>
323 : mortar_sizes,
324 : const gsl::not_null<::dg::MortarMap<Dim, tnsr::I<DataVector, Dim>>*>
325 : all_mortar_inertial_coords,
326 : const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
327 : mortar_jacobians,
328 : const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
329 : penalty_factors,
330 : const Mesh<Dim>& mesh, const Element<Dim>& element,
331 : const DirectionalIdMap<Dim, Mesh<Dim>>& neighbor_meshes,
332 : const ElementMap<Dim, Frame::Inertial>& element_map,
333 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
334 : Frame::Inertial>& inv_jacobian,
335 : const Domain<Dim>& domain,
336 : const domain::FunctionsOfTimeMap& functions_of_time,
337 : const double penalty_parameter, const AmrData&... amr_data) {
338 : apply(face_directions, faces_inertial_coords, face_normals,
339 : face_normal_vectors, face_normal_magnitudes, face_jacobians,
340 : face_jacobian_times_inv_jacobian, deriv_unnormalized_face_normals,
341 : mortar_meshes, mortar_sizes, all_mortar_inertial_coords,
342 : mortar_jacobians, penalty_factors, mesh, element, neighbor_meshes,
343 : element_map, inv_jacobian, domain, functions_of_time,
344 : penalty_parameter, nullptr, nullptr, amr_data...);
345 : }
346 : template <typename Background, typename Metavariables, typename... AmrData>
347 0 : static void apply(
348 : const gsl::not_null<DirectionMap<Dim, Direction<Dim>>*> face_directions,
349 : const gsl::not_null<DirectionMap<Dim, tnsr::I<DataVector, Dim>>*>
350 : faces_inertial_coords,
351 : const gsl::not_null<DirectionMap<Dim, tnsr::i<DataVector, Dim>>*>
352 : face_normals,
353 : const gsl::not_null<DirectionMap<Dim, tnsr::I<DataVector, Dim>>*>
354 : face_normal_vectors,
355 : const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
356 : face_normal_magnitudes,
357 : const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
358 : face_jacobians,
359 : const gsl::not_null<DirectionMap<
360 : Dim, InverseJacobian<DataVector, Dim, Frame::ElementLogical,
361 : Frame::Inertial>>*>
362 : face_jacobian_times_inv_jacobian,
363 : const gsl::not_null<DirectionMap<Dim, tnsr::ij<DataVector, Dim>>*>
364 : deriv_unnormalized_face_normals,
365 : const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_meshes,
366 : const gsl::not_null<::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>*>
367 : mortar_sizes,
368 : const gsl::not_null<::dg::MortarMap<Dim, tnsr::I<DataVector, Dim>>*>
369 : all_mortar_inertial_coords,
370 : const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
371 : mortar_jacobians,
372 : const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
373 : penalty_factors,
374 : const Mesh<Dim>& mesh, const Element<Dim>& element,
375 : const DirectionalIdMap<Dim, Mesh<Dim>>& neighbor_meshes,
376 : const ElementMap<Dim, Frame::Inertial>& element_map,
377 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
378 : Frame::Inertial>& inv_jacobian,
379 : const Domain<Dim>& domain,
380 : const domain::FunctionsOfTimeMap& functions_of_time,
381 : const double penalty_parameter, const Background& background,
382 : const Metavariables& /*meta*/, const AmrData&... /*amr_data*/) {
383 : static_assert(std::is_same_v<InvMetricTag, void> or
384 : not(std::is_same_v<Background, std::nullptr_t>),
385 : "Supply an analytic background from which the 'InvMetricTag' "
386 : "can be retrieved");
387 : [[maybe_unused]] const auto get_inv_metric =
388 : [&background]([[maybe_unused]] const tnsr::I<DataVector, Dim>&
389 : local_inertial_coords)
390 : -> std::optional<tnsr::II<DataVector, Dim>> {
391 : if constexpr (not std::is_same_v<InvMetricTag, void>) {
392 : using factory_classes = typename std::decay_t<
393 : Metavariables>::factory_creation::factory_classes;
394 : return call_with_dynamic_type<tnsr::II<DataVector, Dim>,
395 : tmpl::at<factory_classes, Background>>(
396 : &background, [&local_inertial_coords](const auto* const derived) {
397 : return get<InvMetricTag>(derived->variables(
398 : local_inertial_coords, tmpl::list<InvMetricTag>{}));
399 : });
400 : } else {
401 : (void)background;
402 : return std::nullopt;
403 : }
404 : };
405 : ASSERT(
406 : alg::all_of(mesh.quadrature(),
407 : [&mesh](const auto t) { return t == mesh.quadrature(0); }),
408 : "This function is implemented assuming the quadrature is isotropic");
409 : // Faces
410 : for (const auto& direction : Direction<Dim>::all_directions()) {
411 : const auto face_mesh = mesh.slice_away(direction.dimension());
412 : (*face_directions)[direction] = direction;
413 : // Possible optimization: Not all systems need the coordinates on internal
414 : // faces.
415 : const auto face_logical_coords =
416 : interface_logical_coordinates(face_mesh, direction);
417 : auto& face_inertial_coords = (*faces_inertial_coords)[direction];
418 : face_inertial_coords =
419 : element_map(face_logical_coords, 0., functions_of_time);
420 : auto& face_normal = (*face_normals)[direction];
421 : auto& face_normal_vector = (*face_normal_vectors)[direction];
422 : auto& face_normal_magnitude = (*face_normal_magnitudes)[direction];
423 : // Buffer the inv Jacobian on the face here, then multiply by the face
424 : // Jacobian below
425 : auto& inv_jacobian_on_face =
426 : (*face_jacobian_times_inv_jacobian)[direction];
427 : inv_jacobian_on_face =
428 : element_map.inv_jacobian(face_logical_coords, 0., functions_of_time);
429 : unnormalized_face_normal(make_not_null(&face_normal), face_mesh,
430 : inv_jacobian_on_face, direction);
431 : if constexpr (std::is_same_v<InvMetricTag, void>) {
432 : magnitude(make_not_null(&face_normal_magnitude), face_normal);
433 : for (size_t d = 0; d < Dim; ++d) {
434 : face_normal.get(d) /= get(face_normal_magnitude);
435 : face_normal_vector.get(d) = face_normal.get(d);
436 : }
437 : } else {
438 : const auto inv_metric_on_face = *get_inv_metric(face_inertial_coords);
439 : magnitude(make_not_null(&face_normal_magnitude), face_normal,
440 : inv_metric_on_face);
441 : for (size_t d = 0; d < Dim; ++d) {
442 : face_normal.get(d) /= get(face_normal_magnitude);
443 : }
444 : raise_or_lower_index(make_not_null(&face_normal_vector), face_normal,
445 : inv_metric_on_face);
446 : }
447 : auto& face_jacobian = (*face_jacobians)[direction];
448 : get(face_jacobian) =
449 : get(face_normal_magnitude) / get(determinant(inv_jacobian_on_face));
450 : for (auto& component : inv_jacobian_on_face) {
451 : component *= get(face_jacobian);
452 : }
453 : }
454 : // Compute the Jacobian derivative numerically, because our coordinate maps
455 : // currently don't provide it analytically.
456 : detail::deriv_unnormalized_face_normals_impl(
457 : deriv_unnormalized_face_normals, mesh, element, inv_jacobian);
458 : // Mortars (internal directions)
459 : mortar_meshes->clear();
460 : mortar_sizes->clear();
461 : all_mortar_inertial_coords->clear();
462 : mortar_jacobians->clear();
463 : penalty_factors->clear();
464 : const auto& element_id = element.id();
465 : for (const auto& [direction, neighbors] : element.neighbors()) {
466 : const auto face_mesh = mesh.slice_away(direction.dimension());
467 : ASSERT(element.face_types().at(direction) !=
468 : domain::FaceType::MultipleNonconforming,
469 : "This code needs updating to handle nonconforming blocks.\n");
470 : for (const auto& neighbor_id : neighbors) {
471 : const auto& orientation = neighbors.orientation(neighbor_id);
472 : const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
473 : const auto& neighbor_mesh = neighbor_meshes.at(mortar_id);
474 : mortar_meshes->emplace(
475 : mortar_id,
476 : ::dg::mortar_mesh(face_mesh,
477 : neighbor_mesh.slice_away(direction.dimension())));
478 : mortar_sizes->emplace(
479 : mortar_id, ::dg::mortar_size(element_id, neighbor_id,
480 : direction.dimension(), orientation));
481 : // Mortar Jacobian
482 : const auto& mortar_mesh = mortar_meshes->at(mortar_id);
483 : const auto& mortar_size = mortar_sizes->at(mortar_id);
484 : const auto mortar_logical_coords = detail::mortar_logical_coordinates(
485 : mortar_mesh, mortar_size, direction);
486 : const auto& mortar_inertial_coords =
487 : all_mortar_inertial_coords
488 : ->emplace(mortar_id, element_map(mortar_logical_coords, 0.,
489 : functions_of_time))
490 : .first->second;
491 : Scalar<DataVector> perpendicular_element_size{};
492 : if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
493 : auto& mortar_jacobian = (*mortar_jacobians)[mortar_id];
494 : detail::mortar_jacobian(make_not_null(&mortar_jacobian),
495 : make_not_null(&perpendicular_element_size),
496 : mortar_mesh, mortar_size, direction,
497 : mortar_logical_coords,
498 : get_inv_metric(mortar_inertial_coords),
499 : element_map, functions_of_time);
500 : } else {
501 : // Mortar is identical to face, and we have computed the face normal
502 : // magnitude already above
503 : get(perpendicular_element_size) =
504 : 2. / get(face_normal_magnitudes->at(direction));
505 : }
506 : // Penalty factor
507 : // The penalty factor (like all quantities on mortars) must agree when
508 : // calculated on both sides of the mortar. So we switch perspective to
509 : // the neighbor here.
510 : const auto direction_in_neighbor = orientation(direction).opposite();
511 : const auto reoriented_mortar_mesh = orient_mesh_on_slice(
512 : mortar_mesh, direction.dimension(), orientation);
513 : const auto mortar_size_in_neighbor = ::dg::mortar_size(
514 : neighbor_id, element_id, direction_in_neighbor.dimension(),
515 : orientation.inverse_map());
516 : const auto mortar_logical_coords_in_neighbor =
517 : detail::mortar_logical_coordinates(reoriented_mortar_mesh,
518 : mortar_size_in_neighbor,
519 : direction_in_neighbor);
520 : const ElementMap<Dim, Frame::Inertial> neighbor_element_map{
521 : neighbor_id, domain.blocks()[neighbor_id.block_id()]};
522 : const auto reoriented_mortar_inertial_coords = neighbor_element_map(
523 : mortar_logical_coords_in_neighbor, 0., functions_of_time);
524 : Scalar<DataVector> buffer{};
525 : Scalar<DataVector> reoriented_neighbor_element_size{};
526 : detail::mortar_jacobian(
527 : make_not_null(&buffer),
528 : make_not_null(&reoriented_neighbor_element_size),
529 : reoriented_mortar_mesh, mortar_size_in_neighbor,
530 : direction_in_neighbor, mortar_logical_coords_in_neighbor,
531 : get_inv_metric(reoriented_mortar_inertial_coords),
532 : neighbor_element_map, functions_of_time);
533 : // Orient the result back to the perspective of this element
534 : const auto neighbor_element_size = orient_variables_on_slice(
535 : get(reoriented_neighbor_element_size),
536 : reoriented_mortar_mesh.extents(), direction_in_neighbor.dimension(),
537 : orientation.inverse_map());
538 : penalty_factors->emplace(
539 : mortar_id,
540 : elliptic::dg::penalty(
541 : blaze::min(get(perpendicular_element_size),
542 : neighbor_element_size),
543 : std::max(mesh.extents(direction.dimension()),
544 : neighbor_mesh.extents(direction.dimension())),
545 : penalty_parameter));
546 : } // neighbors
547 : } // internal directions
548 : // Mortars (external directions)
549 : for (const auto& direction : element.external_boundaries()) {
550 : const auto face_mesh = mesh.slice_away(direction.dimension());
551 : const auto mortar_id =
552 : DirectionalId<Dim>{direction, ElementId<Dim>::external_boundary_id()};
553 : mortar_meshes->emplace(mortar_id, face_mesh);
554 : mortar_sizes->emplace(mortar_id,
555 : make_array<Dim - 1>(Spectral::SegmentSize::Full));
556 : penalty_factors->emplace(
557 : mortar_id,
558 : elliptic::dg::penalty(2. / get(face_normal_magnitudes->at(direction)),
559 : mesh.extents(direction.dimension()),
560 : penalty_parameter));
561 : } // external directions
562 : }
563 : };
564 :
565 : /// Initialize background quantities for the elliptic DG operator, possibly
566 : /// including the metric necessary for normalizing face normals
567 : template <size_t Dim, typename BackgroundFields, typename BackgroundTag>
568 1 : struct InitializeBackground : tt::ConformsTo<::amr::protocols::Projector> {
569 0 : using return_tags =
570 : tmpl::list<::Tags::Variables<BackgroundFields>,
571 : domain::Tags::Faces<Dim, ::Tags::Variables<BackgroundFields>>>;
572 0 : using argument_tags =
573 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
574 : domain::Tags::Mesh<Dim>,
575 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
576 : Frame::Inertial>,
577 : BackgroundTag, Parallel::Tags::Metavariables>;
578 :
579 : template <typename BackgroundBase, typename Metavariables,
580 : typename... AmrData>
581 0 : static void apply(
582 : const gsl::not_null<Variables<BackgroundFields>*> background_fields,
583 : const gsl::not_null<DirectionMap<Dim, Variables<BackgroundFields>>*>
584 : face_background_fields,
585 : const tnsr::I<DataVector, Dim>& inertial_coords, const Mesh<Dim>& mesh,
586 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
587 : Frame::Inertial>& inv_jacobian,
588 : const BackgroundBase& background, const Metavariables& /*meta*/,
589 : const AmrData&... amr_data) {
590 : if constexpr (sizeof...(AmrData) == 1) {
591 : if constexpr (std::is_same_v<AmrData...,
592 : std::pair<Mesh<Dim>, Element<Dim>>>) {
593 : if (((mesh == amr_data.first) and ...)) {
594 : // This element hasn't changed during AMR. Nothing to do.
595 : return;
596 : }
597 : }
598 : }
599 :
600 : // Background fields in the volume
601 : using factory_classes =
602 : typename std::decay_t<Metavariables>::factory_creation::factory_classes;
603 : *background_fields =
604 : call_with_dynamic_type<Variables<BackgroundFields>,
605 : tmpl::at<factory_classes, BackgroundBase>>(
606 : &background, [&inertial_coords, &mesh,
607 : &inv_jacobian](const auto* const derived) {
608 : return variables_from_tagged_tuple(derived->variables(
609 : inertial_coords, mesh, inv_jacobian, BackgroundFields{}));
610 : });
611 : // Background fields on faces
612 : for (const auto& direction : Direction<Dim>::all_directions()) {
613 : // Possible optimization: Only the background fields in the
614 : // System::fluxes_computer::argument_tags are needed on internal faces.
615 : // On Gauss grids we could evaluate the background directly on the faces
616 : // instead of projecting the data. However, we need to take derivatives of
617 : // the background fields, so we could evaluate them on a Gauss-Lobatto
618 : // grid in the volume. We could even evaluate the background fields on a
619 : // higher-order grid and project down to get more accurate derivatives if
620 : // needed.
621 : (*face_background_fields)[direction].initialize(
622 : mesh.slice_away(direction.dimension()).number_of_grid_points());
623 : ::dg::project_contiguous_data_to_boundary(
624 : make_not_null(&(*face_background_fields)[direction]),
625 : *background_fields, mesh, direction);
626 : }
627 : }
628 : };
629 :
630 : } // namespace elliptic::dg
|
