|
FhSim
3.1.0
Marine systems simulation
|
|
|
FhSim
3.1.0
Marine systems simulation
|
FhSim assembles the global Jacobian from individual SimObject Jacobians in a two-part structure:
Since SimObjects are independent in the ODE function (each computes its own state derivatives), the on-diagonal blocks form a block-diagonal structure:
\[ J_{diag} = \begin{bmatrix} J_1 & 0 & \cdots & 0 \\ 0 & J_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & J_n \end{bmatrix} \]
where \( J_k \) is the Jacobian of the k-th SimObject.
When stateful SimObjects are connected via ports, the output of one SimObject becomes an input to another, creating off-diagonal Jacobian blocks. These are computed via the chain rule:
\[ \frac{\partial F_{consumer}}{\partial X_{source}} = \frac{\partial F_{consumer}}{\partial u} \cdot \frac{\partial u}{\partial X_{source}} \]
where \( u \) is the port signal. This requires the SimObjects to implement the port Jacobian API (see Cross-Object Port Jacobians below).
In hybrid mode, the global Jacobian is assembled as follows:
The hybrid Jacobian mode is designed for scenarios where some SimObjects provide analytical Jacobians while others do not. This is particularly useful when:
OdeJacobian() implementation) with newer analytical modelsWhen <Jacobian Policy="hybrid"> is set (or Policy="auto" and not all SimObjects provide Jacobians), FhSim:
HasJacobians() = true.HasPortJacobians() = true.Consider a system with a legacy OldHydraulicModel (no Jacobian) connected to a new NewMooringModel (with analytical Jacobian):
The NewMooringModel contributes its analytical Jacobian to the on-diagonal block, while the OldHydraulicModel block is filled numerically. Port connections between them are also handled numerically unless both objects implement port Jacobians.
When two stateful SimObjects are connected via ports, the receiving object's state derivatives depend on the sending object's states through the port signal. This creates off-diagonal blocks in the global Jacobian that must be computed for correct implicit integration.
To provide analytical port Jacobians, implement three additional methods:
Return true to indicate that OutputPortJacobian() and InputPortJacobian() are implemented:
Computes the Jacobian of an output port with respect to local states: \( \frac{\partial (\text{port\_value})}{\partial X_{local}} \).
The matrix is stored in row-major order with dimensions [portSize x nStates].
CommonComputations, ensure they are executed before computing the port Jacobian.Computes the sensitivity of ODE derivatives to an input port value: \( \frac{\partial F}{\partial (\text{input\_port})} \).
The matrix is stored in row-major order with dimensions [nStates x portSize].
For a port connection from SourceObject to ConsumerObject, the off-diagonal Jacobian block is:
\[ J_{block} = \frac{\partial F_{consumer}}{\partial u} \cdot \frac{\partial u}{\partial X_{source}} \]
FhSim computes this using the two methods above and assembles the result into the global Jacobian at the appropriate row/column offsets.
The SpringDamper outputs a force port that depends on its 2 states. The Vessel receives this force and its 6 state derivatives depend on it. The off-diagonal block is a 6x2 matrix computed as:
Both objects must implement HasPortJacobians(), OutputPortJacobian(), and InputPortJacobian() for this block to be computed analytically.
For efficient numerical Jacobian assembly, FhSim uses structural sparsity information to determine which Jacobian entries are potentially non-zero.
The PortDependencyPolicy enum controls how FhSim determines which states an output port depends on:
| Policy | Behavior |
|---|---|
All (default) | Assumes the port depends on all local states. Conservative but safe. May result in more numerical perturbations than necessary. |
UserSpecified | The SimObject author provides an exact list of state indices via GetPortStateDependency(). Enables tighter sparsity patterns and fewer perturbations. |
When GetPortDependencyPolicy() returns UserSpecified, this method provides the exact dependency list:
This affects the sparsity pattern for off-diagonal Jacobian blocks: only the specified state columns are marked as potentially non-zero, reducing the number of numerical perturbations needed.
When numerical finite-difference assembly is required, FhSim uses a greedy largest-first graph coloring algorithm on the column conflict graph. States that do not appear in the same row can be perturbed simultaneously, reducing the number of ODE function evaluations from O(n) to O(number of colors), where the number of colors is typically much smaller than n for sparse systems.
This is an internal optimization; users do not need to configure it directly. Providing accurate GetPortStateDependency() information improves the coloring efficiency.
The Engine Backward Euler path (BackwardEuler_i, legacy alias BackwardEuler_i) uses an Armijo backtracking line search during Newton iteration. This improves convergence for difficult nonlinear systems by adaptively reducing the Newton step size when the residual does not decrease sufficiently.
Parameters: maximum 5 backtracking iterations, sufficient decrease constant c = 1e-4.
FhSim detects singular Jacobians by checking for NaN or Inf values in the Newton step solution. When detected, the simulation aborts the current step with an error message:
This can occur when:
Self-connections (a SimObject's output port connected to its own input port) create algebraic loops that are incompatible with implicit integration. When implicit integration is active (HasJacobians() or HasHybridJacobians() returns true), self-connections are detected and reported as an error:
The on-diagonal Jacobian blocks are assumed to be block-diagonal — that is, each SimObject's state derivatives depend only on its own states, not on other SimObjects' states. This is inherent to the FhSim architecture where OdeFcn() is called per-object. Cross-object dependencies are handled exclusively through port connections (off-diagonal blocks).
Within each off-diagonal block (from a port connection), the Jacobian is treated as dense in the block bounds. Even if GetPortStateDependency() specifies a sparse dependency pattern, the entire submatrix from consumer rows to source columns is considered potentially non-zero. Structural sparsity within the block is only used to reduce numerical perturbations, not to skip zero entries during assembly.
Analytical port Jacobian entries with absolute value below a compression threshold are set to zero during assembly. This prevents near-zero numerical noise from filling the Jacobian. The threshold is defined internally and is typically on the order of 1e-15.
As noted in How analytical Jacobians are assembled, numerical Jacobian entries only fill positions that are still zero after analytical assembly. This means:
| System Size | Recommended Approach |
|---|---|
| < 50 states | Dense Jacobian (default) |
| 50-200 states, limited bandwidth | Banded Jacobian (<Jacobian Type="band"> + <LinearSolver Type="BAND">) |
| 50-200 states, chain topology | Block-tridiagonal (<LinearSolver Type="block_tridiagonal">) |
| 50-5000 states, general sparse | Sparse Jacobian with CSR format (<LinearSolver Type="sparse">, METIS ordering) |
| 50-5000 states, 2D/3D mesh (single SimObject) | <LinearSolver Type="sparse"> — METIS nested dissection is near-optimal for mesh sparsity |
| Multi-body, per-SimObject Jacobians | Component Block GS (<LinearSolver Type="component">) |
| > 5000 states | Sparse Jacobian, consider iterative solver |
For most marine simulation SimObjects, dense Jacobians are appropriate since individual objects typically have fewer than 20 states. Sparse Jacobians become beneficial when a single SimObject has many states (e.g., distributed parameter systems like cables or flexible bodies with many nodes).
For hybrid Jacobian mode, performance depends on the ratio of analytical to numerical entries. Maximizing the number of SimObjects with analytical Jacobians and port Jacobian support will reduce the number of numerical perturbations needed.
FhSim supports computing directional derivatives of outputs with respect to inputs, which is essential for sensitivity analysis, optimization, and FMI 3.0 co-simulation.
Two modes are available for computing directional derivatives:
| Mode | Description | Best When |
|---|---|---|
| Forward (mode=0) | Seed the inputs and solve \((I - \Delta t\,J)\,v = \Delta t\,\text{seed}\) | Few inputs, many outputs |
| Adjoint (mode=1) | Seed the outputs and solve the transposed system \((I - \Delta t\,J)^{T}\mu = \lambda\) | Few outputs, many inputs |
Both modes apply the same operator, one backward-Euler step of size \(\Delta t\), so they are exact transposes of one another: seeding one end and reading the other gives the same number whichever mode is used. Each mode costs one assembly of \(J\) and one LU factorisation regardless of how many inputs or outputs are named, so the choice between them is about which end is convenient to seed, not about cost.
GetDirectionalDerivative locates a port by its SimObject's state offset and does not read the port name, so its \(C\) and \(B\) are plain state selections. Use GetAnalyticalDirectionalDerivative for the derivative that follows port connections through the real OutputPortJacobian and InputPortJacobian factors.
The directional derivative is exposed through the FhSim DLL:
For a single time step with implicit integration, the input-output sensitivity is:
\[ \frac{dy}{du} = C \cdot (I - \Delta t \cdot J)^{-1} \cdot \Delta t \cdot B \]
where:
OutputPortJacobian()OdeJacobian()InputPortJacobian()Forward mode computes dy/du · dvKnown by solving (I - Δt·J) · v = Δt · B · dvKnown and reading C · v. Adjoint mode computes dvKnownᵀ · dy/du by solving (I - Δt·J)ᵀ · μ = Cᵀ · λ and then λ_input = Δt · Bᵀ · μ.
Neither mode adds port-coupling terms of its own. Every coupling — direct, composed through a stateless intermediate, and self-coupling — is already in J, which JacobianAssembler::OdeJacobian assembles, so taking them from J alone counts each one exactly once.
The directional derivative computation uses per-instance workspace buffers allocated during model initialization. Each FhSim DLL context has its own isolated workspaces, making the function safe for concurrent use across multiple FMU instances (FMI 3.0).
FhSim's analytical Jacobian assembler handles SimObjects that have no internal states but mediate forces or signals between stateful objects — so-called stateless intermediates. A typical example is a spring element that receives a position from object A, computes a force, and feeds it back to A and/or forward to object B.
The assembler recognises two one-hop patterns automatically:
| Pattern | Description |
|---|---|
| Cross-object chain (A→S→B) | Stateful object A drives stateless S, whose output drives stateful object B. The composed Jacobian contribution dF_B/dX_A = dF_B/dInput_B · dOutput_S/dInput_S · dPort_A/dX_A is added to the off-diagonal block J[B, A]. |
| Self-coupling (A→S→A) | Stateful object A drives stateless S, whose output feeds back into A itself. The composed contribution augments the diagonal block J[A, A]. |
Both patterns are discovered at model-build time in ModelAssemblyService::BuildPortJacobianLinks (src/engine/model/ModelAssemblyWiring.cpp) and stored on the assembly result:
portJacobianLinks — cross-object composed links (hasIntermediate = true).selfPortCouplingLinks — self-coupling links.Multi-hop chains (A→S1→S2→B) are not handled analytically. A warning is emitted to the log and the affected objects fall back to numerical Jacobian evaluation.
A stateless SimObject must implement InputOutputJacobian to participate in analytical chain assembly:
dY_dU[i + j*outSize] must contain dOutput[i] / dInput[j].
For a scalar spring with Force = -K * Position:
Return true from HasPortJacobians() to advertise the capability (there is no separate HasInputOutputJacobians() — the assembler checks HasPortJacobians() on intermediate objects).
The assembler tracks which stateless intermediates are fully covered by the chain-rule logic in m_coveredStatelessObjects (src/engine/model/JacobianAssembler.h). When deciding whether a stateful object needs a numerical Jacobian (needNumerical), only stateless objects not in this set trigger the fallback. This prevents spurious numerical evaluations when all stateless intermediates implement InputOutputJacobian.
GetJacobianSparsity adds the composed off-diagonal entries (A→S→B) and the extra diagonal columns implied by self-coupling (A→S→A) to the declared sparsity pattern. No additional configuration is required.
The test scenario tests/in/ImplicitIntegration/ImplicitIntegration_BDF_StatelessChain_in.xml demonstrates both patterns:
Analytical solutions for the test:
The test ImplicitIntegrationSim.BDF_StatelessChainCoupling in tests/integrator/methods/ImplicitIntegrationSim_test.cpp verifies both trajectories against these closed-form solutions with BDF integration and diagnostics-enabled Jacobian verification (<Diagnostics><Triggers Once="true"/><Settings><JacobianVerify/></Settings></Diagnostics>).
1.9.8