Merge pull request #5654 from nilsvu/transform_tensor

Move frame transformations to Tensor lib, add function to transform flux tensors to logical frame

src/DataStructures/Tensor/EagerMath/CMakeLists.txt

@@ -4,6 +4,7 @@
 spectre_target_sources(
   ${LIBRARY}
   PRIVATE
+  FrameTransform.cpp
   OrthonormalOneform.cpp
   RaiseOrLowerIndex.cpp
   Trace.cpp
@@ -17,6 +18,7 @@ spectre_target_headers(
   Determinant.hpp
   DeterminantAndInverse.hpp
   DotProduct.hpp
+  FrameTransform.hpp
   Magnitude.hpp
   Norms.hpp
   OrthonormalOneform.hpp

src/PointwiseFunctions/GeneralRelativity/Transform.cpp → src/DataStructures/Tensor/EagerMath/FrameTransform.cpp

@@ -1,7 +1,7 @@
 // Distributed under the MIT License.
 // See LICENSE.txt for details.
 
-#include "PointwiseFunctions/GeneralRelativity/Transform.hpp"
+#include "DataStructures/Tensor/EagerMath/FrameTransform.hpp"
 
 #include <limits>
 
@@ -275,38 +275,48 @@ auto to_different_frame(
   return dest;
 }
 
-template <size_t VolumeDim, typename SrcFrame, typename DestFrame>
+template <typename ResultTensor, typename InputTensor, typename DataType,
+          size_t Dim, typename SourceFrame, typename TargetFrame>
 void first_index_to_different_frame(
-    const gsl::not_null<tnsr::ijj<DataVector, VolumeDim, DestFrame>*> dest,
-    const Tensor<DataVector, tmpl::integral_list<std::int32_t, 2, 1, 1>,
-                 index_list<SpatialIndex<VolumeDim, UpLo::Lo, SrcFrame>,
-                            SpatialIndex<VolumeDim, UpLo::Lo, DestFrame>,
-                            SpatialIndex<VolumeDim, UpLo::Lo, DestFrame>>>& src,
-    const Jacobian<DataVector, VolumeDim, DestFrame, SrcFrame>& jacobian) {
-  for (size_t i = 0; i < VolumeDim; ++i) {
-    for (size_t j = i; j < VolumeDim; ++j) {  // symmetry
-      for (size_t k = 0; k < VolumeDim; ++k) {
-        dest->get(k, i, j) = jacobian.get(0, k) * src.get(0, i, j);
-        for (size_t p = 1; p < VolumeDim; ++p) {
-          dest->get(k, i, j) += jacobian.get(p, k) * src.get(p, i, j);
-        }
-      }
+    const gsl::not_null<ResultTensor*> result, const InputTensor& input,
+    const InverseJacobian<DataType, Dim, SourceFrame, TargetFrame>&
+        inv_jacobian) {
+  static_assert(InputTensor::rank() > 0,
+                "Cannot transform scalar to different frame");
+  static_assert(
+      std::is_same_v<TensorMetafunctions::remove_first_index<ResultTensor>,
+                     TensorMetafunctions::remove_first_index<InputTensor>>,
+      "The input and result tensors must be the same except for the first "
+      "index");
+  using first_index = tmpl::front<typename InputTensor::index_list>;
+  static_assert(
+      std::is_same_v<first_index, SpatialIndex<Dim, UpLo::Up, TargetFrame>>,
+      "This function is currently only tested for transforming an upper "
+      "spatial index but can be generalized.");
+  for (size_t storage_index = 0; storage_index < ResultTensor::size();
+       ++storage_index) {
+    const auto result_index = ResultTensor::get_tensor_index(storage_index);
+    auto input_index = result_index;
+    input_index[0] = 0;
+    result->get(result_index) =
+        input.get(input_index) * inv_jacobian.get(result_index[0], 0);
+    for (size_t d = 1; d < Dim; ++d) {
+      input_index[0] = d;
+      result->get(result_index) +=
+          input.get(input_index) * inv_jacobian.get(result_index[0], d);
     }
   }
 }
 
-template <size_t VolumeDim, typename SrcFrame, typename DestFrame>
-auto first_index_to_different_frame(
-    const Tensor<DataVector, tmpl::integral_list<std::int32_t, 2, 1, 1>,
-                 index_list<SpatialIndex<VolumeDim, UpLo::Lo, SrcFrame>,
-                            SpatialIndex<VolumeDim, UpLo::Lo, DestFrame>,
-                            SpatialIndex<VolumeDim, UpLo::Lo, DestFrame>>>& src,
-    const Jacobian<DataVector, VolumeDim, DestFrame, SrcFrame>& jacobian)
-    -> tnsr::ijj<DataVector, VolumeDim, DestFrame> {
-  auto dest = make_with_value<tnsr::ijj<DataVector, VolumeDim, DestFrame>>(
-      src, std::numeric_limits<double>::signaling_NaN());
-  first_index_to_different_frame(make_not_null(&dest), src, jacobian);
-  return dest;
+template <typename InputTensor, typename DataType, size_t Dim,
+          typename SourceFrame, typename TargetFrame, typename ResultTensor>
+ResultTensor first_index_to_different_frame(
+    const InputTensor& input,
+    const InverseJacobian<DataType, Dim, SourceFrame, TargetFrame>&
+        inv_jacobian) {
+  ResultTensor result{};
+  first_index_to_different_frame(make_not_null(&result), input, inv_jacobian);
+  return result;
 }
 
 }  // namespace transform
@@ -385,36 +395,32 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (Frame::Inertial),
 
 #undef INSTANTIATE
 
-#undef TENSOR
-#undef DTYPE
-
-#define INSTANTIATE(_, data)                                                  \
-  template void transform::first_index_to_different_frame(                    \
-      const gsl::not_null<tnsr::ijj<DataVector, DIM(data), DESTFRAME(data)>*> \
-          dest,                                                               \
-      const Tensor<                                                           \
-          DataVector, tmpl::integral_list<std::int32_t, 2, 1, 1>,             \
-          index_list<SpatialIndex<DIM(data), UpLo::Lo, SRCFRAME(data)>,       \
-                     SpatialIndex<DIM(data), UpLo::Lo, DESTFRAME(data)>,      \
-                     SpatialIndex<DIM(data), UpLo::Lo, DESTFRAME(data)>>>&    \
-          src,                                                                \
-      const Jacobian<DataVector, DIM(data), DESTFRAME(data), SRCFRAME(data)>& \
-          jacobian);                                                          \
-  template auto transform::first_index_to_different_frame(                    \
-      const Tensor<                                                           \
-          DataVector, tmpl::integral_list<std::int32_t, 2, 1, 1>,             \
-          index_list<SpatialIndex<DIM(data), UpLo::Lo, SRCFRAME(data)>,       \
-                     SpatialIndex<DIM(data), UpLo::Lo, DESTFRAME(data)>,      \
-                     SpatialIndex<DIM(data), UpLo::Lo, DESTFRAME(data)>>>&    \
-          src,                                                                \
-      const Jacobian<DataVector, DIM(data), DESTFRAME(data), SRCFRAME(data)>& \
-          jacobian)                                                           \
-      ->tnsr::ijj<DataVector, DIM(data), DESTFRAME(data)>;
-GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3),
-                        (Frame::BlockLogical, Frame::ElementLogical),
-                        (Frame::Grid, Frame::Distorted))
+#define INSTANTIATE(_, data)                                                   \
+  template void transform::first_index_to_different_frame(                     \
+      const gsl::not_null<TensorMetafunctions::prepend_spatial_index<          \
+                              TensorMetafunctions::remove_first_index<         \
+                                  tnsr::TENSOR(data) < DTYPE(data), DIM(data), \
+                                  DESTFRAME(data)>>,                           \
+                          DIM(data), UpLo::Up, SRCFRAME(data)>* > result,      \
+      const tnsr::TENSOR(data) < DTYPE(data), DIM(data),                       \
+      DESTFRAME(data) > &input,                                                \
+      const InverseJacobian<DTYPE(data), DIM(data), SRCFRAME(data),            \
+                            DESTFRAME(data)>& inv_jacobian);                   \
+  template auto transform::first_index_to_different_frame(                     \
+      const tnsr::TENSOR(data) < DTYPE(data), DIM(data),                       \
+      DESTFRAME(data) > &input,                                                \
+      const InverseJacobian<DTYPE(data), DIM(data), SRCFRAME(data),            \
+                            DESTFRAME(data)>& inv_jacobian)                    \
+      ->TensorMetafunctions::prepend_spatial_index<                            \
+          TensorMetafunctions::remove_first_index<                             \
+              tnsr::TENSOR(data) < DTYPE(data), DIM(data), DESTFRAME(data)>>,  \
+      DIM(data), UpLo::Up, SRCFRAME(data) > ;
+GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (Frame::ElementLogical),
+                        (Frame::Inertial), (double, DataVector), (I, II, Ij))
 
 #undef DIM
 #undef SRCFRAME
 #undef DESTFRAME
+#undef TENSOR
+#undef DTYPE
 #undef INSTANTIATE

src/PointwiseFunctions/GeneralRelativity/Transform.hpp → src/DataStructures/Tensor/EagerMath/FrameTransform.hpp

@@ -6,6 +6,7 @@
 #include <cstddef>
 
 #include "DataStructures/DataBox/Tag.hpp"
+#include "DataStructures/Tensor/Metafunctions.hpp"
 #include "DataStructures/Tensor/TypeAliases.hpp"
 #include "Utilities/Gsl.hpp"
 
@@ -183,45 +184,40 @@ auto to_different_frame(
  * \ingroup GeneralRelativityGroup
  * Transforms only the first index to different frame.
  *
- * Often used for derivatives: When representing derivatives as
- * tensors, the first index is typically the derivative index.
- * Numerical derivatives must be computed in the logical frame or
- * sometimes the grid frame (independent of the frame of the tensor
- * being differentiated), and then that derivative index must later
- * be transformed into the same frame as the other indices of the
- * tensor.
+ * ## Examples for transforming some tensors
+ *
+ * ### Flux vector to logical coordinates
+ *
+ * A common example is taking the divergence of a flux vector $F^i$, where the
+ * index $i$ is in inertial coordinates $x^i$ and the divergence is taken in
+ * logical coordinates $\xi^\hat{i}$. So to transform the flux vector to logical
+ * coordinates before taking the divergence we do this:
  *
- * The formula for transforming \f$T_{i\bar{\jmath}\bar{k}}\f$ is
  * \f{align}
- *   T_{\bar{\imath}\bar{\jmath}\bar{k}} &= T_{i\bar{\jmath}\bar{k}}
- *      \frac{\partial x^i}{\partial x^{\bar{\imath}}},
+ *   F^\hat{i} &= F^i \frac{\partial x^\hat{i}}{\partial x^i}
  * \f}
- * where \f$x^i\f$ are the source coordinates and
- * \f$x^{\bar{\imath}}\f$ are the destination coordinates.
  *
- * Note that `Jacobian<DestFrame,SrcFrame>` is the same type as
- * `InverseJacobian<SrcFrame,DestFrame>` and represents
- * \f$\partial x^i/\partial x^{\bar{\jmath}}\f$.
+ * Here, $\frac{\partial x^\hat{i}}{\partial x^i}$ is the inverse Jacobian (see
+ * definitions in TypeAliases.hpp).
  *
- * In principle `first_index_to_different_frame` can be
- * extended/generalized to other tensor types if needed.
+ * Currently, this function is tested for any tensor with a first upper spatial
+ * index. It can be extended/generalized to other tensor types if needed.
  */
-template <size_t VolumeDim, typename SrcFrame, typename DestFrame>
+template <typename ResultTensor, typename InputTensor, typename DataType,
+          size_t Dim, typename SourceFrame, typename TargetFrame>
 void first_index_to_different_frame(
-    const gsl::not_null<tnsr::ijj<DataVector, VolumeDim, DestFrame>*> dest,
-    const Tensor<DataVector, tmpl::integral_list<std::int32_t, 2, 1, 1>,
-                 index_list<SpatialIndex<VolumeDim, UpLo::Lo, SrcFrame>,
-                            SpatialIndex<VolumeDim, UpLo::Lo, DestFrame>,
-                            SpatialIndex<VolumeDim, UpLo::Lo, DestFrame>>>& src,
-    const Jacobian<DataVector, VolumeDim, DestFrame, SrcFrame>& jacobian);
-
-template <size_t VolumeDim, typename SrcFrame, typename DestFrame>
-auto first_index_to_different_frame(
-    const Tensor<DataVector, tmpl::integral_list<std::int32_t, 2, 1, 1>,
-                 index_list<SpatialIndex<VolumeDim, UpLo::Lo, SrcFrame>,
-                            SpatialIndex<VolumeDim, UpLo::Lo, DestFrame>,
-                            SpatialIndex<VolumeDim, UpLo::Lo, DestFrame>>>& src,
-    const Jacobian<DataVector, VolumeDim, DestFrame, SrcFrame>& jacobian)
-    -> tnsr::ijj<DataVector, VolumeDim, DestFrame>;
+    gsl::not_null<ResultTensor*> result, const InputTensor& input,
+    const InverseJacobian<DataType, Dim, SourceFrame, TargetFrame>&
+        inv_jacobian);
+
+template <typename InputTensor, typename DataType, size_t Dim,
+          typename SourceFrame, typename TargetFrame,
+          typename ResultTensor = TensorMetafunctions::prepend_spatial_index<
+              TensorMetafunctions::remove_first_index<InputTensor>, Dim,
+              UpLo::Up, SourceFrame>>
+ResultTensor first_index_to_different_frame(
+    const InputTensor& input,
+    const InverseJacobian<DataType, Dim, SourceFrame, TargetFrame>&
+        inv_jacobian);
 /// @}
 }  // namespace transform

