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 0x7fec7c048380>)
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")
-> modified loglevel to 0 enabled log types :
log status :
- Loglevel: 1, enabled log types :
[xxx] Info: xxx ( logger::info )
[xxx] : xxx ( logger::normal )
[xxx] Warning: xxx ( logger::warn )
[xxx] Error: xxx ( logger::err )
Open the phantom dump
35 dump = shamrock.load_phantom_dump(filename)
36 dump.print_state()
--- dump state ---
table_header_fort_int len = 22
nparttot 576
ntypes 8
npartoftype 576
npartoftype 0
npartoftype 0
npartoftype 0
npartoftype 0
npartoftype 0
npartoftype 0
npartoftype 0
nblocks 1
nptmass 0
ndustlarge 0
ndustsmall 0
idust 7
idtmax_n 1
idtmax_frac 0
idumpfile 10
majorv 2023
minorv 0
microv 0
isink 0
table_header_i8 len = 0
table_header_i16 len = 0
table_header_i32 len = 2
iexternalforce 0
ieos 2
table_header_i64 len = 10
nparttot 576
ntypes 8
npartoftype 576
npartoftype 0
npartoftype 0
npartoftype 0
npartoftype 0
npartoftype 0
npartoftype 0
npartoftype 0
table_header_fort_real len = 37
gamma 1.6666666666666667
RK2 0
polyk2 0
qfacdisc 0.75
qfacdisc2 0.75
time 0.050000000000000225
dtmax 0.005
rhozero 1
hfact 1.2
tolh 0.0001
C_cour 0.3
C_force 0.25
alpha 0
alphau 1
alphaB 1
massoftype 0.001726334915006219
massoftype 0.001726334915006219
massoftype 0.001726334915006219
massoftype 0.001726334915006219
massoftype 0.001726334915006219
massoftype 0.001726334915006219
massoftype 0.001726334915006219
massoftype 0.001726334915006219
Bextx 0
Bexty 0
Bextz 0
dum 0
xmin -0.5
xmax 0.5
ymin -0.4871392896287467
ymax 0.4871392896287467
zmin -0.5103103630798287
zmax 0.5103103630798287
get_conserv -1
etot_in 1.0000000000000013
angtot_in 0
totmom_in 0
table_header_f32 len = 0
table_header_f64 len = 4
udist 1
umass 1
utime 1
umagfd 3.5449077018167907
block 0 :
--blocks_fort_int --
--blocks_i8 --
--blocks_i16 --
--blocks_i32 --
--blocks_i64 --
tag = iorig size = 576
--blocks_fort_real--
tag = x size = 576
tag = y size = 576
tag = z size = 576
tag = vx size = 576
tag = vy size = 576
tag = vz size = 576
tag = u size = 576
--blocks_f32 --
tag = h size = 576
tag = alpha size = 576
tag = divv size = 576
--blocks_f64 --
block 1 :
--blocks_fort_int --
--blocks_i8 --
--blocks_i16 --
--blocks_i32 --
--blocks_i64 --
--blocks_fort_real--
--blocks_f32 --
--blocks_f64 --
------------------
Start a SPH simulation from the phantom dump
40 ctx = shamrock.Context()
41 ctx.pdata_layout_new()
42 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
43
44 cfg = model.gen_config_from_phantom_dump(dump)
45 # Set the solver config to be the one stored in cfg
46 model.set_solver_config(cfg)
47 # Print the solver config
48 model.get_current_config().print_status()
49
50 model.init_scheduler(int(1e8), 1)
51
52 model.init_from_phantom_dump(dump)
adiabatic eos: gamma = 1.6666666666666667
Setting periodic boundaries from phantmdump
----- SPH Solver configuration -----
units :
unit_length : 1
unit_mass : 1
unit_current : 1
unit_temperature : 1
unit_qte : 1
unit_lumint : 1
part mass 0.001726334915006219 ( can be changed using .set_part_mass() )
cfl force 0.25
cfl courant 0.3
--- 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: 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 : 32.39 us (86.5%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 576 min = 576 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 576 min = 576 [LoadBalance][rank=0]
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 : 1042.00 ns (0.1%)
patch tree reduce : 661.00 ns (0.1%)
gen split merge : 591.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1122.00 ns (0.1%)
LB compute : 832.03 us (98.9%)
LB move op cnt : 0
LB apply : 2.01 us (0.2%)
Info: current particle counts : min = 576 max = 576 [Model][rank=0]
Run a simple timestep just for wasting some computing time :)
58 model.timestep()
---------------- t = 0.050000000000000225, dt = 0 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 576 min = 576 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 576 min = 576 [LoadBalance][rank=0]
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 : 15.69 us (1.9%)
patch tree reduce : 2.63 us (0.3%)
gen split merge : 831.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1142.00 ns (0.1%)
LB compute : 781.25 us (96.3%)
LB move op cnt : 0
LB apply : 2.35 us (0.3%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 3.14 us (75.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: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 9.128736121123468e-14, max = 9.545810285423902e-08
iterations = 0
Info: conservation infos : [sph::Model][rank=0]
sum v = (3.3129020355585587e-19,2.8749251562813254e-18,2.6832634791959036e-18)
sum a = (-3.222311945998652e-17,-4.354553303380781e-16,7.359808400080193e-17)
sum e = 0.9998775176108137
sum de = -3.608224830031759e-16
Info: cfl dt = 0.00010701881380824262 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 | 1.6798e+04 | 576 | 3.429e-02 | 20 % | 2 % | 1.10 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
<shamrock.model_sph.TimestepLog object at 0x7fec7756cf70>
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
Plot the result
68 import matplotlib.pyplot as plt
69
70 pixel_x = 1200
71 pixel_y = 1080
72 radius = 0.7
73 center = (0.0, 0.0, 0.0)
74
75 aspect = pixel_x / pixel_y
76 pic_range = [-radius * aspect, radius * aspect, -radius, radius]
77 delta_x = (radius * 2 * aspect, 0.0, 0.0)
78 delta_y = (0.0, radius * 2, 0.0)
79
80 arr_rho = model.render_cartesian_column_integ(
81 "rho", "f64", center=(0.0, 0.0, 0.0), delta_x=delta_x, delta_y=delta_y, nx=pixel_x, ny=pixel_y
82 )
83
84 import copy
85
86 import matplotlib
87
88 my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat")) # copy the default cmap
89 my_cmap.set_bad(color="black")
90
91 fig_width = 6
92 fig_height = fig_width / aspect
93 plt.figure(figsize=(fig_width, fig_height))
94 res = plt.imshow(arr_rho, cmap=my_cmap, origin="lower", extent=pic_range)
95
96 cbar = plt.colorbar(res, extend="both")
97 cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [code unit]")
98 # or r"$\rho$ [code unit]" for slices
99
100 plt.title("t = {:0.3f} [code unit]".format(model.get_time()))
101 plt.xlabel("x")
102 plt.ylabel("z")
103 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, center: (0,0,0), delta_x: (1.5555555555555556,0,0), delta_y: (0,1.4,0), nx: 1200, ny: 1080 [sph::CartesianRender][rank=0]
Info: compute_column_integ took 1889.58 ms [sph::CartesianRender][rank=0]
Total running time of the script: (0 minutes 3.190 seconds)
Estimated memory usage: 139 MB