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