Skip to content

Commit 1b9277f

Browse files
at0m741at0m741
committed
spectral 3d
1 parent c542ebf commit 1b9277f

File tree

5 files changed

+322
-2
lines changed

5 files changed

+322
-2
lines changed

‎Tests/Derivatives/Derivatives.cpp‎

Lines changed: 138 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,142 @@
11
#include "../test.hpp"
2+
#include <fstream>
3+
int deriv_test_spectral_fft() {
4+
using T = float;
5+
using C = std::complex<T>;
6+
using Tensor3D = tensorium::Tensor<C, 3>;
7+
8+
const size_t N = 128;
9+
const T L = 1.0;
10+
const T dx = L / N;
11+
const T pi = static_cast<T>(3.14159265358979323846);
12+
const T TWO_PI = 2 * pi;
13+
14+
Tensor3D f({N, N, N});
15+
Tensor3D df_dx_exact({N, N, N});
16+
Tensor3D df_dx_numeric({N, N, N});
17+
18+
for (size_t i = 0; i < N; ++i) {
19+
for (size_t j = 0; j < N; ++j) {
20+
for (size_t k = 0; k < N; ++k) {
21+
T x = i * dx;
22+
T y = j * dx;
23+
T z = k * dx;
24+
f(i, j, k) = std::sin(TWO_PI * x) * std::sin(TWO_PI * y) * std::sin(TWO_PI * z);
25+
df_dx_exact(i, j, k) = TWO_PI * std::cos(TWO_PI * x) * std::sin(TWO_PI * y) * std::sin(TWO_PI * z);
26+
}
27+
}
28+
}
29+
30+
Tensor3D F = f;
31+
tensorium::SpectralFFT<T>::forward_3D(F);
32+
33+
for (size_t i = 0; i < N; ++i) {
34+
int ki = (i <= N / 2) ? i : i - N;
35+
T kx = TWO_PI * ki / L;
36+
for (size_t j = 0; j < N; ++j) {
37+
for (size_t k = 0; k < N; ++k) {
38+
F(i,j,k) *= C(0, kx);
39+
}
40+
}
41+
}
42+
43+
tensorium::SpectralFFT<T>::backward_3D(F);
44+
df_dx_numeric = F;
45+
46+
T max_err = 0;
47+
for (size_t i = 0; i < N; ++i)
48+
for (size_t j = 0; j < N; ++j)
49+
for (size_t k = 0; k < N; ++k)
50+
max_err = std::max(max_err, std::abs(df_dx_exact(i,j,k).real() - df_dx_numeric(i,j,k).real()));
51+
52+
std::cout << "\n=== Spectral Derivative Test (3D FFT) ===\n";
53+
std::cout << "Erreur max sur d/dx (FFT vs exact): " << max_err << "\n";
54+
55+
std::cout << "\nValeurs (x=0 to 1):\n";
56+
for (size_t i = 0; i < 4; ++i) {
57+
for (size_t j = 0; j < 1; ++j) {
58+
for (size_t k = 0; k < 1; ++k) {
59+
std::cout << "x=" << i * dx << "\tf=" << f(i,j,k).real()
60+
<< "\texact=" << df_dx_exact(i,j,k).real()
61+
<< "\tnumeric=" << df_dx_numeric(i,j,k).real() << "\n";
62+
T err = std::abs(df_dx_numeric(i,j,k).real() - df_dx_exact(i,j,k).real());
63+
std::cout << "\terreur=" << err << "\n";
64+
65+
}
66+
}
67+
}
68+
69+
70+
std::cout << "\n=== Fonction test cos*sin*cos ===\n";
71+
72+
const T A = 2 * pi; // freq_x
73+
const T B = 4 * pi; // freq_y
74+
const T Ce = 6 * pi; // freq_z
75+
76+
for (size_t i = 0; i < N; ++i) {
77+
for (size_t j = 0; j < N; ++j) {
78+
for (size_t k = 0; k < N; ++k) {
79+
T x = i * dx;
80+
T y = j * dx;
81+
T z = k * dx;
82+
f(i, j, k) = std::sin(A * x) * std::cos(B * y) * std::sin(Ce * z);
83+
df_dx_exact(i, j, k) = A * std::cos(A * x) * std::cos(B * y) * std::sin(Ce * z);
84+
}
85+
}
86+
}
87+
88+
F = f;
89+
tensorium::SpectralFFT<T>::forward_3D(F);
90+
91+
for (size_t i = 0; i < N; ++i) {
92+
int ki = (i <= N / 2) ? i : i - N;
93+
T kx = TWO_PI * ki / L;
94+
for (size_t j = 0; j < N; ++j) {
95+
for (size_t k = 0; k < N; ++k) {
96+
F(i,j,k) *= C(0, kx);
97+
}
98+
}
99+
}
100+
101+
tensorium::SpectralFFT<T>::backward_3D(F);
102+
df_dx_numeric = F;
103+
104+
max_err = 0;
105+
for (size_t i = 0; i < N; ++i)
106+
for (size_t j = 0; j < N; ++j)
107+
for (size_t k = 0; k < N; ++k)
108+
max_err = std::max(max_err, std::abs(df_dx_exact(i,j,k).real() - df_dx_numeric(i,j,k).real()));
109+
110+
std::cout << "Erreur max (cos*sin*cos): " << max_err << "\n";
111+
for (size_t i = 0; i < 4; ++i) {
112+
size_t j = 1, k = 1;
113+
T x = i * dx;
114+
std::cout << "x=" << x
115+
<< "\tf=" << f(i,j,k).real()
116+
<< "\texact=" << df_dx_exact(i,j,k).real()
117+
<< "\tnumeric=" << df_dx_numeric(i,j,k).real()
118+
<< "\terreur=" << std::abs(df_dx_exact(i,j,k).real() - df_dx_numeric(i,j,k).real()) << "\n";
119+
}
120+
121+
std::ofstream file("spectral_deriv_output.csv");
122+
file << "x,y,z,f,dfdx_exact,dfdx_numeric\n";
123+
for (size_t i = 0; i < N; ++i) {
124+
for (size_t j = 0; j < N; ++j) {
125+
for (size_t k = 0; k < N; ++k) {
126+
T x = i * dx;
127+
T y = j * dx;
128+
T z = k * dx;
129+
file << x << "," << y << "," << z << ","
130+
<< f(i,j,k).real() << ","
131+
<< df_dx_exact(i,j,k).real() << ","
132+
<< df_dx_numeric(i,j,k).real() << "\n";
133+
}
134+
}
135+
}
136+
file.close();
137+
138+
return 0;
139+
}
2140

