Skip to content

Commit 256fee5

Browse files
at0m741at0m741
committed
derivatives + Christoffel init (need to be computed inside kernels and not precomputed in grid init (-27 memory read at each grid cell))
1 parent 9db805f commit 256fee5

6 files changed

Lines changed: 462 additions & 50 deletions

File tree

‎includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/BSSNCHristoffelTilde.hpp‎

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
#include "BSSNGridDerivatives.hpp"
2+
#include "BSSNGridSoA.hpp"
3+
4+
namespace tensorium_RG::bssn {
5+
6+
template <typename T> inline void compute_tildeGamma_contracted(BSSNGridSoA<T> &G) {
7+
using tensorium_RG::fd::Dx;
8+
using tensorium_RG::fd::Dy;
9+
using tensorium_RG::fd::Dz;
10+
11+
size_t I0, I1, J0, J1, K0, K1;
12+
G.domain_bounds(I0, I1, J0, J1, K0, K1);
13+
14+
const size_t i0 = I0 + 2, i1 = I1 - 2;
15+
const size_t j0 = J0 + 2, j1 = J1 - 2;
16+
const size_t k0 = K0 + 2, k1 = K1 - 2;
17+
18+
#pragma omp parallel for collapse(2)
19+
for (size_t i = i0; i < i1; ++i)
20+
for (size_t j = j0; j < j1; ++j)
21+
for (size_t k = k0; k < k1; ++k) {
22+
23+
const size_t id = G.gamma_tilde[XX].idx(i, j, k);
24+
25+
auto d_g = [&](int m, int a, int b) -> T {
26+
if (a > b) {
27+
int t = a;
28+
a = b;
29+
b = t;
30+
}
31+
const Field3D<T> &F = (a == 0 && b == 0) ? G.gamma_tilde[XX]
32+
: (a == 0 && b == 1) ? G.gamma_tilde[XY]
33+
: (a == 0 && b == 2) ? G.gamma_tilde[XZ]
34+
: (a == 1 && b == 1) ? G.gamma_tilde[YY]
35+
: (a == 1 && b == 2) ? G.gamma_tilde[YZ]
36+
: G.gamma_tilde[ZZ];
37+
38+
if (m == 0)
39+
return Dx(F, i, j, k, G.dx);
40+
if (m == 1)
41+
return Dy(F, i, j, k, G.dy);
42+
return Dz(F, i, j, k, G.dz);
43+
};
44+
45+
T GammaContr[3] = {T(0), T(0), T(0)};
46+
47+
// Γ̃^i = g̃^{jk} Γ̃^i_{jk}
48+
for (int ii = 0; ii < 3; ++ii) {
49+
T acc = T(0);
50+
51+
for (int jj = 0; jj < 3; ++jj)
52+
for (int kk = 0; kk < 3; ++kk) {
53+
T sum_l = T(0);
54+
for (int ll = 0; ll < 3; ++ll) {
55+
const T ginv_il = sym6_inv_get(G.gamma_tilde_inv, id, ii, ll);
56+
57+
const T term = d_g(jj, kk, ll) + d_g(kk, jj, ll) - d_g(ll, jj, kk);
58+
59+
sum_l += ginv_il * term;
60+
}
61+
const T Gamma_i_jk = T(0.5) * sum_l;
62+
63+
const T ginv_jk = sym6_inv_get(G.gamma_tilde_inv, id, jj, kk);
64+
acc += ginv_jk * Gamma_i_jk;
65+
}
66+
67+
GammaContr[ii] = acc;
68+
}
69+
70+
G.tildeGamma[0].ptr()[id] = GammaContr[0];
71+
G.tildeGamma[1].ptr()[id] = GammaContr[1];
72+
G.tildeGamma[2].ptr()[id] = GammaContr[2];
73+
}
74+
}
75+
76+
template <typename T> inline void compute_tildeGamma_full(BSSNGridSoA<T> &G, Field3D<T> *Gamma27) {
77+
using tensorium_RG::fd::Dx;
78+
using tensorium_RG::fd::Dy;
79+
using tensorium_RG::fd::Dz;
80+
81+
size_t I0, I1, J0, J1, K0, K1;
82+
G.domain_bounds(I0, I1, J0, J1, K0, K1);
83+
84+
const size_t i0 = I0 + 2, i1 = I1 - 2;
85+
const size_t j0 = J0 + 2, j1 = J1 - 2;
86+
const size_t k0 = K0 + 2, k1 = K1 - 2;
87+
88+
#pragma omp parallel for collapse(2)
89+
for (size_t i = i0; i < i1; ++i)
90+
for (size_t j = j0; j < j1; ++j)
91+
for (size_t k = k0; k < k1; ++k) {
92+
93+
const size_t id = G.gamma_tilde[XX].idx(i, j, k);
94+
95+
auto d_g = [&](int m, int a, int b) -> T {
96+
if (a > b) {
97+
int t = a;
98+
a = b;
99+
b = t;
100+
}
101+
const Field3D<T> &F = (a == 0 && b == 0) ? G.gamma_tilde[XX]
102+
: (a == 0 && b == 1) ? G.gamma_tilde[XY]
103+
: (a == 0 && b == 2) ? G.gamma_tilde[XZ]
104+
: (a == 1 && b == 1) ? G.gamma_tilde[YY]
105+
: (a == 1 && b == 2) ? G.gamma_tilde[YZ]
106+
: G.gamma_tilde[ZZ];
107+
108+
if (m == 0)
109+
return Dx(F, i, j, k, G.dx);
110+
if (m == 1)
111+
return Dy(F, i, j, k, G.dy);
112+
return Dz(F, i, j, k, G.dz);
113+
};
114+
115+
for (int kk = 0; kk < 3; ++kk)
116+
for (int ii = 0; ii < 3; ++ii)
117+
for (int jj = 0; jj < 3; ++jj) {
118+
119+
T sum = T(0);
120+
for (int ll = 0; ll < 3; ++ll) {
121+
const T ginv_kl = sym6_inv_get(G.gamma_tilde_inv, id, kk, ll);
122+
const T term = d_g(ii, jj, ll) + d_g(jj, ii, ll) - d_g(ll, ii, jj);
123+
sum += ginv_kl * term;
124+
}
125+
const T Gkij = T(0.5) * sum;
126+
Gamma27[kk * 9 + ii * 3 + jj].ptr()[id] = Gkij;
127+
}
128+
}
129+
}
130+
131+
} // namespace tensorium_RG::bssn

