Adom Fields

//! Steady 3D solid conduction. Prescribed convection is NOT a CFD fan model.
use crate::solve::System;
use serde::{Deserialize, Serialize};
use serde_json::{json, Value};
use sha2::{Digest, Sha256};
use std::{
    collections::{HashMap, HashSet},
    path::Path,
};

#[derive(Clone, Deserialize, Serialize)]
#[serde(rename_all = "camelCase", deny_unknown_fields)]
pub struct Body {
    pub name: String,
    /// Union of axis-aligned boxes [xmin,ymin,zmin,xmax,ymax,zmax], mm.
    pub boxes_mm: Vec<[f64; 6]>,
    pub conductivity_w_mk: f64,
    /// Exposed-face h, ordered -x,+x,-y,+y,-z,+z. Zero means insulated.
    pub convection_w_m2k: [f64; 6],
    pub watts: f64,
    pub provenance: String,
}
#[derive(Clone, Deserialize, Serialize)]
#[serde(rename_all = "camelCase", deny_unknown_fields)]
pub struct Contact {
    pub bodies: [String; 2],
    pub resistance_m2k_w: f64,
    pub provenance: String,
}
#[derive(Clone, Deserialize, Serialize)]
#[serde(rename_all = "camelCase", deny_unknown_fields)]
pub struct Spec {
    pub name: String,
    pub ambient_c: f64,
    pub pitch_mm: f64,
    pub bodies: Vec<Body>,
    pub contacts: Vec<Contact>,
    pub assumptions: Vec<String>,
}
#[derive(Clone, Serialize)]
#[serde(rename_all = "camelCase")]
struct Cell {
    body: usize,
    bounds_mm: [f64; 6],
    temperature_c: f64,
}
fn positive(x: f64) -> bool {
    x.is_finite() && x > 0.0
}
fn nonnegative(x: f64) -> bool {
    x.is_finite() && x >= 0.0
}
fn axis(mut edges: Vec<f64>, pitch: f64) -> Result<Vec<f64>, String> {
    edges.sort_by(f64::total_cmp);
    edges.dedup_by(|a, b| (*a - *b).abs() < 1e-8);
    let mut out = Vec::new();
    for p in edges.windows(2) {
        let n = ((p[1] - p[0]) / pitch).ceil() as usize;
        if n > 10000 || out.len() + n > 10000 {
            return Err("axis exceeds 10000 cells".into());
        }
        for i in 0..n {
            out.push(p[0] + (p[1] - p[0]) * i as f64 / n as f64);
        }
    }
    out.push(*edges.last().ok_or("empty geometry")?);
    Ok(out)
}
fn near_index(a: &[f64], x: f64) -> usize {
    a.iter().position(|v| (*v - x).abs() < 1e-7).unwrap()
}

