Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <string>
8 : #include <tuple>
9 : #include <unordered_map>
10 : #include <utility>
11 :
12 : #include "DataStructures/DataBox/Tag.hpp"
13 : #include "DataStructures/DataVector.hpp"
14 : #include "DataStructures/SliceVariables.hpp"
15 : #include "DataStructures/TaggedTuple.hpp"
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "DataStructures/Variables/FrameTransform.hpp"
19 : #include "Domain/Structure/Direction.hpp"
20 : #include "Domain/Structure/DirectionMap.hpp"
21 : #include "Domain/Structure/Element.hpp"
22 : #include "Domain/Structure/ElementId.hpp"
23 : #include "Domain/Structure/IndexToSliceAt.hpp"
24 : #include "Elliptic/Protocols/FirstOrderSystem.hpp"
25 : #include "Elliptic/Systems/GetFluxesComputer.hpp"
26 : #include "Elliptic/Systems/GetSourcesComputer.hpp"
27 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp"
28 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
29 : #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFlux.hpp"
30 : #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFromBoundary.hpp"
31 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
32 : #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
33 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
34 : #include "NumericalAlgorithms/DiscontinuousGalerkin/SimpleBoundaryData.hpp"
35 : #include "NumericalAlgorithms/DiscontinuousGalerkin/SimpleMortarData.hpp"
36 : #include "NumericalAlgorithms/LinearOperators/Divergence.hpp"
37 : #include "NumericalAlgorithms/LinearOperators/Divergence.tpp"
38 : #include "NumericalAlgorithms/LinearOperators/WeakDivergence.hpp"
39 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
40 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
41 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
42 : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
43 : #include "Utilities/ErrorHandling/Assert.hpp"
44 : #include "Utilities/Gsl.hpp"
45 : #include "Utilities/ProtocolHelpers.hpp"
46 : #include "Utilities/TMPL.hpp"
47 :
48 : /*!
49 : * \brief Functionality related to discontinuous Galerkin discretizations of
50 : * elliptic equations
51 : *
52 : * The following is a brief overview of the elliptic DG schemes that are
53 : * implemented here. The scheme is described in detail in \cite Fischer2021voj.
54 : *
55 : * The DG schemes apply to any elliptic PDE that can be formulated in
56 : * first-order flux-form, as detailed by
57 : * `elliptic::protocols::FirstOrderSystem`.
58 : * The DG discretization of equations in this first-order form amounts to
59 : * projecting the equations on the set of basis functions that we also use to
60 : * represent the fields on the computational grid. The currently implemented DG
61 : * operator uses Lagrange interpolating polynomials w.r.t. Legendre-Gauss or
62 : * Legendre-Gauss-Lobatto collocation points as basis functions. Skipping all
63 : * further details here, the discretization results in a linear equation
64 : * \f$A(u)=b\f$ over all grid points and primal variables. Solving the elliptic
65 : * equations amounts to numerically inverting the DG operator \f$A\f$, typically
66 : * without ever constructing the full matrix but by employing an iterative
67 : * linear solver that repeatedly applies the DG operator to "test data". Note
68 : * that the DG operator applies directly to the primal variables. Auxiliary
69 : * variables are only computed temporarily and don't inflate the size of the
70 : * operator. This means the DG operator essentially computes second derivatives
71 : * of the primal variables, modified by the fluxes and sources of the system
72 : * as well as by DG boundary corrections that couple grid points across element
73 : * boundaries.
74 : *
75 : * \par Boundary corrections:
76 : * In this implementation we employ the "internal penalty" DG scheme that
77 : * couples grid points across nearest-neighbor elements through the fluxes:
78 : *
79 : * \f{align}
80 : * \label{eq:internal_penalty_auxiliary}
81 : * u^* &= \frac{1}{2} \left(u^\mathrm{int} + u^\mathrm{ext}\right) \\
82 : * \label{eq:internal_penalty_primal}
83 : * (n_i F^i)^* &= \frac{1}{2} n_i \left(
84 : * F^i_\mathrm{int} + F^i_\mathrm{ext} \right)
85 : * - \sigma n_i F^i(n_j (u^\mathrm{int} - u^\mathrm{ext}))
86 : * \f}
87 : *
88 : * Note that \f$n_i\f$ denotes the face normal on the "interior" side of the
89 : * element under consideration. We assume \f$n^\mathrm{ext}_i=-n_i\f$ in the
90 : * implementation, i.e. face normals don't depend on the dynamic variables
91 : * (which may be discontinuous on element faces). This is the case for the
92 : * problems we are expecting to solve, because those will be on fixed background
93 : * metrics (e.g. a conformal metric for the XCTS system). Numerically, the face
94 : * normals on either side of a mortar may nonetheless be different because the
95 : * two faces adjacent to the mortar may resolve them at different resolutions.
96 : *
97 : * Also note that the numerical fluxes intentionally don't depend on the
98 : * auxiliary field values \f$v\f$. This property allows us to communicate data
99 : * for both the primal and auxiliary boundary corrections together, instead of
100 : * communicating them in two steps. If we were to resort to a two-step
101 : * communication we could replace the derivatives in \f$(n_i F^i)^*\f$ with
102 : * \f$v\f$, which would result in a generalized "stabilized central flux" that
103 : * is slightly less sparse than the internal penalty flux (see e.g.
104 : * \cite HesthavenWarburton, section 7.2). We could also choose to ignore the
105 : * fluxes in the penalty term, but preliminary tests suggest that this may hurt
106 : * convergence.
107 : *
108 : * For a Poisson system (see `Poisson::FirstOrderSystem`) this numerical flux
109 : * reduces to the standard internal penalty flux (see e.g.
110 : * \cite HesthavenWarburton, section 7.2, or \cite Arnold2002):
111 : *
112 : * \f{align}
113 : * u^* &= \frac{1}{2} \left(u^\mathrm{int} + u^\mathrm{ext}\right) \\
114 : * (n_i F^i)^* &= n_i v_i^* = \frac{1}{2} n_i \left(
115 : * \partial_i u^\mathrm{int} + \partial_i u^\mathrm{ext}\right)
116 : * - \sigma \left(u^\mathrm{int} - u^\mathrm{ext}\right)
117 : * \f}
118 : *
119 : * where a sum over repeated indices is assumed, since the equation is
120 : * formulated on a Euclidean geometry.
121 : *
122 : * The penalty factor \f$\sigma\f$ is responsible for removing zero eigenmodes
123 : * and impacts the conditioning of the linear operator to be solved. See
124 : * `elliptic::dg::penalty` for details. For the element size that goes into
125 : * computing the penalty we choose
126 : * \f$h=\frac{J_\mathrm{volume}}{J_\mathrm{face}}\f$, i.e. the ratio of Jacobi
127 : * determinants from logical to inertial coordinates in the element volume and
128 : * on the element face, both evaluated on the face (see \cite Vincent2019qpd).
129 : * Since both \f$N_\mathrm{points}\f$ and \f$h\f$ can be different on either
130 : * side of the element boundary we take the maximum of \f$N_\mathrm{points}\f$
131 : * and the pointwise minimum of \f$h\f$ across the element boundary as is done
132 : * in \cite Vincent2019qpd. Note that we use the number of points
133 : * \f$N_\mathrm{points}\f$ where \cite Vincent2019qpd uses the polynomial degree
134 : * \f$N_\mathrm{points} - 1\f$ because we found unstable configurations on
135 : * curved meshes when using the polynomial degree. Optimizing the penalty on
136 : * curved meshes is subject to further investigation.
137 : *
138 : * \par Discontinuous fluxes:
139 : * The DG operator also supports systems with potentially discontinuous fluxes,
140 : * such as elasticity with layered materials. The way to handle the
141 : * discontinuous fluxes in the DG scheme is described in \cite Vu2023thn.
142 : * Essentially, we evaluate the penalty term in
143 : * Eq. $\ref{eq:internal_penalty_primal}$ on both sides of an element boundary
144 : * and take the average. The other terms in the numerical flux remain unchanged.
145 : */
146 : namespace elliptic::dg {
147 :
148 : /// Data that is projected to mortars and communicated across element
149 : /// boundaries
150 : template <typename PrimalFields, typename PrimalFluxes>
151 1 : using BoundaryData = ::dg::SimpleBoundaryData<
152 : tmpl::append<PrimalFields,
153 : db::wrap_tags_in<::Tags::NormalDotFlux, PrimalFields>>,
154 : tmpl::list<>>;
155 :
156 : /// Boundary data on both sides of a mortar.
157 : ///
158 : /// \note This is a struct (as opposed to a type alias) so it can be used to
159 : /// deduce the template parameters
160 : template <typename TemporalId, typename PrimalFields, typename PrimalFluxes>
161 1 : struct MortarData
162 : : ::dg::SimpleMortarData<TemporalId,
163 : BoundaryData<PrimalFields, PrimalFluxes>,
164 : BoundaryData<PrimalFields, PrimalFluxes>> {};
165 :
166 1 : namespace Tags {
167 : /// Holds `elliptic::dg::MortarData`, i.e. boundary data on both sides of a
168 : /// mortar
169 : template <typename TemporalId, typename PrimalFields, typename PrimalFluxes>
170 1 : struct MortarData : db::SimpleTag {
171 0 : using type = elliptic::dg::MortarData<TemporalId, PrimalFields, PrimalFluxes>;
172 : };
173 : } // namespace Tags
174 :
175 : namespace detail {
176 :
177 : // Mass-conservative restriction: R = M^{-1}_face P^T M_mortar
178 : //
179 : // Note that projecting the mortar data times the Jacobian using
180 : // `Spectral::projection_matrix_child_to_parent(operand_is_massive=false)` is
181 : // equivalent to this implementation on Gauss grids. However, on Gauss-Lobatto
182 : // grids we usually diagonally approximate the mass matrix ("mass lumping") but
183 : // `projection_matrix_child_to_parent(operand_is_massive=false)` uses the full
184 : // mass matrix. Therefore, the two implementations differ slightly on
185 : // Gauss-Lobatto grids. Furthermore, note that
186 : // `projection_matrix_child_to_parent(operand_is_massive=false)` already
187 : // includes the factors of two that account for the mortar size, so they must be
188 : // omitted from the mortar Jacobian when using that approach.
189 : template <typename TagsList, size_t FaceDim>
190 : Variables<TagsList> mass_conservative_restriction(
191 : Variables<TagsList> mortar_vars,
192 : [[maybe_unused]] const Mesh<FaceDim>& mortar_mesh,
193 : [[maybe_unused]] const ::dg::MortarSize<FaceDim>& mortar_size,
194 : [[maybe_unused]] const Scalar<DataVector>& mortar_jacobian,
195 : [[maybe_unused]] const Mesh<FaceDim> face_mesh,
196 : [[maybe_unused]] const Scalar<DataVector>& face_jacobian) {
197 : if constexpr (FaceDim == 0) {
198 : return mortar_vars;
199 : } else {
200 : const auto projection_matrices =
201 : Spectral::projection_matrix_child_to_parent(mortar_mesh, face_mesh,
202 : mortar_size, true);
203 : mortar_vars *= get(mortar_jacobian);
204 : ::dg::apply_mass_matrix(make_not_null(&mortar_vars), mortar_mesh);
205 : auto face_vars =
206 : apply_matrices(projection_matrices, mortar_vars, mortar_mesh.extents());
207 : face_vars /= get(face_jacobian);
208 : ::dg::apply_inverse_mass_matrix(make_not_null(&face_vars), face_mesh);
209 : return face_vars;
210 : }
211 : }
212 :
213 : template <typename System, bool Linearized,
214 : typename PrimalFields = typename System::primal_fields,
215 : typename PrimalFluxes = typename System::primal_fluxes>
216 : struct DgOperatorImpl;
217 :
218 : template <typename System, bool Linearized, typename... PrimalFields,
219 : typename... PrimalFluxes>
220 : struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
221 : tmpl::list<PrimalFluxes...>> {
222 : static_assert(
223 : tt::assert_conforms_to_v<System, elliptic::protocols::FirstOrderSystem>);
224 :
225 : static constexpr size_t Dim = System::volume_dim;
226 : using FluxesComputer = elliptic::get_fluxes_computer<System, Linearized>;
227 : using SourcesComputer = elliptic::get_sources_computer<System, Linearized>;
228 :
229 : struct AllDirections {
230 : bool operator()(const Direction<Dim>& /*unused*/) const { return true; }
231 : };
232 :
233 : struct NoDataIsZero {
234 : bool operator()(const ElementId<Dim>& /*unused*/) const { return false; }
235 : };
236 :
237 : static constexpr auto full_mortar_size =
238 : make_array<Dim - 1>(Spectral::SegmentSize::Full);
239 :
240 : template <bool AllDataIsZero, typename... DerivVars, typename... PrimalVars,
241 : typename... PrimalFluxesVars, typename... PrimalMortarVars,
242 : typename... PrimalMortarFluxes, typename TemporalId,
243 : typename ApplyBoundaryCondition, typename... FluxesArgs,
244 : typename DataIsZero = NoDataIsZero,
245 : typename DirectionsPredicate = AllDirections>
246 : static void prepare_mortar_data(
247 : const gsl::not_null<Variables<tmpl::list<DerivVars...>>*> deriv_vars,
248 : const gsl::not_null<Variables<tmpl::list<PrimalFluxesVars...>>*>
249 : primal_fluxes,
250 : const gsl::not_null<::dg::MortarMap<
251 : Dim, MortarData<TemporalId, tmpl::list<PrimalMortarVars...>,
252 : tmpl::list<PrimalMortarFluxes...>>>*>
253 : all_mortar_data,
254 : const Variables<tmpl::list<PrimalVars...>>& primal_vars,
255 : const Element<Dim>& element, const Mesh<Dim>& mesh,
256 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
257 : Frame::Inertial>& inv_jacobian,
258 : const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
259 : const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
260 : const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
261 : const TemporalId& temporal_id,
262 : const ApplyBoundaryCondition& apply_boundary_condition,
263 : const std::tuple<FluxesArgs...>& fluxes_args,
264 : const DataIsZero& data_is_zero = NoDataIsZero{},
265 : const DirectionsPredicate& directions_predicate = AllDirections{}) {
266 : static_assert(
267 : sizeof...(PrimalVars) == sizeof...(PrimalFields) and
268 : sizeof...(PrimalFluxesVars) == sizeof...(PrimalFluxes),
269 : "The number of variables must match the number of system fields.");
270 : static_assert(
271 : (std::is_same_v<typename PrimalVars::type,
272 : typename PrimalFields::type> and
273 : ...) and
274 : (std::is_same_v<typename PrimalFluxesVars::type,
275 : typename PrimalFluxes::type> and
276 : ...),
277 : "The variables must have the same tensor types as the system fields.");
278 : #ifdef SPECTRE_DEBUG
279 : for (size_t d = 0; d < Dim; ++d) {
280 : ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
281 : (mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto or
282 : mesh.quadrature(d) == Spectral::Quadrature::Gauss),
283 : "The elliptic DG operator is currently only implemented for "
284 : "Legendre-Gauss(-Lobatto) grids. Found basis '"
285 : << mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
286 : << "' in dimension " << d << ".");
287 : }
288 : #endif // SPECTRE_DEBUG
289 : const auto& element_id = element.id();
290 : const bool local_data_is_zero = data_is_zero(element_id);
291 : ASSERT(Linearized or not local_data_is_zero,
292 : "Only a linear operator can take advantage of the knowledge that "
293 : "the operand is zero. Don't return 'true' in 'data_is_zero' unless "
294 : "you also set 'Linearized' to 'true'.");
295 : const size_t num_points = mesh.number_of_grid_points();
296 :
297 : // This function and the one below allocate various Variables to compute
298 : // intermediate quantities. It could be a performance optimization to reduce
299 : // the number of these allocations and/or move some of the memory buffers
300 : // into the DataBox to keep them around permanently. The latter should be
301 : // informed by profiling.
302 :
303 : // Compute partial derivatives grad(u) of the system variables, and from
304 : // those the fluxes F^i(grad(u)) in the volume. We will take the divergence
305 : // of the fluxes in the `apply_operator` function below to compute the full
306 : // elliptic equation -div(F) + S = f(x).
307 : if (AllDataIsZero or local_data_is_zero) {
308 : deriv_vars->initialize(num_points, 0.);
309 : primal_fluxes->initialize(num_points, 0.);
310 : } else {
311 : // Compute partial derivatives of the variables
312 : partial_derivatives(deriv_vars, primal_vars, mesh, inv_jacobian);
313 : // Compute the fluxes
314 : primal_fluxes->initialize(num_points);
315 : std::apply(
316 : [&primal_fluxes, &primal_vars, &deriv_vars,
317 : &element_id](const auto&... expanded_fluxes_args) {
318 : if constexpr (FluxesComputer::is_discontinuous) {
319 : FluxesComputer::apply(
320 : make_not_null(&get<PrimalFluxesVars>(*primal_fluxes))...,
321 : expanded_fluxes_args..., element_id,
322 : get<PrimalVars>(primal_vars)...,
323 : get<DerivVars>(*deriv_vars)...);
324 : } else {
325 : (void)element_id;
326 : FluxesComputer::apply(
327 : make_not_null(&get<PrimalFluxesVars>(*primal_fluxes))...,
328 : expanded_fluxes_args..., get<PrimalVars>(primal_vars)...,
329 : get<DerivVars>(*deriv_vars)...);
330 : }
331 : },
332 : fluxes_args);
333 : }
334 :
335 : // Populate the mortar data on this element's side of the boundary so it's
336 : // ready to be sent to neighbors.
337 : for (const auto& direction : [&element]() -> const auto& {
338 : if constexpr (AllDataIsZero) {
339 : // Skipping internal boundaries for all-zero data because they
340 : // won't contribute boundary corrections anyway (data on both sides
341 : // of the boundary is the same). For all-zero data we are
342 : // interested in external boundaries, to extract inhomogeneous
343 : // boundary conditions from a non-linear operator.
344 : return element.external_boundaries();
345 : } else {
346 : (void)element;
347 : return Direction<Dim>::all_directions();
348 : };
349 : }()) {
350 : if (not directions_predicate(direction)) {
351 : continue;
352 : }
353 : const bool is_internal = element.neighbors().contains(direction);
354 : // Skip directions altogether when both this element and all neighbors in
355 : // the direction have zero data. These boundaries won't contribute
356 : // corrections, because the data is the same on both sides. External
357 : // boundaries also count as zero data here, because they are linearized
358 : // (see assert above).
359 : if (local_data_is_zero and
360 : (not is_internal or
361 : alg::all_of(element.neighbors().at(direction), data_is_zero))) {
362 : continue;
363 : }
364 : const auto face_mesh = mesh.slice_away(direction.dimension());
365 : const size_t face_num_points = face_mesh.number_of_grid_points();
366 : const auto& face_normal = face_normals.at(direction);
367 : Variables<tmpl::list<PrimalFluxesVars...>> primal_fluxes_on_face{};
368 : BoundaryData<tmpl::list<PrimalMortarVars...>,
369 : tmpl::list<PrimalMortarFluxes...>>
370 : boundary_data{};
371 : if (AllDataIsZero or local_data_is_zero) {
372 : if (is_internal) {
373 : // We manufacture zero boundary data directly on the mortars below.
374 : // Nothing to do here.
375 : } else {
376 : boundary_data.field_data.initialize(face_num_points, 0.);
377 : }
378 : } else {
379 : boundary_data.field_data.initialize(face_num_points);
380 : primal_fluxes_on_face.initialize(face_num_points);
381 : // Project fields to faces
382 : // Note: need to convert tags of `Variables` because
383 : // `project_contiguous_data_to_boundary` requires that the face and
384 : // volume tags are subsets.
385 : Variables<
386 : tmpl::list<PrimalVars..., ::Tags::NormalDotFlux<PrimalVars>...>>
387 : boundary_data_ref{};
388 : boundary_data_ref.set_data_ref(boundary_data.field_data.data(),
389 : boundary_data.field_data.size());
390 : ::dg::project_contiguous_data_to_boundary(
391 : make_not_null(&boundary_data_ref), primal_vars, mesh, direction);
392 : // Compute n_i F^i on faces
393 : ::dg::project_contiguous_data_to_boundary(
394 : make_not_null(&primal_fluxes_on_face), *primal_fluxes, mesh,
395 : direction);
396 : EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
397 : make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
398 : boundary_data.field_data)),
399 : face_normal, get<PrimalFluxesVars>(primal_fluxes_on_face)));
400 : }
401 :
402 : if (is_internal) {
403 : if constexpr (not AllDataIsZero) {
404 : // Project boundary data on internal faces to mortars
405 : for (const auto& neighbor_id : element.neighbors().at(direction)) {
406 : if (local_data_is_zero and data_is_zero(neighbor_id)) {
407 : continue;
408 : }
409 : const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
410 : const auto& mortar_mesh = all_mortar_meshes.at(mortar_id);
411 : const auto& mortar_size = all_mortar_sizes.at(mortar_id);
412 : if (local_data_is_zero) {
413 : // No need to project anything. We just manufacture zero boundary
414 : // data on the mortar.
415 : BoundaryData<tmpl::list<PrimalMortarVars...>,
416 : tmpl::list<PrimalMortarFluxes...>>
417 : zero_boundary_data{};
418 : zero_boundary_data.field_data.initialize(
419 : mortar_mesh.number_of_grid_points(), 0.);
420 : (*all_mortar_data)[mortar_id].local_insert(
421 : temporal_id, std::move(zero_boundary_data));
422 : continue;
423 : }
424 : // When no projection is necessary we can safely move the boundary
425 : // data from the face as there is only a single neighbor in this
426 : // direction
427 : auto projected_boundary_data =
428 : Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)
429 : // NOLINTNEXTLINE
430 : ? boundary_data.project_to_mortar(face_mesh, mortar_mesh,
431 : mortar_size)
432 : : std::move(boundary_data); // NOLINT
433 : (*all_mortar_data)[mortar_id].local_insert(
434 : temporal_id, std::move(projected_boundary_data));
435 : }
436 : }
437 : } else {
438 : // No need to do projections on external boundaries
439 : const ::dg::MortarId<Dim> mortar_id{
440 : direction, ElementId<Dim>::external_boundary_id()};
441 : (*all_mortar_data)[mortar_id].local_insert(temporal_id, boundary_data);
442 :
443 : // -------------------------
444 : // Apply boundary conditions
445 : // -------------------------
446 : //
447 : // To apply boundary conditions we fill the boundary data with
448 : // "exterior" or "ghost" data and set it as remote mortar data, so
449 : // external boundaries behave just like internal boundaries when
450 : // applying boundary corrections.
451 : //
452 : // The `apply_boundary_conditions` invocable is expected to impose
453 : // boundary conditions by modifying the fields and fluxes that are
454 : // passed by reference. Dirichlet-type boundary conditions are imposed
455 : // by modifying the fields, and Neumann-type boundary conditions are
456 : // imposed by modifying the interior n.F. Note that all data passed to
457 : // the boundary conditions is taken from the "interior" side of the
458 : // boundary, i.e. with a normal vector that points _out_ of the
459 : // computational domain.
460 : Variables<tmpl::list<DerivVars...>> deriv_vars_on_boundary{};
461 : if (AllDataIsZero or local_data_is_zero) {
462 : deriv_vars_on_boundary.initialize(face_num_points, 0.);
463 : } else {
464 : deriv_vars_on_boundary.initialize(face_num_points);
465 : ::dg::project_contiguous_data_to_boundary(
466 : make_not_null(&deriv_vars_on_boundary), *deriv_vars, mesh,
467 : direction);
468 : }
469 : apply_boundary_condition(
470 : direction,
471 : make_not_null(&get<PrimalMortarVars>(boundary_data.field_data))...,
472 : make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
473 : boundary_data.field_data))...,
474 : get<DerivVars>(deriv_vars_on_boundary)...);
475 :
476 : // Invert the sign of the fluxes to account for the inverted normal on
477 : // exterior faces. Also multiply by 2 and add the interior fluxes to
478 : // impose the boundary conditions on the _average_ instead of just
479 : // setting the fields on the exterior:
480 : // (Dirichlet) u_D = avg(u) = 1/2 (u_int + u_ext)
481 : // => u_ext = 2 u_D - u_int
482 : // (Neumann) (n.F)_N = avg(n.F) = 1/2 [(n.F)_int - (n.F)_ext]
483 : // => (n.F)_ext = -2 (n.F)_N + (n.F)_int]
484 : const auto impose_on_average = [](const auto exterior_field,
485 : const auto& interior_field) {
486 : for (size_t i = 0; i < interior_field.size(); ++i) {
487 : (*exterior_field)[i] *= 2.;
488 : (*exterior_field)[i] -= interior_field[i];
489 : }
490 : };
491 : EXPAND_PACK_LEFT_TO_RIGHT(impose_on_average(
492 : make_not_null(&get<PrimalMortarVars>(boundary_data.field_data)),
493 : get<PrimalMortarVars>(all_mortar_data->at(mortar_id)
494 : .local_data(temporal_id)
495 : .field_data)));
496 : const auto invert_sign_and_impose_on_average =
497 : [](const auto exterior_n_dot_flux,
498 : const auto& interior_n_dot_flux) {
499 : for (size_t i = 0; i < interior_n_dot_flux.size(); ++i) {
500 : (*exterior_n_dot_flux)[i] *= -2.;
501 : (*exterior_n_dot_flux)[i] += interior_n_dot_flux[i];
502 : }
503 : };
504 : EXPAND_PACK_LEFT_TO_RIGHT(invert_sign_and_impose_on_average(
505 : make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
506 : boundary_data.field_data)),
507 : get<::Tags::NormalDotFlux<PrimalMortarVars>>(
508 : all_mortar_data->at(mortar_id)
509 : .local_data(temporal_id)
510 : .field_data)));
511 :
512 : // Store the exterior boundary data on the mortar
513 : all_mortar_data->at(mortar_id).remote_insert(temporal_id,
514 : std::move(boundary_data));
515 : } // if (is_internal)
516 : } // loop directions
517 : }
518 :
519 : // --- This is essentially a break to communicate the mortar data ---
520 :
521 : template <typename... OperatorTags, typename... PrimalVars,
522 : typename... DerivVars, typename... PrimalFluxesVars,
523 : typename... PrimalMortarVars, typename... PrimalMortarFluxes,
524 : typename TemporalId, typename... FluxesArgs,
525 : typename... SourcesArgs, typename... ModifyBoundaryDataArgs,
526 : typename DataIsZero = NoDataIsZero,
527 : typename DirectionsPredicate = AllDirections>
528 : static void apply_operator(
529 : const gsl::not_null<Variables<tmpl::list<OperatorTags...>>*>
530 : operator_applied_to_vars,
531 : const gsl::not_null<::dg::MortarMap<
532 : Dim, MortarData<TemporalId, tmpl::list<PrimalMortarVars...>,
533 : tmpl::list<PrimalMortarFluxes...>>>*>
534 : all_mortar_data,
535 : const Variables<tmpl::list<PrimalVars...>>& primal_vars,
536 : const Variables<tmpl::list<DerivVars...>>& deriv_vars,
537 : // Taking the primal fluxes computed in the `prepare_mortar_data` function
538 : // by const-ref here because other code might use them and so we don't
539 : // want to modify them by adding boundary corrections. E.g. linearized
540 : // sources use the nonlinear fields and fluxes as background fields.
541 : const Variables<tmpl::list<PrimalFluxesVars...>>& primal_fluxes,
542 : const Element<Dim>& element, const Mesh<Dim>& mesh,
543 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
544 : Frame::Inertial>& inv_jacobian,
545 : const Scalar<DataVector>& det_inv_jacobian,
546 : const Scalar<DataVector>& det_jacobian,
547 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
548 : Frame::Inertial>& det_times_inv_jacobian,
549 : const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
550 : const DirectionMap<Dim, tnsr::I<DataVector, Dim>>& face_normal_vectors,
551 : const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
552 : const DirectionMap<Dim, Scalar<DataVector>>& face_jacobians,
553 : const DirectionMap<Dim,
554 : InverseJacobian<DataVector, Dim, Frame::ElementLogical,
555 : Frame::Inertial>>&
556 : face_jacobian_times_inv_jacobians,
557 : const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
558 : const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
559 : const ::dg::MortarMap<Dim, Scalar<DataVector>>& mortar_jacobians,
560 : const ::dg::MortarMap<Dim, Scalar<DataVector>>& penalty_factors,
561 : const bool massive, const ::dg::Formulation formulation,
562 : const TemporalId& /*temporal_id*/,
563 : const DirectionMap<Dim, std::tuple<FluxesArgs...>>& fluxes_args_on_faces,
564 : const std::tuple<SourcesArgs...>& sources_args,
565 : const std::tuple<ModifyBoundaryDataArgs...>& modify_boundary_data_args,
566 : const DataIsZero& data_is_zero = NoDataIsZero{},
567 : const DirectionsPredicate& directions_predicate = AllDirections{}) {
568 : static_assert(
569 : sizeof...(PrimalVars) == sizeof...(PrimalFields) and
570 : sizeof...(PrimalFluxesVars) == sizeof...(PrimalFluxes) and
571 : sizeof...(PrimalMortarVars) == sizeof...(PrimalFields) and
572 : sizeof...(PrimalMortarFluxes) == sizeof...(PrimalFluxes) and
573 : sizeof...(OperatorTags) == sizeof...(PrimalFields),
574 : "The number of variables must match the number of system fields.");
575 : static_assert(
576 : (std::is_same_v<typename PrimalVars::type,
577 : typename PrimalFields::type> and
578 : ...) and
579 : (std::is_same_v<typename PrimalFluxesVars::type,
580 : typename PrimalFluxes::type> and
581 : ...) and
582 : (std::is_same_v<typename PrimalMortarVars::type,
583 : typename PrimalFields::type> and
584 : ...) and
585 : (std::is_same_v<typename PrimalMortarFluxes::type,
586 : typename PrimalFluxes::type> and
587 : ...) and
588 : (std::is_same_v<typename OperatorTags::type,
589 : typename PrimalFields::type> and
590 : ...),
591 : "The variables must have the same tensor types as the system fields.");
592 : #ifdef SPECTRE_DEBUG
593 : for (size_t d = 0; d < Dim; ++d) {
594 : ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
595 : (mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto or
596 : mesh.quadrature(d) == Spectral::Quadrature::Gauss),
597 : "The elliptic DG operator is currently only implemented for "
598 : "Legendre-Gauss(-Lobatto) grids. Found basis '"
599 : << mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
600 : << "' in dimension " << d << ".");
601 : }
602 : #endif // SPECTRE_DEBUG
603 : const auto& element_id = element.id();
604 : const bool local_data_is_zero = data_is_zero(element_id);
605 : ASSERT(Linearized or not local_data_is_zero,
606 : "Only a linear operator can take advantage of the knowledge that "
607 : "the operand is zero. Don't return 'true' in 'data_is_zero' unless "
608 : "you also set 'Linearized' to 'true'.");
609 : const size_t num_points = mesh.number_of_grid_points();
610 :
611 : // This function and the one above allocate various Variables to compute
612 : // intermediate quantities. It could be a performance optimization to reduce
613 : // the number of these allocations and/or move some of the memory buffers
614 : // into the DataBox to keep them around permanently. The latter should be
615 : // informed by profiling.
616 :
617 : // Compute volume terms: -div(F) + S
618 : if (local_data_is_zero) {
619 : operator_applied_to_vars->initialize(num_points, 0.);
620 : } else {
621 : // "Massive" operators retain the factors from the volume integral:
622 : // \int_volume div(F) \phi_p = w_p det(J)_p div(F)_p
623 : // Here, `w` are the quadrature weights (the diagonal logical mass matrix
624 : // with mass-lumping) and det(J) is the Jacobian determinant. The
625 : // quantities are evaluated at the grid point `p`.
626 : if (formulation == ::dg::Formulation::StrongInertial) {
627 : // Compute strong divergence:
628 : // div(F) = (J^\hat{i}_i)_p \sum_q (D_\hat{i})_pq (F^i)_q.
629 : divergence(operator_applied_to_vars, primal_fluxes, mesh,
630 : massive ? det_times_inv_jacobian : inv_jacobian);
631 : // This is the sign flip that makes the operator _minus_ the Laplacian
632 : // for a Poisson system
633 : *operator_applied_to_vars *= -1.;
634 : } else if (formulation == ::dg::Formulation::StrongLogical) {
635 : // Strong divergence but with the Jacobian moved into the divergence:
636 : // div(F) = 1/J_p \sum_q (D_\hat{i})_pq J_q (J^\hat{i}_i)_q (F^i)_q.
637 : const auto logical_fluxes = transform::first_index_to_different_frame(
638 : primal_fluxes, det_times_inv_jacobian);
639 : logical_divergence(operator_applied_to_vars, logical_fluxes, mesh);
640 : if (massive) {
641 : *operator_applied_to_vars *= -1.;
642 : } else {
643 : *operator_applied_to_vars *= -1. * get(det_inv_jacobian);
644 : }
645 : } else if (formulation == ::dg::Formulation::WeakInertial) {
646 : // Compute weak divergence:
647 : // F^i \partial_i \phi = 1/w_p \sum_q
648 : // (D^T_\hat{i})_pq (w det(J) J^\hat{i}_i F^i)_q
649 : weak_divergence(operator_applied_to_vars, primal_fluxes, mesh,
650 : det_times_inv_jacobian);
651 : if (not massive) {
652 : *operator_applied_to_vars *= get(det_inv_jacobian);
653 : }
654 : } else {
655 : ERROR("Unsupported DG formulation: "
656 : << formulation
657 : << "\nSupported formulations are: StrongInertial, WeakInertial, "
658 : "StrongLogical.");
659 : }
660 : if constexpr (not std::is_same_v<SourcesComputer, void>) {
661 : Variables<tmpl::list<OperatorTags...>> sources{num_points, 0.};
662 : std::apply(
663 : [&sources, &primal_vars, &deriv_vars,
664 : &primal_fluxes](const auto&... expanded_sources_args) {
665 : SourcesComputer::apply(
666 : make_not_null(&get<OperatorTags>(sources))...,
667 : expanded_sources_args..., get<PrimalVars>(primal_vars)...,
668 : get<DerivVars>(deriv_vars)...,
669 : get<PrimalFluxesVars>(primal_fluxes)...);
670 : },
671 : sources_args);
672 : if (massive) {
673 : sources *= get(det_jacobian);
674 : }
675 : *operator_applied_to_vars += sources;
676 : }
677 : }
678 : if (massive) {
679 : ::dg::apply_mass_matrix(operator_applied_to_vars, mesh);
680 : }
681 :
682 : // Add boundary corrections
683 : // Keeping track if any corrections were applied here, for an optimization
684 : // below
685 : bool has_any_boundary_corrections = false;
686 : Variables<tmpl::list<transform::Tags::TransformedFirstIndex<
687 : PrimalFluxesVars, Frame::ElementLogical>...>>
688 : lifted_logical_aux_boundary_corrections{num_points, 0.};
689 : for (auto& [mortar_id, mortar_data] : *all_mortar_data) {
690 : const auto& direction = mortar_id.direction();
691 : const auto& neighbor_id = mortar_id.id();
692 : const bool is_internal =
693 : (neighbor_id != ElementId<Dim>::external_boundary_id());
694 : if (not directions_predicate(direction)) {
695 : continue;
696 : }
697 : // When the data on both sides of the mortar is zero then we don't need to
698 : // handle this mortar at all.
699 : if (local_data_is_zero and
700 : (not is_internal or data_is_zero(neighbor_id))) {
701 : continue;
702 : }
703 : has_any_boundary_corrections = true;
704 :
705 : const auto face_mesh = mesh.slice_away(direction.dimension());
706 : auto [local_data, remote_data] = mortar_data.extract();
707 :
708 : if (is_internal) {
709 : if constexpr (Linearized and
710 : not std::is_same_v<typename System::modify_boundary_data,
711 : void>) {
712 : // Apply a linearized modification to received boundary data.
713 : // This allows modifications to depend linearly on the variables.
714 : std::apply(
715 : [&rem_data = remote_data, &loc_data = local_data,
716 : m_id = mortar_id](const auto&... args) {
717 : System::modify_boundary_data::apply_linearized(
718 : make_not_null(
719 : &get<PrimalMortarVars>(rem_data.field_data))...,
720 : make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
721 : rem_data.field_data))...,
722 : get<PrimalMortarVars>(loc_data.field_data)...,
723 : get<::Tags::NormalDotFlux<PrimalMortarVars>>(
724 : loc_data.field_data)...,
725 : m_id, args...);
726 : },
727 : modify_boundary_data_args);
728 : }
729 : }
730 :
731 : const size_t face_num_points = face_mesh.number_of_grid_points();
732 : const auto& face_normal = face_normals.at(direction);
733 : const auto& face_normal_vector = face_normal_vectors.at(direction);
734 : const auto& fluxes_args_on_face = fluxes_args_on_faces.at(direction);
735 : const auto& face_normal_magnitude = face_normal_magnitudes.at(direction);
736 : const auto& face_jacobian = face_jacobians.at(direction);
737 : const auto& face_jacobian_times_inv_jacobian =
738 : face_jacobian_times_inv_jacobians.at(direction);
739 : const auto& mortar_mesh =
740 : is_internal ? all_mortar_meshes.at(mortar_id) : face_mesh;
741 : const auto& mortar_size =
742 : is_internal ? all_mortar_sizes.at(mortar_id) : full_mortar_size;
743 :
744 : // This is the strong auxiliary boundary correction:
745 : // G^i = F^i(n_j (avg(u) - u))
746 : // where
747 : // avg(u) - u = -0.5 * (u_int - u_ext)
748 : auto avg_vars_on_mortar = Variables<tmpl::list<PrimalMortarVars...>>(
749 : local_data.field_data
750 : .template extract_subset<tmpl::list<PrimalMortarVars...>>());
751 : const auto add_remote_contribution = [](auto& lhs, const auto& rhs) {
752 : for (size_t i = 0; i < lhs.size(); ++i) {
753 : lhs[i] -= rhs[i];
754 : }
755 : };
756 : EXPAND_PACK_LEFT_TO_RIGHT(add_remote_contribution(
757 : get<PrimalMortarVars>(avg_vars_on_mortar),
758 : get<PrimalMortarVars>(remote_data.field_data)));
759 : avg_vars_on_mortar *= -0.5;
760 :
761 : // Project from the mortar back down to the face if needed
762 : const auto avg_vars_on_face =
763 : Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)
764 : ? mass_conservative_restriction(
765 : std::move(avg_vars_on_mortar), mortar_mesh, mortar_size,
766 : mortar_jacobians.at(mortar_id), face_mesh, face_jacobian)
767 : : std::move(avg_vars_on_mortar);
768 :
769 : // Apply fluxes to get G^i
770 : Variables<tmpl::list<PrimalFluxesVars...>> auxiliary_boundary_corrections{
771 : face_num_points};
772 : std::apply(
773 : [&auxiliary_boundary_corrections, &face_normal, &face_normal_vector,
774 : &avg_vars_on_face,
775 : &element_id](const auto&... expanded_fluxes_args_on_face) {
776 : if constexpr (FluxesComputer::is_discontinuous) {
777 : FluxesComputer::apply(make_not_null(&get<PrimalFluxesVars>(
778 : auxiliary_boundary_corrections))...,
779 : expanded_fluxes_args_on_face..., element_id,
780 : face_normal, face_normal_vector,
781 : get<PrimalMortarVars>(avg_vars_on_face)...);
782 : } else {
783 : (void)element_id;
784 : FluxesComputer::apply(make_not_null(&get<PrimalFluxesVars>(
785 : auxiliary_boundary_corrections))...,
786 : expanded_fluxes_args_on_face...,
787 : face_normal, face_normal_vector,
788 : get<PrimalMortarVars>(avg_vars_on_face)...);
789 : }
790 : },
791 : fluxes_args_on_face);
792 :
793 : // Lifting for the auxiliary boundary correction:
794 : // \int_face G^i \partial_i \phi
795 : // We first transform the flux index to the logical frame, apply the
796 : // quadrature weights and the Jacobian for the face integral, then take
797 : // the logical weak divergence in the volume after lifting (below the loop
798 : // over faces).
799 : auto logical_aux_boundary_corrections =
800 : transform::first_index_to_different_frame(
801 : auxiliary_boundary_corrections, face_jacobian_times_inv_jacobian);
802 : ::dg::apply_mass_matrix(make_not_null(&logical_aux_boundary_corrections),
803 : face_mesh);
804 : if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
805 : add_slice_to_data(
806 : make_not_null(&lifted_logical_aux_boundary_corrections),
807 : logical_aux_boundary_corrections, mesh.extents(),
808 : direction.dimension(),
809 : index_to_slice_at(mesh.extents(), direction));
810 : } else {
811 : ::dg::lift_boundary_terms_gauss_points(
812 : make_not_null(&lifted_logical_aux_boundary_corrections),
813 : logical_aux_boundary_corrections, mesh, direction);
814 : }
815 :
816 : // This is the strong primal boundary correction:
817 : // -n.H = -avg(n.F) + n.F + penalty * n.F(n_j jump(u))
818 : // Note that the "internal penalty" numerical flux
819 : // (as opposed to the LLF flux) uses the raw field derivatives without
820 : // boundary corrections in the average, which is why we can communicate
821 : // the data so early together with the auxiliary boundary data. In this
822 : // case the penalty needs to include a factor N_points^2 / h (see the
823 : // `penalty` function).
824 : const auto& penalty_factor = penalty_factors.at(mortar_id);
825 : // Compute jump on mortar:
826 : // penalty * jump(u) = penalty * (u_int - u_ext)
827 : const auto add_remote_jump_contribution =
828 : [&penalty_factor](auto& lhs, const auto& rhs) {
829 : for (size_t i = 0; i < lhs.size(); ++i) {
830 : lhs[i] -= rhs[i];
831 : lhs[i] *= get(penalty_factor);
832 : }
833 : };
834 : EXPAND_PACK_LEFT_TO_RIGHT(add_remote_jump_contribution(
835 : get<PrimalMortarVars>(local_data.field_data),
836 : get<PrimalMortarVars>(remote_data.field_data)));
837 : // Compute average on mortar:
838 : // (strong) -avg(n.F) + n.F = 0.5 * (n.F)_int + 0.5 * (n.F)_ext
839 : // (weak) -avg(n.F) = -0.5 * (n.F)_int + 0.5 * (n.F)_ext
840 : const auto add_avg_contribution = [](auto& lhs, const auto& rhs,
841 : const double factor) {
842 : for (size_t i = 0; i < lhs.size(); ++i) {
843 : lhs[i] *= factor;
844 : lhs[i] += 0.5 * rhs[i];
845 : }
846 : };
847 : const bool is_strong_formulation =
848 : formulation == ::dg::Formulation::StrongInertial or
849 : formulation == ::dg::Formulation::StrongLogical;
850 : EXPAND_PACK_LEFT_TO_RIGHT(add_avg_contribution(
851 : get<::Tags::NormalDotFlux<PrimalMortarVars>>(local_data.field_data),
852 : get<::Tags::NormalDotFlux<PrimalMortarVars>>(remote_data.field_data),
853 : is_strong_formulation ? 0.5 : -0.5));
854 :
855 : // Project from the mortar back down to the face if needed, lift and add
856 : // to operator. See auxiliary boundary corrections above for details.
857 : auto primal_boundary_corrections_on_face =
858 : Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)
859 : ? mass_conservative_restriction(
860 : std::move(local_data.field_data), mortar_mesh, mortar_size,
861 : mortar_jacobians.at(mortar_id), face_mesh, face_jacobian)
862 : : std::move(local_data.field_data);
863 :
864 : // Compute fluxes for jump term: n.F(n_j jump(u))
865 : // If the fluxes are trivial (just the spatial metric), we can skip this
866 : // step because the face normal is normalized.
867 : if constexpr (not FluxesComputer::is_trivial) {
868 : // We reuse the memory buffer from above for the result.
869 : std::apply(
870 : [&auxiliary_boundary_corrections, &face_normal, &face_normal_vector,
871 : &primal_boundary_corrections_on_face,
872 : &element_id](const auto&... expanded_fluxes_args_on_face) {
873 : if constexpr (FluxesComputer::is_discontinuous) {
874 : FluxesComputer::apply(
875 : make_not_null(&get<PrimalFluxesVars>(
876 : auxiliary_boundary_corrections))...,
877 : expanded_fluxes_args_on_face..., element_id, face_normal,
878 : face_normal_vector,
879 : get<PrimalMortarVars>(
880 : primal_boundary_corrections_on_face)...);
881 : } else {
882 : (void)element_id;
883 : FluxesComputer::apply(
884 : make_not_null(&get<PrimalFluxesVars>(
885 : auxiliary_boundary_corrections))...,
886 : expanded_fluxes_args_on_face..., face_normal,
887 : face_normal_vector,
888 : get<PrimalMortarVars>(
889 : primal_boundary_corrections_on_face)...);
890 : }
891 : },
892 : fluxes_args_on_face);
893 : if constexpr (FluxesComputer::is_discontinuous) {
894 : if (is_internal) {
895 : // For penalty term with discontinuous fluxes: evaluate the fluxes
896 : // on the other side of the boundary as well and take average
897 : Variables<tmpl::list<PrimalFluxesVars...>> fluxes_other_side{
898 : face_num_points};
899 : std::apply(
900 : [&fluxes_other_side, &face_normal, &face_normal_vector,
901 : &primal_boundary_corrections_on_face,
902 : &local_neighbor_id =
903 : neighbor_id](const auto&... expanded_fluxes_args_on_face) {
904 : FluxesComputer::apply(
905 : make_not_null(
906 : &get<PrimalFluxesVars>(fluxes_other_side))...,
907 : expanded_fluxes_args_on_face..., local_neighbor_id,
908 : face_normal, face_normal_vector,
909 : get<PrimalMortarVars>(
910 : primal_boundary_corrections_on_face)...);
911 : },
912 : fluxes_args_on_face);
913 : auxiliary_boundary_corrections += fluxes_other_side;
914 : auxiliary_boundary_corrections *= 0.5;
915 : }
916 : }
917 : EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
918 : make_not_null(
919 : &get<PrimalMortarVars>(primal_boundary_corrections_on_face)),
920 : face_normal,
921 : get<PrimalFluxesVars>(auxiliary_boundary_corrections)));
922 : }
923 :
924 : // Add penalty term to average term
925 : Variables<tmpl::list<PrimalMortarVars...>> primal_boundary_corrections{};
926 : // First half of the memory allocated above is filled with the penalty
927 : // term, so just use that memory here.
928 : primal_boundary_corrections.set_data_ref(
929 : primal_boundary_corrections_on_face.data(), avg_vars_on_face.size());
930 : // Second half of the memory is filled with the average term. Add that to
931 : // the penalty term.
932 : const auto add_avg_term = [](auto& lhs, const auto& rhs) {
933 : for (size_t i = 0; i < lhs.size(); ++i) {
934 : lhs[i] += rhs[i];
935 : }
936 : };
937 : EXPAND_PACK_LEFT_TO_RIGHT(
938 : add_avg_term(get<PrimalMortarVars>(primal_boundary_corrections),
939 : get<::Tags::NormalDotFlux<PrimalMortarVars>>(
940 : primal_boundary_corrections_on_face)));
941 :
942 : // Lifting for the primal boundary correction:
943 : // \int_face n.H \phi
944 : if (massive) {
945 : // We apply the quadrature weights and Jacobian for the face integral,
946 : // then lift to the volume
947 : primal_boundary_corrections *= get(face_jacobian);
948 : ::dg::apply_mass_matrix(make_not_null(&primal_boundary_corrections),
949 : face_mesh);
950 : if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
951 : add_slice_to_data(operator_applied_to_vars,
952 : primal_boundary_corrections, mesh.extents(),
953 : direction.dimension(),
954 : index_to_slice_at(mesh.extents(), direction));
955 : } else {
956 : ::dg::lift_boundary_terms_gauss_points(operator_applied_to_vars,
957 : primal_boundary_corrections,
958 : mesh, direction);
959 : }
960 : } else {
961 : // Apply an extra inverse mass matrix to the boundary corrections (with
962 : // mass lumping, so it's diagonal).
963 : // For Gauss-Lobatto grids this divides out the quadrature weights and
964 : // Jacobian on the face, leaving only factors perpendicular to the face.
965 : // Those are handled by `dg::lift_flux` (with an extra minus sign since
966 : // the function was written for evolution systems).
967 : // For Gauss grids the quadrature weights and Jacobians are handled
968 : // by `::dg::lift_boundary_terms_gauss_points` (which was also written
969 : // for evolution systems, hence the extra minus sign).
970 : primal_boundary_corrections *= -1.;
971 : if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
972 : ::dg::lift_flux(make_not_null(&primal_boundary_corrections),
973 : mesh.extents(direction.dimension()),
974 : face_normal_magnitude);
975 : add_slice_to_data(operator_applied_to_vars,
976 : primal_boundary_corrections, mesh.extents(),
977 : direction.dimension(),
978 : index_to_slice_at(mesh.extents(), direction));
979 : } else {
980 : // We already have the `face_jacobian = det(J) * magnitude(n)` here,
981 : // so just pass a constant 1 for `magnitude(n)`. This could be
982 : // optimized to avoid allocating the vector of ones.
983 : ::dg::lift_boundary_terms_gauss_points(
984 : operator_applied_to_vars, det_inv_jacobian, mesh, direction,
985 : primal_boundary_corrections,
986 : Scalar<DataVector>{face_mesh.number_of_grid_points(), 1.},
987 : face_jacobian);
988 : }
989 : }
990 : } // loop over all mortars
991 :
992 : if (not has_any_boundary_corrections) {
993 : // No need to handle auxiliary boundary corrections; return early
994 : return;
995 : }
996 :
997 : // Apply weak divergence to lifted auxiliary boundary corrections and add to
998 : // operator
999 : if (massive) {
1000 : logical_weak_divergence(operator_applied_to_vars,
1001 : lifted_logical_aux_boundary_corrections, mesh,
1002 : true);
1003 : } else {
1004 : // Possible optimization: eliminate this allocation by building the
1005 : // inverse mass matrix into `logical_weak_divergence`
1006 : Variables<tmpl::list<OperatorTags...>> massless_aux_boundary_corrections{
1007 : num_points};
1008 : logical_weak_divergence(make_not_null(&massless_aux_boundary_corrections),
1009 : lifted_logical_aux_boundary_corrections, mesh);
1010 : massless_aux_boundary_corrections *= get(det_inv_jacobian);
1011 : ::dg::apply_inverse_mass_matrix(
1012 : make_not_null(&massless_aux_boundary_corrections), mesh);
1013 : *operator_applied_to_vars += massless_aux_boundary_corrections;
1014 : }
1015 : }
1016 :
1017 : template <typename... FixedSourcesTags, typename ApplyBoundaryCondition,
1018 : typename... FluxesArgs, typename... SourcesArgs,
1019 : typename... ModifyBoundaryDataArgs>
1020 : static void impose_inhomogeneous_boundary_conditions_on_source(
1021 : const gsl::not_null<Variables<tmpl::list<FixedSourcesTags...>>*>
1022 : fixed_sources,
1023 : const Element<Dim>& element, const Mesh<Dim>& mesh,
1024 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
1025 : Frame::Inertial>& inv_jacobian,
1026 : const Scalar<DataVector>& det_inv_jacobian,
1027 : const Scalar<DataVector>& det_jacobian,
1028 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
1029 : Frame::Inertial>& det_times_inv_jacobian,
1030 : const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
1031 : const DirectionMap<Dim, tnsr::I<DataVector, Dim>>& face_normal_vectors,
1032 : const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
1033 : const DirectionMap<Dim, Scalar<DataVector>>& face_jacobians,
1034 : const DirectionMap<Dim,
1035 : InverseJacobian<DataVector, Dim, Frame::ElementLogical,
1036 : Frame::Inertial>>&
1037 : face_jacobian_times_inv_jacobians,
1038 : const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
1039 : const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
1040 : const ::dg::MortarMap<Dim, Scalar<DataVector>>& mortar_jacobians,
1041 : const ::dg::MortarMap<Dim, Scalar<DataVector>>& penalty_factors,
1042 : const bool massive, const ::dg::Formulation formulation,
1043 : const ApplyBoundaryCondition& apply_boundary_condition,
1044 : const std::tuple<FluxesArgs...>& fluxes_args,
1045 : const std::tuple<SourcesArgs...>& sources_args,
1046 : const DirectionMap<Dim, std::tuple<FluxesArgs...>>& fluxes_args_on_faces,
1047 : const std::tuple<ModifyBoundaryDataArgs...>& modify_boundary_data_args)
1048 : // This function adds nothing to the fixed sources if the operator is
1049 : // linearized, so it shouldn't be used in that case
1050 : requires(not Linearized)
1051 : {
1052 : // We just feed zero variables through the nonlinear operator to extract the
1053 : // constant contribution at external boundaries. Since the variables are
1054 : // zero the operator simplifies quite a lot. The simplification is probably
1055 : // not very important for performance because this function will only be
1056 : // called when solving a linear elliptic system and only once during
1057 : // initialization, but we specialize the operator for zero data nonetheless
1058 : // just so we can ignore internal boundaries. Only when the system modifies
1059 : // boundary data do we need to handle internal boundaries (see below),
1060 : // otherwise we can skip them.
1061 : const size_t num_points = mesh.number_of_grid_points();
1062 : const Variables<tmpl::list<PrimalFields...>> zero_primal_vars{num_points,
1063 : 0.};
1064 : Variables<tmpl::list<PrimalFluxes...>> primal_fluxes_buffer{num_points};
1065 : Variables<tmpl::list<
1066 : ::Tags::deriv<PrimalFields, tmpl::size_t<Dim>, Frame::Inertial>...>>
1067 : zero_deriv_vars{num_points, 0.};
1068 : Variables<tmpl::list<FixedSourcesTags...>> operator_applied_to_zero_vars{
1069 : num_points};
1070 : // Set up data on mortars
1071 : ::dg::MortarMap<Dim, MortarData<size_t, tmpl::list<PrimalFields...>,
1072 : tmpl::list<PrimalFluxes...>>>
1073 : all_mortar_data{};
1074 : constexpr size_t temporal_id = std::numeric_limits<size_t>::max();
1075 : // Apply the operator to the zero variables, skipping internal boundaries
1076 : prepare_mortar_data<true>(
1077 : make_not_null(&zero_deriv_vars), make_not_null(&primal_fluxes_buffer),
1078 : make_not_null(&all_mortar_data), zero_primal_vars, element, mesh,
1079 : inv_jacobian, face_normals, all_mortar_meshes, all_mortar_sizes,
1080 : temporal_id, apply_boundary_condition, fluxes_args);
1081 : // Modify internal boundary data if needed, e.g. to transform from one
1082 : // variable to another when crossing element boundaries. This is a nonlinear
1083 : // operation in the sense that feeding zero through the operator is nonzero,
1084 : // so we evaluate it here and add the contribution to the fixed sources,
1085 : // just like inhomogeneous external boundary conditions.
1086 : if constexpr (not std::is_same_v<typename System::modify_boundary_data,
1087 : void>) {
1088 : for (const auto& [direction, neighbors] : element.neighbors()) {
1089 : for (const auto& neighbor_id : neighbors) {
1090 : const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
1091 : ASSERT(not all_mortar_data.contains(mortar_id),
1092 : "Mortar data with ID " << mortar_id << " already exists.");
1093 : auto& mortar_data = all_mortar_data[mortar_id];
1094 : // Manufacture zero mortar data, store as local data, then treat as
1095 : // received from neighbor and apply modifications
1096 : BoundaryData<tmpl::list<PrimalFields...>, tmpl::list<PrimalFluxes...>>
1097 : remote_boundary_data_on_mortar{};
1098 : remote_boundary_data_on_mortar.field_data.initialize(
1099 : all_mortar_meshes.at(mortar_id).number_of_grid_points(), 0.);
1100 : mortar_data.local_insert(temporal_id, remote_boundary_data_on_mortar);
1101 : // Modify "received" mortar data
1102 : std::apply(
1103 : [&remote_boundary_data_on_mortar,
1104 : &mortar_id](const auto&... args) {
1105 : System::modify_boundary_data::apply(
1106 : make_not_null(&get<PrimalFields>(
1107 : remote_boundary_data_on_mortar.field_data))...,
1108 : make_not_null(&get<::Tags::NormalDotFlux<PrimalFields>>(
1109 : remote_boundary_data_on_mortar.field_data))...,
1110 : mortar_id, args...);
1111 : },
1112 : modify_boundary_data_args);
1113 : // Insert as remote mortar data
1114 : mortar_data.remote_insert(temporal_id,
1115 : std::move(remote_boundary_data_on_mortar));
1116 : }
1117 : }
1118 : }
1119 : apply_operator(
1120 : make_not_null(&operator_applied_to_zero_vars),
1121 : make_not_null(&all_mortar_data), zero_primal_vars, zero_deriv_vars,
1122 : primal_fluxes_buffer, element, mesh, inv_jacobian, det_inv_jacobian,
1123 : det_jacobian, det_times_inv_jacobian, face_normals, face_normal_vectors,
1124 : face_normal_magnitudes, face_jacobians,
1125 : face_jacobian_times_inv_jacobians, all_mortar_meshes, all_mortar_sizes,
1126 : mortar_jacobians, penalty_factors, massive, formulation, temporal_id,
1127 : fluxes_args_on_faces, sources_args, modify_boundary_data_args);
1128 : // Impose the nonlinear (constant) boundary contribution as fixed sources on
1129 : // the RHS of the equations
1130 : *fixed_sources -= operator_applied_to_zero_vars;
1131 : }
1132 : };
1133 :
1134 : } // namespace detail
1135 :
1136 : /*!
1137 : * \brief Prepare data on mortars so they can be communicated to neighbors
1138 : *
1139 : * Call this function on all elements and communicate the mortar data, then call
1140 : * `elliptic::dg::apply_operator`.
1141 : */
1142 : template <typename System, bool Linearized, typename... Args>
1143 1 : void prepare_mortar_data(Args&&... args) {
1144 : detail::DgOperatorImpl<System, Linearized>::template prepare_mortar_data<
1145 : false>(std::forward<Args>(args)...);
1146 : }
1147 :
1148 : /*!
1149 : * \brief Apply the elliptic DG operator
1150 : *
1151 : * This function applies the elliptic DG operator on an element, assuming all
1152 : * data on mortars is already available. Use the
1153 : * `elliptic::dg::prepare_mortar_data` function to prepare mortar data on
1154 : * neighboring elements, then communicate the data and insert them on the
1155 : * "remote" side of the mortars before calling this function.
1156 : */
1157 : template <typename System, bool Linearized, typename... Args>
1158 1 : void apply_operator(Args&&... args) {
1159 : detail::DgOperatorImpl<System, Linearized>::apply_operator(
1160 : std::forward<Args>(args)...);
1161 : }
1162 :
1163 : /*!
1164 : * \brief For linear systems, impose inhomogeneous boundary conditions as
1165 : * contributions to the fixed sources (i.e. the RHS of the equations).
1166 : *
1167 : * This function exists because the DG operator must typically be linear, but
1168 : * even for linear elliptic equations we typically apply boundary conditions
1169 : * with a constant, and therefore nonlinear, contribution. Standard examples are
1170 : * inhomogeneous (i.e. non-zero) Dirichlet or Neumann boundary conditions. This
1171 : * nonlinear contribution can be added to the fixed sources, leaving only the
1172 : * linearized boundary conditions in the DG operator. For standard constant
1173 : * Dirichlet or Neumann boundary conditions the linearization is of course just
1174 : * zero.
1175 : *
1176 : * This function essentially feeds zero variables through the nonlinear operator
1177 : * and subtracts the result from the fixed sources: `b -= A(x=0)`.
1178 : */
1179 : template <typename System, typename... Args>
1180 1 : void impose_inhomogeneous_boundary_conditions_on_source(Args&&... args) {
1181 : detail::DgOperatorImpl<System, false>::
1182 : impose_inhomogeneous_boundary_conditions_on_source(
1183 : std::forward<Args>(args)...);
1184 : }
1185 :
1186 : } // namespace elliptic::dg
|
