diff --git a/examples/sph/run_dusty.py b/examples/sph/run_dusty.py
new file mode 100644
index 000000000..96e0a57a5
--- /dev/null
+++ b/examples/sph/run_dusty.py
@@ -0,0 +1,170 @@
+"""
+Testing Sod tube with SPH
+=========================
+
+CI test for Sod tube with SPH
+"""
+
+import matplotlib.pyplot as plt
+
+import shamrock
+
+shamrock.enable_experimental_features()
+import numpy as np
+
+gamma = 1.4
+rho = 1
+
+
+def func_rho_t(r):
+    return rho
+
+
+def func_rho_d_j(r, idust):
+    r_ = np.sqrt(r[0] ** 2 + r[1] ** 2 + r[2] ** 2)
+    return (0.1 / ndust) * max(1 - (r_ / rc) ** 2, 0)
+
+
+def func_rho_g(r):
+    return rho - sum([func_rho_d_j(r, i) for i in range(ndust)])
+
+
+def func_s_j(r, idust):
+    rho_t = func_rho_t(r)
+    rho_d_j = [func_rho_d_j(r, i) for i in range(ndust)]
+    eps_j = rho_d_j[idust] / rho_t
+    return np.sqrt(rho_t * eps_j)
+
+
+cs_g = 1
+
+
+def uint_g(r):
+    rho_g = func_rho_g(r)
+    P = rho_g * cs_g * cs_g / gamma
+    return P / ((gamma - 1) * rho_g)
+
+
+ndust = 1
+rc = 0.25
+stopping_times = [1e-1]
+
+
+bmin = (-0.6, -0.6, -0.6)
+bmax = (0.6, 0.6, 0.6)
+
+N_target = 1e4
+scheduler_split_val = int(2e7)
+scheduler_merge_val = int(1)
+
+xm, ym, zm = bmin
+xM, yM, zM = bmax
+vol_b = (xM - xm) * (yM - ym) * (zM - zm)
+
+part_vol = vol_b / N_target
+
+# lattice volume
+HCP_PACKING_DENSITY = 0.74
+part_vol_lattice = HCP_PACKING_DENSITY * part_vol
+
+dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0)
+
+pmass = -1
+
+ctx = shamrock.Context()
+ctx.pdata_layout_new()
+
+model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6")
+
+cfg = model.gen_default_config()
+# cfg.set_artif_viscosity_Constant(alpha_u = 1, alpha_AV = 1, beta_AV = 2)
+# cfg.set_artif_viscosity_VaryingMM97(alpha_min = 0.1,alpha_max = 1,sigma_decay = 0.1, alpha_u = 1, beta_AV = 2)
+cfg.set_artif_viscosity_VaryingCD10(
+    alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
+)
+cfg.set_dust_mode_monofluid_tvi(ndust)
+cfg.set_dust_stopping_times(stopping_times)
+cfg.set_boundary_periodic()
+cfg.set_eos_isothermal(1)
+cfg.print_status()
+model.set_solver_config(cfg)
+
+model.init_scheduler(int(1e8), 1)
+
+
+bmin, bmax = model.get_ideal_hcp_box(dr, bmin, bmax)
+xm, ym, zm = bmin
+xM, yM, zM = bmax
+
+model.resize_simulation_box(bmin, bmax)
+
+setup = model.get_setup()
+gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
+setup.apply_setup(gen, insert_step=scheduler_split_val)
+
+for i in range(ndust):
+
+    def func_s(r):
+        return func_s_j(r, i)
+
+    model.set_field_value_lambda_f64("s_j", func_s, i)
+
+model.set_field_value_lambda_f64("uint", uint_g)
+
+vol_b = (xM - xm) * (yM - ym) * (zM - zm)
+totmass = rho * vol_b
+print("Total mass :", totmass)
+
+pmass = model.total_mass_to_part_mass(totmass)
+model.set_particle_mass(pmass)
+print("Current part mass :", pmass)
+
+
+model.set_cfl_cour(0.1)
+model.set_cfl_force(0.1)
+
+model.timestep()
+
+tnext = 0
+for j in range(10):
+    if j > 0:
+        tnext += 0.1
+        model.evolve_until(tnext)
+
+    dic = ctx.collect_data()
+    print(dic["s_j"])
+
+    print(dic["xyz"].shape)
+
+    x = dic["xyz"][:, 0]
+    y = dic["xyz"][:, 1]
+    z = dic["xyz"][:, 2]
+    s_j = dic["s_j"].reshape(-1, ndust)
+    ds_j_dt = dic["ds_j_dt"].reshape(-1, ndust)
+    cs = dic["soundspeed"]
+
+    print(s_j)
+
+    r = np.sqrt(x * x + y * y + z * z)
+
+    hpart = dic["hpart"]
+    rho = pmass * (model.get_hfact() / np.array(hpart)) ** 3
+
+    fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(15, 5))
+    axs[0].scatter(r, rho, label="rho")
+    axs[0].legend()
+    for i in range(ndust):
+        axs[1].scatter(r, s_j[:, i] ** 2 / rho, label=f"eps_j_{i}")
+    axs[1].legend()
+
+    axs[2].scatter(r, cs, label="soundspeed")
+    axs[2].legend()
+
+    for i in range(ndust):
+        axs[3].scatter(r, ds_j_dt[:, i], label=f"ds_j_dt_{i}")
+    axs[3].legend()
+
+    for k in range(4):
+        axs[k].set_xlim(-0.1, 1.1)
+    plt.savefig(f"mono_{j}.png")
+    plt.close()
diff --git a/examples/sph/run_dusty_settle.py b/examples/sph/run_dusty_settle.py
new file mode 100644
index 000000000..e92741c99
--- /dev/null
+++ b/examples/sph/run_dusty_settle.py
@@ -0,0 +1,257 @@
+"""
+Testing Sod tube with SPH
+=========================
+
+CI test for Sod tube with SPH
+"""
+
+import matplotlib.pyplot as plt
+
+import shamrock
+
+shamrock.enable_experimental_features()
+import numpy as np
+
+gamma = 1.4
+rho_i = 1
+central_mass = 1
+R0 = 1
+H_r_0 = 0.05
+
+si = shamrock.UnitSystem()
+sicte = shamrock.Constants(si)
+codeu = shamrock.UnitSystem(
+    unit_time=3600 * 24 * 365,
+    unit_length=sicte.au(),
+    unit_mass=sicte.sol_mass(),
+)
+ucte = shamrock.Constants(codeu)
+G = ucte.G()
+
+
+def kep_profile(r):
+    return (G * central_mass / r) ** 0.5
+
+
+def omega_k(r):
+    return kep_profile(r) / r
+
+
+H = H_r_0 * R0
+cs = H * omega_k(R0)
+box = 8 * H
+
+print(f"cs = {cs}")
+print(f"H = {H}")
+
+
+def scaling_rho(r):
+    x, y, z = r
+
+    loc_h = H / (2**0.5)
+    gaussian = np.exp(-(y**2) / (2 * loc_h * loc_h)) / (loc_h * np.sqrt(2 * np.pi))
+    return gaussian
+
+
+def func_rho_t(r):
+    return rho_i * scaling_rho(r)
+
+
+def func_rho_d_j(r, idust):
+    return (0.1 / ndust) * rho_i * scaling_rho(r)
+
+
+def func_rho_g(r):
+    return rho_i * scaling_rho(r) - sum([func_rho_d_j(r, i) for i in range(ndust)])
+
+
+cs_g = 1
+
+
+def uint_g(r):
+    rho_g = func_rho_g(r)
+    P = rho_g * cs_g * cs_g / gamma
+    return P / ((gamma - 1) * rho_g)
+
+
+ndust = 4
+rc = 0.25
+stokes = np.logspace(-3, 0, ndust)
+stopping_times = stokes / omega_k(R0) 
+print(stopping_times, omega_k(R0))
+from scipy.special import erfinv
+
+bmin = (-box / 4, -box, -box / 4)
+bmax = (box / 4, box, box / 4)
+
+N_target = 1e5
+scheduler_split_val = int(2e7)
+scheduler_merge_val = int(1)
+
+xm, ym, zm = bmin
+xM, yM, zM = bmax
+vol_b = (xM - xm) * (yM - ym) * (zM - zm)
+
+part_vol = vol_b / N_target
+
+# lattice volume
+HCP_PACKING_DENSITY = 0.74
+part_vol_lattice = HCP_PACKING_DENSITY * part_vol
+
+dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0)
+
+pmass = -1
+
+ctx = shamrock.Context()
+ctx.pdata_layout_new()
+
+model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6")
+
+cfg = model.gen_default_config()
+# cfg.set_artif_viscosity_Constant(alpha_u = 1, alpha_AV = 1, beta_AV = 2)
+# cfg.set_artif_viscosity_VaryingMM97(alpha_min = 0.1,alpha_max = 1,sigma_decay = 0.1, alpha_u = 1, beta_AV = 2)
+cfg.set_artif_viscosity_VaryingCD10(
+    alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
+)
+cfg.set_dust_mode_monofluid_tvi(ndust)
+cfg.set_dust_stopping_times(stopping_times)
+cfg.add_ext_force_vertical_disc_potential(central_mass=1, R0=1)
+cfg.add_ext_force_velocity_dissipation(eta=10)
+cfg.set_boundary_periodic()
+cfg.set_units(codeu)
+cfg.set_eos_isothermal(cs)
+cfg.print_status()
+model.set_solver_config(cfg)
+
+model.init_scheduler(int(1e8), 1)
+
+
+bmin, bmax = model.get_ideal_hcp_box(dr, bmin, bmax)
+xm, ym, zm = bmin
+xM, yM, zM = bmax
+
+model.resize_simulation_box(bmin, bmax)
+
+setup = model.get_setup()
+gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
+setup.apply_setup(gen, insert_step=scheduler_split_val)
+
+vol_b = (xM - xm) * (yM - ym) * (zM - zm)
+totmass = rho_i * vol_b
+print("Total mass :", totmass)
+
+pmass = model.total_mass_to_part_mass(totmass)
+model.set_particle_mass(pmass)
+print("Current part mass :", pmass)
+
+# Correct the barycenter
+analysis_barycenter = shamrock.model_sph.analysisBarycenter(model=model)
+barycenter, disc_mass = analysis_barycenter.get_barycenter()
+
+if shamrock.sys.world_rank() == 0:
+    print(f"disc barycenter = {barycenter}")
+
+model.apply_position_offset((-barycenter[0], -barycenter[1], -barycenter[2]))
+
+
+def f_remap(r):
+    x, y, z = r
+
+    rn = max(abs(yM), abs(ym))
+    print(y, H, H * erfinv(y / rn))
+    y = H * erfinv(y / rn)
+
+    y = min(y, yM)
+    y = max(y, ym)
+    return (x, y, z)
+
+
+model.remap_positions(f_remap)
+model.set_field_value_lambda_f64("uint", uint_g)
+
+
+model.set_cfl_cour(0.1)
+model.set_cfl_force(0.1)
+
+model.timestep()
+
+def compute_sj_new(patchdata):
+    hpart = patchdata["hpart"]
+    rho = pmass * (model.get_hfact() / np.array(hpart)) ** 3
+
+    epsilon_target = 0.1 / ndust
+    s = np.sqrt(rho * epsilon_target)
+
+    return s
+
+
+
+
+
+
+# TODO: add function to modify fields e.g. get rho and do stuff according to it
+
+tnext = 0
+for j in range(1000):
+    if j > 0:
+        tnext += 0.1
+        model.evolve_until(tnext)
+
+        if(j == 20):
+            for k in range(ndust):
+                model.overwrite_field_value_f64("s_j",compute_sj_new,k)
+
+    dic = ctx.collect_data()
+    print(dic["s_j"])
+
+    print(dic["xyz"].shape)
+
+    x = dic["xyz"][:, 0]
+    y = dic["xyz"][:, 1]
+    z = dic["xyz"][:, 2]
+    s_j = dic["s_j"].reshape(-1, ndust)
+    ds_j_dt = dic["ds_j_dt"].reshape(-1, ndust)
+    cs = dic["soundspeed"]
+
+    print(s_j)
+
+    r = y
+
+    hpart = dic["hpart"]
+    rho = pmass * (model.get_hfact() / np.array(hpart)) ** 3
+
+    print("compute original rho")
+    estimated_rho = [func_rho_t(dic["xyz"][kk]) for kk in range(len(dic["xyz"]))]
+
+    sz = 2
+
+    fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(12, 5))
+    time = model.get_time()
+    fig.suptitle(f"t = {time:.2f}")
+
+    fig.subplots_adjust(left=0.07, right=0.98, wspace=0.4)
+
+    axs[0].scatter(y, rho, label="gas", s=sz)
+    for i in range(ndust):
+        axs[0].scatter(y, s_j[:, i] ** 2, label=f"dust {i} ts = {stopping_times[i]:.2f}", s=sz)
+    # axs[0].scatter(y,estimated_rho)
+    axs[0].set_ylabel(r"$\rho$")
+    axs[0].set_xlabel(r"$y$")
+
+    axs[0].set_yscale('log')
+    axs[0].legend()
+    for i in range(ndust):
+        axs[1].scatter(y, s_j[:, i] ** 2 / rho, label=f"dust {i} ts = {stopping_times[i]:.2f}", s=sz)
+    axs[1].set_ylabel(r"$\epsilon_j$")
+    axs[1].set_xlabel(r"$y$")
+    axs[1].legend()
+
+    for i in range(ndust):
+        axs[2].scatter(y, ds_j_dt[:, i], label=f"dust {i} ts = {stopping_times[i]:.2f}", s=sz)
+    axs[2].set_ylabel(r"$\frac{d s_j}{dt}$")
+    axs[2].set_xlabel(r"$y$")
+    axs[2].legend()
+
+    plt.savefig(f"mono_{j}.png")
+    model.do_vtk_dump(f"dump_stratif_{j}.vtk", True)
+    plt.close()
diff --git a/src/shammodels/common/include/shammodels/common/modules/ForwardEuler.hpp b/src/shammodels/common/include/shammodels/common/modules/ForwardEuler.hpp
index b0f90caa1..6450cdeef 100644
--- a/src/shammodels/common/include/shammodels/common/modules/ForwardEuler.hpp
+++ b/src/shammodels/common/include/shammodels/common/modules/ForwardEuler.hpp
@@ -31,8 +31,10 @@ namespace shammodels::common::modules {
 
         using Tscal = shambase::VecComponent<T>;
 
+        u32 nvar;
+
         public:
-        ForwardEuler() = default;
+        ForwardEuler(u32 nvar = 1) : nvar(nvar) {}
 
         struct Edges {
             const shamrock::solvergraph::IDataEdge<Tscal> &dt;
@@ -68,14 +70,32 @@ namespace shammodels::common::modules {
 
             Tscal dt = edges.dt.data;
 
-            sham::distributed_data_kernel_call(
-                shamsys::instance::get_compute_scheduler_ptr(),
-                sham::DDMultiRef{edges.time_derivative.get_spans()},
-                sham::DDMultiRef{edges.field.get_spans()},
-                edges.sizes.indexes,
-                [dt](u32 gid, const T *time_derivative, T *field) {
-                    field[gid] = field[gid] + dt * time_derivative[gid];
+            if (nvar == 1) {
+
+                sham::distributed_data_kernel_call(
+                    shamsys::instance::get_compute_scheduler_ptr(),
+                    sham::DDMultiRef{edges.time_derivative.get_spans()},
+                    sham::DDMultiRef{edges.field.get_spans()},
+                    edges.sizes.indexes,
+                    [dt](u32 gid, const T *time_derivative, T *field) {
+                        field[gid] = field[gid] + dt * time_derivative[gid];
+                    });
+
+            } else {
+
+                auto var_count = edges.sizes.indexes.template map<u32>([&](u64 id, u32 count) {
+                    return count * nvar;
                 });
+
+                sham::distributed_data_kernel_call(
+                    shamsys::instance::get_compute_scheduler_ptr(),
+                    sham::DDMultiRef{edges.time_derivative.get_spans()},
+                    sham::DDMultiRef{edges.field.get_spans()},
+                    var_count,
+                    [dt](u32 gid, const T *time_derivative, T *field) {
+                        field[gid] = field[gid] + dt * time_derivative[gid];
+                    });
+            }
         }
 
         inline virtual std::string _impl_get_label() const { return "ForwardEuler"; };
diff --git a/src/shammodels/sph/CMakeLists.txt b/src/shammodels/sph/CMakeLists.txt
index 295ff990c..795e09301 100644
--- a/src/shammodels/sph/CMakeLists.txt
+++ b/src/shammodels/sph/CMakeLists.txt
@@ -51,6 +51,7 @@ set(Sources
     src/modules/LoopSmoothingLengthIter.cpp
     src/modules/ComputeNeighStats.cpp
     src/modules/NodeUpdateDerivsVaryingAlphaAV.cpp
+    src/modules/NodeUpdateDerivsMonofluidTVI.cpp
     src/modules/setup/GeneratorMCDisc.cpp
     src/modules/setup/ModifierApplyDiscWarp.cpp
     src/modules/setup/ModifierApplyCustomWarp.cpp
diff --git a/src/shammodels/sph/include/shammodels/sph/Model.hpp b/src/shammodels/sph/include/shammodels/sph/Model.hpp
index 41dbc517a..6e43a7bb5 100644
--- a/src/shammodels/sph/include/shammodels/sph/Model.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/Model.hpp
@@ -31,6 +31,7 @@
 #include "shammodels/sph/math/density.hpp"
 #include "shammodels/sph/modules/ComputeLoadBalanceValue.hpp"
 #include "shammodels/sph/modules/SPHSetup.hpp"
+#include "shampylib/PatchDataToPy.hpp"
 #include "shamrock/io/ShamrockDump.hpp"
 #include "shamrock/patch/PatchDataLayer.hpp"
 #include "shamrock/scheduler/ReattributeDataUtility.hpp"
@@ -186,38 +187,70 @@ namespace shammodels::sph {
 
         template<class T>
         inline void set_field_value_lambda(
-            std::string field_name, const std::function<T(Tvec)> pos_to_val) {
+            std::string field_name, const std::function<T(Tvec)> pos_to_val, const i32 offset) {
 
             StackEntry stack_loc{};
+
             PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
+
+            u32 ixyz   = sched.pdl_old().get_field_idx<Tvec>("xyz");
+            u32 ifield = sched.pdl_old().get_field_idx<T>(field_name);
+
             sched.patch_data.for_each_patchdata(
                 [&](u64 patch_id, shamrock::patch::PatchDataLayer &pdat) {
-                    PatchDataField<Tvec> &xyz
-                        = pdat.template get_field<Tvec>(sched.pdl_old().get_field_idx<Tvec>("xyz"));
+                    PatchDataField<Tvec> &xyz = pdat.template get_field<Tvec>(ixyz);
+                    PatchDataField<T> &f      = pdat.template get_field<T>(ifield);
 
-                    PatchDataField<T> &f
-                        = pdat.template get_field<T>(sched.pdl_old().get_field_idx<T>(field_name));
+                    auto f_nvar = f.get_nvar();
 
-                    if (f.get_nvar() != 1) {
-                        shambase::throw_unimplemented();
+                    auto acc     = f.get_buf().copy_to_stdvec();
+                    auto acc_xyz = xyz.get_buf().copy_to_stdvec();
+
+                    for (u32 i = 0; i < pdat.get_obj_cnt(); i++) {
+                        acc[i * f_nvar + offset] = pos_to_val(acc_xyz[i]);
                     }
 
-                    {
-                        auto &buf = f.get_buf();
-                        auto acc  = buf.copy_to_stdvec();
+                    logger::raw_ln(acc);
 
-                        auto &buf_xyz = xyz.get_buf();
-                        auto acc_xyz  = buf_xyz.copy_to_stdvec();
+                    f.get_buf().copy_from_stdvec(acc);
+                });
+        }
 
-                        for (u32 i = 0; i < f.get_obj_cnt(); i++) {
-                            Tvec r = acc_xyz[i];
+        template<class T>
+        inline void overwrite_field_value(
+            std::string field_name,
+            const std::function<std::vector<T>(py::dict)> field_compute,
+            const i32 offset) {
 
-                            acc[i] = pos_to_val(r);
-                        }
+            StackEntry stack_loc{};
+
+            PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
+
+            u32 ifield = sched.pdl_old().get_field_idx<T>(field_name);
 
-                        buf.copy_from_stdvec(acc);
-                        buf_xyz.copy_from_stdvec(acc_xyz);
+            sched.patch_data.for_each_patchdata(
+                [&](u64 patch_id, shamrock::patch::PatchDataLayer &pdat) {
+                    PatchDataField<T> &f = pdat.template get_field<T>(ifield);
+
+                    auto result = field_compute(shamrock::pdat_to_dic(pdat));
+
+                    if (result.size() != f.get_obj_cnt()) {
+                        throw shambase::make_except_with_loc<std::runtime_error>(shambase::format(
+                            "result.size() != f.get_obj_cnt() ({} != {})",
+                            result.size(),
+                            f.get_obj_cnt()));
                     }
+
+                    auto f_nvar = f.get_nvar();
+                    auto acc    = f.get_buf().copy_to_stdvec();
+
+                    for (u32 i = 0; i < pdat.get_obj_cnt(); i++) {
+                        acc[i * f_nvar + offset] = result[i];
+                    }
+
+                    logger::raw_ln(acc);
+
+                    f.get_buf().copy_from_stdvec(acc);
                 });
         }
 
@@ -948,9 +981,9 @@ namespace shammodels::sph {
                         * (1.
                            + sycl::sin(shambase::constants::pi<Tscal> / (2. * Hwarp) * (R - Rwarp)))
                         * sycl::sin(incl_rad));
-                    psi          = shambase::constants::pi<Tscal>
-                                   * Rwarp / (4. * Hwarp) * sycl::sin(incl_rad)
-                                   / sycl::sqrt(1. - (0.5 * sycl::pow(sycl::sin(incl_rad), 2)));
+                    psi = shambase::constants::pi<Tscal>
+                          * Rwarp / (4. * Hwarp) * sycl::sin(incl_rad)
+                          / sycl::sqrt(1. - (0.5 * sycl::pow(sycl::sin(incl_rad), 2)));
                     Tscal psimax = sycl::max(psimax, psi);
                     Tscal x      = pos[i].x();
                     Tscal y      = pos[i].y();
diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
index c6454f00a..6cf93b69d 100644
--- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp
@@ -122,6 +122,8 @@ namespace shammodels::sph {
             u32 ndust;
         };
 
+        std::vector<Tscal> stopping_times;
+
         /// Variant type to store the EOS configuration
         using Variant = std::variant<None, MonofluidTVI, MonofluidComplete>;
 
@@ -131,9 +133,17 @@ namespace shammodels::sph {
         inline void set_monofluid_tvi(u32 nvar) { current_mode = MonofluidTVI{nvar}; }
         inline void set_monofluid_complete(u32 nvar) { current_mode = MonofluidComplete{nvar}; }
 
+        inline bool is_monofluid_tvi() { return bool(std::get_if<MonofluidTVI>(&current_mode)); }
+        inline bool is_monofluid_complete() {
+            return bool(std::get_if<MonofluidComplete>(&current_mode));
+        }
+
+        inline bool has_s_j_field() {
+            return is_monofluid_tvi(); // S_j = sqrt(\rho \epsilon_j)
+        }
+
         inline bool has_epsilon_field() {
-            return bool(std::get_if<MonofluidTVI>(&current_mode))
-                   || bool(std::get_if<MonofluidComplete>(&current_mode));
+            return bool(std::get_if<MonofluidComplete>(&current_mode));
         }
 
         inline bool has_deltav_field() {
@@ -159,10 +169,19 @@ namespace shammodels::sph {
             bool is_not_none = bool(std::get_if<MonofluidTVI>(&current_mode))
                                || bool(std::get_if<MonofluidComplete>(&current_mode));
             if (is_not_none) {
-                ON_RANK_0(
-                    logger::warn_ln(
-                        "SPH::config",
-                        "Dust config != None is work in progress, use it at your own risk"));
+                if (!shamrock::are_experimental_features_allowed()) {
+                    shambase::throw_with_loc<std::runtime_error>(
+                        "Dust config != None is experimental");
+                } else {
+                    ON_RANK_0(
+                        logger::warn_ln(
+                            "SPH::config",
+                            "Dust config != None is work in progress, use it at your own risk"));
+                }
+
+                if (get_dust_nvar() != stopping_times.size()) {
+                    throw "";
+                }
             }
         }
     };
diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp
new file mode 100644
index 000000000..085dcd665
--- /dev/null
+++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp
@@ -0,0 +1,136 @@
+// -------------------------------------------------------//
+//
+// 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 ComputeDustTtilde.hpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief
+ *
+ */
+
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shambackends/vec.hpp"
+#include "shamcomm/logs.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_COMPUTE_DUST_TTILDE_EDGES(X_RO, X_RW)                                                 \
+    /* scalars */                                                                                  \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, gpart_mass)                                     \
+                                                                                                   \
+    /* counts */                                                                                   \
+    X_RO(shamrock::solvergraph::Indexes<u32>, part_counts)                                         \
+                                                                                                   \
+    /* fields */                                                                                   \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, hpart)                                          \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, s_j)                                            \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, t_j)                                            \
+                                                                                                   \
+    /* outputs */                                                                                  \
+    X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, Ttilde_sj)
+
+namespace shammodels::sph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class ComputeDustTtilde : public shamrock::solvergraph::INode {
+
+        using Tscal  = shambase::VecComponent<Tvec>;
+        using Kernel = SPHKernel<Tscal>;
+
+        static constexpr Tscal kernel_radius = SPHKernel<Tscal>::Rkern;
+
+        u32 ndust;
+
+        public:
+        ComputeDustTtilde(u32 ndust) : ndust(ndust) {}
+
+        EXPAND_NODE_EDGES(NODE_COMPUTE_DUST_TTILDE_EDGES)
+
+        inline void _impl_evaluate_internal() {
+
+            __shamrock_stack_entry();
+
+            auto edges = get_edges();
+
+            auto &part_counts = edges.part_counts.indexes;
+
+            // check that all input edges have the particles with ghosts zones
+            edges.hpart.check_sizes(part_counts);
+            edges.s_j.check_sizes(part_counts);
+            edges.t_j.check_sizes(part_counts);
+
+            // ensure that the output edges are of size part_counts (output without ghosts zones)
+            edges.Ttilde_sj.ensure_sizes(part_counts);
+
+            const Tscal pmass = edges.gpart_mass.value;
+
+            auto total_specie_count = part_counts.template map<u32>([&](u64 id, u32 count) {
+                return count * ndust;
+            });
+
+            // call the kernel for each patches with part_counts.get(id_patch) threads of patch
+            // id_patch
+            sham::distributed_data_kernel_call(
+                shamsys::instance::get_compute_scheduler_ptr(),
+                sham::DDMultiRef{
+                    edges.hpart.get_spans(), edges.s_j.get_spans(), edges.t_j.get_spans()},
+                sham::DDMultiRef{edges.Ttilde_sj.get_spans()},
+                total_specie_count,
+                [pmass, ndust = ndust](
+                    u32 thread_id,
+                    const Tscal *__restrict hpart,
+                    const Tscal *__restrict s_j,
+                    const Tscal *__restrict t_j,
+                    Tscal *__restrict Ttilde_sj) {
+                    u32 id_a  = thread_id / ndust;
+                    u32 jdust = thread_id % ndust;
+
+                    Tscal h_a  = hpart[id_a];
+                    Tscal sj_a = s_j[thread_id];
+                    Tscal tj_a = t_j[thread_id];
+
+                    using namespace shamrock::sph;
+                    Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
+
+                    auto epsilon = [&](Tscal sj) {
+                        return sj * sj / rho_a;
+                    };
+                    /*
+                     * Hutchison 2018 eq 15
+                     * \tilde{T}_{sj} = \epsilon_j t_j - \sum_{k} \epsilon_k^2 t_k
+                     */
+
+                    Tscal eps_j_a     = epsilon(sj_a);
+                    Tscal Ttilde_sj_a = eps_j_a * tj_a;
+
+                    for (u32 k = 0; k < ndust; k++) {
+                        Tscal sk_a = s_j[id_a * ndust + k];
+                        Tscal tk_a = t_j[id_a * ndust + k];
+
+                        Tscal eps_k_a = epsilon(sk_a);
+                        Ttilde_sj_a -= eps_k_a * eps_k_a * tk_a;
+                    }
+
+                    // logger::raw_ln("Ttilde_sj_a", jdust, Ttilde_sj_a, eps_j_a, tj_a);
+
+                    Ttilde_sj[thread_id] = Ttilde_sj_a;
+                });
+        }
+
+        inline virtual std::string _impl_get_label() const { return "ComputeDustTtilde"; };
+
+        inline virtual std::string _impl_get_tex() const { return "TODO"; };
+    };
+} // namespace shammodels::sph::modules
diff --git a/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp b/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp
new file mode 100644
index 000000000..0fbc75410
--- /dev/null
+++ b/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp
@@ -0,0 +1,73 @@
+// -------------------------------------------------------//
+//
+// 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 NodeUpdateDerivsMonofluidTVI.hpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief
+ *
+ */
+
+#include "shambackends/vec.hpp"
+#include "shammodels/sph/solvergraph/NeighCache.hpp"
+#include "shamrock/solvergraph/IFieldSpan.hpp"
+#include "shamrock/solvergraph/INode.hpp"
+#include "shamrock/solvergraph/Indexes.hpp"
+#include "shamrock/solvergraph/ScalarEdge.hpp"
+
+#define NODE_UPDATE_DERIVS_MONOFLUID_TVI_EDGES(X_RO, X_RW)                                         \
+    /* scalars */                                                                                  \
+    X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, gpart_mass)                                     \
+                                                                                                   \
+    /* counts */                                                                                   \
+    X_RO(shamrock::solvergraph::Indexes<u32>, part_counts)                                         \
+    X_RO(shamrock::solvergraph::Indexes<u32>, part_counts_with_ghost)                              \
+                                                                                                   \
+    /* fields */                                                                                   \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tvec>, xyz)                                             \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, hpart)                                          \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tvec>, vxyz)                                            \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, omega)                                          \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, pressure)                                       \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, s_j)                                            \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, Ttilde_sj)                                      \
+                                                                                                   \
+    /* neigh */                                                                                    \
+    X_RO(shammodels::sph::solvergraph::NeighCache, neigh_cache)                                    \
+                                                                                                   \
+    /* outputs */                                                                                  \
+    X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, ds_j_dt)
+
+namespace shammodels::sph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class NodeUpdateDerivsMonofluidTVI : public shamrock::solvergraph::INode {
+
+        using Tscal = shambase::VecComponent<Tvec>;
+
+        static constexpr Tscal kernel_radius = SPHKernel<Tscal>::Rkern;
+
+        u32 ndust;
+
+        public:
+        NodeUpdateDerivsMonofluidTVI(u32 ndust) : ndust(ndust) {}
+
+        EXPAND_NODE_EDGES(NODE_UPDATE_DERIVS_MONOFLUID_TVI_EDGES)
+
+        void _impl_evaluate_internal();
+
+        inline virtual std::string _impl_get_label() const {
+            return "NodeUpdateDerivsMonofluidTVI";
+        };
+
+        inline virtual std::string _impl_get_tex() const { return "TODO"; };
+    };
+} // namespace shammodels::sph::modules
diff --git a/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp b/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp
new file mode 100644
index 000000000..cc5e70c86
--- /dev/null
+++ b/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp
@@ -0,0 +1,94 @@
+// -------------------------------------------------------//
+//
+// 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 SetDustStoppingTimeConstant.hpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief
+ *
+ */
+
+#include "shambackends/DeviceBuffer.hpp"
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shambackends/vec.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"
+#include <vector>
+
+#define NODE_SET_DUST_STOPPING_TIME_CTE_EDGES(X_RO, X_RW)                                          \
+    /* scalars */                                                                                  \
+    X_RO(shamrock::solvergraph::ScalarEdge<std::vector<Tscal>>, t_j_0)                             \
+                                                                                                   \
+    /* counts */                                                                                   \
+    X_RO(shamrock::solvergraph::Indexes<u32>, part_counts)                                         \
+                                                                                                   \
+    /* fields */                                                                                   \
+    X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, t_j)
+
+namespace shammodels::sph::modules {
+
+    template<class Tvec, template<class> class SPHKernel>
+    class SetDustStoppingTimeConstant : public shamrock::solvergraph::INode {
+
+        using Tscal  = shambase::VecComponent<Tvec>;
+        using Kernel = SPHKernel<Tscal>;
+
+        static constexpr Tscal kernel_radius = SPHKernel<Tscal>::Rkern;
+
+        u32 ndust;
+
+        public:
+        SetDustStoppingTimeConstant(u32 ndust) : ndust(ndust) {}
+
+        EXPAND_NODE_EDGES(NODE_SET_DUST_STOPPING_TIME_CTE_EDGES)
+
+        inline void _impl_evaluate_internal() {
+
+            __shamrock_stack_entry();
+
+            auto edges = get_edges();
+
+            auto &part_counts                   = edges.part_counts.indexes;
+            const std::vector<Tscal> &inputs_tj = edges.t_j_0.value;
+
+            // ensure that the output edges are of size part_counts
+            edges.t_j.ensure_sizes(part_counts);
+
+            sham::DeviceBuffer<Tscal> t_j_0(ndust, shamsys::instance::get_compute_scheduler_ptr());
+            t_j_0.copy_from_stdvec(inputs_tj);
+
+            auto &q = shamsys::instance::get_compute_scheduler().get_queue();
+
+            part_counts.for_each([&](u64 id, u32 count) {
+                // call the kernel for each patches with part_counts.get(id_patch) threads of patch
+                // id_patch
+                sham::kernel_call(
+                    q,
+                    sham::MultiRef{t_j_0},
+                    sham::MultiRef{edges.t_j.get_spans().get(id)},
+                    part_counts.get(id) * ndust,
+                    [ndust
+                     = ndust](u32 thread_id, const Tscal *__restrict t_j_0, Tscal *__restrict t_j) {
+                        u32 jdust      = thread_id % ndust;
+                        t_j[thread_id] = t_j_0[jdust];
+                    });
+            });
+        }
+
+        inline virtual std::string _impl_get_label() const { return "ComputeDustTtilde"; };
+
+        inline virtual std::string _impl_get_tex() const { return "TODO"; };
+    };
+} // namespace shammodels::sph::modules
diff --git a/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp b/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp
index 3075248eb..78fc5a54b 100644
--- a/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp
@@ -60,6 +60,10 @@ namespace shammodels::sph::modules {
         void update_derivs_cd10(VaryingCD10 cfg);
         void update_derivs_disc_visco(ConstantDisc cfg);
 
+        using DustConfig = typename Config::DustConfig;
+
+        void update_derivs_dust_monofluid_tvi_Sj(DustConfig cfg);
+
         using Cfg_MHD = typename Config::MHDConfig;
 
         using NoneMHD     = typename Cfg_MHD::None;
diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp
index 0ef3a5ad0..55b7386cd 100644
--- a/src/shammodels/sph/src/Solver.cpp
+++ b/src/shammodels/sph/src/Solver.cpp
@@ -112,6 +112,7 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
     bool has_psi_field                         = solver_config.has_field_psi_on_ch();
     bool has_epsilon_field                     = solver_config.dust_config.has_epsilon_field();
     bool has_deltav_field                      = solver_config.dust_config.has_deltav_field();
+    bool has_s_j_field                         = solver_config.dust_config.has_s_j_field();
 
     using namespace shamrock::solvergraph;
 
@@ -148,6 +149,10 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
         solver_graph.register_edge("deltav", FieldRefs<Tvec>("deltav", "\\Delta v"));
         solver_graph.register_edge("dtdeltav", FieldRefs<Tvec>("dtdeltav", "d\\Delta v"));
     }
