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)	ComputePanelJacobianFD	Internal FD of AddNodeForces
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: Uses internal forward finite differences on AddNodeForces. This is correct-by-construction (matches the actual force function derivatives) and avoids the extreme complexity of differentiating through RMatrixFromThreePoints, CalcLocalUVComponents, Reynolds-dependent Cn, and coordinate rotations analytically. The 18-evaluation FD is computed once per panel per Jacobian assembly, which is still cheaper than the integrator's full system-level FD (which would perturb all N states independently).
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 FD Approach

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

Calls panel->AddNodeForces(...) at the base state, then 18 times with forward-perturbed components (eps = 1e-7 * max(1, |x_k|)). Each call zeros the force output buffers first (AddNodeForces adds to them).