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")

Use shamrock documentation style for matplotlib

32 shamrock.matplotlib.set_shamrock_mpl_style()

Setup units

38 si = shamrock.UnitSystem()
39 sicte = shamrock.Constants(si)
40 codeu = shamrock.UnitSystem(
41     unit_time=sicte.second(),
42     unit_length=sicte.au(),
43     unit_mass=sicte.sol_mass(),
44 )
45 ucte = shamrock.Constants(codeu)
46 G = ucte.G()
47 c = ucte.c()

List parameters

 52 # Resolution
 53 Npart = 100000
 54
 55 # Domain decomposition parameters
 56 scheduler_split_val = int(1.0e7)  # split patches with more than 1e7 particles
 57 scheduler_merge_val = scheduler_split_val // 16
 58
 59 # Disc parameter
 60 center_mass = 1e6  # [sol mass]
 61 disc_mass = 0.001  # [sol mass]
 62 Rg = G * center_mass / (c * c)  # [au]
 63 rin = 4.0 * Rg  # [au]
 64 rout = 10 * rin  # [au]
 65 r0 = rin  # [au]
 66
 67 H_r_0 = 0.05
 68 q = 0.75
 69 p = 3.0 / 2.0
 70
 71 Tin = 2 * np.pi * np.sqrt(rin * rin * rin / (G * center_mass))
 72 if shamrock.sys.world_rank() == 0:
 73     print(" Orbital period : ", Tin, " [seconds]")
 74
 75 # Sink parameters
 76 center_racc = rin / 2.0  # [au]
 77 inclination = 30.0 * np.pi / 180.0
 78
 79
 80 # Viscosity parameter
 81 alpha_AV = 1.0e-3 / 0.08
 82 alpha_u = 1.0
 83 beta_AV = 2.0
 84
 85 # Integrator parameters
 86 C_cour = 0.3
 87 C_force = 0.25
 88
 89
 90 # Disc profiles
 91 def sigma_profile(r):
 92     sigma_0 = 1.0  # We do not care as it will be renormalized
 93     return sigma_0 * (r / r0) ** (-p)
 94
 95
 96 def kep_profile(r):
 97     return (G * center_mass / r) ** 0.5
 98
 99
100 def omega_k(r):
101     return kep_profile(r) / r
102
103
104 def cs_profile(r):
105     cs_in = (H_r_0 * r0) * omega_k(r0)
106     return ((r / r0) ** (-q)) * cs_in

Orbital period :  247.58972132551145  [seconds]

Utility functions and quantities deduced from the base one

112 # Deduced quantities
113 pmass = disc_mass / Npart
114
115 bsize = rout * 2
116 bmin = (-bsize, -bsize, -bsize)
117 bmax = (bsize, bsize, bsize)
118
119 cs0 = cs_profile(r0)
120
121
122 def rot_profile(r):
123     return ((kep_profile(r) ** 2) - (2 * p + q) * cs_profile(r) ** 2) ** 0.5
124
125
126 def H_profile(r):
127     H = cs_profile(r) / omega_k(r)
128     # fact = (2.**0.5) * 3. # factor taken from phantom, to fasten thermalizing
129     fact = 1.0
130     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)

138 ctx = shamrock.Context()
139 ctx.pdata_layout_new()

Attach a SPH model to the context