+    if (has_s_j_field) {
+        solver_graph.register_edge("s_j", FieldRefs<Tscal>("s_j", "S_j"));
+        solver_graph.register_edge("ds_j_dt", FieldRefs<Tscal>("ds_j_dt", "dS_j/dt"));
+    }
 
     {
         auto set_gpart_mass = solver_graph.register_node(
@@ -330,6 +335,25 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
             attach_field_sequence.push_back(attach_dtdeltav);
         }
 
+        if (has_s_j_field) {
+            auto attach_s_j
+                = solver_graph.register_node("attach_s_j", GetFieldRefFromLayer<Tscal>(pdl, "s_j"));
+            shambase::get_check_ref(attach_s_j)
+                .set_edges(
+                    solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
+                    solver_graph.get_edge_ptr<FieldRefs<Tscal>>("s_j"));
+            attach_field_sequence.push_back(attach_s_j);
+        }
+
+        if (has_s_j_field) {
+            auto attach_ds_j_dt = solver_graph.register_node(
+                "attach_ds_j_dt", GetFieldRefFromLayer<Tscal>(pdl, "ds_j_dt"));
+            shambase::get_check_ref(attach_ds_j_dt)
+                .set_edges(
+                    solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
+                    solver_graph.get_edge_ptr<FieldRefs<Tscal>>("ds_j_dt"));
+            attach_field_sequence.push_back(attach_ds_j_dt);
+        }
         solver_graph.register_node(
             "attach fields to scheduler",
             OperationSequence("attach fields", std::move(attach_field_sequence)));
@@ -416,6 +440,19 @@ void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
                 half_step_sequence.push_back(half_step_deltav);
             }
 
