noe/mass_inflow_gravity
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

add final scripts

Browse Source
This commit is contained in:
Noé Brucy
2025-12-18 16:51:11 +01:00
commit 49a790ed8c
5 changed files with 3237 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

118
extract_tracers.py Normal file
View File

@@ -0,0 +1,118 @@
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import osyris
from scipy.io import FortranFile
from scipy.spatial import KDTree
def link_tracer(data):
mask_tracer = data["part"]["family"].values == 0
pos_cells = data["mesh"]["position"].to("pc").values
pos_tracers = data["part"]["position"].to("pc").values[mask_tracer]
kdtree = KDTree(pos_cells)
dist, ind_cell = kdtree.query(pos_tracers, k=1, p=np.inf)
return ind_cell
def get_cell_data(data, ind_cell, key, unit):
cell_data = data["mesh"][key].to(unit).values[ind_cell]
return cell_data
def get_tracer_info_osyris(path, snap):
data = osyris.RamsesDataset(snap, path=os.path.expanduser(path)).load()
mesh = data["mesh"]
mesh["temperature"] = ((mesh["pressure"] / mesh["density"] ) * (1.4 *osyris.units("proton_mass") / osyris.units("boltzmann_constant"))).to("K")
mesh["n"] = (mesh["density"] / (1.4 *osyris.units("proton_mass"))).to("cm**-3")
mesh["mass"] = (mesh["density"] * mesh["dx"]**3).to("M_sun")
ind_cell = link_tracer(data)
mask_tracer = data["part"]["family"].values == 0
id_tracers = data["part"]["identity"].values[mask_tracer]
sort = np.argsort(id_tracers)
pos_tracers = data["part"]["position"].to("pc").values[mask_tracer][sort]
T_tracers = get_cell_data(data, ind_cell, "temperature", "K")[sort]
n_tracers = get_cell_data(data, ind_cell, "n", "cm**-3")[sort]
tracer_properties = dict(
id=id_tracers[sort],
position=pos_tracers,
temperature=T_tracers,
density=n_tracers
)
time = data.meta["time"].to("Myr").magnitude
tracer_properties["time"] = time
del data, mesh
return tracer_properties
def get_tracer_fast(path, snap, reextract=False, delete=True):
print("Extracting tracers from snapshot ", snap)
UTILS = os.path.expanduser("~/src/ramses/utils/f90")
out = f"{path}/tracers_{snap:05d}.dat"
if reextract or not os.path.exists(out):
command = f"{UTILS}/part2list -inp ~/simus/mass_inflow/n0.66_rms00000_tracers_memory_70/output_{snap:05d} -out {out}"
os.system(command)
f = FortranFile(out, 'r')
n = f.read_ints(np.int32)[0]
t = f.read_reals(np.float64)[0]
id_tracers = f.read_ints(np.int32)
x = f.read_reals(np.float64)
y = f.read_reals(np.float64)
z = f.read_reals(np.float64)
f.close()
sort = np.argsort(id_tracers)
#pos_tracers = np.array([x[sort], y[sort], z[sort]]).T
tracer_properties = dict(
id=id_tracers[sort],
z=z[sort],
time = t
)
if delete:
os.remove(out)
return tracer_properties
path = os.path.expanduser("~/simus/mass_inflow/n0.66_rms00000_tracers_memory_70")
snaps = np.arange(173, 185, 1)
extractor = "fast"
extractors = {"fast": get_tracer_fast,
"osyris": get_tracer_info_osyris}
tracers_all = []
for i in snaps:
tracers_all.append(extractors[extractor](path, i))
traj_tracers = {}
for key in tracers_all[0]:
if key not in ["id", "time"]:
traj_tracers[key] = np.array([tracers_all[i][key] for i in range(len(tracers_all))])
if traj_tracers[key].ndim == 2:
traj_tracers[key] = traj_tracers[key].transpose(1, 0)
elif traj_tracers[key].ndim == 3:
traj_tracers[key] = traj_tracers[key].transpose(1, 0, 2)
traj_tracers["id"] = tracers_all[0]["id"]
traj_tracers["time"] = [tracers_all[i]["time"] for i in range(len(tracers_all))]
name = os.path.basename(path)
np.savez(f"{path}/{name}_tracers.npz", **traj_tracers)

445
mass_inflow.py Normal file
View File

