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 <string>
       8              : 
       9              : #include "DataStructures/Tensor/TypeAliases.hpp"
      10              : #include "Domain/Tags.hpp"
      11              : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      12              : #include "Utilities/Gsl.hpp"
      13              : #include "Utilities/TMPL.hpp"
      14              : 
      15              : /// \cond
      16              : class DataVector;
      17              : /// \endcond
      18              : 
      19              : namespace domain {
      20              : /*!
      21              :  * \ingroup ComputationalDomainGroup
      22              :  * \brief A diagnostic comparing the analytic and numerical Jacobians for a
      23              :  * map.
      24              :  *
      25              :  * Specifically, returns
      26              :  * \f[
      27              :  * C_{\hat{i}} = 1 -
      28              :  * \frac{\sum_i |\partial_{\hat{i}} x^i|}{\sum_i |D_{\hat{i}} x^i|}
      29              :  * \f], where \f$x^{\hat{i}}\f$ are the logical
      30              :  * coordinates, \f$x^i\f$ are the coordinates in the target frame,
      31              :  * \f$\partial_{\hat{i}}x^i\f$ is the analytic Jacobian, and \f$D_{\hat{i}}
      32              :  * x^i\f$ is the numerical Jacobian.
      33              :  *
      34              :  * \note This function accepts the transpose of the numeric Jacobian as a
      35              :  * parameter, since the numeric Jacobian will typically be computed via
      36              :  * logical_partial_derivative(), which prepends the logical (source frame)
      37              :  * derivative index. Tensors of type Jacobian, in contrast, have the derivative
      38              :  * index second.
      39              :  */
      40              : template <size_t Dim, typename Fr>
      41            1 : void jacobian_diagnostic(
      42              :     const gsl::not_null<
      43              :         tnsr::i<DataVector, Dim, typename Frame::ElementLogical>*>
      44              :         jacobian_diag,
      45              :     const Jacobian<DataVector, Dim, typename Frame::ElementLogical, Fr>&
      46              :         analytic_jacobian,
      47              :     const TensorMetafunctions::prepend_spatial_index<
      48              :         tnsr::I<DataVector, Dim, Fr>, Dim, UpLo::Lo,
      49              :         typename Frame::ElementLogical>& numeric_jacobian_transpose);
      50              : 
      51              : /// @{
      52              : /*!
      53              :  * \ingroup ComputationalDomainGroup
      54              :  * \brief A diagnostic comparing the analytic and numerical Jacobians for a
      55              :  * map.
      56              :  *
      57              :  * Specifically, returns
      58              :  * \f[
      59              :  * C_{\hat{i}} = 1 -
      60              :  * \frac{\sum_i |\partial_{\hat{i}} x^i|}{\sum_i |D_{\hat{i}} x^i|}
      61              :  * \f],
      62              :  * where \f$x^{\hat{i}}\f$ are the logical coordinates, \f$x^i\f$ are the
      63              :  * coordinates in the target frame, \f$\partial_{\hat{i}}x^i\f$ is the analytic
      64              :  * Jacobian, and \f$D_{\hat{i}} x^i\f$ is the numerical Jacobian.
      65              :  *
      66              :  * \note This function accepts the analytic jacobian, mapped coordinates, and
      67              :  * mesh as a parameter. The numeric jacobian is computed
      68              :  * internally by differentiating the mapped coordinates with respect to the
      69              :  * logical coordinates.
      70              :  */
      71              : template <size_t Dim, typename Fr>
      72            1 : void jacobian_diagnostic(
      73              :     const gsl::not_null<tnsr::i<DataVector, Dim, Frame::ElementLogical>*>
      74              :         jacobian_diag,
      75              :     const ::Jacobian<DataVector, Dim, Frame::ElementLogical, Fr>&
      76              :         analytic_jacobian,
      77              :     const tnsr::I<DataVector, Dim, Fr>& mapped_coords, const ::Mesh<Dim>& mesh);
      78              : 
      79              : template <size_t Dim, typename Fr>
      80            1 : tnsr::i<DataVector, Dim, Frame::ElementLogical> jacobian_diagnostic(
      81              :     const ::Jacobian<DataVector, Dim, Frame::ElementLogical, Fr>&
      82              :         analytic_jacobian,
      83              :     const tnsr::I<DataVector, Dim, Fr>& mapped_coords, const ::Mesh<Dim>& mesh);
      84              : /// @}
      85              : 
      86              : namespace Tags {
      87              : /// \ingroup DataBoxTagsGroup
      88              : /// \ingroup ComputationalDomainGroup
      89              : /// \brief A diagnostic comparing the analytic and numerical Jacobians for a
      90              : /// map. See `domain::jacobian_diagnostic` for details.
      91              : template <size_t Dim>
      92            1 : struct JacobianDiagnostic : db::SimpleTag {
      93            0 :   using type = tnsr::i<DataVector, Dim, typename Frame::ElementLogical>;
      94              : };
      95              : 
      96              : /// \ingroup DataBoxTagsGroup
      97              : /// \ingroup ComputationalDomainGroup
      98              : /// \brief Computes the Jacobian diagnostic, which compares the analytic
      99              : /// Jacobian (provided by some coordinate map) to a numerical Jacobian computed
     100              : /// using numerical partial derivatives. The coordinates must be in the target
     101              : /// frame of the map. See `domain::jacobian_diagnostic` for details of the
     102              : /// calculation.
     103              : template <size_t Dim, typename TargetFrame>
     104            1 : struct JacobianDiagnosticCompute : JacobianDiagnostic<Dim>, db::ComputeTag {
     105            0 :   using base = JacobianDiagnostic<Dim>;
     106            0 :   using return_type = typename base::type;
     107            0 :   using argument_tags = tmpl::list<
     108              :       ::domain::Tags::Jacobian<Dim, Frame::ElementLogical, TargetFrame>,
     109              :       ::domain::Tags::Coordinates<Dim, TargetFrame>, Mesh<Dim>>;
     110            0 :   static constexpr auto function = static_cast<void (*)(
     111              :       const gsl::not_null<
     112              :           tnsr::i<DataVector, Dim, typename Frame::ElementLogical>*>,
     113              :       const ::Jacobian<DataVector, Dim, Frame::ElementLogical, TargetFrame>&,
     114              :       const tnsr::I<DataVector, Dim, TargetFrame>&, const ::Mesh<Dim>&)>(
     115              :       &::domain::jacobian_diagnostic<Dim, TargetFrame>);
     116              : };
     117              : }  // namespace Tags
     118              : }  // namespace domain