User-defined root finding algorithm in {Field}ODEEventDetector

[afossa]

Prior to the implementation of #223 it was possible to define a custom root finding algorithm for each event via {Field}ODEIntegrator#addEventHandler().

Would be possible to preserve this functionality, for example by adding a getSolver() method to the new ODEEventDetector interface?

[maisonobe]

There is a redundancy problem with getThreshold in the ODEEventDetector interface and withThreshold in the AbstractODEDetector base class (and the corresponding field variants) and a potentially new withSolver method. From an event detection point of view, the threshold is related to the event itself, as various events accuracy depend on their meaning. From a solver point of view, the threshold depends on the solver. If user specifies a solver and specifies a threshold separately, how do we handle that?

[maisonobe]

Is this need for a solver related to #220? If #220 is solved without allowing user to specify a custom solver, would this suit your needs?

[afossa]

No, this is related to an old question on the forum: https://forum.orekit.org/t/extract-second-to-last-spacecraftstate-after-event-triggering/1165?u=alberto

I have an event whose g function only returns -1/+1, and is quite expensive to evaluate. At the same time, I am not interested in accurately locating the event, but just stop the integration at the step in which the function switches sign. My workaround was to pass a solver to addEventHandler() that literally does nothing, so that the refinement of the event location stops immediately (differently from what I wrote on the forum, I have dropped Orekit and I am using Hiparchus integrators directly). I am open to other options to achieve this, my only constraint is to evaluate the g() function only once per step.

Regarding the first question, the redundancy was already present before the last changes. One option would be let the user to specify either the threshold or the solver, and building a default solver if the user chooses to input the threshold.

[maisonobe]

Yes the redundancy was present and it was not managed properly. What happens if you put a maxCheckInterval larger than the largest possible step and a convergence threshold also larger than the largest possible step? I guess this could drastically reduce the number of calls to the g function, but I also guess the minimum number of evaluations will be 3, which may not suit your needs.

The settings of the numerous parameters in event detectors has been intentionally split in several calls (hence the withXxxxmethods and fluent API) because some detectors were really cumbersome with an ever increasing number of optional parameters. I remember a family of detectors that were merged together and ended up with a dozen different constructors before we switched to the current design.

Perhaps we could assume that the last call to withSolver overrides previous calls to withThreshold and also the other way round, with proper documentation so users refrain from using both, or at least are aware that if they use both settings, the first one will be overridden. This way, the last call would be the only relevant one.

[afossa]

With maxCheckInterval and threshold larger than the largest integrator step and a maxIterationCount strictly greater than 1 it misses the root.

If I set maxIterationCount to 0 or 1, I get the following exception:

Exception in thread "main" org.hipparchus.exception.MathIllegalStateException: ... failed to find root between 61,237.988 (g=-1.0E0) and 62,339.789 (g=1.0E0)
Last iteration at 61,788.889 (g=1.0E0)
    at org.hipparchus.ode.events.FieldEventState.findRoot(FieldEventState.java:384)
    at org.hipparchus.ode.events.FieldEventState.tryAdvance(FieldEventState.java:478)
    at org.hipparchus.ode.AbstractFieldIntegrator.acceptStep(AbstractFieldIntegrator.java:405)
    at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaFieldIntegrator.integrate(EmbeddedRungeKuttaFieldIntegrator.java:309)
    ...
Caused by: org.hipparchus.exception.MathIllegalStateException: maximal count (1) exceeded
    at org.hipparchus.util.Incrementor.lambda$static$0(Incrementor.java:41)
    at org.hipparchus.util.Incrementor.increment(Incrementor.java:237)
    at org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver.solveInterval(FieldBracketingNthOrderBrentSolver.java:240)
    at org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver.solveInterval(BracketedRealFieldUnivariateSolver.java:177)
    at org.hipparchus.ode.events.FieldEventState.findRoot(FieldEventState.java:373)
    ... 7 more

[afossa]

I found a workaround for this:

ODEEventDetector side (assuming forward integration):

• g() evaluates to either +1 or -1 (the expensive evaluation) if it is called at time t greater than the last two times at which it has been called. This ensures we are actually in a new integration step.
• g() returns immediately 0 if it is called at a time t such that t1 < t < t2, with t1 and t2 last two times at which it has been evaluated. This corresponds to a function call within the root refinement machinery.
• maxCheckInterval and tolerance are set to plus infinity
• maxIterationCount is set to 1

ODEStepHandler side: I am logging the output of interpolator.getPreviousState() only if this state is not interpolated, to ensure I am not storing states that correspond to intermediate times (within two integrator steps) computed by the root finding algorithm.

In this way, I evaluate g() and store the "previous state" exactly once per step, without the need of a custom root finding algorithm.

[maisonobe]

I would not rely on this workaround to work if make some unforeseen changes in EventState. I'll add the solver back.

[afossa]

Ok, thank you.

[maisonobe]

Please, tell me if the API suits you. It adds a withSolver method to event detectors (both double and field).

[afossa]

Yes, this works for me. Thanks again.