From ce38f7a36c106e1b42a6751ae02c9c3602f584a8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Thu, 14 Dec 2023 22:13:39 +0100
Subject: [PATCH 1/7] Fix linear solver - UMFPACK as default for sparse - Fix
matirxtype=:banded
---
src/vfvm_linsolve.jl | 15 +++++++++------
src/vfvm_solver.jl | 12 +++---------
2 files changed, 12 insertions(+), 15 deletions(-)
diff --git a/src/vfvm_linsolve.jl b/src/vfvm_linsolve.jl
index 4898a9d79..85eeea9a1 100644
--- a/src/vfvm_linsolve.jl
+++ b/src/vfvm_linsolve.jl
@@ -228,25 +228,28 @@ function LinearAlgebra.ldiv!(u, cache::LinearSolve.LinearCache, b)
copyto!(u, sol.u)
end
+canonical_matrix(A)=A
+canonical_matrix(A::ExtendableSparseMatrix)=SparseMatrixCSC(A)
+
function _solve_linear!(u, system, nlhistory, control, method_linear, A, b)
if isnothing(system.linear_cache)
- Pl = control.precon_linear(SparseMatrixCSC(A))
+ Pl = control.precon_linear(canonical_matrix(A))
nlhistory.nlu += 1
- p = LinearProblem(SparseMatrixCSC(A), b)
+ p = LinearProblem(canonical_matrix(A), b)
system.linear_cache = init(
- p,
- method_linear;
+ p;
+ method_linear
abstol = control.abstol_linear,
reltol = control.reltol_linear,
verbose = doprint(control, 'l'),
Pl,
)
else
- system.linear_cache.A=SparseMatrixCSC(A)
+ system.linear_cache.A=canonical_matrix(A)
system.linear_cache.b=b
if control.keepcurrent_linear
nlhistory.nlu += 1
- system.linear_cache.Pl=control.precon_linear(SparseMatrixCSC(A))
+ system.linear_cache.Pl=control.precon_linear(canonical_matrix(A))
end
end
diff --git a/src/vfvm_solver.jl b/src/vfvm_solver.jl
index bb0177df5..262935074 100644
--- a/src/vfvm_solver.jl
+++ b/src/vfvm_solver.jl
@@ -29,18 +29,12 @@ function _solve_timestep!(solution::AbstractMatrix{Tv}, # old time step solution
update = system.update
_initialize!(solution, system; time, λ = embedparam, params)
- method_linear = control.method_linear
- if isnothing(method_linear)
+ method_linear matrixtype == :sparse ? control.method_linear : nothing;
+ if isnothing(method_linear) && matrixtype == :sparse
if Tv != Float64
method_linear = SparspakFactorization()
else
- if dim_grid(system.grid) == 1
- method_linear = KLUFactorization()
- elseif dim_grid(system.grid) == 2
- method_linear = SparspakFactorization()
- else
- method_linear = UMFPACKFactorization()
- end
+ method_linear = UMFPACKFactorization() # seems to do the best pivoting
end
end
From 7d8972017c563c9aaa66b365c921b9696ca8f1a7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Fri, 15 Dec 2023 23:28:13 +0100
Subject: [PATCH 2/7] Bugfix for assembly of outflow bc
---
src/vfvm_assembly.jl | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/src/vfvm_assembly.jl b/src/vfvm_assembly.jl
index a177905ad..1b3abbf4a 100644
--- a/src/vfvm_assembly.jl
+++ b/src/vfvm_assembly.jl
@@ -273,7 +273,7 @@ function eval_and_assemble(
system.matrix,
idofK,
jdofL,
- -jac_outflow[ispec, jspec+nspecies],
+ jac_outflow[ispec, jspec+nspecies],
edge.fac
)
end
@@ -283,7 +283,7 @@ function eval_and_assemble(
system.matrix,
idofL,
jdofK,
- +jac_outflow[ispec, jspec],
+ -jac_outflow[ispec, jspec],
edge.fac
)
_addnz(
From 3836a39568c48ee88ee7a03b488bd051294aa9ce Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Fri, 15 Dec 2023 23:29:08 +0100
Subject: [PATCH 3/7] fix linsolve bugs - enable banded
---
src/vfvm_linsolve.jl | 4 ++--
src/vfvm_solver.jl | 8 ++++++--
2 files changed, 8 insertions(+), 4 deletions(-)
diff --git a/src/vfvm_linsolve.jl b/src/vfvm_linsolve.jl
index 85eeea9a1..cf2bc13ae 100644
--- a/src/vfvm_linsolve.jl
+++ b/src/vfvm_linsolve.jl
@@ -237,8 +237,8 @@ function _solve_linear!(u, system, nlhistory, control, method_linear, A, b)
nlhistory.nlu += 1
p = LinearProblem(canonical_matrix(A), b)
system.linear_cache = init(
- p;
- method_linear
+ p,
+ method_linear;
abstol = control.abstol_linear,
reltol = control.reltol_linear,
verbose = doprint(control, 'l'),
diff --git a/src/vfvm_solver.jl b/src/vfvm_solver.jl
index 262935074..db1d69743 100644
--- a/src/vfvm_solver.jl
+++ b/src/vfvm_solver.jl
@@ -29,10 +29,14 @@ function _solve_timestep!(solution::AbstractMatrix{Tv}, # old time step solution
update = system.update
_initialize!(solution, system; time, λ = embedparam, params)
- method_linear matrixtype == :sparse ? control.method_linear : nothing;
- if isnothing(method_linear) && matrixtype == :sparse
+ method_linear = system.matrixtype == :sparse ? control.method_linear : nothing;
+ if isnothing(method_linear) && system.matrixtype == :sparse
if Tv != Float64
method_linear = SparspakFactorization()
+ elseif dim_space(system.grid)==1
+ method_linear = KLUFactorization()
+ elseif dim_space(system.grid)==2
+ method_linear = SparspakFactorization()
else
method_linear = UMFPACKFactorization() # seems to do the best pivoting
end
From 3cd1126824336cb839f2b4ed2e1d611e5ad382a7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Fri, 15 Dec 2023 23:30:27 +0100
Subject: [PATCH 4/7] Add plotting of newton updates
---
src/gridvisualize.jl | 41 +++++++++++++++++++++++++++++++++--------
1 file changed, 33 insertions(+), 8 deletions(-)
diff --git a/src/gridvisualize.jl b/src/gridvisualize.jl
index 74f2024da..f5a8c642c 100644
--- a/src/gridvisualize.jl
+++ b/src/gridvisualize.jl
@@ -1,5 +1,6 @@
# Extend GridVisualize methods
import GridVisualize
+import Colors
"""
$(TYPEDSIGNATURES)
@@ -79,14 +80,18 @@ end
"""
plothistory(sys, tsol;
- plots=[:timesteps,:newtonsteps,:updates],
- size=(700,400),
+ plots=[:timestepsizes,
+ :timestepupdates,
+ :newtonsteps,
+ :newtonupdates],
+ size=(700,600),
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...)
+function plothistory(tsol::TransientSolution; plots=[:timestepsizes,:timestepupdates,:newtonsteps,:newtonupdates],
+ size=(700,150*length(plots)), Plotter=GridVisualize.default_plotter(), kwargs...)
hist=history(tsol)
t=hist.times
if length(t)==0
@@ -96,17 +101,37 @@ function plothistory(sys, tsol; plots=[:timesteps,:newtonsteps,:updates], size=(
vis=GridVisualize.GridVisualizer(;layout=(nplots,1), size, Plotter, kwargs...)
for iplot in eachindex(plots)
- if plots[iplot]==:timesteps
+ if plots[iplot]==:timestepsizes
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]==:timestepupdates
+ GridVisualize.scalarplot!(vis[iplot,1],t, hist.updates,xlabel = "t/s", ylabel = "du",
+ title="Time step updates", 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...)
+ title="Newton steps", color=:red, limits=(0,maximum(newtons)+1), kwargs...)
+ elseif plots[iplot]==:newtonupdates
+ nhist=length(hist)
+ dc=1/nhist
+ r=0.0
+ b=1.0
+ g=0.0
+ for h in hist
+ GridVisualize.scalarplot!(vis[iplot,1],
+ 1:length(h),h,
+ xlabel="step",
+ ylabel="du",
+ clear=false,
+ yscale=:log,
+ color=Colors.RGB(r,g,b),
+ linewidth=0.5, title="Newton updates")
+ r+=dc
+ b-=dc
+ end
end
end
GridVisualize.reveal(vis)
end
+
+plothistory(sys,tsol;kwargs...)=plothistory(tsol,kwargs...)
From 9c20dbada37e6ee45567eed3d86ee1b50f872430 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sat, 16 Dec 2023 00:03:04 +0100
Subject: [PATCH 5/7] bump version
---
Project.toml | 4 +++-
docs/src/changes.md | 5 +++++
examples/Example120_ThreeRegions1D.jl | 21 ++++++++++++---------
3 files changed, 20 insertions(+), 10 deletions(-)
diff --git a/Project.toml b/Project.toml
index 164996dd6..09f4ca582 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,10 +1,11 @@
name = "VoronoiFVM"
uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
authors = ["Jürgen Fuhrmann <juergen.fuhrmann@wias-berlin.de>", "Dilara Abdel", "Jan Weidner", "Alexander Seiler", "Patricio Farrell", "Matthias Liero"]
-version = "1.15.1"
+version = "1.16.0"
[deps]
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -28,6 +29,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
[compat]
BandedMatrices = "0.17,1"
+Colors = "0.12"
CommonSolve = "0.2"
DiffResults = "1"
DocStringExtensions = "0.8,0.9"
diff --git a/docs/src/changes.md b/docs/src/changes.md
index 03937ba1e..d001be4b6 100644
--- a/docs/src/changes.md
+++ b/docs/src/changes.md
@@ -1,4 +1,9 @@
# Changes
+## v1.16.0 Dec 15, 2023
+- Bugfix for assembly of outflow bc
+- Bugfix for matrixtype=:banded
+- Updated [`plothistory`](@ref):
+
## 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,
diff --git a/examples/Example120_ThreeRegions1D.jl b/examples/Example120_ThreeRegions1D.jl
index a9d35b604..ba431f743 100644
--- a/examples/Example120_ThreeRegions1D.jl
+++ b/examples/Example120_ThreeRegions1D.jl
@@ -7,6 +7,7 @@ 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,
@@ -124,22 +125,24 @@ function main(; n = 30, Plotter = nothing, plot_grid = false, verbose = false,
clear = false, show = true)
end
- tsol = solve(sys; inival = 0, times = (0, tend), post = plot_timestep, verbose = verbose, Δu_opt = 1.0e-5)
+ tsol = solve(sys; inival = 0, times = (0, tend), post = plot_timestep, verbose = verbose, Δu_opt = 1.0e-5,
+ method_linear=KLUFactorization())
return testval
end
using Test
+
function runtests()
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 &&
- main(; unknown_storage = :dense, rely_on_corrections = true, assembly = :edgewise) ≈ testval &&
- main(; unknown_storage = :sparse, rely_on_corrections = false, assembly = :cellwise) ≈ testval &&
- main(; unknown_storage = :dense, rely_on_corrections = false, assembly = :cellwise) ≈ testval &&
- main(; unknown_storage = :sparse, rely_on_corrections = true, assembly = :cellwise) ≈ testval &&
- main(; unknown_storage = :dense, rely_on_corrections = true, assembly = :cellwise) ≈ testval
+ @test main(; unknown_storage = :sparse, rely_on_corrections = false, assembly = :edgewise) ≈ testval
+ @test main(; unknown_storage = :dense, rely_on_corrections = false, assembly = :edgewise) ≈ testval
+ @test main(; unknown_storage = :sparse, rely_on_corrections = true, assembly = :edgewise) ≈ testval
+ @test main(; unknown_storage = :dense, rely_on_corrections = true, assembly = :edgewise) ≈ testval
+ @test main(; unknown_storage = :sparse, rely_on_corrections = false, assembly = :cellwise) ≈ testval
+ @test main(; unknown_storage = :dense, rely_on_corrections = false, assembly = :cellwise) ≈ testval
+ @test main(; unknown_storage = :sparse, rely_on_corrections = true, assembly = :cellwise) ≈ testval
+ @test main(; unknown_storage = :dense, rely_on_corrections = true, assembly = :cellwise) ≈ testval
end
end
From 987ee07ff5b2462f9a8e1ccbf479382c1354c2aa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sat, 16 Dec 2023 01:10:56 +0100
Subject: [PATCH 6/7] Add compat 3 on RecursiveArrayTools
---
Project.toml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Project.toml b/Project.toml
index 09f4ca582..4c39d987b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -42,7 +42,7 @@ LinearAlgebra = "1.6"
LinearSolve = "2.2"
Printf = "1.6"
Random = "1.6"
-RecursiveArrayTools = "^2.14.2"
+RecursiveArrayTools = "^2.14.2,3"
RecursiveFactorization = "0.2.18"
SparseArrays = "1.6"
SparseDiffTools = "^1.19, 2"
From 3f8e7863185df16a15ad77a6dc2d5669afa6e2be Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sat, 16 Dec 2023 12:49:25 +0100
Subject: [PATCH 7/7] update .JuliaFormatter.toml
---
src/.JuliaFormatter.toml | 3 +++
1 file changed, 3 insertions(+)
diff --git a/src/.JuliaFormatter.toml b/src/.JuliaFormatter.toml
index 905ee5359..5d28f1b62 100644
--- a/src/.JuliaFormatter.toml
+++ b/src/.JuliaFormatter.toml
@@ -1,4 +1,7 @@
style = "sciml"
+pipe_to_function_call=false
always_for_in = false
separate_kwargs_with_semicolon = true
margin = 132
+yas_style_nesting=true
+ignore = "examples/PlutoGridFactory.jl"
\ No newline at end of file