dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Nov 01, 2005 at 10:54:23PM +0100, Johan Jansson wrote:
> 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

Very good!

How do you treat the scaling with the mass matrix?

The mono-adaptive solver handles implicit systems like Mu' = f so we
should be able to include it properly.

Another option would be to introduce a time derivative in FFC and
include it in the form:

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

or alternatively one could throw the Laplacian in the right-hand side
like so
 
  a = v*Dt(u)*dx
  L = (dot(grad(v, grad(u)) + v*f)*dx

but the user can decide where to put it (either assemble matrix or do
action).

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

I've thought about doing it, but we probably need both since
Vector::operator() is not very efficient compared to indexing a plain
array (I think). Can we specify an interface for the ODE solver that
lets the user choose one of two interfaces?

One way would be to overload

    virtual void f(const real u[], real t, real y[]);
    virtual void f(const Vector& u, real t, Vector& y);

and then implement default versions where one of them calls the other
and does conversions with toArray/fromArray.

Since the array one is the default right now, we could add the Vector
version and in the default implementation of the standard version
call the Vector version with proper conversions. (Similar to what we
do in ComplexODE.)

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

Looks great!

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

There can also be a bug in the computation of the Jacobian. I haven't
done any extensive testing.

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

I suggest we add a new class PDE (it's already there) and build some
intelligence into this class. In particular, it should know if the
problem is time-dependent or static. As a start, if FFC sees that a
pair of forms (a, L) is specified it could collect the two forms into
a PDE.

And we should of course make as much use as possible of Garth's new
interface for nonlinear problems, in the ODE solvers to begin with.

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

We should probably try to retire the Octave/MATLAB visualization and
move everything to VTK.

/Anders