PySDM.dynamics.condensation

  1"""
  2bespoke condensational growth solver
  3with implicit-in-particle-size integration and adaptive timestepping
  4as in [Bartman 2020 (MSc thesis, Section 3.3)](https://www.ap.uj.edu.pl/diplomas/attachments/file/download/125485)
  5"""  # pylint: disable=line-too-long
  6
  7from collections import namedtuple
  8
  9import numpy as np
 10
 11from PySDM.physics import si
 12from PySDM.dynamics.impl import register_dynamic
 13
 14DEFAULTS = namedtuple("_", ("rtol_x", "rtol_thd", "cond_range", "schedule"))(
 15    rtol_x=1e-6,
 16    rtol_thd=1e-6,
 17    cond_range=(1e-4 * si.second, 1 * si.second),
 18    schedule="dynamic",
 19)
 20
 21
 22@register_dynamic()
 23class Condensation:  # pylint: disable=too-many-instance-attributes
 24    def __init__(
 25        self,
 26        *,
 27        rtol_x=DEFAULTS.rtol_x,
 28        rtol_thd=DEFAULTS.rtol_thd,
 29        substeps: int = 1,
 30        adaptive: bool = True,
 31        dt_cond_range: tuple = DEFAULTS.cond_range,
 32        schedule: str = DEFAULTS.schedule,
 33        max_iters: int = 16,
 34        update_thd: bool = True,
 35    ):
 36        if adaptive and substeps != 1:
 37            raise ValueError(
 38                "if specifying substeps count manually, adaptivity must be disabled"
 39            )
 40
 41        self.particulator = None
 42        self.enable = True
 43
 44        self.rtol_x = rtol_x
 45        self.rtol_thd = rtol_thd
 46
 47        self.rh_max = None
 48        self.success = None
 49
 50        self.__substeps = substeps
 51        self.adaptive = adaptive
 52        self.counters = {}
 53        self.dt_cond_range = dt_cond_range
 54        self.schedule = schedule
 55        self.max_iters = max_iters
 56
 57        self.cell_order = None
 58
 59        self.update_thd = update_thd
 60
 61    def register(self, builder):
 62        self.particulator = builder.particulator
 63
 64        builder._set_condensation_parameters(
 65            dt_range=self.dt_cond_range,
 66            adaptive=self.adaptive,
 67            fuse=32,
 68            multiplier=2,
 69            RH_rtol=1e-7,
 70            max_iters=self.max_iters,
 71        )
 72        builder.request_attribute("critical volume")
 73        builder.request_attribute("kappa")
 74        builder.request_attribute("dry volume organic fraction")
 75        builder.request_attribute("Reynolds number")
 76
 77        for counter in ("n_substeps", "n_activating", "n_deactivating", "n_ripening"):
 78            self.counters[counter] = self.particulator.Storage.empty(
 79                self.particulator.mesh.n_cell, dtype=int
 80            )
 81            if counter == "n_substeps":
 82                self.counters[counter][:] = self.__substeps if not self.adaptive else -1
 83            else:
 84                self.counters[counter][:] = -1
 85
 86        self.rh_max = self.particulator.Storage.empty(
 87            self.particulator.mesh.n_cell, dtype=float
 88        )
 89        self.rh_max[:] = np.nan
 90        self.success = self.particulator.Storage.empty(
 91            self.particulator.mesh.n_cell, dtype=bool
 92        )
 93        self.success[:] = False
 94        self.cell_order = np.arange(self.particulator.mesh.n_cell)
 95
 96    def __call__(self):
 97        if self.enable:
 98            if self.schedule == "dynamic":
 99                self.cell_order = np.argsort(self.counters["n_substeps"])
100            elif self.schedule == "static":
101                pass
102            else:
103                raise NotImplementedError()
104
105            self.particulator.condensation(
106                rtol_x=self.rtol_x,
107                rtol_thd=self.rtol_thd,
108                counters=self.counters,
109                RH_max=self.rh_max,
110                success=self.success,
111                cell_order=self.cell_order,
112            )
113            if not self.success.all():
114                raise RuntimeError("Condensation failed")
115            if not self.update_thd:
116                self.particulator.environment.get_predicted("thd").ravel(
117                    self.particulator.environment.get_thd()
118                )
119            # note: this makes order of dynamics matter
120            #       (e.g., condensation after chemistry or before)
121            self.particulator.update_TpRH()
122
123            if self.adaptive:
124                self.counters["n_substeps"][:] = np.maximum(
125                    self.counters["n_substeps"][:],
126                    int(self.particulator.dt / self.dt_cond_range[1]),
127                )
128                if self.dt_cond_range[0] != 0:
129                    self.counters["n_substeps"][:] = np.minimum(
130                        self.counters["n_substeps"][:],
131                        int(self.particulator.dt / self.dt_cond_range[0]),
132                    )