Surf/testing

.github/workflows/main.yml

@@ -1,4 +1,4 @@
-name: Build Status
+name: PySDM-examples
 
 defaults:
   run:
@@ -27,11 +27,11 @@ jobs:
           python-version: ${{ matrix.python-version }}
 
       - run: |
-          pip install pytest nbconvert
+          pip install -r test-time-requirements.txt
       - run: |
           pip install -e .
 
-      # TODO #184 (https://github.com/numba/numba/issues/6350#issuecomment-728174860)
+      # https://github.com/numba/numba/issues/6350#issuecomment-728174860
       - if: matrix.platform == 'ubuntu-latest'
         env: 
           NUMBA_THREADING_LAYER: omp

PySDM_examples/Arabas_and_Shima_2017/settings.py

@@ -5,20 +5,24 @@
 from PySDM.physics.constants import si
 from PySDM.physics import constants as const, formulae as phys
 from PySDM.dynamics import condensation
+from PySDM.physics.formulae import Formulae
 import numpy as np
 from pystrict import strict
 
 
 @strict
 class Settings:
-    def __init__(self, w_avg: float, N_STP: float, r_dry: float, mass_of_dry_air: float):
+    def __init__(self, w_avg: float, N_STP: float, r_dry: float, mass_of_dry_air: float, coord: str = 'VolumeLogarithm'):
+        self.formulae = Formulae(saturation_vapour_pressure='AugustRocheMagnus', condensation_coordinate=coord)
+
         self.p0 = 1000 * si.hectopascals
         self.RH0 = .98
         self.kappa = .2  # TODO #441
         self.T0 = 300 * si.kelvin
         self.z_half = 150 * si.metres
 
-        self.q0 = const.eps / (self.p0 / self.RH0 / phys.pvs(self.T0) - 1)
+        pvs = self.formulae.saturation_vapour_pressure.pvs_Celsius(self.T0 - const.T0)
+        self.q0 = const.eps / (self.p0 / self.RH0 / pvs - 1)
         self.w_avg = w_avg
         self.r_dry = r_dry
         self.N_STP = N_STP
