Skip to content

Commit 1a5e7c2

Browse files
at0m741at0m741
committed
Ricci conformal tensor
1 parent bcaee7a commit 1a5e7c2

File tree

5 files changed

+294
-72
lines changed

5 files changed

+294
-72
lines changed

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

Lines changed: 33 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -27,59 +27,68 @@ namespace tensorium_RG {
2727
const tensorium::Tensor<T, 5>& g_tilde_inv,
2828
T dx, T dy, T dz,
2929
int max_iter = 1000,
30-
T tol = 1e-8) {
31-
32-
const size_t NX = Atilde.shape()[0];
33-
const size_t NY = Atilde.shape()[1];
34-
const size_t NZ = Atilde.shape()[2];
35-
36-
tensorium::Tensor<T, 3> psi({NX, NY, NZ});
37-
std::fill(psi.begin(), psi.end(), T(1));;
30+
T tol = T(1e-8))
31+
{
32+
std::array<size_t,5> shp = Atilde.shape();
33+
const size_t NX = shp[0];
34+
const size_t NY = shp[1];
35+
const size_t NZ = shp[2];
36+
37+
tensorium::Tensor<T,3> psi({NX, NY, NZ});
38+
for (size_t i = 0; i < NX; ++i) {
39+
for (size_t j = 0; j < NY; ++j) {
40+
for (size_t k = 0; k < NZ; ++k) {
41+
psi({i,j,k}) = T(1);
42+
}
43+
}
44+
}
3845

3946
const T inv_dx2 = T(1) / (dx * dx);
4047
const T inv_dy2 = T(1) / (dy * dy);
4148
const T inv_dz2 = T(1) / (dz * dz);
42-
const T factor = T(1) / (2 * (inv_dx2 + inv_dy2 + inv_dz2));
49+
const T factor = T(1) / ( T(2) * (inv_dx2 + inv_dy2 + inv_dz2) );
4350

4451
for (int it = 0; it < max_iter; ++it) {
4552
T max_diff = T(0);
4653

47-
#pragma omp parallel for collapse(3) reduction(max:max_diff)
4854
for (size_t i = 1; i < NX - 1; ++i) {
4955
for (size_t j = 1; j < NY - 1; ++j) {
50-
for (size_t k = 1; k < NZ - 1; ++k) {
51-
52-
T lap = inv_dx2 * (psi(i + 1, j, k) + psi(i - 1, j, k) - 2 * psi(i, j, k))
53-
+ inv_dy2 * (psi(i, j + 1, k) + psi(i, j - 1, k) - 2 * psi(i, j, k))
54-
+ inv_dz2 * (psi(i, j, k + 1) + psi(i, j, k - 1) - 2 * psi(i, j, k));
56+
for (size_t k = 1; k < NZ - 1; ++k) {
57+
T lap = inv_dx2 * ( psi({i+1,j,k}) + psi({i-1,j,k}) - T(2)*psi({i,j,k}) )
58+
+ inv_dy2 * ( psi({i,j+1,k}) + psi({i,j-1,k}) - T(2)*psi({i,j,k}) )
59+
+ inv_dz2 * ( psi({i,j,k+1}) + psi({i,j,k-1}) - T(2)*psi({i,j,k}) );
5560

56-
// Compute ò = γ̃^{ia} γ̃^{jb} Ã_{ab} Ã_{ij}
5761
T A2 = T(0);
5862
for (int a = 0; a < 3; ++a) {
5963
for (int b = 0; b < 3; ++b) {
6064
for (int c = 0; c < 3; ++c) {
6165
for (int d = 0; d < 3; ++d) {
62-
A2 += g_tilde_inv(i,j,k,a,c) * g_tilde_inv(i,j,k,b,d)
63-
* Atilde(i,j,k,c,d) * Atilde(i,j,k,a,b);
66+
// accès à g_tilde_inv et Atilde par array d’indices
67+
T g1 = g_tilde_inv({i,j,k,(size_t)a,(size_t)c});
68+
T g2 = g_tilde_inv({i,j,k,(size_t)b,(size_t)d});
69+
T A3 = Atilde({i,j,k,(size_t)c,(size_t)d});
70+
T A4 = Atilde({i,j,k,(size_t)a,(size_t)b});
71+
A2 += g1 * g2 * A3 * A4;
6472
}
6573
}
6674
}
6775
}
6876

69-
T psi_val = std::fmax(psi(i, j, k), T(1e-8));
70-
T rhs = (T(1.0) / 8.0) * A2 / std::pow(psi_val, 7);
77+
T psi_ijk = std::fmax( psi({i,j,k}), T(1e-8) );
78+
T rhs = (T(1)/T(8)) * A2 / std::pow(psi_ijk, T(7));
7179

72-
T new_psi = psi(i, j, k) + T(0.5) * (lap - rhs) * factor;
80+
T new_psi = psi({i,j,k}) + T(0.5)*(lap - rhs)*factor;
7381
new_psi = std::fmax(new_psi, T(0.1));
7482

75-
max_diff = std::fmax(max_diff, std::abs(new_psi - psi(i, j, k)));
76-
psi(i, j, k) = new_psi;
83+
max_diff = std::fmax( max_diff, std::abs(new_psi - psi({i,j,k})) );
84+
psi({i,j,k}) = new_psi;
7785
}
7886
}
7987
}
8088

81-
if (max_diff < tol)
89+
if (max_diff < tol) {
8290
break;
91+
}
8392
}
8493

8594
return psi;

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

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,42 @@ namespace tensorium_RG {
4848

4949
return out;
5050
}
51+
52+
void compute_contracted_christoffel(
53+
const tensorium::Vector<T>& X,
54+
T dx, T dy, T dz,
55+
const tensorium_RG::Metric<T>& metric,
56+
tensorium::Tensor<T, 5>& dGamma_contract
57+
) {
58+
auto shape = dGamma_contract.shape();
59+
const size_t NX = shape[0], NY = shape[1], NZ = shape[2];
60+
61+
for (size_t i = 0; i < NX; ++i)
62+
for (size_t j = 0; j < NY; ++j)
63+
for (size_t k = 0; k < NZ; ++k) {
64+
tensorium::Vector<T> Xs = {
65+
X(0) + dx * (i - NX/2),
66+
X(1) + dy * (j - NY/2),
67+
X(2) + dz * (k - NZ/2)
68+
};
69+
70+
// 1. get metric at Xs
71+
T alpha;
72+
tensorium::Vector<T> beta(3);
73+
tensorium::Tensor<T,2> gamma({3,3});
74+
metric.BSSN(Xs, alpha, beta, gamma);
75+
76+
T chi = compute_conformal_factor(metric, gamma);
77+
78+
tensorium::Vector<T> gamma_contract = tensorium_RG::BSSNContractedGamma<T>::compute(
79+
Xs, metric, dx, dy, dz, chi
80+
);
81+
82+
for (size_t a = 0; a < 3; ++a)
83+
for (size_t b = 0; b < 3; ++b)
84+
dGamma_contract(i,j,k,a,b) = (a == b) ? gamma_contract[a] : T(0);
85+
}
86+
}
5187
};
5288

5389
} // namespace tensorium_RG

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

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,53 @@ namespace tensorium_RG {
4848
}
4949
}
5050

