Implement first demo of Lagrange Restrictions

Author: pjaap

Project.toml

@@ -4,13 +4,16 @@ version = "1.5.0"
 authors = ["Christian Merdon <[email protected]>", "Patrick Jaap <[email protected]>"]
 
 [deps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ExtendableFEMBase = "12fb9182-3d4c-4424-8fd1-727a0899810c"
 ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -26,6 +29,7 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
 Aqua = "0.8"
+BlockArrays = "1.7.0"
 ChunkSplitters = "3.1.2"
 CommonSolve = "0.2"
 DiffResults = "1"
@@ -35,6 +39,7 @@ ExplicitImports = "1"
 ExtendableFEMBase = "1.3.0"
 ExtendableGrids = "1.10.3"
 ExtendableSparse = "1.5.3"
+FillArrays = "1.13.0"
 ForwardDiff = "0.10.35,1"
 GridVisualize = "1.8.1"
 IncompleteLU = "0.2.1"

examples/Example212_PeriodicBoundary2D.jl

@@ -116,6 +116,7 @@ end
 function main(;
         order = 1,
         periodic = true,
+        use_LM_restrictions = true,
         Plotter = nothing,
         force = 10.0,
         h = 1.0e-2,
@@ -154,8 +155,12 @@ function main(;
             return nothing
         end
 
-        @time coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
-        assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        @showtime coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
+        if use_LM_restrictions
+            assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix))
+        else
+            assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        end
     end
 
     sol = solve(PD, FES)
@@ -172,10 +177,15 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example212_PeriodicBoundary2D, "example212.png") #hide
 function runtests()                                                                                  #hide
-    sol, _ = main()                                                                                  #hide
-    @test abs(maximum(view(sol[1])) - 1.3447465095618172) < 1.0e-3                                   #hide
-    sol2, _ = main(threads = 4)                                                                      #hide
-    @test sol.entries ≈ sol2.entries                                                                 #hide
+    sol1, _ = main(use_LM_restrictions = false, threads = 1)                                         #hide
+    @test abs(maximum(view(sol1[1])) - 1.3447465095618172) < 1.0e-3                                  #hide
+
+    sol2, _ = main(use_LM_restrictions = false, threads = 4)                                         #hide
+    @test sol1.entries ≈ sol2.entries                                                                #hide
+
+    sol3, _ = main(use_LM_restrictions = true, threads = 4)                                          #hide
+    @test sol1.entries ≈ sol3.entries                                                                #hide
+
     return nothing                                                                                   #hide
 end                                                                                                  #hide
 

examples/Example312_PeriodicBoundary3D.jl

@@ -130,6 +130,7 @@ end
 function main(;
         order = 1,
         periodic = true,
+        use_LM_restrictions = true,
         Plotter = nothing,
         force = 1.0,
         h = 1.0e-4,
@@ -169,8 +170,12 @@ function main(;
             y[1] = width - x[1]
             return nothing
         end
-        @time coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
-        assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        @showtime coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
+        if use_LM_restrictions
+            assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix))
+        else
+            assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
+        end
     end
 
     ## solve
@@ -185,10 +190,15 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicBoundary3D, "example312.png") #hide
 function runtests()                                                                                  #hide
-    sol, _ = main()                                                                                  #hide
-    @test abs(maximum(view(sol[1])) - 1.8004602502175202) < 2.0e-3                                   #hide
-    sol2, _ = main(threads = 4)                                                                      #hide
-    @test sol.entries ≈ sol2.entries                                                                                 #hide
+    sol1, _ = main(use_LM_restrictions = false, threads = 1)                                         #hide
+    @test abs(maximum(view(sol1[1])) - 1.8004602502175202) < 2.0e-3                                  #hide
+
+    sol2, _ = main(use_LM_restrictions = false, threads = 4)                                         #hide
+    @test sol1.entries ≈ sol2.entries                                                                #hide
+
+    sol3, _ = main(use_LM_restrictions = true, threads = 4)                                          #hide
+    @test sol1.entries ≈ sol3.entries                                                                #hide
+
     return nothing                                                                                   #hide
 end                                                                                                  #hide
 

src/ExtendableFEM.jl

