SpECTRE Documentation Coverage Report
Current view: top level - Elliptic/DiscontinuousGalerkin - DgOperator.hpp Hit Total Coverage
Commit: 1afd040526f6585adf96eb1710263813d4f5ac0e Lines: 7 9 77.8 %
Date: 2024-05-23 14:31:37
Legend: Lines: hit not hit

          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             :         Variables<tmpl::list<DerivTags...>> deriv_vars_on_boundary{};
     459             :         if (AllDataIsZero or local_data_is_zero) {
     460             :           deriv_vars_on_boundary.initialize(face_num_points, 0.);
     461             :         } else {
     462             :           deriv_vars_on_boundary.initialize(face_num_points);
     463             :           ::dg::project_contiguous_data_to_boundary(
     464             :               make_not_null(&deriv_vars_on_boundary), *deriv_vars, mesh,
     465             :               direction);
     466             :         }
     467             :         apply_boundary_condition(
     468             :             direction,
     469             :             make_not_null(&get<PrimalMortarVars>(boundary_data.field_data))...,
     470             :             make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
     471             :                 boundary_data.field_data))...,
     472             :             get<DerivTags>(deriv_vars_on_boundary)...);
     473             : 
     474             :         // Invert the sign of the fluxes to account for the inverted normal on
     475             :         // exterior faces. Also multiply by 2 and add the interior fluxes to
     476             :         // impose the boundary conditions on the _average_ instead of just
     477             :         // setting the fields on the exterior:
     478             :         //   (Dirichlet)  u_D = avg(u) = 1/2 (u_int + u_ext)
     479             :         //                => u_ext = 2 u_D - u_int
     480             :         //   (Neumann)    (n.F)_N = avg(n.F) = 1/2 [(n.F)_int - (n.F)_ext]
     481             :         //                => (n.F)_ext = -2 (n.F)_N + (n.F)_int]
     482             :         const auto impose_on_average = [](const auto exterior_field,
     483             :                                           const auto& interior_field) {
     484             :           for (size_t i = 0; i < interior_field.size(); ++i) {
     485             :             (*exterior_field)[i] *= 2.;
     486             :             (*exterior_field)[i] -= interior_field[i];
     487             :           }
     488             :         };
     489             :         EXPAND_PACK_LEFT_TO_RIGHT(impose_on_average(
     490             :             make_not_null(&get<PrimalMortarVars>(boundary_data.field_data)),
     491             :             get<PrimalMortarVars>(all_mortar_data->at(mortar_id)
     492             :                                       .local_data(temporal_id)
     493             :                                       .field_data)));
     494             :         const auto invert_sign_and_impose_on_average =
     495             :             [](const auto exterior_n_dot_flux,
     496             :                const auto& interior_n_dot_flux) {
     497             :               for (size_t i = 0; i < interior_n_dot_flux.size(); ++i) {
     498             :                 (*exterior_n_dot_flux)[i] *= -2.;
     499             :                 (*exterior_n_dot_flux)[i] += interior_n_dot_flux[i];
     500             :               }
     501             :             };
     502             :         EXPAND_PACK_LEFT_TO_RIGHT(invert_sign_and_impose_on_average(
     503             :             make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
     504             :                 boundary_data.field_data)),
     505             :             get<::Tags::NormalDotFlux<PrimalMortarVars>>(
     506             :                 all_mortar_data->at(mortar_id)
     507             :                     .local_data(temporal_id)
     508             :                     .field_data)));
     509             : 
     510             :         // Store the exterior boundary data on the mortar
     511             :         all_mortar_data->at(mortar_id).remote_insert(temporal_id,
     512             :                                                      std::move(boundary_data));
     513             :       }  // if (is_internal)
     514             :     }    // loop directions
     515             :   }
     516             : 
     517             :   // --- This is essentially a break to communicate the mortar data ---
     518             : 
     519             :   template <bool AllDataIsZero, typename... OperatorTags,
     520             :             typename... PrimalVars, typename... PrimalFluxesVars,
     521             :             typename... PrimalMortarVars, typename... PrimalMortarFluxes,
     522             :             typename TemporalId, typename... FluxesArgs,
     523             :             typename... SourcesArgs, typename DataIsZero = NoDataIsZero,
     524             :             typename DirectionsPredicate = AllDirections>
     525             :   static void apply_operator(
     526             :       const gsl::not_null<Variables<tmpl::list<OperatorTags...>>*>
     527             :           operator_applied_to_vars,
     528             :       const gsl::not_null<::dg::MortarMap<
     529             :           Dim, MortarData<TemporalId, tmpl::list<PrimalMortarVars...>,
     530             :                           tmpl::list<PrimalMortarFluxes...>>>*>
     531             :           all_mortar_data,
     532             :       const Variables<tmpl::list<PrimalVars...>>& primal_vars,
     533             :       // Taking the primal fluxes computed in the `prepare_mortar_data` function
     534             :       // by const-ref here because other code might use them and so we don't
     535             :       // want to modify them by adding boundary corrections. E.g. linearized
     536             :       // sources use the nonlinear fields and fluxes as background fields.
     537             :       const Variables<tmpl::list<PrimalFluxesVars...>>& primal_fluxes,
     538             :       const Element<Dim>& element, const Mesh<Dim>& mesh,
     539             :       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     540             :                             Frame::Inertial>& inv_jacobian,
     541             :       const Scalar<DataVector>& det_inv_jacobian,
     542             :       const Scalar<DataVector>& det_jacobian,
     543             :       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     544             :                             Frame::Inertial>& det_times_inv_jacobian,
     545             :       const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
     546             :       const DirectionMap<Dim, tnsr::I<DataVector, Dim>>& face_normal_vectors,
     547             :       const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
     548             :       const DirectionMap<Dim, Scalar<DataVector>>& face_jacobians,
     549             :       const DirectionMap<Dim,
     550             :                          InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     551             :                                          Frame::Inertial>>&
     552             :           face_jacobian_times_inv_jacobians,
     553             :       const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
     554             :       const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
     555             :       const ::dg::MortarMap<Dim, Scalar<DataVector>>& mortar_jacobians,
     556             :       const ::dg::MortarMap<Dim, Scalar<DataVector>>& penalty_factors,
     557             :       const bool massive, const ::dg::Formulation formulation,
     558             :       const TemporalId& /*temporal_id*/,
     559             :       const DirectionMap<Dim, std::tuple<FluxesArgs...>>& fluxes_args_on_faces,
     560             :       const std::tuple<SourcesArgs...>& sources_args,
     561             :       const DataIsZero& data_is_zero = NoDataIsZero{},
     562             :       const DirectionsPredicate& directions_predicate = AllDirections{}) {
     563             :     static_assert(
     564             :         sizeof...(PrimalVars) == sizeof...(PrimalFields) and
     565             :             sizeof...(PrimalFluxesVars) == sizeof...(PrimalFluxes) and
     566             :             sizeof...(PrimalMortarVars) == sizeof...(PrimalFields) and
     567             :             sizeof...(PrimalMortarFluxes) == sizeof...(PrimalFluxes) and
     568             :             sizeof...(OperatorTags) == sizeof...(PrimalFields),
     569             :         "The number of variables must match the number of system fields.");
     570             :     static_assert(
     571             :         (std::is_same_v<typename PrimalVars::type,
     572             :                         typename PrimalFields::type> and
     573             :          ...) and
     574             :             (std::is_same_v<typename PrimalFluxesVars::type,
     575             :                             typename PrimalFluxes::type> and
     576             :              ...) and
     577             :             (std::is_same_v<typename PrimalMortarVars::type,
     578             :                             typename PrimalFields::type> and
     579             :              ...) and
     580             :             (std::is_same_v<typename PrimalMortarFluxes::type,
     581             :                             typename PrimalFluxes::type> and
     582             :              ...) and
     583             :             (std::is_same_v<typename OperatorTags::type,
     584             :                             typename PrimalFields::type> and
     585             :              ...),
     586             :         "The variables must have the same tensor types as the system fields.");
     587             : #ifdef SPECTRE_DEBUG
     588             :     for (size_t d = 0; d < Dim; ++d) {
     589             :       ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
     590             :                  (mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto or
     591             :                   mesh.quadrature(d) == Spectral::Quadrature::Gauss),
     592             :              "The elliptic DG operator is currently only implemented for "
     593             :              "Legendre-Gauss(-Lobatto) grids. Found basis '"
     594             :                  << mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
     595             :                  << "' in dimension " << d << ".");
     596             :     }
     597             : #endif  // SPECTRE_DEBUG
     598             :     const auto& element_id = element.id();
     599             :     const bool local_data_is_zero = data_is_zero(element_id);
     600             :     ASSERT(Linearized or not local_data_is_zero,
     601             :            "Only a linear operator can take advantage of the knowledge that "
     602             :            "the operand is zero. Don't return 'true' in 'data_is_zero' unless "
     603             :            "you also set 'Linearized' to 'true'.");
     604             :     const size_t num_points = mesh.number_of_grid_points();
     605             : 
     606             :     // This function and the one above allocate various Variables to compute
     607             :     // intermediate quantities. It could be a performance optimization to reduce
     608             :     // the number of these allocations and/or move some of the memory buffers
     609             :     // into the DataBox to keep them around permanently. The latter should be
     610             :     // informed by profiling.
     611             : 
     612             :     // Compute volume terms: -div(F) + S
     613             :     if (local_data_is_zero) {
     614             :       operator_applied_to_vars->initialize(num_points, 0.);
     615             :     } else {
     616             :       // "Massive" operators retain the factors from the volume integral:
     617             :       //   \int_volume div(F) \phi_p = w_p det(J)_p div(F)_p
     618             :       // Here, `w` are the quadrature weights (the diagonal logical mass matrix
     619             :       // with mass-lumping) and det(J) is the Jacobian determinant. The
     620             :       // quantities are evaluated at the grid point `p`.
     621             :       if (formulation == ::dg::Formulation::StrongInertial) {
     622             :         // Compute strong divergence:
     623             :         //   div(F) = (J^\hat{i}_i)_p \sum_q (D_\hat{i})_pq (F^i)_q.
     624             :         divergence(operator_applied_to_vars, primal_fluxes, mesh,
     625             :                    massive ? det_times_inv_jacobian : inv_jacobian);
     626             :         // This is the sign flip that makes the operator _minus_ the Laplacian
     627             :         // for a Poisson system
     628             :         *operator_applied_to_vars *= -1.;
     629             :       } else {
     630             :         // Compute weak divergence:
     631             :         //   F^i \partial_i \phi = 1/w_p \sum_q
     632             :         //     (D^T_\hat{i})_pq (w det(J) J^\hat{i}_i F^i)_q
     633             :         weak_divergence(operator_applied_to_vars, primal_fluxes, mesh,
     634             :                         det_times_inv_jacobian);
     635             :         if (not massive) {
     636             :           *operator_applied_to_vars *= get(det_inv_jacobian);
     637             :         }
     638             :       }
     639             :       if constexpr (not std::is_same_v<SourcesComputer, void>) {
     640             :         Variables<tmpl::list<OperatorTags...>> sources{num_points, 0.};
     641             :         std::apply(
     642             :             [&sources, &primal_vars,
     643             :              &primal_fluxes](const auto&... expanded_sources_args) {
     644             :               SourcesComputer::apply(
     645             :                   make_not_null(&get<OperatorTags>(sources))...,
     646             :                   expanded_sources_args..., get<PrimalVars>(primal_vars)...,
     647             :                   get<PrimalFluxesVars>(primal_fluxes)...);
     648             :             },
     649             :             sources_args);
     650             :         if (massive) {
     651             :           sources *= get(det_jacobian);
     652             :         }
     653             :         *operator_applied_to_vars += sources;
     654             :       }
     655             :     }
     656             :     if (massive) {
     657             :       ::dg::apply_mass_matrix(operator_applied_to_vars, mesh);
     658             :     }
     659             : 
     660             :     // Add boundary corrections
     661             :     // Keeping track if any corrections were applied here, for an optimization
     662             :     // below
     663             :     bool has_any_boundary_corrections = false;
     664             :     Variables<tmpl::list<transform::Tags::TransformedFirstIndex<
     665             :         PrimalFluxesVars, Frame::ElementLogical>...>>
     666             :         lifted_logical_aux_boundary_corrections{num_points, 0.};
     667             :     for (auto& [mortar_id, mortar_data] : *all_mortar_data) {
     668             :       const auto& direction = mortar_id.direction();
     669             :       const auto& neighbor_id = mortar_id.id();
     670             :       const bool is_internal =
     671             :           (neighbor_id != ElementId<Dim>::external_boundary_id());
     672             :       if constexpr (AllDataIsZero) {
     673             :         if (is_internal) {
     674             :           continue;
     675             :         }
     676             :       }
     677             :       if (not directions_predicate(direction)) {
     678             :         continue;
     679             :       }
     680             :       // When the data on both sides of the mortar is zero then we don't need to
     681             :       // handle this mortar at all.
     682             :       if (local_data_is_zero and
     683             :           (not is_internal or data_is_zero(neighbor_id))) {
     684             :         continue;
     685             :       }
     686             :       has_any_boundary_corrections = true;
     687             : 
     688             :       const auto face_mesh = mesh.slice_away(direction.dimension());
     689             :       auto [local_data, remote_data] = mortar_data.extract();
     690             :       const size_t face_num_points = face_mesh.number_of_grid_points();
     691             :       const auto& face_normal = face_normals.at(direction);
     692             :       const auto& face_normal_vector = face_normal_vectors.at(direction);
     693             :       const auto& fluxes_args_on_face = fluxes_args_on_faces.at(direction);
     694             :       const auto& face_normal_magnitude = face_normal_magnitudes.at(direction);
     695             :       const auto& face_jacobian = face_jacobians.at(direction);
     696             :       const auto& face_jacobian_times_inv_jacobian =
     697             :           face_jacobian_times_inv_jacobians.at(direction);
     698             :       const auto& mortar_mesh =
     699             :           is_internal ? all_mortar_meshes.at(mortar_id) : face_mesh;
     700             :       const auto& mortar_size =
     701             :           is_internal ? all_mortar_sizes.at(mortar_id) : full_mortar_size;
     702             : 
     703             :       // This is the strong auxiliary boundary correction:
     704             :       //   G^i = F^i(n_j (avg(u) - u))
     705             :       // where
     706             :       //   avg(u) - u = -0.5 * (u_int - u_ext)
     707             :       auto avg_vars_on_mortar = Variables<tmpl::list<PrimalMortarVars...>>(
     708             :           local_data.field_data
     709             :               .template extract_subset<tmpl::list<PrimalMortarVars...>>());
     710             :       const auto add_remote_contribution = [](auto& lhs, const auto& rhs) {
     711             :         for (size_t i = 0; i < lhs.size(); ++i) {
     712             :           lhs[i] -= rhs[i];
     713             :         }
     714             :       };
     715             :       EXPAND_PACK_LEFT_TO_RIGHT(add_remote_contribution(
     716             :           get<PrimalMortarVars>(avg_vars_on_mortar),
     717             :           get<PrimalMortarVars>(remote_data.field_data)));
     718             :       avg_vars_on_mortar *= -0.5;
     719             : 
     720             :       // Project from the mortar back down to the face if needed
     721             :       const auto avg_vars_on_face =
     722             :           Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)
     723             :               ? mass_conservative_restriction(
     724             :                     std::move(avg_vars_on_mortar), mortar_mesh, mortar_size,
     725             :                     mortar_jacobians.at(mortar_id), face_mesh, face_jacobian)
     726             :               : std::move(avg_vars_on_mortar);
     727             : 
     728             :       // Apply fluxes to get G^i
     729             :       Variables<tmpl::list<PrimalFluxesVars...>> auxiliary_boundary_corrections{
     730             :           face_num_points};
     731             :       std::apply(
     732             :           [&auxiliary_boundary_corrections, &face_normal, &face_normal_vector,
     733             :            &avg_vars_on_face,
     734             :            &element_id](const auto&... expanded_fluxes_args_on_face) {
     735             :             if constexpr (FluxesComputer::is_discontinuous) {
     736             :               FluxesComputer::apply(make_not_null(&get<PrimalFluxesVars>(
     737             :                                         auxiliary_boundary_corrections))...,
     738             :                                     expanded_fluxes_args_on_face..., element_id,
     739             :                                     face_normal, face_normal_vector,
     740             :                                     get<PrimalMortarVars>(avg_vars_on_face)...);
     741             :             } else {
     742             :               (void)element_id;
     743             :               FluxesComputer::apply(make_not_null(&get<PrimalFluxesVars>(
     744             :                                         auxiliary_boundary_corrections))...,
     745             :                                     expanded_fluxes_args_on_face...,
     746             :                                     face_normal, face_normal_vector,
     747             :                                     get<PrimalMortarVars>(avg_vars_on_face)...);
     748             :             }
     749             :           },
     750             :           fluxes_args_on_face);
     751             : 
     752             :       // Lifting for the auxiliary boundary correction:
     753             :       //   \int_face G^i \partial_i \phi
     754             :       // We first transform the flux index to the logical frame, apply the
     755             :       // quadrature weights and the Jacobian for the face integral, then take
     756             :       // the logical weak divergence in the volume after lifting (below the loop
     757             :       // over faces).
     758             :       auto logical_aux_boundary_corrections =
     759             :           transform::first_index_to_different_frame(
     760             :               auxiliary_boundary_corrections, face_jacobian_times_inv_jacobian);
     761             :       ::dg::apply_mass_matrix(make_not_null(&logical_aux_boundary_corrections),
     762             :                               face_mesh);
     763             :       if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
     764             :         add_slice_to_data(
     765             :             make_not_null(&lifted_logical_aux_boundary_corrections),
     766             :             logical_aux_boundary_corrections, mesh.extents(),
     767             :             direction.dimension(),
     768             :             index_to_slice_at(mesh.extents(), direction));
     769             :       } else {
     770             :         ::dg::lift_boundary_terms_gauss_points(
     771             :             make_not_null(&lifted_logical_aux_boundary_corrections),
     772             :             logical_aux_boundary_corrections, mesh, direction);
     773             :       }
     774             : 
     775             :       // This is the strong primal boundary correction:
     776             :       //   -n.H = -avg(n.F) + n.F + penalty * n.F(n_j jump(u))
     777             :       // Note that the "internal penalty" numerical flux
     778             :       // (as opposed to the LLF flux) uses the raw field derivatives without
     779             :       // boundary corrections in the average, which is why we can communicate
     780             :       // the data so early together with the auxiliary boundary data. In this
     781             :       // case the penalty needs to include a factor N_points^2 / h (see the
     782             :       // `penalty` function).
     783             :       const auto& penalty_factor = penalty_factors.at(mortar_id);
     784             :       // Compute jump on mortar:
     785             :       //   penalty * jump(u) = penalty * (u_int - u_ext)
     786             :       const auto add_remote_jump_contribution =
     787             :           [&penalty_factor](auto& lhs, const auto& rhs) {
     788             :             for (size_t i = 0; i < lhs.size(); ++i) {
     789             :               lhs[i] -= rhs[i];
     790             :               lhs[i] *= get(penalty_factor);
     791             :             }
     792             :           };
     793             :       EXPAND_PACK_LEFT_TO_RIGHT(add_remote_jump_contribution(
     794             :           get<PrimalMortarVars>(local_data.field_data),
     795             :           get<PrimalMortarVars>(remote_data.field_data)));
     796             :       // Compute average on mortar:
     797             :       //   (strong)  -avg(n.F) + n.F = 0.5 * (n.F)_int + 0.5 * (n.F)_ext
     798             :       //   (weak)    -avg(n.F) = -0.5 * (n.F)_int + 0.5 * (n.F)_ext
     799             :       const auto add_avg_contribution = [](auto& lhs, const auto& rhs,
     800             :                                            const double factor) {
     801             :         for (size_t i = 0; i < lhs.size(); ++i) {
     802             :           lhs[i] *= factor;
     803             :           lhs[i] += 0.5 * rhs[i];
     804             :         }
     805             :       };
     806             :       EXPAND_PACK_LEFT_TO_RIGHT(add_avg_contribution(
     807             :           get<::Tags::NormalDotFlux<PrimalMortarVars>>(local_data.field_data),
     808             :           get<::Tags::NormalDotFlux<PrimalMortarVars>>(remote_data.field_data),
     809             :           formulation == ::dg::Formulation::StrongInertial ? 0.5 : -0.5));
     810             : 
     811             :       // Project from the mortar back down to the face if needed, lift and add
     812             :       // to operator. See auxiliary boundary corrections above for details.
     813             :       auto primal_boundary_corrections_on_face =
     814             :           Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)
     815             :               ? mass_conservative_restriction(
     816             :                     std::move(local_data.field_data), mortar_mesh, mortar_size,
     817             :                     mortar_jacobians.at(mortar_id), face_mesh, face_jacobian)
     818             :               : std::move(local_data.field_data);
     819             : 
     820             :       // Compute fluxes for jump term: n.F(n_j jump(u))
     821             :       // If the fluxes are trivial (just the spatial metric), we can skip this
     822             :       // step because the face normal is normalized.
     823             :       if constexpr (not FluxesComputer::is_trivial) {
     824             :         // We reuse the memory buffer from above for the result.
     825             :         std::apply(
     826             :             [&auxiliary_boundary_corrections, &face_normal, &face_normal_vector,
     827             :              &primal_boundary_corrections_on_face,
     828             :              &element_id](const auto&... expanded_fluxes_args_on_face) {
     829             :               if constexpr (FluxesComputer::is_discontinuous) {
     830             :                 FluxesComputer::apply(
     831             :                     make_not_null(&get<PrimalFluxesVars>(
     832             :                         auxiliary_boundary_corrections))...,
     833             :                     expanded_fluxes_args_on_face..., element_id, face_normal,
     834             :                     face_normal_vector,
     835             :                     get<PrimalMortarVars>(
     836             :                         primal_boundary_corrections_on_face)...);
     837             :               } else {
     838             :                 (void)element_id;
     839             :                 FluxesComputer::apply(
     840             :                     make_not_null(&get<PrimalFluxesVars>(
     841             :                         auxiliary_boundary_corrections))...,
     842             :                     expanded_fluxes_args_on_face..., face_normal,
     843             :                     face_normal_vector,
     844             :                     get<PrimalMortarVars>(
     845             :                         primal_boundary_corrections_on_face)...);
     846             :               }
     847             :             },
     848             :             fluxes_args_on_face);
     849             :         if constexpr (FluxesComputer::is_discontinuous) {
     850             :           if (is_internal) {
     851             :             // For penalty term with discontinuous fluxes: evaluate the fluxes
     852             :             // on the other side of the boundary as well and take average
     853             :             Variables<tmpl::list<PrimalFluxesVars...>> fluxes_other_side{
     854             :                 face_num_points};
     855             :             std::apply(
     856             :                 [&fluxes_other_side, &face_normal, &face_normal_vector,
     857             :                  &primal_boundary_corrections_on_face,
     858             :                  &local_neighbor_id =
     859             :                      neighbor_id](const auto&... expanded_fluxes_args_on_face) {
     860             :                   FluxesComputer::apply(
     861             :                       make_not_null(
     862             :                           &get<PrimalFluxesVars>(fluxes_other_side))...,
     863             :                       expanded_fluxes_args_on_face..., local_neighbor_id,
     864             :                       face_normal, face_normal_vector,
     865             :                       get<PrimalMortarVars>(
     866             :                           primal_boundary_corrections_on_face)...);
     867             :                 },
     868             :                 fluxes_args_on_face);
     869             :             auxiliary_boundary_corrections += fluxes_other_side;
     870             :             auxiliary_boundary_corrections *= 0.5;
     871             :           }
     872             :         }
     873             :         EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
     874             :             make_not_null(
     875             :                 &get<PrimalMortarVars>(primal_boundary_corrections_on_face)),
     876             :             face_normal,
     877             :             get<PrimalFluxesVars>(auxiliary_boundary_corrections)));
     878             :       }
     879             : 
     880             :       // Add penalty term to average term
     881             :       Variables<tmpl::list<PrimalMortarVars...>> primal_boundary_corrections{};
     882             :       // First half of the memory allocated above is filled with the penalty
     883             :       // term, so just use that memory here.
     884             :       primal_boundary_corrections.set_data_ref(
     885             :           primal_boundary_corrections_on_face.data(), avg_vars_on_face.size());
     886             :       // Second half of the memory is filled with the average term. Add that to
     887             :       // the penalty term.
     888             :       const auto add_avg_term = [](auto& lhs, const auto& rhs) {
     889             :         for (size_t i = 0; i < lhs.size(); ++i) {
     890             :           lhs[i] += rhs[i];
     891             :         }
     892             :       };
     893             :       EXPAND_PACK_LEFT_TO_RIGHT(
     894             :           add_avg_term(get<PrimalMortarVars>(primal_boundary_corrections),
     895             :                        get<::Tags::NormalDotFlux<PrimalMortarVars>>(
     896             :                            primal_boundary_corrections_on_face)));
     897             : 
     898             :       // Lifting for the primal boundary correction:
     899             :       //   \int_face n.H \phi
     900             :       if (massive) {
     901             :         // We apply the quadrature weights and Jacobian for the face integral,
     902             :         // then lift to the volume
     903             :         primal_boundary_corrections *= get(face_jacobian);
     904             :         ::dg::apply_mass_matrix(make_not_null(&primal_boundary_corrections),
     905             :                                 face_mesh);
     906             :         if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
     907             :           add_slice_to_data(operator_applied_to_vars,
     908             :                             primal_boundary_corrections, mesh.extents(),
     909             :                             direction.dimension(),
     910             :                             index_to_slice_at(mesh.extents(), direction));
     911             :         } else {
     912             :           ::dg::lift_boundary_terms_gauss_points(operator_applied_to_vars,
     913             :                                                  primal_boundary_corrections,
     914             :                                                  mesh, direction);
     915             :         }
     916             :       } else {
     917             :         // Apply an extra inverse mass matrix to the boundary corrections (with
     918             :         // mass lumping, so it's diagonal).
     919             :         // For Gauss-Lobatto grids this divides out the quadrature weights and
     920             :         // Jacobian on the face, leaving only factors perpendicular to the face.
     921             :         // Those are handled by `dg::lift_flux` (with an extra minus sign since
     922             :         // the function was written for evolution systems).
     923             :         // For Gauss grids the quadrature weights and Jacobians are handled
     924             :         // by `::dg::lift_boundary_terms_gauss_points` (which was also written
     925             :         // for evolution systems, hence the extra minus sign).
     926             :         primal_boundary_corrections *= -1.;
     927             :         if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
     928             :           ::dg::lift_flux(make_not_null(&primal_boundary_corrections),
     929             :                           mesh.extents(direction.dimension()),
     930             :                           face_normal_magnitude);
     931             :           add_slice_to_data(operator_applied_to_vars,
     932             :                             primal_boundary_corrections, mesh.extents(),
     933             :                             direction.dimension(),
     934             :                             index_to_slice_at(mesh.extents(), direction));
     935             :         } else {
     936             :           // We already have the `face_jacobian = det(J) * magnitude(n)` here,
     937             :           // so just pass a constant 1 for `magnitude(n)`. This could be
     938             :           // optimized to avoid allocating the vector of ones.
     939             :           ::dg::lift_boundary_terms_gauss_points(
     940             :               operator_applied_to_vars, det_inv_jacobian, mesh, direction,
     941             :               primal_boundary_corrections,
     942             :               Scalar<DataVector>{face_mesh.number_of_grid_points(), 1.},
     943             :               face_jacobian);
     944             :         }
     945             :       }
     946             :     }  // loop over all mortars
     947             : 
     948             :     if (not has_any_boundary_corrections) {
     949             :       // No need to handle auxiliary boundary corrections; return early
     950             :       return;
     951             :     }
     952             : 
     953             :     // Apply weak divergence to lifted auxiliary boundary corrections and add to
     954             :     // operator
     955             :     if (massive) {
     956             :       logical_weak_divergence(operator_applied_to_vars,
     957             :                               lifted_logical_aux_boundary_corrections, mesh,
     958             :                               true);
     959             :     } else {
     960             :       // Possible optimization: eliminate this allocation by building the
     961             :       // inverse mass matrix into `logical_weak_divergence`
     962             :       Variables<tmpl::list<OperatorTags...>> massless_aux_boundary_corrections{
     963             :           num_points};
     964             :       logical_weak_divergence(make_not_null(&massless_aux_boundary_corrections),
     965             :                               lifted_logical_aux_boundary_corrections, mesh);
     966             :       massless_aux_boundary_corrections *= get(det_inv_jacobian);
     967             :       ::dg::apply_inverse_mass_matrix(
     968             :           make_not_null(&massless_aux_boundary_corrections), mesh);
     969             :       *operator_applied_to_vars += massless_aux_boundary_corrections;
     970             :     }
     971             :   }
     972             : 
     973             :   template <typename... FixedSourcesTags, typename ApplyBoundaryCondition,
     974             :             typename... FluxesArgs, typename... SourcesArgs,
     975             :             bool LocalLinearized = Linearized,
     976             :             // This function adds nothing to the fixed sources if the operator
     977             :             // is linearized, so it shouldn't be used in that case
     978             :             Requires<not LocalLinearized> = nullptr>
     979             :   static void impose_inhomogeneous_boundary_conditions_on_source(
     980             :       const gsl::not_null<Variables<tmpl::list<FixedSourcesTags...>>*>
     981             :           fixed_sources,
     982             :       const Element<Dim>& element, const Mesh<Dim>& mesh,
     983             :       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     984             :                             Frame::Inertial>& inv_jacobian,
     985             :       const Scalar<DataVector>& det_inv_jacobian,
     986             :       const Scalar<DataVector>& det_jacobian,
     987             :       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     988             :                             Frame::Inertial>& det_times_inv_jacobian,
     989             :       const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
     990             :       const DirectionMap<Dim, tnsr::I<DataVector, Dim>>& face_normal_vectors,
     991             :       const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
     992             :       const DirectionMap<Dim, Scalar<DataVector>>& face_jacobians,
     993             :       const DirectionMap<Dim,
     994             :                          InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     995             :                                          Frame::Inertial>>&
     996             :           face_jacobian_times_inv_jacobians,
     997             :       const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
     998             :       const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
     999             :       const ::dg::MortarMap<Dim, Scalar<DataVector>>& penalty_factors,
    1000             :       const bool massive, const ::dg::Formulation formulation,
    1001             :       const ApplyBoundaryCondition& apply_boundary_condition,
    1002             :       const std::tuple<FluxesArgs...>& fluxes_args,
    1003             :       const std::tuple<SourcesArgs...>& sources_args,
    1004             :       const DirectionMap<Dim, std::tuple<FluxesArgs...>>&
    1005             :           fluxes_args_on_faces) {
    1006             :     // We just feed zero variables through the nonlinear operator to extract the
    1007             :     // constant contribution at external boundaries. Since the variables are
    1008             :     // zero the operator simplifies quite a lot. The simplification is probably
    1009             :     // not very important for performance because this function will only be
    1010             :     // called when solving a linear elliptic system and only once during
    1011             :     // initialization, but we specialize the operator for zero data nonetheless
    1012             :     // just so we can ignore internal boundaries. For internal boundaries we
    1013             :     // would unnecessarily have to copy mortar data around to emulate the
    1014             :     // communication step, so by just skipping internal boundaries we avoid
    1015             :     // that.
    1016             :     const size_t num_points = mesh.number_of_grid_points();
    1017             :     const Variables<tmpl::list<PrimalFields...>> zero_primal_vars{num_points,
    1018             :                                                                   0.};
    1019             :     Variables<tmpl::list<PrimalFluxes...>> primal_fluxes_buffer{num_points};
    1020             :     Variables<tmpl::list<
    1021             :         ::Tags::deriv<PrimalFields, tmpl::size_t<Dim>, Frame::Inertial>...>>
    1022             :         unused_deriv_vars_buffer{};
    1023             :     Variables<tmpl::list<FixedSourcesTags...>> operator_applied_to_zero_vars{
    1024             :         num_points};
    1025             :     // Set up data on mortars. We only need them at external boundaries.
    1026             :     ::dg::MortarMap<Dim, MortarData<size_t, tmpl::list<PrimalFields...>,
    1027             :                                     tmpl::list<PrimalFluxes...>>>
    1028             :         all_mortar_data{};
    1029             :     constexpr size_t temporal_id = std::numeric_limits<size_t>::max();
    1030             :     // Apply the operator to the zero variables, skipping internal boundaries
    1031             :     prepare_mortar_data<true>(make_not_null(&unused_deriv_vars_buffer),
    1032             :                               make_not_null(&primal_fluxes_buffer),
    1033             :                               make_not_null(&all_mortar_data), zero_primal_vars,
    1034             :                               element, mesh, inv_jacobian, face_normals,
    1035             :                               all_mortar_meshes, all_mortar_sizes, temporal_id,
    1036             :                               apply_boundary_condition, fluxes_args);
    1037             :     apply_operator<true>(
    1038             :         make_not_null(&operator_applied_to_zero_vars),
    1039             :         make_not_null(&all_mortar_data), zero_primal_vars, primal_fluxes_buffer,
    1040             :         element, mesh, inv_jacobian, det_inv_jacobian, det_jacobian,
    1041             :         det_times_inv_jacobian, face_normals, face_normal_vectors,
    1042             :         face_normal_magnitudes, face_jacobians,
    1043             :         face_jacobian_times_inv_jacobians, all_mortar_meshes, all_mortar_sizes,
    1044             :         {}, penalty_factors, massive, formulation, temporal_id,
    1045             :         fluxes_args_on_faces, sources_args);
    1046             :     // Impose the nonlinear (constant) boundary contribution as fixed sources on
    1047             :     // the RHS of the equations
    1048             :     *fixed_sources -= operator_applied_to_zero_vars;
    1049             :   }
    1050             : };
    1051             : 
    1052             : }  // namespace detail
    1053             : 
    1054             : /*!
    1055             :  * \brief Prepare data on mortars so they can be communicated to neighbors
    1056             :  *
    1057             :  * Call this function on all elements and communicate the mortar data, then call
    1058             :  * `elliptic::dg::apply_operator`.
    1059             :  */
    1060             : template <typename System, bool Linearized, typename... Args>
    1061           1 : void prepare_mortar_data(Args&&... args) {
    1062             :   detail::DgOperatorImpl<System, Linearized>::template prepare_mortar_data<
    1063             :       false>(std::forward<Args>(args)...);
    1064             : }
    1065             : 
    1066             : /*!
    1067             :  * \brief Apply the elliptic DG operator
    1068             :  *
    1069             :  * This function applies the elliptic DG operator on an element, assuming all
    1070             :  * data on mortars is already available. Use the
    1071             :  * `elliptic::dg::prepare_mortar_data` function to prepare mortar data on
    1072             :  * neighboring elements, then communicate the data and insert them on the
    1073             :  * "remote" side of the mortars before calling this function.
    1074             :  */
    1075             : template <typename System, bool Linearized, typename... Args>
    1076           1 : void apply_operator(Args&&... args) {
    1077             :   detail::DgOperatorImpl<System, Linearized>::template apply_operator<false>(
    1078             :       std::forward<Args>(args)...);
    1079             : }
    1080             : 
    1081             : /*!
    1082             :  * \brief For linear systems, impose inhomogeneous boundary conditions as
    1083             :  * contributions to the fixed sources (i.e. the RHS of the equations).
    1084             :  *
    1085             :  * This function exists because the DG operator must typically be linear, but
    1086             :  * even for linear elliptic equations we typically apply boundary conditions
    1087             :  * with a constant, and therefore nonlinear, contribution. Standard examples are
    1088             :  * inhomogeneous (i.e. non-zero) Dirichlet or Neumann boundary conditions. This
    1089             :  * nonlinear contribution can be added to the fixed sources, leaving only the
    1090             :  * linearized boundary conditions in the DG operator. For standard constant
    1091             :  * Dirichlet or Neumann boundary conditions the linearization is of course just
    1092             :  * zero.
    1093             :  *
    1094             :  * This function essentially feeds zero variables through the nonlinear operator
    1095             :  * and subtracts the result from the fixed sources: `b -= A(x=0)`.
    1096             :  */
    1097             : template <typename System, typename... Args>
    1098           1 : void impose_inhomogeneous_boundary_conditions_on_source(Args&&... args) {
    1099             :   detail::DgOperatorImpl<System, false>::
    1100             :       impose_inhomogeneous_boundary_conditions_on_source(
    1101             :           std::forward<Args>(args)...);
    1102             : }
    1103             : 
    1104             : }  // namespace elliptic::dg

Generated by: LCOV version 1.14