Skip to content

Commit

Permalink
Ruff-ifies all top-level code
Browse files Browse the repository at this point in the history
  • Loading branch information
darothen committed Sep 26, 2023
1 parent 784938e commit 1a8fdb1
Show file tree
Hide file tree
Showing 12 changed files with 161 additions and 274 deletions.
33 changes: 11 additions & 22 deletions pyrcel/_parcel_aux_numba.py
Expand Up @@ -20,16 +20,14 @@
@nb.njit()
@auxcc.export("sigma_w", "f8(f8)")
def sigma_w(T):
"""See :func:`pyrcel.thermo.sigma_w` for full documentation
"""
"""See :func:`pyrcel.thermo.sigma_w` for full documentation"""
return 0.0761 - (1.55e-4) * (T - 273.15)


@nb.njit()
@auxcc.export("ka", "f8(f8, f8, f8)")
def ka(T, r, rho):
"""See :func:`pyrcel.thermo.ka` for full documentation
"""
"""See :func:`pyrcel.thermo.ka` for full documentation"""
ka_cont = 1e-3 * (4.39 + 0.071 * T)
denom = 1.0 + (ka_cont / (c.at * r * rho * c.Cp)) * np.sqrt(
(2 * PI * c.Ma) / (c.R * T)
Expand All @@ -40,43 +38,36 @@ def ka(T, r, rho):
@nb.njit()
@auxcc.export("dv", "f8(f8, f8, f8, f8)")
def dv(T, r, P, accom):
"""See :func:`pyrcel.thermo.dv` for full documentation
"""
"""See :func:`pyrcel.thermo.dv` for full documentation"""
P_atm = P * 1.01325e-5 # Pa -> atm
dv_cont = 1e-4 * (0.211 / P_atm) * ((T / 273.0) ** 1.94)
denom = 1.0 + (dv_cont / (accom * r)) * np.sqrt(
(2 * PI * c.Mw) / (c.R * T)
)
denom = 1.0 + (dv_cont / (accom * r)) * np.sqrt((2 * PI * c.Mw) / (c.R * T))
return dv_cont / denom


@nb.njit()
@auxcc.export("es", "f8(f8)")
def es(T):
"""See :func:`pyrcel.thermo.es` for full documentation
"""
"""See :func:`pyrcel.thermo.es` for full documentation"""
return 611.2 * np.exp(17.67 * T / (T + 243.5))


@nb.njit()
@auxcc.export("Seq", "f8(f8, f8, f8)")
def Seq(r, r_dry, T, kappa):
"""See :func:`pyrcel.thermo.Seq` for full documentation.
"""
"""See :func:`pyrcel.thermo.Seq` for full documentation."""
A = (2.0 * c.Mw * sigma_w(T)) / (c.R * T * c.rho_w * r)
B = 1.0
if kappa > 0.0:
B = (r ** 3 - (r_dry ** 3)) / (r ** 3 - (r_dry ** 3) * (1.0 - kappa))
B = (r**3 - (r_dry**3)) / (r**3 - (r_dry**3) * (1.0 - kappa))
return np.exp(A) * B - 1.0


## RHS Derivative callback function
@nb.njit(parallel=True)
@auxcc.export(
"parcel_ode_sys", "f8[:](f8[:], f8, i4, f8[:], f8[:], f8, f8[:], f8)"
)
@auxcc.export("parcel_ode_sys", "f8[:](f8[:], f8, i4, f8[:], f8[:], f8, f8[:], f8)")
def parcel_ode_sys(y, t, nr, r_drys, Nis, V, kappas, accom):
""" Calculates the instantaneous time-derivative of the parcel model system.
"""Calculates the instantaneous time-derivative of the parcel model system.
Given a current state vector `y` of the parcel model, computes the tendency
of each term including thermodynamic (pressure, temperature, etc) and aerosol
Expand Down Expand Up @@ -177,9 +168,7 @@ def parcel_ode_sys(y, t, nr, r_drys, Nis, V, kappas, accom):
) # Contribution to liq. water tendency due to growth
drs_dt[i] = dr_dt

