MetricIdentityJacobian.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 
9 
10 /// \cond
11 class DataVector;
12 template <size_t Dim>
13 class Mesh;
14 namespace gsl {
15 template <typename T>
16 class not_null;
17 } // namespace gsl
18 /// \endcond
19 
20 namespace dg {
21 /*!
22  * \ingroup DiscontinuousGalerkinGroup
23  * \brief Compute the Jacobian determinant times the inverse Jacobian so that
24  * the result is divergence-free.
25  *
26  * The metric identities are given by
27  *
28  * \f{align*}{
29  *
30  * \partial_{\hat{\imath}}\left(J\frac{\partial\xi^{\hat{\imath}}}
31  * {\partial x^i}\right)=0.
32  *
33  * \f}
34  *
35  * We want to compute \f$J\partial\xi^{\hat{\imath}}/\partial x^i\f$ in such a
36  * way that the above metric identity is satisfied numerically/discretely. We
37  * refer to the inverse Jacobian computed this way as the "metric
38  * indentity-satisfying inverse Jacobian".
39  *
40  * The discretized form, with the Jacobian determinant \f$J\f$ expanded, is
41  * given by
42  *
43  * \f{align*}{
44  *
45  * 2\left(J\frac{\partial \xi^{\hat{\imath}}}{\partial x^i}\right)_{s}
46  * &=\epsilon_{ijk}
47  * \sum_{t_{\hat{1}}} \epsilon^{\hat{\imath}\hat{1}\hat{k}}
48  * D^{(\hat{1})}_{s t_1}
49  * \left(x^j\frac{\partial x^k}{\partial \xi^{\hat{k}}}\right)_{t}
50  * \notag \\
51  * &+\epsilon_{ijk}
52  * \sum_{t_{\hat{2}}} \epsilon^{\hat{\imath}\hat{2}\hat{k}}
53  * D^{(\hat{2})}_{st_2}
54  * \left(x^j\frac{\partial x^k}{\partial \xi^{\hat{k}}}\right)_{t}
55  * \notag \\
56  * &+\epsilon_{ijk}
57  * \sum_{t_{\hat{3}}}\epsilon^{\hat{\imath}\hat{3}\hat{k}}
58  * D^{(\hat{3})}_{st_3}
59  * \left(x^j\frac{\partial x^k}{\partial \xi^{\hat{k}}}\right)_{t}.
60  *
61  * \f}
62  *
63  * where indices \f$s,t,t_1,t_2\f$ and \f$t_3\f$ are over grid points. \f$t_i\f$
64  * are the grid points in the particular logical direction.
65  *
66  * In 1d we have:
67  *
68  * \f{align*}{
69  *
70  * J\frac{\partial \xi^{\hat{\imath}}}{\partial x^i}=
71  * = J\frac{\partial \xi^{\hat{i}}}{\partial x^i} = 1
72  *
73  * \f}
74  *
75  * In 2d we have:
76  *
77  * \f{align*}{
78  *
79  * J\frac{\partial \xi^{\hat{1}}}{\partial x^1}&=D^{\hat{(2)}}x^2 &
80  * J\frac{\partial \xi^{\hat{2}}}{\partial x^1}&=-D^{\hat{(1)}}x^2 \\
81  * J\frac{\partial \xi^{\hat{1}}}{\partial x^2}&=-D^{\hat{(2)}}x^1 &
82  * J\frac{\partial \xi^{\hat{2}}}{\partial x^2}&=D^{\hat{(1)}}x^1 \\
83  *
84  * \f}
85  *
86  * In 3d we have:
87  *
88  * \f{align*}{
89  *
90  * 2J\frac{\partial \xi^{\hat{1}}}{\partial x^1}&=
91  * D^{(\hat{2})}\left(x^2\frac{\partial x^3}{\partial \xi^{\hat{3}}}\right)
92  * -D^{(\hat{2})}\left(x^3\frac{\partial x^2}{\partial \xi^{\hat{3}}}\right)
93  * +D^{(\hat{3})}\left(x^3\frac{\partial x^2}{\partial \xi^{\hat{2}}}\right)
94  * -D^{(\hat{3})}\left(x^2\frac{\partial x^3}{\partial \xi^{\hat{2}}}\right)\\
95  * 2J\frac{\partial \xi^{\hat{2}}}{\partial x^1}&=
96  * D^{(\hat{1})}\left(x^3\frac{\partial x^2}{\partial\xi^{\hat{3}}}\right)
97  * -D^{(\hat{1})}\left(x^2\frac{\partial x^3}{\partial\xi^{\hat{3}}}\right)
98  * +D^{(\hat{3})}\left(x^2\frac{\partial x^3}{\partial\xi^{\hat{1}}}\right)
99  * -D^{(\hat{3})}\left(x^3\frac{\partial x^2}{\partial\xi^{\hat{1}}}\right) \\
100  * 2J\frac{\partial \xi^{\hat{3}}}{\partial x^1}&=
101  * D^{(\hat{1})}\left(x^2\frac{\partial x^3}{\partial\xi^{\hat{2}}}\right)
102  * -D^{(\hat{1})}\left(x^3\frac{\partial x^2}{\partial\xi^{\hat{2}}}\right)
103  * +D^{(\hat{2})}\left(x^3\frac{\partial x^2}{\partial\xi^{\hat{1}}}\right)
104  * -D^{(\hat{2})}\left(x^2\frac{\partial x^3}{\partial\xi^{\hat{1}}}\right)\\
105  * 2J\frac{\partial \xi^{\hat{1}}}{\partial x^2}&=
106  * D^{(\hat{2})}\left(x^3\frac{\partial x^1}{\partial \xi^{\hat{3}}}\right)
107  * -D^{(\hat{2})}\left(x^1\frac{\partial x^3}{\partial \xi^{\hat{3}}}\right)
108  * +D^{(\hat{3})}\left(x^1\frac{\partial x^3}{\partial \xi^{\hat{2}}}\right)
109  * -D^{(\hat{3})}\left(x^3\frac{\partial x^1}{\partial \xi^{\hat{2}}}\right) \\
110  * 2J\frac{\partial \xi^{\hat{2}}}{\partial x^2}&=
111  * D^{(\hat{1})}\left(x^1\frac{\partial x^3}{\partial \xi^{\hat{3}}}\right)
112  * -D^{(\hat{1})}\left(x^3\frac{\partial x^1}{\partial \xi^{\hat{3}}}\right)
113  * +D^{(\hat{3})}\left(x^3\frac{\partial x^1}{\partial \xi^{\hat{1}}}\right)
114  * -D^{(\hat{3})}\left(x^1\frac{\partial x^3}{\partial \xi^{\hat{1}}}\right)\\
115  * 2J\frac{\partial \xi^{\hat{3}}}{\partial x^2}&=
116  * D^{(\hat{1})}\left(x^3\frac{\partial x^1}{\partial \xi^{\hat{2}}}\right)
117  * -D^{(\hat{1})}\left(x^1\frac{\partial x^3}{\partial \xi^{\hat{2}}}\right)
118  * +D^{(\hat{2})}\left(x^1\frac{\partial x^3}{\partial \xi^{\hat{1}}}\right)
119  * -D^{(\hat{2})}\left(x^3\frac{\partial x^1}{\partial \xi^{\hat{1}}}\right) \\
120  * 2J\frac{\partial \xi^{\hat{1}}}{\partial x^3}&=
121  * D^{(\hat{2})}\left(x^1\frac{\partial x^2}{\partial \xi^{\hat{3}}}\right)
122  * -D^{(\hat{2})}\left(x^2\frac{\partial x^1}{\partial \xi^{\hat{3}}}\right)
123  * +D^{(\hat{3})}\left(x^2\frac{\partial x^1}{\partial \xi^{\hat{2}}}\right)
124  * -D^{(\hat{3})}\left(x^1\frac{\partial x^2}{\partial \xi^{\hat{2}}}\right) \\
125  * 2J\frac{\partial \xi^{\hat{2}}}{\partial x^3}&=
126  * D^{(\hat{1})}\left(x^2\frac{\partial x^1}{\partial \xi^{\hat{3}}}\right)
127  * -D^{(\hat{1})}\left(x^1\frac{\partial x^2}{\partial \xi^{\hat{3}}}\right)
128  * +D^{(\hat{3})}\left(x^1\frac{\partial x^2}{\partial \xi^{\hat{1}}}\right)
129  * -D^{(\hat{3})}\left(x^2\frac{\partial x^1}{\partial \xi^{\hat{1}}}\right) \\
130  * 2J\frac{\partial \xi^{\hat{3}}}{\partial x^3}&=
131  * D^{(\hat{1})}\left(x^1\frac{\partial x^2}{\partial \xi^{\hat{2}}}\right)
132  * -D^{(\hat{1})}\left(x^2\frac{\partial x^1}{\partial \xi^{\hat{2}}}\right)
133  * +D^{(\hat{2})}\left(x^2\frac{\partial x^1}{\partial \xi^{\hat{1}}}\right)
134  * -D^{(\hat{2})}\left(x^1\frac{\partial x^2}{\partial \xi^{\hat{1}}}\right)
135  *
136  * \f}
137  *
138  * Again, this ensures that the metric identities are satisfied discretely. That
139  * is,
140  *
141  * \f{align*}{
142  *
143  * \partial_{\hat{\imath}}\left(J\frac{\partial\xi^{\hat{\imath}}}
144  * {\partial x^i}\right)=0
145  *
146  * \f}
147  *
148  * numerically.
149  *
150  * The reason for calculating \f$J\partial\xi^{\hat{\imath}}/\partial x^i\f$ in
151  * this manner is because in the weak form of DG (and most spectral-type methods
152  * can be recast as DG) we effectively evaluate
153  *
154  * \f{align*}{
155  *
156  * \partial_{\hat{\imath}}\left(J\frac{\partial\xi^{\hat{\imath}}}
157  * {\partial x^i} F^i\right),
158  *
159  * \f}
160  *
161  * which should be identically zero if \f$F^i\f$ is constant. This feature of a
162  * scheme is referred to as free-stream preserving. Note that another way to
163  * achieve free-stream preservation is to subtract off the metric identity
164  * violations. That is,
165  *
166  * \f{align*}{
167  *
168  * \partial_{\hat{\imath}}\left(J\frac{\partial\xi^{\hat{\imath}}}
169  * {\partial x^i} F^i\right) -
170  * F^i \partial_{\hat{\imath}}\left(J\frac{\partial\xi^{\hat{\imath}}}
171  * {\partial x^i}\right).
172  *
173  * \f}
174  *
175  * The subtraction technique is most commonly used in finite difference codes.
176  */
177 template <size_t Dim>
178 void metric_identity_jacobian(
180  InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>*>
181  det_jac_times_inverse_jacobian,
182  const Mesh<Dim>& mesh,
183  const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords,
184  const Jacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
185  jacobian) noexcept;
186 } // namespace dg
dg
Functionality related to discontinuous Galerkin schemes.
Definition: ComputeNonconservativeBoundaryFluxes.hpp:23
cstddef
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:47
TypeAliases.hpp
gsl
Implementations from the Guideline Support Library.
Definition: Gsl.hpp:80
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183