Skip to content

Commit 9631c92

Browse files
at0m741at0m741
committed
meh 3D grid init skill issue
1 parent 33ab810 commit 9631c92

File tree

9 files changed

+203
-69
lines changed

9 files changed

+203
-69
lines changed

‎Tests/Matrix/Matrix.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -343,8 +343,9 @@ int matrix_tests() {
343343
gammaj_inv.print();
344344

345345
std::cout << "--- SETUP BSSN TEST---\n";
346-
auto bssn = tensorium::setup_BSSN_grid(X, metric, dx, dy, dz);
346+
auto bssn = tensorium::setup_BSSN(X, metric, dx, dy, dz);
347347

348+
auto bssn3D = tensorium::setup_BSSN_grid(X, metric, dx, dy, dz);
348349
return 0;
349350
}
350351

‎includes/Tensorium/Core/Spectral.hpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -139,6 +139,13 @@ namespace tensorium {
139139
a(i,j,k) = sliceZ(k);
140140
}
141141
}
142+
const T norm = 1.0 / (NX * NY * NZ);
143+
#pragma omp parallel for collapse(3)
144+
for (size_t i = 0; i < NX; ++i)
145+
for (size_t j = 0; j < NY; ++j)
146+
for (size_t k = 0; k < NZ; ++k)
147+
a(i, j, k) *= norm;
148+
142149
}
143150

144151
private:

