dolfin team mailing list archive

[Jed Brown <jed@xxxxxxxx>]

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?
  u = V.function();
  A = V.matrix(a);
  P = V.matrix(p);     // preconditioning matrix, optional [1]
  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);

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.

Attachment: pgpu_Qef7JfS3.pgp
Description: PGP signature