+            if (has_s_j_field) {
+                u32 ndust          = solver_config.dust_config.get_dust_nvar();
+                auto half_step_s_j = solver_graph.register_node(
+                    prefix + "_s_j", shammodels::common::modules::ForwardEuler<Tscal>(ndust));
+                shambase::get_check_ref(half_step_s_j)
+                    .set_edges(
+                        solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
+                        solver_graph.get_edge_ptr<FieldRefs<Tscal>>("ds_j_dt"),
+                        solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
+                        solver_graph.get_edge_ptr<FieldRefs<Tscal>>("s_j"));
+                half_step_sequence.push_back(half_step_s_j);
+            }
+
             return OperationSequence("half step", std::move(half_step_sequence));
         };
 
@@ -1269,6 +1306,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
     bool has_curlB_field   = solver_config.has_field_curlB();
     bool has_epsilon_field = solver_config.dust_config.has_epsilon_field();
     bool has_deltav_field  = solver_config.dust_config.has_deltav_field();
+    bool has_s_j_field     = solver_config.dust_config.has_s_j_field();
 
     PatchDataLayerLayout &pdl = scheduler().pdl_old();
     const u32 ixyz            = pdl.get_field_idx<Tvec>("xyz");
@@ -1299,6 +1337,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
 
     const u32 iepsilon = (has_epsilon_field) ? pdl.get_field_idx<Tscal>("epsilon") : 0;
     const u32 ideltav  = (has_deltav_field) ? pdl.get_field_idx<Tvec>("deltav") : 0;
