Skip to content

Commit 8ace2e4

Browse files
pjaappjaap
authored
Merge pull request #46 from WIAS-PDELib/fix/solve-for-residual
Use LinearSolve routines for the residual equation rather than the original equation
2 parents 988eb90 + ced3fc7 commit 8ace2e4

File tree

3 files changed

+74
-44
lines changed

3 files changed

+74
-44
lines changed

‎CHANGELOG.md‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,12 @@
22

33
## next version
44

5+
## v1.0.0 April 7, 2025
6+
7+
### Changed
8+
9+
- `solve` uses now the residual equation for the linear systems.
10+
511
## v0.9.0 March 22, 2025
612

713
### Added

‎Project.toml‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "ExtendableFEM"
22
uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
33
authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
4-
version = "0.9.1"
4+
version = "1.0.0"
55

66
[deps]
77
CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
@@ -28,7 +28,7 @@ DiffResults = "1"
2828
DocStringExtensions = "0.8,0.9"
2929
ExampleJuggler = "2.2.1"
3030
ExplicitImports = "1"
31-
ExtendableFEMBase = "0.8"
31+
ExtendableFEMBase = "1"
3232
ExtendableGrids = "1.10.3"
3333
ExtendableSparse = "1.5.3"
3434
ForwardDiff = "0.10.35,1"

‎src/solvers.jl‎

Lines changed: 66 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -258,31 +258,29 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
258258
end
259259

260260
## compute nonlinear residual
261-
if !is_linear
262-
fill!(residual.entries, 0)
263-
for j in 1:length(b), k in 1:length(b)
264-
addblock_matmul!(residual[j], A[j, k], sol[unknowns[k]])
265-
end
266-
residual.entries .-= b.entries
267-
#res = A.entries * sol.entries - b.entries
268-
for op in PD.operators
269-
residual.entries[fixed_dofs(op)] .= 0
270-
end
271-
for u_off in SC.parameters[:inactive]
272-
j = get_unknown_id(SC, u_off)
273-
if j > 0
274-
fill!(residual[j], 0)
275-
end
276-
end
277-
if length(freedofs) > 0
278-
nlres = norm(residual.entries[freedofs])
279-
else
280-
nlres = norm(residual.entries)
281-
end
282-
if SC.parameters[:verbosity] > 0 && length(residual) > 1
283-
@info "sub-residuals = $(norms(residual))"
261+
fill!(residual.entries, 0)
262+
for j in 1:length(b), k in 1:length(b)
263+
addblock_matmul!(residual[j], A[j, k], sol[unknowns[k]])
264+
end
265+
residual.entries .-= b.entries
266+
#res = A.entries * sol.entries - b.entries
267+
for op in PD.operators
268+
residual.entries[fixed_dofs(op)] .= 0
269+
end
270+
for u_off in SC.parameters[:inactive]
271+
j = get_unknown_id(SC, u_off)
272+
if j > 0
273+
fill!(residual[j], 0)
284274
end
285275
end
276+
if length(freedofs) > 0
277+
nlres = norm(residual.entries[freedofs])
278+
else
279+
nlres = norm(residual.entries)
280+
end
281+
if SC.parameters[:verbosity] > 0 && length(residual) > 1
282+
@info "sub-residuals = $(norms(residual))"
283+
end
286284
end
287285
time_final += time_assembly
288286
allocs_final += allocs_assembly
@@ -349,25 +347,40 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
349347
linsolve.A = A.entries.cscmatrix
350348
end
351349
end
352-
if !SC.parameters[:constant_rhs] || !SC.parameters[:initialized]
353-
if length(freedofs) > 0
354-
linsolve.b = b.entries[freedofs]
355-
else
356-
linsolve.b = b.entries
357-
end
350+
351+
# we solve for A Δx = r
352+
# and update x = sol - Δx
353+
if length(freedofs) > 0
354+
linsolve.b = residual.entries[freedofs]
355+
else
356+
linsolve.b = residual.entries
358357
end
359358
SC.parameters[:initialized] = true
360359

