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 partial derivatives of the `primal_fields`
153 : * 5. The `primal_fluxes`
154 : *
155 : * The function is expected to _add_ the sources \f$S_\alpha\f$ to the
156 : * output buffers. It must also have the following alias:
157 : *
158 : * - `const_global_cache_tags`: the subset of `argument_tags` that can be
159 : * retrieved from _any_ element's DataBox, because they are stored in the
160 : * global cache.
161 : *
162 : * The `sources_computer` may also be `void`, in which case \f$S_\alpha=0\f$.
163 : *
164 : * - `boundary_conditions_base`: A base class representing the supported
165 : * boundary conditions. Boundary conditions can be factory-created from this
166 : * base class. Currently this should be a specialization of
167 : * `elliptic::BoundaryConditions::BoundaryCondition`.
168 : *
169 : * - `modify_boundary_data` (optional, can be `void`): A class that can modify
170 : * boundary data received from a neighboring element, e.g. to transform from
171 : * one variable to another across element boundaries. This can be used to
172 : * solve for different variables in some regions of the domain to absorb
173 : * singularities. For example, when solving $-\Delta u = f$ we could define
174 : * $u=u_R + u_P$ in some region of the domain with a known (often singular)
175 : * field $u_P$ and solve only for the regular field $u_R$ in this region. So,
176 : * when receiving boundary data from outside this region, we subtract $u_P$,
177 : * and when receiving boundary data from inside this region, we add $u_P$. We
178 : * also add $\Delta u_P$ to the volume fixed sources $f$ inside the
179 : * regularized region.
180 : *
181 : * The `modify_boundary_data` must have an `argument_tags` type alias and an
182 : * `apply` function that takes these arguments in this order:
183 : *
184 : * 1. The primal fields and the normal-dot-fluxes on the mortar as not-null
185 : * pointer. These hold the received data from the neighbor and can be
186 : * modified. Note: the normal-dot-fluxes have been computed with the
187 : * neighbor's normal, so from the perspective of the receiving element they
188 : * are $-n_i F^i$ where $n_i$ is the face normal of the receiving element.
189 : * 2. The `const DirectionalId<Dim>& mortar_id` identifying the mortar.
190 : * 3. The `argument_tags`.
191 : *
192 : * The `modify_boundary_data` may also be used for modifications that depend
193 : * linearly on the variables. In this case, the class must also provide an
194 : * `argument_tags_linearized` type alias and an `apply_linearized` function
195 : * that takes these arguments in this order:
196 : *
197 : * 1. The remote primal fields and the remote normal-dot-fluxes on the mortar
198 : * as not-null pointers. These hold the received data from the neighbor and
199 : * can be modified in-place.
200 : * 2. The local primal fields and the local normal-dot-fluxes on the mortar
201 : * as const references. These hold the data from the local side of the
202 : * boundary.
203 : * 3. The `const DirectionalId<Dim>& mortar_id` identifying the mortar.
204 : * 4. The `argument_tags_linearized`.
205 : */
206 1 : struct FirstOrderSystem {
207 : template <typename ConformingType>
208 0 : struct test {
209 0 : static constexpr size_t volume_dim = ConformingType::volume_dim;
210 :
211 0 : using primal_fields = typename ConformingType::primal_fields;
212 0 : using primal_fluxes = typename ConformingType::primal_fluxes;
213 : static_assert(FirstOrderSystem_detail::test_fields_and_fluxes<
214 : volume_dim, primal_fields, primal_fluxes>::value);
215 :
216 0 : using background_fields = typename ConformingType::background_fields;
217 : static_assert(tt::is_a_v<tmpl::list, background_fields>);
218 0 : using inv_metric_tag = typename ConformingType::inv_metric_tag;
219 : static_assert(std::is_same_v<inv_metric_tag, void> or
220 : tmpl::list_contains_v<background_fields, inv_metric_tag>);
221 :
222 0 : using fluxes_computer = typename ConformingType::fluxes_computer;
223 0 : using sources_computer = typename ConformingType::sources_computer;
224 : static_assert(
225 : tt::is_a_v<tmpl::list, typename fluxes_computer::argument_tags>);
226 :
227 0 : using boundary_conditions_base =
228 : typename ConformingType::boundary_conditions_base;
229 : static_assert(
230 : std::is_base_of_v<domain::BoundaryConditions::BoundaryCondition,
231 : boundary_conditions_base>);
232 0 : using modify_boundary_data = typename ConformingType::modify_boundary_data;
233 : };
234 : };
235 :
236 : } // namespace elliptic::protocols
|
