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/Tensor/EagerMath/Magnitude.hpp"
16 : #include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "DataStructures/VariablesTag.hpp"
20 : #include "Domain/Creators/Tags/Domain.hpp"
21 : #include "Domain/Creators/Tags/InitialExtents.hpp"
22 : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
23 : #include "Domain/Domain.hpp"
24 : #include "Domain/ElementMap.hpp"
25 : #include "Domain/FaceNormal.hpp"
26 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
27 : #include "Domain/FunctionsOfTime/Tags.hpp"
28 : #include "Domain/InterfaceLogicalCoordinates.hpp"
29 : #include "Domain/Structure/CreateInitialMesh.hpp"
30 : #include "Domain/Structure/Direction.hpp"
31 : #include "Domain/Structure/DirectionMap.hpp"
32 : #include "Domain/Structure/DirectionalId.hpp"
33 : #include "Domain/Structure/Element.hpp"
34 : #include "Domain/Structure/IndexToSliceAt.hpp"
35 : #include "Domain/Structure/OrientationMapHelpers.hpp"
36 : #include "Domain/Tags.hpp"
37 : #include "Domain/Tags/FaceNormal.hpp"
38 : #include "Domain/Tags/Faces.hpp"
39 : #include "Domain/Tags/NeighborMesh.hpp"
40 : #include "Domain/Tags/SurfaceJacobian.hpp"
41 : #include "Elliptic/DiscontinuousGalerkin/Penalty.hpp"
42 : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
43 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
44 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
45 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
46 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
47 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
48 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
49 : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
50 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
51 : #include "Utilities/Algorithm.hpp"
52 : #include "Utilities/CallWithDynamicType.hpp"
53 : #include "Utilities/ErrorHandling/Assert.hpp"
54 : #include "Utilities/Gsl.hpp"
55 : #include "Utilities/ProtocolHelpers.hpp"
56 : #include "Utilities/TMPL.hpp"
57 : #include "Utilities/TaggedTuple.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 : for (const auto& neighbor_id : neighbors) {
468 : const auto& orientation = neighbors.orientation(neighbor_id);
469 : const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
470 : const auto& neighbor_mesh = neighbor_meshes.at(mortar_id);
471 : mortar_meshes->emplace(
472 : mortar_id,
473 : ::dg::mortar_mesh(face_mesh,
474 : neighbor_mesh.slice_away(direction.dimension())));
475 : mortar_sizes->emplace(
476 : mortar_id, ::dg::mortar_size(element_id, neighbor_id,
477 : direction.dimension(), orientation));
478 : // Mortar Jacobian
479 : const auto& mortar_mesh = mortar_meshes->at(mortar_id);
480 : const auto& mortar_size = mortar_sizes->at(mortar_id);
481 : const auto mortar_logical_coords = detail::mortar_logical_coordinates(
482 : mortar_mesh, mortar_size, direction);
483 : const auto& mortar_inertial_coords =
484 : all_mortar_inertial_coords
485 : ->emplace(mortar_id, element_map(mortar_logical_coords, 0.,
486 : functions_of_time))
487 : .first->second;
488 : Scalar<DataVector> perpendicular_element_size{};
489 : if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
490 : auto& mortar_jacobian = (*mortar_jacobians)[mortar_id];
491 : detail::mortar_jacobian(make_not_null(&mortar_jacobian),
492 : make_not_null(&perpendicular_element_size),
493 : mortar_mesh, mortar_size, direction,
494 : mortar_logical_coords,
495 : get_inv_metric(mortar_inertial_coords),
496 : element_map, functions_of_time);
497 : } else {
498 : // Mortar is identical to face, and we have computed the face normal
499 : // magnitude already above
500 : get(perpendicular_element_size) =
501 : 2. / get(face_normal_magnitudes->at(direction));
502 : }
503 : // Penalty factor
504 : // The penalty factor (like all quantities on mortars) must agree when
505 : // calculated on both sides of the mortar. So we switch perspective to
506 : // the neighbor here.
507 : const auto direction_in_neighbor = orientation(direction).opposite();
508 : const auto reoriented_mortar_mesh = orient_mesh_on_slice(
509 : mortar_mesh, direction.dimension(), orientation);
510 : const auto mortar_size_in_neighbor = ::dg::mortar_size(
511 : neighbor_id, element_id, direction_in_neighbor.dimension(),
512 : orientation.inverse_map());
513 : const auto mortar_logical_coords_in_neighbor =
514 : detail::mortar_logical_coordinates(reoriented_mortar_mesh,
515 : mortar_size_in_neighbor,
516 : direction_in_neighbor);
517 : const ElementMap<Dim, Frame::Inertial> neighbor_element_map{
518 : neighbor_id, domain.blocks()[neighbor_id.block_id()]};
519 : const auto reoriented_mortar_inertial_coords = neighbor_element_map(
520 : mortar_logical_coords_in_neighbor, 0., functions_of_time);
521 : Scalar<DataVector> buffer{};
522 : Scalar<DataVector> reoriented_neighbor_element_size{};
523 : detail::mortar_jacobian(
524 : make_not_null(&buffer),
525 : make_not_null(&reoriented_neighbor_element_size),
526 : reoriented_mortar_mesh, mortar_size_in_neighbor,
527 : direction_in_neighbor, mortar_logical_coords_in_neighbor,
528 : get_inv_metric(reoriented_mortar_inertial_coords),
529 : neighbor_element_map, functions_of_time);
530 : // Orient the result back to the perspective of this element
531 : const auto neighbor_element_size = orient_variables_on_slice(
532 : get(reoriented_neighbor_element_size),
533 : reoriented_mortar_mesh.extents(), direction_in_neighbor.dimension(),
534 : orientation.inverse_map());
535 : penalty_factors->emplace(
536 : mortar_id,
537 : elliptic::dg::penalty(
538 : blaze::min(get(perpendicular_element_size),
539 : neighbor_element_size),
540 : std::max(mesh.extents(direction.dimension()),
541 : neighbor_mesh.extents(direction.dimension())),
542 : penalty_parameter));
543 : } // neighbors
544 : } // internal directions
545 : // Mortars (external directions)
546 : for (const auto& direction : element.external_boundaries()) {
547 : const auto face_mesh = mesh.slice_away(direction.dimension());
548 : const auto mortar_id =
549 : DirectionalId<Dim>{direction, ElementId<Dim>::external_boundary_id()};
550 : mortar_meshes->emplace(mortar_id, face_mesh);
551 : mortar_sizes->emplace(mortar_id,
552 : make_array<Dim - 1>(Spectral::SegmentSize::Full));
553 : penalty_factors->emplace(
554 : mortar_id,
555 : elliptic::dg::penalty(2. / get(face_normal_magnitudes->at(direction)),
556 : mesh.extents(direction.dimension()),
557 : penalty_parameter));
558 : } // external directions
559 : }
560 : };
561 :
562 : /// Initialize background quantities for the elliptic DG operator, possibly
563 : /// including the metric necessary for normalizing face normals
564 : template <size_t Dim, typename BackgroundFields, typename BackgroundTag>
565 1 : struct InitializeBackground : tt::ConformsTo<::amr::protocols::Projector> {
566 0 : using return_tags =
567 : tmpl::list<::Tags::Variables<BackgroundFields>,
568 : domain::Tags::Faces<Dim, ::Tags::Variables<BackgroundFields>>>;
569 0 : using argument_tags =
570 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
571 : domain::Tags::Mesh<Dim>,
572 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
573 : Frame::Inertial>,
574 : BackgroundTag, Parallel::Tags::Metavariables>;
575 :
576 : template <typename BackgroundBase, typename Metavariables,
577 : typename... AmrData>
578 0 : static void apply(
579 : const gsl::not_null<Variables<BackgroundFields>*> background_fields,
580 : const gsl::not_null<DirectionMap<Dim, Variables<BackgroundFields>>*>
581 : face_background_fields,
582 : const tnsr::I<DataVector, Dim>& inertial_coords, const Mesh<Dim>& mesh,
583 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
584 : Frame::Inertial>& inv_jacobian,
585 : const BackgroundBase& background, const Metavariables& /*meta*/,
586 : const AmrData&... amr_data) {
587 : if constexpr (sizeof...(AmrData) == 1) {
588 : if constexpr (std::is_same_v<AmrData...,
589 : std::pair<Mesh<Dim>, Element<Dim>>>) {
590 : if (((mesh == amr_data.first) and ...)) {
591 : // This element hasn't changed during AMR. Nothing to do.
592 : return;
593 : }
594 : }
595 : }
596 :
597 : // Background fields in the volume
598 : using factory_classes =
599 : typename std::decay_t<Metavariables>::factory_creation::factory_classes;
600 : *background_fields =
601 : call_with_dynamic_type<Variables<BackgroundFields>,
602 : tmpl::at<factory_classes, BackgroundBase>>(
603 : &background, [&inertial_coords, &mesh,
604 : &inv_jacobian](const auto* const derived) {
605 : return variables_from_tagged_tuple(derived->variables(
606 : inertial_coords, mesh, inv_jacobian, BackgroundFields{}));
607 : });
608 : // Background fields on faces
609 : for (const auto& direction : Direction<Dim>::all_directions()) {
610 : // Possible optimization: Only the background fields in the
611 : // System::fluxes_computer::argument_tags are needed on internal faces.
612 : // On Gauss grids we could evaluate the background directly on the faces
613 : // instead of projecting the data. However, we need to take derivatives of
614 : // the background fields, so we could evaluate them on a Gauss-Lobatto
615 : // grid in the volume. We could even evaluate the background fields on a
616 : // higher-order grid and project down to get more accurate derivatives if
617 : // needed.
618 : (*face_background_fields)[direction].initialize(
619 : mesh.slice_away(direction.dimension()).number_of_grid_points());
620 : ::dg::project_contiguous_data_to_boundary(
621 : make_not_null(&(*face_background_fields)[direction]),
622 : *background_fields, mesh, direction);
623 : }
624 : }
625 : };
626 :
627 : } // namespace elliptic::dg
|
