Create a linear version of the types

[santiagobadia]

Up to now, I am considering general nonlinear solvers and operators.

We should consider also linear ones for performance.

For efficiency, we should explore the case in which the TrialFESpace is constant, not being tested yet.

[santiagobadia]

Linear Transient FE Operators

I am thinking about the definition of linear transient FE operators in GridapODEs.

For linear problems, there is not much to do. I think that the ThetaMethod can be kept as it is, since a LinearSolver is a sub-type of a NonlinearSolver. Due to this fact, for linear problems, we can already use linear solvers, and there is no performance hit for linear problems now.

To construct operators using the a-l form (bilinear and linear forms), I would create a constructor of the current TransientFETerm that takes m(t,ut,v), a(t,u,v) and l(t,v) and generates res and jac. That is easy, we have that res(t,u,ut,v) = m(t,u,v) + a(t,u,v)-l(t,v) and jac(t,du,v) = alpha*m(t,du,v)+a(t,du,v). Still to understand what does it imply in terms of implementation. In any case, I would leave this for a later stage.

Constant Transient FE Operators

Here the problem is that we are creating a new operator each time step, and now we want to preserve it. Probably, the idea is to keep the operator in the ode_cache, check that the time step size has not changed, and if it doesn't, just take the operator from the previous iteration. I am not sure we can even check this.

Additionally, we can consider a constant version of the operator that can be constructed from a(u,v), m(u,v) and l(t,v) (or l(v)); a time-dependent RHS does not impact performance. Now, we can construct our transient operator with res(t,u,ut,v) = m(u,t,v) + a(u,v)-l(t,v). Again, that is a aesthetical issue that can be handled later.

I would try hardly not to create new constant versions of ThetaMethod, ThetaMethodNonlinearOperator, TransientFEOperator, TransientFETerms ... I want to think how to provide this functionality with minor modifications, for the sake of readability of the code.

[fverdugo]

I would try to make the transient code as analogous as possible to the steady state.

Disclaimer: the names below can be improved.

For instance, the terms in the heat equations with Neumann bcs:

m(t,ut,v) = ut*v
a(t,u,v) = ∇(v)⊙∇(u)
l_Ω(t,v) = v*f(t)
l_Γ(t,v) = v*g(t)
t_Ω = TransientAffineFETerm(m,a,l_Ω,trian,quad)
t_Γ = TransientFESource(l_Γ,trian_Γ,quad_Γ)
hasconstcoeffs=true
op = TransientAffineFEOperator{hasconstcoeffs}(U,V,t_Ω,t_Γ)

Due to this fact, for linear problems, we can already use linear solvers, and there is no performance hit for linear problems now.

This is true, but I think it would compute more things than strictly needed, you cannot take advantage of elemental matrices to impose Dirichlet boundary conditions, etc.

To avoid this, I would introduce the same ingredients as in the steady state case:

# C flags if it is constant or not
abstract type AffineODEOperator{C} <: ODEOperator end

funciton matrix!(K,op::AffineODEOperator,t::Real)
  @abstractmethod
end

funciton matrix_t!(M,op::AffineODEOperator,t::Real)
  @abstractmethod
end

funciton vector!(b,op::AffineODEOperator,t::Real)
  @abstractmethod
end

abstract type AffineTransientFEOperator{C} <:TransientFEOperator end

funciton matrix!(K,feop::AffineTransientFEOperator,t::Real)
  @abstractmethod
end

funciton matrix_t!(M,feop::AffineTransientFEOperator,t::Real)
  @abstractmethod
end

funciton vector!(b,feop::AffineTransientFEOperator,t::Real)
  @abstractmethod
end

function get_algebraic_operator(feop::AffineTransientFEOperator{C}) where C
  AffineODEOpFromFEOp{C}(feop)
end

struct AffineODEOpFromFEOp{C} <: AffineODEOperator{C} end

# The actual implementation is here
struct AffineTransientFEOperatorFromTerms <: AffineTransientFEOperator{false} end
struct ConstantAffineTransientFEOperatorFromTerms <: AffineTransientFEOperator{true} end

