Seeding examples for 1D environment of S&H 2012 as well as 2D environment of M&G2007

[jtbuch]

Update toward solving a subset of goals raised in #1387.

Specifically, I included functionality for selecting cell id, cell origin, and position in cell of the seeded super-droplets (SDs) in Shipway and Hill's (SH2012) kinematic driver environment. For this I had to modify seeding.py, particulator.py, as well as the backend seeding_methods.py files in a way that does not seem to be straightforwardly backward compatible with the parcel model example that @slayoo introduced earlier.

I also made another major assumption: rather than simply using a default value of water mass for the seeded SDs, I update it using the volume attribute derived from the wet, rather than dry, radius of seed SDs at the time of injection. If I use the dry radius, the resulting precipitation is significantly lower. This might be a subtle point that is worth discussing in some detail, even if it is for my own personal clarification.

Here's the .ipynb file containing the full example.

[claresinger]

Some initial comments.

Something looks odd to me in the plot of superdroplet count over time for the seeding vs. no seeding case. Why in the seeding case is the sd count constant at 32 during the seeding interval?

[claresinger]

should we use equilibrate_wet_radii rather than prescribing a water_mass for the seeds?

[jtbuch]

Something looks odd to me in the plot of superdroplet count over time for the seeding vs. no seeding case. Why in the seeding case is the sd count constant at 32 during the seeding interval?

That is due to the int type casting while plotting the mean SD count. I have fixed it now and it looks better here

However, this made me wonder: why are the number of SDs decreasing with time initially even in the no seeding case?

[claresinger]

However, this made me wonder: why are the number of SDs decreasing with time initially even in the no seeding case?

This is a good question... I thought it was just precip, but the sd_count is decreasing even before precip onset (and actually faster in the initial time than during the precip period). Is there some spin-up that's happening? I don't know why that wouldn't conserve SDs though.

[slayoo]

However, this made me wonder: why are the number of SDs decreasing with time initially even in the no seeding case?

This is a good question... I thought it was just precip, but the sd_count is decreasing even before precip onset (and actually faster in the initial time than during the precip period). Is there some spin-up that's happening? I don't know why that wouldn't conserve SDs though.

I expect it to be caused by the updraft which is only present during the first settings.t_1 seconds (10 minutes). It pushes super droplets beyond the top of the domain.

[claresinger]

I expect it to be caused by the updraft which is only present during the first settings.t_1 seconds (10 minutes). It pushes super droplets beyond the top of the domain.

Oh, this totally makes sense. Thanks, Sylwester!

[jtbuch]

Hi @slayoo and @claresinger, I have committed examples for seeding in both 1D and 2D environments that are compatible with the latest PySDM version. Any comments or criticisms would be much appreciated!

[claresinger]

Is this the same as the kinematic_1d.py file in the PySDM environments? This seems like a lot of code duplication. What have you added? Can it just extend this environment instead of copy/pasting?

[claresinger]

This may be a useful function to have outside of the seeding example. Should it go somewhere common like here? https://github.com/open-atmos/PySDM/tree/main/PySDM/exporters

[claresinger]

Can you write this so it extends the ShipwayHill settings.py? See how it is done here settings1D.py.

[claresinger]

plot time on the x-axis in last cell figures

[claresinger]

plot time on x-axis in last cell figures

[claresinger]

can/should a simulation_common.py file be used for the parcel (0D), 1D, and 2D cases and have all the common ingredients about Dynamics and Products and initialisation?

[jtbuch]

Hi @slayoo , @claresinger and I were discussing a potential issue with the Seeding dynamic as implemented. As you see in this attached plot for 1D seeding (it's also the case in the parcel model) below, the number of super-droplets does not change as expected. Specifically, when the seeding dynamic is called, instead of adding the seed super-droplets to the existing super-droplet buffer, it adds it to the initial buffer state. When it's done, the number of super-droplets returns to the (decreased) number of super-droplets with the additional seed ones.

I think this is due to the way attributes are initialized in the __call()__ method here:

PySDM/PySDM/dynamics/seeding.py

Line 72 in 07512fe

self.post_register_setup_when_attributes_are_known()

Unclear, how tricky the fix might be.

[jtbuch]

I see your point, @slayoo . I have now removed the update step for cell id etc. but still continue to see the same issue of the whole buffer getting updated when seeding.

[slayoo]

I see your point, @slayoo . I have now removed the update step for cell id etc. but still continue to see the same issue of the whole buffer getting updated when seeding.

I must say, I don't yet understand what you described. How does it show up in the parcel model? In the "hello-world" seeding example (https://github.com/open-atmos/PySDM/blob/main/examples/PySDM_examples/seeding/hello_world.ipynb) we clearly see an increase in SD count.

[jtbuch]

I must say, I don't yet understand what you described. How does it show up in the parcel model? In the "hello-world" seeding example (https://github.com/open-atmos/PySDM/blob/main/examples/PySDM_examples/seeding/hello_world.ipynb) we clearly see an increase in SD count.

Yes, sorry, this is not an issue in the parcel model (I incorrectly wrote that earlier) since the number of super-droplets does not change between the start of the simulation and the seeding timestep. This assumption does not translate to 1D and 2D where the number of super-droplets decreases as they are pushed outside the domain.

[jtbuch]

@slayoo I verified the increase in super-droplet count by calculating the total number of super-droplets for each simulation timestep in the 1D example. The count starts at 2336, reduces to 2059 before increasing to 2337 after adding 1 super-droplet. It keeps increasing to 2395 till the end of the seeding period before sharply reverting to 1881.

[slayoo]

doing

np.random.seed(123)
...
formulae= Formulae(seed= np.random.randint(1000)),

seems misleading - would it be more readable to pass the 123 to Formulae ctor?

[slayoo]

since it is a non-SI unit, suggest renaming the var to cloud water mixing ratio [g/kg]

[slayoo]

ditto

[slayoo]

since the unit is non-SI, suggest renaming to "RH [%]"

[slayoo]

since the unit is non-SI, suggest renaming to "... [%]"

[slayoo]

as we have multiple simulations in one notebook, sharing a single backend instance across simulations could significantly reduce runtime by avoiding JIT-compiling the same code multiple times.

[slayoo]

the realisation spread is definitely worth checking here - how the relevant output fluctuates across simulations with the same settings but different random seeds - sharing backend across simulations will speed things up then as well

[slayoo]

@slayoo I verified the increase in super-droplet count by calculating the total number of super-droplets for each simulation timestep in the 1D example. The count starts at 2336, reduces to 2059 before increasing to 2337 after adding 1 super-droplet. It keeps increasing to 2395 till the end of the seeding period before sharply reverting to 1881.

@jtbuch Just confirming I've reproduced it locally and I'm trying to figure out what is the origin of these super-droplet-count jumps

[slayoo]

@jtbuch, not necessarily related to the sd-count-per-cell issue, but it seems that as of now:

• position attributes of seed droplets are initialised to zeros
• the Displacement dynamic still advects these particles
• as a result these are seeded to the sixth cell (because the first 5 minutes are enough to move them up by six cells)

Likely not intentional. We should probably modify the displacement dynamic to ignore particles with zero multiplicity, and then to set different than zero position attributes to spawn the seed particles where intended.

[slayoo]

as for the SD-count jumps, I'm inclined to blame the incompatibility of seeding with the cell_start-computing logic here:

PySDM/PySDM/backends/impl_numba/methods/collisions_methods.py

Lines 685 to 697 in 36811d9

	def _counting_sort_by_cell_id_and_update_cell_start(
	new_idx, idx, cell_id, cell_idx, length, cell_start
	):
	cell_end = cell_start
	# Warning: Assuming len(cell_end) == n_cell+1
	cell_end[:] = 0
	for i in range(length):
	cell_end[cell_idx[cell_id[idx[i]]]] += 1
	for i in range(1, len(cell_end)):
	cell_end[i] += cell_end[i - 1]
	for i in range(length - 1, -1, -1):
	cell_end[cell_idx[cell_id[idx[i]]]] -= 1
	new_idx[cell_end[cell_idx[cell_id[idx[i]]]]] = idx[i]

We need to start with writing a unit test for multi-cell seeding to isolate the problem within a single step and a simplest possible setup (e.g., two cells, two super droplets, etc) - contributions welcome

[jtbuch]

Thank you, @claresinger and @slayoo for your helpful review! I will address these now and then respond to some of the major comments that @slayoo raised.