dolfin team mailing list archive

[Jed Brown <jed@xxxxxxxx>]

On Thu 2008-08-14 20:42, Anders Logg wrote:
> On Wed, Aug 13, 2008 at 09:20:17PM +0000, Jed Brown wrote:
> > On Wed 2008-08-13 18:31, kent-and@xxxxxxxxx wrote:
> > > 
> > > > Kent-Andre Mardal wrote:
> > > >> BTW: the Assembler code that computes A, b and enforces BC
> > > >> simultaneously in a symmetric way,
> > > >
> > > > Very nice. We can now implement some symmetric linear solvers.
> > > >
> > > > Do you keep the 'Dirichlet' dofs and zero the row and column, or do you
> > > > eliminate them from the global system?
> > > 
> > > The Dirichlet dofs are kept. I thought that was simplest to do..
> > 
> > It is certainly more work at this stage, but there are significant
> > benefits to actually removing the degrees of freedom.  In particular, it
> > allows more preconditioners to be effective and makes related matrices
> > (i.e. the Schur complement) better conditioned.
> > 
> > Jed
> 
> Do you have any numbers/references to support this? I was under the
> impression it didn't matter that much.
> 
> It would be a complication to remove dofs (and then put them back in).
> There may be an elegant solution, but it would break the simple
> 
>   A = assemble(a, mesh)
>   b = assemble(L, mesh)
>   solve(A, u.vector(), b)

It can make a huge difference, but we don't normally notice it since
preconditioners are really good at correcting the perturbation.  Here
are a few numbers considering the velocity system for a 50x50 point
Stokes problem.  I calculated these with

  -ksp_rtol 1e-10 -ksp_gmres_restart 500 -pc_type none
  -ksp_monitor_singular_value

and observing the estimate near convergence.  The column `A' below is
the estimated condition number of the unconditioned velocity operator,
`A/ILU' is when preconditioned with ILU(0), `S' is the condition number
of the Schur complement (-B_1 A^{-1} B_2').

N        A         A/ILU    S
0(0)     20000     25       37.8
2(1)     28900     25       38.3
10(5)    28900     24.9     44
100(50)   30000     24.9     100
500(250) 32000+    24.5     773

As you can see, the conditioning of A degrades significantly by just
fixing the velocity at one point.  The Schur complement does not respond
so dramatically at first, but when we fix an entire edge of the boundary
(N=100), the conditioning is terrible.

So as long as we have explicit representations around and we use a
standard preconditioner (ILU is really good for this, AMG is okay) then
it is okay to enforce boundary conditions without eliminating degrees of
freedom.  On the other hand, the Schur complement is dense and
notoriously difficult to precondition, but it is expensive to apply
since every multiplication involves a solve with A (or at least a
multigrid V-cycle).

Note that for slip boundary conditions, it is also important to remove
degrees of freedom, see Bänsch and Höhn 2000, 'Numerical treatment of
the Navier-Stokes equations with slip boundary conditions.'


The mechanics of removing degrees of freedom is much easier in the
nonlinear setting.  We decompose the state vector into a Homogeneous and
Dirichlet part and consider the function

  F(x^H) = A(x^H) (x^D + x^H) - (b^D + b^H)

which is evaluated elementwise and needs essential boundary conditions
present.  On the other hand, the linear solve involves only the block
A^H corresponding to x^H.  That is, we solve

  A^H x^H = -F_H

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.  At the beginning of
function evaluation, scatter Dirichlet and Homogeneous values into the
Combined vector and proceed as usual.  During assembly, simply apply the
Combined -> Homogeneous map when inserting.  This is also completely
natural for the matrix-free operator application, zero the Combined
vector, Homogeneous -> Combined, evaluate action locally, scatter
Combined -> Homogeneous.

Jed

Attachment: pgp7TfH2qUg34.pgp
Description: PGP signature