Merge pull request #6730 from michaeldmurphy1/cartoon_partial_deriv

Implement DG Cartoon partial derivatives

Author: nilsdeppe

src/DataStructures/Tags/TempTensor.hpp

@@ -16,6 +16,9 @@ struct Inertial;
 }  // namespace Frame
 /// \endcond
 
+/*!
+ * \brief Contains objects related to Tags
+ */
 namespace Tags {
 template <size_t N, typename T>
 struct TempTensor : db::SimpleTag {
@@ -25,6 +28,31 @@ struct TempTensor : db::SimpleTag {
   }
 };
 
+/*!
+ * \brief A struct that constructs a `TempTensor` given a label \p N and type
+ * \p TensorType. Call with `make_temp_tensor<N>::apply<Type>`
+ *
+ * \see convert_to_temp_tensors
+ */
+template <size_t N>
+struct make_temp_tensor {
+  template <typename TensorType>
+  struct apply {
+    using type = ::Tags::TempTensor<N, TensorType>;
+  };
+};
+
+/*!
+ * \brief Takes a `tmpl::list` of `Tensor` types and a label \p N which is
+ * applied to all `TempTensor`s in the resulting list.
+ *
+ * \see make_temp_tensor
+ */
+template <typename List, size_t N>
+using convert_to_temp_tensors =
+    tmpl::transform<List,
+                    typename make_temp_tensor<N>::template apply<tmpl::_1>>;
+
 /// @{
 /// \ingroup PeoGroup
 /// Variables Tags for temporary tensors inside a function.

src/NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp

@@ -111,6 +111,79 @@ auto logical_partial_derivative(const Tensor<DataType, SymmList, IndexList>& u,
         Frame::ElementLogical>;
 /// @}
 
+/// @{
+/*!
+ * \ingroup NumericalAlgorithmsGroup
+ * \brief Compute the partial derivatives of each variable, using the Cartoon
+ * method for directions not in the computational domain.
+ *
+ * For a spacetime with some Killing vector \f$\xi^a\f$, not only does the
+ * metric respect the isometry via \f$\mathcal{L}_\xi g_{ab} = 0\f$, but the
+ * fields must as well.
+ *
+ * In the case of axial symmetry about the \f$y\f$-axis, the Killing vector is
+ * \f$\xi^a = (0, -z, 0, x)\f$, so we have that for a vector \f$V^a\f$,
+ *
+ * \f{align*}{
+ *     \mathcal{L}_\xi V^a &= \xi^b \partial_b V^a - V^b \partial_b \xi^a
+ *       = -z \partial_x V^a + x \partial_z V^a - V^b \partial_b \xi^a = 0.
+ * \f}
+ *
+ * Choosing the computational domain to be the \f$z=0\f$ plane, get the result:
+ *
+ * \f{align*}{
+ *     \partial_z V^a &= \frac{V^b \partial_b \xi^a}{x}.
+ * \f}
+ *
+ * In the case that $x=0$ is in the computational domain, L'H&ocirc;pital's
+ * rule is used. For a general tensor,
+ * \f$T^{a_0,\ldots,a_n}{}_{b_0,\ldots,b_m}\f$ we have
+ *
+ * \f{align*}{
+ *   \partial_z T^{a_0,\ldots,a_n}{}_{b_0,\ldots,b_m} &=
+ *   \left\{
+ *   \begin{array}{ll}
+ *     \partial_x(\sum_k
+ *     T^{a_0,\ldots,c_k,\ldots,a_n}{}_{b_0,\ldots,b_m}\partial_{c_k}
+ *     \xi^{a_k} & \\
+ *     \;\;\;\;\;-\sum_k
+ *     T^{a_0,\ldots,a_n}{}_{b_0,\ldots,c_k,\ldots,b_m}\partial_{b_k}
+ *     \xi^{c_k}
+ *     ) & \mathrm{if} \; x=0  \\
+ *     (\sum_k
+ *     T^{a_0,\ldots,c_k,\ldots,a_n}{}_{b_0,\ldots,b_m}\partial_{c_k}
+ *     \xi^{a_k} & \\
+ *     \;-\sum_k
+ *     T^{a_0,\ldots,a_n}{}_{b_0,\ldots,c_k,\ldots,b_m}\partial_{b_k}
+ *     \xi^{c_k})\frac{1}{x} & \mathrm{otherwise}
+ *   \end{array}\right.
+ * \f}
+ *
+ * In the case of spherical symmetry, another Killing vector
+ * \f$\xi'^a = (0, -y, x, 0) \f$ is used as well, with the above equation
+ * easily generalizing. In that case, the computational domain is only the
+ * \f$x\f$-axis.
+ *
+ * This function computes the partial derivatives of _all_ variables in
+ * `VariableTags`. The additional parameter `inertial_coords` is used for
+ * division by the \f$x\f$ coordinates. If \f$x=0\f$ is included in the domain,
+ * it is assumed to be present only at the first index and is handled by
+ * L'H&ocirc;pital's rule.
+ *
+ * The mesh is required to have the Cartoon basis in the last and potentially
+ * second-to-last coordinates and the inverse Jacobian is accordingly used only
+ * in the first and potentially second dimensions.
+ */
+template <typename ResultTags, typename VariableTags, size_t Dim,
+          typename DerivativeFrame, Requires<Dim == 3> = nullptr>
+void cartoon_partial_derivatives(
+    gsl::not_null<Variables<ResultTags>*> du, const Variables<VariableTags>& u,
+    const Mesh<Dim>& mesh,
+    const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                          DerivativeFrame>& inverse_jacobian_3d,
+    const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
+/// @}
+
 /// @{
 /// \ingroup NumericalAlgorithmsGroup
 /// \brief Compute the partial derivatives of each variable with respect to