from __future__ import annotations
import os.path
import shutil
import numpy as np
import phenigraph as g
from scipy.integrate import trapezoid, cumulative_trapezoid, solve_ivp
from scipy.interpolate import interp1d, make_smoothing_spline, RegularGridInterpolator
from scipy.signal import medfilt2d, savgol_filter
from diffenix.Constantes import *
from diffenix.solve_asteroid_belt import SolrhoBelt
from diffenix.numericals_methods import gradient, gradient2, grid_refine_inner_edge
import radmc3dPy.image as radmc_im
import radmc3dPy.analyze as radmc_a
# noinspection PyUnresolvedReferences,PyTypeChecker
[docs]
def merge_y(y: np.ndarray) -> np.ndarray:
"""
Do the inverse operation of split_y
Parameters
----------
y : np.ndarray
The array to be reshaped
Returns
-------
np.ndarray
The new array
"""
return y.flatten()
[docs]
def new_space_grid(old_grid: np.ndarray[np.float64], last_masse_frac: np.ndarray[np.float64]) -> np.ndarray[np.float64]:
if np.any(abs(1-last_masse_frac.sum(axis=0)) > 1e-8):
idx: np.ndarray = np.argwhere(abs(1-last_masse_frac.sum(axis=0)) > 1e-8)[:, 0]
res: list[np.float64] = []
for i in range(len(idx)):
if idx[i] > 0:
if i == 0:
res.extend(list(old_grid[:idx[i]]))
else:
res.extend(list(old_grid[idx[i-1]:idx[i]]))
if abs(1 - last_masse_frac[:, i].sum(axis=0)) < 1e-6:
res.append(old_grid[idx[i] -1] + (old_grid[idx[i]] - old_grid[idx[i] -1]) / 2.)
else:
res.extend(list(np.linspace(old_grid[idx[i] -1], old_grid[idx[i]], 3)[1:-1]))
if idx[-1] < len(old_grid) - 1:
res.extend(list(old_grid[idx[-1]:]))
return np.array(res)
else:
return old_grid
[docs]
def split_y(y: np.ndarray, len_radius: int, nb_species: int) -> np.ndarray:
"""
Split y the argument of the time integration system in n lists, n the number of chemical elements
Parameters
----------
y : np.ndarray
The array to be reshaped
len_radius : int
The lenght of the radius array
nb_species : int
The number of chemicals species
Returns
-------
np.ndarray
masses_fractions(y)
"""
if len(y) != len_radius * nb_species:
raise UserWarning("split_y : the dimension of y doesn't mach with the number of element and the size of the"
" radius list for this Sol", len(y), len_radius * nb_species)
return y.reshape((nb_species, len_radius))
mdot_st: np.float64 =(1e-20 * Mearth / Myr)
# mdot_st: np.float64 = m_dot
[docs]
class Sol:
"""
The object contains all the information needed tu run and exploits a viscous diffusion simulation.
It can be saved into a compressed .npz file.
This object can generate phenigraph's Graphique to analyze the results.
"""
compteur: int = 0
def __init__(self, filename: str = "", consts: dict = None, directory="",
mdot_st: np.float64 = mdot_st,
mus0: np.ndarray = None, mass_fracs0: np.ndarray = None,
noms_mus0: np.ndarray[str] = None):
"""
Initialization of a new Sol object
Parameters
----------
filename: str, optional, default=None
The name of the .npz containing the object.
If not provide a new object is created
consts: dict, optional, default=None
The dictionary containing all the constants and parameters.
If not provided, use the globals constantes
directory: str, optional, default=None
The path top the directory containing the solution relatively to the 'valeur' directory
mdot_st : np.ndarray, optional, default=10**-10 Mearth/Myr
The initial state is assumed to be a steady state generated by a punctual belt with a mass generation belt
mdot_st
mus0 : np.ndarray, optional,
The molar mass (in code unit) of the elements initially in the gas disc,
default CO if there is a kuiper belt, H in the other cases
noms_mus0 : np.ndarray, optional,
The name of the initials elements, default ["CO"] if there is a kuiper belt, ["H"] in others cases
mass_fracs0 : np.ndarray, optional, default=[1.]
The mass fraction of initials elements can be an array of the size of mus0, then the mass fractions will
be equals for each radius or an array of shape (len(mus0), IT_s)
"""
# This calcul is made once and for all and is used to be normalized the surface density if needed
# This object is used to get the mass generation rate by the belt
self.sol_ast = None
self.sol_kuip = None
if filename == "":
if noms_mus0 is not None and len(noms_mus0) != len(mus0):
raise UserWarning("The size of the initial element array mus0 (",
len(mus0), ") must be equal to the size of the "
"initial element names array nom_mus0 (", len(noms_mus0), ")")
if "directory" in consts.keys():
self.directory = str(consts["directory"])
else:
self.directory: str = results
if consts is None or consts == {}:
self.C_L: np.float64 = C_L
self.C_M: np.float64 = C_M
self.C_t: np.float64 = C_t
self.C_E: np.float64 = C_E
self.C_F: np.float64 = C_F # N-> : force
self.C_rho = np.double(self.C_M / (self.C_L * self.C_L * self.C_L)) # kg/M**6 -> : density
self.Ms: np.float64 = Ms # Masse de l'étoile
self.Ts: np.float64 = Ts # Température de l'étoile
self.Rs: np.float64 = Rs # Rayon de l'étoile
self.Ls: np.float64 = Ls # Rayon de l'étoile
self.T0: np.float64 = T0 # Central temperature of the belt
self.pls_temp: np.float64 = pls_temp # Radial power slope of temperature
self.mH_init: np.float64 = zero # Initial primordial hydrogen mass
self.a0: np.float64 = a0 # Central radius of the belt
self.a_in: np.float64 = a_in # Inner radius
self.a_out: np.float64 = a_out # Outer radius
self.alpha: np.float64 = alpha # Viscous parameter
self.IT_s: int = IT_s # Number of spatial steps
self.dip: int = dip # Half-number of sink cells to model the accretion of a planet
self.rtol: np.float64 = rtol # précision relative attendue pour la résolution temporelle
self.atol: np.float64 = atol # précision relative attendue pour la résolution temporelle
self.gamma: np.float64 = gamma
self.method: str = method # Time integrator solveur (RK45-DOP853-LSODA-BDF)
# cf documentation scipy.integrate.solve_ivp
self.a_planets: np.ndarray = a_planets.copy() # The semi-majors axis of planets
self.m_planets: np.ndarray = m_planets.copy() # The mass of planets
self.f_accr: np.float64 = f_accr # hydrodynamic accretion efficiency
self.t0: np.float64 = t0 # Protoplanetary disc lifetime
self.sol_ast_params = dict(sublimation_model=AB_sublimation_model,
mdot_st_init=np.double(1e-10) * Mearth / Myr,
Mbelt=AB_Mbelt, Ms=Ms,
a0=AB_a0, delta_a=AB_delta_r,
tinit=t0, tmax=Gyr,
rho_refr=AB_rho_refr, rho_ice=AB_rho_ice, f_ice=AB_f_ice,
A_ast=AB_Abelt, K=AB_K,
Its=100, Phi=AB_Phi,
dmin=AB_SD_amin, sb=AB_sb, dmax=AB_SD_amax,
e=AB_e, i=AB_i, As=AB_As, bs=AB_bs, bg=AB_bg,
t0_diss=AB_t0_diss, t1_diss=AB_t1_diss,
Its_p=1000, const_mdot=m_dot,
rp=AB_rp, Itt=10000)
self.sol_ast_kwargs_N = dict(amax=AB_SD_amax,
abig=AB_SD_abig,
amed=AB_SD_amed, amin=AB_SD_amin,
qh=AB_SD_qh, qm=AB_SD_qm, ql=AB_SD_ql,
rho=AB_rho_refr)
self.sol_kuip_params = dict(sublimation_model=KB_sublimation_model,
mdot_st_init=np.double(1e-10) * Mearth / Myr,
Mbelt=KB_Mbelt, Ms=Ms,
a0=KB_a0, delta_a=KB_delta_r,
tinit=t0, tmax=Gyr,
rho_refr=KB_rho_refr, rho_ice=KB_rho_ice, f_ice=KB_f_ice,
A_ast=KB_Abelt, K=KB_K,
Its=60, Phi=KB_Phi,
dmin=KB_SD_amin, sb=KB_sb, dmax=KB_SD_amax,
e=KB_e, i=KB_i, As=KB_As, bs=KB_bs, bg=KB_bg,
t0_diss=KB_t0_diss, t1_diss=KB_t1_diss,
Its_p=1000, const_mdot=m_dot,
rp=KB_rp, Itt=10000)
self.sol_kuip_kwargs_N = dict(amax=KB_SD_amax,
abig=KB_SD_abig,
amed=KB_SD_amed, amin=KB_SD_amin,
qh=KB_SD_qh, qm=KB_SD_qm, ql=KB_SD_ql,
rho=KB_rho_refr)
# Constantes globales
else:
self.update_constantes(consts)
self.name: str = "Simulation_diffusion"
self.r: np.ndarray | None = None # r the radius
self.x: np.ndarray | None = None # r ** (1/2)
self.f: list[np.ndarray] = [] # the logarithm of surface density shape=(len(tps),len(r))
self.mass_fractions: list[np.ndarray] = [] # Sigmas_i/Sigmas, shape=(len(tps_mf),len(mu),len(r)
self.idxa: np.ndarray | None = None # The index of 'r' where the asteroid belt is
self.idxk: np.ndarray | None = None # The index of 'r' where the kuiper belt is
self.tps: list[np.float64] = [] # Times for the surface density integration
self.tps_mf: list[np.float64] = [] # Times for the mass fractions integration
self.idx_mu_ast: np.ndarray[int] = np.array([], dtype=int) # index of the chemicals elements
# associated with the asteroid belt ice
# noinspection PyTypeChecker
self.idx_mu_kuip: np.ndarray[int] = np.array([], dtype=int) # index of the chemicals elements
# associated with the kuiper belt ice
mu: list[np.float64] = [] # molar mass for each chemical element
nom_mu: list[str] = [] # name each chemical element
if self.mH_init > zero:
mu.append(np.double(1. * 1e-3 * C_M))
nom_mu.append("H")
minit_asts: np.float64 = mass_steady_state(sol=self, mdot=self.sol_ast_params["mdot_st_init"],
a0=self.sol_ast_params["a0"], mu=np.double(18. * 1e-3 * C_M))
minit_kuips: np.float64 = mass_steady_state(sol=self, mdot=self.sol_kuip_params["mdot_st_init"],
a0=self.sol_kuip_params["a0"], mu=np.double(28. * 1e-3 * C_M))
minit: np.float64 = self.mH_init + minit_asts + minit_kuips
if self.sol_ast_params["sublimation_model"] != "none":
self.idx_mu_ast = len(mu) + np.arange(3)
mu.extend([np.double(18. * 1e-3 * C_M), np.double(16. * 1e-3 * C_M)]) # molar mass
nom_mu.extend(["H2O", "O"])
if not "H" in nom_mu: # Only if there is no water in the system
mu.append(np.double(1. * 1e-3 * C_M))
nom_mu.append("H")
self.idx_mu_ast = np.append(self.idx_mu_ast, len(mu) - 1)
else:
self.idx_mu_ast = np.append(self.idx_mu_ast, np.argwhere(np.array(nom_mu) == "H")[0, 0])
if self.sol_kuip_params["sublimation_model"] != "none":
self.idx_mu_kuip = len(mu) + np.arange(2)
# mu.extend([np.double(28. * 1e-3 * C_M), np.double(12. * 1e-3 * C_M), np.double(16. * 1e-3 * C_M)]
# ) # molar mass
# nom_mu.extend(["CO", "C", "O"])
mu.extend([np.double(28. * 1e-3 * C_M), np.double(12. * 1e-3 * C_M)]
) # molar mass
nom_mu.extend(["CO", "C"])
if not "O" in nom_mu: # Only if there is no water in the system
mu.append(np.double(16. * 1e-3 * C_M))
nom_mu.append("O")
self.idx_mu_kuip = np.append(self.idx_mu_kuip, len(mu) - 1)
else:
self.idx_mu_kuip = np.append(self.idx_mu_kuip, np.argwhere(np.array(nom_mu) == "O")[0, 0])
if mus0 is None and "CO" in nom_mu and "H2O" not in nom_mu:
mus0 = np.array([np.double(28.) * 1e-3 * C_M,
np.double(12.) * 1e-3 * C_M,
np.double(16.) * 1e-3 * C_M])
if self.mH_init == zero:
noms_mus0 = np.array(["CO", "C", "O"])
mass_fracs0 = np.array([0.5, 0.5 * 12 / 28, 0.5 * 16 / 28])
else:
noms_mus0 = np.array(["CO", "C", "O", "H"])
mass_fracs0 = np.array([(minit_kuips / minit) * 0.5,
(minit_kuips / minit) * 0.5 * 12 / 28,
(minit_kuips / minit) * 0.5 * 16 / 28,
self.mH_init / minit])
self.a0 = self.sol_kuip_params["a0"]
elif mus0 is None and "H2O" in nom_mu and "CO" not in nom_mu:
mus0 = np.array([18e-3 * C_M, 1e-3 * C_M, 16e-3 * C_M])
noms_mus0 = np.array(["H2O", "H", "O"])
mass_fracs0 = np.array([(minit_asts / minit) * 0.5,
(minit_asts / minit) * 0.5 * 2. / 18. + self.mH_init / minit,
(minit_asts / minit) * 0.5 * 16. / 18.])
self.a0 = self.sol_ast_params["a0"]
elif mus0 is None and "H2O" in nom_mu and "CO" in nom_mu:
mus0 = np.array([np.double(28.) * 1e-3 * C_M,
np.double(12.) * 1e-3 * C_M,
np.double(16.) * 1e-3 * C_M,
np.double(18.) * 1e-3 * C_M,
np.double(1.) * 1e-3 * C_M,])
noms_mus0 = np.array(["CO", "C", "O", "H2O", "H"])
mass_fracs0 = np.array([(minit_kuips / minit) * 0.5,
(minit_kuips / minit) * 0.5 * 12. / 28.,
(minit_asts / minit) * 0.5 * 16. / 18.
+ (minit_kuips / minit) * 0.5 * 16. / 28.,
(minit_asts / minit) * 0.5,
(minit_asts / minit) * 0.5 * 2. / 18.
+ self.mH_init / minit])
self.a0 = (self.sol_ast_params["a0"] + self.sol_kuip_params["a0"]) / 2.
elif mus0 is None:
mus0 = np.array([1e-3 * C_M])
noms_mus0 = np.array(["H"])
mass_fracs0 = np.array([1.])
for (mu0, n_mu0) in zip(mus0, noms_mus0):
if n_mu0 not in nom_mu:
nom_mu.append(n_mu0)
mu.append(mu0)
self.mu: np.ndarray[np.float64] = np.array(mu, dtype=np.float64)
self.nom_mu: np.ndarray[str] = np.array(nom_mu, dtype=str)
self.diag: dict[str, list[np.ndarray]] = dict(time=[], total_momentum=[], total_mass=[],
phi=[], mdot=[], phi_L=[], Ldot=[])
self.flux_planets: list[np.ndarray] = [
] # mass fluxes at ip+self.dip and ip-self.dip for each planet for each time
self.flux_radmc: list[np.ndarray[np.float64]] = [] # Luminosity flux fo 3 differents wavelenghth integrated for the whole disc.
# Initialisation
self.initialization(mus0=mus0, mass_fracs0=mass_fracs0, nom_mus0=noms_mus0)
self.interp_f = None
else:
# Loading the data from a previously integrated solution saved at directory/filename
if ".npz" not in filename:
filename += ".npz"
if directory == "":
i: int = filename.find("/")
j: int = i
while j != -1:
i = j
j = filename.find("/", i + 1)
directory = filename[:i]
filename = filename[i + 1:]
self.name: str = filename[:-4]
self.directory: str = directory
dic: dict = dict(np.load(directory + "/" + filename))
# if np.any(["sol_ast_" in k for k in dic.keys()]):
# # Extraction of the asteroid belt solution
# self.sol_ast: SolrhoBelt = SolrhoBelt(consts_filename=filename, dossier=directory,
# dic=dic, prefix="sol_ast_")
# else:
# # The default asteroid solution to use if not save with the rest of the solution
# self.sol_ast: SolrhoBelt = SolrhoBelt("Simulation_rho_3.npz")
self.mu = np.array(dic['mu']) * C_M / dic["C_M"]
self.nom_mu = np.array(dic['nom_mu'])
self.r = np.array(dic['r']) * (C_L / dic["C_L"])
self.x = np.sqrt(self.r)
self.tps = list(dic["tps"] * C_t / dic["C_t"])
self.tps_mf = list(dic["tps_mf"] * C_t / dic["C_t"])
self.f = (list(dic["f"] * (C_M * (C_L ** np.float64(dic["pls_temp"])))
/ (dic["C_M"] * (dic["C_L"] ** np.float64(dic["pls_temp"])))))
self.mass_fractions = list(dic["mass_fractions"])
self.flux_planets = list(dic["flux_planets"] * (C_M / C_t) / (dic["C_M"] / dic["C_t"]))
self.idxa = np.array(dic["idxa"])
self.update_constantes(dic)
self.diag: dict[str, list[np.float64]] = {}
for k in dic.keys():
# Datas save at each time-steps to control the mass conservation
if "diag_" in k:
if len(dic[k]) > 0:
self.diag[k[5:]] = list(dic[k])
if type(self.diag[k[5:]][0]) is np.ndarray:
self.diag[k[5:]] = [list(liste) for liste in dic[k]]
if len(self.diag) == 0:
self.diag: dict[str, list[np.ndarray]] = dict(time=[], total_momentum=[], total_mass=[],
phi=[], mdot=[], phi_L=[], Ldot=[])
self.idx_mu_ast = np.array(dic["idx_mu_ast"], dtype=int)
self.idx_mu_kuip = np.array(dic["idx_mu_kuip"], dtype=int)
if "a_planets" in dic.keys():
self.a_planets = np.array(dic["a_planets"])
else:
self.a_planets = np.array([])
if "m_planets" in dic.keys():
self.m_planets = np.array(dic["m_planets"])
else:
self.m_planets = np.array([])
self.idxa = np.intersect1d(np.argwhere(self.r > dic["sol_ast_params_a0"]
- dic["sol_ast_params_delta_a"])[:, 0],
np.argwhere(self.r < dic["sol_ast_params_a0"]
+ dic["sol_ast_params_delta_a"])[:, 0])
self.idxk = np.intersect1d(np.argwhere(self.r > dic["sol_kuip_params_a0"]
- dic["sol_kuip_params_delta_a"])[:, 0],
np.argwhere(self.r < dic["sol_kuip_params_a0"]
+ dic["sol_kuip_params_delta_a"])[:, 0])
if "sol_ast_sig_dots" in dic.keys():
self.sol_ast = SolrhoBelt(kwargs_N=self.sol_ast_kwargs_N, sig_dots=dic["sol_ast_sig_dots"],
tps=dic["sol_ast_tps"], radius=self.r[self.idxa], **self.sol_ast_params)
else:
self.sol_ast = SolrhoBelt(kwargs_N=self.sol_ast_kwargs_N, radius=self.r[self.idxa],
**self.sol_ast_params)
if "sol_kuip_sig_dots" in dic.keys():
self.sol_kuip = SolrhoBelt(kwargs_N=self.sol_kuip_kwargs_N, sig_dots=dic["sol_kuip_sig_dots"],
tps=dic["sol_kuip_tps"], radius=self.r[self.idxk], **self.sol_kuip_params)
else:
self.sol_kuip = SolrhoBelt(kwargs_N=self.sol_kuip_kwargs_N, radius=self.r[self.idxk], **self.sol_kuip_params)
if "flux_radmc" in dic.keys():
self.flux_radmc = list(dic["flux_radmc"])
else:
self.flux_radmc: list[np.ndarray[
np.float64]] = [] # Luminosity flux fo 3 differents wavelenghth integrated for the whole disc.
# Since the equations of mass fractions and surface density are uncoupled, we use an interpolation
# function to get the surface density at each time (radii are the same for both)
if len(self.f) > 2:
tps, idx = np.unique(self.tps, return_index=True)
self.tps = list(tps)
self.f = list(np.array(self.f)[idx])
def linear_interp(t_new):
if isinstance(t_new, list | np.ndarray):
return np.array([linear_interp(t) for t in t_new])
i = np.argmin(np.abs(t_new - np.array(self.tps)))
if i == 0 or (i < len(self.tps) - 1 and t_new > self.tps[i]) :
t1 = self.tps[i]
t2 = self.tps[i + 1]
alpha = (t_new - t1) / (t2 - t1)
return (1 - alpha) * self.f[i] + alpha * self.f[i + 1]
else:
t1 = self.tps[i - 1]
t2 = self.tps[i]
alpha = (t_new - t1) / (t2 - t1)
return (1 - alpha) * self.f[i - 1] + alpha * self.f[i]
self.interp_f = linear_interp
# self.interp_f = CubicSpline(self.tps, self.f, extrapolate=True)
# sigmad0: np.float64 = (1e-3 * Mearth / Myr
# ) / (4.
# * Pi * self.a0 * self.a0 * (np.sqrt(1. + self.delta_r / self.a0)
# - np.sqrt(1. - self.delta_r / self.a0)))
# tps: np.ndarray = np.geomspace(yr, Gyr, 1000)
# sigmasd: np.ndarray = np.array([sigmad0 * ((self.r[self.idxa] / self.a0) ** (-3. / 2.)) for t in tps])
#
# self.interpsigdot = Akima1DInterpolator(tps - yr, 0.1 * sigmasd)
[docs]
def radial_speed(self, i: int | np.ndarray = ii_max) -> np.ndarray:
"""
The radial speed for the temporal index(s) i of tps/f
Parameters
----------
i : int | np.ndarray, optional, default: the radial speed is calculated for every index
The temporal index(s) from which calculate the radial speed
Returns
-------
np.ndarray
The radial speed
"""
if type(i) is int and i == ii.max:
return np.array([radial_speed(sol=self, ys=self.f[i], t=self.tps[i])
for i in range(len(self.tps))])
else:
return radial_speed(sol=self, ys=self.f[i], t=self.tps[i])
[docs]
def mass_flux(self, i: int | np.ndarray = ii_max) -> np.ndarray:
"""
The mass flux for the temporal index(s) i of tps/f
Parameters
----------
i : int | np.ndarray, optional, default: the radial speed is calculated for every index
The temporal index(s) from which calculate the radial speed
Returns
-------
np.ndarray
The mass flux
"""
if type(i) is int and i == ii.max:
return np.array([mass_flux(sol=self, ys=self.f[i], t=self.tps[i])
for i in range(len(self.tps))])
else:
return mass_flux(sol=self, ys=self.f[i], t=self.tps[i])
[docs]
def turbulent_speed(self, i: int | np.ndarray = ii_max) -> np.ndarray:
"""
The turbulent speed for the temporal index(s) i of mass_fractions
Parameters
----------
i : int | np.ndarray, optional, default:the turbulence speed is calculated for every index
The temporal index(s) from which calculate the radial speed
Returns
-------
np.ndarray
The turbulent speeds
"""
if type(i) is int and i == ii.max:
return np.array([turbulent_speed(sol=self,
y=np.array(self.mass_fractions)[i], t=self.tps[i])
for i in range(len(self.tps))])
else:
return turbulent_speed(sol=self,
y=np.array(self.mass_fractions)[i], t=self.tps[i])
[docs]
def constantes(self) -> dict:
"""
All the constants needed for the integration
Returns
-------
dict : dict
A dictionary with all the constantes :
- C_L : np.float64, 1/6371009 m -> au: length conversion factor
- C_M : np.float64, 1/(5.9722*1e24) kg -> Mearth : mass conversion factor
- C_t : np.float64, 1 / (3600 * 24 * 365.25 * 1e6) : s-> 1 Myr : time conversion factor
- Ts : np.float64, Surface temperature of the central star
- Ms : np.float64, Mass of the central star
- Rs : np.float64, Radius of the central star
- Ls : np.float64, Luminosity of the central star
- a_in : np.float64, Inner radius of the integration domain
- a0 : np.float64,
- a_out : np.float64, Outer radius of the integration domain
- T0 : np.float64, Temperature at a0
- mH_init : np.float64, The initial hydrogen mass
- AB_mdot_st_init : np.float64, The gas mass generation rate that generate the initial steady state
- AB_Mbelt : np.float64, Asteroid belt's size
- AB_Abelt : np.float64, Asteroids albedo
- AB_a0 : np.float64, Asteroid belt's central radius
- AB_delta_r : np.float64, Half of asteroid belt size
- AB_K : np.float64, Thermal diffusion coefficient in asteroids
- AB_Phi : np.float64, Porosity
- AB_rp : np.float64, Asteroid's pore radius
- AB_f_ice : np.float64, Ice to total mass ratio
- AB_rho_refr : np.float64, Refractory's materials density in asteroids
- AB_rho_ice : np.float64 Ice density
- AB_sublimation_model : str, {"none", "thermal_full", "constant_rate"} The model use to estimate
the gas generation rate for the asteroid belt :
- none : No gas is produced
- thermal_full : the gas is produced by sublimation : all sublimation, thermal diffusion and molecular
diffusion are take into acompte
- constant_rate : The gas is produced at a constant rate of const_mdot
- AB_const_mdot : np.float64, The gass mass generation rate for a `constant_rate` sublimation model
- AB_SD_amax : np.double, Maximum size of asteroids
- AB_SD_abig : np.double, Intermediate size
- AB_SD_amed : np.double, Medium size for size distribution
- AB_SD_amin : np.double, Minimal size
- AB_SD_qh : np.double, Power slope of the distribution between abig and amax
- AB_SD_qm : np.double, Power slope of the distribution between amed and abig
- AB_SD_ql : np.double, Power slope of the distribution between amin and amed
- AB_s_min: np.float64 = np.double(316) * C_L # Minimum size for collision model
- AB_sb: np.float64 = np.double(316) * C_L # Intermediate size for collision model
- AB_s_max: np.float64 = np.double(316) * C_L # Maximum size for collision model
- AB_e: np.float64 = np.double(0.075) # Eccentricity
- AB_i: np.float64 = AB_e / np.double(2.) # Inclination
- AB_As: np.float64 = np.double(5.) * C_E / C_M # Massique energy to disrupt the asteroids
- AB_bs: np.float64 = np.double(-0.1)
- AB_bg: np.float64 = np.double(0.5)
- AB_t0_diss: np.float64 = 10 * Gyr # Belt's lifetime before dissipation
- AB_t1_diss: np.float64 = 10 * Gyr # Dissipation timescale
- KB_mdot_st_init : np.float64, The gas mass generation rate that generate the initial steady state
- KB_Mbelt : np.float64, Kuiper belt's mass
- KB_Abelt : np.float64, KBO's albedo
- KB_a0 : np.float64, Kuiper belt central radius
- KB_delta_r : np.float64, Half of Kuiper belt size
- KB_K : np.float64, Thermal diffusion coefficient in KBO
- KB_Phi : np.float64, KBO's porosity
- KB_rp : np.float64, KBO's pore's radius
- KB_f_ice : np.float64, Initial ice to total mass ratio in KBO
- KB_rho_refr : np.float64, Refractory's materials density in KBO
- KB_rho_ice : np.float64, Ice density
- KB_sublimation_model : str, {"none", "thermal_full", "constant_rate"} The model use to estimate the
gas generation rate for the Kuiper belt:
- none : No gas is produced
- thermal_full : the gas is produced by sublimation : all sublimation, thermal diffusion and molecular
diffusion are take into acompte
- constant_rate : The gas is produced at a constant rate of const_mdot
- KB_const_mdot : np.float64, The gass mass generation rate for a `constant_rate` sublimation model
- KB_SD_amax : np.double, Maximum size of KBO
- KB_SD_abig : np.double, Intermediate size
- KB_SD_amed : np.double, Medium size
- KB_SD_amin : np.double, Minimal size
- KB_SD_qh : np.double, Power slope of the distribution between abig and amax
- KB_SD_qm : np.double, Power slope of the distribution between amed and abig
- KB_SD_ql : np.double, Power slope of the distribution between amin and amed
- KB_s_min: np.float64 = np.double(316) * C_L # Minimum size for collision model
- KB_sb: np.float64 = np.double(316) * C_L # Intermediate size for collision model
- KB_s_max: np.float64 = np.double(316) * C_L # Maximum size for collision model
- KB_e: np.float64 = np.double(0.075) # Eccentricity
- KB_i: np.float64 = KB_e / np.double(2.) # Inclination
- KB_As: np.float64 = np.double(5.) * C_E / C_M # Massique energy to disrupt the asteroids
- KB_bs: np.float64 = np.double(-0.1)
- KB_bg: np.float64 = np.double(0.5)
- KB_t0_diss: np.float64 = 10 * Gyr # Belt's lifetime before dissipation
- KB_t1_diss: np.float64 = 10 * Gyr # Dissipation timescale
- a_planet : list[np.double], planet's semi-major axis
- m_planet : list[np.double], planet's masses
- f_accr : np.float64, Hydrodynamic accretion efficiency
- t0 : The dissipation time of the protoplanetary disc
- IT_s : int, number of spatial steps
- dip : int, half of the number of sink cells for a given planets to model the accretion
- rtol : np.float64, Relative tolerance for temporal integration
- atol : np.float64, Absolute tolerance for temporal integration (by default too small hence to
relative tolerance will always be the limiting factor)
- method_root : str, 'LSODA','BDF','DOP853','RK45' : Temporal integration method_root
"""
return {"C_L": self.C_L, "C_M": self.C_M, "C_t": self.C_t,
"Ms": self.Ms, "Ts": self.Ts, "Rs": self.Rs, "Ls": self.Ls,
"a_in": self.a_in, "a_out": self.a_out, "mH_init": self.mH_init, "a0": self.a0,
"T0": self.T0, "pls_temp": self.pls_temp,
"alpha": self.alpha, "gamma": self.gamma,
"AB_mdot_st_init" : self.sol_ast_params["mdot_st_init"],
"AB_Mbelt": self.sol_ast_params["Mbelt"], "AB_Abelt": self.sol_ast_params["A_ast"],
"AB_a0": self.sol_ast_params["a0"],
"AB_delta_r": self.sol_ast_params["delta_a"], "AB_K": self.sol_ast_params["K"],
"AB_Phi": self.sol_ast_params["Phi"], "AB_rp": self.sol_ast_params["rp"],
"AB_f_ice": self.sol_ast_params["f_ice"], "AB_rho_refr": self.sol_ast_params["rho_refr"],
"AB_rho_ice": self.sol_ast_params["rho_ice"], "AB_m_init": self.sol_ast.m_init,
"AB_sublimation_model": self.sol_ast_params["sublimation_model"],
"AB_const_mdot": self.sol_ast_params["const_mdot"],
"AB_SD_amax": self.sol_ast_kwargs_N["amax"],
"AB_SD_abig": self.sol_ast_kwargs_N["abig"], "AB_SD_amed": self.sol_ast_kwargs_N["amed"],
"AB_SD_amin": self.sol_ast_kwargs_N["amin"], "AB_SD_qh": self.sol_ast_kwargs_N["qh"],
"AB_SD_qm": self.sol_ast_kwargs_N["qm"], "AB_SD_ql": self.sol_ast_kwargs_N["ql"],
"AB_sb": self.sol_ast_params["sb"],
"AB_e": self.sol_ast_params["e"], "AB_i": self.sol_ast_params["i"],
"AB_As": self.sol_ast_params["As"], "AB_bs": self.sol_ast_params["bs"],
"AB_bg": self.sol_ast_params["bg"],
"AB_t0_diss": self.sol_ast_params["t0_diss"], "AB_t1_diss": self.sol_ast_params["t1_diss"],
"KB_mdot_st_init" : self.sol_kuip_params["mdot_st_init"],
"KB_Mbelt": self.sol_kuip_params["Mbelt"], "KB_Abelt": self.sol_kuip_params["A_ast"],
"KB_a0": self.sol_kuip_params["a0"],
"KB_delta_r": self.sol_kuip_params["delta_a"], "KB_K": self.sol_kuip_params["K"],
"KB_Phi": self.sol_kuip_params["Phi"], "KB_rp": self.sol_kuip_params["rp"],
"KB_f_ice": self.sol_kuip_params["f_ice"], "KB_rho_refr": self.sol_kuip_params["rho_refr"],
"KB_rho_ice": self.sol_kuip_params["rho_ice"], "KB_m_init": self.sol_kuip.m_init,
"KB_sublimation_model": self.sol_kuip_params["sublimation_model"],
"KB_SD_amax": self.sol_kuip_kwargs_N["amax"],
"KB_SD_abig": self.sol_kuip_kwargs_N["abig"], "KB_SD_amed": self.sol_kuip_kwargs_N["amed"],
"KB_SD_amin": self.sol_kuip_kwargs_N["amin"], "KB_SD_qh": self.sol_kuip_kwargs_N["qh"],
"KB_SD_qm": self.sol_kuip_kwargs_N["qm"], "KB_SD_ql": self.sol_kuip_kwargs_N["ql"],
"KB_sb": self.sol_kuip_params["sb"],
"KB_e": self.sol_kuip_params["e"], "KB_i": self.sol_kuip_params["i"],
"KB_As": self.sol_kuip_params["As"], "KB_bs": self.sol_kuip_params["bs"],
"KB_bg": self.sol_kuip_params["bg"],
"KB_t0_diss": self.sol_kuip_params["t0_diss"], "KB_t1_diss": self.sol_kuip_params["t1_diss"],
"a_planets": self.a_planets, "m_planets": self.m_planets, "f_accr": self.f_accr,
"t0": self.t0,
"IT_s": self.IT_s, "dip": self.dip, "rtol": self.rtol, "atol": self.atol, "method_root": self.method,
}
[docs]
def update_constantes(self, dic: dict) -> None:
"""
Update all the constants needed for the integration
Parameters
----------
dic: dict
A dictionary with all the constantes :
- C_L : np.float64, 1/6371009 m -> au: length conversion factor
- C_M : np.float64, 1/(5.9722*1e24) kg -> Mearth : mass conversion factor
- C_t : np.float64, 1 / (3600 * 24 * 365.25 * 1e6) : s-> 1 Myr : time conversion factor
- Ts : np.float64, Surface temperature of the central star
- Ms : np.float64, Mass of the central star
- Rs : np.float64, Radius of the central star
- Ls : np.float64, Luminosity of the central star
- a_in : np.float64, Inner radius of the integration domain
- a0 : np.float64,
- a_out : np.float64, Outer radius of the integration domain
- T0 : np.float64, Temperature at a0
- mH_init : np.float64, The initial hydrogen mass
- AB_mdot_st_init : np.float64, The gas mass generation rate that generate the initial steady state
- AB_Mbelt : np.float64, Asteroid belt's size
- AB_Abelt : np.float64, Asteroids albedo
- AB_a0 : np.float64, Asteroid belt's central radius
- AB_delta_r : np.float64, Half of asteroid belt size
- AB_K : np.float64, Thermal diffusion coefficient in asteroids
- AB_Phi : np.float64, Porosity
- AB_rp : np.float64, Asteroid's pore radius
- AB_f_ice : np.float64, Ice to total mass ratio
- AB_rho_refr : np.float64, Refractory's materials density in asteroids
- AB_rho_ice : np.float64 Ice density
- AB_sublimation_model : str, {"none", "thermal_full", "constant_rate"} The model use to estimate
the gas generation rate for the asteroid belt :
- none : No gas is produced
- thermal_full : the gas is produced by sublimation : all sublimation, thermal diffusion and molecular
diffusion are take into acompte
- constant_rate : The gas is produced at a constant rate of const_mdot
- AB_const_mdot : np.float64, The gass mass generation rate for a `constant_rate` sublimation model
- AB_SD_amax : np.double, Maximum size of asteroids
- AB_SD_abig : np.double, Intermediate size
- AB_SD_amed : np.double, Medium size for size distribution
- AB_SD_amin : np.double, Minimal size
- AB_SD_qh : np.double, Power slope of the distribution between abig and amax
- AB_SD_qm : np.double, Power slope of the distribution between amed and abig
- AB_SD_ql : np.double, Power slope of the distribution between amin and amed
- AB_s_min: np.float64 = np.double(316) * C_L # Minimum size for collision model
- AB_sb: np.float64 = np.double(316) * C_L # Intermediate size for collision model
- AB_s_max: np.float64 = np.double(316) * C_L # Maximum size for collision model
- AB_e: np.float64 = np.double(0.075) # Eccentricity
- AB_i: np.float64 = AB_e / np.double(2.) # Inclination
- AB_As: np.float64 = np.double(5.) * C_E / C_M # Massique energy to disrupt the asteroids
- AB_bs: np.float64 = np.double(-0.1)
- AB_bg: np.float64 = np.double(0.5)
- AB_t0_diss: np.float64 = 10 * Gyr # Belt's lifetime before dissipation
- AB_t1_diss: np.float64 = 10 * Gyr # Dissipation timescale
- KB_mdot_st_init : np.float64, The gas mass generation rate that generate the initial steady state
- KB_Mbelt : np.float64, Kuiper belt's mass
- KB_Abelt : np.float64, KBO's albedo
- KB_a0 : np.float64, Kuiper belt central radius
- KB_delta_r : np.float64, Half of Kuiper belt size
- KB_K : np.float64, Thermal diffusion coefficient in KBO
- KB_Phi : np.float64, KBO's porosity
- KB_rp : np.float64, KBO's pore's radius
- KB_f_ice : np.float64, Initial ice to total mass ratio in KBO
- KB_rho_refr : np.float64, Refractory's materials density in KBO
- KB_rho_ice : np.float64, Ice density
- KB_sublimation_model : str, {"none", "thermal_full", "constant_rate"} The model use to estimate the
gas generation rate for the Kuiper belt:
- none : No gas is produced
- thermal_full : the gas is produced by sublimation : all sublimation, thermal diffusion and molecular
diffusion are take into acompte
- constant_rate : The gas is produced at a constant rate of const_mdot
- KB_const_mdot : np.float64, The gass mass generation rate for a `constant_rate` sublimation model
- KB_SD_amax : np.double, Maximum size of KBO
- KB_SD_abig : np.double, Intermediate size
- KB_SD_amed : np.double, Medium size
- KB_SD_amin : np.double, Minimal size
- KB_SD_qh : np.double, Power slope of the distribution between abig and amax
- KB_SD_qm : np.double, Power slope of the distribution between amed and abig
- KB_SD_ql : np.double, Power slope of the distribution between amin and amed
- KB_s_min: np.float64 = np.double(316) * C_L # Minimum size for collision model
- KB_sb: np.float64 = np.double(316) * C_L # Intermediate size for collision model
- KB_s_max: np.float64 = np.double(316) * C_L # Maximum size for collision model
- KB_e: np.float64 = np.double(0.075) # Eccentricity
- KB_i: np.float64 = KB_e / np.double(2.) # Inclination
- KB_As: np.float64 = np.double(5.) * C_E / C_M # Massique energy to disrupt the asteroids
- KB_bs: np.float64 = np.double(-0.1)
- KB_bg: np.float64 = np.double(0.5)
- KB_t0_diss: np.float64 = 10 * Gyr # Belt's lifetime before dissipation
- KB_t1_diss: np.float64 = 10 * Gyr # Dissipation timescale
- a_planet : list[np.double], planet's semi-major axis
- m_planet : list[np.double], planet's masses
- f_accr : np.float64, Hydrodynamic accretion efficiency
- t0 : The dissipation time of the protoplanetary disc
- IT_s : int, number of spatial steps
- dip : int, half of the number of sink cells for a given planets to model the accretion
- rtol : np.float64, Relative tolerance for temporal integration
- atol : np.float64, Absolute tolerance for temporal integration (by default too small hence to
relative tolerance will always be the limiting factor)
- method_root : str, 'LSODA','BDF','DOP853','RK45' : Temporal integration method_root
Returns
-------
None
"""
# if "directory" in dic.keys():
# self.directory = str(dic["directory"])
# else:
# self.directory = str(dic["results"])
self.method = str(dic["method_root"])
self.rtol = np.double(dic["rtol"])
self.atol = np.double(dic["atol"])
self.IT_s = int(dic["IT_s"])
self.dip = int(dic["dip"])
self.C_L = np.double(dic["C_L"])
self.C_M = np.double(dic["C_M"])
self.C_t = np.double(dic["C_t"])
self.C_E = np.double(dic["C_L"] * dic["C_t"])
self.C_F = np.double(self.C_M * self.C_L / (self.C_t * self.C_t)) # N-> : force
self.C_rho = np.double(self.C_M / (self.C_L * self.C_L * self.C_L)) # kg/M**6 -> : density
self.Ms = np.double(dic["Ms"] * C_M / self.C_M)
self.Ts = np.double(dic["Ts"])
self.Rs = np.double(dic["Rs"] * C_L / self.C_L)
self.Ls = np.double(dic["Ls"] * C_L / self.C_L)
self.T0 = np.double(dic["T0"])
self.pls_temp = np.double(dic["pls_temp"])
if "mH_init" in dic.keys():
self.mH_init = np.double(dic["mH_init"])
else:
self.mH_init = zero
self.a0 = np.double(dic["a0"])
self.a_in = np.double(dic["a_in"]) * C_L / self.C_L
self.a_out = np.double(dic["a_out"]) * C_L / self.C_L
self.alpha = np.double(dic["alpha"]) # viscous parameter
self.gamma = np.double(dic["gamma"])
self.t0 = np.double(dic["t0"]) * C_t / self.C_t
m_init_asts = zero
if "AB_m_init" in dic.keys():
m_init_asts = dic["AB_m_init"] * C_M / self.C_M
self.sol_ast_params = dict(sublimation_model=str(dic["AB_sublimation_model"]),
Mbelt=np.double(dic["AB_Mbelt"]) * C_M / self.C_M,
Ms=np.double(dic["Ms"]) * C_M / self.C_M,
a0=np.double(dic["AB_a0"]) * C_L / self.C_L,
delta_a=np.double(dic["AB_delta_r"]) * C_L / self.C_L,
tinit=self.t0, tmax=101 * Myr,
rho_refr=np.double(dic["AB_rho_refr"]) * C_rho / self.C_rho,
rho_ice=np.double(dic["AB_rho_ice"]) * C_rho / self.C_rho,
m_init=m_init_asts,
f_ice=np.double(dic["AB_f_ice"]), A_ast=np.double(dic["AB_Abelt"]),
K=np.double(dic["AB_K"]) * (C_L * C_L / C_t) / (self.C_L * self.C_L / self.C_t),
Its=90, Phi=np.double(dic["AB_Phi"]),
dmin=np.double(dic["AB_SD_amin"]) * C_L / self.C_L,
sb=np.double(dic["AB_sb"]) * C_L / self.C_L,
dmax=np.double(dic["AB_SD_amax"]) * C_L / self.C_L,
e=np.double(dic["AB_e"]), i=np.double(dic["AB_i"]),
As=np.double(dic["AB_As"]) * (C_E / C_M) / (self.C_E / self.C_M),
bs=np.double(dic["AB_bs"]), bg=np.double(dic["AB_bg"]),
t0_diss=np.double(dic["AB_t0_diss"]) * C_t / self.C_t,
t1_diss=np.double(dic["AB_t1_diss"]) * C_t / self.C_t, Its_p=500,
rp=np.double(dic["AB_rp"]) * C_t / self.C_t, Itt=5000)
if "AB_const_mdot" in dic.keys():
self.sol_ast_params["const_mdot"] = dic["AB_const_mdot"]
else:
self.sol_ast_params["const_mdot"] = m_dot
if "AB_mdot_st_init" in dic.keys():
self.sol_ast_params["mdot_st_init"] = dic["AB_mdot_st_init"]
else:
self.sol_ast_params["mdot_st_init"] = np.double(1e-10) * Mearth / Myr
f_rho_ice: np.float64 = (dic["AB_rho_ice"] - dic["AB_f_ice"] * dic["AB_rho_refr"]
) / (dic["AB_rho_ice"] - dic["AB_rho_refr"])
self.sol_ast_kwargs_N = dict(amax=np.double(dic["AB_SD_amax"]) * C_t / self.C_t,
abig=np.double(dic["AB_SD_abig"]) * C_t / self.C_t,
amed=np.double(dic["AB_SD_amed"]) * C_t / self.C_t,
amin=np.double(dic["AB_SD_amin"]) * C_t / self.C_t,
qh=np.double(dic["AB_SD_qh"]), qm=np.double(dic["AB_SD_qm"]),
ql=np.double(dic["AB_SD_ql"]),
rho=(np.double(dic["AB_rho_ice"]) / dic["AB_f_ice"]
) * C_rho / self.C_rho)
m_init_kuip = zero
if "KB_m_init" in dic.keys():
m_init_kuip = dic["KB_m_init"] * C_M / self.C_M
self.sol_kuip_params = dict(sublimation_model=str(dic["KB_sublimation_model"]),
Mbelt=np.double(dic["KB_Mbelt"]) * C_M / self.C_M,
Ms=np.double(dic["Ms"]) * C_M / self.C_M,
a0=np.double(dic["KB_a0"]) * C_L / self.C_L,
delta_a=np.double(dic["KB_delta_r"]) * C_L / self.C_L,
tinit=self.t0, tmax=101 * Myr,
rho_refr=np.double(dic["KB_rho_refr"]) * C_rho / self.C_rho,
rho_ice=np.double(dic["KB_rho_ice"]) * C_rho / self.C_rho,
m_init=m_init_kuip,
f_ice=np.double(dic["KB_f_ice"]), A_ast=np.double(dic["KB_Abelt"]),
K=np.double(dic["KB_K"]) * (C_L * C_L / C_t) / (self.C_L * self.C_L / self.C_t),
Its=60, Phi=np.double(dic["KB_Phi"]),
dmin=np.double(dic["KB_SD_amin"]) * C_L / self.C_L,
sb=np.double(dic["KB_sb"]) * C_L / self.C_L,
dmax=np.double(dic["KB_SD_amax"]) * C_L / self.C_L,
e=np.double(dic["KB_e"]), i=np.double(dic["KB_i"]),
As=np.double(dic["KB_As"]) * (C_E / C_M) / (self.C_E / self.C_M),
bs=np.double(dic["KB_bs"]), bg=np.double(dic["KB_bg"]),
t0_diss=np.double(dic["KB_t0_diss"]) * C_t / self.C_t,
t1_diss=np.double(dic["KB_t1_diss"]) * C_t / self.C_t, Its_p=100,
rp=np.double(dic["KB_rp"]) * C_t / self.C_t, Itt=500000)
if "KB_const_mdot" in dic.keys():
self.sol_kuip_params["const_mdot"] = dic["KB_const_mdot"]
else:
self.sol_kuip_params["const_mdot"] = m_dot
if "KB_mdot_st_init" in dic.keys():
self.sol_kuip_params["mdot_st_init"] = dic["KB_mdot_st_init"]
else:
self.sol_kuip_params["mdot_st_init"] = np.double(1e-10) * Mearth / Myr
self.sol_kuip_kwargs_N = dict(amax=dic["KB_SD_amax"] * C_L / self.C_L,
abig=dic["KB_SD_abig"] * C_L / self.C_L,
amed=dic["KB_SD_amed"] * C_L / self.C_L,
amin=dic["KB_SD_amin"] * C_L / self.C_L,
qh=dic["KB_SD_qh"], qm=dic["KB_SD_qm"], ql=dic["KB_SD_ql"],
rho=(dic["KB_rho_refr"] / dic["KB_f_ice"]) * C_rho / self.C_rho)
self.a_planets = list(np.array(dic["a_planets"], dtype="double") * C_L / self.C_L)
self.m_planets = list(np.array(dic["m_planets"], dtype="double") * C_L / self.C_L)
self.f_accr = np.double(dic["f_accr"])
self.C_L = np.double(C_L)
self.C_M = np.double(C_M)
self.C_t = np.double(C_t)
self.C_F = np.double(C_F)
[docs]
def initialization(self, mus0: np.ndarray, mass_fracs0: np.ndarray, nom_mus0) -> None:
"""
Initialize the solution with a steady state of the first element
generated by a 10**-10 Mearth/Myr gas generation rate
The radii are linearly spaced between a_min and a_max with IT_s points
Returns
-------
None
"""
self.x = np.linspace(np.sqrt(self.a_in), np.sqrt(self.a_out), self.IT_s)
self.r = self.x * self.x
f0: np.ndarray = (steady_state(sol=self, mdot=self.sol_ast_params["mdot_st_init"], a0=self.sol_ast_params["a0"])
+ steady_state(sol=self, mdot=self.sol_kuip_params["mdot_st_init"], a0=self.sol_kuip_params["a0"])
+ steady_state(sol=self, mtot=self.mH_init, a0=self.a0))
xs0: np.ndarray = np.zeros((len(self.mu), len(self.r)))
if len(mass_fracs0.shape) == 1 and len(mass_fracs0) == len(mus0):
# Identical initial mass fraction for each radius
if abs(mass_fracs0.sum() - 1) > 1e-5:
raise UserWarning("Sol.initialization : The sum of initial mass fractions must be equal to one, not ",
mass_fracs0.sum())
for i in range(len(mus0)):
xs0[np.argwhere(self.nom_mu == nom_mus0[i])[0, 0]] = mass_fracs0[i]
elif len(mass_fracs0.shape) == 1:
raise UserWarning("Sol.initialization : The initial mass fractions array must have the same size of "
"the initials elements arrays mus0: len(mus0)=", len(mus0), "but len(mass_fracs0)=",
len(mass_fracs0))
elif len(mass_fracs0.shape) == 2 and mass_fracs0.shape[0] == len(mus0) and mass_fracs0.shape[1] == len(self.r):
# Identical initial mass fraction for each radius
if np.sum(abs(mass_fracs0.sum(axis=1) - 1)) > 1e-5 / len(self.r):
raise UserWarning("Sol.initialization : The sum of initial mass fractions must be equal to one, not ",
mass_fracs0.sum())
for i in range(len(mus0)):
xs0[np.argwhere(self.mu == mus0[i])[0, 0]] = mass_fracs0[i]
elif len(mass_fracs0.shape) == 2:
raise UserWarning("Sol.initialization : The initial mass fractions array must have the same size of "
"the initials elements arrays mus0 on its first dimension and the size of self.r "
"(self.IT_s) on its second: len(mus0)=", len(mus0), "len(self.r)", len(self.r),
"and mass_fracs0.shape", mass_fracs0.shape)
else:
raise UserWarning("Sol.initialization : The dimension of mass_fracs0 must be 1 or 2, not ",
len(mass_fracs0.shape))
self.idxa = np.intersect1d(np.argwhere(self.r > self.sol_ast_params["a0"]
- self.sol_ast_params["delta_a"])[:, 0],
np.argwhere(self.r < self.sol_ast_params["a0"]
+ self.sol_ast_params["delta_a"])[:, 0])
self.idxk = np.intersect1d(np.argwhere(self.r > self.sol_kuip_params["a0"]
- self.sol_kuip_params["delta_a"])[:, 0],
np.argwhere(self.r < self.sol_kuip_params["a0"]
+ self.sol_kuip_params["delta_a"])[:, 0])
self.f = [f0]
self.mass_fractions = [xs0]
self.tps = [zero]
self.tps_mf = [0.2 * Myr]
# self.tps_mf = [10 * kyr]
self.sol_ast_params["radius"] = self.r[self.idxa]
self.sol_ast_params["Its"] = len(self.idxa)
self.sol_ast = SolrhoBelt(kwargs_N=self.sol_ast_kwargs_N, **self.sol_ast_params)
self.sol_ast_params.update(self.sol_ast.const())
if self.sol_ast.sublimation_model == "thermal_full":
sigma_init = (self.sol_ast.m_init / (4. * Pi * self.sol_ast.a0 * self.sol_ast.a0
* (np.sqrt(1. + self.sol_ast.delta_a / self.sol_ast.a0)
- np.sqrt(1. - self.sol_ast.delta_a / self.sol_ast.a0)))
* ((self.sol_ast.r / self.sol_ast.a0) ** (-3. / 2.)))
# self.mass_fractions[-1][:, self.idxa] = np.maximum(self.mass_fractions[-1][:, self.idxa]
# - sigma_init / (self.f[0][self.idxa] * (len(self.mu) - 1)),
# 0.)
# self.mass_fractions[-1][self.idx_mu_ast[0], self.idxa] = np.minimum(self.mass_fractions[-1][self.idx_mu_ast[0], self.idxa]
# + sigma_init / (self.f[0][self.idxa]
# * (len(self.mu) - 1))
# + sigma_init / (self.f[0][self.idxa]), 1)
self.f[0][self.idxa] += sigma_init / (self.r[self.idxa] ** (-2 - self.pls_temp))
self.sol_kuip_params["radius"] = self.r[self.idxk]
self.sol_kuip_params["Its"] = len(self.idxk)
self.sol_kuip = SolrhoBelt(kwargs_N=self.sol_kuip_kwargs_N, **self.sol_kuip_params)
self.sol_kuip_params.update(self.sol_kuip.const())
self.save()
[docs]
def sigma(self, i: int = ii_max) -> np.float64 | np.ndarray[np.float64]:
"""
Return surface density
Parameters
----------
i : int, optional
The index at wich calculate sigma
Returns
-------
np.float64 | np.ndarray
"""
if type(i) is int and i == ii.max:
return np.array(self.f) * (self.r ** (-2 - self.pls_temp))
else:
return np.array(self.f[i]) * (self.r ** (-2 - self.pls_temp))
[docs]
def interp_sigma(self, t: np.float64) -> np.float64 | np.ndarray[np.float64]:
"""
Return surface density
Parameters
----------
i : int, optional
The index at wich calculate sigma
Returns
-------
np.float64 | np.ndarray
"""
return self.interp_f(t) * (self.r ** (-2 - self.pls_temp))
[docs]
def y(self, i: int) -> np.ndarray:
"""
Return the masses fractions flatten into a one dimensional array for a given time (argument
of the time integration system)
Parameters
----------
i : int
The index to which return the mass fraction
Returns
-------
"""""
return self.mass_fractions[i].flatten()
[docs]
def split_y(self, y: np.ndarray) -> np.ndarray:
"""
Split y the argument of the time integration system in n lists, n the number of chemical elements
Parameters
----------
y : np.ndarray
The array to be reshaped
Returns
-------
np.ndarray
masses_fractions(y)
"""
if len(y) != len(self.r) * (len(self.mu)):
raise UserWarning("split_y : the dimension of y doesn't mach with the number of element and the size of the"
" radius list for this Sol", len(y), len(self.r) * len(self.mu))
return y.reshape((len(self.mu), len(self.r)))
[docs]
def get_final_surface_densitys(self, y: np.ndarray, tps: np.ndarray, step: int = 1) -> None:
"""
Save into this Sol the result of a temporal integration contain in ys. This save a fraction of
y in self.f (fraction given by step) and calculate some test variable such as the flux at the planets levels or
the variables needed to calculate the mass conservation (inner and outer flux, total mass produced...)
Parameters
----------
y : np.ndarray
result of solve_ivp = sigma * (r ** (2 + pls_temp))
tps : np.ndarray
times associated with y
step : int, optional, default=
Only the tps[::step] will be saved
Returns
-------
None
"""
# for i in range(step, len(y), step):
ips = []
for r_planet in self.a_planets:
ips.append(np.argwhere(self.r >= r_planet)[0, 0])
ips = np.array(ips)
for i in range(0, len(y)):
y[i] = boundary_conditions(y[i], self)
self.diag["time"].append(tps[i])
self.diag["total_mass"].append(trapezoid(2. * Pi * self.r[2:-2] * y[i][2:-2] * (self.r[2:-2] ** (-2. - self.pls_temp)),
x=self.r[2:-2]))
self.diag["total_momentum"].append(trapezoid(4. * Pi * np.sqrt(G * self.Ms) * (self.x[2:-2] ** (- 2. * self.pls_temp)) * y[i][2:-2],
x=self.x[2:-2]))
phi = mass_flux(self, y[i], tps[i])[2:-2]
self.diag["Ldot"].append(trapezoid(4 * Pi * np.sqrt(G * self.Ms) * self.r[2:-2] * self.r[2:-2]
* sigma_dot(t=tps[i], ys=y[i], phi=phi, sol=self)[2:-2], x=self.x[2:-2]))
self.diag["mdot"].append(trapezoid(2 * Pi * self.r[2:-2]
* sigma_dot(t=tps[i], ys=y[i], phi=phi, sol=self)[2:-2], x=self.r[2:-2]))
self.diag['phi'].append(- phi[0] + phi[-1])
nu0: np.float64 = (Rgp * T0(tps[i], self.a0, self.Ms) * self.alpha * (self.a0 ** (3 / 2))
/ (self.mu.mean() * np.sqrt(G * self.Ms)))
D: np.float64 = (3 * nu0 / (4. * (self.a0 ** (3 / 2 + self.pls_temp))))
self.diag['phi_L'].append(- 4. * Pi * D * np.sqrt(G * self.Ms) * (
self.x[0] ** (1 + pls_temp) * (y[i][1] - y[i][0]) / (self.x[1] - self.x[0])
- self.x[-1] ** (1 + pls_temp) * (y[i][-1] - y[i][-2]) / (self.x[-1] - self.x[-2]))
)
if len(ips) > 0:
self.flux_planets.append(
np.array([phi[ips + dip], phi[ips - dip]]))
if i % step == 0:
self.tps.append(tps[i])
self.f.append(y[i])
[docs]
def get_final_mass_fractions(self, y: np.ndarray, tps: np.float64, step: int = 1, x: np.float64 = None
) -> np.ndarray:
"""
Transform the result get by solve_ivp for the integration of the mass fractions
to the format of self.sigmas
Parameters
----------
y : np.ndarray
result of solve_ivp
tps : np.ndarray
times associated with y
step : int, optional, default=
Only the tps[::step] will be saved
x : np.ndarray, optional, default: self.x
The radial coordinate
Returns
-------
None
"""
if x is None:
x = self.x
for i in range(0, len(y), step):
self.tps_mf.append(tps[i])
self.mass_fractions.append(boundary_conditions_mf(split_y(y[i], len(x), len(self.mu)), self))
return boundary_conditions_mf(split_y(y[-1], len(x), len(self.mu)), self)
[docs]
def total_mass(self) -> np.ndarray:
"""
Calculate the total mass as a function of time
Returns
-------
"""
return trapezoid(2 * Pi * self.r * self.f * (self.r ** (-2 - self.pls_temp)), x=self.r, axis=1)
[docs]
def mass_conservation(self) -> np.ndarray:
"""
Calculate the mass conservation in the whole system
Returns
-------
np.ndarray
An array fo size self.diag["time"], equal to zero only if the mass is conserved
"""
return abs(np.array(self.diag["total_mass"])
- (self.diag["total_mass"][0]
+ cumulative_trapezoid(-np.array(self.diag["phi"]) + np.array(self.diag["mdot"]),
self.diag["time"], initial=zero))
) / np.array(self.diag["total_mass"])
[docs]
def momentum_conservation(self) -> np.ndarray:
"""
Calculate the momentum conservation in the whole system
Returns
-------
np.ndarray
An array fo size self.diag["time"], equal to zero only if the mass is conserved
"""
return abs(np.array(self.diag["total_momentum"])
- (self.diag["total_momentum"][0]
+ cumulative_trapezoid(-np.array(self.diag["phi_L"]) + np.array(self.diag["Ldot"]),
self.diag["time"], initial=zero))
) / np.array(self.diag["total_momentum"])
[docs]
def time_integration(self, tf: np.float64, dtmax: np.float64 = inf, build_graphs: bool = False,
integrate_only_surface_density: bool = False) -> None:
"""
Time integration of surface density and mass fraction until tf
The surface density is integrated first
Steps of 0.3 Myr are made to save the result
Parameters
----------
tf : np.float64
The final integration time
dtmax : np.float64, optional, default=the maximum timestep to keep a stable spatial integration of surface
density
The maximum timestep
build_graphs : bool, optional, default=True
To build and save surface density and accretions Graphiques
integrate_only_surface_density : bool, optional, default=False
To integrate only the evolution of surfaces densitys
Returns
-------
None
"""
if dtmax == inf:
nu0: np.float64 = (Rgp * T0(self.t0, self.a0, self.Ms) * self.alpha * (self.a0 ** (3 / 2))
/ (self.mu.mean() * np.sqrt(G * self.Ms)))
dtmax: np.float64 = min((self.r[1] - self.r[0]) * (self.r[1] - self.r[0])
/ (nu0 * ((self.r / self.a0) ** (3/2 + self.pls_temp))),
60. * yr)
if self.tps[-1] < tf:
print("\n Integration surfaces density \n")
dtmax = min(dtmax, (tf - self.tps[-1]) / 20)
sol_ivp = solve_ivp(system_sigma_tot, t_span=[self.tps[-1], tf],
y0=self.f[-1],
max_step=dtmax, method=self.method, args=[self],
rtol=self.rtol, atol=self.atol, first_step = dtmax * 1e-6)
try:
if len(sol_ivp.t) > 2:
step: int = max(len(sol_ivp.t) // 150, 1)
self.get_final_surface_densitys(sol_ivp.y.T[0:], tps=sol_ivp.t[0:], step=step)
print("t=", sol_ivp.t[-1] / yr, self.tps[-1] / yr)
self.save()
else:
print(sol_ivp.t)
self.save()
raise RuntimeError("Convergence issue")
except KeyboardInterrupt:
self.save()
raise UserWarning("Keyboard interrupt")
tps, idx = np.unique(self.tps, self.f)
self.tps = list(tps)
self.f = list(np.array(self.f)[idx])
def linear_interp(t_new):
if isinstance(t_new, list | np.ndarray):
return np.array([linear_interp(t) for t in t_new])
i = np.argmin(np.abs(t_new - np.array(self.tps)))
if i == 0 or (i < len(self.tps) - 1 and t_new > self.tps[i]):
t1 = self.tps[i]
t2 = self.tps[i + 1]
alpha = (t_new - t1) / (t2 - t1)
return (1 - alpha) * self.f[i] + alpha * self.f[i + 1]
else:
t1 = self.tps[i - 1]
t2 = self.tps[i]
alpha = (t_new - t1) / (t2 - t1)
return (1 - alpha) * self.f[i - 1] + alpha * self.f[i]
self.interp_f = linear_interp
# self.interp_f = CubicSpline(self.tps, self.f)
if self.tps_mf[-1] < tf and len(self.mu) > 1 and not integrate_only_surface_density:
print("\n Integration masses fractions \n")
dtmax_mf = min(dtmax, yr)
# dtmax_mf = dtmax
tfi: np.float64 = min(tf, self.tps_mf[-1] + dtmax_mf * 1e5)
results = []
tps = []
ti: np.float64 = self.tps_mf[-1]
compt: int = 0
while tfi <= tf:
print("dtmax=", dtmax_mf / yr, "yr", "tfi=", tfi)
if len(results) == 0:
sol_ivp = solve_ivp(system_mass_fraction, t_span=[ti, tfi],
y0=self.mass_fractions[-1].flatten(),
max_step=dtmax_mf, method=self.method, args=[self],
rtol=self.rtol, atol=self.atol)
else:
sol_ivp = solve_ivp(system_mass_fraction, t_span=[ti, tfi],
y0=results[-1].T[-1],
max_step=dtmax_mf, method=self.method, args=[self],
rtol=self.rtol, atol=self.atol)
results.append(sol_ivp.y)
tps.append(sol_ivp.t)
ti = sol_ivp.t[-1]
if tfi < tf and compt < max(2, int(tf / tfi)+1):
dtmax_mf = max(dtmax_mf, yr)
tfi: np.float64 = min(tf, tps[-1][-1] + dtmax_mf * 1e5)
else:
tfi += 2 * tf
try:
if len(results[-1]) > 2:
for (y, t) in zip(results, tps):
step: int = max(len(t) // 150, 1)
self.get_final_mass_fractions(y.T[0:], tps=t[0:], step=step)
print("t=", t[-1] / yr, tps[-1] / yr)
if build_graphs:
self.graph_surface_density(times=np.geomspace(yr, self.tps_mf[-1], 10), idx=None,
species=ii_max,
plot_mean_volumique_density=False, include_last_index=False,
plot_sigma_crit=False, name_planets=None, show=False, save=True,
directory=self.directory, lw=2.)
if len(self.a_planets) > 0:
self.graphs_accretion(name_planets=None, show=False, directory=self.directory,
save=False, times=np.geomspace(yr, self.tps_mf[-1], 10),
absolute=False, include_last_index=False)
self.save()
else:
print(tps[-1])
if build_graphs:
self.graph_surface_density(times=np.geomspace(yr, self.tps_mf[-1], 10), idx=None,
species=ii_max,
plot_mean_volumique_density=False, include_last_index=False,
plot_sigma_crit=False, name_planets=None, show=False, save=True,
directory=self.directory, lw=2.)
if len(self.a_planets) > 0:
self.graphs_accretion(name_planets=None, show=False, directory=self.directory,
save=False, times=np.geomspace(yr, self.tps_mf[-1], 10),
absolute=False, include_last_index=False)
self.save()
raise RuntimeError("Convergence issue", sol_ivp.message)
except KeyboardInterrupt:
if build_graphs:
self.graph_surface_density(times=np.geomspace(yr, self.tps_mf[-1], 10), idx=None, species=ii_max,
plot_mean_volumique_density=False, include_last_index=False,
plot_sigma_crit=False, name_planets=None, show=False, save=True,
directory=self.directory, lw=2.)
if len(self.a_planets) > 0:
self.graphs_accretion(name_planets=None, show=False, directory=self.directory,
save=False, times=np.geomspace(yr, self.tps_mf[-1], 10),
absolute=False, include_last_index=False)
self.save()
raise UserWarning("Keyboard interrupt")
[docs]
def save(self, filename: str = "", directory: str = "./") -> None:
"""
Save the object in the compressed format .npz
Parameters
----------
filename : str, optional, default=self.filename
The name of the .npz file
directory : str, optional, default=self.filename
The name of the directory of the .npz file. If there is any "/" in
filename, then it will be assumed that filename contains the whole path
to the file and directory will be ignored
Returns
-------
None
"""
if "/" in filename:
i: int = filename.find("/")
j: int = i
while j != -1:
i = j
j = filename.find("/", i + 1)
directory = filename[:i]
filename = filename[i + 1:]
if filename != "":
self.name: str = filename
if directory != "./":
self.directory = directory
if os.path.exists("self.directory" + "/" + self.name):
tab = np.load("self.directory" + "/" + self.name)
if tab["tps"] >= self.tps or tab["tps_mf"] >= self.tps_mf:
raise UserWarning("The solution has already been integrated by another process, the Sol cannot be saved")
tab.close()
dic: dict = self.constantes()
dic["name"] = self.name
dic["IT_s"] = np.array(self.IT_s, dtype='double')
dic["r"] = np.array(self.r)
dic["f"] = np.array(self.f)
dic["mass_fractions"] = np.array(self.mass_fractions)
dic["flux_planets"] = np.array(self.flux_planets)
dic["idxa"] = np.array(self.idxa)
dic["idxk"] = np.array(self.idxk)
dic["tps"] = np.array(self.tps)
dic["tps_mf"] = np.array(self.tps_mf)
dic["mu"] = np.array(self.mu)
dic["nom_mu"] = np.array(self.nom_mu)
for k in self.diag.keys():
dic["diag_" + k] = np.array(self.diag[k])
for k in self.sol_ast_params.keys():
dic["sol_ast_params_" + k] = np.array(self.sol_ast_params[k])
for k in self.sol_ast_kwargs_N.keys():
dic["sol_ast_kwargs_N_" + k] = np.array(self.sol_ast_kwargs_N[k])
if self.sol_ast.sig_dots is not None:
dic["sol_ast_sig_dots"] = self.sol_ast.sig_dots
dic["sol_ast_tps"] = self.sol_ast.tps
dic["sol_ast_m_init"] = self.sol_ast.m_init
for k in self.sol_kuip_params.keys():
dic["sol_kuip_params_" + k] = np.array(self.sol_kuip_params[k])
for k in self.sol_kuip_kwargs_N.keys():
dic["sol_kuip_kwargs_N_" + k] = np.array(self.sol_kuip_kwargs_N[k])
if self.sol_kuip.sig_dots is not None:
dic["sol_kuip_sig_dots"] = self.sol_kuip.sig_dots
dic["sol_kuip_tps"] = self.sol_kuip.tps
dic["sol_kuip_m_init"] = self.sol_kuip.m_init
dic["idx_mu_ast"] = self.idx_mu_ast
dic["idx_mu_kuip"] = self.idx_mu_kuip
dic["a_planets"] = self.a_planets
dic["m_planets"] = self.m_planets
dic["flux_radmc"] = self.flux_radmc
np.savez_compressed(self.directory + "/" + self.name, **dic)
[docs]
def add_planets_to_graph(self, graph: g.Graphique, name_planets: list | np.ndarray = None,
show: bool = False) -> None:
"""
Add a gray-shaded area for each planet representing there influence zone (the zone where
they accrete) to a given Graphique (assuming the x-axis is in code unit)
Parameters
----------
graph : g.Graphique
The Graphique to which represent the planets
name_planets : list[str] | np.ndarray[str], optional, default=None
A list with the planet's name, if not None, the gray shaded area will be replaced by
colored shaded area (one for each planet) and the label will be name_planets
show : bool, default=False
To show the Graphique after the operation
Returns
-------
None
"""
colors_planets: list[str] = [g.C14, g.C9, g.C10, g.C11, g.C8, g.C20]
val_min: np.float64 = np.nanmin(graph.lines_y)
val_max: np.float64 = np.nanmax(graph.lines_y) * 10
if name_planets is None:
name_planets = [None for p in self.a_planets]
for (r_planet, m_planet, j) in \
zip(self.a_planets, self.m_planets, range(len(name_planets))):
ip = np.argwhere(self.r >= r_planet)[0, 0]
ind: np.ndarray = np.array([[self.r[ip - self.dip] / au, val_min],
[self.r[ip - self.dip] / au, val_max],
[self.r[ip + self.dip] / au, val_max],
[self.r[ip + self.dip] / au, val_min], ])
if name_planets[j] is None:
graph.polygon(ind=ind, facecolor="#929292",
plot_borders=False)
else:
name: str = name_planets[j]
graph.polygon(ind=ind, label=name, facecolor=colors_planets[j % len(colors_planets)],
plot_borders=False)
if show:
graph.show()
[docs]
def graphs_accretion(self, name_planets: list[str] | np.ndarray[str] | None = None,
show: bool = False, directory: str = None, save: bool = False,
times: np.ndarray[np.float64] = None, absolute: bool = False,
include_last_index: bool = True) -> list[g.Graphique]:
if directory is None:
directory = self.directory
# Accretion_rates
colors_planets: list[str] = [g.C1, g.C2, g.C6, g.C5, g.C3, g.C4]
colors_p: list[str] = []
mdots: np.ndarray = abs(np.array(self.flux_planets)[::max(len(self.flux_planets) // 1000, 1), 1]
- np.array(self.flux_planets)[::max(len(self.flux_planets) // 1000, 1), 0])
for i in range(len(self.a_planets)):
colors_p.append(colors_planets[i % len(colors_planets)])
time = np.array(self.diag["time"][::max(len(self.flux_planets) // 1000, 1)])
print(mdots.shape, time.shape)
# Accretion rates
if name_planets is not None and len(name_planets) == len(self.m_planets):
gr1: g.Graphique = g.loglog(time, mdots.T, linestyle="-",
label=name_planets, color=colors_p, show=False)
else:
gr1: g.Graphique = g.loglog(time, mdots.T, linestyle="-", color=colors_p, show=False)
gr1.config_ax(xlabel="Time [Myr]", ylabel="Accretion rate [M$_\oplus$. Myr$^{-1}$]",
xlim=[1e-3 * Myr, self.tps[-1]], ylim=[1e-11, 1e-1])
gr1.filename = self.name + "_accretion_rates"
if show:
gr1.show()
if save:
gr1.directory = directory
gr1.ext = ".png"
gr1.save_figure(dpi=400)
gr1.save()
# Accreted mass
if name_planets is not None and len(name_planets) == len(self.m_planets):
gr2: g.Graphique = g.loglog(time,
cumulative_trapezoid(mdots, time, initial=zero, axis=0).T,
label=name_planets, color=colors_p, show=False)
else:
gr2: g.Graphique = g.loglog(time,
cumulative_trapezoid(mdots, time, initial=zero, axis=0).T,
color=colors_p, show=False)
gr2.config_ax(xlabel="Time [Myr]", ylabel="Masses accreted [M$_\oplus$]",
xlim=[1e-3 * Myr, self.tps[-1]])
gr2.filename = self.name + "_accreted_masses"
if show:
gr2.show()
if save:
gr2.directory = directory
gr2.ext = ".png"
gr2.save_figure(dpi=400)
gr2.save()
return [gr1, gr2]
[docs]
def graph_surface_density(self, times: np.ndarray = None, idx: np.ndarray = None,
species: int | np.ndarray | list[int] = ii_max,
plot_mean_volumique_density: bool = False,
include_last_index: bool = False,
plot_sigma_crit: bool = True, name_planets: np.ndarray[str] | list[str] = None,
show: bool = True, save: bool = False, directory: str = None,
lw: np.float64 = 2.,
**args_ax
) -> g.Graphique:
"""
Create a Graphique of surface density as a function of radius for different time for given species
Parameters
----------
times : np.ndarray, optional,
table of the times we want to show the surface density.
If None and idx is None, an arbitrary subset of self.tps will be taken
idx : np.ndarray, optional
index table of self.tps we want to show the surface density.
If None, an arbitrary subset of self.tps will be taken
species : int | str | np.ndarray | list[int] | list[str], optional, default: plot all species
the species we want to show (if -1 plot the total surface density),
if < number of species plot the species associated with self.mu[species] else plot all the species
if str, the specie associated with this name is plotted
plot_mean_volumique_density : bool, optional, default=False
If True plot the volumique density instead of the surface density
include_last_index : bool, optional, default=False
Systematically include the last surface density
(even if self.tps[-1] is not in time)
plot_sigma_crit : bool, optional, default=True
To plot a line at the critical surface density to shield the main molecule
name_planets
The planet's names. If not None (default), plot a vertical line at the planet's positions
show : bool, optional, default=True
To show the Graphique before returning it
save : bool, optional, default=False
To save the Graphique
directory : str, optional, default=self.directory
The directory of the Graphique
lw : int, optional, default=2
The line-weight of surface density
args_ax : dict
additional keywords argument for the generation of the axis associated with the graphique (ex
x_lim=[1,100])
Returns
-------
g.Graphique
the requested Graphique
"""
H: np.ndarray = np.sqrt(
self.r * self.r * self.r * Rgp * T(self.r) / (G * self.Ms * self.mu.mean()))
if directory is None:
directory: str = self.directory
if len(self.nom_mu) == 1:
species = -1
if isinstance(species, list):
species = np.array(species)
if isinstance(species, int) and species >= len(self.mu):
species = np.arange(len(self.mu))
if isinstance(species, str) and species in self.nom_mu:
species = int(np.argwhere(self.nom_mu == species)[0, 0])
elif isinstance(species, str):
raise UserWarning("There is no ", species, "in this simulation, only ", self.nom_mu)
if isinstance(species, np.ndarray) and isinstance(species[0], str):
species_new = np.ones(len(species), dtype=int)
for i in range(len(species)):
if isinstance(species[i], str) and species[i] in self.nom_mu:
species_new[i] = int(np.argwhere(self.nom_mu == species[i])[0, 0])
elif isinstance(species[i], str):
raise UserWarning("There is no ", species[i], "in this simulation, only ", self.nom_mu)
species = species_new
if times is None and isinstance(species, int | np.int64) and species < 0:
times: np.ndarray = self.tps[::len(self.tps) // 5]
elif times is None:
times: np.ndarray = self.tps_mf[::len(self.tps_mf) // 5]
if idx is None:
idx: list = [0]
for t in times:
if isinstance(species, int | np.int64) and species < 0 and np.any(self.tps >= t):
idx.append(np.argwhere(self.tps >= t)[0, 0])
elif np.any(self.tps_mf >= t):
idx.append(np.argwhere(self.tps_mf >= t)[0, 0])
if include_last_index and -1 not in idx:
idx.append(-1)
gr = g.Graphique()
gr.config_ax(facecolor=(180 / 255, 180 / 255, 200 / 255, 0.1))
gr.config_spines(spine=["top", "right"], visible=False) # Remouving the top and right limit of the plot
gr.config_spines(spine=["bottom", "left"],
lw=2) # Setting the lineweight of the remaing borders of the figure to 2
gr.config_ticks(width=2, length=6) # Setting the ticks width and length
gr.config_fig(facecolor=g.Ctrensp) # Setting a trensparent background for the rest of the figure
linestyles_aw: list[str] = ["solid", "dashed", "dotted"]
colors_sigmas: list[list[str]] = [[g.C1, g.C2], [g.C17, g.C13], [g.C19, g.C3], [g.C15, g.C11]
, [g.C18, g.C14], [g.C9, g.C7]]
for i in range(len(idx)):
labels: list[str] = None
linestyles: list[str] | str = linestyles_aw[0]
if (type(species) is int and species < 0) and plot_mean_volumique_density:
values: np.ndarray = ((self.sigma(i) / (self.mu.mean() / Na))
* (C_L * C_L * C_L * 1e-6) / H)
colors = g.linear_color_interpolation(np.double(i), 0, np.double(len(idx)), colors_sigmas[0][0],
colors_sigmas[0][1])
elif isinstance(species, int | np.int64) and species < 0:
values: np.ndarray = (self.sigma(i) * au * au / Mearth)
colors = g.linear_color_interpolation(np.double(i), 0, np.double(len(idx)), colors_sigmas[0][0],
colors_sigmas[0][1])
elif ((isinstance(species, np.ndarray)) or species < len(self.mu)) and plot_mean_volumique_density:
values: np.ndarray = ((self.interp_sigma(self.tps_mf[idx[i]])
* self.mass_fractions[idx[i]][species] / (self.mu.mean() / Na))
* (C_L * C_L * C_L * 1e-6) / H)
labels = self.nom_mu[species]
if type(species) is int:
colors = g.linear_color_interpolation(np.double(i), zero, np.double(len(idx)),
colors_sigmas[0][0], colors_sigmas[0][1])
# linestyles = linestyles_aw[0]
else:
# linestyles = np.array([linestyles_aw[j % len(linestyles_aw)] for j in range(len(species))])
colors = [g.linear_color_interpolation(np.double(i), 0, np.double(len(idx)),
colors_sigmas[j % len(colors_sigmas)][0],
colors_sigmas[j % len(colors_sigmas)][1])
for j in range(len(species))]
elif (isinstance(species, np.ndarray)) or species < len(self.mu):
values: np.ndarray = (
self.interp_sigma(self.tps_mf[idx[i]])
* self.mass_fractions[idx[i]][species]
* au * au / Mearth)
labels = self.nom_mu[species]
if type(species) is int:
colors = g.linear_color_interpolation(np.double(i), 0, np.double(len(idx)), g.C1, g.C2)
# linestyles = linestyles_aw[0]
else:
# linestyles = np.array([linestyles_aw[j % len(linestyles_aw)] for j in range(len(species))])
colors = [g.linear_color_interpolation(np.double(i), 0, np.double(len(idx)),
colors_sigmas[j % len(colors_sigmas)][0],
colors_sigmas[j % len(colors_sigmas)][1])
for j in range(len(species))]
elif plot_mean_volumique_density:
values: np.ndarray = ((self.interp_sigma(self.tps_mf[idx[i]])
* self.mass_fractions[idx[i]][species]
/ (self.mu.mean() / Na))
* (C_L * C_L * C_L * 1e-6) / H)
colors = [g.linear_color_interpolation(np.double(i), 0, np.double(len(idx)),
colors_sigmas[j % len(colors_sigmas)][0],
colors_sigmas[j % len(colors_sigmas)][1])
for j in range(len(self.mu))]
# linestyles = np.array([linestyles_aw[j % len(linestyles_aw)] for j in range(len(self.mu))])
labels = self.nom_mu
else:
values: np.ndarray = (self.interp_sigma(self.tps_mf[idx[i]])
* self.mass_fractions[idx[i]] * au * au / Mearth)
colors = [g.linear_color_interpolation(np.double(i), 0, np.double(len(idx)),
colors_sigmas[j % len(colors_sigmas)][0],
colors_sigmas[j % len(colors_sigmas)][1])
for j in range(len(self.mu))]
# linestyles = np.array([linestyles_aw[j % len(linestyles_aw)] for j in range(len(self.mu))])
labels = self.nom_mu
if i == len(idx) - 1 and labels is not None:
gr.loglog(self.r / au, values, color=colors,
label=labels, lw=lw)
else:
gr.loglog(self.r / au, values, color=colors, lw=lw)
self.add_planets_to_graph(gr, name_planets)
if type(species) is int and species < 0:
gr.customized_cmap(np.array(self.tps)[idx],
g.linear_color_interpolation(np.arange(len(idx)),
col_min=colors_sigmas[0][0],
col_max=colors_sigmas[0][1]),
label="time [Myr]",
ticks=np.array(self.tps)[idx][::max(1, len(idx) // 10)],
format="{:.2g}", size=0.02)
elif isinstance(species, int | np.int64):
gr.customized_cmap(np.array(self.tps_mf)[idx],
g.linear_color_interpolation(np.arange(len(idx)),
col_min=colors_sigmas[0][0],
col_max=colors_sigmas[0][1]),
label=self.nom_mu[species] + " : time [Myr]",
ticks=np.array(self.tps_mf)[idx][::max(1, len(idx) // 10)],
format="{:.2g}", size=0.02)
else:
locs: list[str] = ['right']
fracs: list[np.double] = [1.]
for j in range(len(species)):
ticks: np.ndarray = np.array([])
if j == 0:
ticks = np.array(self.tps_mf)[idx][::max(1, len(idx) // 10)]
gr.customized_cmap(np.array(self.tps_mf)[idx],
g.linear_color_interpolation(np.arange(len(idx)),
col_min=colors_sigmas[j % len(colors_sigmas)][0],
col_max=colors_sigmas[j % len(colors_sigmas)][1]),
label="Time [Myr]", ticks=ticks,
location=locs[j % len(locs)], fraction=fracs[j % len(locs)]
, format="{:.2g}", size=0.02)
else:
gr.customized_cmap(np.array(self.tps_mf)[idx],
g.linear_color_interpolation(np.arange(len(idx)),
col_min=colors_sigmas[j % len(colors_sigmas)][0],
col_max=colors_sigmas[j % len(colors_sigmas)][1]),
location=locs[j % len(locs)], fraction=fracs[j % len(locs)], format="{:.2g}",
ticks=ticks, share_axis=True, size=0.02, space_between=0.005)
sigma_crit_co: np.float64 = zero
sigma_crit_h2o: np.float64 = zero
if "CO" in self.nom_mu:
sigma_crit_c0: np.float64 = np.double(1e-7) * Mearth / (au * au)
elif "H2O" in self.nom_mu:
cross_section_h2o: np.float64 = np.double(5.e-22) * C_L * C_L # mean cross-section
sigma_crit_h2o: np.float64 = np.double(1. / ((Na / (18 * 1e-2 * C_M)) * cross_section_h2o))
elif plot_sigma_crit:
print("There neither water or CO in the species, there is no critical surface density to plot")
if plot_sigma_crit and plot_mean_volumique_density:
if sigma_crit_co > zero:
gr.loglog(self.r, (sigma_crit_co / (self.mu.mean() / Na))
* (C_L * C_L * C_L * 1e-6) / H, label="$\sigma_{crit}$(CO)", color=g.C19)
if sigma_crit_h2o > zero:
gr.loglog(self.r, (sigma_crit_h2o / (self.mu.mean() / Na))
* (C_L * C_L * C_L * 1e-6) / H, label="$\sigma_{crit}(H_2O)$", color=g.C19)
elif plot_sigma_crit:
if sigma_crit_co > zero:
gr.loglog([self.r[0], self.r[-1]],
[sigma_crit_co * au * au / Mearth, sigma_crit_co * au * au / Mearth],
label="$\sigma_{crit}$(co)",
color=g.C19)
if sigma_crit_h2o > zero:
gr.loglog([self.r[0], self.r[-1]],
[sigma_crit_h2o * au * au / Mearth, sigma_crit_h2o * au * au / Mearth],
label="$\sigma_{crit}(H_2O)$",
color=g.C19)
gr.title = r"$\alpha$=" + str(self.alpha)
# gr.titre = r"$\alpha$=" + str(self.alpha) + r" $\dot{m}=$" + str(
# self.m_dot * (C_t * 3600 * 24 * 365.25 * 1e6) / Mearth) + r" $M$_dot$.Myr$^{-1}$"
if "xlabel" not in args_ax.keys():
args_ax["xlabel"] = 'Distance to star [au]'
if "ylabel" not in args_ax.keys():
if plot_mean_volumique_density:
args_ax["ylabel"] = r"C$^0$ number density [cm$^{-3}$]"
else:
args_ax["ylabel"] = r"Surface density [M$_\oplus$ . au$^{-2}$]"
gr.config_ax(**args_ax)
gr.filename = "Alpha_" + str(self.alpha) + "_Sigmas"
if show:
gr.show()
if save:
gr.directory = directory
gr.ext = ".png"
gr.save_figure(dpi=400)
gr.save()
return gr
[docs]
def graph_image_surface_density(self, species: int | np.ndarray | list[int] = ii_max,
plot_sigma_crit: bool = True,
vmin: np.float64 | str = 1e-10 * Mearth / Myr, vmax: np.float64 | str= None,
size_time: int = 500, size_radius: int = 1000,
color_sigma_crit_h2o: str | tuple = "w",
color_sigma_crit_co: str | tuple = "g",
cmap: str = "inferno", nb_levels: int = 10,
color_min: str | tuple = None, color_max: str | tuple = None,
plot_planets: bool = True, planets_names: list[str] = None,
show: bool = True, save: bool = False, directory: str = None,
plot_label_levels: bool = True,
**kwargs_ax
) -> g.Graphique:
"""
Build a Graphique that represant the evolution of surface density as function of time and radius
Parameters
----------
species : str | int, optional
The specie for which the surface density is represanted, default the total surface density is represanted
plot_sigma_crit : bool, optional, default = True
To plot a specific level line at the critical surface density
vmin : np.float64 | str , optional, {"auto", np.loat64}, default = 1e-10 * Mearth / Myr
The minimum value of surface density to represant in the result (below it will be considered as saturated)
If "auto", the minimal value will be the minimal value of the surface density represanted
vmax : np.float64 | str, optional, {"auto", np.loat64},
The maximum value of surface density to represant in the result (above it will be considered as saturated)
Default, v_max is not considered like for "auto" where the maximal value will be the maximal value of the surface density represanted
size_time : int, optional, default=1000
The image's size for the time coordinate (x-axis)
size_radius : int, optional, default=1000
The image's size for the radius coordinate (y-axis)
color_sigma_crit_h2o : str, optional, default = "w"
The color of the level line associated with the critical surface density of water
color_sigma_crit_co : str, optional, default = "g"
The color of the level line associated with the critical surface density of CO
cmap : str, optional, default = "inferno"
The colormap used for the image
nb_levels : int, optional, default=10
The number of levels in the image
The number of levels to plot in addition to the image
color_min : str | tuple, optional
The color associated with the minimum value for a `default` cmap
color_max : str | tuple, optional
The color associated with the maximum value for a `default` cmap
plot_planets : bool, optional, default = True
Whether or not to plot a grey shaded era that delimitate the planet's sink cell's era
planets_names : list[str], optional
A list of names for the planets to write close to the grey shadeds eras
plot_label_levels : bool, optional, default=False
To plot (or not) labels for the contours (if True, labels are only plot every two levels)
show : bool, optional, default = True
Whether to show the image
save : bool, optional, default = False
Whether to save the image (if True, both the Graphique and the png image will be saved)
directory : str, optional
The directory to save the Graphique, default=self.directory
Returns
-------
g.Graphique
"""
if isinstance(species, str) and species != "all" and species not in self.nom_mu:
raise UserWarning(f"The specie {species} doesn't exist for this simulation,"
f" Please use 'all' for the total surface density of one of{self.nom_mu}")
if isinstance(species, int) and len(self.nom_mu) <= species < ii_max:
print(f"Graphique.graph_image_surface_density : the specie number {species} cannot be reffered to any"
f" specie of this solution. The total surface density will be represanted")
if isinstance(vmin, str) and vmin !="auto":
raise UserWarning(f"The vmin={vmin} is not a valid value for this simulation")
elif isinstance(vmin, str):
vmin = -inf
if vmax is None:
vmax = np.inf
if isinstance(vmax, str) and vmax !="auto":
raise UserWarning(f"The vmax={vmax} is not a valid value for this simulation")
elif isinstance(vmax, str):
vmax = -inf
if planets_names is not None and len(planets_names) != len(self.a_planets):
raise UserWarning(f"The number of planets {len(planets_names)} in the name_planet parameter does not match"
f"to the actual number of planets in this simulation {len(self.a_planets)}")
sigma_crit_co: np.float64 = zero
sigma_crit_h2o: np.float64 = zero
if "CO" in self.nom_mu:
sigma_crit_co: np.float64 = np.double(1e-7) * Mearth / (au * au)
elif "H2O" in self.nom_mu:
cross_section_h2o: np.float64 = np.double(5.e-22) * C_L * C_L # mean cross-section
sigma_crit_h2o: np.float64 = np.double(1. / ((Na / (18 * 1e-2 * C_M)) * cross_section_h2o))
if ((isinstance(species, int) and (species > len(self.nom_mu) or species < 0))
or (isinstance(species, str) and species == "all")):
tps = np.geomspace(np.min(np.array(self.tps)[np.array(self.tps) > 0]), np.max(self.tps), size_time)
gr: g.Graphique = g.image(self.interp_sigma(tps)[:,
::max(len(self.r) // size_radius, 1)].T,
tps / Myr,
np.array(self.r)[::max(len(self.r) // size_radius, 1)] / au,
cmap=cmap, colorscale="log", show=False, vmin=vmin, vmax=vmax,
color_min=color_min, color_max=color_max)
if vmin == -inf:
vmin = np.min(gr.array_image)
if vmax < inf:
vmax_levels = vmax
else:
vmax_levels: np.float64 = np.max(gr.array_image)
if nb_levels > 0:
levels = np.geomspace(vmin, vmax_levels, nb_levels)
if "CO" in self.nom_mu and "H2O" in self.nom_mu and plot_sigma_crit:
gr.contours(np.append(levels, [sigma_crit_co, sigma_crit_h2o]),
color=np.append(["k" for l in levels], color_sigma_crit_co, color_sigma_crit_h2o))
elif "CO" in self.nom_mu and plot_sigma_crit:
gr.contours(np.append(levels, sigma_crit_co),
color=np.append(["k" for l in levels], color_sigma_crit_co))
elif "H2O" in self.nom_mu and plot_sigma_crit:
gr.contours(np.append(levels, sigma_crit_h2o),
color=np.append(["k" for l in levels], color_sigma_crit_h2o))
else:
gr.contours(levels, color="k")
gr.config_colorbar(ticks=levels,
ticks_labels=["{:.2e}".format(t) for t in
levels],
label="Total surface density [M$_\oplus$ . au$^{-2}]$")
else:
if "CO" in self.nom_mu and "H2O" in self.nom_mu and plot_sigma_crit:
gr.contours([sigma_crit_co, sigma_crit_h2o],
colors=[color_sigma_crit_co, color_sigma_crit_h2o])
elif "CO" in self.nom_mu and plot_sigma_crit:
gr.contours(sigma_crit_co, colors=color_sigma_crit_co)
elif "H2O" in self.nom_mu and plot_sigma_crit:
gr.contours(sigma_crit_h2o, colors=color_sigma_crit_h2o)
gr.config_colorbar(ticks=np.geomspace(vmin, vmax_levels, 10),
ticks_labels=["{:.2e}".format(t) for t in
np.geomspace(vmin, vmax_levels, 10)],
label="Total surface density [M$_\oplus$ . au$^{-2}]$")
gr.config_ax(xlabel="Time [Myr]", ylabel="Distance to the central star [au] ",
xscale="log", yscale="log", xlim=[self.tps[2], self.tps[-1]], ylim=[self.r[0], self.r[-1]])
if vmax < np.inf and np.any(self.sigma() > vmax):
gr.config_colorbar(0, extend="max")
if np.any(self.sigma() < vmin):
if 'extend' in gr.param_colorbar[0].keys():
gr.config_colorbar(0, extend="both")
else:
gr.config_colorbar(0, extend="min")
gr.filename = "Image_total_surface_density"
else:
if isinstance(species, str):
species = np.argwhere(self.nom_mu == species)[0, 0]
idx = np.arange(0, len(self.tps_mf), max(len(self.tps_mf) // size_time, 1))
gr: g.Graphique = g.image(self.interp_sigma(np.array(self.tps_mf)[idx]
)[:, ::max(len(self.r) // size_radius, 1)].T
* np.array(self.mass_fractions)[idx, species,
::max(len(self.r) // size_radius,1)].T,
np.array(self.tps_mf)[idx] / Myr,
np.array(self.r)[::max(len(self.r) // size_radius, 1)] / au,
cmap=cmap, colorscale="log", show=False, vmin=vmin, vmax=vmax,
color_min=color_min, color_max=color_max)
if vmin == -inf:
vmin = np.min(gr.array_image)
if vmax < inf:
vmax_levels = vmax
else:
vmax_levels: np.float64 = np.max(gr.array_image)
if nb_levels > 0:
levels = np.geomspace(vmin, vmax_levels, nb_levels)
if "CO" in self.nom_mu and "H2O" in self.nom_mu and plot_sigma_crit:
gr.contours(np.append(levels, [sigma_crit_co, sigma_crit_h2o]),
color=np.append(["k" for l in levels], color_sigma_crit_co, color_sigma_crit_h2o))
elif "CO" in self.nom_mu and plot_sigma_crit:
gr.contours(np.append(levels, sigma_crit_co),
color=np.append(["k" for l in levels], color_sigma_crit_co))
elif "H2O" in self.nom_mu and plot_sigma_crit:
gr.contours(np.append(levels, sigma_crit_h2o),
color=np.append(["k" for l in levels], color_sigma_crit_h2o))
else:
gr.contours(levels, color="k")
gr.config_colorbar(ticks=levels,
ticks_labels=["{:.2e}".format(t) for t in
levels],
label=self.nom_mu[species] + " surface density [M$_\oplus$ . au$^{-2}]$")
else:
if "CO" in self.nom_mu and "H2O" in self.nom_mu and plot_sigma_crit:
gr.contours([sigma_crit_co, sigma_crit_h2o],
colors=[color_sigma_crit_co, color_sigma_crit_h2o])
elif "CO" in self.nom_mu and plot_sigma_crit:
gr.contours(sigma_crit_co, colors=color_sigma_crit_co)
elif "H2O" in self.nom_mu and plot_sigma_crit:
gr.contours(sigma_crit_h2o, colors=color_sigma_crit_h2o)
gr.config_colorbar(ticks=np.geomspace(vmin, vmax_levels, 10),
ticks_labels=["{:.2e}".format(t) for t in
np.geomspace(vmin, vmax_levels, 10)],
label=self.nom_mu[species] + " surface density [M$_\oplus$ . au$^{-2}]$")
gr.config_ax(xlabel="Time [Myr]", ylabel="Distance to the central star [au] ",
xscale="log", yscale="log", xlim=[self.tps_mf[1], self.tps_mf[-1]], ylim=[self.r[0], self.r[-1]])
if vmax < np.inf and np.any(self.sigma() > vmax):
gr.config_colorbar(0, extend="max")
if np.any(self.sigma() < vmin):
if 'extend' in gr.param_colorbar[0].keys():
gr.config_colorbar(0, extend="both")
else:
gr.config_colorbar(0, extend="min")
gr.filename = "Image_" + self.nom_mu[species] + "_surface_density"
gr.config_ax(**kwargs_ax)
if plot_planets:
for (r_planet, m_planet, j) in \
zip(self.a_planets, self.m_planets, range(len(self.a_planets))):
ip = np.argwhere(self.r >= r_planet)[0, 0]
ind: np.ndarray = np.array([[gr.param_ax["xlim"][0], self.r[ip - self.dip] / au],
[gr.param_ax["xlim"][1], self.r[ip - self.dip] / au],
[gr.param_ax["xlim"][1], self.r[ip + self.dip] / au],
[gr.param_ax["xlim"][0], self.r[ip + self.dip] / au], ])
gr.polygon(ind=ind, facecolor="#929292", plot_borders=False)
if planets_names is not None:
gr.text(gr.param_ax["xlim"][0], self.r[ip - self.dip] / au, planets_names[j],
horizontalalignment="left", verticalalignment="top", rotation="horizontal",
color="#929292")
if plot_label_levels:
gr.config_labels_contours(levels=gr.levels[::2], fmt="%.0e")
else:
gr.config_labels_contours(levels=[], fmt="%.0e")
if directory is not None:
gr.directory = directory
else:
gr.directory = self.directory
if show:
gr.show()
if save:
gr.ext = ".png"
gr.save_figure(dpi=400)
gr.save()
return gr
[docs]
def graphs(self, plot_planets=True, planets_names: list[str] = None,
show: bool = True, save: bool = False, directory: str = None) -> list[g.Graphique]:
print("graph accretion : ")
res: list[g.Graphique] = self.graphs_accretion(name_planets=planets_names, show=show,
save=save, directory=directory)
if directory is None:
directory = self.directory
print("graph image sig_dot :")
res.append(self.sol_ast.graph_image_surface_density(show=show, save=save, directory=directory))
print("graph image surface density:")
res.append(self.graph_image_surface_density(plot_planets=plot_planets, planets_names=planets_names,
show=show, save=save, directory=directory))
for mu_plot in self.nom_mu:
print("mu = ", mu_plot)
res.append(self.graph_image_surface_density(species=mu_plot,
plot_planets=plot_planets, planets_names=planets_names,
show=show, save=save, directory=directory))
return res
[docs]
def radmc_setup(self, i: int, directory: str = None, size_grid: int = 128) -> None:
"""
Generate the configuration file for radmc3d for the surface density at time self.tps[i]
Parameters
----------
i: int
The index at which the surface density is plotted
directory : str, default=self.directory+'/'+i
The directory where the configuration are stored
size_grid : int, default=128
THe size of the spatial grid (the water surface density is interpolated on it)
Returns
-------
None
"""
if directory is not None:
if not os.path.exists(directory):
os.mkdir(directory)
doss: str = directory
else:
if not os.path.exists(self.directory + "/" + str(i)):
os.mkdir(self.directory + "/" + str(i))
doss: str = self.directory + "/" + str(i)
nphot: int = 1000000 # Monte Carlo parameters
# Grid parameters
nr: int = size_grid
nlev_rin: int = 12 # Grid refinement at the inner edge: nr of cycles
nspan_rin: int = 3 # Grid refinement at the inner edge: nr of cells each cycle
ntheta: int = 32
nphi: int = 1
rin: np.float64 = self.r.min()
rout: np.float64 = self.r.max()
thetaup: np.float64 = np.pi * 0.5 - 1.
# Disk parameters
pstar: np.ndarray[np.float64] = np.array([0., 0., 0.])
# Make the coordinates
# ri: np.ndarray[np.float64] = self.r[::max(1, len(self.r) // size_grid)]
ri: np.ndarray[np.float64] = np.geomspace(rin, rout, nr + 1)
# ri: np.ndarray[np.float64] = grid_refine_inner_edge(ri, nlev_rin, nspan_rin) # Refinement at inner edge
thetai: np.ndarray[np.float64] = np.linspace(thetaup, 0.5 * Pi, ntheta + 1)
phii: np.ndarray[np.float64] = np.linspace(0., Pi * 2, nphi + 1)
rc: np.ndarray[np.float64] = 0.5 * (ri[:-1] + ri[1:])
nr: int = len(rc) # Recompute nr, because of refinement at inner edge
thetac: np.ndarray[np.float64] = 0.5 * (thetai[0:ntheta] + thetai[1:ntheta + 1])
phic: np.ndarray[np.float64] = 0.5 * (phii[0:nphi] + phii[1:nphi + 1])
# Make the grid
qq: np.ndarray[np.float64] = np.meshgrid(rc, thetac, phic, indexing='ij')
rr: np.ndarray[np.float64] = qq[0]
tt: np.ndarray[np.float64] = qq[1]
zr: np.ndarray[np.float64] = Pi / 2. - qq[1]
# Make the temperature model
temp: np.ndarray[np.float64] = T(rr, self.tps_mf[i])
# Make the density model
# g.loglog(rr.flatten(), temp.flatten(), ".")
h_h2o: np.ndarray[np.float64] = np.sqrt(kb * temp * rr * rr * rr / (18. * mp * G * self.Ms))
# g.loglog([self.r, rr.flatten()], [ self.interp_sigma(self.tps_mf[i]) * self.mass_fractions[i][self.idx_mu_ast[0]], np.interp(rr, self.r, self.interp_sigma(self.tps_mf[i]) * self.mass_fractions[i][self.idx_mu_ast[0]]).flatten()], ".")
sigma_H2O: np.ndarray[np.float64] = np.interp(rr, self.r, self.interp_sigma(self.tps_mf[i]) * self.mass_fractions[i][self.idx_mu_ast[0]])
rhog_h2o: np.ndarray[np.float64] = (sigma_H2O / (np.sqrt(2. * np.pi) * h_h2o)
) * np.exp(-zr * zr * rr * rr / (h_h2o * h_h2o * 2.))
# Make the velocity model
vr: np.ndarray[np.float64] = np.zeros_like(rr)
vtheta: np.ndarray[np.float64] = np.zeros_like(rr)
vphi: np.ndarray[np.float64] = np.sqrt(G * self.Ms / rr)
vturb: np.ndarray[np.float64] = 0.001 * vphi
# Write the wavelength_micron.inp file
lam1: np.float64 = 1.
lam2: np.float64 = 9.e2
lam3: np.float64 = 4.e3
lam4: np.float64 = 1.0e4
n12: np.float64 = 200
n23: np.float64 = 500
n34: np.float64 = 30
lam12: np.ndarray[np.float64] = np.logspace(np.log10(lam1), np.log10(lam2), n12, endpoint=False)
lam23: np.ndarray[np.float64] = np.logspace(np.log10(lam2), np.log10(lam3), n23, endpoint=False)
lam34: np.ndarray[np.float64] = np.logspace(np.log10(lam3), np.log10(lam4), n34, endpoint=True)
lam: np.ndarray[np.float64] = np.concatenate([lam12, lam23, lam34])
nlam: int = lam.size
# Write the wavelength file
with open(doss + '/wavelength_micron.inp', 'w+') as f:
f.write('%d\n' % nlam)
for value in lam:
f.write('%13.6e\n' % value)
# NB All datas needs to be saved in cgs system....
# Write the stars.inp file
with open(doss + '/stars.inp', 'w+') as f:
f.write('2\n')
f.write('1 %d\n\n' % nlam)
f.write('%13.6e %13.6e %13.6e %13.6e %13.6e\n\n' % (interpRs(self.tps_mf[i], self.Ms) * Rsun * 1e2 / C_L,
self.Ms * 1e3 / C_M,
pstar[0], pstar[1], pstar[2]))
for value in lam:
f.write('%13.6e\n' % value)
f.write('\n%13.6e\n' % (- ((interpLbol(self.tps_mf[i], self.Ms)
/ (4. * Pi * (interpRs(self.tps_mf[i], self.Ms) ** 2) * sigma))) ** (1. / 4.)))
# Write the grid file
with open(doss + '/amr_grid.inp', 'w+') as f:
f.write('1\n') # iformat
f.write('0\n') # AMR grid style (0=regular grid, no AMR)
f.write('100\n') # Coordinate system: spherical
f.write('0\n') # gridinfo
f.write('1 1 0\n') # Include r,theta coordinates
f.write('%d %d %d\n' % (nr, ntheta, 1)) # Size of grid
for value in ri * 1e-2 / C_L:
f.write('%13.6e\n' % value) # X coordinates (cell walls)
for value in thetai:
f.write('%13.6e\n' % value) # Y coordinates (cell walls)
for value in phii:
f.write('%13.6e\n' % value) # Z coordinates (cell walls)
# Write the density file
with open(doss + '/dust_density.inp', 'w+') as f:
f.write('1\n') # Format number
f.write('%d\n' % (nr * ntheta * nphi)) # Nr of cells
f.write('1\n') # Nr of dust species
data = np.zeros_like(rhog_h2o).ravel(order='F') # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%13.6e")
f.write('\n')
# with open(doss + '/dust_density.inp', 'w+') as f:
# f.write('1\n') # Format number
# f.write('%d\n' % (nr * ntheta * nphi)) # Nr of cells
# f.write('0\n') # Nr of dust species
# # data = np.zeros_like(rhog_h2o).ravel(order='F') # Create a 1-D view, fortran-style indexing
# # data.tofile(f, sep='\n', format="%13.6e")
# # f.write('\n')
# Dust opacity control file
with open(doss + '/dustopac.inp', 'w+') as f:
f.write('2 Format number of this file\n')
f.write('1 Nr of dust species\n')
f.write('============================================================================\n')
f.write('1 Way in which this dust species is read\n')
f.write('0 0=Thermal grain\n')
f.write('silicate Extension of name of dustkappa_***.inp file\n')
f.write('----------------------------------------------------------------------------\n')
# with open(doss + '/dustopac.inp', 'w+') as f:
# f.write('2 Format number of this file\n')
# f.write('0 Nr of dust species\n')
# f.write('============================================================================\n')
# f.write('----------------------------------------------------------------------------\n')
# Write the molecule number density file.
mu_h2o = 18 # Mass of a H2O molecule in proton mass unity
nh2o = rhog_h2o / (mu_h2o * mp)
print("test nh2o", nh2o.sum(), rhog_h2o.sum(), sigma_H2O.sum())
with open(doss + '/numberdens_h2op.inp', 'w+') as f:
f.write('1\n') # Format number
f.write('%d\n' % (nr * ntheta * nphi)) # Nr of cells
data = nh2o.ravel(order='F') * 1e-6 / (C_L * C_L * C_L) # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%13.6e")
f.write('\n')
with open(doss + '/numberdens_h2oo.inp', 'w+') as f:
f.write('1\n') # Format number
f.write('%d\n' % (nr * ntheta * nphi)) # Nr of cells
data = nh2o.ravel(order='F') * 1e-6 / (C_L * C_L * C_L) # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%13.6e")
f.write('\n')
# Write the gas velocity field
with open(doss + '/gas_velocity.inp', 'w+') as f:
f.write('1\n') # Format number
f.write('%d\n' % (nr * ntheta * nphi)) # Nr of cells
for iphi in range(nphi):
for itheta in range(ntheta):
for ir in range(nr):
f.write(
'%13.6e %13.6e %13.6e\n' % (vr[ir, itheta, iphi] * 1e-2 * C_t / C_L,
vtheta[ir, itheta, iphi] * 1e-2 * C_t / C_L,
vphi[ir, itheta, iphi] * 1e-2 * C_t / C_L))
# Write the temperature file
with open(doss + '/gas_temperature.inp', 'w+') as f:
f.write('1\n') # Format number
f.write('%d\n' % (nr * ntheta * nphi)) # Nr of cells
data = temp.ravel(order='F') # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%13.6e")
f.write('\n')
# Write the microturbulence file
with open(doss + '/microturbulence.inp', 'w+') as f:
f.write('1\n') # Format number
f.write('%d\n' % (nr * ntheta * nphi)) # Nr of cells
data = vturb.ravel(order='F') * 1e-2 * C_t / C_L # Create a 1-D view, fortran-style indexing
data.tofile(f, sep='\n', format="%13.6e")
f.write('\n')
# Write the lines.inp control file
with open(doss + '/lines.inp', 'w') as f:
f.write('2\n')
f.write('2\n')
f.write('h2oo leiden 0 0 0\n')
f.write('h2op leiden 0 0 0\n')
# Write the radmc3d.inp control file
with open(doss + '/radmc3d.inp', 'w+') as f:
f.write('nphot = %d\n' % nphot)
f.write('scattering_mode_max = 0\n')
f.write('iranfreqmode = 1\n')
# f.write('tgas_eq_tdust = 1\n') # Use the dust temperature of dust species 1 as gas temperature
shutil.copyfile(data_radmc + "/molecule_h2oo.inp", doss + "/molecule_h2oo.inp" )
shutil.copyfile(data_radmc + "/molecule_h2op.inp", doss + "/molecule_h2op.inp" )
shutil.copyfile(data_radmc + "/dustkappa_silicate.inp", doss + "/dustkappa_silicate.inp" )
[docs]
def calc_single_flux_radmc(self, i: int, directory: str = None, size_grid: int = 128,
incl: np.float64 = 13, phi: np.float64 = 59.) -> np.ndarray[np.float64]:
"""
Parameters
----------
i: int
The index at which the surface density is plotted
directory : str, default=self.directory+'/'+i
The directory where the configuration are stored
size_grid : int, default=128
The size of the spatial grid (the water surface density is interpolated on it)
Returns
-------
np.ndarray[np.float64]
The spatially integrated luminosity flux
"""
if len(self.flux_radmc) == 0 or i not in np.array(self.flux_radmc)[:, 0]:
self.radmc_setup(i=i, directory=directory, size_grid=size_grid)
if directory is None:
if not os.path.exists(self.directory + "/" + str(i)):
os.mkdir(self.directory + "/" + str(i))
directory: str = self.directory + "/" + str(i)
cwd: str = os.getcwd()
os.chdir(directory)
vkms: np.float64 = 0.
widthkms: np.float64 = 20.0
npix: np.float64 = 400
linenlam: np.float64 = 40
os.system("radmc3d mctherm")
# H20p_5
print("H2O p-5")
command = 'radmc3d image imolspec 2 iline 6 vkms {} widthkms {} incl {} phi {}' \
' npix {} linenlam {}'.format(vkms, widthkms, incl, phi, npix, linenlam)
os.system(command)
im = radmc_im.readImage()
tab_vkms = np.linspace(vkms - widthkms, vkms + widthkms, linenlam)
im.image = trapezoid(im.image, x=tab_vkms, axis=2).reshape((npix, npix, 1))
im.imageJyppix = trapezoid(im.imageJyppix, x=tab_vkms, axis=2).reshape((npix, npix, 1))
im.nfreq = 1
im.freq = np.array([np.mean(im.freq)])
im.nwav = 1
im.wav = np.array([np.mean(im.freq)])
res_H2Op5: np.float64 = im.imageJyppix.sum()
# H20p_7
command = 'radmc3d image imolspec 2 iline 12 vkms {} widthkms {} incl {} phi {}' \
' npix {} linenlam {}'.format(vkms, widthkms, incl, phi, npix, linenlam)
os.system(command)
im = radmc_im.readImage()
tab_vkms = np.linspace(vkms - widthkms, vkms + widthkms, linenlam)
im.image = trapezoid(im.image, x=tab_vkms, axis=2).reshape((npix, npix, 1))
im.imageJyppix = trapezoid(im.imageJyppix, x=tab_vkms, axis=2).reshape((npix, npix, 1))
im.nfreq = 1
im.freq = np.array([np.mean(im.freq)])
im.nwav = 1
im.wav = np.array([np.mean(im.freq)])
res_H2Op7: np.float64 = im.imageJyppix.sum()
# H20o_7
command = 'radmc3d image imolspec 2 iline 42 vkms {} widthkms {} incl {} phi {}' \
' npix {} linenlam {}'.format(vkms, widthkms, incl, phi, npix, linenlam)
os.system(command)
im = radmc_im.readImage()
tab_vkms = np.linspace(vkms - widthkms, vkms + widthkms, linenlam)
im.image = trapezoid(im.image, x=tab_vkms, axis=2).reshape((npix, npix, 1))
im.imageJyppix = trapezoid(im.imageJyppix, x=tab_vkms, axis=2).reshape((npix, npix, 1))
im.nfreq = 1
im.freq = np.array([np.mean(im.freq)])
im.nwav = 1
im.wav = np.array([np.mean(im.freq)])
res_H2Oo7: np.float64 = im.imageJyppix.sum()
os.system('rm image.out')
os.chdir(cwd)
self.flux_radmc.append(np.array([i, res_H2Op5, res_H2Op7, res_H2Oo7]))
return res_H2Op5, res_H2Op7, res_H2Oo7
else:
return self.flux_radmc[np.argwhere(np.array(self.flux_radmc) == i)[0, 0]][1:]
[docs]
def calc_flux_radmc(self, temporal_size: int = 50, size_grid: int = 20,
incl: np.float64 = 13, phi: np.float64 = 59.) -> None:
indexs: np.ndarray[int] = np.arange(0, len(self.tps_mf), max(1, len(self.tps_mf) // temporal_size))
for i in indexs:
self.calc_single_flux_radmc(i=i, size_grid=size_grid, incl=incl, phi=phi)
[docs]
def mass_steady_state(sol: Sol, mdot: np.float64=None, a0: np.float64 = None, mu: np.float64 = None
) -> np.float64:
"""
Return the mass of a steady state disc where gas is produced at radius a0 with a mass generation rate mdot
Parameters
----------
sol : Sol
The associated Sol
mdot : np.float64
Mass input rate at a0, if mtot is not None, this parameter is ignored
a0 : np.float64, optional, default=sol.a0
The radius at wich the gas is produced
mu : np.float64
Mean molar mass, default sol.mu.mean()
Returns
-------
np.ndarray
the surface density of the steady state
"""
if mu is None:
mu = np.mean(sol.mu) # Mean molar mass
if a0 is None:
a0 = sol.a0
nu0: np.float64 = (Rgp * T0(np.double(100.) * Myr + sol.t0, a0, sol.Ms) * sol.alpha * (a0 ** (3 / 2))
/ (mu * np.sqrt(G * sol.Ms)))
return ((mdot / (3. * Pi * nu0))
* (2. * Pi * (
((a0 ** (-pls_temp + 1 / 2)) - (sol.a_in ** (-pls_temp + 1 / 2)))
/ (a0 ** (-pls_temp - 3 / 2) * (-pls_temp + 1 / 2))
+
((sol.a_out ** (-pls_temp)) - (a0 ** (-pls_temp)))
/ (a0 ** (-pls_temp - 2) * (-pls_temp)))))
[docs]
def steady_state(sol: Sol, mdot: np.float64=None, mtot: np.float64 = None, a0: np.float64 = None, mu: np.float64 = None ):
"""
Return the log of steady state for a given mdot
Parameters
----------
sol : Sol
The associated Sol
mdot : np.float64
Mass input rate at a0, if mtot is not None, this parameter is ignored
mtot : np.float64, optional
The disc's total mass,
a0 : np.float64, optional, default=sol.a0
The radius at wich the gas is produced
mu : np.float64
Mean molar mass, default sol.mu.mean()
Returns
-------
np.ndarray
the surface density of the steady state
"""
if mu is None:
mu = np.mean(sol.mu) # Mean molar mass
if a0 is None:
a0 = sol.a0
if mtot is None:
nu0: np.float64 = (Rgp * T0(np.double(100.) * Myr + sol.t0, a0, sol.Ms) * sol.alpha * (a0 ** (3 / 2))
/ (sol.mu.mean() * np.sqrt(G * sol.Ms)))
f_st: np.ndarray = (mdot * ((1. / a0) ** (-pls_temp - 3 / 2)) / (3. * Pi * nu0)) * sol.x
idx: np.ndarray = np.argwhere(sol.r > a0)[:, 0]
f_st[idx] = mdot * ((1. / a0) ** (-pls_temp - 2.)) / (3. * Pi * nu0)
return f_st
else:
sigma0: np.float64 = mtot / (2. * Pi * (
((a0 ** (-pls_temp + 1 / 2)) - (sol.a_in ** (-pls_temp + 1 / 2)))
/ (a0 ** (-pls_temp - 3 / 2) * (-pls_temp + 1 / 2))
+
((sol.a_out ** (-pls_temp)) - (a0 ** (-pls_temp)))
/ (a0 ** (-pls_temp - 2) * (-pls_temp))))
f_st: np.ndarray = sigma0 * ((1. / a0) ** (-pls_temp - 3 / 2)) * sol.x
idx: np.ndarray = np.argwhere(sol.r > a0)[:, 0]
f_st[idx] = sigma0 * ((1. / a0) ** (-pls_temp - 2.))
return f_st
[docs]
def radial_speed(sol: Sol, ys: np.ndarray, t: np.float64) -> np.ndarray:
"""
Calculate the radial speed for a given sigma, sol and m_out
Parameters
----------
sol : Sol
Associated Sol
ys : np.ndarray
sigma * r ** (2 + pls_temp)
t : np.float64
time
Returns
-------
np.ndarray
the radial speed
"""
nu0: np.float64 = (Rgp * T0(t, sol.a0, sol.Ms) * sol.alpha * (sol.a0 ** (3 / 2))
/ (sol.mu.mean() * np.sqrt(G * sol.Ms)))
D: np.float64 = (3 * nu0 / (4. * (sol.a0 ** (3 / 2 + sol.pls_temp))))
return -2. * D * (sol.x ** (1 + sol.pls_temp)) * np.gradient(np.log(ys), np.log(sol.x))
[docs]
def mass_flux(sol: Sol, ys: np.ndarray, t: np.float64) -> np.ndarray:
"""
Calculate the total mass flux for a given sigma, sol and m_out
Parameters
----------
sol : Sol
Associated Sol
ys : np.ndarray
sigma * r ** (2 + pls_temp)
t : np.float64
time
Returns
-------
np.ndarray
the radial speed
"""
nu0: np.float64 = (Rgp * T0(t, sol.a0, sol.Ms) * sol.alpha * (sol.a0 ** (3 / 2))
/ (sol.mu.mean() * np.sqrt(G * sol.Ms)))
D: np.float64 = (3 * nu0 / (4. * (sol.a0 ** (3 / 2 + sol.pls_temp))))
return - 4. * Pi * D * np.gradient(ys, sol.x)
# return -4. * Pi * D * gradient(ys, sol.x[1]-sol.x[0])
[docs]
def turbulent_speed(y: np.ndarray, sol: Sol, t: np.float64) -> np.ndarray:
"""
Calculate the turbulent speed for a given sigma, sol and m_out
Parameters
----------
y : np.ndarray
mass fractions at time t
sol : Sol
Associated Sol
t : np.float64
time
Returns
-------
np.ndarray
the radial speed
"""
nu0: np.float64 = (Rgp * T0(t, sol.a0, sol.Ms) * sol.alpha * (sol.a0 ** (3 / 2))
/ (sol.mu.mean() * np.sqrt(G * sol.Ms)))
return - (nu0 / (sol.a0 ** (3/2 + sol.pls_temp))) * (sol.x ** (4 + 2 * sol.pls_temp)) * np.gradient(y, sol.x)
[docs]
def T0(t: np.float64 = -1, a0: np.float64 = a0, Ms: np.float64 = Ms
) -> np.ndarray | np.float64:
"""
Calculate and return the temperature at a0
Parameters
----------
t : np.ndarray
The current time
a0 : np.float64
The distance to the central star at wich the temperature is calculated
Ms : np.float64
The mass of the central star
Returns
-------
np.ndarray | np.float
The temperature
"""
if t == -1:
return 278 * ((Ls / Lsun) ** (1 / 4)) * ((a0 / au) ** (-1 / 2))
else:
return 278 * ((interpLbol(t, Ms) / Lsun) ** (1 / 4)) * ((a0 / (au)) ** (-1 / 2))
[docs]
def T(r: np.ndarray | np.float64, t: np.float64 = -1, Ms: np.float64 = Ms) -> np.ndarray | np.float64:
"""
Calculate and return the temperature at radius r
Parameters
----------
r : np.ndarray
radius
t : np.ndarray
The current time
Returns
-------
np.ndarray | np.float
The temperature
"""
if t == -1 or interpLbol is None:
T0 = (278 * ((Ls / Lsun) ** (1 / 4)) * ((a0 / au) ** (-1 / 2)))
else:
T0 = (278 * ((interpLbol(t, Ms) / Lsun) ** (1 / 4)) * ((a0 / au) ** (-1 / 2)))
return T0 * (abs(r) / a0) ** pls_temp
# Shielding functions for CO
tab_Theta: np.ndarray = np.array([1, 1., 9.405e-1, 7.046e-1, 4.015e-1, 9.964e-2, 1.567e-2, 3.162e-3, 4.839e-4, 1e-100])
tab_log_nCO: np.ndarray = np.array([-500, 1., 13, 14, 15, 16, 17, 18, 19, 500]) - np.log10(1e-4 * C_L * C_L)
extrapol_Theta = interp1d(x=tab_log_nCO, y=tab_Theta, fill_value="extrapolate")
sigma_crit_CO: np.float64 = np.double(1e-7) * Mearth / (au * au)
tph_co = np.double(120.) * yr
# Photodissociation timescale for H2O
tph_h2o: np.float64 = 1 / trapezoid(photo_diss["I"] * photo_diss["cs"], photo_diss["lambdas_nm"])
cross_section_h2o: np.float64 = np.double(5.e-22) * C_L * C_L # mean cross-section
sigma_crit_h2o: np.float64 = np.double(1. / ((Na / (18 * 1e-2 * C_M)) * cross_section_h2o))
# tph_h2o: np.float64 = yr
[docs]
def sigma_dot(t: np.float64, ys: np.ndarray, phi: np.ndarray, sol: Sol) \
-> np.ndarray:
"""
Calculate the gas generation and destruction rate
Parameters
----------
t : np.float64
the current time
ys : np.ndarray
sigma * (r**(2 + alpha_t))
phi : np.ndarray
mass flux
sol : Sol
Object Sol associated with the resolution
Returns
-------
np.ndarray
sigma_dot
"""
ydot = np.zeros_like(ys)
# Belt
if "H2O" in sol.nom_mu:
ydot[sol.idxa] += np.maximum(sol.sol_ast.interp_sig_dot(t + sol.t0), zero)
if "CO" in sol.nom_mu:
ydot[sol.idxk] += np.maximum(sol.sol_kuip.interp_sig_dot(t + sol.t0), zero)
# Accretion by planets
for (r_planet, m_planet) in zip(sol.a_planets, sol.m_planets):
ip = np.argwhere(sol.r >= r_planet)[0, 0] # the index of the planet's position (self.r[ip]=a_planet)
cd: np.float64 = np.double(np.sqrt(Rgp * T(r_planet, t + sol.t0) / np.mean(sol.mu))) # sound speed in the disc
H: np.float64 = np.sqrt(r_planet * r_planet * r_planet * Rgp * T(r_planet + sol.t0) / (G * Ms * sol.mu.mean()))
# Vertical scale high
RH: np.float64 = r_planet * ((m_planet / (3. * sol.Ms)) ** (1 / 3.)) # Hill radius
alpha_accr: np.float64 = sol.f_accr * min(np.double(1.), RH / H) # Final accretion efficiency
# alpha_accr: np.float64 = sol.f_accr
idx: np.ndarray = ip + np.arange(-sol.dip, sol.dip)
mdot_in: np.float64 = phi[ip + sol.dip]
mdot_out: np.float64 = (phi[ip - sol.dip - 1] + phi[ip - sol.dip - 2]) / 2.
idx_b: int = idx.max()
idx_b2: int = idx.min()
delta_r: np.float64 = sol.r[idx_b] - sol.r[idx_b2]
sigmar: np.float64 = np.double(delta_r / 2.)
fmod = np.exp(- (sol.r[idx] - r_planet) * (sol.r[idx] - r_planet) / (
2. * sigmar * sigmar))
norm_fmod = trapezoid(2. * Pi * sol.r[idx] * fmod, sol.r[idx])
if mdot_out * mdot_in > 0:
# ydot[idx] -= ((fmod / norm_fmod) * alpha_accr * (min(max((1. / (1 - alpha_accr)) * abs(mdot_out),
# abs(mdot_in)), zero))
# * np.exp(-ys[idx]))
# ydot[idx] -= ((fmod / norm_fmod) * alpha_accr * abs(mdot_in)
# * np.exp(-ys[idx]))
ydot[idx] -= (fmod / norm_fmod) * (alpha_accr / (1 - alpha_accr)) * abs(mdot_out)
return ydot
[docs]
def mass_frac_dot(t: np.float64, ys: np.ndarray, sol: Sol) \
-> np.ndarray:
"""
Calculate the gas generation and destruction rate relatively to each species,
The accretion onto planet is not expected to have any influence onto the mass fraction
Parameters
----------
t : np.float64
the current time
ys : np.ndarray
The mass fractions shape=(n_species, len(self.r))
sol : Sol
Object Sol associated with the resolution
Returns
-------
np.ndarray
"""
ydot = np.zeros_like(ys)
# Belt
surf_dens: np.ndarray = sol.interp_f(t) * (sol.x ** (-4 - 2 * sol.pls_temp))
if "H2O" in sol.nom_mu:
sigmap: np.ndarray = np.maximum(sol.sol_ast.interp_sig_dot(t + sol.t0), zero)
idxa = np.intersect1d(np.argwhere(sol.x > np.sqrt(sol.sol_ast_params["a0"]
- sol.sol_ast_params["delta_a"]))[:, 0],
np.argwhere(sol.x < np.sqrt(sol.sol_ast_params["a0"]
+ sol.sol_ast_params["delta_a"]))[:, 0])
ydot[:, idxa] = -sigmap * ys[:, idxa] / surf_dens[idxa]
ydot[sol.idx_mu_ast[0], idxa] = sigmap * (1 - ys[sol.idx_mu_ast[0], idxa]) / surf_dens[idxa]
if "CO" in sol.nom_mu:
idxk = np.intersect1d(np.argwhere(sol.x > np.sqrt(sol.sol_kuip_params["a0"]
- sol.sol_kuip_params["delta_a"]))[:, 0],
np.argwhere(sol.x < np.sqrt(sol.sol_kuip_params["a0"]
+ sol.sol_kuip_params["delta_a"]))[:, 0])
sigmap: np.ndarray = np.maximum(sol.sol_kuip.interp_sig_dot(t + sol.t0), zero)
ydot[:, idxk] -= sigmap * ys[:, idxk] / surf_dens[idxk]
ydot[sol.idx_mu_kuip[0], idxk] = sigmap * (1 - ys[sol.idx_mu_kuip[0], idxk]) / surf_dens[idxk]
# Photo-dissociation
if interpLbol is not None:
Luv: np.float64 = interpLbol(t + sol.t0, Ms) * (interpReuv_3692(t + sol.t0, Ms) + interpReuv_lya(t + sol.t0, Ms)) / 2.
else:
Luv: np.float64 = zero
if "H2O" in sol.nom_mu:
idx = np.argwhere(ys[sol.idx_mu_ast[0]] > 0)[:, 0]
H: np.ndarray = np.sqrt(
(Rgp * T(sol.x[idx] * sol.x[idx], t + sol.t0)
/ sol.mu[0]) / (G * Ms / (sol.x[idx] ** 6)))
# H2O
if len(idx) > 2:
ys0_mean = ys[sol.idx_mu_ast[0], idx]
sig_dot = (ys0_mean * np.exp(- surf_dens[idx] * ys0_mean / sigma_crit_h2o) / tph_h2o)
sig_dot += (ys0_mean
* (1 / ((Luv * H / (2. * sol.x[idx] * sol.x[idx]))
/ (h * c / (100 * 1e-9 * C_L)) * cross_section_h2o))
* np.exp(-((cumulative_trapezoid(surf_dens[idx] * ys0_mean,
sol.x[idx] * sol.x[idx], initial=zero)
/ (2. * Pi * sol.x[idx] * sol.x[idx] * H)) / sigma_crit_h2o)))
sig_dot = np.maximum(np.minimum(sig_dot, ys[sol.idx_mu_ast[0], idx] / (6000. * yr)),
ys[sol.idx_mu_ast[0], idx] / (1e3 * Myr))
# sig_dot = ys[0, idx] / (100 * kyr)
ydot[sol.idx_mu_ast[0], idx] -= sig_dot
ydot[sol.idx_mu_ast[1], idx] += sig_dot * 2 * sol.mu[sol.idx_mu_ast[1]] / sol.mu[sol.idx_mu_ast[0]]
# H2O -> 2H + O
ydot[sol.idx_mu_ast[2], idx] += sig_dot * sol.mu[sol.idx_mu_ast[2]] / sol.mu[sol.idx_mu_ast[0]]
if "CO" in sol.nom_mu:
idx = np.argwhere(ys[sol.idx_mu_kuip[0]] > -80)[:, 0]
H: np.ndarray = np.sqrt(
(Rgp * T(sol.r[idx], t + sol.t0) / np.mean(sol.mu)) / (G * Ms / (sol.r[idx] ** 6)))
# cross_section_co: np.float64 = np.double(5.e-20) * C_L * C_L # mean cross section
# sigma_crit_CO: np.float64 = 1 / ((Na / sol.mu[sol.idx_mu_kuip[0]]) * cross_section_co)
# sigma_crit_CO: np.float64 = Mearth / (au * au)
if len(idx) > 2:
sig_dot = (ys[sol.idx_mu_kuip[0], idx] * np.exp(- surf_dens[idx] * ys[sol.idx_mu_kuip[1], idx]
/ sigma_crit_CO)
* extrapol_Theta(np.log10(surf_dens[idx] * ys[sol.idx_mu_kuip[0], idx] * Na / sol.mu[0]))
/ tph_co)
# sig_dot += (ys[sol.idx_mu_kuip[0], idx]
# * (1 / (Luv / (h * c / (100 * 1e-9 * C_L)) * cross_section_co))
# * np.exp(-((cumulative_trapezoid(2. * Pi * sol.r[idx]
# * surf_dens[idx] * ys[sol.idx_mu_kuip[0], idx],
# sol.r[idx], initial=zero)
# / (2. * Pi * sol.r[idx] * H)) / sigma_crit_CO)))
sig_dot = np.maximum(np.minimum(sig_dot, ys[sol.idx_mu_kuip[0], idx] / (6000. * yr)),
ys[sol.idx_mu_kuip[0], idx] / (1e3 * Myr))
# sig_dot = savgol_filter(sig_dot, 10,4)
# sig_dot = ys[0, idx] / (100 * kyr)
ydot[sol.idx_mu_kuip[0], idx] -= sig_dot
ydot[sol.idx_mu_kuip[1], idx] += sig_dot * sol.mu[sol.idx_mu_kuip[1]] / sol.mu[sol.idx_mu_kuip[0]]
ydot[sol.idx_mu_kuip[2], idx] += sig_dot * sol.mu[sol.idx_mu_kuip[2]] / sol.mu[sol.idx_mu_kuip[0]]
return ydot
[docs]
def boundary_conditions(ys: np.ndarray, sol: Sol) -> np.ndarray:
"""
Define the inner and outer boundary conditions for the surface density
Parameters
----------
ys : np.ndarray
log of surface density
Nu : np.ndarray
Viscosity
sol : Sol
Associated Sol
Returns
-------
np.ndarray
"""
# ys[0] = ys[1] + np.log(Nu[1] / Nu[0])
ys = np.maximum(ys, 1e-200)
# ys[0] = max(ys[1] * (sol.x[0] / sol.x[1]) ** (
# np.log(ys[1] / ys[2]) / np.log(sol.x[1] / sol.x[2])), 1e-200)
ys[0] = max(min(ys[1] * sol.x[0] / sol.x[1], ys[1] * (sol.x[0] / sol.x[1]) ** (
np.log(ys[1] / ys[2]) / np.log(sol.x[1] / sol.x[2])), ys[1]), 1e-200)
# ys[-1] = max(min(ys[-2] * sol.x[-1] / sol.x[-2], ys[-2] * (sol.x[-1] / sol.x[-2]) ** (
# np.log(ys[-2] / ys[-3]) / np.log(sol.x[-2] / sol.x[-3]))), 1e-200)
ys[-1] = max(min(ys[-2] * (sol.x[-1] / sol.x[-2]) ** (
np.log(ys[-2] / ys[-4]) / np.log(sol.x[-2] / sol.x[-4])), ys[-2]),
1e-200)
return ys
[docs]
def boundary_conditions_mf(ys: np.ndarray, sol: Sol) -> np.ndarray:
"""
Define the inner and outer boundary conditions for the mass fractions
Parameters
----------
ys : np.ndarray
Mass fractions
Nu : np.ndarray
Viscosity
sol : Sol
Associated Sol
Returns
-------
np.ndarray
"""
# ys = savgol_filter(ys, 5, 3)
# ys[1:] = (ys[1:] + ys[:-1]) / 2.
ys = np.minimum(np.maximum(1e-50, ys), 1.)
# if sol.sol_ast.sublimation_model != "none":
# ys[0, sol.idx_mu_ast[0]] = max(ys[0,sol.idx_mu_ast[0]], ys[1, sol.idx_mu_ast[0]]) # no inner H2O flux
# ys[0, sol.idx_mu_ast[1]] = - (2. * sol.mu[sol.idx_mu_ast[1]] / sol.idx_mu_ast[0]) * ys[0, sol.idx_mu_ast[0]]
# ys[0, sol.idx_mu_ast[2]] = - (sol.mu[sol.idx_mu_ast[2]] / sol.idx_mu_ast[0]) * ys[0, sol.idx_mu_ast[0]]
# ys[-1, sol.idx_mu_ast[0]] = min(ys[-1,sol.idx_mu_ast[0]], ys[-2, sol.idx_mu_ast[0]]) # no inner H2O flux
# ys[-1, sol.idx_mu_ast[1]] = - (2. * sol.mu[sol.idx_mu_ast[1]] / sol.idx_mu_ast[0]) * ys[-1, sol.idx_mu_ast[0]]
# ys[-1, sol.idx_mu_ast[2]] = - (sol.mu[sol.idx_mu_ast[2]] / sol.idx_mu_ast[0]) * ys[-1, sol.idx_mu_ast[0]]
# if sol.sol_kuip.sublimation_model != "none":
# ys[0, sol.idx_mu_kuip[0]] = max(ys[0,sol.idx_mu_kuip[0]], ys[1, sol.idx_mu_kuip[0]]) # no inner CO flux
# ys[0, sol.idx_mu_kuip[1]] = - (sol.mu[sol.idx_mu_kuip[1]] / sol.idx_mu_kuip[0]) * ys[0, sol.idx_mu_kuip[0]]
# ys[0, sol.idx_mu_kuip[2]] = - (sol.mu[sol.idx_mu_kuip[2]] / sol.idx_mu_kuip[0]) * ys[0, sol.idx_mu_kuip[0]]
# ys[-1, sol.idx_mu_kuip[0]] = min(ys[-1, sol.idx_mu_kuip[0]], ys[-2, sol.idx_mu_kuip[0]]) # no inner H2O flux
# ys[-1, sol.idx_mu_kuip[1]] = - (2. * sol.mu[sol.idx_mu_kuip[1]] / sol.idx_mu_kuip[0]) * ys[-1, sol.idx_mu_kuip[0]]
# ys[-1, sol.idx_mu_kuip[2]] = - (sol.mu[sol.idx_mu_kuip[2]] / sol.idx_mu_kuip[0]) * ys[-1, sol.idx_mu_kuip[0]]
# ys[:, 0] = np.minimum(np.maximum(
# np.minimum(ys[:, 1] * (sol.x[0] / sol.x[1]) ** (
# np.log(ys[:, 1] / ys[:, 3]) / np.log(sol.x[1] / sol.x[3])),
# ys[:, 1] * sol.x[0] / sol.x[1]), 1e-140),
# 1)
ys[:, 0] = ys[:, 1].copy()
ys[:, -1] = ys[:, -2].copy()
# ys[:, -1] = np.minimum(np.maximum(ys[:, -2] * (sol.x[-1] / sol.x[-2]) ** (
# np.log(ys[:, -2] / ys[:, -4]) / np.log(sol.x[-2] / sol.x[-4])), 1e-140),
# 1).copy()
return ys
[docs]
def system_sigma_tot(t: np.float64, y: np.array, sol: Sol, verbose: bool = True) -> np.ndarray:
"""
Define the system to integrate in time to get the evolution of surface density.
The spatial differentiation is made here to estimate the spatial derivatives
Parameters
----------
t : np.float64
The current time
y : np.ndarray
The log of surface density
sol : Sol
The associated Sol
verbose : bool
To print some information
Returns
-------
np.ndarray
d(log(sigma)/dt
"""
nu0: np.float64 = (Rgp * T0(t, sol.a0, sol.Ms) * sol.alpha * (sol.a0 ** (3 / 2))
/ (sol.mu.mean() * np.sqrt(G * sol.Ms)))
ys = boundary_conditions(y, sol)
D: np.float64 = (3 * nu0 / (4. * (sol.a0 ** (3 / 2 + sol.pls_temp))))
phi = mass_flux(sol, ys, t)
if sol.pls_temp == -0.5:
# yp: np.ndarray = (D * gradient(gradient(ys, sol.x[1]-sol.x[0]), sol.x[1]-sol.x[0])
# + (sol.x ** (4. + 2. * sol.pls_temp)) * sigma_dot(t=t, ys=ys, phi=phi, sol=sol))
yp: np.ndarray = (D * gradient2(ys, sol.x[1]-sol.x[0])
+ (sol.x ** (4. + 2. * sol.pls_temp)) * sigma_dot(t=t, ys=ys, phi=phi, sol=sol))
else:
# yp: np.ndarray = (D * (sol.r ** (-2 - sol.pls_temp))
# * gradient(gradient(ys, sol.x[1]-sol.x[0]), sol.x[1]-sol.x[0])
# + (sol.x ** (4. + 2. * sol.pls_temp)) * sigma_dot(t=t, ys=ys, phi=phi, sol=sol))
yp: np.ndarray = (D * (sol.r ** (-2 - sol.pls_temp))
* gradient2(ys, sol.x[1]-sol.x[0])
+ (sol.x ** (4. + 2. * sol.pls_temp)) * sigma_dot(t=t, ys=ys, phi=phi, sol=sol))
# if sol.pls_temp == -0.5:
# yp: np.ndarray = (D * np.gradient(np.gradient(ys, sol.x), sol.x)
# + (sol.x ** (4. + 2. * sol.pls_temp)) * sigma_dot(t=t, ys=ys, phi=phi, sol=sol))
# else:
# yp: np.ndarray = (D * (sol.r ** (-2 - sol.pls_temp)) * np.gradient(np.gradient(ys, sol.x), sol.x)
# + (sol.x ** (4. + 2. * sol.pls_temp)) * sigma_dot(t=t, ys=ys, phi=phi, sol=sol))
if verbose and sol.compteur % 1000 == 0:
# g.loglog(sol.r, [sigma_dot(t=t, ys=ys, phi=phi, sol=sol)],".")
if t < Myr:
print(sol.compteur, " t=", t / yr, " yr",
"Mbelt=", trapezoid(2 * Pi * sol.r[2:-2] * ys[2:-2] * (sol.r[2:-2] ** (-2. - sol.pls_temp)),
x=sol.r[2:-2]) / Mearth, " Mearth")
elif t < Gyr:
print(sol.compteur, " t=", t / Myr, " Myr",
"Mbelt=", trapezoid(2 * Pi * sol.r[2:-2] * ys[2:-2] * (sol.r[2:-2] ** (-2. - sol.pls_temp)),
x=sol.r[2:-2]) / Mearth, " Mearth")
else:
print(sol.compteur, " t=", t / Gyr, " Gyr",
"Mbelt=", trapezoid(2 * Pi * sol.r[2:-2] * ys[2:-2] * (sol.r[2:-2] ** (-2. - sol.pls_temp)),
x=sol.r[2:-2]) / Mearth, " Mearth")
print("mdot",
trapezoid(2 * Pi * sol.r[2:-2] * yp[2:-2] * (sol.r[2:-2] ** (-2. - sol.pls_temp)),
sol.r[2:-2]) * (Myr / Mearth),
"Mearth / Myr")
print("mdot integrated sigmap : ",
trapezoid(2. * Pi * sol.r[2:-2] * sigma_dot(t=t, ys=ys, phi=phi, sol=sol)[2:-2],
sol.r[2:-2]) * Myr / Mearth, "Mearth / Myr")
print("flux int : ", phi[0] * Myr / Mearth, "Mearth / Myr")
print("flux ext : ", phi[-1] * Myr / Mearth,"Mearth / Myr")
if verbose:
sol.compteur += 1
return yp
[docs]
def system_mass_fraction(t: np.float64, y: np.array, sol: Sol, verbose: bool = True) -> np.ndarray:
"""
Define the system to integrate in time to get the evolution of mass fractions,
the spatial differentiation is made here to estimate the spatial derivatives
Parameters
----------
t : np.float64
The current time
y : np.ndarray
The one dimensional array of mass fractions (shape=(n_species*len(self.r)))
sol : Sol
The associated Sol
verbose : bool, optional, default=True
To print some information
Returns
-------
np.ndarray
d(log(sigma)/dt
"""
# ys = boundary_conditions_mf(sol.split_y(y), sol)
ys = boundary_conditions_mf(split_y(y, len(sol.x), len(sol.mu)), sol)
nu0: np.float64 = (Rgp * T0(t, sol.a0, sol.Ms) * sol.alpha * (sol.a0 ** (3 / 2))
/ (sol.mu.mean() * np.sqrt(G * sol.Ms)))
D: np.float64 = (3 * nu0 / (4. * (sol.a0 ** (3 / 2 + sol.pls_temp))))
f = savgol_filter(sol.interp_f(t), 10, 3)
# d2: np.ndarray = np.gradient(np.gradient(ys, sol.x, axis=1), sol.x, axis=1)
# d2: np.ndarray = gradient(gradient(ys, sol.x[1]-sol.x[0], axis=1), sol.x[1]-sol.x[0], axis=1)
d2: np.ndarray = gradient2(ys, sol.x[1]-sol.x[0], axis=1)
# d2[:5] = np.minimum(d2[:5], zero)
# d2[-5:] = np.maximum(d2[-5:], zero)
# d2[:5] = np.maximum(d2[:5], zero)
# d2[-5:] = np.minimum(d2[-5:], zero)
d2[:5] = zero
d2[-5:] = zero
yp: np.ndarray[np.float64] = ((D / 3.) * (sol.x ** (1 + 2 * sol.pls_temp))
* (d2
+ 4. * (gradient(f, sol.x[1]-sol.x[0]) / f)
* gradient(ys, sol.x[1]-sol.x[0], axis=1, radial_speed=gradient(f, sol.x[1]-sol.x[0])))
+ mass_frac_dot(t=t, ys=ys, sol=sol))
# yp: np.ndarray[np.float64] = ((D / 3.) * (sol.x ** (1 + 2 * sol.pls_temp))
# * (d2
# + 4. * (np.gradient(f, sol.x) / f)
# * np.gradient(ys, sol.x, axis=1))
# + mass_frac_dot(t=t, ys=ys, sol=sol))
yp[:, 0] = zero
yp[:, -1] = zero
if verbose and sol.compteur % 1000 == 0:
print("yp=", trapezoid(2. * Pi * (sol.x ** 2) * yp, sol.x ** 2, axis=1))
if t < Myr:
print(sol.compteur, " t=", t / yr, " yr")
elif t < Gyr:
print(sol.compteur, " t=", t / Myr, " Myr")
else:
print(sol.compteur, " t=", t / Gyr, " Gyr")
# print("Masses fractions = ", trapezoid(2. * Pi * sol.r * sigma * ys, sol.r, axis=1)
# / trapezoid(2. * Pi * sol.r * sigma, sol.r))
sol.compteur += 1
return merge_y(yp)