dolfin team mailing list archive

[Johan Jansson <johanjan@xxxxxxxxxxxxxxxx>]

On Tue, Nov 01, 2005 at 04:31:15PM -0600, Anders Logg wrote:
> On Tue, Nov 01, 2005 at 10:54:23PM +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.

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

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

  Johan