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