From 2af27f4655cd109801082f5313a7f6ce4a73dd46 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Tue, 13 Jan 2026 22:29:46 +0900 Subject: [PATCH 01/18] add Kerr basis --- src/shamphys/include/shamphys/metrics.hpp | 432 ++++++++++++++++++++++ 1 file changed, 432 insertions(+) create mode 100644 src/shamphys/include/shamphys/metrics.hpp diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp new file mode 100644 index 0000000000..0e8dbffa92 --- /dev/null +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -0,0 +1,432 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2025 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +#pragma once + +/** + * @file metrics.hpp + * @brief + */ + +#include "shamunits/Constants.hpp" +#include + +namespace shamphys { + + // The Boyer-Lindquist metric tensor in CARTESIAN-like form + template + inline Tvec get_cartesian_cov_metric(Tvec pos) { + + // metric params + Tscal a = 1; + Tscal mass1 = 1.; + Tscal rs = 2. * mass1; + + // pos coords + Tscal x, y, z = pos; + Tscal x2 = x * x; + Tscal y2 = y * y; + Tscal z2 = z * z; + Tscal a2 = a * a; + Tscal r2_spher = x2 + y2 + z2; + Tscal r2 = 0.5 * (r2_spher - a2) + + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); + Tscal r = sycl::sqrt(r2); + Tscal inv_r2 = 1. / r2; + Tscal rho2 = r2 + a2 * (z2 * inv_r2); + Tscal inv_rho2 = 1. / rho2; + Tscal a2pr2 = a2 + r2; + + // metric components + Tscal delta = a2pr2 - rs * r; + Tscal sintheta2 = 1. - z2 * inv_r2; + Tscal gtt = -1. + (r * rs) * inv_rho2; + Tscal gphiphi = sintheta2 * (a2pr2 + (a2 * r * rs * sintheta2) * inv_rho2); + Tscal gtphi = -((a * r * rs * sintheta2) * inv_rho2); + Tscal omega = a * r * rs / (rho2 * a2pr2 + rs * r * a2 * sintheta2); + + Tscal inv_rho2delta = inv_rho2 / delta; + Tscal gtphi_on_x2py2 = -rs * r * a / (rho2 * a2pr2); + + // let the fun begin + Tvec g(4, 4); + g(0, 0) = gtt; + g(1, 0) = -y * gtphi_on_x2py2; + g(2, 0) = x * gtphi_on_x2py2; + g(3, 0) = 0.; + g(0, 1) = g(1, 0); + g(1, 1) = 1. + + (r * r * r * (a2pr2) *rs * x2 + a2 * delta * r * rs * y2) + / (delta * (a2pr2) * (a2pr2) *rho2); + g(2, 1) = (r * (r * r * r * r + a2 * (-delta + r2)) * rs * x * y) + / (delta * (a2pr2) * (a2pr2) *rho2); + g(3, 1) = (r * rs * x * z) / (delta * rho2); + g(0, 2) = g(2, 0); + g(1, 2) = g(2, 1); + g(2, 2) = 1. + + (a2 * delta * r * rs * x2 + r * r * r * (a2pr2) *rs * y2) + / (delta * (a2pr2) * (a2pr2) *rho2); + g(3, 2) = (r * rs * y * z) / (delta * rho2); + g(0, 3) = g(3, 0); + g(1, 3) = g(3, 1); + g(2, 3) = g(3, 2); + g(3, 3) = 1. + ((a2 + r2) * rs * z2) / (delta * r * rho2); + + return g; + } + + template + inline Tvec get_cartesian_contrav_metric(Tvec pos) { + + // metric params + Tscal a = 1; + Tscal mass1 = 1.; + Tscal rs = 2. * mass1; + + // pos coords + Tscal x, y, z = pos; + Tscal x2 = x * x; + Tscal y2 = y * y; + Tscal z2 = z * z; + Tscal a2 = a * a; + Tscal r2_spher = x2 + y2 + z2; + Tscal r2 = 0.5 * (r2_spher - a2) + + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); + Tscal r = sycl::sqrt(r2); + Tscal inv_r2 = 1. / r2; + Tscal rho2 = r2 + a2 * (z2 * inv_r2); + Tscal inv_rho2 = 1. / rho2; + Tscal a2pr2 = a2 + r2; + + // metric components + Tscal delta = a2pr2 - rs * r; + Tscal sintheta2 = 1. - z2 * inv_r2; + Tscal gtt = -1. + (r * rs) * inv_rho2; + Tscal gphiphi = sintheta2 * (a2pr2 + (a2 * r * rs * sintheta2) * inv_rho2); + Tscal gtphi = -((a * r * rs * sintheta2) * inv_rho2); + Tscal omega = a * r * rs / (rho2 * a2pr2 + rs * r * a2 * sintheta2); + + Tscal inv_rho2delta = inv_rho2 / delta; + Tscal gtphi_on_x2py2 = -rs * r * a / (rho2 * a2pr2); + + // let the fun begin + Tvec g(4, 4); + Tscal domegaterm = 1. / (omega * gtphi + gtt); + g(0, 0) = domegaterm; + g(1, 0) = -y * omega * domegaterm; + g(2, 0) = x * omega * domegaterm; + g(3, 0) = 0.; + g(0, 1) = g(1, 0); + // NOTE from Phantom: the expressions below are NOT regular at x=y=0 and z=r. Needs fixing! + // Or use inv4x4 instead... + g(1, 1) = (delta * r2 * x2) / (a2pr2 * a2pr2 * rho2) + + (gtt * y2) / (-gtphi * gtphi + gphiphi * gtt) + (x2 * z2) / (rho2 * (r2 - z2)); + g(2, 1) = -((gtt * x * y) / (-gtphi * gtphi + gphiphi * gtt)) + + (delta * r2 * x * y) / (a2pr2 * a2pr2 * rho2) + + (x * y * z2) / (rho2 * (r2 - z2)); + g(3, 1) = -((x * z) * inv_rho2) + (delta * x * z) / (a2pr2 * rho2); + g(0, 2) = g(2, 0); + g(1, 2) = g(2, 1); + g(2, 2) = (gtt * x2) / (-gtphi * gtphi + gphiphi * gtt) + + (delta * r2 * y2) / (a2pr2 * a2pr2 * rho2) + (y2 * z2) / (rho2 * (r2 - z2)); + g(3, 2) = -((y * z) * inv_rho2) + (delta * y * z) / (a2pr2 * rho2); + g(0, 3) = g(3, 0); + g(1, 3) = g(3, 1); + g(2, 3) = g(3, 2); + g(3, 3) = (r2 - z2) * inv_rho2 + (delta * z2) / (r2 * rho2); + + return g; + } + + template + inline void metric_cartesian_derivatives( + Tvec pos, Tvec &dgcovdx, Tvec &dgcovdy, Tvec &dgcovdz) { + + // metric params + Tscal a = 1; + Tscal mass1 = 1.; + Tscal rs = 2. * mass1; + + // position coords + Tscal x = pos(0); + Tscal y = pos(1); + Tscal z = pos(2); + + Tscal x2 = x * x; + Tscal y2 = y * y; + Tscal z2 = z * z; + Tscal a2 = a * a; + + Tscal r2_spher = x2 + y2 + z2; + Tscal r2 = 0.5 * (r2_spher - a2) + + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); + Tscal r = sycl::sqrt(r2); + + Tscal rho2 = r2 + a2 * (z2 / r2); + Tscal a2pr2 = a2 + r2; + Tscal delta = a2pr2 - rs * r; + + // Metric components + Tscal sintheta2 = 1. - z2 / r2; + Tscal gtt = -1. + (r * rs) / rho2; + Tscal gphiphi = sintheta2 * (a2 + r2 + (a2 * r * rs * sintheta2) / rho2); + Tscal gtphi = -((a * r * rs * sintheta2) / rho2); + + // inverses + Tscal r21 = 1. / r2; + Tscal rho21 = 1. / rho2; + Tscal delta1 = 1. / delta; + Tscal r2mz2 = r2 - z2; + Tscal r2mz21 = 1. / r2mz2; + Tscal x2py2 = x2 + y2; + Tscal x2py21 = 1. / x2py2; + + // terms + Tscal term1 = a2pr2 * rho21; + Tscal term2 = (2. * a2pr2 * a2pr2 * rho21 * rho21 * z * z * z) / (r2 * r2); + + // Derivatives + Tscal dr2dx = (2. * r2 * x) / rho2; + Tscal dr2dy = (2. * r2 * y) / rho2; + Tscal dr2dz = (2. * a2pr2 * z) / rho2; + + Tscal drho2dx = 2. * r21 * rho21 * x * (r2 * r2 - a2 * z2); + Tscal drho2dy = 2. * r21 * rho21 * y * (r2 * r2 - a2 * z2); + Tscal drho2dz = 2. * r21 * r21 * z * (a2 * r2 + a2pr2 * rho21 * (r2 * r2 - a2 * z2)); + + Tscal ddeltadx = r * rho21 * (2. * r - rs) * x; + Tscal ddeltady = r * rho21 * (2. * r - rs) * y; + Tscal ddeltadz = (a2pr2 * rho21 * (2. * r - rs) * z) / r; + + // metric components + Tscal dgtphidx = (a * rho21 * rho21 * rs * (drho2dx * r2mz2 - x * (r2 + z2))) / r; + Tscal dgtphidy = (a * rho21 * rho21 * rs * (drho2dy * r2mz2 - y * (r2 + z2))) / r; + Tscal dgtphidz + = a * drho2dz * r * rho21 * rho21 * rs * sintheta2 + - (a * a2pr2 * rho21 * rho21 * rs * sintheta2 * z) / r + - a * r * rho21 * rs * (-2. * r21 * z + 2. * a2pr2 * r21 * r21 * rho21 * z * z * z); + + Tscal dsintheta2dx = 2. * r21 * rho21 * x * z2; + Tscal dsintheta2dy = 2. * r21 * rho21 * y * z2; + Tscal dsintheta2dz = -2. * r21 * z + 2. * a2pr2 * r21 * r21 * rho21 * z * z * z; + + Tscal dgphiphidx + = dsintheta2dx * r2 + dr2dx * sintheta2 + + a2 * rho21 * rho21 + * (dsintheta2dx * rho2 * rho2 - drho2dx * r * rs * sintheta2 * sintheta2 + + r * rs * sintheta2 * (2. * dsintheta2dx * rho2 + sintheta2 * x)); + + Tscal dgphiphidy + = dsintheta2dy * r2 + dr2dy * sintheta2 + + a2 * rho21 * rho21 + * (dsintheta2dy * rho2 * rho2 - drho2dy * r * rs * sintheta2 * sintheta2 + + r * rs * sintheta2 * (2. * dsintheta2dy * rho2 + sintheta2 * y)); + + Tscal dgphiphidz + = dsintheta2dz * r2 + dr2dz * sintheta2 + + (a2 * a2 * rho21 * rho21 * rs * sintheta2 * sintheta2 * z) / r + + a2 * rho21 * rho21 + * (dsintheta2dz * rho2 * rho2 - drho2dz * r * rs * sintheta2 * sintheta2 + + r * rs * sintheta2 * (2. * dsintheta2dz * rho2 + sintheta2 * z)); + + // let the fun begin again + dgcovdx = Tvec(4, 4, 0.0); + dgcovdy = Tvec(4, 4, 0.0); + dgcovdz = Tvec(4, 4, 0.0); + + // X + dgcovdx(0, 0) = r * rho21 * rho21 * rs * (-drho2dx + x); + dgcovdx(0, 1) = x2py21 * (-dgtphidx + 2. * gtphi * x * x2py21) * y; + dgcovdx(1, 1) + = -(ddeltadx * delta1 * delta1 * r2 * rho21 * x2) + + delta1 * r2 * rho21 * (2. * x + 2. * rho21 * x * x2 - drho2dx * rho21 * x2) + + x2py21 * x2py21 * (dgphiphidx - 4. * gphiphi * x * x2py21) * y2 + - r2mz21 * rho21 + * (-2. * x + 2. * r2 * r2mz21 * rho21 * x * x2 + drho2dx * rho21 * x2) * z2; + + dgcovdx(0, 2) = x2py21 * (gtphi + dgtphidx * x - 2. * gtphi * x2 * x2py21); + dgcovdx(1, 2) + = -(y + * (ddeltadx * delta1 * delta1 * r2 * rho21 * x + + delta1 * r2 * rho21 * (-1. + drho2dx * rho21 * x - 2. * rho21 * x2) + + x2py21 * x2py21 * (gphiphi + dgphiphidx * x - 4. * gphiphi * x2 * x2py21) + + r2mz21 * rho21 * (-1. + drho2dx * rho21 * x + 2. * r2 * r2mz21 * rho21 * x2) + * z2)); + + dgcovdx(2, 2) = dgphiphidx * x2 * x2py21 * x2py21 + + 2. * gphiphi * x * x2py21 * x2py21 * (1. - 2. * x2 * x2py21) + - rho21 * y2 + * (delta1 * r2 * (ddeltadx * delta1 + rho21 * (drho2dx - 2. * x)) + + r2mz21 * rho21 * (drho2dx + 2. * r2 * r2mz21 * x) * z2); + + dgcovdx(0, 3) = 0.; + dgcovdx(1, 3) + = -(z + * (r2mz21 + delta1 * term1 * (-1. + ddeltadx * delta1 * x + drho2dx * rho21 * x) + - 2. * delta1 * r2 * rho21 * rho21 * x2 + + r2mz21 + * (-2. * rho21 * rho21 * x2 + + r21 * term1 * (-1. + rho21 * x * (drho2dx + 2. * x))) + * z2 + + 2. * r2 * r2mz21 * r2mz21 * rho21 * x2 * (-1. + r21 * term1 * z2))); + + dgcovdx(2, 3) + = -(y * z + * (ddeltadx * delta1 * delta1 * term1 + + delta1 * rho21 * (drho2dx * term1 - 2. * r2 * rho21 * x) + + r2mz21 * rho21 + * (-2. * r2 * r2mz21 * x - 2. * rho21 * x * z2 + + r21 * term1 * (drho2dx + 2. * (x + r2 * r2mz21 * x)) * z2))); + + dgcovdx(3, 3) + = -(a2pr2 * a2pr2 * delta1 * drho2dx * r21 * rho21 * rho21 * z2) + - a2pr2 * a2pr2 * delta1 * r21 * rho21 * (ddeltadx * delta1 + 2. * rho21 * x) * z2 + - 2. * r2 * r2mz21 * r2mz21 * x * (-1. + r21 * term1 * z2) * (-1. + r21 * term1 * z2) + - drho2dx * r2mz21 * (-1. + r21 * term1 * z2) + * (1. + r21 * (-1. + 2. * rho2 * rho21) * term1 * z2) + + 4. * rho21 * x * z2 + * (delta1 * term1 + + r2mz21 * rho2 * (rho21 - r21 * term1) * (-1. + r21 * term1 * z2)); + + // Y + dgcovdy(0, 0) = r * rho21 * rho21 * rs * (-drho2dy + y); + dgcovdy(0, 1) = x2py21 * (-(dgtphidy * y) + gtphi * (-1. + 2. * x2py21 * y2)); + dgcovdy(1, 1) + = -(ddeltady * delta1 * delta1 * r2 * rho21 * x2) + - delta1 * r2 * rho21 * rho21 * x2 * (drho2dy - 2. * y) + + x2py21 * x2py21 + * (2. * gphiphi * y - 4. * gphiphi * x2py21 * y * y2 + dgphiphidy * y2) + - r2mz21 * rho21 * rho21 * x2 * (drho2dy + 2. * r2 * r2mz21 * y) * z2; + + dgcovdy(0, 2) = x * x2py21 * (dgtphidy - 2. * gtphi * x2py21 * y); + dgcovdy(1, 2) + = -(x + * (ddeltady * delta1 * delta1 * r2 * rho21 * y + + delta1 * r2 * rho21 * (-1. + drho2dy * rho21 * y - 2. * rho21 * y2) + + x2py21 * x2py21 * (gphiphi + dgphiphidy * y - 4. * gphiphi * x2py21 * y2) + + r2mz21 * rho21 * (-1. + drho2dy * rho21 * y + 2. * r2 * r2mz21 * rho21 * y2) + * z2)); + + dgcovdy(2, 2) + = dgphiphidy * x2 * x2py21 * x2py21 + 2. * delta1 * r2 * rho21 * y + - 4. * gphiphi * x2 * x2py21 * x2py21 * x2py21 * y + - delta1 * r2 * rho21 + * (-2. * rho21 * y * y2 + ddeltady * delta1 * y2 + drho2dy * rho21 * y2) + - r2mz21 * rho21 + * (-2. * y + 2. * r2 * r2mz21 * rho21 * y * y2 + drho2dy * rho21 * y2) * z2; + + dgcovdy(0, 3) = 0.; + dgcovdy(1, 3) + = -(x * z + * (ddeltady * delta1 * delta1 * term1 + + delta1 * rho21 * (drho2dy * term1 - 2. * r2 * rho21 * y) + + r2mz21 * rho21 + * (-2. * r2 * r2mz21 * y - 2. * rho21 * y * z2 + + r21 * term1 * (drho2dy + 2. * (y + r2 * r2mz21 * y)) * z2))); + + dgcovdy(2, 3) + = -(z + * (r2mz21 + delta1 * term1 * (-1. + ddeltady * delta1 * y + drho2dy * rho21 * y) + - 2. * delta1 * r2 * rho21 * rho21 * y2 + + r2mz21 + * (-2. * rho21 * rho21 * y2 + + r21 * term1 * (-1. + rho21 * y * (drho2dy + 2. * y))) + * z2 + + 2. * r2 * r2mz21 * r2mz21 * rho21 * y2 * (-1. + r21 * term1 * z2))); + + dgcovdy(3, 3) + = -(a2pr2 * a2pr2 * delta1 * drho2dy * r21 * rho21 * rho21 * z2) + - a2pr2 * a2pr2 * delta1 * r21 * rho21 * (ddeltady * delta1 + 2. * rho21 * y) * z2 + - 2. * r2 * r2mz21 * r2mz21 * y * (-1. + r21 * term1 * z2) * (-1. + r21 * term1 * z2) + - drho2dy * r2mz21 * (-1. + r21 * term1 * z2) + * (1. + r21 * (-1. + 2. * rho2 * rho21) * term1 * z2) + + 4. * rho21 * y * z2 + * (delta1 * term1 + + r2mz21 * rho2 * (rho21 - r21 * term1) * (-1. + r21 * term1 * z2)); + + // Z + dgcovdz(0, 0) = -(drho2dz * r * rho21 * rho21 * rs) + (rho21 * rs * term1 * z) / r; + dgcovdz(0, 1) = -(dgtphidz * x2py21 * y); + dgcovdz(1, 1) = -(ddeltadz * delta1 * delta1 * r2 * rho21 * x2) + - delta1 * drho2dz * r2 * rho21 * rho21 * x2 + + dgphiphidz * x2py21 * x2py21 * y2 + 2. * r2mz21 * rho21 * x2 * z + + 2. * delta1 * rho21 * term1 * x2 * z + - drho2dz * r2mz21 * rho21 * rho21 * x2 * z2 + - r2mz21 * r2mz21 * rho21 * x2 * (-2. * z + 2. * term1 * z) * z2; + + dgcovdz(0, 2) = dgtphidz * x * x2py21; + dgcovdz(1, 2) = -(ddeltadz * delta1 * delta1 * r2 * rho21 * x * y) + - delta1 * drho2dz * r2 * rho21 * rho21 * x * y + - dgphiphidz * x * x2py21 * x2py21 * y + 2. * r2mz21 * rho21 * x * y * z + + 2. * delta1 * rho21 * term1 * x * y * z + - drho2dz * r2mz21 * rho21 * rho21 * x * y * z2 + - r2mz21 * r2mz21 * rho21 * x * y * (-2. * z + 2. * term1 * z) * z2; + + dgcovdz(2, 2) + = dgphiphidz * x2 * x2py21 * x2py21 - ddeltadz * delta1 * delta1 * r2 * rho21 * y2 + - delta1 * drho2dz * r2 * rho21 * rho21 * y2 + 2. * r2mz21 * rho21 * y2 * z + + 2. * delta1 * rho21 * term1 * y2 * z - drho2dz * r2mz21 * rho21 * rho21 * y2 * z2 + - r2mz21 * r2mz21 * rho21 * y2 * (-2. * z + 2. * term1 * z) * z2; + + dgcovdz(0, 3) = 0.; + dgcovdz(1, 3) + = delta1 * term1 * x - ddeltadz * delta1 * delta1 * term1 * x * z + - delta1 * drho2dz * rho21 * term1 * x * z + 2. * delta1 * rho21 * term1 * x * z2 + - r2mz21 * x * (1. - r21 * term1 * z2) + + r2mz21 * r2mz21 * x * z * (-2. * z + 2. * term1 * z) * (1. - r21 * term1 * z2) + - r2mz21 * x * z + * (term2 - 2. * r21 * term1 * z - 2. * r21 * rho21 * term1 * z * z2 + + drho2dz * r21 * rho21 * term1 * z2); + + dgcovdz(2, 3) + = delta1 * term1 * y - ddeltadz * delta1 * delta1 * term1 * y * z + - delta1 * drho2dz * rho21 * term1 * y * z + 2. * delta1 * rho21 * term1 * y * z2 + - r2mz21 * y * (1. - r21 * term1 * z2) + + r2mz21 * r2mz21 * y * z * (-2. * z + 2. * term1 * z) * (1. - r21 * term1 * z2) + - r2mz21 * y * z + * (term2 - 2. * r21 * term1 * z - 2. * r21 * rho21 * term1 * z * z2 + + drho2dz * r21 * rho21 * term1 * z2); + + dgcovdz(3, 3) = 2. * a2pr2 * a2pr2 * delta1 * r21 * rho21 * z + - (2. * a2pr2 * a2pr2 * a2pr2 * delta1 * rho21 * rho21 * z * z2) / (r2 * r2) + + 4. * a2pr2 * a2pr2 * delta1 * r21 * rho21 * rho21 * z * z2 + - a2pr2 * a2pr2 * ddeltadz * delta1 * delta1 * r21 * rho21 * z2 + - a2pr2 * a2pr2 * delta1 * drho2dz * r21 * rho21 * rho21 * z2 + + drho2dz * r2mz21 * (1. - r21 * term1 * z2) * (1. - r21 * term1 * z2) + - r2mz21 * r2mz21 * rho2 * (-2. * z + 2. * term1 * z) + * (1. - r21 * term1 * z2) * (1. - r21 * term1 * z2) + + 2. * r2mz21 * rho2 * (1. - r21 * term1 * z2) + * (term2 - 2. * r21 * term1 * z - 2. * r21 * rho21 * term1 * z * z2 + + drho2dz * r21 * rho21 * term1 * z2); + + // Hopefully that was alright + + dgcovdx(1, 0) = dgcovdx(0, 1); + dgcovdx(2, 0) = dgcovdx(0, 2); + dgcovdx(2, 1) = dgcovdx(1, 2); + dgcovdx(3, 0) = dgcovdx(0, 3); + dgcovdx(3, 1) = dgcovdx(1, 3); + dgcovdx(3, 2) = dgcovdx(2, 3); + + dgcovdy(1, 0) = dgcovdy(0, 1); + dgcovdy(2, 0) = dgcovdy(0, 2); + dgcovdy(2, 1) = dgcovdy(1, 2); + dgcovdy(3, 0) = dgcovdy(0, 3); + dgcovdy(3, 1) = dgcovdy(1, 3); + dgcovdy(3, 2) = dgcovdy(2, 3); + + dgcovdz(1, 0) = dgcovdz(0, 1); + dgcovdz(2, 0) = dgcovdz(0, 2); + dgcovdz(2, 1) = dgcovdz(1, 2); + dgcovdz(3, 0) = dgcovdz(0, 3); + dgcovdz(3, 1) = dgcovdz(1, 3); + dgcovdz(3, 2) = dgcovdz(2, 3); + } + +} // namespace shamphys From b09bf765e8bc35b4ed1d525fe625a74dee49abe5 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Tue, 13 Jan 2026 22:30:52 +0900 Subject: [PATCH 02/18] authors --- src/shamphys/include/shamphys/metrics.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp index 0e8dbffa92..7a0ebb8190 100644 --- a/src/shamphys/include/shamphys/metrics.hpp +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -11,6 +11,7 @@ /** * @file metrics.hpp + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) * @brief */ From 872650ed471a5ce016870411a1002adf53b752d1 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 14 Jan 2026 00:15:31 +0900 Subject: [PATCH 03/18] add a tag dispatch --- src/shamphys/include/shamphys/metrics.hpp | 39 ++++++++++++++++++++--- 1 file changed, 34 insertions(+), 5 deletions(-) diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp index 7a0ebb8190..78307b847a 100644 --- a/src/shamphys/include/shamphys/metrics.hpp +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -20,9 +20,13 @@ namespace shamphys { - // The Boyer-Lindquist metric tensor in CARTESIAN-like form + struct Kerr {}; + struct Schwarzschild {}; + + /******************** The Kerr metric ***********************/ + // ################in CARTESIAN-like form ################## template - inline Tvec get_cartesian_cov_metric(Tvec pos) { + inline Tvec get_cartesian_cov_metric_impl(Kerr, Tvec pos) { // metric params Tscal a = 1; @@ -83,7 +87,7 @@ namespace shamphys { } template - inline Tvec get_cartesian_contrav_metric(Tvec pos) { + inline Tvec get_cartesian_contrav_metric_impl(Kerr, Tvec pos) { // metric params Tscal a = 1; @@ -146,8 +150,8 @@ namespace shamphys { } template - inline void metric_cartesian_derivatives( - Tvec pos, Tvec &dgcovdx, Tvec &dgcovdy, Tvec &dgcovdz) { + inline void metric_cartesian_derivatives_impl( + Kerr, Tvec pos, Tvec &dgcovdx, Tvec &dgcovdy, Tvec &dgcovdz) { // metric params Tscal a = 1; @@ -430,4 +434,29 @@ namespace shamphys { dgcovdz(3, 2) = dgcovdz(2, 3); } + // ################in spherical form ################## + // TODO + + /******************** The Kerr metric ***********************/ + // ################in CARTESIAN-like form ################## + // ################in spherical form ################## + // TODO + + /******************** frontend ***********************/ + template + inline Tvec get_cartesian_cov_metric(Tvec pos) { + return get_cartesian_cov_metric_impl(MetricTag{}, pos); + } + + template + inline Tvec get_cartesian_contrav_metric(Tvec pos) { + return get_cartesian_contrav_metric_impl(MetricTag{}, pos); + } + + template + inline void metric_cartesian_derivatives( + Tvec pos, Tvec &dgcovdx, Tvec &dgcovdy, Tvec &dgcovdz) { + metric_cartesian_derivatives_impl(MetricTag{}, pos, dgcovdx, dgcovdy, dgcovdz); + } + } // namespace shamphys From ab2cc36347952a85c37b0e4d534f339448063084 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 14 Jan 2026 00:36:56 +0900 Subject: [PATCH 04/18] add metric params in struct --- src/shamphys/include/shamphys/metrics.hpp | 37 +++++++++++++---------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp index 78307b847a..7677958e48 100644 --- a/src/shamphys/include/shamphys/metrics.hpp +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -20,21 +20,29 @@ namespace shamphys { - struct Kerr {}; + template + struct Kerr { + + Tscal a; + Tscal bh_mass; + Tscal rs; + + Kerr(Tscal spin = 1, Tscal mass = 1) : a(spin), bh_mass(bh_mass), rs(2. * bh_mass) {} + }; + struct Schwarzschild {}; /******************** The Kerr metric ***********************/ // ################in CARTESIAN-like form ################## template - inline Tvec get_cartesian_cov_metric_impl(Kerr, Tvec pos) { + inline Tvec get_cartesian_cov_metric_impl(const Kerr &kerr, Tvec pos) { // metric params - Tscal a = 1; - Tscal mass1 = 1.; - Tscal rs = 2. * mass1; + static constexpr Tscal a = kerr.a; + static constexpr Tscal rs = kerr.rs; // pos coords - Tscal x, y, z = pos; + Tscal x = pos(0), y = pos(1), z = pos(2); Tscal x2 = x * x; Tscal y2 = y * y; Tscal z2 = z * z; @@ -87,15 +95,14 @@ namespace shamphys { } template - inline Tvec get_cartesian_contrav_metric_impl(Kerr, Tvec pos) { + inline Tvec get_cartesian_contrav_metric_impl(const Kerr &kerr, Tvec pos) { // metric params - Tscal a = 1; - Tscal mass1 = 1.; - Tscal rs = 2. * mass1; + static constexpr Tscal a = kerr.a; + static constexpr Tscal rs = kerr.rs; // pos coords - Tscal x, y, z = pos; + Tscal x = pos(0), y = pos(1), z = pos(2); Tscal x2 = x * x; Tscal y2 = y * y; Tscal z2 = z * z; @@ -151,13 +158,11 @@ namespace shamphys { template inline void metric_cartesian_derivatives_impl( - Kerr, Tvec pos, Tvec &dgcovdx, Tvec &dgcovdy, Tvec &dgcovdz) { + const Kerr &kerr, Tvec pos, Tvec &dgcovdx, Tvec &dgcovdy, Tvec &dgcovdz) { // metric params - Tscal a = 1; - Tscal mass1 = 1.; - Tscal rs = 2. * mass1; - + static constexpr Tscal a = kerr.a; + static constexpr Tscal rs = kerr.rs; // position coords Tscal x = pos(0); Tscal y = pos(1); From 80e26b8503e53e79ed19ff72eb4343e14eab45b1 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Wed, 14 Jan 2026 06:57:07 +0000 Subject: [PATCH 05/18] [gh-action] trigger CI with empty commit From a6c2863d5bd6e3e995aebd029a34cf22be4cd4eb Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Sat, 24 Jan 2026 19:36:46 +0900 Subject: [PATCH 06/18] move to mdspan type for the metric --- src/shamphys/include/shamphys/metrics.hpp | 114 +++++++++++++++++----- 1 file changed, 88 insertions(+), 26 deletions(-) diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp index 7677958e48..7bb3a95a8e 100644 --- a/src/shamphys/include/shamphys/metrics.hpp +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -1,7 +1,7 @@ // -------------------------------------------------------// // // SHAMROCK code for hydrodynamics -// Copyright (c) 2021-2025 Timothée David--Cléris +// Copyright (c) 2021-2026 Timothée David--Cléris // SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 // Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information // @@ -15,7 +15,10 @@ * @brief */ +#include "shambase/assert.hpp" +#include "shambackends/math.hpp" #include "shamunits/Constants.hpp" +#include #include namespace shamphys { @@ -34,8 +37,17 @@ namespace shamphys { /******************** The Kerr metric ***********************/ // ################in CARTESIAN-like form ################## - template - inline Tvec get_cartesian_cov_metric_impl(const Kerr &kerr, Tvec pos) { + template< + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> + inline void get_cartesian_covariant_metric_impl( + const Kerr &kerr, + const std::mdspan, Layout1, Accessor1> &pos, + std::mdspan, Layout2, Accessor2> &g) { // metric params static constexpr Tscal a = kerr.a; @@ -68,7 +80,6 @@ namespace shamphys { Tscal gtphi_on_x2py2 = -rs * r * a / (rho2 * a2pr2); // let the fun begin - Tvec g(4, 4); g(0, 0) = gtt; g(1, 0) = -y * gtphi_on_x2py2; g(2, 0) = x * gtphi_on_x2py2; @@ -90,12 +101,19 @@ namespace shamphys { g(1, 3) = g(3, 1); g(2, 3) = g(3, 2); g(3, 3) = 1. + ((a2 + r2) * rs * z2) / (delta * r * rho2); - - return g; } - template - inline Tvec get_cartesian_contrav_metric_impl(const Kerr &kerr, Tvec pos) { + template< + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> + inline void get_cartesian_contravariant_metric_impl( + const Kerr &kerr, + const std::mdspan, Layout1, Accessor1> &pos, + std::mdspan, Layout2, Accessor2> &g) { // metric params static constexpr Tscal a = kerr.a; @@ -128,7 +146,6 @@ namespace shamphys { Tscal gtphi_on_x2py2 = -rs * r * a / (rho2 * a2pr2); // let the fun begin - Tvec g(4, 4); Tscal domegaterm = 1. / (omega * gtphi + gtt); g(0, 0) = domegaterm; g(1, 0) = -y * omega * domegaterm; @@ -152,13 +169,21 @@ namespace shamphys { g(1, 3) = g(3, 1); g(2, 3) = g(3, 2); g(3, 3) = (r2 - z2) * inv_rho2 + (delta * z2) / (r2 * rho2); - - return g; } - template + template< + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> inline void metric_cartesian_derivatives_impl( - const Kerr &kerr, Tvec pos, Tvec &dgcovdx, Tvec &dgcovdy, Tvec &dgcovdz) { + const Kerr &kerr, + const std::mdspan, Layout1, Accessor1> &pos, + std::mdspan, Layout2, Accessor2> &dgcovdx, + std::mdspan, Layout2, Accessor2> &dgcovdy, + std::mdspan, Layout2, Accessor2> &dgcovdz) { // metric params static constexpr Tscal a = kerr.a; @@ -246,10 +271,6 @@ namespace shamphys { + r * rs * sintheta2 * (2. * dsintheta2dz * rho2 + sintheta2 * z)); // let the fun begin again - dgcovdx = Tvec(4, 4, 0.0); - dgcovdy = Tvec(4, 4, 0.0); - dgcovdz = Tvec(4, 4, 0.0); - // X dgcovdx(0, 0) = r * rho21 * rho21 * rs * (-drho2dx + x); dgcovdx(0, 1) = x2py21 * (-dgtphidx + 2. * gtphi * x * x2py21) * y; @@ -448,20 +469,61 @@ namespace shamphys { // TODO /******************** frontend ***********************/ - template - inline Tvec get_cartesian_cov_metric(Tvec pos) { - return get_cartesian_cov_metric_impl(MetricTag{}, pos); + template< + class MetricTag, + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> + inline void get_cartesian_covariant_metric( + const std::mdspan, Layout1, Accessor1> pos, + std::mdspan, Layout2, Accessor2> &g) { + return get_cartesian_covariant_metric_impl< + Tscal, + SizeType, + Layout1, + Layout2, + Accessor1, + Accessor2>(MetricTag{}, pos, g); } - template - inline Tvec get_cartesian_contrav_metric(Tvec pos) { - return get_cartesian_contrav_metric_impl(MetricTag{}, pos); + template< + class MetricTag, + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> + inline void get_cartesian_contravariant_metric( + const std::mdspan, Layout1, Accessor1> pos, + std::mdspan, Layout2, Accessor2> &g) { + return get_cartesian_contravariant_metric_impl< + Tscal, + SizeType, + Layout1, + Layout2, + Accessor1, + Accessor2>(MetricTag{}, pos, g); } - template + template< + class MetricTag, + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> inline void metric_cartesian_derivatives( - Tvec pos, Tvec &dgcovdx, Tvec &dgcovdy, Tvec &dgcovdz) { - metric_cartesian_derivatives_impl(MetricTag{}, pos, dgcovdx, dgcovdy, dgcovdz); + const std::mdspan, Layout1, Accessor1> pos, + std::mdspan, Layout2, Accessor2> &dgcovdx, + std::mdspan, Layout2, Accessor2> &dgcovdy, + std::mdspan, Layout2, Accessor2> &dgcovdz) { + metric_cartesian_derivatives_impl( + MetricTag{}, pos, dgcovdx, dgcovdy, dgcovdz); } } // namespace shamphys From 777133d8cb857555feabf6269f5cabcf0a5e4af7 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 18 Feb 2026 16:28:56 +0900 Subject: [PATCH 07/18] add GR utils --- src/shamphys/include/shamphys/GRUtils.hpp | 57 +++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 src/shamphys/include/shamphys/GRUtils.hpp diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp new file mode 100644 index 0000000000..58197c5880 --- /dev/null +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -0,0 +1,57 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +#pragma once + +/** + * @file GRUtils.hpp + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + */ + +#include "shambase/assert.hpp" +#include "shambackends/math.hpp" +#include "shamunits/Constants.hpp" +#include +#include + +namespace shamphys { + + template + struct GR_physics { + + inline static constexpr Tscal GR_dot(Tvec a, Tvec b, std::mdspan, Layout2, Accessor2> gcov) + { return 0;} + + inline static constexpr Tscal get_U0(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) + { return 0; } + + inline static constexpr Tvec get_V(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) + { return 0; } + + inline static constexpr Tscal get_sqrt_g(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) + { return 0; } + + inline static constexpr Tscal get_sqrt_gamma(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) + { return 0; } + + }; + + + + + + + + +} // namespace shamphys From 7c6286642a7dd06725f409da026af5c4dd025de2 Mon Sep 17 00:00:00 2001 From: "autofix-ci[bot]" <114827586+autofix-ci[bot]@users.noreply.github.com> Date: Wed, 18 Feb 2026 07:29:29 +0000 Subject: [PATCH 08/18] [autofix.ci] automatic fix: pre-commit hooks --- src/shamphys/include/shamphys/GRUtils.hpp | 57 +++++++++++++---------- 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index 58197c5880..b608acf955 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -23,35 +23,42 @@ namespace shamphys { - template - struct GR_physics { - - inline static constexpr Tscal GR_dot(Tvec a, Tvec b, std::mdspan, Layout2, Accessor2> gcov) - { return 0;} - - inline static constexpr Tscal get_U0(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) - { return 0; } - - inline static constexpr Tvec get_V(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) - { return 0; } - - inline static constexpr Tscal get_sqrt_g(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) - { return 0; } - - inline static constexpr Tscal get_sqrt_gamma(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) - { return 0; } - + struct GR_physics { + + inline static constexpr Tscal GR_dot( + Tvec a, + Tvec b, + std::mdspan, Layout2, Accessor2> gcov) { + return 0; + } + + inline static constexpr Tscal get_U0( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + return 0; + } + + inline static constexpr Tvec get_V( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + return 0; + } + + inline static constexpr Tscal get_sqrt_g( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + return 0; + } + + inline static constexpr Tscal get_sqrt_gamma( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + return 0; + } }; - - - - - - - } // namespace shamphys From d8087e2d73d54e8facf5c5fe2c18e75a1d8802da Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 18 Feb 2026 17:32:39 +0900 Subject: [PATCH 09/18] implemented utils --- src/shamphys/include/shamphys/GRUtils.hpp | 85 ++++++++++++++++------- 1 file changed, 60 insertions(+), 25 deletions(-) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index 58197c5880..c318779f45 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -15,6 +15,7 @@ * @brief */ +#include "shambase/aliases_int.hpp" #include "shambase/assert.hpp" #include "shambackends/math.hpp" #include "shamunits/Constants.hpp" @@ -23,35 +24,69 @@ namespace shamphys { - template - struct GR_physics { - - inline static constexpr Tscal GR_dot(Tvec a, Tvec b, std::mdspan, Layout2, Accessor2> gcov) - { return 0;} - - inline static constexpr Tscal get_U0(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) - { return 0; } - - inline static constexpr Tvec get_V(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) - { return 0; } - - inline static constexpr Tscal get_sqrt_g(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) - { return 0; } - - inline static constexpr Tscal get_sqrt_gamma(Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) - { return 0; } - + struct GR_physics { + + inline static constexpr Tscal get_alpha( + std::mdspan, Layout2, Accessor2> gcov) { + // alpha = 1 / sqrt(-g00) + return 1. / sycl::sqrt(-gcov(0, 0)); + } + + inline static constexpr Tvec get_betaUP( + std::mdspan, Layout2, Accessor2> gcov) { + Tscal alpha = get_alpha(gcov); + Tscal alpha2 = alpha * alpha; + Tvec betaUP = {0., 0., 0.}; + + betaUP[0] = gcov(0, 1) * alpha2; + betaUP[1] = gcov(0, 1) * alpha2; + betaUP[2] = gcov(0, 1) * alpha2; + + return betaUP; + } + + inline static constexpr Tscal GR_dot( + Tvec a, + Tvec b, + std::mdspan, Layout2, Accessor2> gcov) { + return 0; + } + + inline static constexpr Tscal get_U0( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + return 0; + } + + inline static constexpr Tvec get_V( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + Tscal alpha = get_alpha(gcov); + Tvec betaUP = get_betaUP(gcov); + Tvec V = (vxyz + betaUP) / alpha; + + return V; + } + + inline static constexpr Tscal get_sqrt_g( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + // for unusual metric, need to compute the determinant + // for now (Kerr), this is enough + return 1; + } + + inline static constexpr Tscal get_sqrt_gamma( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + // for unusual metric, need to compute the determinant of the spatial metric + // for now (Kerr), this is enough + return 1; + } }; - - - - - - - } // namespace shamphys From a5cda8580facf0058b828d5ce9524f9125751676 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 18 Feb 2026 17:48:35 +0900 Subject: [PATCH 10/18] get ADM components --- src/shamphys/include/shamphys/GRUtils.hpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index e1759e1394..857439f24e 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -24,10 +24,6 @@ namespace shamphys { - template< - class Tvec, - class Tscal, - class SizeType, template< class Tvec, class Tscal, @@ -62,14 +58,17 @@ namespace shamphys { Tvec b, std::mdspan, Layout2, Accessor2> gcov) { return 0; + // TODO } inline static constexpr Tscal get_U0( Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { - return 0; - } - inline static constexpr Tscal get_U0( - Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + + // U0 = gamma / alpha + Tscal alpha = get_alpha(gcov); + + // Tscal gamma = mat_inv_33 ... submdspan would be nice ... it's in the 26 standard + // though ...; return 0; } From 0af5a279e8e0770f79cee9451fd1fad8cc287fa9 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 18 Feb 2026 18:11:22 +0900 Subject: [PATCH 11/18] fix typos and bad calls --- src/shamphys/include/shamphys/metrics.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp index 7bb3a95a8e..4623290878 100644 --- a/src/shamphys/include/shamphys/metrics.hpp +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -30,7 +30,7 @@ namespace shamphys { Tscal bh_mass; Tscal rs; - Kerr(Tscal spin = 1, Tscal mass = 1) : a(spin), bh_mass(bh_mass), rs(2. * bh_mass) {} + Kerr(Tscal spin = 1, Tscal mass = 1) : a(spin), bh_mass(mass), rs(2. * mass) {} }; struct Schwarzschild {}; @@ -50,8 +50,8 @@ namespace shamphys { std::mdspan, Layout2, Accessor2> &g) { // metric params - static constexpr Tscal a = kerr.a; - static constexpr Tscal rs = kerr.rs; + const Tscal a = kerr.a; + const Tscal rs = kerr.rs; // pos coords Tscal x = pos(0), y = pos(1), z = pos(2); @@ -116,8 +116,8 @@ namespace shamphys { std::mdspan, Layout2, Accessor2> &g) { // metric params - static constexpr Tscal a = kerr.a; - static constexpr Tscal rs = kerr.rs; + const Tscal a = kerr.a; + const Tscal rs = kerr.rs; // pos coords Tscal x = pos(0), y = pos(1), z = pos(2); @@ -186,8 +186,8 @@ namespace shamphys { std::mdspan, Layout2, Accessor2> &dgcovdz) { // metric params - static constexpr Tscal a = kerr.a; - static constexpr Tscal rs = kerr.rs; + const Tscal a = kerr.a; + const Tscal rs = kerr.rs; // position coords Tscal x = pos(0); Tscal y = pos(1); @@ -480,7 +480,7 @@ namespace shamphys { inline void get_cartesian_covariant_metric( const std::mdspan, Layout1, Accessor1> pos, std::mdspan, Layout2, Accessor2> &g) { - return get_cartesian_covariant_metric_impl< + get_cartesian_covariant_metric_impl< Tscal, SizeType, Layout1, @@ -500,7 +500,7 @@ namespace shamphys { inline void get_cartesian_contravariant_metric( const std::mdspan, Layout1, Accessor1> pos, std::mdspan, Layout2, Accessor2> &g) { - return get_cartesian_contravariant_metric_impl< + get_cartesian_contravariant_metric_impl< Tscal, SizeType, Layout1, From d794403f7d23c029dcb4ab1c7bb63bc2c982dcfa Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Thu, 19 Feb 2026 11:48:17 +0900 Subject: [PATCH 12/18] fix typo --- src/shamphys/include/shamphys/GRUtils.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index 857439f24e..a9d631894c 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -47,8 +47,8 @@ namespace shamphys { Tvec betaUP = {0., 0., 0.}; betaUP[0] = gcov(0, 1) * alpha2; - betaUP[1] = gcov(0, 1) * alpha2; - betaUP[2] = gcov(0, 1) * alpha2; + betaUP[1] = gcov(0, 2) * alpha2; + betaUP[2] = gcov(0, 3) * alpha2; return betaUP; } From 6b4825e7f7431ad7edc1c232a27f472640775d3b Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 11 Mar 2026 09:31:58 +0900 Subject: [PATCH 13/18] incorporate reviews --- src/shammath/include/shammath/matrix_op.hpp | 31 +++++++++++++++++++++ src/shamphys/include/shamphys/GRUtils.hpp | 21 +++++++++++--- src/shamphys/include/shamphys/metrics.hpp | 10 +------ 3 files changed, 49 insertions(+), 13 deletions(-) diff --git a/src/shammath/include/shammath/matrix_op.hpp b/src/shammath/include/shammath/matrix_op.hpp index dde113f365..4a96829f42 100644 --- a/src/shammath/include/shammath/matrix_op.hpp +++ b/src/shammath/include/shammath/matrix_op.hpp @@ -334,6 +334,37 @@ namespace shammath { } } + /** + * @brief Compute the determinant of a 3x3 matrix. + * + * @param input The input matrix to invert. + * + * @details This function computes the determinant of a 3x3 matrix and returns it. + * + */ + template + inline T mat_det_33( + const std::mdspan, Layout, Accessor> &input) { + + T &a00 = input(0, 0); + T &a10 = input(1, 0); + T &a20 = input(2, 0); + + T &a01 = input(0, 1); + T &a11 = input(1, 1); + T &a21 = input(2, 1); + + T &a02 = input(0, 2); + T &a12 = input(1, 2); + T &a22 = input(2, 2); + + T det + = (-a02 * a11 * a20 + a01 * a12 * a20 + a02 * a10 * a21 - a00 * a12 * a21 + - a01 * a10 * a22 + a00 * a11 * a22); + + return det; + } + /** * @brief Compute the inverse of a 3x3 matrix. * diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index a9d631894c..25a9c43d4e 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -61,15 +61,26 @@ namespace shamphys { // TODO } + inline static constexpr Tscal get_gamma( + std::mdspan, Layout2, Accessor2> gcov) { + + std::mdspan, Layout2, Accessor2> gamma_ij; + for (int i = 1; i < 4; i++) { + for (int j = 1; j < 4; j++) { + gamma_ij(i, j) = gcov(i, j); + } + } + + return mat_det_33(gamma_ij); + } + inline static constexpr Tscal get_U0( Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { // U0 = gamma / alpha Tscal alpha = get_alpha(gcov); - - // Tscal gamma = mat_inv_33 ... submdspan would be nice ... it's in the 26 standard - // though ...; - return 0; + Tscal gamma = get_gamma(gcov); + return gamma / alpha; } inline static constexpr Tvec get_V( @@ -81,6 +92,7 @@ namespace shamphys { return V; } + /// placeholder function for more complex metrics inline static constexpr Tscal get_sqrt_g( Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { // for unusual metric, need to compute the determinant @@ -88,6 +100,7 @@ namespace shamphys { return 1; } + /// placeholder function for more complex metrics inline static constexpr Tscal get_sqrt_gamma( Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { // for unusual metric, need to compute the determinant of the spatial metric diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp index 4623290878..9e399df834 100644 --- a/src/shamphys/include/shamphys/metrics.hpp +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -30,7 +30,7 @@ namespace shamphys { Tscal bh_mass; Tscal rs; - Kerr(Tscal spin = 1, Tscal mass = 1) : a(spin), bh_mass(mass), rs(2. * mass) {} + Kerr(Tscal spin, Tscal mass) : a(spin), bh_mass(mass), rs(2. * mass) {} }; struct Schwarzschild {}; @@ -460,14 +460,6 @@ namespace shamphys { dgcovdz(3, 2) = dgcovdz(2, 3); } - // ################in spherical form ################## - // TODO - - /******************** The Kerr metric ***********************/ - // ################in CARTESIAN-like form ################## - // ################in spherical form ################## - // TODO - /******************** frontend ***********************/ template< class MetricTag, From b641a1659269e5de3937f7c65c1d792cf6353c9d Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 11 Mar 2026 15:31:56 +0900 Subject: [PATCH 14/18] add beta, gamma_ij --- src/shamphys/include/shamphys/GRUtils.hpp | 37 +++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index 25a9c43d4e..01652eae0c 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -53,6 +53,19 @@ namespace shamphys { return betaUP; } + inline static constexpr Tvec get_betaDOWN( + std::mdspan, Layout2, Accessor2> gcov) { + Tscal alpha = get_alpha(gcov); + Tscal alpha2 = alpha * alpha; + Tvec betaDOWN = {0., 0., 0.}; + + betaDOWN[0] = gcov(0, 1); + betaDOWN[1] = gcov(0, 2); + betaDOWN[2] = gcov(0, 3); + + return betaDOWN; + } + inline static constexpr Tscal GR_dot( Tvec a, Tvec b, @@ -61,6 +74,30 @@ namespace shamphys { // TODO } + inline static constexpr Tscal get_gammaijUP( + std::mdspan, Layout2, Accessor2> gcov) { + + std::mdspan, Layout2, Accessor2> gamma_ij; + for (int i = 1; i < 4; i++) { + for (int j = 1; j < 4; j++) { + gamma_ij(i, j) = gcov(i, j); + } + } + } + + inline static constexpr Tscal get_gammaijDOWN( + std::mdspan, Layout2, Accessor2> gcov) { + + Tscal alpha = get_alpha(gcov); + Tvec betaUP = get_betaUP(gcov); + std::mdspan, Layout2, Accessor2> gamma_ij; + for (int i = 1; i < 4; i++) { + for (int j = 1; j < 4; j++) { + gamma_ij(i, j) = gcov(i, j) + betaUP[i] * betaUP[j] / (-gcov(0, 0)); + } + } + } + inline static constexpr Tscal get_gamma( std::mdspan, Layout2, Accessor2> gcov) { From c264e22009da9c5451a09e9e2e88ee043341331a Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 11 Mar 2026 16:18:16 +0900 Subject: [PATCH 15/18] fix typos and wrong indexes --- src/shamphys/include/shamphys/GRUtils.hpp | 66 ++++++++++++++--------- src/shamphys/include/shamphys/metrics.hpp | 8 +-- 2 files changed, 45 insertions(+), 29 deletions(-) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index 01652eae0c..9b9cdafad9 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -35,28 +35,26 @@ namespace shamphys { struct GR_physics { inline static constexpr Tscal get_alpha( - std::mdspan, Layout2, Accessor2> gcov) { + std::mdspan, Layout2, Accessor2> gcon) { // alpha = 1 / sqrt(-g00) - return 1. / sycl::sqrt(-gcov(0, 0)); + return 1. / sycl::sqrt(-gcon(0, 0)); } inline static constexpr Tvec get_betaUP( - std::mdspan, Layout2, Accessor2> gcov) { - Tscal alpha = get_alpha(gcov); + std::mdspan, Layout2, Accessor2> gcon) { + Tscal alpha = get_alpha(gcon); Tscal alpha2 = alpha * alpha; Tvec betaUP = {0., 0., 0.}; - betaUP[0] = gcov(0, 1) * alpha2; - betaUP[1] = gcov(0, 2) * alpha2; - betaUP[2] = gcov(0, 3) * alpha2; + betaUP[0] = gcon(0, 1) * alpha2; + betaUP[1] = gcon(0, 2) * alpha2; + betaUP[2] = gcon(0, 3) * alpha2; return betaUP; } inline static constexpr Tvec get_betaDOWN( std::mdspan, Layout2, Accessor2> gcov) { - Tscal alpha = get_alpha(gcov); - Tscal alpha2 = alpha * alpha; Tvec betaDOWN = {0., 0., 0.}; betaDOWN[0] = gcov(0, 1); @@ -66,36 +64,48 @@ namespace shamphys { return betaDOWN; } + /// placeholder inline static constexpr Tscal GR_dot( Tvec a, Tvec b, std::mdspan, Layout2, Accessor2> gcov) { return 0; - // TODO + Tscal result = 0; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + result += gcov(i+1, j+1) * a[i] * b[j]; + } + } + return result; } - inline static constexpr Tscal get_gammaijUP( + inline static constexpr std::mdspan, Layout2, Accessor2> get_gammaijDOWN( std::mdspan, Layout2, Accessor2> gcov) { std::mdspan, Layout2, Accessor2> gamma_ij; for (int i = 1; i < 4; i++) { for (int j = 1; j < 4; j++) { - gamma_ij(i, j) = gcov(i, j); + gamma_ij(i-1, j-1) = gcov(i, j); } } + + return gamma_ij; } - inline static constexpr Tscal get_gammaijDOWN( - std::mdspan, Layout2, Accessor2> gcov) { + inline static constexpr std::mdspan, Layout2, Accessor2> get_gammaijUP( + std::mdspan, Layout2, Accessor2> gcon) { - Tscal alpha = get_alpha(gcov); - Tvec betaUP = get_betaUP(gcov); + Tscal alpha = get_alpha(gcon); + Tvec betaUP = get_betaUP(gcon); std::mdspan, Layout2, Accessor2> gamma_ij; for (int i = 1; i < 4; i++) { for (int j = 1; j < 4; j++) { - gamma_ij(i, j) = gcov(i, j) + betaUP[i] * betaUP[j] / (-gcov(0, 0)); + gamma_ij(i-1, j-1) = gcon(i, j) + betaUP[i-1] * betaUP[j-1] / (-gcon(0, 0)); } } + + return gamma_ij; + } inline static constexpr Tscal get_gamma( @@ -111,19 +121,25 @@ namespace shamphys { return mat_det_33(gamma_ij); } + inline static constexpr Tscal get_lorentz_factor( + Tvec momentum, Tscal enthalpy) { + + return sycl::sqrt(1. + GR_dot(momentum, momentum) / (enthalpy * enthalpy)); + } + inline static constexpr Tscal get_U0( - Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + Tvec momentum, Tscal enthalpy, std::mdspan, Layout2, Accessor2> gcon) { - // U0 = gamma / alpha - Tscal alpha = get_alpha(gcov); - Tscal gamma = get_gamma(gcov); - return gamma / alpha; + // U0 = gamma / alpha wrong + Tscal alpha = get_alpha(gcon); + Tscal lorentz_factor = get_lorentz_factor(momentum, enthalpy); + return lorentz_factor / alpha; } inline static constexpr Tvec get_V( - Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { - Tscal alpha = get_alpha(gcov); - Tvec betaUP = get_betaUP(gcov); + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcon) { + Tscal alpha = get_alpha(gcon); + Tvec betaUP = get_betaUP(gcon); Tvec V = (vxyz + betaUP) / alpha; return V; diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp index 9e399df834..604ed61ac1 100644 --- a/src/shamphys/include/shamphys/metrics.hpp +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -61,7 +61,7 @@ namespace shamphys { Tscal a2 = a * a; Tscal r2_spher = x2 + y2 + z2; Tscal r2 = 0.5 * (r2_spher - a2) - + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); + + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); Tscal r = sycl::sqrt(r2); Tscal inv_r2 = 1. / r2; Tscal rho2 = r2 + a2 * (z2 * inv_r2); @@ -127,7 +127,7 @@ namespace shamphys { Tscal a2 = a * a; Tscal r2_spher = x2 + y2 + z2; Tscal r2 = 0.5 * (r2_spher - a2) - + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); + + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); Tscal r = sycl::sqrt(r2); Tscal inv_r2 = 1. / r2; Tscal rho2 = r2 + a2 * (z2 * inv_r2); @@ -200,8 +200,8 @@ namespace shamphys { Tscal r2_spher = x2 + y2 + z2; Tscal r2 = 0.5 * (r2_spher - a2) - + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); - Tscal r = sycl::sqrt(r2); + + 0.5 * sycl::sqrt((r2_spher - a2) * (r2_spher - a2) + 4. * a2 * z2); + Tscal r = sycl::sqrt(r2); Tscal rho2 = r2 + a2 * (z2 / r2); Tscal a2pr2 = a2 + r2; From 51ad55fa4285a1af43045b62a04f3cb9a60b5bca Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Wed, 11 Mar 2026 07:23:29 +0000 Subject: [PATCH 16/18] [gh-action] trigger CI with empty commit From 35f0c621ef0b0d4869efad9e471cb50d6210d9c0 Mon Sep 17 00:00:00 2001 From: "autofix-ci[bot]" <114827586+autofix-ci[bot]@users.noreply.github.com> Date: Wed, 11 Mar 2026 07:24:04 +0000 Subject: [PATCH 17/18] [autofix.ci] automatic fix: pre-commit hooks --- src/shamphys/include/shamphys/GRUtils.hpp | 25 ++++++++++++----------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index 9b9cdafad9..c45db13249 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -73,39 +73,39 @@ namespace shamphys { Tscal result = 0; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { - result += gcov(i+1, j+1) * a[i] * b[j]; + result += gcov(i + 1, j + 1) * a[i] * b[j]; } } return result; } - inline static constexpr std::mdspan, Layout2, Accessor2> get_gammaijDOWN( - std::mdspan, Layout2, Accessor2> gcov) { + inline static constexpr std::mdspan, Layout2, Accessor2> + get_gammaijDOWN(std::mdspan, Layout2, Accessor2> gcov) { std::mdspan, Layout2, Accessor2> gamma_ij; for (int i = 1; i < 4; i++) { for (int j = 1; j < 4; j++) { - gamma_ij(i-1, j-1) = gcov(i, j); + gamma_ij(i - 1, j - 1) = gcov(i, j); } } return gamma_ij; } - inline static constexpr std::mdspan, Layout2, Accessor2> get_gammaijUP( - std::mdspan, Layout2, Accessor2> gcon) { + inline static constexpr std::mdspan, Layout2, Accessor2> + get_gammaijUP(std::mdspan, Layout2, Accessor2> gcon) { Tscal alpha = get_alpha(gcon); Tvec betaUP = get_betaUP(gcon); std::mdspan, Layout2, Accessor2> gamma_ij; for (int i = 1; i < 4; i++) { for (int j = 1; j < 4; j++) { - gamma_ij(i-1, j-1) = gcon(i, j) + betaUP[i-1] * betaUP[j-1] / (-gcon(0, 0)); + gamma_ij(i - 1, j - 1) + = gcon(i, j) + betaUP[i - 1] * betaUP[j - 1] / (-gcon(0, 0)); } } return gamma_ij; - } inline static constexpr Tscal get_gamma( @@ -121,17 +121,18 @@ namespace shamphys { return mat_det_33(gamma_ij); } - inline static constexpr Tscal get_lorentz_factor( - Tvec momentum, Tscal enthalpy) { + inline static constexpr Tscal get_lorentz_factor(Tvec momentum, Tscal enthalpy) { return sycl::sqrt(1. + GR_dot(momentum, momentum) / (enthalpy * enthalpy)); } inline static constexpr Tscal get_U0( - Tvec momentum, Tscal enthalpy, std::mdspan, Layout2, Accessor2> gcon) { + Tvec momentum, + Tscal enthalpy, + std::mdspan, Layout2, Accessor2> gcon) { // U0 = gamma / alpha wrong - Tscal alpha = get_alpha(gcon); + Tscal alpha = get_alpha(gcon); Tscal lorentz_factor = get_lorentz_factor(momentum, enthalpy); return lorentz_factor / alpha; } From 851b9d4851ba25b8fa19ba4499a89a1c2341817a Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Thu, 12 Mar 2026 17:53:12 +0900 Subject: [PATCH 18/18] fix typos + add doc --- src/shamphys/include/shamphys/GRUtils.hpp | 140 +++++++++++++++++++--- 1 file changed, 124 insertions(+), 16 deletions(-) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index c45db13249..de6f76f051 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -34,12 +34,31 @@ namespace shamphys { class Accessor2> struct GR_physics { + /** + * @brief Compute the alpha factor from the spacetime metric. + * + * Alpha is defined as 1 / sqrt(-g00) where g00 is the + * 0,0 component of the spacetime metric. + * + * @param gcon The spacetime metric as a 4x4 matrix. + * @return The alpha factor. + */ inline static constexpr Tscal get_alpha( std::mdspan, Layout2, Accessor2> gcon) { // alpha = 1 / sqrt(-g00) return 1. / sycl::sqrt(-gcon(0, 0)); } + /** + * @brief Compute the betaUP factor from the spacetime metric. + * + * BetaUP is defined as (g01, g02, g03) / (alpha^2) where + * g0i is the i-th component of the first row of the spacetime + * metric and alpha is the Lorentz factor. + * + * @param gcon The spacetime metric as a 4x4 matrix. + * @return The betaUP factor. + */ inline static constexpr Tvec get_betaUP( std::mdspan, Layout2, Accessor2> gcon) { Tscal alpha = get_alpha(gcon); @@ -53,6 +72,15 @@ namespace shamphys { return betaUP; } + /** + * @brief Compute the betaDOWN factor from the spacetime metric. + * + * BetaDOWN is defined as (g01, g02, g03) where g0i is the i-th + * component of the first row of the spacetime metric. + * + * @param gcov The spacetime metric as a 4x4 matrix. + * @return The betaDOWN factor. + */ inline static constexpr Tvec get_betaDOWN( std::mdspan, Layout2, Accessor2> gcov) { Tvec betaDOWN = {0., 0., 0.}; @@ -64,21 +92,56 @@ namespace shamphys { return betaDOWN; } - /// placeholder + /** + * @brief Compute the dot product of two vectors in the spacetime metric. + * + * @param a The first vector. + * @param b The second vector. + * @param gcov The spacetime metric as a 4x4 matrix. + * @return The dot product of the two vectors in the spacetime metric. + */ inline static constexpr Tscal GR_dot( Tvec a, Tvec b, std::mdspan, Layout2, Accessor2> gcov) { - return 0; + + Tscal result = 0; + for (int i = 0; i < 4; i++) { + for (int j = 0; j < 4; j++) { + result += gcov(i + 1, j + 1) * a[i] * b[j]; + } + } + return result; + } + + /** + * @brief Compute the dot product of two spatial vectors in the spacetime metric. + * + * @param a The first spatial vector. + * @param b The second spatial vector. + * @param gamma_ij The spatial part of the spacetime metric as a 3x3 matrix. + * @return The dot product of the two spatial vectors in the spacetime metric. + */ + inline static constexpr Tscal GR_dot_spatial( + Tvec a, + Tvec b, + std::mdspan, Layout2, Accessor2> gamma_ij) { + Tscal result = 0; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { - result += gcov(i + 1, j + 1) * a[i] * b[j]; + result += gamma_ij(i + 1, j + 1) * a[i] * b[j]; } } return result; } + /** + * @brief Compute the spatial part of the spacetime metric from the full 4x4 metric. + * + * @param gcov The spacetime metric as a 4x4 matrix. + * @return The spatial part of the spacetime metric as a 3x3 matrix. + */ inline static constexpr std::mdspan, Layout2, Accessor2> get_gammaijDOWN(std::mdspan, Layout2, Accessor2> gcov) { @@ -92,6 +155,13 @@ namespace shamphys { return gamma_ij; } + /** + * @brief Compute the contravariant spatial part of the spacetime metric from the full 4x4 + * metric. + * + * @param gcon The spacetime metric as a 4x4 matrix. + * @return The contravariant spatial part of the spacetime metric as a 3x3 matrix. + */ inline static constexpr std::mdspan, Layout2, Accessor2> get_gammaijUP(std::mdspan, Layout2, Accessor2> gcon) { @@ -108,35 +178,73 @@ namespace shamphys { return gamma_ij; } + /** + * @brief Compute the determinant of the spatial part of the spacetime metric. + * + * This function takes a 3x3 matrix representing the spatial part of the + * spacetime metric and returns its determinant. + * + * @param gamma_ij The spatial part of the spacetime metric as a 3x3 matrix. + * @return The determinant of the spatial part of the spacetime metric. + */ inline static constexpr Tscal get_gamma( - std::mdspan, Layout2, Accessor2> gcov) { - - std::mdspan, Layout2, Accessor2> gamma_ij; - for (int i = 1; i < 4; i++) { - for (int j = 1; j < 4; j++) { - gamma_ij(i, j) = gcov(i, j); - } - } + std::mdspan, Layout2, Accessor2> gamma_ij) { return mat_det_33(gamma_ij); } - inline static constexpr Tscal get_lorentz_factor(Tvec momentum, Tscal enthalpy) { + /** + * @brief Compute the Lorentz factor of a fluid element given its momentum and the spacetime + * metric. + * + * This function takes the momentum of a fluid element, its enthalpy, and the spacetime + * metric as inputs and returns the Lorentz factor of the fluid element. + * + * @param momentum The momentum of the fluid element. + * @param enthalpy The enthalpy of the fluid element. + * @param gcov The spacetime metric as a 4x4 matrix. + * @return The Lorentz factor of the fluid element. + */ + inline static constexpr Tscal get_lorentz_factor( + Tvec momentum, + Tscal enthalpy, + std::mdspan, Layout2, Accessor2> gcov) { - return sycl::sqrt(1. + GR_dot(momentum, momentum) / (enthalpy * enthalpy)); + return sycl::sqrt(1. + GR_dot(momentum, momentum, gcov) / (enthalpy * enthalpy)); } + /// @brief Compute the time component of the 4-velocity of a fluid element given its + /// momentum, enthalpy, and the spacetime metric. + /// + /// This function takes the momentum of a fluid element, its enthalpy, and the spacetime + /// metric as inputs and returns the time component of the 4-velocity of the fluid element. + /// + /// @param momentum The momentum of the fluid element. + /// @param enthalpy The enthalpy of the fluid element. + /// @param gcon The spacetime metric as a 4x4 matrix. + /// @param gcov The spacetime metric as a 4x4 matrix. + /// @return The time component of the 4-velocity of the fluid element. inline static constexpr Tscal get_U0( Tvec momentum, Tscal enthalpy, - std::mdspan, Layout2, Accessor2> gcon) { + std::mdspan, Layout2, Accessor2> gcon, + std::mdspan, Layout2, Accessor2> gcov) { - // U0 = gamma / alpha wrong Tscal alpha = get_alpha(gcon); - Tscal lorentz_factor = get_lorentz_factor(momentum, enthalpy); + Tscal lorentz_factor = get_lorentz_factor(momentum, enthalpy, gcov); return lorentz_factor / alpha; } + /** + * @brief Compute the fluid velocity in the local rest frame. + * + * This function takes the fluid velocity in the global frame and the spacetime metric as + * inputs and returns the fluid velocity in the local rest frame. + * + * @param vxyz The fluid velocity in the global frame. + * @param gcon The spacetime metric as a 4x4 matrix. + * @return The fluid velocity in the local rest frame. + */ inline static constexpr Tvec get_V( Tvec vxyz, std::mdspan, Layout2, Accessor2> gcon) { Tscal alpha = get_alpha(gcon);