51+
template<typename T, typename TensorFunc>
52+
void compute_second_derivatives_tensor2D(
53+
const tensorium::Vector<T>& X,
54+
T dx, T dy, T dz,
55+
TensorFunc&& func,
56+
tensorium::Tensor<T, 4>& out
57+
) {
58+
out = tensorium::Tensor<T, 4>({3, 3, 3, 3});
59+
60+
auto shifted = [&](T dx_, T dy_, T dz_) {
61+
tensorium::Vector<T> Xs = X;
62+
Xs(0) += dx_; Xs(1) += dy_; Xs(2) += dz_;
63+
tensorium::Tensor<T, 2> out_tensor({3, 3});
64+
func(Xs, out_tensor);
65+
return out_tensor;
66+
};
67+
68+
for (size_t a = 0; a < 3; ++a) {
69+
for (size_t b = 0; b < 3; ++b) {
70+
// ∂²/∂x²
71+
auto gm2 = shifted(-2 * dx, 0, 0);
72+
auto gm1 = shifted(-dx, 0, 0);
73+
auto g0 = shifted(0, 0, 0);
74+
auto gp1 = shifted(dx, 0, 0);
75+
auto gp2 = shifted(2 * dx, 0, 0);
76+
out(0, a, b, 0) = (-gp2(a,b) + 16*gp1(a,b) - 30*g0(a,b) + 16*gm1(a,b) - gm2(a,b)) / (12 * dx * dx);
77+
78+
// ∂²/∂y²
79+
gm2 = shifted(0, -2 * dy, 0);
80+
gm1 = shifted(0, -dy, 0);
81+
g0 = shifted(0, 0, 0);
82+
gp1 = shifted(0, dy, 0);
83+
gp2 = shifted(0, 2 * dy, 0);
84+
out(1, a, b, 1) = (-gp2(a,b) + 16*gp1(a,b) - 30*g0(a,b) + 16*gm1(a,b) - gm2(a,b)) / (12 * dy * dy);
85+
86+
// ∂²/∂z²
87+
gm2 = shifted(0, 0, -2 * dz);
88+
gm1 = shifted(0, 0, -dz);
89+
g0 = shifted(0, 0, 0);
90+
gp1 = shifted(0, 0, dz);
91+
gp2 = shifted(0, 0, 2 * dz);
92+
out(2, a, b, 2) = (-gp2(a,b) + 16*gp1(a,b) - 30*g0(a,b) + 16*gm1(a,b) - gm2(a,b)) / (12 * dz * dz);
93+
}
94+
}
95+
}
96+
97+
5198
template<typename T, typename VectorFunc>
5299
void compute_partial_derivatives_vector(
53100
const tensorium::Vector<T>& X,
@@ -140,5 +187,46 @@ namespace tensorium_RG {
140187
}
141188
return dtg;
142189
}
190+
191+
template<typename T>
192+
inline void compute_partial_derivatives_vector3D(
193+
const tensorium::Tensor<T, 4>& vec_field,
194+
size_t i, size_t j, size_t k,
195+
T dx, T dy, T dz,
196+
tensorium::Tensor<T, 2>& dvec_out
197+
) {
198+
dvec_out.resize(3, 3);
199+
200+
auto get = [&](int ii, int jj, int kk, int a) -> T {
201+
return vec_field(ii, jj, kk, a);
202+
};
203+
204+
const auto& shape = vec_field.shape();
205+
const size_t NX = shape[0], NY = shape[1], NZ = shape[2];
206+
207+
for (size_t a = 0; a < 3; ++a) {
208+
if (i >= 2 && i + 2 < NX)
209+
dvec_out(0, a) = (-get(i+2,j,k,a) + 8*get(i+1,j,k,a) - 8*get(i-1,j,k,a) + get(i-2,j,k,a)) / (12 * dx);
210+
else if (i > 0 && i + 1 < NX)
211+
dvec_out(0, a) = (get(i+1,j,k,a) - get(i-1,j,k,a)) / (2 * dx);
212+
else
213+
dvec_out(0, a) = T(0);
214+
215+
if (j >= 2 && j + 2 < NY)
216+
dvec_out(1, a) = (-get(i,j+2,k,a) + 8*get(i,j+1,k,a) - 8*get(i,j-1,k,a) + get(i,j-2,k,a)) / (12 * dy);
217+
else if (j > 0 && j + 1 < NY)
218+
dvec_out(1, a) = (get(i,j+1,k,a) - get(i,j-1,k,a)) / (2 * dy);
219+
else
220+
dvec_out(1, a) = T(0);
221+
222+
if (k >= 2 && k + 2 < NZ)
223+
dvec_out(2, a) = (-get(i,j,k+2,a) + 8*get(i,j,k+1,a) - 8*get(i,j,k-1,a) + get(i,j,k-2,a)) / (12 * dz);
224+
else if (k > 0 && k + 1 < NZ)
225+
dvec_out(2, a) = (get(i,j,k+1,a) - get(i,j,k-1,a)) / (2 * dz);
226+
else
227+
dvec_out(2, a) = T(0);
228+
}
229+
}
230+
143231
}
144232

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

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
#pragma once
2+
3+
#include "../Metric.hpp"
4+
#include "BSSNDerivatives.hpp"
5+
#include "../../Core/Tensor.hpp"
6+
#include "../../Core/Matrix.hpp"
7+
#include "../../Core/Vector.hpp"
8+
#include "BSSNChristoffel.hpp"
9+
#include "BSSNDerivatives.hpp"
10+
11+
namespace tensorium_RG {
12+
13+
template<typename T>
14+
class RicciTensor3D {
15+
public:
16+
17+
static tensorium::Tensor<T, 5> compute_ricci_conformal(
18+
const tensorium::Tensor<T, 5>& gamma_tilde,
19+
const tensorium::Tensor<T, 5>& gamma_tilde_inv,
20+
const tensorium::Tensor<T, 6>& christoffel_tilde,
21+
const tensorium::Tensor<T, 5>& dGamma_contract,
22+
const tensorium::Tensor<T, 7>& d2gamma_tilde) {
23+
24+
using tensorium::Tensor;
25+
26+
const auto& shape = gamma_tilde.shape();
27+
const size_t NX = shape[0], NY = shape[1], NZ = shape[2];
28+
29+
Tensor<T, 5> Ricci_tilde({NX, NY, NZ, 3, 3});
30+
Ricci_tilde.fill(T(0));
31+
32+
for (size_t x = 1; x + 1 < NX; ++x) {
33+
for (size_t y = 1; y + 1 < NY; ++y) {
34+
for (size_t z = 1; z + 1 < NZ; ++z) {
35+
36+
for (size_t i = 0; i < 3; ++i) {
37+
for (size_t j = 0; j < 3; ++j) {
38+
39+
T term_laplacian = T(0);
40+
for (size_t k = 0; k < 3; ++k) {
41+
for (size_t l = 0; l < 3; ++l) {
42+
T ginv_kl = gamma_tilde_inv({x, y, z, k, l});
43+
T d2_kl = d2gamma_tilde({x, y, z, i, j, k, l});
44+
term_laplacian += ginv_kl * d2_kl;
45+
}
46+
}
47+
term_laplacian *= T(-0.5);
48+
49+
T term_divGamma = T(0);
50+
for (size_t m = 0; m < 3; ++m) {
51+
T dG_j = dGamma_contract({x, y, z, m, j});
52+
T dG_i = dGamma_contract({x, y, z, m, i});
53+
T g_mi = gamma_tilde({x, y, z, m, i});
54+
T g_mj = gamma_tilde({x, y, z, m, j});
55+
term_divGamma += g_mi * dG_j + g_mj * dG_i;
56+
}
57+
58+
T term_gamma_gamma = T(0);
59+
for (size_t m = 0; m < 3; ++m) {
60+
for (size_t l = 0; l < 3; ++l) {
61+
T g1 = christoffel_tilde({x, y, z, m, l, i});
62+
T g2 = christoffel_tilde({x, y, z, l, m, j});
63+
term_gamma_gamma += g1 * g2;
64+
}
65+
}
66+
67+
Ricci_tilde({x, y, z, i, j}) = term_laplacian + term_divGamma + term_gamma_gamma;
68+
}
69+
}
70+
}
71+
}
72+
}
73+
74+
return Ricci_tilde;
75+
}
76+
77+
};
78+
79+
}

0 commit comments

Comments
 (0)