Add Implicit Trapezoidal method

[Steven-Roberts]

Description

Adds a fixed step implicit trapezoidal method for comparing against adaptive and partitioned integrators. It uses Newton's method with line search for solving nonlinear systems which was found to be key at large time steps.

Proposed changes

• Adds new VectorHandler operations for norms
• Removes internal state management from InfNorm which had an initialization bug and instead use new VectorHandler operations
• Adds ImplicitTrapezoidal class

Checklist

• All tests pass.
• Code compiles cleanly with flags -Wall -Wpedantic -Wconversion -Wextra.
• The new code follows GridKit™ style guidelines.
• There are unit tests for the new code.
• The new code is documented.
• The feature branch is rebased with respect to the target branch.
• The CHANGELOG.md has been updated to reflect the changes. If this is a minor PR that is part of a larger fix already included in the file, state so.

Further comments

Currently, this is not used outside tests. In a separate PR, we'll need to add support for swapping between IDA, Rosenbrock, and this for phasor dynamics application and other "drivers".

[lukelowry]

Will be invaluable for validating EMT. Important!

[superwhiskers]

instead of this scope-guard looking thing, why can't we use something like the experimental scope guards? there are similar things implemented in third-party libraries which could also be used.

[superwhiskers]

this seems indicative of an argument for more use of exceptions to avoid this sort of code being used all over the place. not very relevant to this pr, though, outside of wanting to remove these file-local error handling helpers.

[Steven-Roberts]

I followed the convention in Rosenbrock.cpp, but I agree exceptions would be more idiomatic than int return values and macros like this.

[pelesh]

@alexander-novo, could we think of a better way to handle errors?

@Steven-Roberts, how are similar instances handled in SUNDIALS?

[Steven-Roberts]

@Steven-Roberts, how are similar instances handled in SUNDIALS?

With a similar macro or just writing out the if statement to check the return value. But SUNDIALS does not have access to exceptions, std::expected, and other C++ error handling mechanisms.

[alexander-novo]

@alexander-novo, could we think of a better way to handle errors?

The returning an int for an error code is a very C-like way of handling errors, and using these kind of quality-of-life macros isn't so uncommon in C I think.

A more C++-like approach might be more natural. std::expected or std::variant as suggested by @superwhiskers work very well in Rust, but won't help us with this specific boilerplate issue. Throwing exceptions will reduce the boilerplate heavily, but we'll need to make sure to heavily document which kinds of exceptions can be thrown by functions.

[pelesh]

I would suggest separating the interface to the ImplicitTrapezoidal integrator from the integrator itself. This is on to-do list for the Rosenbrock, as well.

For example Ida class is just an interface to IDA solver. Imho ImplicitTrapezoidal interface should follow the same API and the integrator implementation should be within a directory GridKit/Solver/Dynamics/ImplicitTrapezoidal.

[pelesh]

Is this really the intent here to throw an exception upon failed step? Why not simply clean-up and return with error message? I think exception should be reserved for events that should not happen. For a fixed step solver it is not unexpected that for certain step selection the computation would not converge.

[pelesh]

I would avoid Matlab approach of storing output times a priori. I would especially avoid using std::vector in this code.

[pelesh]

I would avoid these callbacks. We are trying to remove them altogether from the code. I don't think they should be in the core numerical integrator code.