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 -----
[
{
"artif_viscosity": {
"alpha_AV": 0.0125,
"alpha_u": 1.0,
"av_type": "constant_disc",
"beta_AV": 2.0
},
"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,
"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,
"self_grav_config": {
"softening_length": 1e-09,
"softening_mode": "plummer",
"type": "none"
},
"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
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
SPH setup: generating particles ...
SPH setup: Nstep = 1000000 ( 1.0e+06 ) Ntotal = 1000000 ( 1.0e+06 ) rate = 2.776953e+05 N.s^-1
SPH setup: the generation step took : 3.7298867390000003 s
SPH setup: final particle count = 1000000 begining 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 : 18.13 us (51.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 : 1954.00 ns (0.3%)
patch tree reduce : 2.21 us (0.4%)
gen split merge : 942.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1062.00 ns (0.2%)
LB compute : 564.79 us (95.8%)
LB move op cnt : 0
LB apply : 13.51 us (2.3%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 3.44 us (69.4%)
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 : 1302.00 ns (0.3%)
patch tree reduce : 952.00 ns (0.2%)
gen split merge : 501.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 421.00 ns (0.1%)
LB compute : 374.34 us (97.8%)
LB move op cnt : 0
LB apply : 2.27 us (0.6%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.96 us (68.6%)
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 : 1283.00 ns (0.3%)
patch tree reduce : 461.00 ns (0.1%)
gen split merge : 391.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 310.00 ns (0.1%)
LB compute : 422.49 us (98.3%)
LB move op cnt : 0
LB apply : 1914.00 ns (0.4%)
Info: --------------------------------------------- [DataInserterUtility][rank=0]
SPH setup: injected 1000000 / 1000000 => 100.0% | ranks with patchs = 1 / 1 <- global loop ->
SPH setup: the injection step took : 0.082804497 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 | 0.3% 0.0% | 1.01 GB | 5.04 MB |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 3.9875119910000003 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 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)
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 : [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 : 59.70 us (3.5%)
patch tree reduce : 1733.00 ns (0.1%)
gen split merge : 961.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1031.00 ns (0.1%)
LB compute : 1630.34 us (94.4%)
LB move op cnt : 0
LB apply : 12.86 us (0.7%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.83 us (70.3%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.3301323360955917 unconverged cnt = 1000000
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.42917203692426925 unconverged cnt = 1000000
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.55792364800155 unconverged cnt = 999998
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.725300742402015 unconverged cnt = 999992
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9037993151319754 unconverged cnt = 999978
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9037993151318011 unconverged cnt = 999932
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9037993151318012 unconverged cnt = 999768
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9037993151318012 unconverged cnt = 998188
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9037993151318012 unconverged cnt = 939155
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9037993151318012 unconverged cnt = 467894
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9037993151318012 unconverged cnt = 10399
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9037993151318012 unconverged cnt = 13
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4720599584510802e-05,-2.9970972851707412e-05,0)
sum a = (-2.0735366548785272e-18,1.666960840196463e-18,-2.864559548495637e-19)
sum e = 0.04998439673769593
sum de = 3.9367179712623693e-20
Info: CFL hydro = 4.476469486023067e-05 sink sink = inf [SPH][rank=0]
Info: cfl dt = 4.476469486023067e-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 | 2.9052e+04 | 1000000 | 1 | 3.442e+01 | 0.0% | 0.1% 0.0% | 1.16 GB | 5.04 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)
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 203.71 ms, bandwidth = 262.17 MB/s
Info: Dumping state to _to_trash/init_disc.sham [SPH][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 23.88 us (50.4%)
Info: dump to _to_trash/init_disc.sham [Shamrock Dump][rank=0]
- took 275.56 ms, bandwidth = 415.32 MB/s
Single timestep
205 model.evolve_once()
---------------- t = 0, dt = 4.476469486023067e-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 : 20.10 us (1.9%)
patch tree reduce : 4.34 us (0.4%)
gen split merge : 712.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 902.00 ns (0.1%)
LB compute : 986.39 us (94.1%)
LB move op cnt : 0
LB apply : 13.39 us (1.3%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.90 us (71.0%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4722215866676583e-05,-2.9969522809083586e-05,-1.0341650922953199e-11)
sum a = (7.724940478959219e-19,1.8431436932253575e-18,5.0769037276054627e-20)
sum e = 0.04998439732487808
sum de = 1.8894119296603306e-07
Info: CFL hydro = 0.00152208076838718 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.00152208076838718 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.5239e+05 | 1000000 | 1 | 6.562e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.024557601898055503 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 2 [SPH][rank=0]
Info: time since start : 2206.5433387860003 (s) [SPH][rank=0]
Evolve until a given time (code units)
209 model.evolve_until(0.001)
---------------- t = 4.476469486023067e-05, dt = 0.0009552353051397694 ----------------
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.94 us (1.9%)
patch tree reduce : 4.23 us (0.4%)
gen split merge : 752.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1102.00 ns (0.1%)
LB compute : 983.39 us (94.6%)
LB move op cnt : 0
LB apply : 13.74 us (1.3%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.94 us (71.5%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4755084417719843e-05,-2.9940022157721663e-05,-2.2068082113797195e-10)
sum a = (7.928228386300251e-18,-3.4558944247975454e-18,-4.489274620447792e-20)
sum e = 0.04998466429562942
sum de = 4.216635051622411e-06
Info: CFL hydro = 0.002509732122394582 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.002509732122394582 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.4641e+05 | 1000000 | 1 | 6.830e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5034798916470491 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 3 [SPH][rank=0]
Info: time since start : 2213.87810743 (s) [SPH][rank=0]
True
Get the sinks positions
213 print(model.get_sinks())
[{'pos': (-3.6101115374534386e-11, -3.2400736619587124e-11, 2.310224720609737e-13), 'velocity': (-7.209339750784241e-08, -6.497631486829308e-08, 4.620398839606087e-10), 'sph_acceleration': (-7.197438612373167e-05, -6.516750756591994e-05, 4.620343495805267e-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)
223 print(ctx.collect_data())
Info: collected : 1 patches [PatchScheduler][rank=0]
{'xyz': array([[-7.51193244, -6.45496337, -1.73444686],
[-7.50762886, -6.58717578, -0.88969697],
[-7.63272302, -6.4579532 , -0.94683626],
...,
[ 7.69149542, 6.31900019, 0.72686646],
[ 7.71861214, 6.28726935, 0.72867766],
[ 7.64032675, 6.43537698, 0.94725832]], shape=(1000000, 3)), 'vxyz': array([[ 1.29455994e+00, -1.50658397e+00, 3.42125994e-05],
[ 1.30443652e+00, -1.48694427e+00, 3.26449010e-05],
[ 1.27689782e+00, -1.50927741e+00, -1.49121307e-05],
...,
[-1.25774374e+00, 1.53113665e+00, 1.14824075e-04],
[-1.25120684e+00, 1.53627883e+00, 1.66517999e-04],
[-1.27405045e+00, 1.51281387e+00, -2.17082271e-05]],
shape=(1000000, 3)), 'axyz': array([[ 0.28210387, 0.23687329, 0.03421215],
[ 0.15451762, 0.18440075, 0.03264271],
[ 0.2108608 , 0.25759702, -0.01491171],
...,
[-0.27554389, -0.07663619, 0.11481269],
[-0.22472138, -0.12290354, 0.16650197],
[-0.15588396, -0.2105929 , -0.02170517]], shape=(1000000, 3)), 'axyz_ext': array([[ 0.29131382, 0.25032441, 0.0672621 ],
[ 0.2935698 , 0.25757744, 0.0347897 ],
[ 0.29707228, 0.25134921, 0.0368517 ],
...,
[-0.30497984, -0.25055825, -0.02882139],
[-0.30596061, -0.24922314, -0.02888429],
[-0.29814007, -0.25112064, -0.03696382]], shape=(1000000, 3)), 'hpart': array([0.56978443, 0.17686347, 0.24119288, ..., 0.15313134, 0.1440687 ,
0.21701926], shape=(1000000,)), 'uint': array([ 1.28156560e-07, -4.62942584e-07, -8.22843176e-07, ...,
2.03145077e-07, 2.98545507e-08, -5.58217322e-07],
shape=(1000000,)), 'duint': array([ 1.28662390e-04, -4.71658389e-04, -8.24382663e-04, ...,
1.84440340e-04, 5.76388533e-06, -5.63935655e-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 : 2215.922664916 (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 : 20.95 us (61.1%)
Info: dump to _to_trash/dump_0001.sham [Shamrock Dump][rank=0]
- took 272.61 ms, bandwidth = 419.82 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 : 7.06 us (0.8%)
patch tree reduce : 5.71 us (0.7%)
gen split merge : 601.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1263.00 ns (0.1%)
LB compute : 806.81 us (95.6%)
LB move op cnt : 0
LB apply : 16.46 us (1.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 3.72 us (71.6%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.475658677756969e-05,-2.993838909793012e-05,-2.310171747903361e-10)
sum a = (5.285485590866834e-19,-8.809142651444724e-19,-6.569799297141167e-20)
sum e = 0.049984695997130825
sum de = 8.434673094564203e-06
Info: CFL hydro = 0.003171108665537212 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.003171108665537212 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.5053e+05 | 1000000 | 1 | 6.643e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5419210792007849 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 4 [SPH][rank=0]
Info: time since start : 2222.9298788710003 (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.13 us (57.8%)
Info: dump to _to_trash/dump_0002.sham [Shamrock Dump][rank=0]
- took 105.02 ms, bandwidth = 1.06 GB/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 : 8.79 us (1.7%)
patch tree reduce : 2.22 us (0.4%)
gen split merge : 751.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1313.00 ns (0.3%)
LB compute : 489.52 us (95.3%)
LB move op cnt : 0
LB apply : 4.49 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 3.45 us (73.7%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4756466685072136e-05,-2.993819840663045e-05,-2.3100125652963938e-10)
sum a = (-2.100641709190665e-18,9.012430558785756e-18,-1.4097804615863761e-19)
sum e = 0.04998470646505142
sum de = 1.2653392893681442e-05
Info: CFL hydro = 0.003616692252102676 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.003616692252102676 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.5704e+05 | 1000000 | 1 | 6.368e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5653432759024074 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 5 [SPH][rank=0]
Info: time since start : 2229.4410987250003 (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.27 us (58.7%)
Info: dump to _to_trash/dump_0003.sham [Shamrock Dump][rank=0]
- took 106.98 ms, bandwidth = 1.04 GB/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 : 7.42 us (1.4%)
patch tree reduce : 1814.00 ns (0.3%)
gen split merge : 771.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1001.00 ns (0.2%)
LB compute : 516.12 us (95.9%)
LB move op cnt : 0
LB apply : 4.75 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.65 us (67.8%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4756345490273128e-05,-2.993800823390489e-05,-2.309747456286761e-10)
sum a = (-2.927345865710862e-18,4.87890977618477e-18,-1.938117262436246e-19)
sum e = 0.049984721149974654
sum de = 1.687250895006369e-05
Info: CFL hydro = 0.003919751565480149 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.003919751565480149 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.5919e+05 | 1000000 | 1 | 6.282e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5731014667999377 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 6 [SPH][rank=0]
Info: time since start : 2235.8687744020003 (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 : 6.97 us (55.8%)
Info: dump to _to_trash/dump_0004.sham [Shamrock Dump][rank=0]
- took 106.44 ms, bandwidth = 1.05 GB/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 : 8.20 us (1.8%)
patch tree reduce : 2.11 us (0.5%)
gen split merge : 801.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 921.00 ns (0.2%)
LB compute : 438.10 us (95.1%)
LB move op cnt : 0
LB apply : 3.96 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.34 us (70.5%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4756223195670914e-05,-2.9937818585934038e-05,-2.309376666322365e-10)
sum a = (-1.3417001884508117e-18,5.1906179007743525e-18,1.5415999640028266e-19)
sum e = 0.049984740048780285
sum de = 2.1092031617369558e-05
Info: CFL hydro = 0.004128824248681664 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.004128824248681664 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.5847e+05 | 1000000 | 1 | 6.310e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5704927388524883 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 7 [SPH][rank=0]
Info: time since start : 2242.324977003 (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.79 us (58.7%)
Info: dump to _to_trash/dump_0005.sham [Shamrock Dump][rank=0]
- took 105.88 ms, bandwidth = 1.06 GB/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 : 7.70 us (1.5%)
patch tree reduce : 2.07 us (0.4%)
gen split merge : 742.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1012.00 ns (0.2%)
LB compute : 483.42 us (95.6%)
LB move op cnt : 0
LB apply : 3.78 us (0.7%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.46 us (70.8%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4756099803799663e-05,-2.9937629468929124e-05,-2.3089005091855867e-10)
sum a = (4.336808689942018e-18,-1.8973538018496328e-19,-5.0186702124817295e-20)
sum e = 0.04998476316059502
sum de = 2.5311964764406573e-05
Info: CFL hydro = 0.004276156303829419 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.004276156303829419 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.5878e+05 | 1000000 | 1 | 6.298e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5715917940925462 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 8 [SPH][rank=0]
Info: time since start : 2248.767500923 (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.38 us (58.1%)
Info: dump to _to_trash/dump_0006.sham [Shamrock Dump][rank=0]
- took 105.40 ms, bandwidth = 1.06 GB/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 : 7.18 us (1.5%)
patch tree reduce : 2.19 us (0.5%)
gen split merge : 761.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1002.00 ns (0.2%)
LB compute : 448.20 us (95.4%)
LB move op cnt : 0
LB apply : 4.26 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.67 us (72.6%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4755975317245585e-05,-2.9937440889060213e-05,-2.3083193655094437e-10)
sum a = (-3.7947076036992655e-19,4.106415728288848e-18,3.5310685988663645e-20)
sum e = 0.049984790485087
sum de = 2.95323055354425e-05
Info: CFL hydro = 0.004383163238081959 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.004383163238081959 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.5778e+05 | 1000000 | 1 | 6.338e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5679983745186606 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 9 [SPH][rank=0]
Info: time since start : 2255.249423775 (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 : 6.42 us (55.0%)
Info: dump to _to_trash/dump_0007.sham [Shamrock Dump][rank=0]
- took 106.04 ms, bandwidth = 1.05 GB/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 : 7.88 us (1.8%)
patch tree reduce : 2.26 us (0.5%)
gen split merge : 701.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1042.00 ns (0.2%)
LB compute : 425.34 us (95.1%)
LB move op cnt : 0
LB apply : 4.05 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.54 us (68.3%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4755849738620548e-05,-2.9937252852497073e-05,-2.3076336697205662e-10)
sum a = (5.312590645178972e-18,2.6291902682773483e-18,-6.225692162319107e-20)
sum e = 0.0499848220223755
sum de = 3.375304807027477e-05
Info: CFL hydro = 0.004464056200393199 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.004464056200393199 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.5799e+05 | 1000000 | 1 | 6.330e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5687595877077484 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 10 [SPH][rank=0]
Info: time since start : 2261.723394776 (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 : 6.76 us (54.3%)
Info: dump to _to_trash/dump_0008.sham [Shamrock Dump][rank=0]
- took 109.45 ms, bandwidth = 1.02 GB/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 : 8.03 us (1.7%)
patch tree reduce : 1974.00 ns (0.4%)
gen split merge : 772.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 932.00 ns (0.2%)
LB compute : 447.22 us (95.3%)
LB move op cnt : 0
LB apply : 3.96 us (0.8%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.46 us (70.6%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4755723070605262e-05,-2.9937065365406272e-05,-2.3068439026203612e-10)
sum a = (-5.44811591673966e-18,1.328147661294743e-18,-2.2324612116071153e-19)
sum e = 0.049984857772888094
sum de = 3.7974183586963e-05
Info: CFL hydro = 0.004528262027107737 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.004528262027107737 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.5789e+05 | 1000000 | 1 | 6.334e+00 | 0.0% | 0.0% 0.0% | 1.16 GB | 5.04 MB |
+------+------------+---------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0.5683888140667736 (tsim/hr) [sph::Model][rank=0]
Info: iteration since start : 11 [SPH][rank=0]
Info: time since start : 2268.204827471 (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 : 6.30 us (55.7%)
Info: dump to _to_trash/dump_0009.sham [Shamrock Dump][rank=0]
- took 106.28 ms, bandwidth = 1.05 GB/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, rays count: 1296000 [sph::CartesianRender][rank=0]
Info: compute_column_integ took 10.77 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 z = []
287 h = []
288 vz = []
289 az = []
290
291 delta_r = 0.01
292
293 for i in range(len(dat["xyz"])):
294 r = (dat["xyz"][i][0] ** 2 + dat["xyz"][i][1] ** 2) ** 0.5
295 if r < rcenter + delta_r and r > rcenter - delta_r:
296 z.append(dat["xyz"][i][2])
297 h.append(dat["hpart"][i])
298 vz.append(dat["vxyz"][i][2])
299 az.append(dat["axyz"][i][2])
300
301 rho = pmass * (model.get_hfact() / np.array(h)) ** 3
302
303 fig, axs = plt.subplots(nrows=3, ncols=1, sharex=True)
304
305 from scipy.optimize import curve_fit
306
307 def func(x, a, c):
308 return a * np.exp(-((x / c) ** 2) / 2)
309
310 rho_0 = 0.001
311 p0 = [rho_0, H_profile(rcenter)] # a, b, c
312 popt, pcov = curve_fit(func, z, rho, p0=p0)
313
314 z_ana = np.linspace(-5.0 * H_profile(rcenter), 5.0 * H_profile(rcenter), 100)
315 rho_fit = func(z_ana, *popt)
316
317 axs[0].scatter(z, rho, label="rho")
318
319 axs[0].plot(z_ana, rho_fit, c="black", label="gaussian fit")
320 stddev = abs(popt[1])
321 axs[0].annotate(
322 f"Stddev: {stddev:.5f}",
323 xy=(0.05, 0.95),
324 xycoords="axes fraction",
325 fontsize=10,
326 verticalalignment="top",
327 bbox=dict(boxstyle="round", fc="w"),
328 )
329
330 axs[0].set_ylabel("rho")
331 axs[0].legend()
332
333 axs[1].scatter(z, vz, label="vz")
334
335 vz_fit = np.polyfit(z, vz, 1)
336 vz_fit_fn = np.poly1d(vz_fit)
337 axs[1].plot(z_ana, vz_fit_fn(z_ana), c="red", label="linear fit")
338
339 axs[1].set_ylabel("vz")
340 axs[1].legend()
341
342 axs[2].scatter(z, az, label="az")
343
344 az_fit = np.polyfit(z, az, 1)
345 az_fit_fn = np.poly1d(az_fit)
346 print(f"r={rcenter} az_fit={az_fit}")
347 axs[2].plot(z_ana, az_fit_fn(z_ana), c="red", label="linear fit")
348
349 axs[2].set_ylabel("az")
350 axs[2].set_xlabel("z")
351 axs[2].legend()
352
353 plt.tight_layout()
354 plt.show()
Info: collected : 1 patches [PatchScheduler][rank=0]
r=1.0 az_fit=[1.23422503 0.02427204]
r=2.0 az_fit=[0.15016639 0.03196043]
r=3.0 az_fit=[ 0.09963025 -0.0031214 ]
Total running time of the script: (2 minutes 9.423 seconds)
Estimated memory usage: 700 MB