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