BSSNRicciTildeTensor.hpp

Go to the documentation of this file.

 
#pragma once
 
#include "../../Core/Matrix.hpp"
#include "../../Core/Tensor.hpp"
#include "../../Core/Vector.hpp"
#include "../Metric.hpp"
#include "BSSNAutoDiff.hpp"
#include "BSSNChristoffel.hpp"
#include "BSSNContractedChristoffel.hpp"
#include "BSSNChiContext.hpp"
#include "BSSNDerivatives.hpp"
 
namespace tensorium_RG {
template <typename T> class RicciTildeTensor {
  public:
    static tensorium::Tensor<T, 2>
    compute_laplacian_term(const tensorium::Vector<T> &X, T dx, T dy, T dz, const Metric<T> &metric,
                           const tensorium::Tensor<T, 2> &gamma_tilde_inv) {
        using namespace tensorium;
 
        Tensor<T, 2> result({3, 3});
 
        for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
 
                auto tensor_func = [&](const Vector<T> &Xs, Tensor<T, 2> &out_tensor) {
                    T            a;
                    Vector<T>    b(3);
                    Tensor<T, 2> g_phys({3, 3});
                    metric.BSSN(Xs, a, b, g_phys);
                    T            chi = compute_conformal_factor(metric, g_phys);
                    Tensor<T, 2> gtilde = compute_conformal_metric(metric, g_phys, chi);
                    out_tensor(i, j) = gtilde(i, j);
                };
 
                Tensor<T, 4> d2_full({1, 1, 3, 3});
                compute_second_derivatives_tensor2D<T>(X, dx, dy, dz, tensor_func, d2_full);
 
                Tensor<T, 2> d2gamma({3, 3});
                for (int k = 0; k < 3; ++k)
                    for (int l = 0; l < 3; ++l)
                        d2gamma(k, l) = d2_full(0, 0, k, l);
 
                T laplace = 0.0;
                for (int k = 0; k < 3; ++k)
                    for (int l = 0; l < 3; ++l)
                        laplace += gamma_tilde_inv(k, l) * d2gamma(k, l);
 
                result(i, j) = -0.5 * laplace;
            }
        }
 
        for (size_t i = 0; i < 3; ++i) {
            for (size_t j = i + 1; j < 3; ++j) {
                T avg = 0.5 * (result(i, j) + result(j, i));
                result(i, j) = result(j, i) = avg;
            }
        }
 
        return result;
    }
 
    static tensorium::Tensor<T, 2>
    compute_dGamma_term(const tensorium::Vector<T> &X, T dx, T dy, T dz,
                        const tensorium::Vector<T>    &tilde_Gamma,
                        const tensorium::Tensor<T, 2> &tilde_gamma_contract) {
        using namespace tensorium;
 
        Tensor<T, 2> dGamma({3, 3});
        for (int k = 0; k < 3; ++k) {
            auto      scalar_func = [&](const Vector<T> &Xs) -> T { return tilde_Gamma(k); };
            Vector<T> grad_k = partial_scalar(X, dx, dy, dz, scalar_func);
            for (int j = 0; j < 3; ++j) {
                dGamma(j, k) = grad_k(j);
            }
        }
 
        Tensor<T, 2> R2({3, 3});
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                T sum = 0;
                for (int k = 0; k < 3; ++k) {
                    sum += tilde_gamma_contract(k, i) * dGamma(j, k) +
                           tilde_gamma_contract(k, j) * dGamma(i, k);
                }
                R2(i, j) = static_cast<T>(0.5) * sum;
            }
        }
 
        return R2;
    }
 
    static tensorium::Tensor<T, 2>
    compute_GammaGamma_term(const tensorium::Vector<T>    &tilde_Gamma,
                            const tensorium::Tensor<T, 3> &christoffel_tilde,
                            const tensorium::Tensor<T, 2> &gamma_tilde) {
        using namespace tensorium;
        Tensor<T, 2> R3({3, 3});
 
        for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
                T sum_ij = 0;
                for (size_t k = 0; k < 3; ++k) {
                    T G_ijk = 0, G_jik = 0;
                    for (size_t m = 0; m < 3; ++m) {
                        G_ijk += gamma_tilde(i, m) * christoffel_tilde(m, j, k);
                        G_jik += gamma_tilde(j, m) * christoffel_tilde(m, i, k);
                    }
                    sum_ij += tilde_Gamma(k) * (G_ijk + G_jik);
                }
                R3(i, j) = T(0.5) * sum_ij;
            }
        }
 
        return R3;
    }
 
    static tensorium::Tensor<T, 2>
    compute_GammaProduct_term(const tensorium::Tensor<T, 2> &gamma_tilde_inv,
                              const tensorium::Tensor<T, 3> &christoffel_tilde) {
        using namespace tensorium;
        Tensor<T, 2> R4({3, 3});
 
        for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
                T sum = 0;
                for (size_t ell = 0; ell < 3; ++ell) {
                    for (size_t m = 0; m < 3; ++m) {
                        T term1 = 0;
                        for (size_t k = 0; k < 3; ++k) {
                            term1 += christoffel_tilde(k, ell, i) * christoffel_tilde(j, k, m);
                        }
                        for (size_t k = 0; k < 3; ++k) {
                            term1 += christoffel_tilde(k, ell, j) * christoffel_tilde(i, k, m);
                        }
                        T term2 = 0;
                        for (size_t k = 0; k < 3; ++k) {
                            term2 += christoffel_tilde(k, i, m) * christoffel_tilde(j, k, ell);
                        }
                        sum += gamma_tilde_inv(ell, m) * (term1 + term2);
                    }
                }
                R4(i, j) = sum;
            }
        }
        return R4;
    }
    static tensorium::Tensor<T, 2> compute_Ricci_Tilde_tensor(
        const ChiContext<T> &chi_context, const tensorium::Tensor<T, 2> &gamma_tilde_inv,
        const tensorium::Vector<T> &tilde_Gamma, const tensorium::Tensor<T, 3> &christoffel_tilde,
        const tensorium::Tensor<T, 2> &gamma_tilde) {
        using namespace tensorium;
        const auto& X = chi_context.X;
        const T dx = chi_context.dx;
        const T dy = chi_context.dy;
        const T dz = chi_context.dz;
        const auto& metric = chi_context.metric;
        auto R1 = compute_laplacian_term(X, dx, dy, dz, metric, gamma_tilde_inv);
        auto R2 = compute_dGamma_term(X, dx, dy, dz, tilde_Gamma, gamma_tilde);
        auto R3 = compute_GammaGamma_term(tilde_Gamma, christoffel_tilde, gamma_tilde);
        auto R4 = compute_GammaProduct_term(gamma_tilde_inv, christoffel_tilde);
        tensorium::Tensor<T, 2> Ricci({3, 3});
        for (size_t i = 0; i < 3; ++i) {
            for (size_t j = 0; j < 3; ++j) {
                Ricci(i, j) = R1(i, j) + R2(i, j) + R3(i, j) + R4(i, j);
            }
        }
 
        return Ricci;
    }
};
 
} // namespace tensorium_RG