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 from shamrock import NeighCacheStrategy
23
24 # If we use the shamrock executable to run this script instead of the python interpreter,
25 # we should not initialize the system as the shamrock executable needs to handle specific MPI logic
26 if not shamrock.sys.is_initialized():
27     shamrock.change_loglevel(1)
28     shamrock.sys.init("0:0")

Use shamrock documentation style for matplotlib

33 shamrock.matplotlib.set_shamrock_mpl_style()

Setup units

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

List parameters

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

Orbital period :  247.58972132551145  [seconds]

Utility functions and quantities deduced from the base one

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

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

Attach a SPH model to the context

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

----- SPH Solver configuration -----
[
    {
        "artif_viscosity": {
            "alpha_AV": 0.0125,
            "alpha_u": 1.0,
            "beta_AV": 2.0,
            "type": "constant_disc"
        },
        "boundary_config": {
            "bc_type": "free"
        },
        "cfl_config": {
            "cfl_cour": 0.3,
            "cfl_force": 0.25,
            "cfl_multiplier_stiffness": 2.0,
            "eta_sink": 0.05
        },
        "combined_dtdiv_divcurlv_compute": false,
        "debug_dump_filename": "",
        "do_debug_dump": false,
        "dust_config": {
            "ballabio_ts_limiter": false,
            "drag_mode": {
                "type": "none"
            },
            "evol_mode": {
                "type": "none"
            },
            "mode": {
                "type": "none"
            }
        },
        "enable_particle_reordering": false,
        "eos_config": {
            "Tvec": "f64_3",
            "cs0": 5.009972010250009e-05,
            "eos_type": "locally_isothermal_lp07",
            "q": 0.75,
            "r0": 0.03948371767577914
        },
        "epsilon_h": 1e-06,
        "ext_force_config": {
            "force_list": []
        },
        "gpart_mass": 1e-08,
        "h_iter_per_subcycles": 50,
        "h_max_subcycles_count": 100,
        "htol_up_coarse_cycle": 1.1,
        "htol_up_fine_cycle": 1.1,
        "kernel_id": "M4<f64>",
        "mhd_config": {
            "mhd_type": "none"
        },
        "neigh_cache_strategy": "single_stage",
        "particle_killing": [
            {
                "center": [
                    0.0,
                    0.0,
                    0.0
                ],
                "radius": 0.7896743535155828,
                "type": "sphere"
            }
        ],
        "particle_reordering_step_freq": 1000,
        "save_dt_to_fields": false,
        "scheduler_config": {
            "merge_load_value": 625000,
            "split_load_value": 10000000
        },
        "self_grav_config": {
            "softening_length": 1e-09,
            "softening_mode": "plummer",
            "type": "none"
        },
        "show_cfl_detail": false,
        "show_ghost_zone_graph": false,
        "show_neigh_stats": false,
        "smoothing_length_config": {
            "type": "density_based"
        },
        "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
        }
    }
]
------------------------------------
Warning: make_generator_disc_mc: with the current EOS, cs_profile is ignored     [SPHSetup][rank=0]

Show the dot graph of the setup

215 shamrock.utils.plot.DotGraph(gen_disc.get_dot())
G node_0 GeneratorMCDisc node_2 Simulation node_0->node_2


Apply the setup

219 setup.apply_setup(gen_disc)
220
221 model.do_vtk_dump("init_disc.vtk", True)
222
223 model.change_htolerances(coarse=1.3, fine=1.1)
224 model.timestep()
225 model.change_htolerances(coarse=1.1, fine=1.1)
226
227 for i in range(5):
228     model.timestep()

SPH setup: generating particles ...
SPH setup: Nstep = 100000 ( 1.0e+05 ) Ntotal = 100000 ( 1.0e+05 rank min = 2.8e+05 max = 1.0e+05) rate = 1.000000e+05 N.s^-1
SPH setup: the generation step took : 0.37184818000000003 s
SPH setup: final particle count = 100000 beginning injection ...
Info: ---------------------------------------------                   [DataInserterUtility][rank=0]
Info: Compute load ...                                                [DataInserterUtility][rank=0]
Info: run scheduler step ...                                          [DataInserterUtility][rank=0]
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 16.67 us   (67.3%)
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
Info: Loadbalance stats :                                                     [LoadBalance][rank=0]
    npatch = 1
    min = 100000
    max = 100000
    avg = 100000
    efficiency = 100.00%
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 1.70 us    (0.3%)
   patch tree reduce : 2.75 us    (0.5%)
   gen split merge   : 722.00 ns  (0.1%)
   split / merge op  : 0/0
   apply split merge : 832.00 ns  (0.1%)
   LB compute        : 574.85 us  (96.2%)
   LB move op cnt    : 0
   LB apply          : 12.63 us   (2.1%)
