Condensation failed in parcel environment

[mikhailmints]

For some input parameters, the condensation solver is failing in the parcel environment. Here is a minimal example using PySDM_examples.Pyrcel:

from PySDM import Formulae
from PySDM.initialisation.spectra import Lognormal
from PySDM.physics import si

from PySDM_examples.Pyrcel import Settings, Simulation

settings = Settings(
    dz = 1 * si.m,
    n_sd_per_mode = (1000,),
    aerosol_modes_by_kappa = {
        0.1885611585360898: Lognormal(
            norm_factor=40233789.92354 / si.m ** 3,
            m_mode=1.1669322087438114e-09 * si.m,
            s_geom=1.7322676373747328
        ),
    },
    vertical_velocity = 0.0143756504108188 * si.m / si.s,
    initial_pressure = 76309.09460546204 * si.Pa,
    initial_temperature = 296.50590490254416 * si.K,
    initial_relative_humidity = 1,
    displacement = 1000 * si.m,
    formulae = Formulae(constants={'MAC': 1})
)

simulation = Simulation(settings, products=[])

results = simulation.run()

Produces the following output (first and last part shown):

failed to find interval
    file: /home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/backends/impl_numba/methods/condensation_methods.py
    context:
         T
         296.50096771736526
         p
         76304.6452879693
         RH
         1.0002393891789032
         a
         -65.6888252946611
         b
         -70.29399548064919
         fa
         -5445.374951706045
         fb
         -146020.70967949124
failed to find interval
    file: /home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/backends/impl_numba/methods/condensation_methods.py
    context:
         T
         296.50096771736173
         p
         76304.64528796851
         RH
         1.0002393891791526
         a
         -65.6888252946611
         b
         -70.29399548064919
         fa
         -5445.374951705585
         fb
         -146020.7096794852
failed to find interval
    file: /home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/backends/impl_numba/methods/condensation_methods.py
    context:
         T
         296.50096771736173
         p
         76304.64528796851
         RH
         1.0002393891791526
         a
         -65.6888252946611
         b
         -70.29399548064919
         fa
         -5445.374951705585
         fb
         -146020.7096794852

...

failed to find interval
    file: /home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/backends/impl_numba/methods/condensation_methods.py
    context:
         T
         296.50096771736173
         p
         76304.64528796851
         RH
         1.0002393891791526
         a
         -65.6888252946611
         b
         -70.29399548064919
         fa
         -5445.374951705585
         fb
         -146020.7096794852
failed to find interval
    file: /home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/backends/impl_numba/methods/condensation_methods.py
    context:
         T
         296.50096771736173
         p
         76304.64528796851
         RH
         1.0002393891791526
         a
         -65.6888252946611
         b
         -70.29399548064919
         fa
         -5445.374951705585
         fb
         -146020.7096794852
Traceback (most recent call last):
  File "/central/home/mmints/CliMA/AerosolActivationEmulation/CondFailedIssueMinExample.py", line 27, in <module>
    results = simulation.run()
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM_examples/Pyrcel/simulation.py", line 85, in run
    output_products = super()._run(
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM_examples/utils/basic_simulation.py", line 19, in _run
    self.particulator.run(steps=steps_per_output_interval)
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/particulator.py", line 48, in run
    dynamic()
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/dynamics/condensation.py", line 104, in __call__
    raise RuntimeError("Condensation failed")
RuntimeError: Condensation failed

I tried changing rtol_thd, rtol_x, and dt_cond_range to different values, but the error still appeared. In this particular example, changing the number of superdroplets to 100 instead of 1000 makes the error disappear. However, if Logarithmic is used for the spectral sampling instead of ConstantMultiplicity, the issue appears even for low numbers of superdroplets.

[slayoo]

@mikhailmints thank you for reporting it!

It is very likely that #868 solves the problem - it is a PR in which we switch from theta to supersaturation as a driver for timestep adaptivity, what helps to stabilise the solver behaviour further away from the activation point (and in your logs the RH is almost =1, which hints that this might be the case). Help welcome in reviewing and making that PR mergeable.

As a temporary workaround, you can either try disabling timestep adaptivity when instantiating the Condensation dynamic or try using SciPy ODE solver instead of the PySDM internal one using the scipy_solver argument for Pyrcel.Simulation.__init__.

HTH

[mikhailmints]

Thanks for the response! Disabling timestep adaptivity still seems to give the same error. If I use the SciPy solver, I get the following:

 lsoda--  at t (=r1), too much accuracy requested  
       for precision of machine..  see tolsf (=r2) 
      in above,  r1 =  0.2523173154976D-04   r2 =                  NaN
/home/mmints/anaconda3/lib/python3.10/site-packages/scipy/integrate/_ode.py:1348: UserWarning: lsoda: Excess accuracy requested (tolerances too small).
  warnings.warn('{:s}: {:s}'.format(self.__class__.__name__,
Traceback (most recent call last):
  File "/central/home/mmints/CliMA/AerosolActivationEmulation/CondFailedIssueMinExample.py", line 30, in <module>
    results = simulation.run()
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM_examples/Pyrcel/simulation.py", line 85, in run
    output_products = super()._run(
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM_examples/utils/basic_simulation.py", line 19, in _run
    self.particulator.run(steps=steps_per_output_interval)
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/particulator.py", line 48, in run
    dynamic()
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/dynamics/condensation.py", line 95, in __call__
    self.particulator.condensation(
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py", line 32, in _condensation
    func(
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/backends/impl_numba/methods/condensation_methods.py", line 137, in _condensation
    ) = solver(
  File "/home/mmints/anaconda3/lib/python3.10/site-packages/PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py", line 217, in solve
    assert integ.success, integ.message
AssertionError: Unexpected istate in LSODA.

Changing rtol_x and rtol_thd again does not seem to affect anything.

[mikhailmints]

I also tried this with the current version of #868, and I am still getting exactly the same errors.

[slayoo]

Thank you for the follow ups. Indeed seems it is not related in any way with condensation numerics in PySDM (as the same problem arises with SciPy solver).

I confirm I've reproduced it locally.

Seems the problem arises at the very beginning (temperature in the warning logs barely differs from the initial temperature). Setting a lower initial RH could help here, in principle the aerosol initialisation (searching for equilibrium wet radii) is well posed only for RH<1.

[mikhailmints]

It looks like in the SciPy solver specifically, decreasing the value of rtol does help - for example, I changed the value to 1e-8 in #1106, and the test failures that looked similar to this issue disappeared (except now some tests now time out). The example I provided above now also works with the SciPy solver, although if the sampling is changed to Logarithmic, it still fails (in the middle of the simulation, not in the very beginning like before). With the normal condensation solver, however, I still couldn't find any way to make it work other than reducing the number of superdroplets. Changing the initial RH also does not help.

[github-actions[bot]]

Stale issue message

[github-actions[bot]]

Stale issue message

[github-actions[bot]]

Stale issue message