From 216df00961210482ab3ea540572629a1aa36da5b Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 23 May 2026 22:00:46 -0700 Subject: [PATCH 1/2] Port SIM_csr_compare gravity-acceleration octant sweep (#207) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit JEOD's SIM_csr_compare (the name is historical — JEOD 5.4 ships no CSR file) loads earth_GGM05C at degree/order 70, places a non-integrating vehicle, and teleports it through sign-permutations of an ISS-like position at t=1..5s, logging grav_pot/grav_accel/position. Because the body never integrates, this is a per-step gravity-evaluation cross-check, not a trajectory propagation. Adds: - A csr_compare_run01 entry + CSR_COMPARE_SNIPPET (DRAscii logger) to trick/generate_references.sh, and the regenerated reference CSV. - tier3_sim_csr_compare.rs evaluating our GGM05C 70x70 gravity via the production GravityControl::evaluate_accel_only path at each input-deck position, compared against JEOD's logged grav_accel. The sweep positions are JEOD source constants (not output), so computational independence holds. Earth orientation is the planet-generic default (no RNP object in the S_define), so the planet-fixed rotation is identity. - A PERMANENT_GAPS entry (csr_compare) in parity_coverage.rs: the non-integrating gravity-eval has no Bevy-parity counterpart. Our kernel reproduces JEOD's grav_accel to machine precision (observed per-component max <= 4.4e-16 m/s^2 on a ~7 m/s^2 signal); tolerances use a 1e-14 m/s^2 machine-epsilon floor. Closes #207. Co-Authored-By: Claude Opus 4.7 (1M context) --- .../csr_compare_run01_csr_compare.csv | 7 + .../tests/tier3_sim_csr_compare.rs | 173 ++++++++++++++++++ .../tests/parity_coverage.rs | 7 + trick/generate_references.sh | 21 +++ 4 files changed, 208 insertions(+) create mode 100644 crates/astrodyn_verif_jeod/test_data/csr_compare_run01_csr_compare.csv create mode 100644 crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs diff --git a/crates/astrodyn_verif_jeod/test_data/csr_compare_run01_csr_compare.csv b/crates/astrodyn_verif_jeod/test_data/csr_compare_run01_csr_compare.csv new file mode 100644 index 00000000..0bf020b8 --- /dev/null +++ b/crates/astrodyn_verif_jeod/test_data/csr_compare_run01_csr_compare.csv @@ -0,0 +1,7 @@ +sys.exec.out.time {s},vehicle.dyn_body.grav_interaction.grav_pot {m2/s2},vehicle.dyn_body.grav_interaction.grav_accel[0] {m/s2},vehicle.dyn_body.grav_interaction.grav_accel[1] {m/s2},vehicle.dyn_body.grav_interaction.grav_accel[2] {m/s2},vehicle.dyn_body.composite_body.state.trans.position[0] {m},vehicle.dyn_body.composite_body.state.trans.position[1] {m},vehicle.dyn_body.composite_body.state.trans.position[2] {m} + 0, 58908597.80559717, 5.513978457816486, -1.226857697152695, -6.620557498749319, -4292653.41, 955168.47, 5139356.57 + 1, 53346272.25214996, -6.896405135975418,-0.01816912198929226, -1.855999280727087, 7218634.798289895, 18998.64159785956, 1938152.473366886 + 2, 53346195.80388758, 6.896336787430112,-0.01812170195363711, -1.85605387366793, -7218634.798289895, 18998.64159785956, 1938152.473366886 + 3, 53346194.72486015, 6.896336669305437, 0.01817849766546527, -1.856053428009526, -7218634.798289895, -18998.64159785956, 1938152.473366886 + 4, 53346348.38993748, 6.896461906674234, 0.01821455600489722, 1.856021387246724, -7218634.798289895, -18998.64159785956, -1938152.473366886 + 5, 53346350.80331156, 6.896463663957428,-0.01808753348648805, 1.85602180947678, -7218634.798289895, 18998.64159785956, -1938152.473366886 diff --git a/crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs b/crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs new file mode 100644 index 00000000..c3d3e81b --- /dev/null +++ b/crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs @@ -0,0 +1,173 @@ +//! Tier 3: SIM_csr_compare — GGM05C 70×70 gravity-acceleration octant sweep. +//! +//! JEOD's `models/environment/gravity/verif/SIM_csr_compare/RUN_01` places a +//! **non-integrating** vehicle (translational & rotational dynamics off) and +//! teleports it through sign-permutations of the ISS-like position +//! `[x, y, z] = [7218634.798…, 18998.642…, 1938152.473…] m` at t = 1…5 s +//! (t = 0 logs the `iss_typical_state` position), logging the gravity +//! potential and acceleration each second. Because the body never integrates, +//! this is a **per-step gravity-evaluation cross-check, not a trajectory +//! propagation** — so rather than driving `Simulation::step()`, it evaluates +//! the same production gravity path the pipeline's gravity stage uses +//! (`GravityControl::evaluate_accel_only` on a GGM05C source) at each logged +//! position and compares against JEOD's `grav_accel`. (Sibling of the Tier 2 +//! `grav_geospherical` reference-vector test, but at degree/order 70 and +//! against a full sim's logged output rather than a static fixture.) +//! +//! The sweep positions are JEOD *source* constants (from the RUN_01 input +//! deck / `Modified_data/iss_typical_state.py`), not JEOD output — only the +//! `grav_accel` being compared against is JEOD output, so computational +//! independence holds. The Earth orientation in this sim is the default +//! planet-generic (no RNP/rotation object in the S_define), so the +//! planet-fixed frame coincides with the integration frame and no rotation is +//! applied to the evaluation (`t_inertial_pfix = None`). +//! +//! Reference CSV regenerated via `cargo xtask regenerate-tier3` (the +//! `csr_compare_run01` entry in `trick/generate_references.sh`). + +#![allow( + clippy::float_cmp, + reason = "position sanity check asserts exact recovery of the literal JEOD input-deck constants logged back into the CSV" +)] + +use astrodyn::{gravity_fixtures, GravityControl, GravityGradient, GravityModel, GravitySource}; +use astrodyn_verif_jeod::crossval::CrossvalReport; +use astrodyn_verif_jeod::tier3_csv::test_data_path; +use glam::{DMat3, DVec3}; + +const DEGREE: usize = 70; +const ORDER: usize = 70; + +/// One logged record from `csr_compare_run01_csr_compare.csv`. +struct Record { + grav_accel: DVec3, + position: DVec3, +} + +/// Parse the DRAscii CSV (header row + 6 data rows). Columns: +/// time, grav_pot, grav_accel[0..2], position[0..2]. +fn load_csr_compare_csv() -> Vec { + let path = test_data_path("csr_compare_run01_csr_compare.csv"); + let text = std::fs::read_to_string(&path).unwrap_or_else(|e| { + panic!( + "csr_compare reference not found at {}: {e}.\n\ + Regenerate with `JEOD_HOME=… cargo xtask regenerate-tier3`.", + path.display() + ) + }); + text.lines() + .skip(1) // header + .filter(|l| !l.trim().is_empty()) + .map(|line| { + let f: Vec = line + .split(',') + .map(|s| { + s.trim() + .parse::() + .unwrap_or_else(|e| panic!("bad CSV field {s:?}: {e}")) + }) + .collect(); + assert_eq!(f.len(), 8, "expected 8 columns, got {}", f.len()); + Record { + grav_accel: DVec3::new(f[2], f[3], f[4]), + position: DVec3::new(f[5], f[6], f[7]), + } + }) + .collect() +} + +/// The six sweep positions, reconstructed from the RUN_01 input deck +/// (JEOD source — `Modified_data/iss_typical_state.py` for t=0, then the +/// `trick.add_read` octant teleports for t=1…5). +fn expected_positions() -> [DVec3; 6] { + // iss_typical_state.py + let t0 = DVec3::new(-4_292_653.41, 955_168.47, 5_139_356.57); + // input.py octant constants + let x = 7_218_634.798_289_895; + let y = 18_998.641_597_859_56; + let z = 1_938_152.473_366_886; + [ + t0, + DVec3::new(x, y, z), + DVec3::new(-x, y, z), + DVec3::new(-x, -y, z), + DVec3::new(-x, -y, -z), + DVec3::new(-x, y, -z), + ] +} + +#[test] +fn tier3_csr_compare_gravity_octants() { + let records = load_csr_compare_csv(); + assert_eq!(records.len(), 6, "expected 6 logged records (t=0..5)"); + + // Build a GGM05C spherical-harmonics source through the `astrodyn` + // facade (the verif crate reads physics through `astrodyn` only). + let data = gravity_fixtures::load_ggm05c(); + assert!( + data.degree >= DEGREE && data.order >= ORDER, + "GGM05C fixture must cover {DEGREE}×{ORDER}" + ); + let mu = data.mu; + let source = GravitySource { + mu, + model: GravityModel::SphericalHarmonics(Box::new(data)), + }; + // Earth orientation is the planet-generic default (no RNP object in + // the S_define), so the planet-fixed frame coincides with the + // integration frame: the required planet-fixed rotation is identity. + // The source-id type is irrelevant here (evaluate takes the source + // explicitly), so use the usize-id form the kernel tests use. + let control = + GravityControl::::new_nonspherical(0_usize, DEGREE, ORDER, GravityGradient::Skip); + let t_inertial_pfix = DMat3::IDENTITY; + + let positions = expected_positions(); + let mut max_accel_err = [0.0_f64; 3]; + + for (i, rec) in records.iter().enumerate() { + // Sanity: the position we feed our kernel is the same one JEOD + // logged (confirms row alignment + that our input constants match). + let dp = (rec.position - positions[i]).abs(); + assert!( + dp.max_element() < 1e-3, + "record {i}: logged position {:?} != input-deck position {:?}", + rec.position, + positions[i] + ); + + let ga = + control.evaluate_accel_only(&source, positions[i], Some(&t_inertial_pfix), 0.0, false); + let ae = (ga.grav_accel - rec.grav_accel).abs(); + max_accel_err[0] = max_accel_err[0].max(ae.x); + max_accel_err[1] = max_accel_err[1].max(ae.y); + max_accel_err[2] = max_accel_err[2].max(ae.z); + } + + let mut report = CrossvalReport::compute("tier3_csr_compare_gravity_octants", &[], &[]); + report.add_extra("accel_err_x", max_accel_err[0], "m/s^2"); + report.add_extra("accel_err_y", max_accel_err[1], "m/s^2"); + report.add_extra("accel_err_z", max_accel_err[2], "m/s^2"); + report.write(); + + // Tolerances set to 1.05× observed max per component (CLAUDE.md policy). + const ACCEL_TOL: [f64; 3] = [ACCEL_TOL_X, ACCEL_TOL_Y, ACCEL_TOL_Z]; + for k in 0..3 { + assert!( + max_accel_err[k] < ACCEL_TOL[k], + "grav_accel component {k} error {:.6e} m/s^2 exceeds tolerance {:.6e}", + max_accel_err[k], + ACCEL_TOL[k] + ); + } +} + +// Our GGM05C 70×70 gravity kernel reproduces JEOD's logged `grav_accel` +// to machine precision (observed per-component max: x = 0, y ≈ 3.5e-18, +// z ≈ 4.4e-16 m/s² — i.e. ~15 significant figures on a ~7 m/s² signal). +// 1.05× observed would be a sub-ULP / zero tolerance, so these use a +// uniform 1e-14 m/s² machine-epsilon floor that proves bit-level +// agreement while absorbing last-ULP variation in the 70×70 harmonic sum. +const ACCEL_TOL_X: f64 = 1.0e-14; +const ACCEL_TOL_Y: f64 = 1.0e-14; +const ACCEL_TOL_Z: f64 = 1.0e-14; diff --git a/crates/astrodyn_verif_parity/tests/parity_coverage.rs b/crates/astrodyn_verif_parity/tests/parity_coverage.rs index db6d5a59..426df867 100644 --- a/crates/astrodyn_verif_parity/tests/parity_coverage.rs +++ b/crates/astrodyn_verif_parity/tests/parity_coverage.rs @@ -203,6 +203,13 @@ const PERMANENT_GAPS: &[(&str, &str)] = &[ "drag_analytical", "analytical drag verification — out of trait scope (no propagation)", ), + ( + "csr_compare", + "GGM05C 70×70 gravity-acceleration octant cross-check (#207) — the \ + JEOD vehicle is non-integrating (teleported through octant \ + positions), so this is a per-step gravity-evaluation comparison, \ + not a `Simulation::step()` trajectory; no Bevy parity counterpart.", + ), ]; /// Union of the two gap arrays. The coverage check unions both into the diff --git a/trick/generate_references.sh b/trick/generate_references.sh index 969b374f..53aa9c54 100755 --- a/trick/generate_references.sh +++ b/trick/generate_references.sh @@ -523,6 +523,23 @@ for ii in range(3): trick.add_data_record_group(dr) ' +# ── ASCII logging snippet for SIM_csr_compare (gravity-accel octant sweep) ── +# Non-integrating vehicle teleported through octant positions at t=1..5; +# logs gravity potential + acceleration + position at 1 Hz (matches the +# DRBinary "gravity_compare" group the input.py defines). The cross-check +# (#207) evaluates our GGM05C 70x70 accel at these positions vs grav_accel. +CSR_COMPARE_SNIPPET=' +dr = trick.sim_services.DRAscii("csr_compare_ASCII") +dr.set_cycle(1.0) +dr.freq = trick.sim_services.DR_Always +dr.add_variable("vehicle.dyn_body.grav_interaction.grav_pot") +for ii in range(3): + dr.add_variable("vehicle.dyn_body.grav_interaction.grav_accel[" + str(ii) + "]") +for ii in range(3): + dr.add_variable("vehicle.dyn_body.composite_body.state.trans.position[" + str(ii) + "]") +trick.add_data_record_group(dr) +' + # ── ASCII logging snippet for SIM_2_SHADOW_CALC (eclipse geometry) ── # Logs vehicle position, flux magnitude, and radiation force/torque. SHADOW_CALC_SNIPPET=' @@ -906,6 +923,10 @@ PID_INTEG=$LAST_BG_PID throttled_bg run_sim_with_ascii "models/interactions/radiation_pressure/verif/SIM_3_ORBIT" "SET_test/RUN_radiation" "srp_orbit_radiation" "$SRP_ORBIT_SNIPPET" PID_SRP_ORBIT=$LAST_BG_PID +# Group 9b: SIM_csr_compare (GGM05C 70x70 gravity-accel octant sweep, #207) +throttled_bg run_sim_with_ascii "models/environment/gravity/verif/SIM_csr_compare" "SET_test/RUN_01" "csr_compare_run01" "$CSR_COMPARE_SNIPPET" +PID_CSR_COMPARE=$LAST_BG_PID + # Group 10: SIM_torque_compare_simple (high-resolution gravity torque, 6 runs) run_torque_compare_simple_group() { local sim_dir="models/interactions/gravity_torque/verif/SIM_torque_compare_simple" From 07cd6684c0ecfdd0da80ea924f5db08c0111fcd0 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 23 May 2026 23:14:37 -0700 Subject: [PATCH 2/2] Fix SIM_csr_compare regen wait + identity-rotation doc Address Copilot review on #613, verified against JEOD source (models/environment/gravity/verif/SIM_csr_compare): - generate_references.sh: PID_CSR_COMPARE was launched but never waited on, so the script could report success before SIM_csr_compare finished (risking incomplete/missing CSV). Add the wait alongside PID_DYNCOMP. - tier3_sim_csr_compare.rs: module doc claimed `t_inertial_pfix = None`, but the code passes Some(&DMat3::IDENTITY) (evaluate_accel_only panics on None for non-spherical gravity). JEOD evaluates gravity with pfix->state.rot.T_parent_this left at its uninitialized identity since the S_define has no RNP object (spherical_harmonics_calc_nonspherical.cc :109/440), so identity is the correct rotation. Fix the doc to match. Co-Authored-By: Claude Opus 4.7 (1M context) --- crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs | 7 +++++-- trick/generate_references.sh | 1 + 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs b/crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs index c3d3e81b..bf2af210 100644 --- a/crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs +++ b/crates/astrodyn_verif_jeod/tests/tier3_sim_csr_compare.rs @@ -19,8 +19,11 @@ //! `grav_accel` being compared against is JEOD output, so computational //! independence holds. The Earth orientation in this sim is the default //! planet-generic (no RNP/rotation object in the S_define), so the -//! planet-fixed frame coincides with the integration frame and no rotation is -//! applied to the evaluation (`t_inertial_pfix = None`). +//! planet-fixed frame coincides with the integration frame: JEOD evaluates +//! gravity with `pfix->state.rot.T_parent_this` left at its uninitialized +//! identity (see `spherical_harmonics_calc_nonspherical.cc`), so we pass the +//! identity rotation `t_inertial_pfix = Some(&DMat3::IDENTITY)` — required +//! because `evaluate_accel_only` panics on `None` for non-spherical gravity. //! //! Reference CSV regenerated via `cargo xtask regenerate-tier3` (the //! `csr_compare_run01` entry in `trick/generate_references.sh`). diff --git a/trick/generate_references.sh b/trick/generate_references.sh index 53aa9c54..1cdcc463 100755 --- a/trick/generate_references.sh +++ b/trick/generate_references.sh @@ -2611,6 +2611,7 @@ echo "=== Waiting for all sim groups to complete ===" FAIL=0 wait $PID_DYNCOMP || { echo "WARN: SIM_dyncomp group had failures"; FAIL=1; } +wait $PID_CSR_COMPARE || { echo "WARN: SIM_csr_compare group had failures"; FAIL=1; } wait $PID_ORBINIT || { echo "WARN: SIM_orbinit group had failures"; FAIL=1; } wait $PID_ORBELEM || { echo "WARN: SIM_OrbElem failed"; FAIL=1; } wait $PID_LVLH || { echo "WARN: SIM_LVLH group had failures"; FAIL=1; }