‎includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/BSSNGridDerivatives.hpp‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
2+
#pragma once
3+
4+
#include <cstddef>
5+
6+
namespace tensorium_RG::fd {
7+
8+
enum Axis { X = 0, Y = 1, Z = 2 };
9+
10+
template <typename Field> inline double get(const Field &f, size_t i, size_t j, size_t k) {
11+
return f.ptr()[f.idx(i, j, k)];
12+
}
13+
14+
template <typename Field>
15+
inline double Dx(const Field &f, size_t i, size_t j, size_t k, double dx) {
16+
return (-get(f, i + 2, j, k) + 8.0 * get(f, i + 1, j, k) - 8.0 * get(f, i - 1, j, k) +
17+
get(f, i - 2, j, k)) /
18+
(12.0 * dx);
19+
}
20+
21+
template <typename Field>
22+
inline double Dy(const Field &f, size_t i, size_t j, size_t k, double dy) {
23+
return (-get(f, i, j + 2, k) + 8.0 * get(f, i, j + 1, k) - 8.0 * get(f, i, j - 1, k) +
24+
get(f, i, j - 2, k)) /
25+
(12.0 * dy);
26+
}
27+
28+
template <typename Field>
29+
inline double Dz(const Field &f, size_t i, size_t j, size_t k, double dz) {
30+
return (-get(f, i, j, k + 2) + 8.0 * get(f, i, j, k + 1) - 8.0 * get(f, i, j, k - 1) +
31+
get(f, i, j, k - 2)) /
32+
(12.0 * dz);
33+
}
34+
35+
template <typename Field>
36+
inline double D(const Field &f, size_t i, size_t j, size_t k, Axis a, double dx, double dy,
37+
double dz) {
38+
if (a == X)
39+
return Dx(f, i, j, k, dx);
40+
if (a == Y)
41+
return Dy(f, i, j, k, dy);
42+
return Dz(f, i, j, k, dz);
43+
}
44+
45+
} // namespace tensorium_RG::fd

‎includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/BSSNGridSoA.hpp‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,8 @@ template <typename T> class BSSNGridSoA {
2121
Field3D<T> gamma_tilde_inv[6];
2222
Field3D<T> A_tilde[6];
2323

24-
T dx, dy, dz;
24+
Field3D<T> Gamma_tilde[27];
25+
T dx, dy, dz;
2526

2627
BSSNGridSoA(size_t nx, size_t ny, size_t nz, size_t ng, T dx_, T dy_, T dz_)
2728
: dims{nx, ny, nz, ng},
@@ -61,6 +62,8 @@ template <typename T> class BSSNGridSoA {
6162
alloc_field(gamma_tilde_inv[s]);
6263
alloc_field(A_tilde[s]);
6364
}
65+
for (int q = 0; q < 27; ++q)
66+
alloc_field(Gamma_tilde[q]);
6467
}
6568

6669
inline void domain_bounds(size_t &i0, size_t &i1, size_t &j0, size_t &j1, size_t &k0,
@@ -99,6 +102,8 @@ template <typename Boundary, typename T> inline void apply_halos_grid(BSSNGridSo
99102
Boundary::apply(G.gamma_tilde_inv[s], D);
100103
Boundary::apply(G.A_tilde[s], D);
101104
}
105+
for (int q = 0; q < 27; ++q)
106+
Boundary::apply(G.Gamma_tilde[q], D);
102107
}
103108

104109
} // namespace tensorium_RG

0 commit comments

Comments
 (0)