|
FhSim
3.1.0
Marine systems simulation
|
|
|
FhSim
3.1.0
Marine systems simulation
|
src/net/NetStructure.hsrc/net/NetStructure.cppReference implementations:
fhsim_base/src/cable/Cable.cpp — HasJacobians, OdeJacobian, GetJacobianSparsity, computeSegmentStructuralJacobianfhsim_insights/architecture/jacobian-implementation-guide.md| Block | Method | Notes |
|---|---|---|
| Position/velocity identity | OdeJacobian | dXdot_pos/dX_vel = I |
| Cable structural (pos) | ComputeCableStructuralJacobian | Axial + geometric stiffness, from the axial law (cable/subroutines/CableAxialLaw.h: AxialTension, EndBlocks) |
| Cable structural (vel) | ComputeCableStructuralJacobian | Damping: c·√eFactor·e⊗e, 0 past the compression-speed clamp (MARE-0385) |
| Cable drag (vel) | ComputeCableDragVelocityJacobian | Tangential + normal |
| Cable drag (pos) | ComputeCableDragPositionJacobian | Direction rotation effect |
| Panel forces (all) | ComputePanelJacobianAD | Exact: AddNodeForces over sfh::ad::Dual<12> (NetElement3N::NodeForceJacobian, MARE-0205) |
| External force ports | InputPortJacobian | 1/m_inertia[node] diagonal |
| Sparsity | GetJacobianSparsity | CSR: the states of a node and of the nodes it shares a panel or cable with; the velocity rows of a force-input node dense (MARE-0209) |
| Runtime guard | HasJacobians | Returns false for clamps/spheres/disks, and while a node is in bottom contact at the last OdeFcn state (MARE-0134) |
m_NodeWeight is constant)cable/subroutines/CableAxialLaw.h, shared with InternalCable, plan E5). The law holds the compression-speed clamp, the length bounds and the bound of T, and returns T with its two slopes; EndBlocks turns them into the 3x3 blocks.Panel Jacobian: Exact, by forward-mode algorithmic differentiation rather than by a closed form. The panel force function lives in src/net/NetElement3NForces.h and the panel load laws in fhsim_environment/hydrodynamics/, all templated on the scalar type; ComputePanelJacobianAD evaluates it once over sfh::ad::Dual<12> (NetElement3N::NodeForceJacobian), whose nine position slots carry the three position components of the three nodes and whose three velocity slots carry all three node velocities (they enter only through their mean and linearly, see Panel AD Approach below), and builds the whole 18x9 block from the result's gradients (MARE-0205).
This is correct by construction, which the finite difference it replaced was not: there is no step size, so no truncation error and no dependence on the state's magnitude. It avoids the extreme difficulty of differentiating RMatrixFromThreePoints, CalcLocalUVComponents, the Reynolds-dependent Cn and the coordinate rotations by hand, while giving the same answer a correct hand derivation would. At a kink - the Max in the twine tension, the drag saturation floors, the piecewise Ersdal-Faltinsen branch - it returns the derivative of the branch actually taken, which is a valid subgradient.
Cost, measured on 2026-09-28 (WSL2, Release, one 20 x 20-mesh panel, the four laws): about 4.1 µs per panel over Dual<12>, 20 to 27 times the 140-210 ns of one double force evaluation. The eighteen-slot form it replaced (MARE-0205) took about 7.5 µs, 37 to 49 times (this file said fifteen to twenty times until then, which the deep review of 2026-09-28 measured to be wrong). The gain over the 36 evaluations of the central difference is exactness more than speed. See issues MARE-0004 (which this resolves) and MARE-0026 (the enabling work).
GetJacobianSparsity returns the CSR pattern of the connectivity graph: node n's rows couple to the six states of n and of every node it shares a panel or a cable element with, and the three velocity rows of a node with a force input are dense, since a connected object with direct feedthrough may depend on any state (MARE-0209). The solvers keep only the declared entries, so tests/NetStructure_Sparsity_Test.cpp differences OdeFcn and requires every non-negligible entry to be in the pattern.SimpleNet.xml): 4-node rhombus, 6 cables, no panels. Verified with relTol=1e-4, absTol=0.01 across 3 states (tension, compressed, mixed). Max absolute error: 1.1e-3.LowStrainPanelNet.xml): 4-node square, 2 triangular panels of 20 meshes per edge, no cables, twines stretched about 1 %. Verified with relTol=1e-3, absTol=1e-3 at step 1e-6 (MARE-0027). Max absolute error over 40 seeds: 5.1e-5. It also passes at the checker's default step 1e-8 with absTol=0.01 (max 0.0040).tests/NetStructure_Jacobian_Test.cpp; diagnosed on the branch diag/jacobian-twine d4a60da, never merged, with no AD defect found). A test that runs fhsim::test::RunTest with jacobianMode Auto on a net fixture therefore either passes a jacobianConfigOverride with perturbationEps = 1e-6 (the PanelNet_Jacobian settings above) or sets jacobianMode = fhsim::test::TestJacobianMode::No. A Jacobian check that fails at step 1e-8 alone is not evidence of a defect.bool into the arithmetic, and keep kinks as branches. The full list is in DEVELOPING.md, "Automatic-differentiation rules".JacobianChecker cannot be used because the integrator consumes initial conditions during BuildSimulation. Custom tests supply known states directly via IModelStructure.h interface.cable/subroutines/CableAxialLaw.h (cable_axial_law), the same law as CableElDynStiff::AddStructuralForces:
AxialTension returns T, ∂T/∂L (0 outside the length bounds) and ∂T/∂s (0 past the clamp, MARE-0385); both slopes are 0 where T is at a bound. EndBlocks gives dFA/dpB = T (I - e eᵀ)/L + ∂T/∂L e eᵀ + ∂T/∂s e (v_rel - s e)ᵀ/L and dFA/dvA = ∂T/∂s e eᵀ. The mean tension is held through the step, so it has no derivative (MARE-0359). Where T is at maxTension only the turning term T (I - e eᵀ)/L remains (MARE-0387); an element without tension adds nothing.
Calls panel->NodeForceJacobian, which seeds the positions of A, B and C as sfh::ad::Dual<12> variables in slots 0-8 and all three node velocities in the shared slots 9-11 (MARE-0205; the eighteen-slot form until then), and calls AddNodeForces<Dual> once, with the structure's panel load law and fluid. That computes the panel frame and mesh opening angle, the solidity from the panel's rule (NetSolidity.h; a mesh-opening solidity carries its derivative), and net_element_forces::AddNodeForces, whose hydrodynamic part is one evaluation of the load law (MARE-0094). The position blocks are force[node][j].Grad(group*3+k). The velocities enter only through their mean (the law) and linearly (the damping and the added drag, CalcDampingForces), so the shared slots give node n the sum G_n = L'/3 - (a/3) I over the three velocities, with L' the law's derivative with respect to the mean velocity and a the added drag; each velocity block is then dF_n/dv_m = (G_n + (a/3) I)/3 + c (1/3 - delta_nm) I - (a/3) delta_nm I, with c the damping coefficient. NetElement3NForces.TwelveSlotJacobianIsTheEighteenSlotOne* pins it to the eighteen-slot evaluation within 1e-12 of the block's largest entry, for every law, both solidity rules and a positive added drag.
The water velocity (the waked flow at the centroid) and the added drag enter as constants and contribute a zero gradient, as does the constant nodal weight that AddNodeForces adds; the finite difference this replaced relied on that weight cancelling in the difference instead. The force accumulators are default-constructed, so they start at zero as before.
1.9.8