Refactoring of GridapODEs

[AlexandreMagueresse]

As the ODE module of Gridap is expanding, its structure is becoming a little unclear and redundant. I propose the following major changes to the API. The general framework of the module is described here. I made other minor changes internally.

ODE operators

• Added a type hierarchy for ODE operators: NonlinearODE, QuasilinearODE, SemilinearODE and LinearODE
• An ODE operator is now represented in terms of its linear forms (if any) and (sub)residual, so that the full residual is split in the following way, where $N$ be the order of the ODE.

1. NonlinearODE: $\mathrm{residual}(t, \partial_{t}^{0} u, \ldots, \partial_{t}^{N} u) = \mathrm{res}(t, \partial_{t}^{0} u, \ldots, \partial_{t}^{N} u)$
2. QuasilinearODE: $\mathrm{residual}(t, \partial_{t}^{0} u, \ldots, \partial_{t}^{N} u) = \mathrm{mass}(t, \partial_{t}^{0} u, \ldots, \partial_{t}^{N-1} u) \partial_{t}^{N} u + \mathrm{res}(t, \partial_{t}^{0} u, \ldots, \partial_{t}^{N-1} u)$
3. SemilinearODE: $\mathrm{residual}(t, \partial_{t}^{0} u, \ldots, \partial_{t}^{N} u) = \mathrm{mass}(t) \partial_{t}^{N} u + \mathrm{res}(t, \partial_{t}^{0} u, \ldots, \partial_{t}^{N-1} u)$
4. LinearODE: $\mathrm{residual}(t, \partial_{t}^{0} u, \ldots, \partial_{t}^{N} u) = \sum_{0 \leq k \leq N} \mathrm{form}_{k}(t) \partial_{t}^{k} u + \mathrm{res}(t)$

• Implicit-explicit operators are represented as a sum of two operators: an implicit operator of order $N$ and an explicit operator of order $N-1$.

Transient FE spaces, cell fields and operators

