dolfin team mailing list archive

[kent-and@xxxxxxxxx]

> 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

I have never had any trouble with the additional diagonal matrix due to
Dirichlet conditions,  but I have mostly used multigrid. In K.-A. Mardal,
and R. Winther. Uniform Preconditioners for the Time Dependent Stokes
Problem, Numerische Mathematik 98(2):305--327, 2004.
we tested multigrid methods for Stokes type problems with Taylor-Hood,
P2-P0, Mini
and Crouzeix-Raviart elements (the CR experiments did not end up in the
final paper
but the number can be found in my thesis.) The condition number was below
20 for all
methods independent of mesh size. I have done similar experiments with AMG.

Kent