Skip to content

Commit d37331b

Browse files
pjaappjaap
committed
Use LinearSolve routines for the residual equation rather than the original equation
This enables us to use initial values for itearative solvers
1 parent 988eb90 commit d37331b

File tree

1 file changed

+46
-42
lines changed

1 file changed

+46
-42
lines changed

‎src/solvers.jl‎

Lines changed: 46 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,28 @@ 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
364365

365366
## check linear residual with full matrix
366367
if length(freedofs) > 0
367-
soltemp.entries[freedofs] .= x.u
368+
soltemp.entries[freedofs] .= x
368369
residual.entries .= A.entries.cscmatrix * soltemp.entries
369370
else
370-
residual.entries .= A.entries.cscmatrix * x.u
371+
residual.entries .= A.entries.cscmatrix * x
371372
end
372373
residual.entries .-= b.entries
373374
for op in PD.operators
@@ -386,14 +387,14 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
386387
## update solution (incl. damping etc.)
387388
offset = 0
388389
if length(freedofs) > 0
389-
sol.entries[freedofs] .= x.u
390+
sol.entries[freedofs] .= x
390391
else
391392
for u in unknowns
392393
ndofs_u = length(view(sol[u]))
393394
if damping > 0
394-
view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x.u, (offset + 1):(offset + ndofs_u))
395+
view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x, (offset + 1):(offset + ndofs_u))
395396
else
396-
view(sol[u]) .= view(x.u, (offset + 1):(offset + ndofs_u))
397+
view(sol[u]) .= view(x, (offset + 1):(offset + ndofs_u))
397398
end
398399
offset += ndofs_u
399400
end
@@ -645,17 +646,20 @@ function iterate_until_stationarity(
645646
if !SC.parameters[:constant_matrix] || !SC.parameters[:initialized]
646647
linsolve.A = A.entries.cscmatrix
647648
end
648-
if !SC.parameters[:constant_rhs] || !SC.parameters[:initialized]
649-
linsolve.b = b.entries
650-
end
649+
650+
# we solve for A Δx = r
651+
# and update x = sol - Δx
652+
linsolve.b = residual.entries
651653
SC.parameters[:initialized] = true
652654

653655

654656
## solve
655-
x = LinearSolve.solve!(linsolve)
657+
Δx = LinearSolve.solve!(linsolve)
658+
659+
x = sol.entries - Δx.u
656660

657661
fill!(residual.entries, 0)
658-
mul!(residual.entries, A.entries.cscmatrix, x.u)
662+
mul!(residual.entries, A.entries.cscmatrix, x)
659663
residual.entries .-= b.entries
660664
for op in PD.operators
661665
for dof in fixed_dofs(op)
@@ -670,9 +674,9 @@ function iterate_until_stationarity(
670674
for u in unknowns[p]
671675
ndofs_u = length(view(sol[u]))
672676
if damping > 0
673-
view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x.u, (offset + 1):(offset + ndofs_u))
677+
view(sol[u]) .= damping * view(sol[u]) + (1 - damping) * view(x, (offset + 1):(offset + ndofs_u))
674678
else
675-
view(sol[u]) .= view(x.u, (offset + 1):(offset + ndofs_u))
679+
view(sol[u]) .= view(x, (offset + 1):(offset + ndofs_u))
676680
end
677681
offset += ndofs_u
678682
end

0 commit comments

Comments
 (0)