diff --git a/PySDM/environments/_moist.py b/PySDM/environments/_moist.py index 82b2163c1..192d9a464 100644 --- a/PySDM/environments/_moist.py +++ b/PySDM/environments/_moist.py @@ -46,6 +46,14 @@ def sync(self): target['rhod'], target['thd'], target['qv'], target['T'], target['p'], target['RH'] ) + if 'a_w_ice' in self.variables: + self.particulator.backend.a_w_ice( + T=target['T'], + p=target['p'], + RH=target['RH'], + qv=target['qv'], + a_w_ice=target['a_w_ice'] + ) self._values["predicted"] = target @abstractmethod diff --git a/PySDM/exporters/netcdf_exporter.py b/PySDM/exporters/netcdf_exporter.py index 4e4bac741..c0a20e4d5 100644 --- a/PySDM/exporters/netcdf_exporter.py +++ b/PySDM/exporters/netcdf_exporter.py @@ -38,7 +38,8 @@ def _create_variables(self, ncdf): # TODO #340 ParticleVolume var for name, instance in self.simulator.products.items(): - assert name not in self.vars + if name in self.vars: + raise AssertionError(f"product ({name}) has same name as one of netCDF dimensions") n_dimensions = len(instance.shape) if n_dimensions == 3: @@ -51,7 +52,6 @@ def _create_variables(self, ncdf): raise NotImplementedError() self.vars[name] = ncdf.createVariable(name, "f", dimensions) self.vars[name].units = instance.unit - self.vars[name].long_name = instance.description def _write_variables(self, i): self.vars["T"][i] = self.settings.output_steps[i] * self.settings.dt diff --git a/PySDM/impl/product.py b/PySDM/impl/product.py deleted file mode 100644 index 781797725..000000000 --- a/PySDM/impl/product.py +++ /dev/null @@ -1,108 +0,0 @@ -import numpy as np -from PySDM.environments._moist import _Moist - - -class Product: - - def __init__(self, name, unit=None, description=None): - self.name = name - self.unit = unit - self.description = description - self.shape = None - self.buffer = None - self.particulator = None - self.formulae = None - - def register(self, builder): - self.particulator = builder.particulator - self.formulae = self.particulator.formulae - self.shape = self.particulator.mesh.grid - self.buffer = np.empty(self.particulator.mesh.grid) - - def download_to_buffer(self, storage): - storage.download(self.buffer.ravel()) - - -class MomentProduct(Product): - - def __init__(self, name, unit, description): - super().__init__(name, unit, description) - self.moment_0 = None - self.moments = None - - def register(self, builder): - super().register(builder) - self.moment_0 = self.particulator.Storage.empty( - self.particulator.mesh.n_cell, dtype=float) - self.moments = self.particulator.Storage.empty( - (1, self.particulator.mesh.n_cell), dtype=float) - - def download_moment_to_buffer( - self, attr, rank, - filter_attr='volume', filter_range=(-np.inf, np.inf), - weighting_attribute='volume', weighting_rank=0 - ): - self.particulator.moments( - self.moment_0, self.moments, {attr: (rank,)}, - attr_name=filter_attr, attr_range=filter_range, - weighting_attribute=weighting_attribute, weighting_rank=weighting_rank - ) - if rank == 0: # TODO #217 - self.download_to_buffer(self.moment_0) - else: - self.download_to_buffer(self.moments[0, :]) - - -class SpectrumMomentProduct(Product): - def __init__(self, name, unit, description): - super().__init__(name, unit, description) - self.attr_bins_edges = None - self.moment_0 = None - self.moments = None - - def register(self, builder): - super().register(builder) - self.moment_0 = self.particulator.Storage.empty( - (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) - - def recalculate_spectrum_moment( - self, attr, - rank, filter_attr='volume', - weighting_attribute='volume', weighting_rank=0 - ): - self.particulator.spectrum_moments( - self.moment_0, self.moments, attr, rank, self.attr_bins_edges, - attr_name=filter_attr, - weighting_attribute=weighting_attribute, weighting_rank=weighting_rank - ) - - def download_spectrum_moment_to_buffer(self, rank, bin_number): - if rank == 0: # TODO #217 - self.download_to_buffer(self.moment_0[bin_number, :]) - else: - self.download_to_buffer(self.moments[bin_number, :]) - - -class MoistEnvironmentProduct(Product): - def __init__(self, **args): - super().__init__(**args) - self._name = self.name - self.name = self._name + '_env' - self.environment = None - self.source = None - - def register(self, builder): - super().register(builder) - self.particulator.observers.append(self) - self.environment = builder.particulator.environment - self.source = self.environment[self._name] - - def notify(self): - if isinstance(self.environment, _Moist): - self.source = self.environment.get_predicted(self._name) - - def get(self): - self.download_to_buffer(self.source) - return self.buffer diff --git a/PySDM/products/PartMC/VolumeFractalDimension.py b/PySDM/products/PartMC/VolumeFractalDimension.py deleted file mode 100644 index e650d5be4..000000000 --- a/PySDM/products/PartMC/VolumeFractalDimension.py +++ /dev/null @@ -1,9 +0,0 @@ -from PySDM.impl.product import Product - -class VolumeFractalDimension(Product): - def __init__(self): - super().__init__( - name='PMC_fractal_surface_frac_dim', - unit='1', - description='PartMC: volume fractal dimension' - ) diff --git a/PySDM/products/PartMC/__init__.py b/PySDM/products/PartMC/__init__.py deleted file mode 100644 index e65bcf9d8..000000000 --- a/PySDM/products/PartMC/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from .VolumeFractalDimension import VolumeFractalDimension diff --git a/PySDM/products/__init__.py b/PySDM/products/__init__.py index 00d8d0061..72c3f53f3 100644 --- a/PySDM/products/__init__.py +++ b/PySDM/products/__init__.py @@ -8,30 +8,6 @@ from .condensation import * from .displacement import * from .freezing import * -from .PartMC import * - -from .aerosol_specific_concentration import AerosolSpecificConcentration -from .cloud_droplet_effective_radius import CloudDropletEffectiveRadius -from .dry_air_density import DryAirDensity -from .dry_air_potential_temperature import DryAirPotentialTemperature -from .dynamic_wall_time import DynamicWallTime -from .parcel_displacement import ParcelDisplacement -from .particle_mean_radius import ParticleMeanRadius -from .particles_concentration import ( - ParticlesConcentration, AerosolConcentration, CloudDropletConcentration, DrizzleConcentration -) -from .particles_size_spectrum import ( - ParticlesSizeSpectrum, ParticlesDrySizeSpectrum, ParticlesWetSizeSpectrum -) -from .particles_volume_spectrum import ParticlesVolumeSpectrum -from .pressure import Pressure -from .relative_humidity import RelativeHumidity -from .super_droplet_count import SuperDropletCount -from .temperature import Temperature -from .time import Time -from .timers import CPUTime, WallTime -from .total_dry_mass_mixing_ratio import TotalDryMassMixingRatio -from .total_particle_concentration import TotalParticleConcentration -from .total_particle_specific_concentration import TotalParticleSpecificConcentration -from .water_mixing_ratio import WaterMixingRatio -from .water_vapour_mixing_ratio import WaterVapourMixingRatio +from .ambient_thermodynamics import * +from .housekeeping import * +from .size_spectral import * diff --git a/PySDM/products/aerosol_specific_concentration.py b/PySDM/products/aerosol_specific_concentration.py deleted file mode 100644 index b33e99c59..000000000 --- a/PySDM/products/aerosol_specific_concentration.py +++ /dev/null @@ -1,25 +0,0 @@ -from PySDM.physics import constants as const -from PySDM.impl.product import MomentProduct - - -class AerosolSpecificConcentration(MomentProduct): - - def __init__(self, radius_threshold): - self.radius_threshold = radius_threshold - super().__init__( - name='n_a_mg', - unit='mg-1', - description='Aerosol specific concentration' - ) - - def get(self): - self.download_moment_to_buffer( - 'volume', rank=0, - filter_range=[0, self.formulae.trivia.volume(self.radius_threshold)] - ) - result = self.buffer.copy() # TODO #217 - self.download_to_buffer(self.particulator.environment['rhod']) - result[:] /= self.particulator.mesh.dv - result[:] /= self.buffer - const.convert_to(result, const.si.milligram**-1) - return result diff --git a/PySDM/products/ambient_thermodynamics/__init__.py b/PySDM/products/ambient_thermodynamics/__init__.py new file mode 100644 index 000000000..8a7133c58 --- /dev/null +++ b/PySDM/products/ambient_thermodynamics/__init__.py @@ -0,0 +1,6 @@ +from .ambient_pressure import AmbientPressure +from .ambient_temperature import AmbientTemperature +from .ambient_relative_humidity import AmbientRelativeHumidity +from .ambient_dry_air_density import AmbientDryAirDensity +from .ambient_dry_air_potential_temperature import AmbientDryAirPotentialTemperature +from .ambient_water_vapour_mixing_ratio import AmbientWaterVapourMixingRatio diff --git a/PySDM/products/ambient_thermodynamics/ambient_dry_air_density.py b/PySDM/products/ambient_thermodynamics/ambient_dry_air_density.py new file mode 100644 index 000000000..3159c87be --- /dev/null +++ b/PySDM/products/ambient_thermodynamics/ambient_dry_air_density.py @@ -0,0 +1,6 @@ +from PySDM.products.impl.moist_environment_product import MoistEnvironmentProduct + + +class AmbientDryAirDensity(MoistEnvironmentProduct): + def __init__(self, name="rhod", unit="kg/m^3", var=None): + super().__init__(name=name, unit=unit, var=var) diff --git a/PySDM/products/ambient_thermodynamics/ambient_dry_air_potential_temperature.py b/PySDM/products/ambient_thermodynamics/ambient_dry_air_potential_temperature.py new file mode 100644 index 000000000..d7a289b85 --- /dev/null +++ b/PySDM/products/ambient_thermodynamics/ambient_dry_air_potential_temperature.py @@ -0,0 +1,6 @@ +from PySDM.products.impl.moist_environment_product import MoistEnvironmentProduct + + +class AmbientDryAirPotentialTemperature(MoistEnvironmentProduct): + def __init__(self, unit='K', name=None, var=None): + super().__init__(unit=unit, name=name, var=var) diff --git a/PySDM/products/ambient_thermodynamics/ambient_pressure.py b/PySDM/products/ambient_thermodynamics/ambient_pressure.py new file mode 100644 index 000000000..90a69aafa --- /dev/null +++ b/PySDM/products/ambient_thermodynamics/ambient_pressure.py @@ -0,0 +1,6 @@ +from PySDM.products.impl.moist_environment_product import MoistEnvironmentProduct + + +class AmbientPressure(MoistEnvironmentProduct): + def __init__(self, name=None, unit="Pa", var=None): + super().__init__(name=name, unit=unit, var=var) diff --git a/PySDM/products/ambient_thermodynamics/ambient_relative_humidity.py b/PySDM/products/ambient_thermodynamics/ambient_relative_humidity.py new file mode 100644 index 000000000..da9eaa74d --- /dev/null +++ b/PySDM/products/ambient_thermodynamics/ambient_relative_humidity.py @@ -0,0 +1,6 @@ +from PySDM.products.impl.moist_environment_product import MoistEnvironmentProduct + + +class AmbientRelativeHumidity(MoistEnvironmentProduct): + def __init__(self, name=None, unit='dimensionless', var=None): + super().__init__(name=name, unit=unit, var=var) diff --git a/PySDM/products/ambient_thermodynamics/ambient_temperature.py b/PySDM/products/ambient_thermodynamics/ambient_temperature.py new file mode 100644 index 000000000..6e2f6ff75 --- /dev/null +++ b/PySDM/products/ambient_thermodynamics/ambient_temperature.py @@ -0,0 +1,6 @@ +from PySDM.products.impl.moist_environment_product import MoistEnvironmentProduct + + +class AmbientTemperature(MoistEnvironmentProduct): + def __init__(self, name=None, unit="K", var=None): + super().__init__(name=name, unit=unit, var=var) diff --git a/PySDM/products/ambient_thermodynamics/ambient_water_vapour_mixing_ratio.py b/PySDM/products/ambient_thermodynamics/ambient_water_vapour_mixing_ratio.py new file mode 100644 index 000000000..4acce86a0 --- /dev/null +++ b/PySDM/products/ambient_thermodynamics/ambient_water_vapour_mixing_ratio.py @@ -0,0 +1,6 @@ +from PySDM.products.impl.moist_environment_product import MoistEnvironmentProduct + + +class AmbientWaterVapourMixingRatio(MoistEnvironmentProduct): + def __init__(self, name=None, unit="dimensionless", var=None): + super().__init__(unit=unit, name=name, var=var) diff --git a/PySDM/products/aqueous_chemistry/__init__.py b/PySDM/products/aqueous_chemistry/__init__.py index 6859c5a0d..88c7647ca 100644 --- a/PySDM/products/aqueous_chemistry/__init__.py +++ b/PySDM/products/aqueous_chemistry/__init__.py @@ -1,4 +1,5 @@ from .gaseous_mole_fraction import GaseousMoleFraction from .aqueous_mole_fraction import AqueousMoleFraction -from .pH import pH +from .acidity import Acidity from .aqueous_mass_spectrum import AqueousMassSpectrum +from .total_dry_mass_mixing_ratio import TotalDryMassMixingRatio diff --git a/PySDM/products/aqueous_chemistry/pH.py b/PySDM/products/aqueous_chemistry/acidity.py similarity index 70% rename from PySDM/products/aqueous_chemistry/pH.py rename to PySDM/products/aqueous_chemistry/acidity.py index 3fff7fb5c..a3e078e2a 100644 --- a/PySDM/products/aqueous_chemistry/pH.py +++ b/PySDM/products/aqueous_chemistry/acidity.py @@ -2,11 +2,17 @@ average pH (averaging after or before taking the logarithm in pH definition) with weighting either by number or volume """ -from PySDM.impl.product import MomentProduct +import numpy as np +from PySDM.products.impl.moment_product import MomentProduct -class pH(MomentProduct): - def __init__(self, radius_range, weighting='volume', attr='conc_H'): +class Acidity(MomentProduct): + def __init__(self, *, + radius_range=(0, np.inf), + weighting='volume', + attr='conc_H', + unit='dimensionless', + name=None): assert attr in ('pH', 'moles_H', 'conc_H') self.attr = attr @@ -18,18 +24,14 @@ def __init__(self, radius_range, weighting='volume', attr='conc_H'): raise NotImplementedError() self.radius_range = radius_range - super().__init__( - name='pH_' + attr + '_' + weighting + '_weighted', - unit='', - description='number-weighted pH' - ) + super().__init__(name=name, unit=unit) def register(self, builder): builder.request_attribute('conc_H') super().register(builder) - def get(self): - self.download_moment_to_buffer( + def _impl(self, **kwargs): + self._download_moment_to_buffer( self.attr, rank=1, filter_range=(self.formulae.trivia.volume(self.radius_range[0]), self.formulae.trivia.volume(self.radius_range[1])), diff --git a/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py b/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py index a2c78277f..225ea1018 100644 --- a/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py +++ b/PySDM/products/aqueous_chemistry/aqueous_mass_spectrum.py @@ -1,18 +1,13 @@ import numpy as np from chempy import Substance -from PySDM.physics.constants import convert_to, si -from PySDM.impl.product import SpectrumMomentProduct +from PySDM.physics.constants import si +from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct from PySDM.physics.aqueous_chemistry.support import AQUEOUS_COMPOUNDS class AqueousMassSpectrum(SpectrumMomentProduct): - - def __init__(self, key, dry_radius_bins_edges, specific=False): - super().__init__( - name=f'dm_{key}/dlog_10(dry diameter){"_spec" if specific else ""}', - unit=f'µg / {"kg" if specific else "m3"} / (unit dD/D)', - description=f'... {key} ...' - ) + def __init__(self, key, dry_radius_bins_edges, specific=False, name=None, unit='kg/m^3'): + super().__init__(name=name, unit=unit) self.key = key self.dry_radius_bins_edges = dry_radius_bins_edges self.molar_mass = Substance.from_formula(AQUEOUS_COMPOUNDS[key][0]).mass * si.g / si.mole @@ -31,24 +26,29 @@ def register(self, builder): self.shape = (*builder.particulator.mesh.grid, len(self.attr_bins_edges) - 1) - def get(self): + def _impl(self, **kwargs): vals = np.empty([self.particulator.mesh.n_cell, len(self.attr_bins_edges) - 1]) - self.recalculate_spectrum_moment( + self._recalculate_spectrum_moment( attr=f'moles_{self.key}', rank=1, filter_attr='dry volume' ) for i in range(vals.shape[1]): - self.download_spectrum_moment_to_buffer(rank=1, bin_number=i) + self._download_spectrum_moment_to_buffer(rank=1, bin_number=i) vals[:, i] = self.buffer.ravel() - self.download_spectrum_moment_to_buffer(rank=0, bin_number=i) + self._download_spectrum_moment_to_buffer(rank=0, bin_number=i) vals[:, i] *= self.buffer.ravel() d_log10_diameter = np.diff(np.log10(2 * self.dry_radius_bins_edges)) vals *= self.molar_mass / d_log10_diameter / self.particulator.mesh.dv - convert_to(vals, si.ug) if self.specific: - self.download_to_buffer(self.particulator.environment['rhod']) + self._download_to_buffer(self.particulator.environment['rhod']) vals[:] /= self.buffer return vals + + +class SpecificAqueousMassSpectrum(AqueousMassSpectrum): + def __init__(self, key, dry_radius_bins_edges, name=None, unit='dimensionless'): + super().__init__(key=key, dry_radius_bins_edges=dry_radius_bins_edges, + name=name, unit=unit, specific=True) diff --git a/PySDM/products/aqueous_chemistry/aqueous_mole_fraction.py b/PySDM/products/aqueous_chemistry/aqueous_mole_fraction.py index 8e43c87fc..04ffcb667 100644 --- a/PySDM/products/aqueous_chemistry/aqueous_mole_fraction.py +++ b/PySDM/products/aqueous_chemistry/aqueous_mole_fraction.py @@ -1,14 +1,12 @@ -from PySDM.impl.product import MomentProduct -from PySDM.physics.constants import convert_to, PPB, Md +from PySDM.products.impl.moment_product import MomentProduct +from PySDM.physics.constants import Md + +DUMMY_SPECIFIC_GRAVITY = 44 class AqueousMoleFraction(MomentProduct): - def __init__(self, key): - super().__init__( - name=f'aq_{key}_ppb', - unit='ppb', - description=f'aqueous {key} mole fraction' - ) + def __init__(self, key, unit='dimensionless', name=None): + super().__init__(unit=unit, name=name) self.aqueous_chemistry = None self.key = key @@ -16,20 +14,22 @@ def register(self, builder): super().register(builder) self.aqueous_chemistry = self.particulator.dynamics['AqueousChemistry'] - def get(self): + def _impl(self, **kwargs): attr = 'moles_' + self.key - self.download_moment_to_buffer(attr, rank=0) + self._download_moment_to_buffer(attr, rank=0) conc = self.buffer.copy() - self.download_moment_to_buffer(attr, rank=1) + self._download_moment_to_buffer(attr, rank=1) tmp = self.buffer.copy() tmp[:] *= conc - tmp[:] *= 44 * Md + tmp[:] *= DUMMY_SPECIFIC_GRAVITY * Md - self.download_to_buffer(self.particulator.environment['rhod']) + self._download_to_buffer(self.particulator.environment['rhod']) tmp[:] /= self.particulator.mesh.dv tmp[:] /= self.buffer - tmp[:] = self.formulae.trivia.mixing_ratio_2_mole_fraction(tmp[:], specific_gravity=44) - convert_to(tmp, PPB) + tmp[:] = self.formulae.trivia.mixing_ratio_2_mole_fraction( + tmp[:], + specific_gravity=DUMMY_SPECIFIC_GRAVITY + ) return tmp diff --git a/PySDM/products/aqueous_chemistry/gaseous_mole_fraction.py b/PySDM/products/aqueous_chemistry/gaseous_mole_fraction.py index 21fb60f3c..b4b36b490 100644 --- a/PySDM/products/aqueous_chemistry/gaseous_mole_fraction.py +++ b/PySDM/products/aqueous_chemistry/gaseous_mole_fraction.py @@ -1,15 +1,10 @@ -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product from PySDM.physics.aqueous_chemistry.support import GASEOUS_COMPOUNDS, SPECIFIC_GRAVITY -from PySDM.physics.constants import convert_to, PPB class GaseousMoleFraction(Product): - def __init__(self, key): - super().__init__( - name=f'gas_{key}_ppb', - unit='ppb', - description=f'gaseous {key} mole fraction' - ) + def __init__(self, key, unit='dimensionless', name=None): + super().__init__(name=name, unit=unit) self.aqueous_chemistry = None self.compound = GASEOUS_COMPOUNDS[key] @@ -17,10 +12,9 @@ def register(self, builder): super().register(builder) self.aqueous_chemistry = self.particulator.dynamics['AqueousChemistry'] - def get(self): + 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] ) - convert_to(tmp, PPB) return tmp diff --git a/PySDM/products/aqueous_chemistry/total_dry_mass_mixing_ratio.py b/PySDM/products/aqueous_chemistry/total_dry_mass_mixing_ratio.py new file mode 100644 index 000000000..e0f683770 --- /dev/null +++ b/PySDM/products/aqueous_chemistry/total_dry_mass_mixing_ratio.py @@ -0,0 +1,19 @@ +import numpy as np +from PySDM.products.impl.moment_product import MomentProduct + + +class TotalDryMassMixingRatio(MomentProduct): + def __init__(self, density, name=None, unit='kg/kg'): + super().__init__(unit=unit, name=name) + self.density = density + + def _impl(self, **kwargs): + self._download_moment_to_buffer('dry volume', rank=1) + self.buffer[:] *= self.density + result = np.copy(self.buffer) + self._download_moment_to_buffer('dry volume', rank=0) + result[:] *= self.buffer + self._download_to_buffer(self.particulator.environment['rhod']) + result[:] /= self.particulator.mesh.dv + result[:] /= self.buffer + return result diff --git a/PySDM/products/coalescence/__init__.py b/PySDM/products/coalescence/__init__.py index 91c84d90d..a1f4a247e 100644 --- a/PySDM/products/coalescence/__init__.py +++ b/PySDM/products/coalescence/__init__.py @@ -1,4 +1,3 @@ from .coalescence_timestep_mean import CoalescenceTimestepMean from .coalescence_timestep_min import CoalescenceTimestepMin -from .collision_rate import CollisionRate -from .collision_rate_deficit import CollisionRateDeficit +from .collision_rates import CollisionRatePerGridbox, CollisionRateDeficitPerGridbox diff --git a/PySDM/products/coalescence/coalescence_timestep_mean.py b/PySDM/products/coalescence/coalescence_timestep_mean.py index 68c41bbd2..f980e804a 100644 --- a/PySDM/products/coalescence/coalescence_timestep_mean.py +++ b/PySDM/products/coalescence/coalescence_timestep_mean.py @@ -1,17 +1,13 @@ import numpy as np import numba from PySDM.backends.impl_numba.conf import JIT_FLAGS -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product class CoalescenceTimestepMean(Product): - def __init__(self): - super().__init__( - name='dt_coal_avg', - unit='s', - description='Coalescence timestep (mean)' - ) + def __init__(self, unit='s', name=None): + super().__init__(unit=unit, name=name) self.count = 0 self.coalescence = None self.range = None @@ -27,8 +23,8 @@ def register(self, builder): def __get_impl(buffer, count, dt): buffer[:] = np.where(buffer[:] > 0, count * dt / buffer[:], np.nan) - def get(self): - self.download_to_buffer(self.coalescence.stats_n_substep) + def _impl(self, **kwargs): + self._download_to_buffer(self.coalescence.stats_n_substep) CoalescenceTimestepMean.__get_impl(self.buffer, self.count, self.particulator.dt) self.coalescence.stats_n_substep[:] = 0 self.count = 0 diff --git a/PySDM/products/coalescence/coalescence_timestep_min.py b/PySDM/products/coalescence/coalescence_timestep_min.py index e04bb980e..96b86ad15 100644 --- a/PySDM/products/coalescence/coalescence_timestep_min.py +++ b/PySDM/products/coalescence/coalescence_timestep_min.py @@ -1,14 +1,9 @@ -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product class CoalescenceTimestepMin(Product): - - def __init__(self): - super().__init__( - name='dt_coal_min', - unit='s', - description='Coalescence timestep (minimal)' - ) + def __init__(self, unit='s', name=None): + super().__init__(unit=unit, name=name) self.coalescence = None self.range = None @@ -17,7 +12,7 @@ def register(self, builder): self.coalescence = self.particulator.dynamics['Coalescence'] self.range = self.coalescence.dt_coal_range - def get(self): - self.download_to_buffer(self.coalescence.stats_dt_min) + def _impl(self, **kwargs): + self._download_to_buffer(self.coalescence.stats_dt_min) self.coalescence.stats_dt_min[:] = self.coalescence.particulator.dt return self.buffer diff --git a/PySDM/products/coalescence/collision_rate.py b/PySDM/products/coalescence/collision_rate.py deleted file mode 100644 index 3e64a9219..000000000 --- a/PySDM/products/coalescence/collision_rate.py +++ /dev/null @@ -1,20 +0,0 @@ -from PySDM.impl.product import Product - - -class CollisionRate(Product): - - def __init__(self): - super().__init__( - name='collision_rate', - description='Collision rate' - ) - self.coalescence = None - - def register(self, builder): - super().register(builder) - self.coalescence = self.particulator.dynamics['Coalescence'] - - def get(self): # TODO #345 take into account NUMBER of substeps (?) - self.download_to_buffer(self.coalescence.collision_rate) - self.coalescence.collision_rate[:] = 0 - return self.buffer diff --git a/PySDM/products/coalescence/collision_rate_deficit.py b/PySDM/products/coalescence/collision_rate_deficit.py deleted file mode 100644 index 75276b76a..000000000 --- a/PySDM/products/coalescence/collision_rate_deficit.py +++ /dev/null @@ -1,20 +0,0 @@ -from PySDM.impl.product import Product - - -class CollisionRateDeficit(Product): - - def __init__(self): - super().__init__( - name='collision_rate_deficit', - description='Collision rate deficit' - ) - self.coalescence = None - - def register(self, builder): - super().register(builder) - self.coalescence = self.particulator.dynamics['Coalescence'] - - def get(self): # TODO #345 take into account NUMBER of substeps (?) - self.download_to_buffer(self.coalescence.collision_rate_deficit) - self.coalescence.collision_rate_deficit[:] = 0 - return self.buffer diff --git a/PySDM/products/coalescence/collision_rates.py b/PySDM/products/coalescence/collision_rates.py new file mode 100644 index 000000000..8609a4c64 --- /dev/null +++ b/PySDM/products/coalescence/collision_rates.py @@ -0,0 +1,33 @@ +from PySDM.products.impl.product import Product + + +class _CollisionRateProduct(Product): + def __init__(self, name, unit, counter): + super().__init__(name=name, unit=unit) + self.timestep_count = 0 + self.counter = counter + + def register(self, builder): + super().register(builder) + self.counter = getattr(self.particulator.dynamics['Coalescence'], self.counter) + self.particulator.observers.append(self) + + def notify(self): + self.timestep_count += 1 + + def _impl(self, **kwargs): + self._download_to_buffer(self.counter) + if self.timestep_count != 0: + self.buffer[:] /= self.timestep_count * self.particulator.dt + self.counter[:] = 0 + return self.buffer + + +class CollisionRateDeficitPerGridbox(_CollisionRateProduct): + def __init__(self, name=None, unit='s^-1'): + super().__init__(name=name, unit=unit, counter='collision_rate_deficit') + + +class CollisionRatePerGridbox(_CollisionRateProduct): + def __init__(self, name=None, unit='s^-1'): + super().__init__(name=name, unit=unit, counter='collision_rate') diff --git a/PySDM/products/condensation/activable_fraction.py b/PySDM/products/condensation/activable_fraction.py index 30aa05ac8..f4eeffba4 100644 --- a/PySDM/products/condensation/activable_fraction.py +++ b/PySDM/products/condensation/activable_fraction.py @@ -1,27 +1,24 @@ -from PySDM.impl.product import MomentProduct +from PySDM.products.impl.moment_product import MomentProduct class ActivableFraction(MomentProduct): - def __init__(self): - super().__init__( - name="activable fraction", - unit="1", - description="" - ) + def __init__(self, unit="dimensionless", name=None): + super().__init__(name=name, unit=unit) def register(self, builder): super().register(builder) builder.request_attribute('critical supersaturation') - def get(self, S_max): - self.download_moment_to_buffer( + def _impl(self, **kwargs): + s_max = kwargs['S_max'] + self._download_moment_to_buffer( 'volume', rank=0, - filter_range=(0, 1 + S_max / 100), + filter_range=(0, 1 + s_max / 100), filter_attr='critical supersaturation' ) frac = self.buffer.copy() - self.download_moment_to_buffer( + self._download_moment_to_buffer( 'volume', rank=0 ) diff --git a/PySDM/products/condensation/condensation_timestep.py b/PySDM/products/condensation/condensation_timestep.py index f694dc722..f7d4bbe5c 100644 --- a/PySDM/products/condensation/condensation_timestep.py +++ b/PySDM/products/condensation/condensation_timestep.py @@ -1,16 +1,12 @@ import numpy as np -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product class _CondensationTimestep(Product): - def __init__(self, extremum, label, reset_value): - super().__init__( - name=f'dt_cond_{label}', - unit='s', - description=f'Condensation timestep ({label})' - ) + def __init__(self, name, unit, extremum, reset_value): + super().__init__(name=name, unit=unit,) self.extremum = extremum self.reset_value = reset_value self.value = None @@ -25,21 +21,21 @@ def register(self, builder): self.value = np.full_like(self.buffer, np.nan) def notify(self): - self.download_to_buffer(self.condensation.counters['n_substeps']) + self._download_to_buffer(self.condensation.counters['n_substeps']) self.buffer[:] = self.condensation.particulator.dt / self.buffer self.value = self.extremum(self.buffer, self.value) - def get(self): + def _impl(self, **kwargs): self.buffer[:] = self.value[:] self.value[:] = self.reset_value return self.buffer class CondensationTimestepMin(_CondensationTimestep): - def __init__(self): - super().__init__(np.minimum, 'min', np.inf) + def __init__(self, name=None, unit='s'): + super().__init__(name=name, unit=unit, extremum=np.minimum, reset_value=np.inf) class CondensationTimestepMax(_CondensationTimestep): - def __init__(self): - super().__init__(np.maximum, 'max', -np.inf) + def __init__(self, name=None, unit='s'): + super().__init__(name=name, unit=unit, extremum=np.maximum, reset_value=-np.inf) diff --git a/PySDM/products/condensation/event_rates.py b/PySDM/products/condensation/event_rates.py index 546411a91..97c8b6a9d 100644 --- a/PySDM/products/condensation/event_rates.py +++ b/PySDM/products/condensation/event_rates.py @@ -1,16 +1,11 @@ import numpy as np -from PySDM.impl.product import Product -from PySDM.physics.constants import convert_to, si +from PySDM.products.impl.product import Product class EventRate(Product): - def __init__(self, what): - super().__init__( - name=what+'_rate', - description=what.capitalize() + ' rate', - unit='s-1 mg-1' - ) + def __init__(self, what, name=None, unit=None): + super().__init__(name=name, unit=unit) self.condensation = None self.timestep_count = 0 self.what = what @@ -24,35 +19,34 @@ def register(self, builder): def notify(self): self.timestep_count += 1 - self.download_to_buffer(self.condensation.counters['n_'+self.what]) + self._download_to_buffer(self.condensation.counters['n_' + self.what]) self.event_count[:] += self.buffer[:] - def get(self): + def _impl(self, **kwargs): if self.timestep_count == 0: return self.event_count self.event_count[:] /= ( self.timestep_count * self.particulator.dt * self.particulator.mesh.dv ) - self.download_to_buffer(self.particulator.environment['rhod']) + self._download_to_buffer(self.particulator.environment['rhod']) self.event_count[:] /= self.buffer[:] self.buffer[:] = self.event_count[:] self.timestep_count = 0 self.event_count[:] = 0 - convert_to(self.buffer, 1/si.mg) return self.buffer class RipeningRate(EventRate): - def __init__(self): - super().__init__('ripening') + def __init__(self, name=None, unit='s^-1 kg^-1'): + super().__init__('ripening', name=name, unit=unit) class ActivatingRate(EventRate): - def __init__(self): - super().__init__('activating') + def __init__(self, name=None, unit='s^-1 kg^-1'): + super().__init__('activating', name=name, unit=unit) class DeactivatingRate(EventRate): - def __init__(self): - super().__init__('deactivating') + def __init__(self, name=None, unit='s^-1 kg^-1'): + super().__init__('deactivating', name=name, unit=unit) diff --git a/PySDM/products/condensation/peak_supersaturation.py b/PySDM/products/condensation/peak_supersaturation.py index a30a3af6f..288d90ffe 100644 --- a/PySDM/products/condensation/peak_supersaturation.py +++ b/PySDM/products/condensation/peak_supersaturation.py @@ -1,14 +1,10 @@ import numpy as np -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product class PeakSupersaturation(Product): - def __init__(self): - super().__init__( - name='S_max', - unit='%', - description='Peak supersaturation' - ) + def __init__(self, unit='dimensionless', name=None): + super().__init__(unit=unit, name=name) self.condensation = None self.RH_max = None @@ -18,11 +14,11 @@ def register(self, builder): self.condensation = self.particulator.dynamics['Condensation'] self.RH_max = np.full_like(self.buffer, np.nan) - def get(self): - self.buffer[:] = (self.RH_max[:] - 1) * 100 - self.RH_max[:] = -100 + def _impl(self, **kwargs): + self.buffer[:] = (self.RH_max[:] - 1) + self.RH_max[:] = -1 return self.buffer def notify(self): - self.download_to_buffer(self.condensation.rh_max) + self._download_to_buffer(self.condensation.rh_max) self.RH_max[:] = np.maximum(self.buffer[:], self.RH_max[:]) diff --git a/PySDM/products/displacement/surface_precipitation.py b/PySDM/products/displacement/surface_precipitation.py index 900f5fb05..ab3734d65 100644 --- a/PySDM/products/displacement/surface_precipitation.py +++ b/PySDM/products/displacement/surface_precipitation.py @@ -1,15 +1,11 @@ -from PySDM.impl.product import Product -from PySDM.physics.constants import rho_w, convert_to, si +from PySDM.products.impl.product import Product +from PySDM.physics.constants import rho_w class SurfacePrecipitation(Product): - def __init__(self): - super().__init__( - name='surf_precip', - unit='mm/day', - description='Surface precipitation' - ) + def __init__(self, name=None, unit='m/s'): + super().__init__(unit=unit, name=name) self.displacement = None self.dv = None self.dz = None @@ -27,13 +23,12 @@ def register(self, builder): self.dv = self.particulator.mesh.dv self.dz = self.particulator.mesh.dz - def get(self) -> float: + def _impl(self, **kwargs) -> float: if self.elapsed_time == 0.: return 0. result = rho_w * self.accumulated_rainfall / self.elapsed_time / (self.dv / self.dz) self._reset_counters() - convert_to(result, si.mm / si.day) return result def notify(self): diff --git a/PySDM/products/dry_air_density.py b/PySDM/products/dry_air_density.py deleted file mode 100644 index e162d74c0..000000000 --- a/PySDM/products/dry_air_density.py +++ /dev/null @@ -1,11 +0,0 @@ -from PySDM.impl.product import MoistEnvironmentProduct - - -class DryAirDensity(MoistEnvironmentProduct): - - def __init__(self): - super().__init__( - description="Dry-air density", - name="rhod", - unit="kg/m^3" - ) diff --git a/PySDM/products/dry_air_potential_temperature.py b/PySDM/products/dry_air_potential_temperature.py deleted file mode 100644 index f9c90a2e3..000000000 --- a/PySDM/products/dry_air_potential_temperature.py +++ /dev/null @@ -1,11 +0,0 @@ -from PySDM.impl.product import MoistEnvironmentProduct - - -class DryAirPotentialTemperature(MoistEnvironmentProduct): - - def __init__(self): - super().__init__( - description="Dry-air potential temperature", - name="thd", - unit="K" - ) diff --git a/PySDM/products/freezing/__init__.py b/PySDM/products/freezing/__init__.py index c4f8f953d..3ac3e9c2d 100644 --- a/PySDM/products/freezing/__init__.py +++ b/PySDM/products/freezing/__init__.py @@ -1,3 +1,3 @@ -from .ice_water_content import IceWaterContent +from .ice_water_content import IceWaterContent, SpecificIceWaterContent from .freezable_specific_concentration import FreezableSpecificConcentration from .total_unfrozen_immersed_surface_area import TotalUnfrozenImmersedSurfaceArea diff --git a/PySDM/products/freezing/freezable_specific_concentration.py b/PySDM/products/freezing/freezable_specific_concentration.py index e4b2fc11f..71bb82330 100644 --- a/PySDM/products/freezing/freezable_specific_concentration.py +++ b/PySDM/products/freezing/freezable_specific_concentration.py @@ -1,16 +1,10 @@ import numpy as np -from PySDM.physics import constants as const -from PySDM.impl.product import SpectrumMomentProduct +from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct class FreezableSpecificConcentration(SpectrumMomentProduct): - - def __init__(self, temperature_bins_edges): - super().__init__( - name='Freezable specific concentration', - unit="mg-1 K-1", - description='Freezable specific concentration' - ) + def __init__(self, temperature_bins_edges, name=None, unit='kg^-1 K^-1'): + super().__init__(name=name, unit=unit) self.attr_bins_edges = temperature_bins_edges def register(self, builder): @@ -20,20 +14,18 @@ def register(self, builder): super().register(builder) self.shape = (*particulator.mesh.grid, len(self.attr_bins_edges) - 1) - def get(self): + def _impl(self, **kwargs): vals = np.empty([self.particulator.mesh.n_cell, len(self.attr_bins_edges) - 1]) - self.recalculate_spectrum_moment(attr='volume', filter_attr='freezing temperature', rank=0) + self._recalculate_spectrum_moment(attr='volume', filter_attr='freezing temperature', rank=0) for i in range(vals.shape[1]): - self.download_spectrum_moment_to_buffer(rank=0, bin_number=i) + self._download_spectrum_moment_to_buffer(rank=0, bin_number=i) vals[:, i] = self.buffer.ravel() - self.download_to_buffer(self.particulator.environment['rhod']) + self._download_to_buffer(self.particulator.environment['rhod']) rhod = self.buffer.ravel() for i in range(len(self.attr_bins_edges) - 1): dT = abs(self.attr_bins_edges[i + 1] - self.attr_bins_edges[i]) vals[:, i] /= rhod * dT * self.particulator.mesh.dv - const.convert_to(vals, const.si.milligram**-1 / const.si.K) - return np.squeeze(vals.reshape(self.shape)) diff --git a/PySDM/products/freezing/ice_water_content.py b/PySDM/products/freezing/ice_water_content.py index 6b0f11e3b..7ad5ec453 100644 --- a/PySDM/products/freezing/ice_water_content.py +++ b/PySDM/products/freezing/ice_water_content.py @@ -1,28 +1,27 @@ import numpy as np -from PySDM.impl.product import MomentProduct +from PySDM.products.impl.moment_product import MomentProduct from PySDM.physics import constants as const class IceWaterContent(MomentProduct): + def __init__(self, unit='kg/kg', name=None, __specific=False): + super().__init__(unit=unit, name=name) + self.specific = __specific - def __init__(self, specific=True): - super().__init__( - name='qi', - unit='g/kg' if specific else 'kg/m3', - description='Ice water mixing ratio' - ) - self.specific = specific - - def get(self): - self.download_moment_to_buffer('volume', rank=0, filter_range=(-np.inf, 0)) + def _impl(self, **kwargs): + self._download_moment_to_buffer('volume', rank=0, filter_range=(-np.inf, 0)) conc = self.buffer.copy() - self.download_moment_to_buffer('volume', rank=1, filter_range=(-np.inf, 0)) + self._download_moment_to_buffer('volume', rank=1, filter_range=(-np.inf, 0)) result = self.buffer.copy() - result[:] *= -const.rho_i * conc / self.particulator.mesh.dv + result[:] *= -const.rho_i * conc / self.particulator.mesh.dv if self.specific: - self.download_to_buffer(self.particulator.environment['rhod']) + self._download_to_buffer(self.particulator.environment['rhod']) result[:] /= self.buffer - const.convert_to(result, const.si.gram / const.si.kilogram) return result + + +class SpecificIceWaterContent(IceWaterContent): + def __init__(self, unit='kg/m^3', name=None, __specific=True): + super().__init__(unit=unit, name=name) diff --git a/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py b/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py index 57167391f..24482ea42 100644 --- a/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py +++ b/PySDM/products/freezing/total_unfrozen_immersed_surface_area.py @@ -1,23 +1,19 @@ import numpy as np -from PySDM.impl.product import MomentProduct +from PySDM.products.impl.moment_product import MomentProduct class TotalUnfrozenImmersedSurfaceArea(MomentProduct): - def __init__(self): - super().__init__( - name='A_tot', - description='total unfrozen immersed surface area', - unit='m2' - ) + def __init__(self, unit='m^2', name=None): + super().__init__(unit=unit, name=name) - def get(self): + def _impl(self, **kwargs): params = { 'attr': 'immersed surface area', 'filter_attr': 'volume', 'filter_range': (0, np.inf) } - self.download_moment_to_buffer(**params, rank=1) + self._download_moment_to_buffer(**params, rank=1) result = np.copy(self.buffer) - self.download_moment_to_buffer(**params, rank=0) + self._download_moment_to_buffer(**params, rank=0) result[:] *= self.buffer return result diff --git a/PySDM/products/housekeeping/__init__.py b/PySDM/products/housekeeping/__init__.py new file mode 100644 index 000000000..e2cbae7cb --- /dev/null +++ b/PySDM/products/housekeeping/__init__.py @@ -0,0 +1,5 @@ +from .time import Time +from .timers import CPUTime, WallTime +from .parcel_displacement import ParcelDisplacement +from .super_droplet_count_per_gridbox import SuperDropletCountPerGridbox +from .dynamic_wall_time import DynamicWallTime diff --git a/PySDM/products/dynamic_wall_time.py b/PySDM/products/housekeeping/dynamic_wall_time.py similarity index 61% rename from PySDM/products/dynamic_wall_time.py rename to PySDM/products/housekeeping/dynamic_wall_time.py index 5579acfab..4d51a762f 100644 --- a/PySDM/products/dynamic_wall_time.py +++ b/PySDM/products/housekeeping/dynamic_wall_time.py @@ -1,13 +1,9 @@ -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product class DynamicWallTime(Product): - def __init__(self, dynamic): - super().__init__( - name=f'{dynamic}_wall_time', - unit='s', - description=f'{dynamic} wall time', - ) + def __init__(self, dynamic, name=None, unit='s'): + super().__init__(name=name, unit=unit) self.value = 0 self.dynamic = dynamic @@ -16,7 +12,7 @@ def register(self, builder): self.particulator.observers.append(self) self.shape = () - def get(self): + def _impl(self, **kwargs): tmp = self.value self.value = 0 return tmp diff --git a/PySDM/products/parcel_displacement.py b/PySDM/products/housekeeping/parcel_displacement.py similarity index 55% rename from PySDM/products/parcel_displacement.py rename to PySDM/products/housekeeping/parcel_displacement.py index a5f509d52..076989538 100644 --- a/PySDM/products/parcel_displacement.py +++ b/PySDM/products/housekeeping/parcel_displacement.py @@ -1,15 +1,11 @@ from PySDM.environments import Parcel -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product class ParcelDisplacement(Product): - def __init__(self): - super().__init__( - description="Parcel displacement", - name="z", - unit="m" - ) + def __init__(self, unit="m", name=None): + super().__init__(unit=unit, name=name) self.environment = None def register(self, builder): @@ -17,6 +13,6 @@ def register(self, builder): assert isinstance(builder.particulator.environment, Parcel) self.environment = builder.particulator.environment - def get(self): - self.download_to_buffer(self.environment['z']) + def _impl(self, **kwargs): + self._download_to_buffer(self.environment['z']) return self.buffer diff --git a/PySDM/products/housekeeping/super_droplet_count_per_gridbox.py b/PySDM/products/housekeeping/super_droplet_count_per_gridbox.py new file mode 100644 index 000000000..800274db0 --- /dev/null +++ b/PySDM/products/housekeeping/super_droplet_count_per_gridbox.py @@ -0,0 +1,13 @@ +from PySDM.products.impl.product import Product + + +class SuperDropletCountPerGridbox(Product): + def __init__(self, unit='dimensionless', name=None): + super().__init__(unit=unit, name=name) + + def _impl(self, **kwargs): + cell_start = self.particulator.attributes.cell_start + n_cell = cell_start.shape[0] - 1 + for i in range(n_cell): + self.buffer.ravel()[i] = cell_start[i + 1] - cell_start[i] + return self.buffer diff --git a/PySDM/products/time.py b/PySDM/products/housekeeping/time.py similarity index 55% rename from PySDM/products/time.py rename to PySDM/products/housekeeping/time.py index aaa079b09..53b539fcf 100644 --- a/PySDM/products/time.py +++ b/PySDM/products/housekeeping/time.py @@ -1,21 +1,17 @@ -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product class Time(Product): - def __init__(self): - super().__init__( - name='t', - unit='s', - description='Time' - ) + def __init__(self, name=None, unit='s'): + super().__init__(name=name, unit=unit) self.t = 0 def register(self, builder): super().register(builder) self.particulator.observers.append(self) - def get(self): + def _impl(self, **kwargs): return self.t def notify(self): diff --git a/PySDM/products/timers.py b/PySDM/products/housekeeping/timers.py similarity index 63% rename from PySDM/products/timers.py rename to PySDM/products/housekeeping/timers.py index b0ad8784d..8b7571bba 100644 --- a/PySDM/products/timers.py +++ b/PySDM/products/housekeeping/timers.py @@ -1,16 +1,12 @@ from abc import abstractmethod import time -from PySDM.impl.product import Product +from PySDM.products.impl.product import Product class _Timer(Product): - def __init__(self, name, description): - super().__init__( - name=name, - unit='s', - description=description - ) + def __init__(self, name, unit): + super().__init__(name=name, unit=unit) self._time = -1 self.reset() @@ -21,7 +17,7 @@ def register(self, builder): super().register(builder) self.shape = () - def get(self) -> float: + def _impl(self, **kwargs) -> float: result = -self._time self.reset() result += self._time @@ -34,8 +30,8 @@ def clock(): class CPUTime(_Timer): - def __init__(self): - super().__init__('cpu_time', 'CPU Time') + def __init__(self, name='CPU Time', unit='s'): + super().__init__(unit=unit, name=name) @staticmethod def clock(): @@ -43,8 +39,8 @@ def clock(): class WallTime(_Timer): - def __init__(self): - super().__init__('wall_time', 'Wall Time') + def __init__(self, name=None, unit='s'): + super().__init__(unit=unit, name=name) @staticmethod def clock(): diff --git a/PySDM/products/impl/__init__.py b/PySDM/products/impl/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/PySDM/products/impl/moist_environment_product.py b/PySDM/products/impl/moist_environment_product.py new file mode 100644 index 000000000..c2b36ff73 --- /dev/null +++ b/PySDM/products/impl/moist_environment_product.py @@ -0,0 +1,24 @@ +from PySDM.environments._moist import _Moist +from PySDM.products.impl.product import Product + + +class MoistEnvironmentProduct(Product): + def __init__(self, *, name, unit, var=None): + super().__init__(name=name, unit=unit) + self.environment = None + self.source = None + self.var = var or name + + def register(self, builder): + super().register(builder) + self.particulator.observers.append(self) + self.environment = builder.particulator.environment + self.source = self.environment[self.var] + + def notify(self): + if isinstance(self.environment, _Moist): + self.source = self.environment.get_predicted(self.var) + + def _impl(self, **kwargs): + self._download_to_buffer(self.source) + return self.buffer diff --git a/PySDM/products/impl/moment_product.py b/PySDM/products/impl/moment_product.py new file mode 100644 index 000000000..fe4e648fd --- /dev/null +++ b/PySDM/products/impl/moment_product.py @@ -0,0 +1,35 @@ +from abc import ABC + +import numpy as np + +from PySDM.products.impl.product import Product + + +class MomentProduct(Product, ABC): + + def __init__(self, name, unit): + super().__init__(name=name, unit=unit) + self.moment_0 = None + self.moments = None + + def register(self, builder): + super().register(builder) + self.moment_0 = self.particulator.Storage.empty( + self.particulator.mesh.n_cell, dtype=float) + self.moments = self.particulator.Storage.empty( + (1, self.particulator.mesh.n_cell), dtype=float) + + def _download_moment_to_buffer( + self, attr, rank, + filter_attr='volume', filter_range=(-np.inf, np.inf), + weighting_attribute='volume', weighting_rank=0 + ): + self.particulator.moments( + self.moment_0, self.moments, {attr: (rank,)}, + attr_name=filter_attr, attr_range=filter_range, + weighting_attribute=weighting_attribute, weighting_rank=weighting_rank + ) + if rank == 0: # TODO #217 + self._download_to_buffer(self.moment_0) + else: + self._download_to_buffer(self.moments[0, :]) diff --git a/PySDM/products/impl/product.py b/PySDM/products/impl/product.py new file mode 100644 index 000000000..59ab4f1e4 --- /dev/null +++ b/PySDM/products/impl/product.py @@ -0,0 +1,89 @@ +from abc import abstractmethod +from typing import Optional +import re +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])') + + +class Product: + def __init__(self, *, unit: str, name: Optional[str] = None): + self.name = name or self._camel_case_to_words(self.__class__.__name__) + + self._unit = self.__parse_unit(unit) + self._base_unit = self._unit.to_base_units().magnitude + self.__check_unit() + + self.shape = None + self.buffer = None + self.particulator = None + self.formulae = None + + def register(self, builder): + self.particulator = builder.particulator + self.formulae = self.particulator.formulae + self.shape = self.particulator.mesh.grid + self.buffer = np.empty(self.particulator.mesh.grid) + + def _download_to_buffer(self, storage): + storage.download(self.buffer.ravel()) + + @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 + if unit in ('PPM', 'ppm'): + return const.PPM * _UNIT_REGISTRY.dimensionless + if unit in ('PPT', 'ppt'): + return const.PPT * _UNIT_REGISTRY.dimensionless + return _UNIT_REGISTRY[unit] + + @staticmethod + def _camel_case_to_words(string: str): + words = _CAMEL_CASE_PATTERN.findall(string) + words = (word if word.isupper() else word.lower() for word in words) + return ' '.join(words) + + def __check_unit(self): + init = inspect.signature(self.__init__) + if 'unit' not in init.parameters: + raise AssertionError(f"method __init__ of class {type(self).__name__}" + f" is expected to have a unit parameter") + + default_unit_arg = init.parameters['unit'].default + + if default_unit_arg is None or str(default_unit_arg).strip() == '': + raise AssertionError(f"unit parameter of {type(self).__name__}.__init__" + f" is expected to have a non-empty default value") + + default_unit = self.__parse_unit(default_unit_arg) + + if default_unit.to_base_units().magnitude != 1: + raise AssertionError(f'default value "{default_unit_arg}"' + f' of {type(self).__name__}.__init__() unit parameter' + f' is not a base SI unit') + + if self._unit.dimensionality != default_unit.dimensionality: + raise AssertionError(f'provided unit ({self._unit}) has different dimensionality' + f' ({self._unit.dimensionality}) than the default one' + f' ({default_unit.dimensionality})' + f' for product {type(self).__name__}') + + @property + def unit(self): + return str(self._unit) + + @abstractmethod + def _impl(self, **kwargs): + raise NotImplementedError() + + def get(self, **kwargs): + result = self._impl(**kwargs) + const.convert_to(result, self._base_unit) + return result diff --git a/PySDM/products/impl/spectrum_moment_product.py b/PySDM/products/impl/spectrum_moment_product.py new file mode 100644 index 000000000..51d9321ec --- /dev/null +++ b/PySDM/products/impl/spectrum_moment_product.py @@ -0,0 +1,34 @@ +from abc import ABC +from PySDM.products.impl.product import Product + + +class SpectrumMomentProduct(ABC, Product): + def __init__(self, name, unit): + super().__init__(name=name, unit=unit) + self.attr_bins_edges = None + self.moment_0 = None + self.moments = None + + def register(self, builder): + super().register(builder) + self.moment_0 = self.particulator.Storage.empty( + (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) + + def _recalculate_spectrum_moment( + self, attr, + rank, filter_attr='volume', + weighting_attribute='volume', weighting_rank=0 + ): + self.particulator.spectrum_moments( + self.moment_0, self.moments, attr, rank, self.attr_bins_edges, + attr_name=filter_attr, + weighting_attribute=weighting_attribute, weighting_rank=weighting_rank + ) + + def _download_spectrum_moment_to_buffer(self, rank, bin_number): + if rank == 0: # TODO #217 + self._download_to_buffer(self.moment_0[bin_number, :]) + else: + self._download_to_buffer(self.moments[bin_number, :]) diff --git a/PySDM/products/particle_mean_radius.py b/PySDM/products/particle_mean_radius.py deleted file mode 100644 index a1d08d955..000000000 --- a/PySDM/products/particle_mean_radius.py +++ /dev/null @@ -1,18 +0,0 @@ -from PySDM.physics import constants as const -from PySDM.impl.product import MomentProduct - - -class ParticleMeanRadius(MomentProduct): - - def __init__(self): - super().__init__( - name='radius_m1', - unit='um', - description='Mean radius' - ) - - def get(self, unit=const.si.micrometre): - self.download_moment_to_buffer('volume', rank=1/3) - self.buffer[:] /= const.PI_4_3 ** (1 / 3) - const.convert_to(self.buffer, unit) - return self.buffer diff --git a/PySDM/products/particles_concentration.py b/PySDM/products/particles_concentration.py deleted file mode 100644 index ca49a564e..000000000 --- a/PySDM/products/particles_concentration.py +++ /dev/null @@ -1,54 +0,0 @@ -import numpy as np -from PySDM.physics import constants as const -from PySDM.impl.product import MomentProduct - - -class ParticlesConcentration(MomentProduct): - - def __init__(self, radius_range=(0, np.inf), specific=False): - self.radius_range = radius_range - self.specific = specific - super().__init__( - name=f'n_part_{"mg" if specific else "cm3"}', - unit='mg-1' if specific else 'cm-3', - description=f'Particle {"specific " if specific else ""}concentration' - ) - - def get(self): - self.download_moment_to_buffer( - 'volume', rank=0, - filter_range=(self.formulae.trivia.volume(self.radius_range[0]), - self.formulae.trivia.volume(self.radius_range[1]))) - self.buffer[:] /= self.particulator.mesh.dv - if self.specific: - result = self.buffer.copy() # TODO #217 - self.download_to_buffer(self.particulator.environment['rhod']) - result[:] /= self.buffer - else: - result = self.buffer - const.convert_to(result, const.si.mg**-1 if self.specific else const.si.centimetre**-3) - return result - - -class AerosolConcentration(ParticlesConcentration): - - def __init__(self, radius_threshold): - super().__init__((0, radius_threshold)) - self.name = 'n_a_cm3' - self.description = 'Aerosol particles concentration' - - -class CloudDropletConcentration(ParticlesConcentration): - - def __init__(self, radius_range): - super().__init__(radius_range) - self.name = 'n_c_cm3' - self.description = 'Cloud droplets concentration' - - -class DrizzleConcentration(ParticlesConcentration): - - def __init__(self, radius_threshold): - super().__init__((radius_threshold, np.inf)) - self.name = 'n_d_cm3' - self.description = 'Drizzle droplets concentration' diff --git a/PySDM/products/pressure.py b/PySDM/products/pressure.py deleted file mode 100644 index 9d1c560a1..000000000 --- a/PySDM/products/pressure.py +++ /dev/null @@ -1,11 +0,0 @@ -from PySDM.impl.product import MoistEnvironmentProduct - - -class Pressure(MoistEnvironmentProduct): - - def __init__(self): - super().__init__( - description="Pressure", - name="p", - unit="Pa" - ) diff --git a/PySDM/products/relative_humidity.py b/PySDM/products/relative_humidity.py deleted file mode 100644 index 61d69754a..000000000 --- a/PySDM/products/relative_humidity.py +++ /dev/null @@ -1,16 +0,0 @@ -from PySDM.impl.product import MoistEnvironmentProduct - - -class RelativeHumidity(MoistEnvironmentProduct): - - def __init__(self): - super().__init__( - description="Relative humidity", - name="RH", - unit="%" - ) - - def get(self): - super().get() - self.buffer *= 100 - return self.buffer diff --git a/PySDM/products/size_spectral/__init__.py b/PySDM/products/size_spectral/__init__.py new file mode 100644 index 000000000..0a1ff2327 --- /dev/null +++ b/PySDM/products/size_spectral/__init__.py @@ -0,0 +1,8 @@ +from .effective_radius import EffectiveRadius +from .mean_radius import MeanRadius +from .particle_size_spectrum import ParticleSizeSpectrumPerMass, ParticleSizeSpectrumPerVolume +from .particles_concentration import ParticleConcentration, ParticleSpecificConcentration +from .particles_volume_spectrum import ParticlesVolumeSpectrum +from .total_particle_concentration import TotalParticleConcentration +from .total_particle_specific_concentration import TotalParticleSpecificConcentration +from .water_mixing_ratio import WaterMixingRatio diff --git a/PySDM/products/cloud_droplet_effective_radius.py b/PySDM/products/size_spectral/effective_radius.py similarity index 63% rename from PySDM/products/cloud_droplet_effective_radius.py rename to PySDM/products/size_spectral/effective_radius.py index 9d79dcd3a..78af5a6d8 100644 --- a/PySDM/products/cloud_droplet_effective_radius.py +++ b/PySDM/products/size_spectral/effective_radius.py @@ -1,39 +1,33 @@ import numpy as np import numba from PySDM.physics import constants as const -from PySDM.impl.product import MomentProduct +from PySDM.products.impl.moment_product import MomentProduct from PySDM.backends.impl_numba.conf import JIT_FLAGS GEOM_FACTOR = const.PI_4_3 ** (-1 / 3) -class CloudDropletEffectiveRadius(MomentProduct): +class EffectiveRadius(MomentProduct): - def __init__(self, radius_range): + def __init__(self, *, radius_range=(0, np.inf), unit='m', name=None): self.radius_range = radius_range - - super().__init__( - name='r_eff', - unit='um', - description='Cloud Droplet Effective Radius' - ) + super().__init__(name=name, unit=unit) @staticmethod @numba.njit(**JIT_FLAGS) def __get_impl(buffer, tmp): buffer[:] = np.where(tmp[:] > 0, buffer[:] * GEOM_FACTOR / tmp[:], np.nan) - def get(self): + def _impl(self, **kwargs): tmp = np.empty_like(self.buffer) - self.download_moment_to_buffer( + self._download_moment_to_buffer( 'volume', rank=2/3, filter_range=(self.formulae.trivia.volume(self.radius_range[0]), self.formulae.trivia.volume(self.radius_range[1]))) tmp[:] = self.buffer[:] - self.download_moment_to_buffer( + self._download_moment_to_buffer( 'volume', rank=1, filter_range=(self.formulae.trivia.volume(self.radius_range[0]), self.formulae.trivia.volume(self.radius_range[1]))) - CloudDropletEffectiveRadius.__get_impl(self.buffer, tmp) - const.convert_to(self.buffer, const.si.micrometre) + EffectiveRadius.__get_impl(self.buffer, tmp) return self.buffer diff --git a/PySDM/products/size_spectral/mean_radius.py b/PySDM/products/size_spectral/mean_radius.py new file mode 100644 index 000000000..04b2eb9ec --- /dev/null +++ b/PySDM/products/size_spectral/mean_radius.py @@ -0,0 +1,12 @@ +from PySDM.physics import constants as const +from PySDM.products.impl.moment_product import MomentProduct + + +class MeanRadius(MomentProduct): + def __init__(self, name=None, unit='m'): + super().__init__(name=name, unit=unit) + + def _impl(self, **kwargs): + self._download_moment_to_buffer('volume', rank=1 / 3) + self.buffer[:] /= const.PI_4_3 ** (1 / 3) + return self.buffer diff --git a/PySDM/products/particles_size_spectrum.py b/PySDM/products/size_spectral/particle_size_spectrum.py similarity index 55% rename from PySDM/products/particles_size_spectrum.py rename to PySDM/products/size_spectral/particle_size_spectrum.py index 1ffbd4538..8af28fc32 100644 --- a/PySDM/products/particles_size_spectrum.py +++ b/PySDM/products/size_spectral/particle_size_spectrum.py @@ -1,19 +1,14 @@ +from abc import ABC import numpy as np -from PySDM.physics import constants as const -from PySDM.impl.product import SpectrumMomentProduct +from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct -class ParticlesSizeSpectrum(SpectrumMomentProduct): - - def __init__(self, radius_bins_edges, name, dry=False, normalise_by_dv=False): +class ParticleSizeSpectrum(SpectrumMomentProduct, ABC): + def __init__(self, radius_bins_edges, name, unit, dry=False, normalise_by_dv=False): self.volume_attr = 'dry volume' if dry else 'volume' self.radius_bins_edges = radius_bins_edges self.normalise_by_dv = normalise_by_dv - super().__init__( - name=name, - unit=f"mg-1 um-1{'' if normalise_by_dv else ' m^3'}", - description='Specific concentration density' - ) + super().__init__(name=name, unit=unit) def register(self, builder): builder.request_attribute(self.volume_attr) @@ -25,48 +20,48 @@ def register(self, builder): self.shape = (*builder.particulator.mesh.grid, len(self.attr_bins_edges) - 1) - def get(self): + def _impl(self, **kwargs): vals = np.empty([self.particulator.mesh.n_cell, len(self.attr_bins_edges) - 1]) - self.recalculate_spectrum_moment( + self._recalculate_spectrum_moment( attr=self.volume_attr, rank=1, filter_attr=self.volume_attr ) for i in range(vals.shape[1]): - self.download_spectrum_moment_to_buffer(rank=0, bin_number=i) + self._download_spectrum_moment_to_buffer(rank=0, bin_number=i) vals[:, i] = self.buffer.ravel() if self.normalise_by_dv: vals[:] /= self.particulator.mesh.dv - self.download_to_buffer(self.particulator.environment['rhod']) + self._download_to_buffer(self.particulator.environment['rhod']) rhod = self.buffer.ravel() for i in range(len(self.attr_bins_edges) - 1): dr = self.formulae.trivia.radius(volume=self.attr_bins_edges[i + 1]) - \ self.formulae.trivia.radius(volume=self.attr_bins_edges[i]) vals[:, i] /= rhod * dr - const.convert_to(vals, const.si.micrometre**-1 * const.si.milligram**-1) - return np.squeeze(vals.reshape(self.shape)) -class ParticlesWetSizeSpectrum(ParticlesSizeSpectrum): - def __init__(self, radius_bins_edges, normalise_by_dv=False): +class ParticleSizeSpectrumPerMass(ParticleSizeSpectrum): + def __init__(self, radius_bins_edges, dry=False, name=None, unit='kg^-1 m^-1'): super().__init__( radius_bins_edges, - dry=False, - normalise_by_dv=normalise_by_dv, - name='Particles Wet Size Spectrum' + dry=dry, + normalise_by_dv=True, + name=name, + unit=unit ) -class ParticlesDrySizeSpectrum(ParticlesSizeSpectrum): - def __init__(self, radius_bins_edges, normalise_by_dv=False): +class ParticleSizeSpectrumPerVolume(ParticleSizeSpectrum): + def __init__(self, radius_bins_edges, dry=False, name=None, unit='m^-3 m^-1'): super().__init__( radius_bins_edges, - dry=True, - normalise_by_dv=normalise_by_dv, - name='Particles Dry Size Spectrum' + dry=dry, + normalise_by_dv=False, + name=name, + unit=unit ) diff --git a/PySDM/products/size_spectral/particles_concentration.py b/PySDM/products/size_spectral/particles_concentration.py new file mode 100644 index 000000000..fd9ad1244 --- /dev/null +++ b/PySDM/products/size_spectral/particles_concentration.py @@ -0,0 +1,28 @@ +import numpy as np +from PySDM.products.impl.moment_product import MomentProduct + + +class ParticleConcentration(MomentProduct): + def __init__(self, radius_range=(0, np.inf), specific=False, name=None, unit='m^-3'): + self.radius_range = radius_range + self.specific = specific + super().__init__(name=name, unit=unit) + + def _impl(self, **kwargs): + self._download_moment_to_buffer( + 'volume', rank=0, + filter_range=(self.formulae.trivia.volume(self.radius_range[0]), + self.formulae.trivia.volume(self.radius_range[1]))) + self.buffer[:] /= self.particulator.mesh.dv + if self.specific: + result = self.buffer.copy() # TODO #217 + self._download_to_buffer(self.particulator.environment['rhod']) + result[:] /= self.buffer + else: + result = self.buffer + return result + + +class ParticleSpecificConcentration(ParticleConcentration): + def __init__(self, radius_range=(0, np.inf), name=None, unit='kg^-1'): + super().__init__(radius_range=radius_range, specific=True, name=name, unit=unit) diff --git a/PySDM/products/particles_volume_spectrum.py b/PySDM/products/size_spectral/particles_volume_spectrum.py similarity index 64% rename from PySDM/products/particles_volume_spectrum.py rename to PySDM/products/size_spectral/particles_volume_spectrum.py index 446971afd..38df04bc3 100644 --- a/PySDM/products/particles_volume_spectrum.py +++ b/PySDM/products/size_spectral/particles_volume_spectrum.py @@ -1,15 +1,10 @@ import numpy as np -from PySDM.impl.product import SpectrumMomentProduct +from PySDM.products.impl.spectrum_moment_product import SpectrumMomentProduct class ParticlesVolumeSpectrum(SpectrumMomentProduct): - - def __init__(self, radius_bins_edges): - super().__init__( - name='dv/dlnr', - unit='1/(unit dr/r)', - description='Particles volume distribution' - ) + def __init__(self, radius_bins_edges, name=None, unit='dimensionless'): + super().__init__(name=name, unit=unit) self.radius_bins_edges = radius_bins_edges self.moment_0 = None self.moments = None @@ -24,14 +19,14 @@ def register(self, builder): self.shape = (*builder.particulator.mesh.grid, len(self.attr_bins_edges) - 1) - def get(self): + def _impl(self, **kwargs): vals = np.empty([self.particulator.mesh.n_cell, len(self.attr_bins_edges) - 1]) - self.recalculate_spectrum_moment(attr='volume', rank=1, filter_attr='volume') + self._recalculate_spectrum_moment(attr='volume', rank=1, filter_attr='volume') for i in range(vals.shape[1]): - self.download_spectrum_moment_to_buffer(rank=1, bin_number=i) + self._download_spectrum_moment_to_buffer(rank=1, bin_number=i) vals[:, i] = self.buffer.ravel() - self.download_spectrum_moment_to_buffer(rank=0, bin_number=i) + self._download_spectrum_moment_to_buffer(rank=0, bin_number=i) vals[:, i] *= self.buffer.ravel() vals *= 1 / np.diff(np.log(self.radius_bins_edges)) / self.particulator.mesh.dv diff --git a/PySDM/products/size_spectral/total_particle_concentration.py b/PySDM/products/size_spectral/total_particle_concentration.py new file mode 100644 index 000000000..a66d0bac6 --- /dev/null +++ b/PySDM/products/size_spectral/total_particle_concentration.py @@ -0,0 +1,11 @@ +from PySDM.products.impl.moment_product import MomentProduct + + +class TotalParticleConcentration(MomentProduct): + def __init__(self, name=None, unit='m^-3'): + super().__init__(name=name, unit=unit) + + def _impl(self, **kwargs): + self._download_moment_to_buffer('volume', rank=0) + self.buffer[:] /= self.particulator.mesh.dv + return self.buffer diff --git a/PySDM/products/size_spectral/total_particle_specific_concentration.py b/PySDM/products/size_spectral/total_particle_specific_concentration.py new file mode 100644 index 000000000..e36ab25db --- /dev/null +++ b/PySDM/products/size_spectral/total_particle_specific_concentration.py @@ -0,0 +1,14 @@ +from PySDM.products.impl.moment_product import MomentProduct + + +class TotalParticleSpecificConcentration(MomentProduct): + def __init__(self, name=None, unit='kg^-1'): + super().__init__(name=name, unit=unit) + + def _impl(self, **kwargs): + self._download_moment_to_buffer('volume', rank=0) + result = self.buffer.copy() # TODO #217 + self._download_to_buffer(self.particulator.environment['rhod']) + result[:] /= self.particulator.mesh.dv + result[:] /= self.buffer + return result diff --git a/PySDM/products/size_spectral/water_mixing_ratio.py b/PySDM/products/size_spectral/water_mixing_ratio.py new file mode 100644 index 000000000..3f238002b --- /dev/null +++ b/PySDM/products/size_spectral/water_mixing_ratio.py @@ -0,0 +1,32 @@ +import numpy as np +from PySDM.products.impl.moment_product import MomentProduct +from PySDM.physics import constants as const + + +class WaterMixingRatio(MomentProduct): + + def __init__(self, radius_range=(0, np.inf), name=None, unit='dimensionless'): + self.radius_range = radius_range + self.volume_range = None + super().__init__(unit=unit, name=name) + + def register(self, builder): + super().register(builder) + self.volume_range = self.formulae.trivia.volume(np.asarray(self.radius_range)) + self.radius_range = None + + def _impl(self, **kwargs): # TODO #217 + self._download_moment_to_buffer('volume', rank=0, + filter_range=self.volume_range, filter_attr='volume') + conc = self.buffer.copy() + + 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[:] *= conc + result[:] /= self.particulator.mesh.dv + + self._download_to_buffer(self.particulator.environment['rhod']) + result[:] /= self.buffer + return result diff --git a/PySDM/products/super_droplet_count.py b/PySDM/products/super_droplet_count.py deleted file mode 100644 index 97504fc6b..000000000 --- a/PySDM/products/super_droplet_count.py +++ /dev/null @@ -1,18 +0,0 @@ -from PySDM.impl.product import Product - - -class SuperDropletCount(Product): - - def __init__(self): - super().__init__( - name='n_sd', - unit='#/gridbox', - description='Super droplet count' - ) - - def get(self): - cell_start = self.particulator.attributes.cell_start - n_cell = cell_start.shape[0] - 1 - for i in range(n_cell): - self.buffer.ravel()[i] = cell_start[i + 1] - cell_start[i] - return self.buffer diff --git a/PySDM/products/temperature.py b/PySDM/products/temperature.py deleted file mode 100644 index 7e324e301..000000000 --- a/PySDM/products/temperature.py +++ /dev/null @@ -1,11 +0,0 @@ -from PySDM.impl.product import MoistEnvironmentProduct - - -class Temperature(MoistEnvironmentProduct): - - def __init__(self): - super().__init__( - description="Temperature", - name="T", - unit="K" - ) diff --git a/PySDM/products/total_dry_mass_mixing_ratio.py b/PySDM/products/total_dry_mass_mixing_ratio.py deleted file mode 100644 index fdccbfc1b..000000000 --- a/PySDM/products/total_dry_mass_mixing_ratio.py +++ /dev/null @@ -1,25 +0,0 @@ -import numpy as np -from PySDM.physics.constants import convert_to, si -from PySDM.impl.product import MomentProduct - - -class TotalDryMassMixingRatio(MomentProduct): - def __init__(self, density): - super().__init__( - name='q_dry', - description='total dry mass mixing ratio', - unit='μg/kg' - ) - self.density = density - - def get(self): - self.download_moment_to_buffer('dry volume', rank=1) - self.buffer[:] *= self.density - result = np.copy(self.buffer) - self.download_moment_to_buffer('dry volume', rank=0) - result[:] *= self.buffer - self.download_to_buffer(self.particulator.environment['rhod']) - result[:] /= self.particulator.mesh.dv - result[:] /= self.buffer - convert_to(result, si.ug / si.kg) - return result diff --git a/PySDM/products/total_particle_concentration.py b/PySDM/products/total_particle_concentration.py deleted file mode 100644 index bf2f25ccc..000000000 --- a/PySDM/products/total_particle_concentration.py +++ /dev/null @@ -1,18 +0,0 @@ -from PySDM.physics import constants as const -from PySDM.impl.product import MomentProduct - - -class TotalParticleConcentration(MomentProduct): - - def __init__(self): - super().__init__( - name='n_cm3', - unit='cm-3', - description='Total particle concentration' - ) - - def get(self): - self.download_moment_to_buffer('volume', rank=0) - self.buffer[:] /= self.particulator.mesh.dv - const.convert_to(self.buffer, const.si.centimetre**-3) - return self.buffer diff --git a/PySDM/products/total_particle_specific_concentration.py b/PySDM/products/total_particle_specific_concentration.py deleted file mode 100644 index 7a56c67c7..000000000 --- a/PySDM/products/total_particle_specific_concentration.py +++ /dev/null @@ -1,21 +0,0 @@ -from PySDM.physics import constants as const -from PySDM.impl.product import MomentProduct - - -class TotalParticleSpecificConcentration(MomentProduct): - - def __init__(self): - super().__init__( - name='n_mg', - unit='mg-1', - description='Total particle specific concentration' - ) - - def get(self): - self.download_moment_to_buffer('volume', rank=0) - result = self.buffer.copy() # TODO #217 - self.download_to_buffer(self.particulator.environment['rhod']) - result[:] /= self.particulator.mesh.dv - result[:] /= self.buffer - const.convert_to(result, const.si.milligram**-1) - return result diff --git a/PySDM/products/water_mixing_ratio.py b/PySDM/products/water_mixing_ratio.py deleted file mode 100644 index 3693da176..000000000 --- a/PySDM/products/water_mixing_ratio.py +++ /dev/null @@ -1,37 +0,0 @@ -import numpy as np -from PySDM.impl.product import MomentProduct -from PySDM.physics import constants as const - - -class WaterMixingRatio(MomentProduct): - - def __init__(self, radius_range, name='ql', description_prefix='liquid'): - self.radius_range = radius_range - self.volume_range = None - super().__init__( - name=name, - unit='g/kg', - description=f'{description_prefix} water mixing ratio' - ) - - def register(self, builder): - super().register(builder) - self.volume_range = self.formulae.trivia.volume(np.asarray(self.radius_range)) - self.radius_range = None - - def get(self): # TODO #217 - self.download_moment_to_buffer('volume', rank=0, - filter_range=self.volume_range, filter_attr='volume') - conc = self.buffer.copy() - - 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[:] *= conc - result[:] /= self.particulator.mesh.dv - - self.download_to_buffer(self.particulator.environment['rhod']) - result[:] /= self.buffer - const.convert_to(result, const.si.gram / const.si.kilogram) - return result diff --git a/PySDM/products/water_vapour_mixing_ratio.py b/PySDM/products/water_vapour_mixing_ratio.py deleted file mode 100644 index 492716313..000000000 --- a/PySDM/products/water_vapour_mixing_ratio.py +++ /dev/null @@ -1,17 +0,0 @@ -from PySDM.physics import constants as const -from PySDM.impl.product import MoistEnvironmentProduct - - -class WaterVapourMixingRatio(MoistEnvironmentProduct): - - def __init__(self): - super().__init__( - description="Water vapour mixing ratio", - name="qv", - unit="g/kg" - ) - - def get(self): - super().get() - const.convert_to(self.buffer, const.si.gram / const.si.kilogram) - return self.buffer diff --git a/README.md b/README.md index 4db6ed2de..cd9c6da43 100644 --- a/README.md +++ b/README.md @@ -176,7 +176,7 @@ radius_bins_edges = 10 .^ range(log10(10*si.um), log10(5e3*si.um), length=32) builder = Builder(n_sd=n_sd, backend=CPU()) builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m^3)) builder.add_dynamic(Coalescence(kernel=Golovin(b=1.5e3 / si.s))) -products = [ParticlesVolumeSpectrum(radius_bins_edges)] +products = [ParticlesVolumeSpectrum(radius_bins_edges=radius_bins_edges, name="dv/dlnr")] particulator = builder.build(attributes, products) ``` @@ -196,7 +196,10 @@ radius_bins_edges = logspace(log10(10 * si.um), log10(5e3 * si.um), 32); builder = Builder(pyargs('n_sd', int32(n_sd), 'backend', CPU())); builder.set_environment(Box(pyargs('dt', 1 * si.s, 'dv', 1e6 * si.m ^ 3))); builder.add_dynamic(Coalescence(pyargs('kernel', Golovin(1.5e3 / si.s)))); -products = py.list({ ParticlesVolumeSpectrum(py.numpy.array(radius_bins_edges)) }); +products = py.list({ ParticlesVolumeSpectrum(pyargs( ... + 'radius_bins_edges', py.numpy.array(radius_bins_edges), ... + 'name', 'dv/dlnr' ... +)) }); particulator = builder.build(attributes, products); ``` @@ -217,7 +220,7 @@ radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(5e3 * si.um), num builder = Builder(n_sd=n_sd, backend=CPU()) builder.set_environment(Box(dt=1 * si.s, dv=1e6 * si.m**3)) builder.add_dynamic(Coalescence(kernel=Golovin(b=1.5e3 / si.s))) -products = [ParticlesVolumeSpectrum(radius_bins_edges)] +products = [ParticlesVolumeSpectrum(radius_bins_edges=radius_bins_edges, name='dv/dlnr')] particulator = builder.build(attributes, products) ``` @@ -321,13 +324,13 @@ initial humidity. Subsequent particle growth due to [`Condensation`](https://atmos-cloud-sim-uj.github.io/PySDM/dynamics/condensation.html) of water vapour (coupled with the release of latent heat) causes a subset of particles to activate into cloud droplets. Results of the simulation are plotted against vertical -[`ParcelDisplacement`](https://atmos-cloud-sim-uj.github.io/PySDM/products/environments/parcel_displacement.html) +[`ParcelDisplacement`](https://atmos-cloud-sim-uj.github.io/PySDM/products/housekeeping/parcel_displacement.html) and depict the evolution of -[`Supersaturation`](https://atmos-cloud-sim-uj.github.io/PySDM/products/dynamics/condensation/peak_supersaturation.html), -[`CloudDropletEffectiveRadius`](https://atmos-cloud-sim-uj.github.io/PySDM/products/state/cloud_droplet_effective_radius.html), -[`CloudDropletConcentration`](https://atmos-cloud-sim-uj.github.io/PySDM/products/state/particles_concentration.html#PySDM.products.particles_concentration.CloudDropletConcentration) +[`PeakSupersaturation`](https://atmos-cloud-sim-uj.github.io/PySDM/products/condensation/peak_supersaturation.html), +[`EffectiveRadius`](https://atmos-cloud-sim-uj.github.io/PySDM/products/size_spectral/effective_radius.html), +[`ParticleConcentration`](https://atmos-cloud-sim-uj.github.io/PySDM/products/size_spectral/particle_concentration.html#PySDM.products.particles_concentration.ParticleConcentration) and the -[`WaterMixingRatio `](https://atmos-cloud-sim-uj.github.io/PySDM/products/state/water_mixing_ratio.html). +[`WaterMixingRatio `](https://atmos-cloud-sim-uj.github.io/PySDM/products/size_spectral/water_mixing_ratio.html).
Julia (click to expand) @@ -379,11 +382,11 @@ attributes["kappa times dry volume"] = kappa * v_dry attributes["volume"] = builder.formulae.trivia.volume(radius=r_wet) particulator = builder.build(attributes, products=[ - products.PeakSupersaturation(), - products.CloudDropletEffectiveRadius(radius_range=cloud_range), - products.CloudDropletConcentration(radius_range=cloud_range), - products.WaterMixingRatio(radius_range=cloud_range), - products.ParcelDisplacement() + products.PeakSupersaturation(name="S_max", unit="%"), + products.EffectiveRadius(name="r_eff", unit="um", radius_range=cloud_range), + products.ParticleConcentration(name="n_c_cm3", unit="cm^-3", radius_range=cloud_range), + products.WaterMixingRatio(name="ql", unit="g/kg", radius_range=cloud_range), + products.ParcelDisplacement(name="z") ]) cell_id=1 @@ -462,11 +465,11 @@ attributes = py.dict(pyargs( ... )); particulator = builder.build(attributes, py.list({ ... - products.PeakSupersaturation(), ... - products.CloudDropletEffectiveRadius(pyargs('radius_range', cloud_range)), ... - products.CloudDropletConcentration(pyargs('radius_range', cloud_range)), ... - products.WaterMixingRatio(pyargs('radius_range', cloud_range)) ... - products.ParcelDisplacement() ... + products.PeakSupersaturation(pyargs('name', 'S_max', 'unit', '%')), ... + products.EffectiveRadius(pyargs('name', 'r_eff', 'unit', 'um', 'radius_range', cloud_range)), ... + products.ParticleConcentration(pyargs('name', 'n_c_cm3', 'unit', 'cm^-3', 'radius_range', cloud_range)), ... + products.WaterMixingRatio(pyargs('name', 'ql', 'unit', 'g/kg', 'radius_range', cloud_range)) ... + products.ParcelDisplacement(pyargs('name', 'z')) ... })); cell_id = int32(0); @@ -514,7 +517,6 @@ saveas(gcf, "parcel.png") Python (click to expand) ```Python -import numpy as np from matplotlib import pyplot from PySDM.physics import si, spectra from PySDM.initialisation import spectral_sampling, multiplicities, r_wet_init @@ -524,12 +526,12 @@ from PySDM.environments import Parcel from PySDM import Builder, products env = Parcel( - dt=.25 * si.s, - mass_of_dry_air=1e3 * si.kg, - p0=1122 * si.hPa, - q0=20 * si.g / si.kg, - T0=300 * si.K, - w=2.5 * si.m / si.s + dt=.25 * si.s, + mass_of_dry_air=1e3 * si.kg, + p0=1122 * si.hPa, + q0=20 * si.g / si.kg, + T0=300 * si.K, + w=2.5 * si.m / si.s ) spectrum = spectra.Lognormal(norm_factor=1e4 / si.mg, m_mode=50 * si.nm, s_geom=1.5) kappa = .5 * si.dimensionless @@ -548,35 +550,35 @@ v_dry = builder.formulae.trivia.volume(radius=r_dry) r_wet = r_wet_init(r_dry, env, kappa * v_dry) attributes = { - 'n': multiplicities.discretise_n(specific_concentration * env.mass_of_dry_air), - 'dry volume': v_dry, - 'kappa times dry volume': kappa * v_dry, - 'volume': builder.formulae.trivia.volume(radius=r_wet) + 'n': multiplicities.discretise_n(specific_concentration * env.mass_of_dry_air), + 'dry volume': v_dry, + 'kappa times dry volume': kappa * v_dry, + 'volume': builder.formulae.trivia.volume(radius=r_wet) } particulator = builder.build(attributes, products=[ - products.PeakSupersaturation(), - products.CloudDropletEffectiveRadius(radius_range=cloud_range), - products.CloudDropletConcentration(radius_range=cloud_range), - products.WaterMixingRatio(radius_range=cloud_range), - products.ParcelDisplacement() + products.PeakSupersaturation(name='S_max', unit='%'), + products.EffectiveRadius(name='r_eff', unit='um', radius_range=cloud_range), + products.ParticleConcentration(name='n_c_cm3', unit='cm^-3', radius_range=cloud_range), + products.WaterMixingRatio(name='ql', unit='g/kg', radius_range=cloud_range), + products.ParcelDisplacement(name='z') ]) cell_id = 0 output = {product.name: [product.get()[cell_id]] for product in particulator.products.values()} for step in range(output_points): - particulator.run(steps=output_interval) - for product in particulator.products.values(): - output[product.name].append(product.get()[cell_id]) + particulator.run(steps=output_interval) + for product in particulator.products.values(): + output[product.name].append(product.get()[cell_id]) fig, axs = pyplot.subplots(1, len(particulator.products) - 1, sharey="all") for i, (key, product) in enumerate(particulator.products.items()): - if key != 'z': - axs[i].plot(output[key], output['z'], marker='.') - axs[i].set_title(product.name) - axs[i].set_xlabel(product.unit) - axs[i].grid() + if key != 'z': + axs[i].plot(output[key], output['z'], marker='.') + axs[i].set_title(product.name) + axs[i].set_xlabel(product.unit) + axs[i].grid() axs[0].set_ylabel(particulator.products['z'].unit) pyplot.savefig('parcel.svg') ``` diff --git a/test-time-requirements.txt b/test-time-requirements.txt index d49569ae8..755e71cf8 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@81a9bed#egg=PySDM-examples +PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@cb30220#egg=PySDM-examples 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 7126a99a8..7539debab 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 @@ -19,7 +19,7 @@ def test_at_t_0(): output = simulation.run() # Assert - np.testing.assert_allclose(output['RH_env'][0], 95) + np.testing.assert_allclose(output['RH'][0], 95) np.testing.assert_allclose(output['gas_S_IV_ppb'][0], 0.2) np.testing.assert_allclose(output['gas_N_mIII_ppb'][0], 0.1) np.testing.assert_allclose(output['gas_H2O2_ppb'], 0.5) @@ -36,7 +36,7 @@ def test_at_t_0(): np.testing.assert_allclose(num_conc_NH4p, num_conc_SO4mm, rtol=.005) mass_conc_H = num_conc_NH4p * Substance.from_formula("H").mass * si.gram / si.mole np.testing.assert_allclose( - actual=np.asarray(output['q_dry'])*np.asarray(output['rhod_env']), + actual=np.asarray(output['q_dry'])*np.asarray(output['rhod']), desired=mass_conc_NH4p + mass_conc_SO4mm + mass_conc_H, rtol=rtol ) @@ -54,7 +54,7 @@ def test_at_t_0(): settings.formulae.trivia.mole_fraction_2_mixing_ratio( mole_fraction, specific_gravity=SPECIFIC_GRAVITY[compound] - ) * np.asarray(output['rhod_env']) + ) * np.asarray(output['rhod']) ), desired=expected[key], rtol=rtol @@ -73,18 +73,18 @@ def test_at_cloud_base(): # Assert assert round(output['z'][-1]) == (698 - 600) * si.m - np.testing.assert_allclose(output['p_env'][-1], 939 * si.mbar, rtol=.005) - np.testing.assert_allclose(output['T_env'][-1], 284.2 * si.K, rtol=.005) + np.testing.assert_allclose(output['p'][-1], 939 * si.mbar, rtol=.005) + np.testing.assert_allclose(output['T'][-1], 284.2 * si.K, rtol=.005) np.testing.assert_allclose( settings.formulae.state_variable_triplet.rho_of_rhod_qv( - rhod=output['rhod_env'][-1], qv=output['qv_env'][-1]*si.g/si.kg), + rhod=output['rhod'][-1], qv=output['qv'][-1]*si.g/si.kg), 1.15 * si.kg / si.m**3, rtol=.005 ) assert output['ql'][-2] < .00055 assert output['ql'][-1] > .0004 - assert output['RH_env'][-1] > 100 - assert output['RH_env'][-8] < 100 + assert output['RH'][-1] > 100 + assert output['RH'][-8] < 100 @staticmethod def test_at_1200m_above_cloud_base(): diff --git a/tests/smoke_tests/shipway_and_hill_2012/test_few_steps.py b/tests/smoke_tests/shipway_and_hill_2012/test_few_steps.py index ef791c6d5..59e2db911 100644 --- a/tests/smoke_tests/shipway_and_hill_2012/test_few_steps.py +++ b/tests/smoke_tests/shipway_and_hill_2012/test_few_steps.py @@ -19,8 +19,8 @@ def profile(var): return np.mean(output[var][:, -20:], axis=1) if plot: - for var in ('RH_env', 'S_max', 'T_env', 'qv_env', 'p_env', 'ql', - 'ripening_rate', 'activating_rate', 'deactivating_rate'): + for var in ('RH', 'S_max', 'T', 'qv', 'p', 'ql', + 'ripening rate', 'activating rate', 'deactivating rate'): pyplot.plot(profile(var), output['z'], linestyle='--', marker='o') pyplot.ylabel('Z [m]') pyplot.xlabel(var + ' [' + simulation.particulator.products[var].unit + ']') @@ -31,5 +31,5 @@ def profile(var): assert min(profile('ql')) == 0 assert .1 < max(profile('ql')) < 1 # assert max(profile('ripening_rate')) > 0 # TODO #521 - assert max(profile('activating_rate')) == 0 + assert max(profile('activating rate')) == 0 # assert max(profile('deactivating_rate')) > 0 TODO #521 diff --git a/tests/smoke_tests/shipway_and_hill_2012/test_initial_condition.py b/tests/smoke_tests/shipway_and_hill_2012/test_initial_condition.py index 2e9099aa1..318edf695 100644 --- a/tests/smoke_tests/shipway_and_hill_2012/test_initial_condition.py +++ b/tests/smoke_tests/shipway_and_hill_2012/test_initial_condition.py @@ -17,7 +17,7 @@ def test_initial_condition(plot=False): # Plot if plot: - for var in ('RH_env', 'T_env', 'qv_env', 'p_env'): + for var in ('RH', 'T', 'qv', 'p'): pyplot.plot(output[var][:, 0], output['z'], linestyle='--', marker='o') pyplot.ylabel('Z [m]') pyplot.xlabel(var + ' [' + simulation.particulator.products[var].unit + ']') @@ -25,15 +25,15 @@ def test_initial_condition(plot=False): pyplot.show() # Assert - assert output['RH_env'].shape == (settings.nz, 1) + assert output['RH'].shape == (settings.nz, 1) - assert 35 < np.amin(output['RH_env']) < 40 - assert 110 < np.amax(output['RH_env']) < 115 + assert 35 < np.amin(output['RH']) < 40 + assert 110 < np.amax(output['RH']) < 115 - assert 700 * si.hPa < np.amin(output['p_env']) < 710 * si.hPa - assert (np.diff(output['p_env']) < 0).all() - assert 950 * si.hPa < np.amax(output['p_env']) < 1000 * si.hPa + assert 700 * si.hPa < np.amin(output['p']) < 710 * si.hPa + assert (np.diff(output['p']) < 0).all() + assert 950 * si.hPa < np.amax(output['p']) < 1000 * si.hPa - assert 280 * si.K < np.amin(output['T_env']) < 285 * si.K - assert output['T_env'][0] > np.amin(output['T_env']) - assert 295 * si.K < np.amax(output['T_env']) < 300 * si.K + assert 280 * si.K < np.amin(output['T']) < 285 * si.K + assert output['T'][0] > np.amin(output['T']) + assert 295 * si.K < np.amax(output['T']) < 300 * si.K diff --git a/tests/smoke_tests/yang_et_al_2018/test_just_do_it.py b/tests/smoke_tests/yang_et_al_2018/test_just_do_it.py index fb019230b..efd658439 100644 --- a/tests/smoke_tests/yang_et_al_2018/test_just_do_it.py +++ b/tests/smoke_tests/yang_et_al_2018/test_just_do_it.py @@ -52,7 +52,7 @@ def test_just_do_it(backend_class, scheme, adaptive): # TODO #527 if backend_class is not GPU: - assert max(output['ripening_rate']) > 0 + assert max(output['ripening rate']) > 0 def n_tot(n, condition): diff --git a/tests/unit_tests/attributes/test_critical_supersaturation.py b/tests/unit_tests/attributes/test_critical_supersaturation.py index 246a8a3fd..6d7d15c5e 100644 --- a/tests/unit_tests/attributes/test_critical_supersaturation.py +++ b/tests/unit_tests/attributes/test_critical_supersaturation.py @@ -26,11 +26,11 @@ def test_critical_supersaturation(): 'kappa times dry volume': .9 * vdry, 'dry volume organic': np.zeros(n_sd) }, - products=[ActivableFraction()] + products=(ActivableFraction(),) ) # act - AF = particulator.products['activable fraction'].get(S_max) + AF = particulator.products['activable fraction'].get(S_max=S_max) # assert assert 0 < AF < 1 diff --git a/tests/unit_tests/backends/test_freezing_methods.py b/tests/unit_tests/backends/test_freezing_methods.py index a168714dc..2528c7e3a 100644 --- a/tests/unit_tests/backends/test_freezing_methods.py +++ b/tests/unit_tests/backends/test_freezing_methods.py @@ -66,7 +66,7 @@ def test_freeze_time_dependent(plot=False): 'immersed surface area': np.full(n_sd, immersed_surface_area), 'volume': np.full(n_sd, vol) } - products = (IceWaterContent(specific=False),) + products = (IceWaterContent(name='qi'),) particulator = builder.build(attributes=attributes, products=products) env['a_w_ice'] = np.nan diff --git a/tests/unit_tests/test_products.py b/tests/unit_tests/test_products.py new file mode 100644 index 000000000..1c3503132 --- /dev/null +++ b/tests/unit_tests/test_products.py @@ -0,0 +1,53 @@ +import inspect +from typing import Tuple +import sys +import numpy as np +import pytest +from PySDM import products +from PySDM.products.impl.product import Product +from PySDM.products import (AqueousMassSpectrum, AqueousMoleFraction, TotalDryMassMixingRatio, + ParticleSizeSpectrumPerMass, GaseousMoleFraction, + FreezableSpecificConcentration, DynamicWallTime, + ParticleSizeSpectrumPerVolume, ParticlesVolumeSpectrum) + +_ARGUMENTS = { + AqueousMassSpectrum: {'key': 'S_VI', 'dry_radius_bins_edges': (0, np.inf)}, + AqueousMoleFraction: {'key': 'S_VI'}, + TotalDryMassMixingRatio: {'density': 1}, + ParticleSizeSpectrumPerMass: {'radius_bins_edges': (0, np.inf)}, + GaseousMoleFraction: {'key': 'O3'}, + FreezableSpecificConcentration: {'temperature_bins_edges': (0, 300)}, + DynamicWallTime: {'dynamic': 'Condensation'}, + ParticleSizeSpectrumPerVolume: {'radius_bins_edges': (0, np.inf)}, + ParticlesVolumeSpectrum: {'radius_bins_edges': (0, np.inf)} +} + + +@pytest.fixture(params=( + pytest.param(p[1], id=p[0]) + for p in inspect.getmembers(sys.modules[products.__name__], inspect.isclass) +)) +def product(request): + return request.param + + +class TestProducts: + @staticmethod + # pylint: disable=redefined-outer-name + def test_instantiate_all(product): + product(**(_ARGUMENTS[product] if product in _ARGUMENTS else {})) + + @staticmethod + @pytest.mark.parametrize('in_out_pair', ( + ('CPUTime', 'CPU time'), + ('WallTime', 'wall time') + )) + def test_camel_case_to_words(in_out_pair: Tuple[str, str]): + # arrange + test_input, expected_output = in_out_pair + + # act + actual_output = Product._camel_case_to_words(test_input) + + # assert + assert actual_output == expected_output