Modifications for constraints for problems with multiple unknowns

Example285 now supports periodic boundary conditions

CompressedRestriction needs dummy assemble! function (otherwise crashes when linear problem is solved again with same ProblemDescription but without providing old SolverConfiguration), improved show of Unknowns

Author: chmerdon

CHANGELOG.md

@@ -1,5 +1,15 @@
 # CHANGES
 
+## v.1.10.3
+
+### Fixed
+  - Restriction now also work for problems with more than one unknown
+
+### Added
+  - Restrictions are listed in show of ProblemDescription
+  - CompressedRestriction has (empty) assemble! function to prevent crash in some cases
+  - improved show of Unknowns
+  
 ## v.1.10.0
 
 ### Changed

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.10.2"
+version = "1.10.3"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]

examples/Example285_CahnHilliard.jl

@@ -44,6 +44,7 @@ end
 function main(;
         order = 2,                              # finite element order for c and μ
         nref = 4,                               # refinement level
+        periodic = true,
         M = 1.0,
         λ = 1.0e-2,
         iterations_until_next_plot = 20,
@@ -68,6 +69,13 @@ function main(;
     assign_operator!(PD, BilinearOperator([id(μ)]; store = true))
     assign_operator!(PD, BilinearOperator([grad(μ)], [grad(c)]; factor = -λ, store = true))
 
+    if periodic
+        assign_restriction!(PD, CoupledDofsRestriction(c, 1, 3))
+        assign_restriction!(PD, CoupledDofsRestriction(c, 2, 4))
+        assign_restriction!(PD, CoupledDofsRestriction(μ, 1, 3))
+        assign_restriction!(PD, CoupledDofsRestriction(μ, 2, 4))
+    end
+
     ## add nonlinear reaction part (= -df/dc times test function)
     function kernel_dfdc!(result, input, qpinfo)
         return result[1] = -dfdc(input[1])

src/common_restrictions/boundarydata_restriction.jl

@@ -26,7 +26,12 @@ function BoundaryDataRestriction(u::Unknown, data = nothing; kwargs...)
     else
         boundary_operator = InterpolateBoundaryData(u, data; kwargs...)
     end
-    return BoundaryDataRestriction(u, boundary_operator, Dict{Symbol, Any}(:name => "BoundaryDataRestriction"))
+    return BoundaryDataRestriction(
+        u, boundary_operator, Dict{Symbol, Any}(
+            :name => "BoundaryDataRestriction",
+            :unknown => u
+        )
+    )
 end
 
 
@@ -47,7 +52,7 @@ function assemble!(R::BoundaryDataRestriction, sol, SC; kwargs...)
     nvals = length(fixeddofs)
 
     ## assign fixed dofs and vals to restriction
-    n = length(SC.b.entries)
+    n = length(SC.b[R.parameters[:unknown]])
     R.parameters[:matrix] = sparse(fixeddofs, 1:nvals, ones(nvals), n, nvals)
     R.parameters[:rhs] = fixedvals
     R.parameters[:multiplier] = zeros(nvals)

src/common_restrictions/compressed_restriction.jl

@@ -12,13 +12,20 @@ end
 
     This is 𝑛𝑜𝑡 exported, for internal use only!
 """
-function CompressedRestriction(matrix::AbstractMatrix, rhs::AbstractVector)
+function CompressedRestriction(unknown::Union{Unknown, Int}, matrix::AbstractMatrix, rhs::AbstractVector)
     return CompressedRestriction(
         Dict{Symbol, Any}(
             :name => "CompressedRestriction",
+            :unknown => unknown,
             :matrix => matrix,
             :rhs => rhs,
             :multiplier => zero(rhs)
         )
     )
 end
+
+
+function assemble!(R::CompressedRestriction, sol, SC; kwargs...)
+    # do nothing
+    return nothing
+end

src/common_restrictions/coupled_dofs_restriction.jl

@@ -15,11 +15,12 @@ end
 
     The matrix can be obtained from, e.g., `get_periodic_coupling_matrix`.
 """
-function CoupledDofsRestriction(matrix::AbstractMatrix)
+function CoupledDofsRestriction(matrix::AbstractMatrix; unknown = 1)
     return CoupledDofsRestriction(
         [matrix],
         Dict{Symbol, Any}(
             :name => "CoupledDofsRestriction",
+            :unknown => unknown,
             :reduce_col_space => false
         )
     )
@@ -34,11 +35,12 @@ end
     By default, the column space of the matrices is reduced to be of full rank.
     This can toggled by the `:reduce_col_space` parameter.
 """
-function CoupledDofsRestriction(matrices::Vector{AM}) where {AM <: AbstractMatrix}
+function CoupledDofsRestriction(matrices::Vector{AM}; unknown = 1) where {AM <: AbstractMatrix}
     return CoupledDofsRestriction(
         matrices,
         Dict{Symbol, Any}(
             :name => "CoupledDofsRestriction",
+            :unknown => unknown,
             :reduce_col_space => true
         )
     )
@@ -112,7 +114,7 @@ function assemble!(R::CoupledDofsRestriction, sol, SC; kwargs...)
         coupling_matrix = get_periodic_coupling_matrix(FES, source_region, target_region, give_opposite!; R_kwargs...)
 
         # replace R (in this scope)
-        R = CoupledDofsRestriction(coupling_matrix)
+        R = CoupledDofsRestriction(coupling_matrix, unknown = R.parameters[:unknown])
     end
 
 

src/common_restrictions/linear_functional_restriction.jl

@@ -15,8 +15,8 @@ Construct a linear functional restriction for a finite element unknown from a gi
 - `value::T`: The target value for the linear functional (default: `0`).
 
 """
-function LinearFunctionalRestriction(O::LinearOperator{Tv}; value::T = 0) where {Tv, T}
-    return LinearFunctionalRestriction{T, Tv}(value, O, Dict{Symbol, Any}(:name => "LinearFunctionalRestriction"))
+function LinearFunctionalRestriction(u::Unknown, O::LinearOperator{Tv}; value::T = 0) where {Tv, T}
+    return LinearFunctionalRestriction{T, Tv}(value, O, Dict{Symbol, Any}(:name => "LinearFunctionalRestriction", :unknown => u))
 end
 
 
@@ -37,7 +37,12 @@ Internally, it creates a `LinearOperator` using the provided kernel and operator
 """
 function ZeroMeanValueRestriction(u::Unknown; kernel = ExtendableFEMBase.constant_one_kernel, operator = Identity, Tv = Float64)
     linear_operator = LinearOperator(kernel, [(u, operator)]; store = true, Tv)
-    return LinearFunctionalRestriction{Int64, Tv}(0, linear_operator, Dict{Symbol, Any}(:name => "MeanValueRestriction"))
+    return LinearFunctionalRestriction{Int64, Tv}(
+        0, linear_operator, Dict{Symbol, Any}(
+            :name => "MeanValueRestriction",
+            :unknown => u
+        )
+    )
 end
 
 
@@ -59,20 +64,25 @@ Internally, it creates a `LinearOperator` using the provided kernel and operator
 """
 function MassRestriction(u::Unknown; kernel = ExtendableFEMBase.constant_one_kernel, value::T = 0, operator = Identity, Tv = Float64) where {T}
     linear_operator = LinearOperator(kernel, [(u, operator)]; store = true, Tv)
-    return LinearFunctionalRestriction{T, Tv}(value, linear_operator, Dict{Symbol, Any}(:name => "MeanValueRestriction"))
+    return LinearFunctionalRestriction{T, Tv}(
+        value, linear_operator, Dict{Symbol, Any}(
+            :name => "MeanValueRestriction",
+            :unknown => u
+        )
+    )
 end
 
 
 function assemble!(R::LinearFunctionalRestriction{T, Tv}, sol, SC; kwargs...) where {T, Tv}
 
-    n = length(SC.b.entries)
+    n = length(SC.b[R.parameters[:unknown]])
 
-    b = deepcopy(SC.b)
+    b = copy(SC.b)
     fill!(b.entries, 0.0)
 
     assemble!(nothing, b, sol, R.linear_operator, SC; assemble_rhs = true, kwargs...)
 
-    R.parameters[:matrix] = sparse(reshape(b.entries, n, 1))
+    R.parameters[:matrix] = sparse(reshape(view(b[R.parameters[:unknown]]), n, 1))
     R.parameters[:rhs] = Tv[R.value]
     R.parameters[:multiplier] = zeros(1)
 

src/problemdescription.jl

@@ -153,5 +153,13 @@ function Base.show(io::IO, PD::ProblemDescription)
             println(io, "    [$i] $o")
         end
     end
+    println(io, "  Restrictions ($(length(PD.restrictions))):")
+    if isempty(PD.restrictions)
+        println(io, "    (none)")
+    else
+        for (i, r) in enumerate(PD.restrictions)
+            println(io, "    [$i] $r")
+        end
+    end
     return
 end

src/restrictions.jl

@@ -7,7 +7,14 @@ abstract type AbstractRestriction end
 
 
 function Base.show(io::IO, R::AbstractRestriction)
-    print(io, "AbstractRestriction")
+    if haskey(R.parameters, :name)
+        print(io, "$(R.parameters[:name])")
+        if haskey(R.parameters, :name)
+            print(io, " for $(ansatz_function(R.parameters[:unknown]))")
+        end
+    else
+        print(io, "Restriction of type $(typeof(R))")
+    end
     return nothing
 end
 

src/solvers.jl

@@ -43,7 +43,7 @@ function compute_nonlinear_residual!(residual, A, b, sol, unknowns, PD, SC, free
 
     # add Lagrange residuals
     for rs in PD.restrictions
-        mul!(residual.entries, rs.parameters[:matrix], rs.parameters[:multiplier], -1.0, 1.0)
+        mul!(view(residual[rs.parameters[:unknown]]), rs.parameters[:matrix], rs.parameters[:multiplier], -1.0, 1.0)
     end
 
 
@@ -59,14 +59,14 @@ function compute_nonlinear_residual!(residual, A, b, sol, unknowns, PD, SC, free
     end
 
     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]
+    restriction_residuals = [norm(rs.parameters[:matrix]' * view(sol[rs.parameters[:unknown]]) - rs.parameters[:rhs]) for rs in PD.restrictions]
 
     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)
+        view(residual[rs.parameters[:unknown]]) .+= rs.parameters[:matrix] * rs.parameters[:multiplier]
     end
 
     if SC.parameters[:verbosity] > 0
@@ -278,9 +278,6 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
             b_unrestricted = residual.entries
         end
 
-        ## restrictions only involve the blocks coressponding to the unknowns and not necessarily the full sol.entries
-        sol_freedofs = vcat([view(sol[u]) for u in unknowns]...)
-
         if length(PD.restrictions) == 0
             if linsolve_needs_matrix
                 linsolve_A = A_unrestricted
@@ -306,68 +303,97 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
                 if linsolve_needs_matrix
                     if SC.parameters[:compress_restrictions]
 
+                        # compress only once and for all
+                        SC.parameters[:compress_restrictions] = false
+
                         SC.parameters[:verbosity] > 1 && @info "Compressing the restriction matrices..."
 
-                        # combine all restriction matrices into one:
-                        combined_transposed_restriction_matrix = vcat(restriction_matrices'...)
-                        combined_restriction_rhs = vcat(restriction_rhss...)
+                        compressed_restrictions = []
+                        compressed_matrices = []
+                        compressed_rhss = []
+
+                        # find all restricted unknowns
+                        restricted_unknowns = unique([ re.parameters[:unknown] for re in PD.restrictions])
+
+                        for u in restricted_unknowns
 
-                        # compute rank revealing QR decomposition (of the transposed matrix)
-                        qr_result = qr(combined_transposed_restriction_matrix)
+                            # which restrictions belong to `u`?
+                            indices = filter(i -> PD.restrictions[i].parameters[:unknown] == u, 1:length(PD.restrictions))
 
-                        SC.parameters[:verbosity] > 1 && @info "QR decomposition computed"
+                            # combine all restriction matrices into one:
+                            combined_transposed_restriction_matrix = vcat(restriction_matrices[indices]'...)
+                            combined_restriction_rhs = vcat(restriction_rhss[indices]...)
 
-                        # extract components (from docs: Q*R = combined_transposed_restriction_matrix[prow, pcol])
-                        (; Q, R, prow, pcol) = qr_result
+                            # compute rank revealing QR decomposition (of the transposed matrix)
+                            qr_result = qr(combined_transposed_restriction_matrix)
 
-                        # we need the inverse column permutation
-                        ipcol = invperm(pcol)
+                            SC.parameters[:verbosity] > 1 && @info "QR decomposition computed"
 
-                        # the new combined restriction block (add another transpose and convert into CSC matrix)
-                        combined_restriction_matrix = sparse(R[:, ipcol]')
+                            # extract components (from docs: Q*R = combined_transposed_restriction_matrix[prow, pcol])
+                            (; Q, R, prow, pcol) = qr_result
 
-                        # the new combined restriction rhs
-                        combined_restriction_rhs = Q'combined_restriction_rhs[prow]
+                            # we need the inverse column permutation
+                            ipcol = invperm(pcol)
 
-                        # compress the  column space
-                        qr_rank = rank(qr_result)
-                        SC.parameters[:verbosity] > 1 && @info "QR decomposition has rank $qr_rank"
-                        @assert norm(combined_restriction_rhs[(qr_rank + 1):end]) ≤ 1.0e-12 * norm(combined_restriction_rhs) "the rhs of the restriction is not in the image"
+                            # the new combined restriction block (add another transpose and convert into CSC matrix)
+                            combined_restriction_matrix = sparse(R[:, ipcol]')
 
-                        combined_restriction_matrix = combined_restriction_matrix[:, 1:qr_rank]
-                        combined_restriction_rhs = combined_restriction_rhs[1:qr_rank]
+                            # the new combined restriction rhs
+                            combined_restriction_rhs = Q'combined_restriction_rhs[prow]
 
-                        # replace by single entries
-                        restriction_matrices = [combined_restriction_matrix]
-                        restriction_rhss = [combined_restriction_rhs]
+                            # compress the  column space
+                            qr_rank = rank(qr_result)
+                            SC.parameters[:verbosity] > 1 && @info "QR decomposition has rank $qr_rank"
+                            @assert norm(combined_restriction_rhs[(qr_rank + 1):end]) ≤ 1.0e-12 * norm(combined_restriction_rhs) "the rhs of the restriction is not in the image"
 
-                        # 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))]
+                            combined_restriction_matrix = combined_restriction_matrix[:, 1:qr_rank]
+                            combined_restriction_rhs = combined_restriction_rhs[1:qr_rank]
+
+                            # replace by single entries
+                            push!(compressed_matrices, combined_restriction_matrix)
+                            push!(compressed_rhss, combined_restriction_rhs)
+
+                            # explicit copy of the rhs; it is otherwise modified later
+                            compressed_restriction = CompressedRestriction(u, combined_restriction_matrix, deepcopy(combined_restriction_rhs))
+                            push!(compressed_restrictions, compressed_restriction)
+
+                            SC.parameters[:verbosity] > 1 && @info "Created new CompressedRestriction for unknown $u"
+                        end
+
+                        # remove previously defined restrictions
+                        PD.restrictions = compressed_restrictions
+
+                        restriction_matrices = compressed_matrices
+                        restriction_rhss = compressed_rhss
 
                         # update sizes
-                        block_sizes = [length(b_unrestricted), length(combined_restriction_rhs)]
+                        block_sizes = [length(b_unrestricted), length.(compressed_rhss)...]
                         block_ends = cumsum(block_sizes)
-
-                        SC.parameters[:verbosity] > 1 && @info "Created new CompressedRestriction"
                     end
 
                     linsolve_A = spzeros(Tv, block_ends[end], block_ends[end])
                     linsolve_A[1:block_ends[1], 1:block_ends[1]] = A_unrestricted
                 end
 
+                # corresponding sol entries for each restriction
+                restricted_sol_entries = [view(sol[re.parameters[:unknown]]) for re in PD.restrictions]
+
                 ## 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
-                    mul!(rhs, B', sol_freedofs, -1.0, 1.0)
+                for (B, rhs, sol_entries) in zip(restriction_matrices, restriction_rhss, restricted_sol_entries)
+                    # rhs -= B'sol_entries
+                    mul!(rhs, B', sol_entries, -1.0, 1.0)
                 end
 
                 linsolve_b = zeros(Tv, block_ends[end])
                 linsolve_b[1:block_ends[1]] = b_unrestricted
 
                 for i in eachindex(restriction_matrices)
+                    # index range of the restricted unknown
+                    dof_range = (sol[PD.restrictions[i].parameters[:unknown]].offset + 1):(sol[PD.restrictions[i].parameters[:unknown]].last_index)
+
                     if linsolve_needs_matrix
-                        linsolve_A[1:block_ends[1], (block_ends[i] + 1):block_ends[i + 1]] = restriction_matrices[i]
-                        linsolve_A[(block_ends[i] + 1):block_ends[i + 1], 1:block_ends[1]] = restriction_matrices[i]'
+                        linsolve_A[dof_range, (block_ends[i] + 1):block_ends[i + 1]] = restriction_matrices[i]
+                        linsolve_A[(block_ends[i] + 1):block_ends[i + 1], dof_range] = restriction_matrices[i]'
                     end
                     linsolve_b[(block_ends[i] + 1):block_ends[i + 1]] = restriction_rhss[i]