rustytorch/crates/specialized/rtx-fsi/tests/fsi2_harness/overset_march.rs

//! P5 (§5.12): the coupled march on the overset fluid — the harness's
//! rigid phase, release, and per-step subiterated coupling (predictor,
//! IQN-ILS / Aitken passes each re-marching the fluid from the step's
//! snapshot, the acceptance rule of `march.rs`), without the embedded
//! march's rescue machinery (refuted, retired). Returns the harness's
//! `MarchResult` plus the composite's own counters.

use std::cell::RefCell;
use std::io::Write as _;

use nalgebra::Vector3;
use rtx_fea::analysis::{
    AnalysisConfig, ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis,
};
use rtx_fea::materials::{LinearElastic, MaterialDatabase};
use rtx_fea::mesh::{MaterialId, NodeId};
use rtx_fsi::{IqnIls, Subiterated};

use super::march::MarchResult;
use super::overset::OversetFluid;
use super::{clamp_left, env_or, median, mid_amp, BenchmarkCase};

/// The overset march's knobs (`RTX_<PREFIX>_*`).
#[derive(Debug, Clone)]
pub struct OversetMarchConfig {
    pub ny: usize,
    pub flag_nx: usize,
    pub t_release: f64,
    pub t_end: f64,
    pub subcycle: usize,
    pub tol_floor: f64,
    pub rtol: f64,
    pub stall_accept: f64,
    pub max_subiterations: usize,
    pub coupler: String,
    pub reuse: usize,
    pub initial_relaxation: f64,
    pub c1_interface: bool,
    pub predictor: String,
    /// Winslow sweeps per patch regeneration.
    pub sweeps: usize,
    /// Schwarz rounds per corrector (the P5 budget: 3).
    pub max_rounds: usize,
    pub csv_path: Option<String>,
    pub trace_steps: usize,
}

impl OversetMarchConfig {
    pub fn from_env(prefix: &str, d: OversetMarchConfig) -> OversetMarchConfig {
        let num = |k: &str, v: f64| env_or(&format!("RTX_{prefix}_{k}"), v);
        let text =
            |k: &str, v: &str| std::env::var(format!("RTX_{prefix}_{k}")).unwrap_or(v.into());
        OversetMarchConfig {
            ny: num("NY", d.ny as f64) as usize,
            flag_nx: num("FLAG_NX", d.flag_nx as f64) as usize,
            t_release: num("T_RELEASE", d.t_release),
            t_end: num("T_END", d.t_end),
            subcycle: num("SUBCYCLE", d.subcycle as f64) as usize,
            tol_floor: num("TOL_FLOOR", d.tol_floor),
            rtol: num("RTOL", d.rtol),
            stall_accept: num("STALL_ACCEPT", d.stall_accept),
            max_subiterations: num("MAX_SUBIT", d.max_subiterations as f64) as usize,
            coupler: text("COUPLER", &d.coupler),
            reuse: num("REUSE", d.reuse as f64) as usize,
            initial_relaxation: num("OMEGA0", d.initial_relaxation),
            c1_interface: num("C1", if d.c1_interface { 1.0 } else { 0.0 }) > 0.5,
            predictor: text("PREDICTOR", &d.predictor),
            sweeps: num("SWEEPS", d.sweeps as f64) as usize,
            max_rounds: num("MAX_ROUNDS", d.max_rounds as f64) as usize,
            csv_path: std::env::var(format!("RTX_{prefix}_CSV"))
                .ok()
                .or(d.csv_path),
            trace_steps: num("TRACE", d.trace_steps as f64) as usize,
        }
    }
}

/// The march's outcome: the harness's statistics, the composite's
/// counters, and the death (if any) instead of a panic.
pub struct OversetMarchResult {
    pub result: MarchResult,
    pub death: Option<(usize, f64, String)>,
    pub rounds_mean: f64,
    pub reclassified_mean: f64,
    pub fresh_mean: f64,
    pub regen_count: usize,
    pub regen_seconds: f64,
    pub fluid_seconds: f64,
    pub structure_seconds: f64,
    pub faces_used: usize,
}

