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