3141
int deriv_test() {
4142
std::cout << "\n=== Derivate 2D Test (\u2202/\u2202x) ===\n";

‎Tests/main.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,6 @@ int main() {
2222
matrix_tests();
2323
tensor_test();
2424
vector_tests();
25-
25+
deriv_test_spectral_fft();
2626
return 0;
2727
}

‎Tests/test.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,3 +19,4 @@ int linear_solver_test();
1919
int matrix_tests();
2020
int tensor_test();
2121
int vector_tests();
22+
int deriv_test_spectral_fft();

‎includes/Tensorium/Core/Spectral.hpp‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,13 +37,110 @@ namespace tensorium {
3737
* @param a Input/output complex vector (must have power-of-two size)
3838
*/
3939
static inline void forward(CVectorT& a) { transform_impl(a, false); }
40+
41+
static void forward_3D(Tensor<std::complex<T>, 3>& a) {
42+
const auto shape = a.shape();
43+
const size_t NX = shape[0], NY = shape[1], NZ = shape[2];
44+
using CVector = Vector<std::complex<T>>;
45+
46+
#pragma omp parallel for collapse(2)
47+
for (size_t i = 0; i < NX; ++i) {
48+
for (size_t j = 0; j < NY; ++j) {
49+
CVector sliceZ(NZ);
50+
for (size_t k = 0; k < NZ; ++k)
51+
sliceZ(k) = a(i,j,k);
52+
53+
forward(sliceZ);
54+
55+
for (size_t k = 0; k < NZ; ++k)
56+
a(i,j,k) = sliceZ(k);
57+
}
58+
}
59+
60+
#pragma omp parallel for collapse(2)
61+
for (size_t i = 0; i < NX; ++i) {
62+
for (size_t k = 0; k < NZ; ++k) {
63+
CVector sliceY(NY);
64+
for (size_t j = 0; j < NY; ++j)
65+
sliceY(j) = a(i,j,k);
66+
67+
forward(sliceY);
68+
69+
for (size_t j = 0; j < NY; ++j)
70+
a(i,j,k) = sliceY(j);
71+
}
72+
}
73+
74+
#pragma omp parallel for collapse(2)
75+
for (size_t j = 0; j < NY; ++j) {
76+
for (size_t k = 0; k < NZ; ++k) {
77+
CVector sliceX(NX);
78+
for (size_t i = 0; i < NX; ++i)
79+
sliceX(i) = a(i,j,k);
80+
81+
forward(sliceX);
82+
83+
for (size_t i = 0; i < NX; ++i)
84+
a(i,j,k) = sliceX(i);
85+
}
86+
}
87+
}
88+
4089
/**
4190
* @brief Perform inverse FFT (in-place)
4291
*
4392
* @param a Input/output complex vector (must have power-of-two size)
4493
*/
4594
static inline void backward(CVectorT& a) { transform_impl(a, true ); }
4695

96+
static void backward_3D(Tensor<std::complex<T>, 3>& a) {
97+
const auto shape = a.shape();
98+
const size_t NX = shape[0], NY = shape[1], NZ = shape[2];
99+
using CVector = Vector<std::complex<T>>;
100+
101+
#pragma omp parallel for collapse(2)
102+
for (size_t j = 0; j < NY; ++j) {
103+
for (size_t k = 0; k < NZ; ++k) {
104+
CVector sliceX(NX);
105+
for (size_t i = 0; i < NX; ++i)
106+
sliceX(i) = a(i,j,k);
107+
108+
backward(sliceX);
109+
110+
for (size_t i = 0; i < NX; ++i)
111+
a(i,j,k) = sliceX(i);
112+
}
113+
}
114+
115+
#pragma omp parallel for collapse(2)
116+
for (size_t i = 0; i < NX; ++i) {
117+
for (size_t k = 0; k < NZ; ++k) {
118+
CVector sliceY(NY);
119+
for (size_t j = 0; j < NY; ++j)
120+
sliceY(j) = a(i,j,k);
121+
122+
backward(sliceY);
123+
124+
for (size_t j = 0; j < NY; ++j)
125+
a(i,j,k) = sliceY(j);
126+
}
127+
}
128+
129+
#pragma omp parallel for collapse(2)
130+
for (size_t i = 0; i < NX; ++i) {
131+
for (size_t j = 0; j < NY; ++j) {
132+
CVector sliceZ(NZ);
133+
for (size_t k = 0; k < NZ; ++k)
134+
sliceZ(k) = a(i,j,k);
135+
136+
backward(sliceZ);
137+
138+
for (size_t k = 0; k < NZ; ++k)
139+
a(i,j,k) = sliceZ(k);
140+
}
141+
}
142+
}
143+
47144
private:
48145
/**
49146
* @brief Internal FFT implementation (shared by forward/backward)

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

Lines changed: 85 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,95 @@
11
#pragma once
22

33
#include "../Metric.hpp"
4-
4+
#include "../../Core/Spectral.hpp"
55

66
namespace tensorium_RG {
77

88

9+
template<typename T>
10+
void spectral_derivative_1D(tensorium::Vector<T>& field, tensorium::Vector<T>& dfield, T dx) {
11+
using C = std::complex<T>;
12+
using FFT = tensorium::SpectralFFT<T>;
13+
14+
tensorium::Vector<C> field_fft(field.size());
15+
for (size_t i = 0; i < field.size(); ++i)
16+
field_fft(i) = field(i);
17+
18+
FFT::forward(field_fft);
19+
20+
const T L = dx * field.size();
21+
for (size_t k = 0; k < field.size(); ++k) {
22+
int n = (k <= field.size()/2) ? k : (int)k - (int)field.size();
23+
C factor = C(0, 2.0 * M_PI * n / L);
24+
field_fft(k) *= factor;
25+
}
26+
27+
FFT::backward(field_fft);
28+
for (size_t i = 0; i < field.size(); ++i)
29+
dfield(i) = field_fft(i).real();
30+
}
31+
32+
template<typename T>
33+
void spectral_partial_scalar_3D(
34+
const tensorium::Tensor<T, 3>& scalar_field,
35+
T dx, T dy, T dz,
36+
tensorium::Tensor<T, 4>& grad_out) {
37+
38+
using namespace tensorium;
39+
using C = std::complex<T>;
40+
using FFT = SpectralFFT<T>;
41+
42+
const auto& shape = scalar_field.shape();
43+
const size_t NX = shape[0], NY = shape[1], NZ = shape[2];
44+
const size_t N = NX * NY * NZ;
45+
46+
grad_out.resize(NX, NY, NZ, 3);
47+
48+
Tensor<C, 3> field_cplx({NX, NY, NZ});
49+
for (size_t i = 0; i < NX; ++i)
50+
for (size_t j = 0; j < NY; ++j)
51+
for (size_t k = 0; k < NZ; ++k)
52+
field_cplx(i,j,k) = C(scalar_field(i,j,k), 0.0);
53+
54+
FFT::forward_3D(field_cplx);
55+
56+
for (size_t i = 0; i < NX; ++i) {
57+
int ni = (i <= NX/2) ? i : int(i) - int(NX);
58+
T kx = 2.0 * M_PI * ni / (NX * dx);
59+
60+
for (size_t j = 0; j < NY; ++j) {
61+
int nj = (j <= NY/2) ? j : int(j) - int(NY);
62+
T ky = 2.0 * M_PI * nj / (NY * dy);
63+
64+
for (size_t k = 0; k < NZ; ++k) {
65+
int nk = (k <= NZ/2) ? k : int(k) - int(NZ);
66+
T kz = 2.0 * M_PI * nk / (NZ * dz);
67+
68+
C f_hat = field_cplx(i,j,k);
69+
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();
73+
}
74+
}
75+
}
76+
77+
for (int dim = 0; dim < 3; ++dim) {
78+
Tensor<C, 3> dfield_cplx({NX, NY, NZ});
79+
for (size_t i = 0; i < NX; ++i)
80+
for (size_t j = 0; j < NY; ++j)
81+
for (size_t k = 0; k < NZ; ++k)
82+
dfield_cplx(i,j,k) = C(grad_out(i,j,k,dim), 0.0);
83+
84+
FFT::backward_3D(dfield_cplx);
85+
86+
for (size_t i = 0; i < NX; ++i)
87+
for (size_t j = 0; j < NY; ++j)
88+
for (size_t k = 0; k < NZ; ++k)
89+
grad_out(i,j,k,dim) = dfield_cplx(i,j,k).real();
90+
}
91+
}
92+
993

1094
template<typename T, typename ScalarFunc>
1195
tensorium::Vector<T> partial_scalar(

0 commit comments

Comments
 (0)