From e495ccfa233b16b2797d05b6923cc3fb44ce13c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 28 Nov 2023 13:50:07 +0100 Subject: [PATCH] Add plothistory method Improve error messaging for time loop --- Project.toml | 2 +- docs/make.jl | 3 ++- docs/src/changes.md | 8 +++++++- docs/src/post.md | 1 + src/VoronoiFVM.jl | 7 ++++--- src/gridvisualize.jl | 34 ++++++++++++++++++++++++++++++++++ src/vfvm_exceptions.jl | 18 ++++++++++-------- src/vfvm_solver.jl | 14 ++++++++------ 8 files changed, 67 insertions(+), 20 deletions(-) diff --git a/Project.toml b/Project.toml index a8c125058..1759fbf04 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "VoronoiFVM" uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3" authors = ["Jürgen Fuhrmann ", "Dilara Abdel", "Jan Weidner", "Alexander Seiler", "Patricio Farrell", "Matthias Liero"] -version = "1.13.4" +version = "1.14.0" [deps] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" diff --git a/docs/make.jl b/docs/make.jl index 86c5bf8da..52f5d87ef 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,4 +1,5 @@ -using Documenter, ExampleJuggler, VoronoiFVM, ExtendableGrids +using Documenter, ExampleJuggler, VoronoiFVM +using ExtendableGrids, GridVisualize, LinearAlgebra function make_all(; with_examples = true, with_notebooks = true, example = nothing) ExampleJuggler.verbose!(true) diff --git a/docs/src/changes.md b/docs/src/changes.md index 4eaf79206..1d245917d 100644 --- a/docs/src/changes.md +++ b/docs/src/changes.md @@ -1,6 +1,12 @@ +# Changes +## v1.14.0 Nov 27, 2023 +- Add `Δu_max_factor` and `Δt_decrease` to [`VoronoiFVM.SolverControl`](@ref). Formerly (with values 2.0 and 0.5, respectively) the were fixed. New `Δu_max_factor=1.2` +- Add history to [`VoronoiFVM.TransientSolution`](@ref), prepare deprecation of `system.history` +- Add [`plothistory`](@ref) method + # Changes ## v1.13.0 July 24, 2023 -- Add `nodevolumes(sys)` method +- Add [`nodevolumes`](@ref) method ## v1.12.0 July 22, 2023 - Add functionality for outflow boundary conditions diff --git a/docs/src/post.md b/docs/src/post.md index cd7ec4ab8..4617842f4 100644 --- a/docs/src/post.md +++ b/docs/src/post.md @@ -10,6 +10,7 @@ GridVisualize.gridplot GridVisualize.gridplot! GridVisualize.scalarplot GridVisualize.scalarplot! +VoronoiFVM.plothistory ``` ## Norms & volumes diff --git a/src/VoronoiFVM.jl b/src/VoronoiFVM.jl index da703e45d..ba8b3e2a1 100644 --- a/src/VoronoiFVM.jl +++ b/src/VoronoiFVM.jl @@ -42,11 +42,12 @@ export dof export getdof export setdof! +include("vfvm_history.jl") +export NewtonSolverHistory, TransientSolverHistory, details + include("vfvm_transientsolution.jl") export transient_solution, TransientSolution -include("vfvm_history.jl") -export NewtonSolverHistory, TransientSolverHistory, details include("vfvm_xgrid.jl") export cartesian!, circular_symmetric!, spherical_symmetric! @@ -128,7 +129,7 @@ include("vfvm_diffeq_interface.jl") export eval_rhs!, eval_jacobian!, mass_matrix, prepare_diffeq! include("gridvisualize.jl") - +export plothistory #include("precompile.jl") end diff --git a/src/gridvisualize.jl b/src/gridvisualize.jl index e48b9b58e..74f2024da 100644 --- a/src/gridvisualize.jl +++ b/src/gridvisualize.jl @@ -76,3 +76,37 @@ function GridVisualize.scalarplot!(vis, sys::AbstractSystem, sol::AbstractTransi GridVisualize.scalarplot!(vis, simplexgrid(X, T), f; aspect = 1 / xtaspect, ylabel = tlabel, kwargs...) end end + +""" + plothistory(sys, tsol; + plots=[:timesteps,:newtonsteps,:updates], + size=(700,400), + Plotter=GridVisualize.default_plotter(), + kwargs...) + +Plot solution history stored in `tsol`. The `plots` argument allows to choose which plots are shown. +""" +function plothistory(sys, tsol; plots=[:timesteps,:newtonsteps,:updates], size=(700,400), Plotter=GridVisualize.default_plotter(), kwargs...) + hist=history(tsol) + t=hist.times + if length(t)==0 + error("Empty history. Did you pass `log=true` to the `solve()` method ?") + end + nplots=length(plots) + vis=GridVisualize.GridVisualizer(;layout=(nplots,1), size, Plotter, kwargs...) + + for iplot in eachindex(plots) + if plots[iplot]==:timesteps + GridVisualize.scalarplot!(vis[iplot,1],t[1:(end - 1)], t[2:end] - t[1:(end - 1)]; + title = "Time step sizes", xlabel = "t/s", ylabel = "Δt/s", color=:red, kwargs...) + elseif plots[iplot]==:newtonsteps + newtons=map(h->length(h.updatenorm),hist.histories) + GridVisualize.scalarplot!(vis[iplot,1],t, newtons,xlabel = "t/s", ylabel = "n", + title="Newton iterations", color=:red, limits=(0,maximum(newtons)+1), kwargs...) + elseif plots[iplot]==:updates + GridVisualize.scalarplot!(vis[iplot,1],t, hist.updates,xlabel = "t/s", ylabel = "du", + title="Time step updates", color=:red, kwargs...) + end + end + GridVisualize.reveal(vis) +end diff --git a/src/vfvm_exceptions.jl b/src/vfvm_exceptions.jl index f1eefb8ce..c8775e05e 100644 --- a/src/vfvm_exceptions.jl +++ b/src/vfvm_exceptions.jl @@ -32,22 +32,24 @@ end Print error when catching exceptions """ function _print_error(err, st) - println() - println(err) + io=IOBuffer() + println(io) + println(io,err) nlines = 5 for i = 1:min(nlines, length(st)) line = @sprintf("%s", st[i]) L = length(line) if L < 80 - println(line) + println(io,line) else - print(line[1:35]) - print(" ... ") - println(line[(L-35):L]) + print(io,line[1:35]) + print(io," ... ") + println(io,line[(L-35):L]) end end if length(st) > nlines - println("...") + println(io,"...") end - println() + println(io) + @warn String(take!(io)) end diff --git a/src/vfvm_solver.jl b/src/vfvm_solver.jl index e96c6394f..33cda8144 100644 --- a/src/vfvm_solver.jl +++ b/src/vfvm_solver.jl @@ -301,16 +301,16 @@ function CommonSolve.solve( if !called_from_API && doprint(control, 'd') @warn "[d]eprecated: solve(inival,system,times;kwargs...)" end + rd(x)=round(x,sigdigits=5) λstr = "t" if !transient λstr = "p" end - allhistory = TransientSolverHistory() solution = copy(inival) oldsolution = copy(inival) - _initialize_dirichlet!(inival, system; time, λ = Float64(lambdas[1]), params) + _initialize_dirichlet!(solution, system; time, λ = Float64(lambdas[1]), params) Δλ = Δλ_val(control, transient) t0 = @elapsed if !transient @@ -328,7 +328,7 @@ function CommonSolve.solve( kwargs..., ) control.post(solution, oldsolution, lambdas[1], 0) - if control.log + if control.log push!(allhistory, system.history) push!(allhistory.times, lambdas[1]) Δu = control.delta(system, solution, oldsolution, lambdas[1], 0) @@ -388,6 +388,9 @@ function CommonSolve.solve( ) end catch err + if transient + err="Problem at $(λstr)=$(λ|>rd), Δ$(λstr)=$(Δλ|>rd):\n$(err)" + end if (control.handle_exceptions) _print_error(err, stacktrace(catch_backtrace())) else @@ -404,12 +407,11 @@ function CommonSolve.solve( if !solved # reduce time step and retry solution Δλ = Δλ * control.Δt_decrease - rd(x)=round(x,sigdigits=5) if Δλ < Δλ_min(control, transient) if !(control.force_first_step && istep == 0) - err="At $(λstr)=$(λ|>rd): Δ$(λstr)_min=$(Δλ_min(control,transient)|>rd) reached while Δu=$(Δu|>rd) and Δu_opt=$(control.Δu_opt|>rd) " + err="At $(λstr)=$(λ|>rd): Δ$(λstr)_min=$(Δλ_min(control,transient)|>rd) reached while Δu=$(Δu|>rd) and Δu_opt=$(control.Δu_opt|>rd)\n Returning prematurely before $(λstr)=$(lambdas[end]|>rd) " if control.handle_exceptions - @warn err + @warn err else throw(DomainError(err)) end