From 80ffe56804aff1cf81ad45114a72f9c50bff0c51 Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Wed, 8 Jul 2020 19:18:19 +1000
Subject: [PATCH 01/11] prototype

---
 Project.toml                                  |   1 +
 test/TransientFEsTests/DiffEqsWrapperTests.jl | 160 ++++++++++++++++++
 2 files changed, 161 insertions(+)
 create mode 100644 test/TransientFEsTests/DiffEqsWrapperTests.jl

diff --git a/Project.toml b/Project.toml
index dab8d81..d966ade 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
diff --git a/test/TransientFEsTests/DiffEqsWrapperTests.jl b/test/TransientFEsTests/DiffEqsWrapperTests.jl
new file mode 100644
index 0000000..814689a
--- /dev/null
+++ b/test/TransientFEsTests/DiffEqsWrapperTests.jl
@@ -0,0 +1,160 @@
+# module DiffEqsWrapperTests
+
+using Gridap
+using ForwardDiff
+using LinearAlgebra
+using Test
+using GridapODEs.ODETools
+using GridapODEs.TransientFETools
+using Gridap.FESpaces: get_algebraic_operator
+using GridapODEs
+using DifferentialEquations
+
+import Gridap: ∇
+import GridapODEs.TransientFETools: ∂t
+
+##
+
+θ = 1.0
+
+# @santiagobadia : u(x,t) = t not working with Sundials solver, not right sol?
+u(x,t) = 1.0#(1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
+
+# @santiagobadia : Even though the problem that we will solve is linear, the
+# Sundials solvers seems to have no convergence problems in the nonlinear solver
+
+u(x,t) = t#(1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
+
+
+
+u(t::Real) = x -> u(x,t)
+∂tu = ∂t(u)
+
+f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)
+
+domain = (0,1,0,1)
+partition = (4,4)
+model = CartesianDiscreteModel(domain,partition)
+
+order = 1
+
+V0 = FESpace(
+  reffe=:Lagrangian, order=order, valuetype=Float64,
+  conformity=:H1, model=model, dirichlet_tags="boundary")
+U = TransientTrialFESpace(V0,u)
+
+trian = Triangulation(model)
+degree = 2*order
+quad = CellQuadrature(trian,degree)
+
+#
+a(u,v) = ∇(v)⋅∇(u)
+b(v,t) = v*f(t)
+
+res(t,u,ut,v) = a(u,v) + ut*v - b(v,t)
+jac(t,u,ut,du,v) = a(du,v)
+jac_t(t,u,ut,dut,v) = dut*v
+
+t_Ω = FETerm(res,jac,jac_t,trian,quad)
+op = TransientFEOperator(U,V0,t_Ω)
+
+t0 = 0.0
+tF = 1.0
+dt = 0.1
+
+U0 = U(0.0)
+uh0 = interpolate_everywhere(U0,u(0.0))
+
+ls = LUSolver()
+using Gridap.Algebra: NewtonRaphsonSolver
+nls = NLSolver(ls;show_trace=true,method=:newton)
+odes = ThetaMethod(ls,dt,θ)
+solver = TransientFESolver(odes)
+sol_t = Gridap.solve(solver,op,uh0,t0,tF)
+
+#DiffEq wrapper
+
+ode_op = get_algebraic_operator(op)
+ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
+
+function residual(res,du,u,p,t)
+  # Problem with closure if ode_c is named ode_cache
+  # @santiagobadia : Clarify this point
+  ode_c = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
+  GridapODEs.ODETools.residual!(res,ode_op,t,u,du,ode_c)
+end
+
+function jacobian(jac,du,u,p,gamma,t)
+  ode_c = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
+  z = zero(eltype(jac))
+  Gridap.Algebra.fill_entries!(jac,z)
+  GridapODEs.ODETools.jacobian_t!(jac,ode_op,t,u,du,gamma,ode_c)
+  GridapODEs.ODETools.jacobian!(jac,ode_op,t,u,du,ode_c)
+end
+
+#end wrapper
+
+u0 = get_free_values(uh0)
+tθ = 1.0
+dtθ = dt*θ
+
+r = GridapODEs.ODETools.allocate_residual(ode_op,u0,ode_cache)
+J = GridapODEs.ODETools.allocate_jacobian(ode_op,u0,ode_cache)
+##
+residual(r,u0,u0,nothing,tθ)
+jacobian(J,u0,u0,nothing,(1/dtθ),tθ)
+
+using Sundials
+tspan = (0.0,1.0)
+
+# To explore the Sundials solver options, e.g., BE with time step dt
+f_iip = DAEFunction{true}(residual;jac=jacobian)
+prob_iip = DAEProblem{true}(f_iip,u0,u0,tspan,differential_vars=[true])
+sol_iip = Sundials.solve(prob_iip, IDA(), reltol=1e-8, abstol=1e-8)
+@show sol_iip.u
+
+# or iterator version
+integ = init(prob_iip, IDA(), reltol=1e-8, abstol=1e-8)
+step!(integ)
+step!(integ)
+
+# Show of using integrators as iterators
+using Base.Iterators
+for i in take(integ,100)
+      @show integ.u
+end
+
+# Issue 1: No control over the linear solver, we would like to be able to
+# provide certainly linear and possibly nonlinear solvers.
+# Let us assume that we just want to consider a fixed point algorithm
+# and we consider an implicit time integration of a nonlinear PDE.
+# Our solvers are efficient since they re-use cache among
+# nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization)
+# and for the transient case in our library, we can also reuse all this
+# between time steps. Could we attain something like this using
+# DifferentialEquations?
+
+# Issue 2: We have to call twice to update_cache!, for residual and jacobian
+# at each time step
+
+# Issue 3: Iterator-like version as in GridapODEs.
+# HOWEVER, is this computation lazy? I don't think so, so we need to store
+# all time steps, e.g., for computing time-dependent functionals (drag or lift
+# in CFD, etc), which is going to incur in a lot of memory consumption.
+
+
+l2(w) = w*w
+
+tol = 1.0e-6
+_t_n = t0
+
+for (uh_tn, tn) in sol_t
+  global _t_n
+  _t_n += dt
+  e = u(tn) - uh_tn
+  el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
+  @test el2 < tol
+  @show get_free_values(uh_tn)
+end
+
+# end #module

From d5bba91beb114bd4667f9a022bcee344cef5eb70 Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Wed, 8 Jul 2020 19:19:28 +1000
Subject: [PATCH 02/11] Delete DifEqWrappers.jl

---
 tmp/DiffEqWrappers/DifEqWrappers.jl | 177 ----------------------------
 1 file changed, 177 deletions(-)
 delete mode 100644 tmp/DiffEqWrappers/DifEqWrappers.jl

