fixed control flow and nonlinear residual computation in presence of restrictions (#96)

Author: chmerdon

CHANGELOG.md

@@ -1,5 +1,10 @@
 # CHANGES
 
+## v1.7.2
+
+### Fixed
+  - fixed control flow and nonlinear residual computation in presence of restrictions
+
 ## v1.7.1
 
 ### Fixed

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.7.1"
+version = "1.7.2"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]
@@ -14,7 +14,6 @@ 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"
@@ -42,7 +41,6 @@ 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"

src/ExtendableFEM.jl

@@ -64,7 +64,6 @@ 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!

src/common_restrictions/coupled_dofs_restriction.jl

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

src/solvers.jl

@@ -45,14 +45,19 @@ function compute_nonlinear_residual!(residual, A, b, sol, unknowns, PD, SC, free
         residual.entries[fixed_dofs(op)] .= 0
     end
 
+    mask_nonrestricted = ones(Bool, length(residual.entries))
+    for rs in PD.restrictions
+        mask_nonrestricted[fixed_dofs(rs)] .= false
+    end
+
     for u_off in SC.parameters[:inactive]
         j = get_unknown_id(SC, u_off)
         if j > 0
             fill!(residual[j], 0)
         end
     end
 
-    nlres = length(freedofs) > 0 ? norm(residual.entries[freedofs]) : norm(residual.entries)
+    nlres = length(freedofs) > 0 ? norm(residual.entries[freedofs]) : norm(residual.entries[mask_nonrestricted])
 
     if SC.parameters[:verbosity] > 0 && length(residual) > 1
         @info "sub-residuals = $(norms(residual))"
@@ -112,7 +117,7 @@ function assemble_system!(A, b, sol, PD, SC, timer; kwargs...)
         fill!(A.entries.cscmatrix.nzval, 0)
     end
 
-    # Assemble operators
+    # Assemble operators and restrictions
     if SC.parameters[:initialized]
         for op in PD.operators
             @timeit timer "$(op.parameters[:name])" begin
@@ -139,6 +144,10 @@ function assemble_system!(A, b, sol, PD, SC, timer; kwargs...)
             time_assembly += stats.time
             allocs_assembly += stats.bytes
         end
+        ## assemble restrictions
+        for restriction in PD.restrictions
+            @timeit timer "$(restriction.parameters[:name])" assemble!(restriction, sol, SC; kwargs...)
+        end
     end
     flush!(A.entries)
 
@@ -222,17 +231,9 @@ Solves the linear system and updates the solution vector. This includes:
 """
 function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns, freedofs, damping, PD, SC, stats, is_linear, timer, kwargs...)
 
-    @timeit timer "assembly" begin
-        ## assemble restrctions
-        if !SC.parameters[:initialized]
-            for restriction in PD.restrictions
-                @timeit timer "$(restriction.parameters[:name])" assemble!(restriction, sol, SC; kwargs...)
-            end
-        end
-    end
-
     @timeit timer "linear solver" begin
 
+
         # does the linsolve object need a (new) matrix?
         linsolve_needs_matrix = !SC.parameters[:constant_matrix] || !SC.parameters[:initialized]
 
@@ -243,6 +244,8 @@ 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
@@ -264,10 +267,10 @@ 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 = deepcopy([length(freedofs) > 0 ? view(restriction_rhs(re), freedofs) : restriction_rhs(re) for re in PD.restrictions ])
+            restriction_rhss = [length(freedofs) > 0 ? view(restriction_rhs(re), freedofs) : restriction_rhs(re) for re in PD.restrictions ]
 
             # block sizes for the block matrix
-            block_sizes = [size(A_unrestricted, 2), [ size(B, 2) for B in restriction_matrices ]...]
+            block_sizes = [length(b_unrestricted), [ length(b) for b in restriction_rhss ]...]
 
             # total number of additional LM dofs
             nLMs = @views sum(block_sizes[2:end])
@@ -280,7 +283,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
                 end
 
                 total_size = sum(block_sizes)
-                Tv = eltype(A_unrestricted)
+                Tv = eltype(b_unrestricted)
 
                 ## create block matrix
                 if linsolve_needs_matrix
@@ -506,6 +509,8 @@ function check_problem_linearity!(PDs, SCs, unknowns)
         end
         if SCs[j].parameters[:is_linear] == "auto"
             is_linear[j] = !nonlinear[j]
+        elseif SCs[j].parameters[:is_linear] == true
+            is_linear[j] = true
         end
         if is_linear[j] && nonlinear[j]
             @warn "problem $(PD.name) seems nonlinear, but user set is_linear = true (results may be wrong)!!"