Merge pull request #106 from WIAS-PDELib/fix/store-qr

Improve handling of restriction data; compute nonlinear residual correctly

Author: pjaap

CHANGELOG.md

@@ -4,6 +4,11 @@
 
 ### Fixed
   - The linear residual takes also Lagrange multipliers and Lagrange restrictions into account.
+  - The nonlinear residual takes also Lagrange multipliers and Lagrange restrictions into account.
+
+### Added
+  - The compressed restrictions are now stored in a `CompressedRestriction`; for internal use only
+  - Each restriction stores the current Lagrange multipliers `:multipliers` vector in their parameters.
 
 ## v1.9.0
 

src/ExtendableFEM.jl

@@ -210,6 +210,7 @@ include("common_restrictions/linear_functional_restriction.jl")
 export LinearFunctionalRestriction
 export ZeroMeanValueRestriction
 export MassRestriction
+include("common_restrictions/compressed_restriction.jl")
 
 include("plots.jl")
 export plot_convergencehistory!

src/common_restrictions/boundarydata_restriction.jl

@@ -50,7 +50,7 @@ function assemble!(R::BoundaryDataRestriction, sol, SC; kwargs...)
     n = length(SC.b.entries)
     R.parameters[:matrix] = sparse(fixeddofs, 1:nvals, ones(nvals), n, nvals)
     R.parameters[:rhs] = fixedvals
-    R.parameters[:fixed_dofs] = fixeddofs
+    R.parameters[:multiplier] = zeros(nvals)
 
     return nothing
 end

src/common_restrictions/compressed_restriction.jl

@@ -0,0 +1,24 @@
+struct CompressedRestriction <: AbstractRestriction
+    parameters::Dict{Symbol, Any}
+end
+
+
+"""
+    CompressedRestriction(matrix::AbstractMatrix, rhs::AbstractVector)
+
+    Creates a simple linear restriction as a result of the QR compression of multiple restrictions
+    in the `solve` call. The compressed matrix and the compressed can be handed over directly in
+    the constructor.
+
+    This is 𝑛𝑜𝑡 exported, for internal use only!
+"""
+function CompressedRestriction(matrix::AbstractMatrix, rhs::AbstractVector)
+    return CompressedRestriction(
+        Dict{Symbol, Any}(
+            :name => "CompressedRestriction",
+            :matrix => matrix,
+            :rhs => rhs,
+            :multiplier => zero(rhs)
+        )
+    )
+end

src/common_restrictions/coupled_dofs_restriction.jl

@@ -134,9 +134,7 @@ function assemble!(R::CoupledDofsRestriction, sol, SC; kwargs...)
 
     R_orig.parameters[:matrix] = B
     R_orig.parameters[:rhs] = zeros(size(B, 2))
-
-    # fixed dofs are all active rows of B
-    R_orig.parameters[:fixed_dofs] = unique(findnz(B)[1])
+    R_orig.parameters[:multiplier] = zeros(size(B, 2))
 
     return nothing
 end

src/common_restrictions/linear_functional_restriction.jl

@@ -74,7 +74,7 @@ function assemble!(R::LinearFunctionalRestriction{T, Tv}, sol, SC; kwargs...) wh
 
     R.parameters[:matrix] = sparse(reshape(b.entries, n, 1))
     R.parameters[:rhs] = Tv[R.value]
-    R.parameters[:fixed_dofs] = []
+    R.parameters[:multiplier] = zeros(1)
 
     return nothing
 end

src/restrictions.jl

@@ -18,9 +18,9 @@ end
 
 function assemble!(R::AbstractRestriction, SC; kwargs...)
     ## assembles internal restriction matrix in R
+    @warn "assemble! not implemented for $(typeof(R))"
     return nothing
 end
 
 restriction_matrix(R::AbstractRestriction) = R.parameters[:matrix]
 restriction_rhs(R::AbstractRestriction) = R.parameters[:rhs]
-fixed_dofs(R::AbstractRestriction) = R.parameters[:fixed_dofs]

src/solvers.jl

@@ -38,16 +38,17 @@ Compute the nonlinear residual for the current solution.
 function compute_nonlinear_residual!(residual, A, b, sol, unknowns, PD, SC, freedofs)
     residual.entries .= b.entries
     for j in 1:length(b), k in 1:length(b)
-        addblock_matmul!(residual[j], A[j, k], sol[unknowns[k]], factor = -1)
+        addblock_matmul!(residual[j], A[j, k], sol[unknowns[k]], factor = -1.0)
     end
 
-    for op in PD.operators
-        residual.entries[fixed_dofs(op)] .= 0
+    # add Lagrange residuals
+    for rs in PD.restrictions
+        mul!(residual.entries, rs.parameters[:matrix], rs.parameters[:multiplier], -1.0, 1.0)
     end
 
