Note
Go to the end to download the full example code.
CI test: rendering a precomputed Field matches rendering by name#
Creates a tiny SPH gas setup, then for each render entry point (render_slice, render_column_integ, render_azymuthal_integ, render_cartesian_slice, render_cartesian_column_integ) checks that rendering a Field object obtained from model.compute_field(…) gives the exact same result as rendering by field name directly. Also checks that a Field derived with shamrock.map_fields_f64 renders correctly.
SPH setup: generating particles ...
SPH setup: Nstep = 0 ( 0.0e+00 ) Ntotal = 0 ( 0.0e+00 rank min = 0.0e+00 max = 0.0e+00) rate = 0.000000e+00 N.s^-1
SPH setup: Nstep = 70 ( 7.0e+01 ) Ntotal = 70 ( 7.0e+01 rank min = 2.2e+04 max = 7.0e+01) rate = 7.000000e+01 N.s^-1
SPH setup: Nstep = 106 ( 1.1e+02 ) Ntotal = 176 ( 1.8e+02 rank min = 9.4e+04 max = 1.8e+02) rate = 1.760000e+02 N.s^-1
SPH setup: Nstep = 118 ( 1.2e+02 ) Ntotal = 294 ( 2.9e+02 rank min = 1.1e+05 max = 2.9e+02) rate = 2.940000e+02 N.s^-1
SPH setup: Nstep = 109 ( 1.1e+02 ) Ntotal = 403 ( 4.0e+02 rank min = 3.4e+05 max = 4.0e+02) rate = 4.030000e+02 N.s^-1
SPH setup: Nstep = 110 ( 1.1e+02 ) Ntotal = 513 ( 5.1e+02 rank min = 1.2e+05 max = 5.1e+02) rate = 5.130000e+02 N.s^-1
SPH setup: Nstep = 106 ( 1.1e+02 ) Ntotal = 619 ( 6.2e+02 rank min = 4.8e+05 max = 6.2e+02) rate = 6.190000e+02 N.s^-1
SPH setup: Nstep = 118 ( 1.2e+02 ) Ntotal = 737 ( 7.4e+02 rank min = 6.1e+05 max = 7.4e+02) rate = 7.370000e+02 N.s^-1
SPH setup: Nstep = 110 ( 1.1e+02 ) Ntotal = 847 ( 8.5e+02 rank min = 8.1e+04 max = 8.5e+02) rate = 8.470000e+02 N.s^-1
SPH setup: Nstep = 109 ( 1.1e+02 ) Ntotal = 956 ( 9.6e+02 rank min = 4.7e+05 max = 9.6e+02) rate = 9.560000e+02 N.s^-1
SPH setup: Nstep = 115 ( 1.2e+02 ) Ntotal = 1071 ( 1.1e+03 rank min = 4.6e+05 max = 1.1e+03) rate = 1.071000e+03 N.s^-1
SPH setup: Nstep = 121 ( 1.2e+02 ) Ntotal = 1192 ( 1.2e+03 rank min = 4.2e+05 max = 1.2e+03) rate = 1.192000e+03 N.s^-1
SPH setup: Nstep = 118 ( 1.2e+02 ) Ntotal = 1310 ( 1.3e+03 rank min = 1.1e+05 max = 1.3e+03) rate = 1.310000e+03 N.s^-1
SPH setup: Nstep = 119 ( 1.2e+02 ) Ntotal = 1429 ( 1.4e+03 rank min = 5.9e+05 max = 1.4e+03) rate = 1.429000e+03 N.s^-1
SPH setup: Nstep = 121 ( 1.2e+02 ) Ntotal = 1550 ( 1.6e+03 rank min = 6.6e+05 max = 1.6e+03) rate = 1.550000e+03 N.s^-1
SPH setup: Nstep = 121 ( 1.2e+02 ) Ntotal = 1671 ( 1.7e+03 rank min = 5.8e+05 max = 1.7e+03) rate = 1.671000e+03 N.s^-1
SPH setup: Nstep = 118 ( 1.2e+02 ) Ntotal = 1789 ( 1.8e+03 rank min = 5.8e+05 max = 1.8e+03) rate = 1.789000e+03 N.s^-1
SPH setup: Nstep = 119 ( 1.2e+02 ) Ntotal = 1908 ( 1.9e+03 rank min = 6.8e+05 max = 1.9e+03) rate = 1.908000e+03 N.s^-1
SPH setup: Nstep = 122 ( 1.2e+02 ) Ntotal = 2030 ( 2.0e+03 rank min = 1.0e+05 max = 2.0e+03) rate = 2.030000e+03 N.s^-1
SPH setup: Nstep = 120 ( 1.2e+02 ) Ntotal = 2150 ( 2.2e+03 rank min = 6.4e+05 max = 2.2e+03) rate = 2.150000e+03 N.s^-1
SPH setup: Nstep = 118 ( 1.2e+02 ) Ntotal = 2268 ( 2.3e+03 rank min = 5.9e+05 max = 2.3e+03) rate = 2.268000e+03 N.s^-1
SPH setup: Nstep = 63 ( 6.3e+01 ) Ntotal = 2331 ( 2.3e+03 rank min = 2.5e+05 max = 2.3e+03) rate = 2.331000e+03 N.s^-1
SPH setup: the generation step took : 0.023174966 s
SPH setup: final particle count = 2331 beginning injection ...
SPH setup: injected 2000 / 2331 => 85.8% | ranks with patchs = 1 / 1 -> local loop <-
SPH setup: injected 2331 / 2331 => 100.0% | ranks with patchs = 1 / 1 <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.042585828000000006 s
SPH setup: the setup took : 0.071742233 s
---------------- t = 0, dt = 0 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 1.8154761904761907
patch 1 high interf/patch volume: 1.6326530612244896
patch 2 high interf/patch volume: 1.5320512820512822
patch 3 high interf/patch volume: 1.3695652173913047
patch 4 high interf/patch volume: 1.677631578947368
patch 5 high interf/patch volume: 1.5037593984962407
patch 6 high interf/patch volume: 1.4791666666666665
patch 7 high interf/patch volume: 1.3294117647058825
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.05500000000000001 unconverged cnt = 2331
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 2.071428571428573
patch 1 high interf/patch volume: 1.9251700680272104
patch 2 high interf/patch volume: 1.9230769230769234
patch 3 high interf/patch volume: 1.811594202898551
patch 4 high interf/patch volume: 1.9210526315789467
patch 5 high interf/patch volume: 1.781954887218045
patch 6 high interf/patch volume: 1.8333333333333328
patch 7 high interf/patch volume: 1.729411764705882
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.06050000000000001 unconverged cnt = 2331
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 2.071428571428573
patch 1 high interf/patch volume: 1.9251700680272104
patch 2 high interf/patch volume: 1.9230769230769234
patch 3 high interf/patch volume: 1.811594202898551
patch 4 high interf/patch volume: 1.9210526315789467
patch 5 high interf/patch volume: 1.781954887218045
patch 6 high interf/patch volume: 1.8333333333333328
patch 7 high interf/patch volume: 1.729411764705882
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.06655000000000001 unconverged cnt = 2331
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 2.2142857142857144
patch 1 high interf/patch volume: 2.0884353741496593
patch 2 high interf/patch volume: 2.096153846153846
patch 3 high interf/patch volume: 2.0072463768115942
patch 4 high interf/patch volume: 2.0789473684210527
patch 5 high interf/patch volume: 1.9624060150375937
patch 6 high interf/patch volume: 1.9791666666666663
patch 7 high interf/patch volume: 1.8941176470588235
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.07320500000000002 unconverged cnt = 2331
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 2.535714285714285
patch 1 high interf/patch volume: 2.4013605442176864
patch 2 high interf/patch volume: 2.48076923076923
patch 3 high interf/patch volume: 2.384057971014493
patch 4 high interf/patch volume: 2.4342105263157885
patch 5 high interf/patch volume: 2.3082706766917287
patch 6 high interf/patch volume: 2.395833333333333
patch 7 high interf/patch volume: 2.3019607843137253
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.08052550000000003 unconverged cnt = 2331
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 3.211309523809524
patch 1 high interf/patch volume: 3.0442176870748296
patch 2 high interf/patch volume: 3.1506410256410255
patch 3 high interf/patch volume: 3.003623188405797
patch 4 high interf/patch volume: 3.095394736842104
patch 5 high interf/patch volume: 2.936090225563909
patch 6 high interf/patch volume: 3.0451388888888893
patch 7 high interf/patch volume: 2.901960784313724
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.08857805000000003 unconverged cnt = 2331
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 3.520833333333334
patch 1 high interf/patch volume: 3.397959183673469
patch 2 high interf/patch volume: 3.4871794871794872
patch 3 high interf/patch volume: 3.3840579710144927
patch 4 high interf/patch volume: 3.398026315789473
patch 5 high interf/patch volume: 3.2819548872180446
patch 6 high interf/patch volume: 3.3750000000000004
patch 7 high interf/patch volume: 3.2745098039215677
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.09743585500000004 unconverged cnt = 2331
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 4.291666666666666
patch 1 high interf/patch volume: 4.183673469387755
patch 2 high interf/patch volume: 4.160256410256409
patch 3 high interf/patch volume: 4.072463768115942
patch 4 high interf/patch volume: 4.203947368421052
patch 5 high interf/patch volume: 4.097744360902255
patch 6 high interf/patch volume: 4.104166666666667
patch 7 high interf/patch volume: 4.019607843137255
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.10717944050000006 unconverged cnt = 856
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 4.520833333333333
patch 1 high interf/patch volume: 4.445578231292518
patch 2 high interf/patch volume: 4.384615384615383
patch 3 high interf/patch volume: 4.32608695652174
patch 4 high interf/patch volume: 4.411184210526315
patch 5 high interf/patch volume: 4.334586466165414
patch 6 high interf/patch volume: 4.347222222222222
patch 7 high interf/patch volume: 4.294117647058823
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.11789738455000007 unconverged cnt = 36
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
patch 0 high interf/patch volume: 5.607142857142858
patch 1 high interf/patch volume: 5.544217687074831
patch 2 high interf/patch volume: 5.698717948717951
patch 3 high interf/patch volume: 5.659420289855074
patch 4 high interf/patch volume: 5.5131578947368425
patch 5 high interf/patch volume: 5.451127819548869
patch 6 high interf/patch volume: 5.687500000000001
patch 7 high interf/patch volume: 5.658823529411766
Warning: the unit system is not set [sph::Config][rank=0]
Warning: the unit system is not set [sph::Config][rank=0]
Warning: the unit system is not set [sph::Config][rank=0]
Warning: the unit system is not set [sph::Config][rank=0]
Warning: the unit system is not set [sph::Config][rank=0]
Warning: the unit system is not set [sph::Config][rank=0]
Warning: the unit system is not set [sph::Config][rank=0]
Warning: the unit system is not set [sph::Config][rank=0]
render_slice(rho): OK
render_column_integ(rho): OK
render_azymuthal_integ(rho): OK
render_cartesian_slice(rho): OK
render_cartesian_column_integ(rho): OK
render_slice(vxyz): OK
render_cartesian_slice(vxyz): OK
render_slice(map_fields_f64 derived field): OK
run_sph_render_compute_field: OK
13 import numpy as np
14
15 import shamrock
16
17
18 def build_model():
19 ctx = shamrock.Context()
20 ctx.pdata_layout_new()
21
22 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
23
24 cfg = model.gen_default_config()
25 cfg.set_self_gravity_none()
26 cfg.set_artif_viscosity_Constant(alpha_u=1.0, alpha_AV=1.0, beta_AV=2.0)
27 cfg.set_eos_isothermal(1.0)
28 cfg.set_particle_mass(1e-3)
29 cfg.set_boundary_periodic()
30 model.set_solver_config(cfg)
31
32 model.set_cfl_cour(0.1)
33 model.set_cfl_force(0.1)
34
35 model.init_scheduler(1000, 1)
36
37 dr = 0.05
38 bmin = (-0.6, -0.6, -0.6)
39 bmax = (0.6, 0.6, 0.6)
40 model.resize_simulation_box(bmin, bmax)
41
42 setup = model.get_setup()
43 gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
44 setup.apply_setup(gen)
45
46 def vel_func(r):
47 x, y, z = r
48 return (0.1 * y, -0.1 * x, 0.0)
49
50 model.set_field_value_lambda_f64_3("vxyz", vel_func)
51
52 model.timestep()
53
54 return ctx, model
55
56
57 def check_equal(what, ref, got):
58 ref = np.asarray(ref)
59 got = np.asarray(got)
60 if ref.shape != got.shape:
61 raise RuntimeError(f"{what}: shape mismatch, ref={ref.shape}, got={got.shape}")
62 if not np.array_equal(ref, got):
63 abs_diff = np.abs(ref - got)
64 raise RuntimeError(
65 f"{what}: field-based render does not match name-based render\n"
66 f" max abs diff={np.max(abs_diff)}"
67 )
68 if shamrock.sys.world_rank() == 0:
69 print(f"{what}: OK")
70
71
72 def main():
73 ctx, model = build_model()
74
75 positions = [
76 (0.0, 0.0, 0.0),
77 (0.1, 0.0, 0.0),
78 (0.0, 0.1, 0.0),
79 (0.2, 0.1, 0.0),
80 (-0.15, -0.1, 0.05),
81 ]
82
83 rays = [shamrock.math.Ray_f64_3(pos, (0.0, 0.0, 1.0)) for pos in positions]
84
85 ring_rays = [
86 shamrock.math.RingRay_f64_3((0.0, 0.0, z), r, (1.0, 0.0, 0.0), (0.0, 1.0, 0.0))
87 for z, r in [(-0.1, 0.05), (0.0, 0.1), (0.1, 0.15)]
88 ]
89
90 center = (0.0, 0.0, 0.0)
91 delta_x = (0.5, 0.0, 0.0)
92 delta_y = (0.0, 0.5, 0.0)
93 nx = ny = 8
94
95 # ---- scalar field (rho, Field_f64) ----
96 rho_field = model.compute_field("rho", "f64")
97
98 check_equal(
99 "render_slice(rho)",
100 model.render_slice("rho", "f64", positions),
101 model.render_slice(rho_field, positions),
102 )
103
104 check_equal(
105 "render_column_integ(rho)",
106 model.render_column_integ("rho", "f64", rays),
107 model.render_column_integ(rho_field, rays),
108 )
109
110 check_equal(
111 "render_azymuthal_integ(rho)",
112 model.render_azymuthal_integ("rho", "f64", ring_rays),
113 model.render_azymuthal_integ(rho_field, ring_rays),
114 )
115
116 check_equal(
117 "render_cartesian_slice(rho)",
118 model.render_cartesian_slice(
119 "rho", "f64", center=center, delta_x=delta_x, delta_y=delta_y, nx=nx, ny=ny
120 ),
121 model.render_cartesian_slice(
122 rho_field, center=center, delta_x=delta_x, delta_y=delta_y, nx=nx, ny=ny
123 ),
124 )
125
126 check_equal(
127 "render_cartesian_column_integ(rho)",
128 model.render_cartesian_column_integ(
129 "rho", "f64", center=center, delta_x=delta_x, delta_y=delta_y, nx=nx, ny=ny
130 ),
131 model.render_cartesian_column_integ(
132 rho_field, center=center, delta_x=delta_x, delta_y=delta_y, nx=nx, ny=ny
133 ),
134 )
135
136 # ---- vector field (vxyz, Field_f64_3) ----
137 vxyz_field = model.compute_field("vxyz", "f64_3")
138
139 check_equal(
140 "render_slice(vxyz)",
141 model.render_slice("vxyz", "f64_3", positions),
142 model.render_slice(vxyz_field, positions),
143 )
144
145 check_equal(
146 "render_cartesian_slice(vxyz)",
147 model.render_cartesian_slice(
148 "vxyz", "f64_3", center=center, delta_x=delta_x, delta_y=delta_y, nx=nx, ny=ny
149 ),
150 model.render_cartesian_slice(
151 vxyz_field, center=center, delta_x=delta_x, delta_y=delta_y, nx=nx, ny=ny
152 ),
153 )
154
155 # ---- derived field via shamrock.map_fields_f64 ----
156 def scale_by_two(size, x):
157 return 2.0 * x
158
159 rho_x2_field = shamrock.map_fields_f64(scale_by_two, x=rho_field)
160
161 base = np.asarray(model.render_slice(rho_field, positions))
162 derived = np.asarray(model.render_slice(rho_x2_field, positions))
163 if not np.allclose(derived, 2.0 * base, rtol=1e-12, atol=1e-18):
164 raise RuntimeError(
165 "render_slice(map_fields_f64(2*rho)) does not match 2*render_slice(rho)\n"
166 f" base={base}\n"
167 f" derived={derived}"
168 )
169 if shamrock.sys.world_rank() == 0:
170 print("render_slice(map_fields_f64 derived field): OK")
171
172 if shamrock.sys.world_rank() == 0:
173 print("run_sph_render_compute_field: OK")
174
175
176 main()
Total running time of the script: (0 minutes 2.636 seconds)
Estimated memory usage: 160 MB