Note
Go to the end to download the full example code.
Basic disc simulation#
This simple example shows how to run a basic disc simulation in SPH
9 import shamrock
10
11 # If we use the shamrock executable to run this script instead of the python interpreter,
12 # we should not initialize the system as the shamrock executable needs to handle specific MPI logic
13 if not shamrock.sys.is_initialized():
14 shamrock.change_loglevel(1)
15 shamrock.sys.init("0:0")
Setup units
21 si = shamrock.UnitSystem()
22 sicte = shamrock.Constants(si)
23 codeu = shamrock.UnitSystem(
24 unit_time=3600 * 24 * 365,
25 unit_length=sicte.au(),
26 unit_mass=sicte.sol_mass(),
27 )
28 ucte = shamrock.Constants(codeu)
29 G = ucte.G()
List parameters
35 # Resolution
36 Npart = 1000000
37
38 # Sink parameters
39 center_mass = 1.0
40 center_racc = 0.1
41
42 # Disc parameter
43 disc_mass = 0.01 # sol mass
44 rout = 10.0 # au
45 rin = 1.0 # au
46 H_r_0 = 0.05
47 q = 0.5
48 p = 3.0 / 2.0
49 r0 = 1.0
50
51 # Viscosity parameter
52 alpha_AV = 1.0e-3 / 0.08
53 alpha_u = 1.0
54 beta_AV = 2.0
55
56 # Integrator parameters
57 C_cour = 0.3
58 C_force = 0.25
59
60
61 # Disc profiles
62 def sigma_profile(r):
63 sigma_0 = 1.0 # We do not care as it will be renormalized
64 return sigma_0 * (r / r0) ** (-p)
65
66
67 def kep_profile(r):
68 return (G * center_mass / r) ** 0.5
69
70
71 def omega_k(r):
72 return kep_profile(r) / r
73
74
75 def cs_profile(r):
76 cs_in = (H_r_0 * r0) * omega_k(r0)
77 return ((r / r0) ** (-q)) * cs_in
Utility functions and quantities deduced from the base one
83 # Deduced quantities
84 pmass = disc_mass / Npart
85 bmin = (-rout * 2, -rout * 2, -rout * 2)
86 bmax = (rout * 2, rout * 2, rout * 2)
87
88 cs0 = cs_profile(r0)
89
90
91 def rot_profile(r):
92 return ((kep_profile(r) ** 2) - (2 * p + q) * cs_profile(r) ** 2) ** 0.5
93
94
95 def H_profile(r):
96 H = cs_profile(r) / omega_k(r)
97 # fact = (2.**0.5) * 3. # factor taken from phantom, to fasten thermalizing
98 fact = 1.0
99 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)
107 ctx = shamrock.Context()
108 ctx.pdata_layout_new()
Attach a SPH model to the data and configure it
113 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
114
115 # Generate the default config
116 cfg = model.gen_default_config()
117 # Use disc alpha model viscosity
118 cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV)
119 # use the Lodato Price 2007 equation of state
120 cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q, r0=r0)
121 # Use the given code units
122 cfg.set_units(codeu)
123 # Change particle mass
124 cfg.set_particle_mass(pmass)
125 # Set the CFL
126 cfg.set_cfl_cour(C_cour)
127 cfg.set_cfl_force(C_force)
128
129 # Set the solver config to be the one stored in cfg
130 model.set_solver_config(cfg)
131
132 # Print the solver config
133 model.get_current_config().print_status()
134
135 # We want the patches to split above 10^8 part and merge if smaller than 1 part (e.g. disable patch)
136 model.init_scheduler(int(1e8), 1)
137
138 # Set the simulation box size
139 model.resize_simulation_box(bmin, bmax)
----- SPH Solver configuration -----
units :
unit_length : 149597870700
unit_mass : 1.98847e+30
unit_current : 1
unit_temperature : 1
unit_qte : 1
unit_lumint : 1
part mass 1e-08 ( can be changed using .set_part_mass() )
cfl force 0.25
cfl courant 0.3
--- artificial viscosity config
Config Type : constant disc
alpha_AV = 0.0125
alpha_u = 1
beta_AV = 2
--- artificial viscosity config (deduced)
-------------
EOS config f64_3 :
locally isothermal (Lodato Price 2007) :
--- Bondaries config
Config Type : Free boundaries
--- Bondaries config config (deduced)
-------------
------------------------------------
Add the sink particle
144 # null position and velocity
145 model.add_sink(center_mass, (0, 0, 0), (0, 0, 0), center_racc)
Create the setup
150 setup = model.get_setup()
151 gen_disc = setup.make_generator_disc_mc(
152 part_mass=pmass,
153 disc_mass=disc_mass,
154 r_in=rin,
155 r_out=rout,
156 sigma_profile=sigma_profile,
157 H_profile=H_profile,
158 rot_profile=rot_profile,
159 cs_profile=cs_profile,
160 random_seed=666,
161 )
162
163 # Print the dot graph of the setup
164 print(gen_disc.get_dot())
digraph G {
rankdir=LR;
node_0 [label="GeneratorMCDisc"];
node_2 [label="Simulation"];
node_0 -> node_2;
}
Apply the setup
Info: pushing data in scheduler, N = 1000000 [DataInserterUtility][rank=0]
Info: reattributing data ... [DataInserterUtility][rank=0]
Info: reattributing data done in 50.15 ms [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.15 us (48.5%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 1022.00 ns (0.2%)
patch tree reduce : 9.57 us (1.9%)
gen split merge : 811.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 952.00 ns (0.2%)
LB compute : 497.16 us (96.5%)
LB move op cnt : 0
LB apply : 2.00 us (0.4%)
Info: the setup took : 2.9492781580000003 s [SPH setup][rank=0]
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 increse 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
179 model.change_htolerance(1.3)
180 model.timestep()
181 model.change_htolerance(1.1)
---------------- t = 0, dt = 0 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 24.61 us (2.7%)
patch tree reduce : 7.78 us (0.8%)
gen split merge : 812.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1002.00 ns (0.1%)
LB compute : 848.02 us (92.6%)
LB move op cnt : 0
LB apply : 2.04 us (0.2%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.03 us (64.4%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.30782750590945696 unconverged cnt = 999999
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.3566642885216895 unconverged cnt = 999999
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.4636635750781964 unconverged cnt = 999999
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.6027626476016553 unconverged cnt = 999997
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7251187699241742 unconverged cnt = 999987
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7333109890812106 unconverged cnt = 999937
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7529481567236113 unconverged cnt = 999775
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7529481567236118 unconverged cnt = 998210
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7529481567236117 unconverged cnt = 938876
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7529481567236117 unconverged cnt = 468324
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7529481567236117 unconverged cnt = 10415
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7529481567236117 unconverged cnt = 17
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 0, max = 9.350670822948002e-07
iterations = 3
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.877002876192637e-06,-3.3153642116852934e-05,0)
sum a = (-5.304120315706429e-18,2.846030702774449e-19,-1.2705494208814505e-21)
sum e = 0.05004044613227272
sum de = 6.4983308922165854e-21
Info: cfl dt = 3.907389496802996e-05 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 | 3.1474e+04 | 1000000 | 3.177e+01 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
Manipulating the simulation#
Dump files (path relative to where you have started shamrock)
188 dump_folder = "_to_trash"
189 import os
190
191 os.system("mkdir -p " + dump_folder)
192
193 # VTK dump
194 model.do_vtk_dump(dump_folder + "/init_disc.vtk", True)
195
196 # Shamrock restart dump files
197 model.dump(dump_folder + "/init_disc.sham")
198
199 # Phantom dump
200 dump = model.make_phantom_dump()
201 dump.save_dump(dump_folder + "/init_disc.phdump")
Info: dump to _to_trash/init_disc.vtk [VTK Dump][rank=0]
- took 158.36 ms, bandwidth = 337.24 MB/s
Info: Dumping state to _to_trash/init_disc.sham [SPH][rank=0]
Info: dump to _to_trash/init_disc.sham [Shamrock Dump][rank=0]
- took 126.68 ms, bandwidth = 903.44 MB/s
Single timestep
205 model.evolve_once()
---------------- t = 0, dt = 3.907389496802996e-05 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 12.66 us (1.5%)
patch tree reduce : 3.22 us (0.4%)
gen split merge : 591.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 972.00 ns (0.1%)
LB compute : 828.61 us (95.5%)
LB move op cnt : 0
LB apply : 2.54 us (0.3%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 3.00 us (75.9%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 1.0009746271268322e-14, max = 9.999963360720819e-07
iterations = 1
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.877275295883587e-06,-3.315171760280455e-05,4.520315299363933e-11)
sum a = (-1.5246593050577406e-18,1.0706496453294356e-18,-3.006966962752766e-20)
sum e = 0.05004044658210748
sum de = 1.6623973450838555e-07
Info: cfl dt = 0.0013286849078749257 cfl multiplier : 0.34 [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.7085e+05 | 1000000 | 5.853e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.02403323246049789 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 2 [SPH][rank=0]
Info: time since start : 69.466569154 (s) [SPH][rank=0]
Evolve until a given time (code units)
209 model.evolve_until(0.001)
---------------- t = 3.907389496802996e-05, dt = 0.0009609261050319701 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 19.43 us (2.0%)
patch tree reduce : 4.71 us (0.5%)
gen split merge : 731.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 952.00 ns (0.1%)
LB compute : 928.40 us (95.4%)
LB move op cnt : 0
LB apply : 2.37 us (0.2%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1964.00 ns (69.5%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 3.046203435208697e-14, max = 9.999794806296945e-07
iterations = 2
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.883694937533913e-06,-3.310631640840492e-05,1.111660091448701e-09)
sum a = (8.555032767268433e-19,1.0842021724855044e-19,-3.3881317890172014e-21)
sum e = 0.050040718349661346
sum de = 4.250271582630976e-06
Info: cfl dt = 0.002195483366595384 cfl multiplier : 0.56 [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.6297e+05 | 1000000 | 6.136e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.5637713653587878 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 3 [SPH][rank=0]
Info: time since start : 75.827457729 (s) [SPH][rank=0]
True
Get the sinks positions
213 print(model.get_sinks())
[{'pos': (-6.96448103523327e-12, -4.925022249208709e-11, -1.1568632444423742e-12), 'velocity': (-1.374628968486744e-08, -9.842692125776788e-08, -2.3137250214540352e-09), 'sph_acceleration': (-1.35487859454304e-05, -9.834742799767654e-05, -2.313723434691236e-06), '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)
223 print(ctx.collect_data())
adding -> xyz
adding -> vxyz
adding -> axyz
adding -> axyz_ext
adding -> hpart
adding -> uint
adding -> duint
{'xyz': array([[-7.65869271, -6.39642067, -0.95720353],
[-7.55994437, -6.50503034, -1.17167654],
[-7.51907833, -6.42069176, -1.15259277],
...,
[ 7.55442992, 6.37259754, 1.01837309],
[ 7.57542146, 6.28326004, 1.24365905],
[ 7.65988183, 6.31067646, 0.96346708]], shape=(1000000, 3)), 'vxyz': array([[ 1.26848748e+00, -1.51895096e+00, 3.12306504e-05],
[ 1.29103469e+00, -1.50060221e+00, -1.67518492e-05],
[ 1.29096064e+00, -1.51196362e+00, 1.36791610e-05],
...,
[-1.28208848e+00, 1.51994964e+00, -1.62093533e-07],
[-1.27207321e+00, 1.53374707e+00, 1.55803842e-05],
[-1.26170408e+00, 1.53154863e+00, 1.11059041e-05]],
shape=(1000000, 3)), 'axyz': array([[ 2.20260256e-01, 2.08844814e-01, 3.12257824e-02],
[ 2.35952850e-01, 1.26429677e-01, -1.67494312e-02],
[ 2.61260995e-01, 1.51087699e-01, 1.36744045e-02],
...,
[-2.50475689e-01, -2.32575157e-01, -1.60927066e-04],
[-2.64566804e-01, -2.42771666e-01, 1.55786797e-02],
[-2.51138687e-01, -2.23455574e-01, 1.11066195e-02]],
shape=(1000000, 3)), 'axyz_ext': array([[ 0.29974924, 0.25034588, 0.03746345],
[ 0.29432334, 0.25325348, 0.04561565],
[ 0.30052444, 0.25662384, 0.04606712],
...,
[-0.30365396, -0.25614964, -0.040934 ],
[-0.30590711, -0.25372766, -0.05022086],
[-0.30459965, -0.25094772, -0.03831283]], shape=(1000000, 3)), 'hpart': array([0.20771971, 0.24976131, 0.20235232, ..., 0.22745052, 0.30584142,
0.22650342], shape=(1000000,)), 'uint': array([ 3.40709038e-09, 9.06719995e-07, 6.33136482e-07, ...,
-4.58527033e-07, -2.14375395e-07, -4.80424414e-07],
shape=(1000000,)), 'duint': array([ 1.75537810e-06, 8.98428441e-04, 6.23714789e-04, ...,
-4.59495922e-04, -2.15374360e-04, -4.81095692e-04],
shape=(1000000,))}
Performing a timestep loop
227 dt_stop = 0.001
228 for i in range(10):
229 t_target = i * dt_stop
230 # skip if the model is already past the target
231 if model.get_time() > t_target:
232 continue
233
234 model.evolve_until(i * dt_stop)
235
236 # Dump name is "dump_xxxx.sham" where xxxx is the timestep
237 model.dump(dump_folder + f"/dump_{i:04}.sham")
Info: iteration since start : 3 [SPH][rank=0]
Info: time since start : 76.66359348 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0001.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0001.sham [Shamrock Dump][rank=0]
- took 129.19 ms, bandwidth = 885.91 MB/s
---------------- t = 0.001, dt = 0.001 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 6.31 us (0.5%)
patch tree reduce : 1583.00 ns (0.1%)
gen split merge : 671.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1112.00 ns (0.1%)
LB compute : 1292.79 us (98.6%)
LB move op cnt : 0
LB apply : 2.28 us (0.2%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1864.00 ns (65.7%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 3.59756982210602e-14, max = 9.999919608723482e-07
iterations = 2
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.883777269165263e-06,-3.3104468402850225e-05,1.1568617173455943e-09)
sum a = (3.4728350837426314e-19,-3.781155076543197e-18,8.216219588366713e-20)
sum e = 0.05004074705060449
sum de = 8.502063120608242e-06
Info: cfl dt = 0.002784776400983339 cfl multiplier : 0.7066666666666667 [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.6368e+05 | 1000000 | 6.109e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.5892529663416624 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 4 [SPH][rank=0]
Info: time since start : 82.922714908 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0002.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0002.sham [Shamrock Dump][rank=0]
- took 126.21 ms, bandwidth = 906.76 MB/s
---------------- t = 0.002, dt = 0.001 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 6.55 us (1.2%)
patch tree reduce : 1153.00 ns (0.2%)
gen split merge : 691.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1062.00 ns (0.2%)
LB compute : 511.97 us (96.6%)
LB move op cnt : 0
LB apply : 2.28 us (0.4%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1633.00 ns (66.0%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 3.928401332424843e-14, max = 9.999915952970984e-07
iterations = 2
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.883580481113704e-06,-3.3104548693157906e-05,1.1568571334709395e-09)
sum a = (-2.97816784254612e-18,-9.351243737687476e-19,-1.6517142471458857e-20)
sum e = 0.05004075760279957
sum de = 1.275450403348239e-05
Info: cfl dt = 0.003193940312031096 cfl multiplier : 0.8044444444444444 [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.6425e+05 | 1000000 | 6.088e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.5913172349082063 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 5 [SPH][rank=0]
Info: time since start : 89.154881605 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0003.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0003.sham [Shamrock Dump][rank=0]
- took 126.32 ms, bandwidth = 906.03 MB/s
---------------- t = 0.003, dt = 0.001 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 6.24 us (1.2%)
patch tree reduce : 1172.00 ns (0.2%)
gen split merge : 612.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1322.00 ns (0.3%)
LB compute : 499.43 us (96.5%)
LB move op cnt : 0
LB apply : 2.19 us (0.4%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1703.00 ns (65.4%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 4.314372909994344e-14, max = 9.999938675557479e-07
iterations = 2
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.883384425660589e-06,-3.310462979107857e-05,1.1568495020018086e-09)
sum a = (-1.5670109524204556e-18,-2.2090619264392153e-18,2.541098841762901e-20)
sum e = 0.050040772405409085
sum de = 1.700734113447608e-05
Info: cfl dt = 0.003486138014016258 cfl multiplier : 0.8696296296296296 [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.6719e+05 | 1000000 | 5.981e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.6018941014506963 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 6 [SPH][rank=0]
Info: time since start : 95.277418072 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0004.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0004.sham [Shamrock Dump][rank=0]
- took 126.61 ms, bandwidth = 903.90 MB/s
---------------- t = 0.004, dt = 0.001 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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.06 us (1.5%)
patch tree reduce : 1463.00 ns (0.3%)
gen split merge : 602.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1062.00 ns (0.2%)
LB compute : 451.12 us (96.0%)
LB move op cnt : 0
LB apply : 2.23 us (0.5%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1763.00 ns (68.7%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 4.6865656022720254e-14, max = 9.999306744900812e-07
iterations = 2
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.883189105990373e-06,-3.3104711691901545e-05,1.1568388248212036e-09)
sum a = (-2.715587628897287e-18,-3.0493186101154812e-18,-8.512681119905718e-20)
sum e = 0.0500407914555922
sum de = 2.1260579794904812e-05
Info: cfl dt = 0.003702564278479169 cfl multiplier : 0.9130864197530864 [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.6805e+05 | 1000000 | 5.951e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.6049780424631104 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 7 [SPH][rank=0]
Info: time since start : 101.37135826100001 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0005.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0005.sham [Shamrock Dump][rank=0]
- took 129.47 ms, bandwidth = 883.94 MB/s
---------------- t = 0.005, dt = 0.001 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 6.36 us (1.1%)
patch tree reduce : 1222.00 ns (0.2%)
gen split merge : 671.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1062.00 ns (0.2%)
LB compute : 531.79 us (93.9%)
LB move op cnt : 0
LB apply : 2.04 us (0.4%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1714.00 ns (69.0%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 5.169036272632375e-14, max = 9.999633369052757e-07
iterations = 2
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.882994525255597e-06,-3.310479439090972e-05,1.1568250982440996e-09)
sum a = (-3.4558944247975454e-19,-1.4230153513872246e-18,-2.541098841762901e-20)
sum e = 0.050040814752549914
sum de = 2.551421853648265e-05
Info: cfl dt = 0.0038700311671878286 cfl multiplier : 0.9420576131687243 [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.6793e+05 | 1000000 | 5.955e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.6045627296353858 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 8 [SPH][rank=0]
Info: time since start : 107.474370531 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0006.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0006.sham [Shamrock Dump][rank=0]
- took 125.65 ms, bandwidth = 910.82 MB/s
---------------- t = 0.006, dt = 0.001 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 6.60 us (0.6%)
patch tree reduce : 1393.00 ns (0.1%)
gen split merge : 671.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1112.00 ns (0.1%)
LB compute : 1123.51 us (98.3%)
LB move op cnt : 0
LB apply : 2.13 us (0.2%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1793.00 ns (73.4%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 5.118853328034394e-14, max = 9.9998416943092e-07
iterations = 3
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.882800686588224e-06,-3.310487788339004e-05,1.1568083141664227e-09)
sum a = (-5.793705359219414e-18,2.642742795433417e-18,-1.1434944787933055e-20)
sum e = 0.05004084229606356
sum de = 2.9768252145916814e-05
Info: cfl dt = 0.004005966807028481 cfl multiplier : 0.9613717421124829 [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.6773e+05 | 1000000 | 5.962e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.6038234070965299 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 9 [SPH][rank=0]
Info: time since start : 113.57980882800001 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0007.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0007.sham [Shamrock Dump][rank=0]
- took 126.13 ms, bandwidth = 907.39 MB/s
---------------- t = 0.007, dt = 0.001 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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 : 6.12 us (0.5%)
patch tree reduce : 1162.00 ns (0.1%)
gen split merge : 711.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1122.00 ns (0.1%)
LB compute : 1239.17 us (98.6%)
LB move op cnt : 0
LB apply : 2.28 us (0.2%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1874.00 ns (72.5%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 6.175350975866222e-14, max = 9.999783138459902e-07
iterations = 3
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.882607593081656e-06,-3.310496216459064e-05,1.1567884609152483e-09)
sum a = (1.8295911660692887e-19,3.496552006265752e-18,8.512681119905718e-20)
sum e = 0.050040874086274195
sum de = 3.40226702527547e-05
Info: cfl dt = 0.004116601548113821 cfl multiplier : 0.9742478280749886 [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.6662e+05 | 1000000 | 6.002e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.5998454326287141 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 10 [SPH][rank=0]
Info: time since start : 119.72302468800001 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0008.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0008.sham [Shamrock Dump][rank=0]
- took 149.22 ms, bandwidth = 766.98 MB/s
---------------- t = 0.008, dt = 0.0010000000000000009 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 1000000 min = 1000000 [LoadBalance][rank=0]
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.05 us (0.3%)
patch tree reduce : 1352.00 ns (0.0%)
gen split merge : 631.00 ns (0.0%)
split / merge op : 0/0
apply split merge : 1102.00 ns (0.0%)
LB compute : 2.80 ms (99.3%)
LB move op cnt : 0
LB apply : 2.14 us (0.1%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1934.00 ns (67.5%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: smoothing length iteration converged [Smoothinglength][rank=0]
eps min = 6.76811789489175e-14, max = 9.999720767879178e-07
iterations = 3
Info: conservation infos : [sph::Model][rank=0]
sum v = (9.882415247814464e-06,-3.310504722978855e-05,1.1567655234357156e-09)
sum a = (-7.758821796849391e-19,-1.2739375526704677e-18,-1.5246593050577406e-20)
sum e = 0.05004091012360171
sum de = 3.827745882709571e-05
Info: cfl dt = 0.003921166513495738 cfl multiplier : 0.9828318853833258 [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.6536e+05 | 1000000 | 6.048e+00 | 0 % | 0 % | 950.68 MB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0.5952870001124383 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 11 [SPH][rank=0]
Info: time since start : 125.94164024300001 (s) [SPH][rank=0]
Info: Dumping state to _to_trash/dump_0009.sham [SPH][rank=0]
Info: dump to _to_trash/dump_0009.sham [Shamrock Dump][rank=0]
- took 126.98 ms, bandwidth = 901.28 MB/s
Plot column integrated density
241 import matplotlib.pyplot as plt
242
243 pixel_x = 1200
244 pixel_y = 1080
245 radius = 5
246 center = (0.0, 0.0, 0.0)
247
248 aspect = pixel_x / pixel_y
249 pic_range = [-radius * aspect, radius * aspect, -radius, radius]
250 delta_x = (radius * 2 * aspect, 0.0, 0.0)
251 delta_y = (0.0, radius * 2, 0.0)
252
253 arr_rho = model.render_cartesian_column_integ(
254 "rho", "f64", center=(0.0, 0.0, 0.0), delta_x=delta_x, delta_y=delta_y, nx=pixel_x, ny=pixel_y
255 )
256
257 import copy
258
259 import matplotlib
260
261 my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat")) # copy the default cmap
262 my_cmap.set_bad(color="black")
263
264 fig_width = 6
265 fig_height = fig_width / aspect
266 plt.figure(figsize=(fig_width, fig_height))
267 res = plt.imshow(arr_rho, cmap=my_cmap, origin="lower", extent=pic_range, norm="log", vmin=1e-9)
268
269 cbar = plt.colorbar(res, extend="both")
270 cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [code unit]")
271 # or r"$\rho$ [code unit]" for slices
272
273 plt.title("t = {:0.3f} [code unit]".format(model.get_time()))
274 plt.xlabel("x")
275 plt.ylabel("z")
276 plt.show()
![t = 0.009 [code unit]](../_images/sphx_glr_run_sph_basic_disc_001.png)
Info: compute_column_integ field_name: rho, center: (0,0,0), delta_x: (11.11111111111111,0,0), delta_y: (0,10,0), nx: 1200, ny: 1080 [sph::CartesianRender][rank=0]
Info: compute_column_integ took 9.46 s [sph::CartesianRender][rank=0]
Plot vertical profiles at r=1
281 import numpy as np
282
283 dat = ctx.collect_data()
284
285 for rcenter in [1.0, 2.0, 3.0]:
286
287 z = []
288 h = []
289 vz = []
290 az = []
291
292 delta_r = 0.01
293
294 for i in range(len(dat["xyz"])):
295 r = (dat["xyz"][i][0] ** 2 + dat["xyz"][i][1] ** 2) ** 0.5
296 if r < rcenter + delta_r and r > rcenter - delta_r:
297 z.append(dat["xyz"][i][2])
298 h.append(dat["hpart"][i])
299 vz.append(dat["vxyz"][i][2])
300 az.append(dat["axyz"][i][2])
301
302 rho = pmass * (model.get_hfact() / np.array(h)) ** 3
303
304 fig, axs = plt.subplots(nrows=3, ncols=1, sharex=True)
305
306 from scipy.optimize import curve_fit
307
308 def func(x, a, c):
309 return a * np.exp(-((x / c) ** 2) / 2)
310
311 rho_0 = 0.001
312 p0 = [rho_0, H_profile(rcenter)] # a, b, c
313 popt, pcov = curve_fit(func, z, rho, p0=p0)
314
315 z_ana = np.linspace(-5.0 * H_profile(rcenter), 5.0 * H_profile(rcenter), 100)
316 rho_fit = func(z_ana, *popt)
317
318 axs[0].scatter(z, rho, label="rho")
319
320 axs[0].plot(z_ana, rho_fit, c="black", label="gaussian fit")
321 stddev = abs(popt[1])
322 axs[0].annotate(
323 f"Stddev: {stddev:.5f}",
324 xy=(0.05, 0.95),
325 xycoords="axes fraction",
326 fontsize=10,
327 verticalalignment="top",
328 bbox=dict(boxstyle="round", fc="w"),
329 )
330
331 axs[0].set_ylabel("rho")
332 axs[0].legend()
333
334 axs[1].scatter(z, vz, label="vz")
335
336 vz_fit = np.polyfit(z, vz, 1)
337 vz_fit_fn = np.poly1d(vz_fit)
338 axs[1].plot(z_ana, vz_fit_fn(z_ana), c="red", label="linear fit")
339
340 axs[1].set_ylabel("vz")
341 axs[1].legend()
342
343 axs[2].scatter(z, az, label="az")
344
345 az_fit = np.polyfit(z, az, 1)
346 az_fit_fn = np.poly1d(az_fit)
347 print(f"r={rcenter} az_fit={az_fit}")
348 axs[2].plot(z_ana, az_fit_fn(z_ana), c="red", label="linear fit")
349
350 axs[2].set_ylabel("az")
351 axs[2].set_xlabel("z")
352 axs[2].legend()
353
354 plt.tight_layout()
355 plt.show()
adding -> xyz
adding -> vxyz
adding -> axyz
adding -> axyz_ext
adding -> hpart
adding -> uint
adding -> duint
r=1.0 az_fit=[ 2.6665086 -0.0194775]
r=2.0 az_fit=[ 0.14429902 -0.00043202]
r=3.0 az_fit=[0.01686233 0.00431511]
Total running time of the script: (1 minutes 53.087 seconds)
Estimated memory usage: 1119 MB