Covariant derivative interface for kerr initial data (extrinsic curvature tensor construction)

Author: at0m741

Tests/CMakeLists.txt

@@ -27,8 +27,16 @@ set(TEST_BSSN_EVOLUTION_SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/bssn/evolution/ATildeEvolutionTests.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/bssn/evolution/KEvolutionTests.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/bssn/evolution/GammaEvolutionTests.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/evolution/GaugeEvolutionTests.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/evolution/RK4TimeEvolutionTests.cpp)
+    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/evolution/GaugeEvolutionTests.cpp)
+
+set(TEST_BSSN_STABILITY_SOURCES
+    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/minkowski_stability.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/schwarzschild_stability.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/kerr_basic.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/bowen_york_boost_stability.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/brown_spherical_analytic_constraints.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/brown_hyperbolicity_speeds.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/bssn/boundary_conditions.cpp)
 
 set(BENCHMARK_SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/core/MatrixBenchmark.cpp)
@@ -47,13 +55,15 @@ tensorium_add_object_library(test_bssn_constraints_objs ${TEST_BSSN_CONSTRAINT_S
 tensorium_add_object_library(test_bssn_ricci_objs ${TEST_BSSN_RICCI_SOURCES})
 tensorium_add_object_library(test_bssn_solvers_objs ${TEST_BSSN_SOLVER_SOURCES})
 tensorium_add_object_library(test_bssn_evolution_objs ${TEST_BSSN_EVOLUTION_SOURCES})
+tensorium_add_object_library(test_bssn_stability_objs ${TEST_BSSN_STABILITY_SOURCES})
 
 set(BSSN_OBJECTS
     $<TARGET_OBJECTS:test_bssn_initial_objs>
     $<TARGET_OBJECTS:test_bssn_constraints_objs>
     $<TARGET_OBJECTS:test_bssn_ricci_objs>
     $<TARGET_OBJECTS:test_bssn_solvers_objs>
-    $<TARGET_OBJECTS:test_bssn_evolution_objs>)
+    $<TARGET_OBJECTS:test_bssn_evolution_objs>
+    $<TARGET_OBJECTS:test_bssn_stability_objs>)
 
 add_executable(TensoriumChiEvolutionDemo
     bssn/binaries/BSSNChiEvolutionDemo.cpp)

Tests/bssn/boundary_conditions.cpp

@@ -0,0 +1,169 @@
+#include "../framework/Assertions.hpp"
+#include "../framework/TestRegistry.hpp"
+
+#include "../../includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/Geometry/BSSNInvariants.hpp"
+#include "../../includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/Grid/BSSNGridOperations.hpp"
+#include "../../includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/InitialData/BSSNInitialData.hpp"
+
+#include <algorithm>
+#include <cmath>
+
+namespace {
+
+using Grid = tensorium_RG::BSSNGridSoA<double>;
+
+void center_grid(Grid &grid) {
+    const double hx = 0.5 * (static_cast<double>(grid.dims.nx) - 1.0) * grid.dx;
+    const double hy = 0.5 * (static_cast<double>(grid.dims.ny) - 1.0) * grid.dy;
+    const double hz = 0.5 * (static_cast<double>(grid.dims.nz) - 1.0) * grid.dz;
+    grid.x0 = -hx;
+    grid.y0 = -hy;
+    grid.z0 = -hz;
+}
+
+void seed_chi_shell(Grid &grid, double amplitude, size_t shell) {
+    size_t i0, i1, j0, j1, k0, k1;
+    grid.domain_bounds(i0, i1, j0, j1, k0, k1);
+    if (shell == 0)
+        return;
+
+    for (size_t i = i0; i < i1; ++i) {
+        const size_t di = std::min(i - i0, i1 - 1 - i);
+        for (size_t j = j0; j < j1; ++j) {
+            const size_t dj = std::min(j - j0, j1 - 1 - j);
+            for (size_t k = k0; k < k1; ++k) {
+                const size_t dk = std::min(k - k0, k1 - 1 - k);
+                const size_t dmin = std::min({di, dj, dk});
+                if (dmin >= shell)
+                    continue;
+                const double weight = static_cast<double>(shell - dmin) / shell;
+                const size_t idx = grid.chi.idx(i, j, k);
+                grid.chi.ptr()[idx] = 1.0 + amplitude * weight;
+            }
+        }
+    }
+}
+
+double max_shell_deviation(const Grid &grid, size_t shell) {
+    size_t i0, i1, j0, j1, k0, k1;
+    grid.domain_bounds(i0, i1, j0, j1, k0, k1);
+    double max_dev = 0.0;
+    if (shell == 0)
+        return max_dev;
+
+    for (size_t i = i0; i < i1; ++i) {
+        const size_t di = std::min(i - i0, i1 - 1 - i);
+        for (size_t j = j0; j < j1; ++j) {
+            const size_t dj = std::min(j - j0, j1 - 1 - j);
+            for (size_t k = k0; k < k1; ++k) {
+                const size_t dk = std::min(k - k0, k1 - 1 - k);
+                const size_t dmin = std::min({di, dj, dk});
+                if (dmin >= shell)
+                    continue;
+                const size_t idx = grid.chi.idx(i, j, k);
+                const double dev = std::abs(static_cast<double>(grid.chi.ptr()[idx]) - 1.0);
+                max_dev = std::max(max_dev, dev);
+            }
+        }
+    }
+    return max_dev;
+}
+
+void perturb_metric_shell(Grid &grid, double amplitude, size_t shell) {
+    size_t i0, i1, j0, j1, k0, k1;
+    grid.domain_bounds(i0, i1, j0, j1, k0, k1);
+    if (shell == 0)
+        return;
+
+    using tensorium_RG::XX;
+    using tensorium_RG::YY;
+    using tensorium_RG::ZZ;
+
+    for (size_t i = i0; i < i1; ++i) {
+        const size_t di = std::min(i - i0, i1 - 1 - i);
+        for (size_t j = j0; j < j1; ++j) {
+            const size_t dj = std::min(j - j0, j1 - 1 - j);
+            for (size_t k = k0; k < k1; ++k) {
+                const size_t dk = std::min(k - k0, k1 - 1 - k);
+                const size_t dmin = std::min({di, dj, dk});
+                if (dmin >= shell)
+                    continue;
+                const double weight = static_cast<double>(shell - dmin) / shell;
+                const double pulse = amplitude * weight;
+                const double gxx = 1.0 + pulse;
+                const double gyy = 1.0 + pulse;
+                const double gzz = 1.0 - 2.0 * pulse;
+                const size_t idx = grid.gamma_tilde[XX].idx(i, j, k);
+                grid.gamma_tilde[XX].ptr()[idx] = gxx;
+                grid.gamma_tilde[YY].ptr()[idx] = gyy;
+                grid.gamma_tilde[ZZ].ptr()[idx] = gzz;
+                grid.gamma_tilde_inv[XX].ptr()[idx] = 1.0 / gxx;
+                grid.gamma_tilde_inv[YY].ptr()[idx] = 1.0 / gyy;
+                grid.gamma_tilde_inv[ZZ].ptr()[idx] = 1.0 / gzz;
+            }
+        }
+    }
+}
+
+REGISTER_TEST("bssn.bc.outgoing_wave", "Sponge boundary damps outgoing chi pulses", []() {
+    const size_t nx = 40;
+    const size_t ny = 40;
+    const size_t nz = 40;
+    const size_t ng = 8;
+    const double spacing = 0.5;
+    const size_t shell = 6;
+    const double amplitude = 0.1;
+
+    Grid clamp(nx, ny, nz, ng, spacing, spacing, spacing);
+    Grid sponge(nx, ny, nz, ng, spacing, spacing, spacing);
+
+    auto configure = [&](Grid &grid) {
+        center_grid(grid);
+        tensorium_RG::init::minkowski(grid, 0.0);
+        seed_chi_shell(grid, amplitude, shell);
+    };
+
+    configure(clamp);
+    configure(sponge);
+
+    tensorium_RG::bssn::apply_halos_grid<tensorium_RG::bssn::BoundaryClamp>(clamp);
+    tensorium_RG::bssn::apply_halos_grid<tensorium_RG::bssn::BoundarySponge>(sponge);
+
+    const double clamp_dev = max_shell_deviation(clamp, shell);
+    const double sponge_dev = max_shell_deviation(sponge, shell);
+
+    TENSORIUM_TEST_ASSERT(clamp_dev > 0.5 * amplitude);
+    TENSORIUM_TEST_ASSERT(sponge_dev < 0.35 * clamp_dev);
+});
+
+REGISTER_TEST("bssn.bc.by_boundary_localization",
+               "Determinant drift localizes near the computational boundary", []() {
+                   const size_t nx = 48;
+                   const size_t ny = 48;
+                   const size_t nz = 48;
+                   const size_t ng = 8;
+                   const double spacing = 0.5;
+                   const size_t shell = 4;
+                   const double amplitude = 0.2;
+
+                   Grid grid(nx, ny, nz, ng, spacing, spacing, spacing);
+                   center_grid(grid);
+                   tensorium_RG::init::minkowski(grid, 0.0);
+                   perturb_metric_shell(grid, amplitude, shell);
+
+                   const auto stats = tensorium_RG::bssn::compute_invariant_stats(grid, 1);
+                   TENSORIUM_TEST_ASSERT(stats.max_det_deviation > 0.1 * amplitude);
+
+                   size_t i0, i1, j0, j1, k0, k1;
+                   grid.domain_bounds(i0, i1, j0, j1, k0, k1);
+                   const auto dist = [&](size_t idx, size_t lo, size_t hi) {
+                       return std::min(idx - lo, hi - 1 - idx);
+                   };
+                   const size_t di = dist(stats.det_location.i, i0, i1);
+                   const size_t dj = dist(stats.det_location.j, j0, j1);
+                   const size_t dk = dist(stats.det_location.k, k0, k1);
+                   const size_t min_dist = std::min({di, dj, dk});
+                   TENSORIUM_TEST_ASSERT(min_dist < shell);
+               });
+
+} // namespace

Tests/bssn/bowen_york_boost_stability.cpp

@@ -0,0 +1,192 @@
+#include "../framework/TestRegistry.hpp"
+#include "stability_common.hpp"
+#include <array>
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+
+#include "../../includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/Fields/BSSNGridSoA.hpp"
+#include "../../includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/InitialData/BSSNInitialData.hpp"
+#include "../../includes/Tensorium/Physics/DiffGeometry/BSSN_Grid/TimeIntegration/BSSNRK4.hpp"
+
+namespace {
+
+void export_slice_csv(const tensorium::tests::Grid &grid, size_t step,
+                      const std::string &output_dir) {
+    std::stringstream ss;
+    ss << output_dir << "/slice_" << std::setw(4) << std::setfill('0') << step << ".csv";
+
+    std::ofstream file(ss.str());
+    if (!file.is_open())
+        return;
+
+    file << "x,y,alpha,W,mask\n";
+
+    const size_t nx = grid.dims.nx;
+    const size_t ny = grid.dims.ny;
+    const size_t ng = grid.dims.ng;
+    const size_t k = ng + grid.dims.nz / 2;
+
+    for (size_t i = ng; i < nx + ng; ++i) {
+        for (size_t j = ng; j < ny + ng; ++j) {
+            const double x = grid.x0 + (double(i) - double(ng)) * grid.dx;
+            const double y = grid.y0 + (double(j) - double(ng)) * grid.dy;
+
+            const size_t idx = grid.alpha.idx(i, j, k);
+
+            const double alpha = grid.alpha.ptr()[idx];
+            const double chi = grid.chi.ptr()[idx];
+            const double W = std::sqrt(std::max(chi, 1e-16));
+
+            const double mask = (alpha < 0.1) ? 1.0 : 0.0;
+
+            file << x << "," << y << "," << alpha << "," << W << "," << mask << "\n";
+        }
+    }
+}
+
+using tensorium::tests::Grid;
+
+void initialize_single_boost(Grid &grid) {
+    const double m1 = 1.0;
+    const double m2 = 1.0;
+
+    const double x1 = -3.0, y1 = 0.0, z1 = 0.0;
+    const double x2 = 3.0, y2 = 0.0, z2 = 0.0;
+
+    const double Py = 0.095;
+    const double P1[3] = {0.0, Py, 0.0};
+    const double P2[3] = {0.0, -Py, 0.0};
+
+    const double S1[3] = {0.0, 0.0, 0.1};
+    const double S2[3] = {0.0, 0.0, 0.1};
+
+    tensorium_RG::init::binary_bowen_york_puncture_init(grid, m1, x1, y1, z1, P1, S1, m2, x2, y2,
+                                                        z2, P2, S2, 1e-10);
+}
+
+} // namespace
+
+REGISTER_TEST(
+    "bssn.viz.moving_puncture", "Export CSV slices of a moving spinning black hole", []() {
+        tensorium::tests::StabilityRunConfig cfg;
+        cfg.nx = 64;
+        cfg.ny = 64;
+        cfg.nz = 64;
+        cfg.spacing = 0.1;
+        cfg.ng = 4;
+        cfg.padding = 4;
+        cfg.steps = 600;
+        cfg.cfl = 0.25;
+        cfg.gauge_factor = 1.0;
+
+        tensorium_RG::fd::set_fd_dx(cfg.spacing);
+        (void)system("mkdir -p Output/viz");
+
+        tensorium::tests::Grid grid(cfg.nx, cfg.ny, cfg.nz, cfg.ng, cfg.spacing, cfg.spacing,
+                                    cfg.spacing);
+
+        grid.x0 = -0.5 * cfg.spacing * cfg.nx + 0.5 * cfg.spacing;
+        grid.y0 = -0.5 * cfg.spacing * cfg.ny + 0.5 * cfg.spacing;
+        grid.z0 = -0.5 * cfg.spacing * cfg.nz + 0.5 * cfg.spacing;
+
+        initialize_single_boost(grid);
+
+        tensorium_RG::bssn::ProjectionConfig proj_cfg;
+        proj_cfg.padding = 0;
+        proj_cfg.renormalize_metric = true;
+        proj_cfg.project_A_tilde = true;
+        proj_cfg.recompute_inverse = true;
+
+        tensorium_RG::bssn::GaugeParameters<double> params;
+        params.eta = 3.0;
+        params.beta_B_coeff = 0.5;
+
+        tensorium_RG::bssn::project_bssn_state(grid, proj_cfg);
+
+        tensorium_RG::bssn::BSSNRKStepper<double, tensorium_RG::bssn::BoundaryRadiative> stepper(
+            grid, cfg.padding);
+        stepper.set_gauge_parameters(params);
+
+        stepper.set_snapshot_callback([](const tensorium::tests::Grid &g, size_t step) {
+            if (step % 10 == 0) {
+                std::cout << ">> Exporting slice " << step << "..." << std::endl;
+                export_slice_csv(g, step, "Output/viz");
+            }
+        });
+
+        double t = 0.0;
+
+        for (size_t n = 0; n <= cfg.steps; ++n) {
+            const double dt = tensorium_RG::bssn::compute_dt_cfl(
+                grid, tensorium::tests::make_cfl_control(cfg), cfg.padding);
+
+            stepper.step(grid, dt, n);
+            t += dt;
+
+            std::printf("dt = %.4e\n", dt);
+            std::printf("Step %zu / %zu (t=%.4f)\n", n, cfg.steps, t);
+        }
+    });
+REGISTER_TEST(
+    "bssn.init.kerr_schild", "Print formatted BSSN fields after Kerr–Schild initialization", []() {
+        tensorium::tests::StabilityRunConfig cfg;
+        cfg.nx = 64;
+        cfg.ny = 64;
+        cfg.nz = 64;
+        cfg.spacing = 0.1;
+        cfg.ng = 4;
+        cfg.steps = 60;
+        cfg.cfl = 0.15;
+        cfg.gauge_factor = 1.0;
+
+        tensorium_RG::fd::set_fd_dx(cfg.spacing);
+
+        tensorium::tests::Grid grid(cfg.nx, cfg.ny, cfg.nz, cfg.ng, cfg.spacing, cfg.spacing,
+                                    cfg.spacing);
+
+        grid.x0 = -0.5 * cfg.spacing * cfg.nx + 0.5 * cfg.spacing;
+        grid.y0 = -0.5 * cfg.spacing * cfg.ny + 0.5 * cfg.spacing;
+        grid.z0 = -0.5 * cfg.spacing * cfg.nz + 0.5 * cfg.spacing;
+
+        tensorium_RG::init::kerr_schild_single(grid, 1.0, 0.0);
+
+        size_t I0, I1, J0, J1, K0, K1;
+        grid.domain_bounds(I0, I1, J0, J1, K0, K1);
+
+        const size_t                                ic = I0 + grid.dims.nx / 2;
+        const size_t                                jc = J0 + grid.dims.ny / 2;
+        const size_t                                kc = K0 + grid.dims.nz / 2;
+        tensorium_RG::bssn::GaugeParameters<double> params;
+        params.eta = 3.0;
+        params.beta_B_coeff = 0.5;
+
+        tensorium_RG::init::print_bssn_state_at(grid, ic, jc, kc, "Kerr–Schild single BH");
+        tensorium_RG::bssn::BSSNRKStepper<double, tensorium_RG::bssn::BoundaryRadiative> stepper(
+            grid, cfg.padding);
+        stepper.set_gauge_parameters(params);
+
+        stepper.set_snapshot_callback([](const tensorium::tests::Grid &g, size_t step) {
+            if (step % 1 == 0) {
+                std::cout << ">> Exporting slice " << step << "..." << std::endl;
+                export_slice_csv(g, step, "Output/viz");
+            }
+        });
+
+        double t = 0.0;
+
+        for (size_t n = 0; n <= cfg.steps; ++n) {
+            const double dt = tensorium_RG::bssn::compute_dt_cfl(
+                grid, tensorium::tests::make_cfl_control(cfg), cfg.padding);
+
+            stepper.step(grid, dt, n);
+            t += dt;
+
+            std::printf("dt = %.4e\n", dt);
+            std::printf("Step %zu / %zu (t=%.4f)\n", n, cfg.steps, t);
+        }
+        std::cout << "Kerr–Schild init test completed.\n";
+    });