[SPH] Add luminosity field (as its own module)

src/shammodels/sph/CMakeLists.txt

@@ -38,6 +38,7 @@ set(Sources
     src/modules/UpdateDerivs.cpp
     src/modules/ComputeLoadBalanceValue.cpp
     src/modules/ComputeOmega.cpp
+    src/modules/ComputeLuminosity.cpp
     src/modules/NeighbourCache.cpp
     src/modules/ParticleReordering.cpp
     src/modules/AnalysisSodTube.cpp

src/shammodels/sph/include/shammodels/sph/Solver.hpp

@@ -12,7 +12,7 @@
 /**
  * @file Solver.hpp
  * @author Timothée David--Cléris (tim.shamrock@proton.me)
- * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) --no git blame--
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
  * @brief
  */
 
@@ -22,6 +22,7 @@
 #include "shammodels/sph/BasicSPHGhosts.hpp"
 #include "shammodels/sph/SPHUtilities.hpp"
 #include "shammodels/sph/SolverLog.hpp"
+#include "shammodels/sph/config/AVConfig.hpp"
 #include "shammodels/sph/modules/SolverStorage.hpp"
 #include "shamrock/patch/PatchDataLayerLayout.hpp"
 #include "shamrock/scheduler/ComputeField.hpp"
@@ -61,7 +62,8 @@ namespace shammodels::sph {
         static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
         using Kernel             = SPHKernel<Tscal>;
 
-        using Config = SolverConfig<Tvec, SPHKernel>;
+        using Config   = SolverConfig<Tvec, SPHKernel>;
+        using AVConfig = AVConfig<Tvec>;
 
         using u_morton = typename Config::u_morton;
 
@@ -73,6 +75,7 @@ namespace shammodels::sph {
         SolverStorage<Tvec, u_morton> storage{};
 
         Config solver_config;
+        AVConfig av_config;
         SolverLog solve_logs;
 
         inline void init_required_fields() { solver_config.set_layout(context.get_pdl_write()); }

src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp

@@ -808,6 +808,10 @@ struct shammodels::sph::SolverConfig {
     /// @brief Whether the solver has a field for dt divB
     inline bool has_field_dtdivB() { return mhd_config.has_dtdivB_field(); }
 
+    /// @brief Whether to store luminosity
+    bool compute_luminosity = false;
+    inline void use_luminosity(bool enable) { compute_luminosity = enable; }
+
     /// Print the current status of the solver config
     inline void print_status() {
         if (shamcomm::world_rank() != 0) {

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

@@ -0,0 +1,83 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2025 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 ComputeLuminosity.hpp
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
+ * @brief
+ *
+ */
+
+#include "shambackends/typeAliasVec.hpp"
+#include "shambackends/vec.hpp"
+#include "shammodels/sph/SolverConfig.hpp"
+#include "shammodels/sph/modules/SolverStorage.hpp"
+#include "shamrock/scheduler/ShamrockCtx.hpp"
+
+namespace shammodels::sph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class NodeComputeLuminosity : public shamrock::solvergraph::INode {
+
+        using Tscal = shambase::VecComponent<Tvec>;
+
+        static constexpr Tscal kernel_radius = SPHKernel<Tscal>::Rkern;
+        Tscal part_mass;
+        Tscal alpha_u;
+
+        public:
+        NodeComputeLuminosity(Tscal part_mass, Tscal alpha_u)
+            : part_mass(part_mass), alpha_u(alpha_u) {}
+
+        struct Edges {
+            const shamrock::solvergraph::Indexes<u32> &part_counts;
+            const shammodels::sph::solvergraph::NeighCache &neigh_cache;
+            const shamrock::solvergraph::IFieldSpan<Tvec> &xyz;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &hpart;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &omega;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &uint;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &pressure;
+            shamrock::solvergraph::IFieldSpan<Tscal> &luminosity;
+        };
+
+        inline void set_edges(
+            std::shared_ptr<shamrock::solvergraph::Indexes<u32>> part_counts,
+            std::shared_ptr<shammodels::sph::solvergraph::NeighCache> neigh_cache,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tvec>> xyz,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> hpart,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> omega,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> uint,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> pressure,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> luminosity) {
+            __internal_set_ro_edges({part_counts, neigh_cache, xyz, hpart, omega, uint, pressure});
+            __internal_set_rw_edges({luminosity});
+        }
+
+        inline Edges get_edges() {
+            return Edges{
+                get_ro_edge<shamrock::solvergraph::Indexes<u32>>(0),
+                get_ro_edge<shammodels::sph::solvergraph::NeighCache>(1),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tvec>>(2),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(3),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(4),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(5),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(6),
+                get_rw_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(0)};
+        }
+
+        void _impl_evaluate_internal();
+
+        inline virtual std::string _impl_get_label() const { return "ComputeLuminosity"; };
+
+        virtual std::string _impl_get_tex() const;
+    };
+
+} // namespace shammodels::sph::modules

src/shammodels/sph/src/Solver.cpp

@@ -46,6 +46,7 @@
 #include "shammodels/sph/modules/BuildTrees.hpp"
 #include "shammodels/sph/modules/ComputeEos.hpp"
 #include "shammodels/sph/modules/ComputeLoadBalanceValue.hpp"
+#include "shammodels/sph/modules/ComputeLuminosity.hpp"
 #include "shammodels/sph/modules/ComputeNeighStats.hpp"
 #include "shammodels/sph/modules/ComputeOmega.hpp"
 #include "shammodels/sph/modules/ConservativeCheck.hpp"
@@ -2000,6 +2001,92 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
 
         update_derivs();
 
+        bool has_luminosity = solver_config.compute_luminosity;
+
+        if (has_luminosity) {
+            const u32 iluminosity = pdl.get_field_idx<Tscal>("luminosity");
+
+            shambase::get_check_ref(storage.hpart_with_ghosts)
+                .set_refs(storage.merged_xyzh.get()
+                              .template map<std::reference_wrapper<PatchDataField<Tscal>>>(
+                                  [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
+                                      return std::ref(mpdat.get_field<Tscal>(
+                                          1)); // hpart is at index 1 in merged_xyzh
+                                  }));
+
+            auto uint_with_ghost = shamrock::solvergraph::FieldRefs<Tscal>::make_shared("", "");
+
+            shambase::get_check_ref(storage.hpart_with_ghosts)
+                .set_refs(storage.merged_xyzh.get()
+                              .template map<std::reference_wrapper<PatchDataField<Tscal>>>(
+                                  [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
+                                      return std::ref(mpdat.get_field<Tscal>(1));
+                                  }));
+
+            shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::FieldRefs<Tscal>>
+                set_uint_with_ghost_refs(
+                    [&](shamrock::solvergraph::FieldRefs<Tscal> &field_uint_with_ghost_edge) {
+                        shambase::DistributedData<PatchDataLayer> &mpdats
+                            = storage.merged_patchdata_ghost.get();
+
+                        shamrock::solvergraph::DDPatchDataFieldRef<Tscal> field_uint_with_ghost_refs
+                            = {};
+
+                        scheduler().for_each_patchdata_nonempty(
+                            [&](const Patch p, PatchDataLayer &pdat) {
+                                PatchDataLayer &mpdat = mpdats.get(p.id_patch);
+
+                                auto &field = mpdat.get_field<Tscal>(iuint_interf);
+                                field_uint_with_ghost_refs.add_obj(p.id_patch, std::ref(field));
+                            });
+
+                        field_uint_with_ghost_edge.set_refs(field_uint_with_ghost_refs);
+                    });
+
+            set_uint_with_ghost_refs.set_edges(uint_with_ghost);
+            // @@@@@@@@@@@@@@@@
+
+            auto luminosity = shamrock::solvergraph::FieldRefs<Tscal>::make_shared("", "");
+
+            shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::FieldRefs<Tscal>>
+                set_luminosity_refs(
+                    [&](shamrock::solvergraph::FieldRefs<Tscal> &field_luminosity_edge) {
+                        shambase::DistributedData<PatchDataLayer> &mpdats
+                            = storage.merged_patchdata_ghost.get();
+
+                        shamrock::solvergraph::DDPatchDataFieldRef<Tscal> field_luminosity_refs
+                            = {};
+
+                        scheduler().for_each_patchdata_nonempty(
+                            [&](const Patch p, PatchDataLayer &pdat) {
+                                auto &field = pdat.get_field<Tscal>(iluminosity);
+                                field_luminosity_refs.add_obj(p.id_patch, std::ref(field));
+                            });
+                        field_luminosity_edge.set_refs(field_luminosity_refs);
+                    });
+
+            set_luminosity_refs.set_edges(luminosity);
+
+            set_uint_with_ghost_refs.evaluate();
+            set_luminosity_refs.evaluate();
+
+            Tscal alpha_u = 1;
+
+            modules::NodeComputeLuminosity<Tvec, Kern> compute_luminosity{
+                solver_config.gpart_mass, alpha_u};
+
+            compute_luminosity.set_edges(
+                storage.part_counts_with_ghost,
+                storage.neigh_cache,
+                storage.positions_with_ghosts,
+                storage.hpart_with_ghosts,
+                storage.omega,
+                uint_with_ghost,
+                storage.pressure,
+                luminosity);
+            compute_luminosity.evaluate();
+        }
+
         modules::ConservativeCheck<Tvec, Kern> cv_check(context, solver_config, storage);
         cv_check.check_conservation();
 

src/shammodels/sph/src/SolverConfig.cpp

@@ -91,6 +91,11 @@ namespace shammodels::sph {
             pdl.add_field<Tvec>("deltav", ndust);
             pdl.add_field<Tvec>("dtdeltav", ndust);
         }
+
+        if (compute_luminosity) {
+            pdl.add_field<Tscal>("luminosity", 1);
+        }
+
         if (do_MHD_debug()) {
             pdl.add_field<Tvec>("gas_pressure", 1);
             pdl.add_field<Tvec>("mag_pressure", 1);

src/shammodels/sph/src/modules/ComputeLuminosity.cpp

@@ -0,0 +1,115 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2025 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 ComputeLuminosity.cpp
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
+ * @brief
+ *
+ */
+
+#include "shambase/stacktrace.hpp"
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shammodels/sph/SPHUtilities.hpp"
+#include "shammodels/sph/math/forces.hpp"
+#include "shammodels/sph/modules/ComputeLuminosity.hpp"
+#include "shamrock/scheduler/SchedulerUtility.hpp"
+#include "shamrock/solvergraph/IFieldSpan.hpp"
+
+template<class Tvec, template<class> class SPHKernel>
+void shammodels::sph::modules::NodeComputeLuminosity<Tvec, SPHKernel>::_impl_evaluate_internal() {
+
+    __shamrock_stack_entry();
+
+    auto edges = get_edges();
+
+    auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
+
+    edges.luminosity.ensure_sizes(edges.part_counts.indexes);
+
+    sham::distributed_data_kernel_call(
+        dev_sched,
+        sham::DDMultiRef{
+            edges.xyz.get_spans(),
+            edges.hpart.get_spans(),
+            edges.omega.get_spans(),
+            edges.uint.get_spans(),
+            edges.pressure.get_spans(),
+            edges.neigh_cache.neigh_cache},
+        sham::DDMultiRef{edges.luminosity.get_spans()},
+        edges.part_counts.indexes,
+        [part_mass = this->part_mass, alpha_u = this->alpha_u](
+            u32 id_a,
+            const Tvec *r,
+            const Tscal *hpart,
+            const Tscal *omega,
+            const Tscal *uint,
+            const Tscal *pressure,
+            const auto ploop_ptrs,
+            Tscal *luminosity) {
+            shamrock::tree::ObjectCacheIterator particle_looper(ploop_ptrs);
+
+            using namespace shamrock::sph;
+            logger::raw_ln("first batch of loads");
+            Tscal h_a               = hpart[id_a];
+            Tvec xyz_a              = r[id_a];
+            const Tscal u_a         = uint[id_a];
+            const Tscal omega_a     = omega[id_a];
+            const Tscal rho_a       = rho_h(part_mass, h_a, SPHKernel<Tscal>::hfactd);
+            const Tscal P_a         = pressure[id_a];
+            Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
+            logger::raw_ln("passed first batch of loads");
+            Tscal tmp_luminosity = 0;
+
+            particle_looper.for_each_object(id_a, [&](u32 id_b) {
+                logger::raw_ln("second batch of loads");
+                const Tscal u_b = uint[id_b];
+                logger::raw_ln("PB u");
+                const Tscal h_b = hpart[id_b];
+                logger::raw_ln("PB h");
+                const Tscal omega_b = omega[id_b];
+                logger::raw_ln("PB omega");
+                const Tscal P_b = pressure[id_b];
+                logger::raw_ln("PB pressure");
+                const Tscal rho_b = rho_h(part_mass, h_b, SPHKernel<Tscal>::hfactd);
+                Tvec dr           = xyz_a - r[id_b];
+                Tscal rab2        = sycl::dot(dr, dr);
+                Tscal rab         = sycl::sqrt(rab2);
+                logger::raw_ln("passed second batch of loads");
+                Tscal vsigu = vsig_u(P_a, P_b, rho_a, rho_b);
+                Tscal Fab_a = SPHKernel<Tscal>::dW_3d(rab, h_a);
+                Tscal Fab_b = SPHKernel<Tscal>::dW_3d(rab, h_b);
+
+                tmp_luminosity += lambda_shock_conductivity(
+                    part_mass,
+                    alpha_u,
+                    vsigu,
+                    u_a - u_b,
+                    Fab_a * omega_a_rho_a_inv,
+                    Fab_b / (rho_b * omega_b));
+            });
+
+            luminosity[id_a] = tmp_luminosity;
+        });
+}
+
+template<class Tvec, template<class> class SPHKernel>
+std::string shammodels::sph::modules::NodeComputeLuminosity<Tvec, SPHKernel>::_impl_get_tex()
+    const {
+    return "TODO";
+}
+
+using namespace shammath;
+template class shammodels::sph::modules::NodeComputeLuminosity<f64_3, M4>;
+template class shammodels::sph::modules::NodeComputeLuminosity<f64_3, M6>;
+template class shammodels::sph::modules::NodeComputeLuminosity<f64_3, M8>;
+
+template class shammodels::sph::modules::NodeComputeLuminosity<f64_3, C2>;
+template class shammodels::sph::modules::NodeComputeLuminosity<f64_3, C4>;
+template class shammodels::sph::modules::NodeComputeLuminosity<f64_3, C6>;