src/DataStructures/Variables/FrameTransform.hpp

@@ -0,0 +1,59 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/Tensor/EagerMath/FrameTransform.hpp"
+#include "DataStructures/Tensor/Metafunctions.hpp"
+#include "DataStructures/Tensor/TypeAliases.hpp"
+#include "DataStructures/Variables.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace transform {
+
+/// Tags to represent the result of frame-transforming Variables
+namespace Tags {
+/// The `Tag` with the first index transformed to a different frame
+template <typename Tag, typename FirstIndexFrame>
+struct TransformedFirstIndex : db::SimpleTag {
+  using type = TensorMetafunctions::prepend_spatial_index<
+      TensorMetafunctions::remove_first_index<typename Tag::type>,
+      tmpl::front<typename Tag::type::index_list>::dim, UpLo::Up,
+      FirstIndexFrame>;
+};
+}  // namespace Tags
+
+/// @{
+/// Transforms the first index of all tensors in the Variables to a different
+/// frame
+///
+/// See single-Tensor overload for details.
+template <typename... ResultTags, typename... InputTags, size_t Dim,
+          typename SourceFrame, typename TargetFrame>
+void first_index_to_different_frame(
+    const gsl::not_null<Variables<tmpl::list<ResultTags...>>*> result,
+    const Variables<tmpl::list<InputTags...>>& input,
+    const InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>&
+        inv_jacobian) {
+  EXPAND_PACK_LEFT_TO_RIGHT(
+      first_index_to_different_frame(make_not_null(&get<ResultTags>(*result)),
+                                     get<InputTags>(input), inv_jacobian));
+}
+
+template <typename... InputTags, size_t Dim, typename SourceFrame,
+          typename TargetFrame,
+          typename ResultVars = Variables<tmpl::list<
+              Tags::TransformedFirstIndex<InputTags, SourceFrame>...>>>
+ResultVars first_index_to_different_frame(
+    const Variables<tmpl::list<InputTags...>>& input,
+    const InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>&
+        inv_jacobian) {
+  ResultVars result{input.number_of_grid_points()};
+  first_index_to_different_frame(make_not_null(&result), input, inv_jacobian);
+  return result;
+}
+/// @}
+
+}  // namespace transform

