//! P5 (`docs/overset_metal_campaign.md` §5.12): the FSI2 harness's FLUID
//! SIDE on the overset — the background without a body, the cylinder–flag
//! O-grid regenerated around the deformed flag every time the interface
//! moves (`set_patch_mesh`: the overlap rebuilt, fresh cells refilled, the
//! fringe re-stamped, the flux balance on), the patch's wall velocity from
//! the interface velocities (nearest wetted segment, linear along it), and
//! the load from the patch's wall faces (pressure + full-stress traction
//! per face into `WettedSurface::transfer_load`) — no probes, no spike
//! clamp, no smoothing. The structure side is the harness's, unchanged.

use std::cell::Cell;
use std::sync::{Arc, RwLock};

use nalgebra::Vector3;
use rtx_cfd::mesh::PatchSide;
use rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed_tip;
use rtx_cfd::solvers::incompressible::{
    AleBoundaries, ConvectionScheme, CurvilinearParameters, CurvilinearPisoSolver,
    EmbeddedParameters, EmbeddedPisoSolver, FlowField, MgPrecision, NormalDiffusion, OversetField,
    OversetParameters, OversetPisoSolver, OversetResult, OversetSolverState, PatchConvection,
    PatchField, PoissonSolverKind, RobinWall, SideBoundary,
};
use rtx_cfd::{CfdConfig, CfdResult};
use rtx_fea::mesh::{Mesh, NodeId};
use rtx_fsi::{FluidFace, WettedSurface};

use super::{BenchmarkCase, FLAG_X1, H, Interface, L, NU_F, RHO_F, flag_mesh, inflow_for};

const CYL_CENTRE: [f64; 2] = [0.2, 0.2];
const CYL_R: f64 = 0.05;
const FLAG_T: f64 = 0.01;
/// The junction fillet: a fixed 5 mm at every resolution (§5.11).
const FILLET: f64 = 0.5 * 0.41 / 41.0;
const PATCH_ROWS: usize = 12;
const PATCH_STRETCH: f64 = 4.0;

/// The patch's thickness in units of h (`RTX_FSI2O_PATCH_OFFSET`, 6).
/// At the default the overlap band sits a fixed NUMBER of cells from the
/// wall and moves inward in metres with refinement (ny 41 / 62 / 82 →
/// 60 / 40 / 30 mm); P5-3 holds it in metres across the ladder instead.
fn patch_offset_h() -> f64 {
    std::env::var("RTX_FSI2O_PATCH_OFFSET")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(6.0)
}

/// The tip's corner radius (`RTX_FSI2O_TIP_CORNER`, metres; the recorded
/// outline is the full semicircle, corner = t = 0.01): P5-3 option B — the
/// benchmark's flat tip with rounded corners (2.5 mm) against the
/// semicircle, gated on the lift phase and the per-period amplitude.
pub fn tip_corner() -> f64 {
    std::env::var("RTX_FSI2O_TIP_CORNER")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(FLAG_T)
}

/// The outline's fillet radius (`RTX_FSI2O_FILLET`, metres; 5 mm = t/2 in
/// every recorded run, so the tip is a full half-round and the root
/// carries 5 mm fillets — the reference's flag is a sharp rectangle).
/// P5-3's problem-definition probe: the tip's rounding sets how the flag
/// sheds, i.e. the wake's strength, which the h-ladder says is the
/// excitation that keeps growing.
pub fn fillet() -> f64 {
    std::env::var("RTX_FSI2O_FILLET")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(FILLET)
}

/// The patch's across stretch (`RTX_FSI2O_PATCH_STRETCH`, the geometric
/// ratio of the outer to the wall cell, 4.0 in every recorded run): with
/// 12 rows over 6 h the wall cell is 6 h (r − 1)/(r¹² − 1), r¹¹ = stretch —
/// 0.227 h / 0.144 h / 0.088 h at 4 / 8 / 16 — P5-3's wall-normal
/// resolution knob at fixed h (the rows ladder does not build: the overlap
/// needs outer cells ≈ h/2 or larger).
pub fn patch_stretch() -> f64 {
    std::env::var("RTX_FSI2O_PATCH_STRETCH")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(PATCH_STRETCH)
}

/// The background's convection scheme (`RTX_FSI2O_BG_CONVECTION`: `tvd`
/// (default, van Albada), `upwind`) — P5-3's wake-dissipation knob on the
/// fixed-motion replay (the patch's scheme was excluded as the h-term's
/// carrier: both ladders move +1.5–1.7 W/m from ny 41 to 62).
pub fn bg_convection() -> ConvectionScheme {
    match std::env::var("RTX_FSI2O_BG_CONVECTION").as_deref() {
        Ok("upwind") => ConvectionScheme::Upwind,
        _ => ConvectionScheme::TvdVanAlbada,
    }
}

/// The patch's convection scheme (`RTX_FSI2O_PATCH_CONVECTION`: `tvd`
/// (default, van Albada), `upwind`, `none`) — P5-3's near-wake dispersion
/// knob on the fixed-motion replay.
pub fn patch_convection() -> PatchConvection {
    match std::env::var("RTX_FSI2O_PATCH_CONVECTION").as_deref() {
        Ok("upwind") => PatchConvection::Upwind,
        Ok("none") => PatchConvection::None,
        _ => PatchConvection::TvdVanAlbada,
    }
}

/// The patch's across-rows (`RTX_FSI2O_PATCH_ROWS`, 12; scale it with the
/// offset to keep the wall spacing).
fn patch_rows() -> usize {
    std::env::var("RTX_FSI2O_PATCH_ROWS")
        .ok()
        .and_then(|v| v.parse().ok())
        .unwrap_or(PATCH_ROWS)
}

/// The save tag: `fsi2o_ny{ny}`, plus `_off{offset}` off the default
/// patch thickness so a thicker patch never loads the default's state.
fn save_tag(ny: usize) -> String {
    let offset = patch_offset_h();
    if (offset - 6.0).abs() < 1e-12 {
        format!("fsi2o_ny{ny}")
    } else {
        format!("fsi2o_ny{ny}_off{offset}")
    }
}

