Note
Go to the end to download the full example code.
Uniform box in SPH#
This simple example shows a uniform density box in SPH, it is also used to test that the smoothing length iteration find the correct value
10 import shamrock
11
12 # If we use the shamrock executable to run this script instead of the python interpreter,
13 # we should not initialize the system as the shamrock executable needs to handle specific MPI logic
14 if not shamrock.sys.is_initialized():
15 shamrock.change_loglevel(1)
16 shamrock.sys.init("0:0")
Setup parameters
22 gamma = 5.0 / 3.0
23 rho_g = 1
24
25 bmin = (-0.6, -0.6, -0.6)
26 bmax = (0.6, 0.6, 0.6)
27
28 N_target = 1e4
29 scheduler_split_val = int(2e7)
30 scheduler_merge_val = int(1)
Deduced quantities
34 import numpy as np
35
36 xm, ym, zm = bmin
37 xM, yM, zM = bmax
38 vol_b = (xM - xm) * (yM - ym) * (zM - zm)
39
40 part_vol = vol_b / N_target
41
42 # lattice volume
43 HCP_PACKING_DENSITY = 0.74
44 part_vol_lattice = HCP_PACKING_DENSITY * part_vol
45
46 dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0)
47
48 pmass = -1
Setup
53 ctx = shamrock.Context()
54 ctx.pdata_layout_new()
55
56 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
57
58 cfg = model.gen_default_config()
59 cfg.set_artif_viscosity_VaryingCD10(
60 alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
61 )
62 cfg.set_boundary_periodic()
63 cfg.set_eos_adiabatic(gamma)
64 cfg.print_status()
65 model.set_solver_config(cfg)
66 model.init_scheduler(scheduler_split_val, scheduler_merge_val)
67
68
69 bmin, bmax = model.get_ideal_hcp_box(dr, bmin, bmax)
70 xm, ym, zm = bmin
71 xM, yM, zM = bmax
72
73 model.resize_simulation_box(bmin, bmax)
74
75 setup = model.get_setup()
76 gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
77 setup.apply_setup(gen, insert_step=scheduler_split_val)
78
79
80 xc, yc, zc = model.get_closest_part_to((0, 0, 0))
81
82 if shamrock.sys.world_rank() == 0:
83 print("closest part to (0,0,0) is in :", xc, yc, zc)
84
85
86 vol_b = (xM - xm) * (yM - ym) * (zM - zm)
87
88 totmass = rho_g * vol_b
89
90 pmass = model.total_mass_to_part_mass(totmass)
91
92 model.set_value_in_a_box("uint", "f64", 0, bmin, bmax)
93
94 tot_u = pmass * model.get_sum("uint", "f64")
95 if shamrock.sys.world_rank() == 0:
96 print("total u :", tot_u)
97
98 model.set_particle_mass(pmass)
99
100 model.set_cfl_cour(0.1)
101 model.set_cfl_force(0.1)
----- SPH Solver configuration -----
units :
not set
part mass 0 ( can be changed using .set_part_mass() )
cfl force 0
cfl courant 0
--- 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.6666666666666667
--- Bondaries config
Config Type : Periodic boundaries
--- Bondaries config config (deduced)
-------------
------------------------------------
Info: pushing data in scheduler, N = 11520 [DataInserterUtility][rank=0]
Info: reattributing data ... [DataInserterUtility][rank=0]
Info: reattributing data done in 4.41 ms [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 : 3.52 us (55.7%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 11520 min = 11520 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 11520 min = 11520 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 11520
max = 11520
avg = 11520
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1663.00 ns (0.4%)
patch tree reduce : 581.00 ns (0.1%)
gen split merge : 621.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 871.00 ns (0.2%)
LB compute : 381.08 us (97.4%)
LB move op cnt : 0
LB apply : 2.87 us (0.7%)
Info: the setup took : 0.013129176000000001 s [SPH setup][rank=0]
closest part to (0,0,0) is in : 0.0 0.0 0.0
total u : 0.0
Single timestep to iterate the smoothing length
105 model.timestep()
---------------- t = 0, dt = 0 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 11520 min = 11520 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 11520 min = 11520 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 11520
max = 11520
avg = 11520
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 19.88 us (2.2%)
patch tree reduce : 5.24 us (0.6%)
gen split merge : 631.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 972.00 ns (0.1%)
LB compute : 870.23 us (95.4%)
LB move op cnt : 0
LB apply : 3.81 us (0.4%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 3.83 us (78.0%)
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.39592013888888883
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.03437860905706778 unconverged cnt = 11520
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.4238715277777777
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.03781646996277456 unconverged cnt = 11520
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.4238715277777777
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.041598116959052016 unconverged cnt = 11520
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.4506076388888888
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.045757928654957224 unconverged cnt = 11520
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.5162326388888888
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.05033372152045295 unconverged cnt = 11520
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.6863715277777777
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.05536709367249825 unconverged cnt = 11520
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.6863715277777777
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.06090380303974808 unconverged cnt = 11520
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.7886284722222223
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 2.831433839859919e-07, max = 2.831433851253122e-07
iterations = 3
Info: conservation infos : [sph::Model][rank=0]
sum v = (0,0,0)
sum a = (0,0,0)
sum e = 0
sum de = 0
Info: cfl dt = 1.7976931348623157e+308 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) | MPI | alloc | mem (max) |
---------------------------------------------------------------------------------------
| 0 | 3.7830e+04 | 11520 | 3.045e-01 | 0 % | 1 % | 11.43 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
<shamrock.model_sph.TimestepLog object at 0x7f0ca497c730>
Recover data
109 dat = ctx.collect_data()
Info: collected : 1 patches [PatchScheduler][rank=0]
adding -> xyz
adding -> vxyz
adding -> axyz
adding -> axyz_ext
adding -> hpart
adding -> uint
adding -> duint
adding -> alpha_AV
adding -> divv
adding -> dtdivv
adding -> curlv
adding -> soundspeed
Test h value
113 import numpy as np
114
115 min_hpart = np.min(dat["hpart"])
116 max_hpart = np.max(dat["hpart"])
117 mean_hpart = np.mean(dat["hpart"])
118
119 print(f"hpart min={min_hpart} max={max_hpart} delta={max_hpart-min_hpart}")
120
121 assert np.abs(max_hpart - min_hpart) < 1e-15, "hpart delta is too large"
122
123 expected_h = (0.06688949833400996 + 0.06688949833401027) / 2
124
125 assert np.abs(min_hpart - expected_h) < 1e-15, "hpart is off the expected value"
hpart min=0.06688949833400996 max=0.06688949833401027 delta=3.0531133177191805e-16
Plot particle distrib
129 import matplotlib
130 import matplotlib.pyplot as plt
131 from mpl_toolkits.mplot3d import Axes3D
132
133 fig = plt.figure(dpi=120)
134 ax = fig.add_subplot(111, projection="3d")
135 ax.set_xlim3d(bmin[0], bmax[0])
136 ax.set_ylim3d(bmin[1], bmax[1])
137 ax.set_zlim3d(bmin[2], bmax[2])
138
139 cm = matplotlib.colormaps["viridis"]
140 sc = ax.scatter(
141 dat["xyz"][:, 0],
142 dat["xyz"][:, 1],
143 dat["xyz"][:, 2],
144 s=1,
145 vmin=mean_hpart - 1e-10,
146 vmax=mean_hpart + 1e-10,
147 c=dat["hpart"],
148 cmap=cm,
149 )
150 plt.colorbar(sc)
151 plt.show()
Total running time of the script: (0 minutes 1.462 seconds)
Estimated memory usage: 71 MB