Then the ODESolvers like the ThetaMethod can dispatch on the ODEOperator and implement
efficient particular cases in terms of the interface of AffineODEOperator{true/false} and the LinearSolver interface.

This introduces new types, but It opens the door to implement all the optimizations one would expect for linear problems.

Perhaps it is a good idea to discuss it face-to-face.

[santiagobadia]

Before discussing the FE part, and whether or not we want to start creating transient FE terms for all possible combinations (now we have three kinds of terms, not two), I think it is better to have a clear idea of how to deal with constant operators in ODETools. Note that the comment about the Dirichlet data for linear elements is far more complicated here, since you also have the ones that come from Dirichlet data, times a factor that depends on the ODESolver. I think it is quite a mess. We can create a wrapper as the one you propose but always use the jac, jac_t, res implementation, at least for the affine operator. For the constant operator, since everything is constant, we could consider matrix/vector-type implementation.

Being said that, I think that the main issue is how to deal with constant operators (and constant right-hand side and time step size due to Dirichlet data, otherwise I am not sure how we can do this, since we do not keep the Free-Dirichlet block of the matrix in Gridap). What do we need to accommodate this situation in the ODETools part?

I think we must create a ConstantODEOperator and its corresponding ConstantODEOpFromFEOp with a new update_cache!, residual!, jacobian! and jacobian_t!, a new ThetaMethodOperator constructor and a new ConstantThetaMethodOperator with a re-definition of residual! and jacobian! (or matrix/vector instead of jacobian/residual).

When it will work, we can think about a TransientConstantFEOperator with vector!, matrix!, matrix_t! to be used in ConstantODEOpFromFEOp. Again, I am not sure we can do this for affine operators that are not constant (but we can always create a constructor like a linear method), but keeping the current underlying implementation.

Anyway, let us have a chat on this to speed up the process.

[santiagobadia]

I propose:

1. Create a trait for ODEOperator

"""
Trait for ODEOperator that tells us whether the operator depends on the solution
(including its time derivatives), it is an affine operator that depends on time
or it is a constant operator (affine and time-indepedendent)
"""
abstract type ODEDependence end
struct Nonlinear end
struct Affine end
struct Constant end

const AffineODEOperator = ODEOperator{Affine}
const ConstantODEOperator = ODEOperator{Constant}

2. Use the same trait for TransientFEOperator

"""
A transient version of the `Gridap` `FEOperator` that depends on time. The
parameter `C` determines whether it is a nonlinear, affine or constant operator
"""
abstract type TransientFEOperator{T} <: GridapType end

and TransientFEOperatorFromTerms, and define the corresponding constructor using kwargs

function TransientFEOperator(trial,test,terms...;type="nonlinear")
  assem_t = SparseMatrixAssembler(test,evaluate(trial,nothing))
  if (type == "nonlinear")
    T = Nonlinear
  elseif (type == "affine")
    T = Affine
  elseif (type == "constant")
    T = Constant
  else
    error("Operator type not defined, it must be nonlinear, affine or constant")
  end
  TransientFEOperatorFromTerms{T}(trial,∂t(trial),test,assem_t,terms...)
end

3. Use this trait in the ODEOpFromFEOp by modifying the get_algebraic_operator as follows

"""
Returns a `ODEOperator` wrapper of the `TransientFEOperator` that can be
straightforwardly used with the `ODETools` module.
"""
function get_algebraic_operator(feop::TransientFEOperator{T}) where T
  ODEOpFromFEOp{T}(feop)
end

4. Using this trait for dispatching in solver_step! method of ThetaMethod

function solve_step!(uf::AbstractVector,
                     solver::ThetaMethod,
                     op::AffineODEOperator,
                     u0::AbstractVector,
                     t0::Real,
                     cache) 
# here the body

However, it won't work because ODEOpFromFEEOp{Constant} <: ODEOperator{Constant}

Probably doing this

struct ODEOpFromFEOp{C} <: ODEOperator{C} where C <: ODEDepedence
  feop::TransientFEOperator
end

const AffineODEOpFromFEOp = AffineODEOpFromFEOp{Affine}
const ConstantODEOpFromFEOp = ODEOpFromFEOp{Constant}