@@ -29,7 +33,7 @@ def __init__(self, w_avg: float, N_STP: float, r_dry: float, mass_of_dry_air: fl
         self.rtol_x = condensation.default_rtol_x
         self.rtol_thd = condensation.default_rtol_thd
         self.coord = 'volume logarithm'
-        self.dt_cond_range = (1e-10 * si.s, 1 * si.s)
+        self.dt_cond_range = condensation.default_cond_range
 
     @property
     def dt_max(self):

PySDM_examples/Arabas_and_Shima_2017/simulation.py

@@ -24,7 +24,7 @@ def __init__(self, settings, backend=CPU):
         while dt_output / self.n_substeps >= settings.dt_max:  # TODO #334 dt_max
             self.n_substeps += 1
 
-        builder = Builder(backend=backend, n_sd=1)
+        builder = Builder(backend=backend, n_sd=1, formulae=settings.formulae)
         builder.set_environment(Parcel(
             dt=dt_output / self.n_substeps,
             mass_of_dry_air=settings.mass_of_dry_air,
@@ -43,11 +43,11 @@ def __init__(self, settings, backend=CPU):
         ))
         attributes = {}
         r_dry = np.array([settings.r_dry])
-        attributes['dry volume'] = phys.volume(radius=r_dry)
+        attributes['dry volume'] = settings.formulae.trivia.volume(radius=r_dry)
         attributes['n'] = np.array([settings.n_in_dv], dtype=np.int64)
         environment = builder.core.environment
-        r_wet = r_wet_init(r_dry, environment, np.zeros_like(attributes['n']), settings.kappa)
-        attributes['volume'] = phys.volume(radius=r_wet)
+        r_wet = r_wet_init(r_dry, environment, kappa=settings.kappa)
+        attributes['volume'] = settings.formulae.trivia.volume(radius=r_wet)
         products = [
             PySDM_products.ParticleMeanRadius(),
             PySDM_products.CondensationTimestepMin(),
@@ -56,7 +56,8 @@ def __init__(self, settings, backend=CPU):
             PySDM_products.Time(),
             PySDM_products.ActivatingRate(),
             PySDM_products.DeactivatingRate(),
-            PySDM_products.RipeningRate()
+            PySDM_products.RipeningRate(),
+            PySDM_products.PeakSupersaturation()
         ]
 
         self.core = builder.build(attributes, products)
@@ -67,7 +68,7 @@ def save(self, output):
         cell_id = 0
         output["r"].append(self.core.products['radius_m1'].get(unit=const.si.metre)[cell_id])
         output["dt_cond_min"].append(self.core.products['dt_cond_min'].get()[cell_id])
-        output["z"].append(self.core.products["z"].get())
+        output["z"].append(self.core.products["z"].get()[cell_id])
         output["S"].append(self.core.products["RH_env"].get()[cell_id]/100 - 1)
         output["t"].append(self.core.products["t"].get())
 

PySDM_examples/Arabas_et_al_2015/demo_plots.py

@@ -5,7 +5,7 @@
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
-
+from PySDM import products
 
 class _Plot:
 
@@ -37,9 +37,36 @@ def __init__(self, fig, ax, grid, size, product, show=False, lines=False, cmap='
             self.lines['X'][1] = plt.plot([-1] * 2, zlim, **self.line_args)[0]
             self.lines['Z'][1] = plt.plot(xlim, [-1] * 2, **self.line_args)[0]
 
-        data = np.full_like(self.nans, product.range[0])
+        product_range, product_scale = {
+            products.WaterMixingRatio: ((0, 1), 'linear'),
+            products.TotalParticleSpecificConcentration: ((20, 50), 'linear'),
+            products.TotalParticleConcentration: ((20, 50), 'linear'),
+            products.SuperDropletCount: ((0, 10), 'linear'),
+            products.ParticlesVolumeSpectrum: ((20, 50), 'linear'),
+            products.AerosolConcentration: ((1e0, 1e2), 'linear'),
+            products.CloudDropletConcentration: ((1e0, 1e2), 'linear'),
+            products.DrizzleConcentration: ((1e-3, 1e1), 'log'),
+            products.ParticleMeanRadius: ((1, 25), 'linear'),
+            products.CloudDropletEffectiveRadius: ((0, 25), 'linear'),
+            products.AerosolSpecificConcentration: ((0, 3e2), 'linear'),
+            products.WaterVapourMixingRatio: ((5, 7.5), 'linear'),
+            products.Temperature: ((275, 305), 'linear'),
+            products.RelativeHumidity: ((75, 105), 'linear'),
+            products.Pressure: ((90000, 100000), 'linear'),
+            products.DryAirPotentialTemperature: ((275, 300), 'linear'),
+            products.DryAirDensity: ((0.95, 1.3), 'linear'),
+            products.PeakSupersaturation: ((-1, 1), 'linear'),
+            products.RipeningRate: ((1e-1, 1e1), 'log'),
+            products.ActivatingRate: ((1e-1, 1e1), 'log'),
+            products.DeactivatingRate: ((1e-1, 1e1), 'log'),
+            products.CollisionRateDeficit: ((0, 1e10), 'linear'),
+            products.CollisionRate: ((0, 1e10), 'linear'),
+            products.CondensationTimestepMin: ((.01, 10), 'log'),
+            products.CondensationTimestepMax: ((.01, 10), 'log'),
+            products.CoalescenceTimestepMin: ((.01, 10), 'log'),
+        }[product.__class__]
+        data = np.full_like(self.nans, product_range[0])
         label = f"{product.description} [{product.unit}]"
-        scale = product.scale
 
         self.ax.set_xlabel('X [m]')
         self.ax.set_ylabel('Z [m]')
@@ -48,11 +75,11 @@ def __init__(self, fig, ax, grid, size, product, show=False, lines=False, cmap='
                                  origin='lower',
                                  extent=(*xlim, *zlim),
                                  cmap=cmap,
-                                 norm=matplotlib.colors.LogNorm() if scale == 'log' and np.isfinite(
+                                 norm=matplotlib.colors.LogNorm() if product_scale == 'log' and np.isfinite(
                                      data).all() else None  # TODO #37 this is always None!!!
                                  )
         plt.colorbar(self.im, ax=self.ax).set_label(label)
-        self.im.set_clim(vmin=product.range[0], vmax=product.range[1])
+        self.im.set_clim(vmin=product_range[0], vmax=product_range[1])
         if show:
             plt.show()
 
@@ -117,5 +144,5 @@ def update(self, data):
         if data is not None:
             self.ydata[0:len(data)] = data[:]
         else:
-            self.ydata[0:len(data)] = np.nan
+            self.ydata[:] = np.nan
         self.timeseries.set_ydata(self.ydata)

PySDM_examples/Arabas_et_al_2015/demo_settings.py

@@ -4,16 +4,19 @@
 
 from ..utils.widgets import IntSlider, FloatSlider, VBox, Checkbox, Accordion, Dropdown
 from PySDM_examples.Arabas_et_al_2015.settings import Settings
+from PySDM import physics
 import numpy as np
 import numba
 import os
+import inspect
 
 
 class DemoSettings:
     __dir__ = Settings.__dir__
 
     def __init__(self):
         settings = Settings()
+
         self.ui_th_std0 = FloatSlider(description="$\\theta_0$ [K]", value=settings.th_std0, min=280, max=300)
         self.ui_qv0 = FloatSlider(description="q$_{v0}$ [g/kg]", value=settings.qv0 * 1000, min=5, max=10)
         self.ui_p0 = FloatSlider(description="p$_0$ [hPa]", value=settings.p0 / 100, min=900, max=1100)
@@ -43,10 +46,6 @@ def __init__(self):
         self.ui_coalescence_adaptive = Checkbox(
             value=settings.condensation_adaptive,
             description='coalescence adaptive time-step')
-        self.ui_condensation_coord = Dropdown(
-            options=['volume', 'volume logarithm'],
-            value=settings.condensation_coord,
-            description='condensational variable coordinate')
         self.ui_processes = [Checkbox(value=settings.processes[key], description=key) for key in settings.processes.keys()]
         self.ui_sdpg = IntSlider(value=settings.n_sd_per_gridbox, description="n_sd/gridbox", min=1, max=1000)
         self.fct_description = "MPDATA: flux-corrected transport option"
@@ -59,27 +58,65 @@ def __init__(self):
             Checkbox(value=settings.mpdata_iga, description=self.iga_description),
             IntSlider(value=settings.mpdata_iters, description=self.nit_description, min=1, max=5)
         ]
+        self.ui_formulae_options = [
+            Dropdown(
+                description=k,
+                options=physics.formulae._choices(getattr(physics, k)).keys(),
+                value=v.default
+            )
+            for k, v in inspect.signature(physics.Formulae.__init__).parameters.items()
+            if k not in ('self', 'fastmath', 'seed')
+        ]
+        self.ui_formulae_options.append(
+            Checkbox(
+                value=inspect.signature(physics.Formulae.__init__).parameters['fastmath'].default,
+                description='fastmath'
+            )
+        )
+        self.ui_output_options = {
+            'interval': IntSlider(description='interval [s]', value=settings.output_interval, min=30, max=60*15),
+            'aerosol_radius_threshold': FloatSlider(
+                description='aerosol/cloud threshold [um]',
+                value=settings.aerosol_radius_threshold / physics.si.um,
+                min=.1, max=1, step=.1),
+            'drizzle_radius_threshold': IntSlider(
+                description='cloud/drizzle threshold [um]',
+                value=settings.drizzle_radius_threshold / physics.si.um,
+                min=10, max=100, step=5)
+        }
 
         # TODO #37
-        self.v_bins = settings.v_bins
-        self.n_sd = settings.n_sd
+        self.r_bins_edges = settings.r_bins_edges
         self.size = settings.size
-        self.rhod = settings.rhod
-        self.field_values = settings.field_values
-        self.aerosol_radius_threshold = settings.aerosol_radius_threshold
-        self.drizzle_radius_threshold = settings.drizzle_radius_threshold
-        self.stream_function = settings.stream_function
         self.condensation_substeps = settings.condensation_substeps
         self.condensation_dt_cond_range = settings.condensation_dt_cond_range
         self.condensation_schedule = settings.condensation_schedule
         self.kernel = settings.kernel
         self.spectrum_per_mass_of_dry_air = settings.spectrum_per_mass_of_dry_air
-        self.n_spin_up = settings.n_spin_up
         self.coalescence_dt_coal_range = settings.coalescence_dt_coal_range
         self.coalescence_optimized_random = settings.coalescence_optimized_random
         self.coalescence_substeps = settings.coalescence_substeps
-        self.output_interval = settings.output_interval
-        self.versions = settings.versions
+
+        for attr in ('n_sd', 'rhod', 'field_values', 'versions', 'n_spin_up', 'stream_function'):
+            setattr(self, attr, getattr(settings, attr))
+
+    @property
+    def aerosol_radius_threshold(self):
+        return self.ui_output_options['aerosol_radius_threshold'].value * physics.si.um
+
+    @property
+    def drizzle_radius_threshold(self):
+        return self.ui_output_options['drizzle_radius_threshold'].value * physics.si.um
+
+    @property
+    def output_interval(self):
+        return self.ui_output_options['interval'].value
+
+    @property
+    def formulae(self) -> physics.Formulae:
+        return physics.Formulae(
+            **{widget.description: widget.value for widget in self.ui_formulae_options}
+        )
 
     @property
     def steps_per_output_interval(self) -> int:
@@ -137,10 +174,6 @@ def condensation_adaptive(self):
     def coalescence_adaptive(self):
         return self.ui_coalescence_adaptive.value
 
-    @property
-    def condensation_coord(self):
-        return self.ui_condensation_coord.value
-
     @property
     def processes(self):
         result = {}
@@ -186,10 +219,14 @@ def box(self):
             VBox([*self.ui_processes]),
             VBox([self.ui_nx, self.ui_nz, self.ui_sdpg, self.ui_dt, self.ui_simulation_time,
                   self.ui_condensation_rtol_x, self.ui_condensation_rtol_thd,
-                  self.ui_condensation_adaptive, self.ui_coalescence_adaptive, self.ui_condensation_coord,
+                  self.ui_condensation_adaptive, self.ui_coalescence_adaptive,
                   *self.ui_mpdata_options]),
+            VBox([*self.ui_formulae_options]),
+            VBox([*self.ui_output_options.values()])
         ])