@@ -5,6 +5,7 @@ $(read(joinpath(@__DIR__, "..", "README.md"), String))
 """
 module ExtendableFEM
 
+using BlockArrays: BlockMatrix, BlockVector, Block, BlockedMatrix, BlockedVector
 using ChunkSplitters: chunks
 using CommonSolve: CommonSolve
 using DiffResults: DiffResults
@@ -60,6 +61,7 @@ using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry,
 using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!,
     MTExtendableSparseMatrixCSC, findindex,
     rawupdateindex!
+using FillArrays: Zeros
 using ForwardDiff: ForwardDiff
 using GridVisualize: GridVisualize, GridVisualizer, gridplot!, reveal, save,
     scalarplot!, vectorplot!
@@ -68,8 +70,8 @@ using LinearSolve: LinearSolve, LinearProblem, UMFPACKFactorization, deleteat!,
     init, solve
 using Printf: Printf, @printf, @sprintf
 using SparseArrays: SparseArrays, AbstractSparseArray, SparseMatrixCSC, findnz, nnz,
-    nzrange, rowvals, sparse, SparseVector
-using StaticArrays: @MArray
+    nzrange, rowvals, sparse, SparseVector, spzeros
+    using StaticArrays: @MArray
 using SparseDiffTools: SparseDiffTools, ForwardColorJacCache,
     forwarddiff_color_jacobian!, matrix_colors
 using Symbolics: Symbolics
@@ -123,11 +125,16 @@ include("common_operators/reduction_operator.jl")
 #export AbstractReductionOperator
 #export FixbyInterpolation
 
+include("restrictions.jl")
+include("common_restrictions/coupled_dofs_restriction.jl")
+export CoupledDofsRestriction
+
 include("problemdescription.jl")
 export ProblemDescription
 export assign_unknown!
 export assign_operator!
 export replace_operator!
+export assign_restriction!
 
 include("helper_functions.jl")
 export get_periodic_coupling_info

src/common_restrictions/coupled_dofs_restriction.jl

@@ -0,0 +1,26 @@
+struct CoupledDofsRestriction{TM} <: AbstractRestriction
+    coupling_matrix::TM
+    parameters::Dict{Symbol, Any}
+end
+
+function CoupledDofsRestriction(matrix::AbstractMatrix)
+    return CoupledDofsRestriction(matrix, Dict{Symbol, Any}(:name => "CoupledDofsRestriction"))
+end
+
+
+function assemble!(R::CoupledDofsRestriction, A, b, sol, SC; kwargs...)
+
+    # extract all col indices
+    _, J, _ = findnz(R.coupling_matrix)
+
+    # remove duplicates
+    unique_cols = unique(J)
+
+    # subtract diagonal and shrink matrix to non-empty cols
+    B = (R.coupling_matrix - LinearAlgebra.I)[:, unique_cols]
+
+    R.parameters[:matrix] = B
+    return R.parameters[:rhs] = Zeros(length(unique_cols))
+
+    return nothing
+end

src/problemdescription.jl

@@ -26,6 +26,11 @@ mutable struct ProblemDescription
     """
     operators::Array{AbstractOperator, 1}
     #reduction_operators::Array{AbstractReductionOperator,1}
+
+    """
+    A vector of Lagrange restrictions that are involved in the problem.
+    """
+    restrictions::Vector{AbstractRestriction}
 end
 
 """
@@ -41,7 +46,7 @@ Create an empty `ProblemDescription` with the given name.
 
 """
 function ProblemDescription(name = "My problem")
-    return ProblemDescription(name, Array{Unknown, 1}(undef, 0), Array{AbstractOperator, 1}(undef, 0))
+    return ProblemDescription(name, [], [], [])
 end
 
 
@@ -92,6 +97,12 @@ function assign_operator!(PD::ProblemDescription, o::AbstractOperator)
 end
 
 
+function assign_restriction!(PD::ProblemDescription, r::AbstractRestriction)
+    push!(PD.restrictions, r)
+    return length(PD.restrictions)
+end
+
+
 """
 $(TYPEDSIGNATURES)
 

src/restrictions.jl

@@ -0,0 +1,22 @@
+"""
+	AbstractRestriction
+
+Root type for all restrictions
+"""
+abstract type AbstractRestriction end
+
+
+function Base.show(io::IO, R::AbstractRestriction)
+    print(io, "AbstractRestriction")
+    return nothing
+end
+
+# informs solver when operator needs reassembly in a time dependent setting
+function is_timedependent(R::AbstractRestriction)
+    return false
+end
+
+function assemble!(R::AbstractRestriction, A, b, sol, SC; kwargs...)
+    ## assembles internal restriction matrix in R
+    return nothing
+end

src/solvers.jl

@@ -221,27 +221,92 @@ Solves the linear system and updates the solution vector. This includes:
 - Computing the residual
 - Updating the solution with optional damping
 """
