Skip to content

Commit c4bcec4

Browse files
at0m741at0m741
committed
tilde Ricci tensor, now I need the Ricci tensor that depend on the conformal factor \chi to reconstruct the entire Ricci
1 parent 4d942a1 commit c4bcec4

File tree

4 files changed

+329
-218
lines changed

4 files changed

+329
-218
lines changed

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

Lines changed: 50 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -31,30 +31,49 @@ namespace tensorium_RG {
3131

3232
enum class DiffMode { PARTIAL, COV, COV2 };
3333

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-
}
5534

56-
// ===================== RANK 1 (VECTOR) =====================
57-
// Return type : Tensor<T,2> (i.e. une matrice dérivée covariante ou jacobienne)
35+
template<typename T, typename ScalarFunc>
36+
tensorium::Vector<T> covariant_scalar(
37+
const tensorium::Vector<T>& X,
38+
T dx, T dy, T dz,
39+
ScalarFunc&& func,
40+
const tensorium::Tensor<T, 3>& christoffel)
41+
{
42+
tensorium::Vector<T> grad = partial_scalar(X, dx, dy, dz, std::forward<ScalarFunc>(func));
43+
return grad;
44+
}
45+
46+
template<typename T, typename ScalarFunc>
47+
tensorium::Tensor<T, 2> autodiff_scalar_second(
48+
const tensorium::Vector<T>& X,
49+
T dx, T dy, T dz,
50+
ScalarFunc&& func,
51+
const tensorium::Tensor<T, 3>& christoffel)
52+
{
53+
return covariant_scalar_second(X, dx, dy, dz, std::forward<ScalarFunc>(func), christoffel);
54+
}
55+
56+
template<typename T, typename ScalarFunc>
57+
tensorium::Vector<T> autodiff_rank0(
58+
const tensorium::Vector<T>& X,
59+
T dx, T dy, T dz,
60+
ScalarFunc&& func,
61+
DiffMode mode,
62+
const tensorium::Tensor<T, 3>& christoffel)
63+
{
64+
using namespace tensorium;
65+
if (mode == DiffMode::PARTIAL) {
66+
return partial_scalar(X, dx, dy, dz, std::forward<ScalarFunc>(func));
67+
}
68+
else if (mode == DiffMode::COV) {
69+
return covariant_scalar(X, dx, dy, dz, std::forward<ScalarFunc>(func), christoffel);
70+
}
71+
else {
72+
throw std::invalid_argument("autodiff_rank0: PARTIAL ou COV uniquement. Utiliser `autodiff_scalar_second` pour COV2.");
73+
}
74+
}
75+
76+
5877
template<typename T, typename VectorFunc>
5978
tensorium::Tensor<T, 2> autodiff_rank1(
6079
const tensorium::Vector<T>& X,
@@ -75,9 +94,6 @@ namespace tensorium_RG {
7594
}
7695
}
7796

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})
8197
template<typename T, typename TensorFunc>
8298
tensorium::Tensor<T, 3> autodiff_rank2_first(
8399
const tensorium::Vector<T>& X,
@@ -98,7 +114,6 @@ namespace tensorium_RG {
98114
}
99115
}
100116

