Note
Go to the end to download the full example code.
Taylor green vortex in SPH#
This simple example shows a SPH simulation of a Taylor green vortex
8 import matplotlib.pyplot as plt
9 import numpy as np
10
11 import shamrock
12
13 # If we use the shamrock executable to run this script instead of the python interpreter,
14 # we should not initialize the system as the shamrock executable needs to handle specific MPI logic
15 if not shamrock.sys.is_initialized():
16 shamrock.change_loglevel(1)
17 shamrock.sys.init("0:0")
Use shamrock documentation style for matplotlib
21 shamrock.matplotlib.set_shamrock_mpl_style()
Setup parameters
Mach number : 0.10000000000000002
Setup the context, model & sim
43 ctx = shamrock.Context()
44 ctx.pdata_layout_new()
45
46 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6")
47
48 cfg = model.gen_default_config()
49 cfg.set_artif_viscosity_VaryingCD10(
50 alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
51 )
52 cfg.set_boundary_periodic()
53 cfg.set_eos_adiabatic(gamma)
54 # Set the CFL
55 cfg.set_cfl_cour(0.1)
56 cfg.set_cfl_force(0.1)
57
58 # Set the solver config to be the one stored in cfg
59 model.set_solver_config(cfg)
60
61 # Print the solver config
62 model.get_current_config().print_status()
63
64 # We want the patches to split above 10^8 part and merge if smaller than 1 part (e.g. disable patch)
65 model.init_scheduler(int(1e8), 1)
66
67
68 resol = resol_per_green
69 (xs, ys, zs) = model.get_box_dim_fcc_3d(1, resol, resol, resol)
70 dr = 1 / xs
71
72 print("Box size : ", xs, ys, zs)
73 (xs, ys, zs) = (1.0, 1.0, dr * 12)
74
75 model.resize_simulation_box((-xs / 2, -ys / 2, -zs / 2), (xs / 2, ys / 2, zs / 2))
76
77
78 setup = model.get_setup()
79 gen1 = setup.make_generator_lattice_hcp(dr, (-xs / 2, -ys / 2, -zs / 2), (xs / 2, ys / 2, zs / 2))
80 setup.apply_setup(gen1)
81
82 model.set_value_in_a_box("uint", "f64", u_g, (-xs / 2, -ys / 2, -zs / 2), (xs / 2, ys / 2, zs / 2))
83
84
85 def vel_func(r):
86 x, y, z = r
87
88 x += L_green * np.pi
89
90 vx = np.cos(x / L_green) * np.sin(y / L_green)
91 vy = -np.sin(x / L_green) * np.cos(y / L_green)
92 vz = 0
93
94 return (vx, vy, vz)
95
96
97 model.set_field_value_lambda_f64_3("vxyz", vel_func)
98
99
100 vol_b = xs * ys * zs
101
102 totmass = rho_g * vol_b
103
104 print("Total mass :", totmass)
105
106
107 pmass = model.total_mass_to_part_mass(totmass)
108 model.set_particle_mass(pmass)
109 print("Current part mass :", pmass)
110
111
112 dump_folder = "_to_trash"
113 import os
114
115 os.system("mkdir -p " + dump_folder)
116
117 current_fig = None
118
119 cnt_plot = 0
----- SPH Solver configuration -----
[
{
"artif_viscosity": {
"alpha_max": 1.0,
"alpha_min": 0.0,
"alpha_u": 1.0,
"beta_AV": 2.0,
"sigma_decay": 0.1,
"type": "varying_cd10"
},
"boundary_config": {
"bc_type": "periodic"
},
"cfl_config": {
"cfl_cour": 0.1,
"cfl_force": 0.1,
"cfl_multiplier_stiffness": 2.0,
"eta_sink": 0.05
},
"combined_dtdiv_divcurlv_compute": false,
"debug_dump_filename": "",
"do_debug_dump": false,
"dust_config": {
"drag_mode": {
"type": "none"
},
"mode": {
"type": "none"
}
},
"enable_particle_reordering": false,
"eos_config": {
"Tvec": "f64_3",
"eos_type": "adiabatic",
"gamma": 1.4
},
"epsilon_h": 1e-06,
"ext_force_config": {
"force_list": []
},
"gpart_mass": 0.0,
"h_iter_per_subcycles": 50,
"h_max_subcycles_count": 100,
"htol_up_coarse_cycle": 1.1,
"htol_up_fine_cycle": 1.1,
"kernel_id": "M6<f64>",
"mhd_config": {
"mhd_type": "none"
},
"particle_killing": [],
"particle_reordering_step_freq": 1000,
"save_dt_to_fields": false,
"scheduler_config": {
"merge_load_value": 0,
"split_load_value": 0
},
"self_grav_config": {
"softening_length": 1e-09,
"softening_mode": "plummer",
"type": "none"
},
"show_cfl_detail": false,
"show_ghost_zone_graph": false,
"show_neigh_stats": false,
"smoothing_length_config": {
"type": "density_based"
},
"time_state": {
"cfl_multiplier": 0.01,
"dt_sph": 0.0,
"time": 0.0
},
"tree_reduction_level": 3,
"type_id": "sycl::vec<f64,3>",
"unit_sys": null,
"use_two_stage_search": true
}
]
------------------------------------
Box size : 256.0 221.70250336881628 209.02312471749784
SPH setup: generating particles ...
SPH setup: Nstep = 132224 ( 1.3e+05 ) Ntotal = 132224 ( 1.3e+05 rank min = 7.5e+06 max = 1.3e+05) rate = 1.322240e+05 N.s^-1
SPH setup: the generation step took : 0.027431224 s
SPH setup: final particle count = 132224 beginning injection ...
Info: --------------------------------------------- [DataInserterUtility][rank=0]
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 8.17 us (64.2%)
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 132224.0 min = 132224.0 factor = 1
- strategy "round robin" : max = 125612.8 min = 125612.8 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 132224
max = 132224
avg = 132224
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 912.00 ns (0.2%)
patch tree reduce : 1.36 us (0.3%)
gen split merge : 851.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 992.00 ns (0.2%)
LB compute : 411.90 us (97.5%)
LB move op cnt : 0
LB apply : 3.54 us (0.8%)
Info: patch count stable after 1 runs npatch = 1 [DataInserterUtility][rank=0]
Info: --------------------------------------------- [DataInserterUtility][rank=0]
SPH setup: injected 132224 / 132224 => 100.0% | ranks with patchs = 1 / 1 <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.018049718000000003 s
Info: injection perf report: [SPH setup][rank=0]
+======+====================+=======+=============+=============+=============+
| rank | rank get (sum/max) | MPI | alloc d% h% | mem (max) d | mem (max) h |
+======+====================+=======+=============+=============+=============+
| 0 | 0.00s / 0.00s | 0.00s | 0.8% 0.0% | 134.92 MB | 460.80 kB |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 0.07498199500000001 s
Total mass : 0.046875
Current part mass : 3.5451204017424977e-07
Make plotting functions
125 def plot():
126 global current_fig
127 if current_fig is not None:
128 plt.close(current_fig)
129
130 pixel_x = 1080
131 pixel_y = 1080
132 radius = 0.5
133 center = (0.0, 0.0, 0.0)
134 aspect = pixel_x / pixel_y
135 pic_range = [-radius * aspect, radius * aspect, -radius, radius]
136 delta_x = (radius * 2 * aspect, 0.0, 0.0)
137 delta_y = (0.0, radius * 2, 0.0)
138
139 arr_vel = model.render_cartesian_slice(
140 "vxyz",
141 "f64_3",
142 center=(0.0, 0.0, 0.0),
143 delta_x=delta_x,
144 delta_y=delta_y,
145 nx=pixel_x,
146 ny=pixel_y,
147 )
148
149 v_norm = np.sqrt(arr_vel[:, :, 0] ** 2 + arr_vel[:, :, 1] ** 2 + arr_vel[:, :, 2] ** 2)
150
151 import copy
152
153 import matplotlib
154
155 my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat")) # copy the default cmap
156 my_cmap.set_bad(color="black")
157
158 fig_width = 6
159 fig_height = fig_width / aspect
160 current_fig = plt.figure(figsize=(fig_width, fig_height))
161 res = plt.imshow(v_norm, cmap=my_cmap, origin="lower", extent=pic_range)
162
163 cbar = plt.colorbar(res, extend="both")
164 cbar.set_label(r"$\sqrt{vx^2 + vy^2 + vz^2}$ [code unit]")
165
166 plt.title("t = {:0.3f} [code unit]".format(model.get_time()))
167 plt.xlabel("x")
168 plt.ylabel("y")
169 global cnt_plot
170 plt.savefig(dump_folder + f"/test_{cnt_plot}.png")
171 cnt_plot += 1
Run the simulation
176 model.timestep()
177
178 dt_stop = 0.001
179 for i in range(1):
180 t_target = i * dt_stop
181 # skip if the model is already past the target
182 if model.get_time() > t_target:
183 continue
184
185 model.evolve_until(i * dt_stop)
186
187 # Dump name is "dump_xxxx.sham" where xxxx is the timestep
188 model.do_vtk_dump(dump_folder + "/dump_{:04}.vtk".format(i), True)
189 plot()
![t = 0.000 [code unit]](../../_images/sphx_glr_run_sph_taylor_green_vortex_001.png)
---------------- t = 0, dt = 0 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 132224.0 min = 132224.0 factor = 1
- strategy "round robin" : max = 125612.8 min = 125612.8 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 132224
max = 132224
avg = 132224
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 21.85 us (4.3%)
patch tree reduce : 3.06 us (0.6%)
gen split merge : 580.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1.23 us (0.2%)
LB compute : 460.15 us (90.6%)
LB move op cnt : 0
LB apply : 3.94 us (0.8%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.81 us (73.9%)
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 0.6578911544046468
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.004296875 unconverged cnt = 132224
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 0.6578911544046468
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.004726562500000001 unconverged cnt = 132224
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 0.6638809898354309
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.005199218750000002 unconverged cnt = 132224
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 0.9830817400774442
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.0057191406250000024 unconverged cnt = 132224
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 0.9911286907066795
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.006291054687500003 unconverged cnt = 132224
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 1.0163283518877058
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.006920160156250004 unconverged cnt = 132224
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 1.0243753025169406
Warning: the unit system is not set [sph::Config][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.9647884481053865e-19,4.773428345076549e-19,0)
sum a = (1.367610070459655e-16,-9.401966136903069e-16,-1.780827824668828e-14)
sum e = 8.382254464285767
sum de = 3.6917083973131426e-16
Info: cfl dt = 6.796398329179769e-07 cfl multiplier : 0.01 [sph::Model][rank=0]
Info: Timestep perf report: [sph::Model][rank=0]
+======+============+========+========+===========+======+=============+=============+=============+
| rank | rate (N/s) | Nobj | Npatch | tstep | MPI | alloc d% h% | mem (max) d | mem (max) h |
+======+============+========+========+===========+======+=============+=============+=============+
| 0 | 3.8588e+04 | 132224 | 1 | 3.427e+00 | 0.0% | 0.2% 0.0% | 188.57 MB | 5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
Info: evolve_until (target_time = 0.00s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
Info: iteration since start : 1 [SPH][rank=0]
Info: time since start : 424.41294305900004 (s) [SPH][rank=0]
Info: dump to _to_trash/dump_0000.vtk [VTK Dump][rank=0]
- took 8.13 ms, bandwidth = 1.37 GB/s
Info: compute_slice field_name: vxyz, positions count: 1166400 [sph::CartesianRender][rank=0]
Info: compute_slice took 1.95 s [sph::CartesianRender][rank=0]
Total running time of the script: (0 minutes 6.899 seconds)
Estimated memory usage: 345 MB