Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines functions and tags for taking a divergence.
6 :
7 : #pragma once
8 :
9 : #include <cstddef>
10 : #include <string>
11 :
12 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 : #include "DataStructures/DataBox/Tag.hpp"
14 : #include "DataStructures/DataBox/TagName.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "Utilities/TMPL.hpp"
17 : #include "Utilities/TypeTraits.hpp"
18 : #include "Utilities/TypeTraits/IsA.hpp"
19 :
20 : /// \cond
21 : class DataVector;
22 : template <size_t Dim>
23 : class Mesh;
24 : template <typename TagsList>
25 : class Variables;
26 :
27 : namespace domain {
28 : namespace Tags {
29 : template <size_t Dim>
30 : struct Mesh;
31 : } // namespace Tags
32 : } // namespace domain
33 : /// \endcond
34 :
35 : namespace Tags {
36 : /// \ingroup DataBoxTagsGroup
37 : /// \brief Prefix indicating the divergence
38 : ///
39 : /// Prefix indicating the divergence of a Tensor.
40 : ///
41 : /// \see Tags::DivVectorCompute Tags::DivVariablesCompute
42 : template <typename Tag>
43 : requires(tt::is_a_v<Tensor, typename Tag::type>)
44 1 : struct div : db::PrefixTag, db::SimpleTag {
45 0 : using tag = Tag;
46 0 : using type = TensorMetafunctions::remove_first_index<typename Tag::type>;
47 : };
48 : } // namespace Tags
49 :
50 : /// @{
51 : /// \ingroup NumericalAlgorithmsGroup
52 : /// \brief Compute the (Euclidean) divergence of fluxes
53 : template <typename FluxTags, size_t Dim, typename DerivativeFrame>
54 1 : auto divergence(const Variables<FluxTags>& F, const Mesh<Dim>& mesh,
55 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
56 : DerivativeFrame>& inverse_jacobian)
57 : -> Variables<db::wrap_tags_in<Tags::div, FluxTags>>;
58 :
59 : template <typename... DivTags, typename... FluxTags, size_t Dim,
60 : typename DerivativeFrame>
61 1 : void divergence(
62 : gsl::not_null<Variables<tmpl::list<DivTags...>>*> divergence_of_F,
63 : const Variables<tmpl::list<FluxTags...>>& F, const Mesh<Dim>& mesh,
64 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
65 : DerivativeFrame>& inverse_jacobian);
66 : /// @}
67 :
68 : /// @{
69 : /// \ingroup NumericalAlgorithmsGroup
70 : /// \brief Compute the divergence of the vector `input`
71 : template <typename DataType, size_t Dim, typename DerivativeFrame>
72 1 : Scalar<DataType> divergence(
73 : const tnsr::I<DataType, Dim, DerivativeFrame>& input, const Mesh<Dim>& mesh,
74 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
75 : DerivativeFrame>& inverse_jacobian);
76 :
77 : template <typename DataType, size_t Dim, typename DerivativeFrame>
78 1 : void divergence(gsl::not_null<Scalar<DataType>*> div_input,
79 : const tnsr::I<DataType, Dim, DerivativeFrame>& input,
80 : const Mesh<Dim>& mesh,
81 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
82 : DerivativeFrame>& inverse_jacobian);
83 : /// @}
84 :
85 : /// @{
86 : /*!
87 : * \ingroup NumericalAlgorithmsGroup
88 : * \brief Compute the divergence of fluxes where a Cartoon basis is being
89 : * utilized.
90 : *
91 : * The additional parameter `inertial_coords` is used for division by the
92 : * \f$x\f$ coordinates. If \f$x=0\f$ is included in the domain, it is assumed to
93 : * be present only at the first index and is handled by L'Hôpital's rule.
94 : *
95 : * The mesh is required to have the Cartoon basis in the last and potentially
96 : * second-to-last coordinates and the inverse jacobian is accordingly used only
97 : * in the first and potentially second dimensions.
98 : *
99 : * \see cartoon_partial_derivatives for details on Cartoon derivatives.
100 : */
101 : template <typename... DivTags, typename... FluxTags, size_t Dim,
102 : typename DerivativeFrame, Requires<Dim == 3> = nullptr>
103 1 : void cartoon_divergence(
104 : gsl::not_null<Variables<tmpl::list<DivTags...>>*> divergence_of_F,
105 : const Variables<tmpl::list<FluxTags...>>& F, const Mesh<Dim>& mesh,
106 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
107 : DerivativeFrame>& inverse_jacobian_3d,
108 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
109 : /// @}
110 :
111 : /// @{
112 : /*!
113 : * \ingroup NumericalAlgorithmsGroup
114 : * \brief Calls the correct divergence function, either normal divergence or
115 : * cartoon divergence, as determined by mesh basis.
116 : *
117 : * If you have a `Variables` with several tensors with Cartoon bases you need
118 : * to find the divergence of, you should use the `divergence` function
119 : * that operates on `Variables` since that'll be more efficient.
120 : */
121 : template <typename... DivTags, typename... FluxTags, size_t Dim,
122 : typename DerivativeFrame>
123 1 : void divergence(
124 : gsl::not_null<Variables<tmpl::list<DivTags...>>*> div_fluxes,
125 : const Variables<tmpl::list<FluxTags...>>& fluxes, const Mesh<Dim>& mesh,
126 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
127 : DerivativeFrame>& inverse_jacobian_3d,
128 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
129 : template <typename DataType, size_t Dim, typename DerivativeFrame>
130 1 : void divergence(
131 : gsl::not_null<Scalar<DataType>*> div_input,
132 : const tnsr::I<DataType, Dim, DerivativeFrame>& input, const Mesh<Dim>& mesh,
133 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
134 : DerivativeFrame>& inverse_jacobian,
135 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
136 : template <typename DataType, size_t Dim, typename DerivativeFrame>
137 1 : Scalar<DataType> divergence(
138 : const tnsr::I<DataType, Dim, DerivativeFrame>& input, const Mesh<Dim>& mesh,
139 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
140 : DerivativeFrame>& inverse_jacobian,
141 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
142 : /// @}
143 :
144 : /// @{
145 : /*!
146 : * \brief Compute the divergence of fluxes in logical coordinates
147 : *
148 : * Applies the logical differentiation matrix to the fluxes in each dimension
149 : * and sums over dimensions.
150 : *
151 : * \see divergence
152 : */
153 : template <typename ResultTensor, typename FluxTensor, size_t Dim>
154 1 : void logical_divergence(gsl::not_null<ResultTensor*> div_flux,
155 : const FluxTensor& flux, const Mesh<Dim>& mesh);
156 :
157 : template <typename FluxTags, size_t Dim>
158 1 : auto logical_divergence(const Variables<FluxTags>& flux, const Mesh<Dim>& mesh)
159 : -> Variables<db::wrap_tags_in<Tags::div, FluxTags>>;
160 :
161 : template <typename... DivTags, typename... FluxTags, size_t Dim>
162 1 : void logical_divergence(
163 : gsl::not_null<Variables<tmpl::list<DivTags...>>*> div_flux,
164 : const Variables<tmpl::list<FluxTags...>>& flux, const Mesh<Dim>& mesh);
165 : /// @}
166 :
167 : namespace Tags {
168 : /*!
169 : * \ingroup DataBoxTagsGroup
170 : * \brief Compute the divergence of a Variables
171 : *
172 : * Computes the divergence of the every Tensor in the Variables represented by
173 : * `Tag`. The first index of each Tensor must be an upper spatial index, i.e.,
174 : * the first index must have type
175 : * `TensorIndexType<Dim, UpLo::Up, Frame::TargetFrame, IndexType::Spatial>`.
176 : * The divergence is computed in the frame `TargetFrame`, and
177 : * `InverseJacobianTag` must be associated with a map from
178 : * `Frame::ElementLogical` to `Frame::TargetFrame`.
179 : *
180 : * Note that each tensor may have additional tensor indices - in this case the
181 : * divergence is computed for each additional index. For instance, a tensor
182 : * \f$F^i_{ab}\f$ has divergence
183 : * \f$Div_{ab} = \partial_i F^i_{ab}\f$. This is to accommodate evolution
184 : * equations where the evolved variables \f$u_\alpha\f$ are higher-rank tensors
185 : * and thus their fluxes can be written as \f$F^i_\alpha\f$. A simple example
186 : * would be the fluid velocity in hydro systems, where we would write the flux
187 : * as \f$F^{ij}\f$.
188 : *
189 : * This tag inherits from `db::add_tag_prefix<Tags::div, Tag>`.
190 : */
191 : template <typename Tag, typename MeshTag, typename InverseJacobianTag>
192 1 : struct DivVariablesCompute : db::add_tag_prefix<div, Tag>, db::ComputeTag {
193 : private:
194 0 : using inv_jac_indices = typename InverseJacobianTag::type::index_list;
195 0 : static constexpr auto dim = tmpl::back<inv_jac_indices>::dim;
196 : static_assert(std::is_same_v<typename tmpl::front<inv_jac_indices>::Frame,
197 : Frame::ElementLogical>,
198 : "Must map from the logical frame.");
199 :
200 : public:
201 0 : using base = db::add_tag_prefix<div, Tag>;
202 0 : using return_type = typename base::type;
203 0 : static constexpr void (*function)(
204 : const gsl::not_null<return_type*>, const typename Tag::type&,
205 : const Mesh<dim>&, const typename InverseJacobianTag::type&) = &divergence;
206 0 : using argument_tags =
207 : tmpl::list<Tag, domain::Tags::Mesh<dim>, InverseJacobianTag>;
208 : };
209 :
210 : /// \ingroup DataBoxTagsGroup
211 : /// \brief Compute the divergence of a `tnsr::I` (vector)
212 : ///
213 : /// This tag inherits from `db::add_tag_prefix<Tags::div, Tag>`.
214 : ///
215 : /// For an executable that does not allow a Cartoon basis, the last parameter,
216 : /// `InertialCoordsTag`, should not be passed.
217 : template <typename Tag, typename MeshTag, typename InverseJacobianTag,
218 : typename InertialCoordsTag = void>
219 1 : struct DivVectorCompute : div<Tag>, db::ComputeTag {
220 : private:
221 0 : using inv_jac_indices = typename InverseJacobianTag::type::index_list;
222 0 : static constexpr auto dim = tmpl::back<inv_jac_indices>::dim;
223 : static_assert(std::is_same_v<typename tmpl::front<inv_jac_indices>::Frame,
224 : Frame::ElementLogical>,
225 : "Must map from the logical frame.");
226 :
227 : public:
228 0 : using base = div<Tag>;
229 0 : using return_type = typename base::type;
230 0 : static constexpr void function(
231 : const gsl::not_null<return_type*> div_input,
232 : const typename Tag::type& input, const Mesh<dim>& mesh,
233 : const typename InverseJacobianTag::type& inverse_jacobian) {
234 : divergence(div_input, input, mesh, inverse_jacobian);
235 : }
236 0 : static constexpr void function(
237 : const gsl::not_null<return_type*> div_input,
238 : const typename Tag::type& input, const Mesh<dim>& mesh,
239 : const typename InverseJacobianTag::type& inverse_jacobian,
240 : const tnsr::I<DataVector, dim, Frame::Inertial>& inertial_coords) {
241 : divergence(div_input, input, mesh, inverse_jacobian, inertial_coords);
242 : }
243 0 : using argument_tags = tmpl::conditional_t<
244 : std::is_same_v<void, InertialCoordsTag>,
245 : tmpl::list<Tag, MeshTag, InverseJacobianTag>,
246 : tmpl::list<Tag, MeshTag, InverseJacobianTag, InertialCoordsTag>>;
247 : };
248 : } // namespace Tags
|
