[GSPH] Disc physics in GSPH

src/shammodels/gsph/CMakeLists.txt

@@ -19,6 +19,11 @@ set(Sources
     src/modules/GSPHGhostHandler.cpp
     src/modules/ComputeLoadBalanceValue.cpp
     src/modules/io/VTKDump.cpp
+    src/modules/GSPHSetup.cpp
+    src/modules/GeneratorMCDisc.cpp
+    src/modules/ComputeLoadBalanceValue.cpp
+    src/modules/SinkParticlesUpdate.cpp
+    src/modules/ExternalForces.cpp
 )
 
 if(SHAMROCK_USE_SHARED_LIB)

src/shammodels/gsph/include/shammodels/gsph/Model.hpp

@@ -35,6 +35,7 @@
 #include "shamcomm/logs.hpp"
 #include "shammodels/common/setup/generators.hpp"
 #include "shammodels/gsph/Solver.hpp"
+#include "shammodels/gsph/modules/GSPHSetup.hpp"
 #include "shammodels/sph/math/density.hpp"
 #include "shamrock/io/ShamrockDump.hpp"
 #include "shamrock/patch/PatchDataLayer.hpp"
@@ -132,6 +133,22 @@ namespace shammodels::gsph {
         void add_cube_fcc_3d(Tscal dr, std::pair<Tvec, Tvec> _box);
         void add_cube_hcp_3d(Tscal dr, std::pair<Tvec, Tvec> _box);
 
+        inline std::unique_ptr<modules::GSPHSetup<Tvec, SPHKernel>> get_setup() {
+            return std::make_unique<modules::GSPHSetup<Tvec, SPHKernel>>(
+                ctx, solver.solver_config, solver.storage);
+        }
+
+        inline void add_sink(Tscal mass, Tvec pos, Tvec velocity, Tscal accretion_radius) {
+            if (solver.storage.sinks.is_empty()) {
+                solver.storage.sinks.set({});
+            }
+
+            shamlog_debug_ln("SPH", "add sink :", mass, pos, velocity, accretion_radius);
+
+            solver.storage.sinks.get().push_back(
+                {pos, velocity, {}, {}, mass, {}, accretion_radius});
+        }
+
         ////////////////////////////////////////////////////////////////////////////////////////////
         // Field manipulation
         ////////////////////////////////////////////////////////////////////////////////////////////
@@ -412,6 +429,12 @@ namespace shammodels::gsph {
             nlohmann::json metadata;
             metadata["solver_config"] = solver.solver_config;
 
+            if (solver.storage.sinks.is_empty()) {
+                metadata["sinks"] = nlohmann::json{};
+            } else {
+                metadata["sinks"] = solver.storage.sinks.get();
+            }
+
             shamrock::write_shamrock_dump(
                 fname, metadata.dump(4), shambase::get_check_ref(ctx.sched));
         }

src/shammodels/gsph/include/shammodels/gsph/Solver.hpp

@@ -224,10 +224,17 @@ namespace shammodels::gsph {
          * - Courant condition: dt_cour = C_cour * h / vsig
          * - Force condition: dt_force = C_force * sqrt(h / |a|)
          *
-         * @return Minimum CFL timestep across all particles
+         * @return Minimum CFL timestep across all sinks
          */
         Tscal compute_dt_cfl();
 
+        /**
+         * @brief Compute sink timestep constraint
+         *
+         * @return Minimum CFL timestep across all particles
+         */
+        Tscal compute_sink_cfl();
+
         bool apply_corrector(Tscal dt, u64 Npart_all);
 
         void update_sync_load_values();

src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp

@@ -73,6 +73,7 @@ namespace shammodels::gsph {
     struct CFLConfig {
         Tscal cfl_cour  = 0.3;  ///< CFL condition for the courant factor
         Tscal cfl_force = 0.25; ///< CFL condition for the force
+        Tscal eta_sink  = 0.05;
     };
 
     struct SmoothingLengthConfig {
@@ -133,6 +134,16 @@ struct shammodels::gsph::SolverConfig {
         }
     }
 
+    inline Tscal get_constant_c() const {
+        if (!unit_sys) {
+            ON_RANK_0(logger::warn_ln("gsph::Config", "the unit system is not set"));
+            shamunits::Constants<Tscal> ctes{shamunits::UnitSystem<Tscal>{}};
+            return ctes.c();
+        } else {
+            return shamunits::Constants<Tscal>{*unit_sys}.c();
+        }
+    }
+
     //////////////////////////////////////////////////////////////////////////////////////////////
     // Units Config (END)
     //////////////////////////////////////////////////////////////////////////////////////////////
@@ -234,6 +245,29 @@ struct shammodels::gsph::SolverConfig {
 
     inline void set_eos_isothermal(Tscal cs) { eos_config.set_isothermal(cs); }
 
+    /**
+     * `@brief` Set the EOS configuration to a locally isothermal equation of state from Farris 2014
+     *
+     * `@param` h_over_r Disc aspect ratio
+     */
+    inline void set_eos_locally_isothermalFA2014(Tscal h_over_r) {
+        eos_config.set_locally_isothermalFA2014(h_over_r);
+    }
+
+    /**
+     * @brief Set the EOS configuration to a locally isothermal equation of state from Farris 2014
+     * extended to q != 1/2
+     *
+     * @param cs0 Soundspeed at the reference radius
+     * @param q Power exponent of the soundspeed profile
+     * @param r0 Reference radius
+     * @param n_sinks Number of sinks to consider for the equation of state
+     */
+    inline void set_eos_locally_isothermalFA2014_extended(
+        Tscal cs0, Tscal q, Tscal r0, u32 n_sinks) {
+        eos_config.set_locally_isothermalFA2014_extended(cs0, q, r0, n_sinks);
+    }
+
     //////////////////////////////////////////////////////////////////////////////////////////////
     // EOS Config (END)
     //////////////////////////////////////////////////////////////////////////////////////////////

src/shammodels/gsph/include/shammodels/gsph/config/FieldNames.hpp

@@ -45,6 +45,9 @@ namespace shammodels::gsph::names {
         /// 3-acceleration field
         inline constexpr const char *axyz = "axyz";
 
+        /// 3-acceleration field due to external forces
+        inline constexpr const char *axyz_ext = "axyz_ext";
+
         /// Specific internal energy u
         inline constexpr const char *uint = "uint";
 

src/shammodels/gsph/include/shammodels/gsph/modules/ExternalForces.hpp

@@ -0,0 +1,67 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <[email protected]>
+// 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 ExternalForces.hpp
+ * @author Timothée David--Cléris ([email protected])
+ * @author Yona Lapeyre ([email protected])
+ * @brief
+ *
+ */
+
+#include "shambackends/typeAliasVec.hpp"
+#include "shambackends/vec.hpp"
+#include "shammodels/gsph/SolverConfig.hpp"
+#include "shammodels/gsph/modules/SolverStorage.hpp"
+#include "shamrock/scheduler/ShamrockCtx.hpp"
+
+namespace shammodels::gsph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class ExternalForces {
+        public:
+        using Tscal              = shambase::VecComponent<Tvec>;
+        static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
+        using Kernel             = SPHKernel<Tscal>;
+
+        using Config  = SolverConfig<Tvec, SPHKernel>;
+        using Storage = SolverStorage<Tvec, u32>;
+
+        ShamrockCtx &context;
+        Config &solver_config;
+        Storage &storage;
+
+        ExternalForces(ShamrockCtx &context, Config &solver_config, Storage &storage)
+            : context(context), solver_config(solver_config), storage(storage) {}
+
+        /**
+         * @brief is ran once per timestep, it computes the forces that are independant of velocity
+         *
+         */
+        void compute_ext_forces_indep_v();
+
+        /**
+         * @brief add external forces to the particle acceleration, note that forces dependant on
+         * velocity shlould be added here
+         *
+         */
+        void add_ext_forces();
+
+        void point_mass_accrete_particles();
+
+        private:
+        using SolverConfigExtForce = typename Config::ExtForceConfig;
+        using EF_PointMass         = typename SolverConfigExtForce::PointMass;
+
+        inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); }
+    };
+
+} // namespace shammodels::gsph::modules

src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp

@@ -0,0 +1,89 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <[email protected]>
+// 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 GSPHSetup.hpp
+ * @author Yona Lapeyre ([email protected])
+ * @brief
+ *
+ */
+
+#include "shambackends/typeAliasVec.hpp"
+#include "shambackends/vec.hpp"
+#include "shammodels/gsph/SolverConfig.hpp"
+#include "shammodels/gsph/modules/SolverStorage.hpp"
+#include "shammodels/gsph/modules/setup/IGSPHSetupNode.hpp"
+#include "shamrock/scheduler/ShamrockCtx.hpp"
+#include <memory>
+
+namespace shammodels::gsph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class GSPHSetup {
+        public:
+        using Tscal              = shambase::VecComponent<Tvec>;
+        static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
+        using Kernel             = SPHKernel<Tscal>;
+
+        using Config  = SolverConfig<Tvec, SPHKernel>;
+        using Storage = SolverStorage<Tvec, u32>;
+
+        ShamrockCtx &context;
+        Config &solver_config;
+        Storage &storage;
+
+        GSPHSetup(ShamrockCtx &context, Config &solver_config, Storage &storage)
+            : context(context), solver_config(solver_config), storage(storage) {}
+
+        void apply_setup(SetupNodePtr setup, std::optional<u32> insert_step = std::nullopt);
+
+        std::shared_ptr<IGSPHSetupNode> make_generator_disc_mc(
+            Tscal part_mass,
+            Tscal disc_mass,
+            Tscal r_in,
+            Tscal r_out,
+            std::function<Tscal(Tscal)> sigma_profile,
+            std::function<Tscal(Tscal)> H_profile,
+            std::function<Tvec(Tvec)> vel_profile,
+            std::function<Tscal(Tvec)> cs_profile,
+            std::mt19937_64 eng,
+            Tscal init_h_factor);
+
+        std::shared_ptr<IGSPHSetupNode> make_generator_from_context(ShamrockCtx &context_other);
+
+        std::shared_ptr<IGSPHSetupNode> make_combiner_add(
+            SetupNodePtr parent1, SetupNodePtr parent2);
+
+        std::shared_ptr<IGSPHSetupNode> make_modifier_warp_disc(
+            SetupNodePtr parent, Tscal Rwarp, Tscal Hwarp, Tscal inclination, Tscal posangle);
+
+        std::shared_ptr<IGSPHSetupNode> make_modifier_custom_warp(
+            SetupNodePtr parent,
+            std::function<Tscal(Tscal)> inc_profile,
+            std::function<Tscal(Tscal)> psi_profile,
+            std::function<Tvec(Tscal)> k_profile);
+
+        std::shared_ptr<IGSPHSetupNode> make_modifier_add_offset(
+            SetupNodePtr parent, Tvec offset_postion, Tvec offset_velocity);
+
+        std::shared_ptr<IGSPHSetupNode> make_modifier_filter(
+            SetupNodePtr parent, std::function<bool(Tvec)> filter);
+
+        std::shared_ptr<IGSPHSetupNode> make_modifier_split_part(
+            SetupNodePtr parent, u64 n_split, u64 seed, Tscal h_scaling);
+
+        private:
+        inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); }
+
+        u64 injected_parts = 0;
+    };
+
+} // namespace shammodels::gsph::modules

