dolfin team mailing list archive

[Jed Brown <jed@xxxxxxxx>]

On Tue 2008-08-19 18:32, Anders Logg wrote:
> On Tue, Aug 19, 2008 at 03:40:13PM +0200, Jed Brown wrote:
> > On Tue 2008-08-19 14:06, Anders Logg wrote:
> > > On Tue, Aug 19, 2008 at 01:49:22PM +0200, Jed Brown wrote:
> > > > On Tue 2008-08-19 13:40, Anders Logg wrote:
> > > > > On Tue, Aug 19, 2008 at 12:12:50PM +0200, Jed Brown wrote:
> > > > > > On Tue 2008-08-19 11:59, Anders Logg wrote:
> > > > > > > On Thu, Aug 14, 2008 at 10:10:03PM +0000, Jed Brown wrote:
> > > > > > > > One way to implement this is to allocate a vector for Dirichlet values,
> > > > > > > > a vector for Homogeneous values, and a Combined vector.  The Homogeneous
> > > > > > > > vector is the only one that is externally visible.
> > > > > > > 
> > > > > > > Isn't this problematic? I want the entire vector visible externally
> > > > > > > (and not the homogeneous part). It would make it difficult to plot
> > > > > > > solutions, saving to file etc.
> > > > > > > 
> > > > > > > Maybe the Function class could handle the wrapping but it would involve a
> > > > > > > complication.
> > > > > > 
> > > > > > Right, by `externally visible' I mean to the solution process, that is
> > > > > > time-stepping, nonlinear solver, linear solvers, preconditioners.  The
> > > > > > vector you are concerned about is the post-processed state which you can
> > > > > > get with zero communication.  It is inherently tied to the mesh and
> > > > > > anything you do with it likely needs to know mesh connectivity.  I don't
> > > > > > think it is advantageous to lump this in with the global state vector.
> > > > > > 
> > > > > > Jed
> > > > > 
> > > > > I don't understand. What is the global state vector?
> > > > 
> > > > The global state vector is the vector that the solution process sees.
> > > > Every entry in this vector is a real degree of freedom (Dirichlet
> > > > conditions have been removed).  This is the vector used for computing
> > > > norms, applying matrices, etc.  When writing a state to a file, this
> > > > global vector is scattered to a local vector and boundary conditions are
> > > > also scattered into the local vector.  The local vector is serialized
> > > > according to ownership of the mesh (you have to do this anyway).
> > > > 
> > > > Jed
> > > 
> > > I'm only worried about how to create a simple interface. Now, one may
> > > do
> > > 
> > >   u = Function(...);
> > >   A = assemble(a, mesh)
> > >   b = assemble(L, mesh)
> > >   bc.apply(A, b)
> > >   solve(A, u.x(), b)
> > >   plot(u)
> > > 
> > > How would this look if we were to separate out Dirichlet dofs?
> > 
> > How about a FunctionSpace object which manages this distinction.
> > Something like the following should work.
> > 
> >   V = FunctionSpace(mesh, bcs); // Is this name clearer?
> 
> We've discussed introducing a FunctionSpace concept earlier (on ufl-dev)
> to handle boundary conditions, and to enable sharing of function space
> data like meshes and dof maps. This might be a good idea, but it has
> to be something like
> 
>   V = FunctionSpace(element, mesh, bcs)

Yes, of course.

> >   u = V.function();
> 
> I think this should be
> 
>   u = Function(V)

Sure, but the representation is fairly closely tied to information in V
which might not otherwise need to be public.

> >   A = V.matrix(a);
> >   P = V.matrix(p);     // preconditioning matrix, optional [1]
> 
> What do these accomplish? Return a matrix of appropriate size? I don't
> think that's necessary since the assembler can set the size.

It's no problem to have the assembler check if the matrix is allocated
and create it if necessary.  My example was just being explicit that
such allocation only needs to happen once.

> >   solver = LinearSolver(A, P);
> >   u0 = u.copy();       // if nonlinear, set initial guess
> > 
> >   V.assemble(A, a, u0); // builds Jacobian matrix [2]
> >   V.assemble(P, p, u0); // preconditioning form, optional
> >   V.assemble(b, L, u0); // the ``linear form'' (aka nonlinear residual)
> >   solver.solve(u, b);
> 
> The problem here is that the solver needs to be aware of Functions.
> Now, a solver only knows about linear algebra, which is nice.

Whatever syntax you like.  The solver is just going to grab the global
representation and work with that, the line could just as well read

  solver.solve(u.x(), b.x());

My suggestion (relating to your other message about regarding
simplicity) is that solver components always call (in the language
presented below) .dofs() on their arguments while output components
always call .vector().  Then the user only sees the high lever Function
object.

> > Note that u0 disappears if the problem is linear and P disappears if you
> > use the real Jacobian as the preconditioning matrix.  The key point is
> > that boundary conditions are built into the function space.  It is no
> > more work, it just happens in a different place.
> > 
> > [1] It is essential that the test/trial spaces for the preconditioning
> > matrix are the same as for the Jacobian.
> > 
> > [2] I don't like the implicit wiring of the state vector into
> > assemble().  It is really hard to track dependencies and the current
> > state is not a property of the Bilinear/Linear forms.
> 
> Maybe one could do
> 
>   (A, b) = assemble(a, L, mesh, bc)
>   solve(A, u.dofs(), b)
> 
> This would be in addition to the current
> 
>   (A, b) = assemble(a, L, mesh)
>   bc.apply(A, b)
>   solve(A. u.vector(), b)
>
> In the first version, the assembler knows about the boundary
> conditions and in the second it doesn't. This would require having two
> members in the Function class that return a Vector:
> 
>   u.vector() // Returns entire vector
>   u.dofs()   // Returns non-Dirichlet dofs

The Function object can keep track of which version is current and do
the necessary scatter to give the version you ask for.

We seem to have a slightly different philosophy of `simple'.  I would
rather write a couple extra lines which reveal the control flow and
enable doing more complicated things.  For instance, I'd rather avoid
magic such as updating a special vector, put into the Bilinear form at
some earlier point, to base the new Jacobian.  Similarly, if you solve
more than one problem then the (A,b) = assemble(a,L,mesh) version cannot
reuse the matrix.  You can always support both versions but I think this
level of overloading makes it harder, rather than simpler, for the user
to keep track of what is happening.  When someone says `how do I do ...'
and the answer is to basically refactor their code rather than just
insert the required line, it looks less like simplicity and more like
magic.

Jed

Attachment: pgpcDaJDnxPkZ.pgp
Description: PGP signature