|
FhSim
3.1.0
Marine systems simulation
|
|
|
FhSim
3.1.0
Marine systems simulation
|
The <Integrator> block accepts three child elements that together control how the nonlinear and linear algebra is handled inside SUNDIALS at each time step:
| Parameter | XML element / attribute | Default | Description |
|---|---|---|---|
| Jacobian type | <Jacobian Type="..."> | dense | Storage format of the Jacobian matrix passed to SUNDIALS |
| Linear solver | <LinearSolver Type="..."> | SPGMR | Linear solver used within Newton iterations |
| Nonlinear solver | <LinearSolver NonlinearSolver="..."> | newton | Nonlinear iteration strategy |
<Jacobian> and its attributes apply to the SUNDIALS backend only. On the Engine backend <Jacobian Type> is ignored (the storage follows from the linear solver) and <Jacobian Policy> is ignored with a warning in the log.| Value | Matrix | Best for |
|---|---|---|
auto | Exactly the same as dense | — see the note below |
dense | Full N×N dense matrix | Small to medium systems (up to ~200 states) |
band | Banded matrix with upper/lower bandwidth | Chain-coupled or banded systems where bandwidth << N |
sparse | Sparse CSR matrix (N×N allocated, filled analytically or by differences) | Large sparse systems; enables diagonal Jacobi preconditioner for SPGMR |
none | No matrix allocated | Matrix-free Krylov (SPGMR); SUNDIALS uses internal difference quotients |
Type="auto" decides nothing. It is a parse-time spelling of dense: the two are indistinguishable after the input file is read, and no model capability is inspected. The attribute that does inspect the model is Policy — see Implicit Integration with Analytical Jacobians. The linear solver has its own, genuinely adaptive auto; see Linear Solver Selection.When <Jacobian Type="band">, the bandwidth is auto-detected from the Jacobian sparsity pattern (the maximum |i-j| over all non-zero entries). You can override the detected values with BandUpper and BandLower on <LinearSolver>:
If only BandUpper is specified, BandLower defaults to the same value (symmetric bandwidth).
| Value | Solver | Best for |
|---|---|---|
SPGMR | Krylov iterative solver (Scaled Preconditioned GMRES) | Large or sparse systems; also the only option when <Jacobian Type="none"> |
DENSE | Direct dense LU factorisation | Small systems where the full Jacobian is available |
BAND | Direct banded LU factorisation | Systems with limited bandwidth (bw << N); requires <Jacobian Type="band"> |
| Value | Strategy |
|---|---|
newton | Newton iteration (default, recommended for stiff systems) |
fixed_point | Anderson-accelerated fixed-point iteration. Expert tuning via <SundialsExtra>. |
Not all combinations are meaningful. The table below lists every combination and whether it is valid:
| Jacobian Type | LinearSolver Type | Valid | Notes |
|---|---|---|---|
dense | DENSE | Yes | Direct dense LU solve. Recommended for small systems (< ~50 states). |
dense | SPGMR | Yes | Iterative Krylov solve with a dense Jacobian available. Useful when the system is moderately large but a full Jacobian is still affordable. |
band | BAND | Yes | Direct banded LU solve. Best when bandwidth << N (e.g., chain-coupled systems). Storage O(N·bw), factorization O(N·bw²). |
band | DENSE | No | The dense direct solver cannot operate on a band matrix object. Use <LinearSolver Type="BAND">. |
band | SPGMR | No | Use <Jacobian Type="none"> with <LinearSolver Type="SPGMR" Preconditioner="block_tridiagonal"/> for iterative + block preconditioner instead. |
sparse | SPGMR | Yes | Sparse Jacobian is used to build an ILUT (Incomplete LU with Threshold) preconditioner for GMRES. Best choice for large sparse systems (> 200 states). |
sparse | DENSE | No | The dense direct solver cannot operate on a sparse matrix object. FhSim reports an error at startup. |
none | SPGMR | Yes | Matrix-free: SUNDIALS computes Jacobian-vector products by finite differences internally. No Jacobian storage. |
none | DENSE | No | The dense direct solver requires an explicit Jacobian matrix. FhSim reports an error at startup. |
<Jacobian Type> setting still controls the storage format that SUNDIALS uses for that Jacobian. Use <Jacobian Type="dense"> for full dense storage or <Jacobian Type="sparse"> to pass it as a sparse CSR matrix (required for the ILUT preconditioner with SPGMR).Small system, direct solver (recommended default):
Large sparse system, iterative solver with ILUT preconditioner:
Matrix-free solve (no Jacobian storage):
Banded direct solve (chain-coupled or limited bandwidth systems):
In FhSim v3 the integrator no longer controls output timing. The flat <Integrator> attributes for output scheduling used in earlier versions have been removed. Output timing is now configured through the observer/output pipeline — see OBSERVERS section. Internal step-size bounds are set on <StepControl> via StepMin, StepMax, Step0, and MaxSteps.
During Newton iteration, implicit integration requires solving the linear system
\[ (I - \Delta t \cdot J) \cdot \Delta x = r \]
where \(J\) is the system Jacobian. FhSim supports multiple linear solver strategies, selected automatically by default or manually via configuration.
| Solver keyword | Description | Complexity | Best For |
|---|---|---|---|
auto | Heuristic selection based on topology and sparsity | — | Default; works well for most cases |
dense | Dense LU factorization (Eigen::PartialPivLU) | O(n³) | Small systems (n < 30) |
sparse | Sparse LU with METIS nested dissection ordering | O(nnz·√n) typical | Medium/large sparse systems; especially effective for 2D/3D mesh structures (trawl nets, aquaculture nets) |
band | Direct banded LU with automatic RCM reordering | O(n·bw²) | Narrow-bandwidth chains (bw < n/4) |
iterative | BiCGSTAB with ILUT preconditioner | O(k·nnz)/iter | Very large sparse systems (n ≥ 5000) |
iterative_gmres | GMRES with ILUT preconditioner | O(k·nnz)/iter | Ill-conditioned systems where BiCGSTAB stagnates |
block_tridiagonal | Thomas algorithm for block-tridiagonal systems | O(N·bs³) | Pure chain topologies |
near_tridiagonal | Woodbury-corrected block-Thomas | O(N·bs³ + k³) | Single-SO rings, T-junctions, psi-junctions |
component | Direct solve per connected component of the SimObject coupling graph; each component uses the structured local solver its sparsity allows (block-Thomas, Woodbury, dense LU) and is assembled without a dense N×N matrix | O(Σ n_c³) worst case, O(N·bs³) for chain components | Multi-SimObject systems; the only structure it exploits is within a component |
By default, <LinearSolver Type="auto"> applies a topology-aware hierarchy based on the cached Jacobian sparsity pattern and SimObject graph structure:
| Priority | Condition | Selected Solver |
|---|---|---|
| 1 | n < 30 | DenseLU |
| 2 | Pure block-tridiagonal chain detected | BlockTridiagonal |
| 3 | Single SimObject, density < 0.5, near-tridiagonal structure confirmed | NearTridiagonal |
| 4 | Single SimObject, density < 0.5, near-tridiagonal not viable (e.g. 2D mesh) | falls to 5–7 |
| 5 | ≥ 2 stateful SimObjects, density < 0.5 | Component |
| 6 | density < 0.25 AND bandwidth < n/4 | BandedLU |
| 7 | density < 0.25 AND n ≥ 5000 | Iterative (BiCGSTAB) |
| 8 | density < 0.25 | SparseLU (METIS ordering) |
| 9 | Otherwise (dense or irregular) | DenseLU |
The near-tridiagonal viability check at priority 3 calls NearTridiagonalSolver::DetectStructure at selection time. This means 2D/3D mesh SimObjects (trawl nets, aquaculture nets, flexible surfaces) whose sparsity graph has too many off-tridiagonal connections correctly fall through to SparseLU (priority 8) rather than stalling at a failed NearTridiagonal setup and degrading to DenseLU.
IterativeGMRES and explicit NearTridiagonal / Component overrides are also accepted as manual selections. The sparsity pattern and topology are computed once at model initialization with negligible runtime overhead.
SparseLU or Iterative may now be routed to BlockTridiagonal, NearTridiagonal, or Component for better performance.To override automatic selection, use <LinearSolver> child element:
<LinearSolver> attributes MaxIter, Tol, and Restart configure iterative solver parameters. Jacobian verification and sparsity logging are configured in the unified <Diagnostics> block.Example with Jacobian verification enabled:
For banded systems (e.g., a chain of 6-DOF bodies), <LinearSolver Type="band"> is often the fastest direct solver. Bandwidth is auto-detected from the sparsity pattern; user overrides are available:
If BandUpper and BandLower are omitted, they default to the symmetric bandwidth detected from GetJacobianSparsity (or n-1 if no sparsity info is available).
near_tridiagonal (single SO, auto-detected) or sparse.sparse solver with METIS nested dissection ordering is the recommended direct solver. For very large meshes (n > 2000) consider iterative or iterative_gmres.dense.iterative scales better than direct sparse LU. Tune convergence with <LinearSolver MaxIter="..."> and <LinearSolver Tol="..."> attributes.Iterative fails to converge on ill-conditioned systems, switch to iterative_gmres. GMRES has more robust convergence (monotonic residual decrease) at the cost of higher memory (stores a restart-count of Krylov vectors). Adjust the restart count via <LinearSolver Restart="30"> — higher values improve convergence but use more memory.band (banded LU) is often the fastest direct Engine/BackwardEuler_i solver — O(n·bw²) factorize vs O(n·bw³) for block-tridiagonal and O(nnz·√n) for sparse LU. The auto heuristic selects it when bw < n/4.component for a block Gauss-Seidel preconditioned GMRES that exploits the coupling graph topology.Several internal optimizations reduce solve time and memory cost regardless of which linear solver is selected.
The band solver applies Reverse Cuthill-McKee (RCM) reordering before factorization. RCM permutes the rows and columns to minimize the matrix bandwidth, which reduces the fill-in during LU factorization and lowers the O(n·bw²) cost:
The permutation perm[] and its inverse iperm[] are computed once at Setup time and reused throughout the simulation. No user configuration is required; the banded solver always reorders.
BlockTridiagonalSolver (also when it is the local solver of a ComponentSolver component) detects the block size bs at Setup time and dispatches to a compile-time fixed-size kernel for bs = 1–16. Each fixed-size kernel stores blocks as Eigen::Matrix<double, bs, bs> (stack-allocated, fully unrolled) and uses Eigen::PartialPivLU<Matrix<double,bs,bs>> which the compiler can vectorize and inline completely.
| Block size | Kernel storage | LU solver |
|---|---|---|
| 1–16 | unique_ptr<BlockTriKernel<bs>> via macro dispatch | PartialPivLU<Matrix<double,bs,bs>> (compile-time) |
| ≥ 17 | Heap-allocated MatrixXd blocks | PartialPivLU<MatrixXd> (runtime size) |
The dispatch uses a switch over bs with ALL_BLOCK_SIZES(MACRO) expanding 16 cases. Heap-allocated kernels are stored as unique_ptr<void, KernelDeleter> with a typed function-pointer deleter so that the correct destructor is called without exposing the templated type in the header.
For the common 13-state rigid-body case (6-DOF body in FhSim), the bs=13 kernel is approximately 3–5× faster than the dynamic-size path on a typical x86-64 system.
Several solvers own their Jacobian representation directly and do not require FhSim to assemble a full N×N dense matrix first:
| Solver | NeedsDenseJacobian() | How Jacobian is obtained |
|---|---|---|
DenseLUSolver | true | Receives dense N×N matrix from caller |
BlockTridiagonalAdapter | false | Assembles block-tridiagonal pattern directly from SimObject OdeJacobian() and port links |
ComponentSolver | false | Assembles each connected component from the block Jacobians and the couplings inside it |
SparseLUAdapter | false | Fills CSR values via AssembleSparseJacobianValues() |
BandedLUAdapter | true | Receives dense matrix; extracts band after RCM reordering |
When NeedsDenseJacobian() returns false, the implicit method skips the O(n²) dense assembly step entirely and calls FactorizeFromModel() instead, which lets each solver fill only the entries it needs. For large systems this can reduce the per-Newton-step wall time by 30–70%.
BiCGSTAB (iterative) | GMRES (iterative_gmres) | |
|---|---|---|
| Memory | O(n) — fixed | O(k·n) — grows with restart parameter |
| Convergence | Can stagnate on difficult spectra | Monotonic residual decrease, more robust |
| Tuning | None beyond MaxIter/Tol | <LinearSolver Restart="..."> (default 30) |
| Recommendation | Default for large systems | Fallback when BiCGSTAB fails to converge |
For systems with chain topology — such as a series of rigid bodies connected end-to-end — the Jacobian has block-tridiagonal structure. The block_tridiagonal solver exploits this structure using the block Thomas algorithm, achieving O(N × bs³) complexity instead of O((N·bs)³) for dense LU, where N is the number of blocks and bs is the block size.
This is particularly effective for:
With the <Integrator> integrator, block_tridiagonal is used as a direct linear solver during Newton iteration:
With SUNDIALS (CVODE or ARKode), block_tridiagonal configures SPGMR with a block-Thomas preconditioner. The preconditioner assembles the block-tridiagonal portion of the Jacobian directly from SimObject Jacobians and port connections (without forming the full dense Jacobian), then factorizes P = I - gamma·J using the Thomas algorithm. Use <Jacobian Type="none"> for matrix-free Jacobian-vector products:
The preconditioner supports SUNDIALS jok (Jacobian OK) optimization: when only gamma changes between Newton iterations, the raw Jacobian blocks are re-factorized without re-assembling them from the SimObjects.
Detection: When block_tridiagonal is selected, FhSim automatically:
PortJacobianLink topology for chain structure (topological sort from endpoint)Block size: Auto-detected from SimObject state sizes (multi-object chains) or by scanning for the largest block size that divides the state count and matches the sparsity pattern (intra-object). For the fixed 13-state case, compile-time optimized Eigen types (Matrix<double,13,13>, PartialPivLU<Matrix13>) are used.
Performance comparison (N=10 bodies × 13 states = 130 total states):
| Solver | Complexity | Approximate Flops |
|---|---|---|
| Dense LU | O(n³) | ~2.2M |
| Banded LU (bw=25) | O(n·bw²) | ~81K |
| Sparse LU | O(nnz·√n) | ~50K |
| Block Thomas | O(N·bs³) | ~22K |
For single-SimObject systems that are almost block-tridiagonal — such as rings of bodies, T-junctions, or psi-junctions — the near_tridiagonal solver extends the block-Thomas algorithm using the Woodbury identity (Sherman-Morrison-Woodbury matrix identity).
The idea is to decompose the full system matrix as:
\[ A_{full} = A_{tridiag} + U V^T \]
where \(A_{tridiag}\) is the block-tridiagonal core and \(U\), \(V\) (each \(n \times k\)) encode the extra off-diagonal blocks that break tridiagonality (e.g., the closing edge of a ring). With \(k \ll n\) this adds only a cheap \(k \times k\) correction on top of the Thomas solve.
Solve algorithm (per Newton step):
Factorize cost: k additional block-Thomas solves to build \(A_{tridiag}^{-1} U\), plus one \(k \times k\) LU factorization.
Viability: The correction rank \(k\) = (number of extra blocks) × block_size must satisfy \(k < n/4\). When this condition is not met (the system is too far from tridiagonal), the structure is not detected and the auto-selector falls back to SparseLU or Component.
All of the following are detected automatically from the CSR sparsity pattern:
| Topology | Extra blocks | Correction rank |
|---|---|---|
| Ring (N bodies, closed loop) | 2 (blocks (0,N-1) and (N-1,0)) | 2·bs |
| T-junction (main chain + 1 branch) | 2 (branch connection, symmetric) | 2·bs |
| Psi-junction (main chain + 2 branches) | 4 | 4·bs |
| Any combination with k < n/4 | ≤ n/4 / bs extra block pairs | < n/4 |
Engine/BackwardEuler_i (direct solver):
SUNDIALS (as preconditioner for SPGMR):
When <LinearSolver Type="auto"> is in effect and there is a single stateful SimObject with a near-tridiagonal sparsity pattern and density < 0.5, near_tridiagonal is selected automatically.
The dominant cost is the block-Thomas solve, repeated k+1 times at factorize and once at solve time. For a ring of N bodies (k = 2·bs), the factorize cost is roughly 3× that of a pure chain — still far cheaper than SparseLU or dense LU for large N:
| Solver | Ring N=20 bodies (bs=6, n=120) |
|---|---|
| Dense LU | O(n³) ≈ 1.7M flops |
| SparseLU | O(nnz·√n) ≈ 150K flops |
| Near-tridiagonal | O((k+1)·N·bs³) ≈ 25K flops |
The sparse solver (<LinearSolver Type="sparse">) uses Eigen's SparseLU with METIS nested dissection fill-reducing ordering. METIS is the industry standard for sparse direct solvers on mesh-like sparsity patterns.
Why METIS over COLAMD:
COLAMD (Column Approximate Minimum Degree) is a greedy heuristic that minimises fill-in locally. METIS nested dissection finds global separator sets that partition the graph recursively, producing a much smaller fill footprint for 2D and 3D mesh structures. The asymptotic advantage is well-established:
| Problem type | COLAMD fill | METIS fill |
|---|---|---|
| 1D chain (banded) | O(n·bw) | O(n·bw) — same |
| 2D mesh (M×N grid) | O(n^1.5) | O(n log n) — significantly less |
| 3D mesh | O(n^2) | O(n^1.33) |
For a trawl net or aquaculture net modelled as a single SimObject with a triangulated 2D surface (extruded ring or planar mesh):
The ordering is computed once at AnalyzePattern time (model initialisation); subsequent Factorize calls reuse the permutation at no extra cost.
Usage (the sparse solver is auto-selected when density < 0.25 and no structured topology is detected; it can also be forced manually):
Both the banded solver (<Jacobian Type="band"> + <LinearSolver Type="BAND"> for SUNDIALS, or <LinearSolver Type="band"> for Engine/BackwardEuler_i) and the block-tridiagonal solver (<LinearSolver Type="block_tridiagonal">) target systems with limited coupling range, but they serve different use cases:
| Banded LU | Block-Tridiagonal | |
|---|---|---|
| Structure required | Any banded Jacobian | Strict block-tridiagonal (chain topology) |
| Non-uniform blocks | Handled naturally | Requires uniform or compatible block sizes |
| SUNDIALS role | Direct solver (replaces dense LU) | Preconditioner for SPGMR |
| Jacobian assembly | SUNDIALS handles it (analytical or FD) | Custom assembly from SimObject Jacobians |
| Best for | Banded but not cleanly block-tridiagonal | Clean chains of identical bodies |
| Asymptotic cost | O(N·bw²) where bw = 2·bs-1 | O(N·bs³) — ~4× fewer ops for block-tridiagonal |
Rule of thumb:
<LinearSolver Type="block_tridiagonal">.<Jacobian Type="band"> + <LinearSolver Type="BAND">.<LinearSolver Type="component"> (see below).The component solver (<LinearSolver Type="component">) is a direct solver over the connected components of the SimObject coupling graph. The stateful SimObjects are the nodes and the port-Jacobian links the edges. Two SimObjects that are not connected — directly or through other SimObjects — never share a row of \(I - \gamma J\), so each connected component is an independent linear system:
\[ \begin{bmatrix} A_1 & & \\ & A_2 & \\ & & \ddots \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \end{bmatrix}, \qquad A_k = I - \gamma J_k . \]
At Setup the solver partitions the SimObjects into components and picks, per component, the local solver its sparsity allows:
| Component structure | Local solver |
|---|---|
| A correction-free block partition exists | BlockTridiagonalSolver (block-Thomas, fixed-size kernels for bs ≤ 16) |
| One exists only with a small Woodbury correction | NearTridiagonalSolver (block-Thomas + Woodbury) |
| Anything else, a single SimObject, or ≤ 64 states | Eigen::PartialPivLU (dense LU) |
The partition is detected over the component's states in the order the solver discovers its SimObjects (a breadth-first walk of the coupling graph), which is itself bandwidth-reducing. A ring, T-junction or psi-junction chain is therefore usually relabelled into a band and served by BlockTridiagonalSolver at a coarser block size, not by the Woodbury path.
Each component's \(J_k\) is assembled straight from the model — the SimObjects' OdeJacobian() blocks plus the port couplings inside the component — so no dense N×N Jacobian is ever formed (NeedsDenseJacobian() is false). Re-factorizing at a new \(\gamma\) reuses the assembled blocks.
What the solver does not do is decompose a connected system. A model whose stateful SimObjects are all coupled together is one component, and component is then one structured (or dense) direct solve of the whole system. Its saving over dense is the skipped dense assembly and the structured local solver, not a smaller factorization. If your system is one big component with no chain structure, sparse or band is usually the better choice.
The solve is exact — the local solver is picked from the sparsity of everything the component assembles, couplings included, so the partition it settles on covers every assembled entry — so Newton converges quadratically and, as a SUNDIALS SPGMR preconditioner, component is a perfect preconditioner: SPGMR converges in one iteration.
History. Earlier versions described a "Schur complement" mode with an interface between components and a block Gauss-Seidel + GMRES fallback tuned by
GsSweeps/GsSymmetric. Because components are, by construction, not coupled to each other, that interface was always empty and the fallback was never reached; the code was removed and the two attributes are now rejected with an error.
Engine/BackwardEuler_i:
SUNDIALS CVODE/ARKode with SPGMR:
| Scenario | Recommended solver |
|---|---|
| Pure chain of identical bodies | block_tridiagonal |
| Ring, T-junction, psi-junction (single SO) | near_tridiagonal (auto-detected) |
| Several independent (uncoupled) groups of SimObjects | component — one factorization per group |
| Multi-SO sparse system (auto) | component (selected automatically) |
| One large fully-coupled multi-SO system without chain structure | sparse or band |
| Single large monolithic SimObject, chain/ring/junction structure | near_tridiagonal |
| Single large monolithic SimObject, 2D/3D mesh structure (e.g. trawl net) | sparse (auto-selected via METIS SparseLU) |
| Single large monolithic SimObject, very large mesh (n ≥ 5000) | iterative or iterative_gmres |
When using CVODES or ARKODE/DIRK with <LinearSolver Type="SPGMR">, the Preconditioner attribute of <LinearSolver> selects the SUNDIALS preconditioner:
<LinearSolver Preconditioner> | SPGMR Preconditioner | Description |
|---|---|---|
block_tridiagonal | Block-Thomas | Assembles and factorizes block-tridiagonal portion; O(N·bs³) per setup |
near_tridiagonal | Near-Thomas + Woodbury | Block-Thomas + Woodbury correction for rings/junctions |
component | Component direct solve | Exact solve per connected component — a perfect preconditioner |
jacobi | Jacobi diagonal | P = diag(I − γJ); the Jacobian diagonal is read from the model, so it works with every <Jacobian Type> including none |
none | None | SPGMR runs unpreconditioned |
iterative / iterative_gmres | Jacobi diagonal | These name iterative Newton solvers; as an SPGMR preconditioner they fall back to the Jacobi diagonal |
(default, auto) | Structural, else Jacobi diagonal | The structural preconditioner of the model's resolved solver type when it has one, otherwise the Jacobi diagonal |
The Jacobi diagonal preconditioner needs the Jacobian diagonal. Entries where 1 − γ·J_ii comes out zero, negative or non-finite are left unpreconditioned (P_ii = 1) and reported once as a warning. Preconditioner="jacobi" on a model where no SimObject supplies a Jacobian is an error; with auto the same situation only logs a warning and SPGMR runs unpreconditioned.
Preconditioner only applies to <LinearSolver Type="SPGMR">. Combining it with Type="BAND" or Type="DENSE" is an error — those are direct solvers and use no preconditioning.
When using ARKODE with Method="DIRK", FhSim uses a default DIRK method (typically 4th order). To select a specific DIRK Butcher table, use the ARKodeSetTableNum option:
| Method Name | Stages | Order | Embedding | Notes |
|---|---|---|---|---|
BACKWARD_EULER_1_1 | 1 | 1 | – | Implicit Euler, A-stable, L-stable |
IMPLICIT_MIDPOINT_1_2 | 1 | 2 | – | Implicit midpoint rule |
IMPLICIT_TRAPEZOIDAL_2_2 | 2 | 2 | – | Trapezoidal rule (Crank-Nicolson) |
SDIRK_2_1_2 | 2 | 2 | 1 | Default 2nd order, A-stable |
ARK2_DIRK_3_1_2 | 3 | 2 | 1 | Implicit part of ARK2 pair |
BILLINGTON_3_3_2 | 3 | 2 | – | |
TRBDF2_3_3_2 | 3 | 2 | – | TR-BDF2 composite method |
KVAERNO_4_2_3 | 4 | 3 | 2 | L-stable |
ARK324L2SA_DIRK_4_2_3 | 4 | 3 | 2 | Implicit part of ARK3 pair |
ESDIRK324L2SA_4_2_3 | 4 | 3 | 2 | ESDIRK variant |
ESDIRK325L2SA_5_2_3 | 5 | 3 | 2 | ESDIRK variant |
ESDIRK32I5L2SA_5_2_3 | 5 | 3 | 2 | ESDIRK variant |
CASH_5_2_4 | 5 | 4 | 2 | Cash 4th order |
CASH_5_3_4 | 5 | 4 | 3 | Cash 4th order |
SDIRK_5_3_4 | 5 | 4 | 3 | Standard SDIRK 4th order |
KVAERNO_5_3_4 | 5 | 4 | 3 | L-stable |
ARK436L2SA_DIRK_6_3_4 | 6 | 4 | 3 | Implicit part of ARK4 pair |
ARK437L2SA_DIRK_7_3_4 | 7 | 4 | 3 | Default 4th order additive |
ESDIRK436L2SA_6_3_4 | 6 | 4 | 3 | ESDIRK variant |
ESDIRK43I6L2SA_6_3_4 | 6 | 4 | 3 | ESDIRK variant |
QESDIRK436L2SA_6_3_4 | 6 | 4 | 3 | Quasi-ESDIRK variant |
ESDIRK437L2SA_7_3_4 | 7 | 4 | 3 | ESDIRK variant |
KVAERNO_7_4_5 | 7 | 5 | 4 | L-stable |
ARK548L2SA_DIRK_8_4_5 | 8 | 5 | 4 | Implicit part of ARK5 pair |
ARK548L2SAb_DIRK_8_4_5 | 8 | 5 | 4 | Alternate ARK5 pair |
ESDIRK547L2SA_7_4_5 | 7 | 5 | 4 | ESDIRK variant |
ESDIRK547L2SA2_7_4_5 | 7 | 5 | 4 | ESDIRK variant |
Notes:
--) require fixed-step integration.For maximum stability with a simple 1st-order method:
BACKWARD_EULER_1_1 has no embedding, so adaptive stepping is disabled. Use ARKodeSetFixedStep in <SundialsExtra> to set the step size.For CVODES with Method="BDF", the solver automatically adapts between orders 1-5 based on error estimates. To limit the maximum order (e.g., for stability reasons or to reduce memory), use CVodeSetMaxOrd:
| Max Order | Behavior |
|---|---|
| 1 | Equivalent to backward Euler (most stable, least accurate) |
| 2 | BDF-2, good balance of stability and accuracy |
| 3-5 | Higher accuracy, but may have stability issues for very stiff systems |
The default is order 5, which provides the best accuracy for smooth solutions. Reducing the order can help when:
1.9.8