Info: patch count stable after 1 runs npatch = 1                      [DataInserterUtility][rank=0]
Info: ---------------------------------------------                   [DataInserterUtility][rank=0]
SPH setup: injected       100000 / 100000 => 100.0% | ranks with patchs = 1 / 1  <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.010452877000000001 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.2% 0.0% |     1.26 GB |     5.29 MB |
+------+--------------------+-------+-------------+-------------+-------------+
SPH setup: the setup took : 0.399070464 s
Info: dump to init_disc.vtk                                                      [VTK Dump][rank=0]
              - took 14.93 ms, bandwidth = 375.22 MB/s
---------------- t = 0, dt = 0 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
Info: Loadbalance stats :                                                     [LoadBalance][rank=0]
    npatch = 1
    min = 100000
    max = 100000
    avg = 100000
    efficiency = 100.00%
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 7.78 us    (1.9%)
   patch tree reduce : 2.19 us    (0.5%)
   gen split merge   : 1.04 us    (0.2%)
   split / merge op  : 0/0
   apply split merge : 882.00 ns  (0.2%)
   LB compute        : 396.21 us  (94.6%)
   LB move op cnt    : 0
   LB apply          : 3.58 us    (0.9%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.44 us    (67.8%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.006569301129452108 unconverged cnt = 100000
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.008540091468287742 unconverged cnt = 100000
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.011102118908774064 unconverged cnt = 100000
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.014432754581406283 unconverged cnt = 99999
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.018762580955828168 unconverged cnt = 99999
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02439135524257662 unconverged cnt = 99995
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.024830299239077123 unconverged cnt = 99988
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.024830299239077126 unconverged cnt = 99974
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 99926
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 99754
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 98899
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 86081
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 16466
Warning: smoothing length is not converged, rerunning the iterator ...    [Smoothinglength][rank=0]
     largest h = 0.02483029923907713 unconverged cnt = 82
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919894623e-10,3.112830370992019e-11,0)
    sum a = (-2.858736196983264e-29,-4.0234064994579267e-28,8.440683319389103e-27)
    sum e = 1.275752290638463e-10
    sum de = 1.5424291442111275e-31
Info: cfl dt = 0.019407933351461047 cfl multiplier : 0.01                      [sph::Model][rank=0]
Info: Timestep perf report:                                                    [sph::Model][rank=0]
+======+============+========+========+===========+======+=============+=============+=============+
| rank | rate (N/s) |  Nobj  | Npatch |   tstep   | MPI  | alloc d% h% | mem (max) d | mem (max) h |
+======+============+========+========+===========+======+=============+=============+=============+
| 0    | 2.6128e+04 | 100000 |      1 | 3.827e+00 | 0.0% |   0.1% 0.0% |     1.26 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 0 (tsim/hr)                                             [sph::Model][rank=0]
---------------- t = 0, dt = 0.019407933351461047 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
Info: Loadbalance stats :                                                     [LoadBalance][rank=0]
    npatch = 1
    min = 100000
    max = 100000
    avg = 100000
    efficiency = 100.00%
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 6.87 us    (1.8%)
   patch tree reduce : 1.61 us    (0.4%)
   gen split merge   : 862.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 972.00 ns  (0.3%)
   LB compute        : 354.93 us  (94.5%)
   LB move op cnt    : 0
   LB apply          : 3.99 us    (1.1%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.44 us    (67.0%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919895486e-10,3.112830370990416e-11,-1.64741290796061e-28)
    sum a = (9.200895389549837e-28,4.891615270393584e-28,-1.0590029423046891e-26)
    sum e = 1.2757522909678268e-10
    sum de = 6.77927340424307e-32
Info: cfl dt = 0.6595177856935768 cfl multiplier : 0.34                        [sph::Model][rank=0]
Info: Timestep perf report:                                                    [sph::Model][rank=0]
+======+============+========+========+===========+======+=============+=============+=============+
| rank | rate (N/s) |  Nobj  | Npatch |   tstep   | MPI  | alloc d% h% | mem (max) d | mem (max) h |
+======+============+========+========+===========+======+=============+=============+=============+
| 0    | 1.6068e+05 | 100000 |      1 | 6.223e-01 | 0.0% |   0.0% 0.0% |     1.26 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 112.26594440688564 (tsim/hr)                            [sph::Model][rank=0]
---------------- t = 0.019407933351461047, dt = 0.6595177856935768 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
Info: Loadbalance stats :                                                     [LoadBalance][rank=0]
    npatch = 1
    min = 100000
    max = 100000
    avg = 100000
    efficiency = 100.00%
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 7.36 us    (1.6%)
   patch tree reduce : 1.65 us    (0.4%)
   gen split merge   : 832.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 1.19 us    (0.3%)
   LB compute        : 444.94 us  (95.4%)
   LB move op cnt    : 0
   LB apply          : 4.40 us    (0.9%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.76 us    (68.9%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919894062e-10,3.112830370980896e-11,-4.27116363652981e-27)
    sum a = (8.33268661861418e-28,-7.464477847678522e-28,2.031820282226253e-26)
    sum e = 1.2757526707478419e-10
    sum de = 1.8631566652623644e-31
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.5673e+05 | 100000 |      1 | 6.380e-01 | 0.0% |   0.0% 0.0% |     1.26 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 3721.131250343549 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 0.6789257190450378, dt = 1.067486611711989 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
Info: Loadbalance stats :                                                     [LoadBalance][rank=0]
    npatch = 1
    min = 100000
    max = 100000
    avg = 100000
    efficiency = 100.00%
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 6.64 us    (1.7%)
   patch tree reduce : 1.94 us    (0.5%)
   gen split merge   : 922.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 1.03 us    (0.3%)
   LB compute        : 379.57 us  (94.8%)
   LB move op cnt    : 0
   LB apply          : 4.32 us    (1.1%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.42 us    (66.2%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919893851e-10,3.1128303709784325e-11,-1.6606080930920558e-26)
    sum a = (-2.0011153378882847e-28,-4.6480932980579734e-28,-9.740878893424454e-28)
    sum e = 1.275753261081859e-10
    sum de = 2.6524569232519395e-31
Info: cfl dt = 1.3131733244842363 cfl multiplier : 0.7066666666666667          [sph::Model][rank=0]
Info: Timestep perf report:                                                    [sph::Model][rank=0]
+======+============+========+========+===========+======+=============+=============+=============+
| rank | rate (N/s) |  Nobj  | Npatch |   tstep   | MPI  | alloc d% h% | mem (max) d | mem (max) h |
+======+============+========+========+===========+======+=============+=============+=============+
| 0    | 1.5054e+05 | 100000 |      1 | 6.643e-01 | 0.0% |   0.0% 0.0% |     1.26 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 5785.359060006857 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 1.746412330757027, dt = 1.3131733244842363 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
Info: Loadbalance stats :                                                     [LoadBalance][rank=0]
    npatch = 1
    min = 100000
    max = 100000
    avg = 100000
    efficiency = 100.00%
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 6.85 us    (1.9%)
   patch tree reduce : 1.46 us    (0.4%)
   gen split merge   : 922.00 ns  (0.3%)
   split / merge op  : 0/0
   apply split merge : 1.02 us    (0.3%)
   LB compute        : 339.92 us  (94.4%)
   LB move op cnt    : 0
   LB apply          : 3.62 us    (1.0%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.29 us    (66.8%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.29389391989472e-10,3.112830371005654e-11,-2.816808066094176e-26)
    sum a = (7.697411908173454e-28,-8.364450354136216e-29,-1.5337649092407243e-26)
    sum e = 1.2757537036686392e-10
    sum de = 9.331732433663779e-33
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.5637e+05 | 100000 |      1 | 6.395e-01 | 0.0% |   0.0% 0.0% |     1.26 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 7392.352471006853 (tsim/hr)                             [sph::Model][rank=0]
---------------- t = 3.0595856552412632, dt = 1.4558626079520551 ----------------
Info: Summary (strategy = round robin):                                       [LoadBalance][rank=0]
 - strategy "psweep"      : max = 100000.0 min = 100000.0 factor = 1
 - strategy "round robin" : max = 95000.0 min = 95000.0 factor = 0.95
Info: Loadbalance stats :                                                     [LoadBalance][rank=0]
    npatch = 1
    min = 100000
    max = 100000
    avg = 100000
    efficiency = 100.00%
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 6.65 us    (1.8%)
   patch tree reduce : 1.59 us    (0.4%)
   gen split merge   : 831.00 ns  (0.2%)
   split / merge op  : 0/0
   apply split merge : 902.00 ns  (0.2%)
   LB compute        : 352.38 us  (94.6%)
   LB move op cnt    : 0
   LB apply          : 3.98 us    (1.1%)
Info: Scheduler step timings :                                                  [Scheduler][rank=0]
   metadata sync     : 2.46 us    (66.8%)
Info: free boundaries skipping geometry update                            [PositionUpdated][rank=0]
Info: conservation infos :                                                     [sph::Model][rank=0]
    sum v = (-9.293893919895108e-10,3.11283037099426e-11,-1.2934193104573167e-26)
    sum a = (4.309280119156253e-28,-9.211483301390516e-28,-2.689329607532404e-28)
    sum e = 1.2757539490951678e-10
    sum de = -2.026529691326088e-31
Info: cfl dt = 1.7417916889876892 cfl multiplier : 0.8696296296296296          [sph::Model][rank=0]
Info: Timestep perf report:                                                    [sph::Model][rank=0]
+======+============+========+========+===========+======+=============+=============+=============+
| rank | rate (N/s) |  Nobj  | Npatch |   tstep   | MPI  | alloc d% h% | mem (max) d | mem (max) h |
+======+============+========+========+===========+======+=============+=============+=============+
| 0    | 1.5331e+05 | 100000 |      1 | 6.523e-01 | 0.0% |   0.0% 0.0% |     1.26 GB |     5.29 MB |
+------+------------+--------+--------+-----------+------+-------------+-------------+-------------+
Info: estimated rate : 8034.992845577424 (tsim/hr)                             [sph::Model][rank=0]

Usual cartesian rendering

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

Cylindrical rendering

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

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

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

Azymuthal rendering

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

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

Estimated memory usage: 1675 MB

Gallery generated by Sphinx-Gallery