Multiplicity allocation across the droplet size spectrum.

[jbarr444]

I started a conversation with Clare Singer the other day and was curious if anyone else had any thoughts.

In short, what is the best way to initialize super droplet sizes and multiplicities?

I am using between 2^12 and 2^14 super droplets in the parcel environment, dv=1000 m^3. I have been using constant multiplicity as implemented in PySDM the whole time until recently experimenting with constant (multiplicity) * r^p; i.e. constant pth moment (p=1 constant length, p=2 constant area, etc.). It seems to really change things for p=3 for example when using just a single lognormal mode around ~.12 micron. The collision sampling frequency increases on the large tail of the distribution, decreases for the small tail, and the maximum and minumum resolved radii both increase, all leading I think to faster droplet growth over time.

Am I simply not using enough super droplets?

[slayoo]

There are several papers that explore different sampling strategies for super droplets, e.g.:

• Unterstrasser et al. 2017: https://doi.org/10.5194/gmd-10-1521-2017
• Dziekan & Pawlowska: https://doi.org/10.5194/acp-17-13509-2017
• Matsushima et al. 2023: https://doi.org/10.5194/gmd-16-6211-2023

[jbarr444]

Checking the pairs, this does not seem to be the case. Pairings between the two populations seems fairly common; I suppose they're just not passing the collision check for some reason.

[jbarr444]

Here's something: the problem can be reproduced with the one-liner:
self.particulator.attributes["water mass"].data = (self.particulator.attributes['volume'].data - self.particulator.attributes['dry volume'].data) * 997.

I hope this is illuminating to you but I am still not sure...

[jbarr444]

Fixed. The thing wanted a list, not a numpy array. I am on CPU so I thought numpy was ok, but evidently not.

[slayoo]

Well, it's hard to comment not knowing what the "thing" is :)

[jbarr444]

Sorry, I was pretty tired last night :^) The thing would be the one liner sent above, which actually does work as is. If the right side of the = is a numpy array, it doesn’t like it. Sorry for the confusion.

[emmacware]

Hey Jason, Sylwester pointed me to your discussion thread. We've been doing some work lately looking at a different part of the SDM algorithm, but looking at results across different sampling methods, and he thought I might find this interesting (he's right).

I'm very curious what you mean by constant expected collision rate - how are you flattening that from a value calculated in pairs to define a single superdroplet value? How does kernel choice play into that? How does the resampling handle other attributes (such as dry mass), or other attributes that potentially have interesting affects on the way the size spectrum behaves?

I don't know if you're familiar with PyPartMC, another particle based model developed by Sylwester and the open-atmos community geared towards aerosol representation. It has a different coalescence algorithm (the Weighted Flow Algorithm) as opposed to SDM, which you might find interesting (described in DeVille et al., 2011). Instead of having multiplicities, weights are a function of size and are not an attribute of the droplets. This changes the update process in many ways of course but the key thing that "constant expected collision rate" reminded me of is that in WFA, collision rates of particles of a certain size can be precomputed because the multiplicity/weighting of that is a constant, and not particle specific. It allows for some other interesting optimization as well, because of this they are able to do an optimized linear sampling rather than the random pairs. Whether or not its similar to what you've done here, I think there are some interesting analogies there!

Secondly, I'm interested in your resampling method for other reasons: I'm trying to do analysis on a coagulation scheme in the 2D environment. The model runs for a set "SpinUp" time to equilibrate with just condensation and advection before processes like Coagulation/Sedimentation are turned on. The idea is that this gives a chance for the Kohler stiffness/unrealistic prescribed initial supersaturation to even out, but I would love to sanity check that the differences in sampling methods are only through the Coagulation process alone (as expected), and not from anything happening in the SpinUp.

Love the visualizations! Such a cool way to look at droplet growth evolution.

[jbarr444]

Hi Emma,

I'm glad you think the method might be useful!

To begin answering your first question, start with eq. 84 from Shima et. al. 2020. The only logical step I took from there is to let

p(r) ∝ n(r) ∫K(r, r') n(r') dr'

where K can be any non-negative kernel, but preferably a collision kernel. I like to use notation from functional analysis to write

p(r) = (n (K n)) / <n, K n>

where < , > is the usual inner product on L^2.

Ok, this is a nice continuous equation, but as I said above, we're applying this to an already discretized system of super droplets, so they're going to have to be our best estimate of n(r). Now the notation becomes useful:

p(i) = (m_i (K m)) / <m, K m>

Now K is the kernel matrix evaluated on the super droplet radii, and m is the multiplicity vector. Notice now I'm choosing the index i, not r. This answers your third question. An index is chosen with this probability distribution, and the resulting offspring super droplet inherits all the attributes from its ancestor, the ith super droplet. To make this deterministic you would multiply by N_sd and convert to integers.

The main issue you have to deal with is what happens when N_sd(i) < 1. For this I find it helpful to actually just aggregate all the droplets with identical radii and draw r instead of i, then draw i uniformly from the SDs with that r for the other attributes. Then it gets hairy. I do a process that's like 'drop the rarest radius and redistribute its multiplicity to all the other radii with N_sd(r) < 1, repeat until all have N_sd(r) >=1.' This can be improved, but it works.

I have not heard of PyPartMC; I'd have to take some time to familiarize myself with it, but thanks for pointing it out! I like the idea of your SpinUp framework! Let me know if this is helpful, and if you have any followup questions!

[emmacware]

Very interesting, thanks for the explanation of the algorithm! I'll have to let this sink for a while-- I think I mostly understand what you are doing with creating the collisions probability pdf, but I don't have much intuition for the implications beyond that. How unwieldy does the matrix multiplication get at large system Ns? Seems pretty analogous to the way binned models can optimize how to move mass.

I'm still a little confused by the resampling step -- does it only inherit non radius attributes and the radius is chosen based on equal pdf areas, or is the value of the pdf chosen at radius of superdroplet i, and inheriting everything including size attributes (which might not be equal collision rates)?

New questions are:

1. How long do the superdroplets retain the benefits of equal collision rates before the sampling?
2. I can see how it conserves number concentration, but I can't tell what happens to the mass distribution with the resampling, and how does it work to conserve total mass (or distribution of other attributes) if the inherited attributes like radius no longer match the same population size? (I could be confusing the method still)
3. I'm a little surprised you're getting multiplicities under one -- for a superdroplet with such a small multiplicity to have the same collision "activity" as everything else they must be huge, right? If thats the case my instinct is that these are probably important collector drops for the system (Lucky droplets).
4. Are the benefits of equal collision rates made redundant by still checking linear pairs? Is this an opportunity to check less given equal collision rates (I have an old project on sublinear sampling at a bit of a standstill currently)

The analogy to WFA is pretty strong, now that I understand it a little better: different algorithm, but they are able to choose and scale pairs so that everything tested the same collision rate (within a small range), to maximize the ratio of successful events to pairs tested. I got very excited about this when I learned about it and the overlap it might have with my sublinear project, but havent been able to find a way to apply it to SDM with the uneven multiplicity spread, and the need to know the whole system before making being able to make optimized sampling choices.

Very interested to hear your thoughts on these!

[jbarr444]

0. So the radius probability distribution is fully discrete, and only takes nonzero values at the prior SD radii. You can inherit all non-multiplicity attributes from a single SD, or you can draw radius separately, which I prefer. I get that the distinction might be confusing, but the latter helps boost low probabilities by aggregating them from like-radii SDs, but is slower. If all radii are initially unique then there's no difference and you are just cloning SDs. I think I neglected to emphasize that all attributes except multiplicity are inherited from previous SDs. Multiplicity is recalculated according to the Shima equation 84. This is how n(r) is conserved to a great approximation.

1. I would think (but have not tested) that the explosive nature of collisional growth would make it so that you can't count on each droplet having similar rates beyond a few time steps, unless you resample again. I'm toying with the idea of resampling more frequently, but I don't know anything about that yet.

2. Not only is ∫n(r)dr conserved, but n(r) is the same for each r; therefore all moments are (approximately, within rounding) conserved. In theory anyway. In practice, the N_sd(r) < 1 problem forces me to make decisions about which moments to conserve best (number, mass).

3. Not multiplicity < 1; number of super droplets assigned to a radius value < 1 is the issue.

4. My motivation for coming up with this is not to make any process cheaper to compute, just to make the simulation more accurate for a little bit more of a computational price. If that is what you are thinking of trying with it, I'd be very interested to hear what you come up with!

Would you like to give it a try?

[emmacware]

Ahh I was misunderstanding before, I see what you are saying about the unique set, and point 3 that number of superdroplets assigned is the issue, not multiplicity <1.

I'd love to try it out, I feel like that would give me more intuition with the algorithm!

[jbarr444]

Cool! I am not experienced in sharing code over the internet. How would you prefer I send you the function?

[slayoo]

@jbarr444, git and GitHub support multiple branches for each repository, so one way of sharing it would be to commit these changes to a separate branch and make it public by pushing these commits to your fork of PySDM. The commands could look like:

• create a new branch, e.g. with: git checkout -b resampling
• commit the changes to this branch, e.g. with: git commit -m 'COMMIT DESCRIPTION' file.py
• push these changes to your fork, e.g. with: git push --set-upstream origin resampling