[SPH] add CFLs to dust TVA solver

examples/sph/run_dustysettle_tvi.py

@@ -797,11 +797,14 @@ def setup_model():
         beta_AV=av_beta_AV,
     )
 
-    cfg.set_dust_mode_monofluid_tvi(nvar=ndust)
+    cfg.set_dust_mode_monofluid_tvi(
+        nvar=ndust, C_1_fluid=0.1, C_delta_v=1.0, cfl_density_threshold=1e-50
+    )
     cfg.set_dust_drag_epstein(gamma, mrn_distribution.grain_size, mrn_distribution.rho_grains)
     cfg.add_ext_force_vertical_disc_potential(central_mass=1, R0=1)
     cfg.add_ext_force_velocity_dissipation(eta=vel_dissipation_eta)
     cfg.set_two_stage_search(False)
+    cfg.set_show_cfl_detail(True)
     cfg.set_boundary_periodic()
     cfg.set_units(codeu)
     cfg.set_eos_isothermal(cs)

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

@@ -118,6 +118,10 @@ namespace shammodels::sph {
         struct MonofluidTVI {
             u32 ndust;
             bool pure_diffusion_mode = false;
+
+            Tscal C_1_fluid             = 0.1;
+            Tscal C_delta_v             = 1.0;
+            Tscal cfl_density_threshold = shambase::get_epsilon<Tscal>();
         };
 
         struct MonofluidComplete {
@@ -130,8 +134,14 @@ namespace shammodels::sph {
         Variant current_mode = None{};
 
         inline void set_none() { current_mode = None{}; }
-        inline void set_monofluid_tvi(u32 nvar, bool pure_diffusion_mode = false) {
-            current_mode = MonofluidTVI{nvar, pure_diffusion_mode};
+        inline void set_monofluid_tvi(
+            u32 nvar,
+            bool pure_diffusion_mode    = false,
+            Tscal C_1_fluid             = 0.1,
+            Tscal C_delta_v             = 1.0,
+            Tscal cfl_density_threshold = shambase::get_epsilon<Tscal>()) {
+            current_mode = MonofluidTVI{
+                nvar, pure_diffusion_mode, C_1_fluid, C_delta_v, cfl_density_threshold};
         }
         inline void set_monofluid_complete(u32 nvar) { current_mode = MonofluidComplete{nvar}; }
 
@@ -141,14 +151,21 @@ namespace shammodels::sph {
             return bool(std::get_if<MonofluidComplete>(&current_mode));
         }
 
+        inline MonofluidTVI &get_monofluid_tvi() {
+            return shambase::get_check_ref(std::get_if<MonofluidTVI>(&current_mode));
+        }
+
         inline void mode_to_json(nlohmann::json &j) const {
             if (const None *cfg = std::get_if<None>(&current_mode)) {
                 j = {{"type", "none"}};
             } else if (const MonofluidTVI *cfg = std::get_if<MonofluidTVI>(&current_mode)) {
                 j
                     = {{"type", "monofluid_tvi"},
                        {"ndust", cfg->ndust},
-                       {"pure_diffusion_mode", cfg->pure_diffusion_mode}};
+                       {"pure_diffusion_mode", cfg->pure_diffusion_mode},
+                       {"C_1_fluid", cfg->C_1_fluid},
+                       {"C_delta_v", cfg->C_delta_v},
+                       {"cfl_density_threshold", cfg->cfl_density_threshold}};
             } else if (
                 const MonofluidComplete *cfg = std::get_if<MonofluidComplete>(&current_mode)) {
                 j = {{"type", "monofluid_complete"}, {"ndust", cfg->ndust}};
@@ -163,7 +180,11 @@ namespace shammodels::sph {
                 set_none();
             } else if (type == "monofluid_tvi") {
                 set_monofluid_tvi(
-                    j.at("ndust").get<u32>(), j.at("pure_diffusion_mode").get<bool>());
+                    j.at("ndust").get<u32>(),
+                    j.at("pure_diffusion_mode").get<bool>(),
+                    j.at("C_1_fluid").get<Tscal>(),
+                    j.at("C_delta_v").get<Tscal>(),
+                    j.at("cfl_density_threshold").get<Tscal>());
             } else if (type == "monofluid_complete") {
                 set_monofluid_complete(j.at("ndust").get<u32>());
             } else {

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

@@ -0,0 +1,125 @@
+// -------------------------------------------------------//
+//
+// 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 ComputeCFLDust1Fluid.hpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief
+ *
+ */
+
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shammodels/sph/math/density.hpp"
+#include "shamrock/solvergraph/IFieldSpan.hpp"
+#include "shamrock/solvergraph/INode.hpp"
+#include "shamrock/solvergraph/Indexes.hpp"
+#include "shamrock/solvergraph/ScalarEdge.hpp"
+#include "shamsys/NodeInstance.hpp"
+
+#define NODE_EDGES(X_RO, X_RW)                                                                     \
+    X_RO(shamrock::solvergraph::Indexes<u32>, part_counts)                                         \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, C_1_fluid)                                      \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, pmass)                                          \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, hfactd)                                         \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, hpart)                                          \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, soundspeed)                                     \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, s_j)                                            \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, Ts_j)                                           \
+    X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, cfl_dt)
+
+template<class Tvec>
+class ComputeCFLDust1Fluid : public shamrock::solvergraph::INode {
+
+    using Tscal = shambase::VecComponent<Tvec>;
+
+    u32 nbins;
+
+    public:
+    ComputeCFLDust1Fluid(u32 nbins) : nbins(nbins) {}
+
+    EXPAND_NODE_EDGES(NODE_EDGES)
+
+    inline void _impl_evaluate_internal() {
+        auto edges = get_edges();
+
+        auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
+
+        Tscal C_1_fluid = edges.C_1_fluid.value;
+        Tscal pmass     = edges.pmass.value;
+        Tscal hfactd    = edges.hfactd.value;
+
+        sham::distributed_data_kernel_call(
+            dev_sched,
+            sham::DDMultiRef{
+                edges.hpart.get_spans(),
+                edges.soundspeed.get_spans(),
+                edges.s_j.get_spans(),
+                edges.Ts_j.get_spans()},
+            sham::DDMultiRef{edges.cfl_dt.get_spans()},
+            edges.part_counts.indexes,
+            [C_1_fluid, pmass, hfactd, nbins = this->nbins](
+                u32 id_a,
+                const Tscal *hpart,
+                const Tscal *soundspeed,
+                const Tscal *s_j,
+                const Tscal *Ts_j,
+                Tscal *cfl_dt) {
+                u32 id_a_d = id_a * nbins;
+
+                Tscal h_a   = hpart[id_a];
+                Tscal rho_a = shamrock::sph::rho_h(pmass, h_a, hfactd);
+
+                Tscal cs_a  = soundspeed[id_a];
+                Tscal cs2_a = cs_a * cs_a;
+
+                auto rho_dust = [&](int j) {
+                    auto tmp = s_j[id_a_d + j];
+                    return tmp * tmp;
+                };
+
+                auto epsilon_j = [&](int j) {
+                    return rho_dust(j) / rho_a;
+                };
+
+                Tscal sum_eps = 0;
+                for (int j = 0; j < nbins; j++) {
+                    sum_eps += epsilon_j(j);
+                }
+
+                Tscal cs_tilde_2_a = cs2_a * (1 - sum_eps);
+
+                Tscal cs4_over_h2 = cs2_a * cs2_a / (h_a * h_a);
+
+                Tscal cfl_tmp = std::numeric_limits<Tscal>::infinity();
+
+                for (int j = 0; j < nbins; j++) {
+                    Tscal eps_j_a  = epsilon_j(j);
+                    Tscal eps2_j_a = eps_j_a * eps_j_a;
+
+                    Tscal Ts_j_a  = Ts_j[id_a_d + j];
+                    Tscal Ts2_j_a = Ts_j_a * Ts_j_a;
+
+                    Tscal dt_j = h_a / sycl::sqrt(cs_tilde_2_a + Ts2_j_a * eps2_j_a * cs4_over_h2);
+                    cfl_tmp    = sycl::min(cfl_tmp, dt_j);
+                }
+
+                cfl_tmp *= C_1_fluid;
+
+                cfl_dt[id_a] = sycl::min(cfl_dt[id_a], cfl_tmp);
+            });
+    }
+
+    inline virtual std::string _impl_get_label() const { return "ComputeCFLDust1Fluid"; };
+
+    inline virtual std::string _impl_get_tex() const { return "C_{1,fluid}"; };
+};
+
+#undef NODE_EDGES

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

