NaN ephemeris time passed to SPICE spkezr_c

[gaffarelj]

When I run a Mars ascent propagation with a RKF7(8) integrator with an increasingly small time step, and that my propagation needs for some reason to call the spkezr_c method from SPICE, I get the following error:

Toolkit version: N0066

SPICE(DAFNEGADDR) --

Negative value for BEGIN address: -2147053925

A traceback follows.  The name of the highest level module is first.
spkezr_c --> SPKEZR --> SPKEZ --> SPKGEO --> SPKPVN --> SPKR02 --> DAFGDA

Oh, by the way:  The SPICELIB error handling actions are USER-TAILORABLE.  You
can choose whether the Toolkit aborts or continues when errors occur, which
error messages to output, and where to send the output.  Please read the ERROR
"Required Reading" file, or see the routines ERRACT, ERRDEV, and ERRPRT.

After some testing, I changed the getBodyCartesianStateAtEpoch() function from spiceInterface.cpp to print the ephemerisTime input, and it appears that, after some time into the propagagtion, the ephemeris time that is input to this function is -nan.

Now I still have no idea why tudat decides to change my ephemeris time (from 9.82325e+08) to -nan after some time. I expect that this error is caused by something that breaks in my simulation at small time steps.

Nevetherless, I suggest to add a check to the getBodyCartesianStateAtEpoch() function to trigger an error if the ephemerisTime is negative or nan. We could also print a warning and return a zero vector as the body cartesian position, but that could lead to more issues down the line.

[DominicDirkx]

The issue here is that the time-step control goes haywire, and ends up suggesting a NaN time step. I think it's a good idea to add a check, but let's put it at the level of the integrator (and possibly in the spice function as well), that way it will be caught at the 'source'

[gaffarelj]

I would then suggest the following changes:

• Raise an error in the SPICE function in case of incorrect ephemeris time.
• Print a warning and return the minimum step size for variable step integrator if the computed step size is invalid. (Edit: I would suggest taking the minimum dt if it's computed to be 0 or -inf, and the maximum one if it's computed to be inf).

About the latter: if we then use the minimum step size that is user-defined, this step size could be too high for the problem dynamics. I think that the warning that is printed should then describe this, informing that the propagation results could be slightly off due to this. We could also just use the minimum step size divided by an arbitrary factor, to try to alleviate this problem.

I am currently investigating the implementation for all of this.

[DominicDirkx]

Agreed on the first point. For the second, the issue is not that the time step goes below the permitted time step, but that the state becomes NaN/Inf at some point, which is then used in the time step control algorithm. I propose to catch this at the level of the DynamicsStateDerivativeModel class, checking at the beginning if the input is Inf/NaN.

[DominicDirkx]

Since I would propose to completely refactor the time step control anyway, I don't think it's necessary to put a lot of time into this part of the code now.

[DominicDirkx]

Fix pushed to develop :)

[gaffarelj]

Thanks for the fix :)

I am running a different fix on my machine that sets the integration step as the minimum one in RungeKuttaVariableStepSizeIntegrator.performIntegrationStep() if it deems it invalid (nan or inf). This allows me to run my fixed-step simulation with the RungeKuttaVariableStepSizeIntegrator without throwing an error. But it's probably not an optimum solution to push for everyone.

[gaffarelj]

With the exception that you have implemented to throw in the SPICE interface and DynamicsStateDerivativeModel, I get the following (instead of the error message):

python: /home/jerem/miniconda3/envs/tudat-bundle/include/eigen3/Eigen/src/Core/Block.h:146: Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>::Block(XprType&, Eigen::Index, Eigen::Index, Eigen::Index, Eigen::Index) [with XprType = const Eigen::Matrix<double, -1, -1>; int BlockRows = -1; int BlockCols = -1; bool InnerPanel = false; Eigen::Index = long int]: Assertion `startRow >= 0 && blockRows >= 0 && startRow <= xpr.rows() - blockRows && startCol >= 0 && blockCols >= 0 && startCol <= xpr.cols() - blockCols' failed.
Aborted

I think there is an issue in how integrateEquationsFromIntegrator catches exceptions and ends the termination. I have observed the above error in different cases now, and I start suspecting that, in each of these, I was actually supposed to get a more custom error message (from a throw std::runtime_error)...

[gaffarelj]

Turns out that throwing a std::invalid_argument instead of std::runtime_error allows the error to go trough, instead of the above Eigen abort message.

I changed these in DynamicsStateDerivativeModel, RungeKuttaVariableStepSizeIntegrator, and spiceInterface with 1a1ea80.

With this commit, I have also included a few more cases in which the error should be thrown in spiceInterface (negative or inf ephemeris times).

[DominicDirkx]

I think there is an issue in how integrateEquationsFromIntegrator catches exceptions and ends the termination. I have observed the above error in different cases now, and I start suspecting that, in each of these, I was actually supposed to get a more custom error message (from a throw std::runtime_error)...

That's quite odd.. Could you elaborate a bit on what you mean here? What was the output you got in this case before this update?

Turns out that throwing a std::invalid_argument instead of std::runtime_error allows the error to go trough, instead of the above Eigen abort message.

That's good. What output from this is now communicated to the user (if any)?

[gaffarelj]

That's quite odd.. Could you elaborate a bit on what you mean here? What was the output you got in this case before this update?

What I mean is that, if a std::runtime_error is thrown during the propagation itself (in this case from DynamicsStateDerivativeModel or spiceInterface), the following error is actually raised in the terminal:

python: /home/jerem/miniconda3/envs/tudat-bundle/include/eigen3/Eigen/src/Core/Block.h:146: Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>::Block(XprType&, Eigen::Index, Eigen::Index, Eigen::Index, Eigen::Index) [with XprType = const Eigen::Matrix<double, -1, -1>; int BlockRows = -1; int BlockCols = -1; bool InnerPanel = false; Eigen::Index = long int]: Assertion `startRow >= 0 && blockRows >= 0 && startRow <= xpr.rows() - blockRows && startCol >= 0 && blockCols >= 0 && startCol <= xpr.cols() - blockCols' failed.
Aborted

Now that I changed this so that we throw a std::invalid_argument instead, the proper error is raised.

That's good. What output from this is now communicated to the user (if any)?

With the latest implementation, if the cartesian state of a body at -nan is for instance requested to the spiceInterface, the user will see the following:

Error when retrieving Cartesian state from Spice, input time is -nan
Error, propagation terminated at t=982325358.987522, returning propagation data up to current time.

[DominicDirkx]

Hi Jeremie, I am making some small modifications: the isinf function is not platform-independent, so I'm removing this check from the spice functions for now.

[gaffarelj]

Perfect, I didn't know that this was platform-dependent

[DominicDirkx]

Also, the input time can be <0 (this means the time is before J2000). I'm removing this check as well.

[gaffarelj]

Oh, I thought it couldn't, because I thought the ephemeris input to spice was different to the tudat epochs. I got confused by this issue, since the error message states that a negative time is invalid. I assume that these times are then different.

My apologies if this caused any errors to someone.