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.05
 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 cfg.set_scheduler_config(split_load_value=scheduler_split_val, merge_load_value=scheduler_merge_val)
178
179 # Set the solver config to be the one stored in cfg
180 model.set_solver_config(cfg)
181
182 # Print the solver config
183 model.get_current_config().print_status()
184
185 # Init the scheduler & fields
186 model.init()
187
188 # Set the simulation box size
189 model.resize_simulation_box(bmin, bmax)
190
191 # Create the setup
192
193 setup = model.get_setup()
194 gen_disc = setup.make_generator_disc_mc(
195     part_mass=pmass,
196     disc_mass=disc_mass,
197     r_in=rin,
198     r_out=rout,
199     sigma_profile=sigma_profile,
200     H_profile=H_profile,
201     rot_profile=rot_profile,
202     cs_profile=cs_profile,
203     random_seed=666,
204     init_h_factor=0.03,
205 )
206
207 # Print the dot graph of the setup
208 print(gen_disc.get_dot())
209
210 # Apply the setup
211 setup.apply_setup(gen_disc)
212
213 model.do_vtk_dump("init_disc.vtk", True)
214
215 model.change_htolerances(coarse=1.3, fine=1.1)
216 model.timestep()
217 model.change_htolerances(coarse=1.1, fine=1.1)
218
219 for i in range(5):
220     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,
        "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_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
    }
]
------------------------------------
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 ) rate = 2.420389e+05 N.s^-1
SPH setup: the generation step took : 0.42429465600000005 s
SPH setup: final particle count = 100000 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     : 8.42 us    (63.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     : 1623.00 ns (0.3%)
   patch tree reduce : 1933.00 ns (0.3%)
   gen split merge   : 762.00 ns  (0.1%)
   split / merge op  : 0/0
   apply split merge : 932.00 ns  (0.2%)
   LB compute        : 556.05 us  (97.8%)
   LB move op cnt    : 0
   LB apply          : 3.73 us    (0.7%)
Info: Compute load ...                                                [DataInserterUtility][rank=0]
Info: run scheduler step ...                                          [DataInserterUtility][rank=0]
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.58 us    (63.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     : 1924.00 ns (0.5%)
   patch tree reduce : 521.00 ns  (0.1%)
   gen split merge   : 370.00 ns  (0.1%)
   split / merge op  : 0/0
   apply split merge : 311.00 ns  (0.1%)
   LB compute        : 375.24 us  (98.0%)
   LB move op cnt    : 0
   LB apply          : 1703.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.37 us    (65.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     : 1953.00 ns (0.5%)
   patch tree reduce : 521.00 ns  (0.1%)
   gen split merge   : 361.00 ns  (0.1%)
   split / merge op  : 0/0
   apply split merge : 280.00 ns  (0.1%)
   LB compute        : 371.01 us  (98.0%)
   LB move op cnt    : 0
   LB apply          : 1633.00 ns (0.4%)
Info: ---------------------------------------------                   [DataInserterUtility][rank=0]
SPH setup: injected       100000 / 100000 => 100.0% | ranks with patchs = 1 / 1  <- global loop ->
SPH setup: the injection step took : 0.020877974 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.7% 0.0% |     1.16 GB |     5.04 MB |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 0.49140697800000005 s
Info: dump to init_disc.vtk                                                      [VTK Dump][rank=0]
              - took 24.05 ms, bandwidth = 222.09 MB/s
---------------- 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     : 7.87 us    (0.4%)
   patch tree reduce : 2.18 us    (0.1%)
   gen split merge   : 992.00 ns  (0.0%)
   split / merge op  : 0/0
   apply split merge : 1012.00 ns (0.0%)
   LB compute        : 2.11 ms    (99.0%)
   LB move op cnt    : 0
   LB apply          : 3.51 us    (0.2%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.84 us    (66.4%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.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.293893919894618e-10,3.112830370992019e-11,0)
    sum a = (2.6787416956917247e-28,-3.536362554786704e-28,1.5651051282891335e-26)
    sum e = 1.2757522906384626e-10
    sum de = 1.7517536550417152e-31
Info: CFL hydro = 0.019407933351461047 sink sink = inf                                [SPH][rank=0]
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.0419e+04 | 100000 |      1 | 4.897e+00 | 0.0% |   0.1% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0 (tsim/hr)                                             [sph::Model][rank=0]
---------------- t = 0, dt = 0.019407933351461047 ----------------
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.87 us    (0.1%)
   patch tree reduce : 2.11 us    (0.0%)
   gen split merge   : 1192.00 ns (0.0%)
   split / merge op  : 0/0
   apply split merge : 1062.00 ns (0.0%)
   LB compute        : 7.73 ms    (99.7%)
   LB move op cnt    : 0
   LB apply          : 5.58 us    (0.1%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 5.97 us    (82.2%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919895638e-10,3.112830370990416e-11,-8.804510440014427e-29)
    sum a = (1.1974928291807671e-27,-3.7269449679189217e-28,-8.116693217064333e-27)
    sum e = 1.2757522909678268e-10
    sum de = 2.3875512779325514e-32
Info: CFL hydro = 0.6595177856935768 sink sink = inf                                  [SPH][rank=0]
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.3802e+05 | 100000 |      1 | 7.245e-01 | 0.0% |   0.1% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 96.43246434316102 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 0.019407933351461047, dt = 0.6595177856935768 ----------------
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.01 us    (0.1%)
   patch tree reduce : 1954.00 ns (0.0%)
   gen split merge   : 1312.00 ns (0.0%)
   split / merge op  : 0/0
   apply split merge : 1212.00 ns (0.0%)
   LB compute        : 7.42 ms    (99.7%)
   LB move op cnt    : 0
   LB apply          : 4.20 us    (0.1%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 4.41 us    (76.0%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919894064e-10,3.112830370980896e-11,-2.168404344971009e-27)
    sum a = (8.92560968169219e-28,-8.385626177817574e-28,1.8003685293890153e-26)
    sum e = 1.2757526707478419e-10
    sum de = -2.2918566338208107e-32
Info: CFL hydro = 1.067486611711989 sink sink = inf                                   [SPH][rank=0]
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.3539e+05 | 100000 |      1 | 7.386e-01 | 0.0% |   0.0% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 3214.631552802381 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 0.6789257190450378, dt = 1.067486611711989 ----------------
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.08 us    (0.0%)
   patch tree reduce : 2.14 us    (0.0%)
   gen split merge   : 1192.00 ns (0.0%)
   split / merge op  : 0/0
   apply split merge : 942.00 ns  (0.0%)
   LB compute        : 18.71 ms   (99.9%)
   LB move op cnt    : 0
   LB apply          : 6.87 us    (0.0%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 7.93 us    (84.7%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919893851e-10,3.1128303709784325e-11,-2.1950858828095195e-26)
    sum a = (2.3081647812679683e-28,8.470329472543003e-29,-3.14460981668159e-27)
    sum e = 1.275753261081859e-10
    sum de = 1.2414123124344788e-31
Info: CFL hydro = 1.3131733244842363 sink sink = inf                                  [SPH][rank=0]
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.3614e+05 | 100000 |      1 | 7.345e-01 | 0.0% |   0.0% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 5231.946309993719 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 1.746412330757027, dt = 1.3131733244842363 ----------------
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.37 us    (0.1%)
   patch tree reduce : 1963.00 ns (0.0%)
   gen split merge   : 1343.00 ns (0.0%)
   split / merge op  : 0/0
   apply split merge : 1072.00 ns (0.0%)
   LB compute        : 9.04 ms    (99.7%)
   LB move op cnt    : 0
   LB apply          : 4.25 us    (0.0%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.92 us    (70.6%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919894716e-10,3.112830371004766e-11,-2.5995441151234477e-26)
    sum a = (9.020900888258298e-28,2.8693241088239423e-28,-1.5208476567950963e-26)
    sum e = 1.2757537036686392e-10
    sum de = -5.478066985566444e-32
Info: CFL hydro = 1.4558626079520551 sink sink = inf                                  [SPH][rank=0]
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.3208e+05 | 100000 |      1 | 7.571e-01 | 0.0% |   0.0% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 6244.2027448766 (tsim/hr)                               [sph::Model][rank=0]
---------------- t = 3.0595856552412632, dt = 1.4558626079520551 ----------------
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.06 us    (0.1%)
   patch tree reduce : 2.09 us    (0.0%)
   gen split merge   : 1292.00 ns (0.0%)
   split / merge op  : 0/0
   apply split merge : 1022.00 ns (0.0%)
   LB compute        : 9.12 ms    (99.7%)
   LB move op cnt    : 0
   LB apply          : 4.09 us    (0.0%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 3.12 us    (73.9%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919895108e-10,3.112830370994267e-11,-3.118775311790334e-26)
    sum a = (7.369186641112413e-28,-1.122318655111948e-28,-4.493509785184063e-27)
    sum e = 1.2757539490951678e-10
    sum de = -9.534256004378249e-32
Info: CFL hydro = 1.7417916889876892 sink sink = inf                                  [SPH][rank=0]
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.3095e+05 | 100000 |      1 | 7.637e-01 | 0.0% |   0.0% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 6863.124855575635 (tsim/hr)                             [sph::Model][rank=0]

Usual cartesian rendering

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

Cylindrical rendering

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

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

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

Azymuthal rendering

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

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

Estimated memory usage: 1602 MB

Gallery generated by Sphinx-Gallery