Function geodesics

examples/pde/bhrt.cpp:588–620

* @brief Computes the geodesics from the established metric, 4-vector positions * and velocities * * @param pos4 * @param vel4 * @return af::array */

 * @return af::array
 */
af::array geodesics(const af::array& pos4, const af::array& vel4) {
    auto N = vel4.dims()[1];

    af::array uu = af::matmul(af::moddims(vel4, af::dim4(4, 1, N)),
                              af::moddims(vel4, af::dim4(1, 4, N)));
    uu           = af::moddims(uu, af::dim4(1, 4, 4, N));

    af::array metric    = metric4(pos4);
    af::array invmetric = af::moddims(inv_metric(metric), af::dim4(4, 4, 1, N));

    // Compute the partials of the metric with respect to coordinates indices
    af::array dt = af::constant(0, 4, 4, 1, N, f64);

    auto dr     = partials(pos4, 1, 1e-6, 1e-12);
    auto dtheta = partials(pos4, 2, 1e-6, 1e-12);
    auto dphi   = partials(pos4, 3, 1e-6, 1e-12);

    dr     = af::moddims(dr, af::dim4(4, 4, 1, N));
    dtheta = af::moddims(dtheta, af::dim4(4, 4, 1, N));
    dphi   = af::moddims(dphi, af::dim4(4, 4, 1, N));

    // Compute the einsum for each of the christoffel terms
    af::array partials = af::join(2, dt, dr, dtheta, dphi);
    af::array p1       = af::matmul(invmetric, partials);
    af::array p2       = af::reorder(p1, 0, 2, 1, 3);
    af::array p3 = af::matmul(invmetric, af::reorder(partials, 2, 0, 1, 3));

    auto christoffels = -0.5 * (p1 + p2 - p3);

    // Use the geodesics equation to find the 4-vector acceleration
    return af::moddims(af::sum(af::sum(christoffels * uu, 1), 2),
                       af::dim4(4, N));
}

/**
 * @brief Camera struct