Note
Go to the end to download the full example code.
Testing the MHD rotor problem (3D) with SPMHD#
CI test for the 3D MHD rotor problem (Balsara & Spicer 1999; Toth 2000’s ‘first rotor problem’), following the setup described in the Phantom code paper (Price et al. 2018, Section 5.6.4): a rapidly rotating dense disc of material immersed in a uniform ambient medium threaded by a uniform magnetic field, testing the propagation of rotational discontinuities.
13 import os
14
15 import matplotlib.pyplot as plt
16 import numpy as np
17
18 import shamrock
19
20 # If we use the shamrock executable to run this script instead of the python interpreter,
21 # we should not initialize the system as the shamrock executable needs to handle specific MPI logic
22 if not shamrock.sys.is_initialized():
23 shamrock.change_loglevel(1)
24 shamrock.sys.init("0:0")
27 kernel = "M4"
28
29 C_cour = 0.3
30 C_force = 0.25
31
32 nx = 32
33 nz_thin = 12
34
35 gamma = 1.4
36
37 R_disc = 0.1
38 rho_disc = 10.0
39 rho_ambient = 1.0
40 P0 = 1.0
41 v0 = 2.0
42 Bx0 = 5.0 / np.sqrt(4.0 * np.pi)
43
44 fact = (rho_disc / rho_ambient) ** (1.0 / 3.0)
45
46 u_disc = P0 / ((gamma - 1) * rho_disc)
47 u_ambient = P0 / ((gamma - 1) * rho_ambient)
48
49 t_target = 0.15
53 ctx = shamrock.Context()
54 ctx.pdata_layout_new()
55
56 si = shamrock.UnitSystem()
57 codeu = shamrock.UnitSystem(unit_time=1.0, unit_length=1.0, unit_mass=1.2566370621219e-06)
58 ucte = shamrock.Constants(codeu)
59 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel=kernel)
60
61 cfg = model.gen_default_config()
62 cfg.set_units(codeu)
63
64 mu_0 = ucte.mu_0()
65
66 cfg.set_artif_viscosity_None()
67 cfg.set_IdealMHD(sigma_mhd=1, sigma_u=1)
68 cfg.set_boundary_periodic()
69 cfg.set_eos_adiabatic(gamma)
70 model.set_solver_config(cfg)
71 cfg.print_status()
72
73 crit_split = int(1e7)
74 crit_merge = 1
75 model.init_scheduler(crit_split, crit_merge)
----- SPH Solver configuration -----
[
{
"artif_viscosity": {
"type": "none"
},
"boundary_config": {
"bc_type": "periodic"
},
"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",
"eos_type": "adiabatic",
"gamma": 1.4
},
"epsilon_h": 1e-06,
"ext_force_config": {
"force_list": []
},
"gpart_mass": 0.0,
"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": {
"alpha_AV": 1.0,
"alpha_B": 1.0,
"alpha_u": 1.0,
"beta_AV": 1.0,
"mhd_type": "ideal_mhd_constrained_hyper_para",
"sigma_mhd": 1.0
},
"neigh_cache_strategy": "two_stage",
"particle_killing": [],
"particle_reordering_step_freq": 1000,
"save_dt_to_fields": false,
"scheduler_config": {
"merge_load_value": 0,
"split_load_value": 0
},
"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": 3,
"type_id": "sycl::vec<f64,3>",
"unit_sys": {
"unit_current": 1.0,
"unit_length": 1.0,
"unit_lumint": 1.0,
"unit_mass": 1.2566370621219e-06,
"unit_qte": 1.0,
"unit_temperature": 1.0,
"unit_time": 1.0
}
}
]
------------------------------------
79 (xs, ys, zs) = model.get_box_dim_fcc_3d(1, nx, nx, nz_thin)
80 dr = 1.0 / xs
81 (xs, ys, zs) = model.get_box_dim_fcc_3d(dr, nx, nx, nz_thin)
82
83 box_min = (-xs / 2, -ys / 2, -zs / 2)
84 box_max = (xs / 2, ys / 2, zs / 2)
85
86 model.resize_simulation_box(box_min, box_max)
87
88 dr_disc = dr / fact
89 disc_box_min = (-R_disc, -R_disc, -zs / 2)
90 disc_box_max = (R_disc, R_disc, zs / 2)
91
92
93 def in_disc(r):
94 x, y, z = r
95 return (x * x + y * y) < R_disc * R_disc
96
97
98 def out_disc(r):
99 x, y, z = r
100 return (x * x + y * y) >= R_disc * R_disc
101
102
103 setup = model.get_setup()
104 gen_ambient = setup.make_generator_lattice_hcp(dr, box_min, box_max)
105 gen_ambient_masked = setup.make_modifier_filter(parent=gen_ambient, filter=out_disc)
106 gen_disc = setup.make_generator_lattice_hcp(dr_disc, disc_box_min, disc_box_max)
107 gen_disc_masked = setup.make_modifier_filter(parent=gen_disc, filter=in_disc)
108 comb = setup.make_combiner_add(gen_ambient_masked, gen_disc_masked)
109 setup.apply_setup(comb)
110
111 # Internal energy: ambient value everywhere, then overwrite inside the disc
112 model.set_value_in_a_box("uint", "f64", u_ambient, box_min, box_max)
113 model.set_value_in_sphere("uint", "f64", u_disc, (0, 0, 0), R_disc)
114
115
116 def vel_func(r):
117 x, y, z = r
118 r2 = x * x + y * y
119 if r2 >= R_disc * R_disc:
120 return (0.0, 0.0, 0.0)
121 rcyl = np.sqrt(max(r2, 1e-30))
122 vx = -v0 * y / rcyl
123 vy = v0 * x / rcyl
124 return (vx, vy, 0.0)
125
126
127 def mag_func(r):
128 x, y, z = r
129 rho_local = rho_disc if (x * x + y * y) < R_disc * R_disc else rho_ambient
130 return (Bx0 / rho_local, 0.0, 0.0)
131
132
133 model.set_field_value_lambda_f64_3("vxyz", vel_func)
134 model.set_field_value_lambda_f64_3("B/rho", mag_func)
135
136 vol_b = xs * ys * zs
137 vol_disc = np.pi * R_disc * R_disc * zs
138 totmass = rho_ambient * (vol_b - vol_disc) + rho_disc * vol_disc
139 pmass = model.total_mass_to_part_mass(totmass)
140 model.set_particle_mass(pmass)
141
142 print("Total mass :", totmass)
143 print("Current part mass :", pmass)
144 print("Total particle count :", model.get_total_part_count())
145
146 model.set_cfl_cour(C_cour)
147 model.set_cfl_force(C_force)
148
149 model.timestep()
150
151 model.evolve_until(t_target)
SPH setup: generating particles ...
SPH setup: Nstep = 16057 ( 1.6e+04 ) Ntotal = 16057 ( 1.6e+04 rank min = 7.9e+05 max = 1.6e+04) rate = 1.605700e+04 N.s^-1
SPH setup: the generation step took : 0.023222084 s
SPH setup: final particle count = 16057 beginning injection ...
SPH setup: injected 16057 / 16057 => 100.0% | ranks with patchs = 1 / 1 <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.004735663 s
SPH setup: the setup took : 0.032173955000000004 s
Total mass : 0.3517371564808811
Current part mass : 2.1905533815836154e-05
Total particle count : 16057
---------------- 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
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.030624643436074166 unconverged cnt = 15757
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.033441290467993 unconverged cnt = 1412
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.03344129046799447 unconverged cnt = 612
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.03344129046799447 unconverged cnt = 24
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0, dt = 9.361305039293543e-06 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 9.361305039293543e-06, dt = 0.0003179667331182061 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.00032732803815749967, dt = 0.0005068913408230215 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.02140365637482218
---------------- t = 0.0008342193789805212, dt = 0.0003088471520994329 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.0011430665310799542, dt = 0.0004811609334120151 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.01650282944842808
---------------- t = 0.0016242274644919694, dt = 0.00028852941834372825 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.0019127568828356976, dt = 0.0004513250590017535 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.011857336605177899
---------------- t = 0.002364081941837451, dt = 0.0002855860254006608 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.0026496679672381117, dt = 0.0004472462443549473 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.003096914211593059, dt = 0.0005430295859126931 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.010949868689275964
---------------- t = 0.003639943797505752, dt = 0.000299592007521199 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.0039395358050269514, dt = 0.00044019610690170175 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.004379731911928653, dt = 0.0005272079132851575 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.00490693982521381, dt = 0.0005800282773933735 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.005486968102607184, dt = 0.0006117457754147015 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.006098713878021886, dt = 0.0006306730880885195 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.006729386966110406, dt = 0.0006420545873200215 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.007371441553430427, dt = 0.0006674765828164316 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.008038918136246858, dt = 0.0006881561454382069 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.008727074281685065, dt = 0.0006876285612350837 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.009414702842920149, dt = 0.0006876819242558283 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.010102384767175978, dt = 0.0006884869190354837 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.010790871686211461, dt = 0.0006901904934907175 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.011481062179702179, dt = 0.0007246555609302711 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.012205717740632449, dt = 0.0007380571802197036 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.012943774920852153, dt = 0.0007432961584233591 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.013687071079275512, dt = 0.0007946193052820766 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.014481690384557589, dt = 0.0008043857606333408 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.01528607614519093, dt = 0.000810557241282707 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.016096633386473637, dt = 0.000869798410645941 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.01696643179711958, dt = 0.0008880191723396412 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.017854450969459222, dt = 0.0009223967573505114 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.018776847726809733, dt = 0.0009555110223614668 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.0197323587491712, dt = 0.0009716676481982024 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.0207040263973694, dt = 0.0010310594646010163 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.021735085861970416, dt = 0.0010887068756253 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.022823792737595716, dt = 0.0011425071968699525 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.023966299934465668, dt = 0.0011601819691180492 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.025126481903583717, dt = 0.0011760313594803226 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.02630251326306404, dt = 0.0012554323209557296 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.027557945584019767, dt = 0.0012969452529248983 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.028854890836944664, dt = 0.0013214152009816567 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.030176306037926322, dt = 0.0013496314261051182 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.03152593746403144, dt = 0.0013440212592291147 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.03286995872326055, dt = 0.0013384536616935679 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.03420841238495412, dt = 0.0013349538797200041 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.03554336626467412, dt = 0.001333383600951829 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.03687674986562595, dt = 0.001340660334262886 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.03821741019988884, dt = 0.00134301079756071 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.03956042099744955, dt = 0.0013469611038881644 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.040907382101337715, dt = 0.0013773134135356031 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.04228469551487332, dt = 0.0013902880029702594 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.04367498351784358, dt = 0.0013976001253743586 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.04507258364321794, dt = 0.001439730267489728 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.04651231391070766, dt = 0.0014573101904707842 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.04796962410117844, dt = 0.0014769751924910638 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.049446599293669506, dt = 0.0014977457475042247 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.05094434504117373, dt = 0.0016130734755654937 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.05255741851673923, dt = 0.001560132240544997 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.05411755075728423, dt = 0.0015791089188168035 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.05569665967610103, dt = 0.001604555836432781 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.05730121551253381, dt = 0.001624904882025502 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.058926120394559314, dt = 0.0016661173622373644 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.06059223775679668, dt = 0.001664261983037761 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.062256499739834444, dt = 0.0017138723596756653 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.0639703720995101, dt = 0.001688242091809867 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.06565861419131998, dt = 0.0017368068814297465 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.06739542107274972, dt = 0.0017640110309582917 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.06915943210370801, dt = 0.0017810701446123861 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.0709405022483204, dt = 0.0017802871210402858 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.07272078936936067, dt = 0.0018171687504317034 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.07453795811979237, dt = 0.0018639400204733703 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.07640189814026574, dt = 0.0018354687195079366 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.07823736685977367, dt = 0.0018934198296525436 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.08013078668942622, dt = 0.001888444440242755 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.08201923112966897, dt = 0.001962635009866464 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.08398186613953543, dt = 0.002024510240646751 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.08600637638018217, dt = 0.0019941766544564043 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.08800055303463858, dt = 0.0019826948515308633 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.08998324788616945, dt = 0.0020357798477004117 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.09201902773386986, dt = 0.00199297026940463 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.09401199800327449, dt = 0.0020549440261365464 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.09606694202941103, dt = 0.002134463284204309 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.09820140531361535, dt = 0.002201741591010624 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.10040314690462597, dt = 0.0022040208203038 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.10260716772492977, dt = 0.002239544767372461 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.10484671249230224, dt = 0.0022528099849507535 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.10709952247725299, dt = 0.002293877361630998 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.10939339983888399, dt = 0.0023263480987760704 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.11171974793766007, dt = 0.0023577993869540403 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.11407754732461411, dt = 0.002332722196151508 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.11641026952076562, dt = 0.0023578627919003925 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.010095640916472443
---------------- t = 0.11876813231266602, dt = 0.0012047836661210328 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.11997291597878705, dt = 0.0016243952988796132 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.12159731127766667, dt = 0.0019273486369679076 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.12352465991463457, dt = 0.001960795727008812 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.010385101868818584
---------------- t = 0.12548545564164337, dt = 0.0010247839354388978 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.12651023957708227, dt = 0.0014365352514479684 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.12794677482853023, dt = 0.00170953395210459 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.12965630878063483, dt = 0.0019038935488683586 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.011160061032476161
---------------- t = 0.1315602023295032, dt = 0.0010273980210009173 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.1325876003505041, dt = 0.0014600945047063818 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.13404769485521048, dt = 0.001763836488213375 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.13581153134342386, dt = 0.001967828367139376 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.010707174971430781
---------------- t = 0.13777935971056324, dt = 0.0010659001646381998 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.13884525987520144, dt = 0.0015199181104031286 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.14036517798560458, dt = 0.001834190117654256 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.14219936810325884, dt = 0.002042259790331768 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
Warning: the corrector tolerance are broken the step will be re rerunned [BasicGasSPH][rank=0]
eps_v = 0.010512362264616656
---------------- t = 0.1442416278935906, dt = 0.0011120022534457327 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.14535363014703634, dt = 0.001593370625881224 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.14694700077291756, dt = 0.0019219380479688205 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
---------------- t = 0.14886893882088637, dt = 0.0011310611791136294 ----------------
Warning: High interface/patch volume ratio. [InterfaceGen][rank=0]
This can lead to high mpi overhead, try to increase the patch split crit
EvolveUntilResults(reach_target_time=true, reach_niter_max=false, reach_max_walltime=false, iter_count=112)
Diagnostics
156 kinetic_energy = shamrock.model_sph.analysisEnergyKinetic(model=model).get_kinetic_energy()
157 total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum()
158
159 data = ctx.collect_data()
160 hpart = data["hpart"]
161 hfact = model.get_hfact()
162 rho = pmass * (hfact / hpart) ** 3
163 B_on_rho = data["B/rho"]
164
165 magnetic_energy = 0.5 * pmass * np.sum(rho * np.sum(B_on_rho**2, axis=1)) / mu_0
166
167 rho_max = np.max(rho)
168 rho_min = np.min(rho)
169
170 momentum_norm = np.sqrt(total_momentum[0] ** 2 + total_momentum[1] ** 2 + total_momentum[2] ** 2)
171
172 # nx=256
173 do_plot = False
174 dump_folder = "_to_trash"
175 if do_plot and shamrock.sys.world_rank() == 0:
176 os.makedirs(dump_folder, exist_ok=True)
177
178 if do_plot:
179 render_res = 300
180 center = (0.0, 0.0, 0.0)
181 delta_x = (xs, 0.0, 0.0)
182 delta_y = (0.0, ys, 0.0)
183
184 def render_slice_normalized(field, field_type=None, min_normalization=1e-9):
185 kwargs = dict(center=center, delta_x=delta_x, delta_y=delta_y, nx=render_res, ny=render_res)
186 if field_type is None:
187 raw = np.asarray(model.render_cartesian_slice(field, **kwargs))
188 else:
189 raw = np.asarray(model.render_cartesian_slice(field, field_type, **kwargs))
190
191 unity = np.asarray(model.render_cartesian_slice("unity", "f64", **kwargs))
192 unity_b = unity if raw.ndim == 2 else unity[..., None]
193 out = raw / unity_b
194 mask = np.broadcast_to(unity_b, raw.shape) < min_normalization
195 return np.where(mask, np.nan, out)
196
197 rho_field = model.compute_field("rho", "f64")
198 uint_field = model.compute_field("uint", "f64")
199
200 def compute_pressure(size, rho, uint):
201 return (gamma - 1.0) * rho * uint
202
203 P_field = shamrock.map_fields_f64(compute_pressure, rho=rho_field, uint=uint_field)
204
205 def compute_soundspeed(size, rho, P):
206 return np.sqrt(gamma * P / rho)
207
208 cs_field = shamrock.map_fields_f64(compute_soundspeed, rho=rho_field, P=P_field)
209
210 rho_slice = render_slice_normalized(rho_field)
211 P_slice = render_slice_normalized(P_field)
212 cs_slice = render_slice_normalized(cs_field)
213 vxyz_slice = render_slice_normalized("vxyz", "f64_3")
214 B_on_rho_slice = render_slice_normalized("B/rho", "f64_3")
215
216 vmag_slice = np.linalg.norm(vxyz_slice, axis=-1)
217 mach_slice = vmag_slice / cs_slice
218
219 # physical B = (B/rho) * rho
220 B_slice = B_on_rho_slice * rho_slice[..., None]
221 magnetic_pressure_slice = 0.5 * np.sum(B_slice**2, axis=-1) / mu_0
222
223 if shamrock.sys.world_rank() == 0:
224 # Contour limits identical to Toth (2000), as quoted in the
225 # Phantom paper's Figure 34 caption.
226 panels = [
227 ("Density", rho_slice, 0.483, 12.95),
228 ("Pressure", P_slice, 0.0202, 2.008),
229 ("Mach number $|v|/c_s$", mach_slice, 0.0, 1.09),
230 ("Magnetic pressure $\\frac{1}{2}B^2$", magnetic_pressure_slice, 0.0, 2.642),
231 ]
232
233 ny_render, nx_render = rho_slice.shape
234 xg = np.linspace(box_min[0], box_max[0], nx_render)
235 yg = np.linspace(box_min[1], box_max[1], ny_render)
236 X, Y = np.meshgrid(xg, yg)
237
238 fig, axs = plt.subplots(2, 2, figsize=(10, 9))
239
240 for ax, (title, field_slice, vmin, vmax) in zip(axs.flat, panels):
241 levels = np.linspace(vmin, vmax, 30)
242 contour_set = ax.contour(X, Y, field_slice, levels=levels)
243 ax.set_title(title)
244 ax.set_xlabel("x")
245 ax.set_ylabel("y")
246 ax.set_aspect("equal")
247 fig.colorbar(contour_set, ax=ax)
248
249 fig.suptitle(f"MHD rotor problem, $t={t_target}$, nx={nx}")
250 fig.tight_layout()
251 fig.savefig(os.path.join(dump_folder, "mhd_rotor_contours.png"), dpi=150)
252 plt.close(fig)
253
254 print(f"kinetic_energy = {kinetic_energy}")
255 print(f"magnetic_energy = {magnetic_energy}")
256 print(f"momentum_norm = {momentum_norm}")
257 print(f"rho_max = {rho_max}")
258 print(f"rho_min = {rho_min}")
259
260 test_pass = True
261 err_log = ""
262
263 expect_kinetic_energy = 0.05985165376928185
264 expect_magnetic_energy = 0.33746158787028185
265 expect_rho_max = 6.3060516178274515
266 expect_rho_min = 0.24993176850956453
267
268 tol = 1e-3 # relative tolerance
269 max_momentum_norm = 1e-2
270
271
272 def float_equal(val1, val2, prec):
273 return abs(val1 - val2) < prec
274
275
276 if not float_equal(kinetic_energy, expect_kinetic_energy, tol * expect_kinetic_energy):
277 err_log += "error on the kinetic energy is outside of tolerances:\n"
278 err_log += f" expected error = {expect_kinetic_energy} +- {tol * expect_kinetic_energy}\n"
279 err_log += f" obtained error = {kinetic_energy}\n"
280 test_pass = False
281
282 if not float_equal(magnetic_energy, expect_magnetic_energy, tol * expect_magnetic_energy):
283 err_log += "error on the magnetic energy is outside of tolerances:\n"
284 err_log += f" expected error = {expect_magnetic_energy} +- {tol * expect_magnetic_energy}\n"
285 err_log += f" obtained error = {magnetic_energy}\n"
286 test_pass = False
287
288 if momentum_norm > max_momentum_norm:
289 err_log += "total momentum norm is outside of tolerances:\n"
290 err_log += f" expected |p| < {max_momentum_norm}\n"
291 err_log += f" obtained |p| = {momentum_norm}\n"
292 test_pass = False
293
294 if not float_equal(rho_max, expect_rho_max, tol * expect_rho_max):
295 err_log += "error on rho_max is outside of tolerances:\n"
296 err_log += f" expected error = {expect_rho_max} +- {tol * expect_rho_max}\n"
297 err_log += f" obtained error = {rho_max}\n"
298 test_pass = False
299
300 if not float_equal(rho_min, expect_rho_min, tol * expect_rho_min):
301 err_log += "error on rho_min is outside of tolerances:\n"
302 err_log += f" expected error = {expect_rho_min} +- {tol * expect_rho_min}\n"
303 err_log += f" obtained error = {rho_min}\n"
304 test_pass = False
305
306 if test_pass == False:
307 exit("Test did not pass L2 margins : \n" + err_log)
kinetic_energy = 0.05985165376928191
magnetic_energy = 0.33746158787028185
momentum_norm = 0.00015219398459745176
rho_max = 6.3060516178274515
rho_min = 0.24993176850956453
Total running time of the script: (0 minutes 15.989 seconds)
Estimated memory usage: 162 MB