Subject: [UFC-dev] Thoughts on hp-adaptivity and constraints

From: Jed Brown <[email-masked]>
Date: Fri, 9 May 2008 01:50:07 +0200

The biggest obstacle to hp-adaptivity is removing all the constraints from the
global system.  In the general case, the mesh is both h-nonconforming (hanging
nodes) and p-nonconforming (different element order).  Furthermore, we would
like to remove degrees of freedom on Dirichlet boundaries, even in the case that
only part of a vector component is constrained as in the case of slip
boundaries.

I think this requires some rethinking of the managing of degrees of freedom.  I
find it is useful to distinguish faces as being h-dominant or h-subordinate.
All h-conforming faces are h-dominant, but where the mesh is h-nonconforming,
only the largest face is h-dominant.  We can pick a basis for this face by
making its order equal to the lowest order of all the touching elements.  Such
faces may not agree at their edges, so we must identify h-dominant edges and
pick a basis for them similarly.  In this way, we can consider every element to
have its faces constrained against the associated h-dominant faces.  At
conforming interfaces, this constraint is the trivial one.  The process of
forming the stiffness matrix amounts to mapping the element basis to the global
basis by applying these constraints.

This is also an opportunity to make matrix-free application natural, especially
considering that sparse linear algebra is memory bandwidth limited and
increasing p makes memory usage much worse when the real system matrix is
actually formed.  Doing this requires that managing constraints be fast because
it must be done on every matrix application.  To illustrate this, consider the
following application function where I have assumed that the basis can be
decomposed into interior and boundary modes (this is typically true).

for each dominant face:
 for each edge:
   global, dirichlet -> unconstrained face (apply edge constraints)
 global, dirichlet -> face interior (unconstrained)

for each element:
 for each element face:
   associated dominant face -> unconstrained element (apply face constraints)  (*)
 global -> element interior (unconstrained)
 element operation (can be tensor product)                                     (**)
 element interior -> global
 for each element face:
   unconstrained element -> associated dominant face (apply face constraints)  (*)

for each dominant face:
 for each edge:
   face -> global (apply edge constraints)
 face interior -> global (unconstrained)

The complexity is critical in these places:
 (**)  p^6 space, p^6 time with matrix, p^3 space, p^4 time by tensor product
 (*)   p^4 space, p^4 time with matrix, p^2 space, p^3 time by tensor product

Note that entering the full system matrix generates p^6 entries.  For
preconditioning, we'll want to construct a sparse element matrix (we should be
able to do this with a form built against the existing UFC spec).  Global
assembly will require applying constraints, but this code is pretty much already
written if we are doing matrix-free application.  (An obvious optimization is to
only do matrix-free application for sufficiently large p, say p > 4.)

An HpUFC implementation should provide a way to apply the element constraint
matrix for edges and faces in each of the non-conforming cases.  For instance,
when faces meet like:

|------|-------|
|      |       |
| A    |  B    |
|      |-------|
|      |       |
|      |  C    |
|------|-------|

where A is order (p1,q1), B is order (p2,q2), and C is order (p3,q3), the
central edge will have order q0 = min(q1,q2,q3) and we must be able to apply the
constraint for each of the face.  These constraint matrices depend on the
degree of the edge, the degree of the element, and the geometry.  For quad
faces, there are only 3 cases for geometry (plus rotations).  Since the entries
in the constraint matrix are not important, the interface should just provide
application, thus allowing it to be done matrix-free when appropriate
(constraining 3D elements against faces).

The last component is actually evaluating the element contribution.  A very
general interface for the action of the element operator takes a geometry flag
(tet,hex,prism), a tuple describing the anisotropic element order, private local
data (like functions at quad points), and the element trial function.  Forming
the corresponding element matrix requires the geometry flag, order tuple, and
private local data.  As noted before, this element matrix could be sparse,
arising for example from a low-order discretization of a similar form on the
nodes of the high order element basis.

If the constraints are handled efficiently, the kernel becomes application of
the element operator and the preconditioner.  These have the advantage of being
much less memory bandwidth limited than sparse matrix multiplication.

A benefit of this is that we also get a *real* spectral method.  That is, if the
solution is very smooth, we could do pure p-refinement.  The tensor product
method scales extremely well.  I have done some work with polynomial order
(100,100,100) in 3D and (1000,1000) in 2D.  The time to solution for a Poisson
problem is a little more than double the time for the same number of degrees of
freedom with centered finite differences and geometric multigrid.  This is a
huge win considering the difference in accuracy.

What are your thoughts?  I know these are a lot of changes, but it may be
easier to do it all at once instead of incrementally.  I'm by no means an expert
on the form compilers.  How feasible are these requirements?

