BSSN setup ready (need the constraints solver)

Author: at0m741

Tests/Matrix/Matrix.cpp

@@ -303,7 +303,7 @@ int matrix_tests() {
 
 	std::cout << "X = " << X(0) << " " << X(1) << " " << X(2) << " " << X(3) << "\n";
 
-	tensorium_RG::Metric<double> metric("kerr_schild", 1.0, 0.0);
+	tensorium_RG::Metric<double> metric("kerr_schild", 1.0, 0.9);
 	metric(X, g);
 
 	std::cout << "Metric tensor g at X = (t=0, r=10, θ=π/2, φ=0):\n";
@@ -313,7 +313,7 @@ int matrix_tests() {
 	std::cout << "Christoffel symbols Γ^λ_{μν} at X = (t=0, r=10, θ=π/2, φ=0):\n";
 	auto gamma = tensorium::compute_christoffel(X, 1e-5, g, g_inv, metric);
 	gamma.print();
-	auto R = tensorium::compute_riemann_tensor<double>(X, 1e-5, tensorium_RG::Metric<double>("kerr_schild", 1.0, 0.8));
+	auto R = tensorium::compute_riemann_tensor<double>(X, 1e-5, tensorium_RG::Metric<double>("kerr_schild", 1.0, 0.9));
 	tensorium::print_riemann_tensor(R);
 	tensorium::contract_tensor<0, 1>(R);
 	std::cout << "Riemann tensor contracted:\n";

includes/Tensorium/DiffGeometry/BSSN/BSSNConstraintsSolver.hpp

@@ -0,0 +1,90 @@
+#pragma once
+
+#include "../../Core/Tensor.hpp"
+#include "../../Core/Vector.hpp"
+#include "../Metric.hpp"
+#include <cassert>
+#include <cstddef>
+
+namespace tensorium_RG {
+	template<typename T>
+		class ConstraintSolver {
+			public:
+				/**
+				 * @brief Solve the Lichnerowicz equation to enforce the Hamiltonian constraint in BSSN.
+				 *
+				 * @param Atilde        The trace-free conformal extrinsic curvature tensor Ã_{ij} at each grid point.
+				 *                      Shape: [NX][NY][NZ][3][3]
+				 * @param g_tilde_inv   The inverse conformal metric γ̃^{ij} at each grid point.
+				 *                      Shape: [NX][NY][NZ][3][3]
+				 * @param dx, dy, dz    Grid spacings
+				 * @param max_iter      Maximum number of iterations
+				 * @param tol           Convergence threshold
+				 * @return psi          The conformal factor ψ at each grid point. Shape: [NX][NY][NZ]
+				 */
+				static tensorium::Tensor<T, 3> solveLichnerowicz(
+						const tensorium::Tensor<T, 5>& Atilde,
+						const tensorium::Tensor<T, 5>& g_tilde_inv,
+						T dx, T dy, T dz,
+						int max_iter = 1000,
+						T tol = 1e-8) {
+
+					const size_t NX = Atilde.shape()[0];
+					const size_t NY = Atilde.shape()[1];
+					const size_t NZ = Atilde.shape()[2];
+
+					tensorium::Tensor<T, 3> psi({NX, NY, NZ});
+					std::fill(psi.begin(), psi.end(), T(1));;
+
+					const T inv_dx2 = T(1) / (dx * dx);
+					const T inv_dy2 = T(1) / (dy * dy);
+					const T inv_dz2 = T(1) / (dz * dz);
+					const T factor = T(1) / (2 * (inv_dx2 + inv_dy2 + inv_dz2));
+
+					for (int it = 0; it < max_iter; ++it) {
+						T max_diff = T(0);
+
+#pragma omp parallel for collapse(3) reduction(max:max_diff)
+						for (size_t i = 1; i < NX - 1; ++i) {
+							for (size_t j = 1; j < NY - 1; ++j) {
+								for (size_t k = 1; k < NZ - 1; ++k) {
+
+									T lap = inv_dx2 * (psi(i + 1, j, k) + psi(i - 1, j, k) - 2 * psi(i, j, k))
+										+ inv_dy2 * (psi(i, j + 1, k) + psi(i, j - 1, k) - 2 * psi(i, j, k))
+										+ inv_dz2 * (psi(i, j, k + 1) + psi(i, j, k - 1) - 2 * psi(i, j, k));
+
+									// Compute ò = γ̃^{ia} γ̃^{jb} Ã_{ab} Ã_{ij}
+									T A2 = T(0);
+									for (int a = 0; a < 3; ++a) {
+										for (int b = 0; b < 3; ++b) {
+											for (int c = 0; c < 3; ++c) {
+												for (int d = 0; d < 3; ++d) {
+													A2 += g_tilde_inv(i,j,k,a,c) * g_tilde_inv(i,j,k,b,d)
+														* Atilde(i,j,k,c,d)     * Atilde(i,j,k,a,b);
+												}
+											}
+										}
+									}
+
+									T psi_val = std::fmax(psi(i, j, k), T(1e-8));
+									T rhs = (T(1.0) / 8.0) * A2 / std::pow(psi_val, 7);
+
+									T new_psi = psi(i, j, k) + T(0.5) * (lap - rhs) * factor;
+									new_psi = std::fmax(new_psi, T(0.1));
+
+									max_diff = std::fmax(max_diff, std::abs(new_psi - psi(i, j, k)));
+									psi(i, j, k) = new_psi;
+								}
+							}
+						}
+
+						if (max_diff < tol)
+							break;
+					}
+
+					return psi;
+				}
+
+		};
+
+}

includes/Tensorium/DiffGeometry/BSSN/BSSNDerivatives.hpp

@@ -122,5 +122,23 @@ namespace tensorium_RG {
 			gp2 = shifted(0, 0, 2 * dz);
 			out[2] = (-gp2 + 8 * gp1 - 8 * gm1 + gm2) / (12 * dz);
 		}
+
+	template <typename T>
+		tensorium::Tensor<T, 2> compute_dt_gamma_from_beta(
+				const tensorium::Tensor<T, 2>& gamma,
+				const tensorium::Vector<T>& beta,
+				const tensorium::Tensor<T, 2>& partial_beta,
+				const tensorium::Tensor<T, 3>& christoffel) {
+
+			tensorium::Tensor<T, 2> dtg({3, 3});
+			for (size_t i = 0; i < 3; ++i)
+				for (size_t j = 0; j < 3; ++j) {
+					T val = partial_beta(i, j) + partial_beta(j, i);
+					for (size_t k = 0; k < 3; ++k)
+						val -= 2.0 * christoffel(i, j, k) * beta(k);
+					dtg(i, j) = val;
+				}
+			return dtg;
+		}
 }
 

includes/Tensorium/DiffGeometry/BSSN/BSSNPrintDebug.hpp

@@ -0,0 +1,44 @@
+#pragma once
+
+#include "../Metric.hpp"
+#include "../../Core/Vector.hpp"
+#include "../../Core/Tensor.hpp"
+#include "BSSNSetup.hpp"
+
+namespace tensorium_RG {
+
+	template <typename T>
+		inline void print_tensor2(const std::string& name, const tensorium::Tensor<T, 2>& tensor) {
+			std::cout << "--- " << name << " ---\n";
+			for (size_t i = 0; i < 3; ++i) {
+				for (size_t j = 0; j < 3; ++j)
+					std::cout << std::setw(12) << tensor(i, j) << " ";
+				std::cout << "\n";
+			}
+		}
+
+	template <typename T>
+		inline void print_tensor3(const std::string& name, const tensorium::Tensor<T, 3>& tensor) {
+			std::cout << "--- " << name << " ---\n";
+			for (size_t d = 0; d < 3; ++d) {
+				std::cout << "∂_" << "xyz"[d] << ":\n";
+				for (size_t i = 0; i < 3; ++i) {
+					for (size_t j = 0; j < 3; ++j)
+						std::cout << std::setw(12) << tensor(d, i, j) << " ";
+					std::cout << "\n";
+				}
+			}
+		}
+
+	template <typename T>
+		inline void print_vector(const std::string& name, const tensorium::Vector<T>& vec) {
+			std::cout << "--- " << name << " ---\n";
+			for (size_t i = 0; i < vec.size(); ++i)
+				std::cout << std::setw(12) << vec(i) << " ";
+			std::cout << "\n";
+		}
+
+
+
+} // namespace tensorium_RG
+

includes/Tensorium/DiffGeometry/BSSN/BSSNSetup.hpp

@@ -19,6 +19,8 @@
 #include "BSSNAtildeTensor.hpp"
 #include "BSSNTildeChristoffel.hpp"
 #include "BSSNContractedChristoffel.hpp"
+#include "BSSNPrintDebug.hpp"
+#include "BSSNConstraintsSolver.hpp"
 
 namespace tensorium_RG {
 	/**
@@ -84,26 +86,18 @@ namespace tensorium_RG {
 				void init_BSSN(const tensorium::Vector<T>& X,
 						const tensorium_RG::Metric<T>& metric,
 						T dx, T dy, T dz) {
-
 					T alpha;
-
-					// Extract metric quantities
-
 					tensorium::Vector<T> beta(3);
-					tensorium::Tensor<T, 2> gammaj({3, 3});
-					metric.BSSN(X, alpha, beta, gammaj);
-					tensorium::Tensor<T, 3> christoffel({3, 3, 3});
-					tensorium::Tensor<T, 2> gammaj_inv = inv_mat_tensor(gammaj);
-					tensorium::Tensor<T, 2> dgt({3, 3}); 
-					tensorium::Tensor<T, 2> partial_beta({3, 3}); 
-					tensorium::Tensor<T, 2> d_beta({3, 3});
-					T chi = metric.compute_conformal_factor(gammaj);
-					tensorium::Tensor<T, 2> gamma_tilde = compute_conformal_metric(metric, gammaj, chi);
+					tensorium::Tensor<T, 2> gamma_ij({3, 3});
+					metric.BSSN(X, alpha, beta, gamma_ij);
+
+					tensorium::Tensor<T, 2> dgt({3, 3});
+					tensorium::Tensor<T, 2> gamma_ij_inv = inv_mat_tensor(gamma_ij);
+					T chi = metric.compute_conformal_factor(gamma_ij);
+					tensorium::Tensor<T, 2> gamma_tilde = compute_conformal_metric(metric, gamma_ij, chi);
 					tensorium::Tensor<T, 2> gamma_tilde_inv = inv_mat_tensor(gamma_tilde);
-					tensorium::Tensor<T, 3> dgamma_tilde({3, 3, 3});
-					tensorium_RG::ExtrinsicCurvature<T> extr;
-					tensorium_RG::BSSNAtildeTensor<T> Aij;
 
+					tensorium::Tensor<T, 3> dgamma_tilde({3, 3, 3});
 					tensorium_RG::compute_partial_derivatives_tensor2D<T>(
 							X, dx, dy, dz,
 							[&](const tensorium::Vector<T>& Xs, tensorium::Tensor<T, 2>& out) {
@@ -116,13 +110,18 @@ namespace tensorium_RG {
 							},
 							dgamma_tilde
 							);
-				
-					compute_christoffel_3D(gamma_tilde, dgamma_tilde, gamma_tilde_inv, christoffel);
-					auto contracted_Gamma = tensorium_RG::BSSNContractedGamma<T>::compute(X, metric, dx, dy, dz, chi);
+
+					tensorium::Tensor<T, 3> christoffel_tilde({3, 3, 3});
+					compute_christoffel_3D(gamma_tilde, dgamma_tilde, gamma_tilde_inv, christoffel_tilde);
 
 					tensorium::Vector<T> tilde_Gamma(3);
-					TildeGamma<T>::compute(gamma_tilde_inv, christoffel, tilde_Gamma);
+					TildeGamma<T>::compute(gamma_tilde_inv, christoffel_tilde, tilde_Gamma);
 
+					auto contracted_Gamma = tensorium_RG::BSSNContractedGamma<T>::compute(
+							X, metric, dx, dy, dz, chi
+							);
+
+					tensorium::Tensor<T, 2> d_beta({3, 3});
 					tensorium_RG::compute_partial_derivatives_vector<T>(
 							X, dx, dy, dz,
 							[&](const tensorium::Vector<T>& Xs, tensorium::Vector<T>& out) {
@@ -135,31 +134,73 @@ namespace tensorium_RG {
 							d_beta
 							);
 
-					for (size_t i = 0; i < 3; ++i)
-						for (size_t j = 0; j < 3; ++j) {
-							T val = 0;
-							for (size_t l = 0; l < 3; ++l)
-								val += gammaj(j, l) * d_beta(i, l); 
-							partial_beta(i, j) = val;
-						}
-
-					auto Kij = extr.compute_Kij(dgt, gammaj, beta, partial_beta, christoffel, alpha);
-					auto AtildeTensor = Aij.compute_Atilde_tensor(Kij, gammaj_inv, gammaj, chi);
-					
-
-					grid.tilde_Gamma = {tilde_Gamma};
-					grid.contracted_Gamma = {contracted_Gamma};
-					grid.ExtrinsicTensor = {Kij};
-					grid.A_tildeTensor = {AtildeTensor};
-					grid.alpha = {alpha};
-					grid.beta = {beta};
-					grid.gamma_ij = {gammaj};
-					grid.gamma_ij_inv = {gammaj_inv};
-					grid.chi = {chi};
-					grid.gamma_tilde = {gamma_tilde};
-					grid.gamma_tilde_inv = {gamma_tilde_inv};
-					grid.dgamma_tilde = {dgamma_tilde};
-					grid.christoffel_tilde = {christoffel};
+					tensorium::Tensor<T, 3> dgamma_phys({3, 3, 3});
+					tensorium_RG::compute_partial_derivatives_tensor2D<T>(
+							X, dx, dy, dz,
+							[&](const tensorium::Vector<T>& Xs, tensorium::Tensor<T, 2>& out) {
+							T a_tmp;
+							tensorium::Vector<T> b_tmp(3);
+							tensorium::Tensor<T, 2> g_tmp({3, 3});
+							metric.BSSN(Xs, a_tmp, b_tmp, g_tmp);
+							out = g_tmp;
+							},
+							dgamma_phys
+							);
+
+					tensorium::Tensor<T, 3> christoffel_phys({3, 3, 3});
+					compute_christoffel_3D(gamma_ij, dgamma_phys, gamma_ij_inv, christoffel_phys);
+
+					tensorium_RG::ExtrinsicCurvature<T> extr;
+					auto Kij = extr.compute_Kij(dgt, gamma_ij, beta, d_beta, christoffel_phys, alpha);
+
+					tensorium_RG::BSSNAtildeTensor<T> Aij;
+					auto AtildeTensor = Aij.compute_Atilde_tensor(Kij, gamma_ij_inv, gamma_ij, chi);
+
+					grid.alpha             = {alpha};
+					grid.beta              = {beta};
+					grid.gamma_ij          = {gamma_ij};
+					grid.gamma_ij_inv      = {gamma_ij_inv};
+					grid.chi               = {chi};
+					grid.gamma_tilde       = {gamma_tilde};
+					grid.gamma_tilde_inv   = {gamma_tilde_inv};
+					grid.dgamma_tilde      = {dgamma_tilde};
+					grid.christoffel_tilde = {christoffel_tilde};
+					grid.tilde_Gamma       = {tilde_Gamma};
+					grid.contracted_Gamma  = {contracted_Gamma};
+					grid.ExtrinsicTensor   = {Kij};
+					grid.A_tildeTensor     = {AtildeTensor};
+
+					std::cout << std::setprecision(6) << std::fixed;
+					std::cout << "\n========= BSSN Quantities at X =========\n";
+
+					std::cout << "--- alpha ---\n" << grid.alpha[0] << "\n";
+					print_vector("beta", grid.beta[0]);
+
+					std::cout << "dbeta\n";
+					d_beta.print();
+
+					print_tensor2("gamma_ij", grid.gamma_ij[0]);
+					print_tensor2("gamma_ij_inv", grid.gamma_ij_inv[0]);
+
+					std::cout << "--- chi ---\n" << grid.chi[0] << "\n";
+					print_tensor2("gamma_tilde", grid.gamma_tilde[0]);
+					print_tensor2("gamma_tilde_inv", grid.gamma_tilde_inv[0]);
+
+					print_tensor3("dgamma_tilde", grid.dgamma_tilde[0]);
+					print_tensor3("christoffel_tilde", grid.christoffel_tilde[0]);
+					print_vector("tilde_Gamma", grid.tilde_Gamma[0]);
+
+					print_tensor2("∂_t gamma_ij (dgt)", dgt);
+					print_vector("contracted_Gamma", grid.contracted_Gamma[0]);
+
+					print_tensor3("dgamma_phys", dgamma_phys);
+					print_tensor3("christoffel_phys", christoffel_phys);
+					print_tensor2("Extrinsic curvature K_ij", grid.ExtrinsicTensor[0]);
+					Kij.print();
+
+					print_tensor2("A_tilde_ij", grid.A_tildeTensor[0]);
+
+					std::cout << "========================================\n\n";
 				}
 		};
 

includes/Tensorium/DiffGeometry/BSSN/BSSNextrinTensor.hpp

@@ -83,8 +83,9 @@ namespace tensorium_RG {
 							for (size_t j = 0; j < 3; ++j) {
 								T sym_grad_beta = partial_beta(i, j) + partial_beta(j, i);
 								T gamma_beta = T(0.0);
+							
 								for (size_t k = 0; k < 3; ++k)
-									gamma_beta += 2.0 * christoffel(i, j, k) * beta(k);
+									gamma_beta += 2.0 * christoffel(k, i, j) * beta(k);
 
 								Kij(i, j) = -0.5 / alpha * (dgt(i, j) - (sym_grad_beta - gamma_beta));
 							}