diff --git a/tmp/DiffEqWrappers/DifEqWrappers.jl b/tmp/DiffEqWrappers/DifEqWrappers.jl
deleted file mode 100644
index 843ec4e..0000000
--- a/tmp/DiffEqWrappers/DifEqWrappers.jl
+++ /dev/null
@@ -1,177 +0,0 @@
-using DifferentialEquations
-
-f(u,p,t) = 1.01*u
-u0 = 1/2
-tspan = (0.0,1.0)
-prob = ODEProblem(f,u0,tspan)
-sol = solve(prob, ImplicitEuler(), reltol=1e-8, abstol=1e-8)
-
-using Plots
-plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
-     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
-plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")
-
-
-# DAE problem
-u0 = zeros(1)
-u0[1] = 1/2
-
-function residual(res,du,u,p,t)
-  res[1] = - 1.01u[1] + du[1]
-end
-
-function jacobian(jac,du,u,p,gamma,t)
-  jac[1,1] = gamma - 1.01
-  nothing
-end
-
-function residual_ode(du,u,p,t)
-  du[1] = 1.01u[1]
-end
-
-function update_func(mass,u,p,t)
-  mass = ones(1,1)
-end
-
-dependent_M = DiffEqArrayOperator(ones(1,1),update_func=update_func)
-
-# Using https://docs.sciml.ai/v2.1.0/features/performance_overloads.html#Declaring-Explicit-Jacobians-1
-function residual_ode(::Type{Val{:jac}},t,u,J)
-  jac[1,1] = - 1.01
-  nothing
-end
-
-f_iip = ODEFunction{true}(residual_ode, mass_matrix=dependent_M)
-prob_iip = ODEProblem{true}(f_iip,u0,tspan)
-sol_iip = solve(prob_iip, ImplicitEuler(), reltol=1e-8, abstol=1e-8)
-
-# f_iip = ODEFunction{true}(residual;jac=jacobian)
-# Look here for more options
-#https://docs.sciml.ai/latest/features/performance_overloads/#
-# How can we provide a solver for W = M - gamma*J or W_t = gamma*W
-
-
-# Why is it looking for the OOP version of the Jacobian?
-
-Thanks for the quick reply, I will consider this option.
-
-In any case, following
-
-https://docs.sciml.ai/latest/features/performance_overloads/#
-
-I am creating a DAE problem
-
-```julia
-tspan = (0.0,1.0)
-u0 = zeros(1)
-u0[1] = 1/2
-
-function residual(res,du,u,p,t)
-  res[1] = - 1.01u[1] + du[1]
-end
-
-function jacobian(jac,du,u,p,gamma,t)
-  jac[1,1] = gamma - 1.01
-  nothing
-end
-
-f_iip = DAEFunction{true}(residual;jac=jacobian)
-prob_iip = DAEProblem{true}(f_iip,u0,u0,tspan)
-sol_iip = solve(prob_iip, DImplicitEuler(), reltol=1e-8, abstol=1e-8)
-```
-but it returns an error,
-```
-type ODEIntegrator has no field duprev
-```
-which is weird, since I would expect a `DAEIntegrator` being created.
-
-It would be great if you could tell me what I am doing wrong.
-
-Thanks
-
-# f_iip = ODEFunction{true}(residual)#,mass_matrix=M)
-# prob_iip = ODEProblem{true}(f_iip,u0,tspan)
-# sol_iip = solve(prob_iip, Rodas5(), reltol=1e-8, abstol=1e-8)
-
-# Juno.@enter solve(prob_iip, DImplicitEuler(), reltol=1e-8, abstol=1e-8)
-# I understand the problem, in place overwrites the du array, BUT not in terms of residual
-
-
-
-
-
-
-
-
-
-I have the following questions.
-
-I think that even if my "mass" matrix is non-singular, I would like to use the
-ADE interface you provide, since the mass matrix can be nonlinear, and thus, I
-need to compute its derivative with respect to `du` in the Jacobian. The DAE
-interface in `DifferentialEquations.jl` seems to provide this option and
-perfectly matches the interface in our PDE library `Gridap`.
-
-I would like to know whether using such interface instead of a ODE interface
-has an impact on performance (assuming we can solve the same problem in both
-cases).
-
-######
-
-I am trying to combine the ODE solvers in `DifferentialEquations.jl` with the PDE solvers in the `Gridap.jl` library.
-
-It is mandatory for our applications to use the inplace version of the `DifferentialEquations.jl` package. This way, we can pre-allocate the required arrays in our `Gridap` solvers.
-
-On the other hand, even when considering ODE systems, we want to use an implicit ODE expression of the solver, i.e., `A(du,u,p,t)=0`, since our _mass_ matrices are not just identity matrices, and can be fairly complicated or even nonlinear.
-
-Finally, what we would provide is a `residual(res,du,u,p,t)` and `jacobian(J,du,u,p,t)` functions.
-
-As a result, I think that the way to go is to consider an IIP problem in `DifferentialEquations.jl`. It seems that only the `DAE` constructor allows the syntax before.
-
-However, for implicit ODE systems I cannot use an `ODE` algorithm, e.g., Backward-Euler, even though it is acceptable for our problems when the _mass_ matrix is non-singular (thus not a `DAE`). And it seems we must provide an initial value for `du` that does not have mathematical sense in our problems.
-
-Summarizing, is there a way to create an `ODE` solver that takes `residual!(res,du,u,p,t)` and `jacobian!(J,du,u,p,t)` inplace functions, and allows me to use a `ODE` algorithm?
-
-I have written the following silly linear example with what I would like to have
-
-```julia
-tspan = (0.0,1.0)
-u0 = zeros(1)
-u0[1] = 1/2
-
-function residual(res,du,u,p,t)
-  res[1] = - 1.01u[1] + du[1]
-end
-
-function jacobian(jac,du,u,p,gamma,t)
-  jac[1,1] = gamma - 1.01
-  nothing
-end
-
-f_iip = ODEFunction{true}(residual;jac=jacobian)
-prob_iip = ODEProblem{true}(f_iip,u0,tspan)
-sol_iip = solve(prob_iip, ImplicitEuler(), reltol=1e-8, abstol=1e-8)
-```
-
-but it does not work, because it seems the library is looking for an `ODE`-style inplace version of `residual`, i.e., `residual(du,u,p,t)`, based on the error:
-
-```
-MethodError: no method matching residual(::Array{Float64,1}, ::Array{Float64,1}, ::DiffEqBase.NullParameters, ::Float64)
-Closest candidates are:
-  residual(::Any, ::Any, ::Any, ::Any, !Matched::Any) at /.../DifEqWrappers.jl:20
-...
-```
-
-In fact, even without considering the Jacobian, the same problem appears with the residual when doing this:
-```julia
-prob_iip = ODEProblem{false}(residual,u0,tspan)
-sol_iip = solve(prob_iip, ImplicitEuler(), reltol=1e-8, abstol=1e-8)
-```
-
-Is it possible to do what I want in `DifferentialEquations.jl`? If Yes, What am I doing wrong? I could not find an example in the documentation covering this (probably my fault).
-
-# On the other hand, is there a version of the solver that returns an iterator over the sequence of solutions in time?
-#
-# It would be a good feature, in complex PDEs one probably wants to consider the solution at each time step, perform
-# space, compute a norm, etc, but not to store the solution at all time steps at
-# the same time, since it is a lot of memory for a large PDE problem

From a73a73e47aaaef3bb2805bf53e7cf69f5e031ad8 Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Thu, 9 Jul 2020 19:14:46 +1000
Subject: [PATCH 03/11] sundials test

---
 test/TransientFEsTests/DiffEqsWrapperTests.jl | 98 ++++++++++---------
 1 file changed, 51 insertions(+), 47 deletions(-)

diff --git a/test/TransientFEsTests/DiffEqsWrapperTests.jl b/test/TransientFEsTests/DiffEqsWrapperTests.jl
index 814689a..e1e8083 100644
--- a/test/TransientFEsTests/DiffEqsWrapperTests.jl
+++ b/test/TransientFEsTests/DiffEqsWrapperTests.jl
@@ -1,4 +1,4 @@
-# module DiffEqsWrapperTests
+module DiffEqsWrapperTests
 
 using Gridap
 using ForwardDiff
@@ -13,27 +13,28 @@ using DifferentialEquations
 import Gridap: ∇
 import GridapODEs.TransientFETools: ∂t
 
-##
+# Heat equation
 
-θ = 1.0
-
-# @santiagobadia : u(x,t) = t not working with Sundials solver, not right sol?
-u(x,t) = 1.0#(1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
+# With this manufactured solution everything OK
+u(x,t) = 1.0
 
 # @santiagobadia : Even though the problem that we will solve is linear, the
-# Sundials solvers seems to have no convergence problems in the nonlinear solver
-
-u(x,t) = t#(1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
-
-
+# Sundials solvers seems to have convergence issues in the nonlinear solver (?)
+# It seems to work for partitions of at most (3,3), then it returns errors
+# [IDAS ERROR]  IDACalcIC
+# Newton/Linesearch algorithm failed to converge.
+# That is not surprising, because I guess that the J in jacobian is a
+# *dense* matrix, but not in Gridap (or any FEM solver). To clarify this
+# issue.
+u(x,t) = t
 
 u(t::Real) = x -> u(x,t)
-∂tu = ∂t(u)
+# ∂tu = ∂t(u)
 
 f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)
 
 domain = (0,1,0,1)
-partition = (4,4)
+partition = (3,3)
 model = CartesianDiscreteModel(domain,partition)
 
 order = 1
@@ -47,7 +48,6 @@ trian = Triangulation(model)
 degree = 2*order
 quad = CellQuadrature(trian,degree)
 
-#
 a(u,v) = ∇(v)⋅∇(u)
 b(v,t) = v*f(t)
 
@@ -65,21 +65,14 @@ dt = 0.1
 U0 = U(0.0)
 uh0 = interpolate_everywhere(U0,u(0.0))
 
-ls = LUSolver()
-using Gridap.Algebra: NewtonRaphsonSolver
-nls = NLSolver(ls;show_trace=true,method=:newton)
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
-sol_t = Gridap.solve(solver,op,uh0,t0,tF)
-
 #DiffEq wrapper
 
 ode_op = get_algebraic_operator(op)
 ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
 
 function residual(res,du,u,p,t)
-  # Problem with closure if ode_c is named ode_cache
-  # @santiagobadia : Clarify this point
+  # TO DO: Problem with closure if ode_c is named ode_cache
+  # TO DO: Improve update_cache! st do nothing if same time t as in the cache
   ode_c = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
   GridapODEs.ODETools.residual!(res,ode_op,t,u,du,ode_c)
 end
