diff --git a/CHANGELOG.md b/CHANGELOG.md index 52f7f06b..45b22fad 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # CHANGES + +## v1.2.0 May 28, 2025 + +### Changed + - TimerOutputs for measuring/storing/showing runtime and allocations in solve, now also for separate operators + +### Fixed + - HomogeneousData/InterpolateBoundaryData operator fix when system matrix is of type GenericMTExtendableSparseMatrixCSC + ## v1.1.1 April 29, 2025 ### Fixed diff --git a/Project.toml b/Project.toml index fdd10daf..24224a06 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableFEM" uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed" authors = ["Christian Merdon ", "Patrick Jaap "] -version = "1.1.1" +version = "1.2" [deps] CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" @@ -19,6 +19,7 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] @@ -48,6 +49,7 @@ StaticArrays = "1.9.13" Symbolics = "4.2,5,6" Test = "1" TetGen = "2.0.0" +TimerOutputs = "0.5.29" Triangulate = "2.4.0" UnicodePlots = "3.6.5" julia = "1.9" diff --git a/examples/Example103_BurgersEquation.jl b/examples/Example103_BurgersEquation.jl index 5763181b..985c1870 100644 --- a/examples/Example103_BurgersEquation.jl +++ b/examples/Example103_BurgersEquation.jl @@ -96,8 +96,9 @@ function main(; t = 0 for it in 1:Int(floor(T / τ)) t += τ - ExtendableFEM.solve(PD, FES, SC; time = t) + ExtendableFEM.solve(PD, FES, SC; time = t, verbosity = -1, timeroutputs = :hide) end + show(timeroutputs(SC)) end ## plot final state diff --git a/examples/Example205_HeatEquation.jl b/examples/Example205_HeatEquation.jl index 03c240e5..78b49e85 100644 --- a/examples/Example205_HeatEquation.jl +++ b/examples/Example205_HeatEquation.jl @@ -79,8 +79,9 @@ function main(; t = 0 for it in 1:Int(floor(T / τ)) t += τ - ExtendableFEM.solve(PD, FES, SC; time = t) + ExtendableFEM.solve(PD, FES, SC; time = t, verbosity = -1, timeroutputs = :hide) end + @show timeroutputs(SC) end ## plot final state diff --git a/src/ExtendableFEM.jl b/src/ExtendableFEM.jl index 534f21ed..29592823 100644 --- a/src/ExtendableFEM.jl +++ b/src/ExtendableFEM.jl @@ -72,6 +72,7 @@ using SparseDiffTools: SparseDiffTools, ForwardColorJacCache, forwarddiff_color_jacobian!, matrix_colors using Symbolics: Symbolics using SciMLBase: SciMLBase +using TimerOutputs: TimerOutput, print_timer, @timeit using UnicodePlots: UnicodePlots @@ -145,6 +146,7 @@ export tensor_view include("solver_config.jl") export SolverConfiguration export residual +export timeroutputs include("solvers.jl") export solve diff --git a/src/common_operators/combinedofs.jl b/src/common_operators/combinedofs.jl index 0a8bec3c..415ef006 100644 --- a/src/common_operators/combinedofs.jl +++ b/src/common_operators/combinedofs.jl @@ -13,6 +13,7 @@ mutable struct CombineDofs{UT, CT} <: AbstractOperator end default_combop_kwargs() = Dict{Symbol, Tuple{Any, String}}( + :name => ("CombineDofs", "name for operator used in printouts"), :penalty => (1.0e30, "penalty for fixed degrees of freedom"), :verbosity => (0, "verbosity level"), ) diff --git a/src/solver_config.jl b/src/solver_config.jl index 7ee165a0..efc7f58b 100644 --- a/src/solver_config.jl +++ b/src/solver_config.jl @@ -1,13 +1,8 @@ -default_statistics() = Dict{Symbol, Vector{Real}}( - :assembly_times => [], - :solver_times => [], - :assembly_allocations => [], - :solver_allocations => [], - :linear_residuals => [], - :nonlinear_residuals => [], - :matrix_nnz => [], - :total_times => [], - :total_allocations => [], +default_statistics(Tv = Float64, Ti = Int64) = Dict{Symbol, Any}( + :timeroutputs => TimerOutput(), + :linear_residuals => Tv[], + :nonlinear_residuals => Tv[], + :matrix_nnz => Ti[], ) mutable struct SolverConfiguration{AT <: AbstractMatrix, bT, xT} @@ -19,7 +14,7 @@ mutable struct SolverConfiguration{AT <: AbstractMatrix, bT, xT} res::xT freedofs::Vector{Int} ## stores indices of free dofs LP::LinearProblem - statistics::Dict{Symbol, Vector{Real}} + statistics::Dict{Symbol, Any} linsolver::Any unknown_ids_in_sol::Array{Int, 1} unknowns::Array{Unknown, 1} @@ -28,17 +23,54 @@ mutable struct SolverConfiguration{AT <: AbstractMatrix, bT, xT} parameters::Dict{Symbol, Any} # dictionary with user parameters end +""" +```` +residuals(S::SolverConfiguration) +```` + +returns the vector with the residuals of all iterations +""" +residuals(S::SolverConfiguration) = S.statistics[:nonlinear_residuals] + """ ```` residual(S::SolverConfiguration) ```` returns the residual of the last solve - """ residual(S::SolverConfiguration) = S.statistics[:nonlinear_residuals][end] +""" +```` +timeroutputs(S::SolverConfiguration) +```` + +returns TimerOutputs object that contains detailed information on solving and assembly times +""" +timeroutputs(S::SolverConfiguration) = S.statistics[:timeroutputs] + + +""" +```` +lastmatrix(S::SolverConfiguration) +```` + +returns the currently stored system matrix +""" +lastmatrix(S::SolverConfiguration) = S.A + +""" +```` +lastrhs(S::SolverConfiguration) +```` + +returns the currently stored right-hand side +""" +lastrhs(S::SolverConfiguration) = S.b + + # # Default context information with help info. # @@ -65,6 +97,7 @@ default_solver_kwargs() = Dict{Symbol, Tuple{Any, String}}( :constant_rhs => (false, "right-hand side is constant (skips reassembly)"), :method_linear => (UMFPACKFactorization(), "any solver or custom LinearSolveFunction compatible with LinearSolve.jl (default = UMFPACKFactorization())"), :precon_linear => (nothing, "function that computes preconditioner for method_linear in case an iterative solver is chosen"), + :timeroutputs => (:full, "configures show of timeroutputs (choose between :hide, :full, :compact)"), :initialized => (false, "linear system in solver configuration is already assembled (turns true after first solve)"), :plot => (false, "plot all solved unknowns with a (very rough but fast) unicode plot"), ) @@ -83,8 +116,8 @@ end """ ```` -function iterate_until_stationarity( - SolverConfiguration(Problem::ProblemDescription +function SolverConfiguration( + Problem::ProblemDescription [FES::Union{<:FESpace, Vector{<:FESpace}}]; init = nothing, unknowns = Problem.unknowns, @@ -170,5 +203,5 @@ function SolverConfiguration(Problem::ProblemDescription, unknowns::Array{Unknow else LP = LinearProblem(A.entries.cscmatrix, b.entries) end - return SolverConfiguration{typeof(A), typeof(b), typeof(x)}(Problem, A, b, x, x_temp, res, freedofs, LP, default_statistics(), nothing, unknown_ids_in_sol, unknowns, copy(unknowns), offsets, parameters) + return SolverConfiguration{typeof(A), typeof(b), typeof(x)}(Problem, A, b, x, x_temp, res, freedofs, LP, default_statistics(TvM, TiM), nothing, unknown_ids_in_sol, unknowns, copy(unknowns), offsets, parameters) end diff --git a/src/solvers.jl b/src/solvers.jl index 497859b4..4ff95adc 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -70,19 +70,19 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ if SC.parameters[:verbosity] > 0 @info ".... reusing given solver configuration\n" end - time = 0 - allocs = 0 else - time = @elapsed begin - allocs = @allocated begin - SC = SolverConfiguration(PD, unknowns, FES; kwargs...) - if SC.parameters[:verbosity] > 0 - @info ".... init solver configuration\n" - end - end + SC = SolverConfiguration(PD, unknowns, FES; kwargs...) + if SC.parameters[:verbosity] > 0 + @info ".... init solver configuration\n" end end + ## load TimerOutputs + timer = timeroutputs(SC) + if timer == Any[] + timer = TimerOutput() + end + A = SC.A b = SC.b sol = SC.sol @@ -103,7 +103,6 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ damping = SC.parameters[:damping] freedofs = SC.freedofs - if SC.parameters[:verbosity] > -1 if length(freedofs) > 0 @info "SOLVING $(PD.name) @ time = $(SC.parameters[:time]) @@ -133,7 +132,7 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ end end end - if SC.parameters[:verbosity] > -1 + if SC.parameters[:verbosity] > 0 @info " nonlinear = $(nonlinear ? "true" : "false")\n" end if is_linear == "auto" @@ -146,30 +145,20 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ maxits = 0 end - alloc_factor = 1024^2 - if SC.parameters[:verbosity] > -1 - @printf " #IT\t------- RESIDUALS -------\t---- DURATION (s) ----\t\t---- ALLOCATIONS (MiB) ----\n" - @printf " \tNONLINEAR\tLINEAR\t\tASSEMB\tSOLVE\tTOTAL\t\tASSEMB\tSOLVE\tTOTAL\n" - @printf " INI\t\t\t\t\t\t\t%.2f\t\t\t\t%.2f\n" time allocs / alloc_factor + @printf " #IT\t------- RESIDUALS -------\n" + @printf " \tNONLINEAR\tLINEAR\n" end - time_final = time - allocs_final = allocs nlres = 1.1e30 linres = 1.1e30 linsolve = SC.linsolver reduced = false for j in 1:(maxits + 1) - allocs_assembly = 0 - time_assembly = 0 - time_solve_init = 0 - allocs_solve_init = 0 - time_total = 0 if is_linear && j == 2 nlres = linres else - time_total += @elapsed begin + @timeit timer "assembly" begin ## assemble operators if !SC.parameters[:constant_rhs] @@ -179,19 +168,19 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ fill!(A.entries.cscmatrix.nzval, 0) end if SC.parameters[:initialized] - time_assembly += @elapsed for op in PD.operators - allocs_assembly += @allocated assemble!(A, b, sol, op, SC; time = SC.parameters[:time], assemble_matrix = !SC.parameters[:constant_matrix], assemble_rhs = !SC.parameters[:constant_rhs], kwargs...) + for op in PD.operators + @timeit timer "$(op.parameters[:name])" assemble!(A, b, sol, op, SC; time = SC.parameters[:time], assemble_matrix = !SC.parameters[:constant_matrix], assemble_rhs = !SC.parameters[:constant_rhs], kwargs...) end else - time_assembly += @elapsed for op in PD.operators - allocs_assembly += @allocated assemble!(A, b, sol, op, SC; time = SC.parameters[:time], kwargs...) + for op in PD.operators + @timeit timer "$(op.parameters[:name]) (first)" assemble!(A, b, sol, op, SC; time = SC.parameters[:time], kwargs...) end end flush!(A.entries) ## penalize fixed dofs - time_assembly += @elapsed for op in PD.operators - allocs_assembly += @allocated apply_penalties!(A, b, sol, op, SC; assemble_matrix = !SC.parameters[:initialized] || !SC.parameters[:constant_matrix], assemble_rhs = !SC.parameters[:initialized] || !SC.parameters[:constant_rhs], kwargs...) + for op in PD.operators + @timeit timer "$(op.parameters[:name]) (penalties)" apply_penalties!(A, b, sol, op, SC; assemble_matrix = !SC.parameters[:initialized] || !SC.parameters[:constant_matrix], assemble_rhs = !SC.parameters[:initialized] || !SC.parameters[:constant_rhs], kwargs...) end flush!(A.entries) # end @@ -222,48 +211,50 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ # residual = copy(b) # end # end + end - ## show spy - if SC.parameters[:symmetrize] - A.entries.cscmatrix = (A.entries.cscmatrix + A.entries.cscmatrix') / 2 - elseif SC.parameters[:symmetrize_structure] - symmetrize_structure!(A.entries) - end - if SC.parameters[:show_matrix] - @show A - elseif SC.parameters[:spy] - @info ".... spy plot of system matrix:\n$(A.entries.cscmatrix))" - end - if SC.parameters[:check_matrix] - #λ, ϕ = Arpack.eigs(A.entries.cscmatrix; nev = 5, which = :SM, ritzvec = false) - #@info ".... 5 :SM eigs = $(λ)" - #λ, ϕ = Arpack.eigs(A.entries.cscmatrix; nev = 5, which = :LM, ritzvec = false) - #@info ".... 5 :LM eigs = $(λ)" - @info ".... ||A - A'|| = $(norm(A.entries.cscmatrix - A.entries.cscmatrix', Inf))" - @info ".... isposdef = $(isposdef(A.entries.cscmatrix))" - end + ## show spy + if SC.parameters[:symmetrize] + A.entries.cscmatrix = (A.entries.cscmatrix + A.entries.cscmatrix') / 2 + elseif SC.parameters[:symmetrize_structure] + symmetrize_structure!(A.entries) + end + if SC.parameters[:show_matrix] + @show A + elseif SC.parameters[:spy] + @info ".... spy plot of system matrix:\n$(A.entries.cscmatrix))" + end + if SC.parameters[:check_matrix] + #λ, ϕ = Arpack.eigs(A.entries.cscmatrix; nev = 5, which = :SM, ritzvec = false) + #@info ".... 5 :SM eigs = $(λ)" + #λ, ϕ = Arpack.eigs(A.entries.cscmatrix; nev = 5, which = :LM, ritzvec = false) + #@info ".... 5 :LM eigs = $(λ)" + @info ".... ||A - A'|| = $(norm(A.entries.cscmatrix - A.entries.cscmatrix', Inf))" + @info ".... isposdef = $(isposdef(A.entries.cscmatrix))" + end - ## init solver + ## init solver + @timeit timer "linear solver" begin if linsolve === nothing if SC.parameters[:verbosity] > 0 @info ".... initializing linear solver ($(method_linear))\n" end - time_solve_init += @elapsed begin - allocs_solve_init += @allocated begin - abstol = SC.parameters[:abstol] - reltol = SC.parameters[:reltol] - LP = reduced ? LP_reduced : SC.LP - if precon_linear !== nothing - linsolve = init(LP, method_linear; Pl = precon_linear(A.entries.cscmatrix), abstol = abstol, reltol = reltol) - else - linsolve = init(LP, method_linear; abstol = abstol, reltol = reltol) - end - SC.linsolver = linsolve + @timeit timer "initialization" begin + abstol = SC.parameters[:abstol] + reltol = SC.parameters[:reltol] + LP = reduced ? LP_reduced : SC.LP + if precon_linear !== nothing + linsolve = init(LP, method_linear; Pl = precon_linear(A.entries.cscmatrix), abstol = abstol, reltol = reltol) + else + linsolve = init(LP, method_linear; abstol = abstol, reltol = reltol) end + SC.linsolver = linsolve end end + end - ## compute nonlinear residual + ## compute nonlinear residual + @timeit timer "assembly" @timeit timer "residual vector" begin fill!(residual.entries, 0) for j in 1:length(b), k in 1:length(b) addblock_matmul!(residual[j], A[j, k], sol[unknowns[k]]) @@ -288,11 +279,7 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ @info "sub-residuals = $(norms(residual))" end end - time_final += time_assembly + time_solve_init - allocs_final += allocs_assembly + allocs_solve_init end - push!(stats[:assembly_allocations], allocs_assembly) - push!(stats[:assembly_times], time_assembly) if !is_linear push!(stats[:nonlinear_residuals], nlres) end @@ -300,33 +287,21 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ if SC.parameters[:verbosity] > -1 @printf " END\t" @printf "%.3e\t" nlres - @printf "\t\t%.2f\t\t%.2f\t" time_assembly time_total - @printf "\t%.2f\t\t%.2f\n" allocs_assembly / alloc_factor allocs_assembly / alloc_factor - @printf "\tconverged" - @printf "\t\t\t\tSUM -->\t%.2f" time_final - @printf "\t\t\tSUM -->\t%.2f\n\n" allocs_final / alloc_factor + @printf "converged\n" end break elseif isnan(nlres) if SC.parameters[:verbosity] > -1 @printf " END\t" @printf "%.3e\t" nlres - @printf "\t\t\t%.2f\t\t%.2f\t" time_assembly time_total - @printf "\t%.2f\t\t%.2f\n" allocs_assembly / alloc_factor allocs_assembly / alloc_factor - @printf "\tdid not converge" - @printf "\t\t\tSUM -->\t%.2f" time_final - @printf "\t\t\tSUM -->\t%.2f\n\n" allocs_final / alloc_factor + @printf "not converged\n" end break elseif (j == maxits + 1) && !(is_linear) if SC.parameters[:verbosity] > -1 @printf " END\t" @printf "\t\t%.3e\t" linres - @printf "\t\t%.2f\t\t%.2f\t" time_assembly time_total - @printf "\t%.2f\t\t%.2f\n" allocs_assembly / alloc_factor allocs_assembly / alloc_factor - @printf "\tmaxiterations reached" - @printf "\t\t\tSUM -->\t%.2f" time_final - @printf "\t\t\tSUM -->\t%.2f\n\n" allocs_final / alloc_factor + @printf "maxiterations reached\n" end break else @@ -344,30 +319,30 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ end end - time_solve = @elapsed begin - allocs_solve = @allocated begin - if !SC.parameters[:constant_matrix] || !SC.parameters[:initialized] - if length(freedofs) > 0 - linsolve.A = A.entries.cscmatrix[freedofs, freedofs] - else - linsolve.A = A.entries.cscmatrix - end - end - - # we solve for A Δx = r - # and update x = sol - Δx + @timeit timer "linear solver" begin + if !SC.parameters[:constant_matrix] || !SC.parameters[:initialized] if length(freedofs) > 0 - linsolve.b = residual.entries[freedofs] + linsolve.A = A.entries.cscmatrix[freedofs, freedofs] else - linsolve.b = residual.entries + linsolve.A = A.entries.cscmatrix end - SC.parameters[:initialized] = true + end - ## solve - push!(stats[:matrix_nnz], nnz(linsolve.A)) - Δx = LinearSolve.solve!(linsolve) + # we solve for A Δx = r + # and update x = sol - Δx + if length(freedofs) > 0 + linsolve.b = residual.entries[freedofs] + else + linsolve.b = residual.entries + end + SC.parameters[:initialized] = true + + ## solve + push!(stats[:matrix_nnz], nnz(linsolve.A)) + @timeit timer "solve! call" Δx = LinearSolve.solve!(linsolve) - # x = sol.entries - Δx.u for free dofs or partial solutions + # x = sol.entries - Δx.u for free dofs or partial solutions + @timeit timer "update solution" begin if length(freedofs) > 0 x = sol.entries[freedofs] - Δx.u else @@ -380,8 +355,10 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ offset += ndofs_u end end + end - ## check linear residual with full matrix + ## check linear residual with full matrix + @timeit timer "linear residual computation" begin if length(freedofs) > 0 soltemp.entries[freedofs] .= x residual.entries .= A.entries.cscmatrix * soltemp.entries @@ -396,13 +373,15 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ end end end - linres = norm(residual.entries) - push!(stats[:linear_residuals], linres) - if is_linear - push!(stats[:nonlinear_residuals], linres) - end + end + linres = norm(residual.entries) + push!(stats[:linear_residuals], linres) + if is_linear + push!(stats[:nonlinear_residuals], linres) + end - ## update solution (incl. damping etc.) + ## update solution (incl. damping etc.) + @timeit timer "update solution" begin offset = 0 if length(freedofs) > 0 sol.entries[freedofs] .= x @@ -419,23 +398,12 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ end end end - time_total += time_solve - time_final += time_solve - allocs_final += allocs_solve - time_solve += time_solve_init - allocs_solve += allocs_solve_init - 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 - @printf "\t%.2f\t%.2f\t%.2f\n" allocs_assembly / alloc_factor allocs_solve / alloc_factor (allocs_assembly + allocs_solve) / alloc_factor if is_linear - @printf "\tfinished" - @printf "\t\t\t\tSUM -->\t%.2f" time_final - @printf "\t\t\tSUM -->\t%.2f\n\n" allocs_final / alloc_factor + @printf "%.3e\t" linres + @printf "finished\n" + else + @printf "%.3e\n" linres end end end @@ -446,6 +414,14 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Union{<:FESpace, Vector{ end end + # Print the timings in the default way + stats[:timeroutputs] = timer + if SC.parameters[:timeroutputs] in [:full, :compact] + # TimerOutputs.complement!(to) + @printf "\n" + print_timer(timer, compact = SC.parameters[:timeroutputs] == :compact, sortby = :firstexec) + end + if SC.parameters[:return_config] return sol, SC else