dwc_dt *= (
4.0 * PI * c.rho_w / rho_air_dry
) # Hydrated aerosol size -> water mass
dwc_dt *= 4.0 * PI * c.rho_w / rho_air_dry # Hydrated aerosol size -> water mass
# use rho_air_dry for mixing ratio definition consistency
# No freezing implemented yet
dwi_dt = 0.0
Expand Down Expand Up @@ -215,7 +204,7 @@ def parcel_ode_sys(y, t, nr, r_drys, Nis, V, kappas, accom):
"""

## GHAN (2011)
alpha = (c.g * c.Mw * c.L) / (c.Cp * c.R * (T ** 2))
alpha = (c.g * c.Mw * c.L) / (c.Cp * c.R * (T**2))
alpha -= (c.g * c.Ma) / (c.R * T)
gamma = (P * c.Ma) / (c.Mw * pv_sat)
gamma += (c.Mw * c.L * c.L) / (c.Cp * c.R * T * T)
Expand Down
71 changes: 25 additions & 46 deletions pyrcel/activation.py
Expand Up @@ -5,11 +5,11 @@
from scipy.special import erfc

from . import constants as c
from .thermo import es, ka_cont, dv, dv_cont, sigma_w, kohler_crit
from .thermo import dv, dv_cont, es, ka_cont, kohler_crit, sigma_w


def _unpack_aerosols(aerosols):
""" Convert a list of :class:`AerosolSpecies` into lists of aerosol properties.
"""Convert a list of :class:`AerosolSpecies` into lists of aerosol properties.
Parameters
----------
Expand Down Expand Up @@ -38,10 +38,8 @@ def _unpack_aerosols(aerosols):
return dict(species=species, mus=mus, sigmas=sigmas, Ns=Ns, kappas=kappas)


def lognormal_activation(
smax, mu, sigma, N, kappa, sgi=None, T=None, approx=True
):
""" Compute the activated number/fraction from a lognormal mode
def lognormal_activation(smax, mu, sigma, N, kappa, sgi=None, T=None, approx=True):
"""Compute the activated number/fraction from a lognormal mode
Parameters
----------
Expand Down Expand Up @@ -79,7 +77,7 @@ def lognormal_activation(


def binned_activation(Smax, T, rs, aerosol, approx=False):
""" Compute the activation statistics of a given aerosol, its transient
"""Compute the activation statistics of a given aerosol, its transient
size distribution, and updraft characteristics. Following Nenes et al, 2001
also compute the kinetic limitation statistics for the aerosol.
Expand Down Expand Up @@ -159,7 +157,7 @@ def binned_activation(Smax, T, rs, aerosol, approx=False):

# noinspection PyUnresolvedReferences
def multi_mode_activation(Smax, T, aerosols, rss):
""" Compute the activation statistics of a multi-mode, binned_activation
"""Compute the activation statistics of a multi-mode, binned_activation
aerosol population.
Parameters
Expand Down Expand Up @@ -194,7 +192,7 @@ def multi_mode_activation(Smax, T, aerosols, rss):


def _vpres(T):
""" Polynomial approximation of saturated water vapour pressure as
"""Polynomial approximation of saturated water vapour pressure as
a function of temperature.
Parameters
Expand Down Expand Up @@ -231,7 +229,7 @@ def _vpres(T):


def _erfp(x):
""" Polynomial approximation to error function """
"""Polynomial approximation to error function"""
AA = [0.278393, 0.230389, 0.000972, 0.078108]
y = np.abs(1.0 * x)
axx = 1.0 + y * (AA[0] + y * (AA[1] + y * (AA[2] + y * AA[3])))
Expand Down Expand Up @@ -260,7 +258,7 @@ def mbn2014(
tol=1e-6,
max_iters=100,
):
""" Computes droplet activation using an iterative scheme.
"""Computes droplet activation using an iterative scheme.
This method implements the iterative activation scheme under development by
the Nenes' group at Georgia Tech. It encompasses modifications made over a
Expand Down Expand Up @@ -346,9 +344,7 @@ def mbn2014(
A = 4.0 * Mw * surt / R / T / Rho_w
# There are three different ways to do this:
# 1) original formula from MBN2014
f = lambda T, dpg, kappa: np.sqrt(
(A ** 3.0) * 4.0 / 27.0 / kappa / (dpg ** 3.0)
)
f = lambda T, dpg, kappa: np.sqrt((A**3.0) * 4.0 / 27.0 / kappa / (dpg**3.0))
# 2) detailed kohler calculation
# f = lambda T, dpg, kappa: kohler_crit(T, dpg/2., kappa)
# 3) approximate kohler calculation
Expand Down Expand Up @@ -379,9 +375,7 @@ def mbn2014(
)

# Setup constants used in supersaturation equation
wv_pres_sat = _vpres(T) * (
1e5 / 1e3
) # MBN2014: could also use thermo.es()
wv_pres_sat = _vpres(T) * (1e5 / 1e3) # MBN2014: could also use thermo.es()
alpha = g * Mw * L / Cp / R / T / T - g * Ma / R / T
beta1 = P * Ma / wv_pres_sat / Mw + Mw * L * L / Cp / R / T / T
beta2 = (
Expand All @@ -394,14 +388,14 @@ def mbn2014(
cf2 = A / 3.0

def _sintegral(smax):
""" Integrate the activation equation, using ``spar`` as the population
"""Integrate the activation equation, using ``spar`` as the population
splitting threshold.
Inherits the workspace thermodynamic/constant variables from one level
of scope higher.
"""

