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             : #include <utility>
       8             : 
       9             : #include "DataStructures/DataVector.hpp"
      10             : #include "DataStructures/Tensor/Tensor.hpp"
      11             : #include "DataStructures/Variables.hpp"
      12             : #include "Domain/Structure/Direction.hpp"
      13             : #include "Domain/Structure/Side.hpp"
      14             : #include "NumericalAlgorithms/Spectral/BoundaryInterpolationMatrices.hpp"
      15             : #include "NumericalAlgorithms/Spectral/BoundaryLiftingTerm.hpp"
      16             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      17             : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
      18             : #include "Utilities/Algorithm.hpp"
      19             : #include "Utilities/Gsl.hpp"
      20             : 
      21             : namespace dg {
      22             : namespace detail {
      23             : template <typename ValueType, size_t Dim>
      24             : void lift_boundary_terms_gauss_points_impl(
      25             :     gsl::not_null<ValueType*> volume_dt_vars, size_t num_independent_components,
      26             :     const Mesh<Dim>& volume_mesh, size_t dimension,
      27             :     const Scalar<DataVector>& volume_det_inv_jacobian, size_t num_boundary_pts,
      28             :     const gsl::span<const ValueType>& boundary_corrections,
      29             :     const DataVector& boundary_lifting_term,
      30             :     const Scalar<DataVector>& magnitude_of_face_normal,
      31             :     const Scalar<DataVector>& face_det_jacobian);
      32             : 
      33             : template <typename ValueType, size_t Dim>
      34             : void lift_boundary_terms_gauss_points_impl(
      35             :     gsl::not_null<ValueType*> volume_dt_vars, size_t num_independent_components,
      36             :     const Mesh<Dim>& volume_mesh, size_t dimension,
      37             :     const Scalar<DataVector>& volume_det_inv_jacobian, size_t num_boundary_pts,
      38             :     const gsl::span<const ValueType>& upper_boundary_corrections,
      39             :     const DataVector& upper_boundary_lifting_term,
      40             :     const Scalar<DataVector>& upper_magnitude_of_face_normal,
      41             :     const Scalar<DataVector>& upper_face_det_jacobian,
      42             :     const gsl::span<const ValueType>& lower_boundary_corrections,
      43             :     const DataVector& lower_boundary_lifting_term,
      44             :     const Scalar<DataVector>& lower_magnitude_of_face_normal,
      45             :     const Scalar<DataVector>& lower_face_det_jacobian);
      46             : 
      47             : template <typename ValueType, size_t Dim>
      48             : void lift_boundary_terms_gauss_points_impl(
      49             :     gsl::not_null<ValueType*> volume_data, size_t num_independent_components,
      50             :     const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
      51             :     const gsl::span<const ValueType>& boundary_corrections);
      52             : }  // namespace detail
      53             : 
      54             : /*!
      55             :  * \brief Lift the boundary corrections to the volume time derivatives in the
      56             :  * specified direction.
      57             :  *
      58             :  * The general lifting term (for the \f$\xi\f$-dimension) is:
      59             :  *
      60             :  * \f{align*}{
      61             :  * \partial_t u_{\alpha\breve{\imath}\breve{\jmath}\breve{k}}=\cdots
      62             :  * -\frac{\ell_{\breve{\imath}}\left(\xi=1\right)}
      63             :  *   {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
      64             :  *   \left[J\sqrt{
      65             :  *   \frac{\partial\xi}{\partial x^i} \gamma^{ij}
      66             :  *   \frac{\partial\xi}{\partial x^j}}
      67             :  *   \left(G_{\alpha} + D_{\alpha}\right)
      68             :  *   \right]_{\breve{\jmath}\breve{k}}\left(\xi=1\right),
      69             :  * \f}
      70             :  *
      71             :  * where \f$\breve{\imath}\f$, \f$\breve{\jmath}\f$, and \f$\breve{k}\f$ are
      72             :  * indices in the logical \f$\xi\f$, \f$\eta\f$, and \f$\zeta\f$ dimensions.
      73             :  * The \f$G+D\f$ terms correspond to the `boundary_corrections` function
      74             :  * argument, and the function Spectral::boundary_lifting_term() is used to
      75             :  * compute and cache the terms from the lifting terms
      76             :  * \f$\ell_{\breve{\imath}}(\xi=\pm1)/w_{\breve{\imath}}\f$.
      77             :  *
      78             :  * \note that normal vectors are pointing out of the element.
      79             :  *
      80             :  * \param dt_vars The volume time derivatives
      81             :  * $\partial_t u_{\alpha\breve{\imath}\breve{\jmath}\breve{k}}$ to be updated
      82             :  * \param volume_det_inv_jacobian The inverse Jacobian determinant in the volume
      83             :  * $J^{-1}_{\breve{\imath}\breve{\jmath}\breve{k}}$
      84             :  * \param volume_mesh The mesh of the volume
      85             :  * \param direction The direction of the face on which the boundary corrections
      86             :  * are defined
      87             :  * \param boundary_corrections The boundary corrections
      88             :  * $\left(G_{\alpha} + D_{\alpha}\right)$
      89             :  * \param magnitude_of_face_normal The term $\sqrt{
      90             :  * \frac{\partial\xi}{\partial x^i} \gamma^{ij}
      91             :  * \frac{\partial\xi}{\partial x^j}}$
      92             :  * \param face_det_jacobian The volume Jacobian determinant $J$ evaluated on the
      93             :  * face (not the surface Jacobian determinant)
      94             :  */
      95             : template <size_t Dim, typename DtTagsList, typename BoundaryCorrectionTagsList>
      96           1 : void lift_boundary_terms_gauss_points(
      97             :     const gsl::not_null<Variables<DtTagsList>*> dt_vars,
      98             :     const Scalar<DataVector>& volume_det_inv_jacobian,
      99             :     const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
     100             :     const Variables<BoundaryCorrectionTagsList>& boundary_corrections,
     101             :     const Scalar<DataVector>& magnitude_of_face_normal,
     102             :     const Scalar<DataVector>& face_det_jacobian) {
     103             :   ASSERT(alg::all_of(volume_mesh.quadrature(),
     104             :                      [](const Spectral::Quadrature quadrature) {
     105             :                        return quadrature == Spectral::Quadrature::Gauss;
     106             :                      }),
     107             :          "Must use Gauss points in all directions but got the mesh: "
     108             :              << volume_mesh);
     109             :   const Mesh<Dim - 1> boundary_mesh =
     110             :       volume_mesh.slice_away(direction.dimension());
     111             :   const Mesh<1> volume_stripe_mesh =
     112             :       volume_mesh.slice_through(direction.dimension());
     113             :   const size_t num_boundary_grid_points = boundary_mesh.number_of_grid_points();
     114             :   detail::lift_boundary_terms_gauss_points_impl(
     115             :       make_not_null(dt_vars->data()), dt_vars->number_of_independent_components,
     116             :       volume_mesh, direction.dimension(), volume_det_inv_jacobian,
     117             :       num_boundary_grid_points,
     118             :       gsl::make_span(boundary_corrections.data(), boundary_corrections.size()),
     119             :       direction.side() == Side::Upper
     120             :           ? Spectral::boundary_lifting_term(volume_stripe_mesh).second
     121             :           : Spectral::boundary_lifting_term(volume_stripe_mesh).first,
     122             :       magnitude_of_face_normal, face_det_jacobian);
     123             : }
     124             : 
     125             : /*!
     126             :  * \brief Lift both the upper and lower (in logical coordinates) boundary
     127             :  * corrections to the volume time derivatives in the specified logical
     128             :  * dimension.
     129             :  *
     130             :  * The upper and lower boundary corrections in the logical `dimension` are
     131             :  * lifted together in order to reduce the amount of striding through data that
     132             :  * is needed and to improve cache-friendliness.
     133             :  *
     134             :  * The general lifting term (for the \f$\xi\f$-dimension) is:
     135             :  *
     136             :  * \f{align*}{
     137             :  * \partial_t u_{\alpha\breve{\imath}\breve{\jmath}\breve{k}}=\cdots
     138             :  * -\frac{\ell_{\breve{\imath}}\left(\xi=1\right)}
     139             :  *   {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
     140             :  *   \left[J\sqrt{
     141             :  *   \frac{\partial\xi}{\partial x^i} \gamma^{ij}
     142             :  *   \frac{\partial\xi}{\partial x^j}}
     143             :  *   \left(G_{\alpha} + D_{\alpha}\right)
     144             :  *   \right]_{\breve{\jmath}\breve{k}}\left(\xi=1\right)
     145             :  * - \frac{\ell_{\breve{\imath}}\left(\xi=-1\right)}
     146             :  *   {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
     147             :  *   \left[J\sqrt{
     148             :  *   \frac{\partial\xi}{\partial x^i} \gamma^{ij}
     149             :  *   \frac{\partial\xi}{\partial x^j}}
     150             :  *   \left(G_{\alpha} + D_{\alpha}\right)
     151             :  *   \right]_{\breve{\jmath}\breve{k}}\left(\xi=-1\right),
     152             :  * \f}
     153             :  *
     154             :  * where \f$\breve{\imath}\f$, \f$\breve{\jmath}\f$, and \f$\breve{k}\f$ are
     155             :  * indices in the logical \f$\xi\f$, \f$\eta\f$, and \f$\zeta\f$ dimensions.
     156             :  * The \f$G+D\f$ terms correspond to the `upper_boundary_corrections` and
     157             :  * `lower_boundary_corrections` function arguments, and the function
     158             :  * Spectral::boundary_lifting_term() is used to compute and cache the terms
     159             :  * from the lifting terms
     160             :  * \f$\ell_{\breve{\imath}}(\xi=\pm1)/w_{\breve{\imath}}\f$.
     161             :  *
     162             :  * \note that normal vectors are pointing out of the element and therefore both
     163             :  * terms have the same sign.
     164             :  */
     165             : template <size_t Dim, typename DtTagsList, typename BoundaryCorrectionTagsList>
     166           1 : void lift_boundary_terms_gauss_points(
     167             :     const gsl::not_null<Variables<DtTagsList>*> dt_vars,
     168             :     const Scalar<DataVector>& volume_det_inv_jacobian,
     169             :     const Mesh<Dim>& volume_mesh, const size_t dimension,
     170             :     const Variables<BoundaryCorrectionTagsList>& upper_boundary_corrections,
     171             :     const Scalar<DataVector>& upper_magnitude_of_face_normal,
     172             :     const Scalar<DataVector>& upper_face_det_jacobian,
     173             :     const Variables<BoundaryCorrectionTagsList>& lower_boundary_corrections,
     174             :     const Scalar<DataVector>& lower_magnitude_of_face_normal,
     175             :     const Scalar<DataVector>& lower_face_det_jacobian) {
     176             :   ASSERT(alg::all_of(volume_mesh.quadrature(),
     177             :                      [](const Spectral::Quadrature quadrature) {
     178             :                        return quadrature == Spectral::Quadrature::Gauss;
     179             :                      }),
     180             :          "Must use Gauss points in all directions but got the mesh: "
     181             :              << volume_mesh);
     182             :   const Mesh<Dim - 1> boundary_mesh = volume_mesh.slice_away(dimension);
     183             :   const Mesh<1> volume_stripe_mesh = volume_mesh.slice_through(dimension);
     184             :   const size_t num_boundary_grid_points = boundary_mesh.number_of_grid_points();
     185             :   detail::lift_boundary_terms_gauss_points_impl(
     186             :       make_not_null(dt_vars->data()), dt_vars->number_of_independent_components,
     187             :       volume_mesh, dimension, volume_det_inv_jacobian, num_boundary_grid_points,
     188             :       gsl::make_span(upper_boundary_corrections.data(),
     189             :                      upper_boundary_corrections.size()),
     190             :       Spectral::boundary_lifting_term(volume_stripe_mesh).second,
     191             :       upper_magnitude_of_face_normal, upper_face_det_jacobian,
     192             :       gsl::make_span(lower_boundary_corrections.data(),
     193             :                      lower_boundary_corrections.size()),
     194             :       Spectral::boundary_lifting_term(volume_stripe_mesh).first,
     195             :       lower_magnitude_of_face_normal, lower_face_det_jacobian);
     196             : }
     197             : 
     198             : /*!
     199             :  * \brief Lift the boundary corrections from the face to the volume for Gauss
     200             :  * points in the direction perpendicular to the face
     201             :  *
     202             :  * This function doesn't apply any Jacobians or quadrature weights. It only
     203             :  * applies the lifting matrix $\ell_{\breve{\imath}}(\xi=\pm1)$ to the
     204             :  * `boundary_corrections` and adds the result to the `volume_data`.
     205             :  */
     206             : template <size_t Dim, typename... VolumeTags,
     207             :           typename... BoundaryCorrectionTags>
     208           1 : void lift_boundary_terms_gauss_points(
     209             :     const gsl::not_null<Variables<tmpl::list<VolumeTags...>>*> volume_data,
     210             :     const Variables<tmpl::list<BoundaryCorrectionTags...>>&
     211             :         boundary_corrections,
     212             :     const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
     213             :   detail::lift_boundary_terms_gauss_points_impl(
     214             :       make_not_null(volume_data->data()),
     215             :       volume_data->number_of_independent_components, volume_mesh, direction,
     216             :       gsl::make_span(boundary_corrections.data(), boundary_corrections.size()));
     217             : }
     218             : }  // namespace dg