diff --git a/docs/src/changes.md b/docs/src/changes.md index 1d245917d..2ebd6e83d 100644 --- a/docs/src/changes.md +++ b/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 diff --git a/examples/Example107_NonlinearStorage1D.jl b/examples/Example107_NonlinearStorage1D.jl index 8d9ffe77a..855866755 100644 --- a/examples/Example107_NonlinearStorage1D.jl +++ b/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 && diff --git a/examples/Example115_HeterogeneousCatalysis1D.jl b/examples/Example115_HeterogeneousCatalysis1D.jl index 293083317..79d96b04c 100644 --- a/examples/Example115_HeterogeneousCatalysis1D.jl +++ b/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) diff --git a/examples/Example120_ThreeRegions1D.jl b/examples/Example120_ThreeRegions1D.jl index c0fabcd94..a9d35b604 100644 --- a/examples/Example120_ThreeRegions1D.jl +++ b/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,31 +100,21 @@ 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)) @@ -132,16 +122,16 @@ function main(; n = 30, Plotter = nothing, plot_grid = false, verbose = false, 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 && diff --git a/pluto-examples/api-update.jl b/pluto-examples/api-update.jl index 9918987e2..d6b603582 100644 --- a/pluto-examples/api-update.jl +++ b/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""" diff --git a/pluto-examples/problemcase.jl b/pluto-examples/problemcase.jl index 80149715b..a4cffd4da 100644 --- a/pluto-examples/problemcase.jl +++ b/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 diff --git a/src/vfvm_solver.jl b/src/vfvm_solver.jl index 33cda8144..e29a3a747 100644 --- a/src/vfvm_solver.jl +++ b/src/vfvm_solver.jl @@ -4,26 +4,21 @@ # is_precompiling() = ccall(:jl_generating_output, Cint, ()) == 1 - - - """ $(SIGNATURES) Solve time step problem. This is the core routine -for implicit Euler and stationary solve +for implicit Euler and stationary solve. """ -function _solve_timestep!( - solution::AbstractMatrix{Tv}, # old time step solution resp. initial value - oldsol::AbstractMatrix{Tv}, # old time step solution resp. initial value - system::AbstractSystem{Tv,Tc,Ti,Tm}, # Finite volume system - control::SolverControl, - time, - tstep, - embedparam, - params; - called_from_API = false, -) where {Tv,Tc,Ti,Tm} +function _solve_timestep!(solution::AbstractMatrix{Tv}, # old time step solution resp. initial value + oldsol::AbstractMatrix{Tv}, # old time step solution resp. initial value + system::AbstractSystem{Tv, Tc, Ti, Tm}, # Finite volume system + control::SolverControl, + time, + tstep, + embedparam, + params; + called_from_API = false,) where {Tv, Tc, Ti, Tm} _complete!(system; create_newtonvectors = true) nlhistory = NewtonSolverHistory() tasm = 0.0 @@ -71,17 +66,15 @@ function _solve_timestep!( while niter <= control.maxiters # Create Jacobi matrix and RHS for Newton iteration try - tasm += @elapsed nca, nba, nev = eval_and_assemble( - system, - solution, - oldsol, - residual, - time, - tstep, - embedparam, - params; - edge_cutoff = control.edge_cutoff, - ) + tasm += @elapsed nca, nba, nev = eval_and_assemble(system, + solution, + oldsol, + residual, + time, + tstep, + embedparam, + params; + edge_cutoff = control.edge_cutoff,) ncalloc += nca nballoc += nba neval += nev @@ -94,18 +87,16 @@ function _solve_timestep!( end end - tlinsolve += @elapsed _solve_linear!( - values(update), - system, - nlhistory, - control, - method_linear, - system.matrix, - values(residual), - ) + tlinsolve += @elapsed _solve_linear!(values(update), + system, + nlhistory, + control, + method_linear, + system.matrix, + values(residual)) values(solution) .-= damp * values(update) - + # "incremental collection may only sweep so-called young objects" GC.gc(false) @@ -142,16 +133,14 @@ function _solve_timestep!( itstring = @sprintf("it=% 3d", niter) end if control.max_round > 0 - @printf( - "%s %.3e %.3e %.3e % 2d\n", - itstring, - norm, - norm / oldnorm, - dnorm, - nround - ) + @printf("%s %.3e %.3e %.3e % 2d\n", + itstring, + norm, + norm/oldnorm, + dnorm, + nround) else - @printf("%s %.3e %.3e\n", itstring, norm, norm / oldnorm) + @printf("%s %.3e %.3e\n", itstring, norm, norm/oldnorm) end end if niter > 1 && norm / oldnorm > 1.0 / control.tol_mono @@ -187,15 +176,11 @@ function _solve_timestep!( end if doprint(control, 'n') && !system.is_linear - println( - " [n]ewton: $(round(t,sigdigits=3)) seconds asm: $(round(100*tasm/t,sigdigits=3))%, linsolve: $(round(100*tlinsolve/t,sigdigits=3))%", - ) + println(" [n]ewton: $(round(t,sigdigits=3)) seconds asm: $(round(100*tasm/t,sigdigits=3))%, linsolve: $(round(100*tlinsolve/t,sigdigits=3))%") end if doprint(control, 'l') && system.is_linear - println( - " [l]inear($(nameof(typeof(method_linear)))): $(round(t,sigdigits=3)) seconds", - ) + println(" [l]inear($(nameof(typeof(method_linear)))): $(round(t,sigdigits=3)) seconds") end system.history = nlhistory @@ -210,32 +195,28 @@ solve!(solution, inival, system; ```` Mutating version of [`solve(inival,system)`](@ref) """ -function VoronoiFVM.solve!( - solution, # Solution - inival, # Initial value - system::VoronoiFVM.AbstractSystem; # Finite volume system - control = SolverControl(), # Newton solver control information - time = Inf, - tstep = Inf, # Time step size. Inf means stationary solution - embedparam = 0.0, - params = zeros(0), - called_from_API = false, -) +function VoronoiFVM.solve!(solution, # Solution + inival, # Initial value + system::VoronoiFVM.AbstractSystem; # Finite volume system + control = SolverControl(), # Newton solver control information + time = Inf, + tstep = Inf, # Time step size. Inf means stationary solution + embedparam = 0.0, + params = zeros(0), + called_from_API = false,) fix_deprecations!(control) if !called_from_API && doprint(control, 'd') @warn "[d]eprecated: solve(inival,solution, system; kwargs...)" end - _solve_timestep!( - solution, - inival, - system, - control, - time, - tstep, - embedparam, - params; - called_from_API = true, - ) + _solve_timestep!(solution, + inival, + system, + control, + time, + tstep, + embedparam, + params; + called_from_API = true,) return solution end @@ -249,86 +230,92 @@ Alias for [`solve(system::VoronoiFVM.AbstractSystem; kwargs...)`](@ref) with the Solve stationary problem(if `tstep==Inf`) or one step implicit Euler step using Newton's method with `inival` as initial value. Returns a solution array. """ -function CommonSolve.solve( - inival, # Initial value - system::AbstractSystem; # Finite volume system - control = SolverControl(), # Newton solver control information - time = Inf, - tstep = Inf, # Time step size. Inf means stationary solution - params = zeros(0), - called_from_API = false, -) +function CommonSolve.solve(inival, # Initial value + system::AbstractSystem; # Finite volume system + control = SolverControl(), # Newton solver control information + time = Inf, + tstep = Inf, # Time step size. Inf means stationary solution + params = zeros(0), + called_from_API = false,) fix_deprecations!(control) if !called_from_API && doprint(control, 'd') @warn "[d]eprecated: solve(inival,system; kwargs...)" end - solve!( - unknowns(system), - inival, - system; - control = control, - time = time, - tstep = tstep, - params = params, - called_from_API = true, - ) + solve!(unknowns(system), + inival, + system; + control = control, + time = time, + tstep = tstep, + params = params, + called_from_API = true,) end -Δλ_val(control, transient) = transient ? control.Δt : control.Δp -Δλ_min(control, transient) = transient ? control.Δt_min : control.Δp_min -Δλ_max(control, transient) = transient ? control.Δt_max : control.Δp_max -Δλ_grow(control, transient) = transient ? control.Δt_grow : control.Δp_grow - """ solve(inival, system, times; kwargs...) Alias for [`solve(system::VoronoiFVM.AbstractSystem; kwargs...)`](@ref) with the corresponding keyword arguments. """ -function CommonSolve.solve( - inival, - system::VoronoiFVM.AbstractSystem, - lambdas; - control = SolverControl(), - transient = true, # choose between transient and stationary (embedding) case - time = 0.0, - params = zeros(0), - called_from_API = false, - kwargs..., -) +function CommonSolve.solve(inival, + system::VoronoiFVM.AbstractSystem, + lambdas; + control = SolverControl(), + transient = true, # choose between transient and stationary (embedding) case + time = 0.0, + params = zeros(0), + called_from_API = false, + kwargs...,) fix_deprecations!(control) 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" + # rounding in output + rd(x) = round(x; sigdigits = 5) + + # Set initial value of Δλ etc + if transient # λ is time + λstr = "t" + Δλ = control.Δt + Δλ_min = control.Δt_min + Δλ_max = control.Δt_max + Δλ_grow = control.Δt_grow + else # λ is embedding parameter + Δλ = control.Δp + Δλ_min = control.Δp_min + Δλ_max = control.Δp_max + Δλ_grow = control.Δp_grow end + Δu_opt = control.Δu_opt + Δu_max_factor = control.Δu_max_factor + Δt_decrease = control.Δt_decrease + allhistory = TransientSolverHistory() solution = copy(inival) - oldsolution = copy(inival) + oldsolution = copy(inival) # we need a copy as it is later overwritten + + # Initialize Dirichlet boundary values _initialize_dirichlet!(solution, system; time, λ = Float64(lambdas[1]), params) - Δλ = Δλ_val(control, transient) + # If not transient, solve for first embedding lambdas[1] t0 = @elapsed if !transient control.pre(solution, Float64(lambdas[1])) - solution = solve!( - solution, - oldsolution, - system; - called_from_API = true, - control = control, - time = time, - tstep = Inf, - embedparam = Float64(lambdas[1]), - params = params, - kwargs..., - ) + solution = solve!(solution, + oldsolution, + system; + called_from_API = true, + control = control, + time = time, + tstep = Inf, + embedparam = Float64(lambdas[1]), + params = params, + 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) @@ -337,6 +324,7 @@ function CommonSolve.solve( oldsolution .= solution end + # Initialize transient solution struct tsol = TransientSolution(Float64(lambdas[1]), solution; in_memory = control.in_memory) if doprint(control, 'e') @@ -346,50 +334,49 @@ function CommonSolve.solve( λ0 = 0 istep = 0 solved = false - t1 = @elapsed for i = 1:(length(lambdas)-1) - Δλ = max(Δλ, Δλ_min(control, transient)) + # Outer loop over embedding params/ time values + t1 = @elapsed for i = 1:(length(lambdas) - 1) + Δλ = max(Δλ, Δλ_min) λstart = lambdas[i] - λend = lambdas[i+1] + λend = lambdas[i + 1] λ = Float64(λstart) + # Inner loop between two embedding params/time values while λ < λend solved = false λ0 = λ Δu = 0.0 + + # Try to solve, possibly with stepsize decrease while !solved solved = true - try + try # check for non-converging newton λ = λ0 + Δλ control.pre(solution, λ) - if transient - solution = solve!( - solution, - oldsolution, - system; - called_from_API = true, - control = control, - time = λ, - tstep = Δλ, - params = params, - kwargs..., - ) - else - solution = solve!( - solution, - oldsolution, - system; - called_from_API = true, - control = control, - time = time, - tstep = Inf, - embedparam = λ, - params = params, - kwargs..., - ) + + if transient # λ is time + _time = λ + _tstep = Δλ + _embedparam = 0.0 + else # λ is embedding parameter + _time = time + _tstep = Inf + _embedparam = λ end + + solution = solve!(solution, + oldsolution, + system; + called_from_API = true, + control = control, + time = _time, + tstep = _tstep, + embedparam = _embedparam, + params = params, + kwargs...,) catch err if transient - err="Problem at $(λstr)=$(λ|>rd), Δ$(λstr)=$(Δλ|>rd):\n$(err)" + err = "Problem at $(λstr)=$(λ|>rd), Δ$(λstr)=$(Δλ|>rd):\n$(err)" end if (control.handle_exceptions) _print_error(err, stacktrace(catch_backtrace())) @@ -400,24 +387,30 @@ function CommonSolve.solve( end if solved Δu = control.delta(system, solution, oldsolution, λ, Δλ) - if Δu > control.Δu_max_factor * control.Δu_opt + if Δu > Δu_max_factor * Δu_opt solved = false end end if !solved - # reduce time step and retry solution - Δλ = Δλ * control.Δt_decrease - if Δλ < Δλ_min(control, transient) + # reduce time step + Δλ = Δλ * Δt_decrease + if Δλ < Δλ_min 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)\n Returning prematurely before $(λstr)=$(lambdas[end]|>rd) " + err = """ + At $(λstr)=$(λ|>rd): Δ$(λstr)_min=$(Δλ_min|>rd) reached while Δu=$(Δu|>rd) and Δu_opt=$(control.Δu_opt|>rd). + Returning prematurely before $(λstr)=$(lambdas[end]|>rd) + """ if control.handle_exceptions - @warn err + @warn err else throw(DomainError(err)) end - break + break # give up lowering stepsize, break out if "while !solved" loop else - solved=true + if doprint(control, 'e') + println("[e]volution: forced first timestep: Δ$(λstr)=$(Δλ), Δu=$(Δu|>rd), Δu_opt=$(control.Δu_opt|>rd)") + end + solved = true end end if doprint(control, 'e') @@ -426,18 +419,17 @@ function CommonSolve.solve( end end # while !solved - if solved + # Advance step + if solved istep = istep + 1 if doprint(control, 'e') - @printf( - "[e]volution: step=%d %s=%.3e Δ%s=%.3e Δu=%.3e\n", - istep, - λstr, - λ, - λstr, - Δλ, - Δu - ) + @printf("[e]volution: step=%d %s=%.3e Δ%s=%.3e Δu=%.3e\n", + istep, + λstr, + λ, + λstr, + Δλ, + Δu) end if control.log push!(allhistory, system.history) @@ -449,31 +441,47 @@ function CommonSolve.solve( end control.post(solution, oldsolution, λ, Δλ) oldsolution .= solution + + # Adjust last timesteps, so there will be no accidental + # very small timestep + steps_to_go = ceil((λend - λ) / Δλ) + λ_predict = λend - λ + if steps_to_go < control.num_final_steps + λ_predict = (λend - λ) / steps_to_go + end + + # see fixed_timesteps!() + if Δλ_max ≈ Δλ_min + λ_predict = Δλ_max + end + if λ < λend - Δλ = min( - Δλ_max(control, transient), - Δλ * Δλ_grow(control, transient), - Δλ * control.Δu_opt / (Δu + 1.0e-14), - λend - λ, - ) + # ### account for close last timestep + Δλ = min(Δλ_max, + Δλ * Δλ_grow, + Δλ * Δu_opt / (Δu + 1.0e-14), + λ_predict, + λend - λ) end - else - break + else + break # break out of inner loop overt timestep end # if solved end # while λ<λ_end - if !control.store_all # store last solutionobtained + if !control.store_all # store last solution obtained append!(tsol, λ0, solution) end - control.sample(solution, λ0) + + control.sample(solution, λ0) + if solved - if !(λ≈lambdas[i+1]) + if !(λ ≈ lambdas[i + 1]) # check end of interval has been reached in inner loop @warn "λ=$(λ), lambdas[i+1]=$(lambdas[i+1])" end else - break + break # emergency exit end - end # for i = 1:(length(lambdas)-1) + end # for i = 1:(length(lambdas)-1), end outer loop if doprint(control, 'e') println("[e]volution: $(round(t0+t1,sigdigits=3)) seconds") @@ -549,15 +557,13 @@ Keyword arguments: - `tstep`: time step Returns a [`DenseSolutionArray`](@ref) or [`SparseSolutionArray`](@ref) """ -function CommonSolve.solve( - sys::VoronoiFVM.AbstractSystem; - inival = 0, - params = zeros(0), - control = VoronoiFVM.SolverControl(), - time = 0.0, - tstep = Inf, - kwargs..., -) +function CommonSolve.solve(sys::VoronoiFVM.AbstractSystem; + inival = 0, + params = zeros(0), + control = VoronoiFVM.SolverControl(), + time = 0.0, + tstep = Inf, + kwargs...,) fix_deprecations!(control) if isa(inival, Number) @@ -577,27 +583,23 @@ function CommonSolve.solve( sys.linear_cache = nothing if haskey(kwargs, :times) && !isnothing(kwargs[:times]) - solve( - inival, - sys, - kwargs[:times]; - control, - transient = true, - params, - time = kwargs[:times][1], - called_from_API = true, - ) + solve(inival, + sys, + kwargs[:times]; + control, + transient = true, + params, + time = kwargs[:times][1], + called_from_API = true,) elseif haskey(kwargs, :embed) && !isnothing(kwargs[:embed]) - solve( - inival, - sys, - kwargs[:embed]; - called_from_API = true, - transient = false, - control, - params, - time, - ) + solve(inival, + sys, + kwargs[:embed]; + called_from_API = true, + transient = false, + control, + params, + time,) else solve(inival, sys; called_from_API = true, control, params, time, tstep) end diff --git a/src/vfvm_solvercontrol.jl b/src/vfvm_solvercontrol.jl index 3a34db181..9bb841a7a 100644 --- a/src/vfvm_solvercontrol.jl +++ b/src/vfvm_solvercontrol.jl @@ -1,11 +1,11 @@ ################################################ -const FactorizationStrategy= Union{Nothing, Function, Type, ExtendableSparse.AbstractFactorization, LinearSolve.AbstractFactorization,LinearSolve.SciMLLinearSolveAlgorithm} +const FactorizationStrategy = Union{Nothing, Function, Type, ExtendableSparse.AbstractFactorization, + LinearSolve.AbstractFactorization, LinearSolve.SciMLLinearSolveAlgorithm} struct Identity end -Identity(A)=Identity() -LinearAlgebra.ldiv!(u,I::Identity,v)=u.=v -LinearAlgebra.ldiv!(I::Identity,u)=nothing - +Identity(A) = Identity() +LinearAlgebra.ldiv!(u, I::Identity, v) = u .= v +LinearAlgebra.ldiv!(I::Identity, u) = nothing """ SolverControl @@ -106,7 +106,7 @@ Base.@kwdef mutable struct SolverControl Users should experiment with what works best for their problem. For easy access to this functionality, see see also [`VoronoiFVM.LinearSolverStrategy`](@ref). - + """ method_linear::Union{Nothing, LinearSolve.SciMLLinearSolveAlgorithm} = nothing @@ -136,7 +136,7 @@ Base.@kwdef mutable struct SolverControl For easy access to this functionality, see see also [`VoronoiFVM.LinearSolverStrategy`](@ref). """ - precon_linear::FactorizationStrategy= A -> Identity() + precon_linear::FactorizationStrategy = A -> Identity() """ Update preconditioner in each Newton step ? @@ -206,6 +206,12 @@ Base.@kwdef mutable struct SolverControl """ force_first_step::Bool = false + """ + Number of final steps to adjust at end of time interval in order + to prevent breakdown of step size. + """ + num_final_steps::Int = 5 + """ Handle exceptions during transient solver and parameter embedding. If `true`, exceptions in Newton solves are caught, the embedding resp. time step is lowered, @@ -335,5 +341,3 @@ end Legacy name of SolverControl """ const NewtonControl = SolverControl - -