dolfin team mailing list archive

["Johan Hoffman" <jhoffman@xxxxxxxxxxx>]

> 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
>

This looks very good, I am not worried so much about the overhead. The
reason I am still interested in the ability of forming linear systems to
solve in DOLFIN for some problems, is that I do not know if all problems
can be solved efficiently by using the explicit approach of keeping the
whole f(u) on the right hand side of the equation as in the fixed point
iteration you indicate.
I know you can use damped fixed point to help the problem in case it is
not suited for pure explicit update, and I know about the short time step
stabilization for stiff problems. Does this solve all issues? Say for
example a strong non linearity in f(u), can you still use an explicit
approach?

The alternative is then to use Newton, and I think it is great that the
Newton solver in the ODE solver is in good shape. For large problems,
where you cannot store the whole jacobian, does it have the option of
storing sparse versions, or using a matrix-free representation of the
jacobian?

Is it the PETSc non-linear solver you are using together with the ODE
solver today?

I think it is great to go for a general discrete solver of ODE/Newton type
in DOLFIN, where you just feed the right hand side f(u); it is ideal for
automation. I just wonder if you see any limits to this approach in terms
of size or type of problems? Or can we skip the whole business with
forming linear systems of equations Ax=b?

I believe you Johan estimated the work for a matrix-free approach based on
the action of f(u), to be favourable to forming the matrix A and using
matrix-vector multiplication Ax? In general?

/Johan