noe/pipeline
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

black (still without mortimer)

Browse Source
This commit is contained in:
Noe Brucy
2023-01-30 12:12:14 +01:00
parent 8c3db9b7cb
commit 1a8db42d29
11 changed files with 474 additions and 99 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.flake8
+2
View File
@@ -7,3 +7,5 @@ ignore =
W503
# Invalid escape sequence '\['
W605
# Imported but unused
F401
ismfeed/ismfeed.py
+419 -51
View File
@@ -1,13 +1,32 @@
from scipy.integrate import solve_ivp
from plotter import U
import utils.snapshotselector
from snapshotprocessor import mean_by_bins
from utils import snapshotselector as select_snapshot
import numpy as np
import pandas as pd
import os
import matplotlib as mpl
ssfr_base = 2.5e-10
import matplotlib.ticker as tick
import matplotlib.pyplot as plt
from matplotlib import FuncFormatter
from utils.units import convert_exp
ssfr_sun = 2.5e-3
mp = 1.4 * 1.66 * 10 ** (-24) * U.g
z0 = 150 * U.pc
def get_sinks(self, overwrite=False, stellar=True, convert_units=False):
def convert_coldens_s(n0, z0=z0):
return (np.sqrt(2 * np.pi) * mp * z0 * (n0 * U.cm ** (-3))).express(U.coldens)
convert_coldens = np.vectorize(convert_coldens_s)
def get_sinks(self, overwrite=False, stellar=True, convert_units=False, sk=True):
self.sinks_from_log(overwrite=overwrite)
self.coarse_step_from_log(overwrite=overwrite)
@@ -18,7 +37,9 @@ def get_sinks(self, overwrite=False, stellar=True, convert_units=False):
self.stellar = self.get_value("/dataset/stellar_from_log")
self.stellar["sn_time"] = {}
for run in self.runs:
self.stellar["sn_time"][run] = (self.stellar["time"][run] + self.stellar["lifetime"][run] )*self.info["unit_time"].express(U.year)
self.stellar["sn_time"][run] = (
self.stellar["time"][run] + self.stellar["lifetime"][run]
) * self.info["unit_time"].express(U.year)
ind_sort = self.stellar["sn_time"][run].argsort()
for key in self.stellar:
self.stellar[key][run] = self.stellar[key][run][ind_sort]
@@ -32,32 +53,39 @@ def get_sinks(self, overwrite=False, stellar=True, convert_units=False):
for run in self.runs:
if convert_units:
self.sinks["time"][run] *= self.info["unit_time"].express(U.year)
self.sinks["time"][run] *= self.info["unit_time"].express(U.year) # year
self.sinks["cum_mass"][run] = self.sinks["mass_sink"][run].copy()
self.sinks["cum_mass"][run] = self.sinks["mass_sink"][run].copy() # Msun
if stellar:
self.stellar["cum_mass"][run] = np.cumsum(self.stellar["mass"][run]) * self.info["unit_mass"].express(U.Msun)
sn["time"][run], idx, count = np.unique(self.stellar["sn_time"][run], return_index=True, return_counts=True)
self.stellar["cum_mass"][run] = np.cumsum(
self.stellar["mass"][run]
) * self.info["unit_mass"].express(U.Msun)
sn["time"][run], idx, count = np.unique(
self.stellar["sn_time"][run], return_index=True, return_counts=True
)
idx += count - 1
sn["cum_mass"][run] = self.stellar["cum_mass"][run][idx]
ind_sn = np.searchsorted(self.sinks["time"][run], sn["time"][run])
for i, indi in enumerate(ind_sn[ind_sn < ind_sn[-1]]):
self.sinks["cum_mass"][run][indi:ind_sn[i+1]] += sn["cum_mass"][run][i]
self.sinks["cum_mass"][run][indi : ind_sn[i + 1]] += sn["cum_mass"][
run
][i]
time_gas, total_mass, mass_gas = self.total_mass()
if sk:
time_gas, total_mass, mass_gas = self.total_mass() # year, Msun, Msun
sk_ssfr = {}
for run in self.runs:
time_gas[run] *= self.info["unit_time"].express(U.Myr)
sk_ssfr[run] = ssfr_base * 1e6 * (mass_gas[run]/1e6)**1.4
coldens = mass_gas[run] / 1e6 # in Msun/pc²
time_gas[run] *= self.info["unit_time"].express(U.Myr) # Myr
sk_ssfr[run] = ssfr_sun * (coldens / 10) ** 1.4 # Msun/kpc²/yr-1
def dermass(t, sm, run, fact_sfr=1):
ind = min(np.searchsorted(time_gas[run], t), sk_ssfr[run].size - 1)
return sk_ssfr[run][ind]*fact_sfr
return sk_ssfr[run][ind] * fact_sfr * 1e6 # Msun/Myr-1
tspan = np.linspace(0, 200, 100)
tspan = np.linspace(0, 200, 100) # Myr
self.time_esm = {}
self.expected_stellar_mass = {}
@@ -65,66 +93,406 @@ def get_sinks(self, overwrite=False, stellar=True, convert_units=False):
self.expected_stellar_mass_min = {}
for run in self.runs:
sol = solve_ivp(dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run])
self.time_esm[run], self.expected_stellar_mass[run] = sol["t"], sol["y"][0]
sol_max = solve_ivp(dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run, 3])
self.expected_stellar_mass_max[run] = sol_max["y"][0]
sol_min = solve_ivp(dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run, 1./3])
self.expected_stellar_mass_min[run] = sol_min["y"][0]
sol = solve_ivp(
dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run]
)
self.time_esm[run], self.expected_stellar_mass[run] = (
sol["t"],
sol["y"][0],
) # Myr, Msun
sol_max = solve_ivp(
dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run, 3]
)
self.expected_stellar_mass_max[run] = sol_max["y"][0] # Myr, Msun
sol_min = solve_ivp(
dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run, 1.0 / 3]
)
self.expected_stellar_mass_min[run] = sol_min["y"][0] # Myr, Msun
def compute_sfr(self, target_start=0.05, target_end=0.3):
mass = self.sinks["cum_mass"]
time = self.sinks["time"]
mass = self.sinks["cum_mass"] # Msun
time = self.sinks["time"] # year
sfr = {}
sfr_err = {}
tend = {}
tstart = {}
for run in self.runs:
tend[run] = select_runs.time_mcons(self, run, target=target_end) * 1e6
tstart[run] = select_runs.time_mcons(self, run, target=target_start) * 1e6
tend[run] = (
select_snapshot.time_mcons(self, run, target=target_end) * 1e6
) # year
tstart[run] = (
select_snapshot.time_mcons(self, run, target=target_start) * 1e6
) # year
if tend[run] < tstart[run]:
tend[run] = time[run][-1]
for run in self.runs:
idx1 = time[run].searchsorted(tstart[run])
idx2 = min(time[run].searchsorted(tend[run]), len(time[run]) - 1)
sfr[run] = (mass[run][idx2] - mass[run][idx1]) / (time[run][idx2] - time[run][idx1])
sfr[run] = (mass[run][idx2] - mass[run][idx1]) / (
time[run][idx2] - time[run][idx1]
) # Msun/year
sfr_other = []
for i in range(idx1, idx2 - 10):
sfr_other.append((mass[run][i+10:idx2+1] - mass[run][i]) / (time[run][i+10:idx2+1] - time[run][i]))
sfr_other.append(
(mass[run][i + 10 : idx2 + 1] - mass[run][i])
/ (time[run][i + 10 : idx2 + 1] - time[run][i])
) # Msun/year
sfr_err[run] = np.std(np.concatenate(sfr_other))
self.sfr = np.array(list(sfr.values()))
self.sfr_err = np.array(list(sfr_err.values()))
sfr_err[run] = np.std(np.concatenate(sfr_other)) # Msun/year
self.sfr = np.array(list(sfr.values())) # Msun/year
self.sfr_err = np.array(list(sfr_err.values())) # Msun/year
return self.sfr, self.sfr_err
def compute_power(self, target_start=0.05, target_end=0.3):
mass = self.sinks["cum_mass"]
time = self.sinks["time"]
sfr = {}
sfr_err = {}
def get_nml_array(
pl,
xkey,
cmap=None,
logcmap=False,
put_colorbar=True,
cbar_fmt=tick.FormatStrFormatter("%.2g"),
):
pl.study.nml(xkey, overwrite=True)
x = pl.study.get_value(f"/comp/nml_{xkey}")
xname = os.path.basename(xkey)
if xname in pl.value_convert:
x = np.array([pl.value_convert[xname](x_v) for x_v in x])
if xname in pl.label_convert:
xlabel = pl.label_convert[xname]
else:
xlabel = xname
colors = None
if cmap is not None:
if logcmap:
sm = plt.cm.ScalarMappable(
cmap=cmap, norm=mpl.colors.LogNorm(vmin=min(x), vmax=max(x))
)
else:
sm = plt.cm.ScalarMappable(
cmap=cmap, norm=mpl.colors.Normalize(vmin=min(x), vmax=max(x))
)
print(x)
if put_colorbar:
cb = plt.colorbar(sm)
cb.set_label(xlabel)
cb.ax.yaxis.set_major_formatter(cbar_fmt)
def colors(xi):
return sm.cmap(sm.norm(xi))
return x, xname, xlabel, colors
def histo_speed(pli, redo=False):
fig, axes = plt.subplots(1, 4, figsize=(16, 5), sharey=True)
fig2, axes2 = plt.subplots(1, 4, figsize=(16, 5), sharey=True)
fig3, axes3 = plt.subplots(1, 4, figsize=(16, 5), sharey=True)
direction = ["x", "y", "z"]
for pp, ax, ax2, ax3 in zip(
pli.get_snap_list(
select={"filter_nml": (pli.nml_key[0], "in", [1.5, 3, 6, 12])}
),
axes,
axes2,
axes3,
):
pp.load_cells()
vel = pp.cells["vel"][:] * pp.info["unit_velocity"].express(U.km_s)
pos = pp.cells["pos"][:] * pp.info["unit_length"].express(U.pc)
if redo:
pp.vmean_z = []
pp.vmean_x = []
mass = (
pp.cells["rho"]
* pp.info["unit_density"].express(U.Msun / U.pc**3)
* (pp.cells["dx"] * pp.info["unit_length"].express(U.pc)) ** 3
)
for i in range(3):
ax.hist(
vel[:, i],
range=[-100, 100],
bins=100,
histtype="step",
density=True,
ls="-",
label=f"$v_{direction[i]}$",
lw=2,
)
ax.set_yscale("log")
ax.set_ylim(1e-3, 2e-1)
ax.set_xlabel("$v$ [km/s]")
ax.set_title(pli.get_label_run(nml_key=pli.nml_key[:2], run=pp.run))
ax2.set_title(pli.get_label_run(nml_key=pli.nml_key[:2], run=pp.run))
ax3.set_title(pli.get_label_run(nml_key=pli.nml_key[:2], run=pp.run))
if redo:
zbin, vmean_z = mean_by_bins(pos[:, 2], vel[:, i], weights=mass)
pp.zbin = zbin
pp.vmean_z.append(vmean_z)
xbin, vmean_x = mean_by_bins(pos[:, 0], vel[:, i], weights=mass)
pp.xbin = zbin
pp.vmean_x.append(vmean_x)
else:
zbin = pp.zbin
xbin = pp.xbin
vmean_x = pp.vmean_x[i]
vmean_z = pp.vmean_z[i]
ax2.plot(zbin, vmean_z, label=f"$v_{direction[i]}$", lw=2)
ax2.set_xlabel("z [pc]")
ax3.plot(xbin, vmean_x, label=f"$v_{direction[i]}$", lw=2)
ax3.set_xlabel("x [pc]")
fig.suptitle(pli.name)
fig2.suptitle(pli.name)
fig3.suptitle(pli.name)
axes[0].set_ylabel("Mass weighted PDF")
axes[0].legend()
axes2[0].set_ylabel("v [km/s]")
axes2[0].legend()
axes3[0].set_ylabel("v [km/s]")
axes3[0].legend()
# def make_clump_hop(self, threshold_density=10):
# """
# Apply HOP algorithm
# Args:
# threshold_density (float): select only cells over threshold
# """
# # Selection of cells
# mask = (
# self.cells["rho"] * self.info["unit_density"].express(U.H_cc)
# > threshold_density
# )
# ncells = np.sum(mask)
# # fill the matrice with ID, x,y,z and masses of particles
# cells_group = np.zeros((ncells, 14))
# cells_group[:, 0] = np.arange(ncells) # index
# position = self.cells["pos"][mask] * self.info["unit_length"].express(U.pc)
# cells_group[:, 1:4] = position # position
# density = self.cells["rho"][mask] * self.info["unit_density"].express(U.H_cc)
# size = self.cells["dx"][mask] * self.info["unit_length"].express(U.pc)
# cells_group[:, 4] = (
# density * size ** 3 * (U.H_cc * U.pc ** 3).express(U.Msun)
# ) # mass
# cells_group[:, 6] = density
# velocity = self.cells["vel"][mask] * self.info["unit_velocity"].express(U.km_s)
# pressure = self.cells["P"][mask] * self.info["unit_pressure"].express(U.bar)
# temperature = (pressure / density) * (
# (U.bar / U.H_cc) * (1.4 * U.mH) / U.kB
# ).express(U.K)
# cells_group[:, 7] = temperature
# cells_group[:, 8] = pressure
# cells_group[:, 9:12] = velocity
# cells_group[:, 12] = self.cells["phi"][mask]
# cells_group[:, 13] = size ** 3
# # save file.txt
# head = str(ncells)
# np.savetxt(
# self.filename[:-3] + "_hop.txt",
# cells_group[:, :5],
# fmt="%10d %.10e %.10e %.10e %.10e",
# header=head,
# delimiter=" ",
# comments=" ",
# )
# # save file.den
# f = open(self.filename[:-3] + "_hop.den", "wb")
# f.write(pack("I", ncells))
# self.cells["rho"][mask].astype("f").tofile(f)
# f.close()
# # exec HOP algo
# h = HOP(self.filename[:-3] + "_hop.txt", os.path.dirname(self.filename))
# h.process_hop(force=True)
# # get the igroup array
# group_ids = h.get_group_ids()
# # sort it and apply the sorting to the coordinates
# # this means that the particules of group 1 are written first then of group 2 etc...
# ind_sort = np.argsort(group_ids)
# cells_group = cells_group[ind_sort]
# cells_group[:, 5] = group_ids[ind_sort]
# # Make it a pandas' DataFrame
# cells_group = pd.DataFrame(
# cells_group,
# columns=[
# "id",
# "x",
# "y",
# "z",
# "mass",
# "group",
# "density",
# "temperature",
# "pressure",
# "vx",
# "vy",
# "vz",
# "phi",
# "volume",
# ],
# )
# self.clumps = cells_group
# return cells_group
def plot_mass_sfr(
pl,
xkey="turb_params/comp_frac",
start=0.01,
end=0.4,
cmap=None,
logcmap=False,
redo=False,
fig_height=4,
ax_sfr=None,
logsfr=True,
marker="o",
color_sfr="k",
scale=0.9,
do_sfr=True,
plot_dir=".",
):
plt.figure(figsize=np.array([5, fig_height]) * scale)
get_sinks(pl.study, overwrite=redo)
pl.value_convert["Bx"] = lambda x: x * 7.6189439
pl.label_convert["Bx"] = r"B$_0$ [$\mu$G]"
pl.value_convert["bx_bound"] = lambda x: x * 7.6189439
pl.label_convert["bx_bound"] = r"B$_0$ [$\mu$G]"
pl.label_convert["comp_frac"] = r"$\chi$"
x_fmt = FuncFormatter(lambda x, p: f"{convert_exp(x,2)}")
x, xname, xlabel, colors = get_nml_array(pl, xkey, cmap, logcmap, cbar_fmt=x_fmt)
try:
n0 = pl.study.get_nml("galbox_params/dens0", pl.runs[0])
except KeyError():
n0 = pl.study.get_nml("cloud_params/dens0", pl.runs[0])
mass = pl.study.sinks["cum_mass"]
time = pl.study.sinks["time"]
tend = {}
tstart = {}
for run in self.runs:
tend[run] = select_runs.time_mcons(self, run, target=target_end) * 1e6
tstart[run] = select_runs.time_mcons(self, run, target=target_start) * 1e6
for i, run in enumerate(pl.runs):
tend[run] = (select_snapshot.time_mcons(pl.study, run, target=end) * 1e6,)
tstart[run] = select_snapshot.time_mcons(pl.study, run, target=start) * 1e6
if tend[run] < tstart[run]:
tend[run] = time[run][-1]
idx1 = time[run].searchsorted(tend[run])
idx2 = time[run].searchsorted(tstart[run])
mask = time[run] < tend[run]
for run in self.runs:
idx1 = time[run].searchsorted(tstart[run])
idx2 = min(time[run].searchsorted(tend[run]), len(time[run]) - 1)
if cmap is None:
(p,) = plt.plot(
time[run][mask] * 1e-6,
mass[run][mask] * 1e-6,
label=f"{xlabel} = {x[i]:.2f}",
)
else:
color = colors(x[i])
(p,) = plt.plot(time[run][mask] * 1e-6, mass[run][mask] * 1e-6, color=color)
plt.scatter(
[time[run][idx1] * 1e-6],
[mass[run][idx1] * 1e-6],
marker="<",
color=p.get_color(),
)
plt.scatter(
[time[run][idx2] * 1e-6],
[mass[run][idx2] * 1e-6],
marker=">",
color=p.get_color(),
)
sfr[run] = (mass[run][idx2] - mass[run][idx1]) / (time[run][idx2] - time[run][idx1])
sfr_other = []
for i in range(idx1, idx2-10):
sfr_other.append((mass[run][i+10:idx2+1] - mass[run][i]) / (time[run][i+10:idx2+1] - time[run][i]))
plt.title(f"$\Sigma_0$ = {convert_coldens(n0):.1f}" + " M$_\odot$.pc$^{-1}$")
plt.xlabel("Time [Myr]")
plt.ylabel("Mass accreted [$10^6$ M$_\odot$]")
plt.savefig(f"{plot_dir}/{xname}_n{n0}_mass.pdf")
sfr_err[run] = np.std(np.concatenate(sfr_other))
self.sfr = np.array(list(sfr.values()))
self.sfr_err = np.array(list(sfr_err.values()))
return self.sfr, self.sfr_err
##################################################
if do_sfr:
label = f"$\Sigma_0$ = {convert_coldens(n0):.1f}" + " M$_\odot$.pc$^{-2}$"
if ax_sfr is None:
plt.figure(figsize=(5, fig_height))
else:
plt.sca(ax_sfr)
if redo or not hasattr(pl.study, "sfr"):
sfr, sfr_err = compute_sfr(pl.study, target_start=start, target_end=end)
else:
sfr, sfr_err = pl.study.sfr, pl.study.sfr_err
if logsfr:
plt.errorbar(
x=x,
y=np.log10(sfr),
yerr=sfr_err / sfr,
color=color_sfr,
marker=marker,
ls=":",
lw=0.5,
label=label,
elinewidth=1.5,
)
plt.hlines(
np.log10(
ssfr_sun * (convert_coldens(n0) * (1 - end / 2.0) / 10) ** 1.4
),
xmin=min(x),
xmax=max(x),
color=color_sfr,
ls="--",
)
plt.ylabel("log$(\Sigma_{\mathrm{SFR}})$ " + rf"[${U.ssfrK.latex}$]")
else:
plt.errorbar(
x=x,
y=sfr,
yerr=sfr_err,
color=color_sfr,
marker=marker,
ls=":",
lw=0.5,
label=label,
)
plt.hlines(
ssfr_sun * (convert_coldens(n0) * (1 - end / 2.0) / 10) ** 1.4,
xmin=min(x),
xmax=max(x),
color=color_sfr,
ls="--",
)
plt.ylabel("$\Sigma_{\mathrm{SFR}}$ " + rf"[${U.ssfrK.latex}$]")
plt.xlabel(f"{xlabel}")
plt.gca().xaxis.set_major_formatter(x_fmt)
plt.legend()
plt.savefig(f"{plot_dir}/{xname}_n{n0}_logsfr.pdf")
turbox/__init__.py
+7 -2
View File
@@ -1,2 +1,7 @@
from .turbox import get_pspec, get_pspec_slope, build_suite, apply_rule_pdf, span_resolution
from .turbox import (
get_pspec,
get_pspec_slope,
build_suite,
apply_rule_pdf,
span_resolution,
)
turbox/turbox.py
+1 -1
View File
@@ -155,7 +155,7 @@ def build_suite(pl, redo=False, cs0=0.28834810480560674):
)
)
df["sigma"] = list(map(lambda l: np.mean(l), df["sigma_all"].values))
df["sigma"] = list(map(lambda sig_list: np.mean(sig_list), df["sigma_all"].values))
df["Mach"] = df["sigma"] / cs0
df["turnover"] = (df["L"] * U.pc.express(U.km) / (2 * df["sigma"])) * U.s.express(
U.Myr
utils/params.py
+3 -3
View File
@@ -11,9 +11,9 @@ _dir_path = os.path.dirname(os.path.realpath(__file__))
# Add support for '1e3' kind of float
_loader = yaml.SafeLoader
_loader.add_implicit_resolver(
u"tag:yaml.org,2002:float",
"tag:yaml.org,2002:float",
re.compile(
u"""^(?:
"""^(?:
[-+]?(?:[0-9][0-9_]*)\\.[0-9_]*(?:[eE][-+]?[0-9]+)?
|[-+]?(?:[0-9][0-9_]*)(?:[eE][-+]?[0-9]+)
|\\.[0-9_]+(?:[eE][-+][0-9]+)?
@@ -22,7 +22,7 @@ _loader.add_implicit_resolver(
|\\.(?:nan|NaN|NAN))$""",
re.X,
),
list(u"-+0123456789."),
list("-+0123456789."),
)