My long-term research is modeling ice flow (non-Newtonian, nonlinear slip
boundary conditions, free surface).  I wrote a spectral code does this (minus
the free surface) but of course geometry is important so in the end a pure
spectral code is just a toy.  I want to move to an hp-element method without
sacrificing the good properties of the spectral method (like fast tensor product
application and sparse preconditioning).  This isn't possible in any of the
existing libraries without a significant amount of work.  I think the form
compilers are a great tool to avoid lots of error-prone boilerplate.  Also, it
makes the libraries accessible to many more users.  Obviously the finite element
library needs appropriate data structures for the nonconforming elements.  The
closest I've seen to a natural framework for this method is probably Sieve.  The
core of Dolfin is very simple so it may not be so hard to add the functionality.
Once I have a clearer perspective, I'll get to implementing.  As for HpUFC, it
is probably reasonable to draft a spec and write some reference implementations
by hand for a couple problems.  This is fairly easy to do, using FIAT to
tabulate a 1-D basis for a range of orders.

Comments?

Jed

---

From: Kent-Andre Mardal <[email-masked]>
Date: Thu, 22 May 2008 10:29:02 +0200

> matrix for edges and faces in each of the non-conforming cases.  For instance,

I am sorry about the late response. I've been swamped lately. I think
this sounds very interesting although I don't understand it all yet. 

I have a few basic questions related to this subject. 

Concerning the p refinement couldn't we make this discretization
conforming by constructing elements with variable degree across edges. 
I mean, it is simple to create e.g. a P2 element which is P1 on one
edge. This could be done easily with both Fiat and SyFi. With this
approach the addition to the current UFC would be that different
elements would be used a different places in the mesh and that 
the this is chosen based on a criteria. Does this make sense ? 

Concerning constraining the hanging nodes. Is this done with Lagragian 
multipliers ? If so, is it easy to create preconditioners ?  If not 
this is probably a nice thing to study ? 

Kent

---

From: Jed Brown <[email-masked]>
Date: Thu, 22 May 2008 12:42:22 +0200

> the this is chosen based on a criteria. Does this make sense ? 

My preference is not to do this because it makes hanging node constraints
tricky.  Also, it makes the element operation more expensive because it breaks
up the regular tensor product structure.  For relatively low order, this may be
okay, but it is actually a lot harder to implement (it's really an
``optimization'' of the scheme described below).  Also note that you don't gain
nearly as much from tensor products with triangle/tetrahedral elements as with
quad/hex.

> this is probably a nice thing to study ? 

There are two approaches, one with Lagrange multipliers and one by removing
constrained degrees of freedom.  Lagrange multipliers are probably easier to
implement, but they make the system larger, negatively impact condition number,
and make the problem more difficult to precondition.  There are various
intermediate approaches where most constrained degrees of freedom are removed,
but some are enforced with Lagrange multipliers.  I would rather remove them
all.

For this, I like to distinguish between the element degrees of freedom on a face
(which do not correspond directly to global degrees of freedom) and the
canonical degrees of freedom on the adjoining h-dominant face (which do
correspond to global degrees of freedom).  The continuity constraint means that
we can write the element dofs on a face as a linear combination of the face
dofs.  This linear transformation tells us how to convert an element stiffness
matrix written in terms of element dofs to a contribution to the global matrix.
Note that if the face is geometrically conforming (no hanging nodes) and the
polynomial degrees match, this constraint is just the identity matrix, which can
be a special optimization.  With this procedure, the effective element basis is
actually the reduced space dependent on the basis for it's neighboring faces,
but all element operations take place as though it was the full space.

In fact, the face should also be constrained against it's edges, so the same
procedure works recursively.  In general, a mesh entity with topological
dimension n>0 is constrained with regard to it's (h-dominant) neighbors of
topological dimension n-1.  That is, all interior dofs correspond directly to
global degrees of freedom.  All boundary dofs are written in terms of the dofs
for the neighbors with topological dimension n-1.

I hope that helps.  I'm happy to discuss more.

Jed

---

From: Anders Logg <[email-masked]>
Date: Thu, 22 May 2008 17:28:00 +0200

> matrix for edges and faces in each of the non-conforming cases.  For instance,

It's certainly doable, but it requires a major effort from someone
with a good knowledge of DOLFIN + UFC + the form compilers. If we do
it in smaller steps, we could all learn as we go along.

Last time we did something as big as this was when we moved FFC +
DOLFIN to the UFC interface. This was a big pain and we needed to
rewrite large parts of the code.

One of the steps would be to add support for quad meshes, which should
be fairly simple since the mesh data structure should already support
it (in principle).

Then we could add an experimental hpAssembler next to Assembler that
does assembly for hp elements with a special hpufc form.