361360
## solve
362361
push!(stats[:matrix_nnz], nnz(linsolve.A))
363-
x = LinearSolve.solve!(linsolve)
362+
Δx = LinearSolve.solve!(linsolve)
363+
364+
# x = sol.entries - Δx.u for free dofs or partial solutions
365+
if length(freedofs) > 0
366+
x = sol.entries[freedofs] - Δx.u
367+
else
368+
x = zero(Δx)
369+
offset = 0
370+
for u in unknowns
371+
ndofs_u = length(view(sol[u]))
372+
x_range = (offset + 1):(offset + ndofs_u)
373+
x[x_range] .= view(sol[u]) .- view(Δx, x_range)
374+
offset += ndofs_u
375+
end
376+
end
364377

365378
## check linear residual with full matrix
366379
if length(freedofs) > 0
367-
soltemp.entries[freedofs] .= x.u
380+
soltemp.entries[freedofs] .= x
368381
residual.entries .= A.entries.cscmatrix * soltemp.entries
369382
else
370-
residual.entries .= A.entries.cscmatrix * x.u
383+
residual.entries .= A.entries.cscmatrix * x
371384
end
372385
residual.entries .-= b.entries
373386
for op in PD.operators
@@ -386,14 +399,14 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
386399
## update solution (incl. damping etc.)
387400
offset = 0
388401
if length(freedofs) > 0
389-
sol.entries[freedofs] .= x.u
402+
sol.entries[freedofs] .= x
390403
else
391404
for u in unknowns
392405
ndofs_u = length(view(sol[u]))
393406
if damping > 0
394-
view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x.u, (offset + 1):(offset + ndofs_u))
407+
view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x, (offset + 1):(offset + ndofs_u))
395408
else
396-
view(sol[u]) .= view(x.u, (offset + 1):(offset + ndofs_u))
409+
view(sol[u]) .= view(x, (offset + 1):(offset + ndofs_u))
397410
end
398411
offset += ndofs_u
399412
end
@@ -645,17 +658,28 @@ function iterate_until_stationarity(
645658
if !SC.parameters[:constant_matrix] || !SC.parameters[:initialized]
646659
linsolve.A = A.entries.cscmatrix
647660
end
648-
if !SC.parameters[:constant_rhs] || !SC.parameters[:initialized]
649-
linsolve.b = b.entries
650-
end
661+
662+
# we solve for A Δx = r
663+
# and update x = sol - Δx
664+
linsolve.b = residual.entries
651665
SC.parameters[:initialized] = true
652666

653667

654668
## solve
655-
x = LinearSolve.solve!(linsolve)
669+
Δx = LinearSolve.solve!(linsolve)
670+
671+
# x = sol.entries - Δx.u ... in the entry ranges of the present unknowns
672+
x = zero(Δx.u)
673+
offset = 0
674+
for u in unknowns[p]
675+
ndofs_u = length(view(sol[u]))
676+
x_range = (offset + 1):(offset + ndofs_u)
677+
x[x_range] .= view(sol[u]) .- view(Δx.u, x_range)
678+
offset += ndofs_u
679+
end
656680

657681
fill!(residual.entries, 0)
658-
mul!(residual.entries, A.entries.cscmatrix, x.u)
682+
mul!(residual.entries, A.entries.cscmatrix, x)
659683
residual.entries .-= b.entries
660684
for op in PD.operators
661685
for dof in fixed_dofs(op)
@@ -670,9 +694,9 @@ function iterate_until_stationarity(
670694
for u in unknowns[p]
671695
ndofs_u = length(view(sol[u]))
672696
if damping > 0
673-
view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x.u, (offset + 1):(offset + ndofs_u))
697+
view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x, (offset + 1):(offset + ndofs_u))
674698
else
675-
view(sol[u]) .= view(x.u, (offset + 1):(offset + ndofs_u))
699+
view(sol[u]) .= view(x, (offset + 1):(offset + ndofs_u))
676700
end
677701
offset += ndofs_u
678702
end

0 commit comments

Comments
 (0)