From a1d7e1c17873bc2a07408287989a0c25d90d0b5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 25 Jun 2026 18:04:07 +0200 Subject: [PATCH] [SPH] sheared wave setup --- examples/sph/run_sph_sheared_wave.py | 266 ++++++++++++++++++ .../shammodels/sph/config/BCConfig.hpp | 2 +- 2 files changed, 267 insertions(+), 1 deletion(-) create mode 100644 examples/sph/run_sph_sheared_wave.py diff --git a/examples/sph/run_sph_sheared_wave.py b/examples/sph/run_sph_sheared_wave.py new file mode 100644 index 0000000000..f8fdc9d2ec --- /dev/null +++ b/examples/sph/run_sph_sheared_wave.py @@ -0,0 +1,266 @@ +""" +Shearing box in SPH +======================== + +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 + +# 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 +if not shamrock.sys.is_initialized(): + shamrock.change_loglevel(1) + shamrock.sys.init("0:0") + +# %% +# Use shamrock documentation style for matplotlib +shamrock.matplotlib.set_shamrock_mpl_style() + + +# %% +# Initialize context & attach a SPH model to it +ctx = shamrock.Context() +ctx.pdata_layout_new() + +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. + dump_prefix = dump_folder + "/" + "dump" + + # simulation parameters + gamma = 5.0 / 3.0 + rho = 1 + uint = 1 + + dr = 0.005 + bmin = (-0.5, -0.5, -0.5) + bmax = (0.5, 0.5, 0.5) + pmass = -1 + + bmin, bmax = shamrock.math.get_ideal_hcp_box(dr, bmin, bmax) + xm, ym, zm = bmin + xM, yM, zM = bmax + + Omega_0 = 0 + eta = 0.00 + q = 3.0 / 2.0 + + # wave parameters + nx,ny,nz = 0, 1, 4 + amplitude = 1e-6 + + shear_speed = -q * Omega_0 * (xM - xm) + + 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) + + @callback(walltime_interval=30.0) # Checkpoint the simulation every 30 seconds + def checkpoint(self, icheckpoint): + self.do_checkpoint(icheckpoint, purge_old_dumps=True, keep_first=1, keep_last=3) + + @callback(tsim_interval=0.02) # Do the analysis every 0.02 time units + def vtk_dump(self, idump): + model.do_vtk_dump(self.dump_prefix + f"dump_{idump:04}.vtk", True) + + @simulation_setup + def setup(self): + 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), 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.set_units(shamrock.UnitSystem()) + cfg.print_status() + model.set_solver_config(cfg) + model.init_scheduler(int(1e7), 1) + model.resize_simulation_box(self.bmin, self.bmax) + model.add_cube_fcc_3d(self.dr, self.bmin, self.bmax) + + vol_b = (self.xM - self.xm) * (self.yM - self.ym) * (self.zM - self.zm) + totmass = self.rho * vol_b + + pmass = model.total_mass_to_part_mass(totmass) + + model.set_value_in_a_box("uint", "f64", 1, self.bmin, self.bmax) + + model.set_field_value_lambda_f64_3("vxyz", self.vel_func) + model.set_particle_mass(pmass) + + model.set_cfl_cour(0.3) + model.set_cfl_force(0.25) + + +sim = Simulation(model) +sim.run() + +# %% +# 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() diff --git a/src/shammodels/sph/include/shammodels/sph/config/BCConfig.hpp b/src/shammodels/sph/include/shammodels/sph/config/BCConfig.hpp index 78dca7355f..94854e0220 100644 --- a/src/shammodels/sph/include/shammodels/sph/config/BCConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/config/BCConfig.hpp @@ -212,7 +212,7 @@ namespace shammodels::sph { p.set_shearing_periodic( j.at("shear_base").get(), j.at("shear_dir").get(), - j.at("speed").get()); + j.at("shear_speed").get()); } else { shambase::throw_unimplemented("wtf !"); }