text: // Distributed under the MIT License.
  // See LICENSE.txt for details.
  
  #pragma once
  
  #include <array>
  #include <cstddef>
  
  #include "DataStructures/DataBox/Tag.hpp"
  #include "DataStructures/DataVector.hpp"
  #include "DataStructures/Tensor/Tensor.hpp"
  #include "Domain/FaceNormal.hpp"
  #include "Evolution/Systems/ScalarWave/Tags.hpp"
  #include "Utilities/Gsl.hpp"
  #include "Utilities/MakeWithValue.hpp"
  #include "Utilities/TMPL.hpp"
  
  /// \cond
  template <typename>
  class Variables;
  
  namespace Tags {
  template <typename Tag>
  struct Normalized;
  }  // namespace Tags
  /// \endcond
  
  namespace ScalarWave {
  /// @{
  /*!
   * \brief Compute the characteristic speeds for the scalar wave system.
   *
   * Computes the speeds as described in "Optimal constraint projection for
   * hyperbolic evolution systems" by Holst et al. \cite Holst2004wt
   * [see text following Eq.(32)]. The characteristic fields' names used here
   * differ from this paper:
   *
   * \f{align*}
   * \mathrm{SpECTRE} && \mathrm{Holst} \\
   * v^{\hat \psi} && Z^1 \\
   * v^{\hat 0}_{i} && Z^{2}_{i} \\
   * v^{\hat \pm} && u^{1\pm}
   * \f}
   *
   * The corresponding characteristic speeds \f$\lambda_{\hat \alpha}\f$ are given
   * in the text following Eq.(38) of \cite Holst2004wt :
   *
   * \f{align*}
   * \lambda_{\hat \psi} =& 0 \\
   * \lambda_{\hat 0} =& 0 \\
   * \lambda_{\hat \pm} =& \pm 1.
   * \f}
   */
  template <size_t Dim>
  std::array<DataVector, 4> characteristic_speeds(
      const tnsr::i<DataVector, Dim, Frame::Inertial>&
          unit_normal_one_form) noexcept;
  
  template <size_t Dim>
  void characteristic_speeds(
      gsl::not_null<std::array<DataVector, 4>*> char_speeds,
      const tnsr::i<DataVector, Dim, Frame::Inertial>&
          unit_normal_one_form) noexcept;
  
  namespace Tags {
  template <size_t Dim>
  struct CharacteristicSpeedsCompute : Tags::CharacteristicSpeeds<Dim>,
                                       db::ComputeTag {
    using base = Tags::CharacteristicSpeeds<Dim>;
    using return_type = typename base::type;
    using argument_tags =
        tmpl::list<::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>>;
  
    static void function(gsl::not_null<return_type*> char_speeds,
                         const tnsr::i<DataVector, Dim, Frame::Inertial>&
                             unit_normal_one_form) noexcept {
      characteristic_speeds(char_speeds, unit_normal_one_form);
    }
  };
  }  // namespace Tags
  /// @}
  
  /// @{
  /*!
   * \brief Computes characteristic fields from evolved fields
   *
   * \ref Tags::CharacteristicFieldsCompute and
   * \ref Tags::EvolvedFieldsFromCharacteristicFieldsCompute convert between
   * characteristic and evolved fields for the scalar-wave system.
   *
   * \ref Tags::CharacteristicFieldsCompute computes
   * characteristic fields as described in "Optimal constraint projection for
   * hyperbolic evolution systems" by Holst et al. \cite Holst2004wt .
   * Their names used here differ from this paper:
   *
   * \f{align*}
   * \mathrm{SpECTRE} && \mathrm{Holst} \\
   * v^{\hat \psi} && Z^1 \\
   * v^{\hat 0}_{i} && Z^{2}_{i} \\
   * v^{\hat \pm} && u^{1\pm}
   * \f}
   *
   * The characteristic fields \f${v}^{\hat \alpha}\f$ are given in terms of
   * the evolved fields by Eq.(33) - (35) of \cite Holst2004wt, respectively:
   *
   * \f{align*}
   * v^{\hat \psi} =& \psi \\
   * v^{\hat 0}_{i} =& (\delta^k_i - n_i n^k) \Phi_{k} := P^k_i \Phi_{k} \\
   * v^{\hat \pm} =& \Pi \pm n^i \Phi_{i} - \gamma_2\psi
   * \f}
   *
   * where \f$\psi\f$ is the scalar field, \f$\Phi_{i}=\partial_i \psi\f$ is an
   * auxiliary variable, \f$\Pi\f$ is a conjugate momentum, \f$\gamma_2\f$
   * is a constraint damping parameter, and \f$n_k\f$ is the unit normal to the
   * surface along which the characteristic fields are defined.
   *
   * \ref Tags::EvolvedFieldsFromCharacteristicFieldsCompute computes evolved
   * fields \f$u_\alpha\f$ in terms of the characteristic fields. This uses the
   * inverse of above relations:
   *
   * \f{align*}
   * \psi =& v^{\hat \psi}, \\
   * \Pi =& \frac{1}{2}(v^{\hat +} + v^{\hat -}) + \gamma_2 v^{\hat \psi}, \\
   * \Phi_{i} =& \frac{1}{2}(v^{\hat +} - v^{\hat -}) n_i + v^{\hat 0}_{i}.
   * \f}
   *
   * The corresponding characteristic speeds \f$\lambda_{\hat \alpha}\f$
   * are computed by \ref Tags::CharacteristicSpeedsCompute .
   */
  template <size_t Dim>
  Variables<tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus>>
  characteristic_fields(const Scalar<DataVector>& gamma_2,
                        const Scalar<DataVector>& psi,
                        const Scalar<DataVector>& pi,
                        const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
                        const tnsr::i<DataVector, Dim, Frame::Inertial>&
                            unit_normal_one_form) noexcept;
  
  template <size_t Dim>
  void characteristic_fields(
      gsl::not_null<Variables<
          tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus>>*>
          char_fields,
      const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& psi,
      const Scalar<DataVector>& pi,
      const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
      const tnsr::i<DataVector, Dim, Frame::Inertial>&
          unit_normal_one_form) noexcept;
  
  namespace Tags {
  template <size_t Dim>
  struct CharacteristicFieldsCompute : Tags::CharacteristicFields<Dim>,
                                       db::ComputeTag {
    using base = Tags::CharacteristicFields<Dim>;
    using return_type = typename base::type;
    using argument_tags =
        tmpl::list<Tags::ConstraintGamma2, Psi, Pi, Phi<Dim>,
                   ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>>;
  
    static void function(const gsl::not_null<return_type*> char_fields,
                         const Scalar<DataVector>& gamma_2,
                         const Scalar<DataVector>& psi,
                         const Scalar<DataVector>& pi,
                         const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
                         const tnsr::i<DataVector, Dim, Frame::Inertial>&
                             unit_normal_one_form) noexcept {
      characteristic_fields(char_fields, gamma_2, psi, pi, phi,
                            unit_normal_one_form);
    };
  };
  }  // namespace Tags
  /// @}
  
  /// @{
  /*!
   * \brief Compute evolved fields from characteristic fields.
   *
   * For expressions used here to compute evolved fields from characteristic ones,
   * see \ref Tags::CharacteristicFieldsCompute.
   */
  template <size_t Dim>
  Variables<tmpl::list<Psi, Pi, Phi<Dim>>>
  evolved_fields_from_characteristic_fields(
      const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& v_psi,
      const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero,
      const Scalar<DataVector>& v_plus, const Scalar<DataVector>& v_minus,
      const tnsr::i<DataVector, Dim, Frame::Inertial>&
          unit_normal_one_form) noexcept;
  
  template <size_t Dim>
  void evolved_fields_from_characteristic_fields(
      gsl::not_null<Variables<tmpl::list<Psi, Pi, Phi<Dim>>>*> evolved_fields,
      const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& v_psi,
      const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero,
      const Scalar<DataVector>& v_plus, const Scalar<DataVector>& v_minus,
      const tnsr::i<DataVector, Dim, Frame::Inertial>&
          unit_normal_one_form) noexcept;
  
  namespace Tags {
  template <size_t Dim>
  struct EvolvedFieldsFromCharacteristicFieldsCompute
      : Tags::EvolvedFieldsFromCharacteristicFields<Dim>,
        db::ComputeTag {
    using base = Tags::EvolvedFieldsFromCharacteristicFields<Dim>;
    using return_type = typename base::type;
    using argument_tags =
        tmpl::list<Tags::ConstraintGamma2, Tags::VPsi, Tags::VZero<Dim>,
                   Tags::VPlus, Tags::VMinus,
                   ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>>;
  
    static void function(const gsl::not_null<return_type*> evolved_fields,
                         const Scalar<DataVector>& gamma_2,
                         const Scalar<DataVector>& v_psi,
                         const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero,
                         const Scalar<DataVector>& v_plus,
                         const Scalar<DataVector>& v_minus,
                         const tnsr::i<DataVector, Dim, Frame::Inertial>&
                             unit_normal_one_form) noexcept {
      evolved_fields_from_characteristic_fields(evolved_fields, gamma_2, v_psi,
                                                v_zero, v_plus, v_minus,
                                                unit_normal_one_form);
    };
  };
  
  struct LargestCharacteristicSpeed : db::SimpleTag {
    using type = double;
  };
  
  /// Compute the maximum magnitude of the characteristic speeds.
  struct ComputeLargestCharacteristicSpeed : LargestCharacteristicSpeed,
                                             db::ComputeTag {
    using argument_tags = tmpl::list<>;
    using return_type = double;
    using base = LargestCharacteristicSpeed;
    SPECTRE_ALWAYS_INLINE static constexpr void function(
        const gsl::not_null<double*> speed) noexcept {
      *speed = 1.0;
    }
  };
  }  // namespace Tags
  /// @}
  }  // namespace ScalarWave
  