pub fn run_march_overset(case: BenchmarkCase, config: &OversetMarchConfig) -> OversetMarchResult {
    let cfg = config.clone();
    let mut fluid = OversetFluid::build_case(case, cfg.ny, cfg.flag_nx, cfg.sweeps, cfg.max_rounds)
        .expect("overset fluid");
    let dt_fluid = fluid.dt_fluid;
    let dt = dt_fluid * cfg.subcycle as f64;
    let zero_d = vec![0.0; 2 * fluid.interface.wetted.len()];
    // Every setting the acceptance rule reads, printed once: P5-3 lost a
    // day to a floor of 2e-4 against the overnight marches' 1e-6.
    println!(
        "  coupling: {} (reuse {}, ω0 {}, c1 {}), floor {:.1e}, rtol {:.1e}, stall accept {:.1e}, max subit {}, predictor {}, s = {}, patch offset {} h × {} rows",
        cfg.coupler,
        cfg.reuse,
        cfg.initial_relaxation,
        cfg.c1_interface,
        cfg.tol_floor,
        cfg.rtol,
        cfg.stall_accept,
        cfg.max_subiterations,
        cfg.predictor,
        cfg.subcycle,
        std::env::var("RTX_FSI2O_PATCH_OFFSET").unwrap_or_else(|_| "6".into()),
        std::env::var("RTX_FSI2O_PATCH_ROWS").unwrap_or_else(|_| "12".into()),
    );

    // Phase 1: rigid flag to t_release (`RTX_FSI2O_LOAD=dir` replaces the
    // march with the saved state; `RTX_FSI2O_SAVE=dir` saves it).
    let start = std::time::Instant::now();
    let rigid_steps = (cfg.t_release / dt_fluid).round() as usize;
    if let Ok(dir) = std::env::var("RTX_FSI2O_LOAD") {
        let t = fluid.load(&dir).expect("load");
        assert!(
            (t - cfg.t_release).abs() < dt_fluid,
            "the saved state is at t = {t}, not t_release = {}",
            cfg.t_release
        );
        // The saved state is the FIELDS; the solvers' own warm state (the
        // Schwarz acceptor warm start, the fringe flux correction, the
        // embedded body's old volume fractions) is not in it, and a
        // release straight off the load carries a pressure-level
        // transient that the live march does not (P5-3: drag 27 vs 132
        // twenty steps in, lift ± 1000 vs ± 4 — enough to kick the flag's
        // thickness breathing at ny = 62). `RTX_FSI2O_LOAD_SETTLE=N`
        // (default 20) rigid steps rebuild that state on the steady rigid
        // flow; the clock is put back so the runs stay comparable.
        let settle: usize = std::env::var("RTX_FSI2O_LOAD_SETTLE")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(20);
        let (d0, l0) = fluid.measure_force();
        for _ in 0..settle {
            fluid.step().expect("settle fluid step");
        }
        fluid.solver.set_time(t);
        let (d1, l1) = fluid.measure_force();
        println!(
            "  settled the loaded state with {settle} rigid steps: wall drag {d0:.2} → {d1:.2}, lift {l0:.2} → {l1:.2}; clock back to t = {t:.4}"
        );
    } else {
        for _ in 0..rigid_steps {
            fluid.step().expect("rigid fluid step");
        }
        if let Ok(dir) = std::env::var("RTX_FSI2O_SAVE") {
            fluid.save(&dir).expect("save");
        }
    }
    // The release time is the fluid's clock (the rigid step count rounds).
    let t_release = fluid.time();
    fluid.commit_base();
    let (rigid_drag, rigid_lift) = fluid.measure_force();
    println!(
        "  {} OVERSET rigid phase: {rigid_steps} steps (dt {dt_fluid:.3e}) to t = {:.2} s in {:.0} s wall; wall drag {rigid_drag:.2} (rigid-flag reference {:.1}; the overset's own CFD2 at ny = 41: 137.8), lift {rigid_lift:.2}; rounds mean {:.2}",
        case.name,
        cfg.t_release,
        start.elapsed().as_secs_f64(),
        case.rigid_drag_reference,
        fluid.rounds_total.get() as f64 / fluid.correctors_total.get().max(1) as f64
    );

    // The flag: nonlinear Newmark stepper at the coupled dt (as march.rs).
    let mut db = MaterialDatabase::new();
    db.add_material(
        MaterialId(0),
        LinearElastic::new(case.e_s, case.nu_s).with_density(case.rho_s),
        None,
    );
    let analysis = NonlinearDynamicAnalysis::new(
        fluid.mesh.clone(),
        db,
        clamp_left(&fluid.mesh),
        dt,
        1,
        AnalysisConfig::default(),
    )
    .with_total_lagrangian()
    .with_convergence_criteria(ConvergenceCriteria {
        max_iterations: 60,
        ..ConvergenceCriteria::default()
    });
    // `RTX_FSI2O_NEWMARK_GAMMA=γ` (β = (γ + ½)²/4): numerical dissipation
    // of the flag's high-frequency modes (the thickness breathing the
    // body-fitted wall couples to at s = 1); the benchmark's average
    // acceleration (γ = ½) has none.
    let analysis = match std::env::var("RTX_FSI2O_NEWMARK_GAMMA")
        .ok()
        .and_then(|v| v.parse::<f64>().ok())
    {
        Some(g) => {
            let b = (g + 0.5).powi(2) / 4.0;
            println!("  Newmark γ = {g}, β = {b:.4}");
            analysis.with_newmark_parameters(g, b)
        }
        None => analysis,
    };
    let flag = RefCell::new(analysis.stepper().unwrap());
    let wetted_dofs: Vec<[usize; 2]> = fluid
        .interface
        .wetted
        .iter()
        .map(|&id| {
            let dofs = flag.borrow().node_dofs(id);
            [dofs[0], dofs[1]]
        })
        .collect();
    let a_dofs = flag.borrow().node_dofs(fluid.a_node);
    let extract = |state: &DynamicState| -> Vec<f64> {
        let mut d = vec![0.0; 2 * wetted_dofs.len()];
        for (k, dofs) in wetted_dofs.iter().enumerate() {
            d[2 * k] = state.displacement[dofs[0]];
            d[2 * k + 1] = state.displacement[dofs[1]];
        }
        d
    };
    let extract_velocity = |state: &DynamicState| -> Vec<f64> {
        let mut v = vec![0.0; 2 * wetted_dofs.len()];
        for (k, dofs) in wetted_dofs.iter().enumerate() {
            v[2 * k] = state.velocity[dofs[0]];
            v[2 * k + 1] = state.velocity[dofs[1]];
        }
        v
    };

    // Phase 2: release under the current load.
    let (nodal0, conservation0, faces0) = fluid.sample_load(&zero_d);
    flag.borrow_mut().set_nodal_forces(&nodal0);
    let mut flag_state = flag.borrow_mut().rest_state().unwrap();
    let mut committed_nodal = nodal0;
    let mut worst_conservation = conservation0;
    println!(
        "  release: {faces0} wall faces transferred (conservation defect {conservation0:.2e}); initial tip acceleration |a| = {:.3e}",
        (a_dofs.iter().map(|&k| flag_state.acceleration[k].powi(2)).sum::<f64>()).sqrt()
    );

    let fluid = RefCell::new(fluid);
    let mut iqn = (cfg.coupler == "iqn").then(|| {
        IqnIls::new(cfg.max_subiterations, 1.0)
            .unwrap()
            .with_reuse(cfg.reuse)
            .with_initial_relaxation(cfg.initial_relaxation)
            .unwrap()
    });

    let coupled_steps = ((cfg.t_end - cfg.t_release) / dt).round() as usize;
    let mut times = Vec::with_capacity(coupled_steps);
    let mut ux_series = Vec::with_capacity(coupled_steps);
    let mut uy_series = Vec::with_capacity(coupled_steps);
    let mut force_times = Vec::new();
    let mut drag_series = Vec::new();
    let mut lift_series = Vec::new();
    let (mut interval_drag, mut interval_lift) = (Vec::new(), Vec::new());
    let mut total_subiterations = 0usize;
    let mut max_subiterations = 0usize;
    let mut stalled_steps = 0usize;
    let mut retried_steps = 0usize;
    let mut worst_stall = 0.0_f64;
    let mut faces_used = faces0;
    let mut death: Option<(usize, f64, String)> = None;
    let mut csv = cfg.csv_path.as_ref().map(|p| {
        let mut f = std::fs::File::create(p).expect("csv");
        writeln!(f, "t,ux,uy,drag,lift").unwrap();
        f
    });
    let (t_fluid, t_structure) = (std::cell::Cell::new(0.0_f64), std::cell::Cell::new(0.0_f64));
    let prev_area = std::cell::Cell::new(fluid.borrow().shared.read().unwrap().area());
    let save_every: usize = std::env::var("RTX_FSI2O_SAVE_EVERY")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(0);
    let phase_start = std::time::Instant::now();

    for step in 0..coupled_steps {
        let d_n = extract(&flag_state);
        let v_n: Option<Vec<f64>> = cfg.c1_interface.then(|| extract_velocity(&flag_state));
        let d_predicted = if cfg.predictor == "kinematic" {
            let v = extract_velocity(&flag_state);
            d_n.iter().zip(&v).map(|(d, v)| d + dt * v).collect()
        } else {
            flag.borrow_mut().set_nodal_forces(&committed_nodal);
            let (predicted, _) = flag.borrow_mut().step(&flag_state).unwrap();
            extract(&predicted)
        };
        let saved = fluid.borrow().snapshot();
        type PassResult = (DynamicState, Vec<(NodeId, Vector3<f64>)>, f64, usize);
        let latest: RefCell<Option<PassResult>> = RefCell::new(None);
        let pass = |d_candidate: &[f64]| -> Vec<f64> {
            let fs = std::time::Instant::now();
            let mut fl = fluid.borrow_mut();
            fl.restore(&saved);
            fl.advance_subcycled(&d_n, d_candidate, cfg.subcycle, v_n.as_deref())
                .expect("fluid pass");
            let (nodal, conservation, faces) = fl.sample_load(d_candidate);
            t_fluid.set(t_fluid.get() + fs.elapsed().as_secs_f64());
            let ss = std::time::Instant::now();
            let mut flag_ref = flag.borrow_mut();
            flag_ref.set_nodal_forces(&nodal);
            let (candidate_state, _) = flag_ref.step(&flag_state).unwrap();
            t_structure.set(t_structure.get() + ss.elapsed().as_secs_f64());
            let d_new = extract(&candidate_state);
            if step < cfg.trace_steps {
                let residual: f64 = d_new
                    .iter()
                    .zip(d_candidate)
                    .map(|(a, b)| (a - b) * (a - b))
                    .sum::<f64>()
                    .sqrt();
                let load: f64 = nodal.iter().map(|(_, f)| f.norm()).sum();
                // The fluid's power on the flag at this candidate: Σ F · v
                // with v the candidate's mean interface velocity. Physical
                // damping is negative; a spurious velocity-proportional
                // reaction shows as positive power growing with v².
                let (mut power, mut fx, mut fy) = (0.0, 0.0, 0.0);
                for (k, (_, f)) in nodal.iter().enumerate() {
                    let vx = (d_candidate[2 * k] - d_n[2 * k]) / dt;
                    let vy = (d_candidate[2 * k + 1] - d_n[2 * k + 1]) / dt;
                    power += f.x * vx + f.y * vy;
                    fx += f.x;
                    fy += f.y;
                }
                let (pd, bd, div) = fl.last_defects.get();
                let (rough, flips, tn_max, tn_mean) = fl.wall_roughness();
                let (level, wall_flux, wall_area, poly_flux, poly_area) = fl.level_and_wall_flux();
                let area_rate = (poly_area - prev_area.get()) / dt;
                println!(
                    "    step {step} pass: |d_new − d_candidate| = {residual:.3e}, |d_new| = {:.3e}, nodal Σ|F| {load:.2} ΣF ({fx:+.2}, {fy:+.2}), power {power:+.3e} W/m, wall t_n roughness {rough:.2} ({flips} flips, max {tn_max:.1} mean {tn_mean:.1} Pa), patch p level {level:+.1} Pa, wall net flux {wall_flux:+.3e} m²/s over {wall_area:.3} m (polygon {poly_flux:+.3e}; flag area rate {area_rate:+.3e}), patch mass defect {pd:.2e} bg {bd:.2e} div {div:.1e}, {faces} faces",
                    d_new.iter().map(|v| v * v).sum::<f64>().sqrt()
                );
            }
            *latest.borrow_mut() = Some((candidate_state, nodal, conservation, faces));
            d_new
        };

        let increment: f64 = d_predicted
            .iter()
            .zip(&d_n)
            .map(|(a, b)| (a - b) * (a - b))
            .sum::<f64>()
            .sqrt();
        let tol_step = cfg.tol_floor.max(cfg.rtol * increment);
        let retry_at = (5.0 * tol_step).max(0.1 * increment);
        let acceptable = (cfg.stall_accept * tol_step).max(0.1 * increment);
        let mut outcome = if let Some(iqn) = iqn.as_mut() {
            iqn.set_tolerance(tol_step).unwrap();
            iqn.solve(&d_predicted, pass)
        } else {
            Subiterated::aitken(cfg.max_subiterations, tol_step)
                .unwrap()
                .solve(&d_predicted, pass)
        };
        if let (Err(e), Some(iqn_ref)) = (&outcome, iqn.as_mut()) {
            let recoverable = matches!(
                e,
                rtx_fsi::FsiError::CouplingNotConverged { residual, .. }
                | rtx_fsi::FsiError::CouplingDiverged { residual, .. }
                if *residual >= retry_at
            );
            if recoverable {
                iqn_ref.reset_history();
                retried_steps += 1;
                outcome = iqn_ref.solve(&d_predicted, pass);
            }
        }
        let t_now = t_release + (step + 1) as f64 * dt;
        match outcome {
            Ok(c) => {
                total_subiterations += c.iterations;
                max_subiterations = max_subiterations.max(c.iterations);
            }
            Err(
                rtx_fsi::FsiError::CouplingNotConverged {
                    iterations,
                    residual,
                    ..
                }
                | rtx_fsi::FsiError::CouplingDiverged {
                    iterations,
                    residual,
                },
            ) if residual < acceptable => {
                stalled_steps += 1;
                worst_stall = worst_stall.max(residual);
                total_subiterations += iterations;
                max_subiterations = max_subiterations.max(iterations);
            }
            Err(e) => {
                println!(
                    "  {} OVERSET DEATH at step {step} t = {t_now:.4}: {e:?} (increment {increment:.3e}, tol {tol_step:.3e}, acceptable {acceptable:.3e})",
                    case.name
                );
                death = Some((step, t_now, format!("{e:?}")));
                break;
            }
        }
        let (new_state, nodal, conservation, faces) = latest.borrow_mut().take().expect("pass ran");
        flag_state = new_state;
        committed_nodal = nodal;
        prev_area.set(fluid.borrow().shared.read().unwrap().area());
        fluid.borrow_mut().commit_base();
        // `RTX_FSI2O_SAVE_EVERY=N` (+ `RTX_FSI2O_SAVE`): the composite at
        // every N-th committed step, for the offline chain audit.
        if save_every > 0 && (step + 1) % save_every == 0 {
            if let Ok(dir) = std::env::var("RTX_FSI2O_SAVE") {
                let d_now = extract(&flag_state);
                let dd_now = extract_velocity(&flag_state);
                fluid
                    .borrow()
                    .save_instant(&dir, step + 1, &d_now, &dd_now)
                    .expect("save instant");
            }
        }
        worst_conservation = worst_conservation.max(conservation);
        faces_used = faces;
        let ux = flag_state.displacement[a_dofs[0]];
        let uy = flag_state.displacement[a_dofs[1]];
        times.push(t_now);
        ux_series.push(ux);
        uy_series.push(uy);
        let (drag_now, lift_now) = fluid.borrow().measure_force();
        interval_drag.push(drag_now);
        interval_lift.push(lift_now);
        if (step + 1) % 10 == 0 {
            let drag = median(&mut interval_drag);
            let lift = median(&mut interval_lift);
            interval_drag.clear();
            interval_lift.clear();
            force_times.push(t_now);
            drag_series.push(drag);
            lift_series.push(lift);
            if let Some(f) = csv.as_mut() {
                writeln!(f, "{t_now:.6},{ux:.6e},{uy:.6e},{drag:.6e},{lift:.6e}").unwrap();
            }
        } else if let Some(f) = csv.as_mut() {
            writeln!(f, "{t_now:.6},{ux:.6e},{uy:.6e},,").unwrap();
        }
        if (step + 1) % 500 == 0 {
            let window = &uy_series[uy_series.len().saturating_sub(500)..];
            let (w_mid, w_amp) = mid_amp(window);
            let fl = fluid.borrow();
            println!(
                "  t = {t_now:.3} s ({} steps): uy(A) = {uy:.3e} (window mid {w_mid:.3e} amp {w_amp:.3e}), drag {drag_now:.1} lift {lift_now:.1}, {:.1} subit/step, rounds mean {:.2}, reclassified/step {:.1}, regen {:.0} s of {:.0} s fluid, {:.0} s wall",
                step + 1,
                total_subiterations as f64 / (step + 1) as f64,
                fl.rounds_total.get() as f64 / fl.correctors_total.get().max(1) as f64,
                fl.reclassified_total.get() as f64 / (rigid_steps + (step + 1) * cfg.subcycle * 4).max(1) as f64,
                fl.regen_seconds.get(),
                t_fluid.get(),
                phase_start.elapsed().as_secs_f64()
            );
        }
    }
    let fl = fluid.borrow();
    let final_state_finite = flag_state.displacement.iter().all(|v| v.is_finite());
    let steps_done = times.len();
    OversetMarchResult {
        result: MarchResult {
            dt,
            coupled_steps: steps_done,
            times,
            ux: ux_series,
            uy: uy_series,
            force_times,
            drag: drag_series,
            lift: lift_series,
            rigid_drag,
            rigid_lift,
            mean_subiterations: total_subiterations as f64 / steps_done.max(1) as f64,
            max_subiterations,
            stalled_steps,
            retried_steps,
            worst_stall,
            worst_conservation,
            skipped: 0,
            spiked: 0,
            newton_rescues: flag.borrow().rescue_counts(),
            coupling_rescues: 0,
            coupling_rescue_failures: 0,
            rescue_records: Vec::new(),
            final_state_finite,
            elapsed: start.elapsed().as_secs_f64(),
        },
        death,
        rounds_mean: fl.rounds_total.get() as f64 / fl.correctors_total.get().max(1) as f64,
        reclassified_mean: fl.reclassified_total.get() as f64
            / (rigid_steps + steps_done * cfg.subcycle).max(1) as f64,
        fresh_mean: fl.fresh_total.get() as f64
            / (rigid_steps + steps_done * cfg.subcycle).max(1) as f64,
        regen_count: fl.regen_count.get(),
        regen_seconds: fl.regen_seconds.get(),
        fluid_seconds: t_fluid.get(),
        structure_seconds: t_structure.get(),
        faces_used,
    }
}