MetricIdentityJacobian.hpp
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>
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
187 /*!
188  * \ingroup DiscontinuousGalerkinGroup
189  * \brief Compute the Jacobian, inverse Jacobian, and determinant of the
190  * Jacobian so that they satisfy the metric identities.
191  *
192  * Uses dg::metric_identity_jacobian() to compute the determinant of the
193  * Jacobian times the inverse Jacobian. By taking the determinant of this
194  * product, we can isolate \f$J\f$, the determinant of the Jacobian. In \f$d\f$
195  * dimensions, we have:
196  *
197  * \f{align}{
198  *
199  * \mathrm{det}\left(J\frac{\partial \xi^{\hat{\imath}}}{\partial x^i}\right)
200  * = J^{d-1}.
201  *
202  * \f}
203  *
204  * We assume the determinant of the Jacobian is positive, which means logical
205  * and inertial coordinates have the same handedness. With this assumption, we
206  * have
207  *
208  * \f{align}{
209  *
210  * J = \sqrt[(d-1)]{\mathrm{det}\left(J\frac{\partial
211  * \xi^{\hat{\imath}}}{\partial x^i}\right)}
212  *
213  * \f}
214  *
215  * We can now compute the inverse Jacobian using:
216  *
217  * \f{align}{
218  *
219  * \frac{\partial \xi^{\hat{\imath}}}{\partial x^i}=
220  * \frac{1}{J}\left(J\frac{\partial \xi^{\hat{\imath}}}{\partial x^i}\right)
221  *
222  * \f}
223  *
224  * This guarantees that multiplying the determinant of the Jacobian by the
225  * inverse Jacobian gives a result that satisfies the metric identities. We also
226  * compute the Jacobian by inverting the inverse Jacobian, which guarantees they
227  * are (numerical) inverses of each other.
228  *
229  * \warning on entry jacobian must be the Jacobian to use for computing the
230  * determinant of the Jacobian times the inverse Jacobian so that it satisfies
231  * the metric identities. The jacobian can be computed analytically or
232  * numerically, either is fine. On output the jacobian is the inverse of the
233  * inverse Jacobian that satisfies the metric identities.
234  */
235 template <size_t Dim>
238  InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>*>
239  det_jac_times_inverse_jacobian,
241  InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>*>
242  inverse_jacobian,
243  gsl::not_null<Jacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>*>
244  jacobian,
245  gsl::not_null<Scalar<DataVector>*> det_jacobian, const Mesh<Dim>& mesh,
246  const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords) noexcept;
247 } // namespace dg
dg
Functionality related to discontinuous Galerkin schemes.
Definition: ApplyMassMatrix.hpp:13
cstddef
dg::metric_identity_det_jac_times_inv_jac
void metric_identity_det_jac_times_inv_jac(gsl::not_null< InverseJacobian< DataVector, Dim, Frame::Logical, Frame::Inertial > * > det_jac_times_inverse_jacobian, const Mesh< Dim > &mesh, const tnsr::I< DataVector, Dim, Frame::Inertial > &inertial_coords, const Jacobian< DataVector, Dim, Frame::Logical, Frame::Inertial > &jacobian) noexcept
Compute the Jacobian determinant times the inverse Jacobian so that the result is divergence-free.
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
TypeAliases.hpp
gsl
Implementations from the Guideline Support Library.