[Examples] dustydisc setup#1868
Conversation
|
Thanks @tdavidcl for opening this PR! You can do multiple things directly here: Once the workflow completes a message will appear displaying informations related to the run. Also the PR gets automatically reviewed by gemini, you can: |
Workflow reportworkflow report corresponding to commit 8cd27ba Pre-commit check reportSome failures were detected in base source checks checks. ❌ ruff-formatSuggested changesDetailed changes :diff --git a/examples/sph/run_dustydisc.py b/examples/sph/run_dustydisc.py
index 946607d6..1c5638cd 100644
--- a/examples/sph/run_dustydisc.py
+++ b/examples/sph/run_dustydisc.py
@@ -500,7 +500,6 @@ if ndust > 0:
"extra_title": f"[$s_{{grain}}$ = {grain_size_si[jdust]:.2e} m]",
}
-
dust_slice_epsilon_plot = []
for jdust in range(ndust):
@@ -831,8 +830,6 @@ class radial_profile_plot:
self.profile_plot.render_all(self.plot_func)
-
-
class vert_slices_plots:
def __init__(self):
self.profile_plot = AnalysisHelper(
@@ -1043,6 +1040,7 @@ class vert_slices_plots:
def render_all(self):
self.profile_plot.render_all(self.plot_func)
+
if ndust > 0:
rad_plot = radial_profile_plot()
vert_plot = vert_slices_plots()
|
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request introduces a new dusty SPH disc simulation example (run_dustydisc.py), centralizes matplotlib styling across examples via a new shamrock.matplotlib module, and adds dust-specific analysis and plotting helpers to the library. The review feedback highlights several critical issues: potential NameError and division-by-zero warnings in run_dustydisc.py when ndust = 0, console flooding from unprotected print statements in MPI parallel runs, and a potential division-by-zero in DensityPlots.py when gas density is zero.
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.
| for jdust, p in enumerate(dust_column_density_plot): | ||
| p.render_all( | ||
| **p.render_args, | ||
| ) | ||
|
|
||
|
|
||
| vertical_density_plot.render_all( | ||
| **vertical_density_plot.render_args, | ||
| ) | ||
|
|
||
| for jdust, p in enumerate(dust_slice_density_plot): | ||
| p.render_all( | ||
| **p.render_args, | ||
| ) | ||
|
|
||
| for jdust, p in enumerate(dust_slice_epsilon_plot): | ||
| p.render_all( | ||
| **p.render_args, | ||
| ) |
There was a problem hiding this comment.
These plotting loops will raise a NameError when ndust = 0 (the default value) because dust_column_density_plot, dust_slice_density_plot, and dust_slice_epsilon_plot are only defined inside if ndust > 0: blocks. Wrapping these loops in if ndust > 0: guards prevents this crash.
if ndust > 0:
for jdust, p in enumerate(dust_column_density_plot):
p.render_all(
**p.render_args,
)
vertical_density_plot.render_all(
**vertical_density_plot.render_args,
)
if ndust > 0:
for jdust, p in enumerate(dust_slice_density_plot):
p.render_all(
**p.render_args,
)
for jdust, p in enumerate(dust_slice_epsilon_plot):
p.render_all(
**p.render_args,
)| mrn_weight = grain_size ** (4 - mrn_pow) | ||
| mrn_weight *= grain_size_si < mrn_cutoff_si | ||
| mrn_weight = mrn_weight / np.sum(mrn_weight) | ||
|
|
||
| print(f"mrn_weight = {mrn_weight}") |
There was a problem hiding this comment.
When ndust = 0, grain_size is empty, which causes np.sum(mrn_weight) to be 0.0 and results in a division-by-zero RuntimeWarning. Wrapping this calculation in if ndust > 0: and protecting the print statement for MPI parallel runs prevents this warning and avoids console spam.
| mrn_weight = grain_size ** (4 - mrn_pow) | |
| mrn_weight *= grain_size_si < mrn_cutoff_si | |
| mrn_weight = mrn_weight / np.sum(mrn_weight) | |
| print(f"mrn_weight = {mrn_weight}") | |
| if ndust > 0: | |
| mrn_weight = grain_size ** (4 - mrn_pow) | |
| mrn_weight *= grain_size_si < mrn_cutoff_si | |
| mrn_weight = mrn_weight / np.sum(mrn_weight) | |
| if shamrock.sys.world_rank() == 0: | |
| print(f"mrn_weight = {mrn_weight}") | |
| else: | |
| mrn_weight = np.array([]) |
| epsilon_target = epsilon_base * mrn_weight[j] | ||
| print(f"epsilon_target = {epsilon_target} {j}") | ||
| s = np.sqrt(rho * epsilon_target) | ||
|
|
||
| print( | ||
| f"s = {s} {np.isnan(s).any()} epsilon_target = {epsilon_target} mrn_weight = {mrn_weight[j]}, rho = {rho}" | ||
| ) | ||
|
|
||
| return s |
There was a problem hiding this comment.
These print statements inside compute_sj_new_j are executed for every patch on every MPI rank during initialization. This can severely flood the console and degrade setup performance. They should be removed.
| epsilon_target = epsilon_base * mrn_weight[j] | |
| print(f"epsilon_target = {epsilon_target} {j}") | |
| s = np.sqrt(rho * epsilon_target) | |
| print( | |
| f"s = {s} {np.isnan(s).any()} epsilon_target = {epsilon_target} mrn_weight = {mrn_weight[j]}, rho = {rho}" | |
| ) | |
| return s | |
| epsilon_target = epsilon_base * mrn_weight[j] | |
| s = np.sqrt(rho * epsilon_target) | |
| return s |
| print(f"grains sizes = {grain_size_si_edges} [m]") | ||
| print(f"grains dens = {rho_grains_si_edges} [kg.m^-3]") | ||
|
|
||
| grain_size_edges = grain_size_si_edges * codeu.get("m") | ||
| rho_grains_edges = codeu.get("kg") * codeu.get("m", power=-3) * np.array(rho_grains_si_edges) | ||
|
|
||
| print(f"grains sizes = {grain_size_edges} [code u]") | ||
| print(f"grains dens = {rho_grains_edges} [code u]") | ||
|
|
||
| grain_size = np.sqrt(grain_size_edges[:-1] * grain_size_edges[1:]) | ||
| rho_grains = np.sqrt(rho_grains_edges[:-1] * rho_grains_edges[1:]) | ||
|
|
||
| grain_size_si = np.sqrt(grain_size_si_edges[:-1] * grain_size_si_edges[1:]) | ||
| rho_grains_si = np.sqrt(rho_grains_si_edges[:-1] * rho_grains_si_edges[1:]) | ||
|
|
||
| print(f"grains sizes = {grain_size_si} [m]") | ||
| print(f"grains dens = {rho_grains_si} [kg.m^-3]") | ||
|
|
||
| print(f"grains sizes = {grain_size} [code units]") | ||
| print(f"grains dens = {rho_grains} [code units]") |
There was a problem hiding this comment.
These print statements are not protected by rank checks, which will cause duplicate output lines to be printed by every MPI process in parallel runs. Wrapping them with if shamrock.sys.world_rank() == 0: ensures clean console output.
if shamrock.sys.world_rank() == 0:
print(f"grains sizes = {grain_size_si_edges} [m]")
print(f"grains dens = {rho_grains_si_edges} [kg.m^-3]")
grain_size_edges = grain_size_si_edges * codeu.get("m")
rho_grains_edges = codeu.get("kg") * codeu.get("m", power=-3) * np.array(rho_grains_si_edges)
if shamrock.sys.world_rank() == 0:
print(f"grains sizes = {grain_size_edges} [code u]")
print(f"grains dens = {rho_grains_edges} [code u]")
grain_size = np.sqrt(grain_size_edges[:-1] * grain_size_edges[1:])
rho_grains = np.sqrt(rho_grains_edges[:-1] * rho_grains_edges[1:])
grain_size_si = np.sqrt(grain_size_si_edges[:-1] * grain_size_si_edges[1:])
rho_grains_si = np.sqrt(rho_grains_si_edges[:-1] * rho_grains_si_edges[1:])
if shamrock.sys.world_rank() == 0:
print(f"grains sizes = {grain_size_si} [m]")
print(f"grains dens = {rho_grains_si} [kg.m^-3]")
print(f"grains sizes = {grain_size} [code units]")
print(f"grains dens = {rho_grains} [code units]")| tmp = rhod / rho | ||
|
|
||
| print(tmp.min(), tmp.max()) | ||
| print(tmp) | ||
| return tmp |
There was a problem hiding this comment.
The raw division rhod / rho can trigger division-by-zero warnings/errors in regions where density is zero. Using np.divide with a where guard handles this robustly. Additionally, the leftover debug print statements should be removed to keep the library output clean.
| tmp = rhod / rho | |
| print(tmp.min(), tmp.max()) | |
| print(tmp) | |
| return tmp | |
| tmp = np.divide(rhod, rho, out=np.zeros_like(rhod), where=rho > 0) | |
| return tmp |
1 participant
No description provided.