From 1f90761c84c5f596ef7f413134dac719c6489edb Mon Sep 17 00:00:00 2001 From: Aaron David Schneider Date: Wed, 22 Nov 2023 07:49:06 +0100 Subject: [PATCH] black formating --- chemcomp/accretion_models/__init__.py | 2 +- chemcomp/accretion_models/_accretion_class.py | 8 +- chemcomp/accretion_models/gas.py | 40 +- chemcomp/accretion_models/pebbles.py | 70 +- chemcomp/disks/_0pacity/opacity_models.py | 59 +- chemcomp/disks/_chemistry.py | 358 +++++---- chemcomp/disks/_disk_class.py | 691 ++++++++++++------ chemcomp/disks/kees.py | 73 +- chemcomp/disks/mmsn.py | 53 +- chemcomp/disks/twopop.py | 96 +-- chemcomp/helpers/__init__.py | 14 +- chemcomp/helpers/analysis_helper.py | 637 ++++++++++++---- chemcomp/helpers/import_config.py | 17 +- chemcomp/helpers/main_helpers.py | 85 ++- chemcomp/helpers/solvediffonedee.py | 53 +- chemcomp/helpers/tridag.py | 15 +- chemcomp/helpers/units/chemistry_const.py | 6 +- chemcomp/planets/_planet_class.py | 519 +++++++++---- chemcomp/planets/bertmigration.py | 267 ++++--- chemcomp/planets/migration_1.py | 15 +- chemcomp/planets/no_accretion.py | 20 +- chemcomp/planets/simple_planet.py | 24 +- 22 files changed, 2138 insertions(+), 984 deletions(-) diff --git a/chemcomp/accretion_models/__init__.py b/chemcomp/accretion_models/__init__.py index e594695..ed1cf95 100755 --- a/chemcomp/accretion_models/__init__.py +++ b/chemcomp/accretion_models/__init__.py @@ -1,3 +1,3 @@ from ._accretion_class import Accretion from .gas import GasAccretion -from .pebbles import PebbleAccretion \ No newline at end of file +from .pebbles import PebbleAccretion diff --git a/chemcomp/accretion_models/_accretion_class.py b/chemcomp/accretion_models/_accretion_class.py index af17ae6..b6cd5b1 100755 --- a/chemcomp/accretion_models/_accretion_class.py +++ b/chemcomp/accretion_models/_accretion_class.py @@ -13,11 +13,11 @@ def __init__(self, config=None): self.regime = "" def calc_m_dot(self): - """ function that calculates and updates the accretion model """ + """function that calculates and updates the accretion model""" pass def accretion_domain(self): - """ function that calculates the regimes of the accretion model """ + """function that calculates the regimes of the accretion model""" pass def init_accretion(self, planet): @@ -32,11 +32,11 @@ def init_accretion(self, planet): self._init_params() def _init_params(self): - """ function that inits parameters needed """ + """function that inits parameters needed""" pass def update_z(self): - """update the total mass of heavy elements that have been accreted. """ + """update the total mass of heavy elements that have been accreted.""" self.M_z[0] += np.sum(self.m_dot_c_chem[0, :-7] * self.planet.h) self.M_z[1] += np.sum(self.m_dot_a_chem[0, :-7] * self.planet.h) diff --git a/chemcomp/accretion_models/gas.py b/chemcomp/accretion_models/gas.py index 27e1fbc..9e9039f 100755 --- a/chemcomp/accretion_models/gas.py +++ b/chemcomp/accretion_models/gas.py @@ -10,14 +10,14 @@ class GasAccretion(Accretion): - """ Model for GasAccretion. + """Model for GasAccretion. Minimum of Ikoma2000, Machida2010 and disk gas accretion rate (including horseshoe region). """ def __init__(self, config): super().__init__() self.name = "GasAccretion" - self.kappa_env = config.get("kappa_env", 0.05 * u.cm ** 2 / u.g) + self.kappa_env = config.get("kappa_env", 0.05 * u.cm**2 / u.g) self.kappa_env = self.kappa_env.cgs.value self.f = config.get("f", 0.2) self.f_machida = config.get("f_machida", 1) @@ -25,10 +25,7 @@ def __init__(self, config): def _init_params(self): self.const_terms = ( - 1e3 - * (30*earth_mass_cgs_value)**2.5 - * self.kappa_env/0.05 - * year + 1e3 * (30 * earth_mass_cgs_value) ** 2.5 * self.kappa_env / 0.05 * year ) def accretion_domain(self): @@ -41,23 +38,23 @@ def remove_mass_from_flux(self): return def get_min(self): - """" Get the minimal accretion rate """ + """ " Get the minimal accretion rate""" low = ( - 0.83 - * self.f_machida - * self.planet.omega_k - * self.planet.H_g ** 2 - * self.planet.sigma_g - * (self.planet.r_h / self.planet.H_g) ** (9 / 2) + 0.83 + * self.f_machida + * self.planet.omega_k + * self.planet.H_g**2 + * self.planet.sigma_g + * (self.planet.r_h / self.planet.H_g) ** (9 / 2) ) low = np.maximum(low, 1e-60) high = ( - 0.14 - * self.f_machida - * self.planet.omega_k - * self.planet.H_g ** 2 - * self.planet.sigma_g + 0.14 + * self.f_machida + * self.planet.omega_k + * self.planet.H_g**2 + * self.planet.sigma_g ) high = np.maximum(high, 1e-60) @@ -69,7 +66,7 @@ def get_min(self): mdot_HS = np.maximum(mdot_HS, 1e-60) disk_max = disk_max + mdot_HS - ikoma = self.planet.M/(self.planet.M_c**(-2.5)*self.const_terms) + ikoma = self.planet.M / (self.planet.M_c ** (-2.5) * self.const_terms) ikoma = np.maximum(ikoma, 1e-60) valid_accrates = [low, high, disk_max, ikoma] @@ -80,7 +77,7 @@ def get_min(self): return str(minimum), mdot def calc_m_dot(self): - """ main routine that calculates the accretion rate.""" + """main routine that calculates the accretion rate.""" self.update_parameters() self.accretion_domain() self.m_dot_c = 0 @@ -96,6 +93,5 @@ def calc_m_dot(self): self.m_dot = self.m_dot_a + self.m_dot_c self.m_dot_a_chem = self.chem_mask_array * self.m_dot_a - assert np.all(self.m_dot_a_chem >= 0) - return \ No newline at end of file + return diff --git a/chemcomp/accretion_models/pebbles.py b/chemcomp/accretion_models/pebbles.py index 29b12f1..8646859 100755 --- a/chemcomp/accretion_models/pebbles.py +++ b/chemcomp/accretion_models/pebbles.py @@ -8,7 +8,7 @@ class PebbleAccretion(Accretion): - """ Model for PebbleAccretion. """ + """Model for PebbleAccretion.""" def __init__(self, config): super().__init__() @@ -18,7 +18,7 @@ def __init__(self, config): self.u_frag = config.get("u_frag", None) self.t_f = config.get("STOKES", None) - self.factor_peb_accretion = config.get('factor_peb_accretion', 1.0) + self.factor_peb_accretion = config.get("factor_peb_accretion", 1.0) if self.rho_solid is None: self._calculate_rho_solid = True @@ -28,7 +28,9 @@ def __init__(self, config): if self.u_frag is None: self._calculate_u_frag = True - print("Warning: No fixed fragmentation velocity has been set! The code will calculate the fragmentation velocity according to the Drazkowska 2017 approach: 1m/s if silicate, 10m/s if icy.") + print( + "Warning: No fixed fragmentation velocity has been set! The code will calculate the fragmentation velocity according to the Drazkowska 2017 approach: 1m/s if silicate, 10m/s if icy." + ) else: self.u_frag = self.u_frag.cgs.value self._calculate_u_frag = False @@ -37,34 +39,42 @@ def __init__(self, config): self.twoD = False def _init_params(self): - """ Pass some pebble parameters to the disk. """ - self.planet.disk.init_pebble_parameters(self.epsilon_p, self.rho_solid, self.t_f, self.u_frag, - self._calculate_rho_solid, self._calculate_u_frag) + """Pass some pebble parameters to the disk.""" + self.planet.disk.init_pebble_parameters( + self.epsilon_p, + self.rho_solid, + self.t_f, + self.u_frag, + self._calculate_rho_solid, + self._calculate_u_frag, + ) def calc_stokes(self): - """ Compute the stokes number by interpolating the stokes number of the disk to the planets location. """ - self.t_f = np.interp(self.planet.a_p, self.planet.disk.r, self.planet.disk.stokes_number_pebbles) + """Compute the stokes number by interpolating the stokes number of the disk to the planets location.""" + self.t_f = np.interp( + self.planet.a_p, self.planet.disk.r, self.planet.disk.stokes_number_pebbles + ) self.tau_f = self.t_f / self.planet.omega_k def calc_densities(self): - """ calculate the pebble scalehight and the volume density of pebbles. """ + """calculate the pebble scalehight and the volume density of pebbles.""" self.H_p = np.sqrt(self.planet.disk.alpha_height / self.t_f) * self.planet.H_g - self.rho_peb_components = self.planet.sigma_peb_components / (np.sqrt(2 * np.pi) * self.H_p) + self.rho_peb_components = self.planet.sigma_peb_components / ( + np.sqrt(2 * np.pi) * self.H_p + ) return def accretion_domain(self): # checked - """ decide on the regimes of accretion: bondi vs hill """ + """decide on the regimes of accretion: bondi vs hill""" if self.regime == "Bondi": transition_mass = ( - np.sqrt(1 / 3) - * self.planet.del_v ** 3 - / (G * self.planet.omega_k) + np.sqrt(1 / 3) * self.planet.del_v**3 / (G * self.planet.omega_k) ) if self.planet.M > transition_mass: self.regime = "Hill" def twoD_transition(self): # checked - """ Transition between 2D and 3D pebble accretion. """ + """Transition between 2D and 3D pebble accretion.""" transition = np.pi * self.R_acc / (2 * np.sqrt(2 * np.pi)) if transition > self.H_p: self.twoD = True @@ -72,7 +82,7 @@ def twoD_transition(self): # checked self.twoD = False def calc_R_acc(self): # checked - """Calculate the capture radius of pebble accretion. """ + """Calculate the capture radius of pebble accretion.""" if self.regime == "Bondi": t_b = self.planet.r_b / self.planet.del_v # checked self.R_acc = np.sqrt(4 * self.tau_f / t_b) * self.planet.r_b # checked @@ -83,17 +93,22 @@ def calc_R_acc(self): # checked if self.regime in ["Hill", "Bondi"]: t_p = ( - G - * self.planet.M - / ((np.abs(self.planet.del_v) + self.planet.omega_k * self.planet.r_h) ** 3) + G + * self.planet.M + / ( + (np.abs(self.planet.del_v) + self.planet.omega_k * self.planet.r_h) + ** 3 + ) + ) # checked + self.R_acc = self.R_acc * np.exp( + -0.4 * (self.tau_f / t_p) ** 0.65 ) # checked - self.R_acc = self.R_acc * np.exp(-0.4 * (self.tau_f / t_p) ** 0.65) # checked self.del_v = self.planet.del_v + self.planet.omega_k * self.R_acc # checked return def remove_mass_from_flux(self): - """ Ask the disk to remove the accreted pebbles.""" + """Ask the disk to remove the accreted pebbles.""" self.planet.disk.remove_from_peb(self.m_dot_chem, self.planet.idx) return @@ -107,13 +122,19 @@ def calc_m_dot(self): self.twoD_transition() if self.twoD: - self.m_dot_chem = 2 * self.R_acc * self.planet.sigma_peb_components * self.del_v # checked + self.m_dot_chem = ( + 2 * self.R_acc * self.planet.sigma_peb_components * self.del_v + ) # checked else: - self.m_dot_chem = np.pi * self.R_acc ** 2 * self.rho_peb_components * self.del_v # checked + self.m_dot_chem = ( + np.pi * self.R_acc**2 * self.rho_peb_components * self.del_v + ) # checked self.m_dot_chem = self.factor_peb_accretion * self.m_dot_chem - self.m_dot_c_chem = (1 - self.planet.solid_frac) * self.m_dot_chem # checked + self.m_dot_c_chem = ( + 1 - self.planet.solid_frac + ) * self.m_dot_chem # checked self.m_dot_a_chem = self.planet.solid_frac * self.m_dot_chem # checked self.m_dot = np.sum(self.m_dot_chem[0]) # sum over all elements @@ -121,7 +142,6 @@ def calc_m_dot(self): self.m_dot_c = (1 - self.planet.solid_frac) * self.m_dot # checked self.m_dot_a = self.planet.solid_frac * self.m_dot # checked - else: self.m_dot = 0 self.m_dot_c = 1 * self.m_dot diff --git a/chemcomp/disks/_0pacity/opacity_models.py b/chemcomp/disks/_0pacity/opacity_models.py index 337026b..d4adc9b 100755 --- a/chemcomp/disks/_0pacity/opacity_models.py +++ b/chemcomp/disks/_0pacity/opacity_models.py @@ -2,10 +2,15 @@ from chemcomp.helpers.units.berts_units import * -opacity_array = np.genfromtxt(os.path.join(os.path.dirname(os.path.realpath(__file__)), "meanopac5.dat"), - usecols=(0, 1), skip_header=1).T +opacity_array = np.genfromtxt( + os.path.join(os.path.dirname(os.path.realpath(__file__)), "meanopac5.dat"), + usecols=(0, 1), + skip_header=1, +).T opacity_array_derivative = np.gradient(opacity_array[1], opacity_array[0]) -opacity_array_derivative_derivative = np.gradient(opacity_array_derivative, opacity_array[0]) +opacity_array_derivative_derivative = np.gradient( + opacity_array_derivative, opacity_array[0] +) def belllin(rho, temp, Z=0.01, onlygas=False): @@ -19,26 +24,30 @@ def belllin(rho, temp, Z=0.01, onlygas=False): onlygas If set, then do not include the dust part of the _0pacity """ if onlygas: - ki = np.array([1e-8,1e-36,1.5e20,0.348]) - a = np.array([0.6667,0.3333,1.,0.]) - b = np.array([3.,10.,-2.5,0.]) + ki = np.array([1e-8, 1e-36, 1.5e20, 0.348]) + a = np.array([0.6667, 0.3333, 1.0, 0.0]) + b = np.array([3.0, 10.0, -2.5, 0.0]) else: - ki = np.array([2e-4,2e16,0.1,2e81,1e-8,1e-36,1.5e20,0.348]) - a = np.array([0.,0.,0.,1.,0.6667,0.3333,1.,0.]) - b = np.array([2.,-7.,0.5,-24.,3.,10.,-2.5,0.]) - nn = len(ki) + ki = np.array([2e-4, 2e16, 0.1, 2e81, 1e-8, 1e-36, 1.5e20, 0.348]) + a = np.array([0.0, 0.0, 0.0, 1.0, 0.6667, 0.3333, 1.0, 0.0]) + b = np.array([2.0, -7.0, 0.5, -24.0, 3.0, 10.0, -2.5, 0.0]) + nn = len(ki) - n = len(rho) - dm = np.zeros((nn,n)) + n = len(rho) + dm = np.zeros((nn, n)) for ii in range(nn): - dm[ii,:] = ki[ii]*rho[:]**a[ii] - tc = np.zeros((nn+1,n)) - for ii in range(nn-1): - tc[ii+1,:] = (dm[ii,:]/dm[ii+1,:])**(1./(b[ii+1]-b[ii])) - tc[nn,:] = 1e99 + dm[ii, :] = ki[ii] * rho[:] ** a[ii] + tc = np.zeros((nn + 1, n)) + for ii in range(nn - 1): + tc[ii + 1, :] = (dm[ii, :] / dm[ii + 1, :]) ** (1.0 / (b[ii + 1] - b[ii])) + tc[nn, :] = 1e99 kappa = np.zeros_like(rho) for ii in range(nn): - kappa = np.where(np.logical_and(temp>tc[ii,:],temp tc[ii, :], temp < tc[ii + 1, :]), + dm[ii] * temp ** b[ii], + kappa, + ) return kappa * Z / 0.01 @@ -50,7 +59,7 @@ def oplin(_Rho, _Temp, Z=0.01): """ ts4 = 1.0e-4 * _Temp - rho13 = _Rho ** 0.333333333 + rho13 = _Rho**0.333333333 rho23 = rho13 * rho13 ts42 = ts4 * ts4 ts44 = ts42 * ts42 @@ -59,9 +68,9 @@ def oplin(_Rho, _Temp, Z=0.01): # init with no scattering opacity = bk7 * _Rho / (ts42 * np.sqrt(ts4)) - m_0 = _Temp <= t456 * _Rho ** power2 - m_0_1 = np.logical_and(m_0, _Temp > t234 * _Rho ** power1) - m_1 = np.logical_or(_Temp < t678 * _Rho ** power3, _Rho < 1.0e-10) + m_0 = _Temp <= t456 * _Rho**power2 + m_0_1 = np.logical_and(m_0, _Temp > t234 * _Rho**power1) + m_1 = np.logical_or(_Temp < t678 * _Rho**power3, _Rho < 1.0e-10) # disjoint _0pacity laws for 5, 6, and 7. o5 = bk5 * rho23[m_1] * ts42[m_1] * ts4[m_1] @@ -82,8 +91,8 @@ def oplin(_Rho, _Temp, Z=0.01): o4 = bk4 * rho23[m_0_1] / (ts48[m_0_1] * ts4[m_0_1]) o5 = bk5 * rho23[m_0_1] * ts42[m_0_1] * ts4[m_0_1] # parameters used for smoothing - o4an = o4 ** 4 - o3an = o3 ** 4 + o4an = o4**4 + o3an = o3**4 # smoothed and continuous _0pacity law for regions 3, 4, and 5. opacity[m_0_1] = ( @@ -107,7 +116,7 @@ def oplin(_Rho, _Temp, Z=0.01): opacity[m_0] = ( (o1an * o2an / (o1an + o2an)) ** 2 + (o3 / (1 + 1.0e22 / t10)) ** 4 - ) ** 0.25 + ) ** 0.25 return opacity * Z / 0.01 diff --git a/chemcomp/disks/_chemistry.py b/chemcomp/disks/_chemistry.py index e9304a7..33162d7 100755 --- a/chemcomp/disks/_chemistry.py +++ b/chemcomp/disks/_chemistry.py @@ -74,87 +74,105 @@ def split_molecules(abundances_inp): """ ( - rest_mol, - TotMassCO, - TotMassN2, - TotMassCH4, - TotMassCO2, - TotMassNH3, - TotMassH2S, - TotMassH2O, - TotMassFe3O4, - TotMassC, - TotMassFeS, - TotMassNaAlSi3O8, - TotMassKAlSi3O8, - TotmassMg2SiO4, - TotMassFe2O3, - TotMassVO, - TotMassMgSiO3, - TotMassAl2O3, - TotMassTiO, - ) = abundances_inp + rest_mol, + TotMassCO, + TotMassN2, + TotMassCH4, + TotMassCO2, + TotMassNH3, + TotMassH2S, + TotMassH2O, + TotMassFe3O4, + TotMassC, + TotMassFeS, + TotMassNaAlSi3O8, + TotMassKAlSi3O8, + TotmassMg2SiO4, + TotMassFe2O3, + TotMassVO, + TotMassMgSiO3, + TotMassAl2O3, + TotMassTiO, + ) = abundances_inp TotC = ( - TotMassCO * MassC / MassCO - + TotMassCH4 * MassC / MassCH4 - + TotMassCO2 * MassC / MassCO2 - + TotMassC * MassC / MassC + TotMassCO * MassC / MassCO + + TotMassCH4 * MassC / MassCH4 + + TotMassCO2 * MassC / MassCO2 + + TotMassC * MassC / MassC ) TotO = ( - TotMassCO * MassO / MassCO - + TotMassCO2 * MassO * 2.0 / MassCO2 - + TotMassH2O * MassO / MassH2O - + TotMassFe3O4 * 4.0 * MassO / MassFe3O4 - + TotmassMg2SiO4 * 4.0 * MassO / MassMg2SiO4 - + TotMassFe2O3 * 3 * MassO / MassFe2O3 - + TotMassMgSiO3 * MassO * 3.0 / MassMgSiO3 - + TotMassTiO * MassO / MassTiO - + TotMassAl2O3 * 3 * MassO / MassAl2O3 - + TotMassKAlSi3O8 * 8 * MassO / MassKAlSi3O8 - + TotMassNaAlSi3O8 * 8 * MassO / MassNaAlSi3O8 - + TotMassVO * MassO / MassVO + TotMassCO * MassO / MassCO + + TotMassCO2 * MassO * 2.0 / MassCO2 + + TotMassH2O * MassO / MassH2O + + TotMassFe3O4 * 4.0 * MassO / MassFe3O4 + + TotmassMg2SiO4 * 4.0 * MassO / MassMg2SiO4 + + TotMassFe2O3 * 3 * MassO / MassFe2O3 + + TotMassMgSiO3 * MassO * 3.0 / MassMgSiO3 + + TotMassTiO * MassO / MassTiO + + TotMassAl2O3 * 3 * MassO / MassAl2O3 + + TotMassKAlSi3O8 * 8 * MassO / MassKAlSi3O8 + + TotMassNaAlSi3O8 * 8 * MassO / MassNaAlSi3O8 + + TotMassVO * MassO / MassVO ) TotFe = ( - TotMassFe3O4 * 3.0 * MassFe / MassFe3O4 - + TotMassFe2O3 * 2.0 * MassFe / MassFe2O3 - + TotMassFeS * MassFe / MassFeS + TotMassFe3O4 * 3.0 * MassFe / MassFe3O4 + + TotMassFe2O3 * 2.0 * MassFe / MassFe2O3 + + TotMassFeS * MassFe / MassFeS ) TotS = TotMassFeS * MassS / MassFeS + TotMassH2S * MassS / MassH2S TotMg = ( - TotmassMg2SiO4 * MassMg * 2.0 / MassMg2SiO4 - + TotMassMgSiO3 * MassMg / MassMgSiO3 + TotmassMg2SiO4 * MassMg * 2.0 / MassMg2SiO4 + + TotMassMgSiO3 * MassMg / MassMgSiO3 ) TotSi = ( - TotmassMg2SiO4 * MassSi / MassMg2SiO4 - + TotMassMgSiO3 * MassSi / MassMgSiO3 - + TotMassKAlSi3O8 * 3 * MassSi / MassKAlSi3O8 - + TotMassNaAlSi3O8 * 3 * MassSi / MassNaAlSi3O8 - + TotmassMg2SiO4 * MassSi / MassMg2SiO4 + + TotMassMgSiO3 * MassSi / MassMgSiO3 + + TotMassKAlSi3O8 * 3 * MassSi / MassKAlSi3O8 + + TotMassNaAlSi3O8 * 3 * MassSi / MassNaAlSi3O8 ) TotK = TotMassKAlSi3O8 * MassK / MassKAlSi3O8 TotNa = TotMassNaAlSi3O8 * MassNa / MassNaAlSi3O8 TotAl = ( - TotMassKAlSi3O8 * MassAl / MassKAlSi3O8 - + TotMassNaAlSi3O8 * MassAl / MassNaAlSi3O8 - + TotMassAl2O3 * 2 * MassAl / MassAl2O3 + TotMassKAlSi3O8 * MassAl / MassKAlSi3O8 + + TotMassNaAlSi3O8 * MassAl / MassNaAlSi3O8 + + TotMassAl2O3 * 2 * MassAl / MassAl2O3 ) TotTi = TotMassTiO * MassTi / MassTiO TotV = TotMassVO * MassV / MassVO TotN = TotMassNH3 * MassN / MassNH3 + TotMassN2 * 2 * MassN / MassN2 TotH_mol = ( - TotMassCH4 * 4.0 * MassH / MassCH4 - + TotMassH2O * 2.0 * MassH / MassH2O - + TotMassNH3 * 3 * MassH / MassNH3 - + TotMassH2S * 2 * MassH / MassH2S + TotMassCH4 * 4.0 * MassH / MassCH4 + + TotMassH2O * 2.0 * MassH / MassH2O + + TotMassNH3 * 3 * MassH / MassNH3 + + TotMassH2S * 2 * MassH / MassH2S ) # Note how the last entry of the array is TotH_mol - return np.array([TotC, TotO, TotFe, TotS, TotMg, TotSi, TotNa, TotK, TotN, TotAl, TotTi, TotV, TotH_mol]) + return np.array( + [ + TotC, + TotO, + TotFe, + TotS, + TotMg, + TotSi, + TotNa, + TotK, + TotN, + TotAl, + TotTi, + TotV, + TotH_mol, + ] + ) + -def calc_composition(abundances_inp: np.array, solid: bool, He_Mbg: float) -> Tuple[np.array, np.array]: +def calc_composition( + abundances_inp: np.array, solid: bool, He_Mbg: float +) -> Tuple[np.array, np.array]: """ Parameters @@ -169,21 +187,23 @@ def calc_composition(abundances_inp: np.array, solid: bool, He_Mbg: float) -> Tu Nd numpy array: the total mass of Hydrogen in molecules per grid cell """ - + molecule_elements = split_molecules(abundances_inp) rest_mol = abundances_inp[0] rest = np.sum(molecule_elements[:-1], axis=0) TotH_mol = molecule_elements[-1] - if solid: # These are the solid components, getting the H, He abundance is trivial + if solid: # These are the solid components, getting the H, He abundance is trivial TotH = TotH_mol TotHe = np.zeros_like(rest) - else: # The gas components, where information about the Helium content is required and therefore the supplied He_Mbg argument is used + else: # The gas components, where information about the Helium content is required and therefore the supplied He_Mbg argument is used TotHe = He_Mbg * rest_mol - TotH = (1 - rest - TotHe) - + TotH = 1 - rest - TotHe + # Test about the He+H split - assert np.isclose(TotH + TotHe - TotH_mol, rest_mol).all(), "Background gas after split does not equal the total background gas from the molecular abundances" + assert np.isclose( + TotH + TotHe - TotH_mol, rest_mol + ).all(), "Background gas after split does not equal the total background gas from the molecular abundances" elements = np.array( [ @@ -198,15 +218,18 @@ def calc_composition(abundances_inp: np.array, solid: bool, He_Mbg: float) -> Tu ], ) - molecules = abundances_inp.copy() # Avoid unwanted side-effects when returning the supplied molecular abundances + molecules = ( + abundances_inp.copy() + ) # Avoid unwanted side-effects when returning the supplied molecular abundances return np.transpose(np.stack([elements, molecules]), (2, 0, 1)), TotH_mol class Chemistry: - """ Chemistry class. Please be careful when changing or adding molecular species, + """Chemistry class. Please be careful when changing or adding molecular species, some occurancies of the chemistry rely (hardcoded) on the order of the species! """ + def __init__(self, config): # Adundances in terms of H (Asplund et al. 2009) + variation from Chiara plot if not config.get("use_FeH", True): @@ -263,7 +286,9 @@ def __init__(self, config): self.abuFe2O3 = 0.25 * (self.FeHabu - 0.9 * self.SHabu) self.abuFe3O4 = (1.0 / 6.0) * (self.FeHabu - 0.9 * self.SHabu) self.abuFeS = 0.9 * self.SHabu - self.abuMgSiO3 = self.MgHabu - 2 * (self.MgHabu - (self.SiHabu - 3 * self.KHabu - 3 * self.NaHabu)) + self.abuMgSiO3 = self.MgHabu - 2 * ( + self.MgHabu - (self.SiHabu - 3 * self.KHabu - 3 * self.NaHabu) + ) self.abuMg2SiO4 = self.MgHabu - (self.SiHabu - 3 * self.KHabu - 3 * self.NaHabu) self.abuCO = float(config.get("CO_frac", 0.45)) * self.CHabu @@ -284,21 +309,22 @@ def __init__(self, config): self.abuAl2O3 = 0.5 * (self.AlHabu - (self.KHabu + self.NaHabu)) self.abuH2O = self.OHabu - ( - 3 * self.abuMgSiO3 - + 4 * self.abuMg2SiO4 - + self.abuCO - + 2 * self.abuCO2 - + 3 * self.abuFe2O3 - + 4 * self.abuFe3O4 - + self.abuTiO - + 3 * self.abuAl2O3 - + 8 * self.abuNaAlSi3O8 - + 8 * self.abuKAlSi3O8 - + self.abuVO + 3 * self.abuMgSiO3 + + 4 * self.abuMg2SiO4 + + self.abuCO + + 2 * self.abuCO2 + + 3 * self.abuFe2O3 + + 4 * self.abuFe3O4 + + self.abuTiO + + 3 * self.abuAl2O3 + + 8 * self.abuNaAlSi3O8 + + 8 * self.abuKAlSi3O8 + + self.abuVO ) - assert np.isclose(self.abuC + self.abuCO + self.abuCO2 + self.abuCH4, self.CHabu), \ - "Broken chem model, we loose or gain C" + assert np.isclose( + self.abuC + self.abuCO + self.abuCO2 + self.abuCH4, self.CHabu + ), "Broken chem model, we loose or gain C" self.abu_array = np.array( [ @@ -383,22 +409,28 @@ def get_position_of_ice(self, T: np.array) -> np.array: idx = np.minimum(idx, T.size - 2) return idx - - def calc_He_Mbg(self, abundances_gas: np.array, abundances_solid: np.array, nHe_nH: np.array, T: np.array) -> float: + + def calc_He_Mbg( + self, + abundances_gas: np.array, + abundances_solid: np.array, + nHe_nH: np.array, + T: np.array, + ) -> float: """ Calculates the mass ratio of Helium with respect to the background gas. - + Parameters ---------- abundances_gas: 9xN numpy array containing gas composition. Can be generated with Chemistry.abundances abundances_solid: 9xN numpy array containing solid composition. Can be generated with Chemistry.abundances nHe_nH: nHe/nH number density ratio of Helium and Hydrogen. This should match with your initial condition T: Temperature array - + Returns ------- A float describing the mass ratio of Helium w.r.t. the background gas - + """ TotH_solid = split_molecules(abundances_solid)[-1] molecule_elements = split_molecules(abundances_gas) @@ -406,36 +438,58 @@ def calc_He_Mbg(self, abundances_gas: np.array, abundances_solid: np.array, nHe_ rest = np.sum(molecule_elements[:-1], axis=0) rest_mol = abundances_gas[0] Z = self.get_solid_heavies(T) - Mdust_Mgas = Z*rest_mol # = Mmol_solid / Mbg * Mbg / Mgas = Mmol_solid/Mgas = Mdust/Mgas - TotH_solid_Mgas = TotH_solid*Mdust_Mgas # Total mass of Hydrogen in solids over the total gas mass - + Mdust_Mgas = ( + Z * rest_mol + ) # = Mmol_solid / Mbg * Mbg / Mgas = Mmol_solid/Mgas = Mdust/Mgas + TotH_solid_Mgas = ( + TotH_solid * Mdust_Mgas + ) # Total mass of Hydrogen in solids over the total gas mass + # This gives the total hydrogen mass in both gaseous and solid phase, normed to the gas mass, (Hgas+Hsolid)/gas TotH = (1 - rest + TotH_solid_Mgas) / (1 + nHe_nH * MassHe) # This gives Hgas/gas TotH_gas = TotH - TotH_solid_Mgas # This gives H/gas * He/H = He/gas TotHe = TotH * nHe_nH * MassHe - + # Some tests to check for some errors for the calculation of TotH and TotHe - assert np.isclose(1-rest-molecule_elements[-1], rest_mol).all(), "Background gas from the molecular abundance array does not match background gas after splitting the molecules" - assert np.isclose(TotH_gas + TotHe - TotH_mol_gas, rest_mol).all(), "Splitting the background gas into hydrogen and helium did not work" - assert np.isclose(TotHe / (TotH*MassHe), nHe_nH).all(), "After separating solid and gaseous H, [He/H] is no longer correct" - assert np.isclose(TotH_gas+TotHe+rest, 1).all(), "Adding up all gaseous components does not result in the total gas mass" + assert np.isclose( + 1 - rest - molecule_elements[-1], rest_mol + ).all(), "Background gas from the molecular abundance array does not match background gas after splitting the molecules" + assert np.isclose( + TotH_gas + TotHe - TotH_mol_gas, rest_mol + ).all(), "Splitting the background gas into hydrogen and helium did not work" + assert np.isclose( + TotHe / (TotH * MassHe), nHe_nH + ).all(), "After separating solid and gaseous H, [He/H] is no longer correct" + assert np.isclose( + TotH_gas + TotHe + rest, 1 + ).all(), ( + "Adding up all gaseous components does not result in the total gas mass" + ) - # It is expected that He_Mbg is constant in space, so it's just saved as a float assert np.isclose((TotHe / rest_mol).std(), 0), "He_Mbg is not constant." - He_Mbg = (TotHe / rest_mol)[0] + He_Mbg = (TotHe / rest_mol)[0] return He_Mbg def get_solid_composition(self, sigma_dust_mol: np.array, T: np.array) -> np.array: sigma_dust = np.sum(sigma_dust_mol, axis=1) abundances_solid = np.transpose( - np.divide(sigma_dust_mol, sigma_dust[:, np.newaxis], out=np.zeros_like(sigma_dust_mol), - where=np.logical_and(sigma_dust[:, np.newaxis] > 1.0e-60, sigma_dust_mol > 1.0e-60))) - chem_solid, _ = calc_composition(abundances_solid, solid=True, He_Mbg=self.He_Mbg) + np.divide( + sigma_dust_mol, + sigma_dust[:, np.newaxis], + out=np.zeros_like(sigma_dust_mol), + where=np.logical_and( + sigma_dust[:, np.newaxis] > 1.0e-60, sigma_dust_mol > 1.0e-60 + ), + ) + ) + chem_solid, _ = calc_composition( + abundances_solid, solid=True, He_Mbg=self.He_Mbg + ) if np.max(T) > np.max(T_array): - chem_solid[:np.argmax(T) + 1] = 0 + chem_solid[: np.argmax(T) + 1] = 0 return chem_solid @@ -460,29 +514,36 @@ def get_gas_composition(self, sigma_g_components_mol: np.array = None) -> np.arr sigma = np.sum(sigma_g_components_mol, axis=1) abundances_gas = np.transpose( - np.divide(sigma_g_components_mol, sigma[:, np.newaxis], out=np.zeros_like(sigma_g_components_mol), - where=sigma[:, np.newaxis] != 0)) + np.divide( + sigma_g_components_mol, + sigma[:, np.newaxis], + out=np.zeros_like(sigma_g_components_mol), + where=sigma[:, np.newaxis] != 0, + ) + ) chem_gas, _ = calc_composition(abundances_gas, solid=False, He_Mbg=self.He_Mbg) self.mu = sigma / np.sum(sigma_g_components_mol.T / self.M_array, axis=0) return chem_gas def set_z(self, Z0): - """ set the fraction of heavy elements in the disk """ + """set the fraction of heavy elements in the disk""" # The fraction of heavy !molecules! in either gas or solids with respect to the backgroundgas (H+He): self.dtg0 = Z0 # ratio of heavies from the sum over all molecular species: self._ZH = np.sum(self.abu_array[1:] * np.squeeze(self.M_array)[1:]) - rezHabu = MassH + self.HeHabu*MassHe + (self.elabu_init*elmasses).sum() - self._dtg0_of_chem_model = self._ZH/rezHabu + rezHabu = MassH + self.HeHabu * MassHe + (self.elabu_init * elmasses).sum() + self._dtg0_of_chem_model = self._ZH / rezHabu return def get_solid_heavies(self, T): - """ returns a modified dust to gas ratio according to the temperature of the disk """ + """returns a modified dust to gas ratio according to the temperature of the disk""" abu_array = np.ones_like(T) * self.abu_array[:, np.newaxis] - M_masked_solid = np.where(T < T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array)) + M_masked_solid = np.where( + T < T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array) + ) M_masked_solid[0, :] = 0 # exclude restmol Tot_M_solid = np.sum(M_masked_solid * abu_array, axis=0) + 1e-300 @@ -508,60 +569,81 @@ def get_composition(self, T: np.array) -> np.array: """ abundances_solid, abundances_gas = self.abundances(T) - self.mu = np.sum(abundances_gas, axis=0) / np.sum(abundances_gas / self.M_array, axis=0) + self.mu = np.sum(abundances_gas, axis=0) / np.sum( + abundances_gas / self.M_array, axis=0 + ) # massratio tests for the elemental composition: - assert np.all(np.isclose(np.sum(abundances_gas, axis=0), - np.ones_like(T))), "error in chemistry, total sum of gas abundancies is not 1" - assert np.all(np.isclose(np.sum(abundances_solid, axis=0), - np.sum(abundances_solid, axis=0) != 0, np.ones_like(T), - np.zeros_like(T))), "error in chemistry, total sum of solid abundancies is not 1" + assert np.all( + np.isclose(np.sum(abundances_gas, axis=0), np.ones_like(T)) + ), "error in chemistry, total sum of gas abundancies is not 1" + assert np.all( + np.isclose( + np.sum(abundances_solid, axis=0), + np.sum(abundances_solid, axis=0) != 0, + np.ones_like(T), + np.zeros_like(T), + ) + ), "error in chemistry, total sum of solid abundancies is not 1" self.He_Mbg = self.calc_He_Mbg(abundances_gas, abundances_solid, self.HeHabu, T) - chem_solid, _ = calc_composition(abundances_solid, solid=True, He_Mbg=self.He_Mbg) + chem_solid, _ = calc_composition( + abundances_solid, solid=True, He_Mbg=self.He_Mbg + ) chem_gas, _ = calc_composition(abundances_gas, solid=False, He_Mbg=self.He_Mbg) if np.max(T) > np.max(T_array): - chem_solid[:np.argmax(T) + 1] = 0 + chem_solid[: np.argmax(T) + 1] = 0 else: # massratio tests for the molecular composition of the solids: # only test, when there are solids! - assert np.all(np.isclose(np.sum(chem_solid[:, 0, :], axis=1), - np.ones_like(T))), "error in chemistry, total sum of solid abundancies is not 1" + assert np.all( + np.isclose(np.sum(chem_solid[:, 0, :], axis=1), np.ones_like(T)) + ), "error in chemistry, total sum of solid abundancies is not 1" # massratio tests for the molecular composition of the gas: - assert np.all(np.isclose(np.sum(chem_gas[:, 0, :], axis=1), - np.ones_like(T))), "error in chemistry, total sum of gas abundancies is not 1" + assert np.all( + np.isclose(np.sum(chem_gas[:, 0, :], axis=1), np.ones_like(T)) + ), "error in chemistry, total sum of gas abundancies is not 1" ############### # Sanity checks ############### TotH_mol_test = ( - self.abuCH4 * 4.0 - + self.abuH2O * 2.0 - + self.abuNH3 * 3 - + self.abuH2S * 2 + self.abuCH4 * 4.0 + self.abuH2O * 2.0 + self.abuNH3 * 3 + self.abuH2S * 2 ) abu_mol_test = (self.abu_array * self.M_array.squeeze()).sum() - TotH_mol_test abu_asp_test = (self.elabu_init * elmasses).sum() - assert np.isclose(abu_asp_test, abu_mol_test), "error in chemistry, total sum of molecular abundancies is not equal to the total sum of the X/H stellar vaules. Partioning is wrong." + assert np.isclose( + abu_asp_test, abu_mol_test + ), "error in chemistry, total sum of molecular abundancies is not equal to the total sum of the X/H stellar vaules. Partioning is wrong." # Test if the heavy elements are identical to the heavy molecules (gas): - _, abu_gas_test = self.abundances([np.max(T_array)+100]) - chem_gas_test, restH = calc_composition(abu_gas_test, solid=False, He_Mbg=self.He_Mbg) - elabus_heavy_chem_test = np.sum(chem_gas_test[0,0,:-7]) - assert np.isclose(chem_gas_test[0, 1, 1:].sum(), elabus_heavy_chem_test+restH), "error in chemistry, total sum of elemental mass abundancies is not equal to the total sum of molecular mass abundancies." + _, abu_gas_test = self.abundances([np.max(T_array) + 100]) + chem_gas_test, restH = calc_composition( + abu_gas_test, solid=False, He_Mbg=self.He_Mbg + ) + elabus_heavy_chem_test = np.sum(chem_gas_test[0, 0, :-7]) + assert np.isclose( + chem_gas_test[0, 1, 1:].sum(), elabus_heavy_chem_test + restH + ), "error in chemistry, total sum of elemental mass abundancies is not equal to the total sum of molecular mass abundancies." # Test validity of dust to gas ratio: - assert np.isclose(chem_gas_test[0,1,1:].sum()*(1+self.dtg0),self.dtg0), "Error in dust to gas ratio" + assert np.isclose( + chem_gas_test[0, 1, 1:].sum() * (1 + self.dtg0), self.dtg0 + ), "Error in dust to gas ratio" # Test if the heavy elements are identical to the heavy molecules (solids): if np.max(T) < np.max(T_array): # only test, when there are solids! abu_solid_test, _ = self.abundances([0]) - chem_solid_test, restH = calc_composition(abu_solid_test, solid=True, He_Mbg=self.He_Mbg) - elabus_heavy_chem_test = np.sum(chem_solid_test[0,0,:-7]) - assert np.isclose(chem_solid_test[0, 1, 1:].sum(), elabus_heavy_chem_test+restH), "error in chemistry, total sum of elemental mass abundancies is not equal to the total sum of molecular mass abundancies." + chem_solid_test, restH = calc_composition( + abu_solid_test, solid=True, He_Mbg=self.He_Mbg + ) + elabus_heavy_chem_test = np.sum(chem_solid_test[0, 0, :-7]) + assert np.isclose( + chem_solid_test[0, 1, 1:].sum(), elabus_heavy_chem_test + restH + ), "error in chemistry, total sum of elemental mass abundancies is not equal to the total sum of molecular mass abundancies." return chem_solid, chem_gas @@ -582,19 +664,27 @@ def abundances(self, T: np.array) -> Tuple[np.array, np.array]: """ abu_array = np.ones_like(T) * self.abu_array[:, np.newaxis] - M_masked_solid = np.where(T < T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array)) + M_masked_solid = np.where( + T < T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array) + ) M_masked_solid[0, :] = 0 # exclude restmol Tot_M_solid = np.sum(M_masked_solid * abu_array, axis=0) + 1e-300 - Masses_solid = np.where(T < T_array, abu_array * M_masked_solid / Tot_M_solid, - np.zeros_like(abu_array)) + Masses_solid = np.where( + T < T_array, + abu_array * M_masked_solid / Tot_M_solid, + np.zeros_like(abu_array), + ) Z = self.get_solid_heavies(T) - M_masked_gas = np.where(T >= T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array)) + M_masked_gas = np.where( + T >= T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array) + ) M_masked_gas[0, :] = 0 # exclude restmol Tot_M_gas = np.sum(M_masked_gas * abu_array, axis=0) + 1e-300 - Masses_gas = np.where(T >= T_array, abu_array * M_masked_gas / Tot_M_gas, - np.zeros_like(abu_array)) + Masses_gas = np.where( + T >= T_array, abu_array * M_masked_gas / Tot_M_gas, np.zeros_like(abu_array) + ) Masses_gas[0, :] = 1 / (1 + (self.dtg0 - Z)) Masses_gas[1:, :] = Masses_gas[1:, :] * (self.dtg0 - Z) / (1 + (self.dtg0 - Z)) diff --git a/chemcomp/disks/_disk_class.py b/chemcomp/disks/_disk_class.py index dc45bab..75cd1d9 100755 --- a/chemcomp/disks/_disk_class.py +++ b/chemcomp/disks/_disk_class.py @@ -17,17 +17,19 @@ AU = const.au.cgs.value ME = const.M_earth.cgs.value -FWHM_TO_SIGMA = 2*(2*np.log(2))**0.5 +FWHM_TO_SIGMA = 2 * (2 * np.log(2)) ** 0.5 + def gaussian(x, max, mu, FWHM): - """Gaussian helper function. """ + """Gaussian helper function.""" # prevent underflow by maximizing exponent to -1e2 (np.exp(-1e2)=3e-44) - exp = np.minimum(0.5*((x-mu)/FWHM*FWHM_TO_SIGMA)**2.0, 1e2) - return max*np.exp(-exp) + exp = np.minimum(0.5 * ((x - mu) / FWHM * FWHM_TO_SIGMA) ** 2.0, 1e2) + return max * np.exp(-exp) + def tvisc_fct(T, sigma_g, tirr, omega_k, alpha, Z, mu): """Calculate viscous heating and return the difference between our current guess of the temperature - to the last guess of the temperature """ + to the last guess of the temperature""" c_s = (k_B * T / (mu * m_p)) ** 0.5 h = c_s / omega_k rho_g = sigma_g / (np.sqrt(2 * pi) * h) @@ -36,28 +38,31 @@ def tvisc_fct(T, sigma_g, tirr, omega_k, alpha, Z, mu): opacity = opacity_fct(rho_g, T, Z=Z) * 100.0 * Z - qvisc = (9.0 / 4.0) * sigma_g * nu * omega_k ** 2.0 + qvisc = (9.0 / 4.0) * sigma_g * nu * omega_k**2.0 tau = sigma_g * opacity # both sides of the disk - tvisc4 = qvisc * 3.0 * tau / (16.0 * sr) # factor 0.5 because tau is defined on complete disk + tvisc4 = ( + qvisc * 3.0 * tau / (16.0 * sr) + ) # factor 0.5 because tau is defined on complete disk - T_new = np.clip((tvisc4 + tirr ** 4) ** 0.25, tirr, 1.5e5) + T_new = np.clip((tvisc4 + tirr**4) ** 0.25, tirr, 1.5e5) diff = T - T_new return diff + def solve_viscous_heating_globally( - sigma_g, - tirr, - T, - omega_k, - alpha, - Z, - mu, - r, - xtol=1e-1, - nitermax=100, - interpolation_idx=None + sigma_g, + tirr, + T, + omega_k, + alpha, + Z, + mu, + r, + xtol=1e-1, + nitermax=100, + interpolation_idx=None, ): """ Solve the midplane tempeature due to viscous heating (only) @@ -68,44 +73,65 @@ def solve_viscous_heating_globally( sample = np.linspace(0.1, 300.0, 10000) * AU - T = np.array([optimize.root_scalar( - tvisc_fct, - bracket=(tirr[i], 1.5e4), - x0=T[i], # estimate T from last iteration - args=(sigma_g[i], tirr[i], omega_k[i], alpha[i] if hasattr(alpha, "__iter__") else alpha, Z[i] if hasattr(Z, "__iter__") else Z, mu[i] if hasattr(mu, "__iter__") else mu), - options={"maxiter": nitermax, "xtol": xtol} - ).root for i in range(len(T)) if interpolation_idx[i]]) - T_sample = savgol_filter(interp1d(r[interpolation_idx], T, fill_value="extrapolate", assume_sorted=True)(sample), 5, 2) + T = np.array( + [ + optimize.root_scalar( + tvisc_fct, + bracket=(tirr[i], 1.5e4), + x0=T[i], # estimate T from last iteration + args=( + sigma_g[i], + tirr[i], + omega_k[i], + alpha[i] if hasattr(alpha, "__iter__") else alpha, + Z[i] if hasattr(Z, "__iter__") else Z, + mu[i] if hasattr(mu, "__iter__") else mu, + ), + options={"maxiter": nitermax, "xtol": xtol}, + ).root + for i in range(len(T)) + if interpolation_idx[i] + ] + ) + T_sample = savgol_filter( + interp1d(r[interpolation_idx], T, fill_value="extrapolate", assume_sorted=True)( + sample + ), + 5, + 2, + ) T = interp1d(sample, T_sample, fill_value="extrapolate", assume_sorted=True)(r) T = np.clip(T, tirr, 1.5e5) return T class Disk(object): - """ Baseclass of the disk. """ + """Baseclass of the disk.""" def __init__( - self, - defaults, - M_STAR, - ALPHA, - chemistry, - init_grid=True, - time=0 * u.yr, - **kwargs + self, + defaults, + M_STAR, + ALPHA, + chemistry, + init_grid=True, + time=0 * u.yr, + **kwargs, ): - self.sigmin = 1e-60 # sigma_floor value - self.defaults = defaults # dictionary containing the defaults section of the config + self.sigmin = 1e-60 # sigma_floor value + self.defaults = ( + defaults # dictionary containing the defaults section of the config + ) self.M_s = M_STAR.cgs.value # stellar mass self.alpha = float(ALPHA) # alpha viscosity, fix yaml issue - self._chem = chemistry # reference to chemistry object + self._chem = chemistry # reference to chemistry object self.t_0 = 1 * time.cgs.value # starting time of the disk self.t = 1 * time.cgs.value # time in the disk (initialised with starting time) # init timestep self._timestep_input = defaults.get("DEF_TIMESTEP", None) if self._timestep_input is not None: - print(f'WARNING: We are using a custom timestep of {self._timestep_input}') + print(f"WARNING: We are using a custom timestep of {self._timestep_input}") self._timestep_input = self._timestep_input.cgs.value # init t_end @@ -116,20 +142,27 @@ def __init__( r_in = defaults.get("DEF_R_IN", 0.2 * u.au) # inner grid edge r_out = defaults.get("DEF_R_OUT", 100.0 * u.au) # outer grid edge - self.gridsize = defaults.get("DEF_GRIDSIZE", 1000) # number of grid cells + self.gridsize = defaults.get("DEF_GRIDSIZE", 1000) # number of grid cells if defaults.get("DEF_LIN_SPACING", True): - self.r = np.linspace(r_in, r_out, self.gridsize).to(u.au).value # linear grid + self.r = ( + np.linspace(r_in, r_out, self.gridsize).to(u.au).value + ) # linear grid else: - self.r = ( # log grid - np.logspace(np.log10(r_in.to(u.au).value), np.log10(r_out.to(u.au).value), self.gridsize) + self.r = np.logspace( # log grid + np.log10(r_in.to(u.au).value), + np.log10(r_out.to(u.au).value), + self.gridsize, ) - self.r_i = np.array([self.r[0], *((self.r[1:] + self.r[:-1]) / 2.0), self.r[-1]]) * u.au + self.r_i = ( + np.array([self.r[0], *((self.r[1:] + self.r[:-1]) / 2.0), self.r[-1]]) + * u.au + ) self.r = self.r * u.au self.r_i = self.r_i.cgs.value # interfacepositions of the radial grid - self.r = self.r.cgs.value # cellcenter of the radial grid + self.r = self.r.cgs.value # cellcenter of the radial grid self._init_params(**kwargs) @@ -149,45 +182,79 @@ def _init_params(self, **kwargs): # array holding the dust to gas ratio self.DTG = {"total": eval_kwargs(kwargs.pop("DTG_total", 0.015))} - self.DTG_small_grains = self.DTG["total"] * (1.0-0.75) # DTG small grains is only used for the temperature - self.conf_static = eval_kwargs(kwargs.get("static", False)) # config for static disk - self.conf_static_stokes = eval_kwargs(kwargs.get("static_stokes", False)) # config for static stokes number - self.conf_evaporation = eval_kwargs(kwargs.get("evaporation", True)) # config for static stokes number - self.conf_temp_evol = eval_kwargs(kwargs.get("temp_evol", False)) # use the temperature evolution. Advise: dont do! - - self._evap_width = eval_kwargs(kwargs.get("evap_width", 0.1*u.au)).cgs.value # config for static stokes number - - self.alpha_height = eval_kwargs(kwargs.get("ALPHAHEIGHT", self.alpha)) # vertical mixing alpha - self.alpha_frag = eval_kwargs(kwargs.get("ALPHAFRAG", self.alpha)) # fragmentation alpha, not used in Schneider & Bitsch 2021 - self.alpha_mig = eval_kwargs(kwargs.get("ALPHAMIG", self.alpha)) # migration alpha, not used in Schneider & Bitsch 2021 + self.DTG_small_grains = self.DTG["total"] * ( + 1.0 - 0.75 + ) # DTG small grains is only used for the temperature + self.conf_static = eval_kwargs( + kwargs.get("static", False) + ) # config for static disk + self.conf_static_stokes = eval_kwargs( + kwargs.get("static_stokes", False) + ) # config for static stokes number + self.conf_evaporation = eval_kwargs( + kwargs.get("evaporation", True) + ) # config for static stokes number + self.conf_temp_evol = eval_kwargs( + kwargs.get("temp_evol", False) + ) # use the temperature evolution. Advise: dont do! + + self._evap_width = eval_kwargs( + kwargs.get("evap_width", 0.1 * u.au) + ).cgs.value # config for static stokes number + + self.alpha_height = eval_kwargs( + kwargs.get("ALPHAHEIGHT", self.alpha) + ) # vertical mixing alpha + self.alpha_frag = eval_kwargs( + kwargs.get("ALPHAFRAG", self.alpha) + ) # fragmentation alpha, not used in Schneider & Bitsch 2021 + self.alpha_mig = eval_kwargs( + kwargs.get("ALPHAMIG", self.alpha) + ) # migration alpha, not used in Schneider & Bitsch 2021 self._chem.set_z(self.DTG.get("total", 0.015)) # set the solid to gas ratio if self.conf_static_stokes: self.f_m = np.ones_like(self.r) # Massdistribution factor for twopop model - self.stokes_number_small = np.zeros_like(self.r) # stokes number of the small population for twopop model - self.stokes_number_df = np.zeros_like(self.r) # stokes number of the drift limited fragmentation for twopop model - self.stokes_number_drift = np.zeros_like(self.r) # stokes number of the drift limited case for twopop model - self.stokes_number_frag = np.zeros_like(self.r) # stokes number of the fragmentation limited case for twopop model + self.stokes_number_small = np.zeros_like( + self.r + ) # stokes number of the small population for twopop model + self.stokes_number_df = np.zeros_like( + self.r + ) # stokes number of the drift limited fragmentation for twopop model + self.stokes_number_drift = np.zeros_like( + self.r + ) # stokes number of the drift limited case for twopop model + self.stokes_number_frag = np.zeros_like( + self.r + ) # stokes number of the fragmentation limited case for twopop model self.const_f_f = 0.37 # fudge fitting factors from Birnstiel2012 self.const_f_d = 0.55 # fudge fitting factors from Birnstiel2012 - self.const_N = 0.5 # N of Birnstiel2012 - self.a_0 = eval_kwargs(kwargs.get("a_0", 1e-4 * u.cm)) # grainsize of the small population - self.a_0 = self.a_0.cgs.value # grainsize of the small population - self.a_1 = eval_kwargs(kwargs.get("a_0", 1e-4 * u.cm)) * np.ones_like(self.r) # grainsize of the large population - self.a_1 = self.a_1.cgs.value # grainsize of the large population - - self.omega_k = np.sqrt(G * self.M_s / self.r ** 3) # kepler orbitalperiod - self.v_k = self.omega_k * self.r # kepler velocity - self.alpha_factor = np.ones_like(self.r) # gap factor that applies the gap to the alphaviscosity + self.const_N = 0.5 # N of Birnstiel2012 + self.a_0 = eval_kwargs( + kwargs.get("a_0", 1e-4 * u.cm) + ) # grainsize of the small population + self.a_0 = self.a_0.cgs.value # grainsize of the small population + self.a_1 = eval_kwargs(kwargs.get("a_0", 1e-4 * u.cm)) * np.ones_like( + self.r + ) # grainsize of the large population + self.a_1 = self.a_1.cgs.value # grainsize of the large population + + self.omega_k = np.sqrt(G * self.M_s / self.r**3) # kepler orbitalperiod + self.v_k = self.omega_k * self.r # kepler velocity + self.alpha_factor = np.ones_like( + self.r + ) # gap factor that applies the gap to the alphaviscosity # Init for output self.vr_dust = np.zeros_like(self.r[1:]) # dust velocity - self.vr_gas = np.zeros_like(self.r[1:]) # gas velocity + self.vr_gas = np.zeros_like(self.r[1:]) # gas velocity - tau_disk = eval_kwargs(kwargs.get("tau_disk", None)) # disk dispersal time - self.begin_photevap = eval_kwargs(kwargs.get("begin_photoevap", 0 * u.Myr).cgs.value) # time until disk disappears + tau_disk = eval_kwargs(kwargs.get("tau_disk", None)) # disk dispersal time + self.begin_photevap = eval_kwargs( + kwargs.get("begin_photoevap", 0 * u.Myr).cgs.value + ) # time until disk disappears if tau_disk is not None: self.conf_photoevaporation = True @@ -199,45 +266,59 @@ def _init_params(self, **kwargs): # init physical quantities # calculate the chemical composition of sigma: - self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(self.T) # chemistry vectors + self.chemistry_solid, self.chemistry_gas = self._chem.get_composition( + self.T + ) # chemistry vectors dtg = self._chem.get_solid_heavies(self.T) - self.sigma_g_components = self.chemistry_gas * self.sigma_g[:,np.newaxis,np.newaxis] # gas surface density as a compositional vector - self.sigma_g = np.sum(self.sigma_g_components[:, 0, :], axis=1) # gas surface density - - self.rho_g = self.sigma_g / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r) # gas volume density - self.grad_sig = np.gradient( # gradient of the surfacedensity + self.sigma_g_components = ( + self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] + ) # gas surface density as a compositional vector + self.sigma_g = np.sum( + self.sigma_g_components[:, 0, :], axis=1 + ) # gas surface density + + self.rho_g = self.sigma_g / ( + np.sqrt(2 * np.pi) * self.aspect_ratio * self.r + ) # gas volume density + self.grad_sig = np.gradient( # gradient of the surfacedensity np.log10(self.sigma_g), np.log10(self.r) ) # -1.5 for MMSN - self.c_s = self.aspect_ratio * self.r * self.omega_k # soundspeed - self.P = ( # pressure - self.c_s ** 2 - * self.sigma_g - / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r) + self.c_s = self.aspect_ratio * self.r * self.omega_k # soundspeed + self.P = ( # pressure + self.c_s**2 + * self.sigma_g + / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r) ) - self.grad_p = np.gradient( # pressuregradient + self.grad_p = np.gradient( # pressuregradient np.log10(self.P), np.log10(self.r) ) # -2.75 for MMSN, checked - self.T = 2.35 * self.c_s ** 2.0 * m_p / k_B # Temperature + self.T = 2.35 * self.c_s**2.0 * m_p / k_B # Temperature if not hasattr(self, "T_irr"): - self.T_irr = 1.0 * self.T # effective Temperature from irradiation + self.T_irr = 1.0 * self.T # effective Temperature from irradiation if not hasattr(self, "T_visc"): - self.T_visc = np.zeros_like(self.T) # viscous Temperature from viscous heating + self.T_visc = np.zeros_like( + self.T + ) # viscous Temperature from viscous heating self.compute_cs_and_hp() self.compute_nu() self.calc_opacity() - self.grad_T = np.gradient( # Temperature gradient + self.grad_T = np.gradient( # Temperature gradient np.log10(self.T), np.log10(self.r) ) # -0.5 for MMSN - self.eta = -0.5 * self.aspect_ratio ** 2 * self.grad_p + self.eta = -0.5 * self.aspect_ratio**2 * self.grad_p enclosed_mask = self.r < 200.0 * AU - self.disk_mass = np.trapz(2 * np.pi * self.sigma_g[enclosed_mask] * self.r[enclosed_mask], - self.r[enclosed_mask]) + self.disk_mass = np.trapz( + 2 * np.pi * self.sigma_g[enclosed_mask] * self.r[enclosed_mask], + self.r[enclosed_mask], + ) self.disk_mass_dust = 0 - self.m_dot_components = self.compute_mdot_at_interfaces(self.sigma_g_components) + self.m_dot_components = self.compute_mdot_at_interfaces( + self.sigma_g_components + ) self.m_dot = np.sum(self.m_dot_components[:, 0, :], axis=1) self.pebble_flux = np.zeros_like(self.r) @@ -248,8 +329,10 @@ def _init_params(self, **kwargs): print("caution: Unstable disk: Q_min = {}".format(np.min(self.Q))) def calc_pebble_iso(self, diffusion=False): - """Calculate the pebble isolation mass. """ - f_fit = (self.aspect_ratio / 0.05) ** 3.0 * (0.34 * (np.log10(1e-3) / np.log10(self.alpha)) ** 4.0 + 0.66) + """Calculate the pebble isolation mass.""" + f_fit = (self.aspect_ratio / 0.05) ** 3.0 * ( + 0.34 * (np.log10(1e-3) / np.log10(self.alpha)) ** 4.0 + 0.66 + ) if diffusion: self.peb_iso_pi_crit = self.alpha / (2.0 * self.stokes_number_pebbles) @@ -257,7 +340,9 @@ def calc_pebble_iso(self, diffusion=False): else: self.peb_iso = 25.0 * ME * f_fit - def init_pebble_parameters(self, epsilon_p, rho_solid, t_f, u_frag, calculate_rho_solid, calculate_u_frag): + def init_pebble_parameters( + self, epsilon_p, rho_solid, t_f, u_frag, calculate_rho_solid, calculate_u_frag + ): """ Sets the pebble parameters needed for pebble evaporation model all in cgs @@ -275,7 +360,7 @@ def init_pebble_parameters(self, epsilon_p, rho_solid, t_f, u_frag, calculate_rh self.epsilon_p = epsilon_p self.rho_solid = rho_solid if self.rho_solid is not None: - self.rho_solid = np.ones_like(self.sigma_dust)*self.rho_solid + self.rho_solid = np.ones_like(self.sigma_dust) * self.rho_solid self._calculate_rho_solid = calculate_rho_solid self._calculate_u_frag = calculate_u_frag @@ -313,23 +398,62 @@ def calc_sigma_dust(self): self.solid_fraction = self._chem.get_solid_heavies(self.T) self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(self.T) - sigma_g0 = self.sigma_g / (1.0 + (self._chem.dtg0-self.solid_fraction)) # get the fraction of the background gas - - self.sigma_dust_components = self.chemistry_solid * self.solid_fraction[:, np.newaxis, np.newaxis] * sigma_g0[:, np.newaxis, - np.newaxis] - self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] + sigma_g0 = self.sigma_g / ( + 1.0 + (self._chem.dtg0 - self.solid_fraction) + ) # get the fraction of the background gas + + self.sigma_dust_components = ( + self.chemistry_solid + * self.solid_fraction[:, np.newaxis, np.newaxis] + * sigma_g0[:, np.newaxis, np.newaxis] + ) + self.sigma_g_components = ( + self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] + ) self.sigma_dust = np.sum(self.sigma_dust_components[:, 1, :], axis=1) self.sigma_dust_components[:, 1, 0] = np.zeros_like(self.sigma_g) - assert np.all(np.isclose((self.sigma_g + self.sigma_dust), sigma_g0 * (1.0 + self._chem.dtg0))), "error in chemistry: The total masses do not sum up correctly" - assert np.all(np.isclose(((self.sigma_dust_components[:,1,:]+self.sigma_g_components[:,1,:])/(self.sigma_g+self.sigma_dust)[:,np.newaxis]).std(axis=0),0)), "error in chemistry: The individual masses do not sum up to a constant value" - assert np.all(np.isclose(np.sum(self.sigma_dust_components[:, 1, :] + self.sigma_g_components[:, 1, :], axis=1) / sigma_g0, - (1+self._chem.dtg0))), "error in chemistry: the dust to gas ratio is not correct" - assert np.all(np.isclose( - np.sum(self.sigma_dust_components[:, 1, 1:] + self.sigma_g_components[:, 1, 1:], axis=1) / sigma_g0, - self._chem.dtg0)), "error in chemistry: the dust to gas ratio is not correct" + assert np.all( + np.isclose( + (self.sigma_g + self.sigma_dust), sigma_g0 * (1.0 + self._chem.dtg0) + ) + ), "error in chemistry: The total masses do not sum up correctly" + assert np.all( + np.isclose( + ( + ( + self.sigma_dust_components[:, 1, :] + + self.sigma_g_components[:, 1, :] + ) + / (self.sigma_g + self.sigma_dust)[:, np.newaxis] + ).std(axis=0), + 0, + ) + ), "error in chemistry: The individual masses do not sum up to a constant value" + assert np.all( + np.isclose( + np.sum( + self.sigma_dust_components[:, 1, :] + + self.sigma_g_components[:, 1, :], + axis=1, + ) + / sigma_g0, + (1 + self._chem.dtg0), + ) + ), "error in chemistry: the dust to gas ratio is not correct" + assert np.all( + np.isclose( + np.sum( + self.sigma_dust_components[:, 1, 1:] + + self.sigma_g_components[:, 1, 1:], + axis=1, + ) + / sigma_g0, + self._chem.dtg0, + ) + ), "error in chemistry: the dust to gas ratio is not correct" self.sigmin_peb = 1e-60 @@ -349,18 +473,25 @@ def remove_from_gas(self, mdot, f_red, idx): """ try: sum_f_red = np.sum(f_red) - delta_sigma = f_red / ( - (self.r_i[1:][idx] - self.r_i[:-1][idx]) * sum_f_red) + delta_sigma = f_red / ((self.r_i[1:][idx] - self.r_i[:-1][idx]) * sum_f_red) if hasattr(self, "sigdot"): # remove using visc disk evolution self.sigdot_gas_acc = np.zeros_like(self.sigdot) self.sigdot_gas_acc[idx] = delta_sigma[:, np.newaxis] * mdot[1] else: # there is no visc disk evolution, remove directly instead sum_mdot_el = np.sum(mdot, axis=1)[0] - self.sigma_g[idx] -= sum_mdot_el * delta_sigma / (2 * np.pi * self.r[idx]) * self.dt - self.sigma_g_components[idx] -= delta_sigma[:, np.newaxis, np.newaxis] / ( - 2.0 * np.pi * self.r[idx, np.newaxis, np.newaxis]) * mdot * self.dt + self.sigma_g[idx] -= ( + sum_mdot_el * delta_sigma / (2 * np.pi * self.r[idx]) * self.dt + ) + self.sigma_g_components[idx] -= ( + delta_sigma[:, np.newaxis, np.newaxis] + / (2.0 * np.pi * self.r[idx, np.newaxis, np.newaxis]) + * mdot + * self.dt + ) self.sigma_g[self.sigma_g < self.sigmin] = self.sigmin - self.sigma_g_components[self.sigma_g_components < self.sigmin] = self.sigmin + self.sigma_g_components[ + self.sigma_g_components < self.sigmin + ] = self.sigmin return except IndexError: return @@ -382,23 +513,30 @@ def remove_from_peb(self, mdot, idx): self.sigdot_peb_acc[idx] = delta_sigma * mdot[1] # fill sigdot else: mdot_sum = np.sum(mdot, axis=1)[0] - self.sigma_dust[idx] -= delta_sigma / (2.0 * np.pi * self.r[idx]) * mdot_sum * self.dt - self.sigma_dust_components[idx] -= delta_sigma / ( - 2.0 * np.pi * self.r[idx, np.newaxis, np.newaxis]) * mdot * self.dt + self.sigma_dust[idx] -= ( + delta_sigma / (2.0 * np.pi * self.r[idx]) * mdot_sum * self.dt + ) + self.sigma_dust_components[idx] -= ( + delta_sigma + / (2.0 * np.pi * self.r[idx, np.newaxis, np.newaxis]) + * mdot + * self.dt + ) self.sigma_dust[self.sigma_dust < self.sigmin_peb] = self.sigmin_peb self.sigma_dust_components[ - self.sigma_dust_components < self.sigmin_peb] = self.sigmin_peb + self.sigma_dust_components < self.sigmin_peb + ] = self.sigmin_peb return except IndexError: return def apply_gap_profile(self, a_p, fsigma, FWHM): - #fsigma = np.maximum(fsigma, 1e-7) + # fsigma = np.maximum(fsigma, 1e-7) self.alpha_factor = 1 - gaussian(self.r, 1.0 - fsigma, a_p, FWHM) - self.interpolation_idx = np.abs(self.alpha_factor - 1.0 ) < 0.01 + self.interpolation_idx = np.abs(self.alpha_factor - 1.0) < 0.01 self.interpolation_idx[0] = True - alpha_factor_interface = 0.5*(self.alpha_factor[1:]+self.alpha_factor[:-1]) + alpha_factor_interface = 0.5 * (self.alpha_factor[1:] + self.alpha_factor[:-1]) midinterface = np.abs(alpha_factor_interface - 1) < 0.01 self.interpolation_idx_interface = np.array([True, *midinterface, True]) @@ -428,23 +566,38 @@ def compute_sizes(self): msil = np.sum(self.sigma_dust_components[:, 1, 8:], axis=-1) mice = np.sum(self.sigma_dust_components[:, 1, 1:8], axis=-1) # assert np.all(np.isclose(msil + mice, self.sigma_dust)), 'The internal ice densities do not add up correctly' - self.rho_solid = (msil+mice) / (msil/DENSITY_SILICATES + mice/DENSITY_ICE + 1e-90) + self.rho_solid = (msil + mice) / ( + msil / DENSITY_SILICATES + mice / DENSITY_ICE + 1e-90 + ) self.rho_solid = np.maximum(self.rho_solid, DENSITY_ICE) if self._calculate_u_frag: # Calculate the fragmentation velocity according to the Drazkowska 2017 approach: 1m/s if silicate, 10m/s if icy mice = np.sum(self.sigma_dust_components[:, 1, 1:8], axis=-1) - self.u_frag = np.where(mice > 0.01*self.sigma_dust, 1000, 100) - - a_frag = self.const_f_f * 2.0 / (3.0 * np.pi) * self.sigma_g / ( - self.alpha_frag * self.rho_solid) * self.u_frag ** 2.0 / ( - self.c_s ** 2.0) + self.u_frag = np.where(mice > 0.01 * self.sigma_dust, 1000, 100) + + a_frag = ( + self.const_f_f + * 2.0 + / (3.0 * np.pi) + * self.sigma_g + / (self.alpha_frag * self.rho_solid) + * self.u_frag**2.0 + / (self.c_s**2.0) + ) tau_grow = self.sigma_g / ((self.omega_k * self.sigma_dust) + 1e-300) dlnpdlnr = np.abs(self.grad_p) - a_drift = self.const_f_d * 2 * self.sigma_dust / ( - np.pi * self.rho_solid) * self.v_k ** 2 / self.c_s ** 2 / dlnpdlnr - St_df = self.u_frag * self.v_k / (dlnpdlnr * self.c_s ** 2 * (1 - self.const_N)) + a_drift = ( + self.const_f_d + * 2 + * self.sigma_dust + / (np.pi * self.rho_solid) + * self.v_k**2 + / self.c_s**2 + / dlnpdlnr + ) + St_df = self.u_frag * self.v_k / (dlnpdlnr * self.c_s**2 * (1 - self.const_N)) a_df = St_df * 2 * self.sigma_g / (np.pi * self.rho_solid) a_1 = np.min([a_drift, a_frag, a_df], axis=0) @@ -461,9 +614,11 @@ def compute_sizes(self): self.stokes_number_pebbles = a_1 * stokes_factor self.stokes_number_small = self.a_0 * stokes_factor - self.f_m = np.where(np.argmin([a_drift, a_frag, a_df], axis=0) == 0, - 0.97 * np.ones_like(self.stokes_number_pebbles), - 0.75 * np.ones_like(self.stokes_number_pebbles)) + self.f_m = np.where( + np.argmin([a_drift, a_frag, a_df], axis=0) == 0, + 0.97 * np.ones_like(self.stokes_number_pebbles), + 0.75 * np.ones_like(self.stokes_number_pebbles), + ) self.stokes_number_df = 1.0 * St_df self.stokes_number_drift = stokes_factor * a_drift @@ -494,8 +649,11 @@ def compute_viscous_evolution(self): if self.conf_photoevaporation: self.calc_sigdot_photoevap() - sigma_g_components_mol = self.get_viscous_evolution_next_timestep_components(self.sigma_g_components, - self.dt) + sigma_g_components_mol = ( + self.get_viscous_evolution_next_timestep_components( + self.sigma_g_components, self.dt + ) + ) self.compute_dust_next_timestep(self.dt) self.compute_sigma_g_from_mol(sigma_g_components_mol) @@ -503,19 +661,28 @@ def compute_viscous_evolution(self): else: # init sigma_g_components self.compute_disktmid() dtg = self._chem.get_solid_heavies(self.T) - self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(self.T) + self.chemistry_solid, self.chemistry_gas = self._chem.get_composition( + self.T + ) self.mu = self._chem.mu - self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] - + self.sigma_g_components = ( + self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] + ) self.compute_disktmid() if not self.conf_evaporation: - self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(self.T) + self.chemistry_solid, self.chemistry_gas = self._chem.get_composition( + self.T + ) self.mu = self._chem.mu - self.sigma_dust_components = self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis] - self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] + self.sigma_dust_components = ( + self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis] + ) + self.sigma_g_components = ( + self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] + ) if hasattr(self, "vr_peb"): self.calc_pebble_flux(noupdate=False) @@ -525,32 +692,32 @@ def compute_viscous_evolution(self): self._old_disk_mass = 1.0 * self.disk_mass self._old_disk_mass_dust = 1.0 * self.disk_mass_dust - self.disk_mass = np.trapz(2.0 * np.pi * self.sigma_g[enclosed_mask] * self.r[enclosed_mask], - self.r[enclosed_mask]) - self.disk_mass_dust = np.trapz(2 * np.pi * self.sigma_dust[enclosed_mask] * self.r[enclosed_mask], - self.r[enclosed_mask]) + self.disk_mass = np.trapz( + 2.0 * np.pi * self.sigma_g[enclosed_mask] * self.r[enclosed_mask], + self.r[enclosed_mask], + ) + self.disk_mass_dust = np.trapz( + 2 * np.pi * self.sigma_dust[enclosed_mask] * self.r[enclosed_mask], + self.r[enclosed_mask], + ) - #assert (self.disk_mass + self.disk_mass_dust) <= ( + # assert (self.disk_mass + self.disk_mass_dust) <= ( # self._old_disk_mass + self._old_disk_mass_dust), "Disk mass is growing" self.m_dot_components = self.compute_mdot_at_interfaces(self.sigma_g_components) self.m_dot = np.sum(self.m_dot_components[:, 0, :], axis=1) self.P = ( - self.c_s ** 2.0 - * self.sigma_g - / (np.sqrt(2.0 * np.pi) * self.aspect_ratio * self.r) + self.c_s**2.0 + * self.sigma_g + / (np.sqrt(2.0 * np.pi) * self.aspect_ratio * self.r) ) - self.grad_sig = np.gradient( - np.log10(self.sigma_g), np.log10(self.r) - ) - self.grad_T = np.gradient( - np.log10(self.T), np.log10(self.r) - ) # -0.5 for MMSN + self.grad_sig = np.gradient(np.log10(self.sigma_g), np.log10(self.r)) + self.grad_T = np.gradient(np.log10(self.T), np.log10(self.r)) # -0.5 for MMSN self.grad_p = np.gradient(np.log10(self.P), np.log10(self.r)) - self.eta = -0.5 * self.aspect_ratio ** 2 * self.grad_p + self.eta = -0.5 * self.aspect_ratio**2 * self.grad_p self.solid_fraction = self.sigma_dust / self.sigma_g @@ -571,7 +738,9 @@ def get_viscous_evolution_next_timestep_components(self, sigma_g_components, dt) d = 3.0 * nu[:, np.newaxis] v = np.zeros(y.shape) - s = 1 * self.sigdot # if self.t < self.begin_photevap else np.zeros_like(self.sigdot) + s = ( + 1 * self.sigdot + ) # if self.t < self.begin_photevap else np.zeros_like(self.sigdot) if hasattr(self, "sigdot_gas_acc"): s -= self.sigdot_gas_acc @@ -581,7 +750,15 @@ def get_viscous_evolution_next_timestep_components(self, sigma_g_components, dt) s -= self.sigdot_photoevap # checke das vorzeichen! - clip_val = (y - 2.0 * np.pi * 0.9 * x[:, np.newaxis] * self.sigmin * self.chemistry_gas[:, 1, :]) / dt + clip_val = ( + y + - 2.0 + * np.pi + * 0.9 + * x[:, np.newaxis] + * self.sigmin + * self.chemistry_gas[:, 1, :] + ) / dt s = np.clip(s, -clip_val, clip_val) # # Set boundary conditions @@ -612,8 +789,12 @@ def calc_pebble_flux(self, noupdate=False): ------- """ - vr_peb = interp1d(self.r_i[1:-1], self.vr_peb, assume_sorted=True, fill_value="extrapolate")(self.r) - self.pebble_flux = np.abs(2.0 * np.pi * self.r * vr_peb * self.sigma_dust * self.f_m) + vr_peb = interp1d( + self.r_i[1:-1], self.vr_peb, assume_sorted=True, fill_value="extrapolate" + )(self.r) + self.pebble_flux = np.abs( + 2.0 * np.pi * self.r * vr_peb * self.sigma_dust * self.f_m + ) if not noupdate: self.cum_pebble_flux += self.pebble_flux * self.dt @@ -637,32 +818,56 @@ def calc_sigdot(self, icelines_mol): y = 2.0 * np.pi * x[:, np.newaxis] * self.sigma_g_components[:, 1, :] self.sigdot = np.zeros(y.shape) # interpolate back to cell centers - vr_dust = interp1d(self.r_i[1:-1], self.vr_dust, assume_sorted=True, fill_value="extrapolate")(self.r) + vr_dust = interp1d( + self.r_i[1:-1], self.vr_dust, assume_sorted=True, fill_value="extrapolate" + )(self.r) if not self.conf_evaporation: included_icelines = [] else: # curate icelines: (only use icelines that are well enough inside of the grid) - #included_icelines = np.arange(icelines_mol.shape[0])[ + # included_icelines = np.arange(icelines_mol.shape[0])[ # np.logical_and(self.r[icelines_mol] > self.r[5], self.r[icelines_mol] < self.r[-5])] included_icelines = np.arange(1, icelines_mol.shape[0]) # exclude rest_mol - sigdot_cond_const = -3.0 * self.epsilon_p / ( - 2.0 * np.pi * self.rho_solid[:, np.newaxis]) * y * self.sigma_dust[:, np.newaxis] * self.omega_k[:, - np.newaxis] * np.sqrt( - self._chem.mu[:,np.newaxis]) + sigdot_cond_const = ( + -3.0 + * self.epsilon_p + / (2.0 * np.pi * self.rho_solid[:, np.newaxis]) + * y + * self.sigma_dust[:, np.newaxis] + * self.omega_k[:, np.newaxis] + * np.sqrt(self._chem.mu[:, np.newaxis]) + ) sigdot_cond_const *= ( - self.f_m[:, np.newaxis] / self.a_1[:, np.newaxis] + (1.0 - self.f_m[:, np.newaxis]) / self.a_0) + self.f_m[:, np.newaxis] / self.a_1[:, np.newaxis] + + (1.0 - self.f_m[:, np.newaxis]) / self.a_0 + ) for i in included_icelines: # exclude rest_mol - self.sigdot[icelines_mol[i]:-1, i] = sigdot_cond_const[icelines_mol[i]:-1, i] / np.sqrt(np.squeeze(self._chem.M_array)[i]) - self.sigdot[:icelines_mol[i], i] = 2 * np.pi * x[:icelines_mol[i]] * self.sigma_dust_components[ - :icelines_mol[i], 1, i] * np.abs( - vr_dust[:icelines_mol[i]]) / self._evap_width + self.sigdot[icelines_mol[i] : -1, i] = sigdot_cond_const[ + icelines_mol[i] : -1, i + ] / np.sqrt(np.squeeze(self._chem.M_array)[i]) + self.sigdot[: icelines_mol[i], i] = ( + 2 + * np.pi + * x[: icelines_mol[i]] + * self.sigma_dust_components[: icelines_mol[i], 1, i] + * np.abs(vr_dust[: icelines_mol[i]]) + / self._evap_width + ) - min_density = 2.0 * np.pi * x[:, np.newaxis] * 0.9 * np.minimum(self.sigma_dust_components[:, 1, :], - self.sigma_g_components[:, 1, :]) + min_density = ( + 2.0 + * np.pi + * x[:, np.newaxis] + * 0.9 + * np.minimum( + self.sigma_dust_components[:, 1, :], self.sigma_g_components[:, 1, :] + ) + ) - self.sigdot = np.clip(self.sigdot, -min_density / self.dt, - min_density / self.dt) # do not transfer more then available in one timestep + self.sigdot = np.clip( + self.sigdot, -min_density / self.dt, min_density / self.dt + ) # do not transfer more then available in one timestep return def calc_sigdot_photoevap(self): @@ -675,10 +880,18 @@ def calc_sigdot_photoevap(self): """ if self.t > self.begin_photevap: if hasattr(self, "tau_disk"): - self.sigdot_photoevap = 1.0 / self.tau_disk * self.sigma_g_components[:, 1, :] * 2.0 * np.pi * self.r[:, - np.newaxis] + self.sigdot_photoevap = ( + 1.0 + / self.tau_disk + * self.sigma_g_components[:, 1, :] + * 2.0 + * np.pi + * self.r[:, np.newaxis] + ) else: - self.sigdot_photoevap = np.zeros_like(self.r) * self.chemistry_gas[:, 1, :] + self.sigdot_photoevap = ( + np.zeros_like(self.r) * self.chemistry_gas[:, 1, :] + ) def get_dust_radial_drift_next_timestep(self, dt): """ @@ -706,20 +919,31 @@ def get_dust_radial_drift_next_timestep(self, dt): # lmax = np.ones_like(self.r, dtype="bool") lmax = self.T < 2000.0 x = self.r[lmax] - y = 2.0 * np.pi * x[:, np.newaxis] * self.sigma_dust_components[lmax, 1, 1:] # Dust + y = ( + 2.0 * np.pi * x[:, np.newaxis] * self.sigma_dust_components[lmax, 1, 1:] + ) # Dust g = 2.0 * np.pi * x * self.sigma_g[lmax] # Gas g = g[:, np.newaxis] d = self.nu[lmax] di = 0.5 * (d[1:] + d[:-1])[:, np.newaxis] vi = self.vr_dust[lmax[:-1], np.newaxis] - s = - self.sigdot[lmax, 1:] # if self.t < self.begin_photevap else np.zeros_like(self.sigdot[:, 1:]) + s = -self.sigdot[ + lmax, 1: + ] # if self.t < self.begin_photevap else np.zeros_like(self.sigdot[:, 1:]) if hasattr(self, "sigdot_peb_acc"): s -= self.sigdot_peb_acc[lmax, 1:] - with np.errstate(under='ignore'): - clip_val = (y - 2 * np.pi * x[:, np.newaxis] * self.sigmin_peb * self.chemistry_solid[lmax, 1, 1:]) / dt + with np.errstate(under="ignore"): + clip_val = ( + y + - 2 + * np.pi + * x[:, np.newaxis] + * self.sigmin_peb + * self.chemistry_solid[lmax, 1, 1:] + ) / dt s = np.clip(s, -clip_val, clip_val) # @@ -728,21 +952,39 @@ def get_dust_radial_drift_next_timestep(self, dt): bcl = (1, 0, 0, 1) # bcr = (1, 0, 0, 1) # bcl = (0, 1, 2 * np.pi * x[0] * self.chemistry_solid[0, 1, 1:] * self.sigmin_peb, 0) - bcr = (0, 1, 2 * np.pi * x[-1] * self.chemistry_solid[-1, 1, 1:] * self.sigmin_peb, 0) + bcr = ( + 0, + 1, + 2 * np.pi * x[-1] * self.chemistry_solid[-1, 1, 1:] * self.sigmin_peb, + 0, + ) # # Get the new value of y after one time step dt # - y = solvediffonedee_components(x, y, vi, di, g, s, bcl, bcr, dt=dt, int=True, upwind=True) - y = np.maximum(y, 2 * np.pi * x[:, np.newaxis] * self.chemistry_solid[lmax, 1, 1:] * self.sigmin_peb) + y = solvediffonedee_components( + x, y, vi, di, g, s, bcl, bcr, dt=dt, int=True, upwind=True + ) + y = np.maximum( + y, + 2 + * np.pi + * x[:, np.newaxis] + * self.chemistry_solid[lmax, 1, 1:] + * self.sigmin_peb, + ) # y = np.maximum(y, 2 * np.pi * x[:, np.newaxis] * 1e-60) # # Obtain new sigdust # - new_sigma_dust = y / (2 * np.pi * x[:, np.newaxis]) # add offset for numerical reasons + new_sigma_dust = y / ( + 2 * np.pi * x[:, np.newaxis] + ) # add offset for numerical reasons # # Return # - sigma_dust = np.ones_like(self.sigma_dust_components[:, 1, 1:]) * self.sigmin_peb + sigma_dust = ( + np.ones_like(self.sigma_dust_components[:, 1, 1:]) * self.sigmin_peb + ) sigma_dust[lmax] = new_sigma_dust return sigma_dust @@ -770,10 +1012,15 @@ def compute_sigma_g_from_mol(self, sigma_g_components_mol): def compute_sigma_dust_from_mol(self, sigma_dust_components_mol): if self.conf_evaporation: - self.chemistry_solid = self._chem.get_solid_composition(sigma_dust_components_mol, self.T) + self.chemistry_solid = self._chem.get_solid_composition( + sigma_dust_components_mol, self.T + ) - self.sigma_dust = np.where(np.sum(self.chemistry_solid[:, 0], axis=1) > 0.0, - np.sum(sigma_dust_components_mol, axis=1), np.ones_like(self.r) * self.sigmin_peb) + self.sigma_dust = np.where( + np.sum(self.chemistry_solid[:, 0], axis=1) > 0.0, + np.sum(sigma_dust_components_mol, axis=1), + np.ones_like(self.r) * self.sigmin_peb, + ) assert np.all(self.sigma_dust > 0), "negative surface densities!" def compute_mdot_at_interfaces(self, sigma_g, oned=False): @@ -788,14 +1035,16 @@ def compute_mdot_at_interfaces(self, sigma_g, oned=False): y = 2.0 * np.pi * self.r[:, np.newaxis, np.newaxis] * sigma_g else: y = 2.0 * np.pi * self.r * sigma_g - g = self.r ** 0.5 / (self.nu / self.alpha_factor) + g = self.r**0.5 / (self.nu / self.alpha_factor) d = 3.0 * self.nu / self.alpha_factor v = np.zeros(len(x)) # # Compute the mdot # - flux = -getfluxonedee(x, y, v, d, g, int=False, oned=oned) # returns flux at interfaces + flux = -getfluxonedee( + x, y, v, d, g, int=False, oned=oned + ) # returns flux at interfaces return np.array([flux[0], *flux, flux[-1]]) def compute_vr_at_interfaces(self): @@ -816,7 +1065,7 @@ def compute_vr_at_interfaces(self): # # Now compute the radial velocity from the mdot # - self.vr_gas = - self.m_dot[1:-1] / (q + 1e-90) + self.vr_gas = -self.m_dot[1:-1] / (q + 1e-90) def calc_velocity_from_stokes(self, stokes, dlnpdlnr, csi, vki, vrgas): """ @@ -834,8 +1083,9 @@ def calc_velocity_from_stokes(self, stokes, dlnpdlnr, csi, vki, vrgas): """ stokes_i = (stokes[1:] + stokes[:-1]) / 2.0 + 1e-60 - return vrgas / (1.0 + stokes_i ** 2) + dlnpdlnr * csi ** 2 / vki / ( - stokes_i + 1.0 / stokes_i) + return vrgas / (1.0 + stokes_i**2) + dlnpdlnr * csi**2 / vki / ( + stokes_i + 1.0 / stokes_i + ) def compute_interface_velocities(self): """ @@ -856,9 +1106,9 @@ def compute_interface_velocities(self): # self.P = ( - self.c_s ** 2.0 - * self.sigma_g - / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r) + self.c_s**2.0 + * self.sigma_g + / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r) ) self.grad_p = np.gradient( np.log10(self.P), np.log10(self.r) @@ -871,22 +1121,28 @@ def compute_interface_velocities(self): # if not self.conf_static_stokes: f_m_i = (self.f_m[1:] + self.f_m[:-1]) / 2.0 - self.vr_peb = self.calc_velocity_from_stokes(self.stokes_number_pebbles, dlnpdlnr, csi, vki, vrgas) - self.vr_small = self.calc_velocity_from_stokes(self.stokes_number_small, dlnpdlnr, csi, vki, vrgas) + self.vr_peb = self.calc_velocity_from_stokes( + self.stokes_number_pebbles, dlnpdlnr, csi, vki, vrgas + ) + self.vr_small = self.calc_velocity_from_stokes( + self.stokes_number_small, dlnpdlnr, csi, vki, vrgas + ) self.vr_dust = f_m_i * self.vr_peb + (1.0 - f_m_i) * self.vr_small else: - self.vr_dust = self.calc_velocity_from_stokes(self.stokes_number_pebbles, dlnpdlnr, csi, vki, vrgas) + self.vr_dust = self.calc_velocity_from_stokes( + self.stokes_number_pebbles, dlnpdlnr, csi, vki, vrgas + ) self.vr_small = 0.0 self.vr_peb = 1.0 * self.vr_dust def init_T(self): - flux = self.lstar / (4.0 * np.pi * self.r ** 2.0) + flux = self.lstar / (4.0 * np.pi * self.r**2.0) self.qirr = ( - flux * self.flang + flux * self.flang ) # Factor 2 for two sides, factor 0.5 for half irradiated down self.T_irr = ( - 0.5 * self.qirr / sr - ) ** 0.25 # Factor 0.5 because cooling is two-sided + 0.5 * self.qirr / sr + ) ** 0.25 # Factor 0.5 because cooling is two-sided self.T_irr = np.maximum(self.T_irr, 10.0) self.T = 1.0 * self.T_irr @@ -922,13 +1178,13 @@ def compute_disktmid(self): self.qirr = Irradiative heating rate of the disk [erg/cm^2/s] self.qvisc = (if vischeat==True) Viscous heating rate [erg/cm^2/s] """ - if hasattr(self,"alpha_factor"): + if hasattr(self, "alpha_factor"): alpha = self.alpha / self.alpha_factor else: alpha = self.alpha if hasattr(self, "T_irr"): - if (self.conf_temp_evol or not self._T_init_finished): + if self.conf_temp_evol or not self._T_init_finished: self.T = solve_viscous_heating_globally( sigma_g=self.sigma_g, tirr=self.T_irr, @@ -938,9 +1194,12 @@ def compute_disktmid(self): Z=self.DTG_small_grains, r=self.r, mu=self._chem.mu, - interpolation_idx=None if not hasattr(self,"interpolation_idx") else self.interpolation_idx) + interpolation_idx=None + if not hasattr(self, "interpolation_idx") + else self.interpolation_idx, + ) - self.T_visc = (self.T ** 4.0 - self.T_irr ** 4.0) ** 0.25 + self.T_visc = (self.T**4.0 - self.T_irr**4.0) ** 0.25 else: self.init_T() @@ -957,14 +1216,14 @@ def calc_opacity(self): self.opacity = opacity_fct(self.rho_g, self.T, Z=self.DTG_small_grains) def compute_cs_and_hp(self): - """ update all quantities related to T """ + """update all quantities related to T""" self.c_s = (k_B * self.T / (self._chem.mu * m_p)) ** 0.5 - self.aspect_ratio = self.c_s / ( - self.omega_k * self.r - ) + self.aspect_ratio = self.c_s / (self.omega_k * self.r) if hasattr(self, "sigma_g"): - self.rho_g = self.sigma_g / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r) + self.rho_g = self.sigma_g / ( + np.sqrt(2 * np.pi) * self.aspect_ratio * self.r + ) def compute_nu(self): - """ compute the viscosity. """ + """compute the viscosity.""" self.nu = self.alpha * self.c_s * self.c_s / self.omega_k diff --git a/chemcomp/disks/kees.py b/chemcomp/disks/kees.py index 5d27ebc..04868b4 100755 --- a/chemcomp/disks/kees.py +++ b/chemcomp/disks/kees.py @@ -11,18 +11,18 @@ class DiskLabDisk(Disk): - """ LBP, viscous disk evolution and viscous heating. Disk of Schneider & Bitsch (2021)a,b.""" + """LBP, viscous disk evolution and viscous heating. Disk of Schneider & Bitsch (2021)a,b.""" def __init__( - self, - defaults, - chemistry, - M_STAR=1.0 * u.M_sun, - M0=1e-2 * u.M_sun, - R0=200 * u.au, - time=1.0 * u.yr, - ALPHA=1e-3, - **kwargs + self, + defaults, + chemistry, + M_STAR=1.0 * u.M_sun, + M0=1e-2 * u.M_sun, + R0=200 * u.au, + time=1.0 * u.yr, + ALPHA=1e-3, + **kwargs ): super().__init__(defaults, M_STAR, ALPHA, chemistry, time=time, **kwargs) @@ -35,12 +35,17 @@ def __init__( self.init_T() t_iter = 0 while True: - self.make_disk_from_lbp_alpha(M0=M0.cgs.value, R0=R0.cgs.value, alpha=float(ALPHA), time=time.cgs.value) + self.make_disk_from_lbp_alpha( + M0=M0.cgs.value, + R0=R0.cgs.value, + alpha=float(ALPHA), + time=time.cgs.value, + ) T0 = self.T.copy() self.compute_disktmid() self._chem.get_composition(self.T) # updates mu self._init_params(**kwargs) - if np.all(np.abs((T0-self.T)/T0) < 1e-2): + if np.all(np.abs((T0 - self.T) / T0) < 1e-2): break assert t_iter < 1000, "We have trouble to converge the initial T-profile" @@ -63,33 +68,37 @@ def make_disk_from_simplified_lbp(self, Sigc, Rc, gam): gam = The gamma coefficient of the LBP model """ r = self.r / AU - sigma = Sigc * (Rc / r) ** gam * np.exp(-(r / Rc) ** (2. - gam)) + sigma = Sigc * (Rc / r) ** gam * np.exp(-((r / Rc) ** (2.0 - gam))) sigma[sigma < self.sigmin] = self.sigmin self.sigma_g = sigma def make_disk_from_lyndenbellpringle(self, M0, R0, nu0, gam, time): """ - Make a model of the gas surface density of the disk. - Here the model is set by the Lynden-Bell & Pringle solution. - I follow here Lodato, Scardoni, Manara & Testi (2017) - MNRAS 472, 4700, their Eqs. 2, 3 and 4. - - ARGUMENTS: - M0 = Initial disk mass [g] - R0 = Initial disk radius [cm] - nu0 = Viscosity nu at R0 [cm^2/s] - gam = The gamma coefficient of the LBP model - time = Time after start [s] - """ + Make a model of the gas surface density of the disk. + Here the model is set by the Lynden-Bell & Pringle solution. + I follow here Lodato, Scardoni, Manara & Testi (2017) + MNRAS 472, 4700, their Eqs. 2, 3 and 4. + + ARGUMENTS: + M0 = Initial disk mass [g] + R0 = Initial disk radius [cm] + nu0 = Viscosity nu at R0 [cm^2/s] + gam = The gamma coefficient of the LBP model + time = Time after start [s] + """ r = self.r nu = nu0 * (r / R0) ** gam # Eq. 2 - tnu = R0 ** 2 / (3.0 * (2. - gam) ** 2 * nu0) # Eq. 4 + tnu = R0**2 / (3.0 * (2.0 - gam) ** 2 * nu0) # Eq. 4 T = 1.0 + time / tnu # Above Eq. 4 - eta = (2.5 - gam) / (2. - gam) # Below Eq. 3 - self.sigma_g = (M0 / (2 * np.pi * R0 ** 2.0)) * (2. - gam) * \ - (R0 / r) ** gam * T ** (-eta) * \ - np.exp(-(r / R0) ** (2. - gam) / T) # Eq. 3 + eta = (2.5 - gam) / (2.0 - gam) # Below Eq. 3 + self.sigma_g = ( + (M0 / (2 * np.pi * R0**2.0)) + * (2.0 - gam) + * (R0 / r) ** gam + * T ** (-eta) + * np.exp(-((r / R0) ** (2.0 - gam)) / T) + ) # Eq. 3 self.nu = nu self.gam = gam self.sigma_g[self.sigma_g < self.sigmin] = self.sigmin @@ -110,7 +119,7 @@ def make_disk_from_lbp_alpha(self, M0, R0, alpha, time): """ assert hasattr(self, "aspect_ratio") - self.omega_k = np.sqrt(G * self.M_s / self.r ** 3.0) + self.omega_k = np.sqrt(G * self.M_s / self.r**3.0) self.alpha = alpha r = self.r aspect_ratio = 1 * self.aspect_ratio @@ -123,5 +132,5 @@ def make_disk_from_lbp_alpha(self, M0, R0, alpha, time): self.make_disk_from_lyndenbellpringle(M0, R0, nu0, gam, time) def update_parameters(self): - """ Things that should be called whenever the update_parameters() function is called in the planet class. """ + """Things that should be called whenever the update_parameters() function is called in the planet class.""" self.compute_viscous_evolution() diff --git a/chemcomp/disks/mmsn.py b/chemcomp/disks/mmsn.py index 3cc91b7..e8ec6f8 100755 --- a/chemcomp/disks/mmsn.py +++ b/chemcomp/disks/mmsn.py @@ -18,27 +18,33 @@ class MMSN(Disk): - """ docstring for MMSN. Hayashi 1981, Weidenschilling 1977. Likely outdated""" - - def __init__(self, - defaults, - chemistry, - M_STAR=1 * u.M_sun, - R0=80.0 * u.au, - ALPHA=1e-3, - SIGMA_POWER=-15 / 14, - ASPECT_POWER=2 / 7, - SIGMA_0=1500 * u.g / (u.cm ** 2), - ASPECT_0=0.033, - **kwargs): + """docstring for MMSN. Hayashi 1981, Weidenschilling 1977. Likely outdated""" + + def __init__( + self, + defaults, + chemistry, + M_STAR=1 * u.M_sun, + R0=80.0 * u.au, + ALPHA=1e-3, + SIGMA_POWER=-15 / 14, + ASPECT_POWER=2 / 7, + SIGMA_0=1500 * u.g / (u.cm**2), + ASPECT_0=0.033, + **kwargs + ): super().__init__(defaults, M_STAR, ALPHA, chemistry, **kwargs) - self.sigma_init = SIGMA_0.cgs.value * (self.r / AU) ** SIGMA_POWER * np.exp(- self.r / R0.cgs.value) + self.sigma_init = ( + SIGMA_0.cgs.value + * (self.r / AU) ** SIGMA_POWER + * np.exp(-self.r / R0.cgs.value) + ) self.aspect_ratio = ASPECT_0 * (self.r / AU) ** ASPECT_POWER # Hayashi 1981 self._init_params(**kwargs) self.calc_opacity() def compute_disktmid(self, viscmodel=None): - """calculate the disk temperature. """ + """calculate the disk temperature.""" self.compute_cs_and_hp() self.compute_nu() self.calc_opacity() @@ -49,14 +55,21 @@ def calc_sigma_dust_static(self): This is equal to the pebble surface density, if DTG: peb has been set. Otherwise it contains the whole dust (including also the small grains) """ - vr = 2.0 * self.stokes_number_pebbles / (self.stokes_number_pebbles ** 2. + 1.0) * self.eta * self.v_k + vr = ( + 2.0 + * self.stokes_number_pebbles + / (self.stokes_number_pebbles**2.0 + 1.0) + * self.eta + * self.v_k + ) Pebbleflux = 10e-4 * ME / year - self.sigma_dust = Pebbleflux / (2. * np.pi * self.r * vr) + self.sigma_dust = Pebbleflux / (2.0 * np.pi * self.r * vr) self.chemistry_solid, _ = self._chem.get_composition(self.T) - self.sigma_dust_components = self.chemistry_solid * self.sigma_dust[:, np.newaxis, - np.newaxis] + self.sigma_dust_components = ( + self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis] + ) self.sigma_dust_components[:, 1, 0] = np.zeros_like(self.sigma_g) @@ -64,7 +77,7 @@ def calc_sigma_dust_static(self): return def update_parameters(self): - """ compute the viscous evolution. This function is called whenever the planet calls the update.""" + """compute the viscous evolution. This function is called whenever the planet calls the update.""" self.compute_viscous_evolution() def calc_opacity_oplin(self): diff --git a/chemcomp/disks/twopop.py b/chemcomp/disks/twopop.py index 959db67..be366a1 100755 --- a/chemcomp/disks/twopop.py +++ b/chemcomp/disks/twopop.py @@ -11,20 +11,21 @@ G = const.G.cgs.value + class TwoPopDisk(Disk): - """ LBP, Disk of the twopoppy model.""" + """LBP, Disk of the twopoppy model.""" def __init__( - self, - defaults, - chemistry, - M_STAR=0.7 * u.M_sun, - M0=1e-1 * u.M_sun, - R0=20 * u.au, - time=1 * u.yr, - ALPHA=1e-3, - output={}, - **kwargs + self, + defaults, + chemistry, + M_STAR=0.7 * u.M_sun, + M0=1e-1 * u.M_sun, + R0=20 * u.au, + time=1 * u.yr, + ALPHA=1e-3, + output={}, + **kwargs ): super().__init__(defaults, M_STAR, ALPHA, chemistry, time=time, **kwargs) @@ -34,9 +35,15 @@ def __init__( self.tstar = self.tstar.cgs.value self.compute_disktmid() - self.make_disk_from_lbp_alpha(M0=M0.cgs.value, R0=R0.cgs.value, alpha=ALPHA, time=time.cgs.value) + self.make_disk_from_lbp_alpha( + M0=M0.cgs.value, R0=R0.cgs.value, alpha=ALPHA, time=time.cgs.value + ) - self.sigma_g = self.sigma_g / np.trapz(2 * np.pi * self.r * self.sigma_g, x=self.r) * M0.cgs.value + self.sigma_g = ( + self.sigma_g + / np.trapz(2 * np.pi * self.r * self.sigma_g, x=self.r) + * M0.cgs.value + ) self.compute_disktmid() self.compute_disktmid() @@ -46,16 +53,19 @@ def __init__( def evolve_init(self): pass # DEBUGGING - #mask = np.where(self.r < (5 * u.au).cgs.value)[0] - #self.sigma_g_components[mask, 1, 7] = 1e-60 - #self.sigma_dust_components[mask, :, :] = 1e-60 - #self.compute_sigma_g_from_mol(self.sigma_g_components[:, 1, :]) - #self.compute_sigma_dust_from_mol(self.sigma_dust_components[:, 1, :]) - #self.sigma_dust_components = self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis] - #self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] + # mask = np.where(self.r < (5 * u.au).cgs.value)[0] + # self.sigma_g_components[mask, 1, 7] = 1e-60 + # self.sigma_dust_components[mask, :, :] = 1e-60 + # self.compute_sigma_g_from_mol(self.sigma_g_components[:, 1, :]) + # self.compute_sigma_dust_from_mol(self.sigma_dust_components[:, 1, :]) + # self.sigma_dust_components = self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis] + # self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis] def compute_disktmid(self, viscmodel=None): - self.T = ((0.05 ** 0.25 * self.tstar * (self.r / self.rstar) ** -0.5) ** 4 + (7.) ** 4.0) ** 0.25 + self.T = ( + (0.05**0.25 * self.tstar * (self.r / self.rstar) ** -0.5) ** 4 + + (7.0) ** 4.0 + ) ** 0.25 self.tvisc = np.zeros_like(self.T) self.tirr = np.zeros_like(self.T) self.compute_cs_and_hp() @@ -64,31 +74,34 @@ def compute_disktmid(self, viscmodel=None): def make_disk_from_lyndenbellpringle(self, M0, R0, nu0, gam, time): """ - Make a model of the gas surface density of the disk. - Here the model is set by the Lynden-Bell & Pringle solution. - I follow here Lodato, Scardoni, Manara & Testi (2017) - MNRAS 472, 4700, their Eqs. 2, 3 and 4. - - ARGUMENTS: - M0 = Initial disk mass [g] - R0 = Initial disk radius [cm] - nu0 = Viscosity nu at R0 [cm^2/s] - gam = The gamma coefficient of the LBP model - time = Time after start [s] - """ + Make a model of the gas surface density of the disk. + Here the model is set by the Lynden-Bell & Pringle solution. + I follow here Lodato, Scardoni, Manara & Testi (2017) + MNRAS 472, 4700, their Eqs. 2, 3 and 4. + + ARGUMENTS: + M0 = Initial disk mass [g] + R0 = Initial disk radius [cm] + nu0 = Viscosity nu at R0 [cm^2/s] + gam = The gamma coefficient of the LBP model + time = Time after start [s] + """ r = self.r nu = nu0 * (r / R0) ** gam # Eq. 2 - tnu = R0 ** 2.0 / (3.0 * (2. - gam) ** 2.0 * nu0) # Eq. 4 + tnu = R0**2.0 / (3.0 * (2.0 - gam) ** 2.0 * nu0) # Eq. 4 T = 1.0 + time / tnu # Above Eq. 4 - eta = (2.5 - gam) / (2. - gam) # Below Eq. 3 - self.sigma_g = (M0 / (2.0 * np.pi * R0 ** 2)) * (2. - gam) * \ - (R0 / r) ** gam * T ** (-eta) * \ - np.exp(-(r / R0) ** (2. - gam) / T) # Eq. 3 + eta = (2.5 - gam) / (2.0 - gam) # Below Eq. 3 + self.sigma_g = ( + (M0 / (2.0 * np.pi * R0**2)) + * (2.0 - gam) + * (R0 / r) ** gam + * T ** (-eta) + * np.exp(-((r / R0) ** (2.0 - gam)) / T) + ) # Eq. 3 self.nu = nu self.gam = gam self.sigma_g[self.sigma_g < self.sigmin] = self.sigmin - def make_disk_from_lbp_alpha(self, M0, R0, alpha, time): """ Same as make_disk_from_lyndenbellpringle(), but now assuming a constant @@ -105,7 +118,7 @@ def make_disk_from_lbp_alpha(self, M0, R0, alpha, time): """ assert hasattr(self, "aspect_ratio") - self.omega_k = np.sqrt(G * self.M_s / self.r ** 3.0) + self.omega_k = np.sqrt(G * self.M_s / self.r**3.0) self.alpha = alpha r = self.r aspect_ratio = 1.0 * self.aspect_ratio @@ -114,10 +127,9 @@ def make_disk_from_lbp_alpha(self, M0, R0, alpha, time): gam = 1.0 om1 = self.omega_k[0] cs1 = cs[0] - nu0 = self.alpha * cs1 ** 2.0 / om1 + nu0 = self.alpha * cs1**2.0 / om1 self.make_disk_from_lyndenbellpringle(M0, R0, nu0, gam, time) def update_parameters(self): self.compute_viscous_evolution() - diff --git a/chemcomp/helpers/__init__.py b/chemcomp/helpers/__init__.py index be4bf33..8d4fb3f 100755 --- a/chemcomp/helpers/__init__.py +++ b/chemcomp/helpers/__init__.py @@ -1,10 +1,17 @@ import os from ast import literal_eval -CONFIG_PATH = os.path.join(os.getcwd(), "config") # set default path of config directory +CONFIG_PATH = os.path.join( + os.getcwd(), "config" +) # set default path of config directory JOB_PATH = os.path.join(os.getcwd(), "jobs") # set default path of config directory -OUTPUT_PATH = os.path.join(os.getcwd(), "output") # set default path of output directory -ZIP_OUTPUT_PATH = os.path.join(os.getcwd(), "zip_output") # set default path of output directory +OUTPUT_PATH = os.path.join( + os.getcwd(), "output" +) # set default path of output directory +ZIP_OUTPUT_PATH = os.path.join( + os.getcwd(), "zip_output" +) # set default path of output directory + def eval_kwargs(inp): try: @@ -20,4 +27,3 @@ def eval_kwargs(inp): from .main_helpers import run_setup, run_pipeline from .solvediffonedee import solvediffonedee_components, getfluxonedee from .analysis_helper import GrowthPlot, molecules, elements - diff --git a/chemcomp/helpers/analysis_helper.py b/chemcomp/helpers/analysis_helper.py index b6ba31f..8e60629 100755 --- a/chemcomp/helpers/analysis_helper.py +++ b/chemcomp/helpers/analysis_helper.py @@ -12,9 +12,27 @@ np.seterr(divide="ignore", invalid="ignore") elements = ["C", "O", "Fe", "S", "Mg", "Si", "Na", "K", "N", "Al", "Ti", "V", "H", "He"] -molecules = ["rest", "CO", "N2", "CH4", "CO2", "NH3", "H2S", "H2O", "Fe3O4", "C", "FeS", "NaAlSi3O8", "KAlSi3O8", - "Mg2SiO4", - "Fe2O3", "VO", "MgSiO3", "Al2O3", "TiO"] +molecules = [ + "rest", + "CO", + "N2", + "CH4", + "CO2", + "NH3", + "H2S", + "H2O", + "Fe3O4", + "C", + "FeS", + "NaAlSi3O8", + "KAlSi3O8", + "Mg2SiO4", + "Fe2O3", + "VO", + "MgSiO3", + "Al2O3", + "TiO", +] FeH_dict = { "FeH": [-0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4], @@ -54,8 +72,22 @@ MassC = 12.0 MassV = 50.9415 -el_masses = {"C": MassC, "O": MassO, "H": MassH, "Fe": MassFe, "Mg": MassMg, "Si": MassSi, "N": MassN, "S": MassS, - "He": MassHe, "Al": MassAl, "Ti": MassTi, "K": MassK, "Na": MassNa, "V": MassV} +el_masses = { + "C": MassC, + "O": MassO, + "H": MassH, + "Fe": MassFe, + "Mg": MassMg, + "Si": MassSi, + "N": MassN, + "S": MassS, + "He": MassHe, + "Al": MassAl, + "Ti": MassTi, + "K": MassK, + "Na": MassNa, + "V": MassV, +} DOT_SPACE = 1e5 * u.yr # time in years between dots DOT_SPACE_LARGE = 5e5 * u.yr # time in years between dots @@ -75,31 +107,63 @@ ("densely dashdotdotted", (0, (3, 1, 1, 1, 1, 1))), ] -dict_of_planet_read_names = {"t": r"t", "M": r"M", "M_a": r"M_a", - "M_c": r"M_c", "a_p": r"a_p", "r_p": r"R_p", "tau_m": r"\tau_m", "T": "T", - "comp_a": "Atmospheric Content", "comp_c": "Core Content", "sigma_g": r"\Sigma_g", - "gamma_tot": r"\Gamma", "gamma_norm": r"\Gamma / \Gamma_0", "fSigma": "fSigma"} -dict_of_acc_names_template = {r"m_dot": r"\dot M_{\mathrm{", r"m_dot_a": r"$\dot M_{\mathrm{atm,", - r"m_dot_c": r"\dot M_{\mathrm{core,", - r"m_dot_a_chem": r"\dot M_{\mathrm{atm,", r"m_dot_c_chem": r"\dot M_{\mathrm{core,", - r"M_z": r"M_{Z,\mathrm{", r"regime":"regime_"} -dict_of_acc_names = {f"{k}_{acc}": r"{}{}".format(v, acc) + r"}}" for k, v in dict_of_acc_names_template.items() for - acc in ["peb", "gas"]} - -dict_of_extra_names = {"Z_vs_M": r"\frac{M_{Z}}{M}}", - "Z_vs_M_a": r"\frac{M_{a, Z}}{M_a}}", "Z_vs_M_c": r"\frac{M_{c, Z}}{M_c}}"} +dict_of_planet_read_names = { + "t": r"t", + "M": r"M", + "M_a": r"M_a", + "M_c": r"M_c", + "a_p": r"a_p", + "r_p": r"R_p", + "tau_m": r"\tau_m", + "T": "T", + "comp_a": "Atmospheric Content", + "comp_c": "Core Content", + "sigma_g": r"\Sigma_g", + "gamma_tot": r"\Gamma", + "gamma_norm": r"\Gamma / \Gamma_0", + "fSigma": "fSigma", +} +dict_of_acc_names_template = { + r"m_dot": r"\dot M_{\mathrm{", + r"m_dot_a": r"$\dot M_{\mathrm{atm,", + r"m_dot_c": r"\dot M_{\mathrm{core,", + r"m_dot_a_chem": r"\dot M_{\mathrm{atm,", + r"m_dot_c_chem": r"\dot M_{\mathrm{core,", + r"M_z": r"M_{Z,\mathrm{", + r"regime": "regime_", +} +dict_of_acc_names = { + f"{k}_{acc}": r"{}{}".format(v, acc) + r"}}" + for k, v in dict_of_acc_names_template.items() + for acc in ["peb", "gas"] +} + +dict_of_extra_names = { + "Z_vs_M": r"\frac{M_{Z}}{M}}", + "Z_vs_M_a": r"\frac{M_{a, Z}}{M_a}}", + "Z_vs_M_c": r"\frac{M_{c, Z}}{M_c}}", +} NAME = {**dict_of_planet_read_names, **dict_of_acc_names, **dict_of_extra_names} class GrowthPlot: - def __init__(self, files, out, plotlabels=None, FeH=0, close_plots=False, quantities=None, dry_run=False): + def __init__( + self, + files, + out, + plotlabels=None, + FeH=0, + close_plots=False, + quantities=None, + dry_run=False, + ): self.files = files self.out = out self._dry_run = dry_run if quantities is not None: self._used_quantities = set(quantities) - self._used_quantities.update(["M","t"]) + self._used_quantities.update(["M", "t"]) else: self._used_quantities = quantities self.import_data(files) @@ -130,15 +194,27 @@ def solar_from_FeH(FeH): NaHabu = 1 * NaHabu0 VHabu = 1 * VHabu0 - return {"O": OHabu * MassO, "Si": SiHabu * MassSi, "C": CHabu * MassC, "S": SHabu * MassS, - "Fe": FeHabu * MassFe, "Mg": MgHabu * MassMg, "H": 1, "He": 0.085 * MassHe, "Na": NaHabu * MassNa, - "N": NHabu * MassN, "Al": AlHabu * MassAl, "Ti": TiHabu * MassTi, "K": KHabu * MassK, - "V": MassV * VHabu} + return { + "O": OHabu * MassO, + "Si": SiHabu * MassSi, + "C": CHabu * MassC, + "S": SHabu * MassS, + "Fe": FeHabu * MassFe, + "Mg": MgHabu * MassMg, + "H": 1, + "He": 0.085 * MassHe, + "Na": NaHabu * MassNa, + "N": NHabu * MassN, + "Al": AlHabu * MassAl, + "Ti": TiHabu * MassTi, + "K": KHabu * MassK, + "V": MassV * VHabu, + } self.FeH = {} for o in self.out: if type(FeH) == dict: - """ dict needs to have out as keys """ + """dict needs to have out as keys""" self.data[o]["solar"] = solar_from_FeH(FeH[o]) self.FeH[o] = FeH[o] else: @@ -154,8 +230,11 @@ def import_data(self, files): with h5py.File(files[i], "r") as f: units = dict(f["/planet/units"].attrs) keys = list(f["/planet"].keys()) - keys = [k for k in keys if - k in self._used_quantities] if self._used_quantities is not None else keys + keys = ( + [k for k in keys if k in self._used_quantities] + if self._used_quantities is not None + else keys + ) if "units" in keys: keys.remove("units") break @@ -175,18 +254,25 @@ def import_data(self, files): data[k] = u.Quantity(d, units.get(k)) else: data[k] = d - data = {q: data[q][np.where(~np.isnan(data["M"]))] for q in data.keys()} + data = { + q: data[q][np.where(~np.isnan(data["M"]))] + for q in data.keys() + } self.data[self.out[i]] = data except OSError: self.faulty_list.append(i) if len(self.faulty_list) > 0: - print(f"WARNING: There are in total {len(self.faulty_list)} corrupt files!") + print( + f"WARNING: There are in total {len(self.faulty_list)} corrupt files!" + ) else: with h5py.File(files[0], "r") as f: fill_data = {} for k in keys: - fill_data[k] = np.squeeze(u.Quantity(np.array(f["/planet"][k]), units.get(k))) + fill_data[k] = np.squeeze( + u.Quantity(np.array(f["/planet"][k]), units.get(k)) + ) for i in range(len(self.out)): data = {} @@ -195,8 +281,10 @@ def import_data(self, files): data = {q: data[q] for q in data.keys()} self.data[self.out[i]] = data - self.out = [self.out[i] for i in range(len(self.out)) if i not in self.faulty_list] - self.dt = (self.data[self.out[0]]["t"][-1] - self.data[self.out[0]]["t"][-2]) + self.out = [ + self.out[i] for i in range(len(self.out)) if i not in self.faulty_list + ] + self.dt = self.data[self.out[0]]["t"][-1] - self.data[self.out[0]]["t"][-2] def composite_data(self): def add_to_data(fct, quantity, comp_NAME, **kwargs): @@ -206,23 +294,44 @@ def add_to_data(fct, quantity, comp_NAME, **kwargs): NAME[quantity] = comp_NAME def M_Z_vs_M_a(o): - return np.sum(self.data[o]["comp_a"][:, 0, :-7], axis=1) / self.data[o]["M_a"] * u.dimensionless_unscaled + return ( + np.sum(self.data[o]["comp_a"][:, 0, :-7], axis=1) + / self.data[o]["M_a"] + * u.dimensionless_unscaled + ) def M_Z_vs_M(o): - return np.sum(self.data[o]["comp_a"][:, 0, :-7] + self.data[o]["comp_c"][:, 0, :-7], axis=1) / self.data[o][ - "M"] * u.dimensionless_unscaled + return ( + np.sum( + self.data[o]["comp_a"][:, 0, :-7] + + self.data[o]["comp_c"][:, 0, :-7], + axis=1, + ) + / self.data[o]["M"] + * u.dimensionless_unscaled + ) def M_Z(o): - return np.sum(self.data[o]["comp_a"][:, 0, :-7] + self.data[o]["comp_c"][:, 0, :-7], axis=1) + return np.sum( + self.data[o]["comp_a"][:, 0, :-7] + self.data[o]["comp_c"][:, 0, :-7], + axis=1, + ) def M_Z_vs_M_c(o): - return np.sum(self.data[o]["comp_c"][:, 0, :-7], axis=1) / self.data[o]["M_c"] * u.dimensionless_unscaled + return ( + np.sum(self.data[o]["comp_c"][:, 0, :-7], axis=1) + / self.data[o]["M_c"] + * u.dimensionless_unscaled + ) def elem_a(o): return dict( zip( elements, - (self.data[o]["comp_a"][:, 0, :] / self.data[o]["M_a"][:, np.newaxis]).T + ( + self.data[o]["comp_a"][:, 0, :] + / self.data[o]["M_a"][:, np.newaxis] + ).T, ) ) @@ -230,7 +339,10 @@ def elem_c(o): return dict( zip( elements, - (self.data[o]["comp_c"][:, 0, :] / self.data[o]["M_c"][:, np.newaxis]).T + ( + self.data[o]["comp_c"][:, 0, :] + / self.data[o]["M_c"][:, np.newaxis] + ).T, ) ) @@ -238,8 +350,13 @@ def elem_tot(o): return dict( zip( elements, - ((self.data[o]["comp_c"][:, 0, :] + self.data[o]["comp_a"][:, 0, :]) / self.data[o]["M"][:, - np.newaxis]).T + ( + ( + self.data[o]["comp_c"][:, 0, :] + + self.data[o]["comp_a"][:, 0, :] + ) + / self.data[o]["M"][:, np.newaxis] + ).T, ) ) @@ -247,7 +364,10 @@ def mol_a(o): return dict( zip( molecules, - (self.data[o]["comp_a"][:, 1, :] / self.data[o]["M_a"][:, np.newaxis]).T + ( + self.data[o]["comp_a"][:, 1, :] + / self.data[o]["M_a"][:, np.newaxis] + ).T, ) ) @@ -255,7 +375,10 @@ def mol_c(o): return dict( zip( molecules, - (self.data[o]["comp_c"][:, 1, :] / self.data[o]["M_c"][:, np.newaxis]).T + ( + self.data[o]["comp_c"][:, 1, :] + / self.data[o]["M_c"][:, np.newaxis] + ).T, ) ) @@ -263,28 +386,53 @@ def mol_tot(o): return dict( zip( molecules, - ((self.data[o]["comp_c"][:, 1, :] + self.data[o]["comp_a"][:, 1, :]) / self.data[o]["M"][:, - np.newaxis]).T + ( + ( + self.data[o]["comp_c"][:, 1, :] + + self.data[o]["comp_a"][:, 1, :] + ) + / self.data[o]["M"][:, np.newaxis] + ).T, ) ) def XH(o, X, orig): - return self.data[o][f"elem_{orig}"][X] / self.data[o][f"elem_{orig}"]["H"] / self.data[o]["solar"][X] + return ( + self.data[o][f"elem_{orig}"][X] + / self.data[o][f"elem_{orig}"]["H"] + / self.data[o]["solar"][X] + ) def XY(o, X, Y, orig): - return self.data[o][f"elem_{orig}"][X] / self.data[o][f"elem_{orig}"][Y] / self.data[o]["solar"][X] * \ - self.data[o]["solar"][Y] + return ( + self.data[o][f"elem_{orig}"][X] + / self.data[o][f"elem_{orig}"][Y] + / self.data[o]["solar"][X] + * self.data[o]["solar"][Y] + ) def XH_not_normed(o, X, orig): - return self.data[o][f"elem_{orig}"][X] / self.data[o][f"elem_{orig}"]["H"] / el_masses[X] + return ( + self.data[o][f"elem_{orig}"][X] + / self.data[o][f"elem_{orig}"]["H"] + / el_masses[X] + ) def XY_not_normed(o, X, Y, orig): - return self.data[o][f"elem_{orig}"][X] / self.data[o][f"elem_{orig}"][Y] * el_masses[Y] / el_masses[X] + return ( + self.data[o][f"elem_{orig}"][X] + / self.data[o][f"elem_{orig}"][Y] + * el_masses[Y] + / el_masses[X] + ) def pebiso(o): - iso = np.where(self.data[o]["m_dot_gas"].value > 0, np.ones_like(self.data[o]["t"].value), - np.zeros_like(self.data[o]["t"].value)) - iso[np.argmax(iso>0):] = 1 + iso = np.where( + self.data[o]["m_dot_gas"].value > 0, + np.ones_like(self.data[o]["t"].value), + np.zeros_like(self.data[o]["t"].value), + ) + iso[np.argmax(iso > 0) :] = 1 iso = iso * u.dimensionless_unscaled return iso @@ -341,8 +489,10 @@ def convert_units(self): def get_label(self, n, quantity): if hasattr(quantity, "unit") and quantity.unit != u.dimensionless_unscaled: x = r"${}$ / {:latex}".format(n, quantity.unit) - elif hasattr(self.data[self.out[0]].get(n, None), "unit") and self.data[self.out[0]][ - n].unit != u.dimensionless_unscaled: + elif ( + hasattr(self.data[self.out[0]].get(n, None), "unit") + and self.data[self.out[0]][n].unit != u.dimensionless_unscaled + ): quantity = self.data[self.out[0]][n] x = r"${}$ / {:latex}".format(n, quantity.unit) elif n == "a_p": @@ -360,20 +510,20 @@ def get_label(self, n, quantity): return x def plot_data( - self, - quantities, - out=None, - sub_quant=None, - t_dot=None, - label=None, - ulabel=True, - fig=None, - units=None, - ax=None, - plt=None, - color_dict=None, - use_ls=True, - **kwargs + self, + quantities, + out=None, + sub_quant=None, + t_dot=None, + label=None, + ulabel=True, + fig=None, + units=None, + ax=None, + plt=None, + color_dict=None, + use_ls=True, + **kwargs, ): snapshot = kwargs.get("snapshot", None) out = self.out if out is None else out @@ -395,17 +545,25 @@ def plot_data( use_t_dot = False if snapshot is None: + def get_dots(x_data, y_data, t): - t_at_dot = np.arange(dotspace.cgs.value, np.max(t[:-1]).cgs.value, dotspace.cgs.value) + t_at_dot = np.arange( + dotspace.cgs.value, np.max(t[:-1]).cgs.value, dotspace.cgs.value + ) x = np.interp(t_at_dot, t.cgs.value, x_data) y = np.interp(t_at_dot, t.cgs.value, y_data) - t_at_dot = np.arange(dotspace_large.cgs.value, np.max(t[:-1]).cgs.value, dotspace_large.cgs.value) + t_at_dot = np.arange( + dotspace_large.cgs.value, + np.max(t[:-1]).cgs.value, + dotspace_large.cgs.value, + ) x_large = np.interp(t_at_dot, t.cgs.value, x_data) y_large = np.interp(t_at_dot, t.cgs.value, y_data) - x, y = u.Quantity([xi for xi in x if xi not in x_large], x.unit), u.Quantity( - [yi for yi in y if yi not in y_large], y.unit) + x, y = u.Quantity( + [xi for xi in x if xi not in x_large], x.unit + ), u.Quantity([yi for yi in y if yi not in y_large], y.unit) return x, y, x_large, y_large if not hasattr(self, "_plot_x_data"): @@ -434,63 +592,117 @@ def get_dots(x_data, y_data, t): legend = label if units is not None: - self._plot_x_data[ax][o], self._plot_y_data[ax][o] = planet[quantities[0]].to(units[0]), planet[ - quantities[1]].to(units[1]) + self._plot_x_data[ax][o], self._plot_y_data[ax][o] = planet[ + quantities[0] + ].to(units[0]), planet[quantities[1]].to(units[1]) else: - self._plot_x_data[ax][o], self._plot_y_data[ax][o] = planet[quantities[0]], planet[quantities[1]] + self._plot_x_data[ax][o], self._plot_y_data[ax][o] = ( + planet[quantities[0]], + planet[quantities[1]], + ) if ax is not None: - if use_ls: - ls_mask = self.data[o]["pebiso"]==0 + ls_mask = self.data[o]["pebiso"] == 0 not_ls_mask = np.logical_not(ls_mask) - not_ls_mask[np.argmin(ls_mask)-1] = True - self._plot[ax][o], = ax.loglog(self._plot_x_data[ax][o][not_ls_mask], - self._plot_y_data[ax][o][not_ls_mask], label=legend, ls=":", - **kwargs) + not_ls_mask[np.argmin(ls_mask) - 1] = True + (self._plot[ax][o],) = ax.loglog( + self._plot_x_data[ax][o][not_ls_mask], + self._plot_y_data[ax][o][not_ls_mask], + label=legend, + ls=":", + **kwargs, + ) else: ls_mask = np.arange(len(self._plot_x_data[ax][o])) - self._plot[ax][o], = ax.loglog(self._plot_x_data[ax][o][ls_mask], self._plot_y_data[ax][o][ls_mask], label=legend, - **kwargs) - - + (self._plot[ax][o],) = ax.loglog( + self._plot_x_data[ax][o][ls_mask], + self._plot_y_data[ax][o][ls_mask], + label=legend, + **kwargs, + ) - ax.set_xlabel(self.get_label(NAME.get(quantities[0], quantities[0]), self._plot_x_data[ax][o])) - ax.set_ylabel(self.get_label(NAME.get(quantities[1], quantities[1]), self._plot_y_data[ax][o])) + ax.set_xlabel( + self.get_label( + NAME.get(quantities[0], quantities[0]), + self._plot_x_data[ax][o], + ) + ) + ax.set_ylabel( + self.get_label( + NAME.get(quantities[1], quantities[1]), + self._plot_y_data[ax][o], + ) + ) if use_t_dot: x, y, x_large, y_large = get_dots( planet[quantities[0]], planet[quantities[1]], planet["t"] ) if len(x) > 0: - ax.scatter(x, y, marker="o", s=10, color=color_t_dot_small[:len(x)], alpha=1) + ax.scatter( + x, + y, + marker="o", + s=10, + color=color_t_dot_small[: len(x)], + alpha=1, + ) if len(x_large) > 0: - ax.scatter(x_large, y_large, s=25, marker="o", color=color_t_dot[:len(x_large)], alpha=1) + ax.scatter( + x_large, + y_large, + s=25, + marker="o", + color=color_t_dot[: len(x_large)], + alpha=1, + ) elif plt is not None: - - self._plot[o], = plt.loglog( - self._plot_x_data[ax][o], self._plot_y_data[ax][o], label=legend, **kwargs + (self._plot[o],) = plt.loglog( + self._plot_x_data[ax][o], + self._plot_y_data[ax][o], + label=legend, + **kwargs, ) - plt.xlabel(self.get_label(NAME.get(quantities[0], quantities[0]), self._plot_x_data[ax][o])) - plt.ylabel(self.get_label(NAME.get(quantities[1], quantities[1]), self._plot_y_data[ax][o])) + plt.xlabel( + self.get_label( + NAME.get(quantities[0], quantities[0]), + self._plot_x_data[ax][o], + ) + ) + plt.ylabel( + self.get_label( + NAME.get(quantities[1], quantities[1]), + self._plot_y_data[ax][o], + ) + ) plt.title( - self.get_label(NAME.get(quantities[1], quantities[1]), NAME.get(quantities[0], quantities[0]))) + self.get_label( + NAME.get(quantities[1], quantities[1]), + NAME.get(quantities[0], quantities[0]), + ) + ) if use_t_dot: x, y, x_large, y_large = get_dots( planet[quantities[0]], planet[quantities[1]], planet["t"] ) plt.plot(x, y, "ko", markersize=4, color="gray", alpha=1) - plt.plot(x_large, y_large, "ko", markersize=8, color="gray", alpha=1) + plt.plot( + x_large, y_large, "ko", markersize=8, color="gray", alpha=1 + ) else: for o in out: if o in self.out: snapshot = np.minimum(snapshot, len(self._plot_x_data[ax][o])) - self._plot[ax][o].set_data(self._plot_x_data[ax][o][:snapshot], self._plot_y_data[ax][o][:snapshot]) + self._plot[ax][o].set_data( + self._plot_x_data[ax][o][:snapshot], + self._plot_y_data[ax][o][:snapshot], + ) return self._plot[ax] @@ -515,24 +727,70 @@ def apply_mask(self, mask): self.out = mask self.plotlabels = {k: v for (k, v) in self.plotlabels.items() if k in mask} - def pqplot_video(self, plt, quantity, frames, lines=[], sub_qaunt=None, log_data=False, units=None, **kwargs): - """ not working anymore, needs some more fixes""" + def pqplot_video( + self, + plt, + quantity, + frames, + lines=[], + sub_qaunt=None, + log_data=False, + units=None, + **kwargs, + ): + """not working anymore, needs some more fixes""" fig, ax = plt.subplots(sharex=True, sharey=True) div = make_axes_locatable(ax) - cax = div.append_axes('right', '5%', '5%') + cax = div.append_axes("right", "5%", "5%") def animate(i): [ax.clear() for ax in fig.axes] - self.pqplot(plt=plt, snapshot=i, fig=fig, quantity=quantity, lines=lines, sub_qaunt=sub_qaunt, - log_data=log_data, cax=cax, units=units, **kwargs) + self.pqplot( + plt=plt, + snapshot=i, + fig=fig, + quantity=quantity, + lines=lines, + sub_qaunt=sub_qaunt, + log_data=log_data, + cax=cax, + units=units, + **kwargs, + ) return animation.FuncAnimation(fig, animate, frames=frames) - def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=False, symlog_data=False, snapshot=-1, - units=None, icelines=None, ax=None, shape=None, cbar_location=None, - **kwargs): - data_all, extend_r, p_data, p_name, q_data, q_name, r_list, r_name, x_data_all, y_data_all, z_data, z_data_all = self.extract_r0t0_data( - quantity, snapshot, sub_quant, units) + def pqplot( + self, + fig, + quantity, + lines=[], + cax=None, + sub_quant=None, + log_data=False, + symlog_data=False, + snapshot=-1, + units=None, + icelines=None, + ax=None, + shape=None, + cbar_location=None, + **kwargs, + ): + ( + data_all, + extend_r, + p_data, + p_name, + q_data, + q_name, + r_list, + r_name, + x_data_all, + y_data_all, + z_data, + z_data_all, + ) = self.extract_r0t0_data(quantity, snapshot, sub_quant, units) if shape is None: if len(extend_r) > 2: @@ -569,20 +827,32 @@ def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=Fal for l in lines: if isinstance(l[0], str): Z = np.array( - [d[l[0]].value[snapshot if snapshot < len(d["t"]) else -1] for d in data_all[i].values()]) + [ + d[l[0]].value[snapshot if snapshot < len(d["t"]) else -1] + for d in data_all[i].values() + ] + ) else: Z = l[0] - zi, yi, xi = np.histogram2d(p_data, q_data, bins=(len(np.unique(p_data)), len(np.unique(q_data))), - weights=Z, - normed=False) + zi, yi, xi = np.histogram2d( + p_data, + q_data, + bins=(len(np.unique(p_data)), len(np.unique(q_data))), + weights=Z, + normed=False, + ) # zi = gaussian_filter(zi.T, 0.15) x_step = l[1].pop("x_step", 1) y_step = l[1].pop("y_step", 1) zi = zi.T zi = np.ma.masked_equal(zi, np.inf) - CS = ax_temp.contour(x_data_all[i][::x_step, ::y_step], y_data_all[i][::x_step, ::y_step], - zi[::x_step, ::y_step], **l[1]) + CS = ax_temp.contour( + x_data_all[i][::x_step, ::y_step], + y_data_all[i][::x_step, ::y_step], + zi[::x_step, ::y_step], + **l[1], + ) ax_temp.clabel(CS, l[2]) l[1]["x_step"] = x_step l[1]["y_step"] = y_step @@ -593,11 +863,9 @@ def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=Fal time = 3 * u.Myr if snapshot == -1 else snapshot * self.dt if r_name is not None: if r_name == "case": - ax_temp.set_title( - r"{}".format(r_list[i])) + ax_temp.set_title(r"{}".format(r_list[i])) else: - ax_temp.set_title( - r"{}={}".format(r_name, r_list[i])) + ax_temp.set_title(r"{}={}".format(r_name, r_list[i])) shrink = 0.4 if cbar_location == "bottom" else 0.6 aspect = 6 if cbar_location == "bottom" else 20 @@ -607,7 +875,9 @@ def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=Fal cbar = fig.colorbar(im, ax=ax[:-1], cax=cax, shrink=shrink, aspect=aspect) ax.flat[-1].remove() else: - cbar = fig.colorbar(im, ax=ax, cax=cax, location=cbar_location, shrink=shrink, aspect=aspect) + cbar = fig.colorbar( + im, ax=ax, cax=cax, location=cbar_location, shrink=shrink, aspect=aspect + ) try: for temp_ax in ax.flat: @@ -615,15 +885,17 @@ def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=Fal except Exception as e: print(f"no label outer possible {e}") - cbar.ax.set_title(self.get_label(quantity, z_data), y=1.03) fig.suptitle( - r"${}$-${}$ Plot for {}, t = {}".format(p_name, q_name, quantity, time.to(u.Myr))) + r"${}$-${}$ Plot for {}, t = {}".format( + p_name, q_name, quantity, time.to(u.Myr) + ) + ) return im def extract_r0t0_data(self, quantity, snapshot, sub_quant, units): - """ data in pipline should look like: r0 = [20, 10, 5, 20, 10, 5]), t0 = [1e5, 1e5, 1e5, 1e6, 1e6, 1e6]) - in other words: flattened meshgrids""" + """data in pipline should look like: r0 = [20, 10, 5, 20, 10, 5]), t0 = [1e5, 1e5, 1e5, 1e6, 1e6, 1e6]) - in other words: flattened meshgrids""" format_string = "{}:{}, {}:{}" format_string_long = "{}:{}, {}:{}, {}:{}" if parse.parse(format_string_long, self.plotlabels[self.out[0]]) is not None: @@ -637,7 +909,10 @@ def extract_r0t0_data(self, quantity, snapshot, sub_quant, units): if np.all(unique == np.array(["evaporation", "plain"])): unique = unique[::-1] - extend_r = [[self.out[i] for i in range(len(self.out)) if pa[i] == lsta] for lsta in unique] + extend_r = [ + [self.out[i] for i in range(len(self.out)) if pa[i] == lsta] + for lsta in unique + ] elif parse.parse(format_string, self.plotlabels[self.out[0]]) is not None: # -> case of only par_1, par_2 extend_r = [self.out] @@ -650,9 +925,13 @@ def extract_r0t0_data(self, quantity, snapshot, sub_quant, units): p_data, q_data = [], [] for o in sub: if len(extend_r) > 1: - r_name, r, p_name, p, q_name, q = parse.parse(format_string_long, self.plotlabels[o]) + r_name, r, p_name, p, q_name, q = parse.parse( + format_string_long, self.plotlabels[o] + ) else: - p_name, p, q_name, q = parse.parse(format_string, self.plotlabels[o]) + p_name, p, q_name, q = parse.parse( + format_string, self.plotlabels[o] + ) r_name = None r = 0 p_data.append(float(p)) @@ -665,20 +944,35 @@ def extract_r0t0_data(self, quantity, snapshot, sub_quant, units): q_data = (q_data * units[1]).cgs.to(units[1]) if sub_quant is not None: - z_data = np.array([data[quantity][sub_quant][snapshot if snapshot < len(data["t"]) else -1] for data in - data.values()]) + z_data = np.array( + [ + data[quantity][sub_quant][ + snapshot if snapshot < len(data["t"]) else -1 + ] + for data in data.values() + ] + ) else: z_data = np.array( - [data[quantity].value[snapshot if snapshot < len(data["t"]) else -1] for data in - data.values()]) + [ + data[quantity].value[ + snapshot if snapshot < len(data["t"]) else -1 + ] + for data in data.values() + ] + ) if units is not None: z_data = z_data.to(units[2]) - zi, yi, xi = np.histogram2d(p_data, q_data, bins=(len(np.unique(p_data)), len(np.unique(q_data))), - weights=z_data, - normed=False) + zi, yi, xi = np.histogram2d( + p_data, + q_data, + bins=(len(np.unique(p_data)), len(np.unique(q_data))), + weights=z_data, + normed=False, + ) zi = np.ma.masked_equal(zi, 0) zi = np.ma.masked_equal(zi, np.inf) X, Y = np.meshgrid(np.unique(p_data), np.unique(q_data)) @@ -689,35 +983,86 @@ def extract_r0t0_data(self, quantity, snapshot, sub_quant, units): data_all.append(data) x_data_all = np.array(x_data_all) y_data_all = np.array(y_data_all) - return data_all, extend_r, p_data, p_name, q_data, q_name, r_list, r_name, x_data_all, y_data_all, z_data, z_data_all + return ( + data_all, + extend_r, + p_data, + p_name, + q_data, + q_name, + r_list, + r_name, + x_data_all, + y_data_all, + z_data, + z_data_all, + ) def get_iceline_background(self, file, ax): def get_data_iceline(): data = {} - with tables.open_file(file, mode='r') as f: + with tables.open_file(file, mode="r") as f: data["T"] = np.array(f.root.disk.T) data["r"] = np.squeeze(np.array(f.root.disk.r) / const.au.cgs.value) - data["t"] = np.squeeze(np.array(f.root.disk.t)) / 1e6 / 365.25 / 24 / 3600 + data["t"] = ( + np.squeeze(np.array(f.root.disk.t)) / 1e6 / 365.25 / 24 / 3600 + ) - data["icelines"] = np.array([20, 30, 70, 150, 371, 970, 1354, 1423, 1500, 1653, 2000]) + data["icelines"] = np.array( + [20, 30, 70, 150, 371, 970, 1354, 1423, 1500, 1653, 2000] + ) data["iceline_pos"] = np.array( - [np.searchsorted(-data["T"][i], -data["icelines"]) for i in range(len(data["T"]))]) - data["iceline_r"] = np.array([data["r"][data["iceline_pos"][i]] for i in range(len(data["t"]))]) - data["iceline_name"] = ["CO & N2", "CH4", "CO2", "H2O & H2S", "Fe3O4", "NaAlSi3O8 & KAlSi3O8", - "Mg2SiO4 & Fe2O3", - "VO", "MgSiO3", "Al2O3", "TiO"] + [ + np.searchsorted(-data["T"][i], -data["icelines"]) + for i in range(len(data["T"])) + ] + ) + data["iceline_r"] = np.array( + [data["r"][data["iceline_pos"][i]] for i in range(len(data["t"]))] + ) + data["iceline_name"] = [ + "CO & N2", + "CH4", + "CO2", + "H2O & H2S", + "Fe3O4", + "NaAlSi3O8 & KAlSi3O8", + "Mg2SiO4 & Fe2O3", + "VO", + "MgSiO3", + "Al2O3", + "TiO", + ] return data data = get_data_iceline() xlim, ylim = ax.get_xlim(), ax.get_ylim() texts = [] for i in range(len(data["icelines"])): - if np.any(np.logical_and(data["iceline_r"][:, i] > xlim[0], data["iceline_r"][:, i] < xlim[1])): + if np.any( + np.logical_and( + data["iceline_r"][:, i] > xlim[0], data["iceline_r"][:, i] < xlim[1] + ) + ): ymask = np.logical_and(data["t"] > ylim[0], data["t"] < ylim[1]) - ax.plot(data["iceline_r"][ymask, i], data["t"][ymask], ls="--", color="gray", alpha=0.8) - t = ax.text(np.mean(data["iceline_r"][ymask, i]), 0.5 * (np.mean(ax.get_ylim())), - data["iceline_name"][i], fontsize=10, - rotation=90, ha="center", va="center", color="white", backgroundcolor="gray") + ax.plot( + data["iceline_r"][ymask, i], + data["t"][ymask], + ls="--", + color="gray", + alpha=0.8, + ) + t = ax.text( + np.mean(data["iceline_r"][ymask, i]), + 0.5 * (np.mean(ax.get_ylim())), + data["iceline_name"][i], + fontsize=10, + rotation=90, + ha="center", + va="center", + color="white", + backgroundcolor="gray", + ) t.set_bbox(dict(alpha=0.0)) texts.append(t) diff --git a/chemcomp/helpers/import_config.py b/chemcomp/helpers/import_config.py index e72d76d..bb64703 100755 --- a/chemcomp/helpers/import_config.py +++ b/chemcomp/helpers/import_config.py @@ -1,5 +1,6 @@ import astropy.units as u from yaml import load, dump + try: from yaml import CLoader as Loader, CDumper as Dumper except ImportError: @@ -7,7 +8,7 @@ def import_config(config_file): - """Wrapper that imports the config. """ + """Wrapper that imports the config.""" stream = open(config_file, "r") config = load(stream, Loader=Loader) @@ -17,7 +18,7 @@ def import_config(config_file): def scale(dictionary, quantity, unit=1): - """ Helper function to scale dictionary with unit """ + """Helper function to scale dictionary with unit""" if dictionary is not None: if isinstance(dictionary, list): for item in dictionary: @@ -31,15 +32,15 @@ def scale(dictionary, quantity, unit=1): def scale_units(config): - """Add a unit to all unit based quantities from the configuration file. """ + """Add a unit to all unit based quantities from the configuration file.""" planets = config.get("config_planet", None) scale(planets, "M_c", u.earthMass) scale(planets, "M_a", u.earthMass) scale(planets, "a_p", u.au) scale(planets, "t_0", u.Myr) scale(planets, "r_p", u.jupiterRad) - scale(planets, "rho_c", u.g / (u.cm ** 3)) - scale(planets, "rho_0", u.g / (u.cm ** 3)) + scale(planets, "rho_c", u.g / (u.cm**3)) + scale(planets, "rho_0", u.g / (u.cm**3)) scale(planets, "dt", u.yr) scale(planets, "M_end", u.jupiterMass) scale(planets, "r_in", u.au) @@ -50,7 +51,7 @@ def scale_units(config): scale(disk, "M0", u.M_sun) scale(disk, "R0", u.au) scale(disk, "time", u.Myr) - scale(disk, "SIGMA_0", u.g / (u.cm ** 2)) + scale(disk, "SIGMA_0", u.g / (u.cm**2)) scale(disk, "a_0", u.cm) scale(disk, "FUV_MDOT", u.M_sun / u.yr) scale(disk, "tau_disk", u.yr) @@ -60,11 +61,11 @@ def scale_units(config): pebble_accretion = config.get("config_pebble_accretion", None) scale(pebble_accretion, "M_DOT_P", u.earthMass / (u.Myr)) scale(pebble_accretion, "R_peb", u.cm) - scale(pebble_accretion, "rho_solid", u.g / (u.cm ** 3)) + scale(pebble_accretion, "rho_solid", u.g / (u.cm**3)) scale(pebble_accretion, "u_frag", u.m / u.s) gas_accretion = config.get("config_gas_accretion", None) - scale(gas_accretion, "kappa_env", u.cm ** 2 / u.g) + scale(gas_accretion, "kappa_env", u.cm**2 / u.g) defaults = config.get("defaults", None) scale(defaults, "DEF_R_IN", u.au) diff --git a/chemcomp/helpers/main_helpers.py b/chemcomp/helpers/main_helpers.py index bf408b1..21ca67a 100755 --- a/chemcomp/helpers/main_helpers.py +++ b/chemcomp/helpers/main_helpers.py @@ -4,6 +4,7 @@ import numpy as np from slugify import slugify + try: from yaml import CLoader as Loader, CDumper as Dumper except ImportError: @@ -39,7 +40,7 @@ def file_checks_out(config_file, continue_job): if not continue_job: return True - output = import_config(config_file).get("output",{}) + output = import_config(config_file).get("output", {}) name = output.get("name", "planet") filename = output.get("directory", os.path.join(OUTPUT_PATH, "{}.h5".format(name))) @@ -94,11 +95,23 @@ def create_pipeline(job, config): values = np.array([mesh.flatten() for mesh in parameter_mesh]) if double_disks: - files, names = double_disks_wrapper(template, values, parameter_info, name=name, - save_disk=save_disk, save_interval=save_interval) + files, names = double_disks_wrapper( + template, + values, + parameter_info, + name=name, + save_disk=save_disk, + save_interval=save_interval, + ) else: - files = create_yaml_test_range(template, values, parameter_info, name=name, save_disk=save_disk, - save_interval=save_interval) + files = create_yaml_test_range( + template, + values, + parameter_info, + name=name, + save_disk=save_disk, + save_interval=save_interval, + ) names = [name] return files, names @@ -107,7 +120,7 @@ def create_pipeline(job, config): def run_setup(config_file): """load config and start modelling""" if not os.path.isdir(OUTPUT_PATH): - print(f'creating Folder {OUTPUT_PATH}') + print(f"creating Folder {OUTPUT_PATH}") os.mkdir(OUTPUT_PATH) # Load config @@ -131,7 +144,9 @@ def run_setup(config_file): planet_model = "BertPlanet" print("No planetmodel specified, using BertPlanet") - disk = getattr(chemcomp.disks, disk_model)(defaults, chemistry=chemistry, **config_disk) + disk = getattr(chemcomp.disks, disk_model)( + defaults, chemistry=chemistry, **config_disk + ) accretion_models = get_acc_models( config_pebble_accretion, @@ -152,12 +167,16 @@ def run_setup(config_file): def zip_all(names=[], delete=False): if not os.path.isdir(ZIP_OUTPUT_PATH): - print(f'creating Folder {ZIP_OUTPUT_PATH}') + print(f"creating Folder {ZIP_OUTPUT_PATH}") os.mkdir(ZIP_OUTPUT_PATH) zipf = zipfile.ZipFile( - os.path.join(ZIP_OUTPUT_PATH, - "{}_{}.zip".format(os.path.commonprefix(names), datetime.now().strftime("%Y%m%d-%H%M%S"))), + os.path.join( + ZIP_OUTPUT_PATH, + "{}_{}.zip".format( + os.path.commonprefix(names), datetime.now().strftime("%Y%m%d-%H%M%S") + ), + ), "w", zipfile.ZIP_DEFLATED, ) @@ -198,7 +217,15 @@ def is_defined(conf, acc): ] -def write_yaml(template, parameter_info, values, template_name, name="", save_disk=False, save_interval=None): +def write_yaml( + template, + parameter_info, + values, + template_name, + name="", + save_disk=False, + save_interval=None, +): if isinstance(template, str): with open(os.path.join(CONFIG_PATH, "{}.yaml".format(template))) as f: config_template = load(f, Loader=Loader) @@ -235,7 +262,9 @@ def write_yaml(template, parameter_info, values, template_name, name="", save_di return new_file_name, slug, plot_name -def create_yaml_test_range(template, values, parameter_info, name="", save_disk=False, save_interval=None): +def create_yaml_test_range( + template, values, parameter_info, name="", save_disk=False, save_interval=None +): files, slugs, plot_names = [], [], [] if not isinstance(template, str): @@ -244,7 +273,7 @@ def create_yaml_test_range(template, values, parameter_info, name="", save_disk= template_name = template if not os.path.isdir(CONFIG_PATH): - print(f'creating Folder {CONFIG_PATH}') + print(f"creating Folder {CONFIG_PATH}") os.mkdir(CONFIG_PATH) path = os.path.join(CONFIG_PATH, "p_{}".format(name)) @@ -255,9 +284,15 @@ def create_yaml_test_range(template, values, parameter_info, name="", save_disk= os.mkdir(path) for i in range(len(values.T)): - file, slug, plot_name = write_yaml(template, parameter_info, values.T[i], - template_name=template_name, name=name, - save_disk=save_disk, save_interval=save_interval) + file, slug, plot_name = write_yaml( + template, + parameter_info, + values.T[i], + template_name=template_name, + name=name, + save_disk=save_disk, + save_interval=save_interval, + ) files.append(file) slugs.append(slug) plot_names.append(plot_name) @@ -266,14 +301,16 @@ def create_yaml_test_range(template, values, parameter_info, name="", save_disk= slugfile.write("\n".join(slugs)) with open( - os.path.join(OUTPUT_PATH, "plot_names_{}.dat".format(name)), "w" + os.path.join(OUTPUT_PATH, "plot_names_{}.dat".format(name)), "w" ) as plot_name_file: plot_name_file.write("\n".join(plot_names)) return files -def double_disks_wrapper(template, values, parameter_info, name, save_disk, save_interval): +def double_disks_wrapper( + template, values, parameter_info, name, save_disk, save_interval +): """ change template and create four different disks: With Pla, With evap, With Pla+evap, plain Note: currently only two disks. Thank you referee! @@ -288,14 +325,22 @@ def double_disks_wrapper(template, values, parameter_info, name, save_disk, save template_None = config_template.copy() template_None["config_disk"]["DTG_pla"] = 0.0 template_None["config_disk"]["evaporation"] = False - files.append(create_yaml_test_range(template_None, values, parameter_info, name, save_disk, save_interval)) + files.append( + create_yaml_test_range( + template_None, values, parameter_info, name, save_disk, save_interval + ) + ) # evap: template_evap = config_template.copy() template_evap["config_disk"]["DTG_pla"] = 0.0 template_evap["config_disk"]["evaporation"] = True name_evap = f"{name}_evap" - files.append(create_yaml_test_range(template_evap, values, parameter_info, name_evap, save_disk, save_interval)) + files.append( + create_yaml_test_range( + template_evap, values, parameter_info, name_evap, save_disk, save_interval + ) + ) return [item for sublist in files for item in sublist], [name_evap, name] diff --git a/chemcomp/helpers/solvediffonedee.py b/chemcomp/helpers/solvediffonedee.py index 66b1e8d..2bdf9b3 100755 --- a/chemcomp/helpers/solvediffonedee.py +++ b/chemcomp/helpers/solvediffonedee.py @@ -1,8 +1,21 @@ from .tridag import * -def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False, - extracond=None, retall=False): +def solvediffonedee_components( + x, + y, + v, + d, + g, + s, + bcl, + bcr, + dt=None, + int=False, + upwind=False, + extracond=None, + retall=False, +): """ Vectorised version of solvediffonedee (works with sigma_g_components_mol) @@ -121,7 +134,7 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u # # Apply it to y # - a[0] = 0. + a[0] = 0.0 b[0] = bcl[1] - bcl[0] / dxi[0] c[0] = bcl[0] / dxi[0] r[0] = bcl[2] @@ -129,7 +142,7 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u # # Apply it to y/g # - a[0] = 0. + a[0] = 0.0 b[0] = (bcl[1] - bcl[0] / dxi[0]) / g[0] c[0] = (bcl[0] / dxi[0]) / g[1] r[0] = bcl[2] @@ -142,7 +155,7 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u # a[-1] = -bcr[0] / dxi[-1] b[-1] = bcr[1] + bcr[0] / dxi[-1] - c[-1] = 0. + c[-1] = 0.0 r[-1] = bcr[2] else: # @@ -150,7 +163,7 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u # a[-1] = (-bcr[0] / dxi[-1]) / g[-2] b[-1] = (bcr[1] + bcr[0] / dxi[-1]) / g[-1] - c[-1] = 0. + c[-1] = 0.0 r[-1] = bcr[2] # # Internal conditions, if present @@ -196,7 +209,21 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u return sol -def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False, extracond=None, retall=False): +def solvediffonedee( + x, + y, + v, + d, + g, + s, + bcl, + bcr, + dt=None, + int=False, + upwind=False, + extracond=None, + retall=False, +): """ Solve the 1-D linear advection-diffusion equation with boundary conditions, either in stationary state, or as an implicit time step. @@ -312,7 +339,7 @@ def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False # # Apply it to y # - a[0] = 0. + a[0] = 0.0 b[0] = bcl[1] - bcl[0] / dxi[0] c[0] = bcl[0] / dxi[0] r[0] = bcl[2] @@ -320,7 +347,7 @@ def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False # # Apply it to y/g # - a[0] = 0. + a[0] = 0.0 b[0] = (bcl[1] - bcl[0] / dxi[0]) / g[0] c[0] = (bcl[0] / dxi[0]) / g[1] r[0] = bcl[2] @@ -333,7 +360,7 @@ def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False # a[-1] = -bcr[0] / dxi[-1] b[-1] = bcr[1] + bcr[0] / dxi[-1] - c[-1] = 0. + c[-1] = 0.0 r[-1] = bcr[2] else: # @@ -341,7 +368,7 @@ def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False # a[-1] = (-bcr[0] / dxi[-1]) / g[-2] b[-1] = (bcr[1] + bcr[0] / dxi[-1]) / g[-1] - c[-1] = 0. + c[-1] = 0.0 r[-1] = bcr[2] # # Internal conditions, if present @@ -439,7 +466,9 @@ def getfluxonedee(x, y, v, d, g, int=False, upwind=False, oned=False): # Now compute the flux # if not oned: - flux = b[:, np.newaxis, np.newaxis] * y[:-1] + c[:, np.newaxis, np.newaxis] * y[1:] + flux = ( + b[:, np.newaxis, np.newaxis] * y[:-1] + c[:, np.newaxis, np.newaxis] * y[1:] + ) else: flux = b * y[:-1] + c * y[1:] # diff --git a/chemcomp/helpers/tridag.py b/chemcomp/helpers/tridag.py index 24854fe..caef300 100755 --- a/chemcomp/helpers/tridag.py +++ b/chemcomp/helpers/tridag.py @@ -1,10 +1,11 @@ import numpy as np from scipy.linalg import solve_banded + def tridag_components(a, b, c, r): """ - The Numerical Recipes tridag() routine for Python, using the SciPy library. - This is just a wrapper around the solve_banded function of SciPy, to make it + The Numerical Recipes tridag() routine for Python, using the SciPy library. + This is just a wrapper around the solve_banded function of SciPy, to make it work like the tridag routine. """ n = r.shape @@ -13,8 +14,7 @@ def tridag_components(a, b, c, r): ab[1, :] = b[:] ab[2, :-1] = a[1:] - u = np.transpose( - [solve_banded((1, 1), ab[:, :, i], r[:, i]) for i in range(n[1])]) + u = np.transpose([solve_banded((1, 1), ab[:, :, i], r[:, i]) for i in range(n[1])]) return u @@ -43,8 +43,7 @@ def pentdag(aa, a, b, c, cc, r): ab[0, 2:] = cc[:-2] ab[1, 1:] = c[:-1] ab[2, :] = b[:] - ab[3,:-1] = a[1:] - ab[4,:-2] = aa[2:] - u = solve_banded((2,2),ab,r) + ab[3, :-1] = a[1:] + ab[4, :-2] = aa[2:] + u = solve_banded((2, 2), ab, r) return u - diff --git a/chemcomp/helpers/units/chemistry_const.py b/chemcomp/helpers/units/chemistry_const.py index 62f415e..0653bfb 100755 --- a/chemcomp/helpers/units/chemistry_const.py +++ b/chemcomp/helpers/units/chemistry_const.py @@ -66,9 +66,9 @@ MassC = 12.0 # C in terms of H # Masses in terms of H (atomic unit) -MassCO = MassC + MassO #28.0 # CO mass in terms of H +MassCO = MassC + MassO # 28.0 # CO mass in terms of H MassCH4 = MassC + 4 * MassH # CH4 in terms of H -MassCO2 = MassC + 2 * MassO # CO2 in terms of H +MassCO2 = MassC + 2 * MassO # CO2 in terms of H MassH2O = MassO + 2 * MassH # H20 in terms of H MassFe2O3 = 2 * MassFe + 3 * MassO @@ -78,7 +78,7 @@ MassMg2SiO4 = 2 * MassMg + MassSi + 4 * MassO # MgSiO3 in terms of H MassNH3 = 3 * MassH + MassN MassN2 = 2 * MassN -MassH2S = 2 * MassH + MassS +MassH2S = 2 * MassH + MassS MassNaAlSi3O8 = MassNa + MassAl + 3 * MassSi + 8 * MassO MassKAlSi3O8 = MassK + MassAl + 3 * MassSi + 8 * MassO MassTiO = MassTi + MassO diff --git a/chemcomp/planets/_planet_class.py b/chemcomp/planets/_planet_class.py index 1449cf5..236a6f4 100755 --- a/chemcomp/planets/_planet_class.py +++ b/chemcomp/planets/_planet_class.py @@ -25,43 +25,115 @@ def __init__(self, planet, output): self.conf_save_disk = output.get("save_disk", False) self.conf_save_disk_interval = output.get("save_disk_interval", 100) self._save_counter = 0 - filename = output.get("directory", os.path.join(OUTPUT_PATH, "{}.h5".format(name))) + filename = output.get( + "directory", os.path.join(OUTPUT_PATH, "{}.h5".format(name)) + ) if os.path.exists(filename): os.remove(filename) - self._h5file = tables.open_file(filename, mode='w', title="Output of {}".format(name)) + self._h5file = tables.open_file( + filename, mode="w", title="Output of {}".format(name) + ) self._h5file.set_node_attr("/", "finished", 0) - self.output_disk = self._h5file.create_group("/", 'disk', 'Output of the Disk') - self.output_planet = self._h5file.create_group("/", 'planet', 'Output of the Planet') - self.output_acc = self._h5file.create_group("/", 'acc', 'Output of accretion') + self.output_disk = self._h5file.create_group("/", "disk", "Output of the Disk") + self.output_planet = self._h5file.create_group( + "/", "planet", "Output of the Planet" + ) + self.output_acc = self._h5file.create_group("/", "acc", "Output of accretion") self.planet = planet self.disk = planet.disk - list_of_disk_quantities = output.get("quantities", - ["a_1", "t", "sigma_g", "sigma_dust", "sigma_dust_components", "T", - "sigma_g_components", - "f_m", "stokes_number_pebbles", "stokes_number_df", "stokes_number_drift", - "stokes_number_frag", "stokes_number_small", "pebble_flux", - "cum_pebble_flux", - "peb_iso", "peb_iso_pi_crit", "T_visc", "T_irr", "mu", "m_dot", "vr_dust", "vr_gas", "m_dot_components"]) - - self.dict_of_planet_units = {"t": "u.s", "M": "u.g", "M_a": "u.g", "M_c": "u.g", "a_p": "u.cm", - "T": "u.K", "comp_a": "u.g", "comp_c": "u.g", "tau_m": "u.s", - "sigma_g": "u.g/u.cm**2", - "gamma_tot": "u.g*u.cm**2/u.s**2", "sigma_peb": "u.g/u.cm**2", - "pebble_flux": "u.g/u.s", "peb_iso": "u.g"} - self.dict_of_acc_units = {"m_dot": "u.g/u.s", "m_dot_a": "u.g/u.s", "m_dot_c": "u.g/u.s", - "m_dot_a_chem": "u.g/u.s", "m_dot_c_chem": "u.g/u.s", "M_z": "u.g"} - - self.list_of_planet_quantities = ["t", "M", "M_a", "M_c", "a_p", "T", "comp_a", "comp_c", "tau_m", - "gamma_tot", "sigma_g", "gamma_norm", "fSigma", "sigma_peb", "pebble_flux", - "peb_iso"] - self.list_of_acc_quantities = ["m_dot", "m_dot_a", "m_dot_c", "m_dot_a_chem", "m_dot_c_chem", "M_z", "regime"] + list_of_disk_quantities = output.get( + "quantities", + [ + "a_1", + "t", + "sigma_g", + "sigma_dust", + "sigma_dust_components", + "T", + "sigma_g_components", + "f_m", + "stokes_number_pebbles", + "stokes_number_df", + "stokes_number_drift", + "stokes_number_frag", + "stokes_number_small", + "pebble_flux", + "cum_pebble_flux", + "peb_iso", + "peb_iso_pi_crit", + "T_visc", + "T_irr", + "mu", + "m_dot", + "vr_dust", + "vr_gas", + "m_dot_components", + ], + ) + + self.dict_of_planet_units = { + "t": "u.s", + "M": "u.g", + "M_a": "u.g", + "M_c": "u.g", + "a_p": "u.cm", + "T": "u.K", + "comp_a": "u.g", + "comp_c": "u.g", + "tau_m": "u.s", + "sigma_g": "u.g/u.cm**2", + "gamma_tot": "u.g*u.cm**2/u.s**2", + "sigma_peb": "u.g/u.cm**2", + "pebble_flux": "u.g/u.s", + "peb_iso": "u.g", + } + self.dict_of_acc_units = { + "m_dot": "u.g/u.s", + "m_dot_a": "u.g/u.s", + "m_dot_c": "u.g/u.s", + "m_dot_a_chem": "u.g/u.s", + "m_dot_c_chem": "u.g/u.s", + "M_z": "u.g", + } + + self.list_of_planet_quantities = [ + "t", + "M", + "M_a", + "M_c", + "a_p", + "T", + "comp_a", + "comp_c", + "tau_m", + "gamma_tot", + "sigma_g", + "gamma_norm", + "fSigma", + "sigma_peb", + "pebble_flux", + "peb_iso", + ] + self.list_of_acc_quantities = [ + "m_dot", + "m_dot_a", + "m_dot_c", + "m_dot_a_chem", + "m_dot_c_chem", + "M_z", + "regime", + ] # trim list of quantities to existing: _, self.list_of_disk_quantities = [], [] - [self.list_of_disk_quantities.append(quant) if hasattr(self.disk, f"{quant}") else _.append(quant) for quant in - list_of_disk_quantities] + [ + self.list_of_disk_quantities.append(quant) + if hasattr(self.disk, f"{quant}") + else _.append(quant) + for quant in list_of_disk_quantities + ] self._h5file.create_earray(self.output_disk, "r", obj=self.disk.r) self._h5file.create_earray(self.output_disk, "r_i", obj=self.disk.r_i) @@ -69,12 +141,16 @@ def __init__(self, planet, output): self.list_of_disk_pointer = [] for quantity in self.list_of_disk_quantities: obj = np.array(getattr(self.disk, quantity))[np.newaxis] - self.list_of_disk_pointer.append(self._h5file.create_earray(self.output_disk, quantity, obj=obj)) + self.list_of_disk_pointer.append( + self._h5file.create_earray(self.output_disk, quantity, obj=obj) + ) self.list_of_planet_pointer = [] for quantity in self.list_of_planet_quantities: obj = np.array(getattr(self.planet, quantity))[np.newaxis] - self.list_of_planet_pointer.append(self._h5file.create_earray(self.output_planet, quantity, obj=obj)) + self.list_of_planet_pointer.append( + self._h5file.create_earray(self.output_planet, quantity, obj=obj) + ) acc_name_dict = {"PebbleAccretion": "peb", "GasAccretion": "gas"} # drop acc model if not used: @@ -90,12 +166,18 @@ def __init__(self, planet, output): for i in self.acc_index_list: self.list_of_acc_pointer[i] = [] for quantity in self.list_of_acc_quantities: - obj = np.array(getattr(self.planet.accretion_models[i], quantity))[np.newaxis] + obj = np.array(getattr(self.planet.accretion_models[i], quantity))[ + np.newaxis + ] self.list_of_acc_pointer[i].append( - self._h5file.create_earray(self.output_planet, "{}_{}".format(quantity, self.acc_models[i]), - obj=obj)) - - self._planet_units = self._h5file.create_group("/planet", 'units', 'units') + self._h5file.create_earray( + self.output_planet, + "{}_{}".format(quantity, self.acc_models[i]), + obj=obj, + ) + ) + + self._planet_units = self._h5file.create_group("/planet", "units", "units") for k, v in self.dict_of_planet_units.items(): self._h5file.set_node_attr("/planet/units", k, v) @@ -119,13 +201,19 @@ def save_acc(self): for j in self.acc_index_list: ptr = self.list_of_acc_pointer[j] for i in range(len(ptr)): - array = getattr(self.planet.accretion_models[j], self.list_of_acc_quantities[i]) + array = getattr( + self.planet.accretion_models[j], self.list_of_acc_quantities[i] + ) array = np.array(array) self.list_of_acc_pointer[j][i].append(array[np.newaxis]) def save(self, force=False): if self.planet.t >= self.planet.save_interval * self._save_counter or force: - if self.conf_save_disk and self._save_counter % self.conf_save_disk_interval == 0 or force: + if ( + self.conf_save_disk + and self._save_counter % self.conf_save_disk_interval == 0 + or force + ): self.save_disk() self.save_planet() self.save_acc() @@ -138,10 +226,19 @@ def finish(self, error=None): class Planet(object): - """ docstring for Planet. """ + """docstring for Planet.""" def __init__( - self, output, rho_c, a_p, t_0, disk, accretion_models=[], M_c=None, M_a=None, **kwargs + self, + output, + rho_c, + a_p, + t_0, + disk, + accretion_models=[], + M_c=None, + M_a=None, + **kwargs, ): self.save_interval = output.get("save_interval", 50000 * u.yr).cgs.value self._print_number = 0 @@ -161,53 +258,77 @@ def __init__( if M_c is None: # use Lambrecht2012 transitionmass # -> ignored if pebiso_start - M_c = (1-self.solid_frac)* ( - np.sqrt(1.0 / 3.0) - * (disk.eta * disk.v_k) ** 3.0 - / (const.G.cgs.value * disk.omega_k) + M_c = (1 - self.solid_frac) * ( + np.sqrt(1.0 / 3.0) + * (disk.eta * disk.v_k) ** 3.0 + / (const.G.cgs.value * disk.omega_k) ) M_c = M0_fact * np.interp(a_p.cgs.value, disk.r, M_c) * u.g if M_a is None: M_a = 0.0 * u.g if self.solid_frac != 0.0: - M_a = self.solid_frac / (1.0 - self.solid_frac) * M_c + M_a = self.solid_frac / (1.0 - self.solid_frac) * M_c self.M_c = 1.0 * M_c.cgs.value # core mass self.M_a = 1.0 * M_a.cgs.value # atmospheric mass self.M = self.M_a + self.M_c - self.accretion_models = accretion_models # list containing the accretion objects + self.accretion_models = ( + accretion_models # list containing the accretion objects + ) self._accretion_models_dict = {acc.name: acc for acc in self.accretion_models} self._accrete = eval_kwargs(kwargs.get("accrete", True)) self.matter_removal = eval_kwargs(kwargs.get("matter_removal", True)) self.force_insitu = eval_kwargs(kwargs.get("force_insitu", False)) - self.r_in = eval_kwargs(kwargs.get("r_in", self.disk.r_i[0]*u.cm).cgs.value) # where to stop the planet formation + self.r_in = eval_kwargs( + kwargs.get("r_in", self.disk.r_i[0] * u.cm).cgs.value + ) # where to stop the planet formation self.name = output.get("name", "planet") - self.output_file = output.get("directory", os.path.join(OUTPUT_PATH, "{}.h5").format(self.name)) + self.output_file = output.get( + "directory", os.path.join(OUTPUT_PATH, "{}.h5").format(self.name) + ) self.tau_m = np.inf # migration rate - self.fSigma = 1.0 # gap depth + self.fSigma = 1.0 # gap depth - self.M_end = eval_kwargs(kwargs.get("M_end", 1.0 * u.M_sun).cgs.value) # final mass at which we interrupt the formation of the planet + self.M_end = eval_kwargs( + kwargs.get("M_end", 1.0 * u.M_sun).cgs.value + ) # final mass at which we interrupt the formation of the planet - self.past_pebble_iso = False # boolean flag that stores if we passed the pebiso mass - self._keep_peb_iso = eval_kwargs(kwargs.get("keep_peb_iso", True)) # keep pebble isolation mass constant when pebble isolation mass is reached (otherwise we might have multiple phases of pebble accretion) - self._use_pebiso_diffusion = eval_kwargs(kwargs.get("use_pebiso_diffusion", True)) # use the diffusion part of the pebble isolation mass + self.past_pebble_iso = ( + False # boolean flag that stores if we passed the pebiso mass + ) + self._keep_peb_iso = eval_kwargs( + kwargs.get("keep_peb_iso", True) + ) # keep pebble isolation mass constant when pebble isolation mass is reached (otherwise we might have multiple phases of pebble accretion) + self._use_pebiso_diffusion = eval_kwargs( + kwargs.get("use_pebiso_diffusion", True) + ) # use the diffusion part of the pebble isolation mass self._init_accretion() # initialise the accretion models self.update(second=True) # first update of all quantities - self.chem_mask_array = self.disk._chem.mask_array != 0 # see mask_array in chemistry class (needed since there are less elements than mol. species) - self.comp_a = self.chemistry_solid * self.M_a # molecular and elemental composition of the planetary envelope - self.comp_c = self.chemistry_solid * self.M_c # molecular and elemental composition of the planetary core - - self.output = DataObject(self, output) # Initialize the IO class - self.evolve_disk_init() # start evolving the disk to the point where the planet is placed into the disk + self.chem_mask_array = ( + self.disk._chem.mask_array != 0 + ) # see mask_array in chemistry class (needed since there are less elements than mol. species) + self.comp_a = ( + self.chemistry_solid * self.M_a + ) # molecular and elemental composition of the planetary envelope + self.comp_c = ( + self.chemistry_solid * self.M_c + ) # molecular and elemental composition of the planetary core + + self.output = DataObject(self, output) # Initialize the IO class + self.evolve_disk_init() # start evolving the disk to the point where the planet is placed into the disk if pebiso_start: # if planet is seeded with M=M_iso self.calc_pebble_iso() - print("{}: starting at peb_iso with {:.2e} M_e".format(self.name,(self.peb_iso*u.g).to(u.earthMass))) - self.M_c = (1.0-self.solid_frac) * self.peb_iso + print( + "{}: starting at peb_iso with {:.2e} M_e".format( + self.name, (self.peb_iso * u.g).to(u.earthMass) + ) + ) + self.M_c = (1.0 - self.solid_frac) * self.peb_iso self.M_a = self.solid_frac * self.peb_iso self.M = self.M_a + self.M_c self.past_pebble_iso = True @@ -216,18 +337,17 @@ def __init__( self.comp_a = self.chemistry_solid * self.M_a self.comp_c = self.chemistry_solid * self.M_c - def _init_accretion(self): - """ Initialise the Accretionmodels. """ + """Initialise the Accretionmodels.""" for acc in self.accretion_models: acc.init_accretion(self) def update_torque(self): - """ Calculate the Torques""" + """Calculate the Torques""" pass def update_a_p(self): - """ Migrate the Planet. """ + """Migrate the Planet.""" pass def update(self, second, interpolation_idx=None, interpolation_idx_interface=None): @@ -239,32 +359,57 @@ def update(self, second, interpolation_idx=None, interpolation_idx_interface=Non interpolation_idx can be used if gap should be skipped during the calculation of sigma_gap relevant quantities """ - if interpolation_idx is None: # indexes as which we want to interpolate the disk quantities to the planetary position - interpolation_idx = np.ones_like(self.disk.r,dtype=bool) + if ( + interpolation_idx is None + ): # indexes as which we want to interpolate the disk quantities to the planetary position + interpolation_idx = np.ones_like(self.disk.r, dtype=bool) if interpolation_idx_interface is None: - interpolation_idx_interface = np.ones_like(self.disk.r_i,dtype=bool) - - self.sigma_g = ( # surface density of the gas phase - interp1d(self.disk.r[interpolation_idx], self.disk.sigma_g[interpolation_idx], fill_value="extrapolate", assume_sorted=True)(self.a_p) - ) - self.chemistry_gas = interp1d(self.disk.r, self.disk.chemistry_gas, axis=0, fill_value="extrapolate", - assume_sorted=True)(self.a_p) # composition of the gas phase + interpolation_idx_interface = np.ones_like(self.disk.r_i, dtype=bool) + + self.sigma_g = interp1d( # surface density of the gas phase + self.disk.r[interpolation_idx], + self.disk.sigma_g[interpolation_idx], + fill_value="extrapolate", + assume_sorted=True, + )(self.a_p) + self.chemistry_gas = interp1d( + self.disk.r, + self.disk.chemistry_gas, + axis=0, + fill_value="extrapolate", + assume_sorted=True, + )( + self.a_p + ) # composition of the gas phase if self.t < self.disk.begin_photevap: assert np.all(self.chemistry_gas >= 0.0), "gas chemistry is negative!" elif not np.all(self.chemistry_gas >= 0.0): - self.chemistry_gas = np.maximum(self.chemistry_gas, np.zeros_like(self.chemistry_gas)) # composition of the gas phase + self.chemistry_gas = np.maximum( + self.chemistry_gas, np.zeros_like(self.chemistry_gas) + ) # composition of the gas phase if self.t < self.disk.begin_photevap: assert self.sigma_g >= 0.0, "negative surface density!" else: - self.sigma_g = np.maximum(self.sigma_g, 1e-60) # surface density + self.sigma_g = np.maximum(self.sigma_g, 1e-60) # surface density - self.sigma_g_components = self.sigma_g * self.chemistry_gas # surface density as a compositional vector + self.sigma_g_components = ( + self.sigma_g * self.chemistry_gas + ) # surface density as a compositional vector - self.nu = interp1d(self.disk.r, self.disk.nu, fill_value="extrapolate", assume_sorted=True)(self.a_p) # viscosity + self.nu = interp1d( + self.disk.r, self.disk.nu, fill_value="extrapolate", assume_sorted=True + )( + self.a_p + ) # viscosity self.m_dot_disk = ( # viscous accretion rate of the gas phase of the disk - interp1d(self.disk.r_i[interpolation_idx_interface], np.abs(self.disk.m_dot[interpolation_idx_interface]),fill_value="extrapolate", assume_sorted=True)(self.a_p) + interp1d( + self.disk.r_i[interpolation_idx_interface], + np.abs(self.disk.m_dot[interpolation_idx_interface]), + fill_value="extrapolate", + assume_sorted=True, + )(self.a_p) ) if self.m_dot_disk <= 0.0: @@ -274,55 +419,129 @@ def update(self, second, interpolation_idx=None, interpolation_idx_interface=Non self.m_dot_disk_components = self.chemistry_gas * self.m_dot_disk if self.t < self.disk.begin_photevap: - assert self.m_dot_disk >= 0.0, "the absolute value of m_dot_disk is not positive! Interpolation error" + assert ( + self.m_dot_disk >= 0.0 + ), "the absolute value of m_dot_disk is not positive! Interpolation error" else: - self.m_dot_disk_components = np.maximum(self.m_dot_disk_components, np.ones_like(self.chemistry_gas)*1e-60) + self.m_dot_disk_components = np.maximum( + self.m_dot_disk_components, np.ones_like(self.chemistry_gas) * 1e-60 + ) self.m_dot_disk = np.maximum(1e-60, self.m_dot_disk) - self.sigma_g_pert = ( # perturbed surface density of the gap - interp1d(self.disk.r, self.disk.sigma_g, fill_value="extrapolate", assume_sorted=True)(self.a_p) - ) - self.sigma_g_components_pert = interp1d(self.disk.r, # perturbed surface density of the gap as a compositional vector - self.disk.sigma_g_components, axis=0, - fill_value="extrapolate", assume_sorted=True)( - self.a_p) - self.T = interp1d(self.disk.r, self.disk.T, fill_value="extrapolate", assume_sorted=True)(self.a_p) # Temperature - self.H_g = interp1d(self.disk.r, self.disk.aspect_ratio, fill_value="extrapolate", assume_sorted=True)(self.a_p) * self.a_p # Disk scalehight - self.c_s = interp1d(self.disk.r, self.disk.c_s, fill_value="extrapolate", assume_sorted=True)(self.a_p) # soundspeed - self.rho_g = self.sigma_g / (np.sqrt(2.0 * np.pi) * self.H_g) # gas volumedensity - self.grad_sig = interp1d(self.disk.r[interpolation_idx], self.disk.grad_sig[interpolation_idx], fill_value="extrapolate", assume_sorted=True)(self.a_p) # gradient of the surfacedensity - self.grad_T = interp1d(self.disk.r, self.disk.grad_T, fill_value="extrapolate", assume_sorted=True)(self.a_p) # gradient of the Temperature - self.grad_p = interp1d(self.disk.r, self.disk.grad_p, fill_value="extrapolate", assume_sorted=True)(self.a_p) # gradient of the pressure - self.r_h = self.a_p * (self.M / (3.0 * self.disk.M_s)) ** (1.0 / 3.0) # hill radius + self.sigma_g_pert = interp1d( # perturbed surface density of the gap + self.disk.r, self.disk.sigma_g, fill_value="extrapolate", assume_sorted=True + )(self.a_p) + self.sigma_g_components_pert = interp1d( + self.disk.r, # perturbed surface density of the gap as a compositional vector + self.disk.sigma_g_components, + axis=0, + fill_value="extrapolate", + assume_sorted=True, + )(self.a_p) + self.T = interp1d( + self.disk.r, self.disk.T, fill_value="extrapolate", assume_sorted=True + )( + self.a_p + ) # Temperature + self.H_g = ( + interp1d( + self.disk.r, + self.disk.aspect_ratio, + fill_value="extrapolate", + assume_sorted=True, + )(self.a_p) + * self.a_p + ) # Disk scalehight + self.c_s = interp1d( + self.disk.r, self.disk.c_s, fill_value="extrapolate", assume_sorted=True + )( + self.a_p + ) # soundspeed + self.rho_g = self.sigma_g / ( + np.sqrt(2.0 * np.pi) * self.H_g + ) # gas volumedensity + self.grad_sig = interp1d( + self.disk.r[interpolation_idx], + self.disk.grad_sig[interpolation_idx], + fill_value="extrapolate", + assume_sorted=True, + )( + self.a_p + ) # gradient of the surfacedensity + self.grad_T = interp1d( + self.disk.r, self.disk.grad_T, fill_value="extrapolate", assume_sorted=True + )( + self.a_p + ) # gradient of the Temperature + self.grad_p = interp1d( + self.disk.r, self.disk.grad_p, fill_value="extrapolate", assume_sorted=True + )( + self.a_p + ) # gradient of the pressure + self.r_h = self.a_p * (self.M / (3.0 * self.disk.M_s)) ** ( + 1.0 / 3.0 + ) # hill radius self.v_k = np.sqrt(G * self.disk.M_s / self.a_p) # keppler velocity self.omega_k = self.v_k / self.a_p # keppler orbital velocity - self.del_v = interp1d(self.disk.r, self.disk.eta, fill_value="extrapolate", assume_sorted=True)( - self.a_p) * self.v_k # difference between keppler and real gas orbital speed - self.r_b = (G * self.M) / (self.del_v ** 2.0) # bondi radius + self.del_v = ( + interp1d( + self.disk.r, self.disk.eta, fill_value="extrapolate", assume_sorted=True + )(self.a_p) + * self.v_k + ) # difference between keppler and real gas orbital speed + self.r_b = (G * self.M) / (self.del_v**2.0) # bondi radius self.v_h = self.r_h * self.omega_k # hill velocity - self.P = interp1d(self.disk.r, self.disk.P, fill_value="extrapolate", assume_sorted=True)( - self.a_p) + self.P = interp1d( + self.disk.r, self.disk.P, fill_value="extrapolate", assume_sorted=True + )(self.a_p) if second: # This may only be called after migration has been applied # migration calculates the gap parameter according to Kanagawa/Ndugu -> self.fsigma - self.pebble_flux = interp1d(self.disk.r, self.disk.pebble_flux, fill_value="extrapolate", assume_sorted=True)(self.a_p) # pebble flux + self.pebble_flux = interp1d( + self.disk.r, + self.disk.pebble_flux, + fill_value="extrapolate", + assume_sorted=True, + )( + self.a_p + ) # pebble flux if hasattr(self.disk, "f_m"): - self.sigma_peb = interp1d(self.disk.r, self.disk.sigma_dust * self.disk.f_m, fill_value="extrapolate", assume_sorted=True)(self.a_p) # pebble surface density - self.sigma_peb_components = interp1d(self.disk.r, # pebble surface density as a compositional vector - self.disk.sigma_dust_components * self.disk.f_m[:, np.newaxis, - np.newaxis], axis=0, - fill_value="extrapolate", - assume_sorted=True)(self.a_p) + self.sigma_peb = interp1d( + self.disk.r, + self.disk.sigma_dust * self.disk.f_m, + fill_value="extrapolate", + assume_sorted=True, + )( + self.a_p + ) # pebble surface density + self.sigma_peb_components = interp1d( + self.disk.r, # pebble surface density as a compositional vector + self.disk.sigma_dust_components + * self.disk.f_m[:, np.newaxis, np.newaxis], + axis=0, + fill_value="extrapolate", + assume_sorted=True, + )(self.a_p) else: - self.sigma_peb = interp1d(self.disk.r, self.disk.sigma_dust, fill_value="extrapolate", assume_sorted=True)(self.a_p) # checked - self.sigma_peb_components = interp1d(self.disk.r, - self.disk.sigma_dust_components, axis=0, - fill_value="extrapolate", - assume_sorted=True)(self.a_p) + self.sigma_peb = interp1d( + self.disk.r, + self.disk.sigma_dust, + fill_value="extrapolate", + assume_sorted=True, + )( + self.a_p + ) # checked + self.sigma_peb_components = interp1d( + self.disk.r, + self.disk.sigma_dust_components, + axis=0, + fill_value="extrapolate", + assume_sorted=True, + )(self.a_p) self.calc_pebble_iso() @@ -333,18 +552,30 @@ def update(self, second, interpolation_idx=None, interpolation_idx_interface=Non # idx: cell center # tested -> Check - self.idx_i_r = np.searchsorted(self.disk.r_i, self.a_p, - side="right") # interface i+1: idx of right cell interface in self.disk.r_i - self.idx_i_l = max(self.idx_i_r - 1, 0) # interface i: idx of left cell interface in self.disk.r_i - self.idx = 1 * self.idx_i_l # center i: idx of cell (self.disk.r) planet is located in - - self.chemistry_solid = interp1d(self.disk.r, self.disk.chemistry_solid, axis=0, fill_value="extrapolate", - assume_sorted=True)(self.a_p) # compositional vector of pebbles + self.idx_i_r = np.searchsorted( + self.disk.r_i, self.a_p, side="right" + ) # interface i+1: idx of right cell interface in self.disk.r_i + self.idx_i_l = max( + self.idx_i_r - 1, 0 + ) # interface i: idx of left cell interface in self.disk.r_i + self.idx = ( + 1 * self.idx_i_l + ) # center i: idx of cell (self.disk.r) planet is located in + + self.chemistry_solid = interp1d( + self.disk.r, + self.disk.chemistry_solid, + axis=0, + fill_value="extrapolate", + assume_sorted=True, + )( + self.a_p + ) # compositional vector of pebbles self.update_parameters() def update_parameters(self): - """ Parameters that are dependent on planet Mass need to be kept updated """ + """Parameters that are dependent on planet Mass need to be kept updated""" pass @property @@ -352,32 +583,44 @@ def r_c(self): """ radius of the core """ - return (3.0 / 4.0 * self.M_c / (np.pi * self.rho_c)) ** (1.0 / 3.0) # core radius + return (3.0 / 4.0 * self.M_c / (np.pi * self.rho_c)) ** ( + 1.0 / 3.0 + ) # core radius def calc_pebble_iso(self): - """ calculate the pebble isolation mass.""" + """calculate the pebble isolation mass.""" if self._keep_peb_iso and self.past_pebble_iso: return self.disk.calc_pebble_iso(diffusion=self._use_pebiso_diffusion) - if hasattr(self.disk,"interpolation_idx"): - self.peb_iso = interp1d(self.disk.r[self.disk.interpolation_idx], self.disk.peb_iso[self.disk.interpolation_idx], fill_value="extrapolate", assume_sorted=True)(self.a_p) + if hasattr(self.disk, "interpolation_idx"): + self.peb_iso = interp1d( + self.disk.r[self.disk.interpolation_idx], + self.disk.peb_iso[self.disk.interpolation_idx], + fill_value="extrapolate", + assume_sorted=True, + )(self.a_p) else: - self.peb_iso = interp1d(self.disk.r, self.disk.peb_iso, fill_value="extrapolate", assume_sorted=True)(self.a_p) + self.peb_iso = interp1d( + self.disk.r, + self.disk.peb_iso, + fill_value="extrapolate", + assume_sorted=True, + )(self.a_p) self.past_pebble_iso = True if self.M > self.peb_iso else False return def update_all(self): - """ Wrapper that handles the complete update of disk and planet as well as migration.""" + """Wrapper that handles the complete update of disk and planet as well as migration.""" self.disk.update(self.t, self.h) self.update(second=False) self.update_a_p() self.update(second=True) # update all quantities effected by a_p def evolve_disk_init(self): - """ Evolve the disk before the protoplanetary seed is placed into the disk. """ + """Evolve the disk before the protoplanetary seed is placed into the disk.""" self.disk.evolve_init() # initialise disk and save @@ -394,7 +637,7 @@ def evolve_disk_init(self): self.print_params() def calc_timestep(self): - """ Fix the timestep and update time. Important: This is a fixed timestep currently!""" + """Fix the timestep and update time. Important: This is a fixed timestep currently!""" if self.disk._timestep_input: self.h = self.disk._timestep_input else: @@ -416,14 +659,14 @@ def keep_running(self): return True def grow_mass(self) -> None: - """ Main loop and supervisor of chemcomp. + """Main loop and supervisor of chemcomp. dt: the total time that the model should be evolved. Needs to be a astropy Quantity """ def increment_mass(acc): - """ increment the mass of the planet for a certain accretion model """ + """increment the mass of the planet for a certain accretion model""" acc.calc_m_dot() acc.update_z() self.total_m_dot += acc.m_dot @@ -435,7 +678,7 @@ def increment_mass(acc): acc.remove_mass_from_flux() def _grow(): - """ wrap around increment_mass and update final mass. """ + """wrap around increment_mass and update final mass.""" self.total_m_dot = 0.0 list(map(increment_mass, self.accretion_models)) self.M = self.M_a + self.M_c @@ -481,7 +724,7 @@ def _grow(): return def print_params(self, force=False): - """" Handle the IO using self.output -> DataObject. + """ " Handle the IO using self.output -> DataObject. Saves !after! accretion applied """ self.output.save(force) @@ -492,8 +735,14 @@ def print_params(self, force=False): CO = self.comp_a[0, 0] / self.comp_a[0, 1] * 16.0 / 12.0 except FloatingPointError: CO = 0.0 - print("{}: t={:.1f}Myr, dt={:.1f}yr, M={:.1f}M_e, r={:.1f}AU, CO={:.2f}".format(self.name, self.t / Myr, - self.h / year, - self.M / M_earth, - self.a_p / AU, CO)) + print( + "{}: t={:.1f}Myr, dt={:.1f}yr, M={:.1f}M_e, r={:.1f}AU, CO={:.2f}".format( + self.name, + self.t / Myr, + self.h / year, + self.M / M_earth, + self.a_p / AU, + CO, + ) + ) self._print_number += 1 diff --git a/chemcomp/planets/bertmigration.py b/chemcomp/planets/bertmigration.py index db868b8..d1b7f6f 100755 --- a/chemcomp/planets/bertmigration.py +++ b/chemcomp/planets/bertmigration.py @@ -13,24 +13,32 @@ class BertPlanet(Planet): """Migrating planet. torques adapted from Planettrack.F by Bertram Bitsch. Planet used in Schneider & Bitsch (2021a,b)""" def __init__( - self, - output, - disk, - accretion_models=[], - rho_c=5.5 * u.g / (u.cm ** 3), - a_p=5.2 * u.au, - t_0=5e5 * u.yr, - M_c=None, - M_a=None, - **kwargs + self, + output, + disk, + accretion_models=[], + rho_c=5.5 * u.g / (u.cm**3), + a_p=5.2 * u.au, + t_0=5e5 * u.yr, + M_c=None, + M_a=None, + **kwargs ): self.gamma_tot = 0.0 # total torque self.gamma_norm = 0.0 # torque normalized to gamma_0 - self._use_heat_torque = eval_kwargs(kwargs.get("use_heat_torque", True)) # use masset heating torque - self._use_dynamical_torque = eval_kwargs(kwargs.get("use_dynamical_torque", True)) # use paardekooper dynamical torque - self._apply_gap = eval_kwargs(kwargs.get("apply_gap_profile", True)) # if the gap should be added to the gas surface density - self._migrate = eval_kwargs(kwargs.get("migration", True)) # if the planet should migrate + self._use_heat_torque = eval_kwargs( + kwargs.get("use_heat_torque", True) + ) # use masset heating torque + self._use_dynamical_torque = eval_kwargs( + kwargs.get("use_dynamical_torque", True) + ) # use paardekooper dynamical torque + self._apply_gap = eval_kwargs( + kwargs.get("apply_gap_profile", True) + ) # if the gap should be added to the gas surface density + self._migrate = eval_kwargs( + kwargs.get("migration", True) + ) # if the planet should migrate super().__init__( a_p=a_p, @@ -46,11 +54,11 @@ def __init__( self._old_a_p = 1.0 * self.a_p self.init_paardekooper_units() - print(self.name, (self.M*u.g).to(u.earthMass).value) + print(self.name, (self.M * u.g).to(u.earthMass).value) def update_all(self): - """ Update all quantities including disk evolution, migration and planetary quantities. - Overwrites parent update_all() method """ + """Update all quantities including disk evolution, migration and planetary quantities. + Overwrites parent update_all() method""" if self._apply_gap: self.update(second=False) self.calc_gammaeff(migphase=False) @@ -59,25 +67,43 @@ def update_all(self): if self.fSigma < 1 and self.past_pebble_iso: self.disk.apply_gap_profile(self.a_p, self.fSigma, gap_width) - interpolation_idx = self.disk.interpolation_idx if hasattr(self.disk,"interpolation_idx") and self._apply_gap else None - interpolation_idx_interface = self.disk.interpolation_idx_interface if hasattr(self.disk, "interpolation_idx_interface") and self._apply_gap else None + interpolation_idx = ( + self.disk.interpolation_idx + if hasattr(self.disk, "interpolation_idx") and self._apply_gap + else None + ) + interpolation_idx_interface = ( + self.disk.interpolation_idx_interface + if hasattr(self.disk, "interpolation_idx_interface") and self._apply_gap + else None + ) self.disk.update(self.t, self.h) - self.update(second=False, interpolation_idx=interpolation_idx, interpolation_idx_interface=interpolation_idx_interface) + self.update( + second=False, + interpolation_idx=interpolation_idx, + interpolation_idx_interface=interpolation_idx_interface, + ) self.update_a_p() - self.update(second=True, interpolation_idx=interpolation_idx, interpolation_idx_interface=interpolation_idx_interface) # update all quantities effected by a_p + self.update( + second=True, + interpolation_idx=interpolation_idx, + interpolation_idx_interface=interpolation_idx_interface, + ) # update all quantities effected by a_p def init_paardekooper_units(self): - """ Define the paardekooper scaling factors.""" + """Define the paardekooper scaling factors.""" self._r_0 = 1.0 * self.a_p self._omega_0 = 1.0 * self.omega_k self._sigma_0 = 1.0 * self.sigma_g - self._q_d = np.pi * self._r_0 ** 2.0 * self._sigma_0 / self.disk.M_s + self._q_d = np.pi * self._r_0**2.0 * self._sigma_0 / self.disk.M_s self._nu_0 = 1.0 * self.nu # altough mig beta and alpha will change during runtime, we keep them constant since a gap would manipulate this even more! self.mig_beta = -self.grad_T # checked self.mig_alpha = -self.grad_sig # checked - self.xi = self.mig_beta - (self.disk._chem.gamma - 1.0) * self.mig_alpha # checked, see equation 5 + self.xi = ( + self.mig_beta - (self.disk._chem.gamma - 1.0) * self.mig_alpha + ) # checked, see equation 5 def calc_heating_torque(self): """ @@ -89,13 +115,24 @@ def calc_heating_torque(self): if self.M < thermal_crit: L_c_mod_M = 4.0 * np.pi * G * self.Vchi * self.rho_g / self._gammaeff - L_mod_M = G * (self.accretion_models[0].m_dot + self.accretion_models[1].m_dot) / self.r_c + L_mod_M = ( + G + * (self.accretion_models[0].m_dot + self.accretion_models[1].m_dot) + / self.r_c + ) L_norm = L_mod_M / L_c_mod_M lambda_c = np.sqrt(self.Vchi / (1.5 * self.omega_k * self._gammaeff)) eta = self.mig_alpha / 3.0 + (self.mig_beta + 3.0) / 6.0 - gamma_thermal = 1.61 * (self._gammaeff - 1.0) / self._gammaeff * eta * (self.H_g / lambda_c) * (L_norm - 1.0) + gamma_thermal = ( + 1.61 + * (self._gammaeff - 1.0) + / self._gammaeff + * eta + * (self.H_g / lambda_c) + * (L_norm - 1.0) + ) self.gamma_tot += self.gamma_zero * gamma_thermal self.gamma_norm += gamma_thermal @@ -123,9 +160,15 @@ def update_a_p(self): self.tau_m = self.get_static_tau_m() if self.fSigma <= 0.53: - tau_II = - self.a_p ** 2.0 / self.nu * np.maximum(1.0, self.M/(4*np.pi*self.sigma_g*self.a_p**2.0)) + tau_II = ( + -self.a_p**2.0 + / self.nu + * np.maximum(1.0, self.M / (4 * np.pi * self.sigma_g * self.a_p**2.0)) + ) if self.fSigma > 0.1: - self.tau_m = (self.tau_m - tau_II) / (0.53 - 0.1) * (self.fSigma-0.1) + tau_II + self.tau_m = (self.tau_m - tau_II) / (0.53 - 0.1) * ( + self.fSigma - 0.1 + ) + tau_II else: self.tau_m = tau_II @@ -145,12 +188,23 @@ def get_dynamical_tau_m(self): """ - k = 8.0 / 3.0 / np.pi * (3.0 / 2.0 - self.mig_alpha) * self.gamma_norm * self._q_d ** 2.0 * self._xs ** 3.0 / ( - self.aspect_ratio ** 2.0) * self._r_0 ** 2.0 * self._omega_0 / self._nu_0 * ( - self.a_p / self._r_0) ** (5.0 - 3.0 * self.mig_alpha) + k = ( + 8.0 + / 3.0 + / np.pi + * (3.0 / 2.0 - self.mig_alpha) + * self.gamma_norm + * self._q_d**2.0 + * self._xs**3.0 + / (self.aspect_ratio**2.0) + * self._r_0**2.0 + * self._omega_0 + / self._nu_0 + * (self.a_p / self._r_0) ** (5.0 - 3.0 * self.mig_alpha) + ) k = np.minimum(k, 0.5) theta = (1.0 - np.sqrt(1.0 - 2.0 * k)) / k - J_p = self.M * self.a_p ** 2.0 * self.omega_k + J_p = self.M * self.a_p**2.0 * self.omega_k return 0.5 * J_p / self.gamma_tot / theta def get_static_tau_m(self): @@ -162,81 +216,93 @@ def get_static_tau_m(self): tau_m: Migration timescale [s] """ - J_p = self.M * self.a_p ** 2.0 * self.omega_k + J_p = self.M * self.a_p**2.0 * self.omega_k return 0.5 * J_p / self.gamma_tot @property def _kappa(self): """Interpolate opacity to planetary position. Adds minimal offset (for stability)""" - return ( - np.interp(self.a_p, self.disk.r, self.disk.opacity) - + 1.0e-300 - ) + return np.interp(self.a_p, self.disk.r, self.disk.opacity) + 1.0e-300 def update_torque(self): """Calculate Paardekooper2011 torques""" self.calc_gammaeff(migphase=True) - ang_p = self.a_p ** 2.0 * self.omega_k # spec. Ang.mom of planet, checked - pnu = 2.0 / 3.0 * np.sqrt((ang_p * self._xs ** 3.0) / (2.0 * np.pi * self.nu)) # checked - pxi = np.sqrt((ang_p * self._xs ** 3.0) / (2.0 * np.pi * self.Vchi)) # checked + ang_p = self.a_p**2.0 * self.omega_k # spec. Ang.mom of planet, checked + pnu = ( + 2.0 / 3.0 * np.sqrt((ang_p * self._xs**3.0) / (2.0 * np.pi * self.nu)) + ) # checked + pxi = np.sqrt((ang_p * self._xs**3.0) / (2.0 * np.pi * self.Vchi)) # checked Fpnu = 1.0 / (1.0 + (pnu / 1.3) ** 2.0) # checked Fpxi = 1.0 / (1.0 + (pxi / 1.3) ** 2.0) # checked if pnu < np.sqrt(8.0 / (45.0 * np.pi)): - Gpnu = 16.0 / 25.0 * (45.0 * np.pi / 8.0) ** (0.75) * pnu ** (1.5) # checked + Gpnu = ( + 16.0 / 25.0 * (45.0 * np.pi / 8.0) ** (0.75) * pnu ** (1.5) + ) # checked else: Gpnu = 1.0 - 9.0 / 25.0 * (8.0 / (45.0 * np.pi)) ** (1.33333) * pnu ** ( - -8.0 / 3.0 + -8.0 / 3.0 ) # checked if (pxi) < (np.sqrt(8.0 / (45.0 * np.pi))): - Gpxi = 16.0 / 25.0 * (45.0 * np.pi / 8.0) ** (0.75) * pxi ** (1.5) # checked + Gpxi = ( + 16.0 / 25.0 * (45.0 * np.pi / 8.0) ** (0.75) * pxi ** (1.5) + ) # checked else: Gpxi = 1.0 - 9.0 / 25.0 * (8.0 / (45.0 * np.pi)) ** (1.33333) * pxi ** ( - -8.0 / 3.0 + -8.0 / 3.0 ) # checked if (pnu) < (np.sqrt(28.0 / (45.0 * np.pi))): - Kpnu = 16.0 / 25.0 * (45.0 * np.pi / 28.0) ** 0.75 * pnu ** 1.5 # checked + Kpnu = 16.0 / 25.0 * (45.0 * np.pi / 28.0) ** 0.75 * pnu**1.5 # checked else: Kpnu = 1.0 - 9.0 / 25.0 * (28.0 / (45.0 * np.pi)) ** (1.33333) * pnu ** ( - -8.0 / 3.0 + -8.0 / 3.0 ) # checked if (pxi) < (np.sqrt(28.0 / (45.0 * np.pi))): - Kpxi = 16.0 / 25.0 * (45.0 * np.pi / 28.0) ** (0.75) * pxi ** (1.5) # checked + Kpxi = ( + 16.0 / 25.0 * (45.0 * np.pi / 28.0) ** (0.75) * pxi ** (1.5) + ) # checked else: Kpxi = 1.0 - 9.0 / 25.0 * (28.0 / (45.0 * np.pi)) ** (1.33333) * pxi ** ( - -8.0 / 3.0 + -8.0 / 3.0 ) # checked # Calculate torque contributions self.gamma_zero = ( - (self.massratio / self.aspect_ratio) ** 2.0 - * self.sigma_g - * self.a_p ** 4.0 - * self.omega_k ** 2.0 + (self.massratio / self.aspect_ratio) ** 2.0 + * self.sigma_g + * self.a_p**4.0 + * self.omega_k**2.0 ) # checked if self.fSigma < 1.0 and self.fSigma > 0.1: self.gamma_zero = self.gamma_zero * self.fSigma - self._m_gap = np.sqrt(self.aspect_ratio ** 5.0 * self.disk.alpha_mig * 1.0 / 0.04) * self.disk.M_s + self._m_gap = ( + np.sqrt(self.aspect_ratio**5.0 * self.disk.alpha_mig * 1.0 / 0.04) + * self.disk.M_s + ) self._iso_vs_gap = self.peb_iso / self._m_gap xfact = self.gamma_zero / self._gammaeff # checked - GammaLIN = -(2.5 + 1.7 * self.mig_beta - 0.1 * self.mig_alpha) * xfact # checked + GammaLIN = ( + -(2.5 + 1.7 * self.mig_beta - 0.1 * self.mig_alpha) * xfact + ) # checked GammaHSENT = self.xi / self._gammaeff * 7.9 * xfact # equation 5, checked - GammaCLINENT = (2.2 - 1.4 / self._gammaeff) * self.xi * xfact # equation 7, checked + GammaCLINENT = ( + (2.2 - 1.4 / self._gammaeff) * self.xi * xfact + ) # equation 7, checked GammaCLINBARO = 0.7 * (1.5 - self.mig_alpha) * xfact # equation 6, checked GammaHSBARO = 1.1 * (1.5 - self.mig_alpha) * xfact # equation 4, checked - GammaCBARO = (GammaHSBARO * Fpnu * Gpnu + (1.0 - Kpnu) * GammaCLINBARO) # checked + GammaCBARO = GammaHSBARO * Fpnu * Gpnu + (1.0 - Kpnu) * GammaCLINBARO # checked GammaCENT = ( - GammaHSENT * Fpnu * Fpxi * np.sqrt(Gpnu * Gpxi) - + np.sqrt((1.0 - Kpnu) * (1.0 - Kpxi)) * GammaCLINENT + GammaHSENT * Fpnu * Fpxi * np.sqrt(Gpnu * Gpxi) + + np.sqrt((1.0 - Kpnu) * (1.0 - Kpxi)) * GammaCLINENT ) # checked # Normalization factor to Gamma_0 as defined in Paardekooper 2011 @@ -244,85 +310,95 @@ def update_torque(self): self.gamma_norm = self.gamma_tot / self.gamma_zero # checked def calc_gammaeff(self, migphase): - """Calculate the effective gamma AND calculate the gapdepth """ + """Calculate the effective gamma AND calculate the gapdepth""" self.massratio = self.M / self.disk.M_s # checked # rename gradients to self.mig_alpha and beta as in Paardekooper 2011 # Calculate Vchi self.Vchi = ( # eq. 4 in Paardekooper 2011 - 16.0 # faktor of 4 above paardekooper - * self.disk._chem.gamma - * (self.disk._chem.gamma - 1.0) - * sr - * self.T ** 4.0 - / ( - 3.0 - * self._kappa - * self.rho_g ** 2.0 - * self.H_g ** 2 - * self.omega_k ** 2 - ) + 16.0 # faktor of 4 above paardekooper + * self.disk._chem.gamma + * (self.disk._chem.gamma - 1.0) + * sr + * self.T**4.0 + / ( + 3.0 + * self._kappa + * self.rho_g**2.0 + * self.H_g**2 + * self.omega_k**2 + ) ) # checked # calculate Q Q = ( - 2.0 - * self.Vchi - / (3.0 * (self.H_g / self.a_p) ** 3.0) - / (self.a_p ** 2.0 * self.omega_k) + 2.0 + * self.Vchi + / (3.0 * (self.H_g / self.a_p) ** 3.0) + / (self.a_p**2.0 * self.omega_k) ) # checked self._gammaeff = ( - 2.0 - * Q - * self.disk._chem.gamma - / ( - self.disk._chem.gamma * Q - + 0.5 - * np.sqrt( + 2.0 + * Q + * self.disk._chem.gamma + / ( + self.disk._chem.gamma * Q + + 0.5 + * np.sqrt( 2.0 * np.sqrt( - (self.disk._chem.gamma ** 2.0 * Q ** 2.0 + 1.0) ** 2.0 - 16.0 * Q ** 2.0 * (self.disk._chem.gamma - 1.0) + (self.disk._chem.gamma**2.0 * Q**2.0 + 1.0) ** 2.0 + - 16.0 * Q**2.0 * (self.disk._chem.gamma - 1.0) ) - + 2.0 * self.disk._chem.gamma ** 2.0 * Q ** 2.0 + + 2.0 * self.disk._chem.gamma**2.0 * Q**2.0 - 2.0 ) - ) + ) ) # checked # Calculate functions F,G,K self.aspect_ratio = self.H_g / self.a_p # relative disk thickness self._xs = ( - 1.11 / self._gammaeff ** 0.25 * np.sqrt(self.massratio / self.aspect_ratio) + 1.11 / self._gammaeff**0.25 * np.sqrt(self.massratio / self.aspect_ratio) ) # checked, dimensionless in units of r_p ####### # GAP # ####### if not hasattr(self, "_M_HS"): - self._M_HS = 4 * np.pi * self._xs * self.a_p ** 2 * self.sigma_g + self._M_HS = 4 * np.pi * self._xs * self.a_p**2 * self.sigma_g self._sigma_HS = 1 * self.sigma_g - K = self.massratio ** 2.0 / (self.aspect_ratio ** 5.0 * self.disk.alpha_mig) # checked + K = self.massratio**2.0 / ( + self.aspect_ratio**5.0 * self.disk.alpha_mig + ) # checked self.fSigma_kanagawa = 1.0 / (1.0 + 0.04 * K) # checked - R = self.a_p ** 2.0 * self.omega_k / self.nu + R = self.a_p**2.0 * self.omega_k / self.nu Gap = 3.0 / 4.0 * self.H_g / self.r_h + 50.0 / (self.massratio * R) if Gap <= 2.4646: - f_P = (Gap - 0.541) / 4. + f_P = (Gap - 0.541) / 4.0 else: - f_P = 1. - np.exp(-Gap ** 0.75 / 3.) + f_P = 1.0 - np.exp(-(Gap**0.75) / 3.0) if self.past_pebble_iso: _, m_dot_gas = self._accretion_models_dict["GasAccretion"].get_min() if migphase: - self._sigma_HS = self._sigma_HS * (self._old_a_p/self.a_p)**1.5 + self._sigma_HS = self._sigma_HS * (self._old_a_p / self.a_p) ** 1.5 self._old_a_p = 1.0 * self.a_p - M_HS_check = 4.0 * np.pi * self.a_p ** 2.0 * self._xs * self._sigma_HS + M_HS_check = 4.0 * np.pi * self.a_p**2.0 * self._xs * self._sigma_HS if M_HS_check < self._M_HS: self._M_HS = self._M_HS + (self.m_dot_disk - m_dot_gas) * self.h else: - self._M_HS = self._M_HS + (4.0 * np.pi * self.a_p ** 2.0 * self._xs - self._M_HS / self._sigma_HS) * self.sigma_g + self._M_HS = ( + self._M_HS + + ( + 4.0 * np.pi * self.a_p**2.0 * self._xs + - self._M_HS / self._sigma_HS + ) + * self.sigma_g + ) - self._sigma_HS = self._M_HS / (4.0 * np.pi * self.a_p ** 2.0 * self._xs) + self._sigma_HS = self._M_HS / (4.0 * np.pi * self.a_p**2.0 * self._xs) if self.t < self.disk.begin_photevap: assert self._M_HS >= 0.0, "negative mass in horseshoe region" @@ -335,4 +411,3 @@ def calc_gammaeff(self, migphase): self.fSigma = np.maximum(f_A * f_P, 1e-7) return - diff --git a/chemcomp/planets/migration_1.py b/chemcomp/planets/migration_1.py index 9e71af4..d0de876 100755 --- a/chemcomp/planets/migration_1.py +++ b/chemcomp/planets/migration_1.py @@ -8,18 +8,17 @@ class LindbladPlanet(Planet): """simple migrating planet. torques calculated with viscosity of α = 0.004""" def __init__( - self, - output, - disk, - accretion_models=[], - rho_c=5.5 * u.g / (u.cm ** 3.0), + self, + output, + disk, + accretion_models=[], + rho_c=5.5 * u.g / (u.cm**3.0), a_p=5.2 * u.au, t_0=0.0 * u.yr, M_c=None, M_a=None, **kwargs ): - super().__init__( a_p=a_p, M_c=M_c, @@ -34,7 +33,7 @@ def __init__( def update_a_p(self): self.update_torque() - J_p = self.M * self.a_p ** 2.0 * self.omega_k + J_p = self.M * self.a_p**2.0 * self.omega_k self.tau_m = 0.5 * J_p / self._gamma_tot self.a_p += self.a_p / self.tau_m * self.h @@ -42,7 +41,7 @@ def update_torque(self): adiabatic_index = 1.4 q = self.M / self.disk.M_s h = self.H_g / self.a_p - self._gamma_0 = (q / h) ** 2 * self.sigma_g * self.a_p ** 4 * self.omega_k ** 2 + self._gamma_0 = (q / h) ** 2 * self.sigma_g * self.a_p**4 * self.omega_k**2 tgrad = np.interp(self.a_p, self.disk.r, self.disk.grad_T) # to be defined Siggrad = np.interp(self.a_p, self.disk.r, self.disk.grad_sig) # checked self._gamma_tot = ( diff --git a/chemcomp/planets/no_accretion.py b/chemcomp/planets/no_accretion.py index 3044c87..4ff5f55 100755 --- a/chemcomp/planets/no_accretion.py +++ b/chemcomp/planets/no_accretion.py @@ -7,16 +7,16 @@ class NoAccretion(Planet): """NoAccretion planet. Useful if you are just interested in the disk and not in the planet itself""" def __init__( - self, - output, - disk, - accretion_models=[], - rho_c=5.5 * u.g / (u.cm ** 3), - a_p=5.2 * u.au, - t_0=0.0 * u.yr, - M_c=None, - M_a=None, - **kwargs, + self, + output, + disk, + accretion_models=[], + rho_c=5.5 * u.g / (u.cm**3), + a_p=5.2 * u.au, + t_0=0.0 * u.yr, + M_c=None, + M_a=None, + **kwargs, ): self.gamma_tot = 0 self.gamma_norm = 0 diff --git a/chemcomp/planets/simple_planet.py b/chemcomp/planets/simple_planet.py index 74090a3..ec2a179 100755 --- a/chemcomp/planets/simple_planet.py +++ b/chemcomp/planets/simple_planet.py @@ -7,18 +7,17 @@ class SimplePlanet(Planet): """A very simple Planet without anything.""" def __init__( - self, - output, - disk, - accretion_models=[], - rho_c=5.5 * u.g / (u.cm ** 3), - a_p=5.2 * u.au, - t_0=0.0 * u.yr, - M_c=None, - M_a=None, - **kwargs, + self, + output, + disk, + accretion_models=[], + rho_c=5.5 * u.g / (u.cm**3), + a_p=5.2 * u.au, + t_0=0.0 * u.yr, + M_c=None, + M_a=None, + **kwargs, ): - super().__init__( a_p=a_p, M_c=M_c, @@ -27,6 +26,5 @@ def __init__( rho_c=rho_c, disk=disk, accretion_models=accretion_models, - output=output - ** kwargs, + output=output**kwargs, )