pub fn solve(spec: &Spec) -> Result<Value, String> {
    if !spec.ambient_c.is_finite()
        || spec.ambient_c <= -273.15
        || !positive(spec.pitch_mm)
        || spec.bodies.is_empty()
        || spec.assumptions.is_empty()
    {
        return Err(
            "require ambient above absolute zero, positive pitch, bodies and assumptions".into(),
        );
    }
    let mut names = HashMap::new();
    let mut axes: [Vec<f64>; 3] = Default::default();
    for (i, b) in spec.bodies.iter().enumerate() {
        if b.name.trim().is_empty()
            || names.insert(b.name.clone(), i).is_some()
            || b.provenance.trim().is_empty()
            || b.boxes_mm.is_empty()
            || !positive(b.conductivity_w_mk)
            || !nonnegative(b.watts)
            || b.convection_w_m2k.iter().any(|h| !nonnegative(*h))
        {
            return Err("invalid or duplicate body; require material/geometry provenance".into());
        }
        for bb in &b.boxes_mm {
            for d in 0..3 {
                if !bb[d].is_finite() || !bb[d + 3].is_finite() || bb[d + 3] - bb[d] < 1e-7 {
                    return Err("invalid box bounds".into());
                }
                axes[d].extend([bb[d], bb[d + 3]]);
            }
        }
    }
    let mut contacts = HashMap::new();
    for c in &spec.contacts {
        let a = *names.get(&c.bodies[0]).ok_or("unknown contact body")?;
        let b = *names.get(&c.bodies[1]).ok_or("unknown contact body")?;
        if a == b
            || !nonnegative(c.resistance_m2k_w)
            || c.provenance.trim().is_empty()
            || contacts
                .insert((a.min(b), a.max(b)), c.resistance_m2k_w)
                .is_some()
        {
            return Err("invalid/duplicate contact".into());
        }
    }
    for a in &mut axes {
        *a = axis(std::mem::take(a), spec.pitch_mm)?;
    }
    let dims = [axes[0].len() - 1, axes[1].len() - 1, axes[2].len() - 1];
    let total = dims
        .iter()
        .try_fold(1usize, |a, b| a.checked_mul(*b))
        .ok_or("grid overflow")?;
    if total > 2_000_000 {
        return Err("grid exceeds 2 million cells; simplify geometry or increase pitch".into());
    }
    let flat = |p: [usize; 3]| p[0] + dims[0] * (p[1] + dims[1] * p[2]);
    let mut owner = vec![usize::MAX; total];
    for (bi, b) in spec.bodies.iter().enumerate() {
        for bb in &b.boxes_mm {
            let lo = [
                near_index(&axes[0], bb[0]),
                near_index(&axes[1], bb[1]),
                near_index(&axes[2], bb[2]),
            ];
            let hi = [
                near_index(&axes[0], bb[3]),
                near_index(&axes[1], bb[4]),
                near_index(&axes[2], bb[5]),
            ];
            for z in lo[2]..hi[2] {
                for y in lo[1]..hi[1] {
                    for x in lo[0]..hi[0] {
                        let o = &mut owner[flat([x, y, z])];
                        if *o != usize::MAX && *o != bi {
                            return Err(format!(
                                "overlapping bodies {} and {}",
                                spec.bodies[*o].name, b.name
                            ));
                        }
                        *o = bi;
                    }
                }
            }
        }
    }
    let mut map = vec![usize::MAX; total];
    let mut cells = Vec::new();
    let mut coords = Vec::new();
    let mut volumes = vec![0.0; spec.bodies.len()];
    let mut sys = System::new(0, 1, 1);
    for z in 0..dims[2] {
        for y in 0..dims[1] {
            for x in 0..dims[0] {
                let p = [x, y, z];
                let o = owner[flat(p)];
                if o == usize::MAX {
                    continue;
                }
                let bb = [
                    axes[0][x],
                    axes[1][y],
                    axes[2][z],
                    axes[0][x + 1],
                    axes[1][y + 1],
                    axes[2][z + 1],
                ];
                volumes[o] += (bb[3] - bb[0]) * (bb[4] - bb[1]) * (bb[5] - bb[2]);
                map[flat(p)] = sys.add_node(0.0);
                coords.push(p);
                cells.push(Cell {
                    body: o,
                    bounds_mm: bb,
                    temperature_c: 0.0,
                });
            }
        }
    }
    let mut used = HashSet::new();
    let mut exposed = Vec::new();
    let mut edge_meta = Vec::new();
    for (i, c) in cells.iter().enumerate() {
        let p = coords[i];
        let size = [
            c.bounds_mm[3] - c.bounds_mm[0],
            c.bounds_mm[4] - c.bounds_mm[1],
            c.bounds_mm[5] - c.bounds_mm[2],
        ];
        for d in 0..3 {
            for sign in 0..2 {
                let mut nb = p;
                let j = if (sign == 0 && p[d] == 0) || (sign == 1 && p[d] + 1 == dims[d]) {
                    usize::MAX
                } else {
                    if sign == 0 {
                        nb[d] -= 1;
                    } else {
                        nb[d] += 1;
                    }
                    map[flat(nb)]
                };
                let area = size[(d + 1) % 3] * size[(d + 2) % 3] * 1e-6;
                if j == usize::MAX {
                    let h = spec.bodies[c.body].convection_w_m2k[2 * d + sign];
                    // Includes half-cell conduction to the convecting face.
                    let g = if h > 0.0 {
                        area / (1.0 / h + size[d] * 0.0005 / spec.bodies[c.body].conductivity_w_mk)
                    } else {
                        0.0
                    };
                    sys.extra_leak[i] += g;
                    exposed.push((i, d, sign, area, g));
                } else if sign == 1 {
                    let other = &cells[j];
                    let key = (c.body.min(other.body), c.body.max(other.body));
                    let rc = if c.body == other.body {
                        0.0
                    } else {
                        used.insert(key);
                        *contacts.get(&key).ok_or_else(|| {
                            format!(
                                "touching bodies require explicit contact: {} / {}",
                                spec.bodies[c.body].name, spec.bodies[other.body].name
                            )
                        })?
                    };
                    let g = area
                        / (size[d] * 0.0005 / spec.bodies[c.body].conductivity_w_mk
                            + (other.bounds_mm[d + 3] - other.bounds_mm[d]) * 0.0005
                                / spec.bodies[other.body].conductivity_w_mk
                            + rc);
                    sys.edges.push((i, j, g));
                    edge_meta.push((d, area));
                }
            }
        }
    }
    if contacts.keys().any(|k| !used.contains(k)) {
        return Err("declared contact has no shared face; check placement".into());
    }
    // Every connected solid island must have a path to ambient, even unpowered islands.
    let mut adj = vec![Vec::new(); cells.len()];
    for &(a, b, _) in &sys.edges {
        adj[a].push(b);
        adj[b].push(a);
    }
    let mut seen = vec![false; cells.len()];
    for root in 0..cells.len() {
        if !seen[root] {
            let mut stack = vec![root];
            seen[root] = true;
            let mut sink = false;
            while let Some(a) = stack.pop() {
                sink |= sys.extra_leak[a] > 0.0;
                for &b in &adj[a] {
                    if !seen[b] {
                        seen[b] = true;
                        stack.push(b);
                    }
                }
            }
            if !sink {
                return Err("solid island has no thermal path to ambient".into());
            }
        }
    }
    let rhs: Vec<f64> = cells
        .iter()
        .map(|c| {
            let b = c.bounds_mm;
            spec.bodies[c.body].watts * (b[3] - b[0]) * (b[4] - b[1]) * (b[5] - b[2])
                / volumes[c.body]
        })
        .collect();
    let (rise, iterations, residual) = sys.solve(&rhs, 1e-10, 20000);
    let input: f64 = rhs.iter().sum();
    let rejected: f64 = rise.iter().zip(&sys.extra_leak).map(|(t, g)| t * g).sum();
    let balance = (input - rejected).abs() / input.max(1e-12);
    if !residual.is_finite()
        || residual > 1e-8
        || rise.iter().any(|v| !v.is_finite() || *v < -1e-6)
        || balance > 1e-6
    {
        return Err(format!(
            "failed convergence/energy balance: residual {residual}, balance {balance}"
        ));
    }
    for (c, t) in cells.iter_mut().zip(&rise) {
        c.temperature_c = spec.ambient_c + t;
    }
    let body_results:Vec<Value>=spec.bodies.iter().enumerate().map(|(i,b)|{
        let cs:Vec<&Cell>=cells.iter().filter(|c|c.body==i).collect();
        json!({"name":b.name,"watts":b.watts,"minC":cs.iter().map(|c|c.temperature_c).fold(f64::INFINITY,f64::min),"maxC":cs.iter().map(|c|c.temperature_c).fold(f64::NEG_INFINITY,f64::max),"volumeMm3":volumes[i]})
    }).collect();
    let faces:Vec<Value>=exposed.iter().map(|&(i,d,s,area,g)|json!({"cell":i,"axis":d,"positive":s==1,"areaM2":area,"wattsToAmbient":g*rise[i]})).collect();
    let flux:Vec<Value>=sys.edges.iter().zip(&edge_meta).map(|(&(a,b,g),&(d,area))|json!({"from":a,"to":b,"axis":d,"watts":g*(rise[a]-rise[b]),"fluxWm2":g*(rise[a]-rise[b])/area})).collect();
    Ok(
        json!({"schema":"adom-fields-thermal3d-v1","name":spec.name,"method":"Steady finite-volume solid conduction; prescribed convection; no CFD, radiation or heat-pipe transport","ambientC":spec.ambient_c,"pitchMm":spec.pitch_mm,"assumptions":spec.assumptions,"input":spec,"gridDimensions":dims,"cellCount":cells.len(),"iterations":iterations,"relativeResidual":residual,"inputWatts":input,"rejectedWatts":rejected,"relativeEnergyError":balance,"bodies":body_results,"cells":cells,"exposedFaces":faces,"conductiveFaces":flux}),
    )
}