src/ParallelAlgorithms/ApparentHorizonFinder/ComputeExcisionBoundaryVolumeQuantities.tpp

@@ -13,6 +13,7 @@
 #include "DataStructures/Tags/TempTensor.hpp"
 #include "DataStructures/TempBuffer.hpp"
 #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
+#include "DataStructures/Tensor/EagerMath/FrameTransform.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Tensor/TypeAliases.hpp"
 #include "DataStructures/Variables.hpp"
@@ -28,7 +29,6 @@
 #include "PointwiseFunctions/GeneralRelativity/SpacetimeNormalVector.hpp"
 #include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
-#include "PointwiseFunctions/GeneralRelativity/Transform.hpp"
 #include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/TMPL.hpp"

src/ParallelAlgorithms/ApparentHorizonFinder/ComputeHorizonVolumeQuantities.tpp

@@ -13,6 +13,7 @@
 #include "DataStructures/Tags/TempTensor.hpp"
 #include "DataStructures/TempBuffer.hpp"
 #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
+#include "DataStructures/Tensor/EagerMath/FrameTransform.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Tensor/TypeAliases.hpp"
 #include "DataStructures/Variables.hpp"
@@ -28,7 +29,6 @@
 #include "PointwiseFunctions/GeneralRelativity/SpacetimeNormalVector.hpp"
 #include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
