General Runge Kutta Solver

Registered by Johan Hake on 2013-02-20

A general time integration framework for any time dependent PDE/ODE will be added. A user provides a form representing the RHS together with the prognostic variable. If the form is time dependent a Constant for the time variable will also be needed. The form should be of rank 1 and potentially nonlinear in the prognostic variable.

The time integration methods will at first be limited to lower triangular Butcher Tableau based methods. This means all explicit and diagonally implicit RK methods.

In Python a ButcherTablau class takes a UFL form and generates a set of forms which will be used to for the different stages.
* If a stage is Explicit it will only require a single assemble.
* If the stage is implicit a (potential) nonlinear solve is required.
* The last stage is always a linear combination of the stage answers. To handle the final stage an axpy assign will be added to Function in DOLFIN.

At first this will be implemented as a Python only feature. However, a light weight layer could be added to ufl and for generation of DOLFIN wrappers to make it available in C++ too.

A special case for this solver is when the form only consists of PointIntegrals with CG1 (Vector) test functions. Then the dofs are spatially decoupled and a local solve can be preformed for each vertex.

Blueprint information

Status:
Started
Approver:
Johan Hake
Priority:
Essential
Drafter:
Johan Hake
Direction:
Needs approval
Assignee:
Johan Hake
Definition:
Approved
Series goal:
None
Implementation:
Good progress
Milestone target:
milestone icon 1.1.1
Started by
Johan Hake on 2013-02-20

Sprints

Whiteboard

GNW: This just arrived in Boost: Don't know if it's useful.

http://www.boost.org/doc/libs/1_53_0/libs/numeric/odeint/doc/html/index.html

JH: I have looked at it and it looks pretty neat, but it is only AFAIK for ODEs. Even if my initial need limited itself to decoupled local ODEs, I needed the solver to be adjointable. Unfortunately are not RK solver generally self adjointable, which prompted us (Patrick, me and Marie) to implement this in the UFL layer. But once that is done we get a time integrator for any PDE virtually for free.

SP (Simone Pezzuto): Wouldn't it be nice to implement also an interface to TS framework of PETSc? They provide a consistent way to solve time dependent problems, in parallel, with several features. It's also interesting because it might be a good alternative (not replacement of course) for SNES for continuation problems. Also Trilinos has different capabilities in this respect (Loca, Rythmos, ...), but the interface to Epetra is not always consistent (Vector, MultiVector, ...).

JH: Good points. However, the main reason we add this functionality is so we can easily adjoint the solver. The problem with an ODE solver is that it is not self adjointable. By having control of how the forward solution is done, through some sort of scheme, we can more readily generalize the backward solution.
In FEniCS it is also naturally to write your rhs as a form and then choose the time discretization. I do not know the TS library or the Loca/Rythmos libraries well enough to say how well these could interface this abstraction.

(?)

Work Items

Work items:
* Implement AXPY assignment of Functions: DONE
* Add time dependent compiled expressions with a Constant as time variable: DONE
* Add an RKSolver class in DOLFIN (C++) which will be used to step the system: DONE
* Add a ButcherTableau class in Python layer for automatic generation of RK stages: DONE
* Add demo which use the general RKSolver: DONE
* Add a special PointIntegralSolver: TODO

This blueprint contains Public information 
Everyone can see this information.