101-
// ===================== WRAPPER GÉNÉRIQUE =====================
102117
template<typename T, typename FieldFunc>
103118
auto autodiff(
104119
const tensorium::Vector<T>& X,
@@ -109,18 +124,18 @@ namespace tensorium_RG {
109124
{
110125
using namespace tensorium;
111126
auto sample = func(X);
112-
constexpr size_t rank = TensorTraits<std::decay_t<decltype(sample)>>::rank;
127+
128+
constexpr size_t rank = TensorTraits<std::decay_t<decltype(sample)>>::rank;
113129

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) {
130+
if constexpr (rank == 0) {
131+
if (mode == DiffMode::COV2)
132+
return autodiff_scalar_second(X, dx, dy, dz, std::forward<FieldFunc>(func), christoffel); // <- nouveau chemin
133+
else
134+
return autodiff_rank0(X, dx, dy, dz, std::forward<FieldFunc>(func), mode, christoffel);
135+
} else if constexpr (rank == 1) {
118136
return autodiff_rank1(X, dx, dy, dz, std::forward<FieldFunc>(func), mode, christoffel);
119137
}
120138
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.
124139
return autodiff_rank2_first(X, dx, dy, dz, std::forward<FieldFunc>(func), mode, christoffel);
125140
}
126141
else {

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

Lines changed: 87 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,93 @@
66
namespace tensorium_RG {
77

88

9+
10+
template<typename T, typename ScalarFunc>
11+
tensorium::Vector<T> partial_scalar(
12+
const tensorium::Vector<T>& X,
13+
T dx, T dy, T dz,
14+
ScalarFunc&& func)
15+
{
16+
tensorium::Vector<T> grad(3); // <-- ici 3, pas 4
17+
18+
// dérivée selon x¹
19+
{
20+
tensorium::Vector<T> Xs = X;
21+
T gm2, gm1, gp1, gp2;
22+
Xs(1) = X(1) - 2*dx; gm2 = func(Xs);
23+
Xs(1) = X(1) - dx; gm1 = func(Xs);
24+
Xs(1) = X(1) + dx; gp1 = func(Xs);
25+
Xs(1) = X(1) + 2*dx; gp2 = func(Xs);
26+
grad(0) = (-gp2 + 8*gp1 - 8*gm1 + gm2) / (12*dx);
27+
}
28+
// dérivée selon x²
29+
{
30+
tensorium::Vector<T> Xs = X;
31+
T gm2, gm1, gp1, gp2;
32+
Xs(2) = X(2) - 2*dy; gm2 = func(Xs);
33+
Xs(2) = X(2) - dy; gm1 = func(Xs);
34+
Xs(2) = X(2) + dy; gp1 = func(Xs);
35+
Xs(2) = X(2) + 2*dy; gp2 = func(Xs);
36+
grad(1) = (-gp2 + 8*gp1 - 8*gm1 + gm2) / (12*dy);
37+
}
38+
// dérivée selon x³
39+
{
40+
tensorium::Vector<T> Xs = X;
41+
T gm2, gm1, gp1, gp2;
42+
Xs(3) = X(3) - 2*dz; gm2 = func(Xs);
43+
Xs(3) = X(3) - dz; gm1 = func(Xs);
44+
Xs(3) = X(3) + dz; gp1 = func(Xs);
45+
Xs(3) = X(3) + 2*dz; gp2 = func(Xs);
46+
grad(2) = (-gp2 + 8*gp1 - 8*gm1 + gm2) / (12*dz);
47+
}
48+
49+
return grad;
50+
}
51+
52+
53+
template<typename T, typename VectorFunc>
54+
tensorium::Tensor<T,2> partial_vector(
55+
const tensorium::Vector<T>& X,
56+
T dx, T dy, T dz,
57+
VectorFunc&& func)
58+
{
59+
tensorium::Tensor<T,2> result({3,3});
60+
// dérivée ∂/∂x
61+
{
62+
tensorium::Vector<T> Xs = X;
63+
tensorium::Vector<T> Vm2(3), Vm1(3), Vp1(3), Vp2(3);
64+
Xs(1) = X(1) - 2*dx; Vm2 = func(Xs);
65+
Xs(1) = X(1) - dx; Vm1 = func(Xs);
66+
Xs(1) = X(1) + dx; Vp1 = func(Xs);
67+
Xs(1) = X(1) + 2*dx; Vp2 = func(Xs);
68+
for(int i=0;i<3;++i)
69+
result(i,0) = (-Vp2(i) + 8*Vp1(i) - 8*Vm1(i) + Vm2(i)) / (12*dx);
70+
}
71+
// dérivée ∂/∂y
72+
{
73+
tensorium::Vector<T> Xs = X;
74+
tensorium::Vector<T> Vm2(3), Vm1(3), Vp1(3), Vp2(3);
75+
Xs(2) = X(2) - 2*dy; Vm2 = func(Xs);
76+
Xs(2) = X(2) - dy; Vm1 = func(Xs);
77+
Xs(2) = X(2) + dy; Vp1 = func(Xs);
78+
Xs(2) = X(2) + 2*dy; Vp2 = func(Xs);
79+
for(int i=0;i<3;++i)
80+
result(i,1) = (-Vp2(i) + 8*Vp1(i) - 8*Vm1(i) + Vm2(i)) / (12*dy);
81+
}
82+
// dérivée ∂/∂z
83+
{
84+
tensorium::Vector<T> Xs = X;
85+
tensorium::Vector<T> Vm2(3), Vm1(3), Vp1(3), Vp2(3);
86+
Xs(3) = X(3) - 2*dz; Vm2 = func(Xs);
87+
Xs(3) = X(3) - dz; Vm1 = func(Xs);
88+
Xs(3) = X(3) + dz; Vp1 = func(Xs);
89+
Xs(3) = X(3) + 2*dz; Vp2 = func(Xs);
90+
for(int i=0;i<3;++i)
91+
result(i,2) = (-Vp2(i) + 8*Vp1(i) - 8*Vm1(i) + Vm2(i)) / (12*dz);
92+
}
93+
return result;
94+
}
95+
996
template<typename T, typename TensorFunc>
1097
void compute_partial_derivatives_tensor2D(
1198
const tensorium::Vector<T>& X,
@@ -457,50 +544,6 @@ namespace tensorium_RG {
457544
return result;
458545
}
459546

460-
template<typename T, typename ScalarFunc>
461-
tensorium::Vector<T> partial_scalar(
462-
const tensorium::Vector<T>& X,
463-
T dx, T dy, T dz,
464-
ScalarFunc&& func)
465-
{
466-
tensorium::Vector<T> grad(3);
467-
468-
T fxm = func(X - tensorium::Vector<T>{dx, 0, 0});
469-
T fxp = func(X + tensorium::Vector<T>{dx, 0, 0});
470-
grad(0) = (fxp - fxm) / (2 * dx);
471-
472-
T fym = func(X - tensorium::Vector<T>{0, dy, 0});
473-
T fyp = func(X + tensorium::Vector<T>{0, dy, 0});
474-
grad(1) = (fyp - fym) / (2 * dy);
475-
476-
T fzm = func(X - tensorium::Vector<T>{0, 0, dz});
477-
T fzp = func(X + tensorium::Vector<T>{0, 0, dz});
478-
grad(2) = (fzp - fzm) / (2 * dz);
479-
480-
return grad;
481-
}
482-
483-
template<typename T, typename VectorFunc>
484-
tensorium::Tensor<T, 2> partial_vector(
485-
const tensorium::Vector<T>& X,
486-
T dx, T dy, T dz,
487-
VectorFunc&& func)
488-
{
489-
tensorium::Tensor<T, 2> result({3, 3});
490-
491-
for (int i = 0; i < 3; ++i) {
492-
tensorium::Vector<T> dx_vec = {0, 0, 0};
493-
dx_vec(i) = (i == 0) ? dx : (i == 1 ? dy : dz);
494-
495-
tensorium::Vector<T> fm = func(X - dx_vec);
496-
tensorium::Vector<T> fp = func(X + dx_vec);
497-
498-
for (int j = 0; j < 3; ++j)
499-
result(j, i) = (fp(j) - fm(j)) / (2 * dx_vec(i));
500-
}
501-
return result;
502-
}
503-
504547
template<typename T, typename TensorFunc>
505548
tensorium::Tensor<T, 3> partial_tensor2(
506549
const tensorium::Vector<T>& X,

0 commit comments

Comments
 (0)