Merge branch 'master' of https://github.com/gridap/GridapODEs.jl into…

… Fixing_autodiff

NEWS.md

@@ -4,6 +4,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## Unreleased
+
+### Changed
+
+- Hiding the creation of `TransientFESolver` from user code. Since PR [#59](https://github.com/gridap/GridapODEs.jl/pull/59).
+
 ## [0.7.0] - 2021-09-21
 
 ### Added

Project.toml

@@ -13,7 +13,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 ForwardDiff = "0.10"
-Gridap = "0.16"
+Gridap = "0.17"
 LineSearches = "7.1"
 Sundials = "4.5"
 julia = "1.6"

README.md

@@ -48,9 +48,9 @@ a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ
 b(v,t) = ∫( v*f(t) )dΩ
 m(u,v) = ∫( v*u )dΩ
 
-res(t,u,ut,v) = a(u,v) + m(ut,v) - b(v,t)
-jac(t,u,ut,du,v) = a(du,v)
-jac_t(t,u,ut,dut,v) = m(dut,v)
+res(t,(u,ut),v) = a(u,v) + m(ut,v) - b(v,t)
+jac(t,(u,ut),du,v) = a(du,v)
+jac_t(t,(u,ut),dut,v) = m(dut,v)
 
 op = TransientFEOperator(res,jac,jac_t,U,V0)
 
@@ -62,9 +62,8 @@ U0 = U(0.0)
 uh0 = interpolate_everywhere(u(0.0),U0)
 
 ls = LUSolver()
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
-sol_t = solve(solver,op,uh0,t0,tF)
+ode_solver = ThetaMethod(ls,dt,θ)
+sol_t = solve(ode_solver,op,uh0,t0,tF)
 
 for (uh_tn, tn) in sol_t
   # Here we have the solution uh_tn at tn

src/DiffEqsWrappers/DiffEqsWrappers.jl

@@ -13,8 +13,7 @@ using GridapODEs.ODETools: jacobian!
 using Gridap.Algebra: allocate_jacobian
 
 using Gridap.FESpaces: get_algebraic_operator
-
-using Gridap.Algebra: fill_entries!
+using LinearAlgebra: fillstored!
 
 export prototype_jacobian
 export prototype_mass
@@ -51,21 +50,21 @@ function diffeq_wrappers(op)
   function _jacobian!(jac, du, u, p, gamma, t)
     ode_cache = update_cache!(ode_cache, ode_op, t)
     z = zero(eltype(jac))
-    fill_entries!(jac, z)
+    fillstored!(jac, z)
     jacobians!(jac, ode_op, t, (u, du), (1.0, gamma), ode_cache)
   end
 
   function _mass!(mass, du, u, p, t)
     ode_cache = update_cache!(ode_cache, ode_op, t)
     z = zero(eltype(mass))
-    fill_entries!(mass, z)
+    fillstored!(mass, z)
     jacobian!(mass, ode_op, t, (u, du), 2, 1.0, ode_cache)
   end
 
   function _stiffness!(stif, du, u, p, t)
     ode_cache = update_cache!(ode_cache, ode_op, t)
     z = zero(eltype(stif))
-    fill_entries!(stif, z)
+    fillstored!(stif, z)
     jacobian!(stif, ode_op, t, (u, du), 1, 1.0, ode_cache)
   end
 

src/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_entries!(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
 

src/ODETools/AffineThetaMethod.jl

@@ -103,13 +103,13 @@ end
 
 function _matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ)
   z = zero(eltype(A))
-  fill_entries!(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_entries!(A,z)
+  fillstored!(A,z)
   jacobian!(A,odeop,tθ,(vθ,vθ),2,(1/dtθ),ode_cache)
 end
 
@@ -140,7 +140,7 @@ function ThetaMethodConstantOperator(odeop::ConstantODEOperator,tθ::Float64,dt
   residual!(b,odeop,tθ,(u0,vθ),ode_cache)
   b = -1*b
   z = zero(eltype(A))
-  fill_entries!(A,z)
+  fillstored!(A,z)
   jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache)
   return A, b
 end

src/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_entries!(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
 

src/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_entries!(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_entries!(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_entries!(A,z)
+  fillstored!(A,z)
   jacobian!(A,op.odeop,op.t1,(u1,v1,a1),2,1.0,cache)
 end

src/ODETools/ForwardEuler.jl

@@ -60,7 +60,7 @@ function jacobian!(A::AbstractMatrix,op::ForwardEulerNonlinearOperator,x::Abstra
   vf = op.vf
   vf = (x-op.u0)/op.dt
   z = zero(eltype(A))
-  fill_entries!(A,z)
+  fillstored!(A,z)
   jacobian!(A,op.odeop,op.tf,(op.u0,vf),1,(1/op.dt),op.ode_cache)
 end
 

src/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_entries!(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
 

src/ODETools/ODETools.jl

@@ -3,6 +3,7 @@ module ODETools
 using Test
 
 using ForwardDiff
+using LinearAlgebra: fillstored!
 
 const ϵ = 100*eps()
 export ∂t
@@ -59,7 +60,6 @@ export MidPoint
 export ThetaMethod
 export RungeKutta
 export Newmark
-import Gridap.Algebra: fill_entries!
 
 export ODESolution
 export test_ode_solution

src/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_entries!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.ti,(ui,vi),(1/(op.a[op.i,op.i]*op.dt)),op.ode_cache)
 end
 

src/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_entries!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.tθ,(uF,vθ),(1.0,1/op.dtθ),op.ode_cache)
 end
 

src/TransientFETools/TransientFESolutions.jl

@@ -11,36 +11,51 @@ end
 
 
 function TransientFESolution(
-  solver::TransientFESolver, op::TransientFEOperator, uh0, t0::Real, tF::Real)
+  solver::ODESolver, op::TransientFEOperator, uh0, t0::Real, tF::Real)
 
-  ode_solver = solver.odes
   ode_op = get_algebraic_operator(op)
   u0 = get_free_dof_values(uh0)
-  ode_sol = solve(ode_solver,ode_op,u0,t0,tF)
+  ode_sol = solve(solver,ode_op,u0,t0,tF)
   trial = get_trial(op)
 
   TransientFESolution(ode_sol, trial)
 end
 
 function TransientFESolution(
-  solver::TransientFESolver,
+  solver::ODESolver,
   op::TransientFEOperator,
   xh0::Tuple{Vararg{Any}},
   t0::Real,
   tF::Real)
 
-  ode_solver = solver.odes
   ode_op = get_algebraic_operator(op)
   x0 = ()
   for xhi in xh0
     x0 = (x0...,get_free_dof_values(xhi))
   end
-  ode_sol = solve(ode_solver,ode_op,x0,t0,tF)
+  ode_sol = solve(solver,ode_op,x0,t0,tF)
   trial = get_trial(op)
 
   TransientFESolution(ode_sol, trial)
 end
 
+# Solve functions
+
+function solve(
+  solver::ODESolver,op::TransientFEOperator,u0,t0::Real,tf::Real)
+  TransientFESolution(solver,op,u0,t0,tf)
+end
+
+function solve(
+  solver::ODESolver,op::TransientFEOperator,u0,v0,a0,t0::Real,tf::Real)
+  TransientFESolution(solver,op,u0,v0,a0,t0,tf)
+end
+
+function test_transient_fe_solver(solver::ODESolver,op::TransientFEOperator,u0,t0,tf)
+  solution = solve(solver,op,u0,t0,tf)
+  test_transient_fe_solution(solution)
+end
+
 #@fverdugo this is a general implementation of iterate for TransientFESolution
 # We could also implement another one for the very common case that the
 # underlying ode_op is a ODEOpFromFEOp object

src/TransientFETools/TransientFETools.jl

@@ -67,7 +67,6 @@ import Gridap.FESpaces: get_test
 using GridapODEs.ODETools: test_ode_operator
 export test_transient_fe_operator
 
-export TransientFESolver
 import Gridap.FESpaces: FESolver
 import GridapODEs.ODETools: ODESolver
 import Gridap: solve
@@ -94,8 +93,6 @@ include("TransientFEOperators.jl")
 
 include("ODEOperatorInterfaces.jl")
 
-include("TransientFESolvers.jl")
-
 include("TransientFESolutions.jl")
 
 export FETerm

test/ODEsTests/ODESolverMocks.jl

@@ -1,6 +1,5 @@
 using Gridap.Algebra: residual
 using Gridap.Algebra: jacobian
-using Gridap.Algebra: fill_entries!
 import Gridap.Algebra: NonlinearSolver
 import Gridap.Algebra: NonlinearOperator
 import Gridap.Algebra: solve!
@@ -35,7 +34,7 @@ end
 function jacobian!(A::AbstractMatrix,op::OperatorMock,x::AbstractVector)
   uf = x
   uf_t = (x-op.u0)/op.dt
-  fill_entries!(A,0.0)
+  fillstored!(A,0.0)
   jacobian!(A,op.odeop,op.tf,uf,uf_t,op.cache)
   jacobian_t!(A,op.odeop,op.tf,uf,uf_t,(1/op.dt),op.cache)
 end

test/TransientFEsTests/AffineFEOperatorsTests.jl

@@ -54,9 +54,8 @@ U0 = U(0.0)
 uh0 = interpolate_everywhere(u(0.0),U0)
 
 ls = LUSolver()
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
-sol_t = solve(solver,op,uh0,t0,tF)
+ode_solver = ThetaMethod(ls,dt,θ)
+sol_t = solve(ode_solver,op,uh0,t0,tF)
 
 l2(w) = w*w
 

test/TransientFEsTests/BoundaryHeatEquationTests.jl

@@ -70,10 +70,9 @@ uh0 = interpolate_everywhere(u(0.0),U0)
 ls = LUSolver()
 using Gridap.Algebra: NewtonRaphsonSolver
 # nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
+ode_solver = ThetaMethod(ls,dt,θ)
 
-sol_t = solve(solver,op,uh0,t0,tF)
+sol_t = solve(ode_solver,op,uh0,t0,tF)
 
 # Juno.@enter Base.iterate(sol_t)
 

test/TransientFEsTests/ConstantFEOperatorsTests.jl

@@ -54,9 +54,8 @@ U0 = U(0.0)
 uh0 = interpolate_everywhere(u(0.0),U0)
 
 ls = LUSolver()
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
-sol_t = solve(solver,op,uh0,t0,tF)
+ode_solver = ThetaMethod(ls,dt,θ)
+sol_t = solve(ode_solver,op,uh0,t0,tF)
 
 l2(w) = w*w
 

test/TransientFEsTests/DGHeatEquationTests.jl

@@ -69,10 +69,9 @@ uh0 = interpolate_everywhere(u(0.0),U0)
 
 ls = LUSolver()
 using Gridap.Algebra: NewtonRaphsonSolver
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
+ode_solver = ThetaMethod(ls,dt,θ)
 
-sol_t = solve(solver,op,uh0,t0,tF)
+sol_t = solve(ode_solver,op,uh0,t0,tF)
 
 l2(w) = w*w
 

test/TransientFEsTests/ForwardEulerHeatEquationTests.jl

@@ -64,10 +64,9 @@ uh0 = interpolate_everywhere(u(0.0),U0)
 ls = LUSolver()
 using Gridap.Algebra: NewtonRaphsonSolver
 nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
+ode_solver = ThetaMethod(ls,dt,θ)
 
-sol_t = solve(solver,op,uh0,t0,tF)
+sol_t = solve(ode_solver,op,uh0,t0,tF)
 
 l2(w) = w*w
 

test/TransientFEsTests/HeatEquationAutoDiffTests.jl

@@ -71,10 +71,9 @@ uh0 = interpolate_everywhere(u(0.0),U0)
 ls = LUSolver()
 using Gridap.Algebra: NewtonRaphsonSolver
 nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
+ode_solver = ThetaMethod(ls,dt,θ)
 
-sol_t = solve(solver,op,uh0,t0,tF)
+sol_t = solve(ode_solver,op,uh0,t0,tF)
 
 # Juno.@enter Base.iterate(sol_t)