Note
Go to the end to download the full example code.
Start a SPH simulation from a phantom dump#
This simple example shows how to start a SPH simulation from a phantom dump
Download a phantom dump
10 dump_folder = "_to_trash"
11 import os
12
13 os.system("mkdir -p " + dump_folder)
14
15 url = "https://raw.githubusercontent.com/Shamrock-code/reference-files/refs/heads/main/blast_00010"
16
17 filename = dump_folder + "/blast_00010"
18
19 from urllib.request import urlretrieve
20
21 urlretrieve(url, filename)
('_to_trash/blast_00010', <http.client.HTTPMessage object at 0x7fbc84236e00>)
Init shamrock
25 import shamrock
26
27 # If we use the shamrock executable to run this script instead of the python interpreter,
28 # we should not initialize the system as the shamrock executable needs to handle specific MPI logic
29 if not shamrock.sys.is_initialized():
30 shamrock.change_loglevel(1)
31 shamrock.sys.init("0:0")
32
33 ctx, model, in_params = shamrock.utils.phantom.load_simulation(
34 dump_folder, dump_file_name="blast_00010", in_file_name=None
35 )
-----------------------------------------------------------
---------------- Phantom dump loading -----------------
-----------------------------------------------------------
- Loading phantom dump from: _to_trash/blast_00010
- Generating Shamrock solver config from phantom dump
adiabatic eos: gamma = 1.6666666666666667
Setting periodic boundaries from phantmdump
- Setting Shamrock solver config
- Initializing domain scheduler
- Initializing from phantom dump (setup file: False)
Info: Push particles : [Model][rank=0]
rank = 0 patch id=0, add N=576 particles, coords = (-0.5,-0.4871392896287467,-0.5103103630798287) (0.5,0.4871392896287467,0.5103103630798287)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 11.51 us (79.1%)
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 576.0 min = 576.0 factor = 1
- strategy "round robin" : max = 547.2 min = 547.2 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 576
max = 576
avg = 576
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.81 us (0.7%)
patch tree reduce : 831.00 ns (0.2%)
gen split merge : 962.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 992.00 ns (0.2%)
LB compute : 397.27 us (96.6%)
LB move op cnt : 0
LB apply : 4.17 us (1.0%)
Info: current particle counts : min = 576 max = 576 [Model][rank=0]
- Shamrock solver config:
----- 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.3,
"cfl_force": 0.25,
"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.001726334915006219,
"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": 1,
"split_load_value": 100000000
},
"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.050000000000000225
},
"tree_reduction_level": 3,
"type_id": "sycl::vec<f64,3>",
"unit_sys": {
"unit_current": 1.0,
"unit_length": 1.0,
"unit_lumint": 1.0,
"unit_mass": 1.0,
"unit_qte": 1.0,
"unit_temperature": 1.0,
"unit_time": 1.0
},
"use_two_stage_search": true
}
]
------------------------------------
Run a simple timestep just for wasting some computing time :)
42 model.timestep()
---------------- t = 0.050000000000000225, dt = 0 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 576.0 min = 576.0 factor = 1
- strategy "round robin" : max = 547.2 min = 547.2 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 576
max = 576
avg = 576
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 19.15 us (4.7%)
patch tree reduce : 2.38 us (0.6%)
gen split merge : 550.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1.02 us (0.3%)
LB compute : 355.50 us (88.1%)
LB move op cnt : 0
LB apply : 9.99 us (2.5%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.22 us (71.6%)
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: 4.104166666666667
Info: conservation infos : [sph::Model][rank=0]
sum v = (3.1631663503355727e-19,5.832204939435293e-19,1.2313420489508965e-17)
sum a = (-6.33202265670962e-17,-4.423671295679711e-16,5.716848567539504e-16)
sum e = 0.9998775176108142
sum de = 7.582029080098474e-16
Info: cfl dt = 0.00010701881380824262 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 | 1.5609e+04 | 576 | 1 | 3.690e-02 | 0.1% | 2.8% 0.0% | 110.00 MB | 0.00 B |
+------+------------+------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
<shamrock.model_sph.TimestepLog object at 0x7fbc85e0f8b0>
Note
Note that shamrock has to update some smoothing lengths that were in the phantom dump I think that since smoothing length is single precision in phantom they are slightly off from the shamrock point of view hence the update
Use shamrock documentation style for matplotlib
52 shamrock.matplotlib.set_shamrock_mpl_style()
Plot the result
57 import matplotlib.pyplot as plt
58
59 pixel_x = 1200
60 pixel_y = 1080
61 radius = 0.7
62 center = (0.0, 0.0, 0.0)
63
64 aspect = pixel_x / pixel_y
65 pic_range = [-radius * aspect, radius * aspect, -radius, radius]
66 delta_x = (radius * 2 * aspect, 0.0, 0.0)
67 delta_y = (0.0, radius * 2, 0.0)
68
69 arr_rho = model.render_cartesian_column_integ(
70 "rho", "f64", center=(0.0, 0.0, 0.0), delta_x=delta_x, delta_y=delta_y, nx=pixel_x, ny=pixel_y
71 )
72
73 import copy
74
75 import matplotlib
76
77 my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat")) # copy the default cmap
78 my_cmap.set_bad(color="black")
79
80 fig_width = 6
81 fig_height = fig_width / aspect
82 plt.figure(figsize=(fig_width, fig_height))
83 res = plt.imshow(arr_rho, cmap=my_cmap, origin="lower", extent=pic_range)
84
85 cbar = plt.colorbar(res, extend="both")
86 cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [code unit]")
87 # or r"$\rho$ [code unit]" for slices
88
89 plt.title("t = {:0.3f} [code unit]".format(model.get_time()))
90 plt.xlabel("x")
91 plt.ylabel("z")
92 plt.show()
![t = 0.050 [code unit]](../../_images/sphx_glr_run_start_sph_from_phantom_dump_001.png)
Info: compute_column_integ field_name: rho, rays count: 1296000 [sph::CartesianRender][rank=0]
Info: compute_column_integ took 1.82 s [sph::CartesianRender][rank=0]
Total running time of the script: (0 minutes 3.178 seconds)
Estimated memory usage: 283 MB