/// The deforming wall as the patch sees it: the wetted polygon (the
/// `Interface` walk, anchors included) and the velocity at each vertex.
#[derive(Debug, Clone, Default)]
pub struct WallMotion {
    pub polygon: Vec<[f64; 2]>,
    pub velocity: Vec<[f64; 2]>,
    /// The net normal velocity removed from every segment (the imposed
    /// velocity's volume flux over the wall length) when the wall is
    /// made volume-preserving; zero otherwise.
    pub q: f64,
}

impl WallMotion {
    /// The imposed velocity's net volume flux INTO the fluid over the
    /// polygon (counter-clockwise around the body: outward normal
    /// `(dy, −dx)/L`), and the polygon's length.
    pub fn net_flux(&self) -> (f64, f64) {
        let (mut flux, mut len) = (0.0, 0.0);
        for i in 0..self.polygon.len().saturating_sub(1) {
            let (a, b) = (self.polygon[i], self.polygon[i + 1]);
            let (dx, dy) = (b[0] - a[0], b[1] - a[1]);
            let (va, vb) = (self.velocity[i], self.velocity[i + 1]);
            let (vx, vy) = (0.5 * (va[0] + vb[0]), 0.5 * (va[1] + vb[1]));
            flux += vx * dy - vy * dx;
            len += (dx * dx + dy * dy).sqrt();
        }
        (flux, len)
    }

    /// The polygon's enclosed area (shoelace; the walk is closed across
    /// the cylinder by its anchors, so this is the flag's area up to a
    /// constant).
    pub fn area(&self) -> f64 {
        let n = self.polygon.len();
        let mut a = 0.0;
        for i in 0..n {
            let (p, q) = (self.polygon[i], self.polygon[(i + 1) % n]);
            a += p[0] * q[1] - q[0] * p[1];
        }
        0.5 * a
    }
}

impl WallMotion {
    /// Velocity at `(x, y)`: linear along the nearest polygon segment.
    pub fn velocity_at(&self, x: f64, y: f64) -> (f64, f64) {
        let n = self.polygon.len();
        if n < 2 {
            return (0.0, 0.0);
        }
        let (mut best_d, mut best) = (f64::INFINITY, (0.0, 0.0));
        for i in 0..n - 1 {
            let (a, b) = (self.polygon[i], self.polygon[i + 1]);
            let (dx, dy) = (b[0] - a[0], b[1] - a[1]);
            let l2 = dx * dx + dy * dy;
            let t = if l2 > 0.0 {
                (((x - a[0]) * dx + (y - a[1]) * dy) / l2).clamp(0.0, 1.0)
            } else {
                0.0
            };
            let (px, py) = (a[0] + t * dx, a[1] + t * dy);
            let d = (px - x).powi(2) + (py - y).powi(2);
            if d < best_d {
                best_d = d;
                let (va, vb) = (self.velocity[i], self.velocity[i + 1]);
                let l = l2.sqrt().max(1e-300);
                best = (
                    va[0] + t * (vb[0] - va[0]) - self.q * dy / l,
                    va[1] + t * (vb[1] - va[1]) + self.q * dx / l,
                );
            }
        }
        best
    }
}

/// The overset fluid of the coupled march.
pub struct OversetFluid {
    pub case: BenchmarkCase,
    pub mesh: Mesh,
    pub interface: Interface,
    pub a_node: NodeId,
    pub ny: usize,
    pub nx: usize,
    pub h: f64,
    pub mu: f64,
    pub dt_fluid: f64,
    pub solver: OversetPisoSolver,
    pub field: OversetField,
    pub shared: Arc<RwLock<WallMotion>>,
    pub sweeps: usize,
    /// Patch regenerations and their wall time.
    pub regen_count: Cell<usize>,
    pub regen_seconds: Cell<f64>,
    /// PERF-2 P0 wall-time buckets [s]: the composite step (`advance`), the
    /// per-pass restore, the per-step snapshot, the load sampling and the
    /// force measurement (`docs/perf2_campaign.md`).
    pub t_advance: Cell<f64>,
    pub t_restore: Cell<f64>,
    pub t_snapshot: Cell<f64>,
    pub t_sample: Cell<f64>,
    pub t_force: Cell<f64>,
    /// Background cells reclassified, summed over every fluid step.
    pub reclassified_total: Cell<usize>,
    pub fresh_total: Cell<usize>,
    pub rounds_total: Cell<usize>,
    pub correctors_total: Cell<usize>,
    /// The interface of the last `set_geometry`.
    pub last_d: Vec<f64>,
    /// The committed step's patch, the warm start of every regeneration
    /// in the next step (`RTX_FSI2O_WARM_SWEEPS`); `None` = cold builds.
    pub warm_base: Option<rtx_cfd::mesh::PatchMesh>,
    pub warm_sweeps: usize,
    /// The last fluid step's mass defects (patch acceptor ring, background
    /// fringe) and the patch's max divergence.
    pub last_defects: Cell<(f64, f64, f64)>,
}

impl OversetFluid {
    /// Build the composite at rest around the undeformed flag.
    pub fn build_case(
        case: BenchmarkCase,
        ny: usize,
        flag_nx: usize,
        sweeps: usize,
        max_rounds: usize,
    ) -> CfdResult<Self> {
        Self::build_case_with(case, ny, flag_nx, sweeps, max_rounds, None)
    }

    /// [`Self::build_case`] with an initial patch mesh given (a saved
    /// instant's deformed patch: the overlap is built around it).
    pub fn build_case_with(
        case: BenchmarkCase,
        ny: usize,
        flag_nx: usize,
        sweeps: usize,
        max_rounds: usize,
        initial: Option<rtx_cfd::mesh::PatchMesh>,
    ) -> CfdResult<Self> {
        let h = H / ny as f64;
        let nx = (L / h).round() as usize;
        let mu = RHO_F * NU_F;
        let u_mean = case.u_mean;
        let u_peak = 1.5 * 1.5 * u_mean;

        let mesh = flag_mesh(flag_nx, 2);
        let interface = Interface::build(&mesh);
        let a_node = mesh
            .nodes
            .iter()
            .find(|(_, n)| {
                (n.position().x - 0.6).abs() < 1e-9 && (n.position().y - 0.2).abs() < 1e-9
            })
            .map(|(&id, _)| id)
            .expect("point A");
        let zero_d = vec![0.0; 2 * interface.wetted.len()];

        let config = CfdConfig::new()
            .with_density(RHO_F)
            .with_viscosity(mu)
            .with_reference_velocity(u_mean)
            .with_reference_length(0.1);
        let mut background = EmbeddedPisoSolver::new(
            config.clone(),
            EmbeddedParameters {
                corrector_steps: 2,
                tolerance: 1e-7,
                boundaries: AleBoundaries {
                    left: SideBoundary::Velocity,
                    right: SideBoundary::PressureOutlet,
                    bottom: SideBoundary::Velocity,
                    top: SideBoundary::Velocity,
                },
                poisson_solver: PoissonSolverKind::Multigrid,
                poisson_precision: MgPrecision::F64,
                // PERF-2: `RTX_FSI2O_MG_RB=1` selects the red-black smoother
                // (a regime change, band-gated; default = the recorded regime).
                poisson_smoother: if std::env::var("RTX_FSI2O_MG_RB").is_ok_and(|v| v == "1") {
                    println!(
                        "  multigrid smoother: RED-BLACK symmetric Gauss–Seidel (PERF-2 regime, RTX_FSI2O_MG_RB)"
                    );
                    rtx_cfd::solvers::incompressible::MgSmoother::RedBlack
                } else {
                    rtx_cfd::solvers::incompressible::MgSmoother::Lexicographic
                },
                convection_scheme: bg_convection(),
            },
        )?;
        // PERF-2 P2: `RTX_THREADS=n` — rayon's global pool and the multigrid's
        // red-black maps on n threads (bit-identical to 1; no effect on the
        // lexicographic smoother).
        let threads: usize = std::env::var("RTX_THREADS")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(1)
            .max(1);
        if threads > 1 {
            rtx_cfd::solvers::incompressible::configure_threads(threads);
            background.set_poisson_threads(threads);
            println!("  multigrid threads: {threads} (RTX_THREADS)");
        }
        background.set_boundary_velocity(move |x, y, t| {
            if x <= 0.0 {
                (inflow_for(u_mean, y, t), 0.0)
            } else {
                (0.0, 0.0)
            }
        });

        let (cold_mesh, _) = cylinder_flag_patch_deformed_tip(
            CYL_CENTRE,
            CYL_R,
            FLAG_T,
            &interface.edges(&zero_d),
            FLAG_X1,
            h,
            fillet(),
            patch_offset_h() * h,
            patch_rows(),
            patch_stretch(),
            sweeps,
            tip_corner(),
        )?;
        // With the warm-started regeneration (`RTX_FSI2O_WARM_SWEEPS`), the
        // march starts on the converged fixed point of the warm build's
        // sweep + re-spacing map, so the chain starts where it stays.
        let warm_sweeps: usize = std::env::var("RTX_FSI2O_WARM_SWEEPS")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(0);
        let start_mesh = if warm_sweeps > 0 {
            let t0 = std::time::Instant::now();
            let (m, (it, last)) = rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed_from_tip(
                Some(&cold_mesh),
                CYL_CENTRE,
                CYL_R,
                FLAG_T,
                &interface.edges(&zero_d),
                FLAG_X1,
                h,
                fillet(),
                patch_offset_h() * h,
                patch_rows(),
                patch_stretch(),
                20_000,
                tip_corner(),
            )?;
            println!(
                "  warm base: {it} sweep+respace iterations (last move {last:.1e}) in {:.1} s; {warm_sweeps} per regeneration",
                t0.elapsed().as_secs_f64()
            );
            m
        } else {
            cold_mesh
        };
        // The fluid step is the march's: from the mesh the march started
        // on, whatever patch this composite is built around (a saved
        // deformed instant) — 0.2 % of dt read as 0.3–5 N/m of solved-face
        // residual in the audit before this.
        let mut hs = f64::INFINITY;
        for c in 0..start_mesh.cell_count() {
            for (f, _) in start_mesh.cell_faces(c) {
                if start_mesh.is_sface(f) {
                    let d = start_mesh.faces()[f].d;
                    hs = hs.min((d[0] * d[0] + d[1] * d[1]).sqrt());
                }
            }
        }
        let dt_bg = 0.25 / (2.0 * u_peak / h + 4.0 * NU_F / (h * h));
        let dt_patch = 0.4 * (hs * hs / (4.0 * NU_F)).min(hs / u_peak);
        // `RTX_FSI2O_DT_SCALE` (1): the fluid step scaled off its CFL-
        // derived value — P5-3's dt-convergence probe of the first-order
        // Euler steps (the s = 8 pair moved the amplitude 95 → 60 mm).
        let dt_scale: f64 = std::env::var("RTX_FSI2O_DT_SCALE")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(1.0);
        let dt_fluid = dt_bg.min(dt_patch) * dt_scale;
        let patch_mesh = initial.unwrap_or(start_mesh);

        let shared = Arc::new(RwLock::new(WallMotion {
            polygon: interface
                .polygon(&zero_d)
                .iter()
                .map(|&(x, y)| [x, y])
                .collect(),
            velocity: vec![[0.0, 0.0]; interface.walk.len()],
            q: 0.0,
        }));
        let wall = shared.clone();
        let mut patch = CurvilinearPisoSolver::new(
            config,
            CurvilinearParameters {
                tolerance: 1e-5,
                convection: patch_convection(),
                normal_diffusion: NormalDiffusion::LineImplicit,
                ..CurvilinearParameters::default()
            },
            patch_mesh,
        )?;
        // Discriminators (P5-1 matrix): `RTX_FSI2O_NO_WALL_VEL` keeps the
        // wall at rest in the fluid while the geometry moves.
        let no_wall_vel = std::env::var("RTX_FSI2O_NO_WALL_VEL").is_ok();
        patch.set_side_velocity(PatchSide::Inner, move |x, y, _| {
            if no_wall_vel {
                (0.0, 0.0)
            } else {
                wall.read().unwrap().velocity_at(x, y)
            }
        });
        let mut patch_field = PatchField::new(patch.mesh());
        patch.initialize(&mut patch_field, |_, _| (0.0, 0.0));

        let params = OversetParameters {
            stall_rounds: 2,
            max_rounds,
            // `RTX_FSI2O_ROWS`: the overlap depth (patch rows kept non-hole
            // below the acceptor row; default 4). A deeper overlap widens
            // the band the acceptors' donors need on a bent patch.
            overlap_rows: std::env::var("RTX_FSI2O_ROWS")
                .ok()
                .and_then(|v| v.parse().ok())
                .unwrap_or(OversetParameters::default().overlap_rows),
            ..OversetParameters::default()
        };
        let mut solver = OversetPisoSolver::new(background, patch, (nx, ny, h, h), params)?;
        // PERF-2 P2-c: `RTX_FSI2O_INNER_STOP=m` — the background Poisson's inner
        // stop at m × the Schwarz tolerance (default 0.1; band-gated ladder).
        let mut solver = solver;
        if let Some(m) = std::env::var("RTX_FSI2O_INNER_STOP")
            .ok()
            .and_then(|v| v.parse::<f64>().ok())
        {
            solver.set_inner_stop_multiplier(m);
            println!(
                "  background Poisson inner stop: {m} × the Schwarz tolerance (RTX_FSI2O_INNER_STOP; default 0.1)"
            );
        }
        let mut field = OversetField {
            background: FlowField::new(nx, ny, h, h)?,
            patch: patch_field,
        };
        solver.initialize(&mut field)?;
        Ok(Self {
            case,
            mesh,
            interface,
            a_node,
            ny,
            nx,
            h,
            mu,
            dt_fluid,
            solver,
            field,
            shared,
            sweeps,
            regen_count: Cell::new(0),
            regen_seconds: Cell::new(0.0),
            t_advance: Cell::new(0.0),
            t_restore: Cell::new(0.0),
            t_snapshot: Cell::new(0.0),
            t_sample: Cell::new(0.0),
            t_force: Cell::new(0.0),
            reclassified_total: Cell::new(0),
            fresh_total: Cell::new(0),
            rounds_total: Cell::new(0),
            correctors_total: Cell::new(0),
            last_d: zero_d,
            last_defects: Cell::new((0.0, 0.0, 0.0)),
            warm_base: None,
            warm_sweeps,
        })
    }