[fverdugo]

OK, we can start by considering the affine non-constant scenario as a particular non-linear case.

We can define this constructor

function Gridap.AffineFETerm(
  m::Function,a::Function,b::Function,trian::Triangulation,quad::CellQuadrature)
  res(t,u,ut,v) = m(t,ut,v) + a(t,u,v) - b(t,v)
  jac(t,u,ut,du,v) = a(t,du,v)
  jac_t(t,u,ut,dut,v) = m(t,dut,v)
  FETerm(res,jac,jac_t,trian,quad)
end

In order to write

m(t,ut,v) = ut*v
a(t,u,v) = ∇(v)⋅∇(u)
b(t,v) = v*f(t)
t_Ω = AffineFETerm(m,a,b,trian,quad)
op = TransientFEOperator(U,V0,t_Ω)

instead of

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_Ω)

[fverdugo]

On the other hand, I will create a new abstract type ConstantODEOperator instead of introducing a trait. The idea is to extend the library with this new functionality without touching what is already working.

This is the interface for this abstract type:

abstract type ConstantODEOperator <: ODEOperator

get_matrix(op::ConstantODEOperator) = @abstractmethod
get_matrix_t(op::ConstantODEOperator) = @abstractmethod
vector!(b,op::ConstantODEOperator,t::Real,ode_cache) = @abstractmethod
# + Other methods in the ODEOperator interface

Then, we can implement an efficient version of solve_step! in the theta method for this particular case

function solve_step!(uf::AbstractVector,
                     solver::ThetaMethod,
                     op::ConstantODEOperator,
                     u0::AbstractVector,
                     t0::Real,
                     cache)

  dt = solver.dt
  solver.θ == 0.0 ? dtθ = dt : dtθ = dt*solver.θ
  tθ = t0+dtθ

  if cache === nothing
    ode_cache = allocate_cache(op)
    vθ = similar(u0)
    bθ = similar(u0)
    ls_cache = nothing
    K = get_matrix(op)
    M = get_matrix_t(op)
    A = (1/dtθ)*M + K
    ls_cache = nothing
  else
    ode_cache, vθ, bθ, A, ls_cache = cache
  end

  ode_cache = update_cache!(ode_cache,op,tθ)
  vector!(bθ,op,tθ,ode_cache)

  aop = AffineOperator(A,bθ)
  newmatrix = false
  ls_cache = solve!(uf,solver.nls,aop,ls_cache,newmatrix)

  if 0.0 < solver.θ < 1.0
    uf = uf*(1.0/solver.θ)-u0*((1-solver.θ)/solver.θ)
  end

  cache = (ode_cache, vθ, bθ, A, ls_cache)

  tf = t0+dt
  return (uf,tf,cache)

end

[fverdugo]

The last ingredient is a specialization of ConstantODEOperator from all FE-related machinery. In general, this will require types ConstantTransientFEOperator (better name?) and ConstantODEOpFromFEOp

Here, we could even precompute the dirichlet blocks of M and K using the DirichletFESpace so that we only integrate the source terms at each time step

[santiagobadia]

Hi @fverdugo

I have done an implementation of affine and constant operators. I have finally decided to use traits, I find it is a nice way to provide the required functionality with a minimum number of lines.

I have to test the code, we can talk at some point, but I wanted to ask you about this error

GridapODEs.jl/test/TransientFEsTests/TransientFETests.jl

Line 93 in 81522e5

# santiagobadia : @fverdugo I had this problem before related to terms...

I always have issues with the constructor for TransientFEOperatorFromTerms, due to the terms... last argument, I guess. When I have only one term, we cannot iterate on it, and it returns an error. I did a hack and I remember you fixed it, but I would like to understand the problem.

[fverdugo]

Try to remove the ... after terms in this line:

GridapODEs.jl/src/TransientFETools/TransientFEOperators.jl

Line 90 in 81522e5

new(trial,trial_t,test,assem,optype,terms...)

[santiagobadia]

Even though there is still room for optimisation, I think that the implementation after PR #18 solve this issue