Note
Go to the end to download the full example code.
CI test: sink state stays synchronized across MPI ranks through dump/reload#
Creates a tiny SPH setup with a few sinks and gas particles, checks sink sync, evolves, dumps, reloads into a fresh context, and checks again.
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 = 8.8e+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 = 1.1e+05 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 = 8.4e+04 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 = 2.9e+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 = 5.4e+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 = 5.4e+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 = 1.2e+05 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.8e+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.0e+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 = 3.9e+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.3e+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 = 5.7e+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.7e+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.2e+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 = 5.5e+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 = 9.9e+04 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 = 5.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.8e+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.6e+05 max = 2.3e+03) rate = 2.331000e+03 N.s^-1
SPH setup: the generation step took : 0.016854114 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.040355229 s
SPH setup: the setup took : 0.06288004800000001 s
Sinks are in sync !
---------------- 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.730061349693251
patch 1 high interf/patch volume: 1.5783972125435533
patch 2 high interf/patch volume: 1.5133333333333332
patch 3 high interf/patch volume: 1.3568773234200742
patch 4 high interf/patch volume: 1.5966101694915251
patch 5 high interf/patch volume: 1.4597701149425286
patch 6 high interf/patch volume: 1.4240282685512367
patch 7 high interf/patch volume: 1.2936507936507935
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.05500000000000001 unconverged cnt = 2273
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.9938650306748458
patch 1 high interf/patch volume: 1.8780487804878045
patch 2 high interf/patch volume: 1.893333333333333
patch 3 high interf/patch volume: 1.7806691449814123
patch 4 high interf/patch volume: 1.8474576271186438
patch 5 high interf/patch volume: 1.7432950191570882
patch 6 high interf/patch volume: 1.7703180212014136
patch 7 high interf/patch volume: 1.6825396825396823
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.06050000000000001 unconverged cnt = 2273
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.9938650306748458
patch 1 high interf/patch volume: 1.8780487804878045
patch 2 high interf/patch volume: 1.893333333333333
patch 3 high interf/patch volume: 1.7806691449814123
patch 4 high interf/patch volume: 1.8474576271186438
patch 5 high interf/patch volume: 1.7432950191570882
patch 6 high interf/patch volume: 1.7703180212014136
patch 7 high interf/patch volume: 1.6825396825396823
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.06655000000000001 unconverged cnt = 2273
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.128834355828221
patch 1 high interf/patch volume: 2.0313588850174216
patch 2 high interf/patch volume: 2.0733333333333333
patch 3 high interf/patch volume: 1.981412639405204
patch 4 high interf/patch volume: 1.9966101694915255
patch 5 high interf/patch volume: 1.9118773946360152
patch 6 high interf/patch volume: 1.9187279151943462
patch 7 high interf/patch volume: 1.8492063492063493
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.07320500000000002 unconverged cnt = 2273
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.4601226993865026
patch 1 high interf/patch volume: 2.35191637630662
patch 2 high interf/patch volume: 2.4400000000000004
patch 3 high interf/patch volume: 2.3457249070631976
patch 4 high interf/patch volume: 2.3491525423728814
patch 5 high interf/patch volume: 2.2643678160919545
patch 6 high interf/patch volume: 2.314487632508834
patch 7 high interf/patch volume: 2.246031746031746
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.08052550000000003 unconverged cnt = 2273
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.131901840490797
patch 1 high interf/patch volume: 2.9756097560975623
patch 2 high interf/patch volume: 3.11
patch 3 high interf/patch volume: 2.981412639405205
patch 4 high interf/patch volume: 3.0169491525423733
patch 5 high interf/patch volume: 2.8888888888888884
patch 6 high interf/patch volume: 2.9752650176678443
patch 7 high interf/patch volume: 2.8531746031746024
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.08857805000000003 unconverged cnt = 2273
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.4509202453987724
patch 1 high interf/patch volume: 3.337979094076655
patch 2 high interf/patch volume: 3.4599999999999995
patch 3 high interf/patch volume: 3.3717472118959115
patch 4 high interf/patch volume: 3.328813559322034
patch 5 high interf/patch volume: 3.2413793103448274
patch 6 high interf/patch volume: 3.31095406360424
patch 7 high interf/patch volume: 3.2301587301587302
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.09743585500000004 unconverged cnt = 2273
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.233128834355829
patch 1 high interf/patch volume: 4.142857142857143
patch 2 high interf/patch volume: 4.113333333333332
patch 3 high interf/patch volume: 4.040892193308551
patch 4 high interf/patch volume: 4.13220338983051
patch 5 high interf/patch volume: 4.057471264367816
patch 6 high interf/patch volume: 4.053003533568904
patch 7 high interf/patch volume: 3.9841269841269833
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.10717944050000006 unconverged cnt = 1086
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.469325153374233
patch 1 high interf/patch volume: 4.383275261324042
patch 2 high interf/patch volume: 4.346666666666668
patch 3 high interf/patch volume: 4.301115241635688
patch 4 high interf/patch volume: 4.345762711864408
patch 5 high interf/patch volume: 4.283524904214559
patch 6 high interf/patch volume: 4.300353356890459
patch 7 high interf/patch volume: 4.261904761904761
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.11789738455000007 unconverged cnt = 55
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.588957055214724
patch 1 high interf/patch volume: 5.508710801393732
patch 2 high interf/patch volume: 5.673333333333333
patch 3 high interf/patch volume: 5.639405204460967
patch 4 high interf/patch volume: 5.48135593220339
patch 5 high interf/patch volume: 5.421455938697317
patch 6 high interf/patch volume: 5.664310954063605
patch 7 high interf/patch volume: 5.642857142857144
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.12968712300500007 unconverged cnt = 2
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: 6.226993865030673
patch 1 high interf/patch volume: 6.177700348432057
patch 2 high interf/patch volume: 6.183333333333333
patch 3 high interf/patch volume: 6.130111524163569
patch 4 high interf/patch volume: 6.159322033898306
patch 5 high interf/patch volume: 6.095785440613028
patch 6 high interf/patch volume: 6.183745583038868
patch 7 high interf/patch volume: 6.11904761904762
---------------- t = 0, dt = 6.551773540016008e-06 ----------------
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: 6.2879999999999985
patch 1 high interf/patch volume: 6.2057761732852
patch 2 high interf/patch volume: 6.247148288973386
patch 3 high interf/patch volume: 6.180887372013654
patch 4 high interf/patch volume: 6.295774647887325
patch 5 high interf/patch volume: 6.1875
patch 6 high interf/patch volume: 6.247311827956988
patch 7 high interf/patch volume: 6.166123778501629
---------------- t = 6.551773540016008e-06, dt = 0.00022275667803016137 ----------------
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: 6.2879999999999985
patch 1 high interf/patch volume: 6.2057761732852
patch 2 high interf/patch volume: 6.247148288973386
patch 3 high interf/patch volume: 6.180887372013654
patch 4 high interf/patch volume: 6.394366197183099
patch 5 high interf/patch volume: 6.2875000000000005
patch 6 high interf/patch volume: 6.322580645161289
patch 7 high interf/patch volume: 6.2442996742671015
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.010106726402497119
---------------- t = 0.0002293084515701774, dt = 0.00018328014262681863 ----------------
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: 6.2879999999999985
patch 1 high interf/patch volume: 6.2057761732852
patch 2 high interf/patch volume: 6.247148288973386
patch 3 high interf/patch volume: 6.178082191780823
patch 4 high interf/patch volume: 6.394366197183099
patch 5 high interf/patch volume: 6.2875000000000005
patch 6 high interf/patch volume: 6.322580645161289
patch 7 high interf/patch volume: 6.2442996742671015
---------------- t = 0.00041258859419699605, dt = 0.00034020216026755444 ----------------
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: 6.2879999999999985
patch 1 high interf/patch volume: 6.200000000000001
patch 2 high interf/patch volume: 6.247148288973386
patch 3 high interf/patch volume: 6.178082191780823
patch 4 high interf/patch volume: 6.394366197183099
patch 5 high interf/patch volume: 6.2875000000000005
patch 6 high interf/patch volume: 6.322580645161289
patch 7 high interf/patch volume: 6.2442996742671015
Sinks are in sync !
Sinks are in sync !
---------------- t = 0.0007527907544645504, dt = 0.0004464158410420693 ----------------
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: 6.2799999999999985
patch 1 high interf/patch volume: 6.18909090909091
patch 2 high interf/patch volume: 6.247148288973386
patch 3 high interf/patch volume: 6.178082191780823
patch 4 high interf/patch volume: 6.382456140350876
patch 5 high interf/patch volume: 6.275
patch 6 high interf/patch volume: 6.323741007194244
patch 7 high interf/patch volume: 6.2442996742671015
---------------- t = 0.0011992065955066197, dt = 0.000512654960520611 ----------------
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: 6.191999999999998
patch 1 high interf/patch volume: 6.105454545454546
patch 2 high interf/patch volume: 6.235741444866922
patch 3 high interf/patch volume: 6.157534246575343
patch 4 high interf/patch volume: 6.407017543859648
patch 5 high interf/patch volume: 6.2906249999999995
patch 6 high interf/patch volume: 6.384892086330936
patch 7 high interf/patch volume: 6.302931596091206
---------------- t = 0.0017118615560272306, dt = 0.0005523343280162558 ----------------
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: 6.191999999999998
patch 1 high interf/patch volume: 6.050909090909092
patch 2 high interf/patch volume: 6.2015209125475295
patch 3 high interf/patch volume: 6.0479452054794525
patch 4 high interf/patch volume: 6.5069930069930075
patch 5 high interf/patch volume: 6.340624999999999
patch 6 high interf/patch volume: 6.458483754512637
patch 7 high interf/patch volume: 6.32899022801303
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.011840076446588807
---------------- t = 0.0022641958840434865, dt = 0.0002864278912784332 ----------------
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.992
patch 1 high interf/patch volume: 5.923636363636364
patch 2 high interf/patch volume: 5.935361216730037
patch 3 high interf/patch volume: 5.900343642611686
patch 4 high interf/patch volume: 6.157342657342658
patch 5 high interf/patch volume: 6.06875
patch 6 high interf/patch volume: 6.12274368231047
patch 7 high interf/patch volume: 6.061889250814333
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.02576228461305771
---------------- t = 0.0025506237753219196, dt = 0.00019897300678063906 ----------------
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: 6.016
patch 1 high interf/patch volume: 5.941818181818182
patch 2 high interf/patch volume: 5.984732824427481
patch 3 high interf/patch volume: 5.951557093425607
patch 4 high interf/patch volume: 6.234265734265736
patch 5 high interf/patch volume: 6.153124999999999
patch 6 high interf/patch volume: 6.198555956678701
patch 7 high interf/patch volume: 6.143322475570033
Sinks are in sync !
run_test_sink_synchro: OK
Current sinks:
[{'pos': (0.09877660057755726, -0.0003699820152028824, 0.000292858494405309), 'velocity': (-0.7909594697171893, -0.42835509420057705, 0.21579163220191747), 'sph_acceleration': (-2.788118606064703, 6.194080856496939, -3.386302543973059), 'ext_acceleration': (-371.6829482929875, -190.87247862152324, 85.76563352228862), 'mass': 1.0219999999999998, 'angular_momentum': (1.5355132242207976e-05, -3.4042791169362646e-05, -0.0013917210465080319), 'accretion_radius': 0.15}, {'pos': (-0.19723897477916802, 0.0992597386368148, 5.457513126197257e-05), 'velocity': (1.7840284522444005, -0.646223010107699, 0.03996740374797707), 'sph_acceleration': (40.13267477144062, -16.361891265380407, -1.2376352720798174), 'ext_acceleration': (456.6529453658143, -211.77602552231514, 16.18781052807078), 'mass': 0.523, 'angular_momentum': (1.4588956579124865e-08, -1.6035141770323956e-07, -0.0009268830439052644), 'accretion_radius': 0.15}, {'pos': (0.001919965581750201, -0.1466495904953547, 0.048938432104203365), 'velocity': (2.1365106397204268, 3.022034230508638, -0.9498166179408467), 'sph_acceleration': (4.433172833108257, 16.092200428397565, -6.362794011098259), 'ext_acceleration': (522.3351212189342, 1132.705683330991, -355.99519394799995), 'mass': 0.27, 'angular_momentum': (4.446934644992772e-05, 0.0006693090866850239, 0.0020383710492193246), 'accretion_radius': 0.15}]
check_ref_dataset: OK
Current sums: [np.float64(65.79475628688994), np.float64(-15.508903903707498), np.float64(-1.0412208430674998), np.float64(22613.455824802637), np.float64(-5.195701678016274), np.float64(15.008486573476972), np.float64(-19336.88835500463), np.float64(-2117.975619213131), np.float64(5826.038830234722), np.float64(245.85176557074223)]
9 import numpy as np
10
11 import shamrock
12
13 DUMP_NAME = "sink_sync_test.sham"
14
15
16 def check_sinks_are_in_sync(ctx, model):
17 s = str(model.get_sinks())
18 # Collective: every rank must call this
19 hist = shamrock.algs.all_string_histogram([s], delimiter="\n", hash_based=False)
20 if len(hist) != 1:
21 raise RuntimeError(f"sinks not in sync across ranks: {hist}")
22 key, count = next(iter(hist.items()))
23 if count != shamrock.sys.world_size():
24 raise RuntimeError(
25 f"expected count={shamrock.sys.world_size()}, got {count} for key={key!r}"
26 )
27 shamrock.sys.mpi_barrier()
28 if shamrock.sys.world_rank() == 0:
29 print("Sinks are in sync !")
30
31
32 si = shamrock.UnitSystem()
33 sicte = shamrock.Constants(si)
34 codeu = shamrock.UnitSystem(
35 unit_time=sicte.year(),
36 unit_length=sicte.au(),
37 unit_mass=sicte.sol_mass(),
38 )
39 ucte = shamrock.Constants(codeu)
40 G = ucte.G()
41
42
43 def build_model_with_sinks():
44 ctx = shamrock.Context()
45 ctx.pdata_layout_new()
46
47 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
48
49 cfg = model.gen_default_config()
50 cfg.set_self_gravity_none()
51 cfg.set_artif_viscosity_Constant(alpha_u=1.0, alpha_AV=1.0, beta_AV=2.0)
52 cfg.set_eos_isothermal(1.0)
53 cfg.set_particle_mass(1e-3)
54 cfg.set_boundary_periodic()
55 cfg.set_show_cfl_detail(True)
56 cfg.set_units(codeu)
57 model.set_solver_config(cfg)
58
59 model.set_cfl_cour(0.1)
60 model.set_cfl_force(0.1)
61 model.set_eta_sink(1.0)
62
63 model.init_scheduler(1000, 1)
64
65 # Very coarse HCP cube -> handful of SPH particles
66 dr = 0.05
67 bmin = (-0.6, -0.6, -0.6)
68 bmax = (0.6, 0.6, 0.6)
69 model.resize_simulation_box(bmin, bmax)
70
71 setup = model.get_setup()
72 gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
73 setup.apply_setup(gen)
74
75 eng = shamrock.algs.gen_seed(42)
76
77 def vel_func(r):
78 return (10.0, 0.0, 0.0)
79
80 model.set_field_value_lambda_f64_3("vxyz", vel_func)
81
82 # A few sinks (must be added after init_scheduler, on all ranks)
83 model.add_sink(1.0, (0.1, 0.0, 0.0), (0.0, 0.05, 0.0), 0.15)
84 model.add_sink(0.5, (-0.2, 0.1, 0.0), (0.0, -0.03, 0.0), 0.15)
85 model.add_sink(0.25, (0.0, -0.15, 0.05), (0.02, 0.0, 0.0), 0.15)
86
87 return ctx, model
88
89
90 def check_ref_dataset(sinks):
91 if shamrock.sys.world_rank() == 0:
92 print("Current sinks:")
93 print(sinks)
94
95 ref_dataset = [
96 {
97 "pos": (0.09877660057755726, -0.0003699820152028823, 0.00029285849440530886),
98 "velocity": (-0.7909594697171893, -0.42835509420057705, 0.2157916322019174),
99 "sph_acceleration": (-2.7881186060647103, 6.194080856496935, -3.386302543973038),
100 "ext_acceleration": (-371.6829482929875, -190.87247862152324, 85.76563352228864),
101 "mass": 1.0219999999999998,
102 "angular_momentum": (
103 1.5355132242207996e-05,
104 -3.404279116936264e-05,
105 -0.001391721046508032,
106 ),
107 "accretion_radius": 0.15,
108 },
109 {
110 "pos": (-0.19723897477916802, 0.0992597386368148, 5.457513126197254e-05),
111 "velocity": (1.7840284522444, -0.6462230101076988, 0.039967403747977054),
112 "sph_acceleration": (40.132674771440605, -16.361891265380414, -1.2376352720798245),
113 "ext_acceleration": (456.6529453658143, -211.77602552231514, 16.18781052807078),
114 "mass": 0.523,
115 "angular_momentum": (
116 1.4588956579124888e-08,
117 -1.6035141770323982e-07,
118 -0.0009268830439052644,
119 ),
120 "accretion_radius": 0.15,
121 },
122 {
123 "pos": (0.0019199655817502001, -0.1466495904953547, 0.048938432104203365),
124 "velocity": (2.1365106397204268, 3.0220342305086376, -0.9498166179408467),
125 "sph_acceleration": (4.433172833108257, 16.09220042839757, -6.362794011098254),
126 "ext_acceleration": (522.3351212189343, 1132.705683330991, -355.99519394799995),
127 "mass": 0.27,
128 "angular_momentum": (
129 4.4469346449927724e-05,
130 0.0006693090866850239,
131 0.0020383710492193246,
132 ),
133 "accretion_radius": 0.15,
134 },
135 ]
136
137 errors = []
138
139 if len(sinks) != len(ref_dataset):
140 errors.append(f"sink count mismatch: got {len(sinks)}, expected {len(ref_dataset)}")
141 else:
142 for i, (got_sink, ref_sink) in enumerate(zip(sinks, ref_dataset)):
143 for key, ref_val in ref_sink.items():
144 got_val = got_sink[key]
145 got_arr = np.asarray(got_val, dtype=float)
146 ref_arr = np.asarray(ref_val, dtype=float)
147 rtol = 1e-14 if key == "sph_acceleration" else 1e-15
148 if not np.all(np.isclose(got_arr, ref_arr, rtol=rtol, atol=1e-18)):
149 abs_diff = np.abs(got_arr - ref_arr)
150 with np.errstate(divide="ignore", invalid="ignore"):
151 rel_diff = np.where(ref_arr != 0, abs_diff / np.abs(ref_arr), abs_diff)
152 errors.append(
153 f"sink[{i}].{key} mismatch:\n"
154 f" got={got_val}\n"
155 f" ref={ref_val}\n"
156 f" max abs diff={np.max(abs_diff)}\n"
157 f" max rel diff={np.max(rel_diff)}"
158 )
159
160 for err in errors:
161 print(err)
162
163 if errors:
164 raise RuntimeError(f"check_ref_dataset failed with {len(errors)} error(s)")
165
166 if shamrock.sys.world_rank() == 0:
167 print("check_ref_dataset: OK")
168
169
170 def main():
171 ctx, model = build_model_with_sinks()
172
173 check_sinks_are_in_sync(ctx, model)
174
175 for _ in range(5):
176 model.timestep()
177 check_sinks_are_in_sync(ctx, model)
178
179 sinks_before_dump = str(model.get_sinks())
180 model.dump(DUMP_NAME)
181
182 del model
183 del ctx
184
185 ctx2 = shamrock.Context()
186 ctx2.pdata_layout_new()
187 model2 = shamrock.get_Model_SPH(context=ctx2, vector_type="f64_3", sph_kernel="M4")
188 model2.load_from_dump(DUMP_NAME)
189
190 sinks_after_reload = str(model2.get_sinks())
191 if sinks_before_dump != sinks_after_reload:
192 raise RuntimeError(
193 "sink content changed across dump/reload:\n"
194 f" before={sinks_before_dump!r}\n"
195 f" after ={sinks_after_reload!r}"
196 )
197
198 check_sinks_are_in_sync(ctx2, model2)
199
200 for _ in range(5):
201 model2.timestep()
202 check_sinks_are_in_sync(ctx2, model2)
203
204 if shamrock.sys.world_rank() == 0:
205 print("run_test_sink_synchro: OK")
206
207 check_ref_dataset(model2.get_sinks())
208
209 dic = ctx2.collect_data()
210
211 if shamrock.sys.world_rank() > 0:
212 return
213
214 assert 2266 == len(dic["xyz"])
215
216 sum_pos = np.sum(dic["xyz"], axis=0)
217 sum_vel = np.sum(dic["vxyz"], axis=0)
218 sum_acc = np.sum(dic["axyz"], axis=0)
219 sum_hpart = np.sum(dic["hpart"], axis=0)
220
221 dat = np.concatenate([sum_pos, sum_vel, sum_acc, np.atleast_1d(sum_hpart)])
222 print("Current sums: ", [dat[i] for i in range(len(dat))])
223
224 ref_sums = [
225 65.79475628688992,
226 -15.508903903707498,
227 -1.0412208430674998,
228 22613.455824802637,
229 -5.195701678016285,
230 15.008486573476965,
231 -19336.888355004605,
232 -2117.975619213131,
233 5826.038830234707,
234 245.85176557074223,
235 ]
236
237 mismatch = False
238 for i in range(len(dat)):
239 if not np.isclose(dat[i], ref_sums[i], rtol=1e-12, atol=1e-18):
240 abs_diff = np.abs(dat[i] - ref_sums[i])
241 rel_diff = abs_diff / np.abs(ref_sums[i])
242 print(f"sum[{i}] mismatch: got {dat[i]}, expected {ref_sums[i]}")
243 print(f" max abs diff={np.max(abs_diff)}")
244 print(f" max rel diff={np.max(rel_diff)}")
245 mismatch = True
246 if mismatch:
247 raise RuntimeError("sums mismatch")
248
249
250 main()
Total running time of the script: (0 minutes 4.424 seconds)
Estimated memory usage: 160 MB