[SPH] sheared wave setup

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()

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<i32_3>(),
                 j.at("shear_dir").get<i32_3>(),
-                j.at("speed").get<Tscal>());
+                j.at("shear_speed").get<Tscal>());
         } else {
             shambase::throw_unimplemented("wtf !");
         }