src/shammodels/gsph/include/shammodels/gsph/modules/SinkParticlesUpdate.hpp

@@ -0,0 +1,57 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <[email protected]>
+// 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 SinkParticlesUpdate.hpp
+ * @author Timothée David--Cléris ([email protected])
+ * @brief
+ *
+ */
+
+#include "shambackends/typeAliasVec.hpp"
+#include "shambackends/vec.hpp"
+#include "shammodels/gsph/SolverConfig.hpp"
+#include "shammodels/gsph/modules/SolverStorage.hpp"
+#include "shammodels/sph/SinkPartStruct.hpp"
+#include "shamrock/scheduler/ShamrockCtx.hpp"
+
+namespace shammodels::gsph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class SinkParticlesUpdate {
+        public:
+        using Tscal              = shambase::VecComponent<Tvec>;
+        static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
+        using Kernel             = SPHKernel<Tscal>;
+
+        using Config  = SolverConfig<Tvec, SPHKernel>;
+        using Storage = SolverStorage<Tvec, u32>;
+
+        ShamrockCtx &context;
+        Config &solver_config;
+        Storage &storage;
+
+        using Sink = sph::SinkParticle<Tvec>;
+
+        SinkParticlesUpdate(ShamrockCtx &context, Config &solver_config, Storage &storage)
+            : context(context), solver_config(solver_config), storage(storage) {}
+
+        void accrete_particles(Tscal dt);
+        void predictor_step(Tscal dt);
+        void compute_sph_forces();
+        void compute_ext_forces();
+        void corrector_step(Tscal dt);
+
+        private:
+        inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); }
+    };
+
+} // namespace shammodels::gsph::modules

src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp

@@ -29,6 +29,7 @@
 #include "shambackends/vec.hpp"
 #include "shammodels/gsph/modules/GSPHGhostHandler.hpp"
 #include "shammodels/gsph/solvergraph/GhostHandlerEdge.hpp"
+#include "shammodels/sph/SinkPartStruct.hpp"
 #include "shammodels/sph/solvergraph/NeighCache.hpp"
 #include "shamrock/scheduler/SerialPatchTree.hpp"
 #include "shamrock/scheduler/ShamrockCtx.hpp"
@@ -136,6 +137,8 @@ namespace shammodels::gsph {
         Component<shamrock::ComputeField<Tvec>> old_axyz;
         Component<shamrock::ComputeField<Tscal>> old_duint;
 
+        Component<std::vector<sph::SinkParticle<Tvec>>> sinks;
+
         /// Timing statistics
         struct Timings {
             f64 interface = 0;