@@ -0,0 +1,106 @@
+// -------------------------------------------------------//
+//
+// 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 ComputeCFLDustDeltav.hpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief
+ *
+ */
+
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shammodels/sph/math/density.hpp"
+#include "shamrock/solvergraph/IFieldSpan.hpp"
+#include "shamrock/solvergraph/INode.hpp"
+#include "shamrock/solvergraph/Indexes.hpp"
+#include "shamrock/solvergraph/ScalarEdge.hpp"
+#include "shamsys/NodeInstance.hpp"
+
+#define NODE_EDGES(X_RO, X_RW)                                                                     \
+    X_RO(shamrock::solvergraph::Indexes<u32>, part_counts)                                         \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, C_delta_v)                                      \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, cfl_density_threshold)                          \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, pmass)                                          \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, hfactd)                                         \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, hpart)                                          \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, s_j)                                            \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tvec>, delta_v)                                         \
+    X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, cfl_dt)
+
+template<class Tvec>
+class ComputeCFLDustDeltav : public shamrock::solvergraph::INode {
+
+    using Tscal = shambase::VecComponent<Tvec>;
+
+    u32 nbins;
+
+    public:
+    ComputeCFLDustDeltav(u32 nbins) : nbins(nbins) {}
+
+    EXPAND_NODE_EDGES(NODE_EDGES)
+
+    inline void _impl_evaluate_internal() {
+        auto edges = get_edges();
+
+        auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
+
+        Tscal C_delta_v             = edges.C_delta_v.value;
+        Tscal cfl_density_threshold = edges.cfl_density_threshold.value;
+
+        Tscal pmass  = edges.pmass.value;
+        Tscal hfactd = edges.hfactd.value;
+
+        sham::distributed_data_kernel_call(
+            dev_sched,
+            sham::DDMultiRef{
+                edges.hpart.get_spans(), edges.delta_v.get_spans(), edges.s_j.get_spans()},
+            sham::DDMultiRef{edges.cfl_dt.get_spans()},
+            edges.part_counts.indexes,
+            [C_delta_v, cfl_density_threshold, pmass, hfactd, nbins = this->nbins](
+                u32 id_a,
+                const Tscal *hpart,
+                const Tvec *delta_v,
+                const Tscal *s_j,
+                Tscal *cfl_dt) {
+                u32 id_a_d = id_a * nbins;
+
+                Tscal h_a   = hpart[id_a];
+                Tscal rho_a = shamrock::sph::rho_h(pmass, h_a, hfactd);
+
+                auto rho_dust = [&](int j) {
+                    auto tmp = s_j[id_a_d + j];
+                    return tmp * tmp;
+                };
+
+                Tscal cfl_tmp = std::numeric_limits<Tscal>::infinity();
+
+                for (int j = 0; j < nbins; j++) {
+                    Tvec delta_v_j_a       = delta_v[id_a_d + j];
+                    Tscal delta_v_j_a_norm = sycl::length(delta_v_j_a);
+
+                    Tscal rho_d_j_a = rho_dust(j);
+                    if (rho_d_j_a > cfl_density_threshold && delta_v_j_a_norm > 0) {
+                        cfl_tmp = sycl::min(cfl_tmp, h_a / delta_v_j_a_norm);
+                    }
+                }
+
+                cfl_tmp *= C_delta_v;
+
+                cfl_dt[id_a] = sycl::min(cfl_dt[id_a], cfl_tmp);
+            });
+    }
+
+    inline virtual std::string _impl_get_label() const { return "ComputeCFLDustDeltav"; };
+
+    inline virtual std::string _impl_get_tex() const { return "C_{\\Delta_v}"; };
+};
+
+#undef NODE_EDGES