diff --git a/examples/production/compare_torque_free.py b/examples/production/compare_torque_free.py new file mode 100644 index 0000000000..d46bbdc6fa --- /dev/null +++ b/examples/production/compare_torque_free.py @@ -0,0 +1,857 @@ +""" +Production run: Circular disc & central sink particle +===================================================== + +This example demonstrates how to run a smoothed particle hydrodynamics (SPH) +simulation of a circular disc orbiting around a central sink. + +The simulation models: + +- A central star with a given mass and accretion radius +- A gaseous disc with specified mass, inner/outer radii, and vertical structure +- Artificial viscosity for angular momentum transport +- Locally isothermal equation of state + +Also this simulation feature rolling dumps (see `purge_old_dumps` function) to save disk space. + +This example is the accumulation of 3 files in a single one to showcase the complete workflow. + +- The actual run script (runscript.py) +- Plot generation (make_plots.py) +- Animation from the plots (plot_to_gif.py) + +On a cluster or laptop, one can run the code as follows: + +.. code-block:: bash + + mpirun ./shamrock --sycl-cfg 0:0 --loglevel 1 --rscript runscript.py + + +then after the run is done (or while it is running), one can run the following to generate the plots: + +.. code-block:: bash + + python make_plots.py + + +""" + +# %% +# Runscript (runscript.py) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# The runscript is the actual simulation with on the fly analysis & rolling dumps + +import glob +import json +import os # for makedirs + +import numpy as np + +try: + from numba import njit + + _HAS_NUMBA = True +except ImportError: + _HAS_NUMBA = False + + +import shamrock + +# If we use the shamrock executable to run this script instead of the python interpreter, +# we should not initialize the system as the shamrock executable needs to handle specific MPI logic +if not shamrock.sys.is_initialized(): + shamrock.change_loglevel(1) + shamrock.sys.init("0:0") + +# %% +# Setup units + +si = shamrock.UnitSystem() +sicte = shamrock.Constants(si) +codeu = shamrock.UnitSystem( + unit_time=sicte.year(), + unit_length=sicte.au(), + unit_mass=sicte.sol_mass(), +) +ucte = shamrock.Constants(codeu) +G = ucte.G() + +# %% +# List parameters + +# Resolution +Npart = int(float(os.environ.get("NPART"))) + +# Domain decomposition parameters +scheduler_split_val = int(1.0e7) # split patches with more than 1e7 particles +scheduler_merge_val = scheduler_split_val // 16 + +# Dump and plot frequency and duration of the simulation +dump_freq_stop = 2 +plot_freq_stop = 1 + +dt_stop = 1 +nstop = 100 + +# The list of times at which the simulation will pause for analysis / dumping +t_stop = [i * dt_stop for i in range(nstop + 1)] + +run_params = [] + +for center_is_torque_free in [True, False]: + for center_racc in [0.1, 0.25, 0.5, 0.8]: + if center_is_torque_free: + for center_torque_boost_radius_fact in [1.5, 2.0, 3.0]: + run_params.append( + (center_racc, center_is_torque_free, center_torque_boost_radius_fact) + ) + else: + run_params.append((center_racc, center_is_torque_free, 0.0)) + +print(run_params, len(run_params)) + +run_id = int(os.environ.get("RUN_ID")) + +# Sink parameters +center_mass = 1.0 +center_racc = run_params[run_id][0] +center_is_torque_free = run_params[run_id][1] +center_torque_boost_radius_fact = run_params[run_id][2] + +# Disc parameter +disc_mass = 0.01 # sol mass +rout = 20.0 # au +rin = 1.0 # au +H_r_0 = 0.05 +q = 0.5 +p = 3.0 / 2.0 +r0 = 1.0 + +# Viscosity parameter +alpha_AV = 1.0e-3 / 0.08 +alpha_u = 1.0 +beta_AV = 2.0 + +# Integrator parameters +C_cour = 0.3 +C_force = 0.25 + + +sim_folder = f"_to_trash/circular_disc_sink_{center_is_torque_free}_{center_racc}_{center_torque_boost_radius_fact}_{Npart}/" + +print(sim_folder) + +dump_folder = sim_folder + "dump/" +analysis_folder = sim_folder + "analysis/" +plot_folder = analysis_folder + "plots/" + +dump_prefix = dump_folder + "dump_" + + +# Disc profiles +def sigma_profile(r): + sigma_0 = 1.0 # We do not care as it will be renormalized + return sigma_0 * (r / r0) ** (-p) + + +def kep_profile(r): + return (G * center_mass / r) ** 0.5 + + +def omega_k(r): + return kep_profile(r) / r + + +def cs_profile(r): + cs_in = (H_r_0 * r0) * omega_k(r0) + return ((r / r0) ** (-q)) * cs_in + + +# %% +# Create the dump directory if it does not exist +if shamrock.sys.world_rank() == 0: + os.makedirs(sim_folder, exist_ok=True) + os.makedirs(dump_folder, exist_ok=True) + os.makedirs(analysis_folder, exist_ok=True) + os.makedirs(plot_folder, exist_ok=True) + +# %% +# Utility functions and quantities deduced from the base one + +# Deduced quantities +pmass = disc_mass / Npart + +bsize = rout * 2 +bmin = (-bsize, -bsize, -bsize) +bmax = (bsize, bsize, bsize) + +cs0 = cs_profile(r0) + + +def rot_profile(r): + return ((kep_profile(r) ** 2) - (2 * p + q) * cs_profile(r) ** 2) ** 0.5 + + +def H_profile(r): + H = cs_profile(r) / omega_k(r) + # fact = (2.**0.5) * 3. # factor taken from phantom, to fasten thermalizing + fact = 1.0 + return fact * H + + +# %% +# Start the context +# The context holds the data of the code +# We then init the layout of the field (e.g. the list of fields used by the solver) + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +# %% +# Attach a SPH model to the context + +model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4") + + +# %% +# Dump handling + + +def get_vtk_dump_name(idump): + return dump_prefix + f"{idump:07}" + ".vtk" + + +def get_ph_dump_name(idump): + return dump_prefix + f"{idump:07}" + ".phdump" + + +dump_helper = shamrock.utils.dump.ShamrockDumpHandleHelper(model, dump_prefix) + +# %% +# Load the last dump if it exists, setup otherwise + + +def setup_model(): + global disc_mass + + # Generate the default config + cfg = model.gen_default_config() + cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV) + cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q, r0=r0) + + cfg.add_kill_sphere(center=(0, 0, 0), radius=bsize) # kill particles outside the simulation box + + cfg.set_units(codeu) + cfg.set_particle_mass(pmass) + # Set the CFL + cfg.set_cfl_cour(C_cour) + cfg.set_cfl_force(C_force) + + # Enable this to debug the neighbor counts + # cfg.set_show_neigh_stats(True) + + # Standard way to set the smoothing length (e.g. Price et al. 2018) + cfg.set_smoothing_length_density_based() + + cfg.set_save_dt_to_fields(True) + + # Standard density based smoothing lenght but with a neighbor count limit + # Use it if you have large slowdowns due to giant particles + # I recommend to use it if you have a circumbinary discs as the issue is very likely to happen + # cfg.set_smoothing_length_density_based_neigh_lim(500) + + cfg.set_save_dt_to_fields(True) + + # Set the solver config to be the one stored in cfg + model.set_solver_config(cfg) + + # Print the solver config + model.get_current_config().print_status() + + # Init the scheduler & fields + model.init_scheduler(scheduler_split_val, scheduler_merge_val) + + # Set the simulation box size + model.resize_simulation_box(bmin, bmax) + + # Create the setup + + setup = model.get_setup() + gen_disc = setup.make_generator_disc_mc( + part_mass=pmass, + disc_mass=disc_mass, + r_in=rin, + r_out=rout, + sigma_profile=sigma_profile, + H_profile=H_profile, + rot_profile=rot_profile, + cs_profile=cs_profile, + random_seed=666, + ) + + # Print the dot graph of the setup + print(gen_disc.get_dot()) + + # Apply the setup + setup.apply_setup(gen_disc) + + # correct the momentum and barycenter of the disc to 0 + analysis_momentum = shamrock.model_sph.analysisTotalMomentum(model=model) + total_momentum = analysis_momentum.get_total_momentum() + + if shamrock.sys.world_rank() == 0: + print(f"disc momentum = {total_momentum}") + + model.apply_momentum_offset((-total_momentum[0], -total_momentum[1], -total_momentum[2])) + + # Correct the barycenter before adding the sink + analysis_barycenter = shamrock.model_sph.analysisBarycenter(model=model) + barycenter, disc_mass = analysis_barycenter.get_barycenter() + + if shamrock.sys.world_rank() == 0: + print(f"disc barycenter = {barycenter}") + + model.apply_position_offset((-barycenter[0], -barycenter[1], -barycenter[2])) + + total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + + if shamrock.sys.world_rank() == 0: + print(f"disc momentum after correction = {total_momentum}") + + barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() + + if shamrock.sys.world_rank() == 0: + print(f"disc barycenter after correction = {barycenter}") + + if not np.allclose(total_momentum, 0.0): + raise RuntimeError("disc momentum is not 0") + if not np.allclose(barycenter, 0.0): + raise RuntimeError("disc barycenter is not 0") + + # now that the barycenter & momentum are 0, we can add the sink + model.add_sink( + center_mass, + (0, 0, 0), + (0, 0, 0), + center_racc, + is_torque_free=center_is_torque_free, + torque_boost_radius_fact=center_torque_boost_radius_fact, + ) + + # Run a single step to init the integrator and smoothing length of the particles + # Here the htolerance is the maximum factor of evolution of the smoothing length in each + # Smoothing length iterations, increasing it affect the performance negatively but increse the + # convergence rate of the smoothing length + # this is why we increase it temporely to 1.3 before lowering it back to 1.1 (default value) + # Note that both ``change_htolerances`` can be removed and it will work the same but would converge + # more slowly at the first timestep + + model.change_htolerances(coarse=1.3, fine=1.1) + model.timestep() + model.change_htolerances(coarse=1.1, fine=1.1) + + +dump_helper.load_last_dump_or(setup_model) + + +# %% +# On the fly analysis +def save_rho_integ(ext, arr_rho, iplot): + if shamrock.sys.world_rank() == 0: + metadata = {"extent": [-ext, ext, -ext, ext], "time": model.get_time()} + np.save(plot_folder + f"rho_integ_{iplot:07}.npy", arr_rho) + + with open(plot_folder + f"rho_integ_{iplot:07}.json", "w") as fp: + json.dump(metadata, fp) + + +def save_analysis_data(filename, key, value, ianalysis): + """Helper to save analysis data to a JSON file.""" + if shamrock.sys.world_rank() == 0: + filepath = os.path.join(analysis_folder, filename) + try: + with open(filepath, "r") as fp: + data = json.load(fp) + except (FileNotFoundError, json.JSONDecodeError): + data = {key: []} + data[key] = data[key][:ianalysis] + data[key].append({"t": model.get_time(), key: value}) + with open(filepath, "w") as fp: + json.dump(data, fp, indent=4) + + +from shamrock.utils.analysis import ( + AnalysisHelper, + ColumnDensityPlot, + PerfHistory, + SliceDtPart, +) + +perf_analysis = PerfHistory(model, analysis_folder, "perf_history") + +column_density_plot = ColumnDensityPlot( + model, + ext_r=rout * 1.5, + nx=1024, + ny=1024, + ex=(1, 0, 0), + ey=(0, 1, 0), + center=(0, 0, 0), + analysis_folder=analysis_folder, + analysis_prefix="rho_integ_normal", +) + +dt_part_slice_plot = SliceDtPart( + model, + ext_r=rout * 0.5 / (16.0 / 9.0), # aspect ratio of 16:9 + nx=1920, + ny=1080, + ex=(1, 0, 0), + ey=(0, 0, 1), + center=((rin + rout) / 2, 0, 0), + analysis_folder=analysis_folder, + analysis_prefix="dt_part_slice", +) + +profile_plots = AnalysisHelper( + analysis_folder=os.path.join(analysis_folder, "plots"), + analysis_prefix="density_profile", +) + + +def analysis(ianalysis): + column_density_plot.analysis_save(ianalysis) + dt_part_slice_plot.analysis_save(ianalysis) + + barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() + + total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + angular_momentum = shamrock.model_sph.analysisAngularMomentum( + model=model + ).get_angular_momentum() + + potential_energy = shamrock.model_sph.analysisEnergyPotential( + model=model + ).get_potential_energy() + + kinetic_energy = shamrock.model_sph.analysisEnergyKinetic(model=model).get_kinetic_energy() + + save_analysis_data("barycenter.json", "barycenter", barycenter, ianalysis) + save_analysis_data("disc_mass.json", "disc_mass", disc_mass, ianalysis) + save_analysis_data("total_momentum.json", "total_momentum", total_momentum, ianalysis) + save_analysis_data("angular_momentum.json", "angular_momentum", angular_momentum, ianalysis) + save_analysis_data("potential_energy.json", "potential_energy", potential_energy, ianalysis) + save_analysis_data("kinetic_energy.json", "kinetic_energy", kinetic_energy, ianalysis) + + sinks = model.get_sinks() + save_analysis_data("sinks.json", "sinks", sinks, ianalysis) + + perf_analysis.analysis_save(ianalysis) + + rho_field = model.compute_field("rho", "f64") + hpart_field = model.compute_field("hpart", "f64") + + def compute_r_field(): + + def internal(size: int, x: np.array, y: np.array, z: np.array) -> np.array: + r = np.sqrt(x**2 + y**2) + return r + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + ) + + return model.compute_field("custom", "f64", custom_getter) + + def compute_Lz_field(): + pmass = model.get_particle_mass() + hfact = model.get_hfact() + + def internal( + size: int, + x: np.array, + y: np.array, + z: np.array, + hpart: np.array, + vx: np.array, + vy: np.array, + vz: np.array, + ) -> np.array: + rho = pmass * (hfact / hpart) ** 3 + Lz = rho * (x * vy - y * vx) + return Lz + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + dic_out["hpart"], + dic_out["vxyz"][:, 0], + dic_out["vxyz"][:, 1], + dic_out["vxyz"][:, 2], + ) + + return model.compute_field("custom", "f64", custom_getter) + + r_field = compute_r_field() + Lz_field = compute_Lz_field() + + print(rho_field, r_field, Lz_field) + + x_min = 0.25 + x_max = rout * 1.5 + x_min_log = np.log10(x_min) + x_max_log = np.log10(x_max) + + bin_edges_x1d = np.logspace(x_min_log, x_max_log, 2049) + + rho_profile = shamrock.compute_histogram_convolve_x( + bin_edges=bin_edges_x1d, + x_field=r_field, + y_field=rho_field, + size_field=hpart_field, + do_average=True, + ) + + Lz_profile = shamrock.compute_histogram_convolve_x( + bin_edges=bin_edges_x1d, + x_field=r_field, + y_field=Lz_field, + size_field=hpart_field, + do_average=True, + ) + + data = { + "bin_edges_x1d": bin_edges_x1d, + "rho_profile": rho_profile, + "Lz_profile": Lz_profile, + "time": model.get_time(), + } + + profile_plots.analysis_save(ianalysis, data) + + +# %% +# Evolve the simulation +model.solver_logs_reset_cumulated_step_time() +model.solver_logs_reset_step_count() + +t_start = model.get_time() + +idump = 0 +iplot = 0 +istop = 0 +for ttarg in t_stop: + if ttarg >= t_start: + model.evolve_until(ttarg) + + if istop % dump_freq_stop == 0: + dump_helper.write_dump(idump, purge_old_dumps=True, keep_first=1, keep_last=3) + + # dump = model.make_phantom_dump() + # dump.save_dump(get_ph_dump_name(idump)) + + if istop % plot_freq_stop == 0: + analysis(iplot) + + if istop % dump_freq_stop == 0: + idump += 1 + + if istop % plot_freq_stop == 0: + iplot += 1 + + istop += 1 + +# %% +# Plot generation +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Load the on-the-fly analysis after the run to make the plots +# (everything in this section can be in another file) + +import matplotlib +import matplotlib.pyplot as plt + +face_on_render_kwargs = { + "x_unit": "au", + "y_unit": "au", + "time_unit": "year", + "x_label": "x", + "y_label": "y", +} + +sink_params = { + "sink_scale_factor": 1, + "sink_color": "green", + "sink_linewidth": 1, + "sink_fill": False, +} + +column_density_plot.render_all( + **face_on_render_kwargs, + field_unit="kg.m^-2", + field_label="$\\int \\rho \\, \\mathrm{{d}} z$", + vmin=1, + vmax=1e4, + norm="log", + **sink_params, +) + +dt_part_slice_plot.render_all( + **face_on_render_kwargs, + field_unit="year", + field_label="$\\Delta t$", + vmin=1e-4, + vmax=1, + norm="log", + contour_list=[1e-4, 1e-3, 1e-2, 1e-1, 1], + **sink_params, +) + + +def profile_plots_func(iplot, data): + + data = data.item() + + bin_edges_x1d = data["bin_edges_x1d"] + rho_profile = data["rho_profile"] + Lz_profile = data["Lz_profile"] + time = data["time"] + + bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 + + plt.figure(dpi=150) + plt.plot(bin_center, rho_profile) + plt.xlabel("r") + plt.ylabel("density") + + plt.xscale("log") + plt.yscale("log") + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=2) + plt.gca().add_artist(anchored_text) + + plt.savefig(profile_plots.analysis_prefix + f"rho_{iplot:07}.png") + plt.close() + + plt.figure(dpi=150) + plt.plot(bin_center, Lz_profile) + plt.xlabel("r") + plt.ylabel("Lz") + + plt.xscale("log") + plt.yscale("log") + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=2) + plt.gca().add_artist(anchored_text) + + plt.savefig(profile_plots.analysis_prefix + f"Lz_{iplot:07}.png") + plt.close() + + +profile_plots.render_all(profile_plots_func) + +# %% +# Make gif for the doc (plot_to_gif.py) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Convert PNG sequence to Image sequence in mpl + +# sphinx_gallery_multi_image = "single" + + +render_gif = True + +# %% +# Do it for rho integ +if render_gif: + ani = column_density_plot.render_gif(gif_filename="rho_integ.gif", save_animation=True) + if ani is not None: + plt.show() + + +# %% +# Make a gif from the plots +if render_gif and shamrock.sys.world_rank() == 0: + ani = dt_part_slice_plot.render_gif(gif_filename="dt_part_slice.gif", save_animation=True) + if ani is not None: + plt.show() + + +# %% +# helper function to load data from JSON files +def load_data_from_json(filename, key): + filepath = os.path.join(analysis_folder, filename) + with open(filepath, "r") as fp: + data = json.load(fp)[key] + t = [d["t"] for d in data] + values = [d[key] for d in data] + return t, values + + +# %% +# load the json file for barycenter +t, barycenter = load_data_from_json("barycenter.json", "barycenter") +barycenter_x = [d[0] for d in barycenter] +barycenter_y = [d[1] for d in barycenter] +barycenter_z = [d[2] for d in barycenter] + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, barycenter_x) +plt.plot(t, barycenter_y) +plt.plot(t, barycenter_z) +plt.xlabel("t") +plt.ylabel("barycenter") +plt.legend(["x", "y", "z"]) +plt.savefig(analysis_folder + "barycenter.png") +plt.show() + +# %% +# load the json file for disc_mass +t, disc_mass = load_data_from_json("disc_mass.json", "disc_mass") + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, disc_mass) +plt.xlabel("t") +plt.ylabel("disc_mass") +plt.savefig(analysis_folder + "disc_mass.png") +plt.show() + +# %% +# load the json file for total_momentum +t, total_momentum = load_data_from_json("total_momentum.json", "total_momentum") +total_momentum_x = [d[0] for d in total_momentum] +total_momentum_y = [d[1] for d in total_momentum] +total_momentum_z = [d[2] for d in total_momentum] + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, total_momentum_x) +plt.plot(t, total_momentum_y) +plt.plot(t, total_momentum_z) +plt.xlabel("t") +plt.ylabel("total_momentum") +plt.legend(["x", "y", "z"]) +plt.savefig(analysis_folder + "total_momentum.png") +plt.show() + + +# %% +# load the json file for total_momentum +t, angular_momentum = load_data_from_json("angular_momentum.json", "angular_momentum") +angular_momentum_x = [d[0] - angular_momentum[0][0] for d in angular_momentum] +angular_momentum_y = [d[1] - angular_momentum[0][1] for d in angular_momentum] +angular_momentum_z = [d[2] - angular_momentum[0][2] for d in angular_momentum] + + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, angular_momentum_x) +plt.plot(t, angular_momentum_y) +plt.plot(t, angular_momentum_z) +plt.xlabel("t") +plt.ylabel(r"$\mathrm{L} - \mathrm{L}(t=0)$") +plt.legend(["x", "y", "z"]) +plt.savefig(analysis_folder + "angular_momentum.png") + +# %% +# load the json file for energies +t, potential_energy = load_data_from_json("potential_energy.json", "potential_energy") +_, kinetic_energy = load_data_from_json("kinetic_energy.json", "kinetic_energy") + +total_energy = [p + k for p, k in zip(potential_energy, kinetic_energy)] + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, potential_energy) +plt.plot(t, kinetic_energy) +plt.plot(t, total_energy) +plt.xlabel("t") +plt.ylabel("energy") +plt.legend(["potential_energy", "kinetic_energy", "total_energy"]) +plt.savefig(analysis_folder + "energies.png") +plt.show() + +# %% +# load the json file for sinks +t, sinks = load_data_from_json("sinks.json", "sinks") + +sinks_x = [d[0]["pos"][0] for d in sinks] +sinks_y = [d[0]["pos"][1] for d in sinks] +sinks_z = [d[0]["pos"][2] for d in sinks] + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, sinks_x, label="sink 0 (x)") +plt.plot(t, sinks_y, label="sink 0 (y)") +plt.plot(t, sinks_z, label="sink 0 (z)") +plt.xlabel("t") +plt.ylabel("sink position") +plt.legend() +plt.savefig(analysis_folder + "sinks.png") +plt.show() + +# %% +# Sink angular momentum +t, sinks = load_data_from_json("sinks.json", "sinks") + +sinks_lx = np.array([d[0]["angular_momentum"][0] for d in sinks]) +sinks_ly = np.array([d[0]["angular_momentum"][1] for d in sinks]) +sinks_lz = np.array([d[0]["angular_momentum"][2] for d in sinks]) + + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, sinks_lx, label="sink 0 (l_x)") +plt.plot(t, sinks_ly, label="sink 0 (l_y)") +plt.plot(t, sinks_lz, label="sink 0 (l_z)") +plt.xlabel("t") +plt.ylabel("sink spin") +plt.legend() +plt.savefig(analysis_folder + "sink_angular_momentum.png") +plt.show() + +# %% +# Sink to barycenter distance +t, sinks = load_data_from_json("sinks.json", "sinks") +_, barycenter = load_data_from_json("barycenter.json", "barycenter") + +barycenter_x = np.array([d[0] for d in barycenter]) +barycenter_y = np.array([d[1] for d in barycenter]) +barycenter_z = np.array([d[2] for d in barycenter]) + +sinks_x = np.array([d[0]["pos"][0] for d in sinks]) +sinks_y = np.array([d[0]["pos"][1] for d in sinks]) +sinks_z = np.array([d[0]["pos"][2] for d in sinks]) + + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, sinks_x - barycenter_x, label="sink 0 (x)") +plt.plot(t, sinks_y - barycenter_y, label="sink 0 (y)") +plt.plot(t, sinks_z - barycenter_z, label="sink 0 (z)") +plt.xlabel("t") +plt.ylabel("sink pos - barycenter pos") +plt.legend() +plt.savefig(analysis_folder + "sink_to_barycenter_distance.png") +plt.show() + + +# %% +# Plot the performance history (Switch close_plots to True if doing a long run) +perf_analysis.plot_perf_history(close_plots=False) +plt.show() diff --git a/examples/production/make_plots.py b/examples/production/make_plots.py new file mode 100644 index 0000000000..f4736b6896 --- /dev/null +++ b/examples/production/make_plots.py @@ -0,0 +1,209 @@ +cases = { + "racc=0.8": { + "path": "_to_trash/circular_disc_sink_False_0.8_10000000", + "is_torque_free": False, + "racc": 0.8, + }, + "racc=0.5": { + "path": "_to_trash/circular_disc_sink_False_0.5_10000000", + "is_torque_free": False, + "racc": 0.5, + }, + "racc=0.25": { + "path": "_to_trash/circular_disc_sink_False_0.25_10000000", + "is_torque_free": False, + "racc": 0.25, + }, + "racc=0.8 (torque-free)": { + "path": "_to_trash/circular_disc_sink_True_0.8_10000000", + "is_torque_free": True, + "racc": 0.8, + }, + "racc=0.5 (torque-free)": { + "path": "_to_trash/circular_disc_sink_True_0.5_10000000", + "is_torque_free": True, + "racc": 0.5, + }, + "racc=0.25 (torque-free)": { + "path": "_to_trash/circular_disc_sink_True_0.25_10000000", + "is_torque_free": True, + "racc": 0.25, + }, +} + +import json +import os + +import matplotlib.pyplot as plt +import numpy as np +from shamrock.utils.analysis import ( + AnalysisHelper, +) + + +def render_rho_profiles(cases): + helpers = {} + for key, case in cases.items(): + helpers[key] = {} + helpers[key]["helper"] = AnalysisHelper( + analysis_folder=os.path.join(cases[key]["path"], "analysis", "plots"), + analysis_prefix="density_profile", + ) + + for k in helpers: + helpers[k]["list_analysis_id"] = helpers[k]["helper"].get_list_analysis_id() + + profile_list_analysis_id1 = helpers["racc=0.8"]["list_analysis_id"] + + for iplot in profile_list_analysis_id1: + ref_key = None + data = {} + for k in helpers: + if iplot in helpers[k]["list_analysis_id"]: + data[k] = helpers[k]["helper"].load_analysis(iplot).item() + if ref_key is None: + ref_key = k + + time = data[ref_key]["time"] + bin_edges_x1d = data[ref_key]["bin_edges_x1d"] + bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 + + plt.figure(dpi=150) + for k, v in data.items(): + plt.plot(bin_center, v["rho_profile"], label=k) + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=1) + plt.gca().add_artist(anchored_text) + + plt.xscale("log") + plt.yscale("log") + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) + plt.legend(loc="lower left") + plt.savefig(f"_to_trash/compare_torque_free_sink_rho_profile_{iplot:07}.png") + plt.close() + + +def render_Lz_profiles(cases): + helpers = {} + for key, case in cases.items(): + helpers[key] = {} + helpers[key]["helper"] = AnalysisHelper( + analysis_folder=os.path.join(cases[key]["path"], "analysis", "plots"), + analysis_prefix="density_profile", + ) + + for k in helpers: + helpers[k]["list_analysis_id"] = helpers[k]["helper"].get_list_analysis_id() + + profile_list_analysis_id1 = helpers["racc=0.8"]["list_analysis_id"] + + for iplot in profile_list_analysis_id1: + ref_key = None + data = {} + for k in helpers: + if iplot in helpers[k]["list_analysis_id"]: + data[k] = helpers[k]["helper"].load_analysis(iplot).item() + if ref_key is None: + ref_key = k + + time = data[ref_key]["time"] + bin_edges_x1d = data[ref_key]["bin_edges_x1d"] + bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 + + plt.figure(dpi=150) + for k, v in data.items(): + plt.plot(bin_center, v["Lz_profile"], label=k) + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=1) + plt.gca().add_artist(anchored_text) + + plt.xscale("log") + plt.yscale("log") + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) + plt.legend(loc="lower left") + plt.savefig(f"_to_trash/compare_torque_free_sink_Lz_profile_{iplot:07}.png") + plt.close() + + +def compare_acc_rate(cases): + curves = {} + for key, case in cases.items(): + disc_mass_filename = os.path.join(cases[key]["path"], "analysis", "disc_mass.json") + + ret = {} + with open(disc_mass_filename, "r") as fp: + disc_mass_data = json.load(fp)["disc_mass"] + ret["t"] = [d["t"] for d in disc_mass_data] + ret["disc_mass"] = [d["disc_mass"] for d in disc_mass_data] + + # time derivative of the disc mass + ret["disc_mass_rate"] = np.diff(ret["disc_mass"]) / np.diff(ret["t"]) + curves[key] = ret + + plt.figure(dpi=150) + for k, v in curves.items(): + plt.plot(v["t"], v["disc_mass"], label=k) + + # plt.xscale("log") + # plt.yscale("log") + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) + plt.legend(loc="lower left") + plt.savefig("_to_trash/compare_torque_free_sink_acc_rate_disc_mass.png") + plt.close() + + plt.figure(dpi=150) + for k, v in curves.items(): + plt.plot(v["t"][1:], v["disc_mass_rate"], label=k) + + # plt.xscale("log") + # plt.yscale("log") + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) + plt.legend(loc="lower right") + plt.savefig("_to_trash/compare_torque_free_sink_acc_rate_diff.png") + plt.close() + + +def compare_dsink_rate(cases): + + def get_r(pos): + return np.linalg.norm(pos) + + curves = {} + for key, case in cases.items(): + disc_mass_filename = os.path.join(cases[key]["path"], "analysis", "sinks.json") + + ret = {} + with open(disc_mass_filename, "r") as fp: + sinks_dat = json.load(fp)["sinks"] + ret["t"] = [d["t"] for d in sinks_dat] + ret["r"] = [get_r(d["sinks"][0]["pos"]) for d in sinks_dat] + + curves[key] = ret + + plt.figure(dpi=150) + for k, v in curves.items(): + plt.plot(v["t"], v["r"], label=k) + + # plt.xscale("log") + # plt.yscale("log") + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) + plt.legend(loc="lower left") + plt.savefig("_to_trash/compare_torque_free_sink_r.png") + plt.close() + + +compare_dsink_rate(cases) +compare_acc_rate(cases) +render_rho_profiles(cases) +render_Lz_profiles(cases) diff --git a/examples/sph/run_compare_torque_free_sink.py b/examples/sph/run_compare_torque_free_sink.py new file mode 100644 index 0000000000..a0c952708c --- /dev/null +++ b/examples/sph/run_compare_torque_free_sink.py @@ -0,0 +1,579 @@ +""" +Production run: Circular disc & central sink particle +===================================================== + +This example demonstrates how to run a smoothed particle hydrodynamics (SPH) +simulation of a circular disc orbiting around a central sink. + +The simulation models: + +- A central star with a given mass and accretion radius +- A gaseous disc with specified mass, inner/outer radii, and vertical structure +- Artificial viscosity for angular momentum transport +- Locally isothermal equation of state + +Also this simulation feature rolling dumps (see `purge_old_dumps` function) to save disk space. + +This example is the accumulation of 3 files in a single one to showcase the complete workflow. + +- The actual run script (runscript.py) +- Plot generation (make_plots.py) +- Animation from the plots (plot_to_gif.py) + +On a cluster or laptop, one can run the code as follows: + +.. code-block:: bash + + mpirun ./shamrock --sycl-cfg 0:0 --loglevel 1 --rscript runscript.py + + +then after the run is done (or while it is running), one can run the following to generate the plots: + +.. code-block:: bash + + python make_plots.py + + +""" + +# %% +# Runscript (runscript.py) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# The runscript is the actual simulation with on the fly analysis & rolling dumps + +import glob +import json +import os # for makedirs + +import numpy as np + +try: + from numba import njit + + _HAS_NUMBA = True +except ImportError: + _HAS_NUMBA = False + + +from shamrock.utils.analysis import ( + AnalysisHelper, + ColumnDensityPlot, + ColumnParticleCount, + PerfHistory, + SliceDensityPlot, + SliceDiffVthetaProfile, + SliceDtPart, + SliceVzPlot, + VerticalShearGradient, +) + +import shamrock + +# If we use the shamrock executable to run this script instead of the python interpreter, +# we should not initialize the system as the shamrock executable needs to handle specific MPI logic +if not shamrock.sys.is_initialized(): + shamrock.change_loglevel(1) + shamrock.sys.init("0:0") + +# %% +# Setup units + +si = shamrock.UnitSystem() +sicte = shamrock.Constants(si) +codeu = shamrock.UnitSystem( + unit_time=sicte.year(), + unit_length=sicte.au(), + unit_mass=sicte.sol_mass(), +) +ucte = shamrock.Constants(codeu) +G = ucte.G() + +# %% +# List parameters + + +# Resolution +Npart = 200000 + +# Domain decomposition parameters +scheduler_split_val = int(1.0e7) # split patches with more than 1e7 particles +scheduler_merge_val = scheduler_split_val // 16 + +# Dump and plot frequency and duration of the simulation +dump_freq_stop = 2 + +dt_stop = 1 +nstop = 300 + +# The list of times at which the simulation will pause for analysis / dumping +t_stop = [i * dt_stop for i in range(nstop + 1)] + + +# Disc parameter +disc_mass = 0.01 # sol mass +rout = 20.0 # au +rin = 1.0 # au +H_r_0 = 0.05 +q = 0.5 +p = 3.0 / 2.0 +r0 = 1.0 + +# Viscosity parameter +alpha_AV = 1.0e-3 / 0.08 +alpha_u = 1.0 +beta_AV = 2.0 + +# Integrator parameters +C_cour = 0.3 +C_force = 0.25 + + +# Sink parameters +center_mass = 1.0 + + +# Disc profiles +def sigma_profile(r): + sigma_0 = 1.0 # We do not care as it will be renormalized + return sigma_0 * (r / r0) ** (-p) + + +def kep_profile(r): + return (G * center_mass / r) ** 0.5 + + +def omega_k(r): + return kep_profile(r) / r + + +def cs_profile(r): + cs_in = (H_r_0 * r0) * omega_k(r0) + return ((r / r0) ** (-q)) * cs_in + + +# Deduced quantities +pmass = disc_mass / Npart + +bsize = rout * 2 +bmin = (-bsize, -bsize, -bsize) +bmax = (bsize, bsize, bsize) + +cs0 = cs_profile(r0) + + +def rot_profile(r): + return ((kep_profile(r) ** 2) - (2 * p + q) * cs_profile(r) ** 2) ** 0.5 + + +def H_profile(r): + H = cs_profile(r) / omega_k(r) + # fact = (2.**0.5) * 3. # factor taken from phantom, to fasten thermalizing + fact = 1.0 + return fact * H + + +def simulate_disc(center_racc, is_torque_free, sim_folder): + + dump_folder = sim_folder + "dump/" + analysis_folder = sim_folder + "analysis/" + plot_folder = analysis_folder + "plots/" + + dump_prefix = dump_folder + "dump_" + + # %% + # Create the dump directory if it does not exist + if shamrock.sys.world_rank() == 0: + os.makedirs(sim_folder, exist_ok=True) + os.makedirs(dump_folder, exist_ok=True) + os.makedirs(analysis_folder, exist_ok=True) + os.makedirs(plot_folder, exist_ok=True) + + # %% + # Utility functions and quantities deduced from the base one + + # %% + # Start the context + # The context holds the data of the code + # We then init the layout of the field (e.g. the list of fields used by the solver) + + ctx = shamrock.Context() + ctx.pdata_layout_new() + + # %% + # Attach a SPH model to the context + + model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4") + + # %% + # Dump handling + dump_helper = shamrock.utils.dump.ShamrockDumpHandleHelper(model, dump_prefix) + + # %% + # Load the last dump if it exists, setup otherwise + def setup_model(): + global disc_mass + + # Generate the default config + cfg = model.gen_default_config() + cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV) + cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q, r0=r0) + + cfg.add_kill_sphere( + center=(0, 0, 0), radius=bsize + ) # kill particles outside the simulation box + + cfg.set_units(codeu) + cfg.set_particle_mass(pmass) + # Set the CFL + cfg.set_cfl_cour(C_cour) + cfg.set_cfl_force(C_force) + + # Enable this to debug the neighbor counts + # cfg.set_show_neigh_stats(True) + + # Standard way to set the smoothing length (e.g. Price et al. 2018) + cfg.set_smoothing_length_density_based() + + cfg.set_save_dt_to_fields(True) + + # Standard density based smoothing lenght but with a neighbor count limit + # Use it if you have large slowdowns due to giant particles + # I recommend to use it if you have a circumbinary discs as the issue is very likely to happen + # cfg.set_smoothing_length_density_based_neigh_lim(500) + + cfg.set_save_dt_to_fields(True) + + # Set the solver config to be the one stored in cfg + model.set_solver_config(cfg) + + # Print the solver config + model.get_current_config().print_status() + + # Init the scheduler & fields + model.init_scheduler(scheduler_split_val, scheduler_merge_val) + + # Set the simulation box size + model.resize_simulation_box(bmin, bmax) + + # Create the setup + + setup = model.get_setup() + gen_disc = setup.make_generator_disc_mc( + part_mass=pmass, + disc_mass=disc_mass, + r_in=rin, + r_out=rout, + sigma_profile=sigma_profile, + H_profile=H_profile, + rot_profile=rot_profile, + cs_profile=cs_profile, + random_seed=666, + ) + + # Print the dot graph of the setup + print(gen_disc.get_dot()) + + # Apply the setup + setup.apply_setup(gen_disc) + + # correct the momentum and barycenter of the disc to 0 + analysis_momentum = shamrock.model_sph.analysisTotalMomentum(model=model) + total_momentum = analysis_momentum.get_total_momentum() + + if shamrock.sys.world_rank() == 0: + print(f"disc momentum = {total_momentum}") + + model.apply_momentum_offset((-total_momentum[0], -total_momentum[1], -total_momentum[2])) + + # Correct the barycenter before adding the sink + analysis_barycenter = shamrock.model_sph.analysisBarycenter(model=model) + barycenter, disc_mass = analysis_barycenter.get_barycenter() + + if shamrock.sys.world_rank() == 0: + print(f"disc barycenter = {barycenter}") + + model.apply_position_offset((-barycenter[0], -barycenter[1], -barycenter[2])) + + total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + + if shamrock.sys.world_rank() == 0: + print(f"disc momentum after correction = {total_momentum}") + + barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() + + if shamrock.sys.world_rank() == 0: + print(f"disc barycenter after correction = {barycenter}") + + if not np.allclose(total_momentum, 0.0): + raise RuntimeError("disc momentum is not 0") + if not np.allclose(barycenter, 0.0): + raise RuntimeError("disc barycenter is not 0") + + # now that the barycenter & momentum are 0, we can add the sink + model.add_sink( + center_mass, (0, 0, 0), (0, 0, 0), center_racc, is_torque_free=is_torque_free + ) + + # Run a single step to init the integrator and smoothing length of the particles + # Here the htolerance is the maximum factor of evolution of the smoothing length in each + # Smoothing length iterations, increasing it affect the performance negatively but increse the + # convergence rate of the smoothing length + # this is why we increase it temporely to 1.3 before lowering it back to 1.1 (default value) + # Note that both ``change_htolerances`` can be removed and it will work the same but would converge + # more slowly at the first timestep + + model.change_htolerances(coarse=1.3, fine=1.1) + model.timestep() + model.change_htolerances(coarse=1.1, fine=1.1) + + dump_helper.load_last_dump_or(setup_model) + + def save_analysis_data(filename, key, value, ianalysis): + """Helper to save analysis data to a JSON file.""" + if shamrock.sys.world_rank() == 0: + filepath = os.path.join(analysis_folder, filename) + try: + with open(filepath, "r") as fp: + data = json.load(fp) + except (FileNotFoundError, json.JSONDecodeError): + data = {key: []} + data[key] = data[key][:ianalysis] + data[key].append({"t": model.get_time(), key: value}) + with open(filepath, "w") as fp: + json.dump(data, fp, indent=4) + + perf_analysis = PerfHistory(model, analysis_folder, "perf_history") + + column_density_plot = ColumnDensityPlot( + model, + ext_r=rout * 1.5, + nx=1024, + ny=1024, + ex=(1, 0, 0), + ey=(0, 1, 0), + center=(0, 0, 0), + analysis_folder=analysis_folder, + analysis_prefix="rho_integ_normal", + ) + + profile_plot = AnalysisHelper( + analysis_folder=os.path.join(analysis_folder, "plots"), + analysis_prefix="density_profile", + ) + + def analysis(ianalysis): + column_density_plot.analysis_save(ianalysis) + + barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() + + total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + angular_momentum = shamrock.model_sph.analysisAngularMomentum( + model=model + ).get_angular_momentum() + + potential_energy = shamrock.model_sph.analysisEnergyPotential( + model=model + ).get_potential_energy() + + kinetic_energy = shamrock.model_sph.analysisEnergyKinetic(model=model).get_kinetic_energy() + + save_analysis_data("barycenter.json", "barycenter", barycenter, ianalysis) + save_analysis_data("disc_mass.json", "disc_mass", disc_mass, ianalysis) + save_analysis_data("total_momentum.json", "total_momentum", total_momentum, ianalysis) + save_analysis_data("angular_momentum.json", "angular_momentum", angular_momentum, ianalysis) + save_analysis_data("potential_energy.json", "potential_energy", potential_energy, ianalysis) + save_analysis_data("kinetic_energy.json", "kinetic_energy", kinetic_energy, ianalysis) + + sinks = model.get_sinks() + save_analysis_data("sinks.json", "sinks", sinks, ianalysis) + + perf_analysis.analysis_save(ianalysis) + + #''' + rho_field = model.compute_field("rho", "f64") + hpart_field = model.compute_field("hpart", "f64") + + def internal(size: int, x: np.array, y: np.array, z: np.array) -> np.array: + r = np.sqrt(x**2 + y**2 + z**2) + return r + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + ) + + r_field = model.compute_field("custom", "f64", custom_getter) + + print(rho_field, r_field) + + x_min = 0.25 + x_max = rout * 1.5 + x_min_log = np.log10(x_min) + x_max_log = np.log10(x_max) + + bin_edges_x1d = np.logspace(x_min_log, x_max_log, 2049) + + histo = shamrock.compute_histogram( + bin_edges=bin_edges_x1d, + x_field=r_field, + y_field=rho_field, + do_average=True, + ) + + histo_convolve = shamrock.compute_histogram_convolve_x( + bin_edges=bin_edges_x1d, + x_field=r_field, + y_field=rho_field, + size_field=hpart_field, + do_average=True, + ) + + bin_edges_x = np.logspace(x_min_log, x_max_log, 1025) + bin_edges_y = np.logspace(-6, -3, 1025) + histo_top = shamrock.compute_histogram_2d( + bin_edges_x=bin_edges_x, + bin_edges_y=bin_edges_y, + x_field=r_field, + y_field=rho_field, + ) + histo_2d = np.array(histo_top).reshape(len(bin_edges_x) - 1, len(bin_edges_y) - 1) + + data = { + "bin_edges_x1d": bin_edges_x1d, + "bin_edges_x": bin_edges_x, + "bin_edges_y": bin_edges_y, + "histo": histo, + "histo_convolve": histo_convolve, + "histo_2d": histo_2d, + "time": model.get_time(), + } + + profile_plot.analysis_save(ianalysis, data) + + # %% + # Evolve the simulation + model.solver_logs_reset_cumulated_step_time() + model.solver_logs_reset_step_count() + + t_start = model.get_time() + + idump = 0 + iplot = 0 + istop = 0 + for ttarg in t_stop: + if ttarg >= t_start: + model.evolve_until(ttarg) + + if istop % dump_freq_stop == 0: + dump_helper.write_dump(idump, purge_old_dumps=True, keep_first=1, keep_last=3) + + # dump = model.make_phantom_dump() + # dump.save_dump(get_ph_dump_name(idump)) + + analysis(iplot) + + if istop % dump_freq_stop == 0: + idump += 1 + + iplot += 1 + istop += 1 + + +center_racc = 0.8 +sim_folder1 = f"_to_trash/circular_disc_sink_no_torque_free_{Npart}/" +simulate_disc(center_racc, False, sim_folder1) + +center_racc = 0.8 +sim_folder2 = f"_to_trash/circular_disc_sink_torque_free_{Npart}/" +simulate_disc(center_racc, True, sim_folder2) + +center_racc = 0.25 +sim_folder3 = f"_to_trash/circular_disc_sink_no_torque_free_0.25_{Npart}/" +simulate_disc(center_racc, False, sim_folder3) + +center_racc = 0.25 +sim_folder4 = f"_to_trash/circular_disc_sink_torque_free_0.25_{Npart}/" +simulate_disc(center_racc, True, sim_folder4) + +import matplotlib.pyplot as plt + + +def render_profiles(sim_folder1, sim_folder2, sim_folder3, sim_folder4): + profile_plot1 = AnalysisHelper( + analysis_folder=os.path.join(sim_folder1, "analysis", "plots"), + analysis_prefix="density_profile", + ) + profile_plot2 = AnalysisHelper( + analysis_folder=os.path.join(sim_folder2, "analysis", "plots"), + analysis_prefix="density_profile", + ) + profile_plot3 = AnalysisHelper( + analysis_folder=os.path.join(sim_folder3, "analysis", "plots"), + analysis_prefix="density_profile", + ) + profile_plot4 = AnalysisHelper( + analysis_folder=os.path.join(sim_folder4, "analysis", "plots"), + analysis_prefix="density_profile", + ) + profile_list_analysis_id1 = profile_plot1.get_list_analysis_id() + profile_list_analysis_id2 = profile_plot2.get_list_analysis_id() + profile_list_analysis_id3 = profile_plot3.get_list_analysis_id() + profile_list_analysis_id4 = profile_plot4.get_list_analysis_id() + + for iplot in profile_list_analysis_id1: + has_iplot_data_1 = iplot in profile_list_analysis_id1 + has_iplot_data_2 = iplot in profile_list_analysis_id2 + has_iplot_data_3 = iplot in profile_list_analysis_id3 + has_iplot_data_4 = iplot in profile_list_analysis_id4 + + if ( + not has_iplot_data_1 + or not has_iplot_data_2 + or not has_iplot_data_3 + or not has_iplot_data_4 + ): + continue + + data1 = profile_plot1.load_analysis(iplot).item() + data2 = profile_plot2.load_analysis(iplot).item() + data3 = profile_plot3.load_analysis(iplot).item() + data4 = profile_plot4.load_analysis(iplot).item() + + time = data1["time"] + + histo_convolve_torque_free = data1["histo_convolve"] + histo_convolve_not_torque_free = data2["histo_convolve"] + histo_convolve_torque_free_025 = data3["histo_convolve"] + histo_convolve_not_torque_free_025 = data4["histo_convolve"] + + bin_edges_x1d = data1["bin_edges_x1d"] + bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 + + plt.figure(dpi=150) + plt.plot(bin_center, histo_convolve_torque_free, label="racc=0.8") + plt.plot(bin_center, histo_convolve_not_torque_free, label="racc=0.8 torque free") + plt.plot(bin_center, histo_convolve_torque_free_025, label="racc=0.25") + plt.plot(bin_center, histo_convolve_not_torque_free_025, label="racc=0.25 torque free") + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=1) + plt.gca().add_artist(anchored_text) + + plt.xscale("log") + plt.yscale("log") + plt.ylim(1e-6, 1e-3) + plt.xlim(0.1, 10) + plt.legend(loc="upper left") + plt.savefig(f"_to_trash/compare_torque_free_sink_{iplot:07}.png") + plt.close() + + +render_profiles(sim_folder1, sim_folder2, sim_folder3, sim_folder4) diff --git a/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py b/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py new file mode 100644 index 0000000000..220a47dce4 --- /dev/null +++ b/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py @@ -0,0 +1,176 @@ +import numpy as np + +import shamrock.sys + +from .StandardPlotHelper import StandardPlotHelper + +try: + from numba import njit + + _HAS_NUMBA = True +except ImportError: + _HAS_NUMBA = False + + +def SliceAngularMomentum( + model, + ext_r, + nx, + ny, + ex, + ey, + center, + analysis_folder, + analysis_prefix, + do_normalization=True, + min_normalization=1e-9, + Lprojection=[0.0, 0.0, 1.0], +): + def compute_angular_mom(helper): + if _HAS_NUMBA: + if shamrock.sys.world_rank() == 0: + print("Using numba for angular momentum in SliceAngularMomentum") + + pmass = model.get_particle_mass() + hfact = model.get_hfact() + + def internal( + size: int, + hpart: np.array, + x: np.array, + y: np.array, + z: np.array, + vx: np.array, + vy: np.array, + vz: np.array, + ) -> np.array: + + rho = pmass * (hfact / hpart) ** 3 + + r = np.stack([x, y, z], axis=-1) + v = np.stack([vx, vy, vz], axis=-1) + L = rho[:, None] * np.cross(r, v) + + # project L onto the Lprojection vector + L = np.sum(L * Lprojection, axis=-1) + + return L + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["hpart"], + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + dic_out["vxyz"][:, 0], + dic_out["vxyz"][:, 1], + dic_out["vxyz"][:, 2], + ) + + arr_L = helper.slice_render( + "custom", + "f64", + do_normalization, + min_normalization, + custom_getter=custom_getter, + ) + + return arr_L + + return StandardPlotHelper( + model, + ext_r, + nx, + ny, + ex, + ey, + center, + analysis_folder, + analysis_prefix, + compute_function=compute_angular_mom, + ) + + +def ColumnAverageAngularMomentum( + model, + ext_r, + nx, + ny, + ex, + ey, + center, + analysis_folder, + analysis_prefix, + min_normalization=1e-9, + Lprojection=[0.0, 0.0, 1.0], +): + def compute_angular_mom(helper): + if _HAS_NUMBA: + if shamrock.sys.world_rank() == 0: + print("Using numba for angular momentum in SliceAngularMomentum") + + pmass = model.get_particle_mass() + hfact = model.get_hfact() + + def internal( + size: int, + hpart: np.array, + x: np.array, + y: np.array, + z: np.array, + vx: np.array, + vy: np.array, + vz: np.array, + ) -> np.array: + + rho = pmass * (hfact / hpart) ** 3 + + r = np.stack([x, y, z], axis=-1) + v = np.stack([vx, vy, vz], axis=-1) + L = rho[:, None] * np.cross(r, v) + + # project L onto the Lprojection vector + L = np.sum(L * Lprojection, axis=-1) + + return L + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["hpart"], + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + dic_out["vxyz"][:, 0], + dic_out["vxyz"][:, 1], + dic_out["vxyz"][:, 2], + ) + + arr_L = helper.column_average_render( + "custom", + "f64", + min_normalization, + custom_getter=custom_getter, + ) + + return arr_L + + return StandardPlotHelper( + model, + ext_r, + nx, + ny, + ex, + ey, + center, + analysis_folder, + analysis_prefix, + compute_function=compute_angular_mom, + ) diff --git a/src/pylib/shamrock/utils/analysis/__init__.py b/src/pylib/shamrock/utils/analysis/__init__.py index 3772390f1f..d623593586 100644 --- a/src/pylib/shamrock/utils/analysis/__init__.py +++ b/src/pylib/shamrock/utils/analysis/__init__.py @@ -22,5 +22,10 @@ SliceBVerticalShearGradient, ) +from .AngularMomentumPlots import ( + SliceAngularMomentum, + ColumnAverageAngularMomentum, +) + # Performance analysis from .PerfHistory import PerfHistory diff --git a/src/shammath/include/shammath/matrix_op.hpp b/src/shammath/include/shammath/matrix_op.hpp index dde113f365..1d22ffeeeb 100644 --- a/src/shammath/include/shammath/matrix_op.hpp +++ b/src/shammath/include/shammath/matrix_op.hpp @@ -334,6 +334,35 @@ namespace shammath { } } + /** + * @brief Compute the determinant of a 3x3 matrix. + * + * @param input The input matrix to invert. + * @return the determinant + */ + template + inline T mat_det_33( + const std::mdspan, Layout, Accessor> &input) { + + T &a00 = input(0, 0); + T &a10 = input(1, 0); + T &a20 = input(2, 0); + + T &a01 = input(0, 1); + T &a11 = input(1, 1); + T &a21 = input(2, 1); + + T &a02 = input(0, 2); + T &a12 = input(1, 2); + T &a22 = input(2, 2); + + T det + = (-a02 * a11 * a20 + a01 * a12 * a20 + a02 * a10 * a21 - a00 * a12 * a21 + - a01 * a10 * a22 + a00 * a11 * a22); + + return det; + } + /** * @brief Compute the inverse of a 3x3 matrix. * diff --git a/src/shammodels/sph/include/shammodels/sph/Model.hpp b/src/shammodels/sph/include/shammodels/sph/Model.hpp index eb78e1feee..061c6d0002 100644 --- a/src/shammodels/sph/include/shammodels/sph/Model.hpp +++ b/src/shammodels/sph/include/shammodels/sph/Model.hpp @@ -173,15 +173,37 @@ namespace shammodels::sph { Tscal q, std::mt19937 eng); - inline void add_sink(Tscal mass, Tvec pos, Tvec velocity, Tscal accretion_radius) { + inline void add_sink( + Tscal mass, + Tvec pos, + Tvec velocity, + Tscal accretion_radius, + bool is_torque_free, + Tscal torque_boost_radius_fact) { if (solver.storage.sinks.is_empty()) { solver.storage.sinks.set({}); } - shamlog_debug_ln("SPH", "add sink :", mass, pos, velocity, accretion_radius); + shamlog_debug_ln( + "SPH", + "add sink :", + mass, + pos, + velocity, + accretion_radius, + is_torque_free, + torque_boost_radius_fact); solver.storage.sinks.get().push_back( - {pos, velocity, {}, {}, mass, {}, accretion_radius}); + {pos, + velocity, + {}, + {}, + mass, + {}, + accretion_radius, + is_torque_free, + torque_boost_radius_fact}); } template diff --git a/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp b/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp index 0a58f2232d..9eb0e95be7 100644 --- a/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp @@ -34,6 +34,9 @@ namespace shammodels::sph { Tscal mass; Tvec angular_momentum; Tscal accretion_radius; + + bool is_torque_free; + Tscal torque_boost_radius_fact; }; template @@ -50,6 +53,8 @@ namespace shammodels::sph { {"mass", p.mass}, {"angular_momentum", p.angular_momentum}, {"accretion_radius", p.accretion_radius}, + {"is_torque_free", p.is_torque_free}, + {"torque_boost_radius_fact", p.torque_boost_radius_fact}, }; } @@ -65,6 +70,8 @@ namespace shammodels::sph { j.at("mass").get_to(p.mass); j.at("angular_momentum").get_to(p.angular_momentum); j.at("accretion_radius").get_to(p.accretion_radius); + j.at("is_torque_free").get_to(p.is_torque_free); + j.at("torque_boost_radius_fact").get_to(p.torque_boost_radius_fact); } } // namespace shammodels::sph diff --git a/src/shammodels/sph/src/Model.cpp b/src/shammodels/sph/src/Model.cpp index 80b2473925..589515eb83 100644 --- a/src/shammodels/sph/src/Model.cpp +++ b/src/shammodels/sph/src/Model.cpp @@ -1437,7 +1437,9 @@ void shammodels::sph::Model::init_from_phantom_dump( mass[i], {xsink[i], ysink[i], zsink[i]}, {vxsink[i], vysink[i], vzsink[i]}, - Racc[i]); + Racc[i], + false, + 2.0); } } } diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index 5796661fb8..f5b08d9fe1 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -17,14 +17,22 @@ #include "shambase/DistributedData.hpp" #include "shambase/narrowing.hpp" +#include "shambase/string.hpp" #include "shamalgs/collective/reduction.hpp" #include "shamalgs/details/numeric/numeric.hpp" +#include "shamalgs/numeric.hpp" #include "shamalgs/primitives/reduction.hpp" +#include "shamalgs/reduction.hpp" #include "shambackends/DeviceBuffer.hpp" +#include "shambackends/fmt_bindings/fmt_defs.hpp" #include "shambackends/kernel_call.hpp" +#include "shambackends/math.hpp" #include "shamcomm/logs.hpp" +#include "shammath/matrix.hpp" +#include "shammath/matrix_op.hpp" #include "shammath/sphkernels.hpp" #include "shammodels/sph/modules/SinkParticlesUpdate.hpp" +#include "shamsys/legacy/log.hpp" #include template class SPHKernel> @@ -50,6 +58,84 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par std::vector &sink_parts = storage.sinks.get(); + /* // Uncomment this if you want to compute total momentum to check some stuff + auto get_tot_mom = [&]() -> Tvec { + Tvec total_momentum = {}; + + sham::DeviceBuffer total_momentum_part(0, dev_sched); + + scheduler().for_each_patchdata_nonempty( + [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer + &pdat) { u32 len = pdat.get_obj_cnt(); + + total_momentum_part.resize(len); + + sham::DeviceBuffer &vxyz_buf = + pdat.get_field_buf_ref(ivxyz); + + sham::kernel_call( + q, + sham::MultiRef{vxyz_buf}, + sham::MultiRef{total_momentum_part}, + len, + [gpart_mass]( + u32 i, const Tvec *__restrict vxyz, Tvec *__restrict + total_momentum_part) { total_momentum_part[i] = gpart_mass * vxyz[i]; + }); + + total_momentum += shamalgs::primitives::sum(dev_sched, + total_momentum_part, 0, len); + }); + + Tvec tot_total_momentum = + shamalgs::collective::allreduce_sum(total_momentum); + + for (auto &sink : sink_parts) { + tot_total_momentum += sink.mass * sink.velocity; + } + + return tot_total_momentum; + }; + */ + + auto get_ang_mom = [&]() -> Tvec { + Tvec angular_momentum = {}; + + sham::DeviceBuffer total_momentum_part(0, dev_sched); + + scheduler().for_each_patchdata_nonempty([&](const shamrock::patch::Patch p, + shamrock::patch::PatchDataLayer &pdat) { + u32 len = pdat.get_obj_cnt(); + + total_momentum_part.resize(len); + + sham::DeviceBuffer &xyz_buf = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &vxyz_buf = pdat.get_field_buf_ref(ivxyz); + sham::kernel_call( + q, + sham::MultiRef{xyz_buf, vxyz_buf}, + sham::MultiRef{total_momentum_part}, + len, + [gpart_mass]( + u32 i, + const Tvec *__restrict xyz, + const Tvec *__restrict vxyz, + Tvec *__restrict total_momentum_part) { + total_momentum_part[i] = gpart_mass * sycl::cross(xyz[i], vxyz[i]); + }); + + angular_momentum += shamalgs::primitives::sum(dev_sched, total_momentum_part, 0, len); + }); + + Tvec tot_angular_momentum = shamalgs::collective::allreduce_sum(angular_momentum); + + for (auto &sink : sink_parts) { + tot_angular_momentum + += sink.mass * sycl::cross(sink.pos, sink.velocity) + sink.angular_momentum; + } + return tot_angular_momentum; + }; + u32 sink_id = 0; bool had_accretion = false; std::string log = "sink accretion :"; @@ -162,8 +248,9 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par accretion_mr[id_a] = gpart_mass * r; accretion_ma[id_a] = gpart_mass * a; - // dirty trick to account for the residual acceleration in the spin. This - // allows us to maitain a much better angular momentum conservation. + // dirty trick to account for the residual acceleration in the + // spin. This allows us to maitain a much better angular momentum + // conservation. v += a * dt / 2; accretion_l[id_a] = gpart_mass * sycl::cross(r - r_sink, v - v_sink); }); @@ -241,6 +328,176 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par }); } + // now time for torque free torque restitution + for (Sink &s : sink_parts) { + if (s.is_torque_free && sham::dot(s.angular_momentum, s.angular_momentum) > 0) { + // boost particles in a radius of 4 racc to empty the sink angular_momentum + + struct BoostWeight { + Tscal radius; + Tscal weight(Tscal r) const { + if (r < radius * Kernel::Rkern) { + return Kernel::W_3d(r, radius); + } else { + return 0; + } + } + }; + + BoostWeight weight_func{ + s.accretion_radius * s.torque_boost_radius_fact / Kernel::Rkern}; + + Tvec r_sink = s.pos; + + using mat = shammath::mat; + + mat I_rank{}; + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + + sham::DeviceBuffer weighted_intertia(9 * Nobj, dev_sched); + sham::kernel_call( + q, + sham::MultiRef{buf_xyz}, + sham::MultiRef{weighted_intertia}, + Nobj, + [r_sink, gpart_mass, Nobj, weight_func]( + u32 id_a, const Tvec *__restrict xyz, Tscal *weighted_intertia) { + Tvec r_a = xyz[id_a] - r_sink; + Tscal r_a2 = sham::dot(r_a, r_a); + + Tscal w = gpart_mass * weight_func.weight(sycl::sqrt(r_a2)); + + if (w == 0) { + return; + } + + mat I{}; + + { // I = w [Id r_a^2 - r_a r_a^T] + auto Imd = I.get_mdspan(); + + Imd(0, 0) += w * (r_a.y() * r_a.y() + r_a.z() * r_a.z()); + Imd(1, 1) += w * (r_a.x() * r_a.x() + r_a.z() * r_a.z()); + Imd(2, 2) += w * (r_a.x() * r_a.x() + r_a.y() * r_a.y()); + + Imd(0, 1) -= w * r_a.x() * r_a.y(); + Imd(0, 2) -= w * r_a.x() * r_a.z(); + Imd(1, 2) -= w * r_a.y() * r_a.z(); + + // symmetry time ! + Imd(1, 0) = Imd(0, 1); + Imd(2, 0) = Imd(0, 2); + Imd(2, 1) = Imd(1, 2); + } + +#pragma unroll + for (u32 i = 0; i < 9; i++) { + weighted_intertia[Nobj * i + id_a] = I.data[i]; + } + }); + + for (u32 i = 0; i < 9; i++) { + I_rank.data[i] += shamalgs::primitives::sum( + dev_sched, weighted_intertia, Nobj * i, Nobj * (i + 1)); + } + }); + + mat weighted_intertia_sum{}; + + for (u32 i = 0; i < 9; i++) { + weighted_intertia_sum.data[i] = shamalgs::collective::allreduce_sum(I_rank.data[i]); + } + + mat I_inv{}; + if (shammath::mat_det_33(weighted_intertia_sum.get_mdspan()) != 0) { + shammath::mat_inv_33(weighted_intertia_sum.get_mdspan(), I_inv.get_mdspan()); + } + + shammath::vec delta_l{}; + delta_l.data[0] = s.angular_momentum.x(); + delta_l.data[1] = s.angular_momentum.y(); + delta_l.data[2] = s.angular_momentum.z(); + + shammath::vec delta_w{}; + shammath::mat_prod( + I_inv.get_mdspan(), delta_l.get_mdspan_mat_col(), delta_w.get_mdspan_mat_col()); + + Tvec vec_delta_w = {delta_w.data[0], delta_w.data[1], delta_w[2]}; + + if (sham::has_nan_or_inf(vec_delta_w)) { + logger::warn_ln( + "Sink", + shambase::format( + "Nan or inf detected after inertia matrix inversion\n" + " delta_l = {}\n" + " delta_w = {}\n" + " weighted_intertia_sum = {}\n" + " I_inv = {}\n" + "I will therefore skip the boost step (delta_w = 0)", + delta_l.data, + delta_w.data, + weighted_intertia_sum.data, + I_inv.data)); + + delta_w = {}; + vec_delta_w = {delta_w.data[0], delta_w.data[1], delta_w[2]}; + } + + logger::raw_ln("delta l =", delta_l.data[0], delta_l.data[1], delta_l.data[2]); + logger::raw_ln("delta w =", delta_w.data[0], delta_w.data[1], delta_w.data[2]); + + Tvec lboost_sum_rank{}; + Tvec pboost_sum_rank{}; + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); + + sham::DeviceBuffer l_boost(Nobj, dev_sched); + sham::DeviceBuffer p_boost(Nobj, dev_sched); + + sham::kernel_call( + q, + sham::MultiRef{buf_xyz}, + sham::MultiRef{buf_vxyz, l_boost, p_boost}, + Nobj, + [r_sink, delta_w = vec_delta_w, weight_func, gpart_mass]( + u32 id_a, + const Tvec *__restrict xyz, + Tvec *__restrict vxyz, + Tvec *__restrict lboost, + Tvec *__restrict pboost) { + Tvec r_a = xyz[id_a] - r_sink; + Tscal r_a2 = sham::dot(r_a, r_a); + + Tscal W = weight_func.weight(sycl::sqrt(r_a2)); + + Tvec v_boost = W * sycl::cross(delta_w, r_a); + + vxyz[id_a] += v_boost; + lboost[id_a] += gpart_mass * sycl::cross(r_a, v_boost); + pboost[id_a] += gpart_mass * v_boost; + }); + + lboost_sum_rank += shamalgs::primitives::sum(dev_sched, l_boost, 0, Nobj); + pboost_sum_rank += shamalgs::primitives::sum(dev_sched, p_boost, 0, Nobj); + }); + + Tvec lboost = shamalgs::collective::allreduce_sum(lboost_sum_rank); + Tvec pboost = shamalgs::collective::allreduce_sum(pboost_sum_rank); + + logger::raw_ln("lboost =", lboost, "pboost =", pboost); + + s.angular_momentum -= lboost; + s.velocity -= pboost / s.mass; + } + } + if (shamcomm::world_rank() == 0 && had_accretion) { logger::info_ln("sph::Sink", log); } @@ -352,9 +609,10 @@ void shammodels::sph::modules::SinkParticlesUpdate::compute_sph Tvec force = G * delta / (d * d * d); - // This is a hack to avoid the sink kaboom effect - // when the particle is being advected close to the sink before - // being accreted + // This is a hack to avoid the sink + // kaboom effect when the particle is + // being advected close to the sink + // before being accreted if (d < sink_racc) { force = {0, 0, 0}; } diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index 59f97f762f..6ac27434eb 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -825,7 +825,15 @@ void add_instance(py::module &m, std::string name_config, std::string name_model return self.solver.solver_config; }) .def("set_solver_config", &T::set_solver_config) - .def("add_sink", &T::add_sink) + .def( + "add_sink", + &T::add_sink, + py::arg("mass"), + py::arg("pos"), + py::arg("velocity"), + py::arg("accretion_radius"), + py::arg("is_torque_free") = false, + py::arg("torque_boost_radius_fact") = 2.0) .def( "get_sinks", [](T &self) { @@ -834,13 +842,15 @@ void add_instance(py::module &m, std::string name_config, std::string name_model if (!self.solver.storage.sinks.is_empty()) { for (auto &sink : self.solver.storage.sinks.get()) { py::dict sink_dic; - sink_dic["pos"] = sink.pos; - sink_dic["velocity"] = sink.velocity; - sink_dic["sph_acceleration"] = sink.sph_acceleration; - sink_dic["ext_acceleration"] = sink.ext_acceleration; - sink_dic["mass"] = sink.mass; - sink_dic["angular_momentum"] = sink.angular_momentum; - sink_dic["accretion_radius"] = sink.accretion_radius; + sink_dic["pos"] = sink.pos; + sink_dic["velocity"] = sink.velocity; + sink_dic["sph_acceleration"] = sink.sph_acceleration; + sink_dic["ext_acceleration"] = sink.ext_acceleration; + sink_dic["mass"] = sink.mass; + sink_dic["angular_momentum"] = sink.angular_momentum; + sink_dic["accretion_radius"] = sink.accretion_radius; + sink_dic["is_torque_free"] = sink.is_torque_free; + sink_dic["torque_boost_radius_fact"] = sink.torque_boost_radius_fact; list_out.append(sink_dic); } }