From c6c9ece549b15d786fff86f201feeebc765e77e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 20 Jun 2023 15:11:25 +0200 Subject: [PATCH] Linearsolve2 (#74) * mutable linsolve * linearsolve 2.0.1 * update project.toml --- Project.toml | 4 +- examples/Example205_NonlinearPoisson2D.jl | 2 +- pluto-examples/api-update.jl | 2 +- src/VoronoiFVM.jl | 1 - src/vfvm_linsolve.jl | 47 +++++++---------------- src/vfvm_solvercontrol.jl | 7 ++++ src/vfvm_testfunctions.jl | 2 +- 7 files changed, 25 insertions(+), 40 deletions(-) diff --git a/Project.toml b/Project.toml index 035f54ad2..1c438946b 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.7.1" +version = "1.8.0" [deps] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" @@ -39,7 +39,7 @@ ExtendableSparse = "1.2" ForwardDiff = "0.10.35" GridVisualize = "0.5.2,0.6.1,1" JLD2 = "0.4.29" -LinearSolve = "1.37" +LinearSolve = "2.2" PrecompileTools = "1" RecursiveArrayTools = "^2.14.2" RecursiveFactorization = "0.2.18" diff --git a/examples/Example205_NonlinearPoisson2D.jl b/examples/Example205_NonlinearPoisson2D.jl index 25089dbdd..74f6e397b 100644 --- a/examples/Example205_NonlinearPoisson2D.jl +++ b/examples/Example205_NonlinearPoisson2D.jl @@ -13,7 +13,7 @@ using ILUZero function main(; n = 10, Plotter = nothing, verbose = false, unknown_storage = :sparse, method_linear = nothing, assembly=:edgewise, - precon_linear = A -> LinearSolve.Identity()) + precon_linear = A -> VoronoiFVM.Identity()) h = 1.0 / convert(Float64, n) X = collect(0.0:h:1.0) Y = collect(0.0:h:1.0) diff --git a/pluto-examples/api-update.jl b/pluto-examples/api-update.jl index 6c736d2c9..b3d2ab40a 100644 --- a/pluto-examples/api-update.jl +++ b/pluto-examples/api-update.jl @@ -225,7 +225,7 @@ md""" sys0; inival = 0.1, method_linear = KrylovJL_BICGSTAB(), - precon_linear = A -> factorize(A, SparspakFactorization()), + precon_linear = SparspakFactorization(), verbose="nlad" ) diff --git a/src/VoronoiFVM.jl b/src/VoronoiFVM.jl index 1fb202e19..6beb47ae5 100644 --- a/src/VoronoiFVM.jl +++ b/src/VoronoiFVM.jl @@ -13,7 +13,6 @@ using BandedMatrices # using MultidiagonalMatrices using LinearSolve -using LinearSolve: Identity using Statistics using ForwardDiff diff --git a/src/vfvm_linsolve.jl b/src/vfvm_linsolve.jl index 10a035e0a..4898a9d79 100644 --- a/src/vfvm_linsolve.jl +++ b/src/vfvm_linsolve.jl @@ -204,45 +204,27 @@ function Random.rand( end -mutable struct FactorizationPreconditioner{C} - cache::C -end - """ - factorize(A,method::LinearSolve.AbstractFactorization) - -Calculate an LU factorization of A using one of the methods available in LinearSolve.jl. + Make preconditioner constructors from methods """ -function LinearAlgebra.factorize(A, method::LinearSolve.AbstractFactorization) + +function (method::LinearSolve.AbstractFactorization)(A) pr = LinearProblem(A, zeros(eltype(A), size(A, 1))) - p = FactorizationPreconditioner(init(pr, method)) - p + init(pr, method) end -function LinearAlgebra.factorize(A, method::LinearSolve.SciMLLinearSolveAlgorithm) +function (method::LinearSolve.SciMLLinearSolveAlgorithm)(A) pr = LinearProblem(SparseMatrixCSC(A), zeros(eltype(A), size(A, 1))) - p = FactorizationPreconditioner(init(pr, method)) - p -end - - -function (f::LinearSolve.AbstractFactorization)(A) - factorize(A, f) -end - -function (f::LinearSolve.SciMLLinearSolveAlgorithm)(A) - factorize(SparseMatrixCSC(A), f) + init(pr, method) end function (f::ExtendableSparse.AbstractFactorization)(A) factorize!(f, A) end - -function LinearAlgebra.ldiv!(u, p::FactorizationPreconditioner, b) - p.cache = LinearSolve.set_b(p.cache, b) - sol = solve(p.cache) - p.cache = sol.cache +function LinearAlgebra.ldiv!(u, cache::LinearSolve.LinearCache, b) + cache.b=b + sol = solve!(cache) copyto!(u, sol.u) end @@ -260,19 +242,16 @@ function _solve_linear!(u, system, nlhistory, control, method_linear, A, b) Pl, ) else - system.linear_cache = LinearSolve.set_A(system.linear_cache, SparseMatrixCSC(A)) - system.linear_cache = LinearSolve.set_b(system.linear_cache, b) + system.linear_cache.A=SparseMatrixCSC(A) + system.linear_cache.b=b if control.keepcurrent_linear - Pl = control.precon_linear(SparseMatrixCSC(A)) nlhistory.nlu += 1 - system.linear_cache = - LinearSolve.set_prec(system.linear_cache, Pl, LinearSolve.Identity()) + system.linear_cache.Pl=control.precon_linear(SparseMatrixCSC(A)) end end try - sol = LinearSolve.solve(system.linear_cache) - system.linear_cache = sol.cache + sol = LinearSolve.solve!(system.linear_cache) u .= sol.u nliniter = sol.iters nlhistory.nlin = sol.iters diff --git a/src/vfvm_solvercontrol.jl b/src/vfvm_solvercontrol.jl index 42e6cb33b..013ab5c1f 100644 --- a/src/vfvm_solvercontrol.jl +++ b/src/vfvm_solvercontrol.jl @@ -1,5 +1,12 @@ ################################################ 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 + + """ SolverControl SolverControl(;kwargs...) diff --git a/src/vfvm_testfunctions.jl b/src/vfvm_testfunctions.jl index 67bcdb3bd..d2569618b 100644 --- a/src/vfvm_testfunctions.jl +++ b/src/vfvm_testfunctions.jl @@ -79,7 +79,7 @@ function testfunction(factory::TestFunctionFactory, bc0, bc1) method_linear = UMFPACKFactorization() end - p = LinearProblem(factory.tfsystem.matrix, vec(f)) + p = LinearProblem(SparseMatrixCSC(factory.tfsystem.matrix), vec(f)) sol = solve(p, method_linear) sol.u end