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 # 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.03,
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.do_vtk_dump("init_disc.vtk", True)
212
213 model.change_htolerances(coarse=1.3, fine=1.1)
214 model.timestep()
215 model.change_htolerances(coarse=1.1, fine=1.1)
216
217 for i in range(5):
218     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": 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,
        "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.685470e+05 N.s^-1
SPH setup: the generation step took : 0.37931531700000004 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.51 us    (63.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     : 1503.00 ns (0.3%)
   patch tree reduce : 1323.00 ns (0.2%)
   gen split merge   : 781.00 ns  (0.1%)
   split / merge op  : 0/0
   apply split merge : 972.00 ns  (0.2%)
   LB compute        : 514.32 us  (95.9%)
   LB move op cnt    : 0
   LB apply          : 13.60 us   (2.5%)
Info: Compute load ...                                                [DataInserterUtility][rank=0]
Info: run scheduler step ...                                          [DataInserterUtility][rank=0]
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.81 us    (66.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     : 1623.00 ns (0.4%)
   patch tree reduce : 531.00 ns  (0.1%)
   gen split merge   : 371.00 ns  (0.1%)
   split / merge op  : 0/0
   apply split merge : 320.00 ns  (0.1%)
   LB compute        : 399.09 us  (98.1%)
   LB move op cnt    : 0
   LB apply          : 1933.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    (65.9%)
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     : 2.13 us    (0.6%)
   patch tree reduce : 481.00 ns  (0.1%)
   gen split merge   : 371.00 ns  (0.1%)
   split / merge op  : 0/0
   apply split merge : 291.00 ns  (0.1%)
   LB compute        : 372.04 us  (97.8%)
   LB move op cnt    : 0
   LB apply          : 1874.00 ns (0.5%)
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.011817476 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 |   1.5% 0.0% |     1.16 GB |     5.04 MB |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 0.416288394 s
Info: dump to init_disc.vtk                                                      [VTK Dump][rank=0]
              - took 15.09 ms, bandwidth = 353.93 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     : 6.58 us    (1.6%)
   patch tree reduce : 1893.00 ns (0.5%)
   gen split merge   : 992.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 901.00 ns  (0.2%)
   LB compute        : 396.65 us  (94.7%)
   LB move op cnt    : 0
   LB apply          : 4.18 us    (1.0%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.96 us    (70.9%)
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.8084e+04 | 100000 |      1 | 3.561e+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     : 6.78 us    (1.7%)
   patch tree reduce : 1923.00 ns (0.5%)
   gen split merge   : 1032.00 ns (0.3%)
   split / merge op  : 0/0
   apply split merge : 1042.00 ns (0.3%)
   LB compute        : 378.04 us  (94.7%)
   LB move op cnt    : 0
   LB apply          : 3.54 us    (0.9%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.35 us    (71.8%)
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.5827e+05 | 100000 |      1 | 6.318e-01 | 0.0% |   0.1% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 110.58074372649077 (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     : 6.74 us    (1.8%)
   patch tree reduce : 1983.00 ns (0.5%)
   gen split merge   : 1222.00 ns (0.3%)
   split / merge op  : 0/0
   apply split merge : 1012.00 ns (0.3%)
   LB compute        : 357.96 us  (94.4%)
   LB move op cnt    : 0
   LB apply          : 3.63 us    (1.0%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.35 us    (71.6%)
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.5825e+05 | 100000 |      1 | 6.319e-01 | 0.0% |   0.1% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 3757.182547251058 (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     : 5.81 us    (1.6%)
   patch tree reduce : 1964.00 ns (0.5%)
   gen split merge   : 892.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 1022.00 ns (0.3%)
   LB compute        : 341.49 us  (94.5%)
   LB move op cnt    : 0
   LB apply          : 3.41 us    (0.9%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.43 us    (66.6%)
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.5787e+05 | 100000 |      1 | 6.334e-01 | 0.0% |   0.1% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 6066.83923959143 (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     : 6.57 us    (1.7%)
   patch tree reduce : 1943.00 ns (0.5%)
   gen split merge   : 942.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 1073.00 ns (0.3%)
   LB compute        : 357.42 us  (94.6%)
   LB move op cnt    : 0
   LB apply          : 3.53 us    (0.9%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.40 us    (72.7%)
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.5801e+05 | 100000 |      1 | 6.329e-01 | 0.0% |   0.1% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 7470.029145817702 (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.83 us    (2.0%)
   patch tree reduce : 2.46 us    (0.6%)
   gen split merge   : 951.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 1132.00 ns (0.3%)
   LB compute        : 347.07 us  (89.6%)
   LB move op cnt    : 0
   LB apply          : 3.59 us    (0.9%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 3.07 us    (70.2%)
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.5744e+05 | 100000 |      1 | 6.352e-01 | 0.0% |   0.1% 0.0% |     1.16 GB |     5.04 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 8251.390684623859 (tsim/hr)                             [sph::Model][rank=0]

Usual cartesian rendering

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

Cylindrical rendering

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

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

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

Azymuthal rendering

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

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

Estimated memory usage: 1168 MB

Gallery generated by Sphinx-Gallery