144 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
145
146 # Generate the default config
147 cfg = model.gen_default_config()
148 cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV)
149 cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q, r0=r0)
150
151 # cfg.add_ext_force_point_mass(center_mass, center_racc)
152
153 cfg.add_kill_sphere(center=(0, 0, 0), radius=bsize)  # kill particles outside the simulation box
154 # cfg.add_ext_force_lense_thirring(
155 #     central_mass=center_mass,
156 #     Racc=rin,
157 #     a_spin=0.9,
158 #     dir_spin=(np.sin(inclination), np.cos(inclination), 0.0),
159 # )
160
161 cfg.set_units(codeu)
162 cfg.set_particle_mass(pmass)
163 # Set the CFL
164 cfg.set_cfl_cour(C_cour)
165 cfg.set_cfl_force(C_force)
166
167 # On a chaotic disc, we disable to two stage search to avoid giant leaves
168 cfg.set_tree_reduction_level(6)
169 cfg.set_two_stage_search(False)
170
171 # Enable this to debug the neighbor counts
172 # cfg.set_show_neigh_stats(True)
173
174 # Standard way to set the smoothing length (e.g. Price et al. 2018)
175 cfg.set_smoothing_length_density_based()
176
177 # Standard density based smoothing length but with a neighbor count limit
178 # Use it if you have large slowdowns due to giant particles
179 # I recommend to use it if you have a circumbinary discs as the issue is very likely to happen
180 # cfg.set_smoothing_length_density_based_neigh_lim(500)
181
182 cfg.set_scheduler_config(split_load_value=scheduler_split_val, merge_load_value=scheduler_merge_val)
183
184 # Set the solver config to be the one stored in cfg
185 model.set_solver_config(cfg)
186
187 # Print the solver config
188 model.get_current_config().print_status()
189
190 # Init the scheduler & fields
191 model.init()
192
193 # Set the simulation box size
194 model.resize_simulation_box(bmin, bmax)
195
196 # Create the setup
197
198 setup = model.get_setup()
199 gen_disc = setup.make_generator_disc_mc(
200     part_mass=pmass,
201     disc_mass=disc_mass,
202     r_in=rin,
203     r_out=rout,
204     sigma_profile=sigma_profile,
205     H_profile=H_profile,
206     rot_profile=rot_profile,
207     cs_profile=cs_profile,
208     random_seed=666,
209     init_h_factor=0.03,
210 )
211
212 # Print the dot graph of the setup
213 print(gen_disc.get_dot())
214
215 # Apply the setup
216 setup.apply_setup(gen_disc)
217
218 model.do_vtk_dump("init_disc.vtk", True)
219
220 model.change_htolerances(coarse=1.3, fine=1.1)
221 model.timestep()
222 model.change_htolerances(coarse=1.1, fine=1.1)
223
224 for i in range(5):
225     model.timestep()

----- SPH Solver configuration -----
[
    {
        "artif_viscosity": {
            "alpha_AV": 0.0125,
            "alpha_u": 1.0,
            "beta_AV": 2.0,
            "type": "constant_disc"
        },
        "boundary_config": {
            "bc_type": "free"
        },
        "cfl_config": {
            "cfl_cour": 0.3,
            "cfl_force": 0.25,
            "cfl_multiplier_stiffness": 2.0,
            "eta_sink": 0.05
        },
        "combined_dtdiv_divcurlv_compute": false,
        "debug_dump_filename": "",
        "do_debug_dump": false,
        "dust_config": {
            "drag_mode": {
                "type": "none"
            },
            "mode": {
                "type": "none"
            }
        },
        "enable_particle_reordering": false,
        "eos_config": {
            "Tvec": "f64_3",
            "cs0": 5.009972010250009e-05,
            "eos_type": "locally_isothermal_lp07",
            "q": 0.75,
            "r0": 0.03948371767577914
        },
        "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": [
            {
                "center": [
                    0.0,
                    0.0,
                    0.0
                ],
                "radius": 0.7896743535155828,
                "type": "sphere"
            }
        ],
        "particle_reordering_step_freq": 1000,
        "save_dt_to_fields": false,
        "scheduler_config": {
            "merge_load_value": 625000,
            "split_load_value": 10000000
        },
        "self_grav_config": {
            "softening_length": 1e-09,
            "softening_mode": "plummer",
            "type": "none"
        },
        "show_cfl_detail": false,
        "show_ghost_zone_graph": false,
        "show_neigh_stats": false,
        "smoothing_length_config": {
            "type": "density_based"
        },
        "time_state": {
            "cfl_multiplier": 0.01,
            "dt_sph": 0.0,
            "time": 0.0
        },
        "tree_reduction_level": 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
    }
]
------------------------------------
Warning: make_generator_disc_mc: with the current EOS, cs_profile is ignored     [SPHSetup][rank=0]
digraph G {
rankdir=LR;
node_0 [label="GeneratorMCDisc"];
node_2 [label="Simulation"];
node_0 -> node_2;
}

SPH setup: generating particles ...
SPH setup: Nstep = 100000 ( 1.0e+05 ) Ntotal = 100000 ( 1.0e+05 rank min = 3.4e+05 max = 1.0e+05) rate = 1.000000e+05 N.s^-1
SPH setup: the generation step took : 0.305610054 s
SPH setup: final particle count = 100000 beginning injection ...
Info: ---------------------------------------------                   [DataInserterUtility][rank=0]
Info: Compute load ...                                                [DataInserterUtility][rank=0]
Info: run scheduler step ...                                          [DataInserterUtility][rank=0]
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 18.30 us   (81.0%)
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
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     : 921.00 ns  (0.2%)
   patch tree reduce : 852.00 ns  (0.2%)
   gen split merge   : 931.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 982.00 ns  (0.2%)
   LB compute        : 407.78 us  (97.5%)
   LB move op cnt    : 0
   LB apply          : 4.04 us    (1.0%)