zeta_c = ((16.0 / 9.0) * alpha * V * beta2 * (A ** 2)) ** 0.25
zeta_c = ((16.0 / 9.0) * alpha * V * beta2 * (A**2)) ** 0.25
delta = 1.0 - (zeta_c / smax) ** 4.0 # spar -> smax
critical = delta <= 0.0

Expand Down Expand Up @@ -558,7 +552,7 @@ def arg2000(
kappas=[],
min_smax=False,
):
""" Computes droplet activation using a psuedo-analytical scheme.
"""Computes droplet activation using a psuedo-analytical scheme.
This method implements the psuedo-analytical scheme of [ARG2000] to
calculate droplet activation an an adiabatically ascending parcel. It
Expand Down Expand Up @@ -625,12 +619,8 @@ def arg2000(

# Originally from Abdul-Razzak 1998 w/ Ma. Need kappa formulation
wv_sat = es(T - 273.15)
alpha = (c.g * c.Mw * c.L) / (c.Cp * c.R * (T ** 2)) - (c.g * c.Ma) / (
c.R * T
)
gamma = (c.R * T) / (wv_sat * c.Mw) + (c.Mw * (c.L ** 2)) / (
c.Cp * c.Ma * T * P
)
alpha = (c.g * c.Mw * c.L) / (c.Cp * c.R * (T**2)) - (c.g * c.Ma) / (c.R * T)
gamma = (c.R * T) / (wv_sat * c.Mw) + (c.Mw * (c.L**2)) / (c.Cp * c.Ma * T * P)

# Condensation effects - base calculation
G_a = (c.rho_w * c.R * T) / (wv_sat * dv_cont(T, P) * c.Mw)
Expand All @@ -639,8 +629,7 @@ def arg2000(

Smis = []
Sparts = []
for (mu, sigma, N, kappa) in zip(mus, sigmas, Ns, kappas):

for mu, sigma, N, kappa in zip(mus, sigmas, Ns, kappas):
am = mu * 1e-6
N = N * 1e6

Expand All @@ -657,34 +646,26 @@ def arg2000(
# Scale using the formula from [GHAN2011]
# G_ac - estimate using critical radius of number mode radius,
# and new value for condensation coefficient
G_a = (c.rho_w * c.R * T) / (
wv_sat * dv(T, rc_mode, P, accom) * c.Mw
)
G_b = (c.L * c.rho_w * ((c.L * c.Mw / (c.R * T)) - 1)) / (
ka_cont(T) * T
)
G_a = (c.rho_w * c.R * T) / (wv_sat * dv(T, rc_mode, P, accom) * c.Mw)
G_b = (c.L * c.rho_w * ((c.L * c.Mw / (c.R * T)) - 1)) / (ka_cont(T) * T)
G_ac = 1.0 / (G_a + G_b)

# G_ac1 - estimate using critical radius of number mode radius,
# unity condensation coefficient; re_use G_b (no change)
G_a = (c.rho_w * c.R * T) / (
wv_sat * dv(T, rc_mode, P, accom=1.0) * c.Mw
)
G_a = (c.rho_w * c.R * T) / (wv_sat * dv(T, rc_mode, P, accom=1.0) * c.Mw)
G_ac1 = 1.0 / (G_a + G_b)

# Combine using scaling formula (40) from [GHAN2011]
G = G_0 * G_ac / G_ac1

# Parameterization integral solutions
zeta = (2.0 / 3.0) * A * (np.sqrt(alpha * V / G))
etai = ((alpha * V / G) ** (3.0 / 2.0)) / (
N * gamma * c.rho_w * 2.0 * np.pi
)
etai = ((alpha * V / G) ** (3.0 / 2.0)) / (N * gamma * c.rho_w * 2.0 * np.pi)

# Contributions to maximum supersaturation
Spa = fi * ((zeta / etai) ** (1.5))
Spb = gi * (((Smi2 ** 2) / (etai + 3.0 * zeta)) ** (0.75))
S_part = (1.0 / (Smi2 ** 2)) * (Spa + Spb)
Spb = gi * (((Smi2**2) / (etai + 3.0 * zeta)) ** (0.75))
S_part = (1.0 / (Smi2**2)) * (Spa + Spb)

Smis.append(Smi2)
Sparts.append(S_part)
Expand All @@ -700,17 +681,15 @@ def arg2000(

n_acts, act_fracs = [], []
for mu, sigma, N, kappa, sgi in zip(mus, sigmas, Ns, kappas, Smis):
N_act, act_frac = lognormal_activation(
smax, mu * 1e-6, sigma, N, kappa, sgi
)
N_act, act_frac = lognormal_activation(smax, mu * 1e-6, sigma, N, kappa, sgi)
n_acts.append(N_act)
act_fracs.append(act_frac)

return smax, n_acts, act_fracs


def shipwayabel2010(V, T, P, aerosol):
""" Activation scheme following Shipway and Abel, 2010
"""Activation scheme following Shipway and Abel, 2010
(doi:10.1016/j.atmosres.2009.10.005).
"""
Expand Down
24 changes: 8 additions & 16 deletions pyrcel/aerosol.py
Expand Up @@ -7,7 +7,7 @@


def dist_to_conc(dist, r_min, r_max, rule="trapezoid"):
""" Converts a swath of a size distribution function to an actual number
"""Converts a swath of a size distribution function to an actual number
concentration.
Aerosol size distributions are typically reported by normalizing the
Expand Down Expand Up @@ -53,7 +53,7 @@ def dist_to_conc(dist, r_min, r_max, rule="trapezoid"):


class AerosolSpecies(object):
""" Container class for organizing aerosol metadata.
"""Container class for organizing aerosol metadata.
To allow flexibility with how aerosols are defined in the model, this class is
meant to act as a wrapper to contain metadata about aerosols (their species
Expand Down Expand Up @@ -146,9 +146,7 @@ def __init__(
self.kappa = kappa # Kappa hygroscopicity parameter
self.rho = rho # aerosol density kg/m^3
self.mw = mw
self.bins = (
bins
) # Number of bins for discretizing the size distribution
self.bins = bins # Number of bins for discretizing the size distribution

# Handle the size distribution passed to the constructor
self.distribution = distribution
Expand All @@ -164,7 +162,7 @@ def __init__(
lr = (self.r_drys[0] ** 2.0) / mid1
rs = [lr, mid1]
for r_dry in self.r_drys[1:]:
rs.append(r_dry ** 2.0 / rs[-1])
rs.append(r_dry**2.0 / rs[-1])
self.rs = np.array(rs) * 1e6
else:
# Truly mono-disperse, so no boundaries (we don't actually need
Expand All @@ -185,9 +183,7 @@ def __init__(
self.rs = bins[:]
else:
if not r_min:
lr = np.log10(
distribution.mu / (10.0 * distribution.sigma)
)
lr = np.log10(distribution.mu / (10.0 * distribution.sigma))
else:
lr = np.log10(r_min)
if not r_max:
Expand Down Expand Up @@ -247,8 +243,7 @@ def __init__(

else:
raise ValueError(
"Could not work with size distribution of type %r"
% type(distribution)
"Could not work with size distribution of type %r" % type(distribution)
)

# Correct to SI units
Expand All @@ -258,7 +253,7 @@ def __init__(
self.nr = len(self.r_drys)

def stats(self):
""" Compute useful statistics about this aerosol's size distribution.
"""Compute useful statistics about this aerosol's size distribution.
Returns
-------
Expand All @@ -278,10 +273,7 @@ def stats(self):
stats_dict = self.distribution.stats()

if self.rho:

stats_dict["total_mass"] = (
stats_dict["total_volume"] * self.rho
)
stats_dict["total_mass"] = stats_dict["total_volume"] * self.rho
stats_dict["mean_mass"] = stats_dict["mean_volume"] * self.rho
stats_dict["specific_surface_area"] = (
stats_dict["total_surface_area"] / stats_dict["total_mass"]
Expand Down

0 comments on commit 1a8fdb1

Please sign in to comment.