Merge pull request #715 from slayoo/freezing

constants refactor (allowing overriding any value through Formulae ctor arg)

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,

PySDM/attributes/impl/mapper.py

@@ -21,19 +21,19 @@
     'volume': lambda _: Volume,
     'dry volume organic': lambda dynamics: (
         make_dummy_attribute_factory('dry volume organic')
-        if 'Condensation' in dynamics and isinstance(
-            dynamics['Condensation'].particulator.formulae.surface_tension,
-            Constant
+        if 'Condensation' in dynamics and (
+            dynamics['Condensation'].particulator.formulae.surface_tension.__name__ ==
+            Constant.__name__
         )
         else DryVolumeOrganic
     ),
     'dry volume': lambda dynamics:
     DryVolumeDynamic if 'AqueousChemistry' in dynamics else DryVolume,
     'dry volume organic fraction': lambda dynamics: (
         make_dummy_attribute_factory('dry volume organic fraction')
-        if 'Condensation' in dynamics and isinstance(
-            dynamics['Condensation'].particulator.formulae.surface_tension,
-            Constant
+        if 'Condensation' in dynamics and (
+            dynamics['Condensation'].particulator.formulae.surface_tension.__name__ ==
+            Constant.__name__
         )
         else OrganicFraction
     ),

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
         )

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

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

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,40 +66,46 @@ 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
 
             mole_amount_taken += multiplicity[i] * (
                     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

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(

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):

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):

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

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}); 

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