JEOD finding (for nasa/jeod): SIM_ground_contact initialization-state artifact in BodyRefFrame propagation order

[simnaut]

Summary (for the JEOD team)

While porting verif/SIM_ground_contact to a Rust reimplementation
(this repo, issues #88 / #205), we traced the CSV trajectory of
SET_test/RUN_contact_ground to an initialization-state artifact
in ContactGround::initialize_ground: a pre-propagation
GroundInteraction::initialize() → in_contact() call writes a spurious
impulsive force (~1.117 × 10¹⁰ N per vehicle) to subject->force,
which the integrator then consumes at stage 1 of the first RK4 step
(weight 1/6). The trajectory the simulation actually produces — vehicles
launching outward at ~93 km/s within 50 ms of "rest at the planet
surface" — is therefore not the steady-state contact algorithm at work;
it's an interaction between the timing of BodyRefFrame state
propagation and the order of P_DYN("initialization") jobs.

The artifact is fully reproducible from the source tree as committed;
no patch needed to observe it. We're filing this so the next reader of
SIM_ground_contact (us, the JEOD team, or anyone else porting it)
doesn't re-derive the analysis from scratch.

Code references

All paths relative to nasa/jeod:

• models/interactions/contact/verif/SIM_ground_contact/S_modules/contact.sm:70-71
  — P_DYN("initialization") contact.initialize_ground(...) /
  initialize_contact(...) job order.
• models/interactions/contact/verif/SIM_ground_contact/models/contact_ground/src/contact_ground.cc:146-184
  — ContactGround::initialize_ground iterates (ground_facets × contact_pairs) and calls GroundInteraction::initialize for each
  unique subject.
• models/interactions/contact/verif/SIM_ground_contact/models/contact_ground/src/ground_interaction.cc:37-42
  — GroundInteraction::initialize(...) { … in_contact(); } invokes
  in_contact() once during init.
• models/interactions/contact/verif/SIM_ground_contact/models/contact_ground/src/point_ground_interaction.cc::in_contact
  / line_ground_interaction.cc::in_contact — the algorithm reads
  vehicle_point->state.trans.position for the facet_pos term.
• models/dynamics/dyn_body/src/dyn_body_propagate_state.cc:159-198
  — compute_derived_state_forward(structure, mass_point, derived)
  populates derived.state.trans.position = structure.position + T_struct→inertial · mass_point.position. Until this runs,
  vehicle_point.state.trans.position is the default-constructed
  (0, 0, 0).

What the algorithm does at init vs runtime

point_ground_interaction.cc::in_contact (and the line variant) reads
vehicle_point->state.trans.position in two places — to assemble
vec for the terrain query and to compute facet_pos for the
subject_mag < ground_mag test:

Vector3::sum(vehicle_body->structure.state.trans.position,
             vehicle_point->state.trans.position,
             vec);
// ... terrain query gives ground_point ...
Vector3::transform(vehicle_point->state.rot.T_parent_this,
                   vehicle_point->state.trans.position,
                   facet_pos);
Vector3::sum(subject->contact_point, facet_pos, rel_state);
subject_mag = Vector3::vmag(rel_state);
ground_mag  = Vector3::vmag(ground_point);
if (subject_mag < ground_mag) { /* contact, force = k * (ground_point - rel_state) */ }

For a vehicle initialized via trans_init at inertial position
(R, 0, 0) (where R = 6378137 m is Earth's equatorial radius), with
mass center at the structure origin and identity attitude:

Init-time call (vp.state.trans.position == (0, 0, 0) — default-constructed)

• vec = (R, 0, 0) + (0, 0, 0) = (R, 0, 0) — interpreted as the
  vehicle's inertial position.
• Terrain query at altitude 0 gives ground_point body-frame ≈
  (R, 0, 0).
• contact_point = (1, 0, 0) for the point sphere (radius 1) or
  (2, 0, 0) for the line cylinder (length 2 cap-end + radius).
• facet_pos = T_parent_this · (0, 0, 0) = (0, 0, 0).
• rel_state = contact_point + 0 ≈ (1, 0, 0) →
  subject_mag ≈ 1 << R = ground_mag →
  contact triggers with depth ≈ R, force ≈ k · R ≈ 1.117 × 10¹⁰ N.

The force is added to subject->force — and, because
ContactSurface::collect_forces_torques only zeros facet.force after
its own derivative-class collection, the value persists through the
end of init.

Runtime call (vp.state.trans.position propagated to inertial position)

After compute_derived_state_forward runs (from
dyn_manager.initialize or the first propagate_state_from_structure
during the integrator loop), vehicle_point.state.trans.position
becomes (R, 0, 0). Now:

• vec = (R, 0, 0) + (R, 0, 0) = (2R, 0, 0) — the doubled inertial
  position (consistent because vp.pos is now the vehicle's own
  inertial position, summed with the structure's inertial position).
• Terrain query at altitude R gives the same body-frame
  ground_point ≈ (R, 0, 0).
• facet_pos = T · (R, 0, 0) = (R, 0, 0).
• rel_state = (1, 0, 0) + (R, 0, 0) = (R + 1, 0, 0) →
  subject_mag = R + 1 > R = ground_mag →
  no contact at any altitude (including h = 0).

So the steady-state algorithm reports no contact for any vehicle on or
above the surface — physically correct (a body just resting on the
ground has no penetration, no spring force).

Net trajectory

The init-time impulse F = k · R ≈ 1.117 × 10¹⁰ N is consumed
at stage 1 of the first RK4 step (ContactSurface::collect_forces_torques
zeros facet.force between stages, so stages 2–4 of step 1 see no
force, and the steady-state algorithm reports no contact thereafter).
Standard RK4 with force only at stage 1 yields:

Δv_step1 = (1/6) · F · dt / m
        ≈ (1/6) · 1.117 × 10¹⁰ · 0.01 / 200
        ≈ 9.31 × 10⁴ m/s

This matches the t = 0.05 CSV velocity (93081.39 m/s for the line,
93081.40 m/s for the point) to bit-precision (~picometers/second
in our pure-Rust port once we model the artifact explicitly with the
NIST-exact lbf/in conversion).

CSV contact_force doubling

The contact_surface.contact_force[0] column at t = 0 shows
22 339 651 946.5 N — exactly 2 × F. A plausible explanation:
P_DYN("derivative") jobs run once during dyn_manager.initialize
(after initialize_ground) before the data recorder samples at
t = 0. With vp.state.trans.position still at (0, 0, 0) at that
moment, check_contact_ground re-evaluates in_contact and adds a
second copy of F to subject->force. collect_forces_torques
then sums the two contributions into contact_surface.contact_force,
giving 2F — but only at the first sample, before the integrator's
first step propagates vp.state.trans.position. We haven't
verified this with a JEOD live-run trace; the trajectory match
demonstrates the integrator only consumes a single F.

Reproduction steps

To verify against JEOD itself:

1. Check out nasa/jeod and nasa/trick.
2. Build and run models/interactions/contact/verif/SIM_ground_contact/SET_test/RUN_contact_ground
   with the standard trick-CP flow (no patches).
3. Inspect log_contact_data.csv:
   • Row at t = 0: contact_force[0] ≈ 2.234 × 10¹⁰ N,
     position = (R, 0, 0), velocity = (0, 0, 0).
   • Row at t = 0.05: contact_force = 0, velocity ≈ 93 081 m/s.
4. Add a logging breakpoint (or MessageHandler::send) inside
   PointGroundInteraction::in_contact and
   LineGroundInteraction::in_contact printing
   vehicle_point->state.trans.position and subject_mag,
   ground_mag at every entry. Expected output:
   • First call (init via initialize_ground):
     vp.pos = (0, 0, 0), subject_mag ≈ 1,
     ground_mag ≈ 6 378 137, contact = true.
   • Second call (derivative class during dyn_manager.initialize,
     before the first integrator step):
     vp.pos = (0, 0, 0) still, contact = true again — second F
     added to subject->force.
   • Third call onward (every RK4 stage from step 1):
     vp.pos = (R, 0, 0), subject_mag = R + 1, contact = false.

The printout will confirm the "two pre-step in_contact calls write
into subject->force while vp.pos is still default-constructed"
hypothesis.

What our port did

To reproduce JEOD's CSV trajectory bit-for-bit, our Rust port models
the artifact explicitly:

• A Phase enum on compute_ground_contact_geometry with variants
  Initialization (pre-propagation, facet_pos = (0, 0, 0)) and
  SteadyState (post-propagation,
  facet_pos = T_inertial→body · vehicle_pos_inertial).
• Simulation::register_ground_contact_pair evaluates the
  Initialization phase at registration and stashes the impulse on
  the pair config.
• The contact-coupled RK4 closure consumes the pending impulse at
  stage 1 of the first step (mirroring collect_forces_torques's
  facet.force = 0 clearing) and evaluates Phase::SteadyState
  at every stage thereafter.

This achieves a position match of ~3 nm and velocity match of
~5 × 10⁻¹¹ m/s against SET_test/RUN_contact_ground over 10 s.

Suggested action for nasa/jeod

The verification sim runs without error and the unit tests don't
exercise the init-time in_contact call, so this is invisible to
JEOD's own CI. Possible mitigations the JEOD team may want to consider:

1. Document the behavior in
   models/interactions/contact/docs/tex/guide.tex so the next reader
   doesn't infer that SIM_ground_contact's ~93 km/s launch is
   physically meaningful.
2. Move GroundInteraction::initialize's in_contact() call to
   after dyn_manager.initialize's state-propagation phase (or guard
   it behind a check that vehicle_point has been propagated).
3. Reset subject->force to zero at the end of
   ContactGround::initialize_ground so the init-time impulse doesn't
   leak into the integrator. (Conceptually equivalent to calling
   collect_forces_torques once at the end of init.)

We're happy to share the Rust port's reproduction harness or the
diagnostic instrumentation we added if useful.

References

• This investigation was driven by porting work in our Rust
  reimplementation: this repo issue Port GroundFacet to enable SIM_ground_contact Tier 3 cross-validation #88, branch
  claude/prepare-issue-88-TYI0d. The Rust test
  crates/jeod_runner/tests/tier3_sim_contact.rs::tier3_contact_ground
  validates the bit-for-bit match against the committed reference CSV.
• Investigation summary commits: e698b7e, 051e2bb, 20a01e9 in
  this branch.