[GRSPH] Conservative to primitive variable swap module

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)

src/shammodels/sph/include/shammodels/sph/modules/ConsToPrim.hpp

@@ -0,0 +1,76 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
+// 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 Tvec>
+    class NodeConsToPrim : public shamrock::solvergraph::INode {
+        using Tscal = shambase::VecComponent<Tvec>;
+        u32 block_size;
+        Tscal gamma;
+
+        public:
+        NodeConsToPrim(u32 block_size, Tscal gamma) : block_size(block_size), gamma(gamma) {}
+
+        struct Edges {
+            const shamrock::solvergraph::Indexes<u32> &sizes;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &spans_rhostar;
+            const shamrock::solvergraph::IFieldSpan<Tvec> &spans_momentum;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &spans_K;
+            shamrock::solvergraph::IFieldSpan<Tscal> &spans_rho;
+            shamrock::solvergraph::IFieldSpan<Tvec> &spans_vel;
+            shamrock::solvergraph::IFieldSpan<Tscal> &spans_u;
+            shamrock::solvergraph::IFieldSpan<Tscal> &spans_P;
+        };
+
+        inline void set_edges(
+            std::shared_ptr<shamrock::solvergraph::Indexes<u32>> sizes,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> spans_rhostar,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tvec>> spans_momentum,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> spans_K,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> spans_rho,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tvec>> spans_vel,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> spans_u,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> 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<shamrock::solvergraph::Indexes<u32>>(0),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(1),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tvec>>(2),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(3),
+                get_rw_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(0),
+                get_rw_edge<shamrock::solvergraph::IFieldSpan<Tvec>>(1),
+                get_rw_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(2),
+                get_rw_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(3)};
+        }
+
+        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

src/shammodels/sph/src/modules/ConsToPrim.cpp

@@ -0,0 +1,129 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
+// 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 "shammath/riemann.hpp"
+#include "shamrock/patch/PatchDataField.hpp"
+#include "shamsys/NodeInstance.hpp"
+
+namespace {
+
+    template<class Tvec>
+    struct KernelConsToPrim {
+        using Tscal = shambase::VecComponent<Tvec>;
+    };
+
+} // namespace
+
+namespace shammodels::basegodunov::modules {
+
+    template<class Tvec>
+    void NodeConsToPrim<Tvec>::_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 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](
+                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
+                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];
+                bool converged = false;
+                // 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 = 1;
+
+                    converged = sycl::fabs(new_w - w) < 1e-6;
+                    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<class Tvec>
+    std::string NodeConsToPrim<Tvec>::_impl_get_tex() const {
+
+        return "TODO";
+    }
+
+} // namespace shammodels::basegodunov::modules
+
+template class shammodels::basegodunov::modules::NodeConsToPrim<f64_3>;