+    const u32 is_j     = (has_s_j_field) ? pdl.get_field_idx<Tscal>("s_j") : 0;
 
     auto &ghost_layout_ptr                              = storage.ghost_layout;
     shamrock::patch::PatchDataLayerLayout &ghost_layout = shambase::get_check_ref(ghost_layout_ptr);
@@ -1320,6 +1359,7 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
     const u32 iepsilon_interf
         = (has_epsilon_field) ? ghost_layout.get_field_idx<Tscal>("epsilon") : 0;
     const u32 ideltav_interf = (has_deltav_field) ? ghost_layout.get_field_idx<Tvec>("deltav") : 0;
+    const u32 is_j_interf    = (has_s_j_field) ? ghost_layout.get_field_idx<Tscal>("s_j") : 0;
 
     using InterfaceBuildInfos = typename sph::BasicSPHGhostHandler<Tvec>::InterfaceBuildInfos;
 
@@ -1392,6 +1432,11 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
                 sender_patch.get_field<Tvec>(ideltav).append_subset_to(
                     buf_idx, cnt, pdat.get_field<Tvec>(ideltav_interf));
             }
+
+            if (has_s_j_field) {
+                sender_patch.get_field<Tscal>(is_j).append_subset_to(
+                    buf_idx, cnt, pdat.get_field<Tscal>(is_j_interf));
+            }
         });
 
     ghost_handle.template modify_interface_native<PatchDataLayer>(
@@ -1468,6 +1513,10 @@ void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
                     pdat_new.get_field<Tvec>(ideltav_interf).insert(pdat.get_field<Tvec>(ideltav));
                 }
 
