From 28cdfc5dc95f245074ee06b5875ff7bd9a88e4af Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Fri, 6 Feb 2026 18:11:42 +0100 Subject: [PATCH 1/2] 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. --- src/solver_config.jl | 9 +------ src/solvers.jl | 63 ++++++++++++++++++++++++-------------------- 2 files changed, 35 insertions(+), 37 deletions(-) diff --git a/src/solver_config.jl b/src/solver_config.jl index ec096e88..1977e424 100644 --- a/src/solver_config.jl +++ b/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 diff --git a/src/solvers.jl b/src/solvers.jl index 5b2a2503..22866f06 100644 --- a/src/solvers.jl +++ b/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 From d0aebce0b19b9f6e9964d1013738281f3bc40448 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Thu, 12 Feb 2026 13:25:36 +0100 Subject: [PATCH 2/2] CHANGELOG and bump to 1.10.0 --- CHANGELOG.md | 6 ++++++ Project.toml | 2 +- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c25c100c..84612f73 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # CHANGES +## v.1.10.0 + +### Changed + - `SolverConfiguration` has no `LP` property any more: everything is available via the `linsolver`. + - `initialize_linear_solver` is slightly changed and called much later in `solve`. + ## v.1.9.1 ### Fixed diff --git a/Project.toml b/Project.toml index 321553a1..8cd466f3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ExtendableFEM" uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed" -version = "1.9.1" +version = "1.10.0" authors = ["Christian Merdon ", "Patrick Jaap "] [deps]