Skip to content

Bug: Atmosphere.from_pressure may fail to converge for pressures near the tropopause #15

@hiratch

Description

@hiratch

Environment:

ambiance version: 1.3.1  
Python version: 3.9.12

Problem Description:

When using Atmosphere.from_pressure(), a RuntimeError occurs if the input pressure is a value within the troposphere but very close to the tropopause boundary (22632.0 Pa). The scipy.optimize.newton solver fails to converge within the maximum number of iterations.

Steps to Reproduce:

The following minimal code snippet consistently reproduces the error.

from ambiance import Atmosphere

# This pressure is in the troposphere, but very close to the tropopause boundary.
pressure_pa = 22632.2

# This call fails with a RuntimeError.
inst = Atmosphere.from_pressure(pressure_pa)

Traceback:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/path/to/your/env/lib/python3.9/site-packages/ambiance/ambiance.py", line 264, in from_pressure
    return cls(h=opt.newton(f, x0=81e3 - 16e3*np.log10(p)))
  File "/path/to/your/env/lib/python3.9/site-packages/scipy/optimize/_zeros_py.py", line 368, in newton
    raise RuntimeError(msg)
RuntimeError: Failed to converge after 50 iterations, value is [11019.07348227].

Analysis (My Investigation & Thoughts):

I believe the issue might be related to the initial guess (x0) provided to the scipy.optimize.newton solver. Here is my reasoning:

  1. According to the ICAO standard atmosphere, the tropopause is at 11,000 m, where the pressure is approximately 22632.0 Pa. A pressure of 22632.2 Pa should therefore correspond to an altitude slightly below 11,000 m, placing it in the troposphere.
  2. The from_pressure method appears to use a generic formula for the initial guess: x0=81e3 - 16e3*np.log10(p).
  3. For p = 22632.2, this formula yields an initial guess x0 of approximately 11,324 meters. This altitude is located in the stratosphere.
  4. It seems the solver might be starting its search from an incorrect atmospheric layer. I suspect that because of the "kink" in the pressure-altitude function at the layer boundary, the solver struggles to converge when it has to cross this boundary from an initial guess in the wrong layer.

I hope this analysis is helpful for debugging. Thank you for maintaining this useful package.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions