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