Christoffel contraction

Author: at0m741

includes/Tensorium/DiffGeometry/BSSN/BSSNContractedChristoffel.hpp

@@ -0,0 +1,53 @@
+#pragma once
+
+#include "../Metric.hpp"
+#include "BSSNDerivatives.hpp"
+
+namespace tensorium_RG {
+
+	template<typename T>
+		class BSSNContractedGamma {
+			public:
+				/**
+				 * @brief Compute the contracted Christoffel symbol \f$ \Gamma^i_{ij} = -\frac{3}{2} \partial_j \ln \chi \f$
+				 *
+				 * This uses the relation:
+				 * \f[
+				 * \Gamma^i_{ij} = -\frac{3}{2} \frac{\partial_j \chi}{\chi}
+				 * \f]
+				 *
+				 * @param X     Spatial coordinates
+				 * @param metric Metric object implementing `BSSN(X, alpha, beta, gamma)`
+				 * @param dx,dy,dz Grid spacing
+				 * @param chi   Conformal factor \f$ \chi \f$
+				 * @return Vector containing \f$ \Gamma^i_{ij} \f$ for each j
+				 */
+				static tensorium::Vector<T> compute(
+						const tensorium::Vector<T>& X,
+						const tensorium_RG::Metric<T>& metric,
+						T dx, T dy, T dz,
+						T chi)
+				{
+					tensorium::Vector<T> d_chi(3);
+
+					tensorium_RG::compute_partial_derivatives_scalar<T>(
+							X, dx, dy, dz,
+							[&](const tensorium::Vector<T>& Xs) {
+							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);
+							return compute_conformal_factor(metric, g_tmp);
+							},
+							d_chi
+							);
+
+					tensorium::Vector<T> out(3);
+					for (size_t j = 0; j < 3; ++j)
+						out[j] = -1.5 * d_chi[j] / chi;
+
+					return out;
+				}
+		};
+
+} // namespace tensorium_RG

includes/Tensorium/DiffGeometry/BSSN/BSSNDerivatives.hpp

@@ -86,5 +86,41 @@ namespace tensorium_RG {
 			}
 		}
 
+	template<typename T, typename ScalarFunc>
+		void compute_partial_derivatives_scalar(
+				const tensorium::Vector<T>& X,
+				T dx, T dy, T dz,
+				ScalarFunc&& func,
+				tensorium::Vector<T>& out) {
+
+			out.resize(3);
+
+			auto shifted = [&](T dx_, T dy_, T dz_) {
+				tensorium::Vector<T> Xs = X;
+				Xs(0) += dx_; Xs(1) += dy_; Xs(2) += dz_;
+				return func(Xs);
+			};
+
+			// ∂/∂x
+			T gm2 = shifted(-2 * dx, 0, 0);
+			T gm1 = shifted(-dx, 0, 0);
+			T gp1 = shifted(dx, 0, 0);
+			T gp2 = shifted(2 * dx, 0, 0);
+			out[0] = (-gp2 + 8 * gp1 - 8 * gm1 + gm2) / (12 * dx);
+
+			// ∂/∂y
+			gm2 = shifted(0, -2 * dy, 0);
+			gm1 = shifted(0, -dy, 0);
+			gp1 = shifted(0, dy, 0);
+			gp2 = shifted(0, 2 * dy, 0);
+			out[1] = (-gp2 + 8 * gp1 - 8 * gm1 + gm2) / (12 * dy);
+
+			// ∂/∂z
+			gm2 = shifted(0, 0, -2 * dz);
+			gm1 = shifted(0, 0, -dz);
+			gp1 = shifted(0, 0, dz);
+			gp2 = shifted(0, 0, 2 * dz);
+			out[2] = (-gp2 + 8 * gp1 - 8 * gm1 + gm2) / (12 * dz);
+		}
 }
 

includes/Tensorium/DiffGeometry/BSSN/BSSNSetup.hpp

@@ -17,6 +17,8 @@
 #include "BSSNDerivatives.hpp"
 #include "BSSNextrinTensor.hpp"
 #include "BSSNAtildeTensor.hpp"
+#include "BSSNTildeChristoffel.hpp"
+#include "BSSNContractedChristoffel.hpp"
 
 namespace tensorium_RG {
 	/**
@@ -47,6 +49,8 @@ namespace tensorium_RG {
 		std::vector<tensorium::Tensor<double, 3>> christoffel_tilde; ///< Christoffel symbols \f$\tilde{\Gamma}^k_{ij}\f$
 		std::vector<tensorium::Tensor<double, 2>> ExtrinsicTensor;   ///< Extrinsic curvature \f$K_{ij}\f$
 		std::vector<tensorium::Tensor<double, 2>> A_tildeTensor;     ///< Trace-free extrinsic curvature \f$\tilde{A}_{ij}\f$
+		std::vector<tensorium::Vector<double>> tilde_Gamma;			 ///< Conformal contracted symbols \f$\tilde{\Gamma}^i\f$
+		std::vector<tensorium::Vector<double>> contracted_Gamma;	 ///< Contracted symbols \f$\Gamma^i_{ij} = -\frac{3}{2} \partial_j \ln \chi\f$
 	};
 	/**
 	 * @class BSSN
@@ -112,8 +116,12 @@ 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::Vector<T> tilde_Gamma(3);
+					TildeGamma<T>::compute(gamma_tilde_inv, christoffel, tilde_Gamma);
 
 					tensorium_RG::compute_partial_derivatives_vector<T>(
 							X, dx, dy, dz,
@@ -137,7 +145,10 @@ namespace tensorium_RG {
 
 					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};