diff --git a/src/shammodels/sph/CMakeLists.txt b/src/shammodels/sph/CMakeLists.txt index 295ff990ca..66072df4e5 100644 --- a/src/shammodels/sph/CMakeLists.txt +++ b/src/shammodels/sph/CMakeLists.txt @@ -63,6 +63,8 @@ set(Sources src/modules/self_gravity/SGSFMMPlummer.cpp src/modules/render/CartesianRender.cpp src/modules/render/RenderFieldGetter.cpp + src/modules/ConsToPrim.cpp + src/modules/io/VTKDump.cpp ) diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp new file mode 100644 index 0000000000..d552e067a6 --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp @@ -0,0 +1,63 @@ +// -------------------------------------------------------// +// +// 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 ConsToPrim.hpp + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + @brief get primitive variables (rho, vel, internal energy) from conserved variables (rho*, + momentum, entropy) + */ + +#include "shambackends/vec.hpp" +#include "shamrock/solvergraph/Field.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarsEdge.hpp" +#include + +namespace shammodels::sph::modules { + template + class NodeConsToPrim : public shamrock::solvergraph::INode { + using Tscal = shambase::VecComponent; + Tscal gamma; + + public: + NodeConsToPrim(Tscal gamma) : gamma(gamma) {} + + using GcovEdge = std::mdspan, Layout, Accessor>; + +// X_RO(std::mdspan, Layout, Accessor>, gcov) +#define NODE_CONS_TO_PRIM(X_RO, X_RW) \ + /* inputs */ \ + X_RO(shamrock::solvergraph::Indexes, sizes) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_rhostar) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_momentum) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_K) \ + X_RO(GcovEdge, gcon) \ + X_RO(GcovEdge, gcov) \ + \ + /* outputs */ \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_rho) \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_vel) \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_u) \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_P) + + EXPAND_NODE_EDGES(NODE_CONS_TO_PRIM) +#undef NODE_CONS_TO_PRIM + + void _impl_evaluate_internal(); + + inline virtual std::string _impl_get_label() const { return "ConsToPrim"; }; + + virtual std::string _impl_get_tex() const; + }; +} // namespace shammodels::sph::modules diff --git a/src/shammodels/sph/src/modules/ConsToPrim.cpp b/src/shammodels/sph/src/modules/ConsToPrim.cpp new file mode 100644 index 0000000000..f50aadca92 --- /dev/null +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -0,0 +1,176 @@ +// -------------------------------------------------------// +// +// 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 +// +// -------------------------------------------------------// + +/** + * @file ConsToPrim.cpp + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shammodels/sph/modules/ConsToPrim.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shamcomm/logs.hpp" +#include "shamphys/GRUtils.hpp" +#include "shamphys/metrics.hpp" +#include "shamrock/patch/PatchDataField.hpp" +#include "shamsys/NodeInstance.hpp" + +namespace { + + template + struct KernelConsToPrim { + using Tscal = shambase::VecComponent; + }; + +} // namespace + +namespace shammodels::sph::modules { + + template + void NodeConsToPrim::_impl_evaluate_internal() { + auto edges = get_edges(); + auto &thread_counts = edges.sizes.indexes; + + edges.spans_rhostar.check_sizes(thread_counts); + edges.spans_momentum.check_sizes(thread_counts); + edges.spans_K.check_sizes(thread_counts); + + edges.spans_rho.check_sizes(thread_counts); + edges.spans_vel.check_sizes(thread_counts); + edges.spans_u.check_sizes(thread_counts); + edges.spans_P.check_sizes(thread_counts); + + auto &rhostar = edges.spans_rhostar.get_spans(); + auto &momentum = edges.spans_momentum.get_spans(); + auto &K = edges.spans_K.get_spans(); + + auto &rho = edges.spans_rho.get_spans(); + auto &vel = edges.spans_vel.get_spans(); + auto &u = edges.spans_u.get_spans(); + auto &P = edges.spans_P.get_spans(); + + auto &gcov = edges.gcov; + auto &gcon = edges.gcon; + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + + sham::distributed_data_kernel_call( + dev_sched, + sham::DDMultiRef{rhostar, momentum, K}, + sham::DDMultiRef{rho, vel, u, P}, + thread_counts, + [gamma = this->gamma, gcov = edges.gcov]( + u32 id_a, + Tscal *__restrict rhostar, + Tvec *__restrict momentum, + Tscal *__restrict K, + + Tscal *__restrict rho, + Tvec *__restrict vel, + Tscal *__restrict u, + Tscal *__restrict P) { + // on patch, no need of neighbours + + // get metric quantities + const Tscal sqrt_g = get_sqrt_g(vel[id_a], gcov); + const Tscal inv_sqrt_g = 1. / sqrt_g; + const Tscal alpha = get_alpha(gcon); + const Tscal sqrt_gamma_inv = alpha * inv_sqrt_g; + const Tvec betaUP = get_betaUP(gcon); + const Tvec betaDOWN = get_betaDOWN(gcon); + const std::mdspan, Layout, Accessor> gammaijUP + = get_gammaijUP(gcov); + + Tscal rho_a = rho[id_a]; + const Tscal gamfac = gamma / (gamma - 1.); + Tscal w = 1 + gamfac * P[id_a] / rho[id_a]; // initial guess + Tscal pp = 0; + for (u32 i = 0; i < 3; i++) { + for (u32 j = 0; j < 3; j++) { + pp += momentum[id_a][i] * gammaijUP(i, j) * momentum[id_a][j]; + } + } + + Tscal lorentz_factor = get_lorentz_factor(momentum[id_a], w, gcov); + + // Tscal Kstar = K[id_a] + 0.5 * pp / rhostar[id_a]; + + bool converged = false; + constexpr Tscal tol = Tscal(1e-12); + constexpr u32 Nitermax = 100; + // compute u + // iterate + u32 Niter = 0; + + do { + const Tscal w_old = w; + + lorentz_factor = get_lorentz_factor(momentum[id_a], w, gcov); + + // eq 97 + rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; + + // eq 62 + P[id_a] = K[id_a] * sycl::pow(rho_a, gamma); + + // eq B4 + const Tscal f = (1. + gamfac * P[id_a] / rho_a) - w_old; + + // eq B5 ... maybe should use B6 + const Tscal df = -1 + + gamfac + * (1. + - pp * P[id_a] + / (lorentz_factor * lorentz_factor * w_old + * w_old * w_old * rho_a)); + + w = w_old - f / df; + + converged = (sycl::fabs(w - w_old) / w < tol); + Niter++; + } while (Niter < Nitermax and !converged); + + if (converged) { + // compute P, v, u with the last expression of enthalpy + lorentz_factor = get_lorentz_factor(momentum[id_a], w, gcov); + rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; + rho[id_a] = rho_a; + P[id_a] = K[id_a] * sycl::pow(rho_a, gamma); + + Tvec v3d = alpha * momentum[id_a] / (w * lorentz_factor) - betaDOWN; + + // Raise index from down to up + for (u32 i = 0; i < 3; i++) { + Tscal vi = 0.; + for (u32 j = 0; j < 3; j++) { + vi += gammaijUP(i, j) * v3d[j]; + } + vel[id_a][i] = vi; + } + + u[id_a] = (P[id_a] < Tscal(1e-30)) ? Tscal(0) + : P[id_a] / (rho_a * (gamma - Tscal(1))); + } else { + logger::err_ln( + "GRSPH", + "the enthalpy iterator is not converged after", + Niter, + "iterations"); + } + }); + } + + template + std::string NodeConsToPrim::_impl_get_tex() const { + + return "TODO"; + } + +} // namespace shammodels::sph::modules diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp new file mode 100644 index 0000000000..ea047a1dcb --- /dev/null +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -0,0 +1,183 @@ +// -------------------------------------------------------// +// +// 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/aliases_int.hpp" +#include "shambase/assert.hpp" +#include "shambackends/math.hpp" +#include "shamunits/Constants.hpp" +#include +#include + +namespace shamphys { + + template< + class Tvec, + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> + struct GR_physics { + + inline static constexpr Tscal get_alpha( + std::mdspan, Layout2, Accessor2> gcon) { + // alpha = 1 / sqrt(-g00) + return 1. / sycl::sqrt(-gcon(0, 0)); + } + + inline static constexpr Tvec get_betaUP( + std::mdspan, Layout2, Accessor2> gcon) { + Tscal alpha = get_alpha(gcon); + Tscal alpha2 = alpha * alpha; + Tvec betaUP = {0., 0., 0.}; + + 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) { + 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( + std::mdspan, Layout2, Accessor2> a, + std::mdspan, Layout2, Accessor2> b, + std::mdspan, Layout2, Accessor2> gcov) { + + 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; + } + + 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 += gamma_ij(i, j) * a[i] * b[j]; + } + } + return result; + } + + 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); + } + } + + return gamma_ij; + } + + 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)); + } + } + + return gamma_ij; + } + + 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_lorentz_factor( + Tvec momentum, + Tscal enthalpy, + std::mdspan, Layout2, Accessor2> gcov) { + + return sycl::sqrt(1. + GR_dot(momentum, momentum, gcov) / (enthalpy * enthalpy)); + } + + inline static constexpr Tscal get_U0( + Tvec momentum, + Tscal enthalpy, + 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, gcov); + return lorentz_factor / alpha; + } + + inline static constexpr Tvec get_V( + Tvec vxyz, std::mdspan, Layout2, Accessor2> gcon) { + Tscal alpha = get_alpha(gcon); + Tvec betaUP = get_betaUP(gcon); + Tvec V = (vxyz + betaUP) / alpha; + + 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 + // for now (Kerr), this is enough + 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 + // for now (Kerr), this is enough + return 1; + } + }; + +} // namespace shamphys diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp new file mode 100644 index 0000000000..604ed61ac1 --- /dev/null +++ b/src/shamphys/include/shamphys/metrics.hpp @@ -0,0 +1,521 @@ +// -------------------------------------------------------// +// +// 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 metrics.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 Kerr { + + Tscal a; + Tscal bh_mass; + Tscal rs; + + Kerr(Tscal spin, Tscal mass) : a(spin), bh_mass(mass), rs(2. * mass) {} + }; + + struct Schwarzschild {}; + + /******************** The Kerr metric ***********************/ + // ################in CARTESIAN-like form ################## + 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 + const Tscal a = kerr.a; + const Tscal rs = kerr.rs; + + // pos coords + Tscal x = pos(0), y = pos(1), 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 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 + 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); + } + + 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 + const Tscal a = kerr.a; + const Tscal rs = kerr.rs; + + // pos coords + Tscal x = pos(0), y = pos(1), 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 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 + 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); + } + + template< + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> + inline void metric_cartesian_derivatives_impl( + 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 + const Tscal a = kerr.a; + const Tscal rs = kerr.rs; + // 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 + // 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); + } + + /******************** frontend ***********************/ + 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) { + get_cartesian_covariant_metric_impl< + Tscal, + SizeType, + Layout1, + Layout2, + Accessor1, + Accessor2>(MetricTag{}, pos, g); + } + + 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) { + get_cartesian_contravariant_metric_impl< + Tscal, + SizeType, + Layout1, + Layout2, + Accessor1, + Accessor2>(MetricTag{}, pos, g); + } + + template< + class MetricTag, + class Tscal, + class SizeType, + class Layout1, + class Layout2, + class Accessor1, + class Accessor2> + inline void metric_cartesian_derivatives( + 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