dolfin team mailing list archive

["Johan Hoffman" <jhoffman@xxxxxxxxxxx>]

>> > 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 think the most intuitive way for the user to provide the variational
form is to use a time derivative Dt(u). Then the desciption of the
(variational form of the) equation would be something like what you
describe above, which is very close to mathematical notation.

The next step is to determine how to discretize the equation. Today there
are several options avaliable in DOLFIN:

(1) Semi-discretization time->space: this is what is done today in DOLFIN
for PDEs, although the time discretization is done manually. The natural
next step would be for FFC to do the discretization in time from Dt(u),
given a specification in the form-file, either as a space-time element:

space_discretization = FiniteElement("Lagrange", "tetrahedron", 1)
time_discretization = FiniteElement("discontinuous Lagrange", "interval", 2)
space_time_discretization = space_discretization + time_discretization

for cG(1)dG(2), or as a semidiscretization:

space_discretization = FiniteElement("Lagrange", "tetrahedron", 1)
time_discretization = FD_scheme("Crank-Nicolson")
space_time_discretization = space_discretization + time_discretization

for cG(1) with Crank-Nicolson in time. Functions may then be defined using
the space-time discretizations. The same output files from FFC to DOLFIN
would then be generated as is done today; that is computation of element
matrices and vectors.

(2) Semidiscreization space->time: this would be along the lines what
Johan is doing now; computing f(u) to use in the ODE solver (for dot u =
f(u)). We would then still declare a space element:

space_discretization = FiniteElement("Lagrange", "tetrahedron", 1)

but the time derivative Dt(u) would not be discretized by FFC, but instead
would an output file for the ODE solver in DOLFIN be created that contain
the evaluation of f(u).

(3) We also might like to generate the vector F(u) for using the non
linear solver on F(u)=0, although this can be done already today with FFC
defining an L in this form. For Newton we would also like to generate the
Jacobian automatically, which will demand some new functionality.

(4) A next step is moving meshes. Johan is doing this for elasticity I
guess in a pure Lagrangian framework, but we also want to allow for ALE
methods where the mesh is moving independetly from the solution. This
corresponds to a coordinate transformation so that the derivatives and
integrals (dx,ds etc.) needs to be modified according to a specified map.
This would probably be quite easy to do in FFC by implementing
map-dependent operators, where derivatives and integrals are interpreted
using the chain rule wrt a given mapping.

>> > 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 PDE also know about the discretization? For example, if the forms
are on the format (1) or (2); that is if the forms are prepared for the
ODE solver or if it should be assembled to a matrix? So that the PDE class
corresponds to a discretized PDE?

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

In the first alternative, we would still keep a base class called PDE in
the DOLFIN kernel I guess? As we do for BilinearForm etc. I think this is
the best solution, rather than introducing classes for all different
problems that are generated.

/Johan

>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
>