@@ -0,0 +1,445 @@
import osyris
import numpy as np
from numba import jit, prange
pc = osyris.units("pc")
year = osyris.units("year")
from scipy.spatial import KDTree
import matplotlib.pyplot as plt
import os
def recenter(array, center, periodicity=[0,1, 2]):
radial_vector = (array - center).to(pc)
radial_vector = np.array([radial_vector.x.values, radial_vector.y.values, radial_vector.z.values]).T
for i in periodicity:
mask = np.abs(radial_vector[:,i]) > 500
radial_vector[:,i][mask] = radial_vector[:,i][mask] - np.sign(radial_vector[:,i][mask]) * 1000
radial_vector = osyris.Vector(x=radial_vector[:,0], y=radial_vector[:,1], z=radial_vector[:,2], unit=pc)
return radial_vector
def find_tracer_index(data, tracer_id):
mask_tracer = data["part"]["family"].values == 0
tracer_index = np.nonzero((data["part"]["identity"].values == tracer_id) & mask_tracer)[0]
return int(tracer_index)
def find_closest_tracer(data, position):
mask_tracer = data["part"]["family"].values == 0
distance = (data["part"]["position"] - position).norm.to("pc").values[mask_tracer]
closest_tracer_index = np.argmin(distance)
closest_tracer_id = int(data["part"]["identity"][mask_tracer][closest_tracer_index].values)
return closest_tracer_id
def compute_derived_quantities(data, center=None, tracer=None, find_tracer=False, periodicity=[0, 1, 2]):
mesh = data["mesh"]
mesh["grav_potential"] = mesh["grav_potential"].to("pc^2 / Myr^2")
mesh["grav_force"] = mesh["density"] * mesh["grav_acceleration"] * mesh["dx"]**3
mesh["mass"] = (mesh["density"] * (mesh["dx"] ** 3)).to("M_sun")
ind = np.argmin(mesh["grav_potential"])
pos_minimum = mesh["position"][ind]
vel_minimum = mesh["velocity"][ind]
# Recenter the data
if isinstance(tracer, int):
tracer_index = find_tracer_index(data, tracer)
center = data["part"]["position"][tracer_index]
vel_ref = data["part"]["velocity"][tracer_index]
elif center is None and find_tracer:
tracer = find_closest_tracer(data, pos_minimum)
close_tracer_index = find_tracer_index(data, tracer)
center = data["part"]["position"][close_tracer_index]
vel_ref = data["part"]["velocity"][close_tracer_index]
elif center is None:
center = pos_minimum
vel_ref = vel_minimum
else:
tracer = find_closest_tracer(data, center)
close_tracer_index = find_tracer_index(data, tracer)
vel_ref = data["part"]["velocity"][close_tracer_index]
mesh["rel_velocity"] = mesh["velocity"] - vel_ref
sink = data["sink"]
sink["age"] = -(sink["tform"] - data.meta["time"])
part = data["part"]
mesh["radial_vector"] = recenter(mesh["position"], center, periodicity=periodicity)
sink["radial_vector"] = recenter(sink["position"], center, periodicity=periodicity)
part["radial_vector"] = recenter(part["position"], center, periodicity=periodicity)
# Save old position
if "old_position" not in mesh.keys():
mesh["old_position"] = mesh["position"]
sink["old_position"] = sink["position"]
part["old_position"] = part["position"]
mesh["position"].x._array = mesh["radial_vector"].x.to("cm").values
mesh["position"].y._array = mesh["radial_vector"].y.to("cm").values
mesh["position"].z._array = mesh["radial_vector"].z.to("cm").values
sink["position"].x._array = sink["radial_vector"].x.to("cm").values
sink["position"].y._array = sink["radial_vector"].y.to("cm").values
sink["position"].z._array = sink["radial_vector"].z.to("cm").values
part["position"].x._array = part["radial_vector"].x.to("cm").values
part["position"].y._array = part["radial_vector"].y.to("cm").values
part["position"].z._array = part["radial_vector"].z.to("cm").values
# Compute integrated accelerations
data["part"]["velocity_grav"] = data["part"]["vp_grav"] * data.units["velocity"]
data["part"]["velocity_init"] = data["part"]["vp_init"] * data.units["velocity"]
data["part"]["velocity_diff"] = data["part"]["velocity"] - data["part"]["velocity_init"]
# Change potential gauge
min_pot = np.min(data["mesh"]["grav_potential"].values)
data["mesh"]["grav_potential"].values = data["mesh"]["grav_potential"].values - min_pot
data["mesh"]["grav_potential"] = data["mesh"]["grav_potential"].to("km*m/s^2")
# Compute temperature
mesh["temperature"] = ((mesh["pressure"] / mesh["density"] ) * (1.4 *osyris.units("proton_mass") / osyris.units("boltzmann_constant"))).to("K")
# Compute distance to sink
pos_sinks = data["sink"]["position"].to("pc").values
tree = KDTree(pos_sinks)
dist_sinks, _ = tree.query(data["part"]["position"].to("pc").values )
data["part"]["dist_sinks"] = osyris.Array(dist_sinks, unit=pc)
if isinstance(tracer, int):
return tracer
else:
return center
def load_snapshot(dataset, ind_cell_d, path, snap):
print(f"Load dataset {snap}")
dataset[snap] = osyris.RamsesDataset(snap, path=os.path.expanduser(path)).load(select=["mesh","hydro", "part", "sink"])
print(f"Compute derived quantities {snap}")
compute_derived_quantities(dataset[snap], periodicity=[0,1])
print(f"Build tracer index {snap}")
ind_cell_d[snap] = link_tracer(dataset, snap)
## Link tracers to cells
def link_tracer(dataset, snap):
mask_tracer = dataset[snap]["part"]["family"].values == 0
pos_cells = dataset[snap]["mesh"]["position"].to("pc").values
pos_tracers = dataset[snap]["part"]["position"].to("pc").values[mask_tracer]
kdtree = KDTree(pos_cells)
dist, ind_cell = kdtree.query(pos_tracers, k=1, p=np.inf)
return ind_cell
def build_tracer_index(dataset, snaps):
ind_cell_d = {}
for i in snaps:
print(f"Linking tracers to cells for snapshot {i}...")
ind_cell_d[i] = link_tracer(dataset, i)
return ind_cell_d
def get_cell_data(dataset, ind_cell_d, snap, key, unit):
ind_cell = ind_cell_d[snap]
cell_data = dataset[snap]["mesh"][key].to(unit).values[ind_cell]
return cell_data
@jit(nopython=True, parallel=True)
def compute_fraction_gravdom_helper(mass, dist, aligned, bins, cumulative:bool, masks:np.ndarray, stacked=False):
fraction = np.zeros(shape=(max(1, masks.shape[0]), bins.shape[0]), dtype=np.float64)
for i in prange(len(bins)):
if i > 0 and not cumulative:
mask_distance = (bins[i - 1] < dist) & (dist <= bins[i])
else:
mask_distance = (dist <= bins[i])
if np.sum(mask_distance) == 0:
fraction[:,i] = np.nan
else:
for m in range(masks.shape[0]):
if stacked:
mask_total = mask_distance
else:
mask_total = mask_distance & masks[m]
fraction[m][i] = np.nansum(mass[mask_distance & aligned & masks[m]]) / np.nansum(mass[mask_total])
if masks.shape[0] == 0:
fraction[0][i] = np.nansum(mass[mask_distance & aligned]) / np.nansum(mass[mask_distance])
return fraction
def compute_fraction_tracer(dataset, i, bins, id_grav_dom, cumulative=True, masks={}, mask_id=False, stacked=False):
bins = bins
mask_tracer = dataset[i]["part"]["family"].values == 0
mask_tracer_gravdom = np.isin(dataset[i]["part"]["identity"].values[mask_tracer], id_grav_dom, assume_unique=True)
mass = dataset[i]["part"]["mass"].values[mask_tracer]
dist = dataset[i]["part"]["radial_vector"].norm.to(pc).values[mask_tracer]
if mask_id:
masks = {key: np.isin(dataset[i]["part"]["identity"].values[mask_tracer], value, assume_unique=True) for key, value in masks.items()}
fraction = compute_fraction_gravdom_helper(mass, dist, mask_tracer_gravdom, bins, cumulative, np.array(list(masks.values()), dtype=np.bool_), stacked=stacked)
if len(masks) > 0:
fraction = {key: value for key, value in zip(masks.keys(), fraction)}
else:
fraction = fraction[0]
return bins, fraction
def get_aligment(dataset, ind_cell_d, start, end):
mask_tracer_end = dataset[end]["part"]["family"].values == 0
distance_end = dataset[end]["part"]["position"].norm.to(pc).values[mask_tracer_end]
id_end = dataset[end]["part"]["identity"].values[mask_tracer_end]
sort_end = np.argsort(id_end)
unsort_end = np.argsort(sort_end)
velocity_end = dataset[end]["part"]["velocity"].to("km/s").values[mask_tracer_end]
velocity_grav_end = dataset[end]["part"]["velocity_grav"].to("km/s").values[mask_tracer_end]
density_end = get_cell_data(dataset, ind_cell_d, end, "density", "g/cm**3")
mask_tracer_start = dataset[start]["part"]["family"].values == 0
distance_start = dataset[start]["part"]["position"].norm.to(pc).values[mask_tracer_start]
id_start = dataset[start]["part"]["identity"].values[mask_tracer_start]
sort_start = np.argsort(id_start)
velocity_start = dataset[start]["part"]["velocity"].to("km/s").values[mask_tracer_start]
velocity_grav_start = dataset[start]["part"]["velocity_grav"].to("km/s").values[mask_tracer_start]
density_start = get_cell_data(dataset, ind_cell_d, start, "density", "g/cm**3")
distance_diff_sorted = distance_end[sort_end] - distance_start[sort_start]
velocity_diff_sorted = velocity_end[sort_end] - velocity_start[sort_start]
velocity_grav_diff_sorted = velocity_grav_end[sort_end] - velocity_grav_start[sort_start]
density_diff_sorted = density_end[sort_end] - density_start[sort_start]
distance_diff_end = distance_diff_sorted[unsort_end]
velocity_diff_end = velocity_diff_sorted[unsort_end]
velocity_grav_diff_end = velocity_grav_diff_sorted[unsort_end]
density_diff_end = density_diff_sorted[unsort_end]
sorters = dict(
sort_end=sort_end,
unsort_end=unsort_end,
sort_start=sort_start,
)
alignement = dict(
start = start,
end = end,
velocity_diff_end=velocity_diff_end, # km/s
velocity_grav_diff_end=velocity_grav_diff_end, # km/s
distance_diff_end=distance_diff_end, # pc
density_diff_end=density_diff_end, # g/cm3
sorters=sorters,)
delta_t = dataset[end].meta["time"].to("Myr").magnitude - dataset[start].meta["time"].to("Myr").magnitude
alignement['delta_t'] = delta_t # Myr
acc_grav = alignement['velocity_grav_diff_end'] / delta_t
alignement['acc_grav'] = osyris.Vector(x=acc_grav[:,0], y=acc_grav[:,1], z=acc_grav[:,2], unit="km/s/ Myr")
acc_tot = alignement['velocity_diff_end'] / delta_t
alignement['acc_tot'] = osyris.Vector(x=acc_tot[:,0], y=acc_tot[:,1], z=acc_tot[:,2], unit="km/s/ Myr")
acc_other = acc_tot - acc_grav
alignement['acc_other'] = osyris.Vector(x=acc_other[:,0], y=acc_other[:,1], z=acc_other[:,2], unit="km/s/ Myr")
norm_acc_grav = np.linalg.norm(acc_grav, axis=1)
norm_acc_other = np.linalg.norm(acc_other, axis=1)
mask_grav_dom = (norm_acc_grav > norm_acc_other)
id_end = dataset[end]["part"]["identity"].values[mask_tracer_end]
mask_denser = (alignement['density_diff_end'] > 0)
mask_grav_dom_denser = mask_grav_dom & mask_denser
alignement["mask_gravdom"] = mask_grav_dom
alignement['id_grav_dom_denser'] = id_end[mask_grav_dom_denser]
alignement["id_grav_dom"] = id_end[mask_grav_dom]
return alignement
def compute_mass_inflow_distance(dataset, alignement_d, end, masks=None, key_mask_gravdom="id_grav_dom"):
mask_tracer_end = dataset[end]["part"]["family"].values == 0
distance_end = dataset[end]["part"]["position"].norm.to(pc).values[mask_tracer_end]
distance_diff_end = alignement_d[end]['distance_diff_end']
distance_start = distance_end - distance_diff_end
mass_total = dataset[end]["part"]["mass"].to("M_sun").values[mask_tracer_end]
delta_t = alignement_d[end]['delta_t']
dist = np.logspace(np.log10(np.min(distance_end)), np.log10(np.max(distance_end)), 55)
flow_total = np.zeros_like(dist)
flow_gravdom = np.zeros_like(dist)
enclosed_mass_start = np.zeros_like(dist)
enclosed_mass_end = np.zeros_like(dist)
id_tracer_end = dataset[end]["part"]["identity"].values[mask_tracer_end]
if isinstance(key_mask_gravdom, str):
mask_tracer_gravdom = np.isin(id_tracer_end, alignement_d[end][key_mask_gravdom], assume_unique=True)
elif isinstance(key_mask_gravdom, np.ndarray) and isinstance(key_mask_gravdom[0], bool):
mask_tracer_gravdom = key_mask_gravdom
else:
mask_tracer_gravdom = np.isin(id_tracer_end, key_mask_gravdom, assume_unique=True)
print(delta_t)
for i, d in enumerate(dist):
in_flow = (distance_end > d) & (distance_start < d)
flow_total[i] = np.sum(mass_total[in_flow])/ delta_t
flow_gravdom[i] = np.sum(mass_total[in_flow & mask_tracer_gravdom])/ delta_t
enclosed_mass_start[i] = np.sum(mass_total[distance_start > d])
enclosed_mass_end[i] = np.sum(mass_total[distance_end > d])
flow_mask = {}
enclosed_mass_end_mask = {}
enclosed_mass_start_mask = {}
for key in masks:
flow_mask[key] = np.zeros_like(dist)
flow_mask[key + "_gravdom"] = np.zeros_like(dist)
enclosed_mass_end_mask[key] = np.zeros_like(dist)
enclosed_mass_start_mask[key] = np.zeros_like(dist)
for i, d in enumerate(dist):
in_flow = (distance_end > d) & (distance_start < d)
flow_mask[key][i] = np.sum(mass_total[in_flow & masks[key]])/ delta_t
flow_mask[key + "_gravdom"][i] = np.sum(mass_total[in_flow & masks[key] & mask_tracer_gravdom]) / delta_t
enclosed_mass_end_mask[key][i] = np.sum(mass_total[(distance_end > d) & masks[key]])
enclosed_mass_start_mask[key][i] = np.sum(mass_total[(distance_start> d) & masks[key]])
result = dict(distance = dist,
flow_total = flow_total,
flow_gravdom = flow_gravdom,
enclosed_mass_start = enclosed_mass_start,
enclosed_mass_end = enclosed_mass_end,
flow_mask = flow_mask,
enclosed_mass_start_mask = enclosed_mass_start_mask,
enclosed_mass_end_mask = enclosed_mass_end_mask)
return result
def compute_mass_inflow_density(dataset, ind_cell_d, alignement_gravdom, end, masks=None, key_mask_gravdom="id_grav_dom", alignement_inflow=None):
mask_tracer_end = dataset[end]["part"]["family"].values == 0
if alignement_inflow is None:
alignement_inflow = alignement_gravdom
density_end = get_cell_data(dataset, ind_cell_d, end, "density", "g/cm**3")
density_diff_end = alignement_inflow[end]['density_diff_end']
density_start = density_end - density_diff_end
mass_total = dataset[end]["part"]["mass"].to("M_sun").values[mask_tracer_end]
delta_t = alignement_gravdom[end]['delta_t']
dens = np.logspace(np.log10(np.min(density_end)), np.log10(np.max(density_end)), 55)
flow_total = np.zeros_like(dens)
flow_gravdom = np.zeros_like(dens)
enclosed_mass_start = np.zeros_like(dens)
enclosed_mass_end = np.zeros_like(dens)
id_tracer_end = dataset[end]["part"]["identity"].values[mask_tracer_end]
if isinstance(key_mask_gravdom, str):
mask_tracer_gravdom = np.isin(id_tracer_end, alignement_gravdom[end][key_mask_gravdom], assume_unique=True)
elif isinstance(key_mask_gravdom, np.ndarray) and isinstance(key_mask_gravdom[0], bool):
mask_tracer_gravdom = key_mask_gravdom
else:
mask_tracer_gravdom = np.isin(id_tracer_end, key_mask_gravdom, assume_unique=True)
print(delta_t)
for i, d in enumerate(dens):
in_flow = (density_end > d) & (density_start < d)
flow_total[i] = np.sum(mass_total[in_flow])/ delta_t
flow_gravdom[i] = np.sum(mass_total[in_flow & mask_tracer_gravdom])/ delta_t
enclosed_mass_start[i] = np.sum(mass_total[density_start > d])
enclosed_mass_end[i] = np.sum(mass_total[density_end > d])
flow_mask = {}
enclosed_mass_end_mask = {}
enclosed_mass_start_mask = {}
for key in masks:
flow_mask[key] = np.zeros_like(dens)
flow_mask[key + "_gravdom"] = np.zeros_like(dens)
enclosed_mass_end_mask[key] = np.zeros_like(dens)
enclosed_mass_start_mask[key] = np.zeros_like(dens)
for i, d in enumerate(dens):
in_flow = (density_end > d) & (density_start < d)
flow_mask[key][i] = np.sum(mass_total[in_flow & masks[key]])/ delta_t
flow_mask[key + "_gravdom"][i] = np.sum(mass_total[in_flow & masks[key] & mask_tracer_gravdom]) / delta_t
enclosed_mass_end_mask[key][i] = np.sum(mass_total[(density_end > d) & masks[key]])
enclosed_mass_start_mask[key][i] = np.sum(mass_total[(density_start> d) & masks[key]])
result = dict(density = dens,
flow_total = flow_total,
flow_gravdom = flow_gravdom,
enclosed_mass_start = enclosed_mass_start,
enclosed_mass_end = enclosed_mass_end,
flow_mask = flow_mask,
enclosed_mass_start_mask = enclosed_mass_start_mask,
enclosed_mass_end_mask = enclosed_mass_end_mask)
return result
## Plot helpers
def linehist(x, weights, bins=100, density=False, ax=None, **kwargs):
"""
Draw a histogram as a line plot.
"""
if ax is None:
ax = plt.gca()
hist, bin_edges = np.histogram(x, bins=bins, weights=weights, density=density)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
ax.plot(bin_centers, hist, **kwargs)
return bin_centers, hist
def scatter_hist(x, y, weights, ax, ax_histx, ax_histy):
"""
Scatter plot with marginal histograms.
"""
# no labels
ax_histx.tick_params(axis="x", labelbottom=False)
ax_histy.tick_params(axis="y", labelleft=False)
# the hist2d plot:
im = ax.hist2d(x, y, bins=100, weights=weights, norm=mpl.colors.LogNorm())
# now determine nice limits by hand:
ax_histx.hist(x, bins=100, weights=weights, density=True, histtype="step", lw=2)
ax_histy.hist(y, bins=100, orientation='horizontal', weights=weights, histtype="step", density=True, lw=2)
ax_histx.set_yscale("log")
ax_histy.set_xscale("log")
return im
# Data helpers
def groupby(array, key):
"""
Group an array by a key.
"""
sort = np.argsort(key)
key = key[sort]
array = array[sort]
unique, indices = np.unique(key, return_index=True)
return unique, np.split(array, indices[1:])
def aggregate(array, key, function=np.sum):
"""
Aggregate an array by a key using a function.
"""
unique, groups = groupby(array, key)
return unique, np.array([function(x) for x in groups])

