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 <string>
       8             : #include <utility>
       9             : 
      10             : #include "DataStructures/DataBox/Tag.hpp"
      11             : #include "DataStructures/DataBox/TagName.hpp"
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "DataStructures/Tensor/Tensor.hpp"
      14             : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
      15             : #include "NumericalAlgorithms/SphericalHarmonics/StrahlkorperFunctions.hpp"
      16             : #include "NumericalAlgorithms/SphericalHarmonics/TagsDeclarations.hpp"
      17             : #include "NumericalAlgorithms/SphericalHarmonics/TagsTypeAliases.hpp"
      18             : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
      19             : #include "Utilities/ForceInline.hpp"
      20             : #include "Utilities/Gsl.hpp"
      21             : #include "Utilities/TMPL.hpp"
      22             : 
      23             : /// \ingroup SurfacesGroup
      24             : /// Holds tags and ComputeItems associated with a `ylm::Strahlkorper`.
      25           1 : namespace ylm::Tags {
      26             : 
      27             : /// Tag referring to a `ylm::Strahlkorper`
      28             : template <typename Frame>
      29           1 : struct Strahlkorper : db::SimpleTag {
      30           0 :   using type = ylm::Strahlkorper<Frame>;
      31             : };
      32             : 
      33             : /// @{
      34             : /// \f$(\theta,\phi)\f$ on the grid.
      35             : /// Doesn't depend on the shape of the surface.
      36             : template <typename Frame>
      37           1 : struct ThetaPhi : db::SimpleTag {
      38           0 :   using type = tnsr::i<DataVector, 2, ::Frame::Spherical<Frame>>;
      39             : };
      40             : 
      41             : template <typename Frame>
      42           0 : struct ThetaPhiCompute : ThetaPhi<Frame>, db::ComputeTag {
      43           0 :   using base = ThetaPhi<Frame>;
      44           0 :   using return_type = tnsr::i<DataVector, 2, ::Frame::Spherical<Frame>>;
      45           0 :   static constexpr auto function = static_cast<void (*)(
      46             :       const gsl::not_null<tnsr::i<DataVector, 2, ::Frame::Spherical<Frame>>*>,
      47             :       const ylm::Strahlkorper<Frame>&)>(&::ylm::theta_phi<Frame>);
      48           0 :   using argument_tags = tmpl::list<Strahlkorper<Frame>>;
      49             : };
      50             : /// @}
      51             : 
      52             : /// @{
      53             : /// `Rhat(i)` is \f$\hat{r}^i = x_i/\sqrt{x^2+y^2+z^2}\f$ on the grid.
      54             : /// Doesn't depend on the shape of the surface.
      55             : template <typename Frame>
      56           1 : struct Rhat : db::SimpleTag {
      57           0 :   using type = tnsr::i<DataVector, 3, Frame>;
      58             : };
      59             : 
      60             : template <typename Frame>
      61           0 : struct RhatCompute : Rhat<Frame>, db::ComputeTag {
      62           0 :   using base = Rhat<Frame>;
      63           0 :   using return_type = tnsr::i<DataVector, 3, Frame>;
      64           0 :   static constexpr auto function = static_cast<void (*)(
      65             :       const gsl::not_null<tnsr::i<DataVector, 3, Frame>*>,
      66             :       const tnsr::i<DataVector, 2, ::Frame::Spherical<Frame>>&)>(
      67             :       &::ylm::rhat<Frame>);
      68           0 :   using argument_tags = tmpl::list<ThetaPhi<Frame>>;
      69             : };
      70             : /// @}
      71             : 
      72             : /// @{
      73             : /// `Jacobian(i,0)` is \f$\frac{1}{r}\partial x^i/\partial\theta\f$,
      74             : /// and `Jacobian(i,1)`
      75             : /// is \f$\frac{1}{r\sin\theta}\partial x^i/\partial\phi\f$.
      76             : /// Here \f$r\f$ means \f$\sqrt{x^2+y^2+z^2}\f$.
      77             : /// `Jacobian` doesn't depend on the shape of the surface.
      78             : template <typename Frame>
      79           1 : struct Jacobian : db::SimpleTag {
      80           0 :   using type = aliases::Jacobian<Frame>;
      81             : };
      82             : 
      83             : template <typename Frame>
      84           0 : struct JacobianCompute : Jacobian<Frame>, db::ComputeTag {
      85           0 :   using base = Jacobian<Frame>;
      86           0 :   using return_type = aliases::Jacobian<Frame>;
      87           0 :   static constexpr auto function = static_cast<void (*)(
      88             :       const gsl::not_null<ylm::Tags::aliases::Jacobian<Frame>*>,
      89             :       const tnsr::i<DataVector, 2, ::Frame::Spherical<Frame>>&)>(
      90             :       &::ylm::jacobian<Frame>);
      91           0 :   using argument_tags = tmpl::list<ThetaPhi<Frame>>;
      92             : };
      93             : /// @}
      94             : 
      95             : /// @{
      96             : /// `InvJacobian(0,i)` is \f$r\partial\theta/\partial x^i\f$,
      97             : /// and `InvJacobian(1,i)` is \f$r\sin\theta\partial\phi/\partial x^i\f$.
      98             : /// Here \f$r\f$ means \f$\sqrt{x^2+y^2+z^2}\f$.
      99             : /// `InvJacobian` doesn't depend on the shape of the surface.
     100             : template <typename Frame>
     101           1 : struct InvJacobian : db::SimpleTag {
     102           0 :   using type = aliases::InvJacobian<Frame>;
     103             : };
     104             : 
     105             : template <typename Frame>
     106           0 : struct InvJacobianCompute : InvJacobian<Frame>, db::ComputeTag {
     107           0 :   using base = InvJacobian<Frame>;
     108           0 :   using return_type = aliases::InvJacobian<Frame>;
     109           0 :   static constexpr auto function = static_cast<void (*)(
     110             :       const gsl::not_null<ylm::Tags::aliases::InvJacobian<Frame>*>,
     111             :       const tnsr::i<DataVector, 2, ::Frame::Spherical<Frame>>&)>(
     112             :       &::ylm::inv_jacobian<Frame>);
     113           0 :   using argument_tags = tmpl::list<ThetaPhi<Frame>>;
     114             : };
     115             : /// @}
     116             : 
     117             : /// @{
     118             : /// `InvHessian(k,i,j)` is \f$\partial (J^{-1}){}^k_j/\partial x^i\f$,
     119             : /// where \f$(J^{-1}){}^k_j\f$ is the inverse Jacobian.
     120             : /// `InvHessian` is not symmetric because the Jacobians are Pfaffian.
     121             : /// `InvHessian` doesn't depend on the shape of the surface.
     122             : template <typename Frame>
     123           1 : struct InvHessian : db::SimpleTag {
     124           0 :   using type = aliases::InvHessian<Frame>;
     125             : };
     126             : 
     127             : template <typename Frame>
     128           0 : struct InvHessianCompute : InvHessian<Frame>, db::ComputeTag {
     129           0 :   using base = InvHessian<Frame>;
     130           0 :   using return_type = aliases::InvHessian<Frame>;
     131           0 :   static constexpr auto function = static_cast<void (*)(
     132             :       const gsl::not_null<ylm::Tags::aliases::InvHessian<Frame>*>,
     133             :       const tnsr::i<DataVector, 2, ::Frame::Spherical<Frame>>&)>(
     134             :       &::ylm::inv_hessian<Frame>);
     135           0 :   using argument_tags = tmpl::list<ThetaPhi<Frame>>;
     136             : };
     137             : /// @}
     138             : 
     139             : /// @{
     140             : /// (Euclidean) distance \f$r_{\rm surf}(\theta,\phi)\f$ from the center to each
     141             : /// point of the surface.
     142             : template <typename Frame>
     143           1 : struct Radius : db::SimpleTag {
     144           0 :   using type = Scalar<DataVector>;
     145             : };
     146             : 
     147             : template <typename Frame>
     148           0 : struct RadiusCompute : Radius<Frame>, db::ComputeTag {
     149           0 :   using base = Radius<Frame>;
     150           0 :   using return_type = Scalar<DataVector>;
     151           0 :   static constexpr auto function =
     152             :       static_cast<void (*)(const gsl::not_null<Scalar<DataVector>*>,
     153             :                            const ylm::Strahlkorper<Frame>&)>(
     154             :           &(::ylm::radius<Frame>));
     155           0 :   using argument_tags = tmpl::list<Strahlkorper<Frame>>;
     156             : };
     157             : /// @}
     158             : 
     159             : /// @{
     160             : /// The geometrical center of the surface.  Uses
     161             : /// `ylm::Strahlkorper::physical_center`.
     162             : template <typename Frame>
     163           1 : struct PhysicalCenter : db::SimpleTag {
     164           0 :   using type = std::array<double, 3>;
     165             : };
     166             : 
     167             : template <typename Frame>
     168           0 : struct PhysicalCenterCompute : PhysicalCenter<Frame>, db::ComputeTag {
     169           0 :   using base = PhysicalCenter<Frame>;
     170           0 :   using return_type = std::array<double, 3>;
     171           0 :   static void function(gsl::not_null<std::array<double, 3>*> physical_center,
     172             :                        const ylm::Strahlkorper<Frame>& strahlkorper);
     173           0 :   using argument_tags = tmpl::list<Strahlkorper<Frame>>;
     174             : };
     175             : /// @}
     176             : 
     177             : /// @{
     178             : /// `CartesianCoords(i)` is \f$x_{\rm surf}^i\f$,
     179             : /// the vector of \f$(x,y,z)\f$ coordinates of each point
     180             : /// on the surface.
     181             : template <typename Frame>
     182           1 : struct CartesianCoords : db::SimpleTag {
     183           0 :   using type = tnsr::I<DataVector, 3, Frame>;
     184             : };
     185             : 
     186             : template <typename Frame>
     187           0 : struct CartesianCoordsCompute : CartesianCoords<Frame>, db::ComputeTag {
     188           0 :   using base = CartesianCoords<Frame>;
     189           0 :   using return_type = tnsr::I<DataVector, 3, Frame>;
     190           0 :   static constexpr auto function = static_cast<void (*)(
     191             :       const gsl::not_null<tnsr::I<DataVector, 3, Frame>*> coords,
     192             :       const ylm::Strahlkorper<Frame>& strahlkorper,
     193             :       const Scalar<DataVector>& radius,
     194             :       const tnsr::i<DataVector, 3, Frame>& r_hat)>(&ylm::cartesian_coords);
     195           0 :   using argument_tags =
     196             :       tmpl::list<Strahlkorper<Frame>, Radius<Frame>, Rhat<Frame>>;
     197             : };
     198             : /// @}
     199             : 
     200             : /// @{
     201             : /// `DxRadius(i)` is \f$\partial r_{\rm surf}/\partial x^i\f$.  Here
     202             : /// \f$r_{\rm surf}=r_{\rm surf}(\theta,\phi)\f$ is the function
     203             : /// describing the surface, which is considered a function of
     204             : /// Cartesian coordinates
     205             : /// \f$r_{\rm surf}=r_{\rm surf}(\theta(x,y,z),\phi(x,y,z))\f$
     206             : /// for this operation.
     207             : template <typename Frame>
     208           1 : struct DxRadius : db::SimpleTag {
     209           0 :   using type = tnsr::i<DataVector, 3, Frame>;
     210             : };
     211             : 
     212             : template <typename Frame>
     213           0 : struct DxRadiusCompute : DxRadius<Frame>, db::ComputeTag {
     214           0 :   using base = DxRadius<Frame>;
     215           0 :   using return_type = tnsr::i<DataVector, 3, Frame>;
     216           0 :   static constexpr auto function = static_cast<void (*)(
     217             :       const gsl::not_null<tnsr::i<DataVector, 3, Frame>*> dx_radius,
     218             :       const Scalar<DataVector>& scalar,
     219             :       const ylm::Strahlkorper<Frame>& strahlkorper,
     220             :       const Scalar<DataVector>& radius_of_strahlkorper,
     221             :       const aliases::InvJacobian<Frame>& inv_jac)>(
     222             :       &ylm::cartesian_derivs_of_scalar);
     223           0 :   using argument_tags = tmpl::list<Radius<Frame>, Strahlkorper<Frame>,
     224             :                                    Radius<Frame>, InvJacobian<Frame>>;
     225             : };
     226             : /// @}
     227             : 
     228             : /// @{
     229             : /// `D2xRadius(i,j)` is
     230             : /// \f$\partial^2 r_{\rm surf}/\partial x^i\partial x^j\f$. Here
     231             : /// \f$r_{\rm surf}=r_{\rm surf}(\theta,\phi)\f$ is the function
     232             : /// describing the surface, which is considered a function of
     233             : /// Cartesian coordinates
     234             : /// \f$r_{\rm surf}=r_{\rm surf}(\theta(x,y,z),\phi(x,y,z))\f$
     235             : /// for this operation.
     236             : template <typename Frame>
     237           1 : struct D2xRadius : db::SimpleTag {
     238           0 :   using type = tnsr::ii<DataVector, 3, Frame>;
     239             : };
     240             : 
     241             : template <typename Frame>
     242           0 : struct D2xRadiusCompute : D2xRadius<Frame>, db::ComputeTag {
     243           0 :   using base = D2xRadius<Frame>;
     244           0 :   using return_type = tnsr::ii<DataVector, 3, Frame>;
     245           0 :   static constexpr auto function = static_cast<void (*)(
     246             :       gsl::not_null<tnsr::ii<DataVector, 3, Frame>*> d2x_radius,
     247             :       const Scalar<DataVector>& scalar,
     248             :       const ylm::Strahlkorper<Frame>& strahlkorper,
     249             :       const Scalar<DataVector>& radius_of_strahlkorper,
     250             :       const aliases::InvJacobian<Frame>& inv_jac,
     251             :       const aliases::InvHessian<Frame>& inv_hess)>(
     252             :       &ylm::cartesian_second_derivs_of_scalar);
     253           0 :   using argument_tags =
     254             :       tmpl::list<Radius<Frame>, Strahlkorper<Frame>, Radius<Frame>,
     255             :                  InvJacobian<Frame>, InvHessian<Frame>>;
     256             : };
     257             : /// @}
     258             : 
     259             : /// @{
     260             : /// \f$\nabla^2 r_{\rm surf}\f$, the flat Laplacian of the surface.
     261             : /// This is \f$\eta^{ij}\partial^2 r_{\rm surf}/\partial x^i\partial x^j\f$,
     262             : /// where \f$r_{\rm surf}=r_{\rm surf}(\theta(x,y,z),\phi(x,y,z))\f$.
     263             : template <typename Frame>
     264           1 : struct LaplacianRadius : db::SimpleTag {
     265           0 :   using type = Scalar<DataVector>;
     266             : };
     267             : 
     268             : template <typename Frame>
     269           0 : struct LaplacianRadiusCompute : LaplacianRadius<Frame>, db::ComputeTag {
     270           0 :   using base = LaplacianRadius<Frame>;
     271           0 :   using return_type = Scalar<DataVector>;
     272           0 :   static constexpr auto function = static_cast<void (*)(
     273             :       gsl::not_null<Scalar<DataVector>*> lap_radius,
     274             :       const Scalar<DataVector>& radius,
     275             :       const ylm::Strahlkorper<Frame>& strahlkorper,
     276             :       const tnsr::i<DataVector, 2, ::Frame::Spherical<Frame>>& theta_phi)>(
     277             :       &ylm::laplacian_of_scalar);
     278           0 :   using argument_tags =
     279             :       tmpl::list<Radius<Frame>, Strahlkorper<Frame>, ThetaPhi<Frame>>;
     280             : };
     281             : /// @}
     282             : 
     283             : /// @{
     284             : /// `NormalOneForm(i)` is \f$s_i\f$, the (unnormalized) normal one-form
     285             : /// to the surface, expressed in Cartesian components.
     286             : /// This is computed by \f$x_i/r-\partial r_{\rm surf}/\partial x^i\f$,
     287             : /// where \f$x_i/r\f$ is `Rhat` and
     288             : /// \f$\partial r_{\rm surf}/\partial x^i\f$ is `DxRadius`.
     289             : /// See Eq. (8) of \cite Baumgarte1996hh.
     290             : /// Note on the word "normal": \f$s_i\f$ points in the correct direction
     291             : /// (it is "normal" to the surface), but it does not have unit length
     292             : /// (it is not "normalized"; normalization requires a metric).
     293             : template <typename Frame>
     294           1 : struct NormalOneForm : db::SimpleTag {
     295           0 :   using type = tnsr::i<DataVector, 3, Frame>;
     296             : };
     297             : 
     298             : template <typename Frame>
     299           0 : struct NormalOneFormCompute : NormalOneForm<Frame>, db::ComputeTag {
     300           0 :   using base = NormalOneForm<Frame>;
     301           0 :   using return_type = tnsr::i<DataVector, 3, Frame>;
     302           0 :   static constexpr auto function = static_cast<void (*)(
     303             :       gsl::not_null<tnsr::i<DataVector, 3, Frame>*> one_form,
     304             :       const tnsr::i<DataVector, 3, Frame>& dx_radius,
     305             :       const tnsr::i<DataVector, 3, Frame>& r_hat)>(&ylm::normal_one_form);
     306           0 :   using argument_tags = tmpl::list<DxRadius<Frame>, Rhat<Frame>>;
     307             : };
     308             : /// @}
     309             : 
     310             : /// @{
     311             : /// `Tangents(i,j)` is \f$\partial x_{\rm surf}^i/\partial q^j\f$,
     312             : /// where \f$x_{\rm surf}^i\f$ are the Cartesian coordinates of the
     313             : /// surface (i.e. `CartesianCoords`) and are considered functions of
     314             : /// \f$(\theta,\phi)\f$.
     315             : ///
     316             : /// \f$\partial/\partial q^0\f$ means
     317             : /// \f$\partial/\partial\theta\f$; and \f$\partial/\partial q^1\f$
     318             : /// means \f$\csc\theta\,\,\partial/\partial\phi\f$.  Note that the
     319             : /// vectors `Tangents(i,0)` and `Tangents(i,1)` are orthogonal to the
     320             : /// `NormalOneForm` \f$s_i\f$, i.e.
     321             : /// \f$s_i \partial x_{\rm surf}^i/\partial q^j = 0\f$; this statement
     322             : /// is independent of a metric.  Also, `Tangents(i,0)` and
     323             : /// `Tangents(i,1)` are not necessarily orthogonal to each other,
     324             : /// since orthogonality between 2 vectors (as opposed to a vector and
     325             : /// a one-form) is metric-dependent.
     326             : template <typename Frame>
     327           1 : struct Tangents : db::SimpleTag {
     328           0 :   using type = aliases::Jacobian<Frame>;
     329             : };
     330             : 
     331             : template <typename Frame>
     332           0 : struct TangentsCompute : Tangents<Frame>, db::ComputeTag {
     333           0 :   using base = Tangents<Frame>;
     334           0 :   using return_type = aliases::Jacobian<Frame>;
     335           0 :   static constexpr auto function =
     336             :       static_cast<void (*)(gsl::not_null<aliases::Jacobian<Frame>*> tangents,
     337             :                            const ylm::Strahlkorper<Frame>& strahlkorper,
     338             :                            const Scalar<DataVector>& radius,
     339             :                            const tnsr::i<DataVector, 3, Frame>& r_hat,
     340             :                            const aliases::Jacobian<Frame>& jac)>(
     341             :           &ylm::tangents);
     342           0 :   using argument_tags = tmpl::list<Strahlkorper<Frame>, Radius<Frame>,
     343             :                                    Rhat<Frame>, Jacobian<Frame>>;
     344             : };
     345             : /// @}
     346             : 
     347             : /// Tag for holding the previously-found values of a Strahlkorper,
     348             : /// which are saved for extrapolation for future initial guesses
     349             : /// and for computing the time deriv of a Strahlkorper.
     350             : template <typename Frame>
     351           1 : struct PreviousStrahlkorpers : db::SimpleTag {
     352           0 :   using type = std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>;
     353             : };
     354             : 
     355             : /// @{
     356             : /// Tag to compute the time derivative of the coefficients of a Strahlkorper
     357             : /// from a number of previous Strahlkorpers.
     358             : template <typename Frame>
     359           1 : struct TimeDerivStrahlkorper : db::SimpleTag {
     360           0 :   using type = ylm::Strahlkorper<Frame>;
     361             : };
     362             : 
     363             : template <typename Frame>
     364           0 : struct TimeDerivStrahlkorperCompute : db::ComputeTag,
     365             :                                       TimeDerivStrahlkorper<Frame> {
     366           0 :   using base = TimeDerivStrahlkorper<Frame>;
     367           0 :   using return_type = typename base::type;
     368           0 :   static constexpr auto function = static_cast<void (*)(
     369             :       gsl::not_null<ylm::Strahlkorper<Frame>*>,
     370             :       const std::deque<std::pair<double, ylm::Strahlkorper<Frame>>>&)>(
     371             :       &ylm::time_deriv_of_strahlkorper<Frame>);
     372             : 
     373           0 :   using argument_tags = tmpl::list<PreviousStrahlkorpers<Frame>>;
     374             : };
     375             : /// @}
     376             : 
     377             : template <typename Frame>
     378           0 : using items_tags = tmpl::list<Strahlkorper<Frame>>;
     379             : 
     380             : template <typename Frame>
     381           0 : using compute_items_tags =
     382             :     tmpl::list<ThetaPhiCompute<Frame>, RhatCompute<Frame>,
     383             :                JacobianCompute<Frame>, InvJacobianCompute<Frame>,
     384             :                InvHessianCompute<Frame>, RadiusCompute<Frame>,
     385             :                CartesianCoordsCompute<Frame>, DxRadiusCompute<Frame>,
     386             :                D2xRadiusCompute<Frame>, LaplacianRadiusCompute<Frame>,
     387             :                NormalOneFormCompute<Frame>, TangentsCompute<Frame>>;
     388             : }  // namespace ylm::Tags
     389             : 
     390             : /// Tags related to symmetric trace-free tensors
     391           0 : namespace Stf::Tags {
     392             : 
     393             : /*!
     394             :  * \brief Tag used to hold a symmetric trace-free tensor of a certain rank.
     395             :  * \details The type is a symmetric tensor of the requested rank. A
     396             :  * ScalarBaseTag of type `Scalar<double>` is used to identify the tag of which
     397             :  * the symmetric trace-free expansion is done.
     398             :  */
     399             : template <typename ScalarBaseTag, size_t rank, size_t Dim, typename Frame>
     400           1 : struct StfTensor : db::SimpleTag {
     401             :   static_assert(std::is_same_v<typename ScalarBaseTag::type, Scalar<double>>,
     402             :                 "StfTensor base tags must be a Scalar.");
     403             :   static_assert(rank <= 3, "StfTensor tag is only implemented up to rank 3.");
     404           0 :   static std::string name() {
     405             :     return MakeString{} << "StfTensor(" << db::tag_name<ScalarBaseTag>()
     406             :                         << "," << rank << ")";
     407             :   }
     408           0 :   using type_list =
     409             :       tmpl::list<Scalar<double>, tnsr::i<double, Dim, Frame>,
     410             :                  tnsr::ii<double, Dim, Frame>, tnsr::iii<double, Dim, Frame>>;
     411           0 :   using type = tmpl::at<type_list, tmpl::size_t<rank>>;
     412             : };
     413             : }  // namespace Stf::Tags