Skip to content

Commit

Permalink
Merge pull request #5655 from nilsvu/index_manip
Browse files Browse the repository at this point in the history
Move `raise_or_lower_index` and `trace` to Tensor lib
  • Loading branch information
kidder authored Dec 9, 2023
2 parents ef0642c + de80a17 commit 7b3f85c
Show file tree
Hide file tree
Showing 85 changed files with 364 additions and 262 deletions.
4 changes: 4 additions & 0 deletions src/DataStructures/Tensor/EagerMath/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@ spectre_target_sources(
${LIBRARY}
PRIVATE
OrthonormalOneform.cpp
RaiseOrLowerIndex.cpp
Trace.cpp
)

spectre_target_headers(
Expand All @@ -19,6 +21,8 @@ spectre_target_headers(
Norms.hpp
OrthonormalOneform.hpp
OuterProduct.hpp
RaiseOrLowerIndex.hpp
Trace.hpp
)

add_subdirectory(Python)
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "PointwiseFunctions/GeneralRelativity/IndexManipulation.hpp"
#include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"

#include <cstddef>

Expand Down Expand Up @@ -55,49 +55,6 @@ void raise_or_lower_index(
}
}

template <typename DataType, typename Index0, typename Index1>
void trace_last_indices(
const gsl::not_null<Tensor<DataType, Symmetry<1>, index_list<Index0>>*>
trace_of_tensor,
const Tensor<DataType, Symmetry<2, 1, 1>,
index_list<Index0, Index1, Index1>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index1>,
change_index_up_lo<Index1>>>& metric) {
constexpr auto dimension_of_trace = Index0::dim;
constexpr auto dimension_of_metric = Index1::dim;
for (size_t i = 0; i < dimension_of_trace; ++i) {
trace_of_tensor->get(i) = tensor.get(i, 0, 0) * get<0, 0>(metric);
for (size_t j = 0; j < dimension_of_metric; ++j) {
if (j != 0) {
trace_of_tensor->get(i) += tensor.get(i, j, j) * metric.get(j, j);
}
for (size_t k = j + 1; k < dimension_of_metric; ++k) {
trace_of_tensor->get(i) += 2.0 * tensor.get(i, j, k) * metric.get(j, k);
}
}
}
}

template <typename DataType, typename Index0>
void trace(
const gsl::not_null<Scalar<DataType>*> trace,
const Tensor<DataType, Symmetry<1, 1>, index_list<Index0, Index0>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index0>,
change_index_up_lo<Index0>>>& metric) {
constexpr auto dimension = Index0::dim;
get(*trace) = get<0, 0>(tensor) * get<0, 0>(metric);
for (size_t i = 0; i < dimension; ++i) {
if (i != 0) {
get(*trace) += tensor.get(i, i) * metric.get(i, i);
}
for (size_t j = i + 1; j < dimension; ++j) { // symmetry
get(*trace) += 2.0 * tensor.get(i, j) * metric.get(i, j);
}
}
}

#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
#define DTYPE(data) BOOST_PP_TUPLE_ELEM(1, data)
#define FRAME(data) BOOST_PP_TUPLE_ELEM(2, data)
Expand All @@ -123,17 +80,7 @@ void trace(
tensor, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<change_index_up_lo<INDEX0(data)>, \
change_index_up_lo<INDEX0(data)>>>& metric); \
template void trace_last_indices( \
const gsl::not_null< \
Tensor<DTYPE(data), Symmetry<1>, index_list<INDEX0(data)>>*> \
result, \
const Tensor<DTYPE(data), Symmetry<2, 1, 1>, \
index_list<INDEX0(data), INDEX1(data), INDEX1(data)>>& \
tensor, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<change_index_up_lo<INDEX1(data)>, \
change_index_up_lo<INDEX1(data)>>>& metric);
change_index_up_lo<INDEX0(data)>>>& metric);

GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (double, DataVector),
(Frame::Grid, Frame::Distorted, Frame::Inertial,
Expand All @@ -149,23 +96,16 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (double, DataVector),
#undef INDEX1
#undef INSTANTIATE

#define INSTANTIATE2(_, data) \
template void raise_or_lower_index( \
const gsl::not_null< \
Tensor<DTYPE(data), Symmetry<1>, \
index_list<change_index_up_lo<INDEX0(data)>>>*> \
trace_of_tensor, \
const Tensor<DTYPE(data), Symmetry<1>, index_list<INDEX0(data)>>& \
tensor, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<change_index_up_lo<INDEX0(data)>, \
change_index_up_lo<INDEX0(data)>>>& metric); \
template void trace( \
const gsl::not_null<Scalar<DTYPE(data)>*> trace, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<INDEX0(data), INDEX0(data)>>& tensor, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<change_index_up_lo<INDEX0(data)>, \
#define INSTANTIATE2(_, data) \
template void raise_or_lower_index( \
const gsl::not_null< \
Tensor<DTYPE(data), Symmetry<1>, \
index_list<change_index_up_lo<INDEX0(data)>>>*> \
trace_of_tensor, \
const Tensor<DTYPE(data), Symmetry<1>, index_list<INDEX0(data)>>& \
tensor, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<change_index_up_lo<INDEX0(data)>, \
change_index_up_lo<INDEX0(data)>>>& metric);

GENERATE_INSTANTIATIONS(INSTANTIATE2, (1, 2, 3), (double, DataVector),
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -92,67 +92,3 @@ raise_or_lower_index(
return result;
}
/// @}

/// @{
/*!
* \ingroup GeneralRelativityGroup
* \brief Computes trace of a rank 3 tensor, which is symmetric in its last two
* indices, tracing the symmetric indices.
*
* \details For example, if \f$ T_{abc} \f$ is a tensor such that \f$T_{abc} =
* T_{acb} \f$ then \f$ T_a = g^{bc}T_{abc} \f$ is computed, where \f$ g^{bc}
* \f$ is the inverse metric. Note that indices \f$a,b,c,...\f$ can represent
* either spatial or spacetime indices, and can have either valence. You may
* have to add a new instantiation of this template if you need a new use case.
*/
template <typename DataType, typename Index0, typename Index1>
void trace_last_indices(
gsl::not_null<Tensor<DataType, Symmetry<1>, index_list<Index0>>*>
trace_of_tensor,
const Tensor<DataType, Symmetry<2, 1, 1>,
index_list<Index0, Index1, Index1>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index1>,
change_index_up_lo<Index1>>>& metric);

template <typename DataType, typename Index0, typename Index1>
Tensor<DataType, Symmetry<1>, index_list<Index0>> trace_last_indices(
const Tensor<DataType, Symmetry<2, 1, 1>,
index_list<Index0, Index1, Index1>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index1>,
change_index_up_lo<Index1>>>& metric) {
auto trace_of_tensor =
make_with_value<Tensor<DataType, Symmetry<1>, index_list<Index0>>>(metric,
0.);
trace_last_indices(make_not_null(&trace_of_tensor), tensor, metric);
return trace_of_tensor;
}
/// @}

/// @{
/*!
* \ingroup GeneralRelativityGroup
* \brief Computes trace of a rank-2 symmetric tensor.
* \details Computes \f$g^{ab}T_{ab}\f$ or \f$g_{ab}T^{ab}\f$ where \f$(a,b)\f$
* can be spatial or spacetime indices.
*/
template <typename DataType, typename Index0>
void trace(
gsl::not_null<Scalar<DataType>*> trace,
const Tensor<DataType, Symmetry<1, 1>, index_list<Index0, Index0>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index0>,
change_index_up_lo<Index0>>>& metric);

template <typename DataType, typename Index0>
Scalar<DataType> trace(
const Tensor<DataType, Symmetry<1, 1>, index_list<Index0, Index0>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index0>,
change_index_up_lo<Index0>>>& metric) {
Scalar<DataType> trace{};
::trace(make_not_null(&trace), tensor, metric);
return trace;
}
/// @}
113 changes: 113 additions & 0 deletions src/DataStructures/Tensor/EagerMath/Trace.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "DataStructures/Tensor/EagerMath/Trace.hpp"

#include <cstddef>

#include "DataStructures/Tensor/Tensor.hpp"
#include "Utilities/GenerateInstantiations.hpp"

template <typename DataType, typename Index0, typename Index1>
void trace_last_indices(
const gsl::not_null<Tensor<DataType, Symmetry<1>, index_list<Index0>>*>
trace_of_tensor,
const Tensor<DataType, Symmetry<2, 1, 1>,
index_list<Index0, Index1, Index1>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index1>,
change_index_up_lo<Index1>>>& metric) {
constexpr auto dimension_of_trace = Index0::dim;
constexpr auto dimension_of_metric = Index1::dim;
for (size_t i = 0; i < dimension_of_trace; ++i) {
trace_of_tensor->get(i) = tensor.get(i, 0, 0) * get<0, 0>(metric);
for (size_t j = 0; j < dimension_of_metric; ++j) {
if (j != 0) {
trace_of_tensor->get(i) += tensor.get(i, j, j) * metric.get(j, j);
}
for (size_t k = j + 1; k < dimension_of_metric; ++k) {
trace_of_tensor->get(i) += 2.0 * tensor.get(i, j, k) * metric.get(j, k);
}
}
}
}

template <typename DataType, typename Index0>
void trace(
const gsl::not_null<Scalar<DataType>*> trace,
const Tensor<DataType, Symmetry<1, 1>, index_list<Index0, Index0>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index0>,
change_index_up_lo<Index0>>>& metric) {
constexpr auto dimension = Index0::dim;
get(*trace) = get<0, 0>(tensor) * get<0, 0>(metric);
for (size_t i = 0; i < dimension; ++i) {
if (i != 0) {
get(*trace) += tensor.get(i, i) * metric.get(i, i);
}
for (size_t j = i + 1; j < dimension; ++j) { // symmetry
get(*trace) += 2.0 * tensor.get(i, j) * metric.get(i, j);
}
}
}

#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
#define DTYPE(data) BOOST_PP_TUPLE_ELEM(1, data)
#define FRAME(data) BOOST_PP_TUPLE_ELEM(2, data)
#define INDEXTYPE(data) BOOST_PP_TUPLE_ELEM(3, data)
#define UPORLO0(data) BOOST_PP_TUPLE_ELEM(4, data)
#define UPORLO1(data) BOOST_PP_TUPLE_ELEM(5, data)

// Note: currently only instantiate these functions for the cases
// where Index0 and Index1 have the same spatial dimension, Frame and
// IndexType. The functions above are generic and can handle any of
// them being different, if these cases are needed.
#define INDEX0(data) INDEXTYPE(data)<DIM(data), UPORLO0(data), FRAME(data)>
#define INDEX1(data) INDEXTYPE(data)<DIM(data), UPORLO1(data), FRAME(data)>

#define INSTANTIATE(_, data) \
template void trace_last_indices( \
const gsl::not_null< \
Tensor<DTYPE(data), Symmetry<1>, index_list<INDEX0(data)>>*> \
result, \
const Tensor<DTYPE(data), Symmetry<2, 1, 1>, \
index_list<INDEX0(data), INDEX1(data), INDEX1(data)>>& \
tensor, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<change_index_up_lo<INDEX1(data)>, \
change_index_up_lo<INDEX1(data)>>>& metric);

GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (double, DataVector),
(Frame::Grid, Frame::Distorted, Frame::Inertial,
Frame::Spherical<Frame::Inertial>,
Frame::Spherical<Frame::Distorted>,
Frame::Spherical<Frame::Grid>),
(SpatialIndex, SpacetimeIndex), (UpLo::Lo, UpLo::Up),
(UpLo::Lo, UpLo::Up))

#undef FRAME1
#undef UPORLO1
#undef INDEXTYPE1
#undef INDEX1
#undef INSTANTIATE

#define INSTANTIATE2(_, data) \
template void trace( \
const gsl::not_null<Scalar<DTYPE(data)>*> trace, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<INDEX0(data), INDEX0(data)>>& tensor, \
const Tensor<DTYPE(data), Symmetry<1, 1>, \
index_list<change_index_up_lo<INDEX0(data)>, \
change_index_up_lo<INDEX0(data)>>>& metric);

GENERATE_INSTANTIATIONS(INSTANTIATE2, (1, 2, 3), (double, DataVector),
(Frame::Grid, Frame::Distorted, Frame::Inertial),
(SpatialIndex, SpacetimeIndex), (UpLo::Lo, UpLo::Up))

#undef DIM
#undef DTYPE
#undef FRAME0
#undef UPORLO0
#undef INDEXTYPE0
#undef INDEX0
#undef INSTANTIATE2
72 changes: 72 additions & 0 deletions src/DataStructures/Tensor/EagerMath/Trace.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

#pragma once

#include "DataStructures/Tensor/Tensor.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/MakeWithValue.hpp"

/// @{
/*!
* \ingroup GeneralRelativityGroup
* \brief Computes trace of a rank 3 tensor, which is symmetric in its last two
* indices, tracing the symmetric indices.
*
* \details For example, if \f$ T_{abc} \f$ is a tensor such that \f$T_{abc} =
* T_{acb} \f$ then \f$ T_a = g^{bc}T_{abc} \f$ is computed, where \f$ g^{bc}
* \f$ is the inverse metric. Note that indices \f$a,b,c,...\f$ can represent
* either spatial or spacetime indices, and can have either valence. You may
* have to add a new instantiation of this template if you need a new use case.
*/
template <typename DataType, typename Index0, typename Index1>
void trace_last_indices(
gsl::not_null<Tensor<DataType, Symmetry<1>, index_list<Index0>>*>
trace_of_tensor,
const Tensor<DataType, Symmetry<2, 1, 1>,
index_list<Index0, Index1, Index1>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index1>,
change_index_up_lo<Index1>>>& metric);

template <typename DataType, typename Index0, typename Index1>
Tensor<DataType, Symmetry<1>, index_list<Index0>> trace_last_indices(
const Tensor<DataType, Symmetry<2, 1, 1>,
index_list<Index0, Index1, Index1>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index1>,
change_index_up_lo<Index1>>>& metric) {
auto trace_of_tensor =
make_with_value<Tensor<DataType, Symmetry<1>, index_list<Index0>>>(metric,
0.);
trace_last_indices(make_not_null(&trace_of_tensor), tensor, metric);
return trace_of_tensor;
}
/// @}

/// @{
/*!
* \ingroup GeneralRelativityGroup
* \brief Computes trace of a rank-2 symmetric tensor.
* \details Computes \f$g^{ab}T_{ab}\f$ or \f$g_{ab}T^{ab}\f$ where \f$(a,b)\f$
* can be spatial or spacetime indices.
*/
template <typename DataType, typename Index0>
void trace(
gsl::not_null<Scalar<DataType>*> trace,
const Tensor<DataType, Symmetry<1, 1>, index_list<Index0, Index0>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index0>,
change_index_up_lo<Index0>>>& metric);

template <typename DataType, typename Index0>
Scalar<DataType> trace(
const Tensor<DataType, Symmetry<1, 1>, index_list<Index0, Index0>>& tensor,
const Tensor<DataType, Symmetry<1, 1>,
index_list<change_index_up_lo<Index0>,
change_index_up_lo<Index0>>>& metric) {
Scalar<DataType> trace{};
::trace(make_not_null(&trace), tensor, metric);
return trace;
}
/// @}
2 changes: 0 additions & 2 deletions src/Elliptic/Systems/Elasticity/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,6 @@ target_link_libraries(
Utilities
INTERFACE
LinearOperators
PRIVATE
GeneralRelativity
)

add_subdirectory(BoundaryConditions)
Loading

0 comments on commit 7b3f85c

Please sign in to comment.