Norms.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <string>
7 
9 #include "DataStructures/DataVector.hpp"
11 #include "Utilities/Functional.hpp"
12 #include "Utilities/Numeric.hpp"
13 #include "Utilities/TMPL.hpp"
14 
15 namespace L2Norm_detail {
16 template <typename DataType, typename Symm, typename IndexList>
17 Scalar<DataType> pointwise_l2_norm_square(
18  const Tensor<DataType, Symm, IndexList>& tensor) noexcept {
19  auto pointwise_l2_normsq = make_with_value<Scalar<DataType>>(tensor, 0.);
20  for (auto tensor_element = tensor.begin(); tensor_element != tensor.end();
21  ++tensor_element) {
22  // In order to handle tensory symmetries, we multiply the square of each
23  // stored component with its multiplicity
24  get(pointwise_l2_normsq) +=
25  tensor.multiplicity(tensor_element) * square(*tensor_element);
26  }
27  return pointwise_l2_normsq;
28 }
29 } // namespace L2Norm_detail
30 
31 /*!
32  * \ingroup TensorGroup
33  * \brief Compute point-wise Euclidean \f$L^2\f$-norm of arbitrary Tensors.
34  *
35  * \details
36  * At each grid point \f$p\f$ in the element, this function computes the
37  * point-wise Frobenius norm of a given Tensor with arbitrary rank. If the
38  * Tensor \f$A\f$ has rank \f$n\f$ and dimensionality \f$D\f$, then its
39  * Frobenius norm at point \f$p\f$ is computed as:
40  *
41  * \f{equation}
42  * ||A||_2(p) =
43  * \left(\sum^{D-1}_{i_1=0}\sum^{D-1}_{i_2=0}\cdots \sum^{D-1}_{i_n=0}
44  * |A_{i_1 i_2 \cdots i_n}(p)|^2 \right)^{1/2},
45  * \f}
46  *
47  * where both contra-variant and co-variant indices are shown as lower indices.
48  */
49 template <typename DataType, typename Symm, typename IndexList>
51  const Tensor<DataType, Symm, IndexList>& tensor) noexcept {
52  return Scalar<DataType>{
53  sqrt(get(L2Norm_detail::pointwise_l2_norm_square(tensor)))};
54 }
55 
56 // @{
57 /*!
58  * \ingroup TensorGroup
59  * \brief Compute Euclidean \f$L^2\f$-norm of arbitrary Tensors reduced over an
60  * element.
61  *
62  * \details
63  * Computes the RMS value of the point-wise Frobenius norm of a given Tensor
64  * with arbitrary rank over all grid points in an element. If the Tensor \f$A\f$
65  * has rank \f$n\f$ and dimensionality \f$D\f$, and the element (of order
66  * \f$N\f$) has \f$N+1\f$ points, then its element-reduced Frobenius norm is
67  * computed as:
68  *
69  * \f{equation}
70  * ||A||_2 =
71  * \left(\frac{1}{N+1}\sum^{N}_{p=0}
72  * \left(\sum^{D-1}_{i_1=0}\sum^{D-1}_{i_2=0}\cdots
73  * \sum^{D-1}_{i_n=0} |A^p_{i_1 i_2 \cdots i_n}|^2 \right)
74  * \right)^{1/2},
75  * \f}
76  *
77  * where both contra-variant and co-variant indices are shown as lower indices,
78  * and \f$p\f$ indexes grid points in the element.
79  *
80  * \warning This function reduces the Frobenius norm over the element, not the
81  * whole domain.
82  */
83 template <typename DataType, typename Symm, typename IndexList>
84 double l2_norm(const Tensor<DataType, Symm, IndexList>& tensor) noexcept {
85  const auto pointwise_l2_normsq =
86  L2Norm_detail::pointwise_l2_norm_square(tensor);
87  using Plus = funcl::Plus<funcl::Identity>;
88  return sqrt(alg::accumulate(get(pointwise_l2_normsq), 0., Plus{}) /
89  tensor.begin()->size());
90 }
91 
92 namespace Tags {
93 /*!
94  * \ingroup DataBoxTagsGroup
95  * \ingroup DataStructuresGroup
96  * Point-wise Euclidean \f$L^2\f$-norm of a Tensor.
97  * \see `pointwise_l2_norm()` for details.
98  */
99 template <typename Tag>
102  static std::string name() noexcept {
103  return "PointwiseL2Norm(" + db::tag_name<Tag>() + ")";
104  }
105 };
106 
107 /*!
108  * \ingroup DataBoxTagsGroup
109  * \ingroup DataStructuresGroup
110  * Computes the point-wise Euclidean \f$L^2\f$-norm of a Tensor.
111  * \see `pointwise_l2_norm()` for details.
112  */
113 template <typename Tag>
115  using base = PointwiseL2Norm<Tag>;
116  static constexpr db::const_item_type<PointwiseL2Norm<Tag>> (*function)(
118  using argument_tags = tmpl::list<Tag>;
119 };
120 
121 /*!
122  * \ingroup DataBoxTagsGroup
123  * \ingroup DataStructuresGroup
124  * Euclidean \f$L^2\f$-norm of a Tensor, RMS over grid points in element.
125  * \see `l2_norm()` for details.
126  */
127 template <typename Tag>
129  using type = double;
130  static std::string name() noexcept {
131  return "L2Norm(" + db::tag_name<Tag>() + ")";
132  }
133 };
134 
135 /*!
136  * \ingroup DataBoxTagsGroup
137  * \ingroup DataStructuresGroup
138  * Computes the Euclidean \f$L^2\f$-norm of a Tensor, RMS over grid points in
139  * element. \see `l2_norm()` for details.
140  *
141  * \warning This compute tag reduces the Frobenius norm over the element, not
142  * the whole domain.
143  */
144 template <typename Tag>
146  using base = L2Norm<Tag>;
147  static constexpr db::const_item_type<L2Norm<Tag>> (*function)(
149  using argument_tags = tmpl::list<Tag>;
150 };
151 } // namespace Tags
Marks a DataBoxTag as being a compute item that executes a function.
Definition: DataBoxTag.hpp:154
Tags for the DataBox inherit from this type.
Definition: DataBoxTag.hpp:64
decltype(auto) accumulate(const Container &c, T init)
Convenience wrapper around std::accumulate, returns std::accumulate(begin(c), end(c), init).
Definition: Numeric.hpp:54
double l2_norm(const Tensor< DataType, Symm, IndexList > &tensor) noexcept
Compute Euclidean -norm of arbitrary Tensors reduced over an element.
Definition: Norms.hpp:84
Functional for computing + of two objects.
Definition: Functional.hpp:238
Definition: Norms.hpp:100
Definition: DataBoxTag.hpp:29
Defines classes for Tensor.
Definition: Norms.hpp:128
Wraps the template metaprogramming library used (brigand)
Scalar< DataType > pointwise_l2_norm(const Tensor< DataType, Symm, IndexList > &tensor) noexcept
Compute point-wise Euclidean -norm of arbitrary Tensors.
Definition: Norms.hpp:50
Definition: Norms.hpp:114
Definition: Norms.hpp:145
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:52
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Definition: Norms.hpp:15