@@ -91,23 +84,25 @@ function jacobian(jac,du,u,p,gamma,t)
   GridapODEs.ODETools.jacobian_t!(jac,ode_op,t,u,du,gamma,ode_c)
   GridapODEs.ODETools.jacobian!(jac,ode_op,t,u,du,ode_c)
 end
-
 #end wrapper
 
-u0 = get_free_values(uh0)
+
+# Check
+θ = 1.0
 tθ = 1.0
+u0 = get_free_values(uh0)
 dtθ = dt*θ
-
 r = GridapODEs.ODETools.allocate_residual(ode_op,u0,ode_cache)
 J = GridapODEs.ODETools.allocate_jacobian(ode_op,u0,ode_cache)
-##
 residual(r,u0,u0,nothing,tθ)
 jacobian(J,u0,u0,nothing,(1/dtθ),tθ)
+#
 
 using Sundials
 tspan = (0.0,1.0)
 
-# To explore the Sundials solver options, e.g., BE with time step dt
+# To explore the Sundials solver options, e.g., BE with fixed time step dt
+# for comparison
 f_iip = DAEFunction{true}(residual;jac=jacobian)
 prob_iip = DAEProblem{true}(f_iip,u0,u0,tspan,differential_vars=[true])
 sol_iip = Sundials.solve(prob_iip, IDA(), reltol=1e-8, abstol=1e-8)
@@ -124,24 +119,14 @@ for i in take(integ,100)
       @show integ.u
 end
 
-# Issue 1: No control over the linear solver, we would like to be able to
-# provide certainly linear and possibly nonlinear solvers.
-# Let us assume that we just want to consider a fixed point algorithm
-# and we consider an implicit time integration of a nonlinear PDE.
-# Our solvers are efficient since they re-use cache among
-# nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization)
-# and for the transient case in our library, we can also reuse all this
-# between time steps. Could we attain something like this using
-# DifferentialEquations?
-
-# Issue 2: We have to call twice to update_cache!, for residual and jacobian
-# at each time step
-
-# Issue 3: Iterator-like version as in GridapODEs.
-# HOWEVER, is this computation lazy? I don't think so, so we need to store
-# all time steps, e.g., for computing time-dependent functionals (drag or lift
-# in CFD, etc), which is going to incur in a lot of memory consumption.
-
+# GridapODEs version using Backward Euler
+θ = 1.0
+ls = LUSolver()
+using Gridap.Algebra: NewtonRaphsonSolver
+nls = NLSolver(ls;show_trace=true,method=:newton)
+odes = ThetaMethod(ls,dt,θ)
+solver = TransientFESolver(odes)
+sol_t = Gridap.solve(solver,op,uh0,t0,tF)
 
 l2(w) = w*w
 
@@ -157,4 +142,23 @@ for (uh_tn, tn) in sol_t
   @show get_free_values(uh_tn)
 end
 
-# end #module
+# Issue 1: How do I initialize my residual vector (not that important) and my
+# jacobian (sparse CSR matrix). Dense matrices cannot be used in FE codes.
+
+# Issue 2: No control over the (non)linear solver, we would like to be able to
+# provide certainly linear and possibly nonlinear solvers.
+# Let us assume that we just want to consider a fixed point algorithm
+# and we consider an implicit time integration of a nonlinear PDE.
+# Our solvers are efficient since they re-use cache among
+# nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization)
+# and for the transient case in our library, we can also reuse all this
+# between time steps. Could we attain something like this using
+# DifferentialEquations/Sundials?
+
+
+# Issue 3: Iterator-like version as in GridapODEs.
+# HOWEVER, is this computation lazy? I don't think so, so we need to store
+# all time steps, e.g., for computing time-dependent functionals (drag or lift
+# in CFD, etc), which is going to incur in a lot of memory consumption.
+
+end #module

From 2beb716bd61617b5ddf25237b6adce9d5897423c Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Thu, 9 Jul 2020 19:33:05 +1000
Subject: [PATCH 04/11] fix

---
 test/TransientFEsTests/DiffEqsWrapperTests.jl | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/test/TransientFEsTests/DiffEqsWrapperTests.jl b/test/TransientFEsTests/DiffEqsWrapperTests.jl
index e1e8083..6ec05e3 100644
--- a/test/TransientFEsTests/DiffEqsWrapperTests.jl
+++ b/test/TransientFEsTests/DiffEqsWrapperTests.jl
@@ -71,18 +71,20 @@ ode_op = get_algebraic_operator(op)
 ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
 
 function residual(res,du,u,p,t)
-  # TO DO: Problem with closure if ode_c is named ode_cache
+  global ode_cache
   # TO DO: Improve update_cache! st do nothing if same time t as in the cache
-  ode_c = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
-  GridapODEs.ODETools.residual!(res,ode_op,t,u,du,ode_c)
+  # now it would be done twice (residual and jacobian)
+  ode_cache = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
+  GridapODEs.ODETools.residual!(res,ode_op,t,u,du,ode_cache)
 end
 
 function jacobian(jac,du,u,p,gamma,t)
-  ode_c = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
+  global ode_cache
+  ode_cache = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
   z = zero(eltype(jac))
   Gridap.Algebra.fill_entries!(jac,z)
-  GridapODEs.ODETools.jacobian_t!(jac,ode_op,t,u,du,gamma,ode_c)
-  GridapODEs.ODETools.jacobian!(jac,ode_op,t,u,du,ode_c)
+  GridapODEs.ODETools.jacobian_t!(jac,ode_op,t,u,du,gamma,ode_cache)
+  GridapODEs.ODETools.jacobian!(jac,ode_op,t,u,du,ode_cache)
 end
 #end wrapper
 

From 156e672a815257163cd38b62d04ed9712c655b06 Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Fri, 10 Jul 2020 22:08:10 +1000
Subject: [PATCH 05/11] minor changes

---
 test/TransientFEsTests/DiffEqsWrapperTests.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/TransientFEsTests/DiffEqsWrapperTests.jl b/test/TransientFEsTests/DiffEqsWrapperTests.jl
index 6ec05e3..c115eb2 100644
--- a/test/TransientFEsTests/DiffEqsWrapperTests.jl
+++ b/test/TransientFEsTests/DiffEqsWrapperTests.jl
@@ -83,8 +83,7 @@ function jacobian(jac,du,u,p,gamma,t)
   ode_cache = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
   z = zero(eltype(jac))
   Gridap.Algebra.fill_entries!(jac,z)
-  GridapODEs.ODETools.jacobian_t!(jac,ode_op,t,u,du,gamma,ode_cache)
-  GridapODEs.ODETools.jacobian!(jac,ode_op,t,u,du,ode_cache)
+  GridapODEs.ODETools.jacobian_and_jacobian_t!(jac,ode_op,t,u,du,gamma,ode_cache)
 end
 #end wrapper
 

From 8766cc79e536bfe61c04f6c2f95d0d86e95ec509 Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Fri, 24 Jul 2020 10:56:42 +1000
Subject: [PATCH 06/11] working on diffeq wrapper

---
 test/TransientFEsTests/DiffEqsWrapperTests.jl | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/test/TransientFEsTests/DiffEqsWrapperTests.jl b/test/TransientFEsTests/DiffEqsWrapperTests.jl
index c115eb2..14fde3c 100644
--- a/test/TransientFEsTests/DiffEqsWrapperTests.jl
+++ b/test/TransientFEsTests/DiffEqsWrapperTests.jl
@@ -97,6 +97,13 @@ r = GridapODEs.ODETools.allocate_residual(ode_op,u0,ode_cache)
 J = GridapODEs.ODETools.allocate_jacobian(ode_op,u0,ode_cache)
 residual(r,u0,u0,nothing,tθ)
 jacobian(J,u0,u0,nothing,(1/dtθ),tθ)
+
+
+ls = LUSolver()
+
+ls_cache = nothing
+x = copy(u0)
+solve!(x,J,r,ls_cache)
 #
 
 using Sundials
@@ -104,7 +111,7 @@ tspan = (0.0,1.0)
 
 # To explore the Sundials solver options, e.g., BE with fixed time step dt
 # for comparison
-f_iip = DAEFunction{true}(residual;jac=jacobian)
+f_iip = DAEFunction{true}(residual;jac=jacobian,jac_prototype=J)
 prob_iip = DAEProblem{true}(f_iip,u0,u0,tspan,differential_vars=[true])
 sol_iip = Sundials.solve(prob_iip, IDA(), reltol=1e-8, abstol=1e-8)
 @show sol_iip.u

