dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Nov 02, 2005 at 12:30:38AM +0100, Johan Jansson wrote:

> > 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.
> > 
> 
> Ah, I forgot about it for this problem. It doesn't affect the steady
> state, so I didn't catch it. In the updated elasticity case I lump
> it. I guess it should be lumped for the fixed-point case and a matrix
> for the Newton case.

Sounds like a good plan.

> > 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).
> 
> I think these cases imply a lot more work in the compiler though, to
> catch unreasonable cases as well as to detect how the terms should be
> rearranged to be input to the ODE solver. I think it's also good. But
> if it turns out it's straightforward to implement, then why not.

I think the compiler is a good place to check this. There, we have the
symbolics and it's much less work to implement in Python than in C++.

We could require that at most one term contains a time derivative. The
compiler would just need to check which term that is and treat it
specially. I think it could be done without too much trouble,

But until we have figured out how we want to do it, the approach you
have used is a very good start.

> > 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?
> 
> Yes, operator() is terribly inefficient. However, if you assume the
> uniprocessor case, then you can get the array anyway, which is done in
> several places elsewhere, i.e. Vector::array(). 
> 
> > 
> > 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.)
> 
> Ok, this is a good solution.

ok, let's do it.

> > 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.
> > 
> 
> This is a good start.

Should the PDE class be generated as a class called PDE in the
namespace of the form? For Poisson, that would mean that we have a
class

    Poisson::PDE

which would contain a pair of forms

    Poisson::BilinearForm
    Poisson::LinearForm

Or should we replace the namespace Poisson with a class named Poisson
that is a subclass of PDE and put the two forms within its namespace.
The first alternative would probably work better with SWIG and it's
also more consistent with the current naming:

    Poisson::BilinearForm a;
    Poisson::LinearForm L(f);
    Poisson::PDE pde;

What do you think?

/Anders