Note
Go to the end to download the full example code.
Custom warp disc simulation#
This simple example shows how to run a simulation of a disc with a custom initial warp
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")
Use shamrock documentation style for matplotlib
22 shamrock.matplotlib.set_shamrock_mpl_style()
Setup units
28 si = shamrock.UnitSystem()
29 sicte = shamrock.Constants(si)
30 codeu = shamrock.UnitSystem(
31 unit_time=3600 * 24 * 365,
32 unit_length=sicte.au(),
33 unit_mass=sicte.sol_mass(),
34 )
35 ucte = shamrock.Constants(codeu)
36 G = ucte.G()
List parameters
42 # Resolution
43 Npart = 1000000
44
45 # Sink parameters
46 center_mass = 1.0
47 center_racc = 0.1
48
49 # Disc parameter
50 disc_mass = 0.01 # sol mass
51 rout = 60.0 # au
52 rin = 1.0 # au
53 H_r_0 = 0.05
54 q = 0.5
55 p = 3.0 / 2.0
56 r0 = 1.0
57
58 # Viscosity parameter
59 alpha_AV = 1.0e-3 / 0.08
60 alpha_u = 1.0
61 beta_AV = 2.0
62
63 # Integrator parameters
64 C_cour = 0.3
65 C_force = 0.25
66
67 # Warp parameters
68 inclination = 30.0
69
70
71 # Disc profiles
72 def sigma_profile(r):
73 sigma_0 = 1.0 # We do not care as it will be renormalized
74 return sigma_0 * (r / r0) ** (-p)
75
76
77 def kep_profile(r):
78 return (G * center_mass / r) ** 0.5
79
80
81 def omega_k(r):
82 return kep_profile(r) / r
83
84
85 def cs_profile(r):
86 cs_in = (H_r_0 * r0) * omega_k(r0)
87 return ((r / r0) ** (-q)) * cs_in
88
89
90 def inc_profile(r): # profile of the inclination angle
91 if r < 10.0:
92 effective_inc = 0.0
93 elif r < 30.0 and r > 10:
94 lx = 0.1 * (1.0 + np.sin((r - 20) * np.pi / 20.0))
95 lz = np.sqrt(1 - lx * lx)
96 effective_inc = np.arccos(lz)
97
98 else:
99 lx = 0.2
100 lz = np.sqrt(1 - lx * lx)
101 effective_inc = np.arccos(lz)
102
103 return effective_inc
104
105
106 def psi_profile(r):
107 return 0.0
108
109
110 def k_profile(r): # profile of the warp direction
111 incl_rad = inc_profile(r) * np.pi / 180.0
112 return (0.0, np.cos(incl_rad), np.sin(incl_rad))
Utility functions and quantities deduced from the base one
118 # Deduced quantities
119 pmass = disc_mass / Npart
120 bmin = (-rout * 2, -rout * 2, -rout * 2)
121 bmax = (rout * 2, rout * 2, rout * 2)
122
123 cs0 = cs_profile(r0)
124
125
126 def rot_profile(r):
127 return ((kep_profile(r) ** 2) - (2 * p + q) * cs_profile(r) ** 2) ** 0.5
128
129
130 def H_profile(r):
131 H = cs_profile(r) / omega_k(r)
132 # fact = (2.**0.5) * 3. # factor taken from phantom, to fasten thermalizing
133 fact = 1.0
134 return fact * H
Start the context The context holds the data of the code We then init the layout of the field (e.g. the list of fields used by the solver)
142 ctx = shamrock.Context()
143 ctx.pdata_layout_new()
Attach a SPH model to the data and configure it
148 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
149
150 # Generate the default config
151 cfg = model.gen_default_config()
152 # Use disc alpha model viscosity
153 cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV)
154 # use the Lodato Price 2007 equation of state
155 cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q, r0=r0)
156 # Use the given code units
157 cfg.set_units(codeu)
158 # Change particle mass
159 cfg.set_particle_mass(pmass)
160 # Set the CFL
161 cfg.set_cfl_cour(C_cour)
162 cfg.set_cfl_force(C_force)
163
164 # Set the solver config to be the one stored in cfg
165 model.set_solver_config(cfg)
166
167 # Print the solver config
168 model.get_current_config().print_status()
169
170 # We want the patches to split above 10^8 part and merge if smaller than 1 part (e.g. disable patch)
171 model.init_scheduler(int(1e8), 1)
172
173 # Set the simulation box size
174 model.resize_simulation_box(bmin, bmax)
----- SPH Solver configuration -----
[
{
"artif_viscosity": {
"alpha_AV": 0.0125,
"alpha_u": 1.0,
"beta_AV": 2.0,
"type": "constant_disc"
},
"boundary_config": {
"bc_type": "free"
},
"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",
"cs0": 0.31394308776061347,
"eos_type": "locally_isothermal_lp07",
"q": 0.5,
"r0": 1.0
},
"epsilon_h": 1e-06,
"ext_force_config": {
"force_list": []
},
"gpart_mass": 1e-08,
"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": {
"unit_current": 1.0,
"unit_length": 149597870700.0,
"unit_lumint": 1.0,
"unit_mass": 1.98847e+30,
"unit_qte": 1.0,
"unit_temperature": 1.0,
"unit_time": 31536000.0
},
"use_two_stage_search": true
}
]
------------------------------------
Add the sink particle
179 # null position and velocity
180 model.add_sink(center_mass, (0, 0, 0), (0, 0, 0), center_racc)
Create the setup
185 setup = model.get_setup()
186 gen_disc = setup.make_generator_disc_mc(
187 part_mass=pmass,
188 disc_mass=disc_mass,
189 r_in=rin,
190 r_out=rout,
191 sigma_profile=sigma_profile,
192 H_profile=H_profile,
193 rot_profile=rot_profile,
194 cs_profile=cs_profile,
195 random_seed=666,
196 )
197
198 # apply the custom warp
199 warp = setup.make_modifier_custom_warp(
200 parent=gen_disc,
201 inc_profile=inc_profile,
202 psi_profile=psi_profile,
203 k_profile=k_profile,
204 )
205 # Print the dot graph of the setup
206 print(warp.get_dot())
Warning: make_generator_disc_mc: with the current EOS, cs_profile is ignored [SPHSetup][rank=0]
digraph G {
rankdir=LR;
node_0 [label="GeneratorMCDisc"];
node_1 [label="ApplyCustomWarp"];
node_0 -> node_1;
node_3 [label="Simulation"];
node_1 -> node_3;
}
Apply the setup
210 setup.apply_setup(warp)
SPH setup: generating particles ...
SPH setup: Nstep = 1000000 ( 1.0e+06 ) Ntotal = 1000000 ( 1.0e+06 rank min = 1.5e+05 max = 1.0e+06) rate = 1.000000e+06 N.s^-1
SPH setup: the generation step took : 6.511420137 s
SPH setup: final particle count = 1000000 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 : 8.95 us (54.8%)
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 882.00 ns (0.2%)
patch tree reduce : 2.12 us (0.6%)
gen split merge : 1.04 us (0.3%)
split / merge op : 0/0
apply split merge : 1.02 us (0.3%)
LB compute : 360.99 us (96.8%)
LB move op cnt : 0
LB apply : 4.29 us (1.1%)
Info: patch count stable after 1 runs npatch = 1 [DataInserterUtility][rank=0]
Info: --------------------------------------------- [DataInserterUtility][rank=0]
SPH setup: injected 1000000 / 1000000 => 100.0% | ranks with patchs = 1 / 1 <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.043881721000000005 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 | 2.2% 0.0% | 456.18 MB | 5.29 MB |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 6.752447324 s
Run a single step to init the integrator and smoothing length of the particles
Here the htolerance is the maximum factor of evolution of the smoothing length in each
Smoothing length iterations, increasing it affect the performance negatively but increase the
convergence rate of the smoothing length
this is why we increase it temporely to 1.3 before lowering it back to 1.1 (default value)
Note that both change_htolerance can be removed and it will work the same but would converge
more slowly at the first timestep
221 model.change_htolerance(1.3)
222 model.timestep()
223 model.change_htolerance(1.1)
Warning: .change_htolerance(val) is deprecated, [SPH][rank=0]
-> calling this is replaced internally by .change_htolerances(coarse=val, fine=min(val, 1.1))
see: https://shamrock-code.github.io/Shamrock/mkdocs/models/sph/smoothing_length_tolerance
---------------- t = 0, dt = 0 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 29.43 us (2.5%)
patch tree reduce : 7.34 us (0.6%)
gen split merge : 1.20 us (0.1%)
split / merge op : 0/0
apply split merge : 1.07 us (0.1%)
LB compute : 1.10 ms (93.9%)
LB move op cnt : 0
LB apply : 10.72 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.16 us (76.1%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.8074540479960924 unconverged cnt = 1000000
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0496902623949202 unconverged cnt = 1000000
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.3645973411133965 unconverged cnt = 1000000
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.7739765434474155 unconverged cnt = 999999
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 2.3061695064816403 unconverged cnt = 999999
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 2.9980203584261327 unconverged cnt = 999994
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 3.8974264659539726 unconverged cnt = 999986
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 5.066654405740165 unconverged cnt = 999965
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 5.444654479445718 unconverged cnt = 999894
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 5.444654479445718 unconverged cnt = 999617
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 5.444654479445718 unconverged cnt = 996421
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 5.444654479445718 unconverged cnt = 952609
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 5.444654479445718 unconverged cnt = 724147
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 5.444654479445718 unconverged cnt = 135019
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 5.444654479445718 unconverged cnt = 726
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3934716120709599e-05,-1.783643673409164e-05,-5.995703395873715e-07)
sum a = (9.012430558785756e-19,-6.437450399132683e-20,1.8422966602781032e-20)
sum e = 0.01353624682145266
sum de = -5.045801736573469e-22
Info: cfl dt = 7.19642848135117e-05 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.4388e+04 | 1000000 | 1 | 2.908e+01 | 0.0% | 0.1% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
Warning: .change_htolerance(val) is deprecated, [SPH][rank=0]
-> calling this is replaced internally by .change_htolerances(coarse=val, fine=min(val, 1.1))
see: https://shamrock-code.github.io/Shamrock/mkdocs/models/sph/smoothing_length_tolerance
Manipulating the simulation#
Dump files (path relative to where you have started shamrock)
230 dump_folder = "_to_trash"
231 import os
232
233 os.system("mkdir -p " + dump_folder)
234
235 # VTK dump
236 model.do_vtk_dump(dump_folder + "/init_disc.vtk", True)
237
238 # Shamrock restart dump files
239 model.dump(dump_folder + "/init_disc.sham")
240
241 # Phantom dump
242 dump = model.make_phantom_dump()
243 dump.save_dump(dump_folder + "/init_disc.phdump")
Info: dump to _to_trash/init_disc.vtk [VTK Dump][rank=0]
- took 98.93 ms, bandwidth = 566.04 MB/s
Info: Dumping state to _to_trash/init_disc.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 9.13 us (54.1%)
Info: dump to _to_trash/init_disc.sham [Shamrock Dump][rank=0]
- took 38.26 ms, bandwidth = 3.14 GB/s
Single timestep
247 model.evolve_once()
---------------- t = 0, dt = 7.19642848135117e-05 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 22.00 us (4.4%)
patch tree reduce : 2.59 us (0.5%)
gen split merge : 561.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1.14 us (0.2%)
LB compute : 453.75 us (90.5%)
LB move op cnt : 0
LB apply : 10.45 us (2.1%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.47 us (64.5%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3936202893870782e-05,-1.7835748317164945e-05,-5.995819807005979e-07)
sum a = (-1.0503208545953324e-18,4.641740550953566e-19,-2.2234614865425384e-20)
sum e = 0.013536247045198277
sum de = 4.642963559225607e-08
Info: cfl dt = 0.0024469434003790637 cfl multiplier : 0.34 [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.8827e+05 | 1000000 | 1 | 5.311e+00 | 0.0% | 0.1% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.048776445203471464 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 2 [SPH][rank=0]
Info: time since start : 473.675098246 (s) [SPH][rank=0]
Evolve until a given time (code units)
251 model.evolve_until(0.001)
Info: evolve_until (target_time = 0.00s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 7.19642848135117e-05, dt = 0.0009280357151864883 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 21.22 us (4.6%)
patch tree reduce : 2.80 us (0.6%)
gen split merge : 551.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1.04 us (0.2%)
LB compute : 413.76 us (89.7%)
LB move op cnt : 0
LB apply : 10.35 us (2.2%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.72 us (77.4%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3953886521147126e-05,-1.7827553476723058e-05,-5.997204608818979e-07)
sum a = (1.6263032587282567e-18,5.72594272343907e-19,2.8163845496205486e-20)
sum e = 0.013536284075170503
sum de = 6.449358070268885e-07
Info: cfl dt = 0.004033698305238818 cfl multiplier : 0.56 [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.7796e+05 | 1000000 | 1 | 5.619e+00 | 0.0% | 0.1% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5945368177571249 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 3 [SPH][rank=0]
Info: time since start : 479.54306319100004 (s) [SPH][rank=0]
EvolveUntilResults(reach_target_time=true, reach_niter_max=false, reach_max_walltime=false, iter_count=1)
Get the sinks positions
255 print(model.get_sinks())
[{'pos': (-2.065717359730298e-11, -9.571674293603349e-12, 1.6176240765935853e-13), 'velocity': (-4.1279143632583235e-08, -1.9215685349906975e-08, 3.2352509313430425e-10), 'sph_acceleration': (-4.123851286200948e-05, -1.9299173641009363e-05, 3.2352541355366475e-07), 'ext_acceleration': (0.0, 0.0, 0.0), 'mass': 1.0, 'angular_momentum': (0.0, 0.0, 0.0), 'accretion_radius': 0.1}]
Get the fields as python dictionary of numpy arrays
Warning
Do not do this on a large distributed simulation as this gather all data on MPI rank 0 and will use a lot of memory (and crash if the simulation is too large)
265 print(ctx.collect_data())
Info: collected : 1 patches [PatchScheduler][rank=0]
{'xyz': array([[-46.18081976, -38.71910747, -1.18240629],
[-45.84934467, -36.35079718, -1.0951954 ],
[-45.90089759, -30.16149923, -0.51797479],
...,
[ 34.5319474 , 48.88261936, 0.55061546],
[ 34.58505746, 49.21023716, 0.84423995],
[ 49.91803248, 33.88748886, 1.99642314]], shape=(1000000, 3)), 'vxyz': array([[ 0.51838326, -0.61505401, -0.10576831],
[ 0.50908461, -0.63898036, -0.10386844],
[ 0.46471919, -0.70559895, -0.09480775],
...,
[-0.65422734, 0.46066018, 0.13350951],
[-0.65420858, 0.45749044, 0.13350634],
[-0.45598511, 0.66621004, 0.09302728]], shape=(1000000, 3)), 'axyz': array([[ 0.00788434, 0.00659499, -0.00063158],
[ 0.00894982, 0.00722448, -0.00078459],
[ 0.01078347, 0.00750194, -0.00128493],
...,
[-0.00588592, -0.00841436, 0.00135547],
[-0.00558256, -0.0077336 , 0.00172636],
[-0.00863666, -0.00596705, 0.00034757]], shape=(1000000, 3)), 'axyz_ext': array([[ 0.00831349, 0.00697023, 0.00021286],
[ 0.00901889, 0.00715046, 0.00021543],
[ 0.0109206 , 0.00717593, 0.00012324],
...,
[-0.00634959, -0.00898833, -0.00010124],
[-0.00626412, -0.00891306, -0.00015291],
[-0.00894589, -0.00607303, -0.00035778]], shape=(1000000, 3)), 'hpart': array([3.46067577, 3.174374 , 2.93403838, ..., 1.75076501, 1.97237561,
4.40070917], shape=(1000000,)), 'uint': array([1.85772131e-09, 1.42328086e-09, 4.20630317e-10, ...,
4.39544051e-09, 3.77948668e-09, 8.78998943e-11], shape=(1000000,)), 'duint': array([1.85809089e-06, 1.42386071e-06, 4.21118893e-07, ...,
4.39469850e-06, 3.77840063e-06, 8.83487962e-08], shape=(1000000,))}
Performing a timestep loop
269 dt_stop = 0.001
270 for i in range(10):
271 t_target = i * dt_stop
272 # skip if the model is already past the target
273 if model.get_time() > t_target:
274 continue
275
276 model.evolve_until(i * dt_stop)
277
278 # Dump name is "dump_xxxx.sham" where xxxx is the timestep
279 model.dump(dump_folder + f"/dump_{i:04}.sham")
Info: evolve_until (target_time = 0.00s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
Info: iteration since start : 3 [SPH][rank=0]
Info: time since start : 480.417608405 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0001.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 26.81 us (72.2%)
Info: dump to _to_trash/dump_0001.sham [Shamrock Dump][rank=0]
- took 39.09 ms, bandwidth = 3.07 GB/s
Info: evolve_until (target_time = 0.00s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 0.001, dt = 0.001 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.71 us (1.8%)
patch tree reduce : 3.44 us (0.8%)
gen split merge : 841.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1.20 us (0.3%)
LB compute : 417.05 us (94.8%)
LB move op cnt : 0
LB apply : 4.04 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.52 us (70.0%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3955335377140358e-05,-1.7826787147273352e-05,-5.997321022941923e-07)
sum a = (-7.318364664277155e-19,-2.236166980751353e-19,-3.520480687025686e-20)
sum e = 0.013536290988573481
sum de = 1.2903327692788872e-06
Info: cfl dt = 0.005094995112141007 cfl multiplier : 0.7066666666666667 [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.7790e+05 | 1000000 | 1 | 5.621e+00 | 0.0% | 0.1% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.6404363656727183 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 4 [SPH][rank=0]
Info: time since start : 486.09188001900003 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0002.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.49 us (57.5%)
Info: dump to _to_trash/dump_0002.sham [Shamrock Dump][rank=0]
- took 38.56 ms, bandwidth = 3.11 GB/s
Info: evolve_until (target_time = 0.00s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 0.002, dt = 0.001 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 8.10 us (3.0%)
patch tree reduce : 1.89 us (0.7%)
gen split merge : 912.00 ns (0.3%)
split / merge op : 0/0
apply split merge : 1.21 us (0.4%)
LB compute : 250.35 us (91.9%)
LB move op cnt : 0
LB apply : 3.83 us (1.4%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.49 us (69.9%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3955294406859166e-05,-1.7826703866802762e-05,-5.997321020662395e-07)
sum a = (-2.0328790734103208e-19,5.285485590866834e-19,-5.447480642029219e-20)
sum e = 0.013536292598328123
sum de = 1.9362544215220714e-06
Info: cfl dt = 0.005805736233858829 cfl multiplier : 0.8044444444444444 [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.8625e+05 | 1000000 | 1 | 5.369e+00 | 0.0% | 0.0% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.6704895714784297 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 5 [SPH][rank=0]
Info: time since start : 491.51264808800005 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0003.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.80 us (57.7%)
Info: dump to _to_trash/dump_0003.sham [Shamrock Dump][rank=0]
- took 38.63 ms, bandwidth = 3.11 GB/s
Info: evolve_until (target_time = 0.00s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 0.003, dt = 0.001 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.96 us (2.4%)
patch tree reduce : 1.82 us (0.6%)
gen split merge : 811.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1.04 us (0.3%)
LB compute : 305.25 us (93.4%)
LB move op cnt : 0
LB apply : 4.09 us (1.3%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.41 us (68.8%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3955253086524537e-05,-1.7826620802866655e-05,-5.997321012681458e-07)
sum a = (6.844026213814747e-19,3.7269449679189215e-19,-2.9857911390714087e-20)
sum e = 0.013536294854504872
sum de = 2.582635730546344e-06
Info: cfl dt = 0.006282696726571655 cfl multiplier : 0.8696296296296296 [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.8592e+05 | 1000000 | 1 | 5.379e+00 | 0.0% | 0.0% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.6693252054352555 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 6 [SPH][rank=0]
Info: time since start : 496.94265623100006 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0004.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 8.31 us (59.8%)
Info: dump to _to_trash/dump_0004.sham [Shamrock Dump][rank=0]
- took 38.75 ms, bandwidth = 3.10 GB/s
Info: evolve_until (target_time = 0.01s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 0.004, dt = 0.001 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.71 us (2.5%)
patch tree reduce : 1.87 us (0.6%)
gen split merge : 952.00 ns (0.3%)
split / merge op : 0/0
apply split merge : 1.04 us (0.3%)
LB compute : 281.85 us (92.7%)
LB move op cnt : 0
LB apply : 4.48 us (1.5%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.35 us (67.2%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3955211417328168e-05,-1.7826537956975304e-05,-5.99732099897743e-07)
sum a = (-6.098637220230962e-20,6.606856988583543e-19,-1.937587866844212e-20)
sum e = 0.013536297757443975
sum de = 3.2294620889173226e-06
Info: cfl dt = 0.006603774399651343 cfl multiplier : 0.9130864197530864 [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.8674e+05 | 1000000 | 1 | 5.355e+00 | 0.0% | 0.0% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.672257067427579 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 7 [SPH][rank=0]
Info: time since start : 502.34961368200004 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0005.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.46 us (58.4%)
Info: dump to _to_trash/dump_0005.sham [Shamrock Dump][rank=0]
- took 38.93 ms, bandwidth = 3.08 GB/s
Info: evolve_until (target_time = 0.01s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 0.005, dt = 0.001 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 8.33 us (3.1%)
patch tree reduce : 1.84 us (0.7%)
gen split merge : 821.00 ns (0.3%)
split / merge op : 0/0
apply split merge : 1.11 us (0.4%)
LB compute : 245.36 us (91.8%)
LB move op cnt : 0
LB apply : 3.79 us (1.4%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.41 us (69.1%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.395516940045695e-05,-1.7826455330635906e-05,-5.997320979530626e-07)
sum a = (1.1248597539537109e-18,-1.7245590806097555e-18,-4.0763460586613204e-21)
sum e = 0.013536301307567025
sum de = 3.876719232410471e-06
Info: cfl dt = 0.006820935762349557 cfl multiplier : 0.9420576131687243 [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.8698e+05 | 1000000 | 1 | 5.348e+00 | 0.0% | 0.1% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.6731240632279657 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 8 [SPH][rank=0]
Info: time since start : 507.749302666 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0006.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.81 us (59.5%)
Info: dump to _to_trash/dump_0006.sham [Shamrock Dump][rank=0]
- took 39.00 ms, bandwidth = 3.08 GB/s
Info: evolve_until (target_time = 0.01s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 0.006, dt = 0.001 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.43 us (2.4%)
patch tree reduce : 2.10 us (0.7%)
gen split merge : 841.00 ns (0.3%)
split / merge op : 0/0
apply split merge : 1.06 us (0.3%)
LB compute : 286.12 us (93.2%)
LB move op cnt : 0
LB apply : 3.78 us (1.2%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.33 us (69.6%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3955127037116626e-05,-1.7826372925354343e-05,-5.997320954332403e-07)
sum a = (-9.283481101907132e-19,-7.589415207398531e-19,6.400392707690307e-20)
sum e = 0.013536305505279473
sum de = 4.524393001173605e-06
Info: cfl dt = 0.006968850096436335 cfl multiplier : 0.9613717421124829 [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.8586e+05 | 1000000 | 1 | 5.381e+00 | 0.0% | 0.0% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.6690824677194137 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 9 [SPH][rank=0]
Info: time since start : 513.182713616 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0007.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 8.01 us (60.2%)
Info: dump to _to_trash/dump_0007.sham [Shamrock Dump][rank=0]
- took 38.72 ms, bandwidth = 3.10 GB/s
Info: evolve_until (target_time = 0.01s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 0.007, dt = 0.001 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.95 us (3.0%)
patch tree reduce : 1.99 us (0.8%)
gen split merge : 952.00 ns (0.4%)
split / merge op : 0/0
apply split merge : 1.28 us (0.5%)
LB compute : 241.18 us (91.4%)
LB move op cnt : 0
LB apply : 4.16 us (1.6%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.26 us (66.4%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3955084328509234e-05,-1.782629074263273e-05,-5.997320923379383e-07)
sum a = (-2.507217523872729e-19,-1.0808140406964872e-18,-1.6464202912255463e-20)
sum e = 0.013536310350978223
sum de = 5.1724696773583754e-06
Info: cfl dt = 0.00707064604710029 cfl multiplier : 0.9742478280749886 [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.8503e+05 | 1000000 | 1 | 5.405e+00 | 0.0% | 0.1% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.6660979860045211 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 10 [SPH][rank=0]
Info: time since start : 518.638450295 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0008.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.60 us (59.1%)
Info: dump to _to_trash/dump_0008.sham [Shamrock Dump][rank=0]
- took 38.58 ms, bandwidth = 3.11 GB/s
Info: evolve_until (target_time = 0.01s, niter_max = -1, max_walltime = -1.00s) [SPH][rank=0]
---------------- t = 0.008, dt = 0.0010000000000000009 ----------------
Info: Summary (strategy = round robin): [LoadBalance][rank=0]
- strategy "psweep" : max = 1000000.0 min = 1000000.0 factor = 1
- strategy "round robin" : max = 950000.0 min = 950000.0 factor = 0.95
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 1000000
max = 1000000
avg = 1000000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 8.63 us (3.3%)
patch tree reduce : 2.09 us (0.8%)
gen split merge : 811.00 ns (0.3%)
split / merge op : 0/0
apply split merge : 1.06 us (0.4%)
LB compute : 243.15 us (91.7%)
LB move op cnt : 0
LB apply : 3.61 us (1.4%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1.94 us (68.8%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.3955041275850652e-05,-1.7826208783957416e-05,-5.997320886657341e-07)
sum a = (-1.6195269951502222e-18,3.1170812458958252e-19,-9.476181097407485e-21)
sum e = 0.013536315845049659
sum de = 5.820935866955099e-06
Info: cfl dt = 0.00714178873170889 cfl multiplier : 0.9828318853833258 [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.8621e+05 | 1000000 | 1 | 5.370e+00 | 0.0% | 0.1% 0.0% | 1.08 GB | 5.29 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.6703540133455435 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 11 [SPH][rank=0]
Info: time since start : 524.0597609700001 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0009.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.46 us (59.0%)
Info: dump to _to_trash/dump_0009.sham [Shamrock Dump][rank=0]
- took 39.34 ms, bandwidth = 3.05 GB/s
Plot column integrated density
283 import matplotlib.pyplot as plt
284
285 pixel_x = 1200
286 pixel_y = 1080
287 radius = 30
288 center = (0.0, 0.0, 0.0)
289
290 aspect = pixel_x / pixel_y
291 pic_range = [-radius * aspect, radius * aspect, -radius, radius]
292 delta_x = (radius * 2 * aspect, 0.0, 0.0)
293 delta_y = (0.0, 0.0, radius * 2)
294
295 arr_rho = model.render_cartesian_column_integ(
296 "rho", "f64", center=(0.0, 0.0, 0.0), delta_x=delta_x, delta_y=delta_y, nx=pixel_x, ny=pixel_y
297 )
298
299 import copy
300
301 import matplotlib
302
303 my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat")) # copy the default cmap
304 my_cmap.set_bad(color="black")
305
306 fig_width = 6
307 fig_height = fig_width / aspect
308 plt.figure(figsize=(fig_width, fig_height))
309 res = plt.imshow(arr_rho, cmap=my_cmap, origin="lower", extent=pic_range, norm="log", vmin=1e-9)
310
311 cbar = plt.colorbar(res, extend="both")
312 cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [code unit]")
313 # or r"$\rho$ [code unit]" for slices
314
315 plt.title("t = {:0.3f} [code unit]".format(model.get_time()))
316 plt.xlabel("x")
317 plt.ylabel("z")
318 plt.show()
![t = 0.009 [code unit]](../../_images/sphx_glr_run_sph_custom_warp_profile_001.png)
Info: compute_column_integ field_name: rho, rays count: 1296000 [sph::CartesianRender][rank=0]
Info: compute_column_integ took 26.22 s [sph::CartesianRender][rank=0]
Total running time of the script: (2 minutes 1.155 seconds)
Estimated memory usage: 1173 MB