-        layout.set_title(0, 'environment parameters')
+        layout.set_title(0, 'environment')
         layout.set_title(1, 'processes')
         layout.set_title(2, 'discretisation')
+        layout.set_title(3, 'physics')
+        layout.set_title(4, 'output')
         return layout

PySDM_examples/Arabas_et_al_2015/demo_viewer.py

@@ -43,7 +43,7 @@ def reinit(self, products):
             if len(val.shape) == 2
         ]
 
-        r_bins = phys.radius(volume=self.settings.v_bins)
+        r_bins = self.settings.r_bins_edges.copy()
         const.convert_to(r_bins, const.si.micrometres)
         self.spectrumOutput = Output()
         with self.spectrumOutput:

PySDM_examples/Arabas_et_al_2015/mpdata.py

@@ -63,8 +63,5 @@ def wait(self):
                 self.thread.join()
 
     def step(self):
-        try:  # TODO #417: move this to within PyMPDATA
-            for mpdata in self.mpdatas.values():
-                mpdata.advance(1)
-        except NumbaExperimentalFeatureWarning:
-            pass
+        for mpdata in self.mpdatas.values():
+            mpdata.advance(1)

PySDM_examples/Arabas_et_al_2015/netcdf_exporter.py

@@ -20,15 +20,17 @@ def _create_dimensions(self, ncdf):
         ncdf.createDimension("T", len(self.settings.output_steps))
         for index, label in enumerate(self.XZ):
             ncdf.createDimension(label, self.settings.grid[index])
-        ncdf.createDimension("ParticleVolume", len(self.settings.v_bins) - 1)
+        ncdf.createDimension("ParticleVolume",
+                             len(self.settings.formulae.trivia.volume(self.settings.r_bins_edges)) - 1)
 
     def _create_variables(self, ncdf):
         self.vars["T"] = ncdf.createVariable("T", "f", ["T"])
         self.vars["T"].units = "seconds"
 
         for index, label in enumerate(self.XZ):
             self.vars[label] = ncdf.createVariable(label, "f", (label,))
-            self.vars[label][:] = (self.settings.size[index] / self.settings.grid[index]) * (1 / 2 + np.arange(self.settings.grid[index]))
+            self.vars[label][:] = ((self.settings.size[index] / self.settings.grid[index])
+                                   * (1 / 2 + np.arange(self.settings.grid[index])))
             self.vars[label].units = "metres"
 
         # TODO #340 ParticleVolume var