From 5a5ffb911313b2084d4f128985233ebed299762a Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Tue, 28 Dec 2021 22:07:55 +0100 Subject: [PATCH 01/23] constants refactor (allowing overriding any value through Formulae ctor arg) --- PySDM/formulae.py | 101 ++++++++++-------- .../physics/condensation_coordinate/volume.py | 14 +-- .../volume_logarithm.py | 14 +-- PySDM/physics/constants.py | 16 +++ .../diffusion_kinetics/fuchs_sutugin.py | 16 +-- PySDM/physics/diffusion_kinetics/neglect.py | 9 +- PySDM/physics/diffusion_thermics/negclect.py | 6 +- .../diffusion_thermics/tracy_welch_porter.py | 10 +- PySDM/physics/drop_growth/maxwell_mason.py | 6 +- .../bigg_1953.py | 22 ++-- .../niemand_et_al_2012.py | 28 +++-- .../freezing_temperature_spectrum/null.py | 5 +- .../abifm.py | 16 +-- .../constant.py | 11 +- .../heterogeneous_ice_nucleation_rate/null.py | 5 +- PySDM/physics/hydrostatics/default.py | 14 +-- PySDM/physics/hygroscopicity/kappa_koehler.py | 14 +-- .../kappa_koehler_leading_terms.py | 14 +-- PySDM/physics/latent_heat/constant.py | 6 +- PySDM/physics/latent_heat/kirchhoff.py | 6 +- .../particle_advection/explicit_in_space.py | 5 +- .../particle_advection/implicit_in_space.py | 5 +- .../august_roche_magnus.py | 15 +-- .../flatau_walko_cotton.py | 14 +-- .../state_variable_triplet/rhod_thd_qv.py | 13 +-- .../compressed_film_ovadnevaite.py | 17 ++- .../surface_tension/compressed_film_ruehl.py | 14 +-- PySDM/physics/surface_tension/constant.py | 6 +- PySDM/physics/trivia.py | 38 +++---- PySDM/physics/ventilation/neglect.py | 3 +- .../alpert_and_knopf_2016/test_ak16_fig_1.py | 4 +- .../arabas_et_al_2015/test_demo_settings.py | 7 -- .../backends/test_freezing_methods.py | 8 +- .../initialisation/test_r_wet_init.py | 20 ++-- .../test_spectral_discretisation.py | 8 +- tests/unit_tests/physics/test_constants.py | 1 - tests/unit_tests/physics/test_formulae.py | 6 +- .../test_freezing_temperature_spectra.py | 15 +-- 38 files changed, 287 insertions(+), 245 deletions(-) delete mode 100644 tests/smoke_tests/arabas_et_al_2015/test_demo_settings.py diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 2e7369f75..0e740cd80 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -4,24 +4,26 @@ import inspect import re from functools import lru_cache, partial - +from collections import namedtuple +import numbers import numba +import numpy as np 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} - } - ) +def _formula(func=None, constants=None, **kw): + source = "class _:\n" + for line in inspect.getsourcelines(func)[0]: + source += f"{line}\n" + loc = {} + exec( + source.replace(f'def {func.__name__}(const,', f'def {func.__name__}('), + {'const': constants, 'np': np}, + loc + ) return numba.njit( - func, + getattr(loc['_'], func.__name__), **{ **conf.JIT_FLAGS, **{'parallel': False, 'inline': 'always', 'cache': False, **kw} @@ -29,20 +31,20 @@ def _formula(func=None, **kw): ) -def _boost(obj, fastmath): +def _boost(obj, fastmath, constants): 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) + formula = _formula(attr, constants=constants, fastmath=fastmath) setattr(obj, item, formula) - setattr(getattr(obj, item), 'c_inline', partial(_c_inline, fun=formula)) + setattr(getattr(obj, item), 'c_inline', partial(_c_inline, constants=constants, fun=attr)) return obj -def _c_inline(fun, return_type=None, **args): +def _c_inline(fun, return_type=None, constants=None, **args): real_t = 'real_type' return_type = return_type or real_t prae = r"([,+\-*/( ]|^)" @@ -58,23 +60,24 @@ def _c_inline(fun, return_type=None, **args): if stripped.startswith('def '): continue source += stripped - source = source.replace("power(", "pow(") + # 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) + if arg != '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 + ")({const.\\2:" + real_fmt + "})\\3", + "\\1(" + real_t + ")({constants.\\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): +def _pick(value: str, choices: dict, constants: namedtuple): for name, cls in choices.items(): if name == value: - return cls() + return cls(constants) raise ValueError(f"Unknown setting: '{value}'; choices are: {tuple(choices.keys())}") @@ -83,13 +86,14 @@ def _choices(module): @lru_cache() -def _magick(value, module, fastmath): - return _boost(_pick(value, _choices(module)), fastmath) +def _magick(value, module, fastmath, constants): + return _boost(_pick(value, _choices(module), constants), fastmath, constants) class Formulae: def __init__(self, *, - seed: int = const.default_random_seed, + constants: dict = {}, + seed: int = physics.constants.default_random_seed, fastmath: bool = True, condensation_coordinate: str = 'VolumeLogarithm', saturation_vapour_pressure: str = 'FlatauWalkoCotton', @@ -106,40 +110,51 @@ def __init__(self, *, freezing_temperature_spectrum: str = 'Null', heterogeneous_ice_nucleation_rate: str = 'Null' ): + constants_defaults = { + k:getattr(physics.constants, k) + for k in dir(physics.constants) + if isinstance(getattr(physics.constants, k), numbers.Number) + } + constants = namedtuple( + "Constants", + tuple(constants_defaults.keys()) + )( + **{**constants_defaults, **constants} + ) + self.constants = constants self.seed = seed self.fastmath = fastmath - self.trivia = _magick('Trivia', physics.trivia, fastmath) + self.trivia = _magick('Trivia', physics.trivia, fastmath, constants) self.condensation_coordinate = _magick( - condensation_coordinate, physics.condensation_coordinate, fastmath) + condensation_coordinate, physics.condensation_coordinate, fastmath, constants) self.saturation_vapour_pressure = _magick( - saturation_vapour_pressure, physics.saturation_vapour_pressure, fastmath) + saturation_vapour_pressure, physics.saturation_vapour_pressure, fastmath, constants) self.latent_heat = _magick( - latent_heat, physics.latent_heat, fastmath) + latent_heat, physics.latent_heat, fastmath, constants) self.hygroscopicity = _magick( - hygroscopicity, physics.hygroscopicity, fastmath) + hygroscopicity, physics.hygroscopicity, fastmath, constants) self.drop_growth = _magick( - drop_growth, physics.drop_growth, fastmath) + drop_growth, physics.drop_growth, fastmath, constants) self.surface_tension = _magick( - surface_tension, physics.surface_tension, fastmath) + surface_tension, physics.surface_tension, fastmath, constants) self.diffusion_kinetics = _magick( - diffusion_kinetics, physics.diffusion_kinetics, fastmath) + diffusion_kinetics, physics.diffusion_kinetics, fastmath, constants) self.diffusion_thermics = _magick( - diffusion_thermics, physics.diffusion_thermics, fastmath) + diffusion_thermics, physics.diffusion_thermics, fastmath, constants) self.ventilation = _magick( - ventilation, physics.ventilation, fastmath) + ventilation, physics.ventilation, fastmath, constants) self.state_variable_triplet = _magick( - state_variable_triplet, physics.state_variable_triplet, fastmath) + state_variable_triplet, physics.state_variable_triplet, fastmath, constants) self.particle_advection = _magick( - particle_advection, physics.particle_advection, fastmath) + particle_advection, physics.particle_advection, fastmath, constants) self.hydrostatics = _magick( - hydrostatics, physics.hydrostatics, fastmath) + hydrostatics, physics.hydrostatics, fastmath, constants) self.freezing_temperature_spectrum = _magick( - freezing_temperature_spectrum, physics.freezing_temperature_spectrum, fastmath) + freezing_temperature_spectrum, physics.freezing_temperature_spectrum, fastmath, constants) 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) def __str__(self): description = [] @@ -151,9 +166,3 @@ def __str__(self): value = getattr(self, attr).__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() diff --git a/PySDM/physics/condensation_coordinate/volume.py b/PySDM/physics/condensation_coordinate/volume.py index 252197175..ffd684638 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, const): + 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(const, x): return x @staticmethod - def x(volume): + def x(const, volume): return volume diff --git a/PySDM/physics/condensation_coordinate/volume_logarithm.py b/PySDM/physics/condensation_coordinate/volume_logarithm.py index 9f7d5dc78..135e7b8f0 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, const): + 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..c3f491c39 100644 --- a/PySDM/physics/constants.py +++ b/PySDM/physics/constants.py @@ -5,6 +5,8 @@ """ import os import time + +import numpy as np import pint from scipy import constants as sci from chempy import Substance @@ -107,3 +109,17 @@ def convert_to(value, unit): 44 if 'CI' in os.environ # https://en.wikipedia.org/wiki/44_(number) else time.time_ns() ) + +sgm_org = np.nan +delta_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..10b2ba726 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, const): + 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(const, 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..ade6fc1ce 100644 --- a/PySDM/physics/diffusion_kinetics/neglect.py +++ b/PySDM/physics/diffusion_kinetics/neglect.py @@ -4,14 +4,17 @@ class Neglect: + def __init__(self, const): + pass + @staticmethod - def lambdaD(D, T): # pylint: disable=unused-argument + def lambdaD(const, D, T): # pylint: disable=unused-argument return -1 @staticmethod - def lambdaK(T, p): # pylint: disable=unused-argument + def lambdaK(const, T, p): # pylint: disable=unused-argument return -1 @staticmethod - def DK(DK, r, lmbd): # pylint: disable=unused-argument + def DK(const, 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..7d43e3d9c 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, const): + 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..3e0214567 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, const): + 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/drop_growth/maxwell_mason.py b/PySDM/physics/drop_growth/maxwell_mason.py index 6bfe53c17..a13d05ff7 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, const): + 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..0e83048b8 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,27 @@ [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..19fd94f02 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, const): + 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..b74480965 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, const): + 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..b9450c358 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, const): + 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..8eaec3c83 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, const): + 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..146290834 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, const): + 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/latent_heat/constant.py b/PySDM/physics/latent_heat/constant.py index fcafe3c07..22921b30e 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, const): + 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..b6b7e75be 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, const): + 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..d5ced5f4f 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, const): + pass + @staticmethod - def displacement(omega, c_l, c_r): + def displacement(const, 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..17524a9cf 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, const): + pass + @staticmethod - def displacement(omega, c_l, c_r): + def displacement(const, 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..8d7e46447 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, const): + 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(const, 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..58a05bdd4 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, const): + 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 * ( @@ -35,9 +37,9 @@ def ice_Celsius(T): )))))))) @staticmethod - def a_w_ice(T): + def a_w_ice(const, T): return ( - FlatauWalkoCotton.ice_Celsius(T - const.T0) + FlatauWalkoCotton.ice_Celsius(const, T - const.T0) / - FlatauWalkoCotton.pvs_Celsius(T - const.T0) + FlatauWalkoCotton.pvs_Celsius(const, 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..c916ca4f9 100644 --- a/PySDM/physics/state_variable_triplet/rhod_thd_qv.py +++ b/PySDM/physics/state_variable_triplet/rhod_thd_qv.py @@ -2,15 +2,16 @@ 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, const): + pass + # A14 in libcloudph++ 1.0 paper @staticmethod def T(rhod, thd): - return thd * power( + return thd * np.power( rhod * thd / const.p1000 * const.Rd, const.Rd_over_c_pd / (1 - const.Rd_over_c_pd) ) @@ -29,14 +30,14 @@ def dthd_dt(rhod, thd, T, dqv_dt, lv): @staticmethod def th_dry(th_std, qv): - return th_std * power(1 + qv / const.eps, const.Rd / const.c_pd) + return th_std * np.power(1 + qv / const.eps, const.Rd / const.c_pd) @staticmethod def rho_d(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 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..7f8decc4a 100644 --- a/PySDM/physics/surface_tension/compressed_film_ruehl.py +++ b/PySDM/physics/surface_tension/compressed_film_ruehl.py @@ -27,15 +27,15 @@ class CompressedFilmRuehl: 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.nu_org) + assert np.isfinite(const.A0) + assert np.isfinite(const.C0) + assert np.isfinite(const.m_sigma) + assert np.isfinite(const.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 diff --git a/PySDM/physics/surface_tension/constant.py b/PySDM/physics/surface_tension/constant.py index 3412b7eba..c2a76c18c 100644 --- a/PySDM/physics/surface_tension/constant.py +++ b/PySDM/physics/surface_tension/constant.py @@ -10,6 +10,10 @@ class Constant: droplet surface is composed of pure water with constant surface tension `sgm_w`. """ + + def __init__(self, const): + 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..05022e004 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, const): + 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..f1d5cb544 100644 --- a/PySDM/physics/ventilation/neglect.py +++ b/PySDM/physics/ventilation/neglect.py @@ -4,4 +4,5 @@ class Neglect: - pass + def __init__(self, const): + pass 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..94087d83f 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( + 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_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/unit_tests/backends/test_freezing_methods.py b/tests/unit_tests/backends/test_freezing_methods.py index ddab1cc3f..8b02bbad8 100644 --- a/tests/unit_tests/backends/test_freezing_methods.py +++ b/tests/unit_tests/backends/test_freezing_methods.py @@ -33,7 +33,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 +55,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/initialisation/test_r_wet_init.py b/tests/unit_tests/initialisation/test_r_wet_init.py index d5f2aef9c..3b68759a9 100644 --- a/tests/unit_tests/initialisation/test_r_wet_init.py +++ b/tests/unit_tests/initialisation/test_r_wet_init.py @@ -5,25 +5,15 @@ 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.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..2c1fea8d3 100644 --- a/tests/unit_tests/initialisation/test_spectral_discretisation.py +++ b/tests/unit_tests/initialisation/test_spectral_discretisation.py @@ -5,10 +5,7 @@ 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 +14,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_formulae.py b/tests/unit_tests/physics/test_formulae.py index abcc6397e..83b33811c 100644 --- a/tests/unit_tests/physics/test_formulae.py +++ b/tests/unit_tests/physics/test_formulae.py @@ -16,7 +16,7 @@ def test_pvs(): T = 300 * si.kelvins # Act - pvs = sut(T) + pvs = sut(constants, T) # Assert assert pvs.units == si.hectopascals @@ -35,7 +35,7 @@ def test_r_cr(): sgm = constants.sgm_w # Act - r_cr = sut(kp, rd**3, T, sgm) + r_cr = sut(constants, kp, rd**3, T, sgm) # Assert assert r_cr.to_base_units().units == si.metres @@ -51,7 +51,7 @@ def test_lv(): sut = formulae.latent_heat.lv # Act - latent_heat = sut(T) + latent_heat = sut(constants, T) # Assert assert latent_heat.check('[energy]/[mass]') 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 From dcc65df34859e25d719aa3f82114d280a949d172 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 29 Dec 2021 11:26:57 +0100 Subject: [PATCH 02/23] cleanups, bump examples req --- PySDM/formulae.py | 62 +++++++++++++------ .../physics/condensation_coordinate/volume.py | 6 +- .../volume_logarithm.py | 4 +- .../diffusion_kinetics/fuchs_sutugin.py | 4 +- PySDM/physics/diffusion_kinetics/neglect.py | 8 +-- PySDM/physics/diffusion_thermics/negclect.py | 2 +- .../diffusion_thermics/tracy_welch_porter.py | 2 +- PySDM/physics/drop_growth/maxwell_mason.py | 2 +- .../niemand_et_al_2012.py | 12 +++- .../freezing_temperature_spectrum/null.py | 2 +- .../heterogeneous_ice_nucleation_rate/null.py | 2 +- PySDM/physics/hydrostatics/default.py | 2 +- PySDM/physics/hygroscopicity/kappa_koehler.py | 2 +- .../kappa_koehler_leading_terms.py | 2 +- PySDM/physics/latent_heat/constant.py | 2 +- PySDM/physics/latent_heat/kirchhoff.py | 2 +- .../particle_advection/explicit_in_space.py | 4 +- .../particle_advection/implicit_in_space.py | 4 +- .../august_roche_magnus.py | 4 +- .../flatau_walko_cotton.py | 2 +- .../state_variable_triplet/rhod_thd_qv.py | 17 ++--- .../surface_tension/compressed_film_ruehl.py | 2 - PySDM/physics/surface_tension/constant.py | 3 +- PySDM/physics/trivia.py | 2 +- PySDM/physics/ventilation/neglect.py | 2 +- test-time-requirements.txt | 2 +- .../backends/test_freezing_methods.py | 1 - .../test_spectral_discretisation.py | 1 - 28 files changed, 93 insertions(+), 67 deletions(-) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 0e740cd80..c89a90170 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -2,6 +2,7 @@ Logic for enabling common CPU/GPU physics formulae code """ import inspect +from typing import Optional import re from functools import lru_cache, partial from collections import namedtuple @@ -17,8 +18,10 @@ def _formula(func=None, constants=None, **kw): for line in inspect.getsourcelines(func)[0]: source += f"{line}\n" loc = {} - exec( - source.replace(f'def {func.__name__}(const,', f'def {func.__name__}('), + 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 ) @@ -40,7 +43,11 @@ def _boost(obj, fastmath, constants): if callable(attr): formula = _formula(attr, constants=constants, fastmath=fastmath) setattr(obj, item, formula) - setattr(getattr(obj, item), 'c_inline', partial(_c_inline, constants=constants, fun=attr)) + setattr( + getattr(obj, item), + 'c_inline', + partial(_c_inline, constants=constants, fun=attr) + ) return obj @@ -70,6 +77,7 @@ def _c_inline(fun, return_type=None, constants=None, **args): "\\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})' @@ -92,7 +100,7 @@ def _magick(value, module, fastmath, constants): class Formulae: def __init__(self, *, - constants: dict = {}, + constants: Optional[dict] = None, seed: int = physics.constants.default_random_seed, fastmath: bool = True, condensation_coordinate: str = 'VolumeLogarithm', @@ -111,7 +119,7 @@ def __init__(self, *, heterogeneous_ice_nucleation_rate: str = 'Null' ): constants_defaults = { - k:getattr(physics.constants, k) + k: getattr(physics.constants, k) for k in dir(physics.constants) if isinstance(getattr(physics.constants, k), numbers.Number) } @@ -119,7 +127,7 @@ def __init__(self, *, "Constants", tuple(constants_defaults.keys()) )( - **{**constants_defaults, **constants} + **{**constants_defaults, **(constants or {})} ) self.constants = constants self.seed = seed @@ -128,33 +136,47 @@ def __init__(self, *, self.trivia = _magick('Trivia', physics.trivia, fastmath, constants) self.condensation_coordinate = _magick( - condensation_coordinate, physics.condensation_coordinate, fastmath, constants) + condensation_coordinate, + physics.condensation_coordinate, fastmath, constants) self.saturation_vapour_pressure = _magick( - saturation_vapour_pressure, physics.saturation_vapour_pressure, fastmath, constants) + saturation_vapour_pressure, + physics.saturation_vapour_pressure, fastmath, constants) self.latent_heat = _magick( - latent_heat, physics.latent_heat, fastmath, constants) + latent_heat, + physics.latent_heat, fastmath, constants) self.hygroscopicity = _magick( - hygroscopicity, physics.hygroscopicity, fastmath, constants) + hygroscopicity, + physics.hygroscopicity, fastmath, constants) self.drop_growth = _magick( - drop_growth, physics.drop_growth, fastmath, constants) + drop_growth, + physics.drop_growth, fastmath, constants) self.surface_tension = _magick( - surface_tension, physics.surface_tension, fastmath, constants) + surface_tension, + physics.surface_tension, fastmath, constants) self.diffusion_kinetics = _magick( - diffusion_kinetics, physics.diffusion_kinetics, fastmath, constants) + diffusion_kinetics, + physics.diffusion_kinetics, fastmath, constants) self.diffusion_thermics = _magick( - diffusion_thermics, physics.diffusion_thermics, fastmath, constants) + diffusion_thermics, + physics.diffusion_thermics, fastmath, constants) self.ventilation = _magick( - ventilation, physics.ventilation, fastmath, constants) + ventilation, + physics.ventilation, fastmath, constants) self.state_variable_triplet = _magick( - state_variable_triplet, physics.state_variable_triplet, fastmath, constants) + state_variable_triplet, + physics.state_variable_triplet, fastmath, constants) self.particle_advection = _magick( - particle_advection, physics.particle_advection, fastmath, constants) + particle_advection, + physics.particle_advection, fastmath, constants) self.hydrostatics = _magick( - hydrostatics, physics.hydrostatics, fastmath, constants) + hydrostatics, + physics.hydrostatics, fastmath, constants) self.freezing_temperature_spectrum = _magick( - freezing_temperature_spectrum, physics.freezing_temperature_spectrum, fastmath, constants) + freezing_temperature_spectrum, + physics.freezing_temperature_spectrum, fastmath, constants) self.heterogeneous_ice_nucleation_rate = _magick( - heterogeneous_ice_nucleation_rate, physics.heterogeneous_ice_nucleation_rate, fastmath, constants) + heterogeneous_ice_nucleation_rate, + physics.heterogeneous_ice_nucleation_rate, fastmath, constants) def __str__(self): description = [] diff --git a/PySDM/physics/condensation_coordinate/volume.py b/PySDM/physics/condensation_coordinate/volume.py index ffd684638..1b4e8cec0 100644 --- a/PySDM/physics/condensation_coordinate/volume.py +++ b/PySDM/physics/condensation_coordinate/volume.py @@ -5,7 +5,7 @@ class Volume: - def __init__(self, const): + def __init__(self, _): pass @staticmethod @@ -13,9 +13,9 @@ 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(const, x): + def volume(_, x): return x @staticmethod - def x(const, 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 135e7b8f0..88d16370e 100644 --- a/PySDM/physics/condensation_coordinate/volume_logarithm.py +++ b/PySDM/physics/condensation_coordinate/volume_logarithm.py @@ -5,12 +5,12 @@ class VolumeLogarithm: - def __init__(self, const): + def __init__(self, _): pass @staticmethod 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) + return np.exp(-const.TWO_THIRDS * x) * r_dr_dt * 3*np.power(const.PI_4_3, const.TWO_THIRDS) @staticmethod def volume(x): diff --git a/PySDM/physics/diffusion_kinetics/fuchs_sutugin.py b/PySDM/physics/diffusion_kinetics/fuchs_sutugin.py index 10b2ba726..7c657e7aa 100644 --- a/PySDM/physics/diffusion_kinetics/fuchs_sutugin.py +++ b/PySDM/physics/diffusion_kinetics/fuchs_sutugin.py @@ -6,7 +6,7 @@ class FuchsSutugin: - def __init__(self, const): + def __init__(self, _): pass @staticmethod @@ -18,5 +18,5 @@ def lambdaK(const, T, p): return (4. / 5) * const.K0 * T / p / np.sqrt(2 * const.Rd * T) @staticmethod - def DK(const, 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 ade6fc1ce..48187b046 100644 --- a/PySDM/physics/diffusion_kinetics/neglect.py +++ b/PySDM/physics/diffusion_kinetics/neglect.py @@ -4,17 +4,17 @@ class Neglect: - def __init__(self, const): + def __init__(self, _): pass @staticmethod - def lambdaD(const, D, T): # pylint: disable=unused-argument + def lambdaD(_, D, T): # pylint: disable=unused-argument return -1 @staticmethod - def lambdaK(const, T, p): # pylint: disable=unused-argument + def lambdaK(_, T, p): # pylint: disable=unused-argument return -1 @staticmethod - def DK(const, 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 7d43e3d9c..c3c5cfad0 100644 --- a/PySDM/physics/diffusion_thermics/negclect.py +++ b/PySDM/physics/diffusion_thermics/negclect.py @@ -4,7 +4,7 @@ class Neglect: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/diffusion_thermics/tracy_welch_porter.py b/PySDM/physics/diffusion_thermics/tracy_welch_porter.py index 3e0214567..3c382f140 100644 --- a/PySDM/physics/diffusion_thermics/tracy_welch_porter.py +++ b/PySDM/physics/diffusion_thermics/tracy_welch_porter.py @@ -6,7 +6,7 @@ class TracyWelchPorter: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/drop_growth/maxwell_mason.py b/PySDM/physics/drop_growth/maxwell_mason.py index a13d05ff7..1197592b5 100644 --- a/PySDM/physics/drop_growth/maxwell_mason.py +++ b/PySDM/physics/drop_growth/maxwell_mason.py @@ -4,7 +4,7 @@ class MaxwellMason: - def __init__(self, const): + def __init__(self, _): pass @staticmethod 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 0e83048b8..26042b31a 100644 --- a/PySDM/physics/freezing_temperature_spectrum/niemand_et_al_2012.py +++ b/PySDM/physics/freezing_temperature_spectrum/niemand_et_al_2012.py @@ -21,9 +21,17 @@ def pdf(const, T, A_insol): @staticmethod 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)) + 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(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) + 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 19fd94f02..8082b60e4 100644 --- a/PySDM/physics/freezing_temperature_spectrum/null.py +++ b/PySDM/physics/freezing_temperature_spectrum/null.py @@ -5,7 +5,7 @@ class Null: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/heterogeneous_ice_nucleation_rate/null.py b/PySDM/physics/heterogeneous_ice_nucleation_rate/null.py index b74480965..f94abb576 100644 --- a/PySDM/physics/heterogeneous_ice_nucleation_rate/null.py +++ b/PySDM/physics/heterogeneous_ice_nucleation_rate/null.py @@ -5,7 +5,7 @@ class Null: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/hydrostatics/default.py b/PySDM/physics/hydrostatics/default.py index b9450c358..6a6a7b131 100644 --- a/PySDM/physics/hydrostatics/default.py +++ b/PySDM/physics/hydrostatics/default.py @@ -5,7 +5,7 @@ class Default: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/hygroscopicity/kappa_koehler.py b/PySDM/physics/hygroscopicity/kappa_koehler.py index 8eaec3c83..bc0b76f4d 100644 --- a/PySDM/physics/hygroscopicity/kappa_koehler.py +++ b/PySDM/physics/hygroscopicity/kappa_koehler.py @@ -6,7 +6,7 @@ class KappaKoehler: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/hygroscopicity/kappa_koehler_leading_terms.py b/PySDM/physics/hygroscopicity/kappa_koehler_leading_terms.py index 146290834..ed7fb10d8 100644 --- a/PySDM/physics/hygroscopicity/kappa_koehler_leading_terms.py +++ b/PySDM/physics/hygroscopicity/kappa_koehler_leading_terms.py @@ -7,7 +7,7 @@ class KappaKoehlerLeadingTerms: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/latent_heat/constant.py b/PySDM/physics/latent_heat/constant.py index 22921b30e..c70c91de2 100644 --- a/PySDM/physics/latent_heat/constant.py +++ b/PySDM/physics/latent_heat/constant.py @@ -4,7 +4,7 @@ class Constant: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/latent_heat/kirchhoff.py b/PySDM/physics/latent_heat/kirchhoff.py index b6b7e75be..13929d53d 100644 --- a/PySDM/physics/latent_heat/kirchhoff.py +++ b/PySDM/physics/latent_heat/kirchhoff.py @@ -6,7 +6,7 @@ class Kirchhoff: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/particle_advection/explicit_in_space.py b/PySDM/physics/particle_advection/explicit_in_space.py index d5ced5f4f..108e18744 100644 --- a/PySDM/physics/particle_advection/explicit_in_space.py +++ b/PySDM/physics/particle_advection/explicit_in_space.py @@ -4,9 +4,9 @@ class ExplicitInSpace: - def __init__(self, const): + def __init__(self, _): pass @staticmethod - def displacement(const, 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 17524a9cf..8a025b200 100644 --- a/PySDM/physics/particle_advection/implicit_in_space.py +++ b/PySDM/physics/particle_advection/implicit_in_space.py @@ -4,9 +4,9 @@ class ImplicitInSpace: - def __init__(self, const): + def __init__(self, _): pass @staticmethod - def displacement(const, 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 8d7e46447..9677287bb 100644 --- a/PySDM/physics/saturation_vapour_pressure/august_roche_magnus.py +++ b/PySDM/physics/saturation_vapour_pressure/august_roche_magnus.py @@ -7,7 +7,7 @@ class AugustRocheMagnus: - def __init__(self, const): + def __init__(self, _): pass @staticmethod @@ -15,5 +15,5 @@ def pvs_Celsius(const, T): return const.ARM_C1 * np.exp((const.ARM_C2 * T) / (T + const.ARM_C3)) @staticmethod - def ice_Celsius(const, 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 58a05bdd4..e6a2c01b3 100644 --- a/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py +++ b/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py @@ -5,7 +5,7 @@ class FlatauWalkoCotton: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/state_variable_triplet/rhod_thd_qv.py b/PySDM/physics/state_variable_triplet/rhod_thd_qv.py index c916ca4f9..b74718b4d 100644 --- a/PySDM/physics/state_variable_triplet/rhod_thd_qv.py +++ b/PySDM/physics/state_variable_triplet/rhod_thd_qv.py @@ -4,36 +4,37 @@ """ import numpy as np + class RhodThdQv: - def __init__(self, const): + def __init__(self, _): pass # A14 in libcloudph++ 1.0 paper @staticmethod - def T(rhod, thd): + 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): + 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) ) / ( @@ -45,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_ruehl.py b/PySDM/physics/surface_tension/compressed_film_ruehl.py index 7f8decc4a..5457562cf 100644 --- a/PySDM/physics/surface_tension/compressed_film_ruehl.py +++ b/PySDM/physics/surface_tension/compressed_film_ruehl.py @@ -5,7 +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 @@ -26,7 +25,6 @@ class CompressedFilmRuehl: the surface concentration to the surface tension. For the compressed film model it is linear, with slope `m_sigma`. """ - @staticmethod def __init__(self, const): assert np.isfinite(const.nu_org) assert np.isfinite(const.A0) diff --git a/PySDM/physics/surface_tension/constant.py b/PySDM/physics/surface_tension/constant.py index c2a76c18c..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: @@ -11,7 +10,7 @@ class Constant: tension `sgm_w`. """ - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/trivia.py b/PySDM/physics/trivia.py index 05022e004..f423b7a29 100644 --- a/PySDM/physics/trivia.py +++ b/PySDM/physics/trivia.py @@ -5,7 +5,7 @@ class Trivia: - def __init__(self, const): + def __init__(self, _): pass @staticmethod diff --git a/PySDM/physics/ventilation/neglect.py b/PySDM/physics/ventilation/neglect.py index f1d5cb544..2e779e4a1 100644 --- a/PySDM/physics/ventilation/neglect.py +++ b/PySDM/physics/ventilation/neglect.py @@ -4,5 +4,5 @@ class Neglect: - def __init__(self, const): + def __init__(self, _): pass diff --git a/test-time-requirements.txt b/test-time-requirements.txt index cf49cf5c5..53adfe94f 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@fed7024#egg=PySDM-examples diff --git a/tests/unit_tests/backends/test_freezing_methods.py b/tests/unit_tests/backends/test_freezing_methods.py index 8b02bbad8..49e300d7d 100644 --- a/tests/unit_tests/backends/test_freezing_methods.py +++ b/tests/unit_tests/backends/test_freezing_methods.py @@ -2,7 +2,6 @@ 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 import Builder, Formulae from PySDM.backends import CPU from PySDM.environments import Box diff --git a/tests/unit_tests/initialisation/test_spectral_discretisation.py b/tests/unit_tests/initialisation/test_spectral_discretisation.py index 2c1fea8d3..e674ed2c7 100644 --- a/tests/unit_tests/initialisation/test_spectral_discretisation.py +++ b/tests/unit_tests/initialisation/test_spectral_discretisation.py @@ -4,7 +4,6 @@ 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 m_mode = .5e-5 From 690df1f2162f945ab5cf8c21d59423c3b2498e1f Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 29 Dec 2021 12:17:51 +0100 Subject: [PATCH 03/23] c_inline fix --- PySDM/formulae.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index c89a90170..f56300af8 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -70,7 +70,7 @@ def _c_inline(fun, return_type=None, constants=None, **args): # source = source.replace("power(", "pow(") source = re.sub("^return ", "", source) for arg in inspect.signature(fun).parameters: - if arg != 'const': + 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}", From 6feba1e85df5cdfeccf7cdeded5ac024f045d633 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 29 Dec 2021 14:36:22 +0100 Subject: [PATCH 04/23] fixes around surface tension --- PySDM/physics/constants.py | 6 ++++ .../surface_tension/compressed_film_ruehl.py | 34 ++++++++----------- test-time-requirements.txt | 2 +- .../smoke_tests/lowe_et_al_2019/test_fig_1.py | 21 +++++++++--- 4 files changed, 37 insertions(+), 26 deletions(-) diff --git a/PySDM/physics/constants.py b/PySDM/physics/constants.py index c3f491c39..3801682c4 100644 --- a/PySDM/physics/constants.py +++ b/PySDM/physics/constants.py @@ -113,6 +113,12 @@ def convert_to(value, unit): 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 diff --git a/PySDM/physics/surface_tension/compressed_film_ruehl.py b/PySDM/physics/surface_tension/compressed_film_ruehl.py index 5457562cf..55cec9273 100644 --- a/PySDM/physics/surface_tension/compressed_film_ruehl.py +++ b/PySDM/physics/surface_tension/compressed_film_ruehl.py @@ -6,12 +6,6 @@ from scipy import constants as sci from scipy import optimize -nu_org = np.nan -A0 = np.nan -C0 = np.nan -m_sigma = np.nan -sgm_min = np.nan - class CompressedFilmRuehl: """ @@ -26,30 +20,30 @@ class CompressedFilmRuehl: is linear, with slope `m_sigma`. """ def __init__(self, const): - assert np.isfinite(const.nu_org) - assert np.isfinite(const.A0) - assert np.isfinite(const.C0) - assert np.isfinite(const.m_sigma) - assert np.isfinite(const.sgm_min) + 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(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/test-time-requirements.txt b/test-time-requirements.txt index 53adfe94f..3921eb3ba 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@fed7024#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@08ccf4e#egg=PySDM-examples 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..d85358e6c 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 @@ -5,7 +5,6 @@ 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 .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( From 05e6abcaa04883a02ac861d83ca96cd802a7543d Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 29 Dec 2021 15:01:13 +0100 Subject: [PATCH 05/23] fixes, fixes, fixes --- README.md | 11 +++++------ test-time-requirements.txt | 2 +- .../alpert_and_knopf_2016/test_ak16_fig_1.py | 2 +- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 1be872931..ef2245426 100644 --- a/README.md +++ b/README.md @@ -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({ ... @@ -520,7 +520,6 @@ 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 3921eb3ba..312ab31ea 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@08ccf4e#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@3d9c2ff#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 94087d83f..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 @@ -33,7 +33,7 @@ def test_ak16_fig_1(multiplicity, plot=False): n_sd = int(n_sd) data, _ = simulation( - J_HET=1e3 / si.cm ** 2 / si.s, + 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 From b6bf5be498151b3ae0d16ac458fc653cf10aab1a Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 29 Dec 2021 19:41:00 +0100 Subject: [PATCH 06/23] fixing test_freezing --- .../smoke_tests/arabas_et_al_2015/test_freezing.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) 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..80faeafc7 100644 --- a/tests/smoke_tests/arabas_et_al_2015/test_freezing.py +++ b/tests/smoke_tests/arabas_et_al_2015/test_freezing.py @@ -9,12 +9,6 @@ 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 +22,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) From a6c183ea5b3fe3c8995f65447dcb142f6285450d Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 29 Dec 2021 19:42:41 +0100 Subject: [PATCH 07/23] temporary debug info to help with matlab issue --- PySDM/formulae.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index f56300af8..981b71ff3 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -20,6 +20,8 @@ def _formula(func=None, constants=None, **kw): loc = {} for arg_name in ('_', 'const'): source = source.replace(f'def {func.__name__}({arg_name},', f'def {func.__name__}(') + if func.__name__ == 'volume': + print(" TEMP ", source) exec( # pylint:disable=exec-used source, {'const': constants, 'np': np}, From ecd2bf7ce9d9e1883bf00cac04e3f0a5ac387bb3 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 29 Dec 2021 19:50:58 +0100 Subject: [PATCH 08/23] cleanup FlatauWalkoCotton class --- .../saturation_vapour_pressure/flatau_walko_cotton.py | 7 ------- test-time-requirements.txt | 2 +- 2 files changed, 1 insertion(+), 8 deletions(-) diff --git a/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py b/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py index e6a2c01b3..3371bf8ce 100644 --- a/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py +++ b/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py @@ -36,10 +36,3 @@ def ice_Celsius(const, T): const.FWC_I8 )))))))) - @staticmethod - def a_w_ice(const, T): - return ( - FlatauWalkoCotton.ice_Celsius(const, T - const.T0) - / - FlatauWalkoCotton.pvs_Celsius(const, T - const.T0) - ) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 312ab31ea..cb1686b2c 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@3d9c2ff#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@7171815#egg=PySDM-examples From 4404ff2ab242e27ce65c2b4ab35b323baf7ac977 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Wed, 29 Dec 2021 14:03:52 -0600 Subject: [PATCH 09/23] fixing Formulae boost logic to be usable from Matlab + pyint cleanups --- PySDM/formulae.py | 9 +++++---- .../saturation_vapour_pressure/flatau_walko_cotton.py | 1 - README.md | 4 ++-- tests/smoke_tests/arabas_et_al_2015/test_freezing.py | 2 -- 4 files changed, 7 insertions(+), 9 deletions(-) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 981b71ff3..59aededa0 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -3,6 +3,7 @@ """ import inspect from typing import Optional +from types import SimpleNamespace import re from functools import lru_cache, partial from collections import namedtuple @@ -20,8 +21,6 @@ def _formula(func=None, constants=None, **kw): loc = {} for arg_name in ('_', 'const'): source = source.replace(f'def {func.__name__}({arg_name},', f'def {func.__name__}(') - if func.__name__ == 'volume': - print(" TEMP ", source) exec( # pylint:disable=exec-used source, {'const': constants, 'np': np}, @@ -38,18 +37,20 @@ def _formula(func=None, constants=None, **kw): def _boost(obj, fastmath, constants): if not physics.impl.flag.DIMENSIONAL_ANALYSIS: + formulae = {} for item in dir(obj): if item.startswith('__'): continue attr = getattr(obj, item) if callable(attr): formula = _formula(attr, constants=constants, fastmath=fastmath) - setattr(obj, item, formula) setattr( - getattr(obj, item), + formula, 'c_inline', partial(_c_inline, constants=constants, fun=attr) ) + formulae[attr.__name__] = formula + return SimpleNamespace(**formulae) return obj diff --git a/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py b/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py index 3371bf8ce..f65ce8ecf 100644 --- a/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py +++ b/PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py @@ -35,4 +35,3 @@ def ice_Celsius(const, T): const.FWC_I7 + T * const.FWC_I8 )))))))) - diff --git a/README.md b/README.md index ef2245426..bfa261ea7 100644 --- a/README.md +++ b/README.md @@ -272,7 +272,7 @@ savefig("plot.svg") rho_w = py.importlib.import_module('PySDM.physics.constants').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(... @@ -513,7 +513,7 @@ for pykey = py.list(keys(particulator.products)) end i=i+1; end -saveas(gcf, "parcel.png") +saveas(gcf, "parcel.png"); ```
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 80faeafc7..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,8 +5,6 @@ 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 From b9fc820f1ae62531c9b576bbedd386459f0fb09b Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 30 Dec 2021 01:43:35 -0600 Subject: [PATCH 10/23] saner dimensional analysis handling within Formulae, adding __name__ field to the boosted simple namespace (and use it in attribute mapper) --- PySDM/attributes/impl/mapper.py | 8 +++---- PySDM/formulae.py | 37 +++++++++++++++++++-------------- 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/PySDM/attributes/impl/mapper.py b/PySDM/attributes/impl/mapper.py index 987571160..4ab399c4b 100644 --- a/PySDM/attributes/impl/mapper.py +++ b/PySDM/attributes/impl/mapper.py @@ -21,8 +21,8 @@ '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, + if 'Condensation' in dynamics and ( + dynamics['Condensation'].particulator.formulae.surface_tension.__name__ == Constant ) else DryVolumeOrganic @@ -31,8 +31,8 @@ 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, + if 'Condensation' in dynamics and ( + dynamics['Condensation'].particulator.formulae.surface_tension.__name__ == Constant ) else OrganicFraction diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 59aededa0..306c8e2bf 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -15,6 +15,13 @@ def _formula(func=None, constants=None, **kw): + if physics.impl.flag.DIMENSIONAL_ANALYSIS: + first_param = inspect.signature(func).parameters[0] + if first_param in ('_', 'const'): + return partial(func, **{first_param: constants}) + else: + return func + source = "class _:\n" for line in inspect.getsourcelines(func)[0]: source += f"{line}\n" @@ -36,22 +43,20 @@ def _formula(func=None, constants=None, **kw): def _boost(obj, fastmath, constants): - if not physics.impl.flag.DIMENSIONAL_ANALYSIS: - formulae = {} - for item in dir(obj): - if item.startswith('__'): - continue - attr = getattr(obj, item) - if callable(attr): - formula = _formula(attr, constants=constants, fastmath=fastmath) - setattr( - formula, - 'c_inline', - partial(_c_inline, constants=constants, fun=attr) - ) - formulae[attr.__name__] = formula - return SimpleNamespace(**formulae) - return obj + 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) + 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): From 9fe440e7d676773a117d357a76ca41b2e7340694 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 30 Dec 2021 01:49:51 -0600 Subject: [PATCH 11/23] fixing __name__ comparison in mapper --- PySDM/attributes/impl/mapper.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PySDM/attributes/impl/mapper.py b/PySDM/attributes/impl/mapper.py index 4ab399c4b..e3f9ad6b8 100644 --- a/PySDM/attributes/impl/mapper.py +++ b/PySDM/attributes/impl/mapper.py @@ -23,7 +23,7 @@ make_dummy_attribute_factory('dry volume organic') if 'Condensation' in dynamics and ( dynamics['Condensation'].particulator.formulae.surface_tension.__name__ == - Constant + Constant.__class__.__name__ ) else DryVolumeOrganic ), @@ -33,7 +33,7 @@ make_dummy_attribute_factory('dry volume organic fraction') if 'Condensation' in dynamics and ( dynamics['Condensation'].particulator.formulae.surface_tension.__name__ == - Constant + Constant.__class__.__name__ ) else OrganicFraction ), From 6eb3db2369a2fef9bd9894227296acb52aaa079d Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 30 Dec 2021 03:17:54 -0600 Subject: [PATCH 12/23] properly fixing __name__ comparison in mapper; better Formulae __str__ --- PySDM/attributes/impl/mapper.py | 4 +- PySDM/formulae.py | 191 +++++++++++----------- tests/unit_tests/physics/test_formulae.py | 6 +- 3 files changed, 103 insertions(+), 98 deletions(-) diff --git a/PySDM/attributes/impl/mapper.py b/PySDM/attributes/impl/mapper.py index e3f9ad6b8..da09c6abb 100644 --- a/PySDM/attributes/impl/mapper.py +++ b/PySDM/attributes/impl/mapper.py @@ -23,7 +23,7 @@ make_dummy_attribute_factory('dry volume organic') if 'Condensation' in dynamics and ( dynamics['Condensation'].particulator.formulae.surface_tension.__name__ == - Constant.__class__.__name__ + Constant.__name__ ) else DryVolumeOrganic ), @@ -33,7 +33,7 @@ make_dummy_attribute_factory('dry volume organic fraction') if 'Condensation' in dynamics and ( dynamics['Condensation'].particulator.formulae.surface_tension.__name__ == - Constant.__class__.__name__ + Constant.__name__ ) else OrganicFraction ), diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 306c8e2bf..daafb94df 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -14,9 +14,106 @@ from PySDM.backends.impl_numba import conf +class Formulae: + def __init__(self, *, + constants: Optional[dict] = None, + seed: int = physics.constants.default_random_seed, + fastmath: bool = True, + condensation_coordinate: str = 'VolumeLogarithm', + saturation_vapour_pressure: str = 'FlatauWalkoCotton', + latent_heat: str = 'Kirchhoff', + hygroscopicity: str = 'KappaKoehlerLeadingTerms', + drop_growth: str = 'MaxwellMason', + surface_tension: str = 'Constant', + diffusion_kinetics: str = 'FuchsSutugin', + diffusion_thermics: str = 'Neglect', + ventilation: str = 'Neglect', + state_variable_triplet: str = 'RhodThdQv', + particle_advection: str = 'ImplicitInSpace', + hydrostatics: str = 'Default', + freezing_temperature_spectrum: str = 'Null', + heterogeneous_ice_nucleation_rate: str = 'Null' + ): + constants_defaults = { + k: getattr(physics.constants, k) + for k in dir(physics.constants) + if isinstance(getattr(physics.constants, k), numbers.Number) + } + constants = namedtuple( + "Constants", + tuple(constants_defaults.keys()) + )( + **{**constants_defaults, **(constants or {})} + ) + self.constants = constants + self.seed = seed + self.fastmath = fastmath + + self.trivia = _magick('Trivia', physics.trivia, fastmath, constants) + + self.condensation_coordinate = _magick( + condensation_coordinate, + physics.condensation_coordinate, fastmath, constants) + self.saturation_vapour_pressure = _magick( + saturation_vapour_pressure, + physics.saturation_vapour_pressure, fastmath, constants) + self.latent_heat = _magick( + latent_heat, + physics.latent_heat, fastmath, constants) + self.hygroscopicity = _magick( + hygroscopicity, + physics.hygroscopicity, fastmath, constants) + self.drop_growth = _magick( + drop_growth, + physics.drop_growth, fastmath, constants) + self.surface_tension = _magick( + surface_tension, + physics.surface_tension, fastmath, constants) + self.diffusion_kinetics = _magick( + diffusion_kinetics, + physics.diffusion_kinetics, fastmath, constants) + self.diffusion_thermics = _magick( + diffusion_thermics, + physics.diffusion_thermics, fastmath, constants) + self.ventilation = _magick( + ventilation, + physics.ventilation, fastmath, constants) + self.state_variable_triplet = _magick( + state_variable_triplet, + physics.state_variable_triplet, fastmath, constants) + self.particle_advection = _magick( + particle_advection, + physics.particle_advection, fastmath, constants) + self.hydrostatics = _magick( + hydrostatics, + physics.hydrostatics, fastmath, constants) + self.freezing_temperature_spectrum = _magick( + freezing_temperature_spectrum, + physics.freezing_temperature_spectrum, fastmath, constants) + self.heterogeneous_ice_nucleation_rate = _magick( + heterogeneous_ice_nucleation_rate, + physics.heterogeneous_ice_nucleation_rate, fastmath, constants) + + def __str__(self): + description = [] + for attr in dir(self): + if not attr.startswith('_'): + 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 = attr_value.__class__.__name__ + description.append(f"{attr}: {value}") + return ', '.join(description) + + def _formula(func=None, constants=None, **kw): if physics.impl.flag.DIMENSIONAL_ANALYSIS: - first_param = inspect.signature(func).parameters[0] + first_param = tuple(inspect.signature(func).parameters.keys())[0] if first_param in ('_', 'const'): return partial(func, **{first_param: constants}) else: @@ -104,95 +201,3 @@ def _choices(module): @lru_cache() def _magick(value, module, fastmath, constants): return _boost(_pick(value, _choices(module), constants), fastmath, constants) - - -class Formulae: - def __init__(self, *, - constants: Optional[dict] = None, - seed: int = physics.constants.default_random_seed, - fastmath: bool = True, - condensation_coordinate: str = 'VolumeLogarithm', - saturation_vapour_pressure: str = 'FlatauWalkoCotton', - latent_heat: str = 'Kirchhoff', - hygroscopicity: str = 'KappaKoehlerLeadingTerms', - drop_growth: str = 'MaxwellMason', - surface_tension: str = 'Constant', - diffusion_kinetics: str = 'FuchsSutugin', - diffusion_thermics: str = 'Neglect', - ventilation: str = 'Neglect', - state_variable_triplet: str = 'RhodThdQv', - particle_advection: str = 'ImplicitInSpace', - hydrostatics: str = 'Default', - freezing_temperature_spectrum: str = 'Null', - heterogeneous_ice_nucleation_rate: str = 'Null' - ): - constants_defaults = { - k: getattr(physics.constants, k) - for k in dir(physics.constants) - if isinstance(getattr(physics.constants, k), numbers.Number) - } - constants = namedtuple( - "Constants", - tuple(constants_defaults.keys()) - )( - **{**constants_defaults, **(constants or {})} - ) - self.constants = constants - self.seed = seed - self.fastmath = fastmath - - self.trivia = _magick('Trivia', physics.trivia, fastmath, constants) - - self.condensation_coordinate = _magick( - condensation_coordinate, - physics.condensation_coordinate, fastmath, constants) - self.saturation_vapour_pressure = _magick( - saturation_vapour_pressure, - physics.saturation_vapour_pressure, fastmath, constants) - self.latent_heat = _magick( - latent_heat, - physics.latent_heat, fastmath, constants) - self.hygroscopicity = _magick( - hygroscopicity, - physics.hygroscopicity, fastmath, constants) - self.drop_growth = _magick( - drop_growth, - physics.drop_growth, fastmath, constants) - self.surface_tension = _magick( - surface_tension, - physics.surface_tension, fastmath, constants) - self.diffusion_kinetics = _magick( - diffusion_kinetics, - physics.diffusion_kinetics, fastmath, constants) - self.diffusion_thermics = _magick( - diffusion_thermics, - physics.diffusion_thermics, fastmath, constants) - self.ventilation = _magick( - ventilation, - physics.ventilation, fastmath, constants) - self.state_variable_triplet = _magick( - state_variable_triplet, - physics.state_variable_triplet, fastmath, constants) - self.particle_advection = _magick( - particle_advection, - physics.particle_advection, fastmath, constants) - self.hydrostatics = _magick( - hydrostatics, - physics.hydrostatics, fastmath, constants) - self.freezing_temperature_spectrum = _magick( - freezing_temperature_spectrum, - physics.freezing_temperature_spectrum, fastmath, constants) - self.heterogeneous_ice_nucleation_rate = _magick( - heterogeneous_ice_nucleation_rate, - physics.heterogeneous_ice_nucleation_rate, fastmath, constants) - - 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) - else: - value = getattr(self, attr).__class__.__name__ - description.append(f"{attr}: {value}") - return ', '.join(description) diff --git a/tests/unit_tests/physics/test_formulae.py b/tests/unit_tests/physics/test_formulae.py index 83b33811c..abcc6397e 100644 --- a/tests/unit_tests/physics/test_formulae.py +++ b/tests/unit_tests/physics/test_formulae.py @@ -16,7 +16,7 @@ def test_pvs(): T = 300 * si.kelvins # Act - pvs = sut(constants, T) + pvs = sut(T) # Assert assert pvs.units == si.hectopascals @@ -35,7 +35,7 @@ def test_r_cr(): sgm = constants.sgm_w # Act - r_cr = sut(constants, kp, rd**3, T, sgm) + r_cr = sut(kp, rd**3, T, sgm) # Assert assert r_cr.to_base_units().units == si.metres @@ -51,7 +51,7 @@ def test_lv(): sut = formulae.latent_heat.lv # Act - latent_heat = sut(constants, T) + latent_heat = sut(T) # Assert assert latent_heat.check('[energy]/[mass]') From d15d1038838f332c3f4824d9d4eb3185364cff37 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 30 Dec 2021 03:29:17 -0600 Subject: [PATCH 13/23] pylint fixes --- PySDM/formulae.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index daafb94df..63aae915b 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -116,9 +116,8 @@ def _formula(func=None, constants=None, **kw): first_param = tuple(inspect.signature(func).parameters.keys())[0] if first_param in ('_', 'const'): return partial(func, **{first_param: constants}) - else: - return func - + return func + source = "class _:\n" for line in inspect.getsourcelines(func)[0]: source += f"{line}\n" From 13788a2f4231768b2bd4f3ef9084e3bae85d948c Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Thu, 30 Dec 2021 13:46:54 +0100 Subject: [PATCH 14/23] fixing dimensional analysis handling --- PySDM/formulae.py | 55 ++++++++++++++++++++-------------- PySDM/physics/__init__.py | 1 + PySDM/physics/impl/__init__.py | 1 + 3 files changed, 34 insertions(+), 23 deletions(-) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index 63aae915b..d0c91a9bd 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -10,6 +10,8 @@ import numbers import numba import numpy as np +import pint + from PySDM import physics from PySDM.backends.impl_numba import conf @@ -37,7 +39,7 @@ def __init__(self, *, constants_defaults = { k: getattr(physics.constants, k) for k in dir(physics.constants) - if isinstance(getattr(physics.constants, k), numbers.Number) + if isinstance(getattr(physics.constants, k), (numbers.Number, pint.Quantity)) } constants = namedtuple( "Constants", @@ -48,51 +50,52 @@ def __init__(self, *, self.constants = constants self.seed = seed self.fastmath = fastmath + dimensional_analysis = physics.impl.flag.DIMENSIONAL_ANALYSIS - self.trivia = _magick('Trivia', physics.trivia, fastmath, constants) + self.trivia = _magick('Trivia', physics.trivia, fastmath, constants, dimensional_analysis) self.condensation_coordinate = _magick( condensation_coordinate, - physics.condensation_coordinate, fastmath, constants) + physics.condensation_coordinate, fastmath, constants, dimensional_analysis) self.saturation_vapour_pressure = _magick( saturation_vapour_pressure, - physics.saturation_vapour_pressure, fastmath, constants) + physics.saturation_vapour_pressure, fastmath, constants, dimensional_analysis) self.latent_heat = _magick( latent_heat, - physics.latent_heat, fastmath, constants) + physics.latent_heat, fastmath, constants, dimensional_analysis) self.hygroscopicity = _magick( hygroscopicity, - physics.hygroscopicity, fastmath, constants) + physics.hygroscopicity, fastmath, constants, dimensional_analysis) self.drop_growth = _magick( drop_growth, - physics.drop_growth, fastmath, constants) + physics.drop_growth, fastmath, constants, dimensional_analysis) self.surface_tension = _magick( surface_tension, - physics.surface_tension, fastmath, constants) + physics.surface_tension, fastmath, constants, dimensional_analysis) self.diffusion_kinetics = _magick( diffusion_kinetics, - physics.diffusion_kinetics, fastmath, constants) + physics.diffusion_kinetics, fastmath, constants, dimensional_analysis) self.diffusion_thermics = _magick( diffusion_thermics, - physics.diffusion_thermics, fastmath, constants) + physics.diffusion_thermics, fastmath, constants, dimensional_analysis) self.ventilation = _magick( ventilation, - physics.ventilation, fastmath, constants) + physics.ventilation, fastmath, constants, dimensional_analysis) self.state_variable_triplet = _magick( state_variable_triplet, - physics.state_variable_triplet, fastmath, constants) + physics.state_variable_triplet, fastmath, constants, dimensional_analysis) self.particle_advection = _magick( particle_advection, - physics.particle_advection, fastmath, constants) + physics.particle_advection, fastmath, constants, dimensional_analysis) self.hydrostatics = _magick( hydrostatics, - physics.hydrostatics, fastmath, constants) + physics.hydrostatics, fastmath, constants, dimensional_analysis) self.freezing_temperature_spectrum = _magick( freezing_temperature_spectrum, - physics.freezing_temperature_spectrum, fastmath, constants) + physics.freezing_temperature_spectrum, fastmath, constants, dimensional_analysis) self.heterogeneous_ice_nucleation_rate = _magick( heterogeneous_ice_nucleation_rate, - physics.heterogeneous_ice_nucleation_rate, fastmath, constants) + physics.heterogeneous_ice_nucleation_rate, fastmath, constants, dimensional_analysis) def __str__(self): description = [] @@ -111,11 +114,11 @@ def __str__(self): return ', '.join(description) -def _formula(func=None, constants=None, **kw): - if physics.impl.flag.DIMENSIONAL_ANALYSIS: +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, **{first_param: constants}) + return partial(func, constants) return func source = "class _:\n" @@ -138,14 +141,15 @@ def _formula(func=None, constants=None, **kw): ) -def _boost(obj, fastmath, constants): +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) + formula = _formula(attr, constants=constants, fastmath=fastmath, + dimensional_analysis=dimensional_analysis) setattr( formula, 'c_inline', @@ -198,5 +202,10 @@ def _choices(module): @lru_cache() -def _magick(value, module, fastmath, constants): - return _boost(_pick(value, _choices(module), constants), fastmath, constants) +def _magick(value, module, fastmath, constants, dimensional_analysis): + return _boost( + _pick(value, _choices(module), constants), + fastmath, + constants, + dimensional_analysis + ) diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 7c4400b33..5d507fe1e 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -6,3 +6,4 @@ ventilation, state_variable_triplet, trivia, particle_advection, hydrostatics, freezing_temperature_spectrum, heterogeneous_ice_nucleation_rate) from .constants import si +from . import impl 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 From aeb30266c1577b0559a471d6260e188efc69dc13 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 31 Dec 2021 05:00:13 -0600 Subject: [PATCH 15/23] fixing c_inline to work back on GPU; bump examples req --- PySDM/formulae.py | 3 ++- test-time-requirements.txt | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/PySDM/formulae.py b/PySDM/formulae.py index d0c91a9bd..cfc3253a8 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -175,7 +175,8 @@ def _c_inline(fun, return_type=None, constants=None, **args): if stripped.startswith('def '): continue source += stripped - # source = source.replace("power(", "pow(") + 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'): diff --git a/test-time-requirements.txt b/test-time-requirements.txt index cb1686b2c..2864e75f4 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@7171815#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@a045305#egg=PySDM-examples From 06b36fb8e97d41bad54c4aad421ec85e5561aedf Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 31 Dec 2021 12:21:41 +0100 Subject: [PATCH 16/23] less constants imports --- PySDM/attributes/chemistry/acidity.py | 6 ++---- .../physics/critical_supersaturation.py | 9 +++------ PySDM/attributes/physics/dry_radius.py | 3 +-- PySDM/attributes/physics/radius.py | 3 +-- .../impl_numba/methods/condensation_methods.py | 18 +++++++++--------- .../impl_numba/methods/freezing_methods.py | 2 +- .../impl_numba/methods/physics_methods.py | 2 +- .../methods/condensation_methods.py | 2 +- PySDM/environments/parcel.py | 7 ++----- PySDM/initialisation/equilibrate_wet_radii.py | 2 +- PySDM/products/freezing/ice_water_content.py | 3 +-- PySDM/products/impl/product.py | 5 ++--- PySDM/products/size_spectral/mean_radius.py | 3 +-- .../size_spectral/water_mixing_ratio.py | 3 +-- 14 files changed, 27 insertions(+), 41 deletions(-) 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/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/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_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/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/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/products/freezing/ice_water_content.py b/PySDM/products/freezing/ice_water_content.py index b9266e94a..fa6e9b961 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.const.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..4c2e65caa 100644 --- a/PySDM/products/impl/product.py +++ b/PySDM/products/impl/product.py @@ -7,7 +7,6 @@ import inspect import numpy as np import pint -from PySDM.physics import constants as const _UNIT_REGISTRY = pint.UnitRegistry() _CAMEL_CASE_PATTERN = re.compile(r'[A-Z]?[a-z]+|[A-Z]+(?![^A-Z])') @@ -35,8 +34,8 @@ def register(self, builder): def _download_to_buffer(self, storage): storage.download(self.buffer.ravel()) - @staticmethod - def _parse_unit(unit: str): + def _parse_unit(self, unit: str): + const = self.formulae.constants if unit in ('%', 'percent'): return .01 * _UNIT_REGISTRY.dimensionless if unit in ('PPB', 'ppb'): 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 From 43be296c8a7eec0c74dfcbae4bbbf15d27235522 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 31 Dec 2021 12:27:28 +0100 Subject: [PATCH 17/23] fix attr_unit check --- PySDM/products/impl/spectrum_moment_product.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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, From c4b6850f393b4ba03b06012dcb6f967d1b7167f2 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 31 Dec 2021 15:59:05 +0100 Subject: [PATCH 18/23] introducing constants_defaults. closes #713 --- .../impl_numba/methods/chemistry_methods.py | 29 +++-- PySDM/backends/impl_numba/test_helpers/bdf.py | 14 +-- PySDM/dynamics/aqueous_chemistry.py | 6 +- PySDM/formulae.py | 6 +- PySDM/physics/__init__.py | 1 + PySDM/physics/aqueous_chemistry/support.py | 65 ++++++----- PySDM/physics/constants.py | 103 ++---------------- PySDM/physics/constants_defaults.py | 97 +++++++++++++++++ PySDM/physics/dimensional_analysis.py | 4 +- .../aqueous_mole_fraction.py | 3 +- .../gaseous_mole_fraction.py | 4 +- .../displacement/surface_precipitation.py | 4 +- PySDM/products/freezing/ice_water_content.py | 2 +- PySDM/products/impl/product.py | 11 +- .../test_conservation.py | 4 +- .../arabas_and_shima_2017/test_event_rates.py | 5 +- .../test_initialisation.py | 14 +-- .../test_ionic_strength.py | 39 +++---- .../kreidenweis_et_al_2003/test_table_3.py | 5 +- tests/unit_tests/attributes/test_acidity.py | 10 +- .../backends/test_freezing_methods.py | 2 +- tests/unit_tests/backends/test_oxidation.py | 19 ++-- .../initialisation/test_r_wet_init.py | 2 +- .../physics/test_dimensional_analysis.py | 10 +- tests/unit_tests/physics/test_formulae.py | 11 +- .../test_saturation_vapour_pressure.py | 2 +- 26 files changed, 255 insertions(+), 217 deletions(-) create mode 100644 PySDM/physics/constants_defaults.py 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/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/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/formulae.py b/PySDM/formulae.py index cfc3253a8..07e11b98e 100644 --- a/PySDM/formulae.py +++ b/PySDM/formulae.py @@ -37,9 +37,9 @@ def __init__(self, *, heterogeneous_ice_nucleation_rate: str = 'Null' ): constants_defaults = { - k: getattr(physics.constants, k) - for k in dir(physics.constants) - if isinstance(getattr(physics.constants, k), (numbers.Number, pint.Quantity)) + 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", diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 5d507fe1e..9c99e2733 100644 --- a/PySDM/physics/__init__.py +++ b/PySDM/physics/__init__.py @@ -7,3 +7,4 @@ 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/constants.py b/PySDM/physics/constants.py index 3801682c4..2cd40baca 100644 --- a/PySDM/physics/constants.py +++ b/PySDM/physics/constants.py @@ -1,15 +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 numpy as np 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 @@ -28,104 +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() -) - -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/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/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/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 fa6e9b961..df7f9f540 100644 --- a/PySDM/products/freezing/ice_water_content.py +++ b/PySDM/products/freezing/ice_water_content.py @@ -17,7 +17,7 @@ def _impl(self, **kwargs): self._download_moment_to_buffer('volume', rank=0, filter_range=(-np.inf, 0)) conc = self.buffer - result[:] *= -self.formulae.constants.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 4c2e65caa..e9d9098ae 100644 --- a/PySDM/products/impl/product.py +++ b/PySDM/products/impl/product.py @@ -7,6 +7,7 @@ import inspect import numpy as np import pint +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])') @@ -34,16 +35,16 @@ def register(self, builder): def _download_to_buffer(self, storage): storage.download(self.buffer.ravel()) - def _parse_unit(self, unit: str): - const = self.formulae.constants + @staticmethod + 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/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/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/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 49e300d7d..a7cf92669 100644 --- a/tests/unit_tests/backends/test_freezing_methods.py +++ b/tests/unit_tests/backends/test_freezing_methods.py @@ -1,7 +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 import constants_defaults as const from PySDM import Builder, Formulae from PySDM.backends import CPU from PySDM.environments import Box 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 3b68759a9..e46dd5fcb 100644 --- a/tests/unit_tests/initialisation/test_r_wet_init.py +++ b/tests/unit_tests/initialisation/test_r_wet_init.py @@ -4,7 +4,7 @@ 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 import si, constants_defaults as const from PySDM.backends import CPU 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_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): From 497faec0c16147f209cee469ecd3788fbcd88b4b Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 31 Dec 2021 16:01:35 +0100 Subject: [PATCH 19/23] bump examples req --- test-time-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 2864e75f4..cbd96b460 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@a045305#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@099f5c1#egg=PySDM-examples From 21a1487b845494d14a0189eb548eca3a4d897f13 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Fri, 31 Dec 2021 16:02:35 +0100 Subject: [PATCH 20/23] fix --- tests/smoke_tests/lowe_et_al_2019/test_fig_1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 d85358e6c..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,7 +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 import si, constants_defaults as const from .constants import constants assert hasattr(constants, '_pytestfixturefunction') From 263b77dfd4a656264e39d604f8dc3dfe0cef6ce2 Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 2 Jan 2022 00:29:31 +0100 Subject: [PATCH 21/23] bump examples req --- test-time-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index cbd96b460..73c92b8fe 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@099f5c1#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@be6e56d#egg=PySDM-examples From c8b41c9ed675f6256c71e5d5de05066e68755b5a Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 2 Jan 2022 00:39:24 +0100 Subject: [PATCH 22/23] fix rho_w in README.md --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index bfa261ea7..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,7 +269,7 @@ 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)); @@ -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]: From c743528109aa5d3299fb314a9247eb870e2a0fda Mon Sep 17 00:00:00 2001 From: Sylwester Arabas Date: Sun, 2 Jan 2022 12:38:28 +0100 Subject: [PATCH 23/23] bump examples req --- test-time-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test-time-requirements.txt b/test-time-requirements.txt index 73c92b8fe..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@be6e56d#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@1f08b2f#egg=PySDM-examples