Note
Go to the end to download the full example code.
Using Shamrock SPH rendering module#
This example demonstrates how to use the Shamrock SPH rendering module to render the density field or the velocity field of a SPH simulation.
The test simulation to showcase the rendering module
13 import glob
14 import json
15 import os # for makedirs
16
17 import matplotlib
18 import matplotlib.pyplot as plt
19 import numpy as np
20
21 import shamrock
22
23 # If we use the shamrock executable to run this script instead of the python interpreter,
24 # we should not initialize the system as the shamrock executable needs to handle specific MPI logic
25 if not shamrock.sys.is_initialized():
26 shamrock.change_loglevel(1)
27 shamrock.sys.init("0:0")
Setup units
33 si = shamrock.UnitSystem()
34 sicte = shamrock.Constants(si)
35 codeu = shamrock.UnitSystem(
36 unit_time=sicte.second(),
37 unit_length=sicte.au(),
38 unit_mass=sicte.sol_mass(),
39 )
40 ucte = shamrock.Constants(codeu)
41 G = ucte.G()
42 c = ucte.c()
List parameters
47 # Resolution
48 Npart = 100000
49
50 # Domain decomposition parameters
51 scheduler_split_val = int(1.0e7) # split patches with more than 1e7 particles
52 scheduler_merge_val = scheduler_split_val // 16
53
54 # Disc parameter
55 center_mass = 1e6 # [sol mass]
56 disc_mass = 0.001 # [sol mass]
57 Rg = G * center_mass / (c * c) # [au]
58 rin = 4.0 * Rg # [au]
59 rout = 10 * rin # [au]
60 r0 = rin # [au]
61
62 H_r_0 = 0.01
63 q = 0.75
64 p = 3.0 / 2.0
65
66 Tin = 2 * np.pi * np.sqrt(rin * rin * rin / (G * center_mass))
67 if shamrock.sys.world_rank() == 0:
68 print(" Orbital period : ", Tin, " [seconds]")
69
70 # Sink parameters
71 center_racc = rin / 2.0 # [au]
72 inclination = 30.0 * np.pi / 180.0
73
74
75 # Viscosity parameter
76 alpha_AV = 1.0e-3 / 0.08
77 alpha_u = 1.0
78 beta_AV = 2.0
79
80 # Integrator parameters
81 C_cour = 0.3
82 C_force = 0.25
83
84
85 # Disc profiles
86 def sigma_profile(r):
87 sigma_0 = 1.0 # We do not care as it will be renormalized
88 return sigma_0 * (r / r0) ** (-p)
89
90
91 def kep_profile(r):
92 return (G * center_mass / r) ** 0.5
93
94
95 def omega_k(r):
96 return kep_profile(r) / r
97
98
99 def cs_profile(r):
100 cs_in = (H_r_0 * r0) * omega_k(r0)
101 return ((r / r0) ** (-q)) * cs_in
Orbital period : 247.58972132551145 [seconds]
Utility functions and quantities deduced from the base one
107 # Deduced quantities
108 pmass = disc_mass / Npart
109
110 bsize = rout * 2
111 bmin = (-bsize, -bsize, -bsize)
112 bmax = (bsize, bsize, bsize)
113
114 cs0 = cs_profile(r0)
115
116
117 def rot_profile(r):
118 return ((kep_profile(r) ** 2) - (2 * p + q) * cs_profile(r) ** 2) ** 0.5
119
120
121 def H_profile(r):
122 H = cs_profile(r) / omega_k(r)
123 # fact = (2.**0.5) * 3. # factor taken from phantom, to fasten thermalizing
124 fact = 1.0
125 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)
133 ctx = shamrock.Context()
134 ctx.pdata_layout_new()
Attach a SPH model to the context
139 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
140
141 # Generate the default config
142 cfg = model.gen_default_config()
143 cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV)
144 cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q, r0=r0)
145
146 # cfg.add_ext_force_point_mass(center_mass, center_racc)
147
148 cfg.add_kill_sphere(center=(0, 0, 0), radius=bsize) # kill particles outside the simulation box
149 cfg.add_ext_force_lense_thirring(
150 central_mass=center_mass,
151 Racc=rin,
152 a_spin=0.9,
153 dir_spin=(np.sin(inclination), np.cos(inclination), 0.0),
154 )
155
156 cfg.set_units(codeu)
157 cfg.set_particle_mass(pmass)
158 # Set the CFL
159 cfg.set_cfl_cour(C_cour)
160 cfg.set_cfl_force(C_force)
161
162 # On a chaotic disc, we disable to two stage search to avoid giant leaves
163 cfg.set_tree_reduction_level(6)
164 cfg.set_two_stage_search(False)
165
166 # Enable this to debug the neighbor counts
167 # cfg.set_show_neigh_stats(True)
168
169 # Standard way to set the smoothing length (e.g. Price et al. 2018)
170 cfg.set_smoothing_length_density_based()
171
172 # Standard density based smoothing lenght but with a neighbor count limit
173 # Use it if you have large slowdowns due to giant particles
174 # I recommend to use it if you have a circumbinary discs as the issue is very likely to happen
175 # cfg.set_smoothing_length_density_based_neigh_lim(500)
176
177 # Set the solver config to be the one stored in cfg
178 model.set_solver_config(cfg)
179
180 # Print the solver config
181 model.get_current_config().print_status()
182
183 # Init the scheduler & fields
184 model.init_scheduler(scheduler_split_val, scheduler_merge_val)
185
186 # Set the simulation box size
187 model.resize_simulation_box(bmin, bmax)
188
189 # Create the setup
190
191 setup = model.get_setup()
192 gen_disc = setup.make_generator_disc_mc(
193 part_mass=pmass,
194 disc_mass=disc_mass,
195 r_in=rin,
196 r_out=rout,
197 sigma_profile=sigma_profile,
198 H_profile=H_profile,
199 rot_profile=rot_profile,
200 cs_profile=cs_profile,
201 random_seed=666,
202 init_h_factor=0.06,
203 )
204
205 # Print the dot graph of the setup
206 print(gen_disc.get_dot())
207
208 # Apply the setup
209 setup.apply_setup(gen_disc)
210
211 model.change_htolerances(coarse=1.3, fine=1.1)
212 model.timestep()
213 model.change_htolerances(coarse=1.1, fine=1.1)
214
215 for i in range(5):
216 model.timestep()
----- 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": 1.0019944020500018e-05,
"eos_type": "locally_isothermal_lp07",
"q": 0.75,
"r0": 0.03948371767577914
},
"epsilon_h": 1e-06,
"ext_force_config": {
"force_list": [
{
"Racc": 0.03948371767577914,
"a_spin": 0.9,
"central_mass": 1000000.0,
"dir_spin": [
0.49999999999999994,
0.8660254037844387,
0.0
],
"force_type": "lense_thirring"
}
]
},
"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": [
{
"center": [
0.0,
0.0,
0.0
],
"radius": 0.7896743535155828,
"type": "sphere"
}
],
"particle_reordering_step_freq": 1000,
"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": 6,
"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": 1.0
},
"use_two_stage_search": false
}
]
------------------------------------
digraph G {
rankdir=LR;
node_0 [label="GeneratorMCDisc"];
node_2 [label="Simulation"];
node_0 -> node_2;
}
Info: pushing data in scheduler, N = 100000 [DataInserterUtility][rank=0]
Info: reattributing data ... [DataInserterUtility][rank=0]
Info: reattributing data done in 15.00 ms [DataInserterUtility][rank=0]
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 : 13.37 us (69.6%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1603.00 ns (0.3%)
patch tree reduce : 1263.00 ns (0.2%)
gen split merge : 741.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 892.00 ns (0.2%)
LB compute : 570.80 us (98.0%)
LB move op cnt : 0
LB apply : 3.23 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.65 us (63.8%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1753.00 ns (0.4%)
patch tree reduce : 431.00 ns (0.1%)
gen split merge : 391.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 290.00 ns (0.1%)
LB compute : 416.63 us (98.1%)
LB move op cnt : 0
LB apply : 1983.00 ns (0.5%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.50 us (65.7%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1653.00 ns (0.3%)
patch tree reduce : 461.00 ns (0.1%)
gen split merge : 390.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 270.00 ns (0.1%)
LB compute : 481.12 us (98.3%)
LB move op cnt : 0
LB apply : 1953.00 ns (0.4%)
Info: --------------------------------------------- [DataInserterUtility][rank=0]
Info: Final load balancing step 0 of 3 [SPH setup][rank=0]
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 : 2.75 us (63.1%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1263.00 ns (0.3%)
patch tree reduce : 361.00 ns (0.1%)
gen split merge : 391.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 260.00 ns (0.1%)
LB compute : 444.19 us (98.4%)
LB move op cnt : 0
LB apply : 1904.00 ns (0.4%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.60 us (66.4%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1572.00 ns (0.4%)
patch tree reduce : 441.00 ns (0.1%)
gen split merge : 391.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 300.00 ns (0.1%)
LB compute : 384.43 us (98.0%)
LB move op cnt : 0
LB apply : 1993.00 ns (0.5%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.60 us (68.5%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1212.00 ns (0.3%)
patch tree reduce : 401.00 ns (0.1%)
gen split merge : 370.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 250.00 ns (0.1%)
LB compute : 369.97 us (98.1%)
LB move op cnt : 0
LB apply : 1894.00 ns (0.5%)
Info: --------------------------------------------- [DataInserterUtility][rank=0]
Info: Final load balancing step 1 of 3 [SPH setup][rank=0]
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 : 2.65 us (61.1%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1623.00 ns (0.4%)
patch tree reduce : 421.00 ns (0.1%)
gen split merge : 381.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 271.00 ns (0.1%)
LB compute : 432.94 us (98.2%)
LB move op cnt : 0
LB apply : 2.03 us (0.5%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.75 us (64.6%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1612.00 ns (0.4%)
patch tree reduce : 431.00 ns (0.1%)
gen split merge : 371.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 260.00 ns (0.1%)
LB compute : 383.70 us (98.0%)
LB move op cnt : 0
LB apply : 2.02 us (0.5%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.62 us (67.3%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1252.00 ns (0.3%)
patch tree reduce : 401.00 ns (0.1%)
gen split merge : 360.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 260.00 ns (0.1%)
LB compute : 422.54 us (98.2%)
LB move op cnt : 0
LB apply : 2.11 us (0.5%)
Info: --------------------------------------------- [DataInserterUtility][rank=0]
Info: Final load balancing step 2 of 3 [SPH setup][rank=0]
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 : 2.62 us (65.3%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1363.00 ns (0.3%)
patch tree reduce : 401.00 ns (0.1%)
gen split merge : 391.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 260.00 ns (0.1%)
LB compute : 396.96 us (98.1%)
LB move op cnt : 0
LB apply : 2.03 us (0.5%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.73 us (67.5%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1553.00 ns (0.4%)
patch tree reduce : 401.00 ns (0.1%)
gen split merge : 371.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 260.00 ns (0.1%)
LB compute : 382.24 us (98.0%)
LB move op cnt : 0
LB apply : 2.00 us (0.5%)
Info: Compute load ... [DataInserterUtility][rank=0]
Info: run scheduler step ... [DataInserterUtility][rank=0]
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.75 us (64.2%)
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 1603.00 ns (0.4%)
patch tree reduce : 371.00 ns (0.1%)
gen split merge : 391.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 271.00 ns (0.1%)
LB compute : 386.52 us (98.1%)
LB move op cnt : 0
LB apply : 1803.00 ns (0.5%)
Info: --------------------------------------------- [DataInserterUtility][rank=0]
Info: the setup took : 0.345387108 s [SPH setup][rank=0]
---------------- t = 0, dt = 0 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 6.60 us (1.0%)
patch tree reduce : 752.00 ns (0.1%)
gen split merge : 531.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 721.00 ns (0.1%)
LB compute : 653.94 us (97.4%)
LB move op cnt : 0
LB apply : 2.83 us (0.4%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.48 us (68.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.009009940708617091 unconverged cnt = 99991
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.00916897508625909 unconverged cnt = 99969
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.00916897508625909 unconverged cnt = 99937
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.009173539640934479 unconverged cnt = 99800
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.00917353964093448 unconverged cnt = 99408
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.009173539640934482 unconverged cnt = 98057
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.010149085311974134 unconverged cnt = 92700
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.010149085311974134 unconverged cnt = 65335
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.010149085311974134 unconverged cnt = 11471
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.010149085311974134 unconverged cnt = 138
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.567558866379785e-10,2.204738873133719e-09,0)
sum a = (-2.4852304073329302e-11,-5.2417203944085925e-12,5.679286915979897e-12)
sum e = 1.2874775315064656e-10
sum de = -8.517665920295599e-31
Info: CFL hydro = 0.015241645279522596 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.015241645279522596 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) | MPI | alloc | mem (max) |
---------------------------------------------------------------------------------------
| 0 | 4.3934e+04 | 100000 | 2.276e+00 | 0 % | 0 % | 1.16 GB |
---------------------------------------------------------------------------------------
Info: estimated rate : 0 (tsim/hr) [sph::Model][rank=0]
---------------- t = 0, dt = 0.015241645279522596 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.66 us (1.5%)
patch tree reduce : 1322.00 ns (0.3%)
gen split merge : 681.00 ns (0.1%)
split / merge op : 0/0
apply split merge : 1232.00 ns (0.2%)
LB compute : 480.39 us (95.7%)
LB move op cnt : 0
LB apply : 3.76 us (0.7%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.60 us (68.3%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.563770966349765e-10,2.204658980690752e-09,8.656167661400032e-14)
sum a = (-2.4850517878382777e-11,-5.249044853656845e-12,5.681115368334097e-12)
sum e = 1.2874775609029635e-10
sum de = 1.815434405185107e-19
Info: CFL hydro = 0.518214322205733 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.518214322205733 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) | MPI | alloc | mem (max) |
---------------------------------------------------------------------------------------
| 0 | 2.7490e+05 | 100000 | 3.638e-01 | 0 % | 0 % | 1.16 GB |
---------------------------------------------------------------------------------------
Info: estimated rate : 150.83720071152743 (tsim/hr) [sph::Model][rank=0]
---------------- t = 0.015241645279522596, dt = 0.518214322205733 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.33 us (1.6%)
patch tree reduce : 1362.00 ns (0.3%)
gen split merge : 732.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1042.00 ns (0.2%)
LB compute : 440.22 us (95.5%)
LB move op cnt : 0
LB apply : 3.39 us (0.7%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.48 us (71.4%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.4349921596840273e-10,2.201938794651362e-09,3.0306109608989914e-12)
sum a = (-2.478842612289725e-11,-5.4978614697350675e-12,5.742963070288513e-12)
sum e = 1.287511514591596e-10
sum de = 5.913793218247977e-18
Info: CFL hydro = 0.8534471191106839 sink sink = inf [SPH][rank=0]
Info: cfl dt = 0.8534471191106839 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) | MPI | alloc | mem (max) |
---------------------------------------------------------------------------------------
| 0 | 2.6729e+05 | 100000 | 3.741e-01 | 0 % | 0 % | 1.16 GB |
---------------------------------------------------------------------------------------
Info: estimated rate : 4986.558579207298 (tsim/hr) [sph::Model][rank=0]
---------------- t = 0.5334559674852556, dt = 0.8534471191106839 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 6.66 us (1.5%)
patch tree reduce : 1192.00 ns (0.3%)
gen split merge : 692.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1072.00 ns (0.2%)
LB compute : 432.85 us (95.5%)
LB move op cnt : 0
LB apply : 3.82 us (0.8%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.52 us (71.4%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (2.2235969352497844e-10,2.1971821904516175e-09,7.947951430869919e-12)
sum a = (-2.4680578479084503e-11,-5.90683436974251e-12,5.842099589323029e-12)
sum e = 1.2875697699044605e-10
sum de = 1.4336374390089268e-17
Info: CFL hydro = 1.0764417673903492 sink sink = inf [SPH][rank=0]
Info: cfl dt = 1.0764417673903492 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) | MPI | alloc | mem (max) |
---------------------------------------------------------------------------------------
| 0 | 2.6608e+05 | 100000 | 3.758e-01 | 0 % | 0 % | 1.16 GB |
---------------------------------------------------------------------------------------
Info: estimated rate : 8174.963022266753 (tsim/hr) [sph::Model][rank=0]
---------------- t = 1.3869030865959395, dt = 1.0764417673903492 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 7.15 us (1.8%)
patch tree reduce : 1092.00 ns (0.3%)
gen split merge : 852.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1022.00 ns (0.3%)
LB compute : 372.41 us (94.8%)
LB move op cnt : 0
LB apply : 3.56 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.40 us (70.9%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.9583850913728363e-10,2.1906493088514416e-09,1.427893532635529e-11)
sum a = (-2.4534717930977417e-11,-6.421102616967876e-12,5.961980249192321e-12)
sum e = 1.2876244368328456e-10
sum de = 2.4458302969399853e-17
Info: CFL hydro = 1.2233723099251912 sink sink = inf [SPH][rank=0]
Info: cfl dt = 1.2233723099251912 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) | MPI | alloc | mem (max) |
---------------------------------------------------------------------------------------
| 0 | 2.6502e+05 | 100000 | 3.773e-01 | 0 % | 0 % | 1.16 GB |
---------------------------------------------------------------------------------------
Info: estimated rate : 10269.841260634357 (tsim/hr) [sph::Model][rank=0]
---------------- t = 2.463344853986289, dt = 1.2233723099251912 ----------------
Info: summary : [LoadBalance][rank=0]
Info: - strategy "psweep" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: - strategy "round robin" : max = 100000 min = 100000 [LoadBalance][rank=0]
Info: Loadbalance stats : [LoadBalance][rank=0]
npatch = 1
min = 100000
max = 100000
avg = 100000
efficiency = 100.00%
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 6.09 us (1.6%)
patch tree reduce : 1172.00 ns (0.3%)
gen split merge : 862.00 ns (0.2%)
split / merge op : 0/0
apply split merge : 1032.00 ns (0.3%)
LB compute : 362.04 us (94.9%)
LB move op cnt : 0
LB apply : 3.35 us (0.9%)
Info: Scheduler step timings : [Scheduler][rank=0]
metadata sync : 2.69 us (69.5%)
Info: free boundaries skipping geometry update [PositionUpdated][rank=0]
Info: conservation infos : [sph::Model][rank=0]
sum v = (1.6590191978177015e-10,2.182517119800184e-09,2.163717915023062e-11)
sum a = (-2.4355727946293835e-11,-7.003154772316315e-12,6.091129908666387e-12)
sum e = 1.2876675914297396e-10
sum de = 3.6037595924070245e-17
Info: CFL hydro = 1.321630747799398 sink sink = inf [SPH][rank=0]
Info: cfl dt = 1.321630747799398 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) | MPI | alloc | mem (max) |
---------------------------------------------------------------------------------------
| 0 | 2.6436e+05 | 100000 | 3.783e-01 | 0 % | 0 % | 1.16 GB |
---------------------------------------------------------------------------------------
Info: estimated rate : 11642.670481090108 (tsim/hr) [sph::Model][rank=0]
Usual cartesian rendering
222 ext = rout * 1.5
223 center = (0.0, 0.0, 0.0)
224 delta_x = (ext * 2, 0, 0.0)
225 delta_y = (0.0, ext * 2, 0.0)
226 nx = 1024
227 ny = 1024
228 nr = 1024
229 ntheta = 1024
230
231 arr_rho = model.render_cartesian_column_integ(
232 "rho",
233 "f64",
234 center=center,
235 delta_x=delta_x,
236 delta_y=delta_y,
237 nx=nx,
238 ny=ny,
239 )
240
241 arr_vxyz = model.render_cartesian_column_integ(
242 "vxyz",
243 "f64_3",
244 center=center,
245 delta_x=delta_x,
246 delta_y=delta_y,
247 nx=nx,
248 ny=ny,
249 )
250
251
252 def plot_rho_integ(metadata, arr_rho):
253
254 ext = metadata["extent"]
255
256 my_cmap = matplotlib.colormaps["gist_heat"].copy() # copy the default cmap
257 my_cmap.set_bad(color="black")
258
259 res = plt.imshow(
260 arr_rho, cmap=my_cmap, origin="lower", extent=ext, norm="log", vmin=1e-6, vmax=1e-2
261 )
262
263 plt.xlabel("x")
264 plt.ylabel("y")
265 plt.title(f"t = {metadata['time']:0.3f} [seconds]")
266
267 cbar = plt.colorbar(res, extend="both")
268 cbar.set_label(r"$\int \rho \, \mathrm{d}z$ [code unit]")
269
270
271 def plot_vz_integ(metadata, arr_vz):
272 ext = metadata["extent"]
273
274 # if you want an adaptive colorbar
275 v_ext = np.max(arr_vz)
276 v_ext = max(v_ext, np.abs(np.min(arr_vz)))
277 # v_ext = 1e-6
278
279 res = plt.imshow(arr_vz, cmap="seismic", origin="lower", extent=ext, vmin=-v_ext, vmax=v_ext)
280 plt.xlabel("x")
281 plt.ylabel("y")
282 plt.title(f"t = {metadata['time']:0.3f} [seconds]")
283
284 cbar = plt.colorbar(res, extend="both")
285 cbar.set_label(r"$\int v_z \, \mathrm{d}z$ [code unit]")
286
287
288 metadata = {"extent": [-ext, ext, -ext, ext], "time": model.get_time()}
289
290 dpi = 200
291
292 plt.figure(dpi=dpi)
293 plot_rho_integ(metadata, arr_rho)
294
295 plt.figure(dpi=dpi)
296 plot_vz_integ(metadata, arr_vxyz[:, :, 2])
![t = 3.687 [seconds]](../../_images/sphx_glr_run_sph_rendering_001.png)
![t = 3.687 [seconds]](../../_images/sphx_glr_run_sph_rendering_002.png)
Info: compute_column_integ field_name: rho, rays count: 1048576 [sph::CartesianRender][rank=0]
Info: compute_column_integ took 661.61 ms [sph::CartesianRender][rank=0]
Info: compute_column_integ field_name: vxyz, rays count: 1048576 [sph::CartesianRender][rank=0]
Info: compute_column_integ took 646.94 ms [sph::CartesianRender][rank=0]
Cylindrical rendering
302 def make_cylindrical_coords(nr, ntheta):
303 """
304 Generate a list of positions in cylindrical coordinates (r, theta)
305 spanning [0, ext*2] x [-pi, pi] for use with the rendering module.
306
307 Returns:
308 list: List of [x, y, z] coordinate lists
309 """
310
311 # Create the cylindrical coordinate grid
312 r_vals = np.linspace(0, ext, nr)
313 theta_vals = np.linspace(-np.pi, np.pi, ntheta)
314
315 # Create meshgrid
316 r_grid, theta_grid = np.meshgrid(r_vals, theta_vals)
317
318 # Convert to Cartesian coordinates (z = 0 for a disc in the xy-plane)
319 x_grid = r_grid * np.cos(theta_grid)
320 y_grid = r_grid * np.sin(theta_grid)
321 z_grid = np.zeros_like(r_grid)
322
323 # Flatten and stack to create list of positions
324 positions = np.column_stack([x_grid.ravel(), y_grid.ravel(), z_grid.ravel()])
325
326 return [tuple(pos) for pos in positions]
327
328
329 def positions_to_rays(positions):
330 return [shamrock.math.Ray_f64_3(tuple(position), (0.0, 0.0, 1.0)) for position in positions]
331
332
333 positions_cylindrical = make_cylindrical_coords(nr, ntheta)
334 rays_cylindrical = positions_to_rays(positions_cylindrical)
335
336
337 arr_rho_cylindrical = model.render_column_integ("rho", "f64", rays_cylindrical)
338
339 arr_rho_pos = model.render_slice("rho", "f64", positions_cylindrical)
340
341
342 def plot_rho_integ_cylindrical(metadata, arr_rho_cylindrical):
343 ext = metadata["extent"]
344
345 my_cmap = matplotlib.colormaps["gist_heat"].copy() # copy the default cmap
346 my_cmap.set_bad(color="black")
347
348 arr_rho_cylindrical = np.array(arr_rho_cylindrical).reshape(nr, ntheta)
349
350 res = plt.imshow(
351 arr_rho_cylindrical,
352 cmap=my_cmap,
353 origin="lower",
354 extent=ext,
355 norm="log",
356 vmin=1e-6,
357 vmax=1e-2,
358 aspect="auto",
359 )
360 plt.xlabel("r")
361 plt.ylabel(r"$\theta$")
362 plt.title(f"t = {metadata['time']:0.3f} [seconds]")
363 cbar = plt.colorbar(res, extend="both")
364 cbar.set_label(r"$\int \rho \, \mathrm{d}z$ [code unit]")
365
366
367 def plot_rho_slice_cylindrical(metadata, arr_rho_pos):
368 ext = metadata["extent"]
369
370 my_cmap = matplotlib.colormaps["gist_heat"].copy() # copy the default cmap
371 my_cmap.set_bad(color="black")
372
373 arr_rho_pos = np.array(arr_rho_pos).reshape(nr, ntheta)
374
375 res = plt.imshow(
376 arr_rho_pos, cmap=my_cmap, origin="lower", extent=ext, norm="log", vmin=1e-8, aspect="auto"
377 )
378 plt.xlabel("r")
379 plt.ylabel(r"$\theta$")
380 plt.title(f"t = {metadata['time']:0.3f} [seconds]")
381 cbar = plt.colorbar(res, extend="both")
382 cbar.set_label(r"$\rho$ [code unit]")
383
384
385 metadata = {"extent": [0, ext, -np.pi, np.pi], "time": model.get_time()}
386
387 plt.figure(dpi=dpi)
388 plot_rho_integ_cylindrical(metadata, arr_rho_cylindrical)
389
390 plt.figure(dpi=dpi)
391 plot_rho_slice_cylindrical(metadata, arr_rho_pos)
392
393 plt.show()
![t = 3.687 [seconds]](../../_images/sphx_glr_run_sph_rendering_003.png)
![t = 3.687 [seconds]](../../_images/sphx_glr_run_sph_rendering_004.png)
Info: compute_column_integ field_name: rho, rays count: 1048576 [sph::CartesianRender][rank=0]
Info: compute_column_integ took 833.39 ms [sph::CartesianRender][rank=0]
Info: compute_slice field_name: rho, positions count: 1048576 [sph::CartesianRender][rank=0]
Info: compute_slice took 445.10 ms [sph::CartesianRender][rank=0]
Total running time of the script: (0 minutes 13.206 seconds)
Estimated memory usage: 620 MB