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 :
8 : #include "DataStructures/Tensor/TypeAliases.hpp"
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 : * identity-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 1 : void metric_identity_det_jac_times_inv_jac(
179 : gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
180 : 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::ElementLogical, Frame::Inertial>&
185 : jacobian);
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>
236 1 : void metric_identity_jacobian_quantities(
237 : gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
238 : Frame::Inertial>*>
239 : det_jac_times_inverse_jacobian,
240 : gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
241 : Frame::Inertial>*>
242 : inverse_jacobian,
243 : gsl::not_null<
244 : Jacobian<DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
245 : jacobian,
246 : gsl::not_null<Scalar<DataVector>*> det_jacobian, const Mesh<Dim>& mesh,
247 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
248 : } // namespace dg
|
