Skip to content

Commit 11c7e15

Browse files
at0m741at0m741
committed
THEY ROTATE
1 parent 1311ad7 commit 11c7e15

4 files changed

Lines changed: 151 additions & 130 deletions

File tree

‎.gitignore‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,9 @@ compile_commands.json
1212
*.so
1313
*.a
1414
*.csv
15+
*.png
16+
*.jpeg
17+
*.jpg
1518
# python
1619
.venv/
1720
.make

‎Tests/bssn/bowen_york_boost_stability.cpp‎

Lines changed: 15 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -54,8 +54,8 @@ void initialize_single_boost(Grid &grid) {
5454
const double m1 = 1.0;
5555
const double m2 = 1.0;
5656

57-
const double x1 = -3.0, y1 = 0.0, z1 = 0.0;
58-
const double x2 = 3.0, y2 = 0.0, z2 = 0.0;
57+
const double x1 = -4.0, y1 = 0.0, z1 = 0.0;
58+
const double x2 = 4.0, y2 = 0.0, z2 = 0.0;
5959

6060
const double Py = 0.095;
6161
const double P1[3] = {0.0, Py, 0.0};
@@ -73,10 +73,10 @@ void initialize_single_boost(Grid &grid) {
7373
REGISTER_TEST(
7474
"bssn.viz.moving_puncture", "Export CSV slices of a moving spinning black hole", []() {
7575
tensorium::tests::StabilityRunConfig cfg;
76-
cfg.nx = 64;
77-
cfg.ny = 64;
78-
cfg.nz = 64;
79-
cfg.spacing = 0.1;
76+
cfg.nx = 128;
77+
cfg.ny = 128;
78+
cfg.nz = 128;
79+
cfg.spacing = 0.4;
8080
cfg.ng = 4;
8181
cfg.padding = 4;
8282
cfg.steps = 600;
@@ -102,8 +102,8 @@ REGISTER_TEST(
102102
proj_cfg.recompute_inverse = true;
103103

104104
tensorium_RG::bssn::GaugeParameters<double> params;
105-
params.eta = 3.0;
106-
params.beta_B_coeff = 0.5;
105+
params.eta = 2.0;
106+
params.beta_B_coeff = 0.75;
107107

108108
tensorium_RG::bssn::project_bssn_state(grid, proj_cfg);
109109

@@ -134,10 +134,10 @@ REGISTER_TEST(
134134
REGISTER_TEST(
135135
"bssn.init.kerr_schild", "Print formatted BSSN fields after Kerr–Schild initialization", []() {
136136
tensorium::tests::StabilityRunConfig cfg;
137-
cfg.nx = 64;
138-
cfg.ny = 64;
139-
cfg.nz = 64;
140-
cfg.spacing = 0.1;
137+
cfg.nx = 96;
138+
cfg.ny = 96;
139+
cfg.nz = 96;
140+
cfg.spacing = 0.3;
141141
cfg.ng = 4;
142142
cfg.steps = 60;
143143
cfg.cfl = 0.15;
@@ -152,7 +152,7 @@ REGISTER_TEST(
152152
grid.y0 = -0.5 * cfg.spacing * cfg.ny + 0.5 * cfg.spacing;
153153
grid.z0 = -0.5 * cfg.spacing * cfg.nz + 0.5 * cfg.spacing;
154154

155-
tensorium_RG::init::kerr_schild_single(grid, 1.0, 0.0);
155+
tensorium_RG::init::kerr_schild_single(grid, 1.0, 0.9);
156156

157157
size_t I0, I1, J0, J1, K0, K1;
158158
grid.domain_bounds(I0, I1, J0, J1, K0, K1);
@@ -161,8 +161,8 @@ REGISTER_TEST(
161161
const size_t jc = J0 + grid.dims.ny / 2;
162162
const size_t kc = K0 + grid.dims.nz / 2;
163163
tensorium_RG::bssn::GaugeParameters<double> params;
164-
params.eta = 3.0;
165-
params.beta_B_coeff = 0.5;
164+
params.eta = 2.0;
165+
params.beta_B_coeff = 0.75;
166166

167167
tensorium_RG::init::print_bssn_state_at(grid, ic, jc, kc, "Kerr–Schild single BH");
168168
tensorium_RG::bssn::BSSNRKStepper<double, tensorium_RG::bssn::BoundaryRadiative> stepper(

‎includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/Evolution/BSSNEvolutionGammaTilde.hpp‎

Lines changed: 60 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#pragma once
22

3+
#include "../Derivatives/BSSNGridDerivatives.hpp"
34
#include "../Fields/BSSNGridSoA.hpp"
45
#include "Tensorium_Grid/Grid/GridLayout.hpp"
56

@@ -8,15 +9,18 @@
89
* @brief RHS for the conformal metric \f$\tilde{\gamma}_{ij}\f$.
910
* @details Evaluates
1011
* \f[
11-
* \partial_t\tilde{\gamma}_{ij} = -2\alpha\tilde{A}_{ij} + \mathcal{L}_\beta \tilde{\gamma}_{ij} - \tfrac{2}{3}\tilde{\gamma}_{ij}\partial_k\beta^k + \mathcal{D}_6[\tilde{\gamma}_{ij}].
12+
* \partial_t\tilde{\gamma}_{ij} = -2\alpha\tilde{A}_{ij} + \mathcal{L}_\beta \tilde{\gamma}_{ij} -
13+
* \tfrac{2}{3}\tilde{\gamma}_{ij}\partial_k\beta^k + \mathcal{D}_6[\tilde{\gamma}_{ij}].
1214
* \f]
13-
* Lie derivatives expand into \f$\tilde{\gamma}_{ik}\partial_j\beta^k + \tilde{\gamma}_{jk}\partial_i\beta^k\f$ plus the advection term \f$\beta^k\partial_k\tilde{\gamma}_{ij}\f$.
15+
* Lie derivatives expand into \f$\tilde{\gamma}_{ik}\partial_j\beta^k +
16+
* \tilde{\gamma}_{jk}\partial_i\beta^k\f$ plus the advection term
17+
* \f$\beta^k\partial_k\tilde{\gamma}_{ij}\f$.
1418
*/
1519

1620
namespace tensorium_RG::bssn {
1721

1822
#ifndef TENSORIUM_BSSN_EVOLUTION_CLAMP_HELPERS_DEFINED
19-
#define TENSORIUM_BSSN_EVOLUTION_CLAMP_HELPERS_DEFINED
23+
# define TENSORIUM_BSSN_EVOLUTION_CLAMP_HELPERS_DEFINED
2024
inline size_t clamped_lower(size_t lower, size_t guard, size_t upper) {
2125
return std::min(lower + guard, upper);
2226
}
@@ -32,11 +36,12 @@ inline size_t clamped_upper(size_t upper, size_t guard, size_t lower) {
3236
template <typename T>
3337
inline void compute_rhs_gamma_tilde(const BSSNGridSoA<T> &G, Field3D<T> rhs[6],
3438
size_t padding = 4) {
39+
using namespace tensorium_RG::fd;
3540
size_t I0, I1, J0, J1, K0, K1;
3641
G.domain_bounds(I0, I1, J0, J1, K0, K1);
3742

38-
const size_t total = G.gamma_tilde[0].st.nx_tot * G.gamma_tilde[0].st.ny_tot *
39-
G.gamma_tilde[0].st.nz_tot;
43+
const size_t total =
44+
G.gamma_tilde[0].st.nx_tot * G.gamma_tilde[0].st.ny_tot * G.gamma_tilde[0].st.nz_tot;
4045
for (int s = 0; s < 6; ++s)
4146
std::fill(rhs[s].ptr(), rhs[s].ptr() + total, T(0));
4247

@@ -50,58 +55,65 @@ inline void compute_rhs_gamma_tilde(const BSSNGridSoA<T> &G, Field3D<T> rhs[6],
5055
if (i0 >= i1 || j0 >= j1 || k0 >= k1)
5156
return;
5257

53-
const T ko_sigma = T(5.0);
58+
const T ko_sigma = T(0.6);
5459

5560
#pragma omp parallel for collapse(3)
56-
for (size_t i = i0; i < i1; ++i)
57-
for (size_t j = j0; j < j1; ++j)
61+
for (size_t i = i0; i < i1; ++i) {
62+
for (size_t j = j0; j < j1; ++j) {
5863
for (size_t k = k0; k < k1; ++k) {
5964
const size_t id = G.gamma_tilde[0].idx(i, j, k);
60-
const T alpha = G.alpha.ptr()[id];
61-
62-
const T div_beta = tensorium_RG::fd::Dx(G.beta[0], i, j, k, G.dx) +
63-
tensorium_RG::fd::Dy(G.beta[1], i, j, k, G.dy) +
64-
tensorium_RG::fd::Dz(G.beta[2], i, j, k, G.dz);
65-
66-
auto d_beta = [&](int idx, size_t dir_i, size_t dir_j, size_t dir_k) {
67-
if (idx == 0)
68-
return tensorium_RG::fd::Dx(G.beta[dir_i], i, j, k, dir_j == 0 ? G.dx : G.dy);
69-
if (idx == 1)
70-
return tensorium_RG::fd::Dy(G.beta[dir_i], i, j, k, dir_j == 1 ? G.dy : G.dz);
71-
return tensorium_RG::fd::Dz(G.beta[dir_i], i, j, k, G.dz);
72-
};
73-
74-
const T beta_x = G.beta[0].ptr()[id];
75-
const T beta_y = G.beta[1].ptr()[id];
76-
const T beta_z = G.beta[2].ptr()[id];
65+
const T alpha = G.alpha.ptr()[id];
66+
67+
T d_beta[3][3];
68+
d_beta[0][0] = Dx(G.beta[0], i, j, k, G.dx);
69+
d_beta[0][1] = Dy(G.beta[0], i, j, k, G.dy);
70+
d_beta[0][2] = Dz(G.beta[0], i, j, k, G.dz);
71+
d_beta[1][0] = Dx(G.beta[1], i, j, k, G.dx);
72+
d_beta[1][1] = Dy(G.beta[1], i, j, k, G.dy);
73+
d_beta[1][2] = Dz(G.beta[1], i, j, k, G.dz);
74+
d_beta[2][0] = Dx(G.beta[2], i, j, k, G.dx);
75+
d_beta[2][1] = Dy(G.beta[2], i, j, k, G.dy);
76+
d_beta[2][2] = Dz(G.beta[2], i, j, k, G.dz);
77+
78+
const T div_beta = d_beta[0][0] + d_beta[1][1] + d_beta[2][2];
79+
const T beta_vec[3] = {G.beta[0].ptr()[id], G.beta[1].ptr()[id],
80+
G.beta[2].ptr()[id]};
81+
82+
T gam[3][3];
83+
for (int u = 0; u < 3; ++u)
84+
for (int v = u; v < 3; ++v) {
85+
T val = G.gamma_tilde[tensorium_RG::sym6_index(u, v)].ptr()[id];
86+
gam[u][v] = val;
87+
gam[v][u] = val;
88+
}
7789

78-
for (int a = 0; a < 3; ++a)
90+
for (int a = 0; a < 3; ++a) {
7991
for (int b = a; b < 3; ++b) {
8092
const int s = tensorium_RG::sym6_index(a, b);
81-
const T gamma = G.gamma_tilde[s].ptr()[id];
82-
T adv = beta_x * tensorium_RG::fd::Dx_upwind(G.gamma_tilde[s], i, j, k, G.dx, beta_x) +
83-
beta_y * tensorium_RG::fd::Dy_upwind(G.gamma_tilde[s], i, j, k, G.dy, beta_y) +
84-
beta_z * tensorium_RG::fd::Dz_upwind(G.gamma_tilde[s], i, j, k, G.dz, beta_z);
93+
94+
T adv = beta_vec[0] * tensorium_RG::fd::Dx_upwind(G.gamma_tilde[s], i, j, k,
95+
G.dx, beta_vec[0]) +
96+
beta_vec[1] * tensorium_RG::fd::Dy_upwind(G.gamma_tilde[s], i, j, k,
97+
G.dy, beta_vec[1]) +
98+
beta_vec[2] * tensorium_RG::fd::Dz_upwind(G.gamma_tilde[s], i, j, k,
99+
G.dz, beta_vec[2]);
100+
85101
T lie = T(0);
86-
lie += G.gamma_tilde[tensorium_RG::sym6_index(a, 0)].ptr()[id] *
87-
tensorium_RG::fd::Dx(G.beta[b], i, j, k, G.dx);
88-
lie += G.gamma_tilde[tensorium_RG::sym6_index(b, 0)].ptr()[id] *
89-
tensorium_RG::fd::Dx(G.beta[a], i, j, k, G.dx);
90-
lie += G.gamma_tilde[tensorium_RG::sym6_index(a, 1)].ptr()[id] *
91-
tensorium_RG::fd::Dy(G.beta[b], i, j, k, G.dy);
92-
lie += G.gamma_tilde[tensorium_RG::sym6_index(b, 1)].ptr()[id] *
93-
tensorium_RG::fd::Dy(G.beta[a], i, j, k, G.dy);
94-
lie += G.gamma_tilde[tensorium_RG::sym6_index(a, 2)].ptr()[id] *
95-
tensorium_RG::fd::Dz(G.beta[b], i, j, k, G.dz);
96-
lie += G.gamma_tilde[tensorium_RG::sym6_index(b, 2)].ptr()[id] *
97-
tensorium_RG::fd::Dz(G.beta[a], i, j, k, G.dz);
98-
99-
const T result = adv - T(2) * alpha * G.A_tilde[s].ptr()[id] + lie -
100-
T(2.0 / 3.0) * gamma * div_beta +
101-
T(tensorium_RG::fd::KO6(G.gamma_tilde[s], i, j, k, ko_sigma));
102-
rhs[s].ptr()[id] = result;
102+
for (int m = 0; m < 3; ++m) {
103+
lie += gam[a][m] * d_beta[m][b];
104+
lie += gam[b][m] * d_beta[m][a];
105+
}
106+
107+
const T source = -T(2) * alpha * G.A_tilde[s].ptr()[id];
108+
const T trace_rem = -T(2.0 / 3.0) * gam[a][b] * div_beta;
109+
const T diss =
110+
T(tensorium_RG::fd::KO6(G.gamma_tilde[s], i, j, k, ko_sigma));
111+
112+
rhs[s].ptr()[id] = adv + lie + source + trace_rem + diss;
103113
}
114+
}
104115
}
116+
}
117+
}
105118
}
106-
107119
} // namespace tensorium_RG::bssn

0 commit comments

Comments
 (0)