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