Release 4.0.0 · jacobwilliams/ddeabm

New features

Maximum step size. New optional max_step_size argument to initialize and initialize_event. It limits the magnitude of every integration step. By default, there is no limit.
User-selectable maximum order (#22). New optional max_order argument (1 to 20, default 12). With the default, results are identical to previous versions.
- A lower order can improve stability on some problems.
- A higher order mainly helps with REAL128 at very tight tolerances. For example, at a tolerance of 1e-30, max_order=16 needed 1.6 to 1.9 times fewer function evaluations than 12.
Continuing an interrupted integration. New resume_integration method. Call it after idid = -1 (step limit reached), -2 (tolerances increased by the code), -3, or -4 (problem appears stiff), and then call the integration routine again. Previously, there was no way to continue through the class interface.
interpolate has a new optional istat output, which is -1 if no integration step has been taken yet.
New public constants ddeabm_default_max_order (12) and ddeabm_max_order_limit (20).

A leftover value in the caller's idid variable (e.g., -33 or -2 from a previous call) could make a valid call fail immediately.
initial_step_mode=3 could take a first step past tout or tstop.
rtol=atol=0 (most accurate solution) always failed. It can now be continued with resume_integration.
interpolate divided by zero if called before any integration step.
A mismatched size of rtol or atol was silently ignored, and only the first element was used. A size-1 rtol combined with a size-neq atol (or the reverse) now gives vector tolerances.
Calling integrate on an uninitialized or destroyed object crashed. It now returns idid=-33.

Vector events: if any event function was zero at the start, roots of all the functions were ignored over the first interval.
gval was not set when the integration ended without finding a root.
With tstep, the root location was less accurate than the requested tolerance (it was found by extrapolating over several integration steps).
With continue=.true., the rest of the integration step in which the previous event was found was skipped (#17). Events in that part of the step could be missed, and fixed-step output points were skipped.
Vector events: the returned state was not always computed at the root time.
With continue=.true., a tmax beyond tstop was not rejected.
A tstep too small to advance the time (lost in roundoff) could produce an effectively unlimited number of output points. It now returns idid=-33.
When the input was rejected (idid=-33), the event function was still called.

These may affect existing code:

Invalid input sizes: an rtol or atol whose size is neither 1 nor neq now stops the program (error stop) in initialize. A gval of the wrong size (vector events) or an invalid tstep now returns idid=-33.
Event evaluation: the event function is evaluated at the end of each integration step, as well as every tstep (if given). Fixed-step output points are interpolated within each step.
Continuing after an event (continue=.true.): with tstep, the output points now start at the event and continue every tstep from there.
Relaxed tolerances: tolerances increased by the code (idid=-2) are kept when the integration is resumed. They are reset only at the start of a new problem (first_call).
Stiffness detection: it is not done when max_order is 4 or less.

New unit test (test/ddeabm_unit_test.f90) and regression test (test/ddeabm_regression_test.f90). The tests check results against known solutions, and they pass for REAL32, REAL64 and REAL128 builds.
A pixi environment (pixi.toml) for development, and CI updates.