pub fn run(path: &str, out: &Path) -> Result<(), String> {
    let bytes = std::fs::read(path).map_err(|e| e.to_string())?;
    let spec: Spec = serde_json::from_slice(&bytes).map_err(|e| e.to_string())?;
    let mut result = solve(&spec)?;
    result["inputSha256"] = json!(format!("{:x}", Sha256::digest(&bytes)));
    result["solverVersion"] = json!(env!("CARGO_PKG_VERSION"));
    if std::fs::read(path).map_err(|e| e.to_string())? != bytes {
        return Err("input changed during solve".into());
    }
    std::fs::create_dir_all(out).map_err(|e| e.to_string())?;
    std::fs::write(
        out.join("thermal3d.json"),
        serde_json::to_vec(&result).unwrap(),
    )
    .map_err(|e| e.to_string())?;
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;
    fn slab() -> Spec {
        Spec {
            name: "slab".into(),
            ambient_c: 20.0,
            pitch_mm: 1.0,
            bodies: vec![Body {
                name: "solid".into(),
                boxes_mm: vec![[0., 0., 0., 10., 10., 10.]],
                conductivity_w_mk: 10.,
                convection_w_m2k: [0., 0., 0., 0., 0., 100.],
                watts: 1.,
                provenance: "analytical fixture".into(),
            }],
            contacts: vec![],
            assumptions: vec!["uniform volumetric heat, only top convects".into()],
        }
    }
    #[test]
    fn slab_matches_analytic_and_balance() {
        let s = slab();
        let r = solve(&s).unwrap();
        let max = r["bodies"][0]["maxC"].as_f64().unwrap(); // T(bottom cell center) = Ta+Q/hA+q'''*(L²-z²)/2k
        let expected = 20. + 100. + 1000000. * (0.01f64.powi(2) - 0.0005f64.powi(2)) / 20.;
        assert!((max - expected).abs() < 0.03);
        assert!(r["relativeEnergyError"].as_f64().unwrap() < 1e-8);
    }
    #[test]
    fn insulated_island_rejected() {
        let mut s = slab();
        s.bodies[0].convection_w_m2k = [0.; 6];
        assert!(solve(&s).unwrap_err().contains("island"));
    }
    #[test]
    fn overlap_rejected() {
        let mut s = slab();
        let mut b = s.bodies[0].clone();
        b.name = "other".into();
        s.bodies.push(b);
        assert!(solve(&s).unwrap_err().contains("overlapping"));
    }
    #[test]
    fn contact_drop_and_missing_contact() {
        let mut s = slab();
        s.bodies[0].convection_w_m2k = [0.; 6];
        let mut b = s.bodies[0].clone();
        b.name = "sink".into();
        b.boxes_mm = vec![[0., 0., 10., 10., 10., 20.]];
        b.watts = 0.;
        b.convection_w_m2k[5] = 100.;
        s.bodies.push(b);
        assert!(solve(&s).unwrap_err().contains("explicit contact"));
        s.contacts.push(Contact {
            bodies: ["solid".into(), "sink".into()],
            resistance_m2k_w: 0.,
            provenance: "fixture".into(),
        });
        let a = solve(&s).unwrap();
        s.contacts[0].resistance_m2k_w = 0.001;
        let b = solve(&s).unwrap();
        assert!(
            (b["bodies"][0]["maxC"].as_f64().unwrap()
                - a["bodies"][0]["maxC"].as_f64().unwrap()
                - 10.)
                .abs()
                < 1e-5
        );
    }
    #[test]
    fn wrong_contact_rejected() {
        let mut s = slab();
        s.contacts.push(Contact {
            bodies: ["solid".into(), "missing".into()],
            resistance_m2k_w: 0.,
            provenance: "test".into(),
        });
        assert!(solve(&s).is_err());
    }
    #[test]
    fn zero_power_is_ambient() {
        let mut s = slab();
        s.bodies[0].watts = 0.;
        let r = solve(&s).unwrap();
        assert_eq!(r["bodies"][0]["maxC"], 20.);
    }
    #[test]
    fn refinement_converges() {
        let mut s = slab();
        s.pitch_mm = 2.;
        let a = solve(&s).unwrap();
        s.pitch_mm = 0.5;
        let b = solve(&s).unwrap();
        assert!(
            (a["bodies"][0]["maxC"].as_f64().unwrap() - b["bodies"][0]["maxC"].as_f64().unwrap())
                .abs()
                < 0.1
        );
    }
}