[SPH] sheared wave setup#1897
Draft
tdavidcl wants to merge 1 commit into
Draft
Conversation
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
Contributor
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.
Comment on lines
+116
to
+266
| # %% | ||
| # 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() |
Contributor
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.
Comment on lines
+64
to
+70
| 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) |
Contributor
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.
Suggested change
| 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) |
Comment on lines
+88
to
+91
| #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) |
Contributor
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.
Contributor
|
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: |
Contributor
Workflow reportworkflow report corresponding to commit a1d7e1c Pre-commit check reportSome failures were detected in base source checks checks. ❌ ruff-check❌ ruff-formatSuggested changesDetailed changes :diff --git a/examples/sph/run_sph_sheared_wave.py b/examples/sph/run_sph_sheared_wave.py
index f8fdc9d2..8ba4605e 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,7 +203,6 @@ model.set_cfl_force(0.25)
from math import exp
import matplotlib.pyplot as plt
-import numpy as np
def plot(iplot):
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.