From ae20b5271ac6fbed40375c351fb895a76970bffa Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Fri, 31 Jul 2020 15:53:45 +1000
Subject: [PATCH 07/11] DiffEq wrapper updated

---
 test/TransientFEsTests/DiffEqsWrapperTests.jl | 277 ++++++++++--------
 1 file changed, 155 insertions(+), 122 deletions(-)

diff --git a/test/TransientFEsTests/DiffEqsWrapperTests.jl b/test/TransientFEsTests/DiffEqsWrapperTests.jl
index 14fde3c..d9a57ed 100644
--- a/test/TransientFEsTests/DiffEqsWrapperTests.jl
+++ b/test/TransientFEsTests/DiffEqsWrapperTests.jl
@@ -1,160 +1,191 @@
-module DiffEqsWrapperTests
+module DiffEqsWrappersAux
 
-using Gridap
-using ForwardDiff
-using LinearAlgebra
 using Test
+using Gridap
 using GridapODEs.ODETools
 using GridapODEs.TransientFETools
 using Gridap.FESpaces: get_algebraic_operator
 using GridapODEs
-using DifferentialEquations
 
-import Gridap: ∇
-import GridapODEs.TransientFETools: ∂t
+export fe_problem
+export solve_ode_gridap
+export diffeq_wrappers
 
-# Heat equation
+function fe_problem(u, n)
 
-# With this manufactured solution everything OK
-u(x,t) = 1.0
+  f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x)
 
-# @santiagobadia : Even though the problem that we will solve is linear, the
-# Sundials solvers seems to have convergence issues in the nonlinear solver (?)
-# It seems to work for partitions of at most (3,3), then it returns errors
-# [IDAS ERROR]  IDACalcIC
-# Newton/Linesearch algorithm failed to converge.
-# That is not surprising, because I guess that the J in jacobian is a
-# *dense* matrix, but not in Gridap (or any FEM solver). To clarify this
-# issue.
-u(x,t) = t
+  domain = (0, 1, 0, 1)
+  partition = (n, n)
+  model = CartesianDiscreteModel(domain, partition)
 
-u(t::Real) = x -> u(x,t)
-# ∂tu = ∂t(u)
+  order = 1
 
-f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)
+  V0 = FESpace(
+    reffe = :Lagrangian,
+    order = order,
+    valuetype = Float64,
+    conformity = :H1,
+    model = model,
+    dirichlet_tags = "boundary",
+  )
+  U = TransientTrialFESpace(V0, u)
 
-domain = (0,1,0,1)
-partition = (3,3)
-model = CartesianDiscreteModel(domain,partition)
+  trian = Triangulation(model)
+  degree = 2 * order
+  quad = CellQuadrature(trian, degree)
 
-order = 1
+  a(u, v) = ∇(v) ⋅ ∇(u)
+  b(v, t) = v * f(t)
 
-V0 = FESpace(
-  reffe=:Lagrangian, order=order, valuetype=Float64,
-  conformity=:H1, model=model, dirichlet_tags="boundary")
-U = TransientTrialFESpace(V0,u)
+  res(t, u, ut, v) = a(u, v) + ut * v - b(v, t)
+  jac(t, u, ut, du, v) = a(du, v)
+  jac_t(t, u, ut, dut, v) = dut * v
 
-trian = Triangulation(model)
-degree = 2*order
-quad = CellQuadrature(trian,degree)
+  t_Ω = FETerm(res, jac, jac_t, trian, quad)
+  op = TransientFEOperator(U, V0, t_Ω)
 
-a(u,v) = ∇(v)⋅∇(u)
-b(v,t) = v*f(t)
+  U0 = U(0.0)
+  uh0 = interpolate_everywhere(U0, u(0.0))
 
-res(t,u,ut,v) = a(u,v) + ut*v - b(v,t)
-jac(t,u,ut,du,v) = a(du,v)
-jac_t(t,u,ut,dut,v) = dut*v
+  return op, trian, quad, uh0
 
-t_Ω = FETerm(res,jac,jac_t,trian,quad)
-op = TransientFEOperator(U,V0,t_Ω)
+end
 
-t0 = 0.0
-tF = 1.0
-dt = 0.1
+# ODE solver in GridapODEs
+function solve_ode_gridap(op, trian, quad, uh0, θ, dt, t0, tF, u)
+  # GridapODEs version using Backward Euler
+  ls = LUSolver()
+  nls = NLSolver(ls; show_trace = true, method = :newton)
+  odes = ThetaMethod(ls, dt, θ)
+  solver = TransientFESolver(odes)
+  sol_t = Gridap.solve(solver, op, uh0, t0, tF)
 
-U0 = U(0.0)
-uh0 = interpolate_everywhere(U0,u(0.0))
+  l2(w) = w * w
 
-#DiffEq wrapper
+  tol = 1.0e-6
+  _t_n = t0
 
-ode_op = get_algebraic_operator(op)
-ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
+  for (uh_tn, tn) in sol_t
+    _t_n += dt
+    e = u(tn) - uh_tn
+    el2 = sqrt(sum(integrate(l2(e), trian, quad)))
+    @test el2 < tol
+    @show get_free_values(uh_tn)
+  end
 
-function residual(res,du,u,p,t)
-  global ode_cache
-  # TO DO: Improve update_cache! st do nothing if same time t as in the cache
-  # now it would be done twice (residual and jacobian)
-  ode_cache = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
-  GridapODEs.ODETools.residual!(res,ode_op,t,u,du,ode_cache)
 end
 
-function jacobian(jac,du,u,p,gamma,t)
-  global ode_cache
-  ode_cache = GridapODEs.ODETools.update_cache!(ode_cache,ode_op,tθ)
-  z = zero(eltype(jac))
-  Gridap.Algebra.fill_entries!(jac,z)
-  GridapODEs.ODETools.jacobian_and_jacobian_t!(jac,ode_op,t,u,du,gamma,ode_cache)
+# Wrappers for residual and jacobian for DiffEqs
+function diffeq_wrappers(op)
+
+  ode_op = get_algebraic_operator(op)
+  ode_cache = allocate_cache(ode_op)
+
+  function residual(res, du, u, p, t)
+    # TO DO (minor): Improve update_cache! st do nothing if same time t as in the cache
+    # now it would be done twice (residual and jacobian)
+    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
+    GridapODEs.ODETools.residual!(res, ode_op, t, u, du, ode_cache)
+  end
+
+  function jacobian(jac, du, u, p, gamma, t)
+    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
+    z = zero(eltype(jac))
+    Gridap.Algebra.fill_entries!(jac, z)
+    GridapODEs.ODETools.jacobian_and_jacobian_t!(jac, ode_op, t, u, du, gamma, ode_cache)
+  end
+
+  return residual, jacobian
+
 end
-#end wrapper
+
+end #module
 
 
-# Check
+module DiffEqsWrapperTests
+
+using Gridap
+using ForwardDiff
+using LinearAlgebra
+using Test
+using GridapODEs.ODETools
+using GridapODEs.TransientFETools
+using Gridap.FESpaces: get_algebraic_operator
+using GridapODEs
+using DifferentialEquations
+using Gridap.Algebra: NewtonRaphsonSolver
+using Base.Iterators
+
+# Solving the heat equation using GridapODEs and DiffEqs
+
+u(x, t) = t#1.0
+u(t) = x -> u(x, t)
+
+# FE structs
+op, trian, quad, uh0 = DiffEqsWrappersAux.fe_problem(u,3)
+
 θ = 1.0
-tθ = 1.0
-u0 = get_free_values(uh0)
-dtθ = dt*θ
-r = GridapODEs.ODETools.allocate_residual(ode_op,u0,ode_cache)
-J = GridapODEs.ODETools.allocate_jacobian(ode_op,u0,ode_cache)
-residual(r,u0,u0,nothing,tθ)
-jacobian(J,u0,u0,nothing,(1/dtθ),tθ)
+t0 = 0.0
+tF = 1.0
+dt = 0.1
 
+# Check Gridap ODE solver
+DiffEqsWrappersAux.solve_ode_gridap(op,trian,quad,uh0,θ,dt,t0,tF,u)
 
-ls = LUSolver()
+# DiffEq wrappers
+residual, jacobian = DiffEqsWrappersAux.diffeq_wrappers(op)
 
-ls_cache = nothing
-x = copy(u0)
-solve!(x,J,r,ls_cache)
-#
+# Check wrappers
+u0 = get_free_values(uh0)
+dtθ = dt * θ
+ode_op = get_algebraic_operator(op)
+ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
+r = GridapODEs.ODETools.allocate_residual(ode_op, u0, ode_cache)
+J = GridapODEs.ODETools.allocate_jacobian(ode_op, u0, ode_cache)
 