+                if (has_s_j_field) {
+                    pdat_new.get_field<Tscal>(is_j_interf).insert(pdat.get_field<Tscal>(is_j));
+                }
+
                 pdat_new.check_field_obj_cnt_match();
 
                 return pdat_new;
diff --git a/src/shammodels/sph/src/SolverConfig.cpp b/src/shammodels/sph/src/SolverConfig.cpp
index 4259af596..b5d3c8d95 100644
--- a/src/shammodels/sph/src/SolverConfig.cpp
+++ b/src/shammodels/sph/src/SolverConfig.cpp
@@ -111,6 +111,12 @@ namespace shammodels::sph {
         if (should_save_dt_to_fields()) {
             pdl.add_field<Tscal>("dt_part", 1);
         }
+
+        if (dust_config.has_s_j_field()) {
+            u32 ndust = dust_config.get_dust_nvar();
+            pdl.add_field<Tscal>("s_j", ndust);
+            pdl.add_field<Tscal>("ds_j_dt", ndust);
+        }
     }
 
     template<class Tvec, template<class> class SPHKernel>
@@ -151,6 +157,11 @@ namespace shammodels::sph {
             u32 ndust = dust_config.get_dust_nvar();
             ghost_layout.add_field<Tvec>("deltav", ndust);
         }
