Sundials ERROR: UndefVarError: LS not defined
#431
Comments
|
Lines 70 to 77 in b01e0c3 GMRES. @ChrisRackauckas might be able to say more on this
|
|
On a sidenote, for most of our benchmarks we beat the nonlinear solver in Sundials by quite a margin. |
FWIW, I tried |
Do these benchmarks include KINSOL's GMRES or KLU? |
|
They do not yet, but it's not clear why that would change anything. As a quasi-Newton algorithm it would see less optimizations in the GMRES case, so that's a worst case scenario for it. With KLU it's still a non-globalizing method. See: https://docs.sciml.ai/SciMLBenchmarksOutput/dev/NonlinearProblem/nonlinear_solver_23_tests/ See, the point is that it doesn't have a globalizing part (line search, trust region, etc.) and makes heavy convergence assumptions in a quasi-Newton form. So it's issue isn't speed (though it's not the fastest) as much as convergence. It only successfully solves 15 of the 23 test cases, while some of the native NonlinearSolve.jl algorithms are both faster and converge on all 23. From all of the testing we've been able to do, we have not seen KINSOL as competitive at all, so it effectively only exists as a benchmarking tool and not something we'd ever recommend in practice. |
|
But yes, if someone wants this, you just need to do Sundials.jl/src/common_interface/solve.jl Lines 263 to 295 in b01e0c3 Lines 70 to 77 in b01e0c3 |
|
Below is a very ugly draft implementation of gmres for kinsol with
precondition. kinsol_prec is the entry point.
using Sundials
import Sundials: N_Vector, SUNMatrix, @checkflag, KIN_SUCCESS,
UserFunctionAndData, PREC_RIGHT, PREC_NONE, KINCreate, KINSetUserData,
KINInit, Handle, NVector
using LinearAlgebra
using SparseArrays
using IncompleteLU
abstract type AbstractKinsolData{J} end
abstract type AbstractKinsolFunData{J,P} end
abstract type AbstractKinsolPreconditioner{M} end
abstract type AbstractPreconditionerSettings end
struct ILU end
struct AlgebraicMultiGrid end
struct KinsolSettings
linear_solver::Symbol
y0::Vector
jac_lower::Int64
jac_upper::Int64
prec_side::Int64
krylov_dim::Int64
print_level::Int64
end
# struct KinsolFunData{F,D,J,P} <: AbstractKinsolFunData{J,P}
# func::F
# data::D
# end
struct KinsolData{J,P}
jac::J
preconditioner::P
end
# KinsolFunData(func::F, kinsol_data::KinsolData{J,P}) where {F,J,P} =
KinsolFunData{F,D,J,P}(func, kinsol_data)
# KinsolFunData(func::F, jac::J, precond::P) where {F,J<:Array,P} =
KinsolFunData(func, KinsolData(jac, precond))
struct KinsolJacobian end
ILUKinsolPreconditioner = AbstractKinsolPreconditioner{ILU}
AlgebraicMultiGridKinsolPreconditioner =
AbstractKinsolPreconditioner{AlgebraicMultiGrid}
mutable struct KinsolPreconditioner{A,J,S,M} <:
AbstractKinsolPreconditioner{M}
P::A
compute_jacobian::J
settings::S
method::M
end
@kwdef mutable struct AlgebraicMultiGridSettings <:
AbstractPreconditionerSettings
τ::Float64 = 1e-6
end
@kwdef mutable struct ILUSettings <: AbstractPreconditionerSettings
τ::Float64 = 1e-6
end
function KinsolPreconditioner(jac, method::Symbol; kw...)
if method == :ilu
settings = ILUSettings(; kw...)
P = compute_P_ilu(jac.jac, settings)
method_ = ILU()
KinsolPreconditioner(P, jac, settings, method_)
elseif method == :multigrid
settings = AlgebraicMultiGridSettings(; kw...)
P = compute_P_AlgebraicMultiGrid(jac.jac, settings)
method_ = AlgebraicMultiGrid()
KinsolPreconditioner(P, jac, settings, method_)
else
error()
end
end
function compute_P!(precond::T) where {T<:ILUKinsolPreconditioner}
println("precond.settings.τ = $(precond.settings.τ)")
precond.P = ilu(precond.compute_jacobian.jac, τ=precond.settings.τ)
end
compute_P_ilu(jac, settings) = ilu(jac, τ=settings.τ)
compute_P_AlgebraicMultiGrid(jac, settings) = aspreconditioner(
smoothed_aggregation(jac, presmoother=AlgebraicMultigrid.Jacobi(rand(size(jac,
1))), postsmoother=AlgebraicMultigrid.Jacobi(rand(size(jac, 1)))))
psetup(u, uscale, fval, fscal, data::KinsolData{J,P}) where {J,P<:
KinsolPreconditioner} = setup_P!(u, data.preconditioner)
psolve(u, uscale, fval, fscal, v, data::KinsolData{J,P}) where {J,P<:
KinsolPreconditioner} = solve!(data.preconditioner, u, uscale, fval, fscal,
v)
function _psetup(u::N_Vector,
uscale::N_Vector,
fval::N_Vector,
fscale::N_Vector,
userfun::UserFunctionAndData)
psetup(convert(Vector, u),
convert(Vector, uscale),
convert(Vector, fval),
convert(Vector, fscale),
userfun.data)
return KIN_SUCCESS
end
function _psolve(u::N_Vector,
uscale::N_Vector,
f::N_Vector,
fscale::N_Vector,
v::N_Vector,
userfun::UserFunctionAndData
)
psolve(convert(Vector, u),
convert(Vector, uscale),
convert(Vector, f),
convert(Vector, fscale),
convert(Vector, v),
userfun.data)
return KIN_SUCCESS
end
function compute_jac!(pc::KinsolPreconditioner{A,J,S,M}, u::Vector{Float64})
where {A,J,S,M}
println("computing jacobian")
pc.compute_jacobian(u)
end
import Sundials: kinsolfun
# function kinsolfun(y::N_Vector, fy::N_Vector, userfun::KinsolFunData)
# userfun[].func(convert(Vector, fy), convert(Vector, y), userfun[].data)
# return KIN_SUCCESS
# end
function _kinsoljac(
u::N_Vector,
fu::N_Vector,
J::SUNMatrix,
userfun::UserFunctionAndData,
tmp1::N_Vector,
tmp2::N_Vector)
userfun.data.preconditioner.compute_jacobian(convert(Vector, u))
copyto!(J, userfun.data.preconditioner.compute_jacobian.jac)
return KIN_SUCCESS
end
function kinsol_prec(f,
y0::Vector{Float64};
preconditioner::Union{Nothing,AbstractKinsolPreconditioner,Symbol}=nothing,
jac::Union{KinsolJacobian,Nothing}=nothing,
set_jacobian=false,
linear_solver=:Band, ftol=nothing,
jac_upper=0,
jac_lower=0,
su_scale=nothing,
sf_scale=nothing,
printlevel::Int64=1,
krylov_dim=0,
constraints=zeros(length(y0)),
strategy=:none, max_iter=25
)
uvec = vec(deepcopy(y0))
userdata = KinsolData(jac, preconditioner)
# f, Function to be optimized of the form f(y::Vector{Float64},
fy::Vector{Float64})
# where `y` is the input vector, and `fy` is the result of the function
# y0, Vector of initial values
# return: the solution vector
# --- initialize ---
mem_ptr = KINCreate()
(mem_ptr == C_NULL) && error("Failed to allocate KINSOL solver object")
kmem = Handle(mem_ptr)
# use the user_data field to pass a function
# see: JuliaLang/julia#2554
userfun = UserFunctionAndData{typeof(f)}(f, userdata)
function getcfun(::T) where {T}
@cfunction(kinsolfun, Cint, (N_Vector, N_Vector, Ref{T}))
end
flag = @checkflag KINInit(kmem, getcfun(userfun), NVector(y0)) true
flag = @checkflag KINSetUserData(kmem, userfun) true # this is where the
F(U) is passed to kinsol through userfun pointer
# --- preconditioner ---
if preconditioner isa KinsolPreconditioner
prec_side = PREC_RIGHT
else
prec_side = PREC_NONE
end
@Assert prec_side ∈ [PREC_NONE, PREC_RIGHT]
# --- linear solver ---
# J is template for Jacobian
if linear_solver == :Dense
J = Sundials.SUNDenseMatrix(neq, length(uvec))
LS = Sundials.SUNLinSol_Dense(uvec, L)
elseif linear_solver == :LapackDense
J = Sundials.SUNDenseMatrix(length(uvec), length(uvec))
LS = Sundials.SUNLinSol_LapackDense(uvec, J)
elseif linear_solver == :Band
J = Sundials.SUNBandMatrix(length(uvec), jac_upper, jac_lower)
LS = Sundials.SUNLinSol_Band(uvec, J)
elseif linear_solver == :LapackBand
J = Sundials.SUNBandMatrix(length(uvec), jac_upper, jac_lower)
LS = Sundials.SUNLinSol_LapackBand(uvec, J)
elseif linear_solver == :GMRES
LS = Sundials.SUNLinSol_SPGMR(uvec, prec_side, krylov_dim)
J = Sundials.SUNBandMatrix(length(uvec), jac_upper, jac_lower)
elseif linear_solver == :FGMRES
LS = Sundials.SUNLinSol_SPFGMR(uvec, prec_side, krylov_dim)
J = Sundials.SUNBandMatrix(length(uvec), jac_upper, jac_lower)
elseif linear_solver == :KLU
nnz = length(SparseArrays.nonzeros(preconditioner.compute_jacobian.jac))
J = Sundials.SUNSparseMatrix(length(uvec), length(uvec), nnz, Sundials.
CSC_MAT)
LS = Sundials.SUNLinSol_KLU(uvec, J)
else
error("linear solver available for Kinsol are: [:Dense,:LapackDense,:Band,
:GMRES, :FGMRES, :KLU]")
end
flag = @checkflag Sundials.KINDlsSetLinearSolver(kmem, LS, J === nothing ?
C_NULL : J) true
J === nothing ? jacobian = false : jacobian = true
# --- jacobian ---
if preconditioner isa KinsolPreconditioner && set_jacobian
println("Kinsol: setting up Jacobian ...")
function getjacfun(::T) where {T}
@cfunction(_kinsoljac,
Cint,
(
N_Vector,
N_Vector,
SUNMatrix,
Ref{T},
N_Vector,
N_Vector))
end
jac = getjacfun(userfun)
flag = @checkflag Sundials.KINSetJacFn(kmem, jac) true
end
function getpsolvefun(::T) where {T}
@cfunction(_psolve,
Cint,
(
N_Vector,
N_Vector,
N_Vector,
N_Vector,
N_Vector,
Ref{T}))
end
function getpsetupfun(::T) where {T}
@cfunction(_psetup,
Cint,
(
N_Vector,
N_Vector,
N_Vector,
N_Vector,
Ref{T}))
end
psolvefun = getpsolvefun(userfun)
psetupfun = getpsetupfun(userfun)
# --- preconditionner setting ---
if preconditioner isa AbstractKinsolPreconditioner && linear_solver != :KLU
println("KINSOL: setting preconditioner...")
Sundials.KINSetPreconditioner(kmem, psetupfun, psolvefun)
end
## Solve problem
if isnothing(su_scale)
su_scale = ones(length(y0))
end
if isnothing(sf_scale)
sf_scale = ones(length(y0))
end
@Assert length(sf_scale) == length(y0)
@Assert length(su_scale) == length(y0)
if strategy == :none
strategy = Sundials.KIN_NONE
elseif strategy == :linesearch
strategy = Sundials.KIN_LINESEARCH
elseif strategy == :picard
strategy = Sundials.KIN_PICARD
else
error("strategy $strategy not available for kinsol. Available
strategies:$([:none,
:lineasearch, :picard])")
end
# jacobian
Sundials.KINSetPrintLevel(kmem, printlevel)
if ftol isa Number
Sundials.KINSetFuncNormTol(kmem, ftol)
end
Sundials.KINSetNumMaxIters(kmem, max_iter)
Sundials.KINSetConstraints(kmem, constraints)
etachoice = Sundials.KIN_ETACHOICE2
return kmem
end
function run_kinsol(kmem, y, strategy::Union{Symbol,Int32}, su_scale::
Union{Nothing,Vector{Float64}}, sf_scale::Union{Nothing,Vector{Float64}})
if isnothing(su_scale)
su_scale_ = ones(length(y))
else
su_scale_ = su_scale
end
if isnothing(sf_scale)
sf_scale = ones(length(y))
end
y_ = deepcopy(y)
#Sundials.KINSetEtaForm(kmem, etachoice)
#println("su_scale:", su_scale)
#println("sf_scale:", sf_scale)
flag = @checkflag Sundials.KINSol(kmem, y_, kinsol_strategy(strategy),
su_scale_, sf_scale) true
# Sundials.SUNLinSolFree_Dense(LS)
# Sundials.SUNMatDestroy_Dense(J)
return (y_, flag)
end
kinsol_strategy(strategy::Int32) = strategy
function kinsol_strategy(strategy::Symbol)
if strategy == :none
strategy = Sundials.KIN_NONE
elseif strategy == :linesearch
strategy = Sundials.KIN_LINESEARCH
elseif strategy == :picard
strategy = Sundials.KIN_PICARD
else
error("strategy $strategy not available for kinsol. Available
strategies:$([:none,
:lineasearch, :picard])")
end
return strategy
end
struct KinSOLPreconditioner{MP,J,F,S,M} <: AbstractKinsolPreconditioner{M}
P::MP
jac::J
evaluate_jac!::F
settings::S
method::M
end
function kinsolfun(y::N_Vector, fy::N_Vector, userfun::UserFunctionAndData)
userfun.func(convert(Vector, fy), convert(Vector, y), userfun.data)
return Sundials.KIN_SUCCESS
end
function setup_P!(u, precond)
#println("setup_P")
# Build preconditioner
compute_jac!(precond, u)
compute_P!(precond)
end
function solve!(precond::T, u, uscale, fval, fscal, v) where {T<:
ILUKinsolPreconditioner}
#println("solve Px=c")
#p = fscal .* v
ldiv!(precond.P, v)
#v .= p ./ fscal
end
function solve!(precond::T, u, uscale, fval, fscal, v) where {T<:
AlgebraicMultiGridKinsolPreconditioner}
#println("solve Px=c")
#p = fscal .* v
ldiv!(precond.P, v)
#v .= p ./ fscal
end
using AlgebraicMultigrid
function compute_P!(precond::T) where {T<:
AlgebraicMultiGridKinsolPreconditioner}
precond.P = aspreconditioner(smoothed_aggregation(precond.compute_jacobian.jac,
presmoother=AlgebraicMultigrid.Jacobi(rand(size(precond.compute_jacobian.jac,
1))), postsmoother=AlgebraicMultigrid.Jacobi(rand(size(precond.
compute_jacobian.jac, 1)))))
end
# function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev,
solverdata)
# if newW === nothing || newW
# Pl = aspreconditioner(ruge_stuben(convert(AbstractMatrix, W)))
# else
# Pl = Plprev
# end
# Pl, nothing
# end
############# ENTRY POINT
(out_ks, flag) = kinsol_prec((f, u, data) -> f!(f, u), U;
linear_solver=:GMRES,
printlevel=0,
krylov_dim=floor(Int64, sqrt(length(U))),
strategy=:none,
max_iter=50)
…On Wed, Nov 1, 2023 at 10:33 PM Benoit Pasquier ***@***.***> wrote:
Just tested solve(prob, KINSOL(;linear_solver=:GMRES)) on a simple
problem and it threw:
ERROR: UndefVarError: `LS` not defined
Stacktrace:
[1] ___kinsol(f::Function, y0::Vector{Float64}; userdata::Nothing, linear_solver::Symbol, jac_upper::Int64, jac_lower::Int64)
@ Sundials ~/.julia/packages/Sundials/vGumE/src/simple.jl:15
—
Reply to this email directly, view it on GitHub
<#431>, or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AEESMZDWFJHL4VTU6UZRW2TYCMV7TAVCNFSM6AAAAAA62KHFJWVHI2DSMVQWIX3LMV43ASLTON2WKOZRHE3TGNJTHEYTKNI>
.
You are receiving this because you are subscribed to this thread.Message
ID: ***@***.***>
|
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Just tested
solve(prob, KINSOL(;linear_solver=:GMRES))on a simple problem and it threw:The text was updated successfully, but these errors were encountered: