open-atmos / PySDM

Pythonic particle-based (super-droplet) warm-rain/aqueous-chemistry cloud microphysics package with box, parcel & 1D/2D prescribed-flow examples in Python, Julia and Matlab
https://open-atmos.github.io/PySDM/
GNU General Public License v3.0
56 stars 32 forks source link

G&P example fails after turning collisions on (IndexError: index 199 is out of bounds for axis 0 with size 199) #1227

Closed huangynj closed 3 months ago

huangynj commented 9 months ago

Hello,

I tested the example Grabowski_and_Pawlowska_2023, but encountered the following issue.

In the setting.py, I change displacement: float = 1000 * si.m to displacement: float = 2000 * si.m

In the simulation.py, I add builder.add_dynamic(Coalescence(collision_kernel=Golovin(b=1.5e3 / si.second), adaptive=False))

Then run the initial part of figure_1.ipynb:

import numpy as np
import os
from matplotlib import pyplot
from matplotlib import ticker
from open_atmos_jupyter_utils import show_plot

from PySDM import Formulae
from PySDM.physics import si
from PySDM.products import (
    ParcelDisplacement, AmbientRelativeHumidity
)
TRIVIA = Formulae().trivia

from PySDM_examples.Grabowski_and_Pawlowska_2023 import Settings, Simulation

products=(
    ParcelDisplacement(name='z'),
    AmbientRelativeHumidity(name='S_max', unit='%', var='RH'),
)

vertical_velocity = 0.25 * si.m / si.s

Simulation(
    Settings(
        vertical_velocity=vertical_velocity,
        dt=1*si.s,
        n_sd=200,
        aerosol="pristine"
    ), products=products).run()

However, I encountered the following error.

Traceback (most recent call last):
  File “PySDM/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/Parcel.py", line 103, in <module>
    products=products).run()
                       ^^^^^
  File "PySDM/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/simulation.py", line 118, in run
    output_products = super()._run(
                      ^^^^^^^^^^^^^
  File "PySDM/examples/PySDM_examples/utils/basic_simulation.py", line 20, in _run
    self._save(output)
  File "PySDM/examples/PySDM_examples/Grabowski_and_Pawlowska_2023/simulation.py", line 114, in _save
    attr[drop_id].append(attr_data[drop_id])
                         ~~~~~~~~~^^^^^^^^^
IndexError: index 199 is out of bounds for axis 0 with size 199

Does it mean the coalescence process results in removing super-droplets? So that the total size of super-droplets reduces from 200 to 199? If yes, how to know which super-droplet was removed? Is there unique drop_id for each super-droplet?

Thanks,

Yongjie

AgnieszkaMakulska commented 9 months ago

Hello,

Yes, from what I understand the coalescence process results in removing super-droplets. It was not considered in this example because it concerned only condensation.

The code works when I change the _save function in simulation.py. Instead of

           for drop_id in range(self.particulator.n_sd): 
                attr[drop_id].append(attr_data[drop_id])

I wrote

            for drop_id in range(self.particulator.n_sd):
                if drop_id < len(attr_data):
                    attr[drop_id].append(attr_data[drop_id])
                else:
                    attr[drop_id].append(np.nan)

Then the length of attributes remains the same, but attributes of the "missing" super-droplets are np.nan. Maybe there is a more elegant way to do it. I'm not sure if it is possible to have a unique drop_id for each super-droplet, perhaps @slayoo can answer this.

slayoo commented 9 months ago

Thank you @AgnieszkaMakulska ! Upon coalescence or breakup, in a way the super-droplets inherently loose their "identity", and we do not have (yet) any drop-labeling attribute.

Sorting by one of the attributes before plotting could be a way to identify where, in a given attribute dimension, super-droplet are removed from the system?

slayoo commented 9 months ago

@huangynj, please confirm if that helps in your case, so we can add it to the main branch (and add a test with your sample code). Thanks

huangynj commented 9 months ago

@AgnieszkaMakulska @slayoo yes, filling in with np.nan can temporarily solve the issue. Thanks!

github-actions[bot] commented 6 months ago

Stale issue message

slayoo commented 3 months ago

Let me close this one as addressed. Thanks @huangynj for reporting, thanks @AgnieszkaMakulska for help!