dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Aug 19, 2008 at 07:44:00PM +0200, Jed Brown wrote:
> 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

Yes, I have a tendency towards the "magic" end of the spectrum.
I admit this may cause trouble sometimes.

On the other hand, I don't want things to look like PETSc where one
needs to write large chunks of code for simple things.

-- 
Anders

Attachment: signature.asc
Description: Digital signature