diff --git a/PySDM/attributes/chemistry/acidity.py b/PySDM/attributes/chemistry/acidity.py
index 54dda99fa..6e07d46e5 100644
--- a/PySDM/attributes/chemistry/acidity.py
+++ b/PySDM/attributes/chemistry/acidity.py
@@ -3,7 +3,6 @@
"""
from PySDM.attributes.impl.intensive_attribute import DerivedAttribute
from PySDM.physics.aqueous_chemistry.support import AQUEOUS_COMPOUNDS
-from PySDM.physics import constants as const
from PySDM.backends.impl_numba.methods.chemistry_methods import _conc
@@ -16,14 +15,13 @@ def __init__(self, builder):
super().__init__(builder, name='pH', dependencies=self.conc.values())
self.environment = builder.particulator.environment
self.cell_id = builder.get_attribute('cell id')
- self.particles = builder.particulator
def allocate(self, idx):
super().allocate(idx)
- self.data[:] = const.pH_w
+ self.data[:] = self.formulae.constants.pH_w
def recalculate(self):
- dynamic = self.particles.dynamics['AqueousChemistry']
+ dynamic = self.particulator.dynamics['AqueousChemistry']
self.particulator.backend.equilibrate_H(
dynamic.equilibrium_consts,
diff --git a/PySDM/attributes/impl/mapper.py b/PySDM/attributes/impl/mapper.py
index 987571160..da09c6abb 100644
--- a/PySDM/attributes/impl/mapper.py
+++ b/PySDM/attributes/impl/mapper.py
@@ -21,9 +21,9 @@
'volume': lambda _: Volume,
'dry volume organic': lambda dynamics: (
make_dummy_attribute_factory('dry volume organic')
- if 'Condensation' in dynamics and isinstance(
- dynamics['Condensation'].particulator.formulae.surface_tension,
- Constant
+ if 'Condensation' in dynamics and (
+ dynamics['Condensation'].particulator.formulae.surface_tension.__name__ ==
+ Constant.__name__
)
else DryVolumeOrganic
),
@@ -31,9 +31,9 @@
DryVolumeDynamic if 'AqueousChemistry' in dynamics else DryVolume,
'dry volume organic fraction': lambda dynamics: (
make_dummy_attribute_factory('dry volume organic fraction')
- if 'Condensation' in dynamics and isinstance(
- dynamics['Condensation'].particulator.formulae.surface_tension,
- Constant
+ if 'Condensation' in dynamics and (
+ dynamics['Condensation'].particulator.formulae.surface_tension.__name__ ==
+ Constant.__name__
)
else OrganicFraction
),
diff --git a/PySDM/attributes/physics/critical_supersaturation.py b/PySDM/attributes/physics/critical_supersaturation.py
index 0e2eb6cea..39acb6afc 100644
--- a/PySDM/attributes/physics/critical_supersaturation.py
+++ b/PySDM/attributes/physics/critical_supersaturation.py
@@ -2,7 +2,6 @@
kappa-Koehler critical supersaturation calculated for actual environment temperature
"""
from PySDM.attributes.impl.derived_attribute import DerivedAttribute
-from PySDM.physics import constants as const
class CriticalSupersaturation(DerivedAttribute):
@@ -17,15 +16,13 @@ def __init__(self, builder):
name='critical supersaturation',
dependencies=(self.v_crit, self.kappa, self.v_dry, self.f_org)
)
- self.formulae = builder.particulator.formulae
- self.environment = builder.particulator.environment
def recalculate(self):
- if len(self.environment['T']) != 1:
+ if len(self.particulator.environment['T']) != 1:
raise NotImplementedError()
- temperature = self.environment['T'][0]
+ temperature = self.particulator.environment['T'][0]
r_cr = self.formulae.trivia.radius(self.v_crit.data.data)
- rd3 = self.v_dry.data.data / const.PI_4_3
+ rd3 = self.v_dry.data.data / self.formulae.constants.PI_4_3
sgm = self.formulae.surface_tension.sigma(
temperature, self.v_crit.data.data, self.v_dry.data.data, self.f_org.data.data
)
diff --git a/PySDM/attributes/physics/dry_radius.py b/PySDM/attributes/physics/dry_radius.py
index d73f9432f..04c9fa431 100644
--- a/PySDM/attributes/physics/dry_radius.py
+++ b/PySDM/attributes/physics/dry_radius.py
@@ -2,7 +2,6 @@
particle dry radius computed from dry volume
"""
from PySDM.attributes.impl.derived_attribute import DerivedAttribute
-from PySDM.physics import constants as const
class DryRadius(DerivedAttribute):
@@ -12,5 +11,5 @@ def __init__(self, builder):
super().__init__(builder, name='dry radius', dependencies=dependencies)
def recalculate(self):
- self.data.product(self.volume_dry.get(), 1 / const.PI_4_3)
+ self.data.product(self.volume_dry.get(), 1 / self.formulae.constants.PI_4_3)
self.data **= 1/3
diff --git a/PySDM/attributes/physics/radius.py b/PySDM/attributes/physics/radius.py
index f809b57d2..581510e68 100644
--- a/PySDM/attributes/physics/radius.py
+++ b/PySDM/attributes/physics/radius.py
@@ -2,7 +2,6 @@
particle wet radius (calculated from the volume)
"""
from PySDM.attributes.impl.derived_attribute import DerivedAttribute
-from PySDM.physics import constants as const
class Radius(DerivedAttribute):
@@ -13,5 +12,5 @@ def __init__(self, builder):
def recalculate(self):
self.data.idx = self.volume.data.idx
- self.data.product(self.volume.get(), 1 / const.PI_4_3)
+ self.data.product(self.volume.get(), 1 / self.formulae.constants.PI_4_3)
self.data **= 1/3
diff --git a/PySDM/backends/impl_numba/methods/chemistry_methods.py b/PySDM/backends/impl_numba/methods/chemistry_methods.py
index 1877b02ec..b0249b73c 100644
--- a/PySDM/backends/impl_numba/methods/chemistry_methods.py
+++ b/PySDM/backends/impl_numba/methods/chemistry_methods.py
@@ -6,11 +6,11 @@
import numpy as np
from PySDM.backends.impl_numba import conf
from PySDM.backends.impl_numba.toms748 import toms748_solve
-from PySDM.physics.constants import Md, R_str, Rd, K_H2O
-from PySDM.physics.aqueous_chemistry.support import HenryConsts, SPECIFIC_GRAVITY, \
+from PySDM.physics.aqueous_chemistry.support import HenryConsts, SpecificGravities, \
MASS_ACCOMMODATION_COEFFICIENTS, DIFFUSION_CONST, GASEOUS_COMPOUNDS, DISSOCIATION_FACTORS, \
KineticConsts, EquilibriumConsts, k4
from PySDM.backends.impl_common.backend_methods import BackendMethods
+from PySDM.physics.constants import K_H2O
_MAX_ITER_QUITE_CLOSE = 8
@@ -29,6 +29,7 @@ def __init__(self):
self.HENRY_CONST = HenryConsts(self.formulae)
self.KINETIC_CONST = KineticConsts(self.formulae)
self.EQUILIBRIUM_CONST = EquilibriumConsts(self.formulae)
+ self.specific_gravities = SpecificGravities(self.formulae.constants)
def dissolution(self, *, n_cell, n_threads, cell_order, cell_start_arg, idx, do_chemistry_flag,
mole_amounts, env_mixing_ratio, env_T, env_p, env_rho_d, dissociation_factors,
@@ -65,32 +66,38 @@ def dissolution(self, *, n_cell, n_threads, cell_order, cell_start_arg, idx, do_
droplet_volume=droplet_volume.data,
multiplicity=multiplicity.data,
system_type=system_type,
- specific_gravity=SPECIFIC_GRAVITY[compound],
+ specific_gravity=self.specific_gravities[compound],
alpha=MASS_ACCOMMODATION_COEFFICIENTS[compound],
diffusion_const=DIFFUSION_CONST[compound],
dissociation_factor=dissociation_factors[compound].data,
- radius=self.formulae.trivia.radius
+ radius=self.formulae.trivia.radius,
+ const=self.formulae.constants
)
@staticmethod
@numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}})
def dissolution_body(super_droplet_ids, mole_amounts, env_mixing_ratio, henrysConstant, env_p,
env_T, env_rho_d, timestep, dv, droplet_volume, multiplicity, system_type,
- specific_gravity, alpha, diffusion_const, dissociation_factor, radius):
+ specific_gravity, alpha, diffusion_const, dissociation_factor, radius,
+ const):
mole_amount_taken = 0
for i in super_droplet_ids:
- Mc = specific_gravity * Md
- Rc = R_str / Mc
- cinf = env_p / env_T / (Rd / env_mixing_ratio[0] + Rc) / Mc
+ Mc = specific_gravity * const.Md
+ Rc = const.R_str / Mc
+ cinf = env_p / env_T / (const.Rd / env_mixing_ratio[0] + Rc) / Mc
r_w = radius(volume=droplet_volume[i])
- v_avg = np.sqrt(8 * R_str * env_T / (np.pi * Mc))
+ v_avg = np.sqrt(8 * const.R_str * env_T / (np.pi * Mc))
dt_over_scale = timestep / (
4 * r_w / (3 * v_avg * alpha) +
r_w ** 2 / (3 * diffusion_const)
)
A_old = mole_amounts[i] / droplet_volume[i]
H_eff = henrysConstant * dissociation_factor[i]
- A_new = (A_old + dt_over_scale * cinf) / (1 + dt_over_scale / H_eff / R_str / env_T)
+ A_new = (
+ (A_old + dt_over_scale * cinf)
+ /
+ (1 + dt_over_scale / H_eff / const.R_str / env_T)
+ )
new_mole_amount_per_real_droplet = A_new * droplet_volume[i]
assert new_mole_amount_per_real_droplet >= 0
@@ -98,7 +105,7 @@ def dissolution_body(super_droplet_ids, mole_amounts, env_mixing_ratio, henrysCo
new_mole_amount_per_real_droplet - mole_amounts[i]
)
mole_amounts[i] = new_mole_amount_per_real_droplet
- delta_mr = mole_amount_taken * specific_gravity * Md / (dv * env_rho_d)
+ delta_mr = mole_amount_taken * specific_gravity * const.Md / (dv * env_rho_d)
assert delta_mr <= env_mixing_ratio
if system_type == 'closed':
env_mixing_ratio -= delta_mr
diff --git a/PySDM/backends/impl_numba/methods/condensation_methods.py b/PySDM/backends/impl_numba/methods/condensation_methods.py
index af4875774..f1df15a6d 100644
--- a/PySDM/backends/impl_numba/methods/condensation_methods.py
+++ b/PySDM/backends/impl_numba/methods/condensation_methods.py
@@ -5,7 +5,6 @@
from functools import lru_cache
import numpy as np
import numba
-from PySDM.physics import constants as const
from PySDM.backends.impl_numba import conf
from PySDM.backends.impl_numba.toms748 import toms748_solve
from PySDM.backends.impl_common.backend_methods import BackendMethods
@@ -142,7 +141,7 @@ def step(args, dt, n_substeps):
@staticmethod
def make_step_impl(
jit_flags, phys_pvs_C, phys_lv, calculate_ml_old, calculate_ml_new, phys_T,
- phys_p, phys_pv, phys_dthd_dt, phys_D
+ phys_p, phys_pv, phys_dthd_dt, phys_D, const
):
@numba.njit(**jit_flags)
def step_impl(
@@ -188,7 +187,7 @@ def step_impl(
return step_impl
@staticmethod
- def make_calculate_ml_old(jit_flags):
+ def make_calculate_ml_old(jit_flags, const):
@numba.njit(**jit_flags)
def calculate_ml_old(volume, multiplicity, cell_idx):
result = 0
@@ -202,7 +201,7 @@ def calculate_ml_old(volume, multiplicity, cell_idx):
@staticmethod
def make_calculate_ml_new(
jit_flags, dx_dt, volume_of_x, x, phys_r_dr_dt, phys_RH_eq, phys_sigma, radius,
- phys_lambdaK, phys_lambdaD, phys_DK, within_tolerance, max_iters, RH_rtol
+ phys_lambdaK, phys_lambdaD, phys_DK, within_tolerance, max_iters, RH_rtol, const
):
@numba.njit(**jit_flags)
def minfun(x_new, x_old, timestep, kappa, f_org, rd3, temperature, RH, lv, pvs, D, K):
@@ -339,7 +338,8 @@ def make_condensation_solver(self,
fuse=fuse,
multiplier=multiplier,
RH_rtol=RH_rtol,
- max_iters=max_iters
+ max_iters=max_iters,
+ const=self.formulae.constants
)
@staticmethod
@@ -348,18 +348,18 @@ def make_condensation_solver_impl(
fastmath, phys_pvs_C, phys_lv, phys_r_dr_dt, phys_RH_eq, phys_sigma, radius,
phys_T, phys_p, phys_pv, phys_dthd_dt, phys_lambdaK, phys_lambdaD, phys_DK, phys_D,
within_tolerance, dx_dt, volume, x, timestep, dt_range, adaptive,
- fuse, multiplier, RH_rtol, max_iters
+ fuse, multiplier, RH_rtol, max_iters, const
):
jit_flags = {**conf.JIT_FLAGS, **{'parallel': False, 'cache': False, 'fastmath': fastmath}}
- calculate_ml_old = CondensationMethods.make_calculate_ml_old(jit_flags)
+ calculate_ml_old = CondensationMethods.make_calculate_ml_old(jit_flags, const)
calculate_ml_new = CondensationMethods.make_calculate_ml_new(
jit_flags, dx_dt, volume, x, phys_r_dr_dt, phys_RH_eq, phys_sigma, radius,
- phys_lambdaK, phys_lambdaD, phys_DK, within_tolerance, max_iters, RH_rtol
+ phys_lambdaK, phys_lambdaD, phys_DK, within_tolerance, max_iters, RH_rtol, const
)
step_impl = CondensationMethods.make_step_impl(
jit_flags, phys_pvs_C, phys_lv, calculate_ml_old, calculate_ml_new,
- phys_T, phys_p, phys_pv, phys_dthd_dt, phys_D
+ phys_T, phys_p, phys_pv, phys_dthd_dt, phys_D, const
)
step_fake = CondensationMethods.make_step_fake(jit_flags, step_impl)
adapt_substeps = CondensationMethods.make_adapt_substeps(
diff --git a/PySDM/backends/impl_numba/methods/freezing_methods.py b/PySDM/backends/impl_numba/methods/freezing_methods.py
index 19a718a18..f10a6c294 100644
--- a/PySDM/backends/impl_numba/methods/freezing_methods.py
+++ b/PySDM/backends/impl_numba/methods/freezing_methods.py
@@ -6,12 +6,12 @@
from PySDM.backends.impl_common.backend_methods import BackendMethods
from ...impl_numba import conf
from ...impl_common.freezing_attributes import TimeDependentAttributes, SingularAttributes
-from ....physics import constants as const
class FreezingMethods(BackendMethods):
def __init__(self):
super().__init__()
+ const = self.formulae.constants
@numba.njit(**{**conf.JIT_FLAGS, 'fastmath': self.formulae.fastmath, 'parallel': False})
def _unfrozen(volume, i):
diff --git a/PySDM/backends/impl_numba/methods/physics_methods.py b/PySDM/backends/impl_numba/methods/physics_methods.py
index c924fd545..f88c6af08 100644
--- a/PySDM/backends/impl_numba/methods/physics_methods.py
+++ b/PySDM/backends/impl_numba/methods/physics_methods.py
@@ -4,7 +4,6 @@
import numba
from numba import prange
from PySDM.backends.impl_numba import conf
-from PySDM.physics import constants as const
from PySDM.backends.impl_common.backend_methods import BackendMethods
@@ -20,6 +19,7 @@ def __init__(self):
phys_sigma = self.formulae.surface_tension.sigma
phys_volume = self.formulae.trivia.volume
phys_r_cr = self.formulae.hygroscopicity.r_cr
+ const = self.formulae.constants
@numba.njit(**{**conf.JIT_FLAGS, 'fastmath': self.formulae.fastmath})
def explicit_euler_body(y, dt, dy_dt):
diff --git a/PySDM/backends/impl_numba/test_helpers/bdf.py b/PySDM/backends/impl_numba/test_helpers/bdf.py
index be1ddd000..c6b59db9f 100644
--- a/PySDM/backends/impl_numba/test_helpers/bdf.py
+++ b/PySDM/backends/impl_numba/test_helpers/bdf.py
@@ -8,7 +8,7 @@
import numpy as np
import numba
import scipy.integrate
-import PySDM.physics.constants as const
+from PySDM.physics.constants_defaults import rho_w, K0, T0, PI_4_3
from PySDM.backends import Numba
from PySDM.backends.impl_numba.conf import JIT_FLAGS
@@ -80,7 +80,7 @@ def _make_solve(formulae):
@numba.njit(**{**JIT_FLAGS, **{'parallel': False}})
def _ql(n, x, m_d_mean):
- return np.sum(n * volume(x)) * const.rho_w / m_d_mean
+ return np.sum(n * volume(x)) * rho_w / m_d_mean
@numba.njit(**{**JIT_FLAGS, **{'parallel': False}})
def _impl(dy_dt, x, T, p, n, RH, kappa, f_org, dry_volume, thd,
@@ -92,16 +92,16 @@ def _impl(dy_dt, x, T, p, n, RH, kappa, f_org, dry_volume, thd,
v = volume(x_i)
r = phys_radius(v)
Dr = phys_DK(DTp, r, lambdaD)
- Kr = phys_DK(const.K0, r, lambdaK)
+ Kr = phys_DK(K0, r, lambdaK)
sgm = sigma(T, v, dry_volume[i], f_org[i])
dy_dt[idx_x + i] = dx_dt(
x_i,
r_dr_dt(
- RH_eq(r, T, kappa[i], dry_volume[i] / const.PI_4_3, sgm),
+ RH_eq(r, T, kappa[i], dry_volume[i] / PI_4_3, sgm),
T, RH, lv, pvs, Dr, Kr
)
)
- dqv_dt = dot_qv - np.sum(n * volume(x) * dy_dt[idx_x:]) * const.rho_w / m_d_mean
+ dqv_dt = dot_qv - np.sum(n * volume(x) * dy_dt[idx_x:]) * rho_w / m_d_mean
dy_dt[idx_thd] = dot_thd + phys_dthd_dt(rhod_mean, thd, T, dqv_dt, lv)
@numba.njit(**{**JIT_FLAGS, **{'parallel': False}})
@@ -114,7 +114,7 @@ def _odesys(t, y, kappa, f_org, dry_volume, n,
T = phys_T(rhod_mean, thd)
p = phys_p(rhod_mean, T, qv)
pv = phys_pv(p, qv)
- pvs = pvs_C(T - const.T0)
+ pvs = pvs_C(T - T0)
RH = pv / pvs
dy_dt = np.empty_like(y)
@@ -160,7 +160,7 @@ def solve(
m_new = 0
for i in range(n_sd_in_cell):
v_new = volume(y1[idx_x + i])
- m_new += n[cell_idx[i]] * v_new * const.rho_w
+ m_new += n[cell_idx[i]] * v_new * rho_w
v[cell_idx[i]] = v_new
return integ.success, qt - m_new / m_d_mean, y1[idx_thd], 1, 1, 1, 1, np.nan
diff --git a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py
index 8c58890c0..225bb17f9 100644
--- a/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py
+++ b/PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py
@@ -2,7 +2,6 @@
GPU implementation of backend methods for water condensation/evaporation
"""
from typing import Dict, Optional
-from PySDM.physics import constants as const
from PySDM.backends.impl_thrust_rtc.conf import NICE_THRUST_FLAGS
from PySDM.backends.impl_thrust_rtc.nice_thrust import nice_thrust
from PySDM.backends.impl_thrust_rtc.bisection import BISECTION
@@ -27,6 +26,7 @@ def __init__(self):
self.dqv_dt_pred: Optional[StorageBase] = None
self.rhod_mean: Optional[StorageBase] = None
self.vars: Optional[Dict[str, StorageBase]] = None
+ const = self.formulae.constants
self.__calculate_m_l = trtc.For(("ml", "v", "n", "cell_id"), "i", f'''
atomicAdd((real_type*) &ml[cell_id[i]], n[i] * v[i] * {const.rho_w});
diff --git a/PySDM/dynamics/aqueous_chemistry.py b/PySDM/dynamics/aqueous_chemistry.py
index 0bd05d1d1..0edbfe520 100644
--- a/PySDM/dynamics/aqueous_chemistry.py
+++ b/PySDM/dynamics/aqueous_chemistry.py
@@ -4,7 +4,7 @@
from collections import namedtuple
import numpy as np
from PySDM.physics.aqueous_chemistry.support import DIFFUSION_CONST, AQUEOUS_COMPOUNDS, \
- GASEOUS_COMPOUNDS, SPECIFIC_GRAVITY, M
+ GASEOUS_COMPOUNDS, SpecificGravities, M
DEFAULTS = namedtuple("_", ('pH_min', 'pH_max', 'pH_rtol', 'ionic_strength_threshold'))(
@@ -40,16 +40,18 @@ def __init__(self, environment_mole_fractions, system_type, n_substep, dry_rho,
self.equilibrium_consts = {}
self.dissociation_factors = {}
self.do_chemistry_flag = None
+ self.specific_gravities = None
def register(self, builder):
self.particulator = builder.particulator
+ self.specific_gravities = SpecificGravities(self.particulator.formulae.constants)
for key, compound in GASEOUS_COMPOUNDS.items():
shape = (1,)
self.environment_mixing_ratios[compound] = np.full(
shape,
self.particulator.formulae.trivia.mole_fraction_2_mixing_ratio(
- self.environment_mole_fractions[compound], SPECIFIC_GRAVITY[compound])
+ self.environment_mole_fractions[compound], self.specific_gravities[compound])
)
self.environment_mole_fractions = None
diff --git a/PySDM/environments/parcel.py b/PySDM/environments/parcel.py
index a1cb70772..8a185a04a 100644
--- a/PySDM/environments/parcel.py
+++ b/PySDM/environments/parcel.py
@@ -6,18 +6,15 @@
from PySDM.initialisation.equilibrate_wet_radii import equilibrate_wet_radii, default_rtol
from PySDM.initialisation.discretise_multiplicities import discretise_multiplicities
from PySDM.environments.impl.moist import Moist
-from ..physics import constants as const
class Parcel(Moist):
-
def __init__(
self, dt,
mass_of_dry_air: float,
p0: float, q0: float, T0: float,
w: [float, callable],
z0: float = 0,
- g=const.g_std,
mixed_phase=False
):
super().__init__(
@@ -32,7 +29,6 @@ def __init__(
self.T0 = T0
self.z0 = z0
self.mass_of_dry_air = mass_of_dry_air
- self.g = g
self.w = w if callable(w) else lambda _: w
@@ -99,7 +95,8 @@ def advance_parcel_vars(self):
dql_dz = self.dql / dz_dt / dt
lv = self.formulae.latent_heat.lv(T)
- drho_dz = self.formulae.hydrostatics.drho_dz(self.g, p, T, qv, lv, dql_dz=dql_dz)
+ drho_dz = self.formulae.hydrostatics.drho_dz(
+ self.formulae.constants.g_std, p, T, qv, lv, dql_dz=dql_dz)
drhod_dz = drho_dz
self.particulator.backend.explicit_euler(self._tmp['t'], dt, 1)
diff --git a/PySDM/formulae.py b/PySDM/formulae.py
index 2e7369f75..07e11b98e 100644
--- a/PySDM/formulae.py
+++ b/PySDM/formulae.py
@@ -2,94 +2,24 @@
Logic for enabling common CPU/GPU physics formulae code
"""
import inspect
+from typing import Optional
+from types import SimpleNamespace
import re
from functools import lru_cache, partial
-
+from collections import namedtuple
+import numbers
import numba
+import numpy as np
+import pint
+
from PySDM import physics
from PySDM.backends.impl_numba import conf
-# noinspection PyUnresolvedReferences
-from PySDM.physics import constants as const
-
-
-def _formula(func=None, **kw):
- if func is None:
- return numba.njit(
- **{
- **conf.JIT_FLAGS,
- **{'parallel': False, 'inline': 'always', 'cache': False, **kw}
- }
- )
- return numba.njit(
- func,
- **{
- **conf.JIT_FLAGS,
- **{'parallel': False, 'inline': 'always', 'cache': False, **kw}
- }
- )
-
-
-def _boost(obj, fastmath):
- if not physics.impl.flag.DIMENSIONAL_ANALYSIS:
- for item in dir(obj):
- if item.startswith('__'):
- continue
- attr = getattr(obj, item)
- if callable(attr):
- formula = _formula(attr, fastmath=fastmath)
- setattr(obj, item, formula)
- setattr(getattr(obj, item), 'c_inline', partial(_c_inline, fun=formula))
- return obj
-
-
-def _c_inline(fun, return_type=None, **args):
- real_t = 'real_type'
- return_type = return_type or real_t
- prae = r"([,+\-*/( ]|^)"
- post = r"([ )/*\-+,]|$)"
- real_fmt = ".32g"
- source = ''
- for line in inspect.getsourcelines(fun)[0]:
- stripped = line.strip()
- if stripped.startswith('@'):
- continue
- if stripped.startswith('//'):
- continue
- if stripped.startswith('def '):
- continue
- source += stripped
- source = source.replace("power(", "pow(")
- source = re.sub("^return ", "", source)
- for arg in inspect.signature(fun).parameters:
- source = re.sub(f"{prae}({arg}){post}", f"\\1({real_t})({args[arg]})\\3", source)
- source = re.sub(
- f"{prae}const\\.([^\\d\\W]\\w*]*){post}",
- "\\1(" + real_t + ")({const.\\2:" + real_fmt + "})\\3",
- source
- )
- source = eval(f'f"""{source}"""') # pylint: disable=eval-used
- return f'({return_type})({source})'
-
-
-def _pick(value: str, choices: dict):
- for name, cls in choices.items():
- if name == value:
- return cls()
- raise ValueError(f"Unknown setting: '{value}'; choices are: {tuple(choices.keys())}")
-
-
-def _choices(module):
- return {name: cls for name, cls in module.__dict__.items() if isinstance(cls, type)}
-
-
-@lru_cache()
-def _magick(value, module, fastmath):
- return _boost(_pick(value, _choices(module)), fastmath)
class Formulae:
def __init__(self, *,
- seed: int = const.default_random_seed,
+ constants: Optional[dict] = None,
+ seed: int = physics.constants.default_random_seed,
fastmath: bool = True,
condensation_coordinate: str = 'VolumeLogarithm',
saturation_vapour_pressure: str = 'FlatauWalkoCotton',
@@ -106,54 +36,177 @@ def __init__(self, *,
freezing_temperature_spectrum: str = 'Null',
heterogeneous_ice_nucleation_rate: str = 'Null'
):
+ constants_defaults = {
+ k: getattr(physics.constants_defaults, k)
+ for k in dir(physics.constants_defaults)
+ if isinstance(getattr(physics.constants_defaults, k), (numbers.Number, pint.Quantity))
+ }
+ constants = namedtuple(
+ "Constants",
+ tuple(constants_defaults.keys())
+ )(
+ **{**constants_defaults, **(constants or {})}
+ )
+ self.constants = constants
self.seed = seed
self.fastmath = fastmath
+ dimensional_analysis = physics.impl.flag.DIMENSIONAL_ANALYSIS
- self.trivia = _magick('Trivia', physics.trivia, fastmath)
+ self.trivia = _magick('Trivia', physics.trivia, fastmath, constants, dimensional_analysis)
self.condensation_coordinate = _magick(
- condensation_coordinate, physics.condensation_coordinate, fastmath)
+ condensation_coordinate,
+ physics.condensation_coordinate, fastmath, constants, dimensional_analysis)
self.saturation_vapour_pressure = _magick(
- saturation_vapour_pressure, physics.saturation_vapour_pressure, fastmath)
+ saturation_vapour_pressure,
+ physics.saturation_vapour_pressure, fastmath, constants, dimensional_analysis)
self.latent_heat = _magick(
- latent_heat, physics.latent_heat, fastmath)
+ latent_heat,
+ physics.latent_heat, fastmath, constants, dimensional_analysis)
self.hygroscopicity = _magick(
- hygroscopicity, physics.hygroscopicity, fastmath)
+ hygroscopicity,
+ physics.hygroscopicity, fastmath, constants, dimensional_analysis)
self.drop_growth = _magick(
- drop_growth, physics.drop_growth, fastmath)
+ drop_growth,
+ physics.drop_growth, fastmath, constants, dimensional_analysis)
self.surface_tension = _magick(
- surface_tension, physics.surface_tension, fastmath)
+ surface_tension,
+ physics.surface_tension, fastmath, constants, dimensional_analysis)
self.diffusion_kinetics = _magick(
- diffusion_kinetics, physics.diffusion_kinetics, fastmath)
+ diffusion_kinetics,
+ physics.diffusion_kinetics, fastmath, constants, dimensional_analysis)
self.diffusion_thermics = _magick(
- diffusion_thermics, physics.diffusion_thermics, fastmath)
+ diffusion_thermics,
+ physics.diffusion_thermics, fastmath, constants, dimensional_analysis)
self.ventilation = _magick(
- ventilation, physics.ventilation, fastmath)
+ ventilation,
+ physics.ventilation, fastmath, constants, dimensional_analysis)
self.state_variable_triplet = _magick(
- state_variable_triplet, physics.state_variable_triplet, fastmath)
+ state_variable_triplet,
+ physics.state_variable_triplet, fastmath, constants, dimensional_analysis)
self.particle_advection = _magick(
- particle_advection, physics.particle_advection, fastmath)
+ particle_advection,
+ physics.particle_advection, fastmath, constants, dimensional_analysis)
self.hydrostatics = _magick(
- hydrostatics, physics.hydrostatics, fastmath)
+ hydrostatics,
+ physics.hydrostatics, fastmath, constants, dimensional_analysis)
self.freezing_temperature_spectrum = _magick(
- freezing_temperature_spectrum, physics.freezing_temperature_spectrum, fastmath)
+ freezing_temperature_spectrum,
+ physics.freezing_temperature_spectrum, fastmath, constants, dimensional_analysis)
self.heterogeneous_ice_nucleation_rate = _magick(
- heterogeneous_ice_nucleation_rate, physics.heterogeneous_ice_nucleation_rate, fastmath)
- self.__check()
+ heterogeneous_ice_nucleation_rate,
+ physics.heterogeneous_ice_nucleation_rate, fastmath, constants, dimensional_analysis)
def __str__(self):
description = []
for attr in dir(self):
if not attr.startswith('_'):
- if getattr(self, attr).__class__ in (bool, int, float):
- value = getattr(self, attr)
+ attr_value = getattr(self, attr)
+ if attr_value.__class__ in (bool, int, float):
+ value = attr_value
+ elif attr_value.__class__.__name__ == 'Constants':
+ value = str(attr_value)
+ elif attr_value.__class__ in (SimpleNamespace,):
+ value = attr_value.__name__
else:
- value = getattr(self, attr).__class__.__name__
+ value = attr_value.__class__.__name__
description.append(f"{attr}: {value}")
return ', '.join(description)
- def __check(self):
- for attr in dir(self):
- if not attr.startswith('__'):
- if hasattr(getattr(self, attr), '_check'):
- getattr(self, attr)._check()
+
+def _formula(func, constants, dimensional_analysis, **kw):
+ if dimensional_analysis:
+ first_param = tuple(inspect.signature(func).parameters.keys())[0]
+ if first_param in ('_', 'const'):
+ return partial(func, constants)
+ return func
+
+ source = "class _:\n"
+ for line in inspect.getsourcelines(func)[0]:
+ source += f"{line}\n"
+ loc = {}
+ for arg_name in ('_', 'const'):
+ source = source.replace(f'def {func.__name__}({arg_name},', f'def {func.__name__}(')
+ exec( # pylint:disable=exec-used
+ source,
+ {'const': constants, 'np': np},
+ loc
+ )
+ return numba.njit(
+ getattr(loc['_'], func.__name__),
+ **{
+ **conf.JIT_FLAGS,
+ **{'parallel': False, 'inline': 'always', 'cache': False, **kw}
+ }
+ )
+
+
+def _boost(obj, fastmath, constants, dimensional_analysis):
+ formulae = {'__name__': obj.__class__.__name__}
+ for item in dir(obj):
+ attr = getattr(obj, item)
+ if item.startswith('__') or not callable(attr):
+ pass
+ else:
+ formula = _formula(attr, constants=constants, fastmath=fastmath,
+ dimensional_analysis=dimensional_analysis)
+ setattr(
+ formula,
+ 'c_inline',
+ partial(_c_inline, constants=constants, fun=attr)
+ )
+ formulae[attr.__name__] = formula
+ return SimpleNamespace(**formulae)
+
+
+def _c_inline(fun, return_type=None, constants=None, **args):
+ real_t = 'real_type'
+ return_type = return_type or real_t
+ prae = r"([,+\-*/( ]|^)"
+ post = r"([ )/*\-+,]|$)"
+ real_fmt = ".32g"
+ source = ''
+ for line in inspect.getsourcelines(fun)[0]:
+ stripped = line.strip()
+ if stripped.startswith('@'):
+ continue
+ if stripped.startswith('//'):
+ continue
+ if stripped.startswith('def '):
+ continue
+ source += stripped
+ source = source.replace("np.power(", "np.pow(")
+ source = source.replace("np.", "")
+ source = re.sub("^return ", "", source)
+ for arg in inspect.signature(fun).parameters:
+ if arg not in ('_', 'const'):
+ source = re.sub(f"{prae}({arg}){post}", f"\\1({real_t})({args[arg]})\\3", source)
+ source = re.sub(
+ f"{prae}const\\.([^\\d\\W]\\w*]*){post}",
+ "\\1(" + real_t + ")({constants.\\2:" + real_fmt + "})\\3",
+ source
+ )
+ assert constants
+ source = eval(f'f"""{source}"""') # pylint: disable=eval-used
+ return f'({return_type})({source})'
+
+
+def _pick(value: str, choices: dict, constants: namedtuple):
+ for name, cls in choices.items():
+ if name == value:
+ return cls(constants)
+ raise ValueError(f"Unknown setting: '{value}'; choices are: {tuple(choices.keys())}")
+
+
+def _choices(module):
+ return {name: cls for name, cls in module.__dict__.items() if isinstance(cls, type)}
+
+
+@lru_cache()
+def _magick(value, module, fastmath, constants, dimensional_analysis):
+ return _boost(
+ _pick(value, _choices(module), constants),
+ fastmath,
+ constants,
+ dimensional_analysis
+ )
diff --git a/PySDM/initialisation/equilibrate_wet_radii.py b/PySDM/initialisation/equilibrate_wet_radii.py
index 05bb3621c..1372ae5c6 100644
--- a/PySDM/initialisation/equilibrate_wet_radii.py
+++ b/PySDM/initialisation/equilibrate_wet_radii.py
@@ -4,7 +4,6 @@
import numba
import numpy as np
from ..backends.impl_numba.toms748 import toms748_solve
-from ..physics import constants as const
from ..backends.impl_numba.conf import JIT_FLAGS
from ..backends.impl_numba.warnings import warn
@@ -22,6 +21,7 @@ def equilibrate_wet_radii(r_dry: np.ndarray, environment,
if f_org is None:
f_org = np.zeros_like(r_dry, dtype=float)
+ const = environment.particulator.formulae.constants
T = environment["T"].to_ndarray()
RH = environment["RH"].to_ndarray()
diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py
index 7c4400b33..9c99e2733 100644
--- a/PySDM/physics/__init__.py
+++ b/PySDM/physics/__init__.py
@@ -6,3 +6,5 @@
ventilation, state_variable_triplet, trivia, particle_advection, hydrostatics,
freezing_temperature_spectrum, heterogeneous_ice_nucleation_rate)
from .constants import si
+from . import impl
+from . import constants_defaults
diff --git a/PySDM/physics/aqueous_chemistry/support.py b/PySDM/physics/aqueous_chemistry/support.py
index 1f6dfe407..07480bd3e 100644
--- a/PySDM/physics/aqueous_chemistry/support.py
+++ b/PySDM/physics/aqueous_chemistry/support.py
@@ -5,7 +5,7 @@
from chempy import Substance
import numpy as np
from PySDM.physics import si
-from PySDM.physics.constants import R_str, ROOM_TEMP, H_u, dT_u, M, Md, K_H2O
+from PySDM.physics.constants import K_H2O, M
class EqConst:
@@ -23,7 +23,7 @@ class KinConst:
def __init__(self, formulae, k, dT, T_0):
self.formulae = formulae
self.Ea = formulae.trivia.tdep2enthalpy(dT)
- self.A = k * np.exp(self.Ea / (R_str * T_0))
+ self.A = k * np.exp(self.Ea / (self.formulae.constants.R_str * T_0))
def at(self, T):
return self.formulae.trivia.arrhenius(self.A, self.Ea, T)
@@ -31,27 +31,31 @@ def at(self, T):
class HenryConsts:
def __init__(self, formulae):
+ const = formulae.constants
+ T0 = const.ROOM_TEMP
self.HENRY_CONST = {
- "HNO3": EqConst(formulae, 2.1e5 * H_u, 8700 * dT_u, T_0=ROOM_TEMP),
- "H2O2": EqConst(formulae, 7.45e4 * H_u, 7300 * dT_u, T_0=ROOM_TEMP),
- "NH3": EqConst(formulae, 62 * H_u, 4110 * dT_u, T_0=ROOM_TEMP),
- "SO2": EqConst(formulae, 1.23 * H_u, 3150 * dT_u, T_0=ROOM_TEMP),
- "CO2": EqConst(formulae, 3.4e-2 * H_u, 2440 * dT_u, T_0=ROOM_TEMP),
- "O3": EqConst(formulae, 1.13e-2 * H_u, 2540 * dT_u, T_0=ROOM_TEMP),
+ "HNO3": EqConst(formulae, 2.1e5 * const.H_u, 8700 * const.dT_u, T_0=T0),
+ "H2O2": EqConst(formulae, 7.45e4 * const.H_u, 7300 * const.dT_u, T_0=T0),
+ "NH3": EqConst(formulae, 62 * const.H_u, 4110 * const.dT_u, T_0=T0),
+ "SO2": EqConst(formulae, 1.23 * const.H_u, 3150 * const.dT_u, T_0=T0),
+ "CO2": EqConst(formulae, 3.4e-2 * const.H_u, 2440 * const.dT_u, T_0=T0),
+ "O3": EqConst(formulae, 1.13e-2 * const.H_u, 2540 * const.dT_u, T_0=T0),
}
# Table 4 in Kreidenweis et al. 2003
class EquilibriumConsts:
def __init__(self, formulae):
+ const = formulae.constants
+ T0 = const.ROOM_TEMP
self.EQUILIBRIUM_CONST = { # Reaction Specific units, K
- "K_HNO3": EqConst(formulae, 15.4 * M, 8700 * dT_u, T_0=ROOM_TEMP),
- "K_SO2": EqConst(formulae, 1.3e-2 * M, 1960 * dT_u, T_0=ROOM_TEMP),
- "K_NH3": EqConst(formulae, 1.7e-5 * M, -450 * dT_u, T_0=ROOM_TEMP),
- "K_CO2": EqConst(formulae, 4.3e-7 * M, -1000 * dT_u, T_0=ROOM_TEMP),
- "K_HSO3": EqConst(formulae, 6.6e-8 * M, 1500 * dT_u, T_0=ROOM_TEMP),
- "K_HCO3": EqConst(formulae, 4.68e-11 * M, -1760 * dT_u, T_0=ROOM_TEMP),
- "K_HSO4": EqConst(formulae, 1.2e-2 * M, 2720 * dT_u, T_0=ROOM_TEMP),
+ "K_HNO3": EqConst(formulae, 15.4 * const.M, 8700 * const.dT_u, T_0=T0),
+ "K_SO2": EqConst(formulae, 1.3e-2 * const.M, 1960 * const.dT_u, T_0=T0),
+ "K_NH3": EqConst(formulae, 1.7e-5 * const.M, -450 * const.dT_u, T_0=T0),
+ "K_CO2": EqConst(formulae, 4.3e-7 * const.M, -1000 * const.dT_u, T_0=T0),
+ "K_HSO3": EqConst(formulae, 6.6e-8 * const.M, 1500 * const.dT_u, T_0=T0),
+ "K_HCO3": EqConst(formulae, 4.68e-11 * const.M, -1760 * const.dT_u, T_0=T0),
+ "K_HSO4": EqConst(formulae, 1.2e-2 * const.M, 2720 * const.dT_u, T_0=T0),
}
@@ -110,22 +114,31 @@ def __init__(self, formulae):
class KineticConsts:
def __init__(self, formulae):
+ const = formulae.constants
+ T0 = const.ROOM_TEMP
self.KINETIC_CONST = {
- "k0": KinConst(formulae, k=2.4e4 / si.s / M, dT=0 * dT_u, T_0=ROOM_TEMP),
- "k1": KinConst(formulae, k=3.5e5 / si.s / M, dT=-5530 * dT_u, T_0=ROOM_TEMP),
- "k2": KinConst(formulae, k=1.5e9 / si.s / M, dT=-5280 * dT_u, T_0=ROOM_TEMP),
+ "k0": KinConst(formulae, k=2.4e4 / si.s / M, dT=0 * const.dT_u, T_0=T0),
+ "k1": KinConst(formulae, k=3.5e5 / si.s / M, dT=-5530 * const.dT_u, T_0=T0),
+ "k2": KinConst(formulae, k=1.5e9 / si.s / M, dT=-5280 * const.dT_u, T_0=T0),
# Different unit due to a different pseudo-order of kinetics
- "k3": KinConst(formulae, k=7.45e7 / si.s / M / M, dT=-4430 * dT_u, T_0=ROOM_TEMP),
+ "k3": KinConst(formulae, k=7.45e7 / si.s / M / M, dT=-4430 * const.dT_u, T_0=T0),
}
k4 = 13 / M
-SPECIFIC_GRAVITY = {
- compound: Substance.from_formula(compound).mass * si.gram / si.mole / Md
- for compound in GASEOUS_COMPOUNDS.values()
-}
-for compounds in AQUEOUS_COMPOUNDS.values():
- for compound in compounds:
- SPECIFIC_GRAVITY[compound] = Substance.from_formula(compound).mass * si.gram / si.mole / Md
+class SpecificGravities:
+ def __init__(self, constants):
+ self._values = {
+ compound: Substance.from_formula(compound).mass * si.gram / si.mole / constants.Md
+ for compound in GASEOUS_COMPOUNDS.values()
+ }
+
+ for compounds in AQUEOUS_COMPOUNDS.values():
+ for compound in compounds:
+ self._values[compound] = \
+ Substance.from_formula(compound).mass * si.gram / si.mole / constants.Md
+
+ def __getitem__(self, item):
+ return self._values[item]
diff --git a/PySDM/physics/condensation_coordinate/volume.py b/PySDM/physics/condensation_coordinate/volume.py
index 252197175..1b4e8cec0 100644
--- a/PySDM/physics/condensation_coordinate/volume.py
+++ b/PySDM/physics/condensation_coordinate/volume.py
@@ -1,19 +1,21 @@
"""
particle volume as condensation coordinate (i.e. no transformation)
"""
-from numpy import power
-from PySDM.physics import constants as const
+import numpy as np
class Volume:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def dx_dt(x, r_dr_dt):
- return 4 * const.PI * power(x / const.PI_4_3, const.ONE_THIRD) * r_dr_dt
+ def dx_dt(const, x, r_dr_dt):
+ return 4 * const.PI * np.power(x / const.PI_4_3, const.ONE_THIRD) * r_dr_dt
@staticmethod
- def volume(x):
+ def volume(_, x):
return x
@staticmethod
- def x(volume):
+ def x(_, volume):
return volume
diff --git a/PySDM/physics/condensation_coordinate/volume_logarithm.py b/PySDM/physics/condensation_coordinate/volume_logarithm.py
index 9f7d5dc78..88d16370e 100644
--- a/PySDM/physics/condensation_coordinate/volume_logarithm.py
+++ b/PySDM/physics/condensation_coordinate/volume_logarithm.py
@@ -1,19 +1,21 @@
"""
logarithm of particle volume as coordinate (ensures non-negative values)
"""
-from numpy import log, exp, power
-from PySDM.physics import constants as const
+import numpy as np
class VolumeLogarithm:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def dx_dt(x, r_dr_dt):
- return exp(-const.TWO_THIRDS * x) * r_dr_dt * 3 * power(const.PI_4_3, const.TWO_THIRDS)
+ def dx_dt(const, x, r_dr_dt):
+ return np.exp(-const.TWO_THIRDS * x) * r_dr_dt * 3*np.power(const.PI_4_3, const.TWO_THIRDS)
@staticmethod
def volume(x):
- return exp(x)
+ return np.exp(x)
@staticmethod
def x(volume):
- return log(volume)
+ return np.log(volume)
diff --git a/PySDM/physics/constants.py b/PySDM/physics/constants.py
index 7f182feb4..2cd40baca 100644
--- a/PySDM/physics/constants.py
+++ b/PySDM/physics/constants.py
@@ -1,13 +1,13 @@
"""
-Collection of physical constants with dimensional analysis handled with
+Collection of constants with dimensional analysis handled with
[Pint](https://pypi.org/project/Pint/)'s package `UnitRegistry` for test
purposes and mocked with `PySDM.physics.impl.fake_unit_registry.FakeUnitRegistry` by default.
"""
import os
import time
+
import pint
from scipy import constants as sci
-from chempy import Substance
from .impl.fake_unit_registry import FakeUnitRegistry
from .impl.flag import DIMENSIONAL_ANALYSIS
@@ -26,84 +26,17 @@ def convert_to(value, unit):
ONE_THIRD = 1 / 3
TWO_THIRDS = 2 / 3
-Md = (
- 0.78 * Substance.from_formula('N2').mass * si.gram / si.mole +
- 0.21 * Substance.from_formula('O2').mass * si.gram / si.mole +
- 0.01 * Substance.from_formula('Ar').mass * si.gram / si.mole
+default_random_seed = (
+ 44 if 'CI' in os.environ # https://en.wikipedia.org/wiki/44_(number)
+ else time.time_ns()
)
-Mv = Substance.from_formula('H2O').mass * si.gram / si.mole
-
-R_str = sci.R * si.joule / si.kelvin / si.mole
-eps = Mv / Md
-Rd = R_str / Md
-Rv = R_str / Mv
-
-D0 = 2.26e-5 * si.metre ** 2 / si.second
-D_exp = 1.81
-
-K0 = 2.4e-2 * si.joules / si.metres / si.seconds / si.kelvins
-
-p1000 = 1000 * si.hectopascals
-c_pd = 1005 * si.joule / si.kilogram / si.kelvin
-c_pv = 1850 * si.joule / si.kilogram / si.kelvin
-T0 = sci.zero_Celsius * si.kelvin
-g_std = sci.g * si.metre / si.second ** 2
-
-c_pw = 4218 * si.joule / si.kilogram / si.kelvin
-
-Rd_over_c_pd = Rd / c_pd
-
-ARM_C1 = 6.1094 * si.hectopascal
-ARM_C2 = 17.625 * si.dimensionless
-ARM_C3 = 243.04 * si.kelvin
-
-FWC_C0 = 6.115836990e000 * si.hPa
-FWC_C1 = 0.444606896e000 * si.hPa / si.K
-FWC_C2 = 0.143177157e-01 * si.hPa / si.K**2
-FWC_C3 = 0.264224321e-03 * si.hPa / si.K**3
-FWC_C4 = 0.299291081e-05 * si.hPa / si.K**4
-FWC_C5 = 0.203154182e-07 * si.hPa / si.K**5
-FWC_C6 = 0.702620698e-10 * si.hPa / si.K**6
-FWC_C7 = 0.379534310e-13 * si.hPa / si.K**7
-FWC_C8 = -.321582393e-15 * si.hPa / si.K**8
-
-FWC_I0 = 6.098689930e000 * si.hPa
-FWC_I1 = 0.499320233e000 * si.hPa / si.K
-FWC_I2 = 0.184672631e-01 * si.hPa / si.K**2
-FWC_I3 = 0.402737184e-03 * si.hPa / si.K**3
-FWC_I4 = 0.565392987e-05 * si.hPa / si.K**4
-FWC_I5 = 0.521693933e-07 * si.hPa / si.K**5
-FWC_I6 = 0.307839583e-09 * si.hPa / si.K**6
-FWC_I7 = 0.105785160e-11 * si.hPa / si.K**7
-FWC_I8 = 0.161444444e-14 * si.hPa / si.K**8
-
-rho_w = 1 * si.kilograms / si.litres
-rho_i = 916.8 * si.kg / si.metres**3
-pH_w = 7
-sgm_w = 0.072 * si.joule / si.metre ** 2
-nu_w = Mv / rho_w
-
-p_tri = 611.73 * si.pascal
-T_tri = 273.16 * si.kelvin
-l_tri = 2.5e6 * si.joule / si.kilogram
-
-# standard pressure and temperature (ICAO)
-T_STP = T0 + 15 * si.kelvin
-p_STP = 101325 * si.pascal
-rho_STP = p_STP / Rd / T_STP
PPT = 1e-12
PPB = 1e-9
PPM = 1e-6
-ROOM_TEMP = T_tri + 25 * si.K
-M = si.mole / si.litre
-H_u = M / p_STP
-dT_u = si.K
+
+T0 = sci.zero_Celsius * si.kelvin
# there are so few water ions instead of K we have K [H2O] (see Seinfeld & Pandis p 345)
+M = si.mole / si.litre
K_H2O = 1e-14 * M * M
-
-default_random_seed = (
- 44 if 'CI' in os.environ # https://en.wikipedia.org/wiki/44_(number)
- else time.time_ns()
-)
diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py
new file mode 100644
index 000000000..148c66494
--- /dev/null
+++ b/PySDM/physics/constants_defaults.py
@@ -0,0 +1,97 @@
+"""
+default values for constants which can be altered by providing alternative
+values in a constants dictionary passed to Formulae __init__ method
+"""
+import numpy as np
+from scipy import constants as sci
+from chempy import Substance
+from .constants import si, M, PI_4_3, ONE_THIRD, THREE, T0, TWO_THIRDS, PI # pylint: disable=unused-import
+
+Md = (
+ 0.78 * Substance.from_formula('N2').mass * si.gram / si.mole +
+ 0.21 * Substance.from_formula('O2').mass * si.gram / si.mole +
+ 0.01 * Substance.from_formula('Ar').mass * si.gram / si.mole
+)
+Mv = Substance.from_formula('H2O').mass * si.gram / si.mole
+
+R_str = sci.R * si.joule / si.kelvin / si.mole
+eps = Mv / Md
+Rd = R_str / Md
+Rv = R_str / Mv
+
+D0 = 2.26e-5 * si.metre ** 2 / si.second
+D_exp = 1.81
+
+K0 = 2.4e-2 * si.joules / si.metres / si.seconds / si.kelvins
+
+p1000 = 1000 * si.hectopascals
+c_pd = 1005 * si.joule / si.kilogram / si.kelvin
+c_pv = 1850 * si.joule / si.kilogram / si.kelvin
+g_std = sci.g * si.metre / si.second ** 2
+
+c_pw = 4218 * si.joule / si.kilogram / si.kelvin
+
+Rd_over_c_pd = Rd / c_pd
+
+ARM_C1 = 6.1094 * si.hectopascal
+ARM_C2 = 17.625 * si.dimensionless
+ARM_C3 = 243.04 * si.kelvin
+
+FWC_C0 = 6.115836990e000 * si.hPa
+FWC_C1 = 0.444606896e000 * si.hPa / si.K
+FWC_C2 = 0.143177157e-01 * si.hPa / si.K**2
+FWC_C3 = 0.264224321e-03 * si.hPa / si.K**3
+FWC_C4 = 0.299291081e-05 * si.hPa / si.K**4
+FWC_C5 = 0.203154182e-07 * si.hPa / si.K**5
+FWC_C6 = 0.702620698e-10 * si.hPa / si.K**6
+FWC_C7 = 0.379534310e-13 * si.hPa / si.K**7
+FWC_C8 = -.321582393e-15 * si.hPa / si.K**8
+
+FWC_I0 = 6.098689930e000 * si.hPa
+FWC_I1 = 0.499320233e000 * si.hPa / si.K
+FWC_I2 = 0.184672631e-01 * si.hPa / si.K**2
+FWC_I3 = 0.402737184e-03 * si.hPa / si.K**3
+FWC_I4 = 0.565392987e-05 * si.hPa / si.K**4
+FWC_I5 = 0.521693933e-07 * si.hPa / si.K**5
+FWC_I6 = 0.307839583e-09 * si.hPa / si.K**6
+FWC_I7 = 0.105785160e-11 * si.hPa / si.K**7
+FWC_I8 = 0.161444444e-14 * si.hPa / si.K**8
+
+rho_w = 1 * si.kilograms / si.litres
+rho_i = 916.8 * si.kg / si.metres**3
+pH_w = 7
+sgm_w = 0.072 * si.joule / si.metre ** 2
+nu_w = Mv / rho_w
+
+p_tri = 611.73 * si.pascal
+T_tri = 273.16 * si.kelvin
+l_tri = 2.5e6 * si.joule / si.kilogram
+
+# standard pressure and temperature (ICAO)
+T_STP = (sci.zero_Celsius + 15) * si.kelvin
+p_STP = 101325 * si.pascal
+rho_STP = p_STP / Rd / T_STP
+
+ROOM_TEMP = T_tri + 25 * si.K
+H_u = M / p_STP
+dT_u = si.K
+
+sgm_org = np.nan
+delta_min = np.nan
+
+RUEHL_nu_org = np.nan
+RUEHL_A0 = np.nan
+RUEHL_C0 = np.nan
+RUEHL_m_sigma = np.nan
+RUEHL_sgm_min = np.nan
+
+BIGG_DT_MEDIAN = np.nan
+
+NIEMAND_A = np.nan
+NIEMAND_B = np.nan
+
+ABIFM_M = np.inf
+ABIFM_C = np.inf
+ABIFM_UNIT = 1 / si.cm**2 / si.s
+
+J_HET = np.nan
diff --git a/PySDM/physics/diffusion_kinetics/fuchs_sutugin.py b/PySDM/physics/diffusion_kinetics/fuchs_sutugin.py
index 577518d06..7c657e7aa 100644
--- a/PySDM/physics/diffusion_kinetics/fuchs_sutugin.py
+++ b/PySDM/physics/diffusion_kinetics/fuchs_sutugin.py
@@ -2,19 +2,21 @@
Fuch-Sutugin transition-regime correction as advocated for cloud modelling
in [Laaksonen et al. 2005](https://doi.org/10.5194/acp-5-461-2005)
"""
-from numpy import sqrt
-from PySDM.physics import constants as const
+import numpy as np
class FuchsSutugin:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def lambdaD(D, T):
- return D / sqrt(2 * const.Rv * T)
+ def lambdaD(const, D, T):
+ return D / np.sqrt(2 * const.Rv * T)
@staticmethod
- def lambdaK(T, p):
- return (4. / 5) * const.K0 * T / p / sqrt(2 * const.Rd * T)
+ def lambdaK(const, T, p):
+ return (4. / 5) * const.K0 * T / p / np.sqrt(2 * const.Rd * T)
@staticmethod
- def DK(DK, r, lmbd):
+ def DK(_, DK, r, lmbd):
return DK * (1 + lmbd/r) / (1 + 1.71 * lmbd/r + 1.33 * lmbd/r * lmbd/r)
diff --git a/PySDM/physics/diffusion_kinetics/neglect.py b/PySDM/physics/diffusion_kinetics/neglect.py
index b954fba8c..48187b046 100644
--- a/PySDM/physics/diffusion_kinetics/neglect.py
+++ b/PySDM/physics/diffusion_kinetics/neglect.py
@@ -4,14 +4,17 @@
class Neglect:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def lambdaD(D, T): # pylint: disable=unused-argument
+ def lambdaD(_, D, T): # pylint: disable=unused-argument
return -1
@staticmethod
- def lambdaK(T, p): # pylint: disable=unused-argument
+ def lambdaK(_, T, p): # pylint: disable=unused-argument
return -1
@staticmethod
- def DK(DK, r, lmbd): # pylint: disable=unused-argument
+ def DK(_, DK, r, lmbd): # pylint: disable=unused-argument
return DK
diff --git a/PySDM/physics/diffusion_thermics/negclect.py b/PySDM/physics/diffusion_thermics/negclect.py
index 1fe19717e..c3c5cfad0 100644
--- a/PySDM/physics/diffusion_thermics/negclect.py
+++ b/PySDM/physics/diffusion_thermics/negclect.py
@@ -1,10 +1,12 @@
"""
constant diffusion coefficient formulation
"""
-from PySDM.physics import constants as const
class Neglect:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def D(T, p): # pylint: disable=unused-argument
+ def D(const, T, p): # pylint: disable=unused-argument
return const.D0
diff --git a/PySDM/physics/diffusion_thermics/tracy_welch_porter.py b/PySDM/physics/diffusion_thermics/tracy_welch_porter.py
index 5e5572fcd..3c382f140 100644
--- a/PySDM/physics/diffusion_thermics/tracy_welch_porter.py
+++ b/PySDM/physics/diffusion_thermics/tracy_welch_porter.py
@@ -2,11 +2,13 @@
based on "PROPERTIES OF AIR: A Manual for Use in Biophysical Ecology"
(Fourth Edition - 2010, page 22)
"""
-from numpy import power
-from PySDM.physics import constants as const
+import numpy as np
class TracyWelchPorter:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def D(T, p):
- return const.D0 * power(T / const.T0, const.D_exp) * (const.p1000 / p)
+ def D(const, T, p):
+ return const.D0 * np.power(T / const.T0, const.D_exp) * (const.p1000 / p)
diff --git a/PySDM/physics/dimensional_analysis.py b/PySDM/physics/dimensional_analysis.py
index b28894da0..54847c85d 100644
--- a/PySDM/physics/dimensional_analysis.py
+++ b/PySDM/physics/dimensional_analysis.py
@@ -4,7 +4,7 @@
"""
from importlib import reload
from PySDM import formulae
-from . import constants
+from . import constants, constants_defaults
from .impl import flag
@@ -13,9 +13,11 @@ class DimensionalAnalysis:
def __enter__(*_): # pylint: disable=no-method-argument
flag.DIMENSIONAL_ANALYSIS = True
reload(constants)
+ reload(constants_defaults)
reload(formulae)
def __exit__(*_): # pylint: disable=no-method-argument
flag.DIMENSIONAL_ANALYSIS = False
reload(constants)
+ reload(constants_defaults)
reload(formulae)
diff --git a/PySDM/physics/drop_growth/maxwell_mason.py b/PySDM/physics/drop_growth/maxwell_mason.py
index 6bfe53c17..1197592b5 100644
--- a/PySDM/physics/drop_growth/maxwell_mason.py
+++ b/PySDM/physics/drop_growth/maxwell_mason.py
@@ -1,12 +1,14 @@
"""
single-equation approximation of the vapour and heat diffusion problem
"""
-import PySDM.physics.constants as const
class MaxwellMason:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def r_dr_dt(RH_eq, T, RH, lv, pvs, D, K):
+ def r_dr_dt(const, RH_eq, T, RH, lv, pvs, D, K):
return (RH - RH_eq) / (
const.rho_w * const.Rv * T / D / pvs +
const.rho_w * lv / K / T * (lv / const.Rv / T - 1)
diff --git a/PySDM/physics/freezing_temperature_spectrum/bigg_1953.py b/PySDM/physics/freezing_temperature_spectrum/bigg_1953.py
index 4283a5dad..2c89d2b03 100644
--- a/PySDM/physics/freezing_temperature_spectrum/bigg_1953.py
+++ b/PySDM/physics/freezing_temperature_spectrum/bigg_1953.py
@@ -3,26 +3,22 @@
formulae (i.e. immersed surface independent)
"""
import numpy as np
-from PySDM.physics import constants as const
-
-P_median = .5
-DT_median = np.nan
class Bigg_1953:
- def __init__(self):
- assert np.isfinite(DT_median)
+ def __init__(self, const):
+ assert np.isfinite(const.BIGG_DT_MEDIAN)
@staticmethod
- def pdf(T, A_insol): # pylint: disable=unused-argument
- A = np.log(1 - P_median)
- B = DT_median - const.T0
+ def pdf(const, T, A_insol): # pylint: disable=unused-argument
+ A = np.log(1 - .5)
+ B = const.BIGG_DT_MEDIAN - const.T0
return - A * np.exp(A * np.exp(B + T) + B + T)
@staticmethod
- def cdf(T, A_insol): # pylint: disable=unused-argument
- return np.exp(np.log(1 - P_median) * np.exp(DT_median - (const.T0 - T)))
+ def cdf(const, T, A_insol): # pylint: disable=unused-argument
+ return np.exp(np.log(1 - .5) * np.exp(const.BIGG_DT_MEDIAN - (const.T0 - T)))
@staticmethod
- def median():
- return const.T0 - DT_median
+ def median(const):
+ return const.T0 - const.BIGG_DT_median
diff --git a/PySDM/physics/freezing_temperature_spectrum/niemand_et_al_2012.py b/PySDM/physics/freezing_temperature_spectrum/niemand_et_al_2012.py
index fdc68843c..26042b31a 100644
--- a/PySDM/physics/freezing_temperature_spectrum/niemand_et_al_2012.py
+++ b/PySDM/physics/freezing_temperature_spectrum/niemand_et_al_2012.py
@@ -3,31 +3,35 @@
[Niemand et al. 2012](https://doi.org/10.1175/JAS-D-11-0249.1) INAS density parameterization
"""
import numpy as np
-from PySDM.physics import constants as const
-
-a = np.nan
-b = np.nan
class Niemand_et_al_2012:
def __str__(self):
return 'Niemand et al. 2012'
- def __init__(self):
- assert np.isfinite(a)
- assert np.isfinite(b)
+ def __init__(self, const):
+ assert np.isfinite(const.NIEMAND_A)
+ assert np.isfinite(const.NIEMAND_B)
@staticmethod
- def pdf(T, A_insol):
- ns_T = np.exp(a * (T - const.T0) + b)
- return -A_insol * a * ns_T * np.exp(-A_insol * ns_T)
+ def pdf(const, T, A_insol):
+ ns_T = np.exp(const.NIEMAND_A * (T - const.T0) + const.NIEMAND_B)
+ return -A_insol * const.NIEMAND_A * ns_T * np.exp(-A_insol * ns_T)
@staticmethod
- def cdf(T, A_insol):
- ns_T = np.exp(a * (T - const.T0) + b)
- return 1 - np.exp(-A_insol * ns_T) - np.exp(-A_insol*np.exp(-a * const.T0 + b))
+ def cdf(const, T, A_insol):
+ ns_T = np.exp(const.NIEMAND_A * (T - const.T0) + const.NIEMAND_B)
+ return 1 - np.exp(-A_insol * ns_T) - np.exp(
+ -A_insol*np.exp(-const.NIEMAND_A * const.T0 + const.NIEMAND_B)
+ )
@staticmethod
- def invcdf(cdf, A_insol):
- tmp = np.log((np.log(1 - cdf) + np.exp(-A_insol*np.exp(-a * const.T0 + b))) / -A_insol)
- return const.T0 + (tmp - b) / a
+ def invcdf(const, cdf, A_insol):
+ tmp = np.log(
+ (
+ np.log(1 - cdf)
+ +
+ np.exp(-A_insol*np.exp(-const.NIEMAND_A * const.T0 + const.NIEMAND_B))
+ ) / -A_insol
+ )
+ return const.T0 + (tmp - const.NIEMAND_B) / const.NIEMAND_A
diff --git a/PySDM/physics/freezing_temperature_spectrum/null.py b/PySDM/physics/freezing_temperature_spectrum/null.py
index cce2051f8..8082b60e4 100644
--- a/PySDM/physics/freezing_temperature_spectrum/null.py
+++ b/PySDM/physics/freezing_temperature_spectrum/null.py
@@ -5,6 +5,9 @@
class Null:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def cdf(T):
+ def cdf(const, T):
pass
diff --git a/PySDM/physics/heterogeneous_ice_nucleation_rate/abifm.py b/PySDM/physics/heterogeneous_ice_nucleation_rate/abifm.py
index 576a8119f..458a9567c 100644
--- a/PySDM/physics/heterogeneous_ice_nucleation_rate/abifm.py
+++ b/PySDM/physics/heterogeneous_ice_nucleation_rate/abifm.py
@@ -3,19 +3,13 @@
([Knopf & Alpert 2013](https://doi.org/10.1039/C3FD00035D))
"""
import numpy as np
-from ..constants import si
-
-m = np.inf
-c = np.inf
-unit = 1 / si.cm**2 / si.s
class ABIFM:
- @staticmethod
- def _check():
- assert np.isfinite(m)
- assert np.isfinite(c)
+ def __init__(self, const):
+ assert np.isfinite(const.ABIFM_M)
+ assert np.isfinite(const.ABIFM_C)
@staticmethod
- def j_het(a_w_ice):
- return 10**(m * (1 - a_w_ice) + c) * unit
+ def j_het(const, a_w_ice):
+ return 10**(const.ABIFM_M * (1 - a_w_ice) + const.ABIFM_C) * const.ABIFM_UNIT
diff --git a/PySDM/physics/heterogeneous_ice_nucleation_rate/constant.py b/PySDM/physics/heterogeneous_ice_nucleation_rate/constant.py
index ccd082d70..0b1d6d21b 100644
--- a/PySDM/physics/heterogeneous_ice_nucleation_rate/constant.py
+++ b/PySDM/physics/heterogeneous_ice_nucleation_rate/constant.py
@@ -3,14 +3,11 @@
"""
import numpy as np
-J_HET = np.nan
-
class Constant:
- @staticmethod
- def _check():
- assert np.isfinite(J_HET)
+ def __init__(self, const):
+ assert np.isfinite(const.J_HET)
@staticmethod
- def j_het(a_w_ice): # pylint: disable=unused-argument
- return J_HET
+ def j_het(const, a_w_ice): # pylint: disable=unused-argument
+ return const.J_HET
diff --git a/PySDM/physics/heterogeneous_ice_nucleation_rate/null.py b/PySDM/physics/heterogeneous_ice_nucleation_rate/null.py
index 3085c2bad..f94abb576 100644
--- a/PySDM/physics/heterogeneous_ice_nucleation_rate/null.py
+++ b/PySDM/physics/heterogeneous_ice_nucleation_rate/null.py
@@ -5,6 +5,9 @@
class Null:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def j_het(a_w_ice):
+ def j_het(const, a_w_ice):
pass
diff --git a/PySDM/physics/hydrostatics/default.py b/PySDM/physics/hydrostatics/default.py
index ab0d0eb0e..6a6a7b131 100644
--- a/PySDM/physics/hydrostatics/default.py
+++ b/PySDM/physics/hydrostatics/default.py
@@ -1,24 +1,26 @@
"""
# TODO #407
"""
-from numpy import power
-from PySDM.physics import constants as const
+import numpy as np
class Default:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def drho_dz(g, p, T, qv, lv, dql_dz=0):
+ def drho_dz(const, g, p, T, qv, lv, dql_dz=0):
Rq = const.Rv / (1 / qv + 1) + const.Rd / (1 + qv)
cp = const.c_pv / (1 / qv + 1) + const.c_pd / (1 + qv)
rho = p / Rq / T
return (g / T * rho * (Rq / cp - 1) - p * lv / cp / T**2 * dql_dz) / Rq
@staticmethod
- def p_of_z_assuming_const_th_and_qv(g, p0, thstd, qv, z):
+ def p_of_z_assuming_const_th_and_qv(const, g, p0, thstd, qv, z):
z0 = 0
Rq = const.Rv / (1 / qv + 1) + const.Rd / (1 + qv)
- arg = power(
+ arg = np.power(
p0/const.p1000,
const.Rd_over_c_pd
) - (z-z0) * const.Rd_over_c_pd * g / thstd / Rq
- return const.p1000 * power(arg, 1/const.Rd_over_c_pd)
+ return const.p1000 * np.power(arg, 1/const.Rd_over_c_pd)
diff --git a/PySDM/physics/hygroscopicity/kappa_koehler.py b/PySDM/physics/hygroscopicity/kappa_koehler.py
index fdfd6a035..bc0b76f4d 100644
--- a/PySDM/physics/hygroscopicity/kappa_koehler.py
+++ b/PySDM/physics/hygroscopicity/kappa_koehler.py
@@ -2,18 +2,20 @@
kappa-Koehler parameterization
([Petters & Kreidenweis 2007](https://doi.org/10.5194/acp-7-1961-2007))
"""
-from numpy import sqrt, exp
-import PySDM.physics.constants as const
+import numpy as np
class KappaKoehler:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def RH_eq(r, T, kp, rd3, sgm):
- return exp(
+ def RH_eq(const, r, T, kp, rd3, sgm):
+ return np.exp(
(2 * sgm / const.Rv / T / const.rho_w) / r
) * (r**3 - rd3) / (r**3 - rd3 * (1-kp))
@staticmethod
- def r_cr(kp, rd3, T, sgm):
+ def r_cr(const, kp, rd3, T, sgm):
# TODO #493
- return sqrt(3 * kp * rd3 / (2 * sgm / const.Rv / T / const.rho_w))
+ return np.sqrt(3 * kp * rd3 / (2 * sgm / const.Rv / T / const.rho_w))
diff --git a/PySDM/physics/hygroscopicity/kappa_koehler_leading_terms.py b/PySDM/physics/hygroscopicity/kappa_koehler_leading_terms.py
index 71cd0c1c4..ed7fb10d8 100644
--- a/PySDM/physics/hygroscopicity/kappa_koehler_leading_terms.py
+++ b/PySDM/physics/hygroscopicity/kappa_koehler_leading_terms.py
@@ -3,19 +3,21 @@
two-term formulation with 1/r and 1/r^3 terms corresponding to surface-tension
(Kelvin) and soluble substance (Raoult/Wüllner) effects, respectively
"""
-from numpy import sqrt, power
-import PySDM.physics.constants as const
+import numpy as np
class KappaKoehlerLeadingTerms:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def RH_eq(r, T, kp, rd3, sgm):
+ def RH_eq(const, r, T, kp, rd3, sgm):
return (
1 +
(2 * sgm / const.Rv / T / const.rho_w) / r -
- kp * rd3 / power(r, const.THREE)
+ kp * rd3 / np.power(r, const.THREE)
)
@staticmethod
- def r_cr(kp, rd3, T, sgm):
- return sqrt(3 * kp * rd3 / (2 * sgm / const.Rv / T / const.rho_w))
+ def r_cr(const, kp, rd3, T, sgm):
+ return np.sqrt(3 * kp * rd3 / (2 * sgm / const.Rv / T / const.rho_w))
diff --git a/PySDM/physics/impl/__init__.py b/PySDM/physics/impl/__init__.py
index beda360da..ff470b35c 100644
--- a/PySDM/physics/impl/__init__.py
+++ b/PySDM/physics/impl/__init__.py
@@ -2,3 +2,4 @@
physics stuff not intended to be imported from user code
(incl. the `PySDM.physics.impl.fake_unit_registry.FakeUnitRegistry`)
"""
+from . import flag
diff --git a/PySDM/physics/latent_heat/constant.py b/PySDM/physics/latent_heat/constant.py
index fcafe3c07..c70c91de2 100644
--- a/PySDM/physics/latent_heat/constant.py
+++ b/PySDM/physics/latent_heat/constant.py
@@ -1,10 +1,12 @@
"""
temperature-independent latent heat of vaporization
"""
-from PySDM.physics import constants as const
class Constant:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def lv(T): # pylint: disable=unused-argument
+ def lv(const, T): # pylint: disable=unused-argument
return const.l_tri
diff --git a/PySDM/physics/latent_heat/kirchhoff.py b/PySDM/physics/latent_heat/kirchhoff.py
index 7a7cee530..13929d53d 100644
--- a/PySDM/physics/latent_heat/kirchhoff.py
+++ b/PySDM/physics/latent_heat/kirchhoff.py
@@ -3,10 +3,12 @@
law](https://en.wikipedia.org/wiki/Gustav_Kirchhoff#Kirchhoff's_law_of_thermochemistry)
based temperature-dependent latent heat of vaporization
"""
-from PySDM.physics import constants as const
class Kirchhoff:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def lv(T):
+ def lv(const, T):
return const.l_tri + (const.c_pv - const.c_pw) * (T - const.T_tri)
diff --git a/PySDM/physics/particle_advection/explicit_in_space.py b/PySDM/physics/particle_advection/explicit_in_space.py
index 08f0e67c3..108e18744 100644
--- a/PySDM/physics/particle_advection/explicit_in_space.py
+++ b/PySDM/physics/particle_advection/explicit_in_space.py
@@ -4,6 +4,9 @@
class ExplicitInSpace:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def displacement(omega, c_l, c_r):
+ def displacement(_, omega, c_l, c_r):
return c_l * (1 - omega) + c_r * omega
diff --git a/PySDM/physics/particle_advection/implicit_in_space.py b/PySDM/physics/particle_advection/implicit_in_space.py
index a49bdf622..8a025b200 100644
--- a/PySDM/physics/particle_advection/implicit_in_space.py
+++ b/PySDM/physics/particle_advection/implicit_in_space.py
@@ -4,6 +4,9 @@
class ImplicitInSpace:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def displacement(omega, c_l, c_r):
+ def displacement(_, omega, c_l, c_r):
return (omega * c_r + c_l * (1 - omega)) / (1 - c_r + c_l)
diff --git a/PySDM/physics/saturation_vapour_pressure/august_roche_magnus.py b/PySDM/physics/saturation_vapour_pressure/august_roche_magnus.py
index 7fe16dd32..9677287bb 100644
--- a/PySDM/physics/saturation_vapour_pressure/august_roche_magnus.py
+++ b/PySDM/physics/saturation_vapour_pressure/august_roche_magnus.py
@@ -3,16 +3,17 @@
[Wikipedia](https://en.wikipedia.org/wiki/Clausius–Clapeyron_relation#August–Roche–Magnus_formula)
and references therein)
"""
-
-from numpy import exp, nan
-from PySDM.physics import constants as const
+import numpy as np
class AugustRocheMagnus:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def pvs_Celsius(T):
- return const.ARM_C1 * exp((const.ARM_C2 * T) / (T + const.ARM_C3))
+ def pvs_Celsius(const, T):
+ return const.ARM_C1 * np.exp((const.ARM_C2 * T) / (T + const.ARM_C3))
@staticmethod
- def ice_Celsius(T):
- return nan * T
+ def ice_Celsius(_, T):
+ return np.nan * T
diff --git a/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py b/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py
index 2a2db419f..f65ce8ecf 100644
--- a/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py
+++ b/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py
@@ -2,12 +2,14 @@
polynomial fits from
[Flatau et al. 1992](https://doi.org/10.1175/1520-0450(1992)031%3C1507:PFTSVP%3E2.0.CO;2)
"""
-from PySDM.physics import constants as const
class FlatauWalkoCotton:
+ def __init__(self, _):
+ pass
+
@staticmethod
- def pvs_Celsius(T):
+ def pvs_Celsius(const, T):
return (
const.FWC_C0 + T * (
const.FWC_C1 + T * (
@@ -21,7 +23,7 @@ def pvs_Celsius(T):
))))))))
@staticmethod
- def ice_Celsius(T):
+ def ice_Celsius(const, T):
return (
const.FWC_I0 + T * (
const.FWC_I1 + T * (
@@ -33,11 +35,3 @@ def ice_Celsius(T):
const.FWC_I7 + T *
const.FWC_I8
))))))))
-
- @staticmethod
- def a_w_ice(T):
- return (
- FlatauWalkoCotton.ice_Celsius(T - const.T0)
- /
- FlatauWalkoCotton.pvs_Celsius(T - const.T0)
- )
diff --git a/PySDM/physics/state_variable_triplet/rhod_thd_qv.py b/PySDM/physics/state_variable_triplet/rhod_thd_qv.py
index 90c18d22f..b74718b4d 100644
--- a/PySDM/physics/state_variable_triplet/rhod_thd_qv.py
+++ b/PySDM/physics/state_variable_triplet/rhod_thd_qv.py
@@ -2,41 +2,43 @@
dry-air density / dry-air potential temperature / water vapour mixing ratio triplet
(as in libcloudph++)
"""
-from numpy import power
-from PySDM.physics import constants as const
+import numpy as np
class RhodThdQv:
+ def __init__(self, _):
+ pass
+
# A14 in libcloudph++ 1.0 paper
@staticmethod
- def T(rhod, thd):
- return thd * power(
+ def T(const, rhod, thd):
+ return thd * np.power(
rhod * thd / const.p1000 * const.Rd, const.Rd_over_c_pd / (1 - const.Rd_over_c_pd)
)
# A15 in libcloudph++ 1.0 paper
@staticmethod
- def p(rhod, T, qv):
+ def p(const, rhod, T, qv):
return rhod * (1 + qv) * (const.Rv / (1 / qv + 1) + const.Rd / (1 + qv)) * T
@staticmethod
- def pv(p, qv):
+ def pv(const, p, qv):
return p * qv / (qv + const.eps)
@staticmethod
- def dthd_dt(rhod, thd, T, dqv_dt, lv):
+ def dthd_dt(const, rhod, thd, T, dqv_dt, lv):
return - lv * dqv_dt / const.c_pd / T * thd * rhod
@staticmethod
- def th_dry(th_std, qv):
- return th_std * power(1 + qv / const.eps, const.Rd / const.c_pd)
+ def th_dry(const, th_std, qv):
+ return th_std * np.power(1 + qv / const.eps, const.Rd / const.c_pd)
@staticmethod
- def rho_d(p, qv, theta_std):
+ def rho_d(const, p, qv, theta_std):
return p * (
1 - 1 / (1 + const.eps / qv)
) / (
- power(p / const.p1000, const.Rd_over_c_pd) * const.Rd * theta_std
+ np.power(p / const.p1000, const.Rd_over_c_pd) * const.Rd * theta_std
)
@staticmethod
@@ -44,5 +46,5 @@ def rho_of_rhod_qv(rhod, qv):
return rhod * (1 + qv)
@staticmethod
- def rhod_of_pd_T(pd, T):
+ def rhod_of_pd_T(const, pd, T):
return pd / const.Rd / T
diff --git a/PySDM/physics/surface_tension/compressed_film_ovadnevaite.py b/PySDM/physics/surface_tension/compressed_film_ovadnevaite.py
index c38e7efb6..25a5915df 100644
--- a/PySDM/physics/surface_tension/compressed_film_ovadnevaite.py
+++ b/PySDM/physics/surface_tension/compressed_film_ovadnevaite.py
@@ -3,10 +3,6 @@
as in [Ovadnevaite et al. (2017)](https://doi.org/10.1038/nature22806)
"""
import numpy as np
-from PySDM.physics import constants as const
-
-sgm_org = np.nan
-delta_min = np.nan
class CompressedFilmOvadnevaite:
@@ -19,18 +15,17 @@ class CompressedFilmOvadnevaite:
partitions to the surface of the droplet and that the surface tension is a weighted
average of the surface tension of water and the organic material.
"""
- @staticmethod
- def _check():
- assert np.isfinite(sgm_org)
- assert np.isfinite(delta_min)
+ def __init__(self, const):
+ assert np.isfinite(const.sgm_org)
+ assert np.isfinite(const.delta_min)
@staticmethod
- def sigma(T, v_wet, v_dry, f_org): # pylint: disable=unused-argument
+ def sigma(const, T, v_wet, v_dry, f_org): # pylint: disable=unused-argument
# convert wet volume to wet radius
r_wet = ((3 * v_wet) / (4 * np.pi)) ** (1 / 3)
# calculate the minimum shell volume, v_delta
- v_delta = v_wet - ((4 * np.pi) / 3 * (r_wet - delta_min) ** 3)
+ v_delta = v_wet - ((4 * np.pi) / 3 * (r_wet - const.delta_min) ** 3)
# calculate the total volume of organic, v_beta
v_beta = f_org * v_dry
@@ -39,5 +34,5 @@ def sigma(T, v_wet, v_dry, f_org): # pylint: disable=unused-argument
c_beta = np.minimum(v_beta / v_delta, 1)
# calculate sigma
- sgm = (1-c_beta) * const.sgm_w + c_beta * sgm_org
+ sgm = (1-c_beta) * const.sgm_w + c_beta * const.sgm_org
return sgm
diff --git a/PySDM/physics/surface_tension/compressed_film_ruehl.py b/PySDM/physics/surface_tension/compressed_film_ruehl.py
index b4788ea39..55cec9273 100644
--- a/PySDM/physics/surface_tension/compressed_film_ruehl.py
+++ b/PySDM/physics/surface_tension/compressed_film_ruehl.py
@@ -5,13 +5,6 @@
import numpy as np
from scipy import constants as sci
from scipy import optimize
-from PySDM.physics import constants as const
-
-nu_org = np.nan
-A0 = np.nan
-C0 = np.nan
-m_sigma = np.nan
-sgm_min = np.nan
class CompressedFilmRuehl:
@@ -26,32 +19,31 @@ class CompressedFilmRuehl:
the surface concentration to the surface tension. For the compressed film model it
is linear, with slope `m_sigma`.
"""
- @staticmethod
- def _check():
- assert np.isfinite(nu_org)
- assert np.isfinite(A0)
- assert np.isfinite(C0)
- assert np.isfinite(m_sigma)
- assert np.isfinite(sgm_min)
+ def __init__(self, const):
+ assert np.isfinite(const.RUEHL_nu_org)
+ assert np.isfinite(const.RUEHL_A0)
+ assert np.isfinite(const.RUEHL_C0)
+ assert np.isfinite(const.RUEHL_m_sigma)
+ assert np.isfinite(const.RUEHL_sgm_min)
@staticmethod
- def sigma(T, v_wet, v_dry, f_org):
+ def sigma(const, T, v_wet, v_dry, f_org):
r_wet = ((3 * v_wet) / (4 * np.pi))**(1/3) # m - wet radius
# C_bulk is the concentration of the organic in the bulk phase
- Cb_iso = (f_org*v_dry/nu_org) / (v_wet/const.nu_w) # = C_bulk / (1-f_surf)
+ Cb_iso = (f_org*v_dry/const.RUEHL_nu_org) / (v_wet/const.nu_w) # = C_bulk / (1-f_surf)
- # A is the area that one molecule of organic occupies at the droplet surface
- A_iso = (4 * np.pi * r_wet**2) / (f_org * v_dry * sci.N_A / nu_org) # m^2 = A*f_surf
+ # A is the area one molecule of organic occupies at the droplet surface (m^2 = A*f_surf)
+ A_iso = (4 * np.pi * r_wet**2) / (f_org * v_dry * sci.N_A / const.RUEHL_nu_org)
# solve implicitly for fraction of organic at surface
- f = lambda f_surf: Cb_iso*(1-f_surf) - C0*np.exp(
- ((A0**2 - (A_iso/f_surf)**2)*m_sigma*sci.N_A)/(2*sci.R*T)
+ f = lambda f_surf: Cb_iso*(1-f_surf) - const.RUEHL_C0*np.exp(
+ ((const.RUEHL_A0**2 - (A_iso/f_surf)**2)*const.RUEHL_m_sigma*sci.N_A)/(2*sci.R*T)
)
- sol = optimize.root_scalar(f, bracket=[0,1])
+ sol = optimize.root_scalar(f, bracket=[0, 1])
f_surf = sol.root
- # calculate surface tension
- sgm = const.sgm_w - (A0 - A_iso/f_surf)*m_sigma # m^2 * J/m^2 = J = N*m --> N/m - N*m ?
- sgm = np.minimum(np.maximum(sgm, sgm_min), const.sgm_w)
+ # calculate surface tension (m^2 * J/m^2 = J = N*m --> N/m - N*m ?)
+ sgm = const.sgm_w - (const.RUEHL_A0 - A_iso/f_surf)*const.RUEHL_m_sigma
+ sgm = np.minimum(np.maximum(sgm, const.RUEHL_sgm_min), const.sgm_w)
return sgm
diff --git a/PySDM/physics/surface_tension/constant.py b/PySDM/physics/surface_tension/constant.py
index 3412b7eba..c258a2b6c 100644
--- a/PySDM/physics/surface_tension/constant.py
+++ b/PySDM/physics/surface_tension/constant.py
@@ -1,7 +1,6 @@
"""
constant surface tension coefficient
"""
-import PySDM.physics.constants as const
class Constant:
@@ -10,6 +9,10 @@ class Constant:
droplet surface is composed of pure water with constant surface
tension `sgm_w`.
"""
+
+ def __init__(self, _):
+ pass
+
@staticmethod
- def sigma(T, v_wet, v_dry, f_org): # pylint: disable=unused-argument
+ def sigma(const, T, v_wet, v_dry, f_org): # pylint: disable=unused-argument
return const.sgm_w
diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py
index 04cd86e89..f423b7a29 100644
--- a/PySDM/physics/trivia.py
+++ b/PySDM/physics/trivia.py
@@ -1,25 +1,27 @@
"""
Various (hopefully) undebatable formulae
"""
-from numpy import abs, log10, exp, power # pylint: disable=redefined-builtin
-from . import constants as const
+import numpy as np
class Trivia:
+ def __init__(self, _):
+ pass
+
@staticmethod
def volume_of_density_mass(rho, m):
return m / rho
@staticmethod
- def radius(volume):
- return power(volume / const.PI_4_3, const.ONE_THIRD)
+ def radius(const, volume):
+ return np.power(volume / const.PI_4_3, const.ONE_THIRD)
@staticmethod
- def volume(radius):
- return const.PI_4_3 * power(radius, const.THREE)
+ def volume(const, radius):
+ return const.PI_4_3 * np.power(radius, const.THREE)
@staticmethod
- def sphere_surface(diameter):
+ def sphere_surface(const, diameter):
return const.PI * diameter**2
@staticmethod
@@ -28,27 +30,27 @@ def explicit_euler(y, dt, dy_dt):
@staticmethod
def within_tolerance(error_estimate, value, rtol):
- return error_estimate < rtol * abs(value)
+ return error_estimate < rtol * np.abs(value)
@staticmethod
def H2pH(H):
- return -log10(H * 1e-3)
+ return -np.log10(H * 1e-3)
@staticmethod
def pH2H(pH):
- return power(10, -pH) * 1e3
+ return np.power(10, -pH) * 1e3
@staticmethod
- def vant_hoff(K, dH, T, *, T_0):
- return K * exp(-dH / const.R_str * (1 / T - 1 / T_0))
+ def vant_hoff(const, K, dH, T, *, T_0):
+ return K * np.exp(-dH / const.R_str * (1 / T - 1 / T_0))
@staticmethod
- def tdep2enthalpy(tdep):
+ def tdep2enthalpy(const, tdep):
return -tdep * const.R_str
@staticmethod
- def arrhenius(A, Ea, T):
- return A * exp(-Ea / (const.R_str * T))
+ def arrhenius(const, A, Ea, T):
+ return A * np.exp(-Ea / (const.R_str * T))
@staticmethod
def mole_fraction_2_mixing_ratio(mole_fraction, specific_gravity):
@@ -59,9 +61,9 @@ def mixing_ratio_2_mole_fraction(mixing_ratio, specific_gravity):
return mixing_ratio / (specific_gravity + mixing_ratio)
@staticmethod
- def p_d(p, qv):
+ def p_d(const, p, qv):
return p * (1 - 1 / (1 + const.eps / qv))
@staticmethod
- def th_std(p, T):
- return T * power(const.p1000 / p, const.Rd_over_c_pd)
+ def th_std(const, p, T):
+ return T * np.power(const.p1000 / p, const.Rd_over_c_pd)
diff --git a/PySDM/physics/ventilation/neglect.py b/PySDM/physics/ventilation/neglect.py
index e0e5b19e0..2e779e4a1 100644
--- a/PySDM/physics/ventilation/neglect.py
+++ b/PySDM/physics/ventilation/neglect.py
@@ -4,4 +4,5 @@
class Neglect:
- pass
+ def __init__(self, _):
+ pass
diff --git a/PySDM/products/aqueous_chemistry/aqueous_mole_fraction.py b/PySDM/products/aqueous_chemistry/aqueous_mole_fraction.py
index b2d54e18e..24c40e5b9 100644
--- a/PySDM/products/aqueous_chemistry/aqueous_mole_fraction.py
+++ b/PySDM/products/aqueous_chemistry/aqueous_mole_fraction.py
@@ -2,7 +2,6 @@
mole fractions of aqueous-chemistry relevant compounds
"""
from PySDM.products.impl.moment_product import MomentProduct
-from PySDM.physics.constants import Md
DUMMY_SPECIFIC_GRAVITY = 44
@@ -26,7 +25,7 @@ def _impl(self, **kwargs):
self._download_moment_to_buffer(attr, rank=1)
tmp = self.buffer.copy()
tmp[:] *= conc
- tmp[:] *= DUMMY_SPECIFIC_GRAVITY * Md
+ tmp[:] *= DUMMY_SPECIFIC_GRAVITY * self.formulae.constants.Md
self._download_to_buffer(self.particulator.environment['rhod'])
tmp[:] /= self.particulator.mesh.dv
diff --git a/PySDM/products/aqueous_chemistry/gaseous_mole_fraction.py b/PySDM/products/aqueous_chemistry/gaseous_mole_fraction.py
index 6f6918d37..1c979f09b 100644
--- a/PySDM/products/aqueous_chemistry/gaseous_mole_fraction.py
+++ b/PySDM/products/aqueous_chemistry/gaseous_mole_fraction.py
@@ -2,7 +2,7 @@
mole fractions of gaseous compounds relevant for aqueous chemistry
"""
from PySDM.products.impl.product import Product
-from PySDM.physics.aqueous_chemistry.support import GASEOUS_COMPOUNDS, SPECIFIC_GRAVITY
+from PySDM.physics.aqueous_chemistry.support import GASEOUS_COMPOUNDS
class GaseousMoleFraction(Product):
@@ -18,6 +18,6 @@ def register(self, builder):
def _impl(self, **kwargs):
tmp = self.formulae.trivia.mixing_ratio_2_mole_fraction(
self.aqueous_chemistry.environment_mixing_ratios[self.compound],
- specific_gravity=SPECIFIC_GRAVITY[self.compound]
+ specific_gravity=self.aqueous_chemistry.specific_gravities[self.compound]
)
return tmp
diff --git a/PySDM/products/displacement/surface_precipitation.py b/PySDM/products/displacement/surface_precipitation.py
index 5fb7ca146..be7a2eb18 100644
--- a/PySDM/products/displacement/surface_precipitation.py
+++ b/PySDM/products/displacement/surface_precipitation.py
@@ -2,7 +2,6 @@
water volume flux derived from sizes of particles crossing bottom domain boundary
"""
from PySDM.products.impl.product import Product
-from PySDM.physics.constants import rho_w
class SurfacePrecipitation(Product):
@@ -31,7 +30,8 @@ def _impl(self, **kwargs) -> float:
return 0.
# TODO #708
- result = rho_w * self.accumulated_rainfall / self.elapsed_time / (self.dv / self.dz)
+ result = self.formulae.constants.rho_w \
+ * self.accumulated_rainfall / self.elapsed_time / (self.dv / self.dz)
self._reset_counters()
return result
diff --git a/PySDM/products/freezing/ice_water_content.py b/PySDM/products/freezing/ice_water_content.py
index b9266e94a..df7f9f540 100644
--- a/PySDM/products/freezing/ice_water_content.py
+++ b/PySDM/products/freezing/ice_water_content.py
@@ -3,7 +3,6 @@
"""
import numpy as np
from PySDM.products.impl.moment_product import MomentProduct
-from PySDM.physics import constants as const
class IceWaterContent(MomentProduct):
@@ -18,7 +17,7 @@ def _impl(self, **kwargs):
self._download_moment_to_buffer('volume', rank=0, filter_range=(-np.inf, 0))
conc = self.buffer
- result[:] *= -const.rho_i * conc / self.particulator.mesh.dv
+ result[:] *= -self.formulae.constants.rho_i * conc / self.particulator.mesh.dv
if self.specific:
self._download_to_buffer(self.particulator.environment['rhod'])
diff --git a/PySDM/products/impl/product.py b/PySDM/products/impl/product.py
index 8975728b5..e9d9098ae 100644
--- a/PySDM/products/impl/product.py
+++ b/PySDM/products/impl/product.py
@@ -7,7 +7,7 @@
import inspect
import numpy as np
import pint
-from PySDM.physics import constants as const
+from PySDM.physics.constants import PPB, PPM, PPT
_UNIT_REGISTRY = pint.UnitRegistry()
_CAMEL_CASE_PATTERN = re.compile(r'[A-Z]?[a-z]+|[A-Z]+(?![^A-Z])')
@@ -40,11 +40,11 @@ def _parse_unit(unit: str):
if unit in ('%', 'percent'):
return .01 * _UNIT_REGISTRY.dimensionless
if unit in ('PPB', 'ppb'):
- return const.PPB * _UNIT_REGISTRY.dimensionless
+ return PPB * _UNIT_REGISTRY.dimensionless
if unit in ('PPM', 'ppm'):
- return const.PPM * _UNIT_REGISTRY.dimensionless
+ return PPM * _UNIT_REGISTRY.dimensionless
if unit in ('PPT', 'ppt'):
- return const.PPT * _UNIT_REGISTRY.dimensionless
+ return PPT * _UNIT_REGISTRY.dimensionless
return _UNIT_REGISTRY.parse_expression(unit)
@staticmethod
diff --git a/PySDM/products/impl/spectrum_moment_product.py b/PySDM/products/impl/spectrum_moment_product.py
index 201d49214..cab9c0bc9 100644
--- a/PySDM/products/impl/spectrum_moment_product.py
+++ b/PySDM/products/impl/spectrum_moment_product.py
@@ -9,7 +9,6 @@
class SpectrumMomentProduct(ABC, Product):
def __init__(self, name, unit, attr_unit):
super().__init__(name=name, unit=unit)
- _ = self._parse_unit(unit)
self.attr_bins_edges = None
self.attr_unit = attr_unit
self.moment_0 = None
@@ -21,6 +20,7 @@ def register(self, builder):
(len(self.attr_bins_edges) - 1, self.particulator.mesh.n_cell), dtype=float)
self.moments = self.particulator.Storage.empty(
(len(self.attr_bins_edges) - 1, self.particulator.mesh.n_cell), dtype=float)
+ _ = self._parse_unit(self.attr_unit)
def _recalculate_spectrum_moment(
self, attr,
diff --git a/PySDM/products/size_spectral/mean_radius.py b/PySDM/products/size_spectral/mean_radius.py
index 2476d30b4..2d9df2958 100644
--- a/PySDM/products/size_spectral/mean_radius.py
+++ b/PySDM/products/size_spectral/mean_radius.py
@@ -1,7 +1,6 @@
"""
mean radius of particles within a grid cell (optionally restricted to a given size range)
"""
-from PySDM.physics import constants as const
from PySDM.products.impl.moment_product import MomentProduct
@@ -11,5 +10,5 @@ def __init__(self, name=None, unit='m'):
def _impl(self, **kwargs):
self._download_moment_to_buffer('volume', rank=1 / 3)
- self.buffer[:] /= const.PI_4_3 ** (1 / 3)
+ self.buffer[:] /= self.formulae.constants.PI_4_3 ** (1 / 3)
return self.buffer
diff --git a/PySDM/products/size_spectral/water_mixing_ratio.py b/PySDM/products/size_spectral/water_mixing_ratio.py
index 9baabc5e9..c15cc900a 100644
--- a/PySDM/products/size_spectral/water_mixing_ratio.py
+++ b/PySDM/products/size_spectral/water_mixing_ratio.py
@@ -4,7 +4,6 @@
"""
import numpy as np
from PySDM.products.impl.moment_product import MomentProduct
-from PySDM.physics import constants as const
class WaterMixingRatio(MomentProduct):
@@ -27,7 +26,7 @@ def _impl(self, **kwargs): # TODO #217
self._download_moment_to_buffer('volume', rank=1,
filter_range=self.volume_range, filter_attr='volume')
result = self.buffer.copy()
- result[:] *= const.rho_w
+ result[:] *= self.formulae.constants.rho_w
result[:] *= conc
result[:] /= self.particulator.mesh.dv
diff --git a/README.md b/README.md
index 1be872931..979df67c4 100644
--- a/README.md
+++ b/README.md
@@ -247,7 +247,7 @@ In the listing below, its usage is interleaved with plotting logic
Julia (click to expand)
```Julia
-rho_w = pyimport("PySDM.physics.constants").rho_w
+rho_w = pyimport("PySDM.physics.constants_defaults").rho_w
using Plots; plotlyjs()
for step = 0:1200:3600
@@ -269,10 +269,10 @@ savefig("plot.svg")
Matlab (click to expand)
```Matlab
-rho_w = py.importlib.import_module('PySDM.physics.constants').rho_w;
+rho_w = py.importlib.import_module('PySDM.physics.constants_defaults').rho_w;
for step = 0:1200:3600
- particulator.run(int32(step - particulator.n_steps))
+ particulator.run(int32(step - particulator.n_steps));
x = radius_bins_edges / si.um;
y = particulator.products{"dv/dlnr"}.get() * rho_w / si.g;
stairs(...
@@ -293,7 +293,7 @@ legend()
Python (click to expand)
```Python
-from PySDM.physics.constants import rho_w
+from PySDM.physics.constants_defaults import rho_w
from matplotlib import pyplot
for step in [0, 1200, 2400, 3600]:
@@ -448,15 +448,15 @@ output_interval = 4;
output_points = 40;
n_sd = 256;
-formulae = Formulae()
+formulae = Formulae();
builder = Builder(pyargs('backend', CPU(formulae), 'n_sd', int32(n_sd)));
builder.set_environment(env);
-builder.add_dynamic(AmbientThermodynamics())
-builder.add_dynamic(Condensation())
+builder.add_dynamic(AmbientThermodynamics());
+builder.add_dynamic(Condensation());
tmp = spectral_sampling.Logarithmic(spectrum).sample(int32(n_sd));
r_dry = tmp{1};
-v_dry = formulae.trivia.volume(r_dry);
+v_dry = formulae.trivia.volume(pyargs('radius', r_dry));
specific_concentration = tmp{2};
r_wet = equilibrate_wet_radii(r_dry, env, kappa * v_dry);
@@ -464,7 +464,7 @@ attributes = py.dict(pyargs( ...
'n', discretise_multiplicities(specific_concentration * env.mass_of_dry_air), ...
'dry volume', v_dry, ...
'kappa times dry volume', kappa * v_dry, ...
- 'volume', formulae.trivia.volume(r_wet) ...
+ 'volume', formulae.trivia.volume(pyargs('radius', r_wet)) ...
));
particulator = builder.build(attributes, py.list({ ...
@@ -513,14 +513,13 @@ for pykey = py.list(keys(particulator.products))
end
i=i+1;
end
-saveas(gcf, "parcel.png")
+saveas(gcf, "parcel.png");
```
Python (click to expand)
```Python
-import PySDM.initialisation.spectra.lognormal
from matplotlib import pyplot
from PySDM.physics import si
from PySDM.initialisation import discretise_multiplicities, equilibrate_wet_radii
diff --git a/test-time-requirements.txt b/test-time-requirements.txt
index cf49cf5c5..016aab458 100644
--- a/test-time-requirements.txt
+++ b/test-time-requirements.txt
@@ -4,4 +4,4 @@ ghapi
pytest
# note: if cloning both PySDM and PySDM examples, consider "pip install -e"
-PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@ed59c6e#egg=PySDM-examples
+PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@1f08b2f#egg=PySDM-examples
diff --git a/tests/smoke_tests/alpert_and_knopf_2016/test_ak16_fig_1.py b/tests/smoke_tests/alpert_and_knopf_2016/test_ak16_fig_1.py
index 59c1044bc..1802fedb1 100644
--- a/tests/smoke_tests/alpert_and_knopf_2016/test_ak16_fig_1.py
+++ b/tests/smoke_tests/alpert_and_knopf_2016/test_ak16_fig_1.py
@@ -3,7 +3,6 @@
from matplotlib import pylab
import pytest
from PySDM_examples.Alpert_and_Knopf_2016 import simulation, Table1
-from PySDM.physics.heterogeneous_ice_nucleation_rate import constant
from PySDM.physics import si
n_runs_per_case = 3
@@ -12,8 +11,6 @@
@pytest.mark.parametrize("multiplicity", (1, 2, 10))
def test_ak16_fig_1(multiplicity, plot=False):
# Arrange
- constant.J_HET = 1e3 / si.cm ** 2 / si.s
-
dt = 1 * si.s
total_time = 6 * si.min
@@ -36,6 +33,7 @@ def test_ak16_fig_1(multiplicity, plot=False):
n_sd = int(n_sd)
data, _ = simulation(
+ constants={'J_HET': 1e3 / si.cm ** 2 / si.s},
seed=i, n_sd=n_sd, time_step=dt, volume=dv, spectrum=case['ISA'],
droplet_volume=droplet_volume, multiplicity=multiplicity,
total_time=total_time, number_of_real_droplets=number_of_real_droplets
diff --git a/tests/smoke_tests/arabas_and_shima_2017/test_conservation.py b/tests/smoke_tests/arabas_and_shima_2017/test_conservation.py
index e2d755462..423cb24d4 100644
--- a/tests/smoke_tests/arabas_and_shima_2017/test_conservation.py
+++ b/tests/smoke_tests/arabas_and_shima_2017/test_conservation.py
@@ -5,14 +5,14 @@
from PySDM_examples.Arabas_and_Shima_2017.settings import setups
from PySDM_examples.Arabas_and_Shima_2017.settings import Settings, w_avgs
from PySDM.backends.impl_numba.test_helpers import bdf
-from PySDM.physics import constants as const
+from PySDM.physics.constants_defaults import rho_w
from PySDM.backends import CPU, GPU
def liquid_water_mixing_ratio(simulation: Simulation):
droplet_volume = simulation.particulator.attributes['volume'].to_ndarray()[0]
droplet_number = simulation.particulator.attributes['n'].to_ndarray()[0]
- droplet_mass = droplet_number * droplet_volume * const.rho_w
+ droplet_mass = droplet_number * droplet_volume * rho_w
env = simulation.particulator.environment
return droplet_mass / env.mass_of_dry_air
diff --git a/tests/smoke_tests/arabas_and_shima_2017/test_event_rates.py b/tests/smoke_tests/arabas_and_shima_2017/test_event_rates.py
index e63672750..9fc1ffcf1 100644
--- a/tests/smoke_tests/arabas_and_shima_2017/test_event_rates.py
+++ b/tests/smoke_tests/arabas_and_shima_2017/test_event_rates.py
@@ -4,7 +4,7 @@
from PySDM_examples.Arabas_and_Shima_2017.simulation import Simulation
from PySDM_examples.Arabas_and_Shima_2017.settings import setups
from PySDM_examples.Arabas_and_Shima_2017.settings import Settings, w_avgs
-from PySDM.physics.constants import si, rho_STP, convert_to
+from PySDM.physics.constants import si, convert_to
@pytest.mark.parametrize("settings_idx", range(len(w_avgs)))
@@ -16,6 +16,7 @@ def test_event_rates(settings_idx):
r_dry=setups[settings_idx].r_dry,
mass_of_dry_air=1 * si.kg
)
+ const = settings.formulae.constants
settings.n_output = 50
simulation = Simulation(settings)
@@ -26,7 +27,7 @@ def test_event_rates(settings_idx):
rip = np.asarray(output['ripening_rate'])
act = np.asarray(output['activating_rate'])
dea = np.asarray(output['deactivating_rate'])
- act_max = np.full((1,), settings.n_in_dv / simulation.particulator.dt / rho_STP)
+ act_max = np.full((1,), settings.n_in_dv / simulation.particulator.dt / const.rho_STP)
convert_to(act_max, 1/si.mg)
assert (rip == 0).all()
assert (act > 0).any()
diff --git a/tests/smoke_tests/arabas_and_shima_2017/test_initialisation.py b/tests/smoke_tests/arabas_and_shima_2017/test_initialisation.py
index c4e38f177..46870da2e 100644
--- a/tests/smoke_tests/arabas_and_shima_2017/test_initialisation.py
+++ b/tests/smoke_tests/arabas_and_shima_2017/test_initialisation.py
@@ -4,11 +4,11 @@
from PySDM_examples.Arabas_and_Shima_2017.settings import setups
from PySDM_examples.Arabas_and_Shima_2017.simulation import Simulation
from PySDM import Formulae
-from PySDM.physics import constants as const
+CONST = Formulae().constants
-class TestInitialisation:
+class TestInitialisation:
@staticmethod
def simulation_test(var, expected, setup):
simulation = Simulation(setup)
@@ -27,8 +27,8 @@ def test_T_initialisation(settings_idx):
@pytest.mark.parametrize("settings_idx", range(len(setups)))
def test_RH_initialisation(settings_idx):
setup = setups[settings_idx]
- pv0 = setup.p0 / (1 + const.eps / setup.q0)
- pvs = setup.formulae.saturation_vapour_pressure.pvs_Celsius(setup.T0 - const.T0)
+ pv0 = setup.p0 / (1 + CONST.eps / setup.q0)
+ pvs = setup.formulae.saturation_vapour_pressure.pvs_Celsius(setup.T0 - CONST.T0)
TestInitialisation.simulation_test('RH', pv0 / pvs, setup)
@staticmethod
@@ -47,16 +47,16 @@ def test_qv_initialisation(settings_idx):
@pytest.mark.parametrize("settings_idx", range(len(setups)))
def test_rhod_initialisation(settings_idx):
setup = setups[settings_idx]
- pv0 = setup.p0 / (1 + const.eps / setup.q0)
+ pv0 = setup.p0 / (1 + CONST.eps / setup.q0)
pd0 = setup.p0 - pv0
- rhod0 = pd0 / const.Rd / setup.T0
+ rhod0 = pd0 / CONST.Rd / setup.T0
TestInitialisation.simulation_test('rhod', rhod0, setup)
@staticmethod
@pytest.mark.parametrize("settings_idx", range(len(setups)))
def test_thd_initialisation(settings_idx):
setup = setups[settings_idx]
- pv0 = setup.p0 / (1 + const.eps / setup.q0)
+ pv0 = setup.p0 / (1 + CONST.eps / setup.q0)
pd0 = setup.p0 - pv0
phys = Formulae().trivia
thd0 = phys.th_std(pd0, setup.T0)
diff --git a/tests/smoke_tests/arabas_et_al_2015/test_demo_settings.py b/tests/smoke_tests/arabas_et_al_2015/test_demo_settings.py
deleted file mode 100644
index 92134c17d..000000000
--- a/tests/smoke_tests/arabas_et_al_2015/test_demo_settings.py
+++ /dev/null
@@ -1,7 +0,0 @@
-# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
-from PySDM_examples.Arabas_et_al_2015 import Settings
-from PySDM_examples.Szumowski_et_al_1998.gui_settings import GUISettings
-
-
-def test_gui_settings():
- _ = GUISettings(Settings())
diff --git a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py
index 60f40b9fc..ce6c9a0d2 100644
--- a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py
+++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py
@@ -5,16 +5,8 @@
from PySDM import Formulae
from PySDM.physics import si
from PySDM.backends import CPU
-from PySDM.physics.freezing_temperature_spectrum import niemand_et_al_2012
-from PySDM.physics.heterogeneous_ice_nucleation_rate import abifm
from .dummy_storage import DummyStorage
-# TODO #599
-niemand_et_al_2012.a = -0.517
-niemand_et_al_2012.b = 8.934
-abifm.m = 28.13797
-abifm.c = -2.92414
-
@pytest.mark.parametrize("singular", (
pytest.param(False, id="singular: False"),
@@ -28,7 +20,13 @@ def test_freezing(singular):
condensation_coordinate='VolumeLogarithm',
fastmath=True,
freezing_temperature_spectrum='Niemand_et_al_2012',
- heterogeneous_ice_nucleation_rate='ABIFM'
+ heterogeneous_ice_nucleation_rate='ABIFM',
+ constants={
+ 'NIEMAND_A': -0.517,
+ 'NIEMAND_B': 8.934,
+ 'ABIFM_M': 28.13797,
+ 'ABIFM_C': -2.92414
+ }
))
settings.dt = .5 * si.second
settings.grid = (5, 15)
diff --git a/tests/smoke_tests/kreidenweis_et_al_2003/test_ionic_strength.py b/tests/smoke_tests/kreidenweis_et_al_2003/test_ionic_strength.py
index 44001ffe9..12e5be3b4 100644
--- a/tests/smoke_tests/kreidenweis_et_al_2003/test_ionic_strength.py
+++ b/tests/smoke_tests/kreidenweis_et_al_2003/test_ionic_strength.py
@@ -4,7 +4,7 @@
from chempy.electrolytes import ionic_strength
from PySDM_examples.Kreidenweis_et_al_2003 import Settings, Simulation
from PySDM.backends.impl_numba.methods.chemistry_methods import calc_ionic_strength, _K, _conc
-from PySDM.physics.constants import rho_w, ROOM_TEMP, K_H2O
+from PySDM.physics.constants import K_H2O
from PySDM import Formulae
from PySDM.physics.aqueous_chemistry.support import EquilibriumConsts
@@ -13,16 +13,17 @@
@pytest.mark.parametrize("n_sd", (10, 20))
def test_calc_ionic_strength(nt, n_sd):
formulae = Formulae()
+ const = formulae.constants
EQUILIBRIUM_CONST = EquilibriumConsts(formulae).EQUILIBRIUM_CONST
K = _K(
- NH3=EQUILIBRIUM_CONST["K_NH3"].at(ROOM_TEMP),
- SO2=EQUILIBRIUM_CONST["K_SO2"].at(ROOM_TEMP),
- HSO3=EQUILIBRIUM_CONST["K_HSO3"].at(ROOM_TEMP),
- HSO4=EQUILIBRIUM_CONST["K_HSO4"].at(ROOM_TEMP),
- HCO3=EQUILIBRIUM_CONST["K_HCO3"].at(ROOM_TEMP),
- CO2=EQUILIBRIUM_CONST["K_CO2"].at(ROOM_TEMP),
- HNO3=EQUILIBRIUM_CONST["K_HNO3"].at(ROOM_TEMP)
+ NH3=EQUILIBRIUM_CONST["K_NH3"].at(const.ROOM_TEMP),
+ SO2=EQUILIBRIUM_CONST["K_SO2"].at(const.ROOM_TEMP),
+ HSO3=EQUILIBRIUM_CONST["K_HSO3"].at(const.ROOM_TEMP),
+ HSO4=EQUILIBRIUM_CONST["K_HSO4"].at(const.ROOM_TEMP),
+ HCO3=EQUILIBRIUM_CONST["K_HCO3"].at(const.ROOM_TEMP),
+ CO2=EQUILIBRIUM_CONST["K_CO2"].at(const.ROOM_TEMP),
+ HNO3=EQUILIBRIUM_CONST["K_HNO3"].at(const.ROOM_TEMP)
)
settings = Settings(dt=1, n_sd=n_sd, n_substep=5)
@@ -59,16 +60,16 @@ def test_calc_ionic_strength(nt, n_sd):
)
expected = ionic_strength({
- 'H+': conc['H+'] / rho_w,
- 'HCO3-': K.CO2 / conc['H+'] * conc['C+4'] / alpha_C / rho_w,
- 'CO3-2': K.CO2 / conc['H+'] * K.HCO3 / conc['H+'] * conc['C+4'] / alpha_C / rho_w,
- 'HSO3-': K.SO2 / conc['H+'] * conc['S+4'] / alpha_S / rho_w,
- 'SO3-2': K.SO2 / conc['H+'] * K.HSO3 / conc['H+'] * conc['S+4'] / alpha_S / rho_w,
- 'NH4+': K.NH3 / K_H2O * conc['H+'] * conc['N-3'] / alpha_N3 / rho_w,
- 'NO3-': K.HNO3 / conc['H+'] * conc['N+5'] / alpha_N5 / rho_w,
- 'HSO4-': conc['H+'] * conc['S+6'] / (conc['H+'] + K.HSO4) / rho_w,
- 'SO4-2': K.HSO4 * conc['S+6'] / (conc['H+'] + K.HSO4) / rho_w,
- 'OH-': K_H2O / conc['H+'] / rho_w
- }, warn=False) * rho_w
+ 'H+': conc['H+'] / const.rho_w,
+ 'HCO3-': K.CO2 / conc['H+'] * conc['C+4'] / alpha_C / const.rho_w,
+ 'CO3-2': K.CO2 / conc['H+'] * K.HCO3 / conc['H+'] * conc['C+4'] / alpha_C / const.rho_w,
+ 'HSO3-': K.SO2 / conc['H+'] * conc['S+4'] / alpha_S / const.rho_w,
+ 'SO3-2': K.SO2 / conc['H+'] * K.HSO3 / conc['H+'] * conc['S+4'] / alpha_S / const.rho_w,
+ 'NH4+': K.NH3 / K_H2O * conc['H+'] * conc['N-3'] / alpha_N3 / const.rho_w,
+ 'NO3-': K.HNO3 / conc['H+'] * conc['N+5'] / alpha_N5 / const.rho_w,
+ 'HSO4-': conc['H+'] * conc['S+6'] / (conc['H+'] + K.HSO4) / const.rho_w,
+ 'SO4-2': K.HSO4 * conc['S+6'] / (conc['H+'] + K.HSO4) / const.rho_w,
+ 'OH-': K_H2O / conc['H+'] / const.rho_w
+ }, warn=False) * const.rho_w
np.testing.assert_allclose(actual, expected, rtol=1e-15)
diff --git a/tests/smoke_tests/kreidenweis_et_al_2003/test_table_3.py b/tests/smoke_tests/kreidenweis_et_al_2003/test_table_3.py
index 7539debab..6c740267e 100644
--- a/tests/smoke_tests/kreidenweis_et_al_2003/test_table_3.py
+++ b/tests/smoke_tests/kreidenweis_et_al_2003/test_table_3.py
@@ -4,7 +4,7 @@
from PySDM_examples.Kreidenweis_et_al_2003 import Settings, Simulation
from PySDM.physics import si
from PySDM.physics.constants import convert_to, PPB
-from PySDM.physics.aqueous_chemistry.support import AQUEOUS_COMPOUNDS, SPECIFIC_GRAVITY
+from PySDM.physics.aqueous_chemistry.support import AQUEOUS_COMPOUNDS, SpecificGravities
class TestTable3:
@@ -14,6 +14,7 @@ def test_at_t_0():
settings = Settings(n_sd=100, dt=1 * si.s, n_substep=5)
settings.t_max = 0
simulation = Simulation(settings)
+ specific_gravities = SpecificGravities(simulation.particulator.formulae.constants)
# Act
output = simulation.run()
@@ -53,7 +54,7 @@ def test_at_t_0():
actual=(
settings.formulae.trivia.mole_fraction_2_mixing_ratio(
mole_fraction,
- specific_gravity=SPECIFIC_GRAVITY[compound]
+ specific_gravity=specific_gravities[compound]
) * np.asarray(output['rhod'])
),
desired=expected[key],
diff --git a/tests/smoke_tests/lowe_et_al_2019/test_fig_1.py b/tests/smoke_tests/lowe_et_al_2019/test_fig_1.py
index ca7e83b5a..d9115178a 100644
--- a/tests/smoke_tests/lowe_et_al_2019/test_fig_1.py
+++ b/tests/smoke_tests/lowe_et_al_2019/test_fig_1.py
@@ -4,8 +4,7 @@
import pytest
from PySDM_examples.Lowe_et_al_2019 import aerosol
from PySDM import Formulae
-from PySDM.physics import si, constants as const
-from PySDM.physics.surface_tension import compressed_film_ovadnevaite
+from PySDM.physics import si, constants_defaults as const
from .constants import constants
assert hasattr(constants, '_pytestfixturefunction')
@@ -43,7 +42,13 @@ def test_bulk_surface_tension_is_sgm_w():
# pylint: disable=redefined-outer-name,unused-argument
def test_kink_location(constants, aerosol, cutoff):
# arrange
- formulae = Formulae(surface_tension='CompressedFilmOvadnevaite')
+ formulae = Formulae(
+ surface_tension='CompressedFilmOvadnevaite',
+ constants={
+ 'sgm_org': 40 * si.mN / si.m,
+ 'delta_min': 0.1 * si.nm
+ }
+ )
# act
sigma = formulae.surface_tension.sigma(
@@ -51,8 +56,8 @@ def test_kink_location(constants, aerosol, cutoff):
# assert
cutoff_idx = (np.abs(R_WET - cutoff)).argmin()
- assert (sigma[:cutoff_idx] == compressed_film_ovadnevaite.sgm_org).all()
- assert sigma[cutoff_idx] > compressed_film_ovadnevaite.sgm_org
+ assert (sigma[:cutoff_idx] == formulae.constants.sgm_org).all()
+ assert sigma[cutoff_idx] > formulae.constants.sgm_org
assert .98 * const.sgm_w < sigma[-1] <= const.sgm_w
@staticmethod
@@ -70,7 +75,13 @@ def test_kink_location(constants, aerosol, cutoff):
def test_koehler_maxima(constants, aerosol, surface_tension, maximum_x, maximum_y, bimodal):
# arrange
label = {'CompressedFilmOvadnevaite': 'film', 'Constant': 'bulk'}[surface_tension]
- formulae = Formulae(surface_tension=surface_tension)
+ formulae = Formulae(
+ surface_tension=surface_tension,
+ constants={
+ 'sgm_org': 40 * si.mN / si.m,
+ 'delta_min': 0.1 * si.nm
+ }
+ )
sigma = formulae.surface_tension.sigma(
np.nan, V_WET, V_DRY, aerosol.aerosol_modes_per_cc[0]['f_org'])
RH_eq = formulae.hygroscopicity.RH_eq(
diff --git a/tests/unit_tests/attributes/test_acidity.py b/tests/unit_tests/attributes/test_acidity.py
index cf82e3671..2322916f1 100644
--- a/tests/unit_tests/attributes/test_acidity.py
+++ b/tests/unit_tests/attributes/test_acidity.py
@@ -6,7 +6,7 @@
from chempy.equilibria import EqSystem
from chempy.chemistry import Species
from PySDM.physics.aqueous_chemistry.support import M, EquilibriumConsts
-from PySDM.physics.constants import ROOM_TEMP, K_H2O
+from PySDM.physics.constants import K_H2O
from PySDM.formulae import Formulae
from PySDM.backends.impl_numba.methods.chemistry_methods import ChemistryMethods, _K, _conc
from PySDM.dynamics import aqueous_chemistry
@@ -22,7 +22,7 @@ def test_equilibrate_pH_pure_water():
# Arrange
eqs = {}
for key, const in EQUILIBRIUM_CONST.items():
- eqs[key] = np.full(1, const.at(ROOM_TEMP))
+ eqs[key] = np.full(1, const.at(FORMULAE.constants.ROOM_TEMP))
# Act
result = np.empty(1)
@@ -65,7 +65,11 @@ def test_equilibrate_pH_pure_water():
defaultdict(float, {'H2O': 1, 'NH3': 5e-3, 'H2CO3(aq)': 0.01e-3, 'H2SO3(aq)': 0.005e-3}),
defaultdict(float, {'H2O': 1, 'NH3': .5e-3, 'H2CO3(aq)': 0.1e-3, 'H2SO3(aq)': 0.05e-3}),
))
- @pytest.mark.parametrize('env_T', (ROOM_TEMP, ROOM_TEMP-30, ROOM_TEMP+30))
+ @pytest.mark.parametrize('env_T', (
+ FORMULAE.constants.ROOM_TEMP,
+ FORMULAE.constants.ROOM_TEMP-30,
+ FORMULAE.constants.ROOM_TEMP+30
+ ))
def test_equilibrate_pH_non_trivial(init_conc, env_T):
equilibria = {
diff --git a/tests/unit_tests/backends/test_freezing_methods.py b/tests/unit_tests/backends/test_freezing_methods.py
index ddab1cc3f..a7cf92669 100644
--- a/tests/unit_tests/backends/test_freezing_methods.py
+++ b/tests/unit_tests/backends/test_freezing_methods.py
@@ -1,8 +1,7 @@
# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
from matplotlib import pylab
import numpy as np
-from PySDM.physics import constants as const
-from PySDM.physics.heterogeneous_ice_nucleation_rate import constant
+from PySDM.physics import constants_defaults as const
from PySDM import Builder, Formulae
from PySDM.backends import CPU
from PySDM.environments import Box
@@ -33,7 +32,6 @@ def test_freeze_time_dependent(plot=False):
)
rate = 1e-9
immersed_surface_area = 1
- constant.J_HET = rate / immersed_surface_area
number_of_real_droplets = 1024
total_time = 2e9 # effectively interpretted here as seconds, i.e. cycle = 1 * si.s
@@ -56,7 +54,12 @@ def test_freeze_time_dependent(plot=False):
key = f"{case['dt']}:{case['N']}"
output[key] = {'unfrozen_fraction': [], 'dt': case['dt'], 'N': case['N']}
- formulae = Formulae(heterogeneous_ice_nucleation_rate='Constant')
+ formulae = Formulae(
+ heterogeneous_ice_nucleation_rate='Constant',
+ constants={
+ 'J_HET': rate / immersed_surface_area
+ }
+ )
builder = Builder(n_sd=n_sd, backend=CPU(formulae=formulae))
env = Box(dt=case['dt'], dv=d_v)
builder.set_environment(env)
diff --git a/tests/unit_tests/backends/test_oxidation.py b/tests/unit_tests/backends/test_oxidation.py
index 76a64e7b8..5a64d8416 100644
--- a/tests/unit_tests/backends/test_oxidation.py
+++ b/tests/unit_tests/backends/test_oxidation.py
@@ -7,7 +7,8 @@
from PySDM.physics import si
from PySDM.physics.aqueous_chemistry.support import KineticConsts, EquilibriumConsts, \
DISSOCIATION_FACTORS, k4
-from PySDM.physics.constants import T_STP, PI_4_3
+from PySDM.physics.constants import PI_4_3
+from PySDM.physics.constants_defaults import T_STP
formulae = Formulae()
@@ -22,20 +23,18 @@ def __init__(self):
kinetic_consts = KineticConsts(formulae)
equilibrium_consts = EquilibriumConsts(formulae)
-T = T_STP
-
-k0 = Storage.from_ndarray(np.full(1, kinetic_consts.KINETIC_CONST['k0'].at(T)))
-k1 = Storage.from_ndarray(np.full(1, kinetic_consts.KINETIC_CONST['k1'].at(T)))
-k2 = Storage.from_ndarray(np.full(1, kinetic_consts.KINETIC_CONST['k2'].at(T)))
-k3 = Storage.from_ndarray(np.full(1, kinetic_consts.KINETIC_CONST['k3'].at(T)))
-K_SO2 = Storage.from_ndarray(np.full(1, equilibrium_consts.EQUILIBRIUM_CONST['K_SO2'].at(T)))
-K_HSO3 = Storage.from_ndarray(np.full(1, equilibrium_consts.EQUILIBRIUM_CONST['K_HSO3'].at(T)))
+k0 = Storage.from_ndarray(np.full(1, kinetic_consts.KINETIC_CONST['k0'].at(T_STP)))
+k1 = Storage.from_ndarray(np.full(1, kinetic_consts.KINETIC_CONST['k1'].at(T_STP)))
+k2 = Storage.from_ndarray(np.full(1, kinetic_consts.KINETIC_CONST['k2'].at(T_STP)))
+k3 = Storage.from_ndarray(np.full(1, kinetic_consts.KINETIC_CONST['k3'].at(T_STP)))
+K_SO2 = Storage.from_ndarray(np.full(1, equilibrium_consts.EQUILIBRIUM_CONST['K_SO2'].at(T_STP)))
+K_HSO3 = Storage.from_ndarray(np.full(1, equilibrium_consts.EQUILIBRIUM_CONST['K_HSO3'].at(T_STP)))
volume = PI_4_3 * (1 * si.um) ** 3
pH = 5.
n_sd = 1
eqc = {
- k: Storage.from_ndarray(np.full(n_sd, equilibrium_consts.EQUILIBRIUM_CONST[k].at(T)))
+ k: Storage.from_ndarray(np.full(n_sd, equilibrium_consts.EQUILIBRIUM_CONST[k].at(T_STP)))
for k in ('K_HSO3', 'K_SO2')
}
cell_ids = Storage.from_ndarray(np.zeros(n_sd, dtype=int))
diff --git a/tests/unit_tests/initialisation/test_r_wet_init.py b/tests/unit_tests/initialisation/test_r_wet_init.py
index d5f2aef9c..e46dd5fcb 100644
--- a/tests/unit_tests/initialisation/test_r_wet_init.py
+++ b/tests/unit_tests/initialisation/test_r_wet_init.py
@@ -4,26 +4,16 @@
from matplotlib import pyplot
from PySDM.initialisation import equilibrate_wet_radii
from PySDM import Formulae
-from PySDM.physics import si, constants as const
-from PySDM.physics.surface_tension import compressed_film_ovadnevaite
+from PySDM.physics import si, constants_defaults as const
from PySDM.backends import CPU
-@pytest.fixture()
-def constants():
- compressed_film_ovadnevaite.sgm_org = 40 * si.mN / si.m
- compressed_film_ovadnevaite.delta_min = 0.1 * si.nm
- yield
- compressed_film_ovadnevaite.sgm_org = np.nan
- compressed_film_ovadnevaite.delta_min = np.nan
-
-
@pytest.mark.parametrize('r_dry', [
pytest.param(2.4e-09),
pytest.param(2.5e-09)
])
# pylint: disable=unused-argument,redefined-outer-name
-def test_r_wet_init(constants, r_dry, plot=False):
+def test_r_wet_init(r_dry, plot=False):
# Arrange
T = 280
RH = .9
@@ -31,7 +21,13 @@ def test_r_wet_init(constants, r_dry, plot=False):
kappa = .356
class Particulator:
- formulae = Formulae(surface_tension='CompressedFilmOvadnevaite')
+ formulae = Formulae(
+ surface_tension='CompressedFilmOvadnevaite',
+ constants={
+ 'sgm_org': 40 * si.mN / si.m,
+ 'delta_min': 0.1 * si.nm
+ }
+ )
class Env:
particulator = Particulator()
diff --git a/tests/unit_tests/initialisation/test_spectral_discretisation.py b/tests/unit_tests/initialisation/test_spectral_discretisation.py
index 45e9afba2..e674ed2c7 100644
--- a/tests/unit_tests/initialisation/test_spectral_discretisation.py
+++ b/tests/unit_tests/initialisation/test_spectral_discretisation.py
@@ -4,11 +4,7 @@
from PySDM.initialisation.sampling import spectral_sampling, spectro_glacial_sampling
from PySDM.initialisation.spectra.lognormal import Lognormal
from PySDM import Formulae
-from PySDM.physics import constants as const
-from PySDM.physics.freezing_temperature_spectrum import niemand_et_al_2012
-niemand_et_al_2012.a = -0.517
-niemand_et_al_2012.b = 8.934
m_mode = .5e-5
n_part = 256 * 16
@@ -17,7 +13,10 @@
m_range = (.1e-6, 100e-6)
formulae = Formulae(
freezing_temperature_spectrum='Niemand_et_al_2012',
- seed=const.default_random_seed
+ constants={
+ 'NIEMAND_A': -0.517,
+ 'NIEMAND_B': 8.934
+ }
)
diff --git a/tests/unit_tests/physics/test_constants.py b/tests/unit_tests/physics/test_constants.py
index 74ad7cf11..67210b2c1 100644
--- a/tests/unit_tests/physics/test_constants.py
+++ b/tests/unit_tests/physics/test_constants.py
@@ -10,7 +10,6 @@ def consecutive_seeds():
for _ in range(5):
importlib.reload(constants)
seeds.append(constants.default_random_seed)
- print(seeds)
return np.asarray(seeds)
diff --git a/tests/unit_tests/physics/test_dimensional_analysis.py b/tests/unit_tests/physics/test_dimensional_analysis.py
index f3da783e5..e146d659f 100644
--- a/tests/unit_tests/physics/test_dimensional_analysis.py
+++ b/tests/unit_tests/physics/test_dimensional_analysis.py
@@ -3,7 +3,7 @@
import numba
from PySDM.physics.dimensional_analysis import DimensionalAnalysis
from PySDM.formulae import Formulae
-from PySDM.physics import constants
+from PySDM.physics import constants_defaults
assert numba.config.DISABLE_JIT is not None # pylint: disable=no-member
@@ -17,11 +17,11 @@ def test_fake_units():
sut = DimensionalAnalysis()
# Act & Assert
- assert isinstance(constants.D0, float)
+ assert isinstance(constants_defaults.D0, float)
with sut:
- assert not isinstance(constants.D0, float)
- assert isinstance(constants.D0.magnitude, float)
- assert isinstance(constants.D0, float)
+ assert not isinstance(constants_defaults.D0, float)
+ assert isinstance(constants_defaults.D0.magnitude, float)
+ assert isinstance(constants_defaults.D0, float)
@staticmethod
@pytest.mark.skipif("numba.config.DISABLE_JIT")
diff --git a/tests/unit_tests/physics/test_formulae.py b/tests/unit_tests/physics/test_formulae.py
index abcc6397e..5a435359c 100644
--- a/tests/unit_tests/physics/test_formulae.py
+++ b/tests/unit_tests/physics/test_formulae.py
@@ -1,17 +1,16 @@
# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
-from PySDM.physics import constants
+from PySDM.physics import constants_defaults
from PySDM.formulae import Formulae
from PySDM.physics.dimensional_analysis import DimensionalAnalysis
class TestFormulae:
-
@staticmethod
def test_pvs():
with DimensionalAnalysis():
# Arrange
formulae = Formulae()
- si = constants.si
+ si = constants_defaults.si
sut = formulae.saturation_vapour_pressure.pvs_Celsius
T = 300 * si.kelvins
@@ -25,14 +24,14 @@ def test_pvs():
def test_r_cr():
with DimensionalAnalysis():
# Arrange
- si = constants.si
+ si = constants_defaults.si
formulae = Formulae()
sut = formulae.hygroscopicity.r_cr
kp = .5
rd = .1 * si.micrometre
T = 300 * si.kelvins
- sgm = constants.sgm_w
+ sgm = constants_defaults.sgm_w
# Act
r_cr = sut(kp, rd**3, T, sgm)
@@ -44,7 +43,7 @@ def test_r_cr():
def test_lv():
with DimensionalAnalysis():
# Arrange
- si = constants.si
+ si = constants_defaults.si
T = 300 * si.kelvins
formulae = Formulae()
diff --git a/tests/unit_tests/physics/test_freezing_temperature_spectra.py b/tests/unit_tests/physics/test_freezing_temperature_spectra.py
index b616d39a2..cb756a38a 100644
--- a/tests/unit_tests/physics/test_freezing_temperature_spectra.py
+++ b/tests/unit_tests/physics/test_freezing_temperature_spectra.py
@@ -4,7 +4,6 @@
from matplotlib import pylab
from PySDM import Formulae
from PySDM.physics import si, constants as const
-from PySDM.physics.freezing_temperature_spectrum import niemand_et_al_2012, bigg_1953
A = 1 * si.um**2
@@ -16,12 +15,14 @@
))
def test_freezing_temperature_spectra(model, plot=False):
# Arrange
- bigg_1953.DT_median = 33
- niemand_et_al_2012.a = -0.517
- niemand_et_al_2012.b = 8.934
- niemand_et_al_2012.A_insol = 1 * si.um ** 2
-
- formulae = Formulae(freezing_temperature_spectrum=model)
+ formulae = Formulae(
+ freezing_temperature_spectrum=model,
+ constants={
+ 'NIEMAND_A': -0.517,
+ 'NIEMAND_B': 8.934,
+ 'BIGG_DT_MEDIAN': 33
+ }
+ )
temperature = np.linspace(const.T0, const.T0 - 40, num=100)
# Act
diff --git a/tests/unit_tests/physics/test_saturation_vapour_pressure.py b/tests/unit_tests/physics/test_saturation_vapour_pressure.py
index 3911a9697..b3aeda473 100644
--- a/tests/unit_tests/physics/test_saturation_vapour_pressure.py
+++ b/tests/unit_tests/physics/test_saturation_vapour_pressure.py
@@ -3,7 +3,7 @@
import numpy as np
from matplotlib import pyplot
from PySDM import Formulae
-from PySDM.physics import constants as const
+from PySDM.physics import constants_defaults as const
def test_saturation_vapour_pressures(plot=False):