From 93bc5697c09ee9d25cd7d1a3574a755f6c52ccd9 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Mon, 28 Feb 2022 15:20:04 +0100 Subject: [PATCH 1/3] fill! to fillstored! --- src/GridapODEs/ODETools/AffineNewmark.jl | 2 +- src/GridapODEs/ODETools/AffineThetaMethod.jl | 6 +++--- src/GridapODEs/ODETools/ConstantMatrixNewmark.jl | 2 +- src/GridapODEs/ODETools/ConstantNewmark.jl | 6 +++--- src/GridapODEs/ODETools/ForwardEuler.jl | 4 ++-- src/GridapODEs/ODETools/Newmark.jl | 4 ++-- src/GridapODEs/ODETools/ODETools.jl | 2 +- src/GridapODEs/ODETools/RungeKutta.jl | 4 ++-- src/GridapODEs/ODETools/ThetaMethod.jl | 4 ++-- 9 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/GridapODEs/ODETools/AffineNewmark.jl b/src/GridapODEs/ODETools/AffineNewmark.jl index 905af110d..7e229e71c 100644 --- a/src/GridapODEs/ODETools/AffineNewmark.jl +++ b/src/GridapODEs/ODETools/AffineNewmark.jl @@ -97,7 +97,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkAffineOperator,x::AbstractVector a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end diff --git a/src/GridapODEs/ODETools/AffineThetaMethod.jl b/src/GridapODEs/ODETools/AffineThetaMethod.jl index afa110ce6..15ae7e7f7 100644 --- a/src/GridapODEs/ODETools/AffineThetaMethod.jl +++ b/src/GridapODEs/ODETools/AffineThetaMethod.jl @@ -106,13 +106,13 @@ end function _matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ) z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache) end function _mass_matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ) z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobian!(A,odeop,tθ,(vθ,vθ),2,(1/dtθ),ode_cache) end @@ -143,7 +143,7 @@ function ThetaMethodConstantOperator(odeop::ConstantODEOperator,tθ::Float64,dt residual!(b,odeop,tθ,(u0,vθ),ode_cache) b = -1*b z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache) return A, b end diff --git a/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl b/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl index d40a55bc9..06699bcf0 100644 --- a/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl +++ b/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl @@ -92,7 +92,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantMatrixOperator,x::Abstra a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end diff --git a/src/GridapODEs/ODETools/ConstantNewmark.jl b/src/GridapODEs/ODETools/ConstantNewmark.jl index 1b32f67a8..685a117df 100644 --- a/src/GridapODEs/ODETools/ConstantNewmark.jl +++ b/src/GridapODEs/ODETools/ConstantNewmark.jl @@ -125,7 +125,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantOperator,x::AbstractVect a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end @@ -136,7 +136,7 @@ function _mass_matrix!(A::AbstractMatrix,op::NewmarkConstantOperator,x::Abstract a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobian!(A,op.odeop,op.t1,(u1,v1,a1),3,1.0,cache) end @@ -147,6 +147,6 @@ function _damping_matrix!(A::AbstractMatrix,op::NewmarkConstantOperator,x::Abstr a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobian!(A,op.odeop,op.t1,(u1,v1,a1),2,1.0,cache) end diff --git a/src/GridapODEs/ODETools/ForwardEuler.jl b/src/GridapODEs/ODETools/ForwardEuler.jl index 7b68b2f72..5fefb0a86 100644 --- a/src/GridapODEs/ODETools/ForwardEuler.jl +++ b/src/GridapODEs/ODETools/ForwardEuler.jl @@ -59,7 +59,7 @@ function jacobian!(A::AbstractMatrix,op::ForwardEulerNonlinearOperator,x::Abstra vf = op.vf vf = (x-op.u0)/op.dt z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,op.odeop,op.tf,(op.u0,vf),(0,1/op.dt),op.ode_cache) end @@ -73,6 +73,6 @@ end function zero_initial_guess(op::ForwardEulerNonlinearOperator) x0 = similar(op.u0) - fill!(x0,zero(eltype(x0))) + fillstored!(x0,zero(eltype(x0))) x0 end diff --git a/src/GridapODEs/ODETools/Newmark.jl b/src/GridapODEs/ODETools/Newmark.jl index 8588c0ce4..aa04ad808 100644 --- a/src/GridapODEs/ODETools/Newmark.jl +++ b/src/GridapODEs/ODETools/Newmark.jl @@ -75,7 +75,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkNonlinearOperator,x::AbstractVec a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end @@ -91,6 +91,6 @@ end function zero_initial_guess(op::NewmarkNonlinearOperator) x0 = similar(op.u0) - fill!(x0,zero(eltype(x0))) + fillstored!(x0,zero(eltype(x0))) x0 end diff --git a/src/GridapODEs/ODETools/ODETools.jl b/src/GridapODEs/ODETools/ODETools.jl index 9e4924796..81f88a480 100644 --- a/src/GridapODEs/ODETools/ODETools.jl +++ b/src/GridapODEs/ODETools/ODETools.jl @@ -10,7 +10,7 @@ using Test using DocStringExtensions using ForwardDiff -using LinearAlgebra: fill! +using LinearAlgebra: fillstored! const ϵ = 100*eps() export ∂t diff --git a/src/GridapODEs/ODETools/RungeKutta.jl b/src/GridapODEs/ODETools/RungeKutta.jl index 53a2ad900..151dcf007 100644 --- a/src/GridapODEs/ODETools/RungeKutta.jl +++ b/src/GridapODEs/ODETools/RungeKutta.jl @@ -175,7 +175,7 @@ function jacobian!(A::AbstractMatrix,op::RungeKuttaNonlinearOperator,x::Abstract vi = op.vi vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt) z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,op.odeop,op.ti,(ui,vi),(1.0,1.0/(op.a[op.i,op.i]*op.dt)),op.ode_cache) end @@ -189,7 +189,7 @@ end function zero_initial_guess(op::RungeKuttaNonlinearOperator) x0 = similar(op.u0) - fill!(x0,zero(eltype(x0))) + fillstored!(x0,zero(eltype(x0))) x0 end diff --git a/src/GridapODEs/ODETools/ThetaMethod.jl b/src/GridapODEs/ODETools/ThetaMethod.jl index 2bc772d53..42f265ef2 100644 --- a/src/GridapODEs/ODETools/ThetaMethod.jl +++ b/src/GridapODEs/ODETools/ThetaMethod.jl @@ -79,7 +79,7 @@ function jacobian!(A::AbstractMatrix,op::ThetaMethodNonlinearOperator,x::Abstrac vθ = op.vθ vθ = (x-op.u0)/op.dtθ z = zero(eltype(A)) - fill!(A,z) + fillstored!(A,z) jacobians!(A,op.odeop,op.tθ,(uF,vθ),(1.0,1/op.dtθ),op.ode_cache) end @@ -93,6 +93,6 @@ end function zero_initial_guess(op::ThetaMethodNonlinearOperator) x0 = similar(op.u0) - fill!(x0,zero(eltype(x0))) + fillstored!(x0,zero(eltype(x0))) x0 end From 9efb107ee7bc0601650d9274fa52dad2c1aec378 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Mon, 28 Feb 2022 20:46:59 +0100 Subject: [PATCH 2/3] fill! for dense matrices and fillstored! for sparse matrices --- src/GridapODEs/ODETools/AffineNewmark.jl | 6 +++++- src/GridapODEs/ODETools/AffineThetaMethod.jl | 18 +++++++++++++++--- .../ODETools/ConstantMatrixNewmark.jl | 6 +++++- src/GridapODEs/ODETools/ConstantNewmark.jl | 18 +++++++++++++++--- src/GridapODEs/ODETools/ForwardEuler.jl | 6 +++++- src/GridapODEs/ODETools/Newmark.jl | 6 +++++- src/GridapODEs/ODETools/ODETools.jl | 1 + src/GridapODEs/ODETools/RungeKutta.jl | 12 ++++++++++-- src/GridapODEs/ODETools/ThetaMethod.jl | 12 ++++++++++-- 9 files changed, 71 insertions(+), 14 deletions(-) diff --git a/src/GridapODEs/ODETools/AffineNewmark.jl b/src/GridapODEs/ODETools/AffineNewmark.jl index 7e229e71c..935db7e2a 100644 --- a/src/GridapODEs/ODETools/AffineNewmark.jl +++ b/src/GridapODEs/ODETools/AffineNewmark.jl @@ -97,7 +97,11 @@ function jacobian!(A::AbstractMatrix,op::NewmarkAffineOperator,x::AbstractVector a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end diff --git a/src/GridapODEs/ODETools/AffineThetaMethod.jl b/src/GridapODEs/ODETools/AffineThetaMethod.jl index 15ae7e7f7..9b69d7c9d 100644 --- a/src/GridapODEs/ODETools/AffineThetaMethod.jl +++ b/src/GridapODEs/ODETools/AffineThetaMethod.jl @@ -106,13 +106,21 @@ end function _matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ) z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache) end function _mass_matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ) z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobian!(A,odeop,tθ,(vθ,vθ),2,(1/dtθ),ode_cache) end @@ -143,7 +151,11 @@ function ThetaMethodConstantOperator(odeop::ConstantODEOperator,tθ::Float64,dt residual!(b,odeop,tθ,(u0,vθ),ode_cache) b = -1*b z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache) return A, b end diff --git a/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl b/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl index 06699bcf0..c8755d9e9 100644 --- a/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl +++ b/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl @@ -92,7 +92,11 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantMatrixOperator,x::Abstra a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end diff --git a/src/GridapODEs/ODETools/ConstantNewmark.jl b/src/GridapODEs/ODETools/ConstantNewmark.jl index 685a117df..14a5e767c 100644 --- a/src/GridapODEs/ODETools/ConstantNewmark.jl +++ b/src/GridapODEs/ODETools/ConstantNewmark.jl @@ -125,7 +125,11 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantOperator,x::AbstractVect a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end @@ -136,7 +140,11 @@ function _mass_matrix!(A::AbstractMatrix,op::NewmarkConstantOperator,x::Abstract a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobian!(A,op.odeop,op.t1,(u1,v1,a1),3,1.0,cache) end @@ -147,6 +155,10 @@ function _damping_matrix!(A::AbstractMatrix,op::NewmarkConstantOperator,x::Abstr a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobian!(A,op.odeop,op.t1,(u1,v1,a1),2,1.0,cache) end diff --git a/src/GridapODEs/ODETools/ForwardEuler.jl b/src/GridapODEs/ODETools/ForwardEuler.jl index 5fefb0a86..2e9693466 100644 --- a/src/GridapODEs/ODETools/ForwardEuler.jl +++ b/src/GridapODEs/ODETools/ForwardEuler.jl @@ -73,6 +73,10 @@ end function zero_initial_guess(op::ForwardEulerNonlinearOperator) x0 = similar(op.u0) - fillstored!(x0,zero(eltype(x0))) + if issparse(x0) + fillstored!(x0,zero(eltype(x0))) + else + fill!(x0,zero(eltype(x0))) + end x0 end diff --git a/src/GridapODEs/ODETools/Newmark.jl b/src/GridapODEs/ODETools/Newmark.jl index aa04ad808..3cc273b08 100644 --- a/src/GridapODEs/ODETools/Newmark.jl +++ b/src/GridapODEs/ODETools/Newmark.jl @@ -91,6 +91,10 @@ end function zero_initial_guess(op::NewmarkNonlinearOperator) x0 = similar(op.u0) - fillstored!(x0,zero(eltype(x0))) + if issparse(x0) + fillstored!(x0,zero(eltype(x0))) + else + fill!(x0,zero(eltype(x0))) + end x0 end diff --git a/src/GridapODEs/ODETools/ODETools.jl b/src/GridapODEs/ODETools/ODETools.jl index 81f88a480..bb06d8b1a 100644 --- a/src/GridapODEs/ODETools/ODETools.jl +++ b/src/GridapODEs/ODETools/ODETools.jl @@ -11,6 +11,7 @@ using DocStringExtensions using ForwardDiff using LinearAlgebra: fillstored! +using SparseArrays: issparse const ϵ = 100*eps() export ∂t diff --git a/src/GridapODEs/ODETools/RungeKutta.jl b/src/GridapODEs/ODETools/RungeKutta.jl index 151dcf007..7c435a39f 100644 --- a/src/GridapODEs/ODETools/RungeKutta.jl +++ b/src/GridapODEs/ODETools/RungeKutta.jl @@ -175,7 +175,11 @@ function jacobian!(A::AbstractMatrix,op::RungeKuttaNonlinearOperator,x::Abstract vi = op.vi vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt) z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobians!(A,op.odeop,op.ti,(ui,vi),(1.0,1.0/(op.a[op.i,op.i]*op.dt)),op.ode_cache) end @@ -189,7 +193,11 @@ end function zero_initial_guess(op::RungeKuttaNonlinearOperator) x0 = similar(op.u0) - fillstored!(x0,zero(eltype(x0))) + if issparse(x0) + fillstored!(x0,zero(eltype(x0))) + else + fill!(x0,zero(eltype(x0))) + end x0 end diff --git a/src/GridapODEs/ODETools/ThetaMethod.jl b/src/GridapODEs/ODETools/ThetaMethod.jl index 42f265ef2..6767996fe 100644 --- a/src/GridapODEs/ODETools/ThetaMethod.jl +++ b/src/GridapODEs/ODETools/ThetaMethod.jl @@ -79,7 +79,11 @@ function jacobian!(A::AbstractMatrix,op::ThetaMethodNonlinearOperator,x::Abstrac vθ = op.vθ vθ = (x-op.u0)/op.dtθ z = zero(eltype(A)) - fillstored!(A,z) + if issparse(A) + fillstored!(A,z) + else + fill!(A,z) + end jacobians!(A,op.odeop,op.tθ,(uF,vθ),(1.0,1/op.dtθ),op.ode_cache) end @@ -93,6 +97,10 @@ end function zero_initial_guess(op::ThetaMethodNonlinearOperator) x0 = similar(op.u0) - fillstored!(x0,zero(eltype(x0))) + if issparse(x0) + fillstored!(x0,zero(eltype(x0))) + else + fill!(x0,zero(eltype(x0))) + end x0 end From a02502870aecb5ce2840bffe73a516ed906e5e06 Mon Sep 17 00:00:00 2001 From: Oriol Colomes Date: Wed, 2 Mar 2022 10:33:18 +0100 Subject: [PATCH 3/3] sparse matrix in ODEOperatorMocks --- src/GridapODEs/ODETools/AffineNewmark.jl | 6 +----- src/GridapODEs/ODETools/AffineThetaMethod.jl | 18 +++--------------- .../ODETools/ConstantMatrixNewmark.jl | 6 +----- src/GridapODEs/ODETools/ConstantNewmark.jl | 18 +++--------------- src/GridapODEs/ODETools/ForwardEuler.jl | 6 +----- src/GridapODEs/ODETools/Newmark.jl | 6 +----- src/GridapODEs/ODETools/RungeKutta.jl | 12 ++---------- src/GridapODEs/ODETools/ThetaMethod.jl | 12 ++---------- .../ODEsTests/ODEOperatorMocks.jl | 3 ++- .../ODEsTests/ODESolverMocks.jl | 4 ++-- 10 files changed, 18 insertions(+), 73 deletions(-) diff --git a/src/GridapODEs/ODETools/AffineNewmark.jl b/src/GridapODEs/ODETools/AffineNewmark.jl index 935db7e2a..7e229e71c 100644 --- a/src/GridapODEs/ODETools/AffineNewmark.jl +++ b/src/GridapODEs/ODETools/AffineNewmark.jl @@ -97,11 +97,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkAffineOperator,x::AbstractVector a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end diff --git a/src/GridapODEs/ODETools/AffineThetaMethod.jl b/src/GridapODEs/ODETools/AffineThetaMethod.jl index 9b69d7c9d..15ae7e7f7 100644 --- a/src/GridapODEs/ODETools/AffineThetaMethod.jl +++ b/src/GridapODEs/ODETools/AffineThetaMethod.jl @@ -106,21 +106,13 @@ end function _matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ) z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache) end function _mass_matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ) z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobian!(A,odeop,tθ,(vθ,vθ),2,(1/dtθ),ode_cache) end @@ -151,11 +143,7 @@ function ThetaMethodConstantOperator(odeop::ConstantODEOperator,tθ::Float64,dt residual!(b,odeop,tθ,(u0,vθ),ode_cache) b = -1*b z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache) return A, b end diff --git a/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl b/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl index c8755d9e9..06699bcf0 100644 --- a/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl +++ b/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl @@ -92,11 +92,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantMatrixOperator,x::Abstra a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end diff --git a/src/GridapODEs/ODETools/ConstantNewmark.jl b/src/GridapODEs/ODETools/ConstantNewmark.jl index 14a5e767c..685a117df 100644 --- a/src/GridapODEs/ODETools/ConstantNewmark.jl +++ b/src/GridapODEs/ODETools/ConstantNewmark.jl @@ -125,11 +125,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantOperator,x::AbstractVect a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache) end @@ -140,11 +136,7 @@ function _mass_matrix!(A::AbstractMatrix,op::NewmarkConstantOperator,x::Abstract a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobian!(A,op.odeop,op.t1,(u1,v1,a1),3,1.0,cache) end @@ -155,10 +147,6 @@ function _damping_matrix!(A::AbstractMatrix,op::NewmarkConstantOperator,x::Abstr a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0 v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0 z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobian!(A,op.odeop,op.t1,(u1,v1,a1),2,1.0,cache) end diff --git a/src/GridapODEs/ODETools/ForwardEuler.jl b/src/GridapODEs/ODETools/ForwardEuler.jl index 2e9693466..b177b55b6 100644 --- a/src/GridapODEs/ODETools/ForwardEuler.jl +++ b/src/GridapODEs/ODETools/ForwardEuler.jl @@ -73,10 +73,6 @@ end function zero_initial_guess(op::ForwardEulerNonlinearOperator) x0 = similar(op.u0) - if issparse(x0) - fillstored!(x0,zero(eltype(x0))) - else - fill!(x0,zero(eltype(x0))) - end + fill!(x0,zero(eltype(x0))) x0 end diff --git a/src/GridapODEs/ODETools/Newmark.jl b/src/GridapODEs/ODETools/Newmark.jl index 3cc273b08..eda2fe93f 100644 --- a/src/GridapODEs/ODETools/Newmark.jl +++ b/src/GridapODEs/ODETools/Newmark.jl @@ -91,10 +91,6 @@ end function zero_initial_guess(op::NewmarkNonlinearOperator) x0 = similar(op.u0) - if issparse(x0) - fillstored!(x0,zero(eltype(x0))) - else - fill!(x0,zero(eltype(x0))) - end + fill!(x0,zero(eltype(x0))) x0 end diff --git a/src/GridapODEs/ODETools/RungeKutta.jl b/src/GridapODEs/ODETools/RungeKutta.jl index 7c435a39f..ee2768380 100644 --- a/src/GridapODEs/ODETools/RungeKutta.jl +++ b/src/GridapODEs/ODETools/RungeKutta.jl @@ -175,11 +175,7 @@ function jacobian!(A::AbstractMatrix,op::RungeKuttaNonlinearOperator,x::Abstract vi = op.vi vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt) z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobians!(A,op.odeop,op.ti,(ui,vi),(1.0,1.0/(op.a[op.i,op.i]*op.dt)),op.ode_cache) end @@ -193,11 +189,7 @@ end function zero_initial_guess(op::RungeKuttaNonlinearOperator) x0 = similar(op.u0) - if issparse(x0) - fillstored!(x0,zero(eltype(x0))) - else - fill!(x0,zero(eltype(x0))) - end + fill!(x0,zero(eltype(x0))) x0 end diff --git a/src/GridapODEs/ODETools/ThetaMethod.jl b/src/GridapODEs/ODETools/ThetaMethod.jl index 6767996fe..dae23eb32 100644 --- a/src/GridapODEs/ODETools/ThetaMethod.jl +++ b/src/GridapODEs/ODETools/ThetaMethod.jl @@ -79,11 +79,7 @@ function jacobian!(A::AbstractMatrix,op::ThetaMethodNonlinearOperator,x::Abstrac vθ = op.vθ vθ = (x-op.u0)/op.dtθ z = zero(eltype(A)) - if issparse(A) - fillstored!(A,z) - else - fill!(A,z) - end + fillstored!(A,z) jacobians!(A,op.odeop,op.tθ,(uF,vθ),(1.0,1/op.dtθ),op.ode_cache) end @@ -97,10 +93,6 @@ end function zero_initial_guess(op::ThetaMethodNonlinearOperator) x0 = similar(op.u0) - if issparse(x0) - fillstored!(x0,zero(eltype(x0))) - else - fill!(x0,zero(eltype(x0))) - end + fill!(x0,zero(eltype(x0))) x0 end diff --git a/test/GridapODEsTests/ODEsTests/ODEOperatorMocks.jl b/test/GridapODEsTests/ODEsTests/ODEOperatorMocks.jl index ba7453901..5635aa7c6 100644 --- a/test/GridapODEsTests/ODEsTests/ODEOperatorMocks.jl +++ b/test/GridapODEsTests/ODEsTests/ODEOperatorMocks.jl @@ -16,6 +16,7 @@ import Gridap.GridapODEs.ODETools: jacobian! import Gridap.GridapODEs.ODETools: jacobians! import Gridap.GridapODEs.ODETools: allocate_jacobian import Gridap.GridapODEs.ODETools: residual! +using SparseArrays: spzeros struct ODEOperatorMock{T<:Real,C} <: ODEOperator{C} a::T @@ -104,7 +105,7 @@ function jacobians!( end function allocate_jacobian(op::ODEOperatorMock,u::AbstractVector,cache) - zeros(2,2) + spzeros(2,2) end allocate_cache(op::ODEOperatorMock) = nothing diff --git a/test/GridapODEsTests/ODEsTests/ODESolverMocks.jl b/test/GridapODEsTests/ODEsTests/ODESolverMocks.jl index b65c5dc12..7ae341376 100644 --- a/test/GridapODEsTests/ODEsTests/ODESolverMocks.jl +++ b/test/GridapODEsTests/ODEsTests/ODESolverMocks.jl @@ -60,7 +60,7 @@ end function solve!(x::AbstractVector,nls::NLSolverMock,nlop::NonlinearOperator,cache::Nothing) r = residual(nlop,x) J = jacobian(nlop,x) - dx = inv(J)*(-r) + dx = inv(Matrix(J))*(-r) x.= x.+dx cache = (r,J,dx) end @@ -69,7 +69,7 @@ function solve!(x::AbstractVector,nls::NLSolverMock,nlop::NonlinearOperator,cach r, J, dx = cache residual!(r, nlop, x) jacobian!(J, nlop, x) - dx = inv(J)*(-r) + dx = inv(Matrix(J))*(-r) x.= x.+dx end