[Stagging][SPH] Dust paper plots

[tdavidcl]

No description provided.

[coderabbitai]

Important

Review skipped

Draft detected.

Please check the settings in the CodeRabbit UI or the .coderabbit.yaml file in this repository. To trigger a single review, invoke the @coderabbitai review command.

⚙️ Run configuration

Configuration used: Repository UI

Review profile: ASSERTIVE

Plan: Pro Plus

Run ID: 44acd67f-36ad-4929-aa64-587f46b5b186

You can disable this status message by setting the reviews.review_status to false in the CodeRabbit configuration file.

Use the checkbox below for a quick retry:

• 🔍 Trigger review

✨ Finishing Touches

🧪 Generate unit tests (beta)

• Create PR with unit tests

---

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

• X
• Mastodon
• Reddit
• LinkedIn

_{Comment @coderabbitai help to get the list of available commands.}

[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

[tdavidcl]

@claude review this PR

[claude]

Claude finished @tdavidcl's task in 2m 15s —— View job

---

Dust Paper Plots PR Review

Summary

This PR refines 4 dusty hydrodynamics example scripts to produce publication-quality figures. It adds a new sympy-based dispersion relation plotter, improves lattice setup for dustywave, extends the dustydiffuse simulation time, and fixes the Sod dust initialization.

---

Todo

• Review all 4 modified files
• Check for bugs and correctness issues
• Check style consistency with project conventions
• Identify improvements

---

examples/physics/run_dustywave_sympy.py (NEW FILE — +75 lines)

New file — derives the dustywave TVA dispersion relation via SymPy and produces two figure layouts.

1. get_overroots_DCL26_simple has a dead-term coefficient (line 110)

   _r1 = +cs * sp.sqrt(1 - eps) * k + i * k**2 * cs**2 * eps * ts * (-1 + 1) / 2

   The (-1 + 1) term evaluates to zero, making the imaginary part of _r1 vanish entirely. This looks intentional for a "simple" approximation but is confusing as written. Consider * 0 or a comment like # imaginary part drops out in simple approx.

2. get_overroots_DCL26 prints during computation (lines 126-127)

   print(sp.latex(sp.Abs(_r1)))
   print(sp.latex(sp.Abs(_r2)))

   These prints will fire every time get_overroots_DCL26 is called (twice per plot_case). They're likely leftover from development — consider removing or gating behind a debug flag.

3. Duplicate computation in plot_case (lines 191-198 and 248-256)
   The absolute value computations (r1_vals_abs, r2_vals_abs, etc.) are computed twice — once for the 2×2 figure and again for the 2×1 figure. Extract to a helper function or compute once and reuse.

4. _to_trash/ output directory — the file writes to _to_trash/ but doesn't create it (unlike the other scripts that call os.makedirs("_to_trash", exist_ok=True)). This will crash if the directory doesn't exist. Fix: add os.makedirs("_to_trash", exist_ok=True) before the first plt.savefig.

---

examples/sph/run_dustywave_tva.py (+77 / -32)

5. nwave type change (line 307)
   Changed from nwave = 2 (int) to nwave = 2.0 (float). This is fine since it's used in range(int(Twave_cnt * nwave)) where the int cast handles it. The float is arguably more intentional for a physical quantity, but the int cast is still there. Consider: nwave_steps = int(Twave_cnt * nwave) at the top to make the type clearer.

6. Missing plt.close() after final multi-panel figure (line 508)
   The final combined figure dustywave_tva_scan_all.png/pdf calls plt.close() at the end, but the animation section above (lines 476-491) calls plt.show() after each ani.save() — this is fine for interactive use but note that in Sphinx/gallery contexts these plt.show() calls are no-ops. The existing sphinx_gallery_multi_image setting handles this.

7. Hardcoded LZ env var with default 18 (line 38)

   lx = int(os.environ.get("LZ", 18))

   Inconsistent with the other files that use RESOL as the env var name. Consider standardizing on one env var name across all dusty examples, or documenting the convention.

---

examples/sph/run_dustydiffuse_tva.py (+10 / -4)

8. Extended simulation range — range(20) → range(31) (t=0 to t=3.0)
   Good improvement for capturing the full diffusion behavior.

9. plt.figure(dpi=250) without plt.close() (line 298)
   The final figure uses plt.figure(dpi=250) and ends with plt.close(), which is correct. However, note that plt.figure() creates a new figure on top of any remaining open figures. The animation section above leaves figures in an unknown state. Consider adding plt.clf() or plt.close('all') before the final figure section.

10. rasterized=True on particle plots — good addition for PDF output quality.

---

examples/sph/run_sod_dust_tva.py (+28 / -18)

11. ts changed from 0.001 to 1e-4 (line 38)
    Significant physics change — 10× smaller drag time constant. This is a physics parameter change, so verify this is the intended physical regime. If intentional, add a comment explaining why.

12. RESOL default changed from 128 to 32 (line 52)
    4× fewer particles. Good for CI speed, but verify the L2 convergence still holds at this resolution. The previous PR ([SPH] add dusty TVA Sod example and improve Sod L2 analysis #1939) already improved L2 analysis — this should be fine but worth confirming.

13. s_j initialization refactoring (lines 221-234)
    The old code used set_value_in_a_box for gas and dust regions separately. The new code computes rho from particle masses and kernel densities, then sets s_j = sqrt(rho * eps). This is more physically correct (density varies spatially), but it means the initial s_j is now position-dependent rather than uniform. Verify that the Sod initial condition still produces the expected shock structure with this spatially-varying s_j.

14. Missing plt.close() on L2 and dust mass diagnostic plots (lines 273, 294)
    These plots call plt.savefig() and plt.show() but no plt.close(). Not a crash risk, but could accumulate open figure handles in long-running contexts. Add plt.close() after each for consistency.

15. rasterized=True on all scatter plots — good for publication-quality PDF output.

---

Style Notes

16. # ## title separator pattern — consistent across all files (e.g., line 2-3: ===================================). Good.

17. _to_trash/ directory — used consistently for intermediate outputs. Consider documenting this convention in a shared utility or at least a one-line comment.

18. No if __name__ == "__main__" guard — consistent with the project's existing script style. Good.

---

Critical Items (must fix before merge)

#	File	Issue
4	run_dustywave_sympy.py	Missing os.makedirs("_to_trash") — will crash
11	run_sod_dust_tva.py	ts physics change needs documentation
13	run_sod_dust_tva.py	s_j initialization change may alter Sod solution — verify

Recommendations (nice to have)

#	File	Issue
2	run_dustywave_sympy.py	Remove debug print(sp.latex(...)) calls
3	run_dustywave_sympy.py	Deduplicate absolute value computations
15	run_sod_dust_tva.py	Add plt.close() on diagnostic plots

[github-actions]

Workflow report

workflow report corresponding to commit 416c343
Commiter email is 37929162+mergify[bot]@users.noreply.github.com
You are using github private e-mail. This prevent proper tracing of who contributed what, please disable it (see Keep my email addresses private).

Light CI is enabled (the default for pull requests). This will only run the basic tests and not the full tests.
Full CI runs if the full-ci label is set, or automatically on Mergify merge-queue branches (mergify/merge-queue/*).
The merge gate job "on PR / all" is skipped in this case. Queue entry uses "on PR / all_light"; full CI runs in the merge queue.

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�[91mPLW1508�[0m�[1m Invalid type for environment variable default; expected `str` or `None`�[0m
  �[1m�[94m--> �[0mexamples/sph/run_dustydisc.py:45:37
   �[1m�[94m|�[0m
�[1m�[94m44�[0m �[1m�[94m|�[0m # Resolution
�[1m�[94m45�[0m �[1m�[94m|�[0m Npart = int(os.environ.get("NPART", 100000))
   �[1m�[94m|�[0m                                     �[1m�[91m^^^^^^�[0m
�[1m�[94m46�[0m �[1m�[94m|�[0m print(f"Npart = {Npart} (NPART={os.environ.get('NPART', 'not set')!r})")
   �[1m�[94m|�[0m

�[1m�[91mPLW1508�[0m�[1m Invalid type for environment variable default; expected `str` or `None`�[0m
  �[1m�[94m--> �[0mexamples/sph/run_dustydisc.py:88:37
   �[1m�[94m|�[0m
�[1m�[94m87�[0m �[1m�[94m|�[0m # Dust parameters
�[1m�[94m88�[0m �[1m�[94m|�[0m ndust = int(os.environ.get("NDUST", 5))
   �[1m�[94m|�[0m                                     �[1m�[91m^�[0m
�[1m�[94m89�[0m �[1m�[94m|�[0m gamma = 1.4
   �[1m�[94m|�[0m

�[1m�[91mPLW1508�[0m�[1m Invalid type for environment variable default; expected `str` or `None`�[0m
  �[1m�[94m--> �[0mexamples/sph/run_dustydisc_perftest.py:55:41
   �[1m�[94m|�[0m
�[1m�[94m54�[0m �[1m�[94m|�[0m     # Resolution
�[1m�[94m55�[0m �[1m�[94m|�[0m     Npart = int(os.environ.get("NPART", 100000))
   �[1m�[94m|�[0m                                         �[1m�[91m^^^^^^�[0m
�[1m�[94m56�[0m �[1m�[94m|�[0m     print(f"Npart = {Npart} (NPART={os.environ.get('NPART', 'not set')!r})")
   �[1m�[94m|�[0m

Found 3 errors.

Suggested changes

Detailed changes :