SpECTRE Documentation Coverage Report
Current view: top level - Elliptic/DiscontinuousGalerkin - DgOperator.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 7 9 77.8 %
Date: 2024-03-29 00:33:31
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             :         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

Generated by: LCOV version 1.14