← Back to team overview

dolfin team mailing list archive

Re: PDE<->ODE interface

 

On Fri, Nov 04, 2005 at 09:15:32AM -0600, Anders Logg wrote:

...

> The ODE solvers are starting to come in pretty good shape, at least
> the standard mono-adaptive solvers.
> 
> Johan Jansson has done some benchmarking for elasticity and the
> overhead was very small compared to explictly forming the system by
> hand and solving. I don't remember the exact numbers, maybe Johan can
> give us the details?
> 
> /Anders
> 

Sure. In my elasticity solver I have implemented dG(0) (aka backward
Euler) with fixed-point iteration. This is a very simple algorithm, so
there are not many ways to do it, and we can assume that I have not
introduced any extra overhead. The algorithm looks like this:

while(R > TOL)
  R = -U1 + U0 + k * f(U1)
  U1 = U1 + R

Where "R" is the discrete residual.

I compared my simple implementation to using the general ODE solver,
also with dG(0). I forced both solvers to take 10 iterations for each
time step and both solvers to take a fixed time step of 1e-3 over a
time interval of 5s. The problem I'm solving is a spinning cube which
is represented as 3630 dofs in the ODE solver. Here are the results
("simple" is my implementation, "ode" is the general ODE solver):

               Computation time   # evaluations of f()
simple dG(0):  125s               50000
ode dG(0):     132s               50000
ode cG(1):     144s               55001

Here we can see that the overhead of using the general ODE solver is
6%. I also included cG(1) (aka Crank-Nicolson), which requires an
extra evaluation of f() every time step, but incurs no additional
overhead.

In this test I turned off residual computation for the general ODE
solver (since my simple solver doesn't compute the residual). It would
have added an extra evaluation of f() for every time step for dG(0),
but no extra evaluations for cG(1).

An overhead of 6% is negligible, and it gives the freedom to choose
any cG(q)/dG(q) method with residual computation and adaptivity. So I
cannot see a reason not to use the ODE solver included in DOLFIN for
the fixed-point case.

I have not benchmarked the Newton case yet, but I can't see why it
would incur any extra overhead. The ODE solver in DOLFIN can compute
the Jacobian df/du numerically, and efficiently as far as I can see,
but there's also the option of supplying your own function for
df/du. In that case the overhead should be negligible as well compared
to implementing Newton's method yourself with a known df/du. Perhaps I
can perform a similar benchmark with the heat module and Newton's
method.

  Johan



Follow ups

References