add new methods and update testing for accuracy

[Copilot]

This pull request introduces several updates and improvements to the documentation, examples, and diagnostics for the operator-partitioned multiphysics solver (op_engine). The main changes expand the documentation to cover new stiff ODE methods, add a new diagnostic example for comparing solver splits, and update existing examples for clarity and reproducibility. The most important changes are grouped below.

Documentation and Guidance Improvements

• Expanded method documentation and usage guidance: The README.md now describes additional stiff ODE solvers (implicit-euler, trapezoidal, bdf2, ros2), provides explicit guidance on method selection, operator configuration, and adaptive stepping, and clarifies model shapes and configuration options. [1] [2] [3]

New Diagnostic Example

• Added comprehensive IMEX/implicit diagnostic: A new script, examples/biogeo_split_diagnostic.py, runs and compares various IMEX splits and fully implicit methods in op_engine against SciPy's BDF solver, reporting L2 trajectory and final-state errors as well as timings. This aids in benchmarking and validating solver accuracy and stability.

Example and API Updates

• Updated biogeochemical network example: The main example (examples/biogeochemical_network.py) now covers all available solver methods, clarifies documentation, and recommends preferred IMEX splits for the model. Output descriptions and recommendations have been improved for user clarity.

• Improved reproducibility and parameter defaults: Tolerances for both SciPy and op_engine are tightened for more accurate comparisons, and model parameters are adjusted for consistency (e.g., disabling seasonality by default, changing n_t and sea). [1] [2] [3]

Codebase Simplification

• Removed unused report/manifest utilities: Unused runtime report, sweep plotting, and manifest-writing code has been removed from the biogeochemical example, simplifying the codebase and focusing on canonical runs and plotting. [1] [2]

• Minor usability and logging improvements: Added a simple logging helper and more informative progress messages in the canonical run orchestration. [1] [2] [3]

[MacdonaldJoshuaCaleb]

check CI and fix failing CI, note there have been updates to CI on mian

[Copilot]

check CI and fix failing CI, note there have been updates to CI on mian

I checked CI and fixed the failing job caused by a Ruff lint error (PLR0913) in examples/biogeo_split_diagnostic.py by adding a targeted inline ignore on the parameter-heavy helper. Fix is in commit f162a54.

Screenshot (CI status via MCP/API): latest CI run for this commit is in progress with no failed jobs currently reported.

[TimothyWillard]

I think perhaps the examples should be extracted out to a separate PR, because there's a lot of work to do on that front in terms of documenting and thinking about usage. Such as:

• How are the examples expected to be used? E.g. as reference point, as a starting point, an explanation, a demo, etc.?
• How are users expected to run examples? It seems like it is structured so a user clones the repository and then does uv run python examples/biogeo_split_diagnostic.py but that is not noted/advertised anywhere.
• I think the lack of documentation, especially inline comments explaining what's going on make it hard to approach for new users.
• Also splitting between multiple files is less than ideal.

I think it might be better to focus on the documentation site and adding guides there rather than the code dumps in the examples. Can think about ways to incorporate examples into the documentation as well when focusing on documentation (Demo notebooks like vaxflux or documentation assets like flepimop2).

[TimothyWillard]

Maybe this would be good as a standalone guide at the documentation site: https://accidda.github.io/op_engine/latest/? The current descriptions are informative, but perhaps brief. I can see an argument for this being a follow up issue.

[TimothyWillard]

I think perhaps blanket ignore PLR0914 (https://docs.astral.sh/ruff/rules/too-many-locals/), especially for tests, is less than ideal. I would be curious to know what tests are using more than 15 local vars.

[TimothyWillard]

Why not use from op_engine.core_solver import MethodName?

[TimothyWillard]

Minor: I don't think these test file specific helpers like _exact_* need to be private, these are just test files so no need to worry about curating the API.

[TimothyWillard]

If this constant is only used once by the @pytest.mark.parametrize decorator I'm not sure I see the value in declaring it a constant?

[TimothyWillard]

This seems a bit odd to put this here in the CoreSolver class when this information is available to RunConfig/RunPlan and could be validated there using the __post_init__ hook for dataclasses.

[TimothyWillard]

Similar to another comment this seems odd to me to have this on the CoreSolver class when it should be done at the RunConfig/RunPlan level?

[TimothyWillard]

This seems like a pretty big code smell to me, especially C901.

[TimothyWillard]

Suggested change

	if plan.method == "implicit-euler":
	if plan.jacobian is None:
	raise RuntimeError(_INTERNAL_ERROR_OP_AXIS_MSG)
	return self._step_implicit_euler_linearized(
	rhs_func,
	step=step,
	jacobian=plan.jacobian,
	)
	if plan.method == "trapezoidal":
	if plan.jacobian is None:
	raise RuntimeError(_INTERNAL_ERROR_OP_AXIS_MSG)
	return self._trapezoidal_linearized_step(
	rhs_func,
	step=step,
	jacobian=plan.jacobian,
	)
	if plan.method == "bdf2":
	if plan.jacobian is None:
	raise RuntimeError(_INTERNAL_ERROR_OP_AXIS_MSG)
	return self._step_bdf2_linearized(
	rhs_func,
	step=step,
	jacobian=plan.jacobian,
	)
	if plan.method == "ros2":
	if plan.jacobian is None:
	raise RuntimeError(_INTERNAL_ERROR_OP_AXIS_MSG)
	return self._step_rosenbrock_w2(
	rhs_func,
	step=step,
	jacobian=plan.jacobian,
	)
	if (step_func := _PLAN_TO_STEP_FUNC.get(plan.method)) is not None:
	if plan.jacobian is None
	raise RuntimeError(_INTERNAL_ERROR_OP_AXIS_MSG)
	return step_func(rhs_func, step=step, jacobian=plan.jacobian)

And then have _PLAN_TO_STEP_FUNC be a constant, or I guess since these are bound methods maybe make it a method/classvar? But you get the idea.

[TimothyWillard]

Something that concerns me about the CoreSolver class in particular is there is the confluence of mutating state and returning state that I think makes it a nightmare to unit test the class currently. As is it is very unclear how I would write a unit test for an individual method in this class, even though these are largely just compositions of mathematical operations which should be ideal to unit test.

Perhaps outside the scope of this PR but perhaps need to consider an API refactor that focuses on end user usability and unit testability. Maybe could think about ways to group methods so they implement a parent class or have a scipy style solve_ivp/optimize functional interface. I'm not really sure what's best or pros/cons of both ahead of time.