gridap · amartinhuertas · Apr 5, 2023 · Mar 22, 2023 · Mar 23, 2023 · Mar 23, 2023
diff --git a/NEWS.md b/NEWS.md
@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- ODE operators cache linear system at initial time or the time stored by the operator. Before, the linear system was cached at time `t = 0.0`, which cannot be done if the operator is not well-defined at `t = 0.0`. Since PR [#891](https://github.com/gridap/Gridap.jl/pull/891).
 - Fixed the method `get_normal_vector` for `AdaptedTriangulation`. The method `get_facet_normal`
   was using default, it's now using the spetialized implementation for the underlying triangulation type.
   Since PR [#884](https://github.com/gridap/Gridap.jl/pull/884).

diff --git a/src/ODEs/ODETools/AffineNewmark.jl b/src/ODEs/ODETools/AffineNewmark.jl
@@ -20,7 +20,7 @@ function solve_step!(
     (v,a, ode_cache) = newmark_cache
 
     # Allocate matrices and vectors
-    A, b = _allocate_matrix_and_vector(op,x0,ode_cache)
+    A, b = _allocate_matrix_and_vector(op,t0,x0,ode_cache)
 
     # Create affine operator cache
     affOp_cache = (A,b,nothing)
@@ -101,14 +101,14 @@ function jacobian!(A::AbstractMatrix,op::NewmarkAffineOperator,x::AbstractVector
   jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
 end
 
-function _allocate_matrix(odeop::ODEOperator,x::Tuple{Vararg{AbstractVector}},ode_cache)
-  A = allocate_jacobian(odeop,x[1],ode_cache)
+function _allocate_matrix(odeop::ODEOperator,t0::Real,x::Tuple{Vararg{AbstractVector}},ode_cache)
+  A = allocate_jacobian(odeop,t0,x[1],ode_cache)
   return A
 end
 
-function _allocate_matrix_and_vector(odeop::ODEOperator,x::Tuple{Vararg{AbstractVector}},ode_cache)
-  b = allocate_residual(odeop,x[1],ode_cache)
-  A = allocate_jacobian(odeop,x[1],ode_cache)
+function _allocate_matrix_and_vector(odeop::ODEOperator,t0::Real,x::Tuple{Vararg{AbstractVector}},ode_cache)
+  b = allocate_residual(odeop,t0,x[1],ode_cache)
+  A = allocate_jacobian(odeop,t0,x[1],ode_cache)
   return A, b
 end
 

diff --git a/src/ODEs/ODETools/AffineThetaMethod.jl b/src/ODEs/ODETools/AffineThetaMethod.jl
@@ -14,7 +14,7 @@ function solve_step!(uf::AbstractVector,
     vθ = similar(u0)
     vθ .= 0.0
     l_cache = nothing
-    A, b = _allocate_matrix_and_vector(op,u0,ode_cache)
+    A, b = _allocate_matrix_and_vector(op,t0,u0,ode_cache)
   else
     ode_cache, vθ, A, b, l_cache = cache
   end
@@ -54,10 +54,10 @@ function solve_step!(uf::AbstractVector,
     ode_cache = allocate_cache(op)
     vθ = similar(u0)
     vθ .= 0.0
-    A, b = _allocate_matrix_and_vector(op,u0,ode_cache)
+    A, b = _allocate_matrix_and_vector(op,t0,u0,ode_cache)
     A = _matrix!(A,op,tθ,dtθ,u0,ode_cache,vθ)
     b = _vector!(b,op,tθ,dtθ,vθ,ode_cache,vθ)
-    M = _allocate_matrix(op,u0,ode_cache)
+    M = _allocate_matrix(op,t0,u0,ode_cache)
     M = _mass_matrix!(M,op,tθ,dtθ,u0,ode_cache,vθ)
     _u0 =  similar(u0,(axes(M)[2],)) # Needed for the distributed case
     copy!(_u0,u0)
@@ -94,7 +94,7 @@ given time step, i.e., M(t)(u_n+θ-u_n)/dt + K(t)u_n+θ + b(t)
 function ThetaMethodAffineOperator(odeop::AffineODEOperator,tθ::Float64,dtθ::Float64,
                                    u0::AbstractVector,ode_cache,vθ::AbstractVector)
   # vθ .= 0.0
-  A, b = _allocate_matrix_and_vector(odeop,u0,ode_cache)
+  A, b = _allocate_matrix_and_vector(odeop,t0,u0,ode_cache)
   _matrix_and_vector!(A,b,odeop,tθ,dtθ,u0,ode_cache,vθ)
   afop = AffineOperator(A,b)
 end
@@ -121,14 +121,14 @@ function _vector!(b,odeop,tθ,dtθ,u0,ode_cache,vθ)
   b .*= -1.0
 end
 
-function _allocate_matrix(odeop,u0,ode_cache)
-  A = allocate_jacobian(odeop,u0,ode_cache)
+function _allocate_matrix(odeop,t0,u0,ode_cache)
+  A = allocate_jacobian(odeop,t0,u0,ode_cache)
   return A
 end
 
-function _allocate_matrix_and_vector(odeop,u0,ode_cache)
-  b = allocate_residual(odeop,u0,ode_cache)
-  A = allocate_jacobian(odeop,u0,ode_cache)
+function _allocate_matrix_and_vector(odeop,t0,u0,ode_cache)
+  b = allocate_residual(odeop,t0,u0,ode_cache)
+  A = allocate_jacobian(odeop,t0,u0,ode_cache)
   return A, b
 end
 
@@ -137,9 +137,9 @@ Affine operator that represents the θ-method affine operator at a
 given time step, i.e., M(t)(u_n+θ-u_n)/dt + K(t)u_n+θ + b(t)
 """
 function ThetaMethodConstantOperator(odeop::ConstantODEOperator,tθ::Float64,dtθ::Float64,
-                                   u0::AbstractVector,ode_cache,vθ::AbstractVector)
-  b = allocate_residual(odeop,u0,ode_cache)
-  A = allocate_jacobian(odeop,u0,ode_cache)
+                                     u0::AbstractVector,ode_cache,vθ::AbstractVector)
+  b = allocate_residual(odeop,tθ,u0,ode_cache)
+  A = allocate_jacobian(odeop,tθ,u0,ode_cache)
   residual!(b,odeop,tθ,(u0,vθ),ode_cache)
   @. b = -1.0 * b
   z = zero(eltype(A))

diff --git a/src/ODEs/ODETools/ConstantMatrixNewmark.jl b/src/ODEs/ODETools/ConstantMatrixNewmark.jl
@@ -24,7 +24,7 @@ function solve_step!(
     newmark_affOp = NewmarkConstantMatrixOperator(op,t1,dt,γ,β,x0,newmark_cache)
 
     # Allocate matrices and vectors
-    A, b = _allocate_matrix_and_vector(op,x0,ode_cache)
+    A, b = _allocate_matrix_and_vector(op,t0,x0,ode_cache)
     jacobian!(A,newmark_affOp,u1)
 
     # Create affine operator cache
@@ -96,13 +96,13 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantMatrixOperator,x::Abstra
   jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
 end
 
-function _allocate_matrix(odeop::NewmarkConstantMatrixOperator,x::Tuple{Vararg{AbstractVector}},ode_cache)
-  A = allocate_jacobian(odeop,x[1],ode_cache)
+function _allocate_matrix(odeop::NewmarkConstantMatrixOperator,t0::Real,x::Tuple{Vararg{AbstractVector}},ode_cache)
+  A = allocate_jacobian(odeop,t0,x[1],ode_cache)
   return A
 end
 
-function _allocate_matrix_and_vector(odeop::NewmarkConstantMatrixOperator,x::Tuple{Vararg{AbstractVector}},ode_cache)
-  b = allocate_residual(odeop,x[1],ode_cache)
-  A = allocate_jacobian(odeop,x[1],ode_cache)
+function _allocate_matrix_and_vector(odeop::NewmarkConstantMatrixOperator,t0::Real,x::Tuple{Vararg{AbstractVector}},ode_cache)
+  b = allocate_residual(odeop,t0,x[1],ode_cache)
+  A = allocate_jacobian(odeop,t0,x[1],ode_cache)
   return A, b
 end
diff --git a/src/ODEs/ODETools/ConstantNewmark.jl b/src/ODEs/ODETools/ConstantNewmark.jl
@@ -22,9 +22,9 @@ function solve_step!(
     (v,a, ode_cache) = newmark_cache
 
     # Allocate matrices and vectors
-    A, b = _allocate_matrix_and_vector(op,x0,ode_cache)
-    M = _allocate_matrix(op,x0,ode_cache)
-    C = _allocate_matrix(op,x0,ode_cache)
+    A, b = _allocate_matrix_and_vector(op,t0,x0,ode_cache)
+    M = _allocate_matrix(op,t0,x0,ode_cache)
+    C = _allocate_matrix(op,t0,x0,ode_cache)
     b1 = similar(b)
     b1 .= 0.0
 

diff --git a/src/ODEs/ODETools/ForwardEuler.jl b/src/ODEs/ODETools/ForwardEuler.jl
@@ -64,11 +64,11 @@ function jacobian!(A::AbstractMatrix,op::ForwardEulerNonlinearOperator,x::Abstra
 end
 
 function allocate_residual(op::ForwardEulerNonlinearOperator,x::AbstractVector)
-  allocate_residual(op.odeop,x,op.ode_cache)
+  allocate_residual(op.odeop,op.tf,x,op.ode_cache)
 end
 
 function allocate_jacobian(op::ForwardEulerNonlinearOperator,x::AbstractVector)
-  allocate_jacobian(op.odeop,x,op.ode_cache)
+  allocate_jacobian(op.odeop,op.tf,x,op.ode_cache)
 end
 
 function zero_initial_guess(op::ForwardEulerNonlinearOperator)

diff --git a/src/ODEs/ODETools/GeneralizedAlpha.jl b/src/ODEs/ODETools/GeneralizedAlpha.jl
@@ -135,12 +135,12 @@ end
 
 function allocate_residual(op::GeneralizedAlphaNonlinearOperator,x::AbstractVector)
   vαm, cache = op.ode_cache
-  allocate_residual(op.odeop,x,cache)
+  allocate_residual(op.odeop,op.tαf,x,cache)
 end
 
 function allocate_jacobian(op::GeneralizedAlphaNonlinearOperator,x::AbstractVector)
   vαm, cache = op.ode_cache
-  allocate_jacobian(op.odeop,x,cache)
+  allocate_jacobian(op.odeop,op.tαf,x,cache)
 end
 
 function zero_initial_guess(op::GeneralizedAlphaNonlinearOperator)
@@ -214,13 +214,13 @@ end
 
 function allocate_residual(op::GeneralizedAlphaDttNonlinearOperator,x::AbstractVector)
   vαf, aαm, cache = op.ode_cache
-  allocate_residual(op.odeop,x,cache)
+  allocate_residual(op.odeop,op.tαf,x,cache)
 end
 
 
 function allocate_jacobian(op::GeneralizedAlphaDttNonlinearOperator,x::AbstractVector)
   vαf, aαm, cache = op.ode_cache
-  allocate_jacobian(op.odeop,x,cache)
+  allocate_jacobian(op.odeop,op.tαf,x,cache)
 end
 
 

diff --git a/src/ODEs/ODETools/Newmark.jl b/src/ODEs/ODETools/Newmark.jl
@@ -81,12 +81,12 @@ end
 
 function allocate_residual(op::NewmarkNonlinearOperator,x::AbstractVector)
   v1, a1, cache = op.ode_cache
-  allocate_residual(op.odeop,x,cache)
+  allocate_residual(op.odeop,op.t1,x,cache)
 end
 
 function allocate_jacobian(op::NewmarkNonlinearOperator,x::AbstractVector)
   v1, a1, cache = op.ode_cache
-  allocate_jacobian(op.odeop,x,cache)
+  allocate_jacobian(op.odeop,op.t1,x,cache)
 end
 
 function zero_initial_guess(op::NewmarkNonlinearOperator)

diff --git a/src/ODEs/ODETools/ODEOperators.jl b/src/ODEs/ODETools/ODEOperators.jl
@@ -65,6 +65,7 @@ end
 """
 function allocate_residual(
   op::ODEOperator,
+  t0::Real,
   u::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
   ode_cache)
   @abstractmethod
@@ -106,7 +107,7 @@ end
 
 """
 """
-function allocate_jacobian(op::ODEOperator,u::AbstractVector,ode_cache)
+function allocate_jacobian(op::ODEOperator,t0::Real,u::AbstractVector,ode_cache)
   @abstractmethod
 end
 
@@ -124,9 +125,9 @@ Tests the interface of `ODEOperator` specializations
 function test_ode_operator(op::ODEOperator,t::Real,u::AbstractVector,u_t::AbstractVector)
   cache = allocate_cache(op)
   cache = update_cache!(cache,op,0.0)
-  r = allocate_residual(op,u,cache)
+  r = allocate_residual(op,0.0,u,cache)
   residual!(r,op,t,(u,u_t),cache)
-  J = allocate_jacobian(op,u,cache)
+  J = allocate_jacobian(op,0.0,u,cache)
   jacobian!(J,op,t,(u,u_t),1,1.0,cache)
   jacobian!(J,op,t,(u,u_t),2,1.0,cache)
   jacobians!(J,op,t,(u,u_t),(1.0,1.0),cache)

diff --git a/src/ODEs/ODETools/RungeKutta.jl b/src/ODEs/ODETools/RungeKutta.jl
@@ -180,11 +180,11 @@ function jacobian!(A::AbstractMatrix,op::RungeKuttaNonlinearOperator,x::Abstract
 end
 
 function allocate_residual(op::RungeKuttaNonlinearOperator,x::AbstractVector)
-  allocate_residual(op.odeop,x,op.ode_cache)
+  allocate_residual(op.odeop,op.ti,x,op.ode_cache)
 end
 
 function allocate_jacobian(op::RungeKuttaNonlinearOperator,x::AbstractVector)
-  allocate_jacobian(op.odeop,x,op.ode_cache)
+  allocate_jacobian(op.odeop,op.ti,x,op.ode_cache)
 end
 
 function zero_initial_guess(op::RungeKuttaNonlinearOperator)

diff --git a/src/ODEs/ODETools/ThetaMethod.jl b/src/ODEs/ODETools/ThetaMethod.jl
@@ -84,11 +84,11 @@ function jacobian!(A::AbstractMatrix,op::ThetaMethodNonlinearOperator,x::Abstrac
 end
 
 function allocate_residual(op::ThetaMethodNonlinearOperator,x::AbstractVector)
-  allocate_residual(op.odeop,x,op.ode_cache)
+  allocate_residual(op.odeop,op.tθ,x,op.ode_cache)
 end
 
 function allocate_jacobian(op::ThetaMethodNonlinearOperator,x::AbstractVector)
-  allocate_jacobian(op.odeop,x,op.ode_cache)
+  allocate_jacobian(op.odeop,op.tθ,x,op.ode_cache)
 end
 
 function zero_initial_guess(op::ThetaMethodNonlinearOperator)

diff --git a/src/ODEs/TransientFETools/ODEOperatorInterfaces.jl b/src/ODEs/TransientFETools/ODEOperatorInterfaces.jl
@@ -47,16 +47,16 @@ function update_cache!(ode_cache,op::ODEOpFromFEOp,t::Real)
   (Us,Uts,fecache)
 end
 
-function allocate_residual(op::ODEOpFromFEOp,uhF::AbstractVector,ode_cache)
+function allocate_residual(op::ODEOpFromFEOp,t0::Real,uhF::AbstractVector,ode_cache)
   Us,Uts,fecache = ode_cache
   uh = EvaluationFunction(Us[1],uhF)
-  allocate_residual(op.feop,uh,fecache)
+  allocate_residual(op.feop,t0,uh,fecache)
 end
 
-function allocate_jacobian(op::ODEOpFromFEOp,uhF::AbstractVector,ode_cache)
+function allocate_jacobian(op::ODEOpFromFEOp,t0::Real,uhF::AbstractVector,ode_cache)
   Us,Uts,fecache = ode_cache
   uh = EvaluationFunction(Us[1],uhF)
-  allocate_jacobian(op.feop,uh,fecache)
+  allocate_jacobian(op.feop,t0,uh,fecache)
 end
 
 """

diff --git a/src/ODEs/TransientFETools/TransientFEOperators.jl b/src/ODEs/TransientFETools/TransientFEOperators.jl
@@ -17,11 +17,11 @@ function get_trial(op::TransientFEOperator)
   @abstractmethod # time dependent
 end
 
-function allocate_residual(op::TransientFEOperator,uh,cache)
+function allocate_residual(op::TransientFEOperator,t0,uh,cache)
   @abstractmethod
 end
 
-function allocate_jacobian(op::TransientFEOperator,uh,cache)
+function allocate_jacobian(op::TransientFEOperator,t0,uh,cache)
   @notimplemented
 end
 
@@ -231,6 +231,7 @@ get_order(op::TransientFEOperatorFromWeakForm) = op.order
 
 function allocate_residual(
   op::TransientFEOperatorFromWeakForm,
+  t0::Real,
   uh::T,
   cache) where T
   V = get_test(op)
@@ -240,7 +241,7 @@ function allocate_residual(
     dxh = (dxh...,uh)
   end
   xh = TransientCellField(uh,dxh)
-  vecdata = collect_cell_vector(V,op.res(0.0,xh,v))
+  vecdata = collect_cell_vector(V,op.res(t0,xh,v))
   allocate_vector(op.assem_t,vecdata)
 end
 
@@ -259,9 +260,10 @@ end
 
 function allocate_jacobian(
   op::TransientFEOperatorFromWeakForm,
+  t0::Real,
   uh::CellField,
   cache)
-  _matdata_jacobians = fill_initial_jacobians(op,uh)
+  _matdata_jacobians = fill_initial_jacobians(op,t0,uh)
   matdata = _vcat_matdata(_matdata_jacobians)
   allocate_matrix(op.assem_t,matdata)
 end
@@ -292,15 +294,15 @@ function jacobians!(
   A
 end
 
-function fill_initial_jacobians(op::TransientFEOperatorFromWeakForm,uh)
+function fill_initial_jacobians(op::TransientFEOperatorFromWeakForm,t0::Real,uh)
   dxh = ()
   for i in 1:get_order(op)
     dxh = (dxh...,uh)
   end
   xh = TransientCellField(uh,dxh)
   _matdata = ()
   for i in 1:get_order(op)+1
-    _matdata = (_matdata...,_matdata_jacobian(op,0.0,xh,i,0.0))
+    _matdata = (_matdata...,_matdata_jacobian(op,t0,xh,i,0.0))
   end
   return _matdata
 end
@@ -360,12 +362,12 @@ function test_transient_fe_operator(op::TransientFEOperator,uh)
   U = get_trial(op)
   U0 = U(0.0)
   @test isa(U0,FESpace)
-  r = allocate_residual(op,uh,cache)
+  r = allocate_residual(op,0.0,uh,cache)
   @test isa(r,AbstractVector)
   xh = TransientCellField(uh,(uh,))
   residual!(r,op,0.0,xh,cache)
   @test isa(r,AbstractVector)
-  J = allocate_jacobian(op,uh,cache)
+  J = allocate_jacobian(op,0.0,uh,cache)
   @test isa(J,AbstractMatrix)
   jacobian!(J,op,0.0,xh,1,1.0,cache)
   @test isa(J,AbstractMatrix)

diff --git a/test/ODEsTests/ODEsTests/ODEOperatorMocks.jl b/test/ODEsTests/ODEsTests/ODEOperatorMocks.jl
@@ -43,7 +43,7 @@ function residual!(r::AbstractVector,op::ODEOperatorMock,t::Real,x::NTuple{3,Abs
   r
 end
 
-function allocate_residual(op::ODEOperatorMock,u::AbstractVector,cache)
+function allocate_residual(op::ODEOperatorMock,t0::Real,u::AbstractVector,cache)
   zeros(2)
 end
 
@@ -104,7 +104,7 @@ function jacobians!(
   J
 end
 
-function allocate_jacobian(op::ODEOperatorMock,u::AbstractVector,cache)
+function allocate_jacobian(op::ODEOperatorMock,t0::Real,u::AbstractVector,cache)
   spzeros(2,2)
 end
 

diff --git a/test/ODEsTests/ODEsTests/ODEOperatorsTests.jl b/test/ODEsTests/ODEsTests/ODEOperatorsTests.jl
@@ -18,13 +18,13 @@ u_t = ones(2)*2.0
 cache = allocate_cache(op)
 update_cache!(cache,op,0.0)
 
-r = allocate_residual(op,u,cache)
+t = 0.0
+r = allocate_residual(op,t,u,cache)
 @test r == zeros(2)
 
-J = allocate_jacobian(op,u,cache)
+J = allocate_jacobian(op,t,u,cache)
 @test J == zeros(2,2)
 
-t = 0.0
 residual!(r,op,t,(u,u_t),cache)
 _r = zeros(2)
 _r[1] = u_t[1] - op.a * u[1]