‎includes/Tensorium/Core/Tensor.hpp‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -68,12 +68,18 @@ namespace tensorium {
6868
}
6969
/** @brief Resize 2D tensor */
7070
///@{
71-
void resize(const std::array<size_t, 2>& dims) {
71+
void resize(const std::array<size_t, Rank>& dims) {
7272
dimensions = dims;
7373
update_strides();
74-
total_size = dims[0] * dims[1];
75-
data.resize(total_size);
74+
75+
size_t total = 1;
76+
for (size_t i = 0; i < Rank; ++i)
77+
total *= dims[i];
78+
79+
total_size = total;
80+
data.resize(total);
7681
}
82+
7783
/** @brief Resize 2D tensor */
7884
void resize(size_t d0, size_t d1) {
7985
resize(std::array<size_t, 2>{d0, d1});

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

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ namespace tensorium {
2929

3030
namespace tensorium_RG {
3131

32-
enum class DiffMode { PARTIAL, COV, COV2 };
32+
enum class DiffMode { PARTIAL, COV, COV2, SPEC };
3333

3434

3535
template<typename T, typename ScalarFunc>
@@ -108,8 +108,9 @@ namespace tensorium_RG {
108108
}
109109
else if (mode == DiffMode::COV) {
110110
return covariant_tensor2(X, dx, dy, dz, std::forward<TensorFunc>(func), christoffel);
111-
}
112-
else {
111+
}else if (mode == DiffMode::SPEC) {
112+
return spectral_partial_tensor2<T>(X, dx, dy, dz, std::forward<TensorFunc>(func), 64, 64, 64);
113+
} else {
113114
throw std::invalid_argument("autodiff_rank2_first: uniquement PARTIAL ou COV autorisés pour un tenseur d'ordre 2");
114115
}
115116
}

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

Lines changed: 55 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,46 @@ namespace tensorium_RG {
6565
}
6666
}
6767

68+
static void compute3D(
69+
const tensorium::Tensor<T, 5>& gamma_tilde, // [NX, NY, NZ, 3, 3]
70+
const tensorium::Tensor<T, 6>& dgamma_tilde, // [NX, NY, NZ, 3, 3, 3]
71+
const tensorium::Tensor<T, 5>& gamma_tilde_inv, // [NX, NY, NZ, 3, 3]
72+
tensorium::Tensor<T, 6>& Christoffel // [NX, NY, NZ, 3, 3, 3] (output)
73+
) {
74+
using namespace tensorium;
75+
76+
const size_t NX = gamma_tilde.dim(0);
77+
const size_t NY = gamma_tilde.dim(1);
78+
const size_t NZ = gamma_tilde.dim(2);
79+
80+
Christoffel.resize(NX, NY, NZ, 3, 3, 3);
81+
82+
for (size_t x = 0; x < NX; ++x) {
83+
for (size_t y = 0; y < NY; ++y) {
84+
for (size_t z = 0; z < NZ; ++z) {
85+
Tensor<T, 2> gtilde(3, 3);
86+
Tensor<T, 2> ginv(3, 3);
87+
Tensor<T, 3> dg(3, 3, 3);
88+
Tensor<T, 3> chr(3, 3, 3);
89+
90+
for (size_t i = 0; i < 3; ++i)
91+
for (size_t j = 0; j < 3; ++j) {
92+
gtilde(i, j) = gamma_tilde(x, y, z, i, j);
93+
ginv(i, j) = gamma_tilde_inv(x, y, z, i, j);
94+
for (size_t k = 0; k < 3; ++k)
95+
dg(i, j, k) = dgamma_tilde(x, y, z, i, j, k);
96+
}
97+
98+
BSSNChristoffel<T>::compute(gtilde, dg, ginv, chr);
99+
100+
for (size_t i = 0; i < 3; ++i)
101+
for (size_t j = 0; j < 3; ++j)
102+
for (size_t k = 0; k < 3; ++k)
103+
Christoffel(x, y, z, i, j, k) = chr(i, j, k);
104+
}
105+
}
106+
}
107+
}
68108

69109
};
70110
/**
@@ -91,25 +131,25 @@ namespace tensorium_RG {
91131
const tensorium::Tensor<T, 5>& gamma_field,
92132
size_t i, size_t j, size_t k,
93133
T dx, T dy, T dz,
94-
tensorium::Tensor<T, 3>& dgamma_out) {
134+
tensorium::Tensor<T, 3>& dgamma_out) {
95135

96-
dgamma_out.resize(3, 3, 3);
136+
dgamma_out.resize(3, 3, 3);
97137

98-
auto get = [&](int ii, int jj, int kk, int a, int b) -> T {
99-
return gamma_field(ii, jj, kk, a, b);
100-
};
138+
auto get = [&](int ii, int jj, int kk, int a, int b) -> T {
139+
return gamma_field(ii, jj, kk, a, b);
140+
};
101141

102-
auto dim = gamma_field.shape();
103-
size_t Nx = dim[0], Ny = dim[1], Nz = dim[2];
142+
auto dim = gamma_field.shape();
143+
size_t Nx = dim[0], Ny = dim[1], Nz = dim[2];
104144

105-
for (size_t a = 0; a < 3; ++a) {
106-
for (size_t b = 0; b < 3; ++b) {
107-
if (i >= 2 && i + 2 < Nx)
108-
dgamma_out(0, a, b) = (-get(i+2,j,k,a,b) + 8*get(i+1,j,k,a,b) - 8*get(i-1,j,k,a,b) + get(i-2,j,k,a,b)) / (12 * dx);
109-
else if (i > 0 && i + 1 < Nx)
110-
dgamma_out(0, a, b) = (get(i+1,j,k,a,b) - get(i-1,j,k,a,b)) / (2 * dx);
111-
else
112-
dgamma_out(0, a, b) = T(0);
145+
for (size_t a = 0; a < 3; ++a) {
146+
for (size_t b = 0; b < 3; ++b) {
147+
if (i >= 2 && i + 2 < Nx)
148+
dgamma_out(0, a, b) = (-get(i+2,j,k,a,b) + 8*get(i+1,j,k,a,b) - 8*get(i-1,j,k,a,b) + get(i-2,j,k,a,b)) / (12 * dx);
149+
else if (i > 0 && i + 1 < Nx)
150+
dgamma_out(0, a, b) = (get(i+1,j,k,a,b) - get(i-1,j,k,a,b)) / (2 * dx);
151+
else
152+
dgamma_out(0, a, b) = T(0);
113153

114154
if (j >= 2 && j + 2 < Ny)
115155
dgamma_out(1, a, b) = (-get(i,j+2,k,a,b) + 8*get(i,j+1,k,a,b) - 8*get(i,j-1,k,a,b) + get(i,j-2,k,a,b)) / (12 * dy);

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

Lines changed: 48 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
#include "../Metric.hpp"
44
#include "../../Core/Spectral.hpp"
5-
5+
#include "BSSNGridSetupUtils.hpp"
66
namespace tensorium_RG {
77

88

@@ -33,73 +33,104 @@ namespace tensorium_RG {
3333
void spectral_partial_scalar_3D(
3434
const tensorium::Tensor<T, 3>& scalar_field,
3535
T dx, T dy, T dz,
36-
tensorium::Tensor<T, 4>& grad_out) {
37-
36+
tensorium::Tensor<T, 4>& grad_out)
37+
{
3838
using namespace tensorium;
3939
using C = std::complex<T>;
4040
using FFT = SpectralFFT<T>;
4141

4242
const auto& shape = scalar_field.shape();
4343
const size_t NX = shape[0], NY = shape[1], NZ = shape[2];
44-
const size_t N = NX * NY * NZ;
4544

46-
grad_out.resize(NX, NY, NZ, 3);
45+
grad_out.resize({NX, NY, NZ, 3});
4746

4847
Tensor<C, 3> field_cplx({NX, NY, NZ});
4948
for (size_t i = 0; i < NX; ++i)
5049
for (size_t j = 0; j < NY; ++j)
5150
for (size_t k = 0; k < NZ; ++k)
52-
field_cplx(i,j,k) = C(scalar_field(i,j,k), 0.0);
51+
field_cplx(i, j, k) = C(scalar_field(i, j, k), 0.0);
5352

5453
FFT::forward_3D(field_cplx);
5554

5655
for (size_t i = 0; i < NX; ++i) {
57-
int ni = (i <= NX/2) ? i : int(i) - int(NX);
56+
int ni = (i <= NX / 2) ? int(i) : int(i) - int(NX);
5857
T kx = 2.0 * M_PI * ni / (NX * dx);
5958

6059
for (size_t j = 0; j < NY; ++j) {
61-
int nj = (j <= NY/2) ? j : int(j) - int(NY);
60+
int nj = (j <= NY / 2) ? int(j) : int(j) - int(NY);
6261
T ky = 2.0 * M_PI * nj / (NY * dy);
6362

6463
for (size_t k = 0; k < NZ; ++k) {
65-
int nk = (k <= NZ/2) ? k : int(k) - int(NZ);
64+
int nk = (k <= NZ / 2) ? int(k) : int(k) - int(NZ);
6665
T kz = 2.0 * M_PI * nk / (NZ * dz);
6766

68-
C f_hat = field_cplx(i,j,k);
67+
C f_hat = field_cplx(i, j, k);
6968

70-
grad_out(i,j,k,0) = (f_hat * C(0, kx)).real();
71-
grad_out(i,j,k,1) = (f_hat * C(0, ky)).real();
72-
grad_out(i,j,k,2) = (f_hat * C(0, kz)).real();
69+
field_cplx(i, j, k) = f_hat;
70+
71+
grad_out(i, j, k, 0) = (f_hat * C(0, kx)).real();
72+
grad_out(i, j, k, 1) = (f_hat * C(0, ky)).real();
73+
grad_out(i, j, k, 2) = (f_hat * C(0, kz)).real();
7374
}
7475
}
7576
}
7677

7778
for (int dim = 0; dim < 3; ++dim) {
7879
Tensor<C, 3> dfield_cplx({NX, NY, NZ});
80+
7981
for (size_t i = 0; i < NX; ++i)
8082
for (size_t j = 0; j < NY; ++j)
8183
for (size_t k = 0; k < NZ; ++k)
82-
dfield_cplx(i,j,k) = C(grad_out(i,j,k,dim), 0.0);
84+
dfield_cplx(i, j, k) = C(grad_out(i, j, k, dim), 0.0);
8385

8486
FFT::backward_3D(dfield_cplx);
8587

8688
for (size_t i = 0; i < NX; ++i)
8789
for (size_t j = 0; j < NY; ++j)
8890
for (size_t k = 0; k < NZ; ++k)
89-
grad_out(i,j,k,dim) = dfield_cplx(i,j,k).real();
91+
grad_out(i, j, k, dim) = dfield_cplx(i, j, k).real();
9092
}
9193
}
9294

95+
template<typename T, typename TensorFunc>
96+
tensorium::Tensor<T, 3> spectral_partial_tensor2(
97+
const tensorium::Vector<T>& X,
98+
T dx, T dy, T dz,
99+
TensorFunc&& func,
100+
size_t NX, size_t NY, size_t NZ)
101+
{
102+
using namespace tensorium;
103+
using C = std::complex<T>;
104+
using FFT = SpectralFFT<T>;
105+
106+
Tensor<T, 2> gamma = func(X);
107+
Tensor<T, 3> grad_gamma({3, 3, 3});
108+
109+
for (size_t i = 0; i < 3; ++i) {
110+
for (size_t j = 0; j < 3; ++j) {
111+
Tensor<T, 3> gamma_ij_grid =
112+
tensorium_RG::populate_tensor3D_component<T>(
113+
i, j, func, dx, dy, dz, NX, NY, NZ);
114+
115+
Tensor<T, 4> grad_tmp;
116+
spectral_partial_scalar_3D(gamma_ij_grid, dx, dy, dz, grad_tmp);
117+
118+
for (size_t k = 0; k < 3; ++k)
119+
grad_gamma(i, j, k) = grad_tmp(NX/2, NY/2, NZ/2, k);
120+
}
121+
}
122+
123+
return grad_gamma;
124+
}
93125

94126
template<typename T, typename ScalarFunc>
95127
tensorium::Vector<T> partial_scalar(
96128
const tensorium::Vector<T>& X,
97129
T dx, T dy, T dz,
98130
ScalarFunc&& func)
99131
{
100-
tensorium::Vector<T> grad(3); // <-- ici 3, pas 4
132+
tensorium::Vector<T> grad(3);
101133

102-
// dérivée selon x¹
103134
{
104135
tensorium::Vector<T> Xs = X;
105136
T gm2, gm1, gp1, gp2;

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

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
#pragma once
2+
#include "../../Core/Tensor.hpp"
3+
#include "../../Core/Vector.hpp"
4+
5+
namespace tensorium_RG {
6+
7+
template<typename T, typename TensorFunc>
8+
tensorium::Tensor<T, 3> populate_tensor3D_component(
9+
size_t i, size_t j,
10+
TensorFunc&& func,
11+
T dx, T dy, T dz,
12+
size_t NX, size_t NY, size_t NZ)
13+
{
14+
using namespace tensorium;
15+
Tensor<T, 3> output({NX, NY, NZ});
16+
17+
for (size_t x = 0; x < NX; ++x) {
18+
for (size_t y = 0; y < NY; ++y) {
19+
for (size_t z = 0; z < NZ; ++z) {
20+
Vector<T> X(4);
21+
X(0) = 0.0;
22+
X(1) = x * dx;
23+
X(2) = y * dy;
24+
X(3) = z * dz;
25+
26+
Tensor<T, 2> gamma = func(X);
27+
output(x, y, z) = gamma(i, j);
28+
}
29+
}
30+
}
31+
return output;
32+
}
33+
34+
}

0 commit comments

Comments
 (0)