-- 
Anders

---

From: [email-masked]
Date: Tue, 27 May 2008 00:17:02 +0200 (CEST)

> tricky.

I meant to remove the hanging nodes in this way.

> Also, it makes the element operation more expensive because it
> breaks

Ok, it might be more expensive in this way, but again you do not need to
remove
the constraints. But anyway, it was just a thought. I think the HpUFC
maybe could be general enough to support both apporaches. This approach
would basically
only mean that UFC support different elements throught the mesh, which
is needed anyway.

> but some are enforced with Lagrange multipliers.  I would rather remove

ok.

> but all element operations take place as though it was the full space.

Ok, let me see if I understand.

UFC needs to support different elements througout the mesh.
Furthermore, in an assembly loop a standard element matrix is computed on
an element, then this element matrix is multiplied with a constraint
matrix to enforce the continuity. This constraint matrix can be computed
based on the (different) elements associated  with an edge. Is this
correct ?

Kent

> _______________________________________________
> UFC-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufc-dev
>

---

From: Jed Brown <[email-masked]>
Date: Tue, 27 May 2008 13:11:49 +0200

> I meant to remove the hanging nodes in this way.

Right, the advantage of working on the full element space (apply constraints
separately) is that you can have a highly optimized kernel for evaluating
residuals and applying the Jacobian.  If you build it into the space, then this
kernel needs to know (anisotropic) orders of all it's faces, and the faces'
edges and their relative location (h-nonconforming).  This effectively kills the
tensor product (at least unless you are very clever) which is the key to making
the hp version fast.  This is why I think it's the wrong way to proceed if you
would like to accommodate high order.

> remove the constraints. But anyway, it was just a thought. I think the HpUFC
> maybe could be general enough to support both apporaches. This approach would
> basically only mean that UFC support different elements throught the mesh,
> which is needed anyway.

I think it would not be too hard to support both since adapting the element
basis requires setup code which is somewhere between the externally managed
continuous basis and a discontinuous basis.  But I think managed mode should
come first because it is easier to implement and permits highly optimized
element kernels.

> based on the (different) elements associated  with an edge. Is this

Yes, and this constraint matrix can also be applied by tensor product for
appropriate choices of element shape basis.

