Add plothistory method

Improve error messaging for time loop
WIAS-PDELib · Nov 29, 2023 · e495ccf · e495ccf · j-fu · Nov 29, 2023
1 parent 9b31b96
commit e495ccf
Show file tree

Hide file tree

Showing 8 changed files with 67 additions and 20 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "VoronoiFVM"
 uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
 authors = ["Jürgen Fuhrmann <[email protected]>", "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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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