Note
Go to the end to download the full example code.
Disc Monte-Carlo setup density profiles#
Setup a disc with the Monte-Carlo disc generator, run a single step and compare the resulting density profiles to the setup surface density profile.
----- SPH Solver configuration -----
[
{
"artif_viscosity": {
"alpha_AV": 0.0765404518492,
"alpha_u": 1.0,
"beta_AV": 2.0,
"type": "constant_disc"
},
"boundary_config": {
"bc_type": "free"
},
"cfl_config": {
"cfl_cour": 0.3,
"cfl_force": 0.25,
"cfl_multiplier_stiffness": 2.0,
"eta_sink": 0.05
},
"combined_dtdiv_divcurlv_compute": false,
"debug_dump_filename": "",
"do_debug_dump": false,
"dust_config": {
"ballabio_ts_limiter": false,
"drag_mode": {
"type": "none"
},
"evol_mode": {
"type": "none"
},
"mode": {
"type": "none"
}
},
"enable_particle_reordering": false,
"eos_config": {
"Tvec": "f64_3",
"cs0": 0.04999981727640601,
"eos_type": "locally_isothermal_lp07",
"q": 0.5,
"r0": 1.0
},
"epsilon_h": 1e-06,
"ext_force_config": {
"force_list": []
},
"gpart_mass": 5e-07,
"h_iter_per_subcycles": 50,
"h_max_subcycles_count": 100,
"htol_up_coarse_cycle": 1.1,
"htol_up_fine_cycle": 1.1,
"kernel_id": "M4<f64>",
"mhd_config": {
"mhd_type": "none"
},
"neigh_cache_strategy": "two_stage",
"particle_killing": [
{
"center": [
0.0,
0.0,
0.0
],
"radius": 12.0,
"type": "sphere"
}
],
"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": 149597870700.0,
"unit_lumint": 1.0,
"unit_mass": 1.98847e+30,
"unit_qte": 1.0,
"unit_temperature": 1.0,
"unit_time": 5022548.032116797
}
}
]
------------------------------------
Warning: make_generator_disc_mc: with the current EOS, cs_profile is ignored [SPHSetup][rank=0]
SPH setup: generating particles ...
SPH setup: Nstep = 100000 ( 1.0e+05 ) Ntotal = 100000 ( 1.0e+05 rank min = 4.9e+05 max = 1.0e+05) rate = 1.000000e+05 N.s^-1
SPH setup: the generation step took : 0.211467649 s
SPH setup: final particle count = 100000 beginning injection ...
SPH setup: injected 100000 / 100000 => 100.0% | ranks with patchs = 1 / 1 <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.009937711 s
SPH setup: the setup took : 0.23699259700000003 s
---------------- t = 0, dt = 0 ----------------
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.3780158886363346 unconverged cnt = 99999
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.4158174774999681 unconverged cnt = 99998
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.45739922524996496 unconverged cnt = 99998
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.5031391477749615 unconverged cnt = 99996
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.5534530625524577 unconverged cnt = 99995
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.6087983688077035 unconverged cnt = 99992
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.6696782056884739 unconverged cnt = 99992
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7366460262573213 unconverged cnt = 99989
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.8103106288830535 unconverged cnt = 99986
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.8913416917713589 unconverged cnt = 99983
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757872 unconverged cnt = 99975
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 99964
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 99948
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 99921
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 99874
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 99787
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 99600
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 99081
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 97730
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 94665
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 89429
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 82181
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 72804
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 60804
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 45114
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 27369
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 11761
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 3193
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 714
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757873 unconverged cnt = 147
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9422982425757874 unconverged cnt = 19
----- SPH Solver configuration -----
[
{
"artif_viscosity": {
"alpha_AV": 0.0765404518492,
"alpha_u": 1.0,
"beta_AV": 2.0,
"type": "constant_disc"
},
"boundary_config": {
"bc_type": "free"
},
"cfl_config": {
"cfl_cour": 0.3,
"cfl_force": 0.25,
"cfl_multiplier_stiffness": 2.0,
"eta_sink": 0.05
},
"combined_dtdiv_divcurlv_compute": false,
"debug_dump_filename": "",
"do_debug_dump": false,
"dust_config": {
"ballabio_ts_limiter": false,
"drag_mode": {
"type": "none"
},
"evol_mode": {
"type": "none"
},
"mode": {
"type": "none"
}
},
"enable_particle_reordering": false,
"eos_config": {
"Tvec": "f64_3",
"cs0": 0.04999981727640601,
"eos_type": "locally_isothermal_lp07",
"q": 0.5,
"r0": 1.0
},
"epsilon_h": 1e-06,
"ext_force_config": {
"force_list": []
},
"gpart_mass": 5e-07,
"h_iter_per_subcycles": 50,
"h_max_subcycles_count": 100,
"htol_up_coarse_cycle": 1.1,
"htol_up_fine_cycle": 1.1,
"kernel_id": "M4<f64>",
"mhd_config": {
"mhd_type": "none"
},
"neigh_cache_strategy": "two_stage",
"particle_killing": [
{
"center": [
0.0,
0.0,
0.0
],
"radius": 12.0,
"type": "sphere"
}
],
"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": 149597870700.0,
"unit_lumint": 1.0,
"unit_mass": 1.98847e+30,
"unit_qte": 1.0,
"unit_temperature": 1.0,
"unit_time": 5022548.032116797
}
}
]
------------------------------------
Warning: make_generator_disc_mc: with the current EOS, cs_profile is ignored [SPHSetup][rank=0]
SPH setup: generating particles ...
SPH setup: Nstep = 100000 ( 1.0e+05 ) Ntotal = 100000 ( 1.0e+05 rank min = 4.1e+05 max = 1.0e+05) rate = 1.000000e+05 N.s^-1
SPH setup: the generation step took : 0.24411791600000002 s
SPH setup: final particle count = 100000 beginning injection ...
SPH setup: injected 100000 / 100000 => 100.0% | ranks with patchs = 1 / 1 <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.0041283150000000005 s
SPH setup: the setup took : 0.256689263 s
---------------- t = 0, dt = 0 ----------------
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.6783553284643802 unconverged cnt = 99997
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.6783553284643803 unconverged cnt = 99994
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.6783553284643803 unconverged cnt = 99993
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.6783553284643803 unconverged cnt = 99991
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.720214305328142 unconverged cnt = 99987
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.7922357358609563 unconverged cnt = 99982
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.8330097352644652 unconverged cnt = 99978
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.8330097352644652 unconverged cnt = 99973
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9040452070332429 unconverged cnt = 99961
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9944497277365673 unconverged cnt = 99947
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.093894700510224 unconverged cnt = 99917
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469016082 unconverged cnt = 99866
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 99771
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 99575
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 99063
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 97530
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 93597
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 86509
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 77376
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 66725
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 54334
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 40929
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 26638
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 13009
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 4322
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 922
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 195
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 21
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0959018469017643 unconverged cnt = 2
----- SPH Solver configuration -----
[
{
"artif_viscosity": {
"alpha_AV": 0.0765404518492,
"alpha_u": 1.0,
"beta_AV": 2.0,
"type": "constant_disc"
},
"boundary_config": {
"bc_type": "free"
},
"cfl_config": {
"cfl_cour": 0.3,
"cfl_force": 0.25,
"cfl_multiplier_stiffness": 2.0,
"eta_sink": 0.05
},
"combined_dtdiv_divcurlv_compute": false,
"debug_dump_filename": "",
"do_debug_dump": false,
"dust_config": {
"ballabio_ts_limiter": false,
"drag_mode": {
"type": "none"
},
"evol_mode": {
"type": "none"
},
"mode": {
"type": "none"
}
},
"enable_particle_reordering": false,
"eos_config": {
"Tvec": "f64_3",
"cs0": 0.04999981727640601,
"eos_type": "locally_isothermal_lp07",
"q": 0.5,
"r0": 1.0
},
"epsilon_h": 1e-06,
"ext_force_config": {
"force_list": []
},
"gpart_mass": 5e-07,
"h_iter_per_subcycles": 50,
"h_max_subcycles_count": 100,
"htol_up_coarse_cycle": 1.1,
"htol_up_fine_cycle": 1.1,
"kernel_id": "M4<f64>",
"mhd_config": {
"mhd_type": "none"
},
"neigh_cache_strategy": "two_stage",
"particle_killing": [
{
"center": [
0.0,
0.0,
0.0
],
"radius": 12.0,
"type": "sphere"
}
],
"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": 149597870700.0,
"unit_lumint": 1.0,
"unit_mass": 1.98847e+30,
"unit_qte": 1.0,
"unit_temperature": 1.0,
"unit_time": 5022548.032116797
}
}
]
------------------------------------
Warning: make_generator_disc_mc: with the current EOS, cs_profile is ignored [SPHSetup][rank=0]
SPH setup: generating particles ...
SPH setup: Nstep = 100000 ( 1.0e+05 ) Ntotal = 100000 ( 1.0e+05 rank min = 3.6e+05 max = 1.0e+05) rate = 1.000000e+05 N.s^-1
SPH setup: the generation step took : 0.280932824 s
SPH setup: final particle count = 100000 beginning injection ...
SPH setup: injected 100000 / 100000 => 100.0% | ranks with patchs = 1 / 1 <- global loop -> (msg count : 0)
SPH setup: the injection step took : 0.004177049 s
SPH setup: the setup took : 0.29394506800000003 s
---------------- t = 0, dt = 0 ----------------
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.8182003274206934 unconverged cnt = 99994
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9000203601627628 unconverged cnt = 99992
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 0.9900223961790391 unconverged cnt = 99985
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0022740097793776 unconverged cnt = 99980
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0022740097793779 unconverged cnt = 99976
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0022740097793779 unconverged cnt = 99972
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0022740097793779 unconverged cnt = 99958
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.0022740097793779 unconverged cnt = 99938
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.030546130122973 unconverged cnt = 99898
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1336007431352704 unconverged cnt = 99823
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.170768558067543 unconverged cnt = 99688
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675435 unconverged cnt = 99350
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675438 unconverged cnt = 98317
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675435 unconverged cnt = 95079
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675438 unconverged cnt = 87683
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675435 unconverged cnt = 76395
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675438 unconverged cnt = 63975
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675435 unconverged cnt = 51239
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675438 unconverged cnt = 39078
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675435 unconverged cnt = 27257
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675438 unconverged cnt = 16110
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675435 unconverged cnt = 7295
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675438 unconverged cnt = 2161
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675435 unconverged cnt = 482
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675438 unconverged cnt = 113
Warning: smoothing length is not converged, rerunning the iterator ... [Smoothinglength][rank=0]
largest h = 1.1707685580675435 unconverged cnt = 15
p=1.0: max relative delta sigma (2 R_in < r < 0.8 R_out): 1.2814e-01, fitted slope: -1.0266 (expected -1.0)
p=1.5: max relative delta sigma (2 R_in < r < 0.8 R_out): 1.1869e-01, fitted slope: -1.5255 (expected -1.5)
p=2.0: max relative delta sigma (2 R_in < r < 0.8 R_out): 1.2396e-01, fitted slope: -2.0090 (expected -2.0)
9 import matplotlib.pyplot as plt
10 import numpy as np
11
12 import shamrock
13
14 # If we use the shamrock executable to run this script instead of the python interpreter,
15 # we should not initialize the system as the shamrock executable needs to handle specific MPI logic
16 if not shamrock.sys.is_initialized():
17 shamrock.change_loglevel(1)
18 shamrock.sys.init("0:0")
19
20 si = shamrock.UnitSystem()
21 sicte = shamrock.Constants(si)
22 codeu = shamrock.UnitSystem(
23 unit_time=sicte.year() / (2 * np.pi),
24 unit_length=sicte.au(),
25 unit_mass=sicte.sol_mass(),
26 )
27 ucte = shamrock.Constants(codeu)
28 G = ucte.G()
29
30 R_in = 1.0
31 R_out = 10.0
32 R_ref = 1.0
33 disc_mass = 0.05
34 Npart = 100000
35 p_index_list = [1.0, 1.5, 2.0]
36 p_index = p_index_list[0] # overwritten in the loop below
37 q_index = 0.50
38 H_R = 0.05
39 alpha_SS = 0.005
40 m1 = 1.0
41
42 pmass = disc_mass / Npart
43
44
45 # Disc profiles
46 def sigma_profile(r):
47 sigma_0 = 1.0 # We do not care as it will be renormalized
48 return sigma_0 * (r / R_ref) ** (-p_index)
49
50
51 def kep_profile(r):
52 return (G * m1 / r) ** 0.5
53
54
55 def omega_k(r):
56 return kep_profile(r) / r
57
58
59 def cs_profile(r):
60 cs_in = (H_R * R_ref) * omega_k(R_ref)
61 return ((r / R_ref) ** (-q_index)) * cs_in
62
63
64 def get_sigma_norm():
65 x_list = np.linspace(R_in, R_out, 2048)
66 term = [sigma_profile(x) * 2 * np.pi * x for x in x_list]
67 return disc_mass / (np.sum(term) * (x_list[1] - x_list[0]))
68
69
70 def rot_profile(r):
71 return ((kep_profile(r) ** 2) - (2 * p_index + q_index) * cs_profile(r) ** 2) ** 0.5
72
73
74 def H_profile(r):
75 H = cs_profile(r) / omega_k(r)
76 fact = 1.0
77 return fact * H
78
79
80 def setup_shamrock_disc():
81 ctx = shamrock.Context()
82 ctx.pdata_layout_new()
83 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
84
85 alpha_AV = 0.0765404518492
86 alpha_u = 1.0
87 beta_AV = 2.0
88 C_cour = 0.3
89 C_force = 0.25
90 cs0 = cs_profile(R_ref)
91 bsize = R_out * 1.2
92
93 # Generate the default config
94 cfg = model.gen_default_config()
95 cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV)
96 cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q_index, r0=R_ref)
97
98 cfg.add_kill_sphere(center=(0, 0, 0), radius=bsize) # kill particles outside the simulation box
99
100 cfg.set_units(codeu)
101 cfg.set_particle_mass(pmass)
102 # Set the CFL
103 cfg.set_cfl_cour(C_cour)
104 cfg.set_cfl_force(C_force)
105
106 cfg.set_smoothing_length_density_based()
107
108 # Set the solver config to be the one stored in cfg
109 model.set_solver_config(cfg)
110
111 # Print the solver config
112 model.get_current_config().print_status()
113
114 # Init the scheduler & fields
115 model.init_scheduler(int(1e8), 1)
116
117 # Set the simulation box size
118 ext = R_out * 1.2
119 model.resize_simulation_box((-ext, -ext, -ext), (ext, ext, ext))
120
121 # Create the setup
122 setup = model.get_setup()
123 gen_disc = setup.make_generator_disc_mc(
124 part_mass=pmass,
125 disc_mass=disc_mass,
126 r_in=R_in,
127 r_out=R_out,
128 sigma_profile=sigma_profile,
129 H_profile=H_profile,
130 rot_profile=rot_profile,
131 cs_profile=cs_profile,
132 random_seed=666,
133 )
134
135 # Apply the setup
136 setup.apply_setup(gen_disc, insert_step=1000000)
137
138 # now that the barycenter & momentum are 0, we can add the sink
139 model.add_sink(m1, (0, 0, 0), (0, 0, 0), 1.0)
140
141 return ctx, model
142
143
144 def positions_to_rays(positions):
145 return [shamrock.math.Ray_f64_3(tuple(position), (0.0, 0.0, 1.0)) for position in positions]
146
147
148 def compute_avg_sigma_profile(model, ntheta, r):
149 theta = np.linspace(0, 2 * np.pi, ntheta)
150
151 r_grid, theta_grid = np.meshgrid(r, theta)
152 x_grid = r_grid * np.cos(theta_grid)
153 y_grid = r_grid * np.sin(theta_grid)
154 z_grid = np.zeros_like(r_grid)
155
156 positions = np.column_stack([x_grid.ravel(), y_grid.ravel(), z_grid.ravel()])
157
158 rays = positions_to_rays(positions)
159 arr_sigma = model.render_column_integ("rho", "f64", rays)
160
161 arr_sigma = np.array(arr_sigma).reshape(ntheta, len(r))
162
163 # average over the theta direction
164 arr_sigma = np.mean(arr_sigma, axis=0)
165 return arr_sigma
166
167
168 def get_profiles(model):
169 x_list = np.linspace(0, R_out * 1.2, 2049)[1:]
170 positions = [(x, 0.0, 0.0) for x in x_list.tolist()]
171
172 arr_rho = np.array(model.render_slice("rho", "f64", positions))
173 arr_sigma_avg = compute_avg_sigma_profile(model, 32, x_list)
174
175 arr_sigma = np.array([sigma_profile(x) * sigma_norm for x in x_list])
176 arr_sigma[(x_list < R_in) | (x_list > R_out)] = 0.0
177
178 return x_list, arr_rho, arr_sigma_avg, arr_sigma
179
180
181 dpi = 200
182
183 fig_sigma, ax_sigma = plt.subplots(dpi=dpi)
184 fig_rho, ax_rho = plt.subplots(dpi=dpi)
185 fig_err, ax_err = plt.subplots(dpi=dpi)
186
187 max_delta_sigma = {}
188 fitted_slope = {}
189
190 for i, p_index in enumerate(p_index_list):
191 # the disc profiles read p_index & sigma_norm from the module globals
192 sigma_norm = get_sigma_norm()
193
194 ctx, model = setup_shamrock_disc()
195
196 # Run a single step so that the smoothing length & density are computed
197 model.timestep()
198
199 x_list, arr_rho, arr_sigma_avg, arr_sigma = get_profiles(model)
200
201 del model
202 del ctx
203
204 H = np.array([H_profile(x) for x in x_list])
205 # midplane density expected from the setup surface density profile
206 arr_rho_setup = arr_sigma / (np.sqrt(2 * np.pi) * H)
207
208 color = f"C{i}"
209 lbl = f" (p={p_index})"
210
211 # Surface density profile
212 ax_sigma.plot(x_list, arr_sigma_avg, color=color, label="sigma avg (shamrock)" + lbl)
213 ax_sigma.plot(x_list, arr_sigma, color=color, ls="--", label="sigma (setup)" + lbl)
214
215 # Midplane density profile
216 ax_rho.plot(x_list, arr_rho, color=color, label="rho (shamrock)" + lbl)
217 ax_rho.plot(
218 x_list,
219 arr_sigma_avg / (np.sqrt(2 * np.pi) * H),
220 color=color,
221 ls=":",
222 label="sigma avg/(sqrt(2*pi)*H)" + lbl,
223 )
224 ax_rho.plot(
225 x_list, arr_rho_setup, color=color, ls="--", label="sigma/(sqrt(2*pi)*H) (setup)" + lbl
226 )
227
228 # Relative error on the surface density
229 mask = arr_sigma > 0
230 delta_sigma = np.zeros_like(arr_sigma)
231 delta_sigma[mask] = (arr_sigma_avg[mask] - arr_sigma[mask]) / arr_sigma[mask]
232 ax_err.plot(x_list, np.abs(delta_sigma), color=color, label="|delta sigma| / sigma" + lbl)
233
234 # ignore the disc edges where the SPH smoothing makes the profile deviate
235 mask_inner = (x_list > 2 * R_in) & (x_list < 0.8 * R_out)
236 max_delta_sigma[p_index] = np.max(np.abs(delta_sigma[mask_inner]))
237
238 # power law index of the measured surface density (should be -p_index)
239 fitted_slope[p_index] = np.polyfit(
240 np.log(x_list[mask_inner]), np.log(arr_sigma_avg[mask_inner]), 1
241 )[0]
242
243 ax_sigma.set_xlabel("r")
244 ax_sigma.set_ylabel("sigma")
245 ax_sigma.set_title("Surface density after one step")
246 ax_sigma.legend(fontsize="small")
247 ax_sigma.set_ylim(1e-6, 1e-2)
248 ax_sigma.set_yscale("log")
249
250 ax_rho.set_xlabel("r")
251 ax_rho.set_ylabel("rho")
252 ax_rho.set_title("Midplane density after one step")
253 ax_rho.legend(fontsize="small")
254 ax_rho.set_ylim(1e-5, 1e-1)
255 ax_rho.set_yscale("log")
256
257 ax_err.set_xlabel("r")
258 ax_err.set_ylabel("relative error")
259 ax_err.set_yscale("log")
260 ax_err.set_title("Surface density relative error after one step")
261 ax_err.legend(fontsize="small")
262
263 plt.show()
264
265 # With a biased radial sampling (e.g. the rejection sampling bound being too low, which
266 # makes r uniform) the max relative error is ~0.36 for p=1.5 and ~1.05 for p=2,
267 # and the fitted slope is -1 regardless of p
268 max_rel_err_tol = 0.2
269 slope_tol = 0.1
270
271 errors = []
272 for p_index in p_index_list:
273 err = max_delta_sigma[p_index]
274 slope = fitted_slope[p_index]
275 print(
276 f"p={p_index}: max relative delta sigma (2 R_in < r < 0.8 R_out): {err:.4e}, "
277 f"fitted slope: {slope:.4f} (expected {-p_index})"
278 )
279
280 if err > max_rel_err_tol:
281 errors.append(f"p={p_index}: max relative delta sigma {err:.4e} > {max_rel_err_tol}")
282 if abs(slope + p_index) > slope_tol:
283 errors.append(
284 f"p={p_index}: fitted slope {slope:.4f} differs from {-p_index} by more than {slope_tol}"
285 )
286
287 if errors:
288 raise ValueError("Disc MC profile check failed:\n" + "\n".join(errors))
Total running time of the script: (0 minutes 14.744 seconds)
Estimated memory usage: 209 MB