Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : 
       9             : #include "DataStructures/Tensor/TypeAliases.hpp"
      10             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      11             : #include "Utilities/Gsl.hpp"
      12             : 
      13             : /// \cond
      14             : class DataVector;
      15             : template <size_t Dim, typename T>
      16             : class DirectionMap;
      17             : template <size_t Dim>
      18             : class Mesh;
      19             : template <typename TagsList>
      20             : class Variables;
      21             : /// \endcond
      22             : 
      23             : /*!
      24             : * \brief The finite difference second derivatives originally
      25             : * designed for implementing the CCZ4 system.
      26             : */
      27             : namespace fd {
      28             : /*!
      29             :  * \brief Compute the pure and mixed second logical partial derivatives using
      30             :  * finite difference derivatives. Only 3 dimensions 4th-order fd is supported.
      31             :  *
      32             :  * We use a symmetric stencil in the bulk and slightly asymetric stencil
      33             :  * for the mixed partial derivatives where ghost data from two
      34             :  * neighbors is needed (to avoid communication across diagonal elements).
      35             :  * The `mixed_second_logical_derivatives` return an 3D array corresponding to
      36             :  * xy, yz, xz mixed partial derivatives.
      37             :  *
      38             :  * The stencil for the pure second derivatives is
      39             :  \f[ \frac{\partial^2 f}{\partial x^2} =
      40             :     - \frac{1}{12\delta^2}[f(x+2\delta,y)+f(x-2\delta,y)]
      41             :     + \frac{4}{3\delta^2}[f(x+\delta,y)+f(x-\delta,y)]-\frac{5}{2\delta^2}f(x,y)
      42             :  + O(\delta^4). \f]
      43             :  *
      44             :  * The stencil for the mixed second derivatives in the bulk is
      45             :  \f[ \begin{aligned} \frac{\partial^2 f}{\partial x \partial y} =&
      46             :       \frac{1}{3\delta^2}[f(x+\delta,y+\delta) + f(x-\delta,y-\delta) -
      47             :  f(x+\delta,y-\delta) - f(x-\delta,y+\delta)] \\
      48             :     &+ \frac{1}{48\delta^2}[f(x+2\delta,y-2\delta) + f(x-2\delta,y+2\delta)
      49             :         - f(x+2\delta,y+2\delta) - f(x-2\delta,y-2\delta)] + O(\delta^4).
      50             :  \end{aligned}\f]
      51             :  *
      52             :  * The stencil for the mixed second derivatives in the ghost zone is (with
      53             :  * descending diagonal)
      54             :  \f[ \begin{aligned} \frac{\partial^2 f}{\partial x \partial y} =&
      55             :       \frac{1}{24\delta^2}[f(x+2\delta,y-2\delta) + f(x-2\delta,y+2\delta)
      56             :         - f(x+2\delta,y) - f(x-2\delta,y) - f(x,y+2\delta) - f(x,y-2\delta)] \\
      57             :     &+ \frac{2}{3\delta^2}[f(x+\delta,y) + f(x-\delta,y) + f(x,y+\delta) +
      58             :         f(x,y-\delta) - f(x+\delta,y-\delta) - f(x-\delta,y+\delta)]
      59             :     - \frac{5}{4\delta^2}f(x,y) + O(\delta^4),
      60             :  \end{aligned}\f]
      61             :  *
      62             :  * and the ascending diagonal stencil has the opposite weights.
      63             :  *
      64             :  * \note Currently the stride is always one because we transpose the data before
      65             :  * reconstruction.
      66             :  */
      67             : template <size_t Dim>
      68           1 : void second_logical_partial_derivatives(
      69             :     gsl::not_null<std::array<gsl::span<double>, Dim>*>
      70             :         pure_second_logical_derivatives,
      71             :     gsl::not_null<std::array<gsl::span<double>, Dim>*>
      72             :         mixed_second_logical_derivatives,
      73             :     const gsl::span<const double>& volume_vars,
      74             :     const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
      75             :     const Mesh<Dim>& volume_mesh, size_t number_of_variables, size_t fd_order);
      76             : 
      77             : /*!
      78             :  * \brief Compute the second partial derivatives on the `DerivativeFrame` using
      79             :  * the `inverse_jacobian` and `inverse_hessian`.
      80             :  * Only 3 dimensions 4th-order fd is supported.
      81             :  *
      82             :  * First and second logical partial derivatives are first computed using the
      83             :  * `fd::logical_partial_derivatives()` and
      84             :  * `fd::second_logical_partial_derivatives()` function.
      85             :  *
      86             :  * \note The `inverse_hessian` has not been implemented, so currently inertial
      87             :  * coordinates cannot mix logical coordinates.
      88             :  */
      89             : template <typename DerivativeTags, size_t Dim, typename DerivativeFrame>
      90           1 : void second_partial_derivatives(
      91             :     gsl::not_null<
      92             :         Variables<db::wrap_tags_in<::Tags::second_deriv, DerivativeTags,
      93             :                                    tmpl::size_t<Dim>, DerivativeFrame>>*>
      94             :         second_partial_derivatives,
      95             :     const gsl::span<const double>& volume_vars,
      96             :     const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
      97             :     const Mesh<Dim>& volume_mesh, size_t number_of_variables,
      98             :     size_t fd_order,
      99             :     const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     100             :                           DerivativeFrame>& inverse_jacobian);
     101             : }  // namespace fd