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
![t = 0.000 [code unit]](../../_images/sphx_glr_run_sph_taylor_green_vortex_001.png)
Mach number : 0.10000000000000002
----- SPH Solver configuration -----
[
{
"artif_viscosity": {
"alpha_max": 1.0,
"alpha_min": 0.0,
"alpha_u": 1.0,
"av_type": "varying_cd10",
"beta_AV": 2.0,
"sigma_decay": 0.1
},
"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,
"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,
"self_grav_config": {
"softening_length": 1e-09,
"softening_mode": "plummer",
"type": "none"
},
"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 ) rate = 8.155878e+06 N.s^-1
SPH setup: the generation step took : 0.068307171 s
SPH setup: final particle count = 132224 begining 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.68 us (66.8%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 132224 min = 132224 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 132224 min = 132224 [LoadBalance][rank=0]
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 : 1834.00 ns (0.4%)
patch tree reduce : 1002.00 ns (0.2%)
gen split merge : 892.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 902.00 ns (0.2%)
LB compute : 451.43 us (95.4%)
LB move op cnt : 0
LB apply : 13.27 us (2.8%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.62 us (65.6%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 132224 min = 132224 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 132224 min = 132224 [LoadBalance][rank=0]
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 : 1343.00 ns (0.4%)
patch tree reduce : 401.00 ns (0.1%)
gen split merge : 391.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 291.00 ns (0.1%)
LB compute : 357.13 us (98.0%)
LB move op cnt : 0
LB apply : 1904.00 ns (0.5%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.52 us (63.8%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 132224 min = 132224 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 132224 min = 132224 [LoadBalance][rank=0]
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 : 1383.00 ns (0.4%)
patch tree reduce : 401.00 ns (0.1%)
gen split merge : 381.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 291.00 ns (0.1%)
LB compute : 337.20 us (97.9%)
LB move op cnt : 0
LB apply : 1844.00 ns (0.5%)
Info: --------------------------------------------- [DataInserterUtility][rank=0]
SPH setup: injected 132224 / 132224 => 100.0% | ranks with patchs = 1 / 1 <- global loop ->
SPH setup: the injection step took : 0.017392315000000002 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 | 1.1% 0.0% | 935.14 MB | 450.00 kB |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 0.12992277600000002 s
Total mass : 0.046875
Current part mass : 3.5451204017424977e-07
---------------- t = 0, dt = 0 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 132224 min = 132224 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 132224 min = 132224 [LoadBalance][rank=0]
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 : 7.00 us (1.3%)
patch tree reduce : 1803.00 ns (0.3%)
gen split merge : 692.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 901.00 ns (0.2%)
LB compute : 502.86 us (96.0%)
LB move op cnt : 0
LB apply : 4.05 us (0.8%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.90 us (70.1%)
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 hydro = 6.796398329179769e-07 sink sink = inf [SPH][rank=0]
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.2101e+04 | 132224 | 1 | 4.119e+00 | 0.0% | 0.2% 0.0% | 935.14 MB | 5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 1 [SPH][rank=0]
Info: time since start : 1521.450724727 (s) [SPH][rank=0]
Info: dump to _to_trash/dump_0000.vtk [VTK Dump][rank=0]
- took 19.21 ms, bandwidth = 551.33 MB/s
Info: compute_slice field_name: vxyz, positions count: 1166400 [sph::CartesianRender][rank=0]
Info: compute_slice took 2.32 s [sph::CartesianRender][rank=0]
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")
18
19 gamma = 1.4
20
21 rho_g = 1
22
23 Mach = 0.1
24 P_g = (Mach**-2) * rho_g / gamma
25
26 print(f"Mach number : {1 / np.sqrt(gamma * P_g / rho_g)}")
27
28 u_g = P_g / ((gamma - 1) * rho_g)
29
30 resol_per_green = 128
31 vortex_size = 1
32
33 L_green = vortex_size / (2 * np.pi)
34
35 ctx = shamrock.Context()
36 ctx.pdata_layout_new()
37
38 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6")
39
40 cfg = model.gen_default_config()
41 cfg.set_artif_viscosity_VaryingCD10(
42 alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
43 )
44 cfg.set_boundary_periodic()
45 cfg.set_eos_adiabatic(gamma)
46 # Set the CFL
47 cfg.set_cfl_cour(0.1)
48 cfg.set_cfl_force(0.1)
49
50 # Set the solver config to be the one stored in cfg
51 model.set_solver_config(cfg)
52
53 # Print the solver config
54 model.get_current_config().print_status()
55
56 # We want the patches to split above 10^8 part and merge if smaller than 1 part (e.g. disable patch)
57 model.init_scheduler(int(1e8), 1)
58
59
60 resol = resol_per_green
61 (xs, ys, zs) = model.get_box_dim_fcc_3d(1, resol, resol, resol)
62 dr = 1 / xs
63
64 print("Box size : ", xs, ys, zs)
65 (xs, ys, zs) = (1.0, 1.0, dr * 12)
66
67 model.resize_simulation_box((-xs / 2, -ys / 2, -zs / 2), (xs / 2, ys / 2, zs / 2))
68
69
70 setup = model.get_setup()
71 gen1 = setup.make_generator_lattice_hcp(dr, (-xs / 2, -ys / 2, -zs / 2), (xs / 2, ys / 2, zs / 2))
72 setup.apply_setup(gen1)
73
74 model.set_value_in_a_box("uint", "f64", u_g, (-xs / 2, -ys / 2, -zs / 2), (xs / 2, ys / 2, zs / 2))
75
76
77 def vel_func(r):
78 x, y, z = r
79
80 x += L_green * np.pi
81
82 vx = np.cos(x / L_green) * np.sin(y / L_green)
83 vy = -np.sin(x / L_green) * np.cos(y / L_green)
84 vz = 0
85
86 return (vx, vy, vz)
87
88
89 model.set_field_value_lambda_f64_3("vxyz", vel_func)
90
91
92 vol_b = xs * ys * zs
93
94 totmass = rho_g * vol_b
95
96 print("Total mass :", totmass)
97
98
99 pmass = model.total_mass_to_part_mass(totmass)
100 model.set_particle_mass(pmass)
101 print("Current part mass :", pmass)
102
103
104 dump_folder = "_to_trash"
105 import os
106
107 os.system("mkdir -p " + dump_folder)
108
109 current_fig = None
110
111 cnt_plot = 0
112
113
114 def plot():
115 global current_fig
116 if current_fig is not None:
117 plt.close(current_fig)
118
119 pixel_x = 1080
120 pixel_y = 1080
121 radius = 0.5
122 center = (0.0, 0.0, 0.0)
123 aspect = pixel_x / pixel_y
124 pic_range = [-radius * aspect, radius * aspect, -radius, radius]
125 delta_x = (radius * 2 * aspect, 0.0, 0.0)
126 delta_y = (0.0, radius * 2, 0.0)
127
128 arr_vel = model.render_cartesian_slice(
129 "vxyz",
130 "f64_3",
131 center=(0.0, 0.0, 0.0),
132 delta_x=delta_x,
133 delta_y=delta_y,
134 nx=pixel_x,
135 ny=pixel_y,
136 )
137
138 v_norm = np.sqrt(arr_vel[:, :, 0] ** 2 + arr_vel[:, :, 1] ** 2 + arr_vel[:, :, 2] ** 2)
139
140 import copy
141
142 import matplotlib
143
144 my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat")) # copy the default cmap
145 my_cmap.set_bad(color="black")
146
147 fig_width = 6
148 fig_height = fig_width / aspect
149 current_fig = plt.figure(figsize=(fig_width, fig_height))
150 res = plt.imshow(v_norm, cmap=my_cmap, origin="lower", extent=pic_range)
151
152 cbar = plt.colorbar(res, extend="both")
153 cbar.set_label(r"$\sqrt{vx^2 + vy^2 + vz^2}$ [code unit]")
154
155 plt.title("t = {:0.3f} [code unit]".format(model.get_time()))
156 plt.xlabel("x")
157 plt.ylabel("y")
158 global cnt_plot
159 plt.savefig(dump_folder + f"/test_{cnt_plot}.png")
160 cnt_plot += 1
161
162
163 model.timestep()
164
165 dt_stop = 0.001
166 for i in range(1):
167 t_target = i * dt_stop
168 # skip if the model is already past the target
169 if model.get_time() > t_target:
170 continue
171
172 model.evolve_until(i * dt_stop)
173
174 # Dump name is "dump_xxxx.sham" where xxxx is the timestep
175 model.do_vtk_dump(dump_folder + "/dump_{:04}.vtk".format(i), True)
176 plot()
Total running time of the script: (0 minutes 7.362 seconds)
Estimated memory usage: 161 MB