From d37331b99ce62947e1cb3a1ffd5bbe99cef9e3d8 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Mon, 7 Apr 2025 15:52:42 +0200 Subject: [PATCH 1/3] Use LinearSolve routines for the residual equation rather than the original equation This enables us to use initial values for itearative solvers --- src/solvers.jl | 88 ++++++++++++++++++++++++++------------------------ 1 file changed, 46 insertions(+), 42 deletions(-) diff --git a/src/solvers.jl b/src/solvers.jl index 07981d7c..47272eb4 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -258,31 +258,29 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ end ## compute nonlinear residual - if !is_linear - fill!(residual.entries, 0) - for j in 1:length(b), k in 1:length(b) - addblock_matmul!(residual[j], A[j, k], sol[unknowns[k]]) - end - residual.entries .-= b.entries - #res = A.entries * sol.entries - b.entries - for op in PD.operators - residual.entries[fixed_dofs(op)] .= 0 - end - for u_off in SC.parameters[:inactive] - j = get_unknown_id(SC, u_off) - if j > 0 - fill!(residual[j], 0) - end - end - if length(freedofs) > 0 - nlres = norm(residual.entries[freedofs]) - else - nlres = norm(residual.entries) - end - if SC.parameters[:verbosity] > 0 && length(residual) > 1 - @info "sub-residuals = $(norms(residual))" + fill!(residual.entries, 0) + for j in 1:length(b), k in 1:length(b) + addblock_matmul!(residual[j], A[j, k], sol[unknowns[k]]) + end + residual.entries .-= b.entries + #res = A.entries * sol.entries - b.entries + for op in PD.operators + residual.entries[fixed_dofs(op)] .= 0 + end + for u_off in SC.parameters[:inactive] + j = get_unknown_id(SC, u_off) + if j > 0 + fill!(residual[j], 0) end end + if length(freedofs) > 0 + nlres = norm(residual.entries[freedofs]) + else + nlres = norm(residual.entries) + end + if SC.parameters[:verbosity] > 0 && length(residual) > 1 + @info "sub-residuals = $(norms(residual))" + end end time_final += time_assembly allocs_final += allocs_assembly @@ -349,25 +347,28 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ linsolve.A = A.entries.cscmatrix end end - if !SC.parameters[:constant_rhs] || !SC.parameters[:initialized] - if length(freedofs) > 0 - linsolve.b = b.entries[freedofs] - else - linsolve.b = b.entries - end + + # we solve for A Δx = r + # and update x = sol - Δx + if length(freedofs) > 0 + linsolve.b = residual.entries[freedofs] + else + linsolve.b = residual.entries end SC.parameters[:initialized] = true ## solve push!(stats[:matrix_nnz], nnz(linsolve.A)) - x = LinearSolve.solve!(linsolve) + Δx = LinearSolve.solve!(linsolve) + + x = sol.entries - Δx.u ## check linear residual with full matrix if length(freedofs) > 0 - soltemp.entries[freedofs] .= x.u + soltemp.entries[freedofs] .= x residual.entries .= A.entries.cscmatrix * soltemp.entries else - residual.entries .= A.entries.cscmatrix * x.u + residual.entries .= A.entries.cscmatrix * x end residual.entries .-= b.entries for op in PD.operators @@ -386,14 +387,14 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ ## update solution (incl. damping etc.) offset = 0 if length(freedofs) > 0 - sol.entries[freedofs] .= x.u + sol.entries[freedofs] .= x else for u in unknowns ndofs_u = length(view(sol[u])) if damping > 0 - view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x.u, (offset + 1):(offset + ndofs_u)) + view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x, (offset + 1):(offset + ndofs_u)) else - view(sol[u]) .= view(x.u, (offset + 1):(offset + ndofs_u)) + view(sol[u]) .= view(x, (offset + 1):(offset + ndofs_u)) end offset += ndofs_u end @@ -645,17 +646,20 @@ function iterate_until_stationarity( if !SC.parameters[:constant_matrix] || !SC.parameters[:initialized] linsolve.A = A.entries.cscmatrix end - if !SC.parameters[:constant_rhs] || !SC.parameters[:initialized] - linsolve.b = b.entries - end + + # we solve for A Δx = r + # and update x = sol - Δx + linsolve.b = residual.entries SC.parameters[:initialized] = true ## solve - x = LinearSolve.solve!(linsolve) + Δx = LinearSolve.solve!(linsolve) + + x = sol.entries - Δx.u fill!(residual.entries, 0) - mul!(residual.entries, A.entries.cscmatrix, x.u) + mul!(residual.entries, A.entries.cscmatrix, x) residual.entries .-= b.entries for op in PD.operators for dof in fixed_dofs(op) @@ -670,9 +674,9 @@ function iterate_until_stationarity( for u in unknowns[p] ndofs_u = length(view(sol[u])) if damping > 0 - view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x.u, (offset + 1):(offset + ndofs_u)) + view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x, (offset + 1):(offset + ndofs_u)) else - view(sol[u]) .= view(x.u, (offset + 1):(offset + ndofs_u)) + view(sol[u]) .= view(x, (offset + 1):(offset + ndofs_u)) end offset += ndofs_u end From cd4d3cdd740e7cba2ec09eded5250e46500c02bb Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Mon, 7 Apr 2025 15:57:37 +0200 Subject: [PATCH 2/3] Add CHANGELOG entry and bump compats --- CHANGELOG.md | 6 ++++++ Project.toml | 4 ++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 22396269..d4807f57 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,12 @@ ## next version +## v1.0.0 April 7, 2025 + +### Changed + + - `solve` uses now the residual equation for the linear systems. + ## v0.9.0 March 22, 2025 ### Added diff --git a/Project.toml b/Project.toml index ac50328f..b4a0aa63 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableFEM" uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed" authors = ["Christian Merdon ", "Patrick Jaap "] -version = "0.9.1" +version = "1.0.0" [deps] CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" @@ -28,7 +28,7 @@ DiffResults = "1" DocStringExtensions = "0.8,0.9" ExampleJuggler = "2.2.1" ExplicitImports = "1" -ExtendableFEMBase = "0.8" +ExtendableFEMBase = "1" ExtendableGrids = "1.10.3" ExtendableSparse = "1.5.3" ForwardDiff = "0.10.35,1" From ced3fc79479233f7ba9ec785395afcb0c42479f5 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Mon, 7 Apr 2025 16:33:02 +0200 Subject: [PATCH 3/3] Fix ranges of the solution vector --- src/solvers.jl | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/solvers.jl b/src/solvers.jl index 47272eb4..0bf2ea65 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -361,7 +361,19 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ push!(stats[:matrix_nnz], nnz(linsolve.A)) Δx = LinearSolve.solve!(linsolve) - x = sol.entries - Δx.u + # x = sol.entries - Δx.u for free dofs or partial solutions + if length(freedofs) > 0 + x = sol.entries[freedofs] - Δx.u + else + x = zero(Δx) + offset = 0 + for u in unknowns + ndofs_u = length(view(sol[u])) + x_range = (offset + 1):(offset + ndofs_u) + x[x_range] .= view(sol[u]) .- view(Δx, x_range) + offset += ndofs_u + end + end ## check linear residual with full matrix if length(freedofs) > 0 @@ -656,7 +668,15 @@ function iterate_until_stationarity( ## solve Δx = LinearSolve.solve!(linsolve) - x = sol.entries - Δx.u + # x = sol.entries - Δx.u ... in the entry ranges of the present unknowns + x = zero(Δx.u) + offset = 0 + for u in unknowns[p] + ndofs_u = length(view(sol[u])) + x_range = (offset + 1):(offset + ndofs_u) + x[x_range] .= view(sol[u]) .- view(Δx.u, x_range) + offset += ndofs_u + end fill!(residual.entries, 0) mul!(residual.entries, A.entries.cscmatrix, x)