diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7d62eb3..ce394af 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -54,6 +54,12 @@ Modified LSMR is now the sole iterative solver, replacing CG and GMRES.
   `FePreconditioner::additive_schwarz_diagnostics`, and the Python
   `AdditiveSchwarzDiagnostics` class. Scheduling metrics are now private
   to the `Auto` heuristic (closes #34).
+- **BREAKING:** `Operator::apply` / `apply_adjoint` now return
+  `Result<(), SolveError>`; the infallible pair and the `try_apply` /
+  `try_apply_adjoint` defaults are gone. `ApplyError` is removed; its
+  variants moved onto `SolveError`. `PyFePreconditioner.apply` raises
+  `RuntimeError` on local-solver failure instead of returning NaNs
+  (closes #29).
 
 ## [0.1.0] - 2026-03-12
 
diff --git a/crates/schwarz-precond/examples/additive_schwarz.rs b/crates/schwarz-precond/examples/additive_schwarz.rs
index f6a5141..d02cb6a 100644
--- a/crates/schwarz-precond/examples/additive_schwarz.rs
+++ b/crates/schwarz-precond/examples/additive_schwarz.rs
@@ -4,8 +4,8 @@
 //! subdomains and diagonal local solvers.
 
 use schwarz_precond::{
-    lsmr, mlsmr, LocalSolveError, LocalSolver, Operator, SchwarzPreconditioner, SubdomainCore,
-    SubdomainEntry,
+    lsmr, mlsmr, LocalSolveError, LocalSolver, Operator, SchwarzPreconditioner, SolveError,
+    SubdomainCore, SubdomainEntry,
 };
 
 // ---------------------------------------------------------------------------
@@ -23,7 +23,7 @@ impl Operator for TridiagOperator {
     fn ncols(&self) -> usize {
         self.n
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         for i in 0..self.n {
             y[i] = 3.0 * x[i];
             if i > 0 {
@@ -33,12 +33,12 @@ impl Operator for TridiagOperator {
                 y[i] -= x[i + 1];
             }
         }
+        Ok(())
     }
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-        self.apply(x, y); // symmetric
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
+        self.apply(x, y) // symmetric
     }
 }
-
 // ---------------------------------------------------------------------------
 // Diagonal local solver: y = rhs / diag_val
 // ---------------------------------------------------------------------------
diff --git a/crates/schwarz-precond/examples/custom_local_solver.rs b/crates/schwarz-precond/examples/custom_local_solver.rs
index 10a1218..2ca8966 100644
--- a/crates/schwarz-precond/examples/custom_local_solver.rs
+++ b/crates/schwarz-precond/examples/custom_local_solver.rs
@@ -6,10 +6,9 @@
 use faer::{MatRef, Side};
 
 use schwarz_precond::{
-    mlsmr, LocalSolveError, LocalSolver, Operator, SchwarzPreconditioner, SparseMatrix,
+    mlsmr, LocalSolveError, LocalSolver, Operator, SchwarzPreconditioner, SolveError, SparseMatrix,
     SubdomainCore, SubdomainEntry,
 };
-
 // ---------------------------------------------------------------------------
 // Tridiagonal operator
 // ---------------------------------------------------------------------------
@@ -25,7 +24,7 @@ impl Operator for TridiagOperator {
     fn ncols(&self) -> usize {
         self.n
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         for i in 0..self.n {
             y[i] = 3.0 * x[i];
             if i > 0 {
@@ -35,12 +34,12 @@ impl Operator for TridiagOperator {
                 y[i] -= x[i + 1];
             }
         }
+        Ok(())
     }
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-        self.apply(x, y);
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
+        self.apply(x, y)
     }
 }
-
 // ---------------------------------------------------------------------------
 // Build the same tridiag as a SparseMatrix (CSR) for submatrix extraction
 // ---------------------------------------------------------------------------
diff --git a/crates/schwarz-precond/src/error.rs b/crates/schwarz-precond/src/error.rs
index 902eecc..89f8855 100644
--- a/crates/schwarz-precond/src/error.rs
+++ b/crates/schwarz-precond/src/error.rs
@@ -1,14 +1,13 @@
 //! Error types for the `schwarz-precond` crate.
 //!
-//! Errors are layered to match the build → apply → solve lifecycle:
+//! Errors are layered to match the build → solve lifecycle:
 //!
 //! - **Build errors** ([`SubdomainCoreBuildError`], [`SubdomainEntryBuildError`],
 //!   [`PreconditionerBuildError`]) — caught during construction, before any
 //!   solve begins.
-//! - **Apply errors** ([`ApplyError`]) — runtime failures during a single
-//!   preconditioner or operator application (e.g. a local solver diverges).
-//! - **Solve errors** ([`SolveError`]) — wraps `ApplyError` for the iterative
-//!   solver layer.
+//! - **Solve errors** ([`SolveError`]) — runtime failures during a solve,
+//!   including operator/preconditioner application (e.g. a local solver
+//!   diverges) and iterative-solver input validation.
 //!
 //! Each level chains to its source via [`Error::source`].
 
@@ -145,51 +144,25 @@ impl Display for LocalSolveError {
 
 impl Error for LocalSolveError {}
 
-/// Runtime failure while applying a Schwarz preconditioner/operator.
+/// Runtime failure while executing a solve.
+///
+/// Covers both operator/preconditioner application failures (e.g. a local
+/// solver diverges) and iterative-solver input validation.
 #[derive(Debug, Clone, PartialEq, Eq)]
-pub enum ApplyError {
-    /// A local subdomain solve failed.
+#[non_exhaustive]
+pub enum SolveError {
+    /// A local subdomain solve failed during a preconditioner apply.
     LocalSolveFailed {
         /// Index of the failing subdomain entry in the preconditioner.
         subdomain: usize,
         /// Local solver error.
         source: LocalSolveError,
     },
-    /// Internal synchronization failed (e.g. poisoned mutex).
+    /// Internal synchronization failed (e.g. poisoned mutex) during an apply.
     Synchronization {
         /// Context string identifying the lock/synchronization site.
         context: &'static str,
     },
-}
-
-impl Display for ApplyError {
-    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
-        match self {
-            Self::LocalSolveFailed { subdomain, source } => {
-                write!(f, "subdomain {subdomain} local solve failed: {source}")
-            }
-            Self::Synchronization { context } => {
-                write!(f, "synchronization failure at {context}")
-            }
-        }
-    }
-}
-
-impl Error for ApplyError {
-    fn source(&self) -> Option<&(dyn Error + 'static)> {
-        match self {
-            Self::LocalSolveFailed { source, .. } => Some(source),
-            Self::Synchronization { .. } => None,
-        }
-    }
-}
-
-/// Runtime failure while executing an iterative solver.
-#[derive(Debug, Clone, PartialEq, Eq)]
-#[non_exhaustive]
-pub enum SolveError {
-    /// Operator/preconditioner apply failed.
-    Apply(ApplyError),
     /// Solver input was invalid before any iteration was attempted.
     InvalidInput {
         /// Context string identifying the validation site.
@@ -202,7 +175,12 @@ pub enum SolveError {
 impl Display for SolveError {
     fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
         match self {
-            Self::Apply(err) => write!(f, "operator apply failed: {err}"),
+            Self::LocalSolveFailed { subdomain, source } => {
+                write!(f, "subdomain {subdomain} local solve failed: {source}")
+            }
+            Self::Synchronization { context } => {
+                write!(f, "synchronization failure at {context}")
+            }
             Self::InvalidInput { context, message } => {
                 write!(f, "invalid solver input at {context}: {message}")
             }
@@ -213,18 +191,13 @@ impl Display for SolveError {
 impl Error for SolveError {
     fn source(&self) -> Option<&(dyn Error + 'static)> {
         match self {
-            Self::Apply(err) => Some(err),
+            Self::LocalSolveFailed { source, .. } => Some(source),
+            Self::Synchronization { .. } => None,
             Self::InvalidInput { .. } => None,
         }
     }
 }
 
-impl From<ApplyError> for SolveError {
-    fn from(value: ApplyError) -> Self {
-        Self::Apply(value)
-    }
-}
-
 pub(crate) fn validate_entries<S: LocalSolver>(
     entries: &[SubdomainEntry<S>],
     n_dofs: usize,
diff --git a/crates/schwarz-precond/src/lib.rs b/crates/schwarz-precond/src/lib.rs
index c98d602..a282ac4 100644
--- a/crates/schwarz-precond/src/lib.rs
+++ b/crates/schwarz-precond/src/lib.rs
@@ -92,31 +92,22 @@
 ///
 /// Preconditioners are operators too (M^{-1} is a linear map).
 /// All implementors must be Send + Sync to enable Rayon parallelism.
+///
+/// Both apply methods are fallible: implementors that cannot fail in practice
+/// (matrices, identity operators) still return `Result<(), SolveError>` so
+/// callers can use a uniform `?` propagation path. Symmetric operators
+/// should delegate `apply_adjoint` to `apply`.
 pub trait Operator: Send + Sync {
     /// Number of rows in the operator.
     fn nrows(&self) -> usize;
     /// Number of columns in the operator.
     fn ncols(&self) -> usize;
-    /// y = A*x
-    fn apply(&self, x: &[f64], y: &mut [f64]);
+    /// Computes y = A * x. Returns an error if the apply fails at runtime
+    /// (e.g. a local subdomain solver diverges).
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), error::SolveError>;
     /// Computes y = A^T * x. For symmetric operators, this should delegate to `apply`.
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]);
-
-    /// Fallible version of [`Operator::apply`].
-    ///
-    /// Implementors with runtime failure modes should override this.
-    fn try_apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), error::ApplyError> {
-        self.apply(x, y);
-        Ok(())
-    }
-
-    /// Fallible version of [`Operator::apply_adjoint`].
-    ///
-    /// Implementors with runtime failure modes should override this.
-    fn try_apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), error::ApplyError> {
-        self.apply_adjoint(x, y);
-        Ok(())
-    }
+    /// Returns an error under the same conditions as `apply`.
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), error::SolveError>;
 }
 
 // ============================================================================
@@ -133,7 +124,7 @@ mod sparse_matrix;
 
 pub use domain::{PartitionWeights, SubdomainCore};
 pub use error::{
-    ApplyError, LocalSolveError, PreconditionerBuildError, SolveError, SubdomainCoreBuildError,
+    LocalSolveError, PreconditionerBuildError, SolveError, SubdomainCoreBuildError,
     SubdomainEntryBuildError,
 };
 pub use local_solve::{LocalSolver, SubdomainEntry};
diff --git a/crates/schwarz-precond/src/lsmr/bidiag.rs b/crates/schwarz-precond/src/lsmr/bidiag.rs
index 56f1444..adc26f8 100644
--- a/crates/schwarz-precond/src/lsmr/bidiag.rs
+++ b/crates/schwarz-precond/src/lsmr/bidiag.rs
@@ -278,7 +278,7 @@ pub(super) trait Bidiagonalization {
 impl<A: Operator + ?Sized> Bidiagonalization for GolubKahan<'_, A> {
     fn step(&mut self) -> Result<BidiagStep, SolveError> {
         // β_{k+1} u_{k+1} = A v_k − α_k u_k
-        self.operator.try_apply(&self.bufs.v, &mut self.bufs.av)?;
+        self.operator.apply(&self.bufs.v, &mut self.bufs.av)?;
         let beta_sq = axpy_with_sq_norm(&mut self.bufs.u, &self.bufs.av, -self.alpha);
         let beta = beta_sq.sqrt();
         if beta == 0.0 {
@@ -295,7 +295,7 @@ impl<A: Operator + ?Sized> Bidiagonalization for GolubKahan<'_, A> {
 
         // α_{k+1} v_{k+1} = Aᵀ u_{k+1} − β_{k+1} v_k
         self.operator
-            .try_apply_adjoint(&self.bufs.u, &mut self.bufs.atu)?;
+            .apply_adjoint(&self.bufs.u, &mut self.bufs.atu)?;
         let mut alpha_sq = axpy_with_sq_norm(&mut self.bufs.v, &self.bufs.atu, -beta);
 
         // Windowed MGS in the standard inner product, before normalization.
@@ -330,7 +330,7 @@ impl<A: Operator + ?Sized, M: Operator + ?Sized> Bidiagonalization
         // Phase 1 — ũ_{k+1} unnormalized: scale = −α_k / β_k.
         // Compute ũ_{k+1} = A ṽ_k − (α_k / β_k) ũ_k, then β_{k+1} = ‖ũ_{k+1}‖.
         let scale = -(self.alpha * self.beta_prev_inv);
-        self.operator.try_apply(&self.bufs.v, &mut self.bufs.av)?;
+        self.operator.apply(&self.bufs.v, &mut self.bufs.av)?;
         let beta_sq = axpy_with_sq_norm(&mut self.bufs.u, &self.bufs.av, scale);
         let beta = beta_sq.sqrt();
         if beta == 0.0 {
@@ -353,7 +353,7 @@ impl<A: Operator + ?Sized, M: Operator + ?Sized> Bidiagonalization
         // Precondition: α_k > 0 (the outer loop guard in lsmr_from_bidiag
         // never calls step() once α has collapsed to zero).
         self.operator
-            .try_apply_adjoint(&self.bufs.u, &mut self.bufs.atu)?;
+            .apply_adjoint(&self.bufs.u, &mut self.bufs.atu)?;
         debug_assert!(
             self.alpha > 0.0,
             "self.alpha must be > 0; lsmr_from_bidiag's loop guard prevents step() after alpha=0",
@@ -371,7 +371,7 @@ impl<A: Operator + ?Sized, M: Operator + ?Sized> Bidiagonalization
         // step), and push the normalized v_{k+1} and p̃_{k+1,norm} = p_tilde /
         // α_new into the ring for the next step's MGS.
         self.preconditioner
-            .try_apply(&self.bufs.p_tilde, &mut self.bufs.v)?;
+            .apply(&self.bufs.p_tilde, &mut self.bufs.v)?;
 
         if let Some(reorth) = &self.bufs.local_reorth {
             reorth.reorthogonalize(&mut self.bufs.v, &mut self.bufs.p_tilde);
@@ -464,7 +464,7 @@ impl<'a, A: Operator + ?Sized> GolubKahan<'a, A> {
         }
 
         // α₁ v₁ = Aᵀ u₁
-        operator.try_apply_adjoint(&bufs.u, &mut bufs.v)?;
+        operator.apply_adjoint(&bufs.u, &mut bufs.v)?;
         let alpha = par_dot(&bufs.v, &bufs.v).sqrt();
         if alpha > 0.0 {
             scale_in_place(&mut bufs.v, 1.0 / alpha);
@@ -555,10 +555,10 @@ impl<'a, A: Operator + ?Sized, M: Operator + ?Sized> ModifiedGolubKahan<'a, A, M
         }
 
         // p̃ = Aᵀ u₁
-        operator.try_apply_adjoint(&bufs.u, &mut bufs.p_tilde)?;
+        operator.apply_adjoint(&bufs.u, &mut bufs.p_tilde)?;
 
         // ṽ₁ = M⁻¹ p̃
-        preconditioner.try_apply(&bufs.p_tilde, &mut bufs.v)?;
+        preconditioner.apply(&bufs.p_tilde, &mut bufs.v)?;
 
         // α₁ = √⟨ṽ₁, p̃⟩ via the M-norm dot product trick.
         let vp = par_dot(&bufs.v, &bufs.p_tilde);
diff --git a/crates/schwarz-precond/src/lsmr/tests.rs b/crates/schwarz-precond/src/lsmr/tests.rs
index 44aeced..3a74fc1 100644
--- a/crates/schwarz-precond/src/lsmr/tests.rs
+++ b/crates/schwarz-precond/src/lsmr/tests.rs
@@ -17,11 +17,13 @@ impl Operator for IdentityOp {
     fn ncols(&self) -> usize {
         self.n
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         y.copy_from_slice(x);
+        Ok(())
     }
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         y.copy_from_slice(x);
+        Ok(())
     }
 }
 
@@ -36,16 +38,18 @@ impl Operator for OverdeterminedOp {
     fn ncols(&self) -> usize {
         3
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         y[0] = x[0];
         y[1] = x[1];
         y[2] = x[2];
         y[3] = x[0] + x[1];
+        Ok(())
     }
-    fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) {
+    fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) -> Result<(), SolveError> {
         x[0] = u[0] + u[3];
         x[1] = u[1] + u[3];
         x[2] = u[2];
+        Ok(())
     }
 }
 
@@ -60,13 +64,14 @@ impl Operator for DiagPrecond {
     fn ncols(&self) -> usize {
         3
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         y[0] = x[0] / 2.0;
         y[1] = x[1] / 2.0;
         y[2] = x[2];
+        Ok(())
     }
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-        self.apply(x, y);
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
+        self.apply(x, y)
     }
 }
 
@@ -74,10 +79,11 @@ impl Operator for DiagPrecond {
 /// "did we actually solve the least-squares problem?" check.
 fn normal_equation_residual<O: Operator + ?Sized>(op: &O, x: &[f64], b: &[f64]) -> f64 {
     let mut ax = vec![0.0; op.nrows()];
-    op.apply(x, &mut ax);
+    op.apply(x, &mut ax).expect("apply succeeds");
     let resid: Vec<f64> = b.iter().zip(&ax).map(|(bi, ai)| bi - ai).collect();
     let mut atr = vec![0.0; op.ncols()];
-    op.apply_adjoint(&resid, &mut atr);
+    op.apply_adjoint(&resid, &mut atr)
+        .expect("apply_adjoint succeeds");
     vec_norm(&atr)
 }
 
@@ -137,14 +143,16 @@ fn test_mlsmr_underdetermined_system() {
         fn ncols(&self) -> usize {
             3
         }
-        fn apply(&self, x: &[f64], y: &mut [f64]) {
+        fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
             y[0] = x[0];
             y[1] = x[1];
+            Ok(())
         }
-        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) {
+        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) -> Result<(), SolveError> {
             x[0] = u[0];
             x[1] = u[1];
             x[2] = 0.0;
+            Ok(())
         }
     }
 
@@ -166,15 +174,17 @@ fn test_mlsmr_rank_deficient_system() {
         fn ncols(&self) -> usize {
             2
         }
-        fn apply(&self, x: &[f64], y: &mut [f64]) {
+        fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
             let s = x[0] + x[1];
             y[0] = s;
             y[1] = 2.0 * s;
+            Ok(())
         }
-        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) {
+        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) -> Result<(), SolveError> {
             let s = u[0] + 2.0 * u[1];
             x[0] = s;
             x[1] = s;
+            Ok(())
         }
     }
 
@@ -195,13 +205,15 @@ fn test_mlsmr_zero_column_and_zero_row() {
         fn ncols(&self) -> usize {
             2
         }
-        fn apply(&self, x: &[f64], y: &mut [f64]) {
+        fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
             y[0] = x[0];
             y[1] = 0.0;
+            Ok(())
         }
-        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) {
+        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) -> Result<(), SolveError> {
             x[0] = u[0];
             x[1] = 0.0;
+            Ok(())
         }
     }
 
@@ -227,13 +239,15 @@ fn test_mlsmr_mid_stream_beta_zero_breakdown() {
         fn ncols(&self) -> usize {
             2
         }
-        fn apply(&self, x: &[f64], y: &mut [f64]) {
+        fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
             y[0] = x[0];
             y[1] = 0.0;
+            Ok(())
         }
-        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) {
+        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) -> Result<(), SolveError> {
             x[0] = u[0];
             x[1] = 0.0;
+            Ok(())
         }
     }
 
@@ -394,12 +408,13 @@ impl Operator for DenseOp {
     fn ncols(&self) -> usize {
         self.cols
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         for (yi, row) in y.iter_mut().zip(self.data.chunks_exact(self.cols)) {
             *yi = row.iter().zip(x).map(|(a, b)| a * b).sum();
         }
+        Ok(())
     }
-    fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) {
+    fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) -> Result<(), SolveError> {
         for (j, xj) in x.iter_mut().enumerate() {
             let mut s = 0.0;
             for (ui, row) in u.iter().zip(self.data.chunks_exact(self.cols)) {
@@ -407,6 +422,7 @@ impl Operator for DenseOp {
             }
             *xj = s;
         }
+        Ok(())
     }
 }
 
@@ -495,13 +511,14 @@ fn test_mlsmr_local_reorth_preconditioned() {
         fn ncols(&self) -> usize {
             self.0.len()
         }
-        fn apply(&self, x: &[f64], y: &mut [f64]) {
+        fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
             for ((yi, &xi), &di) in y.iter_mut().zip(x).zip(self.0.iter()) {
                 *yi = di * xi;
             }
+            Ok(())
         }
-        fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-            self.apply(x, y);
+        fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
+            self.apply(x, y)
         }
     }
     let m = DiagOp(diag_inv);
@@ -602,12 +619,14 @@ fn test_mlsmr_step1_alpha_zero_early_exit() {
         fn ncols(&self) -> usize {
             1
         }
-        fn apply(&self, x: &[f64], y: &mut [f64]) {
+        fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
             y[0] = x[0];
             y[1] = 0.0;
+            Ok(())
         }
-        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) {
+        fn apply_adjoint(&self, u: &[f64], x: &mut [f64]) -> Result<(), SolveError> {
             x[0] = u[0];
+            Ok(())
         }
     }
 
@@ -696,11 +715,12 @@ fn test_mlsmr_rejects_bad_preconditioner_shape() {
         fn ncols(&self) -> usize {
             2
         }
-        fn apply(&self, x: &[f64], y: &mut [f64]) {
+        fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
             y.copy_from_slice(x);
+            Ok(())
         }
-        fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-            self.apply(x, y);
+        fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
+            self.apply(x, y)
         }
     }
 
diff --git a/crates/schwarz-precond/src/schwarz/executor.rs b/crates/schwarz-precond/src/schwarz/executor.rs
index b33f2b6..0674593 100644
--- a/crates/schwarz-precond/src/schwarz/executor.rs
+++ b/crates/schwarz-precond/src/schwarz/executor.rs
@@ -1,7 +1,7 @@
 //! Additive Schwarz execution engine.
 //!
 //! [`AdditiveExecutor`] owns the subdomain entries and dispatches
-//! `try_apply` using the reduction plan chosen by the scheduler.
+//! `apply` using the reduction plan chosen by the scheduler.
 //! It manages the [`BufferPool`] for zero-allocation steady-state
 //! operation.
 //!
@@ -26,7 +26,7 @@ use std::sync::{Arc, Mutex};
 use rayon::prelude::*;
 use thread_local::ThreadLocal;
 
-use crate::error::ApplyError;
+use crate::error::SolveError;
 use crate::local_solve::{LocalSolver, SubdomainEntry};
 
 use super::planning::{ReductionPlan, ResolvedReductionStrategy};
@@ -48,8 +48,8 @@ impl BufferPool {
         strategy: ResolvedReductionStrategy,
         n_dofs: usize,
         max_scratch_size: usize,
-    ) -> Result<SchwarzBuffers, ApplyError> {
-        let mut pool = self.inner.lock().map_err(|_| ApplyError::Synchronization {
+    ) -> Result<SchwarzBuffers, SolveError> {
+        let mut pool = self.inner.lock().map_err(|_| SolveError::Synchronization {
             context: "additive.buf_pool.lock.pop",
         })?;
         if let Some(idx) = pool.iter().position(|bufs| bufs.strategy() == strategy) {
@@ -61,15 +61,15 @@ impl BufferPool {
     fn put(
         &self,
         bufs: SchwarzBuffers,
-        apply_result: &Result<(), ApplyError>,
-    ) -> Result<(), ApplyError> {
+        apply_result: &Result<(), SolveError>,
+    ) -> Result<(), SolveError> {
         if let Ok(mut pool) = self.inner.lock() {
             if pool.len() < Self::MAX_POOL_SIZE {
                 pool.push(bufs);
             }
             Ok(())
         } else if apply_result.is_ok() {
-            Err(ApplyError::Synchronization {
+            Err(SolveError::Synchronization {
                 context: "additive.buf_pool.lock.push",
             })
         } else {
@@ -178,8 +178,8 @@ impl WorkerReductionBuffers {
     fn finish_round(
         self,
         z: &mut [f64],
-        apply_result: &Result<(), ApplyError>,
-    ) -> Result<Vec<AdditiveSweepBuffers>, ApplyError> {
+        apply_result: &Result<(), SolveError>,
+    ) -> Result<Vec<AdditiveSweepBuffers>, SolveError> {
         let mut buffers = self.into_buffers()?;
         if apply_result.is_ok() {
             reduce_into(z, &buffers);
@@ -188,11 +188,11 @@ impl WorkerReductionBuffers {
         Ok(buffers)
     }
 
-    fn into_buffers(mut self) -> Result<Vec<AdditiveSweepBuffers>, ApplyError> {
+    fn into_buffers(mut self) -> Result<Vec<AdditiveSweepBuffers>, SolveError> {
         let mut buffers =
             self.shared_pool
                 .into_inner()
-                .map_err(|_| ApplyError::Synchronization {
+                .map_err(|_| SolveError::Synchronization {
                     context: "additive.reduction.pool.into_inner",
                 })?;
         for worker_stack in self.worker_stacks.iter_mut() {
@@ -278,12 +278,12 @@ impl<S: LocalSolver> AdditiveExecutor<S> {
         }
     }
 
-    pub(super) fn try_apply(
+    pub(super) fn apply(
         &self,
         plan: ReductionPlan,
         r: &[f64],
         z: &mut [f64],
-    ) -> Result<(), ApplyError> {
+    ) -> Result<(), SolveError> {
         let mut bufs = self
             .buf_pool
             .take(plan.strategy, self.n_dofs, self.max_scratch_size)?;
@@ -305,7 +305,7 @@ impl<S: LocalSolver> AdditiveExecutor<S> {
         r: &[f64],
         z: &mut [f64],
         accum: &[AtomicU64],
-    ) -> Result<(), ApplyError> {
+    ) -> Result<(), SolveError> {
         self.subdomains.par_iter().enumerate().try_for_each_init(
             || LocalSolveScratch::new(self.max_scratch_size),
             |scratch, (subdomain, entry)| {
@@ -317,7 +317,7 @@ impl<S: LocalSolver> AdditiveExecutor<S> {
                         &mut scratch.z_scratch,
                         allow_inner_parallelism,
                     )
-                    .map_err(|source| ApplyError::LocalSolveFailed { subdomain, source })
+                    .map_err(|source| SolveError::LocalSolveFailed { subdomain, source })
             },
         )?;
 
@@ -340,7 +340,7 @@ impl<S: LocalSolver> AdditiveExecutor<S> {
         r: &[f64],
         z: &mut [f64],
         pool: &mut Vec<AdditiveSweepBuffers>,
-    ) -> Result<(), ApplyError> {
+    ) -> Result<(), SolveError> {
         let worker_buffers =
             WorkerReductionBuffers::new(std::mem::take(pool), self.n_dofs, self.max_scratch_size);
         let apply_result =
@@ -357,7 +357,7 @@ impl<S: LocalSolver> AdditiveExecutor<S> {
                                 &mut buffers.scratch.z_scratch,
                                 allow_inner_parallelism,
                             )
-                            .map_err(|source| ApplyError::LocalSolveFailed { subdomain, source })
+                            .map_err(|source| SolveError::LocalSolveFailed { subdomain, source })
                     })
                 });
 
diff --git a/crates/schwarz-precond/src/schwarz/preconditioner.rs b/crates/schwarz-precond/src/schwarz/preconditioner.rs
index 483eb57..7f90cb2 100644
--- a/crates/schwarz-precond/src/schwarz/preconditioner.rs
+++ b/crates/schwarz-precond/src/schwarz/preconditioner.rs
@@ -6,7 +6,7 @@
 //! the executor. `apply` is lock-free in steady state (buffers are
 //! borrowed from a pool).
 
-use crate::error::{validate_entries, ApplyError, PreconditionerBuildError};
+use crate::error::{validate_entries, PreconditionerBuildError, SolveError};
 use crate::local_solve::{LocalSolver, SubdomainEntry};
 use crate::Operator;
 
@@ -131,10 +131,10 @@ impl<S: LocalSolver> SchwarzPreconditioner<S> {
         }
     }
 
-    /// Fallible operator apply that propagates local-solver failures.
-    pub fn try_apply(&self, r: &[f64], z: &mut [f64]) -> Result<(), ApplyError> {
+    /// Operator apply that propagates local-solver failures.
+    pub fn apply(&self, r: &[f64], z: &mut [f64]) -> Result<(), SolveError> {
         let plan = self.reduction_plan();
-        self.executor.try_apply(plan, r, z)
+        self.executor.apply(plan, r, z)
     }
 
     fn reduction_plan(&self) -> ReductionPlan {
@@ -165,21 +165,11 @@ impl<S: LocalSolver> Operator for SchwarzPreconditioner<S> {
         self.executor.n_dofs
     }
 
-    fn apply(&self, r: &[f64], z: &mut [f64]) {
-        if self.try_apply(r, z).is_err() {
-            z.fill(f64::NAN);
-        }
-    }
-
-    fn apply_adjoint(&self, r: &[f64], z: &mut [f64]) {
-        self.apply(r, z);
-    }
-
-    fn try_apply(&self, r: &[f64], z: &mut [f64]) -> Result<(), ApplyError> {
-        SchwarzPreconditioner::try_apply(self, r, z)
+    fn apply(&self, r: &[f64], z: &mut [f64]) -> Result<(), SolveError> {
+        SchwarzPreconditioner::apply(self, r, z)
     }
 
-    fn try_apply_adjoint(&self, r: &[f64], z: &mut [f64]) -> Result<(), ApplyError> {
-        SchwarzPreconditioner::try_apply(self, r, z)
+    fn apply_adjoint(&self, r: &[f64], z: &mut [f64]) -> Result<(), SolveError> {
+        SchwarzPreconditioner::apply(self, r, z)
     }
 }
diff --git a/crates/schwarz-precond/tests/common/mod.rs b/crates/schwarz-precond/tests/common/mod.rs
index 908ae92..d1caf79 100644
--- a/crates/schwarz-precond/tests/common/mod.rs
+++ b/crates/schwarz-precond/tests/common/mod.rs
@@ -2,7 +2,9 @@
 
 #![allow(dead_code)]
 
-use schwarz_precond::{LocalSolveError, LocalSolver, Operator, SubdomainCore, SubdomainEntry};
+use schwarz_precond::{
+    LocalSolveError, LocalSolver, Operator, SolveError, SubdomainCore, SubdomainEntry,
+};
 
 // ---------------------------------------------------------------------------
 // Operators
@@ -27,7 +29,7 @@ impl Operator for TridiagOperator {
     fn ncols(&self) -> usize {
         self.n
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         for i in 0..self.n {
             y[i] = self.diag_val * x[i];
             if i > 0 {
@@ -37,9 +39,10 @@ impl Operator for TridiagOperator {
                 y[i] -= x[i + 1];
             }
         }
+        Ok(())
     }
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-        self.apply(x, y);
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
+        self.apply(x, y)
     }
 }
 
@@ -142,7 +145,7 @@ impl LocalSolver for FailingLocalSolver {
 pub fn check_residual<A: Operator>(op: &A, x: &[f64], b: &[f64], tol: f64) {
     let n = b.len();
     let mut ax = vec![0.0; n];
-    op.apply(x, &mut ax);
+    op.apply(x, &mut ax).expect("apply succeeds");
     let err: f64 = ax
         .iter()
         .zip(b.iter())
diff --git a/crates/schwarz-precond/tests/schwarz.rs b/crates/schwarz-precond/tests/schwarz.rs
index 62e64e8..83058c0 100644
--- a/crates/schwarz-precond/tests/schwarz.rs
+++ b/crates/schwarz-precond/tests/schwarz.rs
@@ -131,8 +131,7 @@ fn run_nested_parallel_reduction_regression_case() {
 
         for _ in 0..8 {
             let mut z = vec![0.0; n];
-            schwarz.apply(&rhs, &mut z);
-
+            schwarz.apply(&rhs, &mut z).expect("apply succeeds");
             for (i, (&zi, &ri)) in z.iter().zip(&rhs).enumerate() {
                 assert!(
                     (zi - 2.0 * ri).abs() <= 1e-12,
@@ -225,13 +224,18 @@ fn test_clone_produces_independent_preconditioner() {
     let mut z_orig = vec![0.0; n];
     let mut z_clone = vec![0.0; n];
 
-    original.apply(&r1, &mut z_orig);
-    cloned.apply(&r2, &mut z_clone);
+    original
+        .apply(&r1, &mut z_orig)
+        .expect("original apply succeeds");
+    cloned
+        .apply(&r2, &mut z_clone)
+        .expect("cloned apply succeeds");
 
     // Verify independently: apply the original with r2 to check the clone's result.
     let mut z_check = vec![0.0; n];
-    original.apply(&r2, &mut z_check);
-
+    original
+        .apply(&r2, &mut z_check)
+        .expect("check apply succeeds");
     for i in 0..n {
         assert!(
             (z_clone[i] - z_check[i]).abs() < 1e-14,
@@ -244,7 +248,7 @@ fn test_clone_produces_independent_preconditioner() {
 
     // Verify the original was not corrupted by the clone's apply.
     let mut z_orig2 = vec![0.0; n];
-    original.apply(&r1, &mut z_orig2);
+    original.apply(&r1, &mut z_orig2).expect("apply succeeds");
     for i in 0..n {
         assert!(
             (z_orig[i] - z_orig2[i]).abs() < 1e-14,
@@ -266,9 +270,10 @@ fn test_additive_schwarz_operator_dimensions() {
     let r = vec![1.0; n];
     let mut z1 = vec![0.0; n];
     let mut z2 = vec![0.0; n];
-    schwarz.apply(&r, &mut z1);
-    schwarz.apply_adjoint(&r, &mut z2);
-
+    schwarz.apply(&r, &mut z1).expect("apply succeeds");
+    schwarz
+        .apply_adjoint(&r, &mut z2)
+        .expect("apply_adjoint succeeds");
     for i in 0..n {
         assert!(
             (z1[i] - z2[i]).abs() < 1e-14,
@@ -299,12 +304,11 @@ fn test_additive_schwarz_parallel_apply_stress_no_panics() {
         .map(|rhs| {
             let mut z = vec![0.0; n];
             for _ in 0..16 {
-                schwarz.apply(rhs, &mut z);
+                schwarz.apply(rhs, &mut z).expect("apply succeeds");
             }
             z
         })
         .collect();
-
     assert_eq!(outputs.len(), rhs_columns.len());
     assert!(
         outputs.iter().flatten().all(|v| v.is_finite()),
@@ -337,9 +341,12 @@ fn test_additive_backends_match_on_overlapping_subdomains() {
 
         let mut z_atomic = vec![0.0; n];
         let mut z_reduction = vec![0.0; n];
-        atomic.apply(&rhs, &mut z_atomic);
-        reduction.apply(&rhs, &mut z_reduction);
-
+        atomic
+            .apply(&rhs, &mut z_atomic)
+            .expect("atomic apply succeeds");
+        reduction
+            .apply(&rhs, &mut z_reduction)
+            .expect("reduction apply succeeds");
         assert_vec_close(&z_atomic, &z_reduction, 1e-12);
     });
 }
@@ -374,48 +381,22 @@ fn test_additive_auto_matches_resolved_backend() {
 
             let mut z_auto = vec![0.0; n];
             let mut z_explicit = vec![0.0; n];
-            auto.apply(&rhs, &mut z_auto);
-            explicit.apply(&rhs, &mut z_explicit);
-
+            auto.apply(&rhs, &mut z_auto).expect("auto apply succeeds");
+            explicit
+                .apply(&rhs, &mut z_explicit)
+                .expect("explicit apply succeeds");
             assert_vec_close(&z_auto, &z_explicit, 1e-12);
         });
     }
 }
 
 use schwarz_precond::{
-    ApplyError, PreconditionerBuildError, SolveError, SubdomainCoreBuildError,
-    SubdomainEntryBuildError,
+    PreconditionerBuildError, SolveError, SubdomainCoreBuildError, SubdomainEntryBuildError,
 };
 use std::error::Error;
 
 // ============================================================================
-// Default try_apply / try_apply_adjoint tests
-// ============================================================================
-
-#[test]
-fn test_try_apply_default_succeeds() {
-    let a = TridiagOperator::new(5, 3.0);
-    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];
-    let mut y_apply = vec![0.0; 5];
-    let mut y_try = vec![0.0; 5];
-    a.apply(&x, &mut y_apply);
-    let result = a.try_apply(&x, &mut y_try);
-    assert!(result.is_ok());
-    assert_eq!(y_apply, y_try);
-}
-
-#[test]
-fn test_try_apply_adjoint_default_succeeds() {
-    let a = TridiagOperator::new(5, 3.0);
-    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];
-    let mut y_apply = vec![0.0; 5];
-    let mut y_try = vec![0.0; 5];
-    a.apply_adjoint(&x, &mut y_apply);
-    let result = a.try_apply_adjoint(&x, &mut y_try);
-    assert!(result.is_ok());
-    assert_eq!(y_apply, y_try);
-}
-
+// Additive Schwarz edge cases
 // ============================================================================
 // Additive Schwarz edge cases
 // ============================================================================
@@ -444,8 +425,9 @@ fn test_additive_schwarz_apply_subdomain_empty_indices() {
             // Verify apply works with empty subdomain
             let r = vec![1.0; 5];
             let mut z = vec![0.0; 5];
-            schwarz.apply(&r, &mut z);
-            // Empty subdomain contributes nothing
+            schwarz
+                .apply(&r, &mut z)
+                .expect("apply with empty subdomain succeeds");
             for &v in &z {
                 assert!((v - 0.0).abs() < 1e-14);
             }
@@ -522,12 +504,12 @@ fn test_local_solve_error_display() {
 }
 
 #[test]
-fn test_apply_error_display_local_solve_failed() {
+fn test_solve_error_display_local_solve_failed() {
     let local_err = LocalSolveError::ApproxCholSolveFailed {
         context: "test",
         message: "fail".to_string(),
     };
-    let err = ApplyError::LocalSolveFailed {
+    let err = SolveError::LocalSolveFailed {
         subdomain: 7,
         source: local_err,
     };
@@ -540,8 +522,8 @@ fn test_apply_error_display_local_solve_failed() {
 }
 
 #[test]
-fn test_apply_error_display_synchronization() {
-    let err = ApplyError::Synchronization {
+fn test_solve_error_display_synchronization() {
+    let err = SolveError::Synchronization {
         context: "mutex.lock",
     };
     let msg = err.to_string();
@@ -553,12 +535,12 @@ fn test_apply_error_display_synchronization() {
 }
 
 #[test]
-fn test_apply_error_source() {
+fn test_solve_error_source() {
     let local_err = LocalSolveError::ApproxCholSolveFailed {
         context: "test",
         message: "err".to_string(),
     };
-    let err = ApplyError::LocalSolveFailed {
+    let err = SolveError::LocalSolveFailed {
         subdomain: 0,
         source: local_err,
     };
@@ -567,34 +549,13 @@ fn test_apply_error_source() {
         "LocalSolveFailed should have a source"
     );
 
-    let err2 = ApplyError::Synchronization { context: "test" };
+    let err2 = SolveError::Synchronization { context: "test" };
     assert!(
         err2.source().is_none(),
         "Synchronization should have no source"
     );
 }
 
-#[test]
-fn test_solve_error_display() {
-    let apply_err = ApplyError::Synchronization { context: "test" };
-    let err = SolveError::Apply(apply_err);
-    let msg = err.to_string();
-    assert!(
-        msg.contains("operator apply failed"),
-        "missing prefix: {msg}"
-    );
-}
-
-#[test]
-fn test_solve_error_source() {
-    let apply_err = ApplyError::Synchronization { context: "test" };
-    let err = SolveError::Apply(apply_err);
-    assert!(
-        err.source().is_some(),
-        "SolveError::Apply should have a source"
-    );
-}
-
 #[test]
 fn test_solve_error_invalid_input_display_and_source() {
     let err = SolveError::InvalidInput {
@@ -606,19 +567,6 @@ fn test_solve_error_invalid_input_display_and_source() {
     assert!(msg.contains("bad dimension"));
     assert!(err.source().is_none());
 }
-
-#[test]
-fn test_solve_error_from_apply_error() {
-    let apply_err = ApplyError::Synchronization { context: "conv" };
-    let solve_err: SolveError = apply_err.into();
-    match solve_err {
-        SolveError::Apply(ApplyError::Synchronization { context }) => {
-            assert_eq!(context, "conv");
-        }
-        _ => panic!("expected SolveError::Apply(Synchronization)"),
-    }
-}
-
 // ============================================================================
 // Validation error tests
 // ============================================================================
@@ -693,22 +641,23 @@ fn test_additive_schwarz_parallel_readout_large_n() {
 
     let rhs = vec![4.0; n];
     let mut z = vec![0.0; n];
-    let result = schwarz.try_apply(&rhs, &mut z);
-    assert!(result.is_ok(), "try_apply should succeed: {:?}", result);
+    let result = schwarz.apply(&rhs, &mut z);
+    assert!(result.is_ok(), "apply should succeed: {:?}", result);
 
     // Each DOF: output = 1.0 * (1.0 * 4.0 / 2.0) = 2.0
     for (i, &v) in z.iter().enumerate() {
-        assert!((v - 2.0).abs() < 1e-12, "z[{i}] = {v}, expected 2.0",);
+        assert!((v - 2.0_f64).abs() < 1e-12, "z[{i}] = {v}, expected 2.0",);
     }
 }
 
 // ============================================================================
-// NaN-fill on solver error (apply, not try_apply)
+// apply propagates local-solver failure
 // ============================================================================
 
 #[test]
-fn test_additive_schwarz_apply_fills_nan_on_solver_failure() {
-    // FailingLocalSolver always returns Err — apply must fill z with NAN.
+fn test_additive_schwarz_apply_returns_err_on_solver_failure() {
+    // FailingLocalSolver always returns Err — apply must propagate the error
+    // (previously this silently filled z with NaN).
     let n = 4;
     let solver = FailingLocalSolver {
         n_local: 2,
@@ -722,53 +671,12 @@ fn test_additive_schwarz_apply_fills_nan_on_solver_failure() {
 
     let rhs = vec![1.0; n];
     let mut z = vec![0.0; n];
-    schwarz.apply(&rhs, &mut z);
-
-    assert!(
-        z.iter().all(|v| v.is_nan()),
-        "all outputs should be NaN when solver fails, got: {:?}",
-        z,
-    );
-}
-
-// ============================================================================
-// try_apply / try_apply_adjoint direct calls match apply / apply_adjoint
-// ============================================================================
-
-#[test]
-fn test_additive_schwarz_try_apply_matches_apply() {
-    let n = 10;
-    let schwarz =
-        SchwarzPreconditioner::new(make_schwarz_entries(n), n).expect("valid additive schwarz");
-
-    let rhs: Vec<f64> = (0..n).map(|i| (i + 1) as f64).collect();
-
-    // apply path
-    let mut z_apply = vec![0.0; n];
-    schwarz.apply(&rhs, &mut z_apply);
-
-    // try_apply direct call
-    let mut z_try = vec![0.0; n];
-    let res = SchwarzPreconditioner::try_apply(&schwarz, &rhs, &mut z_try);
-    assert!(res.is_ok(), "try_apply should return Ok");
-
-    // try_apply_adjoint direct call (symmetric — same result)
-    let mut z_try_adj = vec![0.0; n];
-    let res_adj = Operator::try_apply_adjoint(&schwarz, &rhs, &mut z_try_adj);
-    assert!(res_adj.is_ok(), "try_apply_adjoint should return Ok");
+    let result = schwarz.apply(&rhs, &mut z);
 
-    for i in 0..n {
-        assert!(
-            (z_apply[i] - z_try[i]).abs() < 1e-14,
-            "try_apply mismatch at {i}: apply={}, try={}",
-            z_apply[i],
-            z_try[i],
-        );
-        assert!(
-            (z_apply[i] - z_try_adj[i]).abs() < 1e-14,
-            "try_apply_adjoint mismatch at {i}: apply={}, try_adj={}",
-            z_apply[i],
-            z_try_adj[i],
-        );
+    match result {
+        Err(SolveError::LocalSolveFailed { subdomain, .. }) => {
+            assert_eq!(subdomain, 0, "failure should be reported for subdomain 0");
+        }
+        other => panic!("expected LocalSolveFailed, got: {:?}", other),
     }
 }
diff --git a/crates/within-py/src/lib.rs b/crates/within-py/src/lib.rs
index 987891f..3c977cd 100644
--- a/crates/within-py/src/lib.rs
+++ b/crates/within-py/src/lib.rs
@@ -602,10 +602,11 @@ impl PyFePreconditioner {
             )));
         }
         let mut y = vec![0.0; self.inner.nrows()];
-        self.inner.apply(x_slice, &mut y);
+        self.inner
+            .apply(x_slice, &mut y)
+            .map_err(|e| pyo3::exceptions::PyRuntimeError::new_err(e.to_string()))?;
         Ok(numpy::PyArray1::from_vec(py, y))
     }
-
     /// Number of rows (DOFs).
     #[getter]
     fn nrows(&self) -> usize {
diff --git a/crates/within/benches/fixest.rs b/crates/within/benches/fixest.rs
index e5149a5..d70936f 100644
--- a/crates/within/benches/fixest.rs
+++ b/crates/within/benches/fixest.rs
@@ -285,7 +285,7 @@ fn bench_matvec(c: &mut Criterion) {
         let x: Vec<f64> = (0..n_dofs).map(|i| (i as f64).sin()).collect();
         let mut y = vec![0.0; n_obs];
         group.bench_function(BenchmarkId::new("apply", &label), |b| {
-            b.iter(|| op.apply(&x, &mut y))
+            b.iter(|| op.apply(&x, &mut y).expect("apply succeeds"))
         });
     }
     group.finish();
diff --git a/crates/within/src/operator.rs b/crates/within/src/operator.rs
index 9bc8a86..69ef810 100644
--- a/crates/within/src/operator.rs
+++ b/crates/within/src/operator.rs
@@ -94,7 +94,6 @@ impl<'a, S: ObservationStore> WeightedDesignOperator<'a, S> {
         }
     }
 }
-
 impl<S: ObservationStore> Operator for WeightedDesignOperator<'_, S> {
     fn nrows(&self) -> usize {
         self.design.n_rows
@@ -104,7 +103,7 @@ impl<S: ObservationStore> Operator for WeightedDesignOperator<'_, S> {
         self.design.n_dofs
     }
 
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), schwarz_precond::SolveError> {
         // y = W^{1/2} (D x)
         self.design.matvec_d(x, y);
         if let Some(sw) = &self.sqrt_weights {
@@ -112,9 +111,10 @@ impl<S: ObservationStore> Operator for WeightedDesignOperator<'_, S> {
                 *yi *= swi;
             }
         }
+        Ok(())
     }
 
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), schwarz_precond::SolveError> {
         // y = D^T (W^{1/2} x)
         match &self.sqrt_weights {
             None => self.design.rmatvec_dt(x, y),
@@ -126,5 +126,6 @@ impl<S: ObservationStore> Operator for WeightedDesignOperator<'_, S> {
                 self.design.rmatvec_dt(&tmp, y);
             }
         }
+        Ok(())
     }
 }
diff --git a/crates/within/src/operator/preconditioner.rs b/crates/within/src/operator/preconditioner.rs
index 93bd9e4..0f1cb4e 100644
--- a/crates/within/src/operator/preconditioner.rs
+++ b/crates/within/src/operator/preconditioner.rs
@@ -9,9 +9,9 @@
 //! # Integration with `schwarz-precond`
 //!
 //! The enum implements the [`Operator`] trait from the `schwarz-precond`
-//! crate, so it can be passed directly to LSMR as a preconditioner. Error
-//! handling flows through `try_apply` for graceful reporting of local-solver
-//! failures.
+//! crate, so it can be passed directly to LSMR as a preconditioner. The
+//! `apply` method is fallible — local-solver failures propagate to the
+//! caller as `SolveError`.
 
 use schwarz_precond::{LocalSolver, Operator, ReductionStrategy};
 use serde::{Deserialize, Serialize};
@@ -80,33 +80,17 @@ impl Operator for FePreconditioner {
         }
     }
 
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), schwarz_precond::SolveError> {
         match self {
             Self::Additive(p) => p.apply(x, y),
         }
     }
 
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), schwarz_precond::SolveError> {
         match self {
             Self::Additive(p) => p.apply_adjoint(x, y),
         }
     }
-
-    fn try_apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), schwarz_precond::ApplyError> {
-        match self {
-            Self::Additive(p) => p.try_apply(x, y),
-        }
-    }
-
-    fn try_apply_adjoint(
-        &self,
-        x: &[f64],
-        y: &mut [f64],
-    ) -> Result<(), schwarz_precond::ApplyError> {
-        match self {
-            Self::Additive(p) => p.try_apply_adjoint(x, y),
-        }
-    }
 }
 
 /// Build a [`FePreconditioner`] from a design and configuration.
diff --git a/crates/within/src/operator/schwarz.rs b/crates/within/src/operator/schwarz.rs
index 93c57e1..df1a8b9 100644
--- a/crates/within/src/operator/schwarz.rs
+++ b/crates/within/src/operator/schwarz.rs
@@ -65,8 +65,8 @@ impl FeSchwarz {
     }
 
     /// Apply the preconditioner, returning an error on local-solver failure.
-    pub fn try_apply(&self, r: &[f64], z: &mut [f64]) -> Result<(), schwarz_precond::ApplyError> {
-        self.0.try_apply(r, z)
+    pub fn apply(&self, r: &[f64], z: &mut [f64]) -> Result<(), schwarz_precond::SolveError> {
+        self.0.apply(r, z)
     }
 
     #[cfg(test)]
@@ -84,24 +84,12 @@ impl schwarz_precond::Operator for FeSchwarz {
         self.0.ncols()
     }
 
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
-        self.0.apply(x, y);
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), schwarz_precond::SolveError> {
+        self.0.apply(x, y)
     }
 
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-        self.0.apply_adjoint(x, y);
-    }
-
-    fn try_apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), schwarz_precond::ApplyError> {
-        self.0.try_apply(x, y)
-    }
-
-    fn try_apply_adjoint(
-        &self,
-        x: &[f64],
-        y: &mut [f64],
-    ) -> Result<(), schwarz_precond::ApplyError> {
-        self.0.try_apply_adjoint(x, y)
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), schwarz_precond::SolveError> {
+        self.0.apply_adjoint(x, y)
     }
 }
 
diff --git a/crates/within/src/operator/tests.rs b/crates/within/src/operator/tests.rs
index ea45285..45d3a1b 100644
--- a/crates/within/src/operator/tests.rs
+++ b/crates/within/src/operator/tests.rs
@@ -35,13 +35,13 @@ mod design_tests {
         let r = vec![0.1, 0.2, -0.3, 0.4, -0.5];
 
         let mut dx = vec![0.0; 5];
-        op.apply(&x, &mut dx);
+        op.apply(&x, &mut dx).expect("apply succeeds");
         let lhs: f64 = dx.iter().zip(r.iter()).map(|(a, b)| a * b).sum();
 
         let mut dtr = vec![0.0; 7];
-        op.apply_adjoint(&r, &mut dtr);
+        op.apply_adjoint(&r, &mut dtr)
+            .expect("apply_adjoint succeeds");
         let rhs: f64 = x.iter().zip(dtr.iter()).map(|(a, b)| a * b).sum();
-
         assert!((lhs - rhs).abs() < 1e-12);
     }
 
@@ -51,7 +51,7 @@ mod design_tests {
         let op = WeightedDesignOperator::new(&schema);
         let x = vec![1.0, 2.0, 3.0, 10.0, 20.0, 30.0, 40.0];
         let mut y = vec![0.0; 5];
-        op.apply(&x, &mut y);
+        op.apply(&x, &mut y).expect("apply succeeds");
         assert_eq!(y, vec![11.0, 22.0, 33.0, 41.0, 12.0]);
     }
 
@@ -61,7 +61,8 @@ mod design_tests {
         let op = WeightedDesignOperator::new(&schema);
         let r = vec![1.0, 2.0, 3.0, 4.0, 5.0];
         let mut x = vec![0.0; 7];
-        op.apply_adjoint(&r, &mut x);
+        op.apply_adjoint(&r, &mut x)
+            .expect("apply_adjoint succeeds");
         assert_eq!(x, vec![5.0, 7.0, 3.0, 6.0, 2.0, 3.0, 4.0]);
     }
 }
@@ -436,7 +437,7 @@ mod schwarz_tests {
     use crate::operator::schwarz::{
         build_additive_with_strategy, build_entry, build_reduced_schur_factor, ReducedSchurConfig,
     };
-    use schwarz_precond::{LocalSolver, Operator, ReductionStrategy};
+    use schwarz_precond::{LocalSolver, ReductionStrategy};
 
     const BLOCK_ELIM_NESTED_RAYON_CHILD_ENV: &str = "WITHIN_TEST_BLOCK_ELIM_NESTED_RAYON_CHILD";
 
@@ -616,9 +617,12 @@ mod schwarz_tests {
             for _ in 0..4 {
                 let mut z_reduction = vec![0.0; n_dofs];
                 let mut z_atomic = vec![0.0; n_dofs];
-                reduction.apply(&rhs, &mut z_reduction);
-                atomic.apply(&rhs, &mut z_atomic);
-
+                reduction
+                    .apply(&rhs, &mut z_reduction)
+                    .expect("reduction apply succeeds");
+                atomic
+                    .apply(&rhs, &mut z_atomic)
+                    .expect("atomic apply succeeds");
                 for (i, (&zr, &za)) in z_reduction.iter().zip(&z_atomic).enumerate() {
                     assert!(
                         zr.is_finite() && za.is_finite(),
@@ -761,8 +765,7 @@ mod schwarz_tests {
 
         let r = vec![1.0; design.n_dofs];
         let mut z = vec![0.0; design.n_dofs];
-        schwarz.apply(&r, &mut z);
-        assert!(z.iter().all(|&v| v.is_finite()));
+        schwarz.apply(&r, &mut z).expect("schwarz apply succeeds");
     }
 
     #[test]
diff --git a/crates/within/tests/error_paths.rs b/crates/within/tests/error_paths.rs
index 7817d36..bf4e81e 100644
--- a/crates/within/tests/error_paths.rs
+++ b/crates/within/tests/error_paths.rs
@@ -1,7 +1,7 @@
 use std::error::Error;
 
 use ndarray::Array2;
-use schwarz_precond::{ApplyError, PreconditionerBuildError, SolveError};
+use schwarz_precond::{PreconditionerBuildError, SolveError};
 use within::observation::{FactorMajorStore, ObservationWeights};
 use within::{solve, Preconditioner, SolverParams, WeightedDesign, WithinError};
 
@@ -121,7 +121,7 @@ fn test_within_error_display_preconditioner_build() {
 
 #[test]
 fn test_within_error_display_iterative_solve() {
-    let inner = SolveError::Apply(ApplyError::Synchronization { context: "test" });
+    let inner = SolveError::Synchronization { context: "test" };
     let e = WithinError::IterativeSolve(inner);
     assert!(e.to_string().contains("test"));
 }
@@ -168,7 +168,7 @@ fn test_within_error_source_preconditioner_build() {
 
 #[test]
 fn test_within_error_source_iterative_solve() {
-    let inner = SolveError::Apply(ApplyError::Synchronization { context: "test" });
+    let inner = SolveError::Synchronization { context: "test" };
     let e = WithinError::IterativeSolve(inner);
     assert!(e.source().is_some());
 }
@@ -194,7 +194,7 @@ fn test_within_error_from_preconditioner_build_error() {
 
 #[test]
 fn test_within_error_from_solve_error() {
-    let inner = SolveError::Apply(ApplyError::Synchronization { context: "test" });
+    let inner = SolveError::Synchronization { context: "test" };
     let e: WithinError = inner.into();
     match e {
         WithinError::IterativeSolve(_) => {}
diff --git a/crates/within/tests/properties.rs b/crates/within/tests/properties.rs
index b4a861b..593d2a7 100644
--- a/crates/within/tests/properties.rs
+++ b/crates/within/tests/properties.rs
@@ -61,8 +61,8 @@ proptest! {
 
         let mut y1 = vec![0.0; n];
         let mut y2 = vec![0.0; n];
-        fe_precond.apply(&x, &mut y1);
-        deserialized.apply(&x, &mut y2);
+        fe_precond.apply(&x, &mut y1).expect("apply succeeds");
+        deserialized.apply(&x, &mut y2).expect("deserialized apply succeeds");
 
         for (a, b) in y1.iter().zip(y2.iter()) {
             prop_assert!((a - b).abs() < 1e-12, "serde roundtrip mismatch: {} vs {}", a, b);