Note
Go to the end to download the full example code.
CI test: external forces point mass accretion#
Creates a tiny SPH setup with a few point mass external forces with an accretion radius, evolves, dumps, reloads into a fresh context, evolves again, and checks the resulting particle count and field sums.
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.3e+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 = 7.2e+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.2e+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.6e+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.1e+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.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 = 4.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 = 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 = 5.1e+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 = 3.7e+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.8e+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.5e+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.2e+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 = 4.6e+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 = 4.7e+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.2e+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.0e+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.0e+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.3e+05 max = 2.3e+03) rate = 2.331000e+03 N.s^-1
SPH setup: the generation step took : 0.017440261000000002 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.043354654 s
SPH setup: the setup took : 0.06660589800000001 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.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 = 8.357077848469307e-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 = 8.357077848469307e-06, dt = 0.0002841410065758001 ----------------
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
---------------- t = 0.0002924980844242694, dt = 0.00046793507433955827 ----------------
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.0007604331587638276, dt = 0.0005899254064286999 ----------------
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.2719999999999985
patch 1 high interf/patch volume: 6.174545454545456
patch 2 high interf/patch volume: 6.247148288973386
patch 3 high interf/patch volume: 6.167808219178084
patch 4 high interf/patch volume: 6.376760563380283
patch 5 high interf/patch volume: 6.250000000000001
patch 6 high interf/patch volume: 6.322580645161289
patch 7 high interf/patch volume: 6.234527687296417
---------------- t = 0.0013503585651925276, dt = 0.000669882067954046 ----------------
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.107999999999998
patch 1 high interf/patch volume: 6.014545454545455
patch 2 high interf/patch volume: 6.152091254752853
patch 3 high interf/patch volume: 6.065068493150685
patch 4 high interf/patch volume: 6.507042253521128
patch 5 high interf/patch volume: 6.378125
patch 6 high interf/patch volume: 6.4587813620071675
patch 7 high interf/patch volume: 6.364820846905538
---------------- t = 0.0020202406331465738, dt = 0.0007209398751302408 ----------------
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.192
patch 1 high interf/patch volume: 5.98909090909091
patch 2 high interf/patch volume: 6.114068441064637
patch 3 high interf/patch volume: 5.931506849315069
patch 4 high interf/patch volume: 6.515789473684208
patch 5 high interf/patch volume: 6.296874999999999
patch 6 high interf/patch volume: 6.420863309352518
patch 7 high interf/patch volume: 6.247557003257329
---------------- t = 0.0027411805082768145, dt = 0.0007299132315052103 ----------------
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.952
patch 1 high interf/patch volume: 5.890510948905109
patch 2 high interf/patch volume: 5.896946564885496
patch 3 high interf/patch volume: 5.855172413793104
patch 4 high interf/patch volume: 6.007017543859647
patch 5 high interf/patch volume: 5.928125
patch 6 high interf/patch volume: 5.928057553956834
patch 7 high interf/patch volume: 5.892508143322476
---------------- t = 0.003471093739782025, dt = 0.0006701883402267496 ----------------
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.112
patch 1 high interf/patch volume: 6.032727272727273
patch 2 high interf/patch volume: 6.125954198473283
patch 3 high interf/patch volume: 6.076124567474049
patch 4 high interf/patch volume: 6.224561403508769
patch 5 high interf/patch volume: 6.153124999999999
patch 6 high interf/patch volume: 6.220216606498196
patch 7 high interf/patch volume: 6.1830065359477135
---------------- t = 0.004141282080008774, dt = 0.0005932073639326702 ----------------
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.092369477911648
patch 1 high interf/patch volume: 6.047101449275361
patch 2 high interf/patch volume: 6.091954022988505
patch 3 high interf/patch volume: 6.097560975609758
patch 4 high interf/patch volume: 6.125435540069687
patch 5 high interf/patch volume: 6.093749999999999
patch 6 high interf/patch volume: 6.144927536231886
patch 7 high interf/patch volume: 6.124183006535948
run_test_ext_forces_accretion: OK
Current sums: [111.13342631354446, -15.790263346174845, -1.0024165122154156, 22523.573252726077, 0.16482394866530714, 25.96527810306144, -25326.18011836384, -69.26477834526008, 6203.398030417404, 245.61542121817263]
10 import numpy as np
11
12 import shamrock
13
14 DUMP_NAME = "ext_forces_accretion_test.sham"
15
16 si = shamrock.UnitSystem()
17 sicte = shamrock.Constants(si)
18 codeu = shamrock.UnitSystem(
19 unit_time=sicte.year(),
20 unit_length=sicte.au(),
21 unit_mass=sicte.sol_mass(),
22 )
23 ucte = shamrock.Constants(codeu)
24 G = ucte.G()
25
26
27 def build_model_with_ext_forces():
28 ctx = shamrock.Context()
29 ctx.pdata_layout_new()
30
31 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
32
33 cfg = model.gen_default_config()
34 cfg.set_self_gravity_none()
35 cfg.set_artif_viscosity_Constant(alpha_u=1.0, alpha_AV=1.0, beta_AV=2.0)
36 cfg.set_eos_isothermal(1.0)
37 cfg.set_particle_mass(1e-3)
38 cfg.set_boundary_periodic()
39 cfg.set_show_cfl_detail(True)
40 cfg.set_units(codeu)
41
42 cfg.add_ext_force_paczynski_wiita(1.0, (0.1, 0.0, 0.0), 0.15)
43 cfg.add_ext_force_paczynski_wiita(0.5, (-0.2, 0.1, 0.0), 0.15)
44 cfg.add_ext_force_paczynski_wiita(0.25, (0.0, -0.15, 0.05), 0.15)
45
46 model.set_solver_config(cfg)
47
48 model.set_cfl_cour(0.1)
49 model.set_cfl_force(0.1)
50
51 model.init_scheduler(1000, 1)
52
53 dr = 0.05
54 bmin = (-0.6, -0.6, -0.6)
55 bmax = (0.6, 0.6, 0.6)
56 model.resize_simulation_box(bmin, bmax)
57
58 setup = model.get_setup()
59 gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
60 setup.apply_setup(gen)
61
62 def vel_func(r):
63 return (10.0, 0.0, 0.0)
64
65 model.set_field_value_lambda_f64_3("vxyz", vel_func)
66
67 return ctx, model
68
69
70 def main():
71 ctx, model = build_model_with_ext_forces()
72
73 for _ in range(5):
74 model.timestep()
75
76 model.dump(DUMP_NAME)
77
78 del model
79 del ctx
80
81 ctx2 = shamrock.Context()
82 ctx2.pdata_layout_new()
83 model2 = shamrock.get_Model_SPH(context=ctx2, vector_type="f64_3", sph_kernel="M4")
84 model2.load_from_dump(DUMP_NAME)
85
86 for _ in range(5):
87 model2.timestep()
88
89 if shamrock.sys.world_rank() == 0:
90 print("run_test_ext_forces_accretion: OK")
91
92 dic = ctx2.collect_data()
93
94 if shamrock.sys.world_rank() > 0:
95 return
96
97 assert 2262 == len(dic["xyz"])
98
99 sum_pos = np.sum(dic["xyz"], axis=0)
100 sum_vel = np.sum(dic["vxyz"], axis=0)
101 sum_acc = np.sum(dic["axyz"], axis=0)
102 sum_hpart = np.sum(dic["hpart"], axis=0)
103
104 dat = np.concatenate([sum_pos, sum_vel, sum_acc, np.atleast_1d(sum_hpart)])
105 print("Current sums: ", [float(dat[i]) for i in range(len(dat))])
106
107 ref_sums = [
108 111.13342631354446,
109 -15.790263346174845,
110 -1.0024165122154156,
111 22523.573252726077,
112 0.16482394866530714,
113 25.96527810306144,
114 -25326.18011836384,
115 -69.26477834526008,
116 6203.398030417404,
117 245.61542121817263,
118 ]
119
120 mismatch = False
121 for i in range(len(dat)):
122 if not np.isclose(dat[i], ref_sums[i], rtol=1e-11, atol=1e-18):
123 abs_diff = np.abs(dat[i] - ref_sums[i])
124 rel_diff = abs_diff / np.abs(ref_sums[i])
125 print(f"sum[{i}] mismatch: got {dat[i]}, expected {ref_sums[i]}")
126 print(f" max abs diff={np.max(abs_diff)}")
127 print(f" max rel diff={np.max(rel_diff)}")
128 mismatch = True
129 if mismatch:
130 raise RuntimeError("sums mismatch")
131
132
133 main()
Total running time of the script: (0 minutes 4.227 seconds)
Estimated memory usage: 161 MB