fixed bug LinearOperatorDG, SolverStatistics replaced by Dict statist…

…ics in SolverConfiguration, now also tracks nnz of sparse matrix, small extensions of Example260, version bump

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
 authors = ["Christian Merdon <[email protected]>"]
-version = "0.0.16"
+version = "0.0.17"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"

examples/Example260_AxisymmetricNavierStokesProblem.jl

@@ -105,7 +105,7 @@ function main(; μ = 0.1, nrefs = 4, nonlinear = false, uniform = false, Plotter
     if nonlinear
         assign_operator!(PD, NonlinearOperator(kernel_convection!, [id(u)], [id(u),grad(u)]; bonus_quadorder = 1, kwargs...))#; jacobian = kernel_jacobian!)) 
     end
-    assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 1:2))
+    assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = [3]))
     assign_operator!(PD, HomogeneousBoundaryData(u; regions = [4], mask = (1,0,1)))
     assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1], mask = (0,1,1)))
 
@@ -114,7 +114,7 @@ function main(; μ = 0.1, nrefs = 4, nonlinear = false, uniform = false, Plotter
         xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs)
     else
         xgrid = simplexgrid(Triangulate;
-        points=[0 0 ; 1 0 ; 1 1 ; 0 1]',
+        points=[0 0 ; 5 0 ; 5 1 ; 0 1]',
         bfaces=[1 2 ; 2 3 ; 3 4 ; 4 1 ]',
         bfaceregions=[1, 2, 3, 4],
         regionpoints=[0.5 0.5;]',
@@ -123,14 +123,24 @@ function main(; μ = 0.1, nrefs = 4, nonlinear = false, uniform = false, Plotter
     end
 
     ## solve
-    FES = [FESpace{H1P2{2,2}}(xgrid), FESpace{H1P1{1}}(xgrid)]
+    FES = [FESpace{H1P2B{2,2}}(xgrid), FESpace{L2P1{1}}(xgrid)]
     sol = ExtendableFEM.solve(PD, FES; kwargs...)
 
     ## compute divergence in cylindrical coordinates by volume integrals
-    DivIntegrator = ItemIntegrator(kernel_l2div, [id(u), div(u)]; quadorder = 2, resultdim = 1)
+    DivIntegrator = ItemIntegrator(kernel_l2div, [id(u), div(u)]; quadorder = 3, resultdim = 1)
     div_error = sqrt(sum(evaluate(DivIntegrator, sol)))
-    @info "||div(u)|| = $div_error"
+    @info "||div(u_h)|| = $div_error"
 
+    ## compute L2error
+    function kernel_l2error(result, u_ops, qpinfo)
+        u!(result, qpinfo)
+        result .= (result - u_ops).^2
+    end
+    ErrorIntegratorExact = ItemIntegrator(kernel_l2error, [id(1)]; entities = ON_BFACES, regions = [3], quadorder = 4, kwargs...)
+    error = evaluate(ErrorIntegratorExact, sol)
+    L2error = sqrt(sum(view(error,1,:)) + sum(view(error,2,:)))
+    @info "||u - u_h|| = $L2error"
+
     ## plot
     plt = plot([id(u)], sol; Plotter = Plotter)
 

src/ExtendableFEM.jl

@@ -79,7 +79,7 @@ export get_periodic_coupling_info
 
 include("solver_config.jl")
 export SolverConfiguration
-export SolverStatistics, residual
+export residual
 
 include("solvers.jl")
 export solve

src/common_operators/linear_operator_dg.jl

@@ -680,11 +680,16 @@ function build_assembler!(b, O::LinearOperatorDG{Tv}, FE_test; time = 0.0) where
 			end
 		end
 		O.assembler = assembler
+	else
+		## update the time
+		for j = 1 : length(O.QP_infos)
+			O.QP_infos[j].time = time
+		end
 	end
 end
 
-function ExtendableFEM.assemble!(A, b, sol, O::LinearOperatorDG{Tv, UT}, SC::SolverConfiguration; assemble_matrix = true, kwargs...) where {Tv, UT}
-	if !assemble_matrix
+function ExtendableFEM.assemble!(A, b, sol, O::LinearOperatorDG{Tv, UT}, SC::SolverConfiguration; assemble_rhs = true, kwargs...) where {Tv, UT}
+	if !assemble_rhs
 		return nothing
 	end
 	if UT <: Integer
@@ -703,8 +708,8 @@ function ExtendableFEM.assemble!(A, b, sol, O::LinearOperatorDG{Tv, UT}, SC::Sol
 	end
 end
 
-function ExtendableFEM.assemble!(b, O::LinearOperatorDG{Tv, UT}, sol = nothing; assemble_matrix = true, kwargs...) where {Tv, UT}
-	if !assemble_matrix
+function ExtendableFEM.assemble!(b, O::LinearOperatorDG{Tv, UT}, sol = nothing; assemble_rhs = true, kwargs...) where {Tv, UT}
+	if !assemble_rhs
 		return nothing
 	end
 	@assert UT <: Integer

src/solver_config.jl

@@ -1,13 +1,14 @@
-mutable struct SolverStatistics{Tv}
-	assembly_times::Vector{Tv}
-	solver_times::Vector{Tv}
-	assembly_allocs::Vector{Tv}
-	solver_allocs::Vector{Tv}
-	linear_residuals::Vector{Tv}
-	nonlinear_residuals::Vector{Tv}
-end
-
-residual(S::SolverStatistics) = S.nonlinear_residuals[end]
+default_statistics() = Dict{Symbol, Vector{Real}}(
+	:assembly_times => [],
+	:solver_times => [],
+	:assembly_allocations => [],
+	:solver_allocations => [],
+	:linear_residuals => [],
+	:nonlinear_residuals => [],
+	:matrix_nnz => [],
+	:total_times => [],
+	:total_allocations => [],
+)
 
 mutable struct SolverConfiguration{AT <: AbstractMatrix, bT, xT}
 	PD::ProblemDescription
@@ -16,7 +17,7 @@ mutable struct SolverConfiguration{AT <: AbstractMatrix, bT, xT}
 	sol::xT					## stores solution
 	res::xT
 	LP::LinearProblem
-	statistics::SolverStatistics{Float64}
+	statistics::Dict{Symbol, Vector{Real}}
 	linsolver::Any
 	unknown_ids_in_sol::Array{Int, 1}
 	unknowns::Array{Unknown, 1}
@@ -25,6 +26,9 @@ mutable struct SolverConfiguration{AT <: AbstractMatrix, bT, xT}
 	parameters::Dict{Symbol, Any} # dictionary with user parameters
 end
 
+residual(S::SolverConfiguration) = S.statistics[:nonlinear_residuals][end]
+
+
 #
 # Default context information with help info.
 #
@@ -125,5 +129,5 @@ function SolverConfiguration(Problem::ProblemDescription, unknowns::Array{Unknow
 	end
 	res = deepcopy(b)
 	LP = LinearProblem(A.entries.cscmatrix, b.entries)
-	return SolverConfiguration{typeof(A), typeof(b), typeof(x)}(Problem, A, b, x, res, LP, SolverStatistics(Float64[],Float64[],Float64[],Float64[],Float64[],Float64[]), nothing, unknown_ids_in_sol, unknowns, copy(unknowns), offsets, parameters)
+	return SolverConfiguration{typeof(A), typeof(b), typeof(x)}(Problem, A, b, x, res, LP, default_statistics(), nothing, unknown_ids_in_sol, unknowns, copy(unknowns), offsets, parameters)
 end

src/solvers.jl

@@ -77,6 +77,9 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace,Vector{<
 	method_linear = SC.parameters[:method_linear]
 	precon_linear = SC.parameters[:precon_linear]
 	stats = SC.statistics
+	for (key, value) in stats
+		stats[key] = []
+	end
 
 	## unpack solver parameters
 	maxits = SC.parameters[:maxiterations]
@@ -245,10 +248,10 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace,Vector{<
 			time_final += time_assembly
 			allocs_final += allocs_assembly
 		end
-		push!(stats.assembly_allocs, allocs_assembly)
-		push!(stats.assembly_times, time_assembly)
+		push!(stats[:assembly_allocations], allocs_assembly)
+		push!(stats[:assembly_times], time_assembly)
 		if !is_linear
-			push!(stats.nonlinear_residuals, nlres)
+			push!(stats[:nonlinear_residuals], nlres)
 		end
 		if nlres < nltol
 			if SC.parameters[:verbosity] > -1
@@ -298,6 +301,7 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace,Vector{<
 				SC.parameters[:initialized] = true
 
 				## solve
+				push!(stats[:matrix_nnz], nnz(A.entries.cscmatrix))
 				x = LinearSolve.solve!(linsolve)
 
 				fill!(residual.entries, 0)
@@ -312,9 +316,9 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace,Vector{<
 				end
 				#@info residual.entries, norms(residual)
 				linres = norm(residual.entries)
-				push!(stats.linear_residuals, linres)
+				push!(stats[:linear_residuals], linres)
 				if is_linear
-					push!(stats.nonlinear_residuals, linres)
+					push!(stats[:nonlinear_residuals], linres)
 				end
 				offset = 0
 				for u in unknowns
@@ -331,8 +335,10 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace,Vector{<
 		time_total += time_solve
 		time_final += time_solve
 		allocs_final += allocs_solve
-		push!(stats.solver_allocs, allocs_solve)
-		push!(stats.solver_times, time_solve)
+		push!(stats[:solver_allocations], allocs_solve)
+		push!(stats[:solver_times], time_solve)
+		push!(stats[:total_times], time_total)
+		push!(stats[:total_allocations], (allocs_assembly + allocs_solve))
 		if SC.parameters[:verbosity] > -1
 			@printf "%.3e\t" linres
 			@printf "%.2f\t%.2f\t%.2f\t" time_assembly time_solve time_total