    /// Save the fluid state (`RTX_FSI2O_SAVE=dir`, tag `fsi2o_ny{ny}`):
    /// the background as `FlowField::save`, the patch vectors raw, the
    /// time.
    pub fn save(&self, dir: &str) -> CfdResult<()> {
        let dir = std::path::Path::new(dir);
        std::fs::create_dir_all(dir).expect("save dir");
        let tag = save_tag(self.ny);
        self.field
            .background
            .save(&dir.join(format!("bg_{tag}.bin")))?;
        for (name, vals) in [
            ("u", &self.field.patch.u),
            ("v", &self.field.patch.v),
            ("p", &self.field.patch.p),
            ("flux", &self.field.patch.flux),
        ] {
            let bytes: Vec<u8> = vals.iter().flat_map(|x| x.to_le_bytes()).collect();
            std::fs::write(dir.join(format!("patch_{tag}_{name}.bin")), bytes).expect("save patch");
        }
        std::fs::write(
            dir.join(format!("time_{tag}.txt")),
            format!("{:.17e}", self.solver.time()),
        )
        .expect("save time");
        println!(
            "  saved the fluid state {tag} at t = {:.4} to {}",
            self.solver.time(),
            dir.display()
        );
        Ok(())
    }

    /// Save the composite at an instant of the coupled march for the
    /// offline chain audit: the background, the patch's node coordinates
    /// (the deformed mesh), its vectors, the time and the interface `d`
    /// (`dir/inst_<tag>_<step>/`).
    pub fn save_instant(&self, dir: &str, step: usize, d: &[f64], ddot: &[f64]) -> CfdResult<()> {
        let tag = save_tag(self.ny);
        let dir = std::path::Path::new(dir).join(format!("inst_{tag}_{step:06}"));
        std::fs::create_dir_all(&dir).expect("instant dir");
        self.field.background.save(&dir.join("bg.bin"))?;
        let mesh = self.solver.patch().mesh();
        let (ns, nn) = (mesh.ns(), mesh.nn());
        let mut nodes: Vec<f64> = Vec::with_capacity(2 * (ns + 1) * (nn + 1));
        for k in 0..=nn {
            for i in 0..=ns {
                let p = mesh.node_xy(mesh.node(k, i));
                nodes.push(p[0]);
                nodes.push(p[1]);
            }
        }
        for (name, vals) in [
            ("nodes", &nodes),
            ("u", &self.field.patch.u),
            ("v", &self.field.patch.v),
            ("p", &self.field.patch.p),
            ("flux", &self.field.patch.flux),
            ("d", &d.to_vec()),
            ("dd", &ddot.to_vec()),
        ] {
            let bytes: Vec<u8> = vals.iter().flat_map(|x| x.to_le_bytes()).collect();
            std::fs::write(dir.join(format!("patch_{name}.bin")), bytes).expect("save");
        }
        std::fs::write(
            dir.join("meta.txt"),
            format!(
                "t {:.17e}\nns {ns}\nnn {nn}\ndt {:.17e}\n",
                self.solver.time(),
                self.dt_fluid
            ),
        )
        .expect("meta");
        Ok(())
    }

    /// Rebuild the composite on a saved instant (`save_instant`) — the
    /// deformed patch from its node coordinates, the fields, the wall
    /// motion — and return it with the interface and the time; the chain
    /// then runs on it as on any settled state.
    pub fn from_instant(
        case: BenchmarkCase,
        ny: usize,
        flag_nx: usize,
        dir: &std::path::Path,
    ) -> CfdResult<(Self, Vec<f64>, f64)> {
        let read = |name: &str| -> Vec<f64> {
            let bytes = std::fs::read(dir.join(format!("patch_{name}.bin")))
                .unwrap_or_else(|e| panic!("instant {name}: {e}"));
            bytes
                .chunks_exact(8)
                .map(|c| f64::from_le_bytes(c.try_into().expect("8 bytes")))
                .collect()
        };
        let meta = std::fs::read_to_string(dir.join("meta.txt")).expect("meta");
        let mut t = 0.0;
        let (mut ns, mut nn) = (0usize, 0usize);
        let mut dt_saved: Option<f64> = None;
        for line in meta.lines() {
            let mut it = line.split_whitespace();
            match (it.next(), it.next()) {
                (Some("t"), Some(v)) => t = v.parse().expect("t"),
                (Some("ns"), Some(v)) => ns = v.parse().expect("ns"),
                (Some("nn"), Some(v)) => nn = v.parse().expect("nn"),
                (Some("dt"), Some(v)) => dt_saved = v.parse().ok(),
                _ => {}
            }
        }
        let nodes = read("nodes");
        let (xs, ys): (Vec<f64>, Vec<f64>) = nodes.chunks_exact(2).map(|c| (c[0], c[1])).unzip();
        let mesh = rtx_cfd::mesh::PatchMesh::from_nodes(ns, nn, xs, ys, Some([0.0, 0.0]))?;
        // The composite around the deformed patch: built from scratch so
        // the overlap is the instant's.
        let mut fluid = Self::build_case_with(case, ny, flag_nx, 100, 3, Some(mesh))?;
        let d = read("d");
        let dd = read("dd");
        {
            let mut w = fluid.shared.write().unwrap();
            w.polygon = fluid
                .interface
                .polygon(&d)
                .iter()
                .map(|&(x, y)| [x, y])
                .collect();
            w.velocity = fluid
                .interface
                .walk_velocities(&dd)
                .iter()
                .map(|&(u, v)| [u, v])
                .collect();
        }
        fluid.field.background = FlowField::load(&dir.join("bg.bin"))?;
        fluid.field.patch.u = read("u");
        fluid.field.patch.v = read("v");
        fluid.field.patch.p = read("p");
        fluid.field.patch.flux = read("flux");
        fluid.solver.set_time(t);
        fluid.last_d = d.clone();
        if let Some(dt) = dt_saved {
            fluid.dt_fluid = dt;
        }
        Ok((fluid, d, t))
    }