489
mass_inflow_time_evo.ipynb Normal file
View File

@@ -0,0 +1,489 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib as mpl\n",
"import hickle as hkl\n",
"import osyris"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"snaps = [174, 175, 177, 178, 179, 180, 181, 182, 184, 185]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"allresults = {}\n",
"for snap in snaps:\n",
" allresults[snap] = hkl.load(f\"/home/nbrucyci/src/mass-inflow/results_sinks_173_173_{snap}.h5\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"time = np.array([allresults[snap]['time'] for snap in snaps])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"from matplotlib import ticker as tick\n",
"import matplotlib.colors as colors\n",
"\n",
"def get_colors(ax, labels, cmap, clabel=\"\", put_colorbar=True,\n",
" cbar_fmt=tick.FormatStrFormatter('%.3g')):\n",
" \n",
" colors = None\n",
" if isinstance(cmap, str):\n",
" cmap = plt.get_cmap(cmap)\n",
" bounds = np.concatenate([labels[0:1] - (labels[1] - labels[0]) / 2, (labels[:-1] + labels[1:]) / 2, labels[-1:] + (labels[-1] - labels[-2])/2 ])\n",
" norm = mpl.colors.BoundaryNorm(bounds, cmap.N) \n",
" sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)\n",
" if put_colorbar:\n",
" cb = plt.colorbar(sm, ticks=labels, ax=ax)\n",
" cb.set_label(clabel)\n",
" cb.ax.minorticks_off()\n",
" cb.ax.yaxis.set_major_formatter(cbar_fmt)\n",
" colors = lambda xi: sm.cmap(sm.norm(xi)) \n",
" return colors\n",
"\n",
"\n",
"def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):\n",
" new_cmap = colors.LinearSegmentedColormap.from_list(\n",
" 'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),\n",
" cmap(np.linspace(minval, maxval, n)))\n",
" return new_cmap"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Hcc = (1.4 * osyris.units(\"m_p/cm**3\")).to(\"g/cm**3\")\n",
"sink_threshold = 1e3 * Hcc\n",
"CNM = 1e-23\n",
"H2 = 2e-22\n",
"\n",
"def density_to_n(density):\n",
" return (density /Hcc).magnitude\n",
"\n",
"fig, axs = plt.subplots(2, 1, figsize=(6,9), constrained_layout=True)\n",
"\n",
"greys2 = truncate_colormap(plt.get_cmap(\"Greys\"), 0.2, 1)\n",
"\n",
"ax = axs[0]\n",
"\n",
"c = get_colors(ax, time, greys2, clabel=\"$t_e$\", put_colorbar=True, cbar_fmt=tick.FormatStrFormatter('%.3g'))\n",
"\n",
"ax_top = ax.twiny()\n",
"ax_top.set_xlabel(r\"$n$ [cm$^{-3}$]\")\n",
"ax_top.set_xscale(\"log\")\n",
"\n",
"def convert_lim(ax):\n",
" x1, x2 = ax.get_xlim()\n",
" ax_top.set_xlim(density_to_n(x1), density_to_n(x2))\n",
" ax_top.figure.canvas.draw()\n",
"\n",
"ax.callbacks.connect(\"xlim_changed\", convert_lim)\n",
"\n",
"ax.vlines(sink_threshold.magnitude, 0, 60,color=\"darkorange\", lw=3, label=\"Sink threshold\")\n",
"ax.axvspan(CNM, H2, 0, 60, color=\"tab:blue\", alpha=0.2, label='\"CNM\"')\n",
"ax.axvspan(H2, 1.2*sink_threshold.magnitude, 0, 60, color=\"tab:blue\", alpha=0.5, label='\"Molecular gas\"')\n",
"\n",
"\n",
"ax.legend(loc=\"upper left\")\n",
"ax.set_title(\"Within 100 pc of the minimum of potential\")\n",
"\n",
"\n",
"ax.set_ylabel(r\"$\\dot{\\mathrm{M}}_\\mathrm{H}$($t_s$, $t_e$) [M$_\\odot$/kyr]\")\n",
"\n",
"ax.set_ylim(0, 30)\n",
"ax.legend(loc=\"upper left\", frameon=False, ncol=1)\n",
"\n",
"ax.set_xscale(\"log\")\n",
"ax.set_xlim(1e-25, 1.2*sink_threshold.magnitude)\n",
"ax.set_xlabel(r\"$\\rho$ [g cm$^{-3}$]\")\n",
"\n",
"\n",
"for i, snap in enumerate(snaps):\n",
" mass_inflow = allresults[snap][\"mass_inflow\"]\n",
" \n",
" t = time[i] \n",
"\n",
" ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist\"]/1e3, label=\"Total\", alpha=0.7, lw=4, ls='-', color=c(t))\n",
" ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/1e3, label=\"Gravity-driven\", alpha=0.7, lw=4, ls='--', color=c(t))\n",
" #ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/mass_inflow[\"flow_mask\"][\"dist\"], \n",
" # label=\"Gravity-driven / Total\", lw=4, ls='-', color=c(t))\n",
"\n",
"\n",
"ax = axs[1]\n",
"\n",
"c = get_colors(ax, time, \"Reds\", clabel=\"$t_e$\", put_colorbar=True, cbar_fmt=tick.FormatStrFormatter('%.3g'))\n",
"\n",
"ax_top = ax.twiny()\n",
"ax_top.set_xlabel(r\"$n$ [cm$^{-3}$]\")\n",
"ax_top.set_xscale(\"log\")\n",
"\n",
"def convert_lim(ax):\n",
" x1, x2 = ax.get_xlim()\n",
" ax_top.set_xlim(density_to_n(x1), density_to_n(x2))\n",
" ax_top.figure.canvas.draw()\n",
"\n",
"ax.callbacks.connect(\"xlim_changed\", convert_lim)\n",
"\n",
"ax.vlines(sink_threshold.magnitude, 0, 60,color=\"darkorange\", lw=3, label=\"Sink threshold\")\n",
"ax.axvspan(CNM, H2, 0, 60, color=\"tab:blue\", alpha=0.2, label='\"CNM\"')\n",
"ax.axvspan(H2, 1.2*sink_threshold.magnitude, 0, 60, color=\"tab:blue\", alpha=0.5, label='\"Molecular gas\"')\n",
"\n",
"\n",
"ax.legend(loc=\"upper left\")\n",
"\n",
"\n",
"ax.set_ylabel(r\"Fraction of gravity-driven inflow\")\n",
"\n",
"ax.set_ylim(0, 0.58)\n",
"ax.legend(loc=\"upper left\", frameon=False, ncol=1)\n",
"\n",
"ax.set_xscale(\"log\")\n",
"ax.set_xlim(1e-25, 1.2*sink_threshold.magnitude)\n",
"ax.set_xlabel(r\"$\\rho$ [g cm$^{-3}$]\")\n",
"\n",
"for i, snap in enumerate(snaps):\n",
" mass_inflow = allresults[snap][\"mass_inflow\"]\n",
" \n",
" t = time[i] \n",
"\n",
" #ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist\"]/1e3, label=\"Total\", alpha=0.7, lw=4, ls='-', color=\"k\")\n",
" #ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/1e3, label=\"Gravity-driven\", alpha=0.7, lw=4, ls='--', color=\"k\")\n",
" ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/mass_inflow[\"flow_mask\"][\"dist\"], \n",
" label=\"Gravity-driven / Total\", lw=4, ls='-', color=c(t))\n",
"\n",
"\n",
"plt.savefig(f\"mass_inflow_density_all.pdf\", bbox_inches=\"tight\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Hcc = (1.4 * osyris.units(\"m_p/cm**3\")).to(\"g/cm**3\")\n",
"sink_threshold = 1e3 * Hcc\n",
"CNM = 1e-23\n",
"H2 = 2e-22\n",
"\n",
"def density_to_n(density):\n",
" return (density /Hcc).magnitude\n",
"\n",
"fig, ax = plt.subplots(1, 1, figsize=(6,4), constrained_layout=True)\n",
"\n",
"#c2 = get_colors(ax, time, \"Purples\", clabel=\"$t_e$ (without sinks' vicinity)\", put_colorbar=True, cbar_fmt=tick.FormatStrFormatter('%.3g'))\n",
"\n",
"c = get_colors(ax, time, \"Reds\", clabel=\"$t_e$\", put_colorbar=True, cbar_fmt=tick.FormatStrFormatter('%.3g'))\n",
"\n",
"\n",
"ax_top = ax.twiny()\n",
"ax_top.set_xlabel(r\"$n$ [cm$^{-3}$]\")\n",
"ax_top.set_xscale(\"log\")\n",
"\n",
"def convert_lim(ax):\n",
" x1, x2 = ax.get_xlim()\n",
" ax_top.set_xlim(density_to_n(x1), density_to_n(x2))\n",
" ax_top.figure.canvas.draw()\n",
"\n",
"ax.callbacks.connect(\"xlim_changed\", convert_lim)\n",
"\n",
"ax.vlines(sink_threshold.magnitude, 0, 60,color=\"darkorange\", lw=3, label=\"Sink threshold\")\n",
"ax.axvspan(CNM, H2, 0, 60, color=\"tab:blue\", alpha=0.2, label='\"CNM\"')\n",
"ax.axvspan(H2, 1.2*sink_threshold.magnitude, 0, 60, color=\"tab:blue\", alpha=0.5, label='\"Molecular gas\"')\n",
"\n",
"\n",
"ax.legend(loc=\"upper left\")\n",
"\n",
"\n",
"ax.set_ylabel(r\"Fraction of gravity-driven inflow\")\n",
"\n",
"ax.set_ylim(0, 0.58)\n",
"\n",
"ax.set_xscale(\"log\")\n",
"\n",
"ax.set_xlim(1e-25, 1.2*sink_threshold.magnitude)\n",
"ax.set_xlabel(r\"$\\rho$ [g cm$^{-3}$]\")\n",
"\n",
"for i, snap in enumerate(snaps):\n",
" mass_inflow = allresults[snap][\"mass_inflow\"]\n",
" \n",
" t = time[i] \n",
" \n",
" frac_with_sinks = mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/mass_inflow[\"flow_mask\"][\"dist\"]\n",
" frac_without_sinks = mass_inflow[\"flow_mask\"][\"dist_sinks_gravdom\"]/mass_inflow[\"flow_mask\"][\"dist_sinks\"]\n",
" \n",
" ax.plot( mass_inflow[\"density\"], frac_with_sinks, lw=2, ls='-', color=c(t))\n",
"\n",
" ax.plot( mass_inflow[\"density\"], frac_without_sinks, lw=2, ls='-.', color=c(t))\n",
" \n",
"ax.plot([], [], lw=2, ls='-', color='k', label=\"with sinks' vicinity\")\n",
"ax.plot([], [], lw=2, ls='-.', color='k', label=\"without sinks' vicinity\")\n",
"ax.legend(loc=\"upper left\", frameon=False, ncol=1)\n",
"\n",
"\n",
"plt.savefig(f\"mass_inflow_density_sink_time_evo.pdf\", bbox_inches=\"tight\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(6,5), constrained_layout=True)\n",
"ax = plt.gca()\n",
"\n",
"ax_top = ax.twiny()\n",
"ax_top.set_xlabel(r\"$n$ [cm$^{-3}$]\")\n",
"ax_top.set_xscale(\"log\")\n",
"\n",
"def convert_lim(ax):\n",
" x1, x2 = ax.get_xlim()\n",
" ax_top.set_xlim(density_to_n(x1), density_to_n(x2))\n",
" ax_top.figure.canvas.draw()\n",
"\n",
"ax.callbacks.connect(\"xlim_changed\", convert_lim)\n",
"c = get_colors(ax, time, greys2, clabel=\"t\", put_colorbar=True, cbar_fmt=tick.FormatStrFormatter('%.3g'))\n",
"\n",
"\n",
"\n",
"ax.vlines(sink_threshold.magnitude, 0, 1e6, color=\"tab:orange\", lw=2, ls=\":\", label=\"Sink threshold\")\n",
"ax.axvspan(CNM, H2, 0, 60, color=\"tab:blue\", alpha=0.2, label=\"CNM\")\n",
"ax.axvspan(H2, 1.2*sink_threshold.magnitude, 0, 60, color=\"tab:blue\", alpha=0.5, label=\"Molecular gas\")\n",
"\n",
"ax.set_xlim(5e-25, 1.2* sink_threshold.magnitude)\n",
"ax.set_ylim(1e2, 1e6)\n",
"\n",
"for i, snap in enumerate(snaps):\n",
" mass_inflow = allresults[snap][\"mass_inflow\"]\n",
"\n",
" ax.loglog(mass_inflow[\"density\"], mass_inflow[\"enclosed_mass_end_mask\"][\"dist\"], lw=4, color=c(time[i]), alpha=0.6)\n",
"\n",
"\n",
"\n",
"ax.set_xlabel(r\"$\\rho$ [g/cm$^3$]\")\n",
"ax.set_ylabel(r\"M($\\rho^\\prime > \\rho$ ) [M$_\\odot$]\")\n",
"ax.set_title(\"Tracer mass within 100 pc of the minimum of the potential\")\n",
"ax.legend()\n",
"plt.savefig(f\"mass_enclosed_density_evo.pdf\", bbox_inches=\"tight\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"plt.figure(figsize=(5,4), constrained_layout=True)\n",
"ax = plt.gca()\n",
"c = get_colors(ax, time, \"Reds\", clabel=\"time [Myr]\", put_colorbar=True, cbar_fmt=tick.FormatStrFormatter('%.3g'))\n",
"\n",
"for i, snap in enumerate(snaps):\n",
" dist, frac_end = allresults[snap][\"mass_fraction\"][\"dist\"], allresults[snap][\"mass_fraction\"][\"frac_end\"]\n",
" t = time[i]\n",
" plt.plot(dist, frac_end, alpha=0.7, lw=4, ls='-', color=c(t), label=\"all\")\n",
"\n",
"plt.xlabel(\"Distance from potential minimum (pc)\")\n",
"plt.ylabel(\"Mass fraction of gravity-driven gas\")\n",
"plt.xscale(\"log\")\n",
"plt.savefig(\"mass_fraction_stratified_tracers.pdf\", bbox_inches=\"tight\")\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Normalized time"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"allresults_normalized = {}\n",
"for snap in snaps[:-2]:\n",
" allresults_normalized[snap] = hkl.load(f\"/home/nbrucyci/src/mass-inflow/results__{snap}_184_185.h5\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Hcc = (1.4 * osyris.units(\"m_p/cm**3\")).to(\"g/cm**3\")\n",
"sink_threshold = 1e3 * Hcc\n",
"CNM = 1e-23\n",
"H2 = 2e-22\n",
"\n",
"def density_to_n(density):\n",
" return (density /Hcc).magnitude\n",
"\n",
"plt.figure(figsize=(6,5), constrained_layout=True)\n",
"ax = plt.gca()\n",
"\n",
"\n",
"c = get_colors(ax, time, \"Reds\", clabel=\"t\", put_colorbar=True, cbar_fmt=tick.FormatStrFormatter('%.3g'))\n",
"\n",
"ax_top = ax.twiny()\n",
"ax_top.set_xlabel(r\"$n$ [cm$^{-3}$]\")\n",
"ax_top.set_xscale(\"log\")\n",
"\n",
"def convert_lim(ax):\n",
" x1, x2 = ax.get_xlim()\n",
" ax_top.set_xlim(density_to_n(x1), density_to_n(x2))\n",
" ax_top.figure.canvas.draw()\n",
"\n",
"ax.callbacks.connect(\"xlim_changed\", convert_lim)\n",
"\n",
"ax.vlines(sink_threshold.magnitude, 0, 60,color=\"darkorange\", lw=3, label=\"Sink threshold\")\n",
"ax.axvspan(CNM, H2, 0, 60, color=\"tab:blue\", alpha=0.2, label='\"CNM\"')\n",
"ax.axvspan(H2, 1.2*sink_threshold.magnitude, 0, 60, color=\"tab:blue\", alpha=0.5, label=\"Molecular gas\")\n",
"\n",
"\n",
"ax.legend(loc=\"upper left\")\n",
"ax.set_title(\"Within 100 pc of the minimum of potential\")\n",
"\n",
"\n",
"ax.set_ylabel(r\"Fraction of gravity-driven inflow\")\n",
"\n",
"ax.set_ylim(0, 0.45)\n",
"ax.legend(loc=\"upper left\", frameon=False, ncol=1)\n",
"\n",
"ax.set_xscale(\"log\")\n",
"ax.set_xlim(1e-25, 1.2*sink_threshold.magnitude)\n",
"\n",
"\n",
"for i, snap in enumerate(snaps):\n",
" mass_inflow = allresults_normalized[snap][\"mass_inflow\"]\n",
" \n",
" t = time[i] \n",
"\n",
" #ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist\"]/1e3, label=\"Total\", alpha=0.7, lw=4, ls='-', color=\"k\")\n",
" #ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/1e3, label=\"Gravity-driven\", alpha=0.7, lw=4, ls='--', color=\"k\")\n",
" ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/mass_inflow[\"flow_mask\"][\"dist\"], \n",
" label=\"Gravity-driven / Total\", lw=4, ls='-', color=c(t))\n",
"\n",
"\n",
"plt.savefig(f\"mass_inflow_density_all_lastmyr.pdf\", bbox_inches=\"tight\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Hcc = (1.4 * osyris.units(\"m_p/cm**3\")).to(\"g/cm**3\")\n",
"sink_threshold = 1e3 * Hcc\n",
"CNM = 1e-23\n",
"H2 = 2e-22\n",
"\n",
"def density_to_n(density):\n",
" return (density /Hcc).magnitude\n",
"\n",
"plt.figure(figsize=(6,5), constrained_layout=True)\n",
"ax = plt.gca()\n",
"\n",
"\n",
"c = get_colors(ax, time, \"Greys\", clabel=\"t\", put_colorbar=True, cbar_fmt=tick.FormatStrFormatter('%.3g'))\n",
"\n",
"ax_top = ax.twiny()\n",
"ax_top.set_xlabel(r\"$n$ [cm$^{-3}$]\")\n",
"ax_top.set_xscale(\"log\")\n",
"\n",
"def convert_lim(ax):\n",
" x1, x2 = ax.get_xlim()\n",
" ax_top.set_xlim(density_to_n(x1), density_to_n(x2))\n",
" ax_top.figure.canvas.draw()\n",
"\n",
"ax.callbacks.connect(\"xlim_changed\", convert_lim)\n",
"\n",
"ax.vlines(sink_threshold.magnitude, 0, 60,color=\"darkorange\", lw=3, label=\"Sink threshold\")\n",
"ax.axvspan(CNM, H2, 0, 60, color=\"tab:blue\", alpha=0.2, label='\"CNM\"')\n",
"ax.axvspan(H2, 1.2*sink_threshold.magnitude, 0, 60, color=\"tab:blue\", alpha=0.5, label=\"Molecular gas\")\n",
"\n",
"\n",
"ax.legend(loc=\"upper left\")\n",
"ax.set_title(\"Within 100 pc of the minimum of potential\")\n",
"\n",
"\n",
"ax.set_ylabel(r\"Mass inflow rate [M$_\\odot$/kyr]\")\n",
"\n",
"ax.set_ylim(0, 100)\n",
"ax.legend(loc=\"upper left\", frameon=False, ncol=1)\n",
"\n",
"ax.set_xscale(\"log\")\n",
"ax.set_xlim(1e-25, 1.2*sink_threshold.magnitude)\n",
"\n",
"\n",
"for i, snap in enumerate(snaps[:-2]):\n",
" mass_inflow = allresults_normalized[snap][\"mass_inflow\"]\n",
" \n",
" t = time[i] \n",
"\n",
" ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist\"]/1e3, label=\"Total\", alpha=0.7, lw=4, ls='-', color=c(t))\n",
" ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/1e3, label=\"Gravity-driven\", alpha=0.7, lw=4, ls='--', color=c(t))\n",
" #ax.plot(mass_inflow[\"density\"], mass_inflow[\"flow_mask\"][\"dist_gravdom\"]/mass_inflow[\"flow_mask\"][\"dist\"], \n",
" # label=\"Gravity-driven / Total\", lw=4, ls='-', color=c(t))\n",
"\n",
"\n",
"plt.savefig(f\"mass_inflow_density_all_flow_lastmyr.pdf\", bbox_inches=\"tight\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

1841
mass_inflow_tracers_classical.ipynb Normal file
View File

File diff suppressed because it is too large Load Diff

344
run_mass_inflow_time.py Normal file
View File

@@ -0,0 +1,344 @@
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import osyris
import numpy as np
from numba import jit, prange
from pysankey import sankey
import hickle as hkl
pc = osyris.units("pc")
year = osyris.units("year")
center = osyris.Vector(0.5, 0.5, 0.5, unit=pc) * 1000
import mass_inflow as mi
## Mass fraction
def mass_fraction(iend, alignement):
print(f"Mass fraction for {iend}")
dist, frac_end = mi.compute_fraction_tracer(dataset, iend, np.logspace(np.log10(5), 3, 30), alignement[iend]['id_grav_dom'], cumulative=False)
plt.figure(figsize=(5,4), constrained_layout=True)
plt.plot(dist, frac_end, alpha=0.7, lw=4, ls='-', color="k", label="all")
plt.xlabel("Distance from potential minimum (pc)")
plt.ylabel("Mass fraction of gravity-driven gas")
plt.xscale("log")
plt.savefig(f"mass_fraction_stratified_tracers_{iend}.pdf", bbox_inches="tight")
results[iend]["mass_fraction"] = {"dist": dist, "frac_end": frac_end}
def comp_acc(iend, alignement):
print(f"Compare accelerations for {iend}")
mask_tracer = dataset[iend]["part"]["family"].values == 0
id_grav_dom = alignement[iend]['id_grav_dom']
mask_tracer_gravdom = np.isin(dataset[iend]["part"]["identity"].values[mask_tracer], id_grav_dom, assume_unique=True)
mass_all = dataset[iend]["part"][mask_tracer]["mass"].to("M_sun").values
atot_tracers_all = alignement[iend]["acc_tot"].to("km/s/Myr").values
plt.figure(figsize=(5, 4))
bins = np.linspace(0,10, 100)
mass_tot = np.sum(mass_all)
anorm = np.sqrt(np.sum(atot_tracers_all**2, axis=1))
weights_all= mass_all/mass_tot
mean_all = np.sum(anorm*weights_all)
mi.linehist(anorm, weights_all, bins=bins, color="k", label="all gas", lw=3)
plt.vlines(mean_all, 0, 0.06, color="k", lw=1, label=f"mean: {mean_all:.1f} km/s/Myr")
mass_gravdom = np.copy(mass_all)
mass_gravdom[np.logical_not(mask_tracer_gravdom)] = 0
weights_grav = mass_gravdom/mass_tot
mi.linehist(anorm, weights=weights_grav, bins=bins, label="gravity-driven gas", color="k", ls='--', lw=3)
mean_grav = np.sum(anorm*mass_gravdom)/np.sum(mass_gravdom)
plt.vlines(mean_grav, 0, 0.009, color="k", lw=1, label=f"mean: {mean_grav:.1f} km/s/Myr", ls="--")
results[iend]["comp_acc"] = {"anorm": anorm,
"weights_all": weights_all,
"weights_grav": weights_grav,
"bins": bins,
"mean_all": mean_all,
"mean_grav": mean_grav}
plt.xlabel(r"$a_\mathrm{tot}$ [km/s/Myr]")
plt.ylabel("Mass weighted distribution function")
plt.legend()
plt.yscale("log")
plt.savefig("distrib_atot.pdf")
def plot_mass_inflow(mass_inflow, iend):
Hcc = (1.4 * osyris.units("m_p/cm**3")).to("g/cm**3")
sink_threshold = 1e3 * Hcc
CNM = 1e-23
H2 = 2e-22
def density_to_n(density):
return (density/Hcc).magnitude
plt.figure(figsize=(6,5), constrained_layout=True)
ax = plt.gca()
ax_top = ax.twiny()
ax_top.set_xlabel(r"$n$ [cm$^{-3}$]")
ax_top.set_xscale("log")
def convert_lim(ax):
x1, x2 = ax.get_xlim()
ax_top.set_xlim(density_to_n(x1), density_to_n(x2))
ax_top.figure.canvas.draw()
ax.callbacks.connect("xlim_changed", convert_lim)
ax.loglog(mass_inflow["density"], mass_inflow["enclosed_mass_start_mask"]["dist"], label="$t_s$", lw=4, ls='-', color="k")
ax.loglog(mass_inflow["density"], mass_inflow["enclosed_mass_end_mask"]["dist"], label="$t_e$", lw=4, ls='--', color="k")
ax.vlines(sink_threshold.magnitude, 0, 1e6, color="tab:orange", lw=2, ls=":", label="Sink threshold")
ax.axvspan(CNM, H2, 0, 60, color="tab:blue", alpha=0.2, label="CNM")
ax.axvspan(H2, 1.2*sink_threshold.magnitude, 0, 60, color="tab:blue", alpha=0.5, label="Molecular gas")
ax.set_xlim(5e-25, 1.2* sink_threshold.magnitude)
ax.set_xlabel(r"$\rho$ [g/cm$^3$]")
ax.set_ylabel(r"M($\rho^\prime > \rho$ ) [M$_\odot$]")
ax.set_title("Tracer mass within 100 pc of the minimum of the potential")
ax.legend()
plt.savefig(f"mass_enclosed_density_{iend}.pdf", bbox_inches="tight")
plt.figure(figsize=(6,5), constrained_layout=True)
ax = plt.gca()
ax_top = plt.gca().twiny()
ax_top.set_xlabel(r"$n$ [cm$^{-3}$]")
ax_top.set_xscale("log")
ax.callbacks.connect("xlim_changed", convert_lim)
ax.plot(mass_inflow["density"], mass_inflow["flow_mask"]["dist"]/1e3, label="Total", alpha=0.7, lw=4, ls='-', color="k")
ax.plot(mass_inflow["density"], mass_inflow["flow_mask"]["dist_gravdom"]/1e3, label="Gravity-driven", alpha=0.7, lw=4, ls='--', color="k")
ax.vlines(sink_threshold.magnitude, 0, 60, color="tab:orange", lw=2, ls=":", label="Sink threshold")
ax.axvspan(CNM, H2, 0, 60, color="tab:blue", alpha=0.2, label="CNM")
ax.axvspan(H2, 1.2*sink_threshold.magnitude, 0, 60, color="tab:blue", alpha=0.5, label="Molecular gas")
ax.set_xlabel("Density [g / cm$^3$]")
ax.set_ylabel(r"Mass inflow rate [ M$_\odot$ / kyr]")
ax.legend(loc="upper left")
ax.set_title("Within 100 pc of the minimum of potential")
ax.set_xscale("log")
ax.set_xlim(1e-25, 1.2*sink_threshold.magnitude)
ax.set_ylim(0, 45)
ax2 = ax.twinx()
ax2.set_ylabel(r"Fraction of gravity-driven inflow", color="tab:red")
ax2.plot(mass_inflow["density"], mass_inflow["flow_mask"]["dist_gravdom"]/mass_inflow["flow_mask"]["dist"],
label="Gravity-driven / Total", alpha=1, lw=4, ls=':', color="tab:red")
ax2.set_ylim(0, 0.45)
ax2.legend(loc="upper right", frameon=False, ncol=1)
ax2.spines['right'].set_color('tab:red')
ax2.yaxis.label.set_color('tab:red')
ax2.tick_params(axis='y', colors='tab:red')
plt.savefig(f"mass_inflow_density_{iend}.pdf", bbox_inches="tight")
def mass_inflow(iend, alignement_gravdom, alignement_inflow=None):
print(f"Mass inflow for {iend}")
mask_tracer_end = dataset[iend]["part"]["family"].values == 0
# Select tracers within 100 pc of the potential minimum
mask_dist = np.abs(dataset[iend]["part"]["position"][mask_tracer_end].to("pc").norm.values) < 100
# Also select gas not within twice the accretion radius of a sink
mask_dist_sinks = mask_dist & (dataset[iend]["part"]["dist_sinks"][mask_tracer_end].to("pc").values > 8)
mass_inflow = mi.compute_mass_inflow_density(dataset, ind_cell_d, alignement_gravdom, iend,
masks={"dist": mask_dist, "dist_sinks": mask_dist_sinks},
alignement_inflow=alignement_inflow)
results[iend]["mass_inflow"] = mass_inflow
def pdf(iend, alignement):
print(f"PDF for {iend}")
mask_tracer_end = dataset[iend]["part"]["family"].values == 0
mass_tracers = dataset[iend]["part"]["mass"][mask_tracer_end].to("Msun").values
mask_tracer_gravdom = np.isin(dataset[iend]["part"]["identity"].values[mask_tracer_end], alignement[iend]['id_grav_dom'], assume_unique=True)
grav_cells, mass_tracers_gravdom_cells = mi.aggregate(mass_tracers[mask_tracer_gravdom], ind_cell_d[iend][mask_tracer_gravdom])
all_cells, mass_tracers_cells = mi.aggregate(mass_tracers, ind_cell_d[iend])
dataset[iend]["mesh"]["mass_tracers"] = osyris.Array(np.zeros_like(dataset[iend]["mesh"]["mass"].values), unit="Msun")
dataset[iend]["mesh"]["mass_tracers"].values[all_cells] = mass_tracers_cells
dataset[iend]["mesh"]["mass_tracers_gravdom"] = osyris.Array(np.zeros_like(dataset[iend]["mesh"]["mass"].values), unit="Msun")
dataset[iend]["mesh"]["mass_tracers_gravdom"].values[grav_cells] = mass_tracers_gravdom_cells
dataset[iend]["mesh"]["frac_grav_dom"] = osyris.Array(np.zeros_like(dataset[iend]["mesh"]["mass"].values))
dataset[iend]["mesh"]["frac_grav_dom"] = dataset[iend]["mesh"]["mass_tracers_gravdom"] / dataset[iend]["mesh"]["mass_tracers"]
mask_zero = dataset[iend]["mesh"]["mass_tracers"].values == 0
dataset[iend]["mesh"]["frac_grav_dom"].values[mask_zero] = 0
mask_dist_mesh = np.abs(dataset[iend]["mesh"]["position"].to("pc").norm.values) < 100
T_cells = dataset[iend]["mesh"]["temperature"].to("K").values[mask_dist_mesh]
masks_cells_end = {"cnm": T_cells < 250, "wnm": T_cells > 4000}
masks_cells_end["lnm"] = ~masks_cells_end["cnm"] & ~masks_cells_end["wnm"]
plt.figure()
mass = dataset[iend]["mesh"]["mass"].values[mask_dist_mesh]
density = dataset[iend]["mesh"]["density"].values[mask_dist_mesh]
mass_gravdom = dataset[iend]["mesh"]["frac_grav_dom"].values[mask_dist_mesh] * mass
bins = np.linspace(-26, -20, 100)
mi.linehist(np.log10(density), weights=mass, bins=bins, lw=3, label="All", color="k", alpha=0.5)
mi.linehist(np.log10(density), weights=mass_gravdom, bins=bins, lw=3, label="gravity-dominated", color="k", ls="--")
mi.linehist(np.log10(density[masks_cells_end["cnm"]]), weights=mass[masks_cells_end["cnm"]], bins=bins, lw=3, alpha=0.5, label="cnm", color="tab:blue")
mi.linehist(np.log10(density[masks_cells_end["lnm"]]), weights=mass[masks_cells_end["lnm"]], bins=bins, lw=3, alpha=0.5, label="lnm", color="tab:green")
mi.linehist(np.log10(density[masks_cells_end["wnm"]]), weights=mass[masks_cells_end["wnm"]], bins=bins, lw=3, alpha=0.5, label="wnm", color="tab:red")
mi.linehist(np.log10(density[masks_cells_end["cnm"]]), weights=mass_gravdom[masks_cells_end["cnm"]], bins=bins, lw=3, ls="--", color="tab:blue")
mi.linehist(np.log10(density[masks_cells_end["lnm"]]), weights=mass_gravdom[masks_cells_end["lnm"]], bins=bins, lw=3, ls="--", color="tab:green")
mi.linehist(np.log10(density[masks_cells_end["wnm"]]), weights=mass_gravdom[masks_cells_end["wnm"]], bins=bins, lw=3, ls="--", color="tab:red")
plt.xlabel("log Density [g / cm$^3$]")
plt.ylabel("Mass [M$_\odot$]")
plt.yscale("log")
plt.ylim(1, 1e5)
plt.xlim(-25.7, -20.5)
plt.legend(loc="upper left")
plt.savefig(f"pdf_dist_{iend}.pdf", bbox_inches="tight")
masks_100K = (T_cells < 30)
plt.figure()
bins = np.linspace(-24, -20, 50)
mi.linehist(np.log10(density[masks_100K]), weights=mass[masks_100K], bins=bins, lw=3, color="#283383", label="10-30 K", density=True)
mi.linehist(np.log10(density[masks_100K]), weights=mass_gravdom[masks_100K], bins=bins, lw=3, color="#283383", label="gravity-dominated", ls="--", density=True)
plt.ylim(1e-2, 1.5)
plt.xlim(-22.7, -20.3)
plt.xlabel("log Density [g / cm$^3$]")
plt.ylabel("Mass [M$_\odot$]")
plt.yscale("log")
plt.legend(loc="upper left")
plt.savefig(f"pdf_100K_{iend}.pdf", bbox_inches="tight")
x_span, y_span, z_span = 50, 50, 50
min_dens = 1e-22
mask_los = ((np.abs(dataset[iend]["mesh"]["position"].to("pc").values[:,0] - x) < x_span)
& (np.abs(dataset[iend]["mesh"]["position"].to("pc").values[:,1] - y) < y_span)
& (np.abs(dataset[iend]["mesh"]["position"].to("pc").values[:,2] - z) < z_span))
velocity_x = dataset[iend]["mesh"]["rel_velocity"].to("km/s").values[:,0][mask_los]
density = dataset[iend]["mesh"]["density"].values[mask_los]
mass = dataset[iend]["mesh"]["mass"].values[mask_los]
mass_gravdom = dataset[iend]["mesh"]["frac_grav_dom"].values[mask_los] * mass
plt.figure(figsize=(6, 4), constrained_layout=True)
mask_dens = (density > min_dens) # & (density < 1e-21)
mass_tot_dens = np.sum(mass[mask_dens])
bins = np.linspace(-11, 11, 60)
bin_width = np.diff(bins)[0]
x,y = mi.linehist(velocity_x[mask_dens], weights=mass[mask_dens], bins=bins, density=True, lw=3, color="k", label="All")
x_g, y_g = mi.linehist(velocity_x[mask_dens], weights=mass_gravdom[mask_dens]/(mass_tot_dens*bin_width), bins=bins, lw=3, ls='--', color="k", label="Gravity-dominated")
x_ng, y_ng = mi.linehist(velocity_x[mask_dens], weights=(mass[mask_dens] - mass_gravdom[mask_dens])/(mass_tot_dens*bin_width), bins=bins, lw=3, ls=':', color="k", label="non gravity-dominated")
results[iend]["pdf_vx"] = {"x": x, "y": y, "x_g": x_g, "y_g": y_g, "x_ng": x_ng, "y_ng": y_ng}
plt.legend()
plt.title(rf"$\rho > 10^"+ "{" + f"{np.log10(min_dens):.0f}" + "}$ g/cm$^3$")
plt.xlabel("$v_x$ [km/s]")
plt.ylabel("Distribution function")
plt.savefig(f"linewidth_{iend}.pdf", bbox_inches="tight")
def process_snapshot(end, start_gravdom, start_inflow=None, free_snap=None):
if start_inflow is None:
start_inflow = start_gravdom
if end not in dataset:
mi.load_snapshot(dataset, ind_cell_d, path, end)
if start_gravdom not in dataset:
mi.load_snapshot(dataset, ind_cell_d, path, start_gravdom)
if start_inflow not in dataset:
mi.load_snapshot(dataset, ind_cell_d, path, start_inflow)
if free_snap is None:
free_snap = end
alignement[start_gravdom] = {}
if end != start_gravdom:
print(f"Compute alignement {end}")
alignement[start_gravdom][end] = mi.get_aligment(dataset, ind_cell_d, start_gravdom, end)
if start_gravdom != start_inflow:
alignement[start_inflow]
print(f"Compute alignement for inflow from {start_inflow} to {end}")
alignement[start_inflow][end]= mi.get_aligment(dataset, ind_cell_d, start_inflow, end)
results[end] = {}
print(f"Compute results for {end}")
mass_fraction(end, alignement[start_gravdom])
comp_acc(end, alignement[start_gravdom])
mass_inflow(end, alignement[start_gravdom], alignement_inflow=alignement[start_inflow])
pdf(end, alignement=alignement[start_gravdom])
time = dataset[end].meta["time"].to("Myr").magnitude
print(f"Save results for {end} [time={time:.2f} Myr]")
results[end]["time"] = time
hkl.dump(results[end], f"results_{start_gravdom}_{start_inflow}_{end}.h5")
if free_snap in dataset:
print(f"Free memory for {free_snap}")
del dataset[free_snap]
del ind_cell_d[free_snap]
if free_snap in alignement:
del alignement[free_snap]
del results[free_snap]
dataset = {}
alignement = {}
ind_cell_d = {}
results = {}
snaps = np.arange(173, 186)
start = 173
end = 185
snaps = [173, 180]
path = "~/simus/mass_inflow/n0.66_rms00000_tracers_memory_70/"
for iend in snaps:
process_snapshot(iend, start_gravdom=start, start_inflow=start, free_snap=iend)