|
FhSim
3.1.0
Marine systems simulation
|
|
|
FhSim
3.1.0
Marine systems simulation
|
FhSim supports implicit integration methods for simulating stiff systems that are challenging or impossible to solve efficiently with explicit methods. Implicit methods require solving nonlinear equations at each time step, which in turn requires Jacobian matrix information. FhSim supports multiple Jacobian computation strategies:
OdeJacobian() methodWhen all SimObjects with states provide analytical Jacobians, FhSim automatically uses them with implicit solvers, resulting in better performance and numerical stability.
HasJacobians() returns false, causing the system to fall back to numerical Jacobian approximation. Your existing simulations will continue to work unchanged.Implicit methods are particularly beneficial for stiff systems — systems where some states evolve much faster than others. Common examples in marine simulation include:
Signs that your system may benefit from implicit integration:
<SIMULATION> + <Integrator Method="...">.| Method | Description | Best For |
|---|---|---|
BackwardEuler_i | Implicit Backward Euler with Newton iteration | Simple stiff systems, fixed-step requirements |
StableSolver_i, despite its name, is a stability-enhanced explicit fixed-step scheme: it runs no Newton iteration and uses no linear solver, so nothing in this chapter applies to it.
Backward Euler uses the formula:
\[ x_{n+1} = x_n + \Delta t \cdot f(t_{n+1}, x_{n+1}) \]
The implicit equation is solved with Newton iteration using analytical, hybrid, or numerical Jacobian information depending on model capabilities and settings.
Use in the input file (v3):
| Method | Description | Best For |
|---|---|---|
BDF | Backward Differentiation Formula (CVODES, order 1-5) | General stiff systems, adaptive stepping |
DIRK | Diagonally implicit Runge-Kutta (ARKODE) | Stiff systems requiring high accuracy |
Use in the input file (v3):
Deriving OdeJacobian() for a non-trivial model is days of work, so decide deliberately.
Leave <Jacobian Policy> on auto unless one of these holds. auto already uses every analytical Jacobian your model provides and fills the rest numerically, so there is nothing to gain by naming a policy:
Policy="analytical" when you must guarantee no numerical entry is used — for example when a regression suite has to fail loudly if someone adds a SimObject without a Jacobian. It is an assertion, not an optimisation: auto already picks the analytical path when the model qualifies.Policy="numeric" to turn every analytical Jacobian off, for A/B comparison while debugging a suspect derivative.Policy="hybrid" only to prevent the fully-analytical path from being chosen while still using the analytical blocks that exist.<Jacobian Policy> is honoured by the SUNDIALS backend only. The Engine backend ignores it and logs a warning; Euler1imp decides purely on whether the model reports analytical Jacobians. Setting a policy on an Engine method changes nothing.Write the analytical Jacobian when the model is stiff and the run is long enough that the finite-difference cost matters — each numerical column costs one extra model evaluation per Newton step, which is why the payoff grows with the number of states and with how often the solver refactorises. It is not worth it for a model that already runs happily under an explicit method.
An analytical Jacobian that is subtly wrong does not fail loudly. It makes the Newton iteration converge slowly, or to the wrong answer, and the symptom looks like a stiffness problem. Run the built-in verification the first time you implement OdeJacobian(), and again after every change to the derivative expressions. It compares the assembled matrix against a forward-difference Jacobian and names the SimObject and port link behind each mismatch. Treat a clean verification run as part of "done", not as an extra.
The SUNDIALS integrator supports a Jacobian policy that controls how Jacobians are computed. Set it with the Policy attribute of <Jacobian>:
| Policy | Behavior | When to Use |
|---|---|---|
auto (default) | Uses analytical Jacobians if all SimObjects provide them (HasJacobians() returns true for all). Falls back to hybrid if some provide them (HasHybridJacobians() returns true). Falls back to SUNDIALS built-in numerical approximation if none provide them. | Recommended for most users. Automatically selects the best available strategy. |
analytical | Requires all SimObjects with states to provide analytical Jacobians. Fails with an error if any SimObject lacks OdeJacobian(). | When you need to guarantee that no numerical Jacobian approximation is used. |
hybrid | Uses the hybrid Jacobian assembly if at least one SimObject provides analytical Jacobians. Falls back to SUNDIALS built-in if none do. | When mixing legacy models (no Jacobians) with newer analytical models. See Hybrid Jacobian Mode below. |
numeric | Always uses SUNDIALS built-in numerical Jacobian approximation, ignoring any analytical Jacobians. | For debugging or comparing analytical vs numerical results. |
All virtual methods on SimObject that participate in implicit integration are listed below. Every method has a safe default — existing SimObjects require no changes. Methods are grouped by concern; implement them in roughly the order shown to build up capability incrementally.
| Method | Default | When to implement | Impl. difficulty | Performance impact |
|---|---|---|---|---|
| `HasJacobians()` | false | Whenever you implement OdeJacobian() | Trivial — 1 line | Gate: must return true before any analytical Jacobian is used |
| `OdeJacobian()` | No-op | Any stateful SO that benefits from analytical J = dF/dX | Medium — derive all partial derivatives | High — eliminates per-SO finite-difference cost; enables all structured solvers |
| `GetJacobianSparsity()` | Dense (-1) | SOs with many states and sparse internal coupling | Low–Medium — enumerate non-zeros in CSR | High for large SOs — unlocks sparse, band, block_tridiagonal, near_tridiagonal solvers and graph-coloring optimisation |
| `HasPortJacobians()` | false | Whenever you implement OutputPortJacobian / InputPortJacobian | Trivial — 1 line | Gate: must return true before off-diagonal port blocks are computed analytically |
| `OutputPortJacobian()` | No-op | Stateful SO whose output ports are consumed by another stateful SO | Low–Medium — derive d(port)/dX | High for coupled systems — fills off-diagonal Jacobian block analytically; avoids O(states) extra finite-difference perturbations per Newton step |
| `InputPortJacobian()` | No-op | Stateful SO that receives input from another stateful SO | Low–Medium — derive dF/d(input) | (same as OutputPortJacobian) |
| `InputOutputJacobian()` | No-op | Stateless SOs only — intermediaries that transform a signal without holding state | Low — single matrix of gains | Moderate — enables analytical chain-rule through stateless intermediaries (A→S→B); without it the link falls back to numerical perturbation |
| `GetPortDependencyPolicy()` | All | When you know an output port depends on only a subset of local states | Trivial — 1 line | — |
| `GetPortStateDependency()` | {} | Only when GetPortDependencyPolicy() returns UserSpecified | Low — return an index list | Moderate for large SOs — tighter sparsity pattern reduces graph-coloring perturbations |
HasHybridJacobians() is a system-level query on ModelStructure, not a SimObject method. It returns true when at least one SO provides analytical Jacobians (enabling the hybrid assembly path described in Hybrid Jacobian Mode).HasInputOutputJacobians() method on SimObject. The assembler detects stateless-chain capability by checking HasPortJacobians() on the intermediate object and calling InputOutputJacobian() directly. Simply override InputOutputJacobian() and return true from HasPortJacobians().To enable analytical Jacobian support in your SimObject, implement three methods:
Return true to indicate that OdeJacobian() is implemented:
Compute the Jacobian matrix J = dF/dX, where F is defined by your OdeFcn():
The matrix is stored in row-major order: element (i, j) is at index i * nStates + j.
For small systems, return -1 to indicate a dense Jacobian:
For large systems with sparse structure, return the sparsity pattern in CSR format to reduce memory and computation. See SimObject methods for details.
The pattern is model-wide, not per SimObject. A single SimObject returning -1, or a single port coupling whose objects do not provide port Jacobians, abandons the pattern for the whole system: the Jacobian is then assembled dense and every sparse linear solver (sparse, band, block_tridiagonal, near_tridiagonal) becomes unavailable, however sparse the rest of the model is. A model of 30 states or more that ends up on the dense path logs one warning at startup naming the SimObjects that declared nothing and counting the couplings without port Jacobians; below 30 states the dense solver is chosen regardless, so nothing is lost and nothing is logged.
The Van der Pol oscillator is a classic stiff system for \( \mu \gg 1 \):
\[ \begin{aligned} \dot{x}_0 &= x_1 \\ \dot{x}_1 &= \mu (1 - x_0^2) x_1 - x_0 \end{aligned} \]
To verify that your analytical Jacobian is correct, compare it against a numerical approximation. A simple finite difference check:
A well-implemented Jacobian should have errors on the order of \( 10^{-6} \) or less.
The rest of the implicit-integration material is split by audience:
1.9.8