-function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns, freedofs, damping, PD, SC, stats, is_linear, timer)
-    # Update system matrix if needed
-    if !SC.parameters[:constant_matrix] || !SC.parameters[:initialized]
+function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns, freedofs, damping, PD, SC, stats, is_linear, timer, kwargs...)
+
+    @timeit timer "Lagrange restrictions" begin
+        ## assemble restrctions
+        if !SC.parameters[:initialized]
+            for restriction in PD.restrictions
+                @timeit timer "$(restriction.parameters[:name])" assemble!(restriction, A, b, sol, SC; kwargs...)
+            end
+        end
+    end
+
+    @timeit timer "linear solver" begin
+        ## start with the assembled matrix containing all assembled operators
+        if !SC.parameters[:constant_matrix] || !SC.parameters[:initialized]
+            if length(freedofs) > 0
+                A_unrestricted = A.entries.cscmatrix[freedofs, freedofs]
+            else
+                A_unrestricted = A.entries.cscmatrix
+            end
+        end
+
+        # we solve for A Δx = r
+        # and update x = sol - Δx
         if length(freedofs) > 0
-            linsolve.A = A.entries.cscmatrix[freedofs, freedofs]
+            b_unrestricted = residual.entries[freedofs]
         else
-            linsolve.A = A.entries.cscmatrix
+            b_unrestricted = residual.entries
         end
-    end
 
-    # Set right-hand side
-    if length(freedofs) > 0
-        linsolve.b = residual.entries[freedofs]
-    else
-        linsolve.b = residual.entries
-    end
-    SC.parameters[:initialized] = true
+        @timeit timer "LM restrictions" begin
+            ## add possible Lagrange restrictions
+            @timeit timer "prepare" begin
+                restriction_matrices = [length(freedofs) > 0 ? re.parameters[:matrix][freedofs, :] : re.parameters[:matrix] for re in PD.restrictions ]
+                restriction_rhs = [length(freedofs) > 0 ? re.parameters[:rhs][freedofs] : re.parameters[:rhs] 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_rhs)
+                    rhs .+= B'sol.entries
+                end
+            end
+
+
+            @timeit timer "compute blocks" begin
+                # block sizes for the block matrix
+                block_sizes = [size(A_unrestricted, 2), [ size(B, 2) for B in restriction_matrices ]...]
 
-    # Solve linear system
-    push!(stats[:matrix_nnz], nnz(linsolve.A))
-    @timeit timer "solve! call" Δx = LinearSolve.solve!(linsolve)
+                total_size = sum(block_sizes)
+                Tv = eltype(A_unrestricted)
+
+                ## create block matrix
+                A_block = BlockMatrix(spzeros(Tv, total_size, total_size), block_sizes, block_sizes)
+                A_block[Block(1, 1)] = A_unrestricted
+
+                b_block = BlockVector(zeros(Tv, total_size), block_sizes)
+                b_block[Block(1)] = b_unrestricted
+
+                u_unrestricted = linsolve.u
+                u_block = BlockVector(zeros(Tv, total_size), block_sizes)
+                u_block[Block(1)] = u_unrestricted
+
+                for i in eachindex(PD.restrictions)
+                    A_block[Block(1, i + 1)] = restriction_matrices[i]
+                    A_block[Block(i + 1, 1)] = transpose(restriction_matrices[i])
+                    b_block[Block(i + 1)] = restriction_rhs[i]
+
+                end
+            end
+
+            @timeit timer "convert" begin
+
+                linsolve.A = sparse(A_block) # convert to CSC Matrix
+                linsolve.b = Vector(b_block) # convert to dense vector
+                linsolve.u = Vector(u_block) # convert to dense vector
+
+            end
+        end
+
+        SC.parameters[:initialized] = true
+
+        # Solve linear system
+        push!(stats[:matrix_nnz], nnz(linsolve.A))
+        @timeit timer "solve! call" blocked_Δx = LinearSolve.solve!(linsolve)
+
+        # extract the solution / dismiss the lagrange multipliers
+        @views Δx = blocked_Δx[1:length(b_unrestricted)]
+    end
 
     # Compute solution update
     @timeit timer "update solution" begin
@@ -570,12 +635,12 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{
 
         linsolve = SC.linsolver
 
-        @timeit timer "linear solver" begin
-            linres = solve_linear_system!(
-                A, b, sol, soltemp, residual, linsolve, unknowns,
-                freedofs, damping, PD, SC, stats, is_linear, timer
-            )
-        end
+
+        linres = solve_linear_system!(
+            A, b, sol, soltemp, residual, linsolve, unknowns,
+            freedofs, damping, PD, SC, stats, is_linear, timer
+        )
+
         if SC.parameters[:verbosity] > -1
             if is_linear
                 @printf "%.3e\t" linres