+
+        if (dust_config.has_s_j_field()) {
+            u32 ndust = dust_config.get_dust_nvar();
+            ghost_layout.add_field<Tscal>("s_j", ndust);
+        }
     }
 
 }; // namespace shammodels::sph
diff --git a/src/shammodels/sph/src/modules/ComputeEos.cpp b/src/shammodels/sph/src/modules/ComputeEos.cpp
index ded3f9b76..81de5b527 100644
--- a/src/shammodels/sph/src/modules/ComputeEos.cpp
+++ b/src/shammodels/sph/src/modules/ComputeEos.cpp
@@ -19,6 +19,7 @@
 #include "shambase/DistributedData.hpp"
 #include "shambase/exception.hpp"
 #include "shambackends/kernel_call.hpp"
+#include "shamcomm/logs.hpp"
 #include "shammath/sphkernels.hpp"
 #include "shammodels/sph/math/density.hpp"
 #include "shammodels/sph/modules/ComputeEos.hpp"
@@ -89,6 +90,48 @@ struct RhoGetterMonofluid {
     void complete_event_state(sycl::event e) { buf_h.complete_event_state(e); }
 };
 
+template<class Tscal>
+struct RhoGetterSJ {
+    sham::DeviceBuffer<Tscal> &buf_h;
+    sham::DeviceBuffer<Tscal> &buf_s_j;
+    u32 nvar_dust;
+    Tscal pmass;
+    Tscal hfact;
+
+    struct accessed {
+        const Tscal *h;
+        const Tscal *buf_s_j;
+        u32 nvar_dust;
+        Tscal pmass;
+        Tscal hfact;
+
+        Tscal operator()(u32 i) const {
+            using namespace shamrock::sph;
+            Tscal rho = rho_h(pmass, h[i], hfact);
+
+            Tscal epsilon_sum = 0;
+            for (u32 j = 0; j < nvar_dust; j++) {
+                Tscal s_j = buf_s_j[i * nvar_dust + j];
+                epsilon_sum += s_j * s_j / rho;
+            }
+
+            return rho * (1 - epsilon_sum);
+        }
+    };
+
+    accessed get_read_access(sham::EventList &depends_list) {
+        auto h   = buf_h.get_read_access(depends_list);
+        auto s_j = buf_s_j.get_read_access(depends_list);
+
+        return accessed{h, s_j, nvar_dust, pmass, hfact};
+    }
+
+    void complete_event_state(sycl::event e) {
+        buf_h.complete_event_state(e);
+        buf_s_j.complete_event_state(e);
+    }
+};
+
 template<class Tvec, template<class> class SPHKernel>
 template<class RhoGetGen>
 void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos_internal(
@@ -561,6 +604,9 @@ void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos() {
     shambase::get_check_ref(storage.pressure).ensure_sizes(counts_with_ghosts);
     shambase::get_check_ref(storage.soundspeed).ensure_sizes(counts_with_ghosts);
 
+    logger::raw_ln(
+        solver_config.dust_config.has_epsilon_field(), solver_config.dust_config.has_s_j_field());
+
     if (solver_config.dust_config.has_epsilon_field()) {
 
         u32 iepsilon_interf = ghost_layout.get_field_idx<Tscal>("epsilon");
@@ -574,6 +620,21 @@ void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos() {
                 gpart_mass,
                 Kernel::hfactd};
         });
+    } else if (solver_config.dust_config.has_s_j_field()) {
+
+        u32 is_j_interf = ghost_layout.get_field_idx<Tscal>("s_j");
+        u32 nvar_dust   = solver_config.dust_config.get_dust_nvar();
+
+        compute_eos_internal([&](PatchDataLayer &mpdat) {
+            return RhoGetterSJ<Tscal>{
+                mpdat.get_field_buf_ref<Tscal>(ihpart_interf),
+                mpdat.get_field_buf_ref<Tscal>(is_j_interf),
+                nvar_dust,
+                gpart_mass,
+                Kernel::hfactd};
+        });
+
+        logger::info_ln("ComputeEos", "SJ field has been computed");
     } else {
         compute_eos_internal([&](PatchDataLayer &mpdat) {
             return RhoGetterBase<Tscal>{
diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp
new file mode 100644
index 000000000..1ba183206
--- /dev/null
+++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp
@@ -0,0 +1,166 @@
+// -------------------------------------------------------//
+//
+// 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 NodeUpdateDerivsMonofluidTVI.cpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief
+ *
+ */
+
+#include "shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp"
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shamcomm/logs.hpp"
+#include "shammath/sphkernels.hpp"
+#include "shammodels/sph/math/density.hpp"
+#include "shamrock/patch/PatchDataField.hpp"
+
+template<class Tvec, template<class> class SPHKernel>
+struct KernelUpdateDerivsMonofluidTVI {
+    using Tscal                   = shambase::VecComponent<Tvec>;
+    using Kernel                  = SPHKernel<Tscal>;
+    static constexpr Tscal hfactd = Kernel::hfactd;
+    static constexpr Tscal Rkern  = Kernel::Rkern;
+    static constexpr Tscal Rker2  = Rkern * Rkern;
+
+    Tscal pmass;
+    u32 ndust;
+
+    inline void operator()(
+        u32 thread_id,
+        // input
+        const Tvec *__restrict xyz,
+        const Tscal *__restrict hpart,
+        const Tvec *__restrict vxyz,
+        const Tscal *__restrict omega,
+        const Tscal *__restrict pressure,
+        const Tscal *__restrict s_j,
+        const Tscal *__restrict Ttilde_sj,
+        shamrock::tree::ObjectCache::ptrs_read ploop_ptrs,
+        // output
+        Tscal *__restrict ds_j_dt) const {
+
+        u32 id_a  = thread_id / ndust;
+        u32 jdust = thread_id % ndust;
+
+        Tscal h_a         = hpart[id_a];
+        Tvec xyz_a        = xyz[id_a];
+        Tvec vxyz_a       = vxyz[id_a];
+        Tscal P_a         = pressure[id_a];
+        Tscal omega_a     = omega[id_a];
+        Tscal s_j_a       = s_j[thread_id];
+        Tscal Ttilde_sj_a = Ttilde_sj[thread_id];
+
+        using namespace shamrock::sph;
+        Tscal rho_a             = rho_h(pmass, h_a, Kernel::hfactd);
+        Tscal rho_a_sq          = rho_a * rho_a;
+        Tscal rho_a_inv         = 1. / rho_a;
+        Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
+
+        Tscal term1 = 0;
+        Tscal term2 = 0;
+
+        shamrock::tree::ObjectCacheIterator particle_looper(ploop_ptrs);
+        particle_looper.for_each_object(id_a, [&](u32 id_b) {
+            Tvec dr    = xyz_a - xyz[id_b];
+            Tscal rab2 = sycl::dot(dr, dr);
+            Tscal h_b  = hpart[id_b];
+
+            if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
+                return;
+            }
+
+            Tvec vxyz_b       = vxyz[id_b];
+            Tscal P_b         = pressure[id_b];
+            Tscal omega_b     = omega[id_b];
+            Tscal s_j_b       = s_j[id_b * ndust + jdust];
+            Tscal Ttilde_sj_b = Ttilde_sj[id_b * ndust + jdust];
+
+            Tscal rab = sycl::sqrt(rab2);
+
+            Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd);
+
+            Tscal Fab_a = Kernel::dW_3d(rab, h_a);
+            Tscal Fab_b = Kernel::dW_3d(rab, h_b);
+
+            Tvec v_ab = vxyz_a - vxyz_b;
+
+            Tvec r_ab_unit = dr * sham::inv_sat_positive(rab);
+
+            Tscal F_ab_bar    = (Fab_a + Fab_b) / 2;
+            Tscal delta_P     = P_a - P_b;
+            Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b);
+
+            // logger::raw_ln("Ts_weighted", Ts_weighted);
+
+            term1 += (pmass * s_j_b / rho_b) * Ts_weighted * delta_P * F_ab_bar;
+            term2 += pmass * sham::dot(v_ab, r_ab_unit * Fab_a);
+        });
+
+        // eq 51, Hutchison 2018
+        ds_j_dt[thread_id] = Tscal{-0.5} * term1 + (s_j_a / (2 * rho_a * omega_a)) * term2;
+    }
+};
+
+template<class Tvec, template<class> class SPHKernel>
+void shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI<Tvec, SPHKernel>::
+    _impl_evaluate_internal() {
+
+    __shamrock_stack_entry();
+
+    auto edges = get_edges();
+
+    auto &part_counts_with_ghost = edges.part_counts_with_ghost.indexes;
+    auto &part_counts            = edges.part_counts.indexes;
+
+    // check that all input edges have the particles with ghosts zones
+    edges.xyz.check_sizes(part_counts_with_ghost);
+    edges.hpart.check_sizes(part_counts_with_ghost);
+    edges.vxyz.check_sizes(part_counts_with_ghost);
+    edges.omega.check_sizes(part_counts_with_ghost);
+    edges.pressure.check_sizes(part_counts_with_ghost);
+    edges.s_j.check_sizes(part_counts_with_ghost);
+    edges.Ttilde_sj.check_sizes(part_counts_with_ghost);
+
+    // ensure that the output edges are of size part_counts (output without ghosts zones)
+    edges.ds_j_dt.ensure_sizes(part_counts);
+
+    const Tscal pmass = edges.gpart_mass.value;
+
+    using ComputeKernel = KernelUpdateDerivsMonofluidTVI<Tvec, SPHKernel>;
+
+    auto total_specie_count = part_counts.template map<u32>([&](u64 id, u32 count) {
+        return count * ndust;
+    });
+
+    // call the kernel for each patches with part_counts.get(id_patch) threads of patch id_patch
+    sham::distributed_data_kernel_call(
+        shamsys::instance::get_compute_scheduler_ptr(),
+        sham::DDMultiRef{
+            edges.xyz.get_spans(),
+            edges.hpart.get_spans(),
+            edges.vxyz.get_spans(),
+            edges.omega.get_spans(),
+            edges.pressure.get_spans(),
+            edges.s_j.get_spans(),
+            edges.Ttilde_sj.get_spans(),
+            edges.neigh_cache},
+        sham::DDMultiRef{edges.ds_j_dt.get_spans()},
+        total_specie_count,
+        ComputeKernel{pmass, ndust});
+}
+
+using namespace shammath;
+template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI<f64_3, M4>;
+template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI<f64_3, M6>;
+template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI<f64_3, M8>;
+
+template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI<f64_3, C2>;
+template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI<f64_3, C4>;
+template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI<f64_3, C6>;
diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp
index 809c5f7e2..6f69865fb 100644
--- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp
+++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp
@@ -26,7 +26,10 @@
 #include "shammodels/sph/math/forces.hpp"
 #include "shammodels/sph/math/mhd.hpp"
 #include "shammodels/sph/math/q_ab.hpp"
+#include "shammodels/sph/modules/ComputeDustTtilde.hpp"
+#include "shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp"
 #include "shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAV.hpp"
+#include "shammodels/sph/modules/SetDustStoppingTimeConstant.hpp"
 #include "shammodels/sph/modules/UpdateDerivs.hpp"
 #include "shamphys/mhd.hpp"
 #include "shamrock/patch/PatchDataFieldSpan.hpp"
@@ -34,12 +37,15 @@
 #include "shamrock/solvergraph/IFieldSpan.hpp"
 #include "shamrock/solvergraph/Indexes.hpp"
 #include "shamrock/solvergraph/ScalarEdge.hpp"
+#include <memory>
+#include <vector>
 
 template<class Tvec, template<class> class SPHKernel>
 void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs() {
 
-    Cfg_AV cfg_av   = solver_config.artif_viscosity;
-    Cfg_MHD cfg_mhd = solver_config.mhd_config;
+    Cfg_AV cfg_av       = solver_config.artif_viscosity;
+    Cfg_MHD cfg_mhd     = solver_config.mhd_config;
+    DustConfig cfg_dust = solver_config.dust_config;
 
     if (Constant *v = std::get_if<Constant>(&cfg_av.config)) {
         update_derivs_constantAV(*v);
@@ -60,6 +66,11 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs() {
     } else {
         shambase::throw_unimplemented();
     }
+
+    if (cfg_dust.has_s_j_field()) {
+        // we can do it separately because the backreaction is done only through the pressure
+        update_derivs_dust_monofluid_tvi_Sj(cfg_dust);
+    }
 }
 
 template<class Tvec, template<class> class SPHKernel>
@@ -1051,6 +1062,134 @@ void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(
     });
 }
 
