Skip to content

Commit e4ddacd

Browse files
emmacwareslayootluettm
authored
super-droplet attribute "alpha" sampling (incl. new example based on Fig 1 from Matsushima et al. 2023) + refactor of the spectral sampling logic to enable deterministic/pseudorandom/quasirandom choice for each sampling type (#1651)
Co-authored-by: Sylwester Arabas <sylwester.arabas@agh.edu.pl> Co-authored-by: Tim Lüttmer <154344756+tluettm@users.noreply.github.com>
1 parent a14c605 commit e4ddacd

File tree

53 files changed

+16504
-186
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

53 files changed

+16504
-186
lines changed

‎PySDM/environments/kinematic_1d.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -60,12 +60,12 @@ def init_attributes(
6060
) = self.mesh.cellular_attributes(positions)
6161

6262
if collisions_only:
63-
v_wet, n_per_kg = spectral_discretisation.sample(
63+
v_wet, n_per_kg = spectral_discretisation.sample_deterministic(
6464
backend=self.particulator.backend, n_sd=self.particulator.n_sd
6565
)
6666
attributes["volume"] = v_wet
6767
else:
68-
r_dry, n_per_kg = spectral_discretisation.sample(
68+
r_dry, n_per_kg = spectral_discretisation.sample_deterministic(
6969
backend=self.particulator.backend, n_sd=self.particulator.n_sd
7070
)
7171
attributes["dry volume"] = self.formulae.trivia.volume(radius=r_dry)

‎PySDM/environments/kinematic_2d.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -60,9 +60,9 @@ def init_attributes(
6060
attributes["position in cell"],
6161
) = self.mesh.cellular_attributes(positions)
6262

63-
r_dry, n_per_kg = spectral_sampling(spectrum=dry_radius_spectrum).sample(
64-
n_sd=n_sd, backend=self.particulator.backend
65-
)
63+
r_dry, n_per_kg = spectral_sampling(
64+
spectrum=dry_radius_spectrum
65+
).sample_deterministic(n_sd=n_sd, backend=self.particulator.backend)
6666

6767
attributes["dry volume"] = self.formulae.trivia.volume(radius=r_dry)
6868
attributes["kappa times dry volume"] = kappa * attributes["dry volume"]

‎PySDM/initialisation/sampling/spectral_sampling.py‎

Lines changed: 108 additions & 76 deletions
Original file line numberDiff line numberDiff line change
@@ -1,58 +1,51 @@
11
"""
2-
spectral sampling logic incl. linear, logarithmic, uniform-random and constant-multiplicity
3-
sampling classes
2+
spectral discretisation logic incl. linear, logarithmic, and constant-multiplicity
3+
layouts with deterministic, pseudorandom and quasirandom sampling
44
"""
55

6+
import abc
67
from typing import Optional, Tuple
78

89
import numpy as np
9-
from scipy import optimize
10+
from scipy.interpolate import interp1d
1011

1112
default_cdf_range = (0.00001, 0.99999)
1213

1314

14-
class SpectralSampling: # pylint: disable=too-few-public-methods
15-
def __init__(self, spectrum, size_range: Optional[Tuple[float, float]] = None):
15+
class SpectralSampling:
16+
def __init__(
17+
self,
18+
spectrum,
19+
*,
20+
size_range: Optional[Tuple[float, float]] = None,
21+
error_threshold: Optional[float] = None,
22+
):
1623
self.spectrum = spectrum
24+
self.error_threshold = error_threshold or 0.01
1725

1826
if size_range is None:
19-
if hasattr(spectrum, "percentiles"):
20-
self.size_range = spectrum.percentiles(default_cdf_range)
21-
else:
22-
self.size_range = [np.nan, np.nan]
23-
for i in (0, 1):
24-
result = optimize.root(
25-
lambda x, value=default_cdf_range[i]: spectrum.cdf(x) - value,
26-
x0=spectrum.median(),
27-
)
28-
assert result.success
29-
self.size_range[i] = result.x
27+
self.cdf_range = default_cdf_range
28+
self.size_range = spectrum.percentiles(self.cdf_range)
3029
else:
3130
assert len(size_range) == 2
3231
assert size_range[0] > 0
3332
assert size_range[1] > size_range[0]
3433
self.size_range = size_range
34+
self.cdf_range = (
35+
spectrum.cdf(size_range[0]),
36+
spectrum.cdf(size_range[1]),
37+
)
3538

39+
@abc.abstractmethod
40+
def _sample(self, frac_values):
41+
pass
3642

37-
class DeterministicSpectralSampling(
38-
SpectralSampling
39-
): # pylint: disable=too-few-public-methods
40-
# TODO #1031 - error_threshold will be also used in non-deterministic sampling
41-
def __init__(
42-
self,
43-
spectrum,
44-
size_range: Optional[Tuple[float, float]] = None,
45-
error_threshold: Optional[float] = None,
46-
):
47-
super().__init__(spectrum, size_range)
48-
self.error_threshold = error_threshold or 0.01
49-
50-
def _sample(self, grid, spectrum):
43+
def _sample_with_grid(self, grid):
5144
x = grid[1:-1:2]
52-
cdf = spectrum.cumulative(grid[0::2])
45+
cdf = self.spectrum.cumulative(grid[0::2])
5346
y_float = cdf[1:] - cdf[0:-1]
5447

55-
diff = abs(1 - np.sum(y_float) / spectrum.norm_factor)
48+
diff = abs(1 - np.sum(y_float) / self.spectrum.norm_factor)
5649
if diff > self.error_threshold:
5750
raise ValueError(
5851
f"{diff * 100:.3g}% error in total real-droplet number due to sampling "
@@ -61,61 +54,100 @@ def _sample(self, grid, spectrum):
6154

6255
return x, y_float
6356

57+
def sample_deterministic(
58+
self, n_sd, *, backend=None
59+
): # pylint: disable=unused-argument
60+
return self._sample(
61+
frac_values=np.linspace(
62+
self.cdf_range[0], self.cdf_range[1], num=2 * n_sd + 1
63+
)
64+
)
6465

65-
class Linear(DeterministicSpectralSampling): # pylint: disable=too-few-public-methods
66-
def sample(self, n_sd, *, backend=None): # pylint: disable=unused-argument
67-
grid = np.linspace(*self.size_range, num=2 * n_sd + 1)
68-
return self._sample(grid, self.spectrum)
66+
def sample_quasirandom(self, n_sd, *, backend):
67+
num_elements = n_sd
68+
storage = backend.Storage.empty(num_elements, dtype=float)
69+
backend.Random(seed=backend.formulae.seed, size=num_elements)(storage)
70+
u01 = storage.to_ndarray()
6971

72+
frac_values = np.linspace(
73+
self.cdf_range[0], self.cdf_range[1], num=2 * n_sd + 1
74+
)
7075

71-
class Logarithmic(
72-
DeterministicSpectralSampling
73-
): # pylint: disable=too-few-public-methods
74-
def __init__(
75-
self,
76-
spectrum,
77-
size_range: [None, Tuple[float, float]] = None,
78-
error_threshold: Optional[float] = None,
79-
):
80-
super().__init__(spectrum, size_range, error_threshold)
81-
self.start = np.log10(self.size_range[0])
82-
self.stop = np.log10(self.size_range[1])
76+
for i in range(1, len(frac_values) - 1, 2):
77+
frac_values[i] = frac_values[i - 1] + u01[i // 2] * (
78+
frac_values[i + 1] - frac_values[i - 1]
79+
)
8380

84-
def sample(self, n_sd, *, backend=None): # pylint: disable=unused-argument
85-
grid = np.logspace(self.start, self.stop, num=2 * n_sd + 1)
86-
return self._sample(grid, self.spectrum)
81+
return self._sample(frac_values=frac_values)
8782

83+
def sample_pseudorandom(self, n_sd, *, backend):
84+
num_elements = 2 * n_sd + 1
85+
storage = backend.Storage.empty(num_elements, dtype=float)
86+
backend.Random(seed=backend.formulae.seed, size=num_elements)(storage)
87+
u01 = storage.to_ndarray()
88+
89+
frac_values = np.sort(
90+
self.cdf_range[0] + u01 * (self.cdf_range[1] - self.cdf_range[0])
91+
)
92+
return self._sample(frac_values=frac_values)
8893

89-
class ConstantMultiplicity(
90-
DeterministicSpectralSampling
91-
): # pylint: disable=too-few-public-methods
92-
def __init__(self, spectrum, size_range=None):
93-
super().__init__(spectrum, size_range)
9494

95-
self.cdf_range = (
96-
spectrum.cumulative(self.size_range[0]),
97-
spectrum.cumulative(self.size_range[1]),
95+
class Logarithmic(SpectralSampling):
96+
def _sample(self, frac_values):
97+
grid = np.exp(
98+
(np.log(self.size_range[1]) - np.log(self.size_range[0])) * frac_values
99+
+ np.log(self.size_range[0])
98100
)
99-
assert 0 < self.cdf_range[0] < self.cdf_range[1]
101+
return self._sample_with_grid(grid)
100102

101-
def sample(self, n_sd, *, backend=None): # pylint: disable=unused-argument
102-
cdf_arg = np.linspace(self.cdf_range[0], self.cdf_range[1], num=2 * n_sd + 1)
103-
cdf_arg /= self.spectrum.norm_factor
104-
percentiles = self.spectrum.percentiles(cdf_arg)
105103

106-
assert np.isfinite(percentiles).all()
104+
class Linear(SpectralSampling):
105+
def _sample(self, frac_values):
106+
grid = self.size_range[0] + frac_values * (
107+
self.size_range[1] - self.size_range[0]
108+
)
109+
return self._sample_with_grid(grid)
107110

108-
return self._sample(percentiles, self.spectrum)
109111

112+
class ConstantMultiplicity(SpectralSampling):
113+
def _sample(self, frac_values):
114+
grid = self.spectrum.percentiles(frac_values)
115+
assert np.isfinite(grid).all()
116+
return self._sample_with_grid(grid)
110117

111-
class UniformRandom(SpectralSampling): # pylint: disable=too-few-public-methods
112-
def sample(self, n_sd, *, backend):
113-
n_elements = n_sd
114-
storage = backend.Storage.empty(n_elements, dtype=float)
115-
backend.Random(seed=backend.formulae.seed, size=n_elements)(storage)
116-
u01 = storage.to_ndarray()
117118

118-
pdf_arg = self.size_range[0] + u01 * (self.size_range[1] - self.size_range[0])
119-
dr = abs(self.size_range[1] - self.size_range[0]) / n_sd
120-
# TODO #1031 - should also handle error_threshold check
121-
return pdf_arg, dr * self.spectrum.size_distribution(pdf_arg)
119+
class AlphaSampling(SpectralSampling):
120+
"""as in [Matsushima et al. 2023](https://doi.org/10.5194/gmd-16-6211-2023)"""
121+
122+
def __init__(
123+
self,
124+
spectrum,
125+
*,
126+
alpha,
127+
dist_1_inv,
128+
interp_points,
129+
size_range=None,
130+
error_threshold: Optional[float] = None,
131+
):
132+
super().__init__(
133+
spectrum, size_range=size_range, error_threshold=error_threshold
134+
)
135+
self.alpha = alpha
136+
self.dist_0_cdf = self.spectrum.cdf
137+
self.dist_1_inv = dist_1_inv
138+
self.x_prime = np.linspace(
139+
self.size_range[0], self.size_range[1], num=interp_points
140+
)
141+
142+
def _sample(self, frac_values):
143+
if self.alpha == 0:
144+
frac_values = self.spectrum.percentiles(frac_values)
145+
elif self.alpha == 1:
146+
frac_values = self.dist_1_inv(frac_values, self.size_range)
147+
else:
148+
sd_cdf = self.dist_0_cdf(self.x_prime)
149+
x_sd_cdf = (1 - self.alpha) * self.x_prime + self.alpha * self.dist_1_inv(
150+
sd_cdf, self.size_range
151+
)
152+
frac_values = interp1d(sd_cdf, x_sd_cdf)(frac_values)
153+
return self._sample_with_grid(frac_values)

‎PySDM/initialisation/spectra/sum.py‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ def __init__(self, spectra: tuple, interpolation_grid=None):
2020
cdf_arg = np.zeros(len(interpolation_grid) * len(self.spectra) + 1)
2121
cdf_arg[1:] = np.concatenate(percentiles)
2222
cdf = self.cumulative(cdf_arg) / self.norm_factor
23-
self.inverse_cdf = interp1d(cdf, cdf_arg)
23+
self.inverse_cdf = interp1d(cdf, cdf_arg, bounds_error=False)
2424

2525
def size_distribution(self, arg):
2626
result = 0.0
@@ -36,3 +36,9 @@ def cumulative(self, arg):
3636

3737
def percentiles(self, cdf_values):
3838
return self.inverse_cdf(cdf_values)
39+
40+
def pdf(self, arg):
41+
return self.size_distribution(arg) / self.norm_factor
42+
43+
def cdf(self, arg):
44+
return self.cumulative(arg) / self.norm_factor

‎docs/bibliography.json‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -951,6 +951,15 @@
951951
"title": "A physically constrained classical description of the homogeneous nucleation of ice in water",
952952
"label": "Koop and Murray 2016 (J. Chem. Phys. 145)"
953953
},
954+
"https://doi.org/10.5194/gmd-16-6211-2023": {
955+
"title": "Overcoming computational challenges to realize meter- to submeter-scale resolution in cloud simulations using the super-droplet method",
956+
"label": "Matsushima et al. 2023 (Geosci. Model Dev. 16)",
957+
"usages": [
958+
"PySDM/initialisation/sampling/spectral_sampling.py",
959+
"examples/PySDM_examples/Matsushima_et_al_2023/figure_1.ipynb",
960+
"examples/PySDM_examples/Matsushima_et_al_2023/__init__.py"
961+
]
962+
},
954963
"https://doi.org/10.48550/arXiv.2509.05536": {
955964
"usages": [
956965
"examples/PySDM_examples/Ware_et_al_2025/__init__.py",

‎docs/markdown/pysdm_landing.md‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -103,7 +103,7 @@ Exponential = pyimport("PySDM.initialisation.spectra").Exponential
103103
n_sd = 2^15
104104
initial_spectrum = Exponential(norm_factor=8.39e12, scale=1.19e5 * si.um^3)
105105
attributes = Dict()
106-
attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(spectrum=initial_spectrum).sample(n_sd)
106+
attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(spectrum=initial_spectrum).sample_deterministic(n_sd)
107107
```
108108
</details>
109109
<details>
@@ -119,7 +119,7 @@ initial_spectrum = Exponential(pyargs(...
119119
'norm_factor', 8.39e12, ...
120120
'scale', 1.19e5 * si.um ^ 3 ...
121121
));
122-
tmp = ConstantMultiplicity(initial_spectrum).sample(int32(n_sd));
122+
tmp = ConstantMultiplicity(initial_spectrum).sample_deterministic(int32(n_sd));
123123
attributes = py.dict(pyargs('volume', tmp{1}, 'multiplicity', tmp{2}));
124124
```
125125
</details>
@@ -133,7 +133,7 @@ from PySDM.initialisation.spectra.exponential import Exponential
133133
n_sd = 2 ** 15
134134
initial_spectrum = Exponential(norm_factor=8.39e12, scale=1.19e5 * si.um ** 3)
135135
attributes = {}
136-
attributes['volume'], attributes['multiplicity'] = ConstantMultiplicity(initial_spectrum).sample(n_sd)
136+
attributes['volume'], attributes['multiplicity'] = ConstantMultiplicity(initial_spectrum).sample_deterministic(n_sd)
137137
```
138138
</details>
139139

@@ -323,7 +323,7 @@ The component submodules used to create this simulation are visualized below:
323323
IS["initial_spectrum :Exponential"] -->|passed as arg to| CM_INIT
324324
CM_INIT(["ConstantMultiplicity.__init__()"]) -->|instantiates| CM_INSTANCE
325325
CM_INSTANCE[":ConstantMultiplicity"] -.-|has a method| SAMPLE
326-
SAMPLE(["ConstantMultiplicity.sample()"]) -->|returns| n
326+
SAMPLE(["ConstantMultiplicity.sample_deterministic()"]) -->|returns| n
327327
SAMPLE -->|returns| volume
328328
n -->|added as element of| ATTRIBUTES
329329
PARTICULATOR_INSTANCE -.-|has a method| PARTICULATOR_RUN(["Particulator.run()"])
@@ -414,7 +414,7 @@ builder = Builder(backend=CPU(formulae), n_sd=n_sd, environment=env)
414414
builder.add_dynamic(AmbientThermodynamics())
415415
builder.add_dynamic(Condensation())
416416

417-
r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample(n_sd)
417+
r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample_deterministic(n_sd)
418418
v_dry = formulae.trivia.volume(radius=r_dry)
419419
r_wet = equilibrate_wet_radii(r_dry=r_dry, environment=builder.particulator.environment, kappa_times_dry_volume=kappa * v_dry)
420420

@@ -495,7 +495,7 @@ builder = Builder(pyargs('backend', CPU(formulae), 'n_sd', int32(n_sd), 'environ
495495
builder.add_dynamic(AmbientThermodynamics());
496496
builder.add_dynamic(Condensation());
497497
498-
tmp = spectral_sampling.Logarithmic(spectrum).sample(int32(n_sd));
498+
tmp = spectral_sampling.Logarithmic(spectrum).sample_deterministic(int32(n_sd));
499499
r_dry = tmp{1};
500500
v_dry = formulae.trivia.volume(pyargs('radius', r_dry));
501501
specific_concentration = tmp{2};
@@ -596,7 +596,7 @@ builder = Builder(backend=CPU(formulae), n_sd=n_sd, environment=env)
596596
builder.add_dynamic(AmbientThermodynamics())
597597
builder.add_dynamic(Condensation())
598598

599-
r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample(n_sd)
599+
r_dry, specific_concentration = spectral_sampling.Logarithmic(spectrum).sample_deterministic(n_sd)
600600
v_dry = formulae.trivia.volume(radius=r_dry)
601601
r_wet = equilibrate_wet_radii(r_dry=r_dry, environment=builder.particulator.environment, kappa_times_dry_volume=kappa * v_dry)
602602

‎examples/PySDM_examples/Abade_and_Albuquerque_2024/simulation.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -49,9 +49,9 @@ def __init__(self, settings):
4949
if settings.enable_vapour_deposition_on_ice:
5050
builder.add_dynamic(VapourDepositionOnIce(adaptive=True))
5151

52-
r_dry, n_in_dv = ConstantMultiplicity(settings.soluble_aerosol).sample(
53-
n_sd=settings.n_sd
54-
)
52+
r_dry, n_in_dv = ConstantMultiplicity(
53+
settings.soluble_aerosol
54+
).sample_deterministic(n_sd=settings.n_sd)
5555
attributes = builder.particulator.environment.init_attributes(
5656
n_in_dv=n_in_dv, kappa=settings.kappa, r_dry=r_dry
5757
)

‎examples/PySDM_examples/Abdul_Razzak_Ghan_2000/run_ARG_parcel.py‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -61,7 +61,9 @@ def run_parcel(
6161
}
6262
for i, mode in enumerate(aerosol.modes):
6363
kappa, spectrum = mode["kappa"]["CompressedFilmOvadnevaite"], mode["spectrum"]
64-
r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd_per_mode)
64+
r_dry, concentration = ConstantMultiplicity(spectrum).sample_deterministic(
65+
n_sd_per_mode
66+
)
6567
v_dry = builder.formulae.trivia.volume(radius=r_dry)
6668
specific_concentration = concentration / builder.formulae.constants.rho_STP
6769
attributes["multiplicity"] = np.append(

‎examples/PySDM_examples/Alpert_and_Knopf_2016/simulation.py‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -227,7 +227,9 @@ def simulation(
227227
n_sd, multiplicity / volume
228228
)
229229
else:
230-
_isa, _conc = spectral_sampling.ConstantMultiplicity(spectrum).sample(n_sd)
230+
_isa, _conc = spectral_sampling.ConstantMultiplicity(
231+
spectrum
232+
).sample_deterministic(n_sd)
231233
attributes = {
232234
"multiplicity": discretise_multiplicities(_conc * volume),
233235
"immersed surface area": _isa,

‎examples/PySDM_examples/Arabas_et_al_2025/aida.ipynb‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -301,7 +301,7 @@
301301
}
302302
},
303303
"source": [
304-
"d_dry, conc = spectral_sampling.ConstantMultiplicity(dNdD).sample(common['n_sd'])\n",
304+
"d_dry, conc = spectral_sampling.ConstantMultiplicity(dNdD).sample_deterministic(common['n_sd'])\n",
305305
"\n",
306306
"ATTRIBUTES = {\n",
307307
" 'n': conc * common['volume'],\n",

0 commit comments

Comments
 (0)