Rework initialize_linear_solver

This removes the `LP` property of the `SolverConfiguration` since all information is
stored in the `linsolver`.

The initializer is now triggered much later in the `solve`-call: after system matrix
and right hand side are known. This is important when restrictions are in use.

Moreover, we do not need to assemble a `u` vector for the linear solver. This is done
in the init-call.

Author: pjaap

src/solver_config.jl

@@ -14,7 +14,6 @@ mutable struct SolverConfiguration{AT <: AbstractMatrix, bT, xT}
     tempsol::xT ## temporary solution
     res::xT
     freedofs::Vector{Int} ## stores indices of free dofs
-    LP::LinearProblem
     statistics::Dict{Symbol, Any}
     linsolver::Any
     unknown_ids_in_sol::Array{Int, 1}
@@ -221,11 +220,5 @@ function SolverConfiguration(Problem::ProblemDescription, unknowns::Array{Unknow
         x_temp = x
     end
 
-    ## construct linear problem
-    if length(freedofs) > 0
-        LP = LinearProblem(A.entries.cscmatrix[freedofs, freedofs], b.entries[freedofs])
-    else
-        LP = LinearProblem(A.entries.cscmatrix, b.entries)
-    end
-    return SolverConfiguration{typeof(A), typeof(b), typeof(x)}(Problem, A, b, x, x_temp, res, freedofs, LP, default_statistics(TvM, TiM), nothing, unknown_ids_in_sol, unknowns, copy(unknowns), offsets, parameters)
+    return SolverConfiguration{typeof(A), typeof(b), typeof(x)}(Problem, A, b, x, x_temp, res, freedofs, default_statistics(TvM, TiM), nothing, unknown_ids_in_sol, unknowns, copy(unknowns), offsets, parameters)
 end

src/solvers.jl

@@ -84,7 +84,7 @@ end
 
 Initialize the linear solver for the given system.
 """
-function init_linear_solver!(SC, A, timer, method_linear, precon_linear)
+function init_linear_solver!(SC, A, b, timer)
 
     # TODO use the timer
     time_assembly = 0.0
@@ -98,11 +98,14 @@ function init_linear_solver!(SC, A, timer, method_linear, precon_linear)
             stats = @timed begin
                 abstol = SC.parameters[:abstol]
                 reltol = SC.parameters[:reltol]
-                LP = SC.LP
+                method_linear = SC.parameters[:method_linear]
+                precon_linear = SC.parameters[:precon_linear]
+
+                LP = LinearProblem(A, b)
                 if precon_linear !== nothing
-                    SC.linsolver = init(LP, method_linear; Pl = precon_linear(A.entries.cscmatrix), abstol = abstol, reltol = reltol)
+                    SC.linsolver = init(LP, method_linear; Pl = precon_linear(A), abstol, reltol)
                 else
-                    SC.linsolver = init(LP, method_linear; abstol = abstol, reltol = reltol)
+                    SC.linsolver = init(LP, method_linear; abstol, reltol)
                 end
             end
             time_assembly += stats.time
@@ -242,11 +245,10 @@ Solves the linear system and updates the solution vector. This includes:
 - Computing the residual
 - Updating the solution with optional damping
 """
-function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns, freedofs, damping, PD, SC, stats, is_linear, timer, kwargs...)
+function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs, damping, PD, SC, stats, is_linear, timer, kwargs...)
 
     @timeit timer "linear solver" begin
 
-
         # does the linsolve object need a (new) matrix?
         linsolve_needs_matrix = !SC.parameters[:constant_matrix] || !SC.parameters[:initialized]
 
@@ -272,9 +274,9 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
 
         if length(PD.restrictions) == 0
             if linsolve_needs_matrix
-                linsolve.A = A_unrestricted
+                linsolve_A = A_unrestricted
             end
-            linsolve.b = b_unrestricted
+            linsolve_b = b_unrestricted
         else
             # add possible Lagrange restrictions
             restriction_matrices = [length(freedofs) > 0 ? view(restriction_matrix(re), freedofs, :) : restriction_matrix(re) for re in PD.restrictions ]
@@ -342,11 +344,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
                 end
 
                 b_block = BlockVector(zeros(Tv, total_size), block_sizes)
-
                 b_block[Block(1)] = b_unrestricted
-                u_unrestricted = @views linsolve.u[1:block_sizes[1]]
-                u_block = BlockVector(zeros(Tv, total_size), block_sizes)
-                u_block[Block(1)] = u_unrestricted
 
                 for i in eachindex(restriction_matrices)
                     if linsolve_needs_matrix
@@ -358,8 +356,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
                 end
 
                 # convert to dense vectors
-                linsolve.b = Vector(b_block)
-                linsolve.u = Vector(u_block)
+                linsolve_b = Vector(b_block)
 
                 if linsolve_needs_matrix
 
@@ -379,12 +376,28 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
                         # write each block directly in the resulting matrix
                         A_flat[range_row, range_col] = A_block[Block(i, j)]
                     end
-                    linsolve.A = A_flat
+                    linsolve_A = A_flat
                 end
             end
         end
 
-        SC.parameters[:initialized] = true
+
+        if SC.parameters[:initialized]
+            # set/update the linear system
+            SC.linsolver.b = linsolve_b
+            if linsolve_needs_matrix
+                SC.linsolver.A = linsolve_A
+            end
+        else
+            # init solver if not done before
+            @timeit timer "linear solver" begin
+                init_linear_solver!(SC, linsolve_A, linsolve_b, timer)
+            end
+            SC.parameters[:initialized] = true
+        end
+
+        # now the linear solver is definitely ready
+        linsolve = SC.linsolver
 
         # Solve linear system
         push!(stats[:matrix_nnz], nnz(linsolve.A))
@@ -725,11 +738,6 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
                 @info "....  isposdef  = $(isposdef(A.entries.cscmatrix))"
             end
 
-            ## init solver
-            @timeit timer "linear solver" begin
-                init_linear_solver!(SC, A, timer, method_linear, precon_linear)
-            end
-
             ## compute nonlinear residual
             @timeit timer "assembly" @timeit timer "residual vector" begin
                 nlres = compute_nonlinear_residual!(residual, A, b, sol, unknowns, PD, SC, freedofs)
@@ -744,11 +752,8 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
             break
         end
 
-        linsolve = SC.linsolver
-
-
         linres = solve_linear_system!(
-            A, b, sol, soltemp, residual, linsolve, unknowns,
+            A, b, sol, soltemp, residual, unknowns,
             freedofs, damping, PD, SC, stats, is_linear, timer
         )
 
@@ -874,11 +879,11 @@ function iterate_until_stationarity(
                     # Initialize linear solver if needed
                     solve_init_time, solve_init_allocs = init_linear_solver!(
                         SC,
-                        A,
-                        TimerOutput(),
-                        SC.parameters[:method_linear],
-                        SC.parameters[:precon_linear]
+                        A.entries.cscmatrix,
+                        b.entries,
+                        TimerOutput()
                     )
+
                     time_solve_init += solve_init_time
                     allocs_solve_init += solve_init_allocs