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 <type_traits>
       8             : 
       9             : #include "DataStructures/Tensor/Tensor.hpp"
      10             : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
      11             : #include "Utilities/TMPL.hpp"
      12             : #include "Utilities/TypeTraits/IsA.hpp"
      13             : 
      14             : /// \ref protocols related to elliptic systems
      15           1 : namespace elliptic::protocols {
      16             : 
      17             : namespace FirstOrderSystem_detail {
      18             : template <size_t Dim, typename PrimalFields, typename PrimalFluxes>
      19             : struct test_fields_and_fluxes;
      20             : template <size_t Dim, typename... PrimalFields, typename... PrimalFluxes>
      21             : struct test_fields_and_fluxes<Dim, tmpl::list<PrimalFields...>,
      22             :                               tmpl::list<PrimalFluxes...>> : std::true_type {
      23             :   static_assert(sizeof...(PrimalFields) == sizeof...(PrimalFluxes),
      24             :                 "The system must have the same number of fields and fluxes.");
      25             :   static_assert(
      26             :       ((tmpl::size<typename PrimalFluxes::type::index_list>::value ==
      27             :         tmpl::size<typename PrimalFields::type::index_list>::value + 1) and
      28             :        ...) and
      29             :           (std::is_same_v<tmpl::front<typename PrimalFluxes::type::index_list>,
      30             :                           SpatialIndex<Dim, UpLo::Up, Frame::Inertial>> and
      31             :            ...),
      32             :       "Primal fluxes and primal fields must correspond to each "
      33             :       "other. In particular, each primal flux must have one "
      34             :       "index more than its corresponding primal field and an upper-spatial "
      35             :       "first index.");
      36             : };
      37             : }  // namespace FirstOrderSystem_detail
      38             : 
      39             : /*!
      40             :  * \brief A system of elliptic equations in first-order "flux" formulation
      41             :  *
      42             :  * Classes conforming to this protocol represent a set of elliptic partial
      43             :  * differential equations in first-order "flux" formulation:
      44             :  *
      45             :  * \f{equation}
      46             :  * -\partial_i F^i_\alpha + S_\alpha = f_\alpha(x)
      47             :  * \f}
      48             :  *
      49             :  * in terms of fluxes \f$F_\alpha^i\f$, sources \f$S_\alpha\f$ and fixed-sources
      50             :  * \f$f_\alpha(x)\f$ \cite Fischer2021voj.
      51             :  * It resembles closely formulations of hyperbolic
      52             :  * conservation laws but allows the fluxes \f$F_\alpha^i\f$ to be higher-rank
      53             :  * tensor fields. The fluxes and sources are functionals of the system
      54             :  * variables \f$u_\alpha(x)\f$ and their derivatives. The fixed-sources
      55             :  * \f$f_\alpha(x)\f$ are independent of the system variables. See the
      56             :  * `Poisson::FirstOrderSystem` and the `Elasticity::FirstOrderSystem` for
      57             :  * examples.
      58             :  *
      59             :  * Note that this formulation has been simplified since \cite Fischer2021voj :
      60             :  * We assume that the fluxes are linear in the fields and their derivatives and
      61             :  * removed the notion of "auxiliary variables" from the formulation altogether.
      62             :  * In the language of \cite Fischer2021voj we always just choose the partial
      63             :  * derivatives of the fields as auxiliary variables.
      64             :  *
      65             :  * Conforming classes must have these static member variables:
      66             :  *
      67             :  * - `size_t volume_dim`: The number of spatial dimensions.
      68             :  *
      69             :  * Conforming classes must have these type aliases:
      70             :  *
      71             :  * - `primal_fields`: A list of tags representing the primal fields. These are
      72             :  *   the fields we solve for, e.g. \f$u\f$ for a Poisson equation.
      73             :  *   (we may rename this to just "fields" since we removed the notion of
      74             :  *   "auxiliary fields")
      75             :  *
      76             :  * - `primal_fluxes`: A list of tags representing the primal fluxes
      77             :  *   \f$F_\alpha^i\f$. These are typically some linear combination of the
      78             :  *   derivatives of the system fields with raised indices, e.g. \f$v^i = g^{ij}
      79             :  *   \partial_j u\f$ for a curved-space Poisson equation on a background metric
      80             :  *   \f$g_{ij}\f$. They must have an upper-spatial first index, because their
      81             :  *   divergence defines the elliptic equation.
      82             :  *
      83             :  * - `background_fields`: A list of tags representing the variable-independent
      84             :  *   background fields in the equations. Examples are a background metric,
      85             :  *   associated fixed geometry quantities such as Christoffel symbols or the
      86             :  *   Ricci scalar, or any other fixed field that determines the problem to be
      87             :  *   solved such as matter sources in the Einstein constraint equations.
      88             :  *
      89             :  * - `inv_metric_tag`: The tag that defines the background geometry, i.e. the
      90             :  *   the geometry that the elliptic equations are formulated on. This is the
      91             :  *   metric responsible for normalizing one-forms, such as face normals.
      92             :  *
      93             :  * - `fluxes_computer`: A class that defines the fluxes \f$F_\alpha^i\f$. Must
      94             :  *   have an `argument_tags` type alias and an `apply` function that takes these
      95             :  *   arguments in this order:
      96             :  *
      97             :  *   1. The `primal_fluxes` as not-null pointer
      98             :  *   2. The `argument_tags`
      99             :  *   3. If `is_discontinuous` is `true` (see below):
     100             :  *      const ElementId<Dim>& element_id
     101             :  *   4. The `primal_fields`
     102             :  *   5. The partial derivatives of the `primal_fields`
     103             :  *
     104             :  *   The function can assume the output buffers are already correctly sized,
     105             :  *   but no guarantee is made on the values that the buffers hold at input.
     106             :  *
     107             :  *   The `fluxes_computer` must also have an `apply` function overload that is
     108             :  *   evaluated on faces of DG elements. It computes the same fluxes
     109             :  *   \f$F_\alpha^i\f$, but with the field derivatives replaced by the the face
     110             :  *   normal times the fields, and with the non-principal (non-derivative) terms
     111             :  *   set to zero. Having this separate function is an optimization to take
     112             :  *   advantage of the face normal remaining constant throughout the solve, so it
     113             :  *   can be "baked in" to the flux. The function takes these arguments in this
     114             :  *   order:
     115             :  *
     116             :  *   1. The `primal_fluxes` as not-null pointer
     117             :  *   2. The `argument_tags`
     118             :  *   3. If `is_discontinuous` is `true` (see below):
     119             :  *      const ElementId<Dim>& element_id
     120             :  *   4. The `const tnsr::i<DataVector, Dim>& face_normal` ($n_i$)
     121             :  *   5. The `const tnsr::I<DataVector, Dim>& face_normal_vector` ($n^i$)
     122             :  *   6. The `primal_fields`
     123             :  *
     124             :  *   The `fluxes_computer` class must also have the following additional
     125             :  *   type aliases and static member variables:
     126             :  *
     127             :  *   - `volume_tags`: the subset of `argument_tags` that will be retrieved
     128             :  *     directly from the DataBox, instead of retrieving it from the face of an
     129             :  *     element, when fluxes are applied on a face.
     130             :  *   - `const_global_cache_tags`: the subset of `argument_tags` that can be
     131             :  *     retrieved from _any_ element's DataBox, because they are stored in the
     132             :  *     global cache.
     133             :  *   - `bool is_trivial`: a boolean indicating whether the fluxes are simply
     134             :  *     the spatial metric, as is the case for the Poisson equation. Some
     135             :  *     computations can be skipped in this case.
     136             :  *   - `bool is_discontinuous`: a boolean indicating whether the fluxes are
     137             :  *     potentially discontinuous across element boundaries. This is `true` for
     138             :  *     systems where the equations on both sides of the boundary are different,
     139             :  *     e.g. elasticity with different materials on either side of the boundary.
     140             :  *     An additional `element_id` argument is passed to the `apply` functions in
     141             :  *     this case, which identifies on which side of an element boundary the
     142             :  *     fluxes are being evaluated.
     143             :  *
     144             :  * - `sources_computer`: A class that defines the sources \f$S_\alpha\f$. Must
     145             :  *   have an `argument_tags` type alias and an `apply` function that adds the
     146             :  *   sources to the equations. It takes these arguments in this order:
     147             :  *
     148             :  *   1. The types of the `primal_fields` as not-null pointer. These are the
     149             :  *      primal equations.
     150             :  *   2. The `argument_tags`
     151             :  *   3. The `primal_fields`
     152             :  *   4. The `primal_fluxes`
     153             :  *
     154             :  *   The function is expected to _add_ the sources \f$S_\alpha\f$ to the
     155             :  *   output buffers.
     156             :  *   The `sources_computer` may also be `void`, in which case \f$S_\alpha=0\f$.
     157             :  *
     158             :  * - `boundary_conditions_base`: A base class representing the supported
     159             :  *   boundary conditions. Boundary conditions can be factory-created from this
     160             :  *   base class. Currently this should be a specialization of
     161             :  *   `elliptic::BoundaryConditions::BoundaryCondition`.
     162             :  */
     163           1 : struct FirstOrderSystem {
     164             :   template <typename ConformingType>
     165           0 :   struct test {
     166           0 :     static constexpr size_t volume_dim = ConformingType::volume_dim;
     167             : 
     168           0 :     using primal_fields = typename ConformingType::primal_fields;
     169           0 :     using primal_fluxes = typename ConformingType::primal_fluxes;
     170             :     static_assert(FirstOrderSystem_detail::test_fields_and_fluxes<
     171             :                   volume_dim, primal_fields, primal_fluxes>::value);
     172             : 
     173           0 :     using background_fields = typename ConformingType::background_fields;
     174             :     static_assert(tt::is_a_v<tmpl::list, background_fields>);
     175           0 :     using inv_metric_tag = typename ConformingType::inv_metric_tag;
     176             :     static_assert(std::is_same_v<inv_metric_tag, void> or
     177             :                   tmpl::list_contains_v<background_fields, inv_metric_tag>);
     178             : 
     179           0 :     using fluxes_computer = typename ConformingType::fluxes_computer;
     180           0 :     using sources_computer = typename ConformingType::sources_computer;
     181           0 :     using fluxes_argument_tags = typename fluxes_computer::argument_tags;
     182             :     static_assert(tt::is_a_v<tmpl::list, fluxes_argument_tags>);
     183             : 
     184           0 :     using boundary_conditions_base =
     185             :         typename ConformingType::boundary_conditions_base;
     186             :     static_assert(
     187             :         std::is_base_of_v<domain::BoundaryConditions::BoundaryCondition,
     188             :                           boundary_conditions_base>);
     189             :   };
     190             : };
     191             : 
     192             : }  // namespace elliptic::protocols