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              : template <typename DerivativeTags, size_t Dim, typename DerivativeFrame>
      87            1 : void second_partial_derivatives(
      88              :     gsl::not_null<
      89              :         Variables<db::wrap_tags_in<::Tags::second_deriv, DerivativeTags,
      90              :                                    tmpl::size_t<Dim>, DerivativeFrame>>*>
      91              :         second_partial_derivatives,
      92              :     const gsl::span<const double>& volume_vars,
      93              :     const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
      94              :     const Mesh<Dim>& volume_mesh, size_t number_of_variables, size_t fd_order,
      95              :     const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
      96              :                           DerivativeFrame>& inverse_jacobian,
      97              :     const InverseHessian<DataVector, Dim, Frame::ElementLogical,
      98              :                          DerivativeFrame>& inverse_hessian);
      99              : }  // namespace fd