+tθ = 1.0
+residual(r, u0, u0, nothing, tθ)
+jacobian(J, u0, u0, nothing, (1 / dtθ), tθ)
+
+J
+
+# Using Sundials to solve the resulting ODE
 using Sundials
-tspan = (0.0,1.0)
+tspan = (0.0, 1.0)
+
+u(x, t) = t
+u(t) = x -> u(x, t)
+
+# ISSUE 1: When I choose n > 2, even though the problem that we will solve is
+# linear, the Sundials solvers seems to have convergence issues in the nonlinear
+# solver (?). Ut returns errors
+# [IDAS ERROR]  IDACalcIC Newton/Linesearch algorithm failed to converge.
+# ISSUE 2: When I pass `jac_prototype` the code gets stuck
+n = 2
+op, trian, quad, uh0 = DiffEqsWrappersAux.fe_problem(u,n)
+residual, jacobian = DiffEqsWrappersAux.diffeq_wrappers(op)
+u0 = get_free_values(uh0)
 
-# To explore the Sundials solver options, e.g., BE with fixed time step dt
-# for comparison
-f_iip = DAEFunction{true}(residual;jac=jacobian,jac_prototype=J)
-prob_iip = DAEProblem{true}(f_iip,u0,u0,tspan,differential_vars=[true])
-sol_iip = Sundials.solve(prob_iip, IDA(), reltol=1e-8, abstol=1e-8)
+# To explore the Sundials solver options, e.g., BE with fixed time step dtd
+f_iip = DAEFunction{true}(residual; jac = jacobian) #,jac_prototype=J)
+# jac_prototype is the way to pass my pre-allocated jacobian matrix
+prob_iip = DAEProblem{true}(f_iip, u0, u0, tspan, differential_vars = [true])
+sol_iip = Sundials.solve(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
 @show sol_iip.u
 
 # or iterator version
-integ = init(prob_iip, IDA(), reltol=1e-8, abstol=1e-8)
-step!(integ)
-step!(integ)
-
-# Show of using integrators as iterators
-using Base.Iterators
-for i in take(integ,100)
-      @show integ.u
-end
+integ = init(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
+# step!(integ)
 
-# GridapODEs version using Backward Euler
-θ = 1.0
-ls = LUSolver()
-using Gridap.Algebra: NewtonRaphsonSolver
-nls = NLSolver(ls;show_trace=true,method=:newton)
-odes = ThetaMethod(ls,dt,θ)
-solver = TransientFESolver(odes)
-sol_t = Gridap.solve(solver,op,uh0,t0,tF)
-
-l2(w) = w*w
-
-tol = 1.0e-6
-_t_n = t0
-
-for (uh_tn, tn) in sol_t
-  global _t_n
-  _t_n += dt
-  e = u(tn) - uh_tn
-  el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
-  @test el2 < tol
-  @show get_free_values(uh_tn)
+# Show using integrators as iterators
+for i in take(integ, 100)
+  @show integ.u
 end
 
-# Issue 1: How do I initialize my residual vector (not that important) and my
-# jacobian (sparse CSR matrix). Dense matrices cannot be used in FE codes.
+end # module
 
-# Issue 2: No control over the (non)linear solver, we would like to be able to
-# provide certainly linear and possibly nonlinear solvers.
+# ISSUE: Future work, add own (non)linear solver.
 # Let us assume that we just want to consider a fixed point algorithm
 # and we consider an implicit time integration of a nonlinear PDE.
 # Our solvers are efficient since they re-use cache among
@@ -162,11 +193,13 @@ end
 # and for the transient case in our library, we can also reuse all this
 # between time steps. Could we attain something like this using
 # DifferentialEquations/Sundials?
-
-
-# Issue 3: Iterator-like version as in GridapODEs.
-# HOWEVER, is this computation lazy? I don't think so, so we need to store
-# all time steps, e.g., for computing time-dependent functionals (drag or lift
-# in CFD, etc), which is going to incur in a lot of memory consumption.
-
-end #module
+# @ChrisRackauckas suggests to take a look at:
+# https://docs.sciml.ai/latest/tutorials/advanced_ode_example/ shows swapping out linear solvers.
+# https://docs.sciml.ai/latest/features/linear_nonlinear/ is all of the extra details.
+
+# Try to pass solver too
+# ls = LUSolver()
+# ls_cache = nothing
+# x = copy(u0)
+# solve!(x,J,r,ls_cache)
+#

From dbd898e7e98d00a00a4aaeb17d848acf6248767f Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Tue, 4 Aug 2020 11:16:28 +1000
Subject: [PATCH 08/11] Added simple test

---
 test/TransientFEsTests/SimpleDiffEqsTest.jl | 161 ++++++++++++++++++++
 1 file changed, 161 insertions(+)
 create mode 100644 test/TransientFEsTests/SimpleDiffEqsTest.jl

diff --git a/test/TransientFEsTests/SimpleDiffEqsTest.jl b/test/TransientFEsTests/SimpleDiffEqsTest.jl
new file mode 100644
index 0000000..12b14b2
--- /dev/null
+++ b/test/TransientFEsTests/SimpleDiffEqsTest.jl
@@ -0,0 +1,161 @@
+module GridapFESolver
+
+using Test
+using Gridap
+using GridapODEs.ODETools
+using GridapODEs.TransientFETools
+using Gridap.FESpaces: get_algebraic_operator
+using GridapODEs
+
+export fe_problem
+export solve_ode_gridap
+export diffeq_wrappers
+
+function fe_problem(u, n)
+
+  f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x)
+
+  domain = (0, 1, 0, 1)
+  partition = (n, n)
+  model = CartesianDiscreteModel(domain, partition)
+
+  order = 1
+
+  V0 = FESpace(
+    reffe = :Lagrangian,
+    order = order,
+    valuetype = Float64,
+    conformity = :H1,
+    model = model,
+    dirichlet_tags = "boundary",
+  )
+  U = TransientTrialFESpace(V0, u)
+
+  trian = Triangulation(model)
+  degree = 2 * order
+  quad = CellQuadrature(trian, degree)
+
+  a(u, v) = ∇(v) ⋅ ∇(u)
+  b(v, t) = v * f(t)
+
+  res(t, u, ut, v) = a(u, v) + ut * v - b(v, t)
+  jac(t, u, ut, du, v) = a(du, v)
+  jac_t(t, u, ut, dut, v) = dut * v
+
+  t_Ω = FETerm(res, jac, jac_t, trian, quad)
+  op = TransientFEOperator(U, V0, t_Ω)
+
+  U0 = U(0.0)
+  uh0 = interpolate_everywhere(U0, u(0.0))
+  u0 = get_free_values(uh0)
+
+  return op, u0
+
+end
+
+# Wrappers for residual and jacobian for DiffEqs
+function diffeq_wrappers(op)
+
+  ode_op = get_algebraic_operator(op)
+  ode_cache = allocate_cache(ode_op)
+
+  function residual(res, du, u, p, t)
+    # TO DO (minor): Improve update_cache! st do nothing if same time t as in the cache
+    # now it would be done twice (residual and jacobian)
+    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
+    GridapODEs.ODETools.residual!(res, ode_op, t, u, du, ode_cache)
+  end
+
+  function jacobian(jac, du, u, p, gamma, t)
+    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
+    z = zero(eltype(jac))
+    Gridap.Algebra.fill_entries!(jac, z)
+    GridapODEs.ODETools.jacobian_and_jacobian_t!(jac, ode_op, t, u, du, gamma, ode_cache)
+  end
+
+  return residual, jacobian
+
+end
+
+# Allocate Jacobian
+function allocate_jac(op,u0)
+  ode_op = get_algebraic_operator(op)
+  ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
+  return GridapODEs.ODETools.allocate_jacobian(ode_op, u0, ode_cache)
+end
+
+end #module
+
+
+module DiffEqsWrapperTests
+
+using Gridap
+using ForwardDiff
+using LinearAlgebra
+using Test
+using GridapODEs.ODETools
+using GridapODEs.TransientFETools
+using Gridap.FESpaces: get_algebraic_operator
+using GridapODEs
+using DifferentialEquations
+using Gridap.Algebra: NewtonRaphsonSolver
+using Base.Iterators
+
+# Solving the heat equation using GridapODEs and DiffEqs
+using Sundials
+tspan = (0.0, 1.0)
+
+u(x, t) = t
+u(t) = x -> u(x, t)
+
+# ISSUE 1: When I choose n > 2, even though the problem that we will solve is
+# linear, the Sundials solvers seems to have convergence issues in the nonlinear
+# solver (?). Ut returns errors
+# [IDAS ERROR]  IDACalcIC Newton/Linesearch algorithm failed to converge.
+# ISSUE 2: When I pass `jac_prototype` the code gets stuck
+S = GridapFESolver
+
+n = 9
+op, u0 = S.fe_problem(u,n)
+residual, jacobian = S.diffeq_wrappers(op)
+J = S.allocate_jac(op,u0)
+
+# To explore the Sundials solver options, e.g., BE with fixed time step dtd
+f_iip = DAEFunction{true}(residual; jac = jacobian)#, jac_prototype=J)
+# jac_prototype is the way to pass my pre-allocated jacobian matrix
+prob_iip = DAEProblem{true}(f_iip, u0, u0, tspan, differential_vars = [true])
+# When I pass `jac_prototype` the code get stuck here:
+sol_iip = Sundials.solve(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
+@show sol_iip.u
+
+# or iterator version
+integ = init(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
+# step!(integ)
+
+# Show using integrators as iterators
+for i in take(integ, 100)
+  @show integ.u
+end
+
+end # module
+
+# FUTURE WORK: Check other options, not only Sundials
+
+# ISSUE: Future work, add own (non)linear solver.
+# Let us assume that we just want to consider a fixed point algorithm
+# and we consider an implicit time integration of a nonlinear PDE.
+# Our solvers are efficient since they re-use cache among
+# nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization)
+# and for the transient case in our library, we can also reuse all this
+# between time steps. Could we attain something like this using
+# DifferentialEquations/Sundials?
+# @ChrisRackauckas suggests to take a look at:
+# https://docs.sciml.ai/latest/tutorials/advanced_ode_example/ shows swapping out linear solvers.
+# https://docs.sciml.ai/latest/features/linear_nonlinear/ is all of the extra details.
+
+# Try to pass solver too
+# ls = LUSolver()
+# ls_cache = nothing
+# x = copy(u0)
+# solve!(x,J,r,ls_cache)
+#

From 768820785095ad9371f9e30fa2febc9e2f39f936 Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Fri, 7 Aug 2020 22:07:22 +1000
Subject: [PATCH 09/11] Mass matrix in DiffEq

---
 test/TransientFEsTests/SimpleDiffEqsTest.jl | 37 ++++++++++++++++++---
 1 file changed, 33 insertions(+), 4 deletions(-)

diff --git a/test/TransientFEsTests/SimpleDiffEqsTest.jl b/test/TransientFEsTests/SimpleDiffEqsTest.jl
index 12b14b2..60e3dd2 100644
--- a/test/TransientFEsTests/SimpleDiffEqsTest.jl
+++ b/test/TransientFEsTests/SimpleDiffEqsTest.jl
@@ -11,6 +11,7 @@ export fe_problem
 export solve_ode_gridap
 export diffeq_wrappers
 
+# FE problem (heat eq) using Gridap
 function fe_problem(u, n)
 
   f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x)
@@ -73,7 +74,21 @@ function diffeq_wrappers(op)
     GridapODEs.ODETools.jacobian_and_jacobian_t!(jac, ode_op, t, u, du, gamma, ode_cache)
   end
 
-  return residual, jacobian
+  function mass(mass, du, u, p, t)
+    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
+    z = zero(eltype(mass))
+    Gridap.Algebra.fill_entries!(mass, z)
+    GridapODEs.ODETools.jacobian_t!(mass, ode_op, t, u, du, 1.0, ode_cache)
+  end
+
+  function stiffness(stif, du, u, p, t)
+    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
+    z = zero(eltype(stif))
+    Gridap.Algebra.fill_entries!(stif, z)
+    GridapODEs.ODETools.jacobian!(stif, ode_op, t, u, du, ode_cache)
+  end
+
+  return residual, jacobian, mass, stiffness
 
 end
 
@@ -113,12 +128,26 @@ u(t) = x -> u(x, t)
 # solver (?). Ut returns errors
 # [IDAS ERROR]  IDACalcIC Newton/Linesearch algorithm failed to converge.
 # ISSUE 2: When I pass `jac_prototype` the code gets stuck
-S = GridapFESolver
+S = Main.GridapFESolver
 
-n = 9
+n = 3 # cells per dim (2D)
 op, u0 = S.fe_problem(u,n)
-residual, jacobian = S.diffeq_wrappers(op)
+
+# Some checks
+residual, jacobian, mass, stiffness = S.diffeq_wrappers(op)
 J = S.allocate_jac(op,u0)
+r = copy(u0)
+θ = 1.0; t0 = 0.0; tF = 1.0; dt = 0.1; tθ = 1.0; dtθ = dt*θ
+residual(r, u0, u0, nothing, tθ)
+jacobian(J, u0, u0, nothing, (1 / dtθ), tθ)
+
+K = S.allocate_jac(op,u0)
+M = S.allocate_jac(op,u0)
+stiffness(K, u0, u0, nothing, tθ)
+mass(M, u0, u0, nothing, tθ)
+# Here you have the mass matrix M
+
+@test (1/dtθ)*M+K ≈ J
 
 # To explore the Sundials solver options, e.g., BE with fixed time step dtd
 f_iip = DAEFunction{true}(residual; jac = jacobian)#, jac_prototype=J)

From f25370468525d625bfddd818952eb1b472ba7df0 Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Thu, 13 Aug 2020 14:47:28 +1000
Subject: [PATCH 10/11] changing solver sundials

---
 test/TransientFEsTests/DiffEqsWrapperTests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/TransientFEsTests/DiffEqsWrapperTests.jl b/test/TransientFEsTests/DiffEqsWrapperTests.jl
index d9a57ed..87416a3 100644
--- a/test/TransientFEsTests/DiffEqsWrapperTests.jl
+++ b/test/TransientFEsTests/DiffEqsWrapperTests.jl
@@ -171,7 +171,7 @@ u0 = get_free_values(uh0)
 f_iip = DAEFunction{true}(residual; jac = jacobian) #,jac_prototype=J)
 # jac_prototype is the way to pass my pre-allocated jacobian matrix
 prob_iip = DAEProblem{true}(f_iip, u0, u0, tspan, differential_vars = [true])
-sol_iip = Sundials.solve(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
+sol_iip = Sundials.solve(prob_iip, IDA(linear_solver=:KLU), reltol = 1e-8, abstol = 1e-8)
 @show sol_iip.u
 
 # or iterator version

From e067ab9ab15935fff820873de01c17b334d4e334 Mon Sep 17 00:00:00 2001
From: Santiago Badia <santiago.badia@monash.edu>
Date: Thu, 13 Aug 2020 22:27:38 +1000
Subject: [PATCH 11/11] added wrappers for sundials

---
 Project.toml                                  |   2 +-
 src/DiffEqsWrappers/DiffEqsWrappers.jl        |  91 ++++++++
 src/GridapODEs.jl                             |   2 +-
 .../DiffEqsTests.jl}                          | 101 ++-------
 test/DiffEqsWrappersTests/runtests.jl         |   7 +
 test/TransientFEsTests/DiffEqsWrapperTests.jl | 205 ------------------
 test/runtests.jl                              |   2 +
 7 files changed, 121 insertions(+), 289 deletions(-)
 create mode 100644 src/DiffEqsWrappers/DiffEqsWrappers.jl
 rename test/{TransientFEsTests/SimpleDiffEqsTest.jl => DiffEqsWrappersTests/DiffEqsTests.jl} (56%)
 create mode 100644 test/DiffEqsWrappersTests/runtests.jl
 delete mode 100644 test/TransientFEsTests/DiffEqsWrapperTests.jl

diff --git a/Project.toml b/Project.toml
index d966ade..f570b18 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 ForwardDiff = "0.10.10"
-Gridap = "0.12"
+Gridap = "0.12, 0.13"
 LineSearches = "7.0.1"
 julia = "1.0"
 
diff --git a/src/DiffEqsWrappers/DiffEqsWrappers.jl b/src/DiffEqsWrappers/DiffEqsWrappers.jl
new file mode 100644
index 0000000..580a778
--- /dev/null
+++ b/src/DiffEqsWrappers/DiffEqsWrappers.jl
@@ -0,0 +1,91 @@
+module DiffEqWrappers
+
+using Test
+
+using GridapODEs.TransientFETools: TransientFEOperator
+
+using GridapODEs.ODETools: allocate_cache
+using GridapODEs.ODETools: update_cache!
+using GridapODEs.ODETools: residual!
+using GridapODEs.ODETools: jacobian_and_jacobian_t!
+using GridapODEs.ODETools: jacobian!
+using GridapODEs.ODETools: jacobian_t!
+
+using Gridap.Algebra: allocate_jacobian
+
+using Gridap.FESpaces: get_algebraic_operator
+
+using Gridap.Algebra: fill_entries!
+
+export prototype_jacobian
+export prototype_mass
+export prototype_stiffness
+
+export diffeq_wrappers
+
+"""
+  This method takes a `FEOperator` and returns some methods that can be used
+  in `DifferentialEquations.jl`. Assuming we want to solve the nonlinear ODE
+  `res(t,u,du) = 0`, we return:
+  1. `residual!(res, du, u, p, t)`: It returns the residual (in `res`) at
+     `(u,du,t)`, following the signature in `DifferentialEquations`.
+     For the moment, we do not support parameters.
+  2. `jacobian!(jac, du, u, p, gamma, t)`: Idem for the Jacobian. It returns
+     `∂res/∂du*γ + ∂res/∂u`
+  3. `mass!(mass, du, u, p, t)`: Idem for the mass matrix. It returns
+     `∂res/∂du`
+  4. `stiffness!(stif, du, u, p, t)`: Idem for the mass matrix. It returns
+     `∂res/∂u`
+"""
+function diffeq_wrappers(op)
+
+  ode_op = get_algebraic_operator(op)
+  ode_cache = allocate_cache(ode_op)
+
+  function _residual!(res, du, u, p, t)
+    # TO DO (minor): Improve update_cache! st do nothing if same time t as in the cache
+    # now it would be done twice (residual and jacobian)
+    ode_cache = update_cache!(ode_cache, ode_op, t)
+    residual!(res, ode_op, t, u, du, ode_cache)
+  end
+
+  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)
+    jacobian_and_jacobian_t!(jac, ode_op, t, u, du, 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)
+    jacobian_t!(mass, ode_op, t, u, du, 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)
+    jacobian!(stif, ode_op, t, u, du, ode_cache)
+  end
+
+  return _residual!, _jacobian!, _mass!, _stiffness!
+
+end
+
+"""
+  It allocates the Jacobian (or mass or stiffness) matrix, given the `FEOperator`
+  and a vector of size total number of unknowns
+"""
+function prototype_jacobian(op::TransientFEOperator,u0)
+  ode_op = get_algebraic_operator(op)
+  ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
+  return allocate_jacobian(ode_op, u0, ode_cache)
+end
+
+const prototype_mass = prototype_jacobian
+
+const prototype_stiffness = prototype_jacobian
+
+end #module
diff --git a/src/GridapODEs.jl b/src/GridapODEs.jl
index df068de..97a096b 100644
--- a/src/GridapODEs.jl
+++ b/src/GridapODEs.jl
@@ -4,6 +4,6 @@ include("ODETools/ODETools.jl")
 
 include("TransientFETools/TransientFETools.jl")
 
-# include("Exports.jl")
+include("DiffEqsWrappers/DiffEqsWrappers.jl")
 
 end #module
diff --git a/test/TransientFEsTests/SimpleDiffEqsTest.jl b/test/DiffEqsWrappersTests/DiffEqsTests.jl
similarity index 56%
rename from test/TransientFEsTests/SimpleDiffEqsTest.jl
rename to test/DiffEqsWrappersTests/DiffEqsTests.jl
index 60e3dd2..36a750f 100644
--- a/test/TransientFEsTests/SimpleDiffEqsTest.jl
+++ b/test/DiffEqsWrappersTests/DiffEqsTests.jl
@@ -1,15 +1,16 @@
-module GridapFESolver
+module DiffEqsWrapperTests
 
 using Test
 using Gridap
+using GridapODEs
 using GridapODEs.ODETools
 using GridapODEs.TransientFETools
-using Gridap.FESpaces: get_algebraic_operator
-using GridapODEs
+using GridapODEs.DiffEqWrappers
 
-export fe_problem
-export solve_ode_gridap
-export diffeq_wrappers
+# using DifferentialEquations
+using Sundials
+using Gridap.Algebra: NewtonRaphsonSolver
+using Base.Iterators
 
 # FE problem (heat eq) using Gridap
 function fe_problem(u, n)
@@ -54,70 +55,7 @@ function fe_problem(u, n)
 
 end
 
-# Wrappers for residual and jacobian for DiffEqs
-function diffeq_wrappers(op)
-
-  ode_op = get_algebraic_operator(op)
-  ode_cache = allocate_cache(ode_op)
-
-  function residual(res, du, u, p, t)
-    # TO DO (minor): Improve update_cache! st do nothing if same time t as in the cache
-    # now it would be done twice (residual and jacobian)
-    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
-    GridapODEs.ODETools.residual!(res, ode_op, t, u, du, ode_cache)
-  end
-
-  function jacobian(jac, du, u, p, gamma, t)
-    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
-    z = zero(eltype(jac))
-    Gridap.Algebra.fill_entries!(jac, z)
-    GridapODEs.ODETools.jacobian_and_jacobian_t!(jac, ode_op, t, u, du, gamma, ode_cache)
-  end
-
-  function mass(mass, du, u, p, t)
-    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
-    z = zero(eltype(mass))
-    Gridap.Algebra.fill_entries!(mass, z)
-    GridapODEs.ODETools.jacobian_t!(mass, ode_op, t, u, du, 1.0, ode_cache)
-  end
-
-  function stiffness(stif, du, u, p, t)
-    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
-    z = zero(eltype(stif))
-    Gridap.Algebra.fill_entries!(stif, z)
-    GridapODEs.ODETools.jacobian!(stif, ode_op, t, u, du, ode_cache)
-  end
-
-  return residual, jacobian, mass, stiffness
-
-end
-
-# Allocate Jacobian
-function allocate_jac(op,u0)
-  ode_op = get_algebraic_operator(op)
-  ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
-  return GridapODEs.ODETools.allocate_jacobian(ode_op, u0, ode_cache)
-end
-
-end #module
-
-
-module DiffEqsWrapperTests
-
-using Gridap
-using ForwardDiff
-using LinearAlgebra
-using Test
-using GridapODEs.ODETools
-using GridapODEs.TransientFETools
-using Gridap.FESpaces: get_algebraic_operator
-using GridapODEs
-using DifferentialEquations
-using Gridap.Algebra: NewtonRaphsonSolver
-using Base.Iterators
-
 # Solving the heat equation using GridapODEs and DiffEqs
-using Sundials
 tspan = (0.0, 1.0)
 
 u(x, t) = t
@@ -128,34 +66,33 @@ u(t) = x -> u(x, t)
 # solver (?). Ut returns errors
 # [IDAS ERROR]  IDACalcIC Newton/Linesearch algorithm failed to converge.
 # ISSUE 2: When I pass `jac_prototype` the code gets stuck
-S = Main.GridapFESolver
 
 n = 3 # cells per dim (2D)
-op, u0 = S.fe_problem(u,n)
+op, u0 = fe_problem(u,n)
 
 # Some checks
-residual, jacobian, mass, stiffness = S.diffeq_wrappers(op)
-J = S.allocate_jac(op,u0)
+res!, jac!, mass!, stif! = diffeq_wrappers(op)
+J = prototype_jacobian(op,u0)
 r = copy(u0)
 θ = 1.0; t0 = 0.0; tF = 1.0; dt = 0.1; tθ = 1.0; dtθ = dt*θ
-residual(r, u0, u0, nothing, tθ)
-jacobian(J, u0, u0, nothing, (1 / dtθ), tθ)
+res!(r, u0, u0, nothing, tθ)
+jac!(J, u0, u0, nothing, (1 / dtθ), tθ)
 
-K = S.allocate_jac(op,u0)
-M = S.allocate_jac(op,u0)
-stiffness(K, u0, u0, nothing, tθ)
-mass(M, u0, u0, nothing, tθ)
+K = prototype_jacobian(op,u0)
+M = prototype_jacobian(op,u0)
+stif!(K, u0, u0, nothing, tθ)
+mass!(M, u0, u0, nothing, tθ)
 # Here you have the mass matrix M
 
 @test (1/dtθ)*M+K ≈ J
 
 # To explore the Sundials solver options, e.g., BE with fixed time step dtd
-f_iip = DAEFunction{true}(residual; jac = jacobian)#, jac_prototype=J)
+f_iip = DAEFunction{true}(res!; jac = jac!)#, jac_prototype=J)
 # jac_prototype is the way to pass my pre-allocated jacobian matrix
 prob_iip = DAEProblem{true}(f_iip, u0, u0, tspan, differential_vars = [true])
 # When I pass `jac_prototype` the code get stuck here:
 sol_iip = Sundials.solve(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
-@show sol_iip.u
+# @show sol_iip.u
 
 # or iterator version
 integ = init(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
@@ -163,7 +100,7 @@ integ = init(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
 
 # Show using integrators as iterators
 for i in take(integ, 100)
-  @show integ.u
+  # @show integ.u
 end
 
 end # module
diff --git a/test/DiffEqsWrappersTests/runtests.jl b/test/DiffEqsWrappersTests/runtests.jl
new file mode 100644
index 0000000..fd0d0a1
--- /dev/null
+++ b/test/DiffEqsWrappersTests/runtests.jl
@@ -0,0 +1,7 @@
+module DiffEqsWrappersTests
+
+using Test
+
+@testset "DiffEqWrappers" begin include("DiffEqsTests.jl") end
+
+end # module
diff --git a/test/TransientFEsTests/DiffEqsWrapperTests.jl b/test/TransientFEsTests/DiffEqsWrapperTests.jl
deleted file mode 100644
index 87416a3..0000000
--- a/test/TransientFEsTests/DiffEqsWrapperTests.jl
+++ /dev/null
@@ -1,205 +0,0 @@
-module DiffEqsWrappersAux
-
-using Test
-using Gridap
-using GridapODEs.ODETools
-using GridapODEs.TransientFETools
-using Gridap.FESpaces: get_algebraic_operator
-using GridapODEs
-
-export fe_problem
-export solve_ode_gridap
-export diffeq_wrappers
-
-function fe_problem(u, n)
-
-  f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x)
-
-  domain = (0, 1, 0, 1)
-  partition = (n, n)
-  model = CartesianDiscreteModel(domain, partition)
-
-  order = 1
-
-  V0 = FESpace(
-    reffe = :Lagrangian,
-    order = order,
-    valuetype = Float64,
-    conformity = :H1,
-    model = model,
-    dirichlet_tags = "boundary",
-  )
-  U = TransientTrialFESpace(V0, u)
-
-  trian = Triangulation(model)
-  degree = 2 * order
-  quad = CellQuadrature(trian, degree)
-
-  a(u, v) = ∇(v) ⋅ ∇(u)
-  b(v, t) = v * f(t)
-
-  res(t, u, ut, v) = a(u, v) + ut * v - b(v, t)
-  jac(t, u, ut, du, v) = a(du, v)
-  jac_t(t, u, ut, dut, v) = dut * v
-
-  t_Ω = FETerm(res, jac, jac_t, trian, quad)
-  op = TransientFEOperator(U, V0, t_Ω)
-
-  U0 = U(0.0)
-  uh0 = interpolate_everywhere(U0, u(0.0))
-
-  return op, trian, quad, uh0
-
-end
-
-# ODE solver in GridapODEs
-function solve_ode_gridap(op, trian, quad, uh0, θ, dt, t0, tF, u)
-  # GridapODEs version using Backward Euler
-  ls = LUSolver()
-  nls = NLSolver(ls; show_trace = true, method = :newton)
-  odes = ThetaMethod(ls, dt, θ)
-  solver = TransientFESolver(odes)
-  sol_t = Gridap.solve(solver, op, uh0, t0, tF)
-
-  l2(w) = w * w
-
-  tol = 1.0e-6
-  _t_n = t0
-
-  for (uh_tn, tn) in sol_t
-    _t_n += dt
-    e = u(tn) - uh_tn
-    el2 = sqrt(sum(integrate(l2(e), trian, quad)))
-    @test el2 < tol
-    @show get_free_values(uh_tn)
-  end
-
-end
-
-# Wrappers for residual and jacobian for DiffEqs
-function diffeq_wrappers(op)
-
-  ode_op = get_algebraic_operator(op)
-  ode_cache = allocate_cache(ode_op)
-
-  function residual(res, du, u, p, t)
-    # TO DO (minor): Improve update_cache! st do nothing if same time t as in the cache
-    # now it would be done twice (residual and jacobian)
-    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
-    GridapODEs.ODETools.residual!(res, ode_op, t, u, du, ode_cache)
-  end
-
-  function jacobian(jac, du, u, p, gamma, t)
-    ode_cache = GridapODEs.ODETools.update_cache!(ode_cache, ode_op, t)
-    z = zero(eltype(jac))
-    Gridap.Algebra.fill_entries!(jac, z)
-    GridapODEs.ODETools.jacobian_and_jacobian_t!(jac, ode_op, t, u, du, gamma, ode_cache)
-  end
-
-  return residual, jacobian
-
-end
-
-end #module
-
-
-module DiffEqsWrapperTests
-
-using Gridap
-using ForwardDiff
-using LinearAlgebra
-using Test
-using GridapODEs.ODETools
-using GridapODEs.TransientFETools
-using Gridap.FESpaces: get_algebraic_operator
-using GridapODEs
-using DifferentialEquations
-using Gridap.Algebra: NewtonRaphsonSolver
-using Base.Iterators
-
-# Solving the heat equation using GridapODEs and DiffEqs
-
-u(x, t) = t#1.0
-u(t) = x -> u(x, t)
-
-# FE structs
-op, trian, quad, uh0 = DiffEqsWrappersAux.fe_problem(u,3)
-
-θ = 1.0
-t0 = 0.0
-tF = 1.0
-dt = 0.1
-
-# Check Gridap ODE solver
-DiffEqsWrappersAux.solve_ode_gridap(op,trian,quad,uh0,θ,dt,t0,tF,u)
-
-# DiffEq wrappers
-residual, jacobian = DiffEqsWrappersAux.diffeq_wrappers(op)
-
-# Check wrappers
-u0 = get_free_values(uh0)
-dtθ = dt * θ
-ode_op = get_algebraic_operator(op)
-ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
-r = GridapODEs.ODETools.allocate_residual(ode_op, u0, ode_cache)
-J = GridapODEs.ODETools.allocate_jacobian(ode_op, u0, ode_cache)
-
-tθ = 1.0
-residual(r, u0, u0, nothing, tθ)
-jacobian(J, u0, u0, nothing, (1 / dtθ), tθ)
-
-J
-
-# Using Sundials to solve the resulting ODE
-using Sundials
-tspan = (0.0, 1.0)
-
-u(x, t) = t
-u(t) = x -> u(x, t)
-
-# ISSUE 1: When I choose n > 2, even though the problem that we will solve is
-# linear, the Sundials solvers seems to have convergence issues in the nonlinear
-# solver (?). Ut returns errors
-# [IDAS ERROR]  IDACalcIC Newton/Linesearch algorithm failed to converge.
-# ISSUE 2: When I pass `jac_prototype` the code gets stuck
-n = 2
-op, trian, quad, uh0 = DiffEqsWrappersAux.fe_problem(u,n)
-residual, jacobian = DiffEqsWrappersAux.diffeq_wrappers(op)
-u0 = get_free_values(uh0)
-
-# To explore the Sundials solver options, e.g., BE with fixed time step dtd
-f_iip = DAEFunction{true}(residual; jac = jacobian) #,jac_prototype=J)
-# jac_prototype is the way to pass my pre-allocated jacobian matrix
-prob_iip = DAEProblem{true}(f_iip, u0, u0, tspan, differential_vars = [true])
-sol_iip = Sundials.solve(prob_iip, IDA(linear_solver=:KLU), reltol = 1e-8, abstol = 1e-8)
-@show sol_iip.u
-
-# or iterator version
-integ = init(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
-# step!(integ)
-
-# Show using integrators as iterators
-for i in take(integ, 100)
-  @show integ.u
-end
-
-end # module
-
-# ISSUE: Future work, add own (non)linear solver.
-# Let us assume that we just want to consider a fixed point algorithm
-# and we consider an implicit time integration of a nonlinear PDE.
-# Our solvers are efficient since they re-use cache among
-# nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization)
-# and for the transient case in our library, we can also reuse all this
-# between time steps. Could we attain something like this using
-# DifferentialEquations/Sundials?
-# @ChrisRackauckas suggests to take a look at:
-# https://docs.sciml.ai/latest/tutorials/advanced_ode_example/ shows swapping out linear solvers.
-# https://docs.sciml.ai/latest/features/linear_nonlinear/ is all of the extra details.
-
-# Try to pass solver too
-# ls = LUSolver()
-# ls_cache = nothing
-# x = copy(u0)
-# solve!(x,J,r,ls_cache)
-#
diff --git a/test/runtests.jl b/test/runtests.jl
index 3a77daf..edc4e51 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,11 +1,13 @@
 module GridapODEsRunTests
 
 using Test
+@time @testset "DiffEqsWrappers" begin include("DiffEqsWrappersTests/runtests.jl") end
 
 @time @testset "ODETools" begin include("ODEsTests/runtests.jl") end
 
 @time @testset "TransientFETools" begin include("TransientFEsTests/runtests.jl") end
 
+
 # include("../bench/runbenchs.jl")
 
 end #module