+template<class Tvec, template<class> class SPHKernel>
+void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_dust_monofluid_tvi_Sj(
+    DustConfig cfg) {
+    StackEntry stack_loc{};
+
+    using namespace shamrock;
+    using namespace shamrock::patch;
+
+    PatchDataLayerLayout &pdl = scheduler().pdl_old();
+
+    const u32 ixyz     = pdl.get_field_idx<Tvec>("xyz");
+    const u32 ivxyz    = pdl.get_field_idx<Tvec>("vxyz");
+    const u32 iaxyz    = pdl.get_field_idx<Tvec>("axyz");
+    const u32 ihpart   = pdl.get_field_idx<Tscal>("hpart");
+    const u32 is_j     = pdl.get_field_idx<Tscal>("s_j");
+    const u32 ids_j_dt = pdl.get_field_idx<Tscal>("ds_j_dt");
+
+    shamrock::patch::PatchDataLayerLayout &ghost_layout
+        = shambase::get_check_ref(storage.ghost_layout.get());
+    u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
+    u32 ivxyz_interf  = ghost_layout.get_field_idx<Tvec>("vxyz");
+    u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
+    u32 is_j_interf   = ghost_layout.get_field_idx<Tscal>("s_j");
+
+    u32 ndust = cfg.get_dust_nvar();
+
+    auto &merged_xyzh                                 = storage.merged_xyzh.get();
+    shamrock::solvergraph::Field<Tscal> &omega        = shambase::get_check_ref(storage.omega);
+    shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
+
+    auto &part_counts            = storage.part_counts;
+    auto &part_counts_with_ghost = storage.part_counts_with_ghost;
+    auto &xyz_refs               = storage.positions_with_ghosts;
+    auto &pressure_field         = storage.pressure;
+
+    std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> vxyz_refs
+        = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("vxyz", "v");
+    {
+        shambase::get_check_ref(vxyz_refs).set_refs(
+            mpdats.map<std::reference_wrapper<PatchDataField<Tvec>>>(
+                [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
+                    return std::ref(mpdat.get_field<Tvec>(ivxyz_interf));
+                }));
+    }
+
+    std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_refs
+        = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("hpart", "h");
+    { // if was just reset before this call
+        shambase::get_check_ref(hpart_refs)
+            .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
+                [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
+                    return std::ref(mpdat.get_field<Tscal>(ihpart_interf));
+                }));
+    }
+
+    std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> omega_refs
+        = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("omega", "omega");
+    {
+        shambase::get_check_ref(omega_refs)
+            .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
+                [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
+                    return std::ref(mpdat.get_field<Tscal>(iomega_interf));
+                }));
+    }
+
+    // s_j_interf
+    std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> s_j_refs
+        = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("s_j", "s_j");
+    {
+        shambase::get_check_ref(s_j_refs).set_refs(
+            mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
+                [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
+                    return std::ref(mpdat.get_field<Tscal>(is_j_interf));
+                }));
+    }
+
+    std::shared_ptr<shamrock::solvergraph::Field<Tscal>> t_j_field
+        = std::make_shared<shamrock::solvergraph::Field<Tscal>>(ndust, "t_j", "t_j");
+
+    std::shared_ptr<shamrock::solvergraph::ScalarEdge<std::vector<Tscal>>> input_t_j
+        = std::make_shared<shamrock::solvergraph::ScalarEdge<std::vector<Tscal>>>("", "");
+    input_t_j->value = cfg.stopping_times;
+
+    std::shared_ptr<SetDustStoppingTimeConstant<Tvec, SPHKernel>> node_set_tj
+        = std::make_shared<SetDustStoppingTimeConstant<Tvec, SPHKernel>>(ndust);
+    {
+        node_set_tj->set_edges(input_t_j, part_counts_with_ghost, t_j_field);
+    }
+    node_set_tj->evaluate();
+
+    std::shared_ptr<shamrock::solvergraph::Field<Tscal>> Ttilde_sj_field
+        = std::make_shared<shamrock::solvergraph::Field<Tscal>>(ndust, "Ttilde_sj", "Ttilde_sj");
+
+    shamrock::solvergraph::SolverGraph &solver_graph = storage.solver_graph;
+    auto ds_j_dt_refs
+        = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("ds_j_dt");
+
+    auto gpart_mass
+        = solver_graph.get_edge_ptr<shamrock::solvergraph::ScalarEdge<Tscal>>("gpart_mass");
+
+    std::shared_ptr<ComputeDustTtilde<Tvec, SPHKernel>> node_tj
+        = std::make_shared<ComputeDustTtilde<Tvec, SPHKernel>>(ndust);
+    {
+        node_tj->set_edges(
+            gpart_mass, part_counts_with_ghost, hpart_refs, s_j_refs, t_j_field, Ttilde_sj_field);
+    }
+    node_tj->evaluate();
+
+    std::shared_ptr<NodeUpdateDerivsMonofluidTVI<Tvec, SPHKernel>> node
+        = std::make_shared<NodeUpdateDerivsMonofluidTVI<Tvec, SPHKernel>>(ndust);
+    {
+        node->set_edges(
+            gpart_mass,
+            part_counts,
+            part_counts_with_ghost,
+            xyz_refs,
+            hpart_refs,
+            vxyz_refs,
+            omega_refs,
+            pressure_field,
+            s_j_refs,
+            Ttilde_sj_field,
+            storage.neigh_cache,
+            ds_j_dt_refs);
+    }
+    node->evaluate();
+}
+
 using namespace shammath;
 template class shammodels::sph::modules::UpdateDerivs<f64_3, M4>;
 template class shammodels::sph::modules::UpdateDerivs<f64_3, M6>;
diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp
index 74da727eb..dafeb57de 100644
--- a/src/shammodels/sph/src/pySPHModel.cpp
+++ b/src/shammodels/sph/src/pySPHModel.cpp
@@ -39,6 +39,7 @@
 #include <pybind11/pytypes.h>
 #include <memory>
 #include <random>
+#include <utility>
 
 template<class Tvec, template<class> class SPHKernel>
 void add_instance(py::module &m, std::string name_config, std::string name_model) {
@@ -270,6 +271,11 @@ void add_instance(py::module &m, std::string name_config, std::string name_model
             },
             py::kw_only(),
             py::arg("ndust"))
+        .def(
+            "set_dust_stopping_times",
+            [](TConfig &self, std::vector<Tscal> ts) {
+                self.dust_config.stopping_times = ts;
+            })
         .def("add_ext_force_point_mass", &TConfig::add_ext_force_point_mass)
         .def(
             "add_ext_force_lense_thirring",
@@ -767,8 +773,32 @@ void add_instance(py::module &m, std::string name_config, std::string name_model
                         "unknown field type");
                 }
             })
-        .def("set_field_value_lambda_f64_3", &T::template set_field_value_lambda<f64_3>)
-        .def("set_field_value_lambda_f64", &T::template set_field_value_lambda<f64>)
+        .def(
+            "set_field_value_lambda_f64",
+            [](T &self,
+               std::string field_name,
+               const std::function<f64(Tvec)> pos_to_val,
+               const i32 offset) {
+                return self.template set_field_value_lambda<f64>(
+                    std::move(field_name), pos_to_val, offset);
+            },
+            py::arg("field_name"),
+            py::arg("pos_to_val"),
+            py::arg("offset") = 0)
+        .def(
+            "set_field_value_lambda_f64_3",
+            [](T &self,
+               std::string field_name,
+               const std::function<f64_3(Tvec)> pos_to_val,
+               const i32 offset) {
+                return self.template set_field_value_lambda<f64_3>(
+                    std::move(field_name), pos_to_val, offset);
+            },
+            py::arg("field_name"),
+            py::arg("pos_to_val"),
+            py::arg("offset") = 0)
+        .def("overwrite_field_value_f64", &T::template overwrite_field_value<f64>)
+        .def("overwrite_field_value_f64_3", &T::template overwrite_field_value<f64_3>)
         .def("remap_positions", &T::remap_positions)
         //.def("set_field_value_lambda_f64_3",[](T&self,std::string field_name, const
         // std::function<f64_3 (Tscal, Tscal , Tscal)> pos_to_val){