• Only minor changes were brought to transient FE spaces: TransientTrialFESpace and TransientMultiFieldFESpace are now subtypes of AbstractTransientTrialFESpace, which is itself a subtype of `FESpace.
• The time derivative operator has been extended to arbitrary order. In particular, transient FE spaces and cell-fields now support arbitrary-order time derivatives.
• The previous constructors for transient FE operators have been replaced by the following ones:

1. TransientFEOperator(res, jacs, trial, test) and TransientFEOperator(res, trial, test; order) for the version with AD. The residual is expected to have the signature residual(t, u, v).
2. TransientQuasilinearFEOperator(mass, res, jacs, trial, test) and TransientQuasilinearFEOperator(mass, res, trial, test; order) for the version with AD. The mass and residual are expected to have the signatures mass(t, u, dtu, v) and residual(t, u, v), i.e. mass is written as a linear form of the highest-order time derivative dtu. Since in this setting the mass matrix is supposed to depend on lower-order time derivatives, u is provided for the nonlinearity of the mass matrix.
3. TransientSemilinearFEOperator(mass, res, jacs, trial, test; constant_mass) and TransientSemilinearFEOperator(mass, res, trial, test; order, constant_mass) for the version with AD (no AD is needed for the mass term). The mass and residual are expected to have the signatures mass(t, dtu, v) and residual(t, u, v), where here again dtu is the highest-order derivative, i.e. mass is specified as a linear form of du.
4. TransientLinearFEOperator(forms, res, jacs, trial, test; constant_forms) and TransientLinearFEOperator(forms, res, trial, test; constant_forms) for the version with AD (in fact no AD is needed at all). The forms and residual are expected to have the signatures form_k(t, dtku, v) and residual(t, v), i.e. form_k is a linear form of the $k$-th order derivative, and the residual does not depend on u.

• I am not sure that my indications of the signatures is clear above, so just to fix ideas, considering the heat equation:

• As a TransientFEOperator: res(t, u, v) = ∫(∂t(u) ⋅ v + ∇(u) ⋅ ∇(v) - f(t) ⋅ v) * dΩ
• As a TransientQuasilinearFEOperator: mass(t, u, dtu, v) = ∫(dtu ⋅ v) * dΩ and res(t, u, v) = ∫(∇(u) ⋅ ∇(v) - f(t) ⋅ v) * dΩ
• As a TransientQuasilinearFEOperator: mass(t, dtu, v) = ∫(dtu ⋅ v) * dΩ and res(t, u, v) = ∫(∇(u) ⋅ ∇(v) - f(t) ⋅ v) * dΩ
• As a TransientLinearFEOperator: mass(t, dtu, v) = ∫(dtu ⋅ v) ) * dΩ, stiffness(t, u, v) = ∫(∇(u) ⋅ ∇(v)) * dΩ and res(t, u, v) = ∫(- f(t) ⋅ v) * dΩ. Important: note the minus sign here. We could change that if it is more intuitive to have res equal to the forcing term rather than its opposite.

• As shown above, it is possible to indicate that linear forms are constant. If this is the case, the will be assembled once, store and only retrieved when calling jacobian!.
• In order to limit the number of conversions between vector of degrees of freedom and FE function, especially when we don't have to convert when the jacobians are stored, the computation of residuals and jacobians is no longer performed on the transient FE operator, but rather on the ODEOpFromFEOp, which is a subtype of ODEOperator. In that sense TransientFEOperator does not fully implement the interface of an FEOperator, it is rather at the intersection of FEOperator and ODEOperator.
• To mimic the ODESolution interface, FESolution is now an abstract type and the equivalent of GenericODESolution is GenericFESolution.
• The order of the arguments to create an FE solution has changed: instead of solve(solver, operator, u0, t0, tF), it is now solve(solver, operator, t0, tF, u0). Similarly, when iterating over an FE solution, the output is (t_n, u_n) rather than (u_n, t_n). These changes follow from using the order (space, time, functional space) in the whole module.

ODE solvers

• When applying an ODE solver to an ODE operator, a DiscreteODEOperator is created (it is a subtype of NonlinearOperator). It represents the discretisation of the time derivative performed by the ODE solver. However, whenever possible, an ODE solver should instantiate a LinearDiscreteODEOperator instead, which behaves like an AffineOperator. This is always possible when an explicit ODE solver is applied to a quasilinear ODE, or when a (diagonally-) implicit ODE solver is applied to a linear ODE.
• By default, the jacobian matrix of the DiscreteODEOperator needs to be updated across consecutive time steps, even when the DiscreteODEOperator is actually a LinearDiscreteODEOperator and the solver for this system of equation is a LinearSolver. In particular, the solve! function has been overwritten to recompute the jacobian matrix and its factorisation. However, under special circumstances, the jacobian matrix may be constant, and it is possible to indicate so at the LinearDiscreteODEOperator level.
• All the previous ODE solvers have been rewritten to solve for the highest-order derivative of the ODE. For example with the $\theta$ scheme, the previous version performed the following update $$\textrm{Find } x \textrm{ such that } \mathrm{residual}(t_{\theta}, x, \frac{1}{\theta h} (x - u_{n})) = 0, \qquad u_{n+1} = u_{n} + \frac{1}{\theta} (x - u_{n}),$$ where $h = t_{n+1} - t_{n}$ is the time step. The new implementation of the $\theta$ scheme now solves for $y = \frac{1}{\theta h} (x - u_{n})$ and performs the equivalent update $$\textrm{Find } x \textrm{ such that } \mathrm{residual}(t_{\theta}, u_{n} + \theta h x, x) = 0, \qquad u_{n+1} = u_{n} + h x.$$

Notes

Additionally to these changes to the ODE module, I made the following changes.

• I added an axpy_entries! function to Algebra that enables adding two matrices, with an efficient specialisation for sparse matrices.
• Together with @JordiManyer, we fixed the behaviour of automatic differentiation on MultiFieldCellField. The issue was that the return type of the jacobian provided by ForwardDiff was BlockMatrix{<:SubArray}, and its corresponding testvalue was wrong (it was Array{T, 0} instead of Array{T, P}, where P is the dimension of the quantity being differentiated, i.e. 2 for jacobians).
• MultiFieldCellField and MultiFieldFEFunction now implement Base.length so that they can be zipped and mapped over.

To discuss

I would like to discuss the following changes that directly impact the user.

• When iterating over an FESolution, the iterator returns the time first and then the FEFunction. It used to return the FEFunction first and then the time.
• When describing a TransientLinearFEOperator, the order is stiffness, then mass, then residual, where residual is the opposite of the forcing term. This has the advantage of keeping the same order as the other constructors (nonlinear, semilinear and quasilinear) but is quite unnatural from a user perspective. Another question is the order in which the jacobians are given, and whether the linear forms are constant. I am not happy about this but I have not been able to come up with a satisfying answer for now.

I have some final minor questions.

• Is const GridapODEs = ODEs useful? See here.
• In keeping with the previous code, (space::FESpace)(t) = evaluate(space, t) is only defined for julia versions above 1.3. Do we need to keep compatibility with julia versions lower than 1.3? Also, why is this not being defined for earlier versions? See here.

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 80.12976% with 490 lines in your changes are missing coverage. Please review.

Project coverage is 87.59%. Comparing base (faacb16) to head (f578f90). Report is 2 commits behind head on master.

Files	Patch %	Lines
src/ODEs/ODESolvers/RungeKuttaEX.jl	0.00%	84 Missing :warning:
src/ODEs/TransientFEOperators.jl	71.25%	71 Missing :warning:
src/ODEs/ODESolvers/ForwardEuler.jl	0.00%	63 Missing :warning:
src/ODEs/TransientCellFields.jl	50.96%	51 Missing :warning:
src/ODEs/ODEOperators.jl	75.88%	41 Missing :warning:
src/ODEs/ODESolvers/Tableaus.jl	61.76%	26 Missing :warning:
src/ODEs/ODESolvers/RungeKuttaIMEX.jl	88.50%	20 Missing :warning:
src/ODEs/ODESolvers/TableausEX.jl	88.60%	18 Missing :warning:
src/ODEs/TransientFESpaces.jl	83.83%	16 Missing :warning:
src/ODEs/TransientFESolutions.jl	71.69%	15 Missing :warning:
... and 12 more

:exclamation: Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #965      +/-   ##
==========================================
- Coverage   88.12%   87.59%   -0.54%     
==========================================
  Files         177      176       -1     
  Lines       21122    21944     +822     
==========================================
+ Hits        18613    19221     +608     
- Misses       2509     2723     +214     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[AlexandreMagueresse]

Hi @santiagobadia @amartinhuertas @JordiManyer @oriolcg @tamaratambyah, I am finally done with all the changes I had in mind. I think the pull request is ready for review. The only thing I would like to add is a markdown file docs/ODEs.md to describe the framework I followed and briefly comment on the numerical methods. I will add this file later today.

[oriolcg]

Hi @AlexandreMagueresse, I'm encountering the following error when using the TransientSemilinearFEOperator:

ERROR: LoadError: AssertionError: It is only possible to efficiently add two sparse matrices that have the same
sparsity pattern.

Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/Gridap/EALIA/src/Helpers/Macros.jl:60 [inlined]
  [2] axpy_entries!(α::Int64, A::SparseArrays.SparseMatrixCSC{Float64, Int64}, B::SparseArrays.SparseMatrixCSC{Float64, Int64}; check::Bool)
    @ Gridap.Algebra ~/.julia/packages/Gridap/EALIA/src/Algebra/AlgebraInterfaces.jl:277
  [3] #62
    @ ~/.julia/packages/GridapDistributed/xLNTs/src/Algebra.jl:57 [inlined]
  [4] map(::GridapDistributed.var"#62#63"{Bool, Int64}, ::PartitionedArrays.MPIArray{SparseArrays.SparseMatrixCSC{Float64, Int64}, 1}, ::PartitionedArrays.MPIArray{SparseArrays.SparseMatrixCSC{Float64, Int64}, 1})
    @ PartitionedArrays ~/.julia/packages/PartitionedArrays/D9r9E/src/mpi_array.jl:229
  [5] #axpy_entries!#61
    @ ~/.julia/packages/GridapDistributed/xLNTs/src/Algebra.jl:56 [inlined]
  [6] axpy_entries!(α::Int64, A::PartitionedArrays.PSparseMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, PartitionedArrays.MPIArray{SparseArrays.SparseMatrixCSC{Float64, Int64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.SparseMatrixAssemblyCache, 1}, Float64}, B::PartitionedArrays.PSparseMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, PartitionedArrays.MPIArray{SparseArrays.SparseMatrixCSC{Float64, Int64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.SparseMatrixAssemblyCache, 1}, Float64})
    @ GridapDistributed ~/.julia/packages/GridapDistributed/xLNTs/src/Algebra.jl:49
  [7] jacobian_add!(J::PartitionedArrays.PSparseMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, PartitionedArrays.MPIArray{SparseArrays.SparseMatrixCSC{Float64, Int64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.SparseMatrixAssemblyCache, 1}, Float64}, odeop::Gridap.ODEs.ODEOpFromTFEOp{Gridap.ODEs.SemilinearODE}, t::Float64, us::Tuple{PartitionedArrays.PVector{Vector{Float64}, PartitionedArrays.MPIArray{Vector{Float64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.LocalIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.VectorAssemblyCache{Float64}, 1}, Float64}, PartitionedArrays.PVector{Vector{Float64}, PartitionedArrays.MPIArray{Vector{Float64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.LocalIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.VectorAssemblyCache{Float64}, 1}, Float64}}, ws::Tuple{Float64, Int64}, odeopcache::Gridap.ODEs.ODEOpFromTFEOpCache)
    @ Gridap.ODEs ~/.julia/packages/Gridap/EALIA/src/ODEs/ODEOpsFromTFEOps.jl:347
  [8] jacobian!(J::PartitionedArrays.PSparseMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, PartitionedArrays.MPIArray{SparseArrays.SparseMatrixCSC{Float64, Int64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.SparseMatrixAssemblyCache, 1}, Float64}, odeop::Gridap.ODEs.ODEOpFromTFEOp{Gridap.ODEs.SemilinearODE}, t::Float64, us::Tuple{PartitionedArrays.PVector{Vector{Float64}, PartitionedArrays.MPIArray{Vector{Float64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.LocalIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.VectorAssemblyCache{Float64}, 1}, Float64}, PartitionedArrays.PVector{Vector{Float64}, PartitionedArrays.MPIArray{Vector{Float64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.LocalIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.VectorAssemblyCache{Float64}, 1}, Float64}}, ws::Tuple{Float64, Int64}, odeopcache::Gridap.ODEs.ODEOpFromTFEOpCache)
    @ Gridap.ODEs ~/.julia/packages/Gridap/EALIA/src/ODEs/ODEOperators.jl:353
  [9] jacobian!(J::PartitionedArrays.PSparseMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, PartitionedArrays.MPIArray{SparseArrays.SparseMatrixCSC{Float64, Int64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.OwnAndGhostIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.SparseMatrixAssemblyCache, 1}, Float64}, nlop::Gridap.ODEs.NonlinearStageOperator, x::PartitionedArrays.PVector{Vector{Float64}, PartitionedArrays.MPIArray{Vector{Float64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.LocalIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.VectorAssemblyCache{Float64}, 1}, Float64})
...

I'm solving the transient Navier-Stokes equations where, indeed, the mass and stiffness matrices don't have the same sparsity pattern. Have you thought about this issue? Is there a way to address this? Or should I just use another operator?

[AlexandreMagueresse]

Hi @oriolcg, thank you for finding that issue. My last commit should fix it.

[oriolcg]

Hi @AlexandreMagueresse, now I'm getting another error when using AD with multifield. I'm not sure if this is an issue only on parallel or if in serial also appears...

ERROR: LoadError: type DistributedMultiFieldCellField has no field cellfield
Stacktrace:
  [1] getproperty(x::GridapDistributed.DistributedMultiFieldCellField{Vector{GridapDistributed.DistributedCellField{PartitionedArrays.MPIArray{Gridap.ODEs.TransientSingleFieldCellField{Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{Gridap.CellData.ReferenceDomain}}}, 1}, GridapDistributed.DistributedTriangulation{2, 2, ...
    @ Gridap.CellData ~/.julia/packages/Gridap/iy2tm/src/CellData/CellFields.jl:674
  [2] (::Gridap.ODEs.var"#jac_0#90"{PerforatedCylinder.var"#res#49"{PerforatedCylinder.var"#c#47", PerforatedCylinder.var"#τc#46"{PerforatedCylinder.var"#τₘ#45"{GridapDistributed.DistributedCellField{PartitionedArrays.MPIArray{Gridap.CellData.GenericCellField{Gridap.CellData.ReferenceDomain}, 1}, GridapDistributed.DistributedTriangulation{2, 2, ...
    @ Gridap.ODEs ~/.julia/packages/Gridap/iy2tm/src/ODEs/TransientFEOperators.jl:558
  [3] allocate_odeopcache(odeop::Gridap.ODEs.ODEOpFromTFEOp{Gridap.ODEs.SemilinearODE}, t::Float64, us::Tuple{PartitionedArrays.PVector{Vector{Float64}, PartitionedArrays.MPIArray{Vector{Float64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.LocalIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.VectorAssemblyCache{Float64}, 1}, Float64}, ...
    @ Gridap.ODEs ~/.julia/packages/Gridap/iy2tm/src/ODEs/ODEOpsFromTFEOps.jl:109
  [4] allocate_odecache(odeslvr::Gridap.ODEs.ThetaMethod, odeop::Gridap.ODEs.ODEOpFromTFEOp{Gridap.ODEs.SemilinearODE}, t0::Float64, us0::Tuple{PartitionedArrays.PVector{Vector{Float64}, PartitionedArrays.MPIArray{Vector{Float64}, 1}, PartitionedArrays.MPIArray{PartitionedArrays.LocalIndices, 1}, PartitionedArrays.MPIArray{PartitionedArrays.VectorAssemblyCache{Float64}, 1}, Float64}})
    @ Gridap.ODEs ~/.julia/packages/Gridap/iy2tm/src/ODEs/ODESolvers/ThetaMethod.jl:50
...

It's failing in the allocation of the jacobians when using a SemilinearOperator, I think it is expecting a TransientCellField and it's receiving just a DistributedCellField. Any clue on this?

[JordiManyer]

@oriolcg AD is not supported by GridapDistributed yet (and never was). We have a working prototype, but it's not part of the library.

Although the error you point out is indeed caused by accessing properties of a TransientCellField, the issue is much more involved (and has to do with how measures are currently embedded in the residual functions).

So long story short, this is an issue in distributed only, and it has nothing to do with this PR.