-    mask_nonrestricted = ones(Bool, length(residual.entries))
-    for rs in PD.restrictions
-        mask_nonrestricted[fixed_dofs(rs)] .= false
+
+    for op in PD.operators
+        residual.entries[fixed_dofs(op)] .= 0
     end
 
     for u_off in SC.parameters[:inactive]
@@ -57,10 +58,22 @@ function compute_nonlinear_residual!(residual, A, b, sol, unknowns, PD, SC, free
         end
     end
 
-    nlres = length(freedofs) > 0 ? norm(residual.entries[freedofs]) : norm(residual.entries[mask_nonrestricted])
+    nlres = length(freedofs) > 0 ? norm(residual.entries[freedofs]) : norm(residual.entries)
+    restriction_residuals = [norm(rs.parameters[:matrix]' * sol.entries - rs.parameters[:rhs]) for rs in PD.restrictions]
 
-    if SC.parameters[:verbosity] > 0 && length(residual) > 1
-        @info "sub-residuals = $(norms(residual))"
+    if length(PD.restrictions) > 0
+        nlres = sqrt(nlres^2 + norm(restriction_residuals)^2)
+    end
+
+    for rs in PD.restrictions
+        mul!(residual.entries, rs.parameters[:matrix], rs.parameters[:multiplier], 1.0, 1.0)
+    end
+
+    if SC.parameters[:verbosity] > 0
+        if length(residual) > 1
+            @info "sub-residuals = $(norms(residual))"
+        end
+        @info "nonlinear residuals of restrictions = $restriction_residuals"
     end
 
     return nlres
@@ -244,8 +257,6 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
             else
                 A_unrestricted = A.entries.cscmatrix
             end
-            #else
-            #    A_unrestricted = A.entries.cscmatrix
         end
 
         # we solve for A Δx = r
@@ -267,7 +278,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
         else
             # add possible Lagrange restrictions
             restriction_matrices = [length(freedofs) > 0 ? view(restriction_matrix(re), freedofs, :) : restriction_matrix(re) for re in PD.restrictions ]
-            restriction_rhss = [length(freedofs) > 0 ? view(restriction_rhs(re), freedofs) : restriction_rhs(re) for re in PD.restrictions ]
+            restriction_rhss = deepcopy([length(freedofs) > 0 ? view(restriction_rhs(re), freedofs) : restriction_rhs(re) for re in PD.restrictions ])
 
             # block sizes for the block matrix
             block_sizes = [length(b_unrestricted), [ length(b) for b in restriction_rhss ]...]
@@ -278,11 +289,6 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
 
             @timeit timer "LM restrictions (nLMs = $nLMs)" begin
 
-                ## we need to add the (initial) solution to the rhs, since we work with the residual equation
-                for (B, rhs) in zip(restriction_matrices, restriction_rhss)
-                    rhs .-= B'sol_freedofs
-                end
-
                 Tv = eltype(b_unrestricted)
 
                 ## create block matrix
@@ -318,6 +324,9 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
                         restriction_matrices = [combined_restriction_matrix]
                         restriction_rhss = [combined_restriction_rhs]
 
+                        # remove previously defined restrictions (with explicit copy of the rhs; it is otherwise modified later)
+                        PD.restrictions = [CompressedRestriction(combined_restriction_matrix, deepcopy(combined_restriction_rhs))]
+
                         # update sizes
                         block_sizes = [length(b_unrestricted), length(combined_restriction_rhs)]
                         total_size = sum(block_sizes)
@@ -327,6 +336,11 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
                     A_block[Block(1, 1)] = A_unrestricted
                 end
 
+                ## we need to add the (initial) solution to the rhs, since we work with the residual equation
+                for (B, rhs) in zip(restriction_matrices, restriction_rhss)
+                    rhs .-= B'sol_freedofs
+                end
+
                 b_block = BlockVector(zeros(Tv, total_size), block_sizes)
 
                 b_block[Block(1)] = b_unrestricted
@@ -406,6 +420,11 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
             block_ends = cumsum(block_sizes)
             restriction_residuals = [norm(residual_flat[(block_ends[i] + 1):block_ends[i + 1]]) for i in 1:(length(block_sizes) - 1) ]
             push!(stats[:restriction_residuals], restriction_residuals)
+
+            # store Lagrange multipliers
+            for (i, rs) in enumerate(PD.restrictions)
+                rs.parameters[:multiplier] = linsolve.u[(block_ends[i] + 1):block_ends[i + 1]]
+            end
         end
 
         linres = norm(residual_flat)