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 -----
units :
not set
part mass 0 ( can be changed using .set_part_mass() )
cfl force 0.1
cfl courant 0.1
--- artificial viscosity config
Config Type : VaryingCD10 (Cullen & Dehnen 2010)
alpha_min = 0
alpha_max = 1
sigma_decay = 0.1
alpha_u = 1
beta_AV = 2
--- artificial viscosity config (deduced)
-------------
EOS config f64_3 :
adiabatic :
gamma 1.4
--- Bondaries config
Config Type : Periodic boundaries
--- Bondaries config config (deduced)
-------------
------------------------------------
Box size : 256.0 221.70250336881628 209.02312471749784
Info: pushing data in scheduler, N = 132224 [DataInserterUtility][rank=0]
Info: reattributing data ... [DataInserterUtility][rank=0]
Info: reattributing data done in 23.55 ms [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 3.61 us (39.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 : 1052.00 ns (0.2%)
patch tree reduce : 1052.00 ns (0.2%)
gen split merge : 1052.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 912.00 ns (0.2%)
LB compute : 446.74 us (97.9%)
LB move op cnt : 0
LB apply : 1883.00 ns (0.4%)
Info: the setup took : 0.07165858400000001 s [SPH setup][rank=0]
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.38 us (1.3%)
patch tree reduce : 1403.00 ns (0.2%)
gen split merge : 1032.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1433.00 ns (0.3%)
LB compute : 548.78 us (96.4%)
LB move op cnt : 0
LB apply : 1974.00 ns (0.3%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.06 us (68.7%)
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
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 1.9302601001271766e-10, max = 9.282277784123224e-07
iterations = 3
Info: conservation infos : [sph::Model][rank=0]
sum v = (-6.045502917247344e-20,1.763271684197142e-20,0)
sum a = (1.0513199723902636e-15,-6.716150707533982e-16,-1.2009754275286882e-14)
sum e = 8.382254464285712
sum de = 7.605678239985814e-17
Info: cfl dt = 6.796398329179769e-07 cfl multiplier : 0.01 [sph::Model][rank=0]
Info: processing rate infos : [sph::Model][rank=0]
---------------------------------------------------------------------------------------
| rank | rate (N.s^-1) | Nobj | t compute (s) | interf | alloc | mem (max) |
---------------------------------------------------------------------------------------
| 0 | 3.5718e+04 | 132224 | 3.702e+00 | 14 % | 0 % | 160.18 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 1 [SPH][rank=0]
Info: time since start : 14.692532593000001 (s) [SPH][rank=0]
Info: dump to _to_trash/dump_0000.vtk [VTK Dump][rank=0]
- took 15.20 ms, bandwidth = 696.95 MB/s
Info: compute_slice field_name: vxyz, center: (0,0,0), delta_x: (1,0,0), delta_y: (0,1,0), nx: 1080, ny: 1080 [sph::CartesianRender][rank=0]
Info: compute_slice took 2.33 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
168 t_target = i * dt_stop
169 # skip if the model is already past the target
170 if model.get_time() > t_target:
171 continue
172
173 model.evolve_until(i * dt_stop)
174
175 # Dump name is "dump_xxxx.sham" where xxxx is the timestep
176 model.do_vtk_dump(dump_folder + "/dump_{:04}.vtk".format(i), True)
177 plot()
Total running time of the script: (0 minutes 7.081 seconds)
Estimated memory usage: 253 MB