[documentation] Update the example for the Poisson equation with halo regions

Author: amontoison

docs/src/custom_workspaces.md

@@ -136,6 +136,12 @@ function LinearAlgebra.mul!(y::HaloVector{Float64}, A::LaplacianOperator, u::Hal
             
             # Update the output vector with the Laplacian result
             y.data[i,j] = dx2 + dy2
+
+            # Add a tiny value ε to the diagonal to ensure the operator is
+            # symmetric positive definite (SPD) and to avoid numerical issues
+            if i == j
+                y.data[i,j] += 1e-8
+            end
         end
     end
 
@@ -330,7 +336,7 @@ kc = KrylovConstructor(b)
 workspace = CgWorkspace(kc)
 
 # Solve the system with CG
-Krylov.cg!(workspace, A, b, atol=1e-12, rtol=0.0, verbose=1)
+Krylov.cg!(workspace, A, b, atol=1e-10, rtol=0.0, verbose=1)
 u_sol = solution(workspace)
 stats = statistics(workspace)
 ```

src/cg.jl

@@ -160,7 +160,7 @@ kwargs_cg = (:M, :ldiv, :radius, :linesearch, :atol, :rtol, :itmax, :timemax, :v
     MisI || mulorldiv!(z, M, r, ldiv)
     kcopy!(n, p, z)  # p ← z
     γ = kdotr(n, r, z)
-    γ ≥ 0 || error("The preconditioner `M` is not symmetric positive definite.")
+    γ ≥ 0 || error("The linear operator `A` or the preconditioner `M` is not symmetric positive definite.")
     rNorm = sqrt(γ)
     history && push!(rNorms, rNorm)
     if γ == 0
@@ -240,7 +240,7 @@ kwargs_cg = (:M, :ldiv, :radius, :linesearch, :atol, :rtol, :itmax, :timemax, :v
       kaxpy!(n, -α, Ap, r)
       MisI || mulorldiv!(z, M, r, ldiv)
       γ_next = kdotr(n, r, z)
-      γ_next ≥ 0 || error("The preconditioner `M` is not symmetric positive definite.")
+      γ_next ≥ 0 || error("The linear operator `A` or the preconditioner `M` is not symmetric positive definite.")
       rNorm = sqrt(γ_next)
       history && push!(rNorms, rNorm)