NaN appears during the integration

Orekit Python Wrapper usage
1

Hi everyone,

I have recently started to work with Orekit wrapped in Python and I am currently facing a problem with a VLEO propagation that also considers air drag effect.
I have already opened an issue on the GitLab repository (NaN appears during the integration (#463) · Issues · Orekit Labs / Orekit Python Wrapper · GitLab) where you can find a precise description of the scenario, the error I am getting and a script you can use to test the case.

To sum it up, at a certain point of the propagation the run stops because NaN values appear. The error happens inside the propagation class, so I cannot debug it from Python.

Has anyone faced a similar problem? Can anyone help me?

Thank you in advance!

2

Hi and welcome,

I haven’t looked at your code yet but if you say it’s in VLEO it’s probably to do with atmospheric model at low altitude.
Does the NaN disappear if you remove drag from the force models or switch to a different atmospheric model?

Cheers,
Romain.

3

Hi Romain, thank you for the answer!

The NaNs disappear if I remove the drag, but they do not if I switch to a different atmospheric model. I tried NRLMSISE00 and DTM2000.

In this .txt file, you can see the altitudes in km along the run. The last value is the last before the NaNs appear (351.01363869109844 km).
altitude.txt (2.5 MB)

I link also here the script I am using to run the test if you want to try it directly.
orekit_free_propagation.py (8.9 KB)

Have a nice day,
Simone

4

Hi @Simo-Dalbo,

I’ve looked a bit into your code and I think it’s a bug, thanks for finding it !!

Here’s a Java class pinpointing the error:
Forum_2024_03_29_Simo_Dalbo_NaN_Appears_During_Integration2.java (5.4 KB)
With this the error occurs between 0.011s and 0.012s of propagation.
The error disappears when the drag model is not used.
Propagation starts at 2005-09-10T00:26:47.22 so about 20.5 days after your own initial date.
The error happens in Hipparchus interface ExplicitRungeKuttaIntegrator, method applyInternalButcherWeights, when calling method equations.computeDerivatives.
The underlying error is that the atmospheric density from NRL MSIS 00 model is NaN.

The error does not disappear when using a fixed-step Runge-Kutta (RK) integrator.
Getting to it with a RK integrator:

  • Put a conditional breakpoint on l.131 of ExplicitRungeKuttaIntegrator (condition: k == 3 && t0 >=0.011)
  • Activate breakpoint at l.77 of DragForce at density computation
  • If you step you will see that rho = NaN
  • If you don’t and put a break point in NRLMSISE00.Output.gts7 at L.1745
  • The error happens when computing the density of ATOMIC_NITROGEN, function densu at L.1745 returns infinity, and from there the NaNs appear.

Looking a bit further, in densu, L.2612:

double expl  = (tt <= 0) ? 50. : FastMath.min(50., FastMath.exp(-s2 * gamma * zg2));
double densu = dlb * FastMath.pow(tlb / tt, 1.0 + alpha + gamma) * expl;

Here expl is 5.26e-299 and FastMath.pow(tlb / tt, 1.0 + alpha + gamma) is 2.16e298.
Multiplying the two of them should give around 1.136 but, for some reason, Java evaluates them to +Infinity… :smiling_face_with_tear:

Could you please open an issue on the forge mentioning the problem, this forum thread, and your example in Python and mine in Java?

Thanks!
Maxime

5

It appears that just splitting the computation in 2 steps will correct the bug.

double expl  = (tt <= 0) ? 50. : FastMath.min(50., FastMath.exp(-s2 * gamma * zg2));
double xTest = FastMath.pow(tlb / tt, 1.0 + alpha + gamma);
double densu = dlb * xTest * expl;

I don’t know why. Does anyone have an idea?
With this @Simo-Dalbo the test crashes after around 2 months and 10 days because the satellite crashed on the ground.

1 Like
6

Hi @MaximeJ,
Thank you very much for your help! I opened an issue on the forge, as you suggested.
You can find it here: NaN appears during the integration (#1365) · Issues · Orekit / Orekit · GitLab.

Have a nice day,
Simone

2 Likes
7

Thanks a lot for investigating Maxime. Simone also reported that the crash happened with another atmospheric model. Were you also able to reproduce that?

Cheers,
Romain.

8

Ah… No, I missed that one :pensive:
I’ll have a look at it when I have time. I guess it must happen for similar reasons.

9

Hi again @Simo-Dalbo, @Serrof,

The issue with DTM2000 is in the same idea and happens here (L. 641).
It is fixed by replacing:

final double fzI = FastMath.pow(t120tz, upapg) * FastMath.exp(-sigzeta * gamma);

With

final double fzI = FastMath.exp(FastMath.log(t120tz) * upapg - sigzeta * gamma);

I’ll add it to the correction of the bug.

Cheers,
Maxime

2 Likes
10

Hi,
I’ve proposed a merge request to fix this issue.
It will be available in patch 12.0.3 or version 12.1, depending on whether the next release is a patch or a minor version.
@Serrof I’ve chosen you as a reviewer for this one :wink:

11

Thank you very much again for the help, glad it has been solved! :smiley:

Greetings,
Simone