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