[SPH] add Dust TVA CFL

[tdavidcl]

No description provided.

[gemini-code-assist]

Code Review

This pull request introduces a new ComputeCFLDust1Fluid module to calculate CFL conditions for dust in a 1-fluid model, registers the Ts_j field in the solver, and updates the run_dustysettle_tvi.py example to use periodic HCP box calculations and compute L2 errors. Feedback on these changes highlights a critical compilation error in ComputeCFLDust1Fluid.hpp caused by variable shadowing of Ts_j, commented-out CFL dust computation code in Solver.cpp that should be enabled, and an unprofessional string literal exception. Additionally, there are several instances of leftover debugging code and print statements in NodeUpdateDerivsMonofluidTVI.cpp and run_dustysettle_tvi.py that should be removed before merging.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

I am having trouble creating individual review comments. Click here to see my feedback.

src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDust1Fluid.hpp (107-108)

The local variable Ts_j shadows the parameter pointer Ts_j of the lambda. This leads to a compilation error because Ts_j on the right-hand side of the assignment is resolved to the local variable itself (which is a scalar Tscal and cannot be subscripted). Rename the local variable to avoid shadowing and fix the compilation error.

                    Tscal Ts_val = Ts_j[id_a_d + j];
                    Tscal Ts2_j  = Ts_val * Ts_val;

src/shammodels/sph/src/Solver.cpp (2703-2712)

The code to instantiate and set up ComputeCFLDust1Fluid is commented out, meaning the dust CFL condition is not being computed or applied. Furthermore, the commented-out code references s_j_refs and Ts_j_refs which are not defined in this scope. To enable this feature, uncomment the block, retrieve the edges from the solver graph, and evaluate the node.

            std::shared_ptr<ComputeCFLDust1Fluid<Tscal>> compute_cfl_dust1_fluid;
            if (solver_config.dust_config.has_s_j_field()) {
                u32 ndust = solver_config.dust_config.get_dust_nvar();

                compute_cfl_dust1_fluid = std::make_shared<ComputeCFLDust1Fluid<Tscal>>(ndust);

                auto s_j_refs = storage.solver_graph.template get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("s_j");
                auto Ts_j_refs = storage.solver_graph.template get_edge_ptr<shamrock::solvergraph::Field<Tscal>>("Ts_j");

                compute_cfl_dust1_fluid->set_edges(
                    storage.part_counts, C_cour_edge, hpart_refs, storage.soundspeed, s_j_refs,
                    Ts_j_refs, cfl_dt);
            }

src/shammodels/sph/src/Solver.cpp (1595-1599)

Throwing a string literal is bad practice in C++ and the message is unprofessional. It should throw a standard exception class like std::runtime_error with a descriptive message.

        if (std::holds_alternative<None>(cfg.dust_drag_mode)) {
            throw std::runtime_error("Dust drag mode 'None' is not supported when s_j field is enabled.");
        }

src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp (102-111)

This is leftover debugging code with hardcoded particle/dust indices and printf statements. It should be removed from the production codebase.

src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp (126-144)

These blocks are leftover debugging code with hardcoded conditions and printf statements. They should be removed from the production codebase.

examples/sph/run_dustysettle_tvi.py (773-774)

Setting reference_dusty_settle = None immediately after instantiating it disables the reference model and the subsequent L2 error analysis. This appears to be leftover debugging code and should be removed.

        reference_dusty_settle = ReferenceDustySettleAll()

[github-actions]

Thanks @tdavidcl for opening this PR!

You can do multiple things directly here:
1 - Comment pre-commit.ci run to run pre-commit checks.
2 - Comment pre-commit.ci autofix to apply fixes.
3 - Add label autofix.ci to fix authorship & pre-commit for every commit made.
4 - Add label light-ci to only trigger a reduced & faster version of the CI (need the full one before merge).
5 - Add label trigger-ci to create an empty commit to trigger the CI.

Once the workflow completes a message will appear displaying informations related to the run.

Also the PR gets automatically reviewed by gemini, you can:
1 - Comment /gemini review to trigger a review
2 - Comment /gemini summary for a summary
3 - Tag it using @gemini-code-assist either in the PR or in review comments on files

[github-actions]

Workflow report

workflow report corresponding to commit 512c7f7
Commiter email is timothee.davidcleris@proton.me

Pre-commit check report

Some failures were detected in base source checks checks.
Check the On PR / Linting / Base source checks (pull_request) job in the tests for more detailed output

❌ ruff-check

�[1m�[91mNPY201 �[0m�[1m`np.trapz` will be removed in NumPy 2.0. Use `numpy.trapezoid` on NumPy 2.0, or ignore this warning on earlier versions.�[0m
   �[1m�[94m-->�[0m examples/sph/run_dustysettle_tvi.py:568:23
    �[1m�[94m|�[0m
�[1m�[94m566 |�[0m         L2_integral = np.trapezoid(L2_func, z_ref)  # NumPy >= 2.0
�[1m�[94m567 |�[0m     else:
�[1m�[94m568 |�[0m         L2_integral = np.trapz(L2_func, z_ref)  # NumPy < 2.0
    �[1m�[94m|�[0m                       �[1m�[91m^^^^^^^^�[0m
�[1m�[94m569 |�[0m
�[1m�[94m570 |�[0m     return np.sqrt(L2_integral)
    �[1m�[94m|�[0m
�[1m�[96mhelp�[0m: �[1mReplace with `numpy.trapezoid` (requires NumPy 2.0 or greater)�[0m

Found 1 error.
No fixes available (1 hidden fix can be enabled with the `--unsafe-fixes` option).

Suggested changes

Detailed changes :

diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDustDeltav.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDustDeltav.hpp
index 82bc3369..31202b3c 100644
--- a/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDustDeltav.hpp
+++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeCFLDustDeltav.hpp
@@ -28,7 +28,7 @@
     X_RO(shamrock::solvergraph::Indexes<u32>, part_counts)                                         \
     X_RO(shamrock::solvergraph::ScalarEdge<Tscal>, C_delta_v)                                      \
     X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, hpart)                                          \
-    X_RO(shamrock::solvergraph::IFieldSpan<Tvec>, delta_v)                                     \
+    X_RO(shamrock::solvergraph::IFieldSpan<Tvec>, delta_v)                                         \
     X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, cfl_dt)
 
 template<class Tvec>
@@ -56,20 +56,17 @@ class ComputeCFLDustDeltav : public shamrock::solvergraph::INode {
             sham::DDMultiRef{edges.cfl_dt.get_spans()},
             edges.part_counts.indexes,
             [C_delta_v, nbins = this->nbins](
-                u32 id_a,
-                const Tscal *hpart,
-                const Tvec *delta_v,
-                Tscal *cfl_dt) {
+                u32 id_a, const Tscal *hpart, const Tvec *delta_v, Tscal *cfl_dt) {
                 u32 id_a_d = id_a * nbins;
 
-                Tscal h_a   = hpart[id_a];
+                Tscal h_a = hpart[id_a];
 
                 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];
+                    Tvec delta_v_j_a       = delta_v[id_a_d + j];
                     Tscal delta_v_j_a_norm = sycl::length(delta_v_j_a);
-                    cfl_tmp = sycl::min(cfl_tmp, h_a / delta_v_j_a_norm);
+                    cfl_tmp                = sycl::min(cfl_tmp, h_a / delta_v_j_a_norm);
                 }
 
                 cfl_tmp *= C_delta_v;