Info: patch count stable after 1 runs npatch = 1                      [DataInserterUtility][rank=0]
Info: ---------------------------------------------                   [DataInserterUtility][rank=0]
SPH setup: injected       100000 / 100000 => 100.0% | ranks with patchs = 1 / 1  <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.012666371000000001 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.8% 0.0% |     1.24 GB |     5.29 MB |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 0.348118641 s
Info: dump to init_disc.vtk                                                      [VTK Dump][rank=0]
              - took 14.73 ms, bandwidth = 380.25 MB/s
---------------- t = 0, dt = 0 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
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.04 us    (2.2%)
   patch tree reduce : 1.67 us    (0.5%)
   gen split merge   : 1.07 us    (0.3%)
   split / merge op  : 0/0
   apply split merge : 1.45 us    (0.4%)
   LB compute        : 303.67 us  (93.2%)
   LB move op cnt    : 0
   LB apply          : 4.35 us    (1.3%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 1.78 us    (73.6%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.006569301129452108 unconverged cnt = 100000
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.008540091468287742 unconverged cnt = 100000
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.011102118908774064 unconverged cnt = 100000
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.014432754581406283 unconverged cnt = 99999
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.018762580955828168 unconverged cnt = 99999
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02439135524257662 unconverged cnt = 99995
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.024830299239077123 unconverged cnt = 99988
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.024830299239077126 unconverged cnt = 99974
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 99926
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 99754
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 98899
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 86081
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 16466
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 82
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919894623e-10,3.112830370992019e-11,0)
    sum a = (-2.858736196983264e-29,-3.8539999100070664e-28,8.440683319389103e-27)
    sum e = 1.275752290638463e-10
    sum de = 1.604058902431519e-31
Info: cfl dt = 0.019407933351461047 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.7797e+04 | 100000 |      1 | 3.598e+00 | 0.0% |   0.1% 0.0% |     1.24 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0 (tsim/hr)                                             [sph::Model][rank=0]
---------------- t = 0, dt = 0.019407933351461047 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
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.38 us    (2.2%)
   patch tree reduce : 1.82 us    (0.5%)
   gen split merge   : 1.06 us    (0.3%)
   split / merge op  : 0/0
   apply split merge : 1.32 us    (0.4%)
   LB compute        : 315.36 us  (93.4%)
   LB move op cnt    : 0
   LB apply          : 4.12 us    (1.2%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 1.71 us    (71.2%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919895486e-10,3.112830370990416e-11,-1.64741290796061e-28)
    sum a = (7.845642673942957e-28,4.891615270393584e-28,-1.0590029423046891e-26)
    sum e = 1.2757522909678268e-10
    sum de = 4.1600086798764294e-32
Info: cfl dt = 0.6595177856935768 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.8576e+05 | 100000 |      1 | 5.383e-01 | 0.0% |   0.1% 0.0% |     1.24 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 129.7896706500745 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 0.019407933351461047, dt = 0.6595177856935768 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
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.09 us    (2.1%)
   patch tree reduce : 1.97 us    (0.6%)
   gen split merge   : 1.03 us    (0.3%)
   split / merge op  : 0/0
   apply split merge : 1.04 us    (0.3%)
   LB compute        : 309.92 us  (93.5%)
   LB move op cnt    : 0
   LB apply          : 4.03 us    (1.2%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 1.95 us    (74.4%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919894062e-10,3.112830370980896e-11,-4.27116363652981e-27)
    sum a = (9.68793933422106e-28,-3.3987197008578803e-28,2.031820282226253e-26)
    sum e = 1.2757526707478419e-10
    sum de = 7.538210172953165e-32
Info: cfl dt = 1.067486611711989 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.8176e+05 | 100000 |      1 | 5.502e-01 | 0.0% |   0.0% 0.0% |     1.24 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 4315.397857278035 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 0.6789257190450378, dt = 1.067486611711989 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
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.30 us    (1.9%)
   patch tree reduce : 1.67 us    (0.5%)
   gen split merge   : 941.00 ns  (0.3%)
   split / merge op  : 0/0
   apply split merge : 1.06 us    (0.3%)
   LB compute        : 318.64 us  (93.9%)
   LB move op cnt    : 0
   LB apply          : 4.35 us    (1.3%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 1.66 us    (54.8%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919893851e-10,3.1128303709784325e-11,-1.6606080930920558e-26)
    sum a = (-2.0011153378882847e-28,-4.6480932980579734e-28,-9.740878893424454e-28)
    sum e = 1.275753261081859e-10
    sum de = 3.0838652307946804e-31
Info: cfl dt = 1.3131733244842363 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.7947e+05 | 100000 |      1 | 5.572e-01 | 0.0% |   0.0% 0.0% |     1.24 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 6896.849780228866 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 1.746412330757027, dt = 1.3131733244842363 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
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.74 us    (2.0%)
   patch tree reduce : 1.73 us    (0.5%)
   gen split merge   : 942.00 ns  (0.3%)
   split / merge op  : 0/0
   apply split merge : 1.06 us    (0.3%)
   LB compute        : 314.47 us  (93.9%)
   LB move op cnt    : 0
   LB apply          : 3.76 us    (1.1%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 1.45 us    (68.4%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.29389391989472e-10,3.112830371005654e-11,-2.816808066094176e-26)
    sum a = (7.697411908173454e-28,-1.514071393217062e-28,-1.3169244747436234e-26)
    sum e = 1.2757537036686392e-10
    sum de = 8.482818625364342e-32
Info: cfl dt = 1.4558626079520551 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.7281e+05 | 100000 |      1 | 5.787e-01 | 0.0% |   0.0% 0.0% |     1.24 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 8169.625082769658 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 3.0595856552412632, dt = 1.4558626079520551 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
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.9%)
   patch tree reduce : 1.86 us    (0.5%)
   gen split merge   : 1.20 us    (0.3%)
   split / merge op  : 0/0
   apply split merge : 1.32 us    (0.4%)
   LB compute        : 323.21 us  (93.8%)
   LB move op cnt    : 0
   LB apply          : 4.38 us    (1.3%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.00 us    (70.4%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919895108e-10,3.11283037099426e-11,-2.1607810484457203e-26)
    sum a = (9.730290981583775e-28,-1.1921988732604278e-27,-2.689329607532404e-28)
    sum e = 1.2757539490951678e-10
    sum de = -3.382384372174702e-31
Info: cfl dt = 1.7417916889876892 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.7931e+05 | 100000 |      1 | 5.577e-01 | 0.0% |   0.0% 0.0% |     1.24 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 9397.649443975697 (tsim/hr)                             [sph::Model][rank=0]

Usual cartesian rendering

231 ext = rout * 1.5
232 center = (0.0, 0.0, 0.0)
233 delta_x = (ext * 2, 0, 0.0)
234 delta_y = (0.0, ext * 2, 0.0)
235 nx = 1024
236 ny = 1024
237 nr = 1024
238 ntheta = 1024
239
240 arr_rho = model.render_cartesian_column_integ(
241     "rho",
242     "f64",
243     center=center,
244     delta_x=delta_x,
245     delta_y=delta_y,
246     nx=nx,
247     ny=ny,
248 )
249
250 arr_vxyz = model.render_cartesian_column_integ(
251     "vxyz",
252     "f64_3",
253     center=center,
254     delta_x=delta_x,
255     delta_y=delta_y,
256     nx=nx,
257     ny=ny,
258 )
259
260
261 def plot_rho_integ(metadata, arr_rho):
262     ext = metadata["extent"]
263
264     my_cmap = matplotlib.colormaps["gist_heat"].copy()  # copy the default cmap
265     my_cmap.set_bad(color="black")
266
267     res = plt.imshow(
268         arr_rho, cmap=my_cmap, origin="lower", extent=ext, norm="log", vmin=1e-6, vmax=1e-2
269     )
270
271     plt.xlabel("x")
272     plt.ylabel("y")
273     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
274
275     cbar = plt.colorbar(res, extend="both")
276     cbar.set_label(r"$\int \rho \, \mathrm{d}z$ [code unit]")
277
278
279 def plot_vz_integ(metadata, arr_vz):
280     ext = metadata["extent"]
281
282     # if you want an adaptive colorbar
283     v_ext = np.max(arr_vz)
284     v_ext = max(v_ext, np.abs(np.min(arr_vz)))
285     # v_ext = 1e-6
286
287     res = plt.imshow(arr_vz, cmap="seismic", origin="lower", extent=ext, vmin=-v_ext, vmax=v_ext)
288     plt.xlabel("x")
289     plt.ylabel("y")
290     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
291
292     cbar = plt.colorbar(res, extend="both")
293     cbar.set_label(r"$\int v_z \, \mathrm{d}z$ [code unit]")
294
295
296 metadata = {"extent": [-ext, ext, -ext, ext], "time": model.get_time()}
297
298 dpi = 200
299
300 plt.figure(dpi=dpi)
301 plot_rho_integ(metadata, arr_rho)
302
303 plt.figure(dpi=dpi)
304 plot_vz_integ(metadata, arr_vxyz[:, :, 2])
  • t = 4.515 [seconds]
  • t = 4.515 [seconds]
Info: compute_column_integ field_name: rho, rays count: 1048576      [sph::CartesianRender][rank=0]
Info: compute_column_integ took 1.32 s                               [sph::CartesianRender][rank=0]
Info: compute_column_integ field_name: vxyz, rays count: 1048576     [sph::CartesianRender][rank=0]
Info: compute_column_integ took 1.35 s                               [sph::CartesianRender][rank=0]

Cylindrical rendering

310 def make_cylindrical_coords(nr, ntheta):
311     """
312     Generate a list of positions in cylindrical coordinates (r, theta)
313     spanning [0, ext*2] x [-pi, pi] for use with the rendering module.
314
315     Returns:
316         list: List of [x, y, z] coordinate lists
317     """
318
319     # Create the cylindrical coordinate grid
320     r_vals = np.linspace(0, ext, nr)
321     theta_vals = np.linspace(-np.pi, np.pi, ntheta)
322
323     # Create meshgrid
324     r_grid, theta_grid = np.meshgrid(r_vals, theta_vals)
325
326     # Convert to Cartesian coordinates (z = 0 for a disc in the xy-plane)
327     x_grid = r_grid * np.cos(theta_grid)
328     y_grid = r_grid * np.sin(theta_grid)
329     z_grid = np.zeros_like(r_grid)
330
331     # Flatten and stack to create list of positions
332     positions = np.column_stack([x_grid.ravel(), y_grid.ravel(), z_grid.ravel()])
333
334     return [tuple(pos) for pos in positions]
335
336
337 def positions_to_rays(positions):
338     return [shamrock.math.Ray_f64_3(tuple(position), (0.0, 0.0, 1.0)) for position in positions]
339
340
341 positions_cylindrical = make_cylindrical_coords(nr, ntheta)
342 rays_cylindrical = positions_to_rays(positions_cylindrical)
343
344
345 arr_rho_cylindrical = model.render_column_integ("rho", "f64", rays_cylindrical)
346
347 arr_rho_pos = model.render_slice("rho", "f64", positions_cylindrical)
348
349
350 def plot_rho_integ_cylindrical(metadata, arr_rho_cylindrical):
351     ext = metadata["extent"]
352
353     my_cmap = matplotlib.colormaps["gist_heat"].copy()  # copy the default cmap
354     my_cmap.set_bad(color="black")
355
356     arr_rho_cylindrical = np.array(arr_rho_cylindrical).reshape(nr, ntheta)
357
358     res = plt.imshow(
359         arr_rho_cylindrical,
360         cmap=my_cmap,
361         origin="lower",
362         extent=ext,
363         norm="log",
364         vmin=1e-6,
365         vmax=1e-2,
366         aspect="auto",
367     )
368     plt.xlabel("r")
369     plt.ylabel(r"$\theta$")
370     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
371     cbar = plt.colorbar(res, extend="both")
372     cbar.set_label(r"$\int \rho \, \mathrm{d}z$ [code unit]")
373
374
375 def plot_rho_slice_cylindrical(metadata, arr_rho_pos):
376     ext = metadata["extent"]
377
378     my_cmap = matplotlib.colormaps["gist_heat"].copy()  # copy the default cmap
379     my_cmap.set_bad(color="black")
380
381     arr_rho_pos = np.array(arr_rho_pos).reshape(nr, ntheta)
382
383     res = plt.imshow(
384         arr_rho_pos, cmap=my_cmap, origin="lower", extent=ext, norm="log", vmin=1e-8, aspect="auto"
385     )
386     plt.xlabel("r")
387     plt.ylabel(r"$\theta$")
388     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
389     cbar = plt.colorbar(res, extend="both")
390     cbar.set_label(r"$\rho$ [code unit]")
391
392
393 metadata = {"extent": [0, ext, -np.pi, np.pi], "time": model.get_time()}
394
395 plt.figure(dpi=dpi)
396 plot_rho_integ_cylindrical(metadata, arr_rho_cylindrical)
397
398 plt.figure(dpi=dpi)
399 plot_rho_slice_cylindrical(metadata, arr_rho_pos)
400
401 plt.show()
  • t = 4.515 [seconds]
  • t = 4.515 [seconds]
Info: compute_column_integ field_name: rho, rays count: 1048576      [sph::CartesianRender][rank=0]
Info: compute_column_integ took 1.71 s                               [sph::CartesianRender][rank=0]
Info: compute_slice field_name: rho, positions count: 1048576        [sph::CartesianRender][rank=0]
Info: compute_slice took 727.29 ms                                   [sph::CartesianRender][rank=0]

Cylindrical rendering with custom getter (vtheta but with a mask)

406 positions_cylindrical = make_cylindrical_coords(nr, ntheta)
407 rays_cylindrical = positions_to_rays(positions_cylindrical)
408
409
410 def custom_getter(size, dic_out):
411     x = dic_out["xyz"][:, 0]
412     y = dic_out["xyz"][:, 1]
413     z = dic_out["xyz"][:, 2]
414     vx = dic_out["vxyz"][:, 0]
415     vy = dic_out["vxyz"][:, 1]
416     vz = dic_out["vxyz"][:, 2]
417
418     v_theta = np.zeros(size)
419     for i in range(size):
420         e_theta = np.array([-y[i], x[i], 0])
421         e_theta /= np.linalg.norm(e_theta) + 1e-9  # Avoid division by zero
422         v_theta[i] = np.dot(e_theta, np.array([vx[i], vy[i], vz[i]]))
423
424         if x[i] > 0.2:
425             v_theta[i] = 0.0  # To show that we have full control on the rendering
426
427     return v_theta
428
429
430 arr_vtheta_cylindrical = model.render_column_integ("custom", "f64", rays_cylindrical, custom_getter)
431 arr_vtheta_pos = model.render_slice("custom", "f64", positions_cylindrical, custom_getter)
432
433
434 def plot_vtheta_integ_cylindrical(metadata, arr_vtheta_cylindrical):
435     ext = metadata["extent"]
436
437     my_cmap = matplotlib.colormaps["gist_heat"].copy()  # copy the default cmap
438     my_cmap.set_bad(color="black")
439
440     arr_vtheta_cylindrical = np.array(arr_vtheta_cylindrical).reshape(nr, ntheta)
441
442     res = plt.imshow(
443         arr_vtheta_cylindrical,
444         cmap=my_cmap,
445         origin="lower",
446         extent=ext,
447         norm="log",
448         vmin=1e-7,
449         vmax=1e-5,
450         aspect="auto",
451     )
452     plt.xlabel("r")
453     plt.ylabel(r"$\theta$")
454     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
455     cbar = plt.colorbar(res, extend="both")
456     cbar.set_label(r"$\int v_\theta \, \mathrm{d}z$ [code unit]")
457
458
459 def plot_vtheta_slice_cylindrical(metadata, arr_vtheta_pos):
460     ext = metadata["extent"]
461
462     my_cmap = matplotlib.colormaps["gist_heat"].copy()  # copy the default cmap
463     my_cmap.set_bad(color="black")
464
465     arr_vtheta_pos = np.array(arr_vtheta_pos).reshape(nr, ntheta)
466
467     res = plt.imshow(
468         arr_vtheta_pos,
469         cmap=my_cmap,
470         origin="lower",
471         extent=ext,
472         norm="log",
473         vmin=1e-8,
474         aspect="auto",
475     )
476     plt.xlabel("r")
477     plt.ylabel(r"$\theta$")
478     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
479     cbar = plt.colorbar(res, extend="both")
480     cbar.set_label(r"$v_\theta$ [code unit]")
481
482
483 metadata = {"extent": [0, ext, -np.pi, np.pi], "time": model.get_time()}
484
485 plt.figure(dpi=dpi)
486 plot_vtheta_integ_cylindrical(metadata, arr_vtheta_cylindrical)
487
488 plt.figure(dpi=dpi)
489 plot_vtheta_slice_cylindrical(metadata, arr_vtheta_pos)
490
491 plt.show()
  • t = 4.515 [seconds]
  • t = 4.515 [seconds]
Info: compute_column_integ field_name: custom, rays count: 1048576   [sph::CartesianRender][rank=0]
sph::RenderFieldGetter compute custom field took :  0.41132996200000005 s
Info: compute_column_integ took 2.12 s                               [sph::CartesianRender][rank=0]
Info: compute_slice field_name: custom, positions count: 1048576     [sph::CartesianRender][rank=0]
sph::RenderFieldGetter compute custom field took :  0.43110169000000004 s
Info: compute_slice took 1.16 s                                      [sph::CartesianRender][rank=0]

Azymuthal rendering

496 H_r_render = H_r_0 * 4
497
498
499 def make_azymuthal_coords(nr, nz):
500     """
501     Generate a list of positions in cylindrical coordinates (r, theta)
502     spanning [0, ext*2] x [-pi, pi] for use with the rendering module.
503
504     Returns:
505         list: List of [x, y, z] coordinate lists
506     """
507
508     # Create the cylindrical coordinate grid
509     r_vals = np.linspace(0, ext, nr)
510     z_vals = np.linspace(-ext * H_r_render, ext * H_r_render, nz)
511
512     # Create meshgrid
513     r_grid, z_grid = np.meshgrid(r_vals, z_vals)
514
515     # Flatten and stack to create list of positions
516     positions = np.column_stack([r_grid.ravel(), z_grid.ravel()])
517
518     return [tuple(pos) for pos in positions]
519
520
521 def make_ring_rays(positions):
522     def position_to_ring_ray(position):
523         r = position[0]
524         z = position[1]
525         e_x = (1.0, 0.0, 0.0)
526         e_y = (0.0, 1.0, 0.0)
527         center = (0.0, 0.0, z)
528         return shamrock.math.RingRay_f64_3(center, r, e_x, e_y)
529
530     return [position_to_ring_ray(position) for position in positions]
531
532
533 def make_slice_coord_for_azymuthal(positions):
534     def position_to_ring_ray(position):
535         r = position[0]
536         z = position[1]
537         e_x = (1.0, 0.0, 0.0)
538         e_y = (0.0, 1.0, 0.0)
539         center = (0.0, 0.0, z)
540         return (r, 0.0, z)
541
542     return [position_to_ring_ray(position) for position in positions]
543
544
545 nr = 1024
546 nz = 1024
547
548 positions_azymuthal = make_azymuthal_coords(nr, nz)
549 ring_rays_azymuthal = make_ring_rays(positions_azymuthal)
550 slice_coords_azymuthal = make_slice_coord_for_azymuthal(positions_azymuthal)
551
552 arr_rho_azymuthal = model.render_azymuthal_integ("rho", "f64", ring_rays_azymuthal)
553 arr_rho_slice_azymuthal = model.render_slice("rho", "f64", slice_coords_azymuthal)
554
555 arr_vxyz_azymuthal = model.render_azymuthal_integ("vxyz", "f64_3", ring_rays_azymuthal)
556 arr_vxyz_slice_azymuthal = model.render_slice("vxyz", "f64_3", slice_coords_azymuthal)
557
558
559 def plot_rho_integ_azymuthal(metadata, arr_rho_azymuthal):
560     ext = metadata["extent"]
561
562     my_cmap = matplotlib.colormaps["gist_heat"].copy()  # copy the default cmap
563     my_cmap.set_bad(color="black")
564
565     arr_rho_azymuthal = np.array(arr_rho_azymuthal).reshape(nr, nz)
566
567     res = plt.imshow(
568         arr_rho_azymuthal, cmap=my_cmap, origin="lower", extent=ext, norm="log", vmin=1e-5, vmax=1
569     )
570     plt.xlabel("r")
571     plt.ylabel("z")
572     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
573     cbar = plt.colorbar(res, extend="both")
574     cbar.set_label(r"$\int \rho \, \mathrm{d}\theta$ [code unit]")
575
576
577 def plot_rho_slice_azymuthal(metadata, arr_rho_slice_azymuthal):
578     ext = metadata["extent"]
579
580     my_cmap = matplotlib.colormaps["gist_heat"].copy()  # copy the default cmap
581     my_cmap.set_bad(color="black")
582
583     arr_rho_slice_azymuthal = np.array(arr_rho_slice_azymuthal).reshape(nr, nz)
584
585     res = plt.imshow(
586         arr_rho_slice_azymuthal,
587         cmap=my_cmap,
588         origin="lower",
589         extent=ext,
590         norm="log",
591         vmin=1e-5,
592         vmax=1,
593     )
594     plt.xlabel("r")
595     plt.ylabel("z")
596     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
597     cbar = plt.colorbar(res, extend="both")
598     cbar.set_label(r"$\rho$ [code unit]")
599
600
601 def plot_vz_integ_azymuthal(metadata, arr_vxyz_azymuthal):
602     ext = metadata["extent"]
603
604     my_cmap = matplotlib.colormaps["seismic"].copy()  # copy the default cmap
605     my_cmap.set_bad(color="black")
606
607     arr_vz_azymuthal = np.array(arr_vxyz_azymuthal).reshape(nr, nz, 3)[:, :, 2]
608
609     res = plt.imshow(
610         arr_vz_azymuthal, cmap=my_cmap, origin="lower", extent=ext, vmin=-1e-6, vmax=1e-6
611     )
612     plt.xlabel("r")
613     plt.ylabel("z")
614     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
615     cbar = plt.colorbar(res, extend="both")
616     cbar.set_label(r"$\int v_z \, \mathrm{d}\theta$ [code unit]")
617
618
619 def plot_vz_slice_azymuthal(metadata, arr_vxyz_slice_azymuthal):
620     ext = metadata["extent"]
621
622     my_cmap = matplotlib.colormaps["seismic"].copy()  # copy the default cmap
623     my_cmap.set_bad(color="black")
624
625     arr_vz_slice_azymuthal = np.array(arr_vxyz_slice_azymuthal).reshape(nr, nz, 3)[:, :, 2]
626
627     res = plt.imshow(
628         arr_vz_slice_azymuthal, cmap=my_cmap, origin="lower", extent=ext, vmin=-5e-6, vmax=5e-6
629     )
630     plt.xlabel("r")
631     plt.ylabel("z")
632     plt.title(f"t = {metadata['time']:0.3f} [seconds]")
633     cbar = plt.colorbar(res, extend="both")
634     cbar.set_label(r"$v_z$ [code unit]")
635
636
637 metadata = {"extent": [0, ext, -ext * H_r_render, ext * H_r_render], "time": model.get_time()}
638 fig_size = (6, 3)
639 plt.figure(dpi=dpi, figsize=fig_size)
640 plot_rho_integ_azymuthal(metadata, arr_rho_azymuthal)
641
642 plt.figure(dpi=dpi, figsize=fig_size)
643 plot_rho_slice_azymuthal(metadata, arr_rho_slice_azymuthal)
644
645 plt.figure(dpi=dpi, figsize=fig_size)
646 plot_vz_integ_azymuthal(metadata, arr_vxyz_azymuthal)
647
648 plt.figure(dpi=dpi, figsize=fig_size)
649 plot_vz_slice_azymuthal(metadata, arr_vxyz_slice_azymuthal)
650
651 plt.show()
  • t = 4.515 [seconds]
  • t = 4.515 [seconds]
  • t = 4.515 [seconds]
  • t = 4.515 [seconds]
Info: compute_azymuthal_integ field_name: rho, ring_rays count: 1048576  [sph::CartesianRender][rank=0]
Info: compute_azymuthal_integ took 31.78 s                           [sph::CartesianRender][rank=0]
Info: compute_slice field_name: rho, positions count: 1048576        [sph::CartesianRender][rank=0]
Info: compute_slice took 345.72 ms                                   [sph::CartesianRender][rank=0]
Info: compute_azymuthal_integ field_name: vxyz, ring_rays count: 1048576  [sph::CartesianRender][rank=0]
Info: compute_azymuthal_integ took 31.75 s                           [sph::CartesianRender][rank=0]
Info: compute_slice field_name: vxyz, positions count: 1048576       [sph::CartesianRender][rank=0]
Info: compute_slice took 308.32 ms                                   [sph::CartesianRender][rank=0]

Total running time of the script: (1 minutes 33.271 seconds)

Estimated memory usage: 1636 MB

Gallery generated by Sphinx-Gallery