Report factor value bytes and preserve overflow-safe reduction limits

crates/solver/src/analysis.rs

@@ -266,9 +266,7 @@ pub(crate) fn factor_symmetric(a: &Csr, limits: &Limits, kind: Kind) -> std::res
 /// a singular pivot.
 pub(crate) fn factor_error(e: FactorError, mesh: &Mesh, reduction: &Reduction) -> Error {
     match e {
-        FactorError::TooLarge { entries, limit } => {
+        FactorError::TooLarge { entries, limit } => Error::factor_limit(entries, limit),
-            Error::Limit(format!("the factorization needs {entries} entries; the limit is {limit} (refine less or raise the limit)"))
-        }
         FactorError::OutOfMemory { entries } => Error::Limit(format!("not enough memory for the factorization ({entries} entries)")),
         FactorError::Internal(what) => Error::Internal(what),
         FactorError::Cancelled => Error::Model(fieldcraft_model::Error::Cancelled),
@@ -733,3 +731,20 @@ impl StepResult {
         .map(Location::Point)
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn raised_factor_budget_keeps_the_lanczos_vector_limit() {
+        let neq = Limits::MAX_DOFS;
+        for (entries, expected) in [(400_000_000, 33), (2_000_000_000, 166)] {
+            let limits = Limits { max_factor_entries: entries, ..Limits::default() };
+            let vectors = lanczos_basis(&limits, neq);
+            assert_eq!(vectors, expected);
+            assert!(2 * neq * vectors <= entries);
+            assert_eq!(lanczos_basis(&Limits { max_lanczos_vectors: 7, ..limits }, neq), 7);
+        }
+    }
+}

crates/solver/src/constraints.rs

@@ -367,7 +367,7 @@ impl Budget {
     }
 
     fn exceeded(&self, entries: usize) -> bool {
-        self.bytes(entries) > self.max_bytes
+        entries > self.max_bytes / (8 * (1 + self.matrices))
     }
 }
 
@@ -770,3 +770,18 @@ impl Reduction {
         out
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn reduced_matrix_budget_refuses_byte_counts_that_overflow() {
+        for matrices in [1, 2] {
+            let budget = Budget { max_bytes: usize::MAX, pattern_entries: 0, matrices };
+            let entries = usize::MAX / (8 * (1 + matrices));
+            assert!(!budget.exceeded(entries));
+            assert!(budget.exceeded(entries + 1));
+        }
+    }
+}

crates/solver/src/heat.rs

@@ -303,9 +303,7 @@ fn prescribed(model: &Model, resolved: &Resolved, mesh: &HeatMesh, step: &Step)
 
 fn factor(a: &Csr, limits: &Limits, mesh: &HeatMesh, eqs: &Reduction) -> Result<Factor> {
     Factor::new(&UpperColumns::from_symmetric(a), limits.max_factor_entries).map_err(|e| match e {
-        FactorError::TooLarge { entries, limit } => {
+        FactorError::TooLarge { entries, limit } => Error::factor_limit(entries, limit),
-            Error::Limit(format!("the factorization needs {entries} entries; the limit is {limit} (refine less or raise the limit)"))
-        }
         FactorError::OutOfMemory { entries } => Error::Limit(format!("not enough memory for the factorization ({entries} entries)")),
         FactorError::Internal(what) => Error::Internal(what),
         FactorError::Cancelled => Error::Model(fieldcraft_model::Error::Cancelled),

crates/solver/src/lib.rs

@@ -76,6 +76,14 @@ pub enum Error {
 pub type Result<T> = std::result::Result<T, Error>;
 
 impl Error {
+    /// Factor-value bytes alone; indices and numeric workspace need additional memory.
+    pub(crate) fn factor_limit(entries: usize, limit: usize) -> Error {
+        let bytes = entries as u128 * std::mem::size_of::<f64>() as u128;
+        Error::Limit(format!(
+            "the factorization needs {entries} entries ({bytes} bytes of factor values alone; indices and workspace are additional); the limit is {limit} entries (refine less or raise max_factor_entries)"
+        ))
+    }
+
     /// A kernel error, naming the element.
     pub fn element(label: u32, source: fieldcraft_fem::Error) -> Error {
         Error::Element { label, source }
@@ -112,3 +120,19 @@ impl Default for Limits {
         }
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn factor_limit_estimates_bytes_without_allocating_the_factor() {
+        for (entries, bytes) in [(467_000_000, "3736000000"), (2_000_000_000, "16000000000")] {
+            let e = Error::factor_limit(entries, 400_000_000).to_string();
+            assert!(e.contains(&format!("{entries} entries ({bytes} bytes of factor values alone")), "{e}");
+            assert!(e.contains("indices and workspace are additional"), "{e}");
+            assert!(e.contains("limit is 400000000 entries"), "{e}");
+            assert!(e.contains("raise max_factor_entries"), "{e}");
+        }
+    }
+}

crates/solver/tests/couplings.rs

@@ -368,16 +368,17 @@ fn coupling_errors_name_the_problem() {
 /// budget stops while counting, before allocating.
 #[test]
 fn oversized_reductions_are_refused_before_allocating() {
-    let mut g = grid([83, 83, 1], ElementKind::Hex8, |q| [83.0 * q[0], 83.0 * q[1], q[2]]);
+    let mut g = grid([110, 110, 1], ElementKind::Hex8, |q| [110.0 * q[0], 110.0 * q[1], q[2]]);
     g.define(iso(1000.0, 0.3), None, full(), Procedure::Static {});
-    g.couple("Base", [41.5, 41.5, 0.0], CouplingKind::Distributing, &ALL, |q| near(q[2], 0.0));
+    g.couple("Base", [55.0, 55.0, 0.0], CouplingKind::Distributing, &ALL, |q| near(q[2], 0.0));
     g.bc("hold", BoundaryCondition::Fixed { target: rp("Base") });
     g.face_set("top", |q| near(q[2], 1.0));
     g.load("p", Load::Pressure { set: "top".into(), magnitude: 1.0 });
     let start = std::time::Instant::now();
-    let e = solve_static(&g.model, Some("Step"), &Limits::default()).unwrap_err().to_string();
+    let raised = Limits { max_factor_entries: 2_000_000_000, ..Limits::default() };
+    let e = solve_static(&g.model, Some("Step"), &raised).unwrap_err().to_string();
     eprintln!("refused in {:?}: {e}", start.elapsed());
-    assert!(e.contains("couples 7056 nodes") && e.contains("GB") && e.contains("kinematic coupling"), "{e}");
+    assert!(e.contains("couples 12321 nodes") && e.contains("GB") && e.contains("kinematic coupling"), "{e}");
 
     let mut g = grid([2, 1, 1], ElementKind::Hex8, |q| [2.0 * q[0], q[1], q[2]]);
     g.define(iso(1000.0, 0.3), None, full(), Procedure::Static {});

crates/solver/tests/thermal.rs

@@ -438,3 +438,17 @@ fn trusses_take_thermal_strain_along_their_axis() {
         assert!(r.elements.iter().flat_map(|el| &el.points).all(|p| p.stress[0].abs() < 1e-9));
     }
 }
+
+#[test]
+fn heat_factor_cap_reports_value_bytes_before_allocation() {
+    let mut g = grid([2, 1, 1], ElementKind::Hex8, |p| [2.0 * p[0], p[1], p[2]]);
+    g.define(iso(1000.0, 0.3), None, full(), Procedure::HeatTransfer { transient: None });
+    thermal(&mut g, iso(1000.0, 0.3), None, conductor(K));
+    hold(&mut g, "cold", 0.0, |p| near(p[0], 0.0));
+    let before = g.model.clone();
+    let tight = Limits { max_factor_entries: 1, ..Limits::default() };
+    let e = solve_heat(&g.model, None, &tight).unwrap_err().to_string();
+    assert!(e.contains("bytes of factor values alone") && e.contains("limit is 1 entries"), "{e}");
+    assert_eq!(g.model, before);
+    assert!(solve_heat(&g.model, None, &Limits::default()).is_ok());
+}

crates/solver/tests/verification.rs

@@ -362,7 +362,8 @@ fn errors_name_the_problem() {
     let tight = Limits { max_dofs: 10, ..Limits::default() };
     assert!(solve_static(&g.model, Some("Step"), &tight).unwrap_err().to_string().contains("limit"));
     let tight = Limits { max_factor_entries: 5, ..Limits::default() };
-    assert!(solve_static(&g.model, Some("Step"), &tight).unwrap_err().to_string().contains("limit"));
+    let e = solve_static(&g.model, Some("Step"), &tight).unwrap_err().to_string();
+    assert!(e.contains("bytes of factor values alone") && e.contains("limit is 5 entries"), "{e}");
     assert!(solve_static(&g.model, Some("Nope"), &Limits::default()).is_err());
     assert!(solve_frequency(&g.model, Some("Step"), &FrequencyOptions::default(), &Limits::default()).is_err(), "a static step");
     assert!(solve_static(&g.model, Some("Step"), &Limits::default()).is_ok());