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")
Use shamrock documentation style for matplotlib
21 shamrock.matplotlib.set_shamrock_mpl_style()
Setup parameters
27 gamma = 5.0 / 3.0
28 rho_g = 1
29
30 bmin = (-0.6, -0.6, -0.6)
31 bmax = (0.6, 0.6, 0.6)
32
33 N_target = 1e4
34 scheduler_split_val = int(2e7)
35 scheduler_merge_val = int(1)
Deduced quantities
39 import numpy as np
40
41 xm, ym, zm = bmin
42 xM, yM, zM = bmax
43 vol_b = (xM - xm) * (yM - ym) * (zM - zm)
44
45 part_vol = vol_b / N_target
46
47 # lattice volume
48 HCP_PACKING_DENSITY = 0.74
49 part_vol_lattice = HCP_PACKING_DENSITY * part_vol
50
51 dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0)
52
53
54 bmin, bmax = shamrock.math.get_ideal_hcp_box(dr, bmin, bmax)
55 xm, ym, zm = bmin
56 xM, yM, zM = bmax
57 vol_b = (xM - xm) * (yM - ym) * (zM - zm)
58
59 totmass = rho_g * vol_b
60
61 pmass = -1
Setup
66 ctx = shamrock.Context()
67 ctx.pdata_layout_new()
68
69 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
70
71 cfg = model.gen_default_config()
72 cfg.set_artif_viscosity_VaryingCD10(
73 alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
74 )
75 cfg.set_boundary_periodic()
76 cfg.set_eos_adiabatic(gamma)
77 cfg.print_status()
78 model.set_solver_config(cfg)
79 model.init_scheduler(scheduler_split_val, scheduler_merge_val)
80
81 model.resize_simulation_box(bmin, bmax)
82
83 setup = model.get_setup()
84 gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
85 setup.apply_setup(gen, insert_step=scheduler_split_val)
86
87
88 xc, yc, zc = model.get_closest_part_to((0, 0, 0))
89
90 if shamrock.sys.world_rank() == 0:
91 print("closest part to (0,0,0) is in :", xc, yc, zc)
92
93
94 pmass = model.total_mass_to_part_mass(totmass)
95
96 model.set_value_in_a_box("uint", "f64", 0, bmin, bmax)
97
98 tot_u = pmass * model.get_sum("uint", "f64")
99 if shamrock.sys.world_rank() == 0:
100 print("total u :", tot_u)
101
102 model.set_particle_mass(pmass)
103
104 model.set_cfl_cour(0.1)
105 model.set_cfl_force(0.1)
----- 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.0,
"cfl_force": 0.0,
"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.6666666666666667
},
"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": "M4<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
}
]
------------------------------------
SPH setup: generating particles ...
SPH setup: Nstep = 11520 ( 1.2e+04 ) Ntotal = 11520 ( 1.2e+04 rank min = 4.9e+06 max = 1.2e+04) rate = 1.152000e+04 N.s^-1
SPH setup: the generation step took : 0.006088026000000001 s
SPH setup: final particle count = 11520 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 : 5.21 us (64.0%)
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 11520.0 min = 11520.0 factor = 1
- strategy "round robin" : max = 10944.0 min = 10944.0 factor = 0.95
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 : 951.00 ns (0.3%)
patch tree reduce : 1.44 us (0.4%)
gen split merge : 1.08 us (0.3%)
split / merge op : 0/0
apply split merge : 842.00 ns (0.2%)
LB compute : 359.30 us (97.2%)
LB move op cnt : 0
LB apply : 3.46 us (0.9%)
Info: patch count stable after 1 runs npatch = 1 [DataInserterUtility][rank=0]
Info: --------------------------------------------- [DataInserterUtility][rank=0]
SPH setup: injected 11520 / 11520 => 100.0% | ranks with patchs = 1 / 1 <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.006079323 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.9% 0.0% | 134.92 MB | 0.00 B |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 0.016436429000000002 s
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
109 model.timestep()
---------------- t = 0, dt = 0 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 11520.0 min = 11520.0 factor = 1
- strategy "round robin" : max = 10944.0 min = 10944.0 factor = 0.95
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 : 23.03 us (5.5%)
patch tree reduce : 2.48 us (0.6%)
gen split merge : 591.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 931.00 ns (0.2%)
LB compute : 369.05 us (88.4%)
LB move op cnt : 0
LB apply : 4.46 us (1.1%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.99 us (74.3%)
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
Warning: the unit system is not set [sph::Config][rank=0]
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: 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.0019e+04 | 11520 | 1 | 3.838e-01 | 0.0% | 0.8% 0.0% | 134.92 MB | 460.80 kB |
+------+------------+-------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
<shamrock.model_sph.TimestepLog object at 0x7fbc60e5fbb0>
Recover data
113 dat = ctx.collect_data()
Info: collected : 1 patches [PatchScheduler][rank=0]
Test h value
117 import numpy as np
118
119 min_hpart = np.min(dat["hpart"])
120 max_hpart = np.max(dat["hpart"])
121 mean_hpart = np.mean(dat["hpart"])
122
123 print(f"hpart min={min_hpart} max={max_hpart} delta={max_hpart - min_hpart}")
124
125 assert np.abs(max_hpart - min_hpart) < 1e-15, "hpart delta is too large"
126
127 expected_h = (0.06688949833400996 + 0.06688949833401027) / 2
128
129 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
133 import matplotlib
134 import matplotlib.pyplot as plt
135 from mpl_toolkits.mplot3d import Axes3D
136
137 fig = plt.figure(dpi=120)
138 ax = fig.add_subplot(111, projection="3d")
139 ax.set_xlim3d(bmin[0], bmax[0])
140 ax.set_ylim3d(bmin[1], bmax[1])
141 ax.set_zlim3d(bmin[2], bmax[2])
142
143 cm = matplotlib.colormaps["viridis"]
144 sc = ax.scatter(
145 dat["xyz"][:, 0],
146 dat["xyz"][:, 1],
147 dat["xyz"][:, 2],
148 s=1,
149 vmin=mean_hpart - 1e-10,
150 vmax=mean_hpart + 1e-10,
151 c=dat["hpart"],
152 cmap=cm,
153 )
154
155 ax.minorticks_off()
156
157 plt.colorbar(sc)
158 plt.show()
Total running time of the script: (0 minutes 2.257 seconds)
Estimated memory usage: 159 MB