Conversation
There was a problem hiding this comment.
Code Review
This pull request introduces a new example script, run_sph_sheared_wave.py, for running an unstratified shearing box simulation in SPH, and updates BCConfig.hpp to parse shear_speed instead of speed from the JSON configuration. The review feedback highlights a large block of redundant, broken procedural code at the end of the example script that should be removed. Additionally, the reviewer suggests scaling the wave phase in vel_func to ensure perfect periodicity across adjusted box boundaries, and enabling the commented-out shearing box boundary conditions and external forces to align with the example's intended purpose.
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.
| # %% | ||
| # Setup parameters | ||
|
|
||
|
|
||
|
|
||
| render_gif = True | ||
|
|
||
|
|
||
| import os | ||
|
|
||
| # Create the dump directory if it does not exist | ||
| if shamrock.sys.world_rank() == 0: | ||
| os.makedirs(dump_folder, exist_ok=True) | ||
|
|
||
| # %% | ||
| # Generate the config & init the scheduler | ||
| 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_boundary_shearing_periodic((1, 0, 0), (0, 1, 0), shear_speed) | ||
| cfg.set_eos_adiabatic(gamma) | ||
| cfg.add_ext_force_shearing_box(Omega_0=Omega_0, eta=eta, q=q) | ||
| cfg.set_units(shamrock.UnitSystem()) | ||
| cfg.print_status() | ||
| model.set_solver_config(cfg) | ||
|
|
||
| model.init_scheduler(int(1e7), 1) | ||
|
|
||
| model.resize_simulation_box(bmin, bmax) | ||
|
|
||
|
|
||
| # %% | ||
| # Add the particles & set fields values | ||
| # Note that every field that are not mentioned are set to zero | ||
| model.add_cube_fcc_3d(dr, bmin, bmax) | ||
|
|
||
| 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_value_in_a_box("uint", "f64", 1, bmin, bmax) | ||
| # model.set_value_in_a_box("vxyz","f64_3", (-10,0,0) , bmin,bmax) | ||
|
|
||
| pen_sz = 0.1 | ||
|
|
||
| mm = 1 | ||
| MM = 0 | ||
|
|
||
|
|
||
| def vel_func(r): | ||
| global mm, MM | ||
| x, y, z = r | ||
|
|
||
| s = (x - (xM + xm) / 2) / (xM - xm) | ||
| vel = (shear_speed) * s | ||
|
|
||
| mm = min(mm, vel) | ||
| MM = max(MM, vel) | ||
|
|
||
| return (0, vel, 0.0) | ||
| # return (1,0,0) | ||
|
|
||
|
|
||
| model.set_field_value_lambda_f64_3("vxyz", vel_func) | ||
| # print("Current part mass :", pmass) | ||
| model.set_particle_mass(pmass) | ||
|
|
||
|
|
||
| tot_u = pmass * model.get_sum("uint", "f64") | ||
| # print("total u :",tot_u) | ||
|
|
||
| print(f"v_shear = {shear_speed} | dv = {MM - mm}") | ||
|
|
||
|
|
||
| model.set_cfl_cour(0.3) | ||
| model.set_cfl_force(0.25) | ||
|
|
||
| # %% | ||
| # Perform the plot | ||
|
|
||
| from math import exp | ||
|
|
||
| import matplotlib.pyplot as plt | ||
| import numpy as np | ||
|
|
||
|
|
||
| def plot(iplot): | ||
| dic = ctx.collect_data() | ||
| fig, axs = plt.subplots(2, 1, figsize=(5, 8), sharex=True) | ||
| fig.suptitle("t = {:.2f}".format(model.get_time())) | ||
| axs[0].scatter(dic["xyz"][:, 0], dic["xyz"][:, 1], s=1) | ||
| axs[1].scatter(dic["xyz"][:, 0], dic["vxyz"][:, 1], s=1) | ||
|
|
||
| axs[0].set_ylabel("y") | ||
| axs[1].set_ylabel("vy") | ||
| axs[1].set_xlabel("x") | ||
|
|
||
| axs[0].set_xlim(xm - 0.1, xM + 0.1) | ||
| axs[0].set_ylim(ym - 0.1, yM + 0.1) | ||
|
|
||
| axs[1].set_xlim(xm - 0.1, xM + 0.1) | ||
| axs[1].set_ylim(shear_speed * 0.7, -shear_speed * 0.7) | ||
|
|
||
| plt.tight_layout() | ||
| plt.savefig(os.path.join(dump_folder, f"{sim_name}_{iplot:04}.png")) | ||
| plt.close(fig) | ||
|
|
||
|
|
||
| # %% | ||
| # Performing the timestep loop | ||
| model.timestep() | ||
|
|
||
| dt_stop = 0.02 | ||
| for i in range(20): | ||
| t_target = i * dt_stop | ||
| # skip if the model is already past the target | ||
| if model.get_time() > t_target: | ||
| continue | ||
|
|
||
| model.evolve_until(i * dt_stop) | ||
|
|
||
| # Dump name is "dump_xxxx.sham" where xxxx is the timestep | ||
| model.do_vtk_dump(os.path.join(dump_folder, f"{sim_name}_{i:04}.vtk"), True) | ||
| plot(i) | ||
|
|
||
| #################################################### | ||
| # Convert PNG sequence to Image sequence in mpl | ||
| #################################################### | ||
|
|
||
| import matplotlib.animation as animation | ||
| from shamrock.utils.plot import show_image_sequence | ||
|
|
||
| # If the animation is not returned only a static image will be shown in the doc | ||
| glob_str = os.path.join(dump_folder, f"{sim_name}_*.png") | ||
| ani = show_image_sequence(glob_str, render_gif=render_gif) | ||
|
|
||
| if render_gif and shamrock.sys.world_rank() == 0: | ||
| # To save the animation using Pillow as a gif | ||
| # writer = animation.PillowWriter(fps=15, | ||
| # metadata=dict(artist='Me'), | ||
| # bitrate=1800) | ||
| # ani.save('scatter.gif', writer=writer) | ||
|
|
||
| # Show the animation | ||
| plt.show() |
There was a problem hiding this comment.
This entire block of procedural code at the bottom of the file appears to be a redundant/leftover version of the simulation setup. It is completely broken and will raise multiple NameErrors if executed because:
- Variables like
shear_speed,gamma,Omega_0,eta,q,bmin,bmax,dr,rho,xm,xM,ym,yM,zm,zMare defined as class attributes ofSimulationand are not available in the global scope. sim_nameis used but never defined anywhere in the script.
Since the simulation is already successfully configured and run using the Simulation class on lines 35-114, this entire section (lines 116-266) should be removed. If you want to keep the plotting or animation functionality, it should be refactored as a callback method inside the Simulation class.
| def vel_func(self,r): | ||
| x, y, z = r | ||
|
|
||
| s = (x - (self.xM + self.xm) / 2) / (self.xM - self.xm) | ||
| vel = (self.shear_speed) * s | ||
|
|
||
| return (self.amplitude*np.sin(2.0*np.pi*(y+4*z)), vel, 0.0) |
There was a problem hiding this comment.
The wave phase in vel_func is currently hardcoded as y + 4*z. Since the box boundaries are adjusted by get_ideal_hcp_box (which can slightly modify the box size to fit the HCP lattice), the actual box lengths in
To ensure perfect periodicity, scale the wave phase using the actual box dimensions and the wave parameters self.ny and self.nz.
| def vel_func(self,r): | |
| x, y, z = r | |
| s = (x - (self.xM + self.xm) / 2) / (self.xM - self.xm) | |
| vel = (self.shear_speed) * s | |
| return (self.amplitude*np.sin(2.0*np.pi*(y+4*z)), vel, 0.0) | |
| def vel_func(self, r): | |
| x, y, z = r | |
| s = (x - (self.xM + self.xm) / 2) / (self.xM - self.xm) | |
| vel = self.shear_speed * s | |
| phase = self.ny * (y - self.ym) / (self.yM - self.ym) + self.nz * (z - self.zm) / (self.zM - self.zm) | |
| return (self.amplitude * np.sin(2.0 * np.pi * phase), vel, 0.0) |
| #cfg.set_boundary_shearing_periodic((1, 0, 0), (0, 1, 0), self.shear_speed) | ||
| cfg.set_boundary_periodic() | ||
| cfg.set_eos_adiabatic(self.gamma) | ||
| #cfg.add_ext_force_shearing_box(Omega_0=self.Omega_0, eta=self.eta, q=self.q) |
There was a problem hiding this comment.
The shearing box boundary conditions and external forces are currently commented out, and standard periodic boundaries are used instead. Since this example is intended to demonstrate a shearing box simulation (as described in the module docstring), these configurations should be enabled. If you want to run a standard periodic wave, consider creating a separate example or parameterizing this script.
|
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: |
|
Important Draft PR not reviewedDraft PRs are not automatically reviewed by default.
To automatically review draft PRs, update your CodeRabbit configuration: reviews:
auto_review:
drafts: trueThanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
Workflow reportworkflow report corresponding to commit d696f99 Light CI is enabled (the default for pull requests). This will only run the basic tests and not the full tests. Pre-commit check reportSome failures were detected in base source checks checks. ❌ ruff-format❌ ruff-checkSuggested changesDetailed changes :diff --git a/examples/sph/run_sph_sheared_wave.py b/examples/sph/run_sph_sheared_wave.py
index f8fdc9d2..fd125410 100644
--- a/examples/sph/run_sph_sheared_wave.py
+++ b/examples/sph/run_sph_sheared_wave.py
@@ -7,9 +7,10 @@ This simple example shows how to run an unstratified shearing box simulation
# sphinx_gallery_multi_image = "single"
-import shamrock
-from shamrock.utils.SimulationRunner import SimulationRunner, callback, simulation_setup
import numpy as np
+from shamrock.utils.SimulationRunner import SimulationRunner, callback, simulation_setup
+
+import shamrock
# If we use the shamrock executable to run this script instead of the python interpreter,
# we should not initialize the system as the shamrock executable needs to handle specific MPI logic
@@ -32,9 +33,10 @@ model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4"
dump_folder = "_to_trash/sph_sheared_wave"
+
class Simulation(SimulationRunner):
# Use the global vars defined at the top of the file
- t_end = 15.
+ t_end = 15.0
dump_prefix = dump_folder + "/" + "dump"
# simulation parameters
@@ -56,18 +58,18 @@ class Simulation(SimulationRunner):
q = 3.0 / 2.0
# wave parameters
- nx,ny,nz = 0, 1, 4
+ nx, ny, nz = 0, 1, 4
amplitude = 1e-6
shear_speed = -q * Omega_0 * (xM - xm)
- def vel_func(self,r):
+ def vel_func(self, r):
x, y, z = r
s = (x - (self.xM + self.xm) / 2) / (self.xM - self.xm)
vel = (self.shear_speed) * s
- return (self.amplitude*np.sin(2.0*np.pi*(y+4*z)), vel, 0.0)
+ return (self.amplitude * np.sin(2.0 * np.pi * (y + 4 * z)), vel, 0.0)
@callback(walltime_interval=30.0) # Checkpoint the simulation every 30 seconds
def checkpoint(self, icheckpoint):
@@ -85,10 +87,10 @@ class Simulation(SimulationRunner):
cfg.set_artif_viscosity_VaryingCD10(
alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
)
- #cfg.set_boundary_shearing_periodic((1, 0, 0), (0, 1, 0), self.shear_speed)
+ # cfg.set_boundary_shearing_periodic((1, 0, 0), (0, 1, 0), self.shear_speed)
cfg.set_boundary_periodic()
cfg.set_eos_adiabatic(self.gamma)
- #cfg.add_ext_force_shearing_box(Omega_0=self.Omega_0, eta=self.eta, q=self.q)
+ # cfg.add_ext_force_shearing_box(Omega_0=self.Omega_0, eta=self.eta, q=self.q)
cfg.set_units(shamrock.UnitSystem())
cfg.print_status()
model.set_solver_config(cfg)
@@ -117,7 +119,6 @@ sim.run()
# Setup parameters
-
render_gif = True
@@ -202,13 +203,12 @@ model.set_cfl_force(0.25)
from math import exp
import matplotlib.pyplot as plt
-import numpy as np
def plot(iplot):
dic = ctx.collect_data()
fig, axs = plt.subplots(2, 1, figsize=(5, 8), sharex=True)
- fig.suptitle("t = {:.2f}".format(model.get_time()))
+ fig.suptitle(f"t = {model.get_time():.2f}")
axs[0].scatter(dic["xyz"][:, 0], dic["xyz"][:, 1], s=1)
axs[1].scatter(dic["xyz"][:, 0], dic["vxyz"][:, 1], s=1)
@@ -248,7 +248,7 @@ for i in range(20):
# Convert PNG sequence to Image sequence in mpl
####################################################
-import matplotlib.animation as animation
+from matplotlib import animation
from shamrock.utils.plot import show_image_sequence
# If the animation is not returned only a static image will be shown in the doc
|
No description provided.