Recipe: Analytical Jacobian for NetStructure

Target file pair:

• src/net/NetStructure.h
• src/net/NetStructure.cpp

Reference implementations:

• fhsim_base/src/cable/Cable.cpp — HasJacobians, OdeJacobian, GetJacobianSparsity, computeSegmentStructuralJacobian
• fhsim_insights/architecture/jacobian-implementation-guide.md

---

Implementation Status

Completed

Block	Method	Notes
Position/velocity identity	OdeJacobian	dXdot_pos/dX_vel = I
Cable structural (pos)	ComputeCableStructuralJacobian	Axial + geometric stiffness
Cable structural (vel)	ComputeCableStructuralJacobian	Damping: c·e⊗e
Cable drag (vel)	ComputeCableDragVelocityJacobian	Tangential + normal
Cable drag (pos)	ComputeCableDragPositionJacobian	Direction rotation effect
Panel forces (all)	ComputePanelJacobianAD	Exact: AddNodeForces over sfh::ad::Dual<18>
External force ports	InputPortJacobian	1/m_inertia[node] diagonal
Sparsity	GetJacobianSparsity	Dense (-1), CSR deferred
Runtime guard	HasJacobians	Returns false for clamps/spheres/disks

Omitted (zero derivative or state-independent)

• Gravity/buoyancy (m_NodeWeight is constant)
• Bottom contact (only active near seabed; not gated since HasJacobians is setup-time only)
• Air/water blending of rhoWater across z=0 surface

Architecture Decisions

1. Cable Jacobian: Fully analytical. Covers Bound(spring + damping, 0, maxTension) saturation by returning zero when force is at either bound.

2. Panel Jacobian: Exact, by forward-mode algorithmic differentiation rather than by a closed form. The panel force function lives in src/net/NetElement3NForces.h templated on the scalar type; ComputePanelJacobianAD evaluates it once over sfh::ad::Dual<18>, whose eighteen gradient slots carry the three position and three velocity components of the three nodes, and reads the whole 18x9 block out of the result's gradients.

   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 and the same convention as the Net/Disk drag Jacobian.

   Cost is roughly fifteen to twenty times one force evaluation, since every operation carries eighteen gradient components. That is the same order as the 36 evaluations the central difference needed, so the gain here is exactness rather than speed. See issues MARE-0004 (which this resolves) and MARE-0026 (the enabling work).

3. Dense Jacobian: GetJacobianSparsity returns -1. For large nets, CSR would exploit the cable+panel connectivity graph, but the connectivity is determined at runtime and the dense form is simpler. Revisit if profiling shows this is the bottleneck.

---

Test Configuration

• Cable test (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.
• Panel test (SimplePanelNet.xml): 4-node square, 2 triangular panels, no cables. Verified with relTol=5e-3, absTol=0.1 across 3 states (nominal, non-planar, compressed). Max absolute error: 0.09 (double-FD round trip noise).
• Framework StealICs bug: The built-in JacobianChecker cannot be used because the integrator consumes initial conditions during BuildSimulation. Custom tests supply known states directly via IModelStructure.h interface.

---

Key Implementation Details

Cable Structural Force Law

r = pB - pA, L = |r|, e = r/L
s = (vA - vB) · e  (compression speed)
spring = (L - L0) * EA / L0  (can be negative)
damp = c * s  (c < 0 for dissipation)
rawTD = spring + damp
TD = Bound(rawTD, 0, maxTension)  ← actual force law
FA = TD * e, FB = -FA

When rawTD ≤ 0 or rawTD ≥ maxTension: force is flat, derivative = 0.

Panel AD Approach

ComputePanelJacobianAD(panel, pA, pB, pC, vA, vB, vC, vWater, addedDrag,
                       dFA_dpA, ..., dFC_dvC)

Seeds the eighteen inputs as sfh::ad::Dual<18> variables - positions of A, B and C in slots 0-8 and velocities in slots 9-17 - takes the panel's fixed properties from panel->PanelForceParams(), and calls net_element_forces::AddNodeForces once. Output block [group][node][j*3+k] is then just force[node][j].Grad(group*3+k).

The water velocity and the two scalar arguments 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.