    /// The wall load by REGION at the current state (P5-3): for the top
    /// face, the bottom face, the tip arc, and the fillets + cylinder —
    /// `(name, faces, length, drag, lift, mean t_n, min t_n, max t_n)`
    /// with `t_n` the normal traction (≈ −p at the wall), and the patch's
    /// pressure level.
    pub fn wall_regions(
        &self,
        d: &[f64],
    ) -> (
        Vec<(&'static str, usize, f64, f64, f64, f64, f64, f64)>,
        f64,
    ) {
        let e = self.interface.edges(d);
        let tip_mid = [
            0.5 * (e.tip[0][0] + e.tip[e.tip.len() - 1][0]),
            0.5 * (e.tip[0][1] + e.tip[e.tip.len() - 1][1]),
        ];
        let nearest = |pts: &[[f64; 2]], q: [f64; 2]| -> f64 {
            pts.iter()
                .map(|p| (p[0] - q[0]).powi(2) + (p[1] - q[1]).powi(2))
                .fold(f64::INFINITY, f64::min)
        };
        let mut acc: Vec<(&'static str, usize, f64, f64, f64, f64, f64, f64)> = vec![
            (
                "top",
                0,
                0.0,
                0.0,
                0.0,
                0.0,
                f64::INFINITY,
                f64::NEG_INFINITY,
            ),
            (
                "bottom",
                0,
                0.0,
                0.0,
                0.0,
                0.0,
                f64::INFINITY,
                f64::NEG_INFINITY,
            ),
            (
                "tip arc",
                0,
                0.0,
                0.0,
                0.0,
                0.0,
                f64::INFINITY,
                f64::NEG_INFINITY,
            ),
            (
                "fillets+cyl",
                0,
                0.0,
                0.0,
                0.0,
                0.0,
                f64::INFINITY,
                f64::NEG_INFINITY,
            ),
        ];
        for (centre, normal, len, traction) in self.solver.patch().wall_tractions(
            &self.field.patch,
            PatchSide::Inner,
            self.solver.time(),
        ) {
            let on_cyl =
                ((centre[0] - CYL_CENTRE[0]).powi(2) + (centre[1] - CYL_CENTRE[1]).powi(2)).sqrt()
                    < CYL_R + 1.5 * fillet();
            let near_tip = ((centre[0] - tip_mid[0]).powi(2) + (centre[1] - tip_mid[1]).powi(2))
                .sqrt()
                < 2.5 * FLAG_T;
            let k = if on_cyl {
                3
            } else if near_tip {
                2
            } else if nearest(&e.top, centre) <= nearest(&e.bottom, centre) {
                0
            } else {
                1
            };
            let tn = traction[0] * normal[0] + traction[1] * normal[1];
            let r = &mut acc[k];
            r.1 += 1;
            r.2 += len;
            r.3 += traction[0] * len;
            r.4 += traction[1] * len;
            r.5 += tn * len;
            r.6 = r.6.min(tn);
            r.7 = r.7.max(tn);
        }
        for r in &mut acc {
            if r.2 > 0.0 {
                r.5 /= r.2;
            }
        }
        let p = &self.field.patch.p;
        let level = p.iter().sum::<f64>() / p.len().max(1) as f64;
        (acc, level)
    }

    /// The solver-metric momentum chain at the current state (§5.11's
    /// instrument on the moving patch): the solved-face pin, the ring, the
    /// box in the solver's flux form, the patch's balance, the wall.
    pub fn chain_line(&self) -> String {
        let dt = self.dt_fluid;
        let h = self.h;
        let (nx, ny) = (self.nx, self.ny);
        let cv = (
            (0.10 / h).round() as usize,
            (0.75 / h).round() as usize,
            (0.05 / h).round() as usize,
            (0.36 / h).round() as usize,
        );
        let _ = (nx, ny);
        let wall = self.measure_force();
        let mr = self.solver.momentum_residual(&self.field, dt);
        let ring = mr.fringe_fringe.fx + mr.fringe_hole.fx;
        let bx = self.solver.solver_metric_force(&self.field, dt, cv).0;
        let bx2 = self
            .solver
            .solver_metric_force(
                &self.field,
                dt,
                (
                    (0.09 / h).round() as usize,
                    (0.70 / h).round() as usize,
                    (0.07 / h).round() as usize,
                    (0.34 / h).round() as usize,
                ),
            )
            .0;
        let pb = self
            .solver
            .patch()
            .momentum_balance(&self.field.patch, self.solver.time());
        let ff = pb.flux_force()[0];
        let fw = pb.wall_force()[0];
        format!(
            "solved far Σ|r| ({:.2e}, {:.2e}) {}/{} | near ring Σr ({:+.2e}, {:+.2e}) | box (solver flux form) {:.3} [2 boxes spread {:.1e}] → ring Σr {:+.3} ({} + {} faces of {} + {}) → hole flux {:.3} → band {:+.3} → patch interface {:.3} → interior {:+.3} (δP {:+.3}, balance residual {:+.3}) → wall, scheme fluxes {:.3} → wall formula {:+.3} → wall ({:.3}, {:.3}); total wall − box {:+.3} ({:+.2} %); reclassified so far {}",
            mr.solved_far.abs_x,
            mr.solved_far.abs_y,
            mr.solved_far.evaluated,
            mr.solved_far.total,
            mr.solved_near.fx,
            mr.solved_near.fy,
            bx,
            (bx - bx2).abs(),
            ring,
            mr.fringe_fringe.evaluated,
            mr.fringe_hole.evaluated,
            mr.fringe_fringe.total,
            mr.fringe_hole.total,
            bx + ring,
            ff - (bx + ring),
            ff,
            fw - ff,
            pb.pressure_defect()[0],
            pb.balance()[0],
            fw,
            wall.0 - fw,
            wall.0,
            wall.1,
            wall.0 - bx,
            100.0 * (wall.0 - bx) / wall.0,
            self.reclassified_total.get()
        )
    }

    /// Load a state saved by [`Self::save`] (same ny, undeformed patch);
    /// returns its time.
    pub fn load(&mut self, dir: &str) -> CfdResult<f64> {
        let dir = std::path::Path::new(dir);
        let tag = save_tag(self.ny);
        self.field.background = FlowField::load(&dir.join(format!("bg_{tag}.bin")))?;
        let read = |name: &str| -> Vec<f64> {
            let bytes = std::fs::read(dir.join(format!("patch_{tag}_{name}.bin")))
                .unwrap_or_else(|e| panic!("load patch {name}: {e}"));
            bytes
                .chunks_exact(8)
                .map(|c| f64::from_le_bytes(c.try_into().expect("8 bytes")))
                .collect()
        };
        self.field.patch.u = read("u");
        self.field.patch.v = read("v");
        self.field.patch.p = read("p");
        self.field.patch.flux = read("flux");
        let t: f64 = std::fs::read_to_string(dir.join(format!("time_{tag}.txt")))
            .expect("time")
            .trim()
            .parse()
            .expect("time value");
        self.solver.set_time(t);
        println!(
            "  loaded the fluid state {tag} at t = {t:.4} from {}",
            dir.display()
        );
        Ok(t)
    }

    /// The patch around the interface `d`.
    pub fn patch_for(&self, d: &[f64]) -> CfdResult<rtx_cfd::mesh::PatchMesh> {
        let start = std::time::Instant::now();
        let (mesh, _) = rtx_cfd::mesh::patch_gen::cylinder_flag_patch_deformed_from_tip(
            if self.warm_sweeps > 0 {
                self.warm_base.as_ref()
            } else {
                None
            },
            CYL_CENTRE,
            CYL_R,
            FLAG_T,
            &self.interface.edges(d),
            FLAG_X1,
            self.h,
            fillet(),
            patch_offset_h() * self.h,
            patch_rows(),
            patch_stretch(),
            if self.warm_sweeps > 0 && self.warm_base.is_some() {
                self.warm_sweeps
            } else {
                self.sweeps
            },
            tip_corner(),
        )?;
        self.regen_count.set(self.regen_count.get() + 1);
        self.regen_seconds
            .set(self.regen_seconds.get() + start.elapsed().as_secs_f64());
        Ok(mesh)
    }

    /// The wall for the next fluid step: geometry `d`, velocity `ddot`.
    pub fn set_geometry(&mut self, d: &[f64], ddot: &[f64]) -> CfdResult<()> {
        {
            let mut w = self.shared.write().unwrap();
            w.polygon = self
                .interface
                .polygon(d)
                .iter()
                .map(|&(x, y)| [x, y])
                .collect();
            w.velocity = self
                .interface
                .walk_velocities(ddot)
                .iter()
                .map(|&(u, v)| [u, v])
                .collect();
            // `RTX_FSI2O_VOLUME_PRESERVE`: the fluid sees a volume-preserving
            // wall — the imposed velocity's net flux (the flag's thickness
            // breathing under the pressure step, a ~300 Hz mode of the
            // ν = 0.4 solid) is removed uniformly along the wall.
            w.q = 0.0;
            if std::env::var("RTX_FSI2O_VOLUME_PRESERVE").is_ok() {
                let (flux, len) = w.net_flux();
                w.q = flux / len.max(1e-300);
            }
        }
        // `RTX_FSI2O_FREEZE_PATCH`: the undeformed patch throughout (the
        // wall velocity still moves) — is the regeneration the driver?
        if std::env::var("RTX_FSI2O_FREEZE_PATCH").is_ok() {
            self.last_d = d.to_vec();
            return Ok(());
        }
        let mesh = match self.patch_for(d) {
            Ok(m) => m,
            Err(e) => {
                self.dump_edges(d, &format!("generator refused ({e:?})"));
                return Err(e);
            }
        };
        self.last_d = d.to_vec();
        self.solver.set_patch_mesh(mesh)
    }

    /// The last geometry, for the offline reproduction
    /// (`patch_cylinder_flag_deformed.rs`, `RTX_CF_EDGES_FILE`): one edge
    /// per block, `x y` per line, when `RTX_FSI2O_DUMP_DIR` is set.
    fn dump_edges(&self, d: &[f64], why: &str) {
        let Ok(dir) = std::env::var("RTX_FSI2O_DUMP_DIR") else {
            return;
        };
        let e = self.interface.edges(d);
        let mut out = String::new();
        for (name, pts) in [("bottom", &e.bottom), ("tip", &e.tip), ("top", &e.top)] {
            out.push_str(&format!("# {name} {}\n", pts.len()));
            for p in pts {
                out.push_str(&format!("{:.17e} {:.17e}\n", p[0], p[1]));
            }
        }
        let path = std::path::Path::new(&dir).join("p5_death_edges.txt");
        std::fs::write(&path, out).expect("dump edges");
        println!("  {why}; edges dumped to {}", path.display());
    }

    /// One fluid step at the current wall.
    pub fn step(&mut self) -> CfdResult<OversetResult> {
        let t0 = std::time::Instant::now();
        let r = match futures::executor::block_on(
            self.solver.advance(&mut self.field, self.dt_fluid),
        ) {
            Ok(r) => r,
            Err(e) => {
                // The overlap is rebuilt inside `advance`: a refused patch
                // surfaces here, with the last geometry.
                let d = self.last_d.clone();
                self.dump_edges(&d, &format!("advance refused ({e:?})"));
                return Err(e);
            }
        };
        self.last_defects.set((
            r.patch_mass_defect,
            r.background_mass_defect,
            r.patch_max_divergence,
        ));
        self.reclassified_total
            .set(self.reclassified_total.get() + r.reclassified_cells);
        self.fresh_total.set(self.fresh_total.get() + r.fresh_cells);
        self.rounds_total
            .set(self.rounds_total.get() + r.rounds.iter().sum::<usize>());
        self.correctors_total
            .set(self.correctors_total.get() + r.rounds.len());
        self.t_advance
            .set(self.t_advance.get() + t0.elapsed().as_secs_f64());
        Ok(r)
    }

    /// `subcycle` fluid substeps from the current state with the interface
    /// interpolated from `d_n` to `d_candidate` (the harness's
    /// `advance_subcycled`, digit for digit in the kinematics).
    pub fn advance_subcycled(
        &mut self,
        d_n: &[f64],
        d_candidate: &[f64],
        subcycle: usize,
        v_n: Option<&[f64]>,
    ) -> CfdResult<()> {
        let dt = self.dt_fluid * subcycle as f64;
        let mean_velocity: Vec<f64> = d_candidate
            .iter()
            .zip(d_n)
            .map(|(new, old)| (new - old) / dt)
            .collect();
        for m in 1..=subcycle {
            let fraction = m as f64 / subcycle as f64;
            let (d_sub, ddot_sub): (Vec<f64>, Vec<f64>) = match v_n {
                None => (
                    d_n.iter()
                        .zip(d_candidate)
                        .map(|(old, new)| old + fraction * (new - old))
                        .collect(),
                    mean_velocity.clone(),
                ),
                Some(v_start) => {
                    let mut d_sub = Vec::with_capacity(d_n.len());
                    let mut ddot_sub = Vec::with_capacity(d_n.len());
                    for k in 0..d_n.len() {
                        let v_end = 2.0 * mean_velocity[k] - v_start[k];
                        let accel = (v_end - v_start[k]) / dt;
                        let tau = fraction * dt;
                        d_sub.push(d_n[k] + v_start[k] * tau + 0.5 * accel * tau * tau);
                        ddot_sub.push(v_start[k] + accel * tau);
                    }
                    (d_sub, ddot_sub)
                }
            };
            self.set_geometry(&d_sub, &ddot_sub)?;
            self.step()?;
        }
        Ok(())
    }

    /// The fluid's traction on every Inner (wall) face, in the patch's
    /// Inner-face order — the Robin wall's datum for the next pass.
    pub fn inner_tractions(&self) -> Vec<[f64; 2]> {
        self.solver
            .patch()
            .wall_tractions(&self.field.patch, PatchSide::Inner, self.solver.time())
            .into_iter()
            .map(|(_, _, _, t)| t)
            .collect()
    }
    /// Put the Robin wall (impedance `alpha`, datum per Inner face) on the
    /// patch, or remove it.
    pub fn set_robin(&mut self, alpha: f64, datum: Option<Vec<[f64; 2]>>) {
        self.solver
            .patch_mut()
            .set_robin_wall(datum.map(|d| RobinWall { alpha, datum: d }));
    }

    /// Drag and lift on cylinder + flag from the patch's wall stress.
    pub fn measure_force(&self) -> (f64, f64) {
        let t0 = std::time::Instant::now();
        let f = self
            .solver
            .patch()
            .surface_force(&self.field.patch, PatchSide::Inner, self.solver.time())
            .total();
        self.t_force
            .set(self.t_force.get() + t0.elapsed().as_secs_f64());
        (f[0], f[1])
    }

    /// The wall faces' tractions transferred to the flag's wetted nodes
    /// at geometry `d`: `(nodal forces, conservation defect, faces used)`.
    /// Faces on the cylinder proper are skipped; the fillets' load goes
    /// to the nearest (clamped) root nodes.
    pub fn sample_load(&self, d: &[f64]) -> (Vec<(NodeId, Vector3<f64>)>, f64, usize) {
        let t0 = std::time::Instant::now();
        let out = self.sample_load_inner(d);
        self.t_sample
            .set(self.t_sample.get() + t0.elapsed().as_secs_f64());
        out
    }

    fn sample_load_inner(&self, d: &[f64]) -> (Vec<(NodeId, Vector3<f64>)>, f64, usize) {
        let mut faces = Vec::new();
        let mut tractions: Vec<Vector3<f64>> = Vec::new();
        for (centre, normal, len, traction) in self.solver.patch().wall_tractions(
            &self.field.patch,
            PatchSide::Inner,
            self.solver.time(),
        ) {
            let on_cylinder =
                ((centre[0] - CYL_CENTRE[0]).powi(2) + (centre[1] - CYL_CENTRE[1]).powi(2)).sqrt()
                    < CYL_R + 1e-9;
            if on_cylinder {
                continue;
            }
            faces.push(FluidFace {
                centroid: Vector3::new(centre[0], centre[1], 0.0),
                normal: Vector3::new(normal[0], normal[1], 0.0),
                area: len,
            });
            tractions.push(Vector3::new(traction[0], traction[1], 0.0));
        }
        let nodes_now = self.interface.deformed_nodes(d);
        let surface = WettedSurface::build(&faces, &nodes_now).expect("transfer build");
        let nodal = surface.transfer_load(&faces, &tractions).unwrap();
        let total_sampled: Vector3<f64> =
            faces.iter().zip(&tractions).map(|(f, t)| t * f.area).sum();
        let total_nodal: Vector3<f64> = nodal.iter().sum();
        let conservation = (total_nodal - total_sampled).norm() / total_sampled.norm().max(1e-30);
        (
            self.interface
                .wetted
                .iter()
                .zip(nodal)
                .map(|(&id, f)| (id, f))
                .collect(),
            conservation,
            faces.len(),
        )
    }

    /// The wall pressure's roughness along the flag: over the flag's
    /// wall faces in order, `Σ|t_n(f+1) − t_n(f)| / Σ|t_n(f)|` with `t_n`
    /// the normal traction, and the count of sign changes of the
    /// face-to-face difference — a checkerboard reads ≫ 1 with a sign
    /// change at every face.
    pub fn wall_roughness(&self) -> (f64, usize, f64, f64) {
        let mut tn: Vec<f64> = Vec::new();
        for (centre, normal, _, traction) in self.solver.patch().wall_tractions(
            &self.field.patch,
            PatchSide::Inner,
            self.solver.time(),
        ) {
            let on_cylinder =
                ((centre[0] - CYL_CENTRE[0]).powi(2) + (centre[1] - CYL_CENTRE[1]).powi(2)).sqrt()
                    < CYL_R + 1e-9;
            if !on_cylinder {
                tn.push(traction[0] * normal[0] + traction[1] * normal[1]);
            }
        }
        let total: f64 = tn.iter().map(|v| v.abs()).sum();
        let mut jumps = 0.0;
        let mut flips = 0usize;
        let mut prev_diff = 0.0;
        for w in tn.windows(2) {
            let d = w[1] - w[0];
            jumps += d.abs();
            if prev_diff * d < 0.0 {
                flips += 1;
            }
            prev_diff = d;
        }
        let max = tn.iter().fold(0.0_f64, |m, v| m.max(v.abs()));
        (
            jumps / total.max(1e-300),
            flips,
            max,
            total / tn.len().max(1) as f64,
        )
    }

    /// The patch's pressure level (mean over its cells) and the wall's
    /// net volume flux from the imposed velocity (`Σ u_wall · S_out` over
    /// the Inner faces, positive = fluid leaving the patch through the
    /// wall), against the inner ring's enclosed-area rate would be the
    /// mesh's own; both must match for a consistent moving wall.
    pub fn level_and_wall_flux(&self) -> (f64, f64, f64, f64, f64) {
        let mesh = self.solver.patch().mesh();
        let p = &self.field.patch.p;
        let level = p.iter().sum::<f64>() / p.len().max(1) as f64;
        let wall = self.shared.read().unwrap();
        let (poly_flux, _) = wall.net_flux();
        let poly_area = wall.area();
        let (mut flux, mut area) = (0.0, 0.0);
        for (f, face) in mesh.faces().iter().enumerate() {
            if mesh.side(f) != Some(PatchSide::Inner) {
                continue;
            }
            // Inner side: S points from the body into the fluid (the
            // neighbour is the cell), so −S is out of the fluid.
            let sign = if face.neigh.is_some() { 1.0 } else { -1.0 };
            let s_in = [sign * face.s[0], sign * face.s[1]];
            let (u, v) = wall.velocity_at(face.centre[0], face.centre[1]);
            // Volume flux INTO the fluid = u_wall · S_into_fluid.
            flux += u * s_in[0] + v * s_in[1];
            area += (s_in[0] * s_in[0] + s_in[1] * s_in[1]).sqrt();
        }
        (level, flux, area, poly_flux, poly_area)
    }

    /// The current patch becomes the warm start of the next step's
    /// regenerations (called after a committed step).
    pub fn commit_base(&mut self) {
        if self.warm_sweeps > 0 {
            self.warm_base = Some(self.solver.patch().mesh().clone());
        }
    }

    pub fn snapshot(&self) -> (OversetSolverState, OversetField) {
        let t0 = std::time::Instant::now();
        let out = (self.solver.snapshot(), self.field.clone());
        self.t_snapshot
            .set(self.t_snapshot.get() + t0.elapsed().as_secs_f64());
        out
    }

    pub fn restore(&mut self, saved: &(OversetSolverState, OversetField)) {
        let t0 = std::time::Instant::now();
        self.solver.restore(&saved.0);
        self.field = saved.1.clone();
        self.t_restore
            .set(self.t_restore.get() + t0.elapsed().as_secs_f64());
    }

    /// PERF-2 P0: the solver's own split of `advance` (when `RTX_PROFILE`
    /// is set) as one printable line, with the background Poisson's
    /// setup / iterate split.
    pub fn profile_line(&self) -> Option<String> {
        let t = self.solver.timers()?;
        let (ps, pi, pc, pk) = self.solver.background().poisson_profile();
        let rg = rtx_cfd::mesh::patch_gen::regen_profile();
        let s = |ns: u64| ns as f64 * 1e-9;
        let adv = s(t.advance_ns).max(1e-300);
        let pct = |ns: u64| 100.0 * s(ns) / adv;
        Some(
            format!(
                "  advance split over {} steps ({} rounds), {:.0} s: overlap build {:.0} s ({:.1}%), predictor bg {:.0} s ({:.1}%), predictor patch {:.0} s ({:.1}%), bg Poisson {:.0} s ({:.1}%) [setup {:.0} s, iterate {:.0} s, {} solves, {} CG iterations, {:.2} ms per V-cycle], patch BiCGSTAB {:.0} s ({:.1}%), round exchange {:.0} s ({:.1}%), apply {:.0} s ({:.1}%), end exchange {:.0} s ({:.1}%), other {:.0} s",
                t.steps,
                t.rounds,
                adv,
                s(t.overlap_build_ns),
                pct(t.overlap_build_ns),
                s(t.predictor_bg_ns),
                pct(t.predictor_bg_ns),
                s(t.predictor_patch_ns),
                pct(t.predictor_patch_ns),
                s(t.bg_solve_ns),
                pct(t.bg_solve_ns),
                s(ps),
                s(pi),
                pc,
                pk,
                1e3 * s(pi) / (pk + pc).max(1) as f64,
                s(t.patch_solve_ns),
                pct(t.patch_solve_ns),
                s(t.round_exchange_ns),
                pct(t.round_exchange_ns),
                s(t.apply_ns),
                pct(t.apply_ns),
                s(t.end_exchange_ns),
                pct(t.end_exchange_ns),
                adv - s(t.overlap_build_ns
                    + t.predictor_bg_ns
                    + t.predictor_patch_ns
                    + t.bg_solve_ns
                    + t.patch_solve_ns
                    + t.round_exchange_ns
                    + t.apply_ns
                    + t.end_exchange_ns),
            ) + &format!(
                "\n  regeneration split ({} builds): outline {:.0} s, hull + offset {:.0} s, ring projection {:.0} s, Winslow sweeps {:.0} s, respace {:.0} s, warm bookkeeping {:.0} s, mesh finalisation {:.0} s",
                rg[7],
                s(rg[0]),
                s(rg[1]),
                s(rg[2]),
                s(rg[3]),
                s(rg[4]),
                s(rg[5]),
                s(rg[6])
            ),
        )
    }

    pub fn time(&self) -> f64 {
        self.solver.time()
    }
}