I think it would be hard to design a HpUFC spec until we have a reference
implementation of something similar to look at.  I'm working on one which is
intended to be fast for tensor product hexahedral elements.  When I get it to a
useful stage, I'll make it available.  At this point, I'm leaning towards using
Maxima and a mangler to produce the needed code since there is a lot of near
duplication if it is done by hand.  That is, suppose we are working with the
Stokes problem for non-Newtonian rheology or fully implicit Navier-Stokes.
This Jacobian has saddle point structure (dF/du)(u) = J = [A B'; B 0].
We want functions which:

* evaluate F
* set up Jacobian for unassembled application (stored values at quad pts)
* apply J, J', B, B', A, A'
* assemble approximations to A, S = -B A^-1 B'  (for the preconditioner)

In a library, we can use these functions to build many saddle-point
preconditioners, but these components are problem specific (though structurally
similar) and they only depend on F.  I can get all the pieces easily and
generically from Maxima and I also get `optimize' to produce a minimal set of
values which need to be stored at quadrature points.  (I think it is not
desirable to project the stored values back into the function space because it
reduces accuracy and requires more work when applying the Jacobian unassembled.
If you just store the value at which the Jacobian is evaluated, it's even more
work, but the accuracy is not affected.)

I'd like to keep the discussion open as to how to do all of this with the right
amount of generality.  I have a vision for one way to put this together, but I
certainly don't have all the details worked out and I'm sure there are things I
have overlooked.

Jed

---

From: Anders Logg <[email-masked]>
Date: Tue, 27 May 2008 20:42:39 +0200

> certainly don't have all the details worked out and I'm sure there are things I

Bangerth and Kayser-Herold recently wrote a paper about how to handle
the constraint matrices. It's a while since I read it but from what I
can remember it sounds similar to what you suggest:

  www.isc.tamu.edu/publications-reports/tr/0704.pdf

I have a good feeling it will appear in TOMS shortly...

I agree it will be beneficial to have a proof of concept
implementation that we can look at before we modify the UFC interface,
in particular if you can try and make your code "UFC-like" from the
start.

For DOLFIN, there are two major things we need to handle (in addition
to the constraint matrices), which are (i) quad and hex meshes and
(ii) varying cell types throughout the domain.

-- 
Anders

---

From: Jed Brown <[email-masked]>
Date: Tue, 27 May 2008 21:55:26 +0200

>   www.isc.tamu.edu/publications-reports/tr/0704.pdf

I have read this paper, but I think we can do better in a couple of ways.

The biggest obstacle to general implementations of p- and hp-element codes is
that the global matrix has a lot of entries.  Clever modal basis choices can
make it sparse for linear problems, but fails for nonlinear problems.  This
means that unassembled evaluation of the Jacobian is desired which means that
constraint enforcement needs to be fast and can't be done by manipulating the
global matrix.  The paper describes eliminating constraints by manipulating the
assembled matrix.  This is a bad thing to do in serial and much worse in
parallel.  They find that removing constrained degrees of freedom costs as much
as global assembly (and of course the solver takes nearly forever because there
are lots of entries).  This means that their method basically won't work for
unassembled application because you have to enforce constraints at every Krylov
iteration.

My idea is to replace all those global operations with purely local optimized
projections implemented using a table of precomputed 1-D projection matrices.
This gives great data locality and very little branching.  It's possible that it
will not work well, but the asymptotics are strongly on my side.  For
sufficiently low order, simple assembly is bound to still win.  The central
question is where the crossover point occurs, but I won't know until I
implement.  (Of course we can split the Jacobian into an assembled part for the
low order elements and an unassembled part for the high order elements.)

Another hp-element implementation to be aware of is Concepts which chooses
enrichment of deficient elements instead of constraining abundant elements.

  http://www.concepts.math.ethz.ch/

and a paper, also at esaim-m2an.org  vol 36, no 5

  http://e-collection.ethbib.ethz.ch/eserv/eth:25713/eth-25713-01.pdf

Also, Philipp Frauenfelder's thesis goes into a lot of detail on the
implementation.

Jed

---

From: [email-masked]
Date: Tue, 27 May 2008 21:59:16 +0200 (CEST)

>> > > up the regular tensor product structure.  For relatively low order,
>> this
>> > > may be
>> > > ``optimization'' of the scheme described below).  Also note that you
>> don't
>> > > gain
>> > > nearly as much from tensor products with triangle/tetrahedral
>> elements as
>> > > with
>> > > quad/hex.
>> > >
>> >
>> > > For this, I like to distinguish between the element degrees of
>> freedom on
>> > > a face
>> > > correspond to global degrees of freedom).  The continuity constraint
>> means
>> > > that
>> > > we can write the element dofs on a face as a linear combination of
>> the
>> > > face
>> > > matrix written in terms of element dofs to a contribution to the
>> global
>> > > matrix.
>> > > Note that if the face is geometrically conforming (no hanging nodes)
>> and
>> > > the
>> > > polynomial degrees match, this constraint is just the identity
>> matrix,
>> > > which can
>> > > be a special optimization.  With this procedure, the effective
>> element
>> > > basis is
>> > > actually the reduced space dependent on the basis for it's
>> neighboring
>> > > faces,

Agree that proof of concept may be good.
Anyway, in UFC we have the following function

    virtual void tabulate_tensor(double* A,
                                 const double * const * w,
                                 const cell& c) const = 0;

Could an addition be a similar function

   virtual void tabulate_constraint_tensor(double* A,
                                 const double * const * w,
                                 const cell& c0, const cell& c1) const = 0;

where c0 and c1 are the cells meeting at edge e.

the assemble loop will the be:

for cell in domain
     tabulate_tensor(..);
     for ncell in cell.neighbourhood()
          tabulate_constraint_tensor(..);
     multiply_constraints(...);

Is this correct ?

Kent

---

From: Jed Brown <[email-masked]>
Date: Tue, 27 May 2008 22:26:18 +0200

>                                  const cell& c0, const cell& c1) const = 0;

I think we need to tabulate on the (h-dominant) facets.  I think it would take a
facet number (i.e. 0-3 for tet, 0-5 for hex) and orientation (right/left handed,
rotation, size/translation).

This is what I have in mind for Jacobian application:

for each dominant face:
 for each edge:
   global -> unconstrained face (apply edge constraints)
 global -> face interior (unconstrained)

for each element:
 for each element face:
   associated dominant face -> unconstrained element (apply face constraints)
 global -> element interior (unconstrained)
 element operation applied by tensor product         (*)
 element interior -> global
 for each element face:
   unconstrained element -> associated dominant face (apply face constraints)

for each dominant face:
 for each edge:

Maybe it's useful to separate out the problem-specific part from that which is
only basis dependent.  In this case, only (*) is problem dependent.  The key
basis-specific part is projecting a function on facet dofs to element dofs (can
be done exactly because of the minimum rule) and projecting a function on
element dofs onto facet dofs (not exact because the facet is lower order, but
still consistent).  For preconditioning matrix assembly, everything below (*) is
applied to the rows and columns of the element matrix with appropriate hoisting
out of the element loop.

Jed