-#include "PointwiseFunctions/GeneralRelativity/Transform.hpp"
 #include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/TMPL.hpp"
@@ -199,11 +199,6 @@ void ComputeHorizonVolumeQuantities::apply(
   using inertial_metric_tag = gr::Tags::SpatialMetric<DataVector, 3>;
   using inertial_inv_metric_tag = gr::Tags::InverseSpatialMetric<DataVector, 3>;
   using inertial_ex_curvature_tag = gr::Tags::ExtrinsicCurvature<DataVector, 3>;
-  using logical_deriv_metric_tag = ::Tags::TempTensor<
-      6, Tensor<DataVector, tmpl::integral_list<std::int32_t, 2, 1, 1>,
-                index_list<SpatialIndex<3, UpLo::Lo, Frame::ElementLogical>,
-                           SpatialIndex<3, UpLo::Lo, TargetFrame>,
-                           SpatialIndex<3, UpLo::Lo, TargetFrame>>>>;
   using deriv_metric_tag = ::Tags::Tempijj<7, 3, TargetFrame, DataVector>;
   using inertial_spatial_ricci_tag = gr::Tags::SpatialRicci<DataVector, 3>;
 
@@ -213,7 +208,7 @@ void ComputeHorizonVolumeQuantities::apply(
       tmpl::list<metric_tag, inv_metric_tag, lapse_tag, shift_tag,
                  spacetime_normal_vector_tag, inertial_metric_tag,
                  inertial_inv_metric_tag, inertial_ex_curvature_tag,
-                 logical_deriv_metric_tag, deriv_metric_tag>;
+                 deriv_metric_tag>;
   using full_temp_tags_list = std::conditional_t<
       tmpl::list_contains_v<DestTagList,
                             gr::Tags::SpatialRicci<DataVector, 3, TargetFrame>>,
@@ -236,7 +231,6 @@ void ComputeHorizonVolumeQuantities::apply(
   auto& lapse = get<lapse_tag>(buffer);
   auto& shift = get<shift_tag>(buffer);
   auto& spacetime_normal_vector = get<spacetime_normal_vector_tag>(buffer);
-  auto& logical_deriv_metric = get<logical_deriv_metric_tag>(buffer);
   auto& deriv_metric = get<deriv_metric_tag>(buffer);
 
   // These may or may not be temporaries, depending on if they are asked for
@@ -278,11 +272,8 @@ void ComputeHorizonVolumeQuantities::apply(
                           metric);
 
   // Differentiate 3-metric.
-  logical_partial_derivative(make_not_null(&logical_deriv_metric), metric,
-                             mesh);
-  transform::first_index_to_different_frame(make_not_null(&deriv_metric),
-                                            logical_deriv_metric,
-                                            invjac_logical_to_target);
+  partial_derivative(make_not_null(&deriv_metric), metric, mesh,
+                     invjac_logical_to_target);
   gr::christoffel_second_kind(make_not_null(&spatial_christoffel_second_kind),
                               deriv_metric, inv_metric);
 

src/PointwiseFunctions/AnalyticData/GrMhd/PolarMagnetizedFmDisk.hpp

@@ -11,14 +11,14 @@
 #include <type_traits>
 
 #include "DataStructures/Tensor/EagerMath/Determinant.hpp"
+#include "DataStructures/Tensor/EagerMath/FrameTransform.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "Options/Options.hpp"
 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
 #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
 #include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp"
 #include "PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp"
-#include "PointwiseFunctions/GeneralRelativity/Transform.hpp"
 #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/Serialization/CharmPupable.hpp"