← Back to team overview

dolfin team mailing list archive

Re: Applying Dirichlet conditions by removing degrees of freedom (was [Fwd: Re: [HG DOLFIN] merge])

 

On Thu, Aug 14, 2008 at 10:10:03PM +0000, Jed Brown wrote:
> 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.

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.

-- 
Anders

Attachment: signature.asc
Description: Digital signature


Follow ups

References