Repository navigation
Radau5ODE / RodasODE: use the PyFMI Jacobian by default, and a 'thet' default to match - #422
Draft
hubertus65 wants to merge 2 commits into
Draft
hubertus65 wants to merge 2 commits into
hubertus65 wants to merge 2 commits into
Conversation
AssimuloFMIAlg turns `with_jacobian` on by default only when the FMU provides directional derivatives, or (for FMUs without them, with at least PYFMI_JACOBIAN_LIMIT states) when the solver is CVode. Radau5ODE and RodasODE were excluded and therefore fell back to their own dense finite-difference Jacobian: nx rhs evaluations per Jacobian, and RodasODE needs a Jacobian every step. FMIODE2's Jacobian uses the dependency structure with CPR colouring, i.e. a handful of rhs evaluations per Jacobian, and both solvers accept it (dense, or CSC for Radau5ODE's SPARSE linear solver). Measured on Modelica.Fluid.Examples.BranchingDynamicPipes (82 states, no directional derivatives, rtol 1e-4): Radau5ODE 5.5 s -> 3.3 s, and RodasODE completes where it previously failed (the failing FMU evaluation happened at one of rodas.f's own finite-difference perturbation points). The automatic SPARSE selection for large, sparse systems stays limited to the solvers that have a `linear_solver` option (CVode, Radau5ODE); RodasODE would reject it. Users can still pass with_jacobian=False. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…e PyFMI Jacobian
Radau5 recomputes the Jacobian after every accepted step in which Newton's
contraction rate exceeded THET (Assimulo option 'thet', default 1e-3), which
on stiff FMU models is nearly every step: 1493 Jacobian evaluations for 1629
steps on a 141-state vehicle model, 1473 for 1783 steps on a 226-state plant.
CVode on the same models needs 124 and 454. With the PyFMI Jacobian each
evaluation costs 100-240 rhs calls, so this dominates the run.
Hairer's documentation recommends raising THET to about 0.1 when Jacobian
evaluations are costly. Measured at rtol 1e-6 with correct event counts and
unchanged accuracy (deviation vs an rtol 1e-8 reference 1e-5 .. 4e-4):
model thet 1e-3 -> 0.1 -> 0.5 -> 0.9 (Jacobians / solver time)
vehicle, 141 st. 1493 / 379 s -> 1242 / 347 s -> 669 / 226 s -> 654 / 219 s
plant, 226 st. 1473 / 444 s -> 1379 / 499 s* -> 1177 / 420 s -> 1189 / 399 s
data center, 86 272 -> 213 (6.8 -> 6.6 s); event-dominated small models unchanged
(* under parallel load; counts are the measure)
'Radau5ODE_options["thet"]' is now "Default": 0.1 when 'with_jacobian' is in
effect, otherwise Assimulo's 1e-3; an explicit value is passed through.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Two defaults for
Radau5ODEandRodasODE, from a benchmark study of PyFMI's stiff solvers on 8MSL and 4 Modelon-library ME FMUs (86–318 states) measured against rtol 1e-8 references.
8f6dd19with_jacobiandefaults to True forRadau5ODEandRodasODE(it was CVode only)a0f8bc4Radau5ODE_options["thet"], the Jacobian-reuse threshold, defaults to 0.1 when the PyFMI Jacobian is in effect (Assimulo's own default is 1e-3)Why
The Jacobian.
PYFMI_JACOBIAN_SOLVERSlisted CVode only, soRadau5ODE— an implicitRunge-Kutta method that needs a Jacobian on every Newton solve — fell back to Assimulo's internal
differences, i.e.
nx+1rhs evaluations per Jacobian with no colouring, while the FMU'sdirectional derivatives or PyFMI's coloured finite differences sat unused.
thet. Hairer'sTHETis the Newton contraction rate above which the Jacobian is recomputedafter an accepted step. Assimulo's 1e-3 recomputes after nearly every step. That is the right
default for a cheap analytic Jacobian and the wrong one for an FMU Jacobian costing 100–240 rhs
evaluations.
Accuracy and event counts are unchanged in every case.
Behaviour change, and the test that encodes it
This changes results for existing users who run
Radau5ODEorRodasODEon an ME FMU: theJacobian source changes, so step sequences change (within tolerance). That is the point of the
PR, but it is not a silent internal cleanup, so it gets its own PR rather than riding along
inside a feature.
tests/test_fmi2.py::test_with_jacobian_optionasserts the old behaviour and is updated here.test_radau5_thet_defaultis added: 0.1 with the PyFMI Jacobian, Assimulo's default without, anexplicit user value always wins.
thetis a new key inRadau5ODE_optionswith the sentinel"Default"; when the PyFMI Jacobianis not in effect the key is deleted before reaching Assimulo, so Assimulo's own default applies
and nothing changes for that path.
🤖 Generated with Claude Code
https://claude.ai/code/session_01PyrxGeWuwCmAH9LishMW6p