Fix SCC residuals transfer for linear problems

Previously, when using the SCC algorithm with linear problems,
residuals were not computed and passed as nothing, resulting in
a vector with nothing values in the solution object.

This fix checks if sol.resid is nothing (which it often is for
LinearSolution) and computes residuals using the formula
`resid = A*u - b` when needed, then properly transfers them
to the solution object.

Fixes #706

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <[email protected]>

Author: ChrisRackauckas

Project.toml

@@ -23,6 +23,7 @@ NonlinearSolveSpectralMethods = "26075421-4e9a-44e1-8bd1-420ed7ad02b2"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SCCNonlinearSolve = "9dfe8606-65a1-4bb3-9748-cb89d1561431"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
@@ -43,23 +44,13 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
-[sources.BracketingNonlinearSolve]
-path = "lib/BracketingNonlinearSolve"
-
-[sources.NonlinearSolveBase]
-path = "lib/NonlinearSolveBase"
-
-[sources.NonlinearSolveFirstOrder]
-path = "lib/NonlinearSolveFirstOrder"
-
-[sources.NonlinearSolveQuasiNewton]
-path = "lib/NonlinearSolveQuasiNewton"
-
-[sources.NonlinearSolveSpectralMethods]
-path = "lib/NonlinearSolveSpectralMethods"
-
-[sources.SimpleNonlinearSolve]
-path = "lib/SimpleNonlinearSolve"
+[sources]
+BracketingNonlinearSolve = {path = "lib/BracketingNonlinearSolve"}
+NonlinearSolveBase = {path = "lib/NonlinearSolveBase"}
+NonlinearSolveFirstOrder = {path = "lib/NonlinearSolveFirstOrder"}
+NonlinearSolveQuasiNewton = {path = "lib/NonlinearSolveQuasiNewton"}
+NonlinearSolveSpectralMethods = {path = "lib/NonlinearSolveSpectralMethods"}
+SimpleNonlinearSolve = {path = "lib/SimpleNonlinearSolve"}
 
 [extensions]
 NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"

lib/SCCNonlinearSolve/src/SCCNonlinearSolve.jl

@@ -47,8 +47,10 @@ function iteratively_build_sols(alg, sols, (prob, explicitfun), args...; kwargs.
         # `remake` to recalculate `A` and `b` based on updated parameters from `explicitfun`.
         # Pass `A` and `b` to avoid unnecessarily copying them.
         sol = SciMLBase.solve(SciMLBase.remake(prob; A, b), alg.linalg; kwargs...)
+        # LinearSolution may have resid=nothing, so compute it: resid = A*u - b
+        resid = isnothing(sol.resid) ? A * sol.u - b : sol.resid
         SciMLBase.build_linear_solution(
-            alg.linalg, sol.u, nothing, nothing, retcode = sol.retcode)
+            alg.linalg, sol.u, resid, nothing, retcode = sol.retcode)
     else
         sol = SciMLBase.solve(prob, alg.nlalg; kwargs...)
         SciMLBase.build_solution(

lib/SCCNonlinearSolve/test/core_tests.jl

@@ -75,3 +75,75 @@ end
     scc_sol_default=solve(sccprob)
     @test sol ≈ manualscc ≈ scc_sol_default
 end
+
+@testitem "SCC Residuals Transfer" setup=[CoreRootfindTesting] tags=[:core] begin
+    using NonlinearSolveFirstOrder
+    using LinearAlgebra
+
+    # Create a simple SCC problem with both nonlinear and linear components
+    # to test that residuals are properly computed and transferred
+
+    # Nonlinear problem
+    function f1(du, u, p)
+        du[1] = u[1]^2 - 2.0
+        du[2] = u[2] - u[1]
+    end
+    explicitfun1(p, sols) = nothing
+    prob1 = NonlinearProblem(
+        NonlinearFunction{true, SciMLBase.NoSpecialize}(f1), [1.0, 1.0], nothing)
+
+    # Linear problem
+    A2 = [1.0 0.5; 0.5 1.0]
+    b2 = [1.0, 2.0]
+    prob2 = LinearProblem(A2, b2)
+    explicitfun2(p, sols) = nothing
+
+    # Another nonlinear problem
+    function f3(du, u, p)
+        du[1] = u[1] + u[2] - 3.0
+        du[2] = u[1] * u[2] - 2.0
+    end
+    explicitfun3(p, sols) = nothing
+    prob3 = NonlinearProblem(
+        NonlinearFunction{true, SciMLBase.NoSpecialize}(f3), [1.0, 2.0], nothing)
+
+    # Create SCC problem
+    sccprob = SciMLBase.SCCNonlinearProblem((prob1, prob2, prob3),
+        SciMLBase.Void{Any}.([explicitfun1, explicitfun2, explicitfun3]))
+
+    # Solve with SCCAlg
+    using SCCNonlinearSolve
+    scc_alg = SCCNonlinearSolve.SCCAlg(nlalg = NewtonRaphson(), linalg = nothing)
+    scc_sol = solve(sccprob, scc_alg)
+
+    # Test that solution was successful
+    @test SciMLBase.successful_retcode(scc_sol)
+
+    # Test that residuals are not nothing
+    @test scc_sol.resid !== nothing
+    @test !any(isnothing, scc_sol.resid)
+
+    # Test that residuals have the correct length
+    expected_length = length(prob1.u0) + length(prob2.b) + length(prob3.u0)
+    @test length(scc_sol.resid) == expected_length
+
+    # Test that residuals are small (near zero for converged solution)
+    @test norm(scc_sol.resid) < 1e-10
+
+    # Manually compute residuals to verify correctness
+    u1 = scc_sol.u[1:2]
+    u2 = scc_sol.u[3:4]
+    u3 = scc_sol.u[5:6]
+
+    # Compute residuals for each component
+    resid1 = zeros(2)
+    f1(resid1, u1, nothing)
+
+    resid2 = A2 * u2 - b2
+
+    resid3 = zeros(2)
+    f3(resid3, u3, nothing)
+
+    expected_resid = vcat(resid1, resid2, resid3)
+    @test scc_sol.resid ≈ expected_resid atol=1e-10
+end