Fix timestep sizes at end of evolution

* clean + comment time evolution code

docs/src/changes.md

@@ -1,6 +1,12 @@
 # Changes
+## v1.15.0 Dec 1, 2023
+- Adjusted time/embedding stepping scheme, added `num_final_steps` to  [`VoronoiFVM.SolverControl`](@ref). This may lead to
+  sligthly different results when solving time dependent problems. Some unit test values have been adapted. Before,
+  accidentally very small time steps at the end of an evolution were possible.
+
 ## 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 `Δu_max_factor` and `Δt_decrease` to [`VoronoiFVM.SolverControl`](@ref). 
+  Before (with values 2.0 and 0.5, respectively) the were fixed. New `Δu_max_factor=1.2` default.
 - Add history to [`VoronoiFVM.TransientSolution`](@ref), prepare deprecation of `system.history`
 - Add [`plothistory`](@ref) method
 

examples/Example107_NonlinearStorage1D.jl

@@ -81,23 +81,24 @@ function main(; n = 20, m = 2.0, Plotter = nothing, verbose = false,
     control.force_first_step = true
     tsol = solve(sys; inival, times = [t0, tend], control)
 
-    p = GridVisualizer(; Plotter = Plotter, layout = (1, 1), fast = true)
-    for i = 1:length(tsol)
-        time = tsol.t[i]
-        scalarplot!(p[1, 1], grid, tsol[1, :, i]; title = @sprintf("t=%.3g", time),
-                    color = :red, label = "numerical")
-        scalarplot!(p[1, 1], grid, map(x -> barenblatt(x, time, m)^m, grid); clear = false,
-                    color = :green, label = "exact")
-        reveal(p)
-        sleep(1.0e-2)
+    if Plotter != nothing
+        p = GridVisualizer(; Plotter = Plotter, layout = (1, 1), fast = true)
+        for i = 1:length(tsol)
+            time = tsol.t[i]
+            scalarplot!(p[1, 1], grid, tsol[1, :, i]; title = @sprintf("t=%.3g", time),
+                        color = :red, label = "numerical")
+            scalarplot!(p[1, 1], grid, map(x -> barenblatt(x, time, m)^m, grid); clear = false,
+                        color = :green, label = "exact")
+            reveal(p)
+            sleep(1.0e-2)
+        end
     end
-
     return sum(tsol[end])
 end
 
 using Test
 function runtests()
-    testval = 175.20261258406686
+    testval = 175.2022110336759
     @test main(; unknown_storage = :sparse, assembly = :edgewise) ≈ testval &&
           main(; unknown_storage = :dense, assembly = :edgewise) ≈ testval &&
           main(; unknown_storage = :sparse, assembly = :cellwise) ≈ testval &&

examples/Example115_HeterogeneousCatalysis1D.jl

@@ -156,7 +156,7 @@ function main(; n = 10, Plotter = nothing, verbose = false, tend = 1,
     p = GridVisualizer(; Plotter = Plotter, layout = (3, 1))
 
     control = fixed_timesteps!(VoronoiFVM.NewtonControl(), tstep)
-    tsol = solve(sys; inival, times = [0, tend], control)
+    tsol = solve(sys; inival, times = [0, tend], control, verbose = verbose)
 
     p = GridVisualizer(; Plotter = Plotter, layout = (3, 1), fast = true)
     for it = 1:length(tsol)

examples/Example120_ThreeRegions1D.jl

@@ -7,10 +7,10 @@ using Printf
 using VoronoiFVM
 using ExtendableGrids
 using GridVisualize
-using LinearSolve
 
 function main(; n = 30, Plotter = nothing, plot_grid = false, verbose = false,
-              unknown_storage = :sparse, tend = 10, rely_on_corrections = false, assembly = :edgewise)
+              unknown_storage = :sparse, tend = 10,
+              rely_on_corrections = false, assembly = :edgewise)
     h = 3.0 / (n - 1)
     X = collect(0:h:3.0)
     grid = VoronoiFVM.Grid(X)
@@ -100,48 +100,38 @@ function main(; n = 30, Plotter = nothing, plot_grid = false, verbose = false,
 
     boundary_dirichlet!(sys, 3, 2, 0.0)
 
-    U = unknowns(sys)
-    U .= 0
-
-    control = VoronoiFVM.NewtonControl()
-    control.verbose = verbose
-    control.method_linear = SparspakFactorization()
-    tstep = 0.01
-    time = 0.0
-    istep = 0
     testval = 0
     p = GridVisualizer(; Plotter = Plotter, layout = (1, 1))
-    while time < tend
-        time = time + tstep
-        U = solve(sys; inival = U, control, tstep)
-        if verbose
-            @printf("time=%g\n", time)
-        end
-        tstep *= 1.1
-        istep = istep + 1
-        testval = U[2, 5]
 
+    testval = 0.0
+    function plot_timestep(U, Uold, time, Δt)
         U1 = view(U[1, :], subgrid1)
         U2 = view(U[2, :], subgrid2)
         U3 = view(U[3, :], subgrid3)
 
+        testval += sum(U2)
+
+        if Plotter == nothing
+            return
+        end
+
         scalarplot!(p[1, 1], subgrid1, U1; label = "spec1", color = (0.5, 0, 0),
                     xlimits = (0, 3), flimits = (0, 1e-3),
                     title = @sprintf("three regions t=%.3g", time))
         scalarplot!(p[1, 1], subgrid2, U2; label = "spec2", color = (0.0, 0.5, 0),
                     clear = false)
         scalarplot!(p[1, 1], subgrid3, U3; label = "spec3", color = (0.0, 0.0, 0.5),
                     clear = false, show = true)
-        if Plotter != nothing
-            sleep(1.0e-2)
-        end
     end
+
+    tsol = solve(sys; inival = 0, times = (0, tend), post = plot_timestep, verbose = verbose, Δu_opt = 1.0e-5)
+
     return testval
 end
 
 using Test
 function runtests()
-    testval = 0.0005967243505359461
+    testval = 0.359448515181824
     @test main(; unknown_storage = :sparse, rely_on_corrections = false, assembly = :edgewise) ≈ testval &&
           main(; unknown_storage = :dense, rely_on_corrections = false, assembly = :edgewise) ≈ testval &&
           main(; unknown_storage = :sparse, rely_on_corrections = true, assembly = :edgewise) ≈ testval &&

pluto-examples/api-update.jl

@@ -420,7 +420,7 @@ let
 end
 
 # ╔═╡ 4cbea340-9c02-4e69-8f5e-62bf45312bdd
-@test isapprox(sum(sol2) / length(sol2), 2.4908109508494247, rtol = 1.0e-14)
+@test isapprox(sum(sol2) / length(sol2), 2.4921650158811794, rtol = 1.0e-14)
 
 # ╔═╡ 1c18b5a0-cca6-46a1-bb9f-b3d65b8043c5
 md"""

pluto-examples/problemcase.jl

@@ -395,7 +395,7 @@ grid_n, sol_n, bt_n = trsolve(grid_2d(; nref = nref); tend = tend);
 sum(bt_n)
 
 # ╔═╡ c52ed973-2250-423a-b427-e91972f7ce74
-@test sum(bt_n) ≈ 17.643110936180495
+@test sum(bt_n) ≈ 18.143158169851787
 
 # ╔═╡ 732e79fa-5b81-4401-974f-37ea3427e770
 scalarplot(grid_n, sol_n(t)[ic, :]; resolution = (500, 200), show = true)
@@ -404,7 +404,7 @@ scalarplot(grid_n, sol_n(t)[ic, :]; resolution = (500, 200), show = true)
 grid_1, sol_1, bt_1 = trsolve(grid_1d(; nref = nref); tend = tend);
 
 # ╔═╡ 02330841-fdf9-4ebe-9da6-cf96529b223c
-@test sum(bt_1) ≈ 20.412099101959157
+@test sum(bt_1) ≈ 20.66209910195916
 
 # ╔═╡ e36d2aef-1b5a-45a7-9289-8d1e544bcedd
 scalarplot(grid_1,
@@ -419,7 +419,7 @@ scalarplot(grid_1,
 grid_f, sol_f, bt_f = trsolve(grid_2d(; nref = nref, ε_fix = ε_fix); tend = tend);
 
 # ╔═╡ d23d6634-266c-43e3-9493-b61fb390bbe7
-@test sum(bt_f) ≈ 20.411131554885404
+@test sum(bt_f) ≈ 20.661131375044135
 
 # ╔═╡ f42d4eb6-3e07-40c9-a8b3-dc772e674222
 scalarplot(grid_f, sol_f(t)[ic, :]; resolution = (500, 200), show = true)
@@ -428,7 +428,7 @@ scalarplot(grid_f, sol_f(t)[ic, :]; resolution = (500, 200), show = true)
 grid_ϕ, sol_ϕ, bt_ϕ = trsolve(grid_2d(; nref = nref); ϕ = [1.0e-3, 1], tend = tend);
 
 # ╔═╡ b260df8a-3721-4203-bc0c-a23bcab9a311
-@test sum(bt_ϕ) ≈ 20.4122562994476
+@test sum(bt_ϕ) ≈ 20.412256299447236
 
 # ╔═╡ ce49bb25-b2d0-4d17-a8fe-d7b62e9b20be
 begin