dolfin team mailing list archive

[Johan Jansson <johanjan@xxxxxxxxxxxxxxxx>]

Hi!

I've done some more work on the PDE<->ODE solver interface, and I
think it looks very promising. I've just added a module for the Heat
equation as a demo and test for the interface.

The idea is that all equations should be put in the format:

dot(u) = f(u)

This is the form for the Heat equation in such format:

element = FiniteElement("Lagrange", "triangle", 1)

v = BasisFunction(element)
u = Function(element)
f = Function(element)

L = (-dot(grad(u), grad(v)) + f * v) * dx

The solver then consists of an ODE representation which simply
assembles the above form as its right-hand side. The only other thing
the solver does is then time-step the ODE. This could be generalized
so that the same solver could be used for all equations.

This method is evaluated by the right-hand side of the ODE:

void HeatSolver::fu()
{
  FEM::assemble(L, dotu, mesh);
  FEM::applyBC(Dummy, dotu, mesh, element, bc);
}

The right-hand side of the ODE evaluates this and copies between
arrays and Vectors (perhaps the ODE interface should use Vectors as
well, then that step would be eliminated).

The ODE solver can then use either fixed-point iteration or Newton's
method (which is the choice in this case). I've verified the computed
Jacobian matrix for this Heat form (it's constant since the problem is
linear, this is not exploited yet though).

I've solved the same problem as in src/demo/poisson, but
time-dependent (it converges to the same steady-state
solution). Here's the solution in VTK-format (with a Python script for
visualizing it with Mayavi 1.4):

http://ct32.elmagn.chalmers.se/transfer/2005-11-01/heatexample.zip

Here's an image of the solution as well:

http://ct32.elmagn.chalmers.se/transfer/2005-11-01/heat0107.png

Here's a plot (ODE-style) for one of the components:

http://ct32.elmagn.chalmers.se/transfer/2005-11-01/heat-20.eps

Here you can see the residual-based adaptivity a bit, it's not so
clear for this type of problem however (the time step just grows).


I've also converted the updated elasticity module to the same
PDE<->ODE solver interface. I've tested both fixed-point iteration and
Newton's method, but Newton's method doesn't converge so well, so some
more testing is needed there (or perhaps Newton's method is simply not
suitable). The fixed-point method works very well though, here's the
solution of a simple problem I've set up:

http://ct32.elmagn.chalmers.se/transfer/2005-11-01/ko_time-adaptivity01.mpg

And here's the ODE-style plot for the position, velocity and stress
for a point in the body:

http://ct32.elmagn.chalmers.se/transfer/2005-11-01/elasticity-241.eps
http://ct32.elmagn.chalmers.se/transfer/2005-11-01/elasticity-928.eps
http://ct32.elmagn.chalmers.se/transfer/2005-11-01/elasticity-1383.eps

Here the adaptivity is much clearer, the impact produces a fast
transient and the time step is reduced to track it.

I've used the cG(1) method (in time) in both examples, but an
arbitrary order cG(q) or dG(q) method can be selected. For example,
I've tested cG(3) and cG(5) for a similar example, and the result is
that the residual is reduced (and thus the time step is larger for the
same accuracy). But this particular problem is not solved more
efficiently by higher-order methods.


I think the way to move forward is to incorporate such a PDE<->ODE
solver interface in the Solver class, which could then be used by all
modules. We could then work on special cases (time-independent,
linear) as well as on the Newton solver, which could perhaps be
re-implemented using Garth's new PETSc interface for non-linear
equations (which I assume has underwent a lot more testing and has
more features).

What do you think?


As a side note, I also added a Python visualization script for the
Convection-Diffusion solver, here's the solution+script:

http://ct32.elmagn.chalmers.se/transfer/2005-11-01/convdiffexample.zip

And an image:

http://ct32.elmagn.chalmers.se/transfer/2005-11-01/convdiff0010.png

The VTK format (with Mayavi or ParaView) is really great, we should
use it also to store time steps and residuals from the ODE solver for
PDE problems. I haven't tested 3d-visualization yet, but it looks to
be supported even better than for the 2d case.

  Johan