From 8e5898b0a174604db74049455505f611db880678 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 18 Feb 2026 15:06:26 +0900 Subject: [PATCH 01/13] baseline --- src/shammodels/ramses/CMakeLists.txt | 1 + .../shammodels/sph/modules/ConsToPrim.hpp | 73 ++++++++++++++ src/shammodels/sph/src/modules/ConsToPrim.cpp | 96 +++++++++++++++++++ 3 files changed, 170 insertions(+) create mode 100644 src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp create mode 100644 src/shammodels/sph/src/modules/ConsToPrim.cpp diff --git a/src/shammodels/ramses/CMakeLists.txt b/src/shammodels/ramses/CMakeLists.txt index 4d5e89756d..696cccde65 100644 --- a/src/shammodels/ramses/CMakeLists.txt +++ b/src/shammodels/ramses/CMakeLists.txt @@ -44,6 +44,7 @@ set(Sources src/modules/TransformGhostLayer.cpp src/modules/FuseGhostLayer.cpp src/modules/render/GridRender.cpp + src/modules/ConsToPrim.cpp ) if(SHAMROCK_USE_SHARED_LIB) 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..3ae2f5b713 --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp @@ -0,0 +1,73 @@ +// -------------------------------------------------------// +// +// 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 conserved variables (rho*, momentum, entropy) from primitive variables (rho, vel, internal energy) + */ + +#include "shambackends/vec.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" + +namespace shammodels::basegodunov::modules { + template + class NodeConsToPrim : public shamrock::solvergraph::INode { + using Tscal = shambase::VecComponent; + u32 block_size; + Tscal gamma; + + public: + NodeConsToPrim(u32 block_size, Tscal gamma) : block_size(block_size), gamma(gamma) {} + + struct Edges { + const shamrock::solvergraph::Indexes &sizes; + const shamrock::solvergraph::IFieldSpan &spans_rhostar; + const shamrock::solvergraph::IFieldSpan &spans_momentum; + const shamrock::solvergraph::IFieldSpan &spans_K; + shamrock::solvergraph::IFieldSpan &spans_vel; + shamrock::solvergraph::IFieldSpan &spans_u; + shamrock::solvergraph::IFieldSpan &spans_P; + }; + + inline void set_edges( + std::shared_ptr> sizes, + std::shared_ptr> spans_rhostar, + std::shared_ptr> spans_momentum, + std::shared_ptr> spans_K, + std::shared_ptr> spans_vel, + std::shared_ptr> spans_u, + std::shared_ptr> spans_P) { + __internal_set_ro_edges({sizes, spans_rhostar, spans_momentum, spans_K}); + __internal_set_rw_edges({spans_vel, spans_u, spans_P}); + } + + inline Edges get_edges() { + return Edges{ + get_ro_edge>(0), + get_ro_edge>(1), + get_ro_edge>(2), + get_ro_edge>(3), + get_rw_edge>(0), + get_rw_edge>(1), + get_rw_edge>(2) + }; + } + + void _impl_evaluate_internal(); + + inline virtual std::string _impl_get_label() const { return "ConsToPrim"; }; + + virtual std::string _impl_get_tex() const; + }; +} // namespace shammodels::basegodunov::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..fa353191f3 --- /dev/null +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -0,0 +1,96 @@ +// -------------------------------------------------------// +// +// 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 "shammath/riemann.hpp" +#include "shamrock/patch/PatchDataField.hpp" +#include "shamsys/NodeInstance.hpp" + +namespace { + + template + struct KernelConsToPrim { + using Tscal = shambase::VecComponent; + }; + +} // namespace + +namespace shammodels::basegodunov::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_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 &u = edges.spans_u.get_spans(); + auto &P = edges.spans_P.get_spans(); + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + + sham::distributed_data_kernel_call( + dev_sched, + sham::DDMultiRef{rhostar, momentum, K}, + sham::DDMultiRef{u, P}, + thread_counts, + [gamma = this->gamma]( + u32 id_a, + Tscal *__restrict rhostar, + Tvec *__restrict momentum, + Tscal *__restrict K, + Tscal *__restrict u, + Tscal *__restrict P) { + // on patch, no need of neighbours + + //get metric + + //guess enthalpy w, with adiabatic EOS and previous values + Tscal w = gamma/(gamma-1) * P[id_a] / rhostar[id_a]; + bool converged = false; + //compute u + //iterate + u32 Niter = 0; + do { + + Tscal new_w = 0; + + converged = sycl::fabs(new_w - w) < 1e-6; + w = new_w; + Niter++; + } while (Niter < 100 or converged); + }); + + + } + + template + std::string NodeConsToPrim::_impl_get_tex() const { + + return "TODO"; + } + +} // namespace shammodels::basegodunov::modules + +template class shammodels::basegodunov::modules::NodeConsToPrim; From 8b6fcf51f9f806faf4fe7b8ea1b9a9d55e0111ea Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 18 Feb 2026 16:01:39 +0900 Subject: [PATCH 02/13] enthalpy iterator --- .../shammodels/sph/modules/ConsToPrim.hpp | 15 ++-- src/shammodels/sph/src/modules/ConsToPrim.cpp | 75 +++++++++++++------ 2 files changed, 63 insertions(+), 27 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp index 3ae2f5b713..d7f7eddb97 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp @@ -12,7 +12,8 @@ /** * @file ConsToPrim.hpp * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) - * @brief get conserved variables (rho*, momentum, entropy) from primitive variables (rho, vel, internal energy) + * @brief get conserved variables (rho*, momentum, entropy) from primitive variables (rho, vel, + * internal energy) */ #include "shambackends/vec.hpp" @@ -35,6 +36,7 @@ namespace shammodels::basegodunov::modules { const shamrock::solvergraph::IFieldSpan &spans_rhostar; const shamrock::solvergraph::IFieldSpan &spans_momentum; const shamrock::solvergraph::IFieldSpan &spans_K; + shamrock::solvergraph::IFieldSpan &spans_rho; shamrock::solvergraph::IFieldSpan &spans_vel; shamrock::solvergraph::IFieldSpan &spans_u; shamrock::solvergraph::IFieldSpan &spans_P; @@ -45,11 +47,12 @@ namespace shammodels::basegodunov::modules { std::shared_ptr> spans_rhostar, std::shared_ptr> spans_momentum, std::shared_ptr> spans_K, + std::shared_ptr> spans_rho, std::shared_ptr> spans_vel, std::shared_ptr> spans_u, std::shared_ptr> spans_P) { __internal_set_ro_edges({sizes, spans_rhostar, spans_momentum, spans_K}); - __internal_set_rw_edges({spans_vel, spans_u, spans_P}); + __internal_set_rw_edges({spans_rho, spans_vel, spans_u, spans_P}); } inline Edges get_edges() { @@ -58,10 +61,10 @@ namespace shammodels::basegodunov::modules { get_ro_edge>(1), get_ro_edge>(2), get_ro_edge>(3), - get_rw_edge>(0), - get_rw_edge>(1), - get_rw_edge>(2) - }; + get_rw_edge>(0), + get_rw_edge>(1), + get_rw_edge>(2), + get_rw_edge>(3)}; } void _impl_evaluate_internal(); diff --git a/src/shammodels/sph/src/modules/ConsToPrim.cpp b/src/shammodels/sph/src/modules/ConsToPrim.cpp index fa353191f3..2c36e2e381 100644 --- a/src/shammodels/sph/src/modules/ConsToPrim.cpp +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -16,6 +16,7 @@ #include "shammodels/sph/modules/ConsToPrim.hpp" #include "shambackends/kernel_call_distrib.hpp" +#include "shamcomm/logs.hpp" #include "shammath/riemann.hpp" #include "shamrock/patch/PatchDataField.hpp" #include "shamsys/NodeInstance.hpp" @@ -33,56 +34,88 @@ namespace shammodels::basegodunov::modules { template void NodeConsToPrim::_impl_evaluate_internal() { - auto edges = get_edges(); + 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 &rhostar = edges.spans_rhostar.get_spans(); auto &momentum = edges.spans_momentum.get_spans(); - auto &K = edges.spans_K.get_spans(); - auto &u = edges.spans_u.get_spans(); - auto &P = edges.spans_P.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 dev_sched = shamsys::instance::get_compute_scheduler_ptr(); sham::distributed_data_kernel_call( dev_sched, sham::DDMultiRef{rhostar, momentum, K}, - sham::DDMultiRef{u, P}, - thread_counts, + sham::DDMultiRef{rho, vel, u, P}, + thread_counts, [gamma = this->gamma]( - u32 id_a, - Tscal *__restrict rhostar, - Tvec *__restrict momentum, - Tscal *__restrict K, - Tscal *__restrict u, + 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 + // get metric + Tscal sqrt_g = get_sqrtg(gcov); + Tscal inv_sqrt_g = 1. / sqrt_g; + Tscal sqrt_gamma = get_sqrt_gamma(gcov); + Tscal sqrt_gamma_inv = alpha * inv_sqrt_g; - //guess enthalpy w, with adiabatic EOS and previous values - Tscal w = gamma/(gamma-1) * P[id_a] / rhostar[id_a]; + // guess enthalpy w, with adiabatic EOS and previous values + Tscal w = gamma / (gamma - 1) * P[id_a] / rhostar[id_a]; bool converged = false; - //compute u - //iterate + // compute u + // iterate u32 Niter = 0; do { + // get values of density and pressure from alod w + Tscal lorentz_factor + = sycl::sqrt(1. + sycl::dot(momentum[id_a], momentum[id_a] / (w * w))); + + rho[id_a] = sqrt_gamma_inv * rhostar[id_a] / lorentz_factor; + Tscal polyk = 1.; + P[id_a] = (gamma - 1.) * rho[id_a] * polyk; - Tscal new_w = 0; + Tscal new_w = 1; converged = sycl::fabs(new_w - w) < 1e-6; - w = new_w; + w = new_w; Niter++; } while (Niter < 100 or converged); - }); - + if (converged) { + Tvec v3d = alpha * momentum[id_a] / (w * lorentz_factor) - betadown; + // Raise index from down to up + for (u32 i = 0; i < 4; i++) { + vel[id_a][i] = sycl::dot(gammaijUP[:, i], v3d); + } + } else { + logger::err_ln( + "GRSPH", + "the enthalpy iterator is not converged after", + Niter, + "iterations"); + } + }); } template From 6a2996ceda2390979c0279cbf9a1528884cc7cd5 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 11 Mar 2026 15:18:04 +0900 Subject: [PATCH 03/13] fix types and typos --- src/shammodels/ramses/CMakeLists.txt | 1 - src/shammodels/sph/CMakeLists.txt | 1 + .../shammodels/sph/modules/ConsToPrim.hpp | 21 +++++++++-------- src/shammodels/sph/src/modules/ConsToPrim.cpp | 23 ++++++++++++------- 4 files changed, 28 insertions(+), 18 deletions(-) diff --git a/src/shammodels/ramses/CMakeLists.txt b/src/shammodels/ramses/CMakeLists.txt index 696cccde65..4d5e89756d 100644 --- a/src/shammodels/ramses/CMakeLists.txt +++ b/src/shammodels/ramses/CMakeLists.txt @@ -44,7 +44,6 @@ set(Sources src/modules/TransformGhostLayer.cpp src/modules/FuseGhostLayer.cpp src/modules/render/GridRender.cpp - src/modules/ConsToPrim.cpp ) if(SHAMROCK_USE_SHARED_LIB) diff --git a/src/shammodels/sph/CMakeLists.txt b/src/shammodels/sph/CMakeLists.txt index 7bf09cb47f..74aed19673 100644 --- a/src/shammodels/sph/CMakeLists.txt +++ b/src/shammodels/sph/CMakeLists.txt @@ -62,6 +62,7 @@ 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 index d7f7eddb97..9bf5527d04 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp @@ -20,18 +20,19 @@ #include "shamrock/solvergraph/IFieldSpan.hpp" #include "shamrock/solvergraph/INode.hpp" #include "shamrock/solvergraph/Indexes.hpp" +#include -namespace shammodels::basegodunov::modules { - template +namespace shammodels::sph::modules { + template class NodeConsToPrim : public shamrock::solvergraph::INode { using Tscal = shambase::VecComponent; - u32 block_size; Tscal gamma; public: - NodeConsToPrim(u32 block_size, Tscal gamma) : block_size(block_size), gamma(gamma) {} + NodeConsToPrim(Tscal gamma) : gamma(gamma) {} struct Edges { + std::mdspan, Layout, Accessor> &gcov; const shamrock::solvergraph::Indexes &sizes; const shamrock::solvergraph::IFieldSpan &spans_rhostar; const shamrock::solvergraph::IFieldSpan &spans_momentum; @@ -43,6 +44,7 @@ namespace shammodels::basegodunov::modules { }; inline void set_edges( + std::mdspan, Layout, Accessor> &gcov, std::shared_ptr> sizes, std::shared_ptr> spans_rhostar, std::shared_ptr> spans_momentum, @@ -57,10 +59,11 @@ namespace shammodels::basegodunov::modules { inline Edges get_edges() { return Edges{ - get_ro_edge>(0), - get_ro_edge>(1), - get_ro_edge>(2), - get_ro_edge>(3), + get_ro_edge>(0), // fix types + get_ro_edge>(1), + get_ro_edge>(2), + get_ro_edge>(3), + get_ro_edge>(4), get_rw_edge>(0), get_rw_edge>(1), get_rw_edge>(2), @@ -73,4 +76,4 @@ namespace shammodels::basegodunov::modules { virtual std::string _impl_get_tex() const; }; -} // namespace shammodels::basegodunov::modules +} // namespace shammodels::sph::modules diff --git a/src/shammodels/sph/src/modules/ConsToPrim.cpp b/src/shammodels/sph/src/modules/ConsToPrim.cpp index 2c36e2e381..215c21927e 100644 --- a/src/shammodels/sph/src/modules/ConsToPrim.cpp +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -18,6 +18,8 @@ #include "shambackends/kernel_call_distrib.hpp" #include "shamcomm/logs.hpp" #include "shammath/riemann.hpp" +#include "shamphys/GRUtils.hpp" +#include "shamphys/metrics.hpp" #include "shamrock/patch/PatchDataField.hpp" #include "shamsys/NodeInstance.hpp" @@ -30,10 +32,10 @@ namespace { } // namespace -namespace shammodels::basegodunov::modules { +namespace shammodels::sph::modules { - template - void NodeConsToPrim::_impl_evaluate_internal() { + template + void NodeConsToPrim::_impl_evaluate_internal() { auto edges = get_edges(); auto &thread_counts = edges.sizes.indexes; @@ -55,6 +57,8 @@ namespace shammodels::basegodunov::modules { auto &u = edges.spans_u.get_spans(); auto &P = edges.spans_P.get_spans(); + auto &gcov = edges.gcov; + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); sham::distributed_data_kernel_call( @@ -78,6 +82,7 @@ namespace shammodels::basegodunov::modules { Tscal sqrt_g = get_sqrtg(gcov); Tscal inv_sqrt_g = 1. / sqrt_g; Tscal sqrt_gamma = get_sqrt_gamma(gcov); + Tscal alpha = get_alpha(gcov); Tscal sqrt_gamma_inv = alpha * inv_sqrt_g; // guess enthalpy w, with adiabatic EOS and previous values @@ -86,9 +91,10 @@ namespace shammodels::basegodunov::modules { // compute u // iterate u32 Niter = 0; + Tscal lorentz_factor; do { // get values of density and pressure from alod w - Tscal lorentz_factor + lorentz_factor = sycl::sqrt(1. + sycl::dot(momentum[id_a], momentum[id_a] / (w * w))); rho[id_a] = sqrt_gamma_inv * rhostar[id_a] / lorentz_factor; @@ -118,12 +124,13 @@ namespace shammodels::basegodunov::modules { }); } - template - std::string NodeConsToPrim::_impl_get_tex() const { + template + std::string NodeConsToPrim::_impl_get_tex() const { return "TODO"; } -} // namespace shammodels::basegodunov::modules +} // namespace shammodels::sph::modules -template class shammodels::basegodunov::modules::NodeConsToPrim; +template class shammodels::sph::modules:: + NodeConsToPrim; From c91b588f2032f5386c7932f4a1684c8d769be3d3 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 11 Mar 2026 15:20:38 +0900 Subject: [PATCH 04/13] from PR1 --- src/shamphys/include/shamphys/GRUtils.hpp | 112 +++++ src/shamphys/include/shamphys/metrics.hpp | 521 ++++++++++++++++++++++ 2 files changed, 633 insertions(+) create mode 100644 src/shamphys/include/shamphys/GRUtils.hpp create mode 100644 src/shamphys/include/shamphys/metrics.hpp diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp new file mode 100644 index 0000000000..c3a83153ac --- /dev/null +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -0,0 +1,112 @@ +// -------------------------------------------------------// +// +// 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> 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, 2) * alpha2; + betaUP[2] = gcov(0, 3) * alpha2; + + return betaUP; + } + + inline static constexpr Tscal GR_dot( + Tvec a, + Tvec b, + std::mdspan, Layout2, Accessor2> gcov) { + return 0; + // 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 = get_gamma(gcov); + return gamma / 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 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 \ No newline at end of file diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp new file mode 100644 index 0000000000..cd72b3071a --- /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 \ No newline at end of file From c71fcc378c68d9ac57ddba3def996f788099b9f3 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 11 Mar 2026 15:51:23 +0900 Subject: [PATCH 05/13] compiling first version --- src/shammodels/sph/src/modules/ConsToPrim.cpp | 11 +++--- src/shamphys/include/shamphys/GRUtils.hpp | 39 ++++++++++++++++++- 2 files changed, 43 insertions(+), 7 deletions(-) diff --git a/src/shammodels/sph/src/modules/ConsToPrim.cpp b/src/shammodels/sph/src/modules/ConsToPrim.cpp index 215c21927e..758926a475 100644 --- a/src/shammodels/sph/src/modules/ConsToPrim.cpp +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -109,10 +109,12 @@ namespace shammodels::sph::modules { } while (Niter < 100 or converged); if (converged) { - Tvec v3d = alpha * momentum[id_a] / (w * lorentz_factor) - betadown; + Tvec betaDOWN = get_betaDOWN(gcov); + Tvec gammaijUP = get_gammaijUP(gcov); + Tvec v3d = alpha * momentum[id_a] / (w * lorentz_factor) - betaDOWN; // Raise index from down to up - for (u32 i = 0; i < 4; i++) { - vel[id_a][i] = sycl::dot(gammaijUP[:, i], v3d); + for (u32 i = 0; i < 3; i++) { + vel[id_a][i] = sycl::dot(gammaijUP[i + 1], v3d); } } else { logger::err_ln( @@ -131,6 +133,3 @@ namespace shammodels::sph::modules { } } // namespace shammodels::sph::modules - -template class shammodels::sph::modules:: - NodeConsToPrim; diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index c3a83153ac..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) { @@ -109,4 +146,4 @@ namespace shamphys { } }; -} // namespace shamphys \ No newline at end of file +} // namespace shamphys From 9d8aab056c0eaaef94ab8a46846842dd0f0ea478 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:50:58 +0000 Subject: [PATCH 06/13] [autofix.ci] automatic fix: pre-commit hooks --- src/shamphys/include/shamphys/metrics.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/shamphys/include/shamphys/metrics.hpp b/src/shamphys/include/shamphys/metrics.hpp index cd72b3071a..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; @@ -518,4 +518,4 @@ namespace shamphys { MetricTag{}, pos, dgcovdx, dgcovdy, dgcovdz); } -} // namespace shamphys \ No newline at end of file +} // namespace shamphys From db029f538db4c4f17395cd1b32bff09cabc154cd 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:51:05 +0000 Subject: [PATCH 07/13] [gh-action] trigger CI with empty commit From 840323c85630c21428473f5e4f8a157247ee4c2c Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Wed, 11 Mar 2026 17:33:11 +0900 Subject: [PATCH 08/13] new solvergraph edges intro --- .../shammodels/sph/modules/ConsToPrim.hpp | 56 ++++++------------- src/shammodels/sph/src/modules/ConsToPrim.cpp | 1 - 2 files changed, 18 insertions(+), 39 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp index 9bf5527d04..33008bf6a7 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp @@ -17,9 +17,11 @@ */ #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 { @@ -31,44 +33,22 @@ namespace shammodels::sph::modules { public: NodeConsToPrim(Tscal gamma) : gamma(gamma) {} - struct Edges { - std::mdspan, Layout, Accessor> &gcov; - const shamrock::solvergraph::Indexes &sizes; - const shamrock::solvergraph::IFieldSpan &spans_rhostar; - const shamrock::solvergraph::IFieldSpan &spans_momentum; - const shamrock::solvergraph::IFieldSpan &spans_K; - shamrock::solvergraph::IFieldSpan &spans_rho; - shamrock::solvergraph::IFieldSpan &spans_vel; - shamrock::solvergraph::IFieldSpan &spans_u; - shamrock::solvergraph::IFieldSpan &spans_P; - }; - - inline void set_edges( - std::mdspan, Layout, Accessor> &gcov, - std::shared_ptr> sizes, - std::shared_ptr> spans_rhostar, - std::shared_ptr> spans_momentum, - std::shared_ptr> spans_K, - std::shared_ptr> spans_rho, - std::shared_ptr> spans_vel, - std::shared_ptr> spans_u, - std::shared_ptr> spans_P) { - __internal_set_ro_edges({sizes, spans_rhostar, spans_momentum, spans_K}); - __internal_set_rw_edges({spans_rho, spans_vel, spans_u, spans_P}); - } - - inline Edges get_edges() { - return Edges{ - get_ro_edge>(0), // fix types - get_ro_edge>(1), - get_ro_edge>(2), - get_ro_edge>(3), - get_ro_edge>(4), - get_rw_edge>(0), - get_rw_edge>(1), - get_rw_edge>(2), - get_rw_edge>(3)}; - } +// 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) \ + \ + /* 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(); diff --git a/src/shammodels/sph/src/modules/ConsToPrim.cpp b/src/shammodels/sph/src/modules/ConsToPrim.cpp index 758926a475..b3294b07c9 100644 --- a/src/shammodels/sph/src/modules/ConsToPrim.cpp +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -17,7 +17,6 @@ #include "shammodels/sph/modules/ConsToPrim.hpp" #include "shambackends/kernel_call_distrib.hpp" #include "shamcomm/logs.hpp" -#include "shammath/riemann.hpp" #include "shamphys/GRUtils.hpp" #include "shamphys/metrics.hpp" #include "shamrock/patch/PatchDataField.hpp" From b42c17d22a692aa282198309d7ffddf36a412002 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Thu, 12 Mar 2026 16:36:28 +0900 Subject: [PATCH 09/13] finished implementation of enthalpy iterator --- src/shammodels/sph/src/modules/ConsToPrim.cpp | 85 +++++++++++++------ 1 file changed, 57 insertions(+), 28 deletions(-) diff --git a/src/shammodels/sph/src/modules/ConsToPrim.cpp b/src/shammodels/sph/src/modules/ConsToPrim.cpp index b3294b07c9..a1a6c9b23f 100644 --- a/src/shammodels/sph/src/modules/ConsToPrim.cpp +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -77,50 +77,79 @@ namespace shammodels::sph::modules { Tscal *__restrict P) { // on patch, no need of neighbours - // get metric - Tscal sqrt_g = get_sqrtg(gcov); - Tscal inv_sqrt_g = 1. / sqrt_g; - Tscal sqrt_gamma = get_sqrt_gamma(gcov); - Tscal alpha = get_alpha(gcov); - Tscal sqrt_gamma_inv = alpha * inv_sqrt_g; - - // guess enthalpy w, with adiabatic EOS and previous values - Tscal w = gamma / (gamma - 1) * P[id_a] / rhostar[id_a]; - bool converged = false; + // get metric quantities + const Tscal sqrt_g = get_sqrtg(gcov); + const Tscal inv_sqrt_g = 1. / sqrt_g; + const Tscal sqrt_gamma = get_sqrt_gamma(gcov); + const Tscal alpha = get_alpha(gcov); + const Tscal sqrt_gamma_inv = alpha * inv_sqrt_g; + const Tvec betaUP = get_betaUP(gcov); + const Tvec betaDOWN = get_betaDOWN(gcov); // careul cov con + const std::mdspan, Layout, Accessor> gammaijUP + = get_gammaijUP(gcov); + + Tscal pp = 0; + for (u32 i = 0; i < 3; i++) { + pp += momentum[id_a][i] * sycl::dot(gammaijUP[i + 1], momentum[id_a]); + } + + Tscal lorentz_factor = get_lorentz_factor(momentum[id_a], gcov); + + Tscal rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; + const Tscal gamfac = gamma / (gamma - 1.); + Tscal w = 1 + gamfac * P[id_a] / rho[id_a]; // initial guess + + 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; - Tscal lorentz_factor; + do { - // get values of density and pressure from alod w - lorentz_factor - = sycl::sqrt(1. + sycl::dot(momentum[id_a], momentum[id_a] / (w * w))); + const Tscal w_old = w; - rho[id_a] = sqrt_gamma_inv * rhostar[id_a] / lorentz_factor; - Tscal polyk = 1.; - P[id_a] = (gamma - 1.) * rho[id_a] * polyk; + lorentz_factor = get_lorentz_factor(momentum[id_a], gcov); - Tscal new_w = 1; + // eq 97 + Tscal rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; - converged = sycl::fabs(new_w - w) < 1e-6; - w = new_w; + // 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 + 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 < 100 or converged); + } while (Niter < Nitermax and !converged); if (converged) { - Tvec betaDOWN = get_betaDOWN(gcov); - Tvec gammaijUP = get_gammaijUP(gcov); - Tvec v3d = alpha * momentum[id_a] / (w * lorentz_factor) - betaDOWN; + lorentz_factor = get_lorentz_factor(momentum[id_a], gcov); + + Tscal rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; + 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++) { vel[id_a][i] = sycl::dot(gammaijUP[i + 1], v3d); } } else { logger::err_ln( - "GRSPH", - "the enthalpy iterator is not converged after", - Niter, - "iterations"); + "GRSPH", "the walpy iterator is not converged after", Niter, "iterations"); } }); } From 6f4bd2f02b4bda21029e835467cd5f3962398634 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Thu, 12 Mar 2026 17:45:22 +0900 Subject: [PATCH 10/13] debugged --- src/shammodels/sph/src/modules/ConsToPrim.cpp | 50 ++++++++++++------- 1 file changed, 31 insertions(+), 19 deletions(-) diff --git a/src/shammodels/sph/src/modules/ConsToPrim.cpp b/src/shammodels/sph/src/modules/ConsToPrim.cpp index a1a6c9b23f..38a294d438 100644 --- a/src/shammodels/sph/src/modules/ConsToPrim.cpp +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -65,7 +65,7 @@ namespace shammodels::sph::modules { sham::DDMultiRef{rhostar, momentum, K}, sham::DDMultiRef{rho, vel, u, P}, thread_counts, - [gamma = this->gamma]( + [gamma = this->gamma, gcov = edges.gcov]( u32 id_a, Tscal *__restrict rhostar, Tvec *__restrict momentum, @@ -80,7 +80,6 @@ namespace shammodels::sph::modules { // get metric quantities const Tscal sqrt_g = get_sqrtg(gcov); const Tscal inv_sqrt_g = 1. / sqrt_g; - const Tscal sqrt_gamma = get_sqrt_gamma(gcov); const Tscal alpha = get_alpha(gcov); const Tscal sqrt_gamma_inv = alpha * inv_sqrt_g; const Tvec betaUP = get_betaUP(gcov); @@ -88,18 +87,19 @@ namespace shammodels::sph::modules { const std::mdspan, Layout, Accessor> gammaijUP = get_gammaijUP(gcov); - Tscal pp = 0; + 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++) { - pp += momentum[id_a][i] * sycl::dot(gammaijUP[i + 1], momentum[id_a]); + 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], gcov); + Tscal lorentz_factor = get_lorentz_factor(momentum[id_a], w, gcov); - Tscal rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; - const Tscal gamfac = gamma / (gamma - 1.); - Tscal w = 1 + gamfac * P[id_a] / rho[id_a]; // initial guess - - Tscal Kstar = K[id_a] + 0.5 * pp / rhostar[id_a]; + // Tscal Kstar = K[id_a] + 0.5 * pp / rhostar[id_a]; bool converged = false; constexpr Tscal tol = Tscal(1e-12); @@ -111,10 +111,10 @@ namespace shammodels::sph::modules { do { const Tscal w_old = w; - lorentz_factor = get_lorentz_factor(momentum[id_a], gcov); + lorentz_factor = get_lorentz_factor(momentum[id_a], w, gcov); // eq 97 - Tscal rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; + rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; // eq 62 P[id_a] = K[id_a] * sycl::pow(rho_a, gamma); @@ -122,7 +122,7 @@ namespace shammodels::sph::modules { // eq B4 const Tscal f = (1. + gamfac * P[id_a] / rho_a) - w_old; - // eq B5 + // eq B5 ... maybe should use B6 const Tscal df = -1 + gamfac * (1. @@ -137,19 +137,31 @@ namespace shammodels::sph::modules { } while (Niter < Nitermax and !converged); if (converged) { - lorentz_factor = get_lorentz_factor(momentum[id_a], gcov); - - Tscal rho_a = rhostar[id_a] * sqrt_gamma_inv / lorentz_factor; - P[id_a] = K[id_a] * sycl::pow(rho_a, gamma); + // 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++) { - vel[id_a][i] = sycl::dot(gammaijUP[i + 1], v3d); + 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 walpy iterator is not converged after", Niter, "iterations"); + "GRSPH", + "the enthalpy iterator is not converged after", + Niter, + "iterations"); } }); } From d38d045a350b9d54ab9a1c6fadd4687982f8025f Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Thu, 12 Mar 2026 17:46:09 +0900 Subject: [PATCH 11/13] utils update --- src/shamphys/include/shamphys/GRUtils.hpp | 89 ++++++++++++++++------- 1 file changed, 61 insertions(+), 28 deletions(-) diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index 01652eae0c..b5f739b3bf 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); @@ -70,32 +68,57 @@ namespace shamphys { Tvec a, Tvec b, std::mdspan, Layout2, Accessor2> gcov) { - return 0; - // TODO + + 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 get_gammaijUP( - std::mdspan, Layout2, Accessor2> gcov) { + 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 + 1, j + 1) * 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, 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 +134,29 @@ namespace shamphys { return mat_det_33(gamma_ij); } - inline static constexpr Tscal get_U0( - Tvec vxyz, std::mdspan, Layout2, Accessor2> gcov) { + 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)); + } - // U0 = gamma / alpha - Tscal alpha = get_alpha(gcov); - Tscal gamma = get_gamma(gcov); - return gamma / alpha; + inline static constexpr Tscal get_U0( + Tvec momentum, + Tscal enthalpy, + std::mdspan, Layout2, Accessor2> gcon) { + + // 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; From 2a4ae2237ed21c559be06d0e2679508ff98c402e Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Fri, 13 Mar 2026 16:30:18 +0900 Subject: [PATCH 12/13] add metric as edge --- .../sph/include/shammodels/sph/modules/ConsToPrim.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp index 33008bf6a7..d667e20e93 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp @@ -33,6 +33,8 @@ namespace shammodels::sph::modules { 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 */ \ @@ -40,6 +42,7 @@ namespace shammodels::sph::modules { X_RO(shamrock::solvergraph::IFieldSpan, spans_rhostar) \ X_RO(shamrock::solvergraph::IFieldSpan, spans_momentum) \ X_RO(shamrock::solvergraph::IFieldSpan, spans_K) \ + X_RO(GcovEdge, gcov) \ \ /* outputs */ \ X_RW(shamrock::solvergraph::IFieldSpan, spans_rho) \ From dcafee7873968680db4057d93ef668ae2b1ca21b Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Fri, 13 Mar 2026 16:46:56 +0900 Subject: [PATCH 13/13] add the contravariant metric to edges --- .../sph/include/shammodels/sph/modules/ConsToPrim.hpp | 5 +++-- src/shammodels/sph/src/modules/ConsToPrim.cpp | 9 +++++---- src/shamphys/include/shamphys/GRUtils.hpp | 11 ++++++----- 3 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp index d667e20e93..d552e067a6 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp @@ -12,8 +12,8 @@ /** * @file ConsToPrim.hpp * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) - * @brief get conserved variables (rho*, momentum, entropy) from primitive variables (rho, vel, - * internal energy) + @brief get primitive variables (rho, vel, internal energy) from conserved variables (rho*, + momentum, entropy) */ #include "shambackends/vec.hpp" @@ -42,6 +42,7 @@ namespace shammodels::sph::modules { 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 */ \ diff --git a/src/shammodels/sph/src/modules/ConsToPrim.cpp b/src/shammodels/sph/src/modules/ConsToPrim.cpp index 38a294d438..f50aadca92 100644 --- a/src/shammodels/sph/src/modules/ConsToPrim.cpp +++ b/src/shammodels/sph/src/modules/ConsToPrim.cpp @@ -57,6 +57,7 @@ namespace shammodels::sph::modules { auto &P = edges.spans_P.get_spans(); auto &gcov = edges.gcov; + auto &gcon = edges.gcon; auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); @@ -78,12 +79,12 @@ namespace shammodels::sph::modules { // on patch, no need of neighbours // get metric quantities - const Tscal sqrt_g = get_sqrtg(gcov); + const Tscal sqrt_g = get_sqrt_g(vel[id_a], gcov); const Tscal inv_sqrt_g = 1. / sqrt_g; - const Tscal alpha = get_alpha(gcov); + const Tscal alpha = get_alpha(gcon); const Tscal sqrt_gamma_inv = alpha * inv_sqrt_g; - const Tvec betaUP = get_betaUP(gcov); - const Tvec betaDOWN = get_betaDOWN(gcov); // careul cov con + const Tvec betaUP = get_betaUP(gcon); + const Tvec betaDOWN = get_betaDOWN(gcon); const std::mdspan, Layout, Accessor> gammaijUP = get_gammaijUP(gcov); diff --git a/src/shamphys/include/shamphys/GRUtils.hpp b/src/shamphys/include/shamphys/GRUtils.hpp index b5f739b3bf..ea047a1dcb 100644 --- a/src/shamphys/include/shamphys/GRUtils.hpp +++ b/src/shamphys/include/shamphys/GRUtils.hpp @@ -65,8 +65,8 @@ namespace shamphys { } inline static constexpr Tscal GR_dot( - Tvec a, - Tvec b, + std::mdspan, Layout2, Accessor2> a, + std::mdspan, Layout2, Accessor2> b, std::mdspan, Layout2, Accessor2> gcov) { Tscal result = 0; @@ -86,7 +86,7 @@ namespace shamphys { Tscal result = 0; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { - result += gamma_ij(i + 1, j + 1) * a[i] * b[j]; + result += gamma_ij(i, j) * a[i] * b[j]; } } return result; @@ -145,11 +145,12 @@ namespace shamphys { 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; }