Skip to content

Commit 1b7d2f5

Browse files
at0m741at0m741
committed
AutoDiff proto working on BSSN
1 parent 1a5e7c2 commit 1b7d2f5

File tree

5 files changed

+652
-65
lines changed

5 files changed

+652
-65
lines changed

‎includes/Tensorium/Core/Vector.hpp‎

Lines changed: 42 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,8 @@ namespace tensorium {
8484
void resize(size_t n) { data.resize(n); }
8585
///@}
8686

87-
87+
88+
8889
/** @name Debug and Utilities */
8990
///@{
9091

@@ -100,20 +101,33 @@ namespace tensorium {
100101
/** @name Basic Operations */
101102
///@{
102103

104+
static Vector<K> canonical(int index, K dx, K dy, K dz) {
105+
Vector<K> out(3, K(0));
106+
if (index == 0) out(0) = dx;
107+
else if (index == 1) out(1) = dy;
108+
else if (index == 2) out(2) = dz;
109+
else throw std::invalid_argument("Index must be 0, 1, or 2");
110+
return out;
111+
}
103112
/**
104113
* @brief Subtract two vectors.
105114
* @param other The vector to subtract.
106115
* @return A new Vector containing the difference.
107116
*/
108117
__attribute__((always_inline, hot, flatten))
118+
119+
__attribute__((always_inline, hot, flatten))
109120
Vector<K> operator-(const Vector<K>& other) const {
110-
assert(data.size() == other.data.size());
111121
Vector<K> result(data.size());
112-
for (size_t i = 0; i < data.size(); ++i)
122+
size_t m = std::min(data.size(), other.data.size());
123+
for (size_t i = 0; i < m; ++i)
113124
result[i] = data[i] - other[i];
125+
for (size_t i = m; i < data.size(); ++i)
126+
result[i] = data[i];
114127
return result;
115128
}
116129

130+
117131
/**
118132
* @brief Add another vector to this one (in-place).
119133
* @param v The vector to add.
@@ -451,4 +465,28 @@ namespace tensorium {
451465
return r;
452466
}
453467
};
454-
}
468+
469+
470+
template<typename K>
471+
inline Vector<K> operator+(const Vector<K>& a, const Vector<K>& b) {
472+
size_t n = a.size();
473+
Vector<K> result(n);
474+
size_t m = std::min(n, b.size());
475+
for (size_t i = 0; i < m; ++i)
476+
result[i] = a[i] + b[i];
477+
for (size_t i = m; i < n; ++i)
478+
result[i] = a[i];
479+
return result;
480+
}
481+
482+
483+
template<typename K>
484+
inline Vector<K> operator-(const Vector<K>& a, const Vector<K>& b) {
485+
assert(a.size() == b.size());
486+
Vector<K> result(a.size());
487+
for (size_t i = 0; i < a.size(); ++i)
488+
result[i] = a[i] - b[i];
489+
return result;
490+
}
491+
492+
}

‎includes/Tensorium/DiffGeometry/BSSN/BSSNAutoDiff.hpp‎

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
#pragma once
2+
3+
#include "../Metric.hpp"
4+
#include "BSSNDerivatives.hpp"
5+
#include "../../Core/Tensor.hpp"
6+
7+
namespace tensorium {
8+
template<typename T>
9+
struct TensorTraits;
10+
11+
template<typename T, size_t Rank>
12+
struct TensorTraits<Tensor<T, Rank>> {
13+
static constexpr size_t rank = Rank;
14+
using value_type = T;
15+
};
16+
17+
template<typename T>
18+
struct TensorTraits {
19+
static constexpr size_t rank = 0;
20+
using value_type = T;
21+
};
22+
23+
template<typename T>
24+
struct TensorTraits<Vector<T>> {
25+
static constexpr size_t rank = 1;
26+
using value_type = T;
27+
};
28+
}
29+
30+
namespace tensorium_RG {
31+
32+
enum class DiffMode { PARTIAL, COV, COV2 };
33+
34+
// ===================== RANK 0 (SCALAR) =====================
35+
// Return type : Tensor<T,1> (i.e. un vecteur gradient)
36+
template<typename T, typename ScalarFunc>
37+
tensorium::Tensor<T, 1> autodiff_rank0(
38+
const tensorium::Vector<T>& X,
39+
T dx, T dy, T dz,
40+
ScalarFunc&& func,
41+
DiffMode mode,
42+
const tensorium::Tensor<T, 3>& christoffel)
43+
{
44+
using namespace tensorium;
45+
if (mode == DiffMode::PARTIAL) {
46+
return partial_scalar(X, dx, dy, dz, std::forward<ScalarFunc>(func));
47+
}
48+
else if (mode == DiffMode::COV2) {
49+
return covariant_scalar_second(X, dx, dy, dz, std::forward<ScalarFunc>(func), christoffel);
50+
}
51+
else {
52+
throw std::invalid_argument("autodiff_rank0: uniquement PARTIAL ou COV2 autorisés pour un champ scalaire");
53+
}
54+
}
55+
56+
// ===================== RANK 1 (VECTOR) =====================
57+
// Return type : Tensor<T,2> (i.e. une matrice dérivée covariante ou jacobienne)
58+
template<typename T, typename VectorFunc>
59+
tensorium::Tensor<T, 2> autodiff_rank1(
60+
const tensorium::Vector<T>& X,
61+
T dx, T dy, T dz,
62+
VectorFunc&& func,
63+
DiffMode mode,
64+
const tensorium::Tensor<T, 3>& christoffel)
65+
{
66+
using namespace tensorium;
67+
if (mode == DiffMode::PARTIAL) {
68+
return partial_vector(X, dx, dy, dz, std::forward<VectorFunc>(func));
69+
}
70+
else if (mode == DiffMode::COV) {
71+
return covariant_vector(X, dx, dy, dz, std::forward<VectorFunc>(func), christoffel);
72+
}
73+
else {
74+
throw std::invalid_argument("autodiff_rank1: uniquement PARTIAL ou COV autorisés pour un champ vectoriel");
75+
}
76+
}
77+
78+
// ============== RANK 2 (TENSOR 2D) : PREMIÈRE DÉRIVÉE SEULEMENT ==============
79+
// Nous n'implémentons ici que PARTIAL et COV. Si on demande COV2, on lève une exception.
80+
// Return type : Tensor<T,3> (i.e. un tenseur à 3 indices : ∂_k T_{ij} ou ∇_k T_{ij})
81+
template<typename T, typename TensorFunc>
82+
tensorium::Tensor<T, 3> autodiff_rank2_first(
83+
const tensorium::Vector<T>& X,
84+
T dx, T dy, T dz,
85+
TensorFunc&& func,
86+
DiffMode mode,
87+
const tensorium::Tensor<T, 3>& christoffel)
88+
{
89+
using namespace tensorium;
90+
if (mode == DiffMode::PARTIAL) {
91+
return partial_tensor2(X, dx, dy, dz, std::forward<TensorFunc>(func));
92+
}
93+
else if (mode == DiffMode::COV) {
94+
return covariant_tensor2(X, dx, dy, dz, std::forward<TensorFunc>(func), christoffel);
95+
}
96+
else {
97+
throw std::invalid_argument("autodiff_rank2_first: uniquement PARTIAL ou COV autorisés pour un tenseur d'ordre 2");
98+
}
99+
}
100+
101+
// ===================== WRAPPER GÉNÉRIQUE =====================
102+
template<typename T, typename FieldFunc>
103+
auto autodiff(
104+
const tensorium::Vector<T>& X,
105+
T dx, T dy, T dz,
106+
FieldFunc&& func,
107+
DiffMode mode,
108+
const tensorium::Tensor<T, 3>& christoffel = {})
109+
{
110+
using namespace tensorium;
111+
auto sample = func(X);
112+
constexpr size_t rank = TensorTraits<std::decay_t<decltype(sample)>>::rank;
113+
114+
if constexpr (rank == 0) {
115+
return autodiff_rank0(X, dx, dy, dz, std::forward<FieldFunc>(func), mode, christoffel);
116+
}
117+
else if constexpr (rank == 1) {
118+
return autodiff_rank1(X, dx, dy, dz, std::forward<FieldFunc>(func), mode, christoffel);
119+
}
120+
else if constexpr (rank == 2) {
121+
// On autorise ici PARTIAL et COV uniquement pour un champ « Tensor, rank2 ».
122+
// Si on veut une dérivée covariante seconde d'un tenseur 2D, l'utilisateur
123+
// devra appeler explicitement la fonction covariant_tensor2_second(...) par lui-même.
124+
return autodiff_rank2_first(X, dx, dy, dz, std::forward<FieldFunc>(func), mode, christoffel);
125+
}
126+
else {
127+
static_assert(rank <= 2, "autodiff: seuls les rangs 0 (scalaire), 1 (vecteur) et 2 (tenseur2) sont supportés.");
128+
}
129+
}
130+
131+
} // namespace tensorium_RG

0 commit comments

Comments
 (0)