Weighting.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cstddef>
8 #include <unordered_set>
9 
10 #include "DataStructures/DataVector.hpp"
13 #include "Utilities/Gsl.hpp"
14 
15 namespace LinearSolver::Schwarz {
16 
17 /*!
18  * \brief Weights for the solution on an element-centered subdomain, decreasing
19  * from 1 to 0.5 towards the `side` over the logical distance `width`, and
20  * further to 0 over the same distance outside the element.
21  *
22  * The weighting function over a full element-centered subdomain is
23  *
24  * \f{equation}
25  * w(\xi) = \frac{1}{2}\left( \phi\left( \frac{\xi + 1}{\delta} \right) -
26  * \phi\left( \frac{\xi - 1}{\delta} \right) \right) \f}
27  *
28  * where \f$\phi(\xi)\f$ is a second-order `::smoothstep`, i.e. the quintic
29  * polynomial
30  *
31  * \f{align*}
32  * \phi(\xi) = \begin{cases} \mathrm{sign}(\xi) \quad \text{for}
33  * \quad |\xi| > 1 \\
34  * \frac{1}{8}\left(15\xi - 10\xi^3 + 3\xi^5\right) \end{cases}
35  * \f}
36  *
37  * (see Eq. (39) in \cite Vincent2019qpd).
38  *
39  * The `LinearSolver::Schwarz::extruding_weight` and
40  * `LinearSolver::Schwarz::intruding_weight` functions each compute one of the
41  * two terms in \f$w(\xi)\f$. For example, consider an element-centered
42  * subdomain `A` that overlaps with a neighboring element-centered subdomain
43  * `B`. To combine solutions on `A` and `B` to a weighted solution on `A`,
44  * multiply the solution on `A` with the `extruding_weight` and the solution on
45  * `B` with the `intruding_weight`, both evaluated at the logical coordinates in
46  * `A` and at the `side` of `A` that faces `B`.
47  */
48 DataVector extruding_weight(const DataVector& logical_coords, double width,
49  const Side& side) noexcept;
50 
51 // @{
52 /*!
53  * \brief Weights for data on the central element of an element-centered
54  * subdomain
55  *
56  * Constructs the weighting field
57  *
58  * \f{equation}
59  * W(\boldsymbol{\xi}) = \prod^d_{i=0} w(\xi^i)
60  * \f}
61  *
62  * where \f$w(\xi^i)\f$ is the one-dimensional weighting function described in
63  * `LinearSolver::Schwarz::extruding_weight` and \f$\xi^i\f$ are the
64  * element-logical coordinates (see Eq. (41) in \cite Vincent2019qpd).
65  */
66 template <size_t Dim>
67 void element_weight(
68  gsl::not_null<Scalar<DataVector>*> element_weight,
69  const tnsr::I<DataVector, Dim, Frame::Logical>& logical_coords,
70  const std::array<double, Dim>& overlap_widths,
71  const std::unordered_set<Direction<Dim>>& external_boundaries) noexcept;
72 
73 template <size_t Dim>
74 Scalar<DataVector> element_weight(
75  const tnsr::I<DataVector, Dim, Frame::Logical>& logical_coords,
76  const std::array<double, Dim>& overlap_widths,
77  const std::unordered_set<Direction<Dim>>& external_boundaries) noexcept;
78 // @}
79 
80 /*!
81  * \brief Weights for the intruding solution of a neighboring element-centered
82  * subdomain, increasing from 0 to 0.5 towards the `side` over the logical
83  * distance `width`, and further to 1 over the same distance outside the
84  * element.
85  *
86  * \see `LinearSolver::Schwarz::extruding_weight`
87  */
88 DataVector intruding_weight(const DataVector& logical_coords, double width,
89  const Side& side) noexcept;
90 
91 // @{
92 /*!
93  * \brief Weights for data on overlap regions intruding into an element-centered
94  * subdomain
95  *
96  * Constructs the weighting field \f$W(\xi)\f$ as described in
97  * `LinearSolver::Schwarz::element_weight` for the data that overlaps with the
98  * central element of an element-centered subdomain. The weights are constructed
99  * in such a way that all weights at a grid point sum to one, i.e. the weight is
100  * conserved. The `logical_coords` are the element-logical coordinates of the
101  * central element.
102  *
103  * This function assumes that corner- and edge-neighbors of the central element
104  * are not part of the subdomain, which means that no contributions from those
105  * neighbors are expected although the weighting field is non-zero in overlap
106  * regions with those neighbors. Therefore, to retain conservation we must
107  * account for this missing weight by adding it to the central element, to the
108  * intruding overlaps from face-neighbors, or split it between the two. We
109  * choose to add the weight to the intruding overlaps, since that's where
110  * information from the corner- and edge-regions propagates through in a DG
111  * context.
112  */
113 template <size_t Dim>
114 void intruding_weight(
116  const tnsr::I<DataVector, Dim, Frame::Logical>& logical_coords,
117  const Direction<Dim>& direction,
118  const std::array<double, Dim>& overlap_widths,
119  size_t num_intruding_overlaps,
120  const std::unordered_set<Direction<Dim>>& external_boundaries) noexcept;
121 
122 template <size_t Dim>
124  const tnsr::I<DataVector, Dim, Frame::Logical>& logical_coords,
125  const Direction<Dim>& direction,
126  const std::array<double, Dim>& overlap_widths,
127  size_t num_intruding_overlaps,
128  const std::unordered_set<Direction<Dim>>& external_boundaries) noexcept;
129 // @}
130 
131 } // namespace LinearSolver::Schwarz
LinearSolver::Schwarz::intruding_weight
DataVector intruding_weight(const DataVector &logical_coords, const double width, const Side &side) noexcept
Weights for the intruding solution of a neighboring element-centered subdomain, increasing from 0 to ...
Definition: Weighting.cpp:79
LinearSolver::Schwarz::extruding_weight
DataVector extruding_weight(const DataVector &logical_coords, const double width, const Side &side) noexcept
Weights for the solution on an element-centered subdomain, decreasing from 1 to 0....
Definition: Weighting.cpp:25
unordered_set
Direction
Definition: Direction.hpp:23
LinearSolver::Schwarz
Items related to the Schwarz linear solver.
Definition: CommunicateOverlapFields.hpp:36
cstddef
array
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
Side
Side
Definition: Side.hpp:17
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
Tensor.hpp
Direction.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183