Subject: [Syfi-dev] Matrix-free form application with tensor products

From: Jed Brown <[email-masked]>
Date: Fri, 2 May 2008 09:42:02 +0200

I'm cross-posting because I don't know which this is more appropriate for,
though everyone probably reads both lists anyway.

How difficult would it be to produce optimized code for the action of a form
instead of computing matrix entries?  In particular, when the elements are high
order tensor products, say order p about 10 in each Cartesian direction, it is quite
cheap (O(p^3) memory, O(p^4) time) to compute the action, but the element matrix
is dense so requires O(p^6) memory.  It is also desirable to have code that is
independent of polynomial order so that in the context of an hp-element method,
the order in each direction may be different and may vary from element to
element.  This forward application should be simpler than computing matrix
entries.

Of course, we'll need some matrix for preconditioning, but we can get this from
a low-order (i.e. Q1) discretization on the nodes of the high order version.
This will produce a sparse matrix which is spectrally equivalent to the
(unformed) high order one.  I've recently done some experiments with this method
in the context of a Chebyshev collocation scheme (here application costs only O(p^3
log p) due to the FFT) and found it to be very effective.  The trouble is that
assembling the sparse approximation requires describing your problem twice, once
as a tensor product and again as a finite difference (better with collocation)
or traditional finite element (better for high order Galerkin) scheme.  If we
can remove this duplication, I suspect it would make a huge difference in the
adoption of hp-element methods.

Any thoughts?

Jed

---

From: [email-masked]
Date: Fri, 2 May 2008 10:20:46 +0200 (CEST)

> element.  This forward application should be simpler than computing matrix

It should be fairly easy to create a matrix that rely on a UFC form, some
BC and
a mesh that computes the matrix-vector product my running over the elements
in the mesh. In Dolfin the Matrix class would be derived from GenericMatrix.
Hence, this approach is independent of whether you use FFC of SyFi/SFC.
But UFC and also most of SFC and FFC assumes (for now) that the element is
the same troughout the mesh. This is something we would like to fix in the
future
and you are welcome, if you want to begin.

> (unformed) high order one.
> I've recently done some experiments with this
> method

A starting point could be a high-order matrix-free system
preconditioned with a low order algebraic multigrid cycle.
The hard part in the implementation might be the mapping of the high-order
elements to the low order elements.

Kent

---

From: Jed Brown <[email-masked]>
Date: Fri, 2 May 2008 11:24:39 +0200

> and you are welcome, if you want to begin.

My understanding is that FFC/SFC generate matrix elements.  I don't want matrix
elements.  I want the application of the form to a vector.  For any particular
form, it is easy to write the tensor product application independently of
polynomial degree.  I'm mostly interested in hexahedral elements.  (Even if the
mesh has some tets and prisms, it should be possible to use them sparingly so that
they don't contribute too many entries to the global matrix.)  Here, one only
needs 1D basis functions and quadrature points for each order, say 2 to 16.
Then the function which applies the form would take the order triplet (p,q,r)
where 2 <= p,q,r <= 16, the element vector, and the result vector.  What FFC/SFC
do is much more complicated.  I don't know the internals, so I don't know which
would be the most appropriate to implement this feature.

I'm happy using PETSc directly, so I'll keep using MatShell for this type of
operation.  My question is really with regard to generating the element
application function.  Of course, we can use the same form file to generate
matrix elements for the low-order discretization, or we could use a different
one which leaves out some terms in exchange for nicer matrix structure.

> preconditioned with a low order algebraic multigrid cycle.

This is exactly what I've been doing.  It works quite well.

> elements to the low order elements.

This is the advantage of using a nodal basis.  There is a natural duality between
high order hexahedral tensor product elements and Q1 elements on the nodes.

Jed

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 2 May 2008 12:34:30 +0200

> high order hexahedral tensor product elements and Q1 elements on the nodes.

Let me try to summarize what's already there.

In FFC, I think the action of a form can be accomplished by simply changing
  u = TrialFunction(element)
to
  u = Function(element)
such that u becomes a coefficient, the vector to multiply with.

In SFC, you can write
  b = Action(a)
if a is a Form. This will simply do the element matrix-vector
product symbolically prior to compilation.
(This hasn't been used or tested very well, but it is short and simple)

SFC supports hexahedral elements, FFC does not.

Adding mappings to the form compilers shouldn't be too hard with quadrature.

UFC doesn't cover high order geometries, varying cell types in a mesh
or varying element orders. This is probably the largest obstacle.

--
Martin

---

From: Anders Logg <[email-masked]>
Date: Fri, 2 May 2008 12:56:01 +0200

> or varying element orders. This is probably the largest obstacle.

Yes, if you have a bilinear form a in FFC that gives you the matrix,
just change the trial function argument to a Function and you will get
a linear form L which you can assemble into a vector.

You can assemble this into a PETScVector and then extract the PETSc
Vec pointer by x.vec() and then wrap that into a MatShell.

Or you can let DOLFIN handle this and just overload the mult()
function in PETScKrylovMatrix (which will handle the communication
with PETSc through a MatShell).

Since FFC and SyFi both generate UFC code, you can choose either of
the two to generate the code for the action.

Then as Martin and Kent point out, the main problem may be that all
elements need to be of the same type. This might change in future
versions of UFC.

-- 
Anders

---

From: Jed Brown <[email-masked]>
Date: Fri, 2 May 2008 12:57:52 +0200

> such that u becomes a coefficient, the vector to multiply with.

Oh, right.  We can use a linear form to compute the action of the bilinear form
without getting matrix entries.  So this part is already implemented and well
tested.  Great.

> (This hasn't been used or tested very well, but it is short and simple)

How do you think this would compare to using a linear form?  For high order
elements, the symbolic expression could be so huge the compiler has trouble
doing the necessary CSE.

> SFC supports hexahedral elements, FFC does not.

It sounds like I can get what I want if I write a form file, specify high order
hex elements, just have a linear form, and specify quadrature.  Can I expect
that applying this form will be O(p^4)?

> or varying element orders. This is probably the largest obstacle.

Indeed.  We can do some proof of concept work without these things, but they are
clearly desirable in the long term.  It would greatly complicate the spec,
though.  Would it be better to create a new HpUFC with these features?  Of
course, if hp-adaptivity is made easy, everyone will be doing it :-)

Jed

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 2 May 2008 14:25:04 +0200

> elements, the symbolic expression could be so huge the compiler has trouble

It shouldn't be very different, but I think those problems will be the
same any way.

Frankly, I don't think SFC will cope well with p=10 anyway.
I plan to fix some bottlenecks to make it more robust.
However, there are some things I need UFL to be able to do,
so it may take a little time.

> that applying this form will be O(p^4)?

Note that SFC doesn't yet have form files in the same sense as FFC does.
I'm working on the unified form language (UFL), but it's not finished.
There are some SFC demos in the newest syfi under examples/solvers/.

> course, if hp-adaptivity is made easy, everyone will be doing it :-)

That is a good idea, and it can be merged back in UFC for UFC 2.0.

--
Martin

---

From: Kent-Andre Mardal <[email-masked]>
Date: Mon, 05 May 2008 11:36:02 +0200

> so it may take a little time.

I have just printed out the basis functions on a unit-cube p=10. 
It looks like: 

fe.N(0)=(1.0/9.3329542656000000e+10)*
( 602595.0*(y*(y*y*y)*(y*y*y))*pow( y-1.0,3.0)+87498.0*y*pow( y-1.0,9.0)+4536.0*pow( y-1.0,10.0)+87498.0*(((y*y)*(y*y))*((y*y)*(y*y)))*pow( y-1.0,2.0)+1996867.0*(y*y*y)*pow( y-1.0,7.0)+602595.0*(y*y)*pow( y-1.0,8.0)+3558504.0*(y*(y*y)*(y*y))*pow( y-1.0,5.0)+4536.0*(y*((y*y)*(y*y))*((y*y)*(y*y)))*( y-1.0)+3558504.0*((y*y)*(y*y))*pow( y-1.0,6.0)+1996867.0*((y*y*y)*(y*y*y))*pow( y-1.0,4.0))*( 602595.0*pow( x-1.0,3.0)*(x*(x*x*x)*(x*x*x))+4536.0*pow( x-1.0,10.0)+602595.0*pow( x-1.0,8.0)*(x*x)+1996867.0*pow( x-1.0,7.0)*(x*x*x)+4536.0*( x-1.0)*(x*((x*x)*(x*x))*((x*x)*(x*x)))+3558504.0*pow( x-1.0,6.0)*((x*x)*(x*x))+1996867.0*pow( x-1.0,4.0)*((x*x*x)*(x*x*x))+87498.0*pow( x-1.0,9.0)*x+3558504.0*pow( x-1.0,5.0)*(x*(x*x)*(x*x))+87498.0*pow( x-1.0,2.0)*(((x*x)*(x*x))*((x*x)*(x*x))))*( 602595.0*(z*(z*z*z)*(z*z*z))*pow( z-1.0,3.0)+87498.0*z*pow( z-1.0,9.0)+87498.0*(((z*z)*(z*z))*((z*z)*(z*z)))*pow( z-1.0,2.0)+1996867.0*(z*z*z)*pow( z-1.0,7.0)+602595.0*(z*z)*pow( z-1.0,8.0)+3558504.0*(z*(z*z)*(z*z))*pow( z-1.0,5.0)+4536.0*(z*((z*z)*(z*z))*((z*z)*(z*z)))*( z-1.0)+3558504.0*((z*z)*(z*z))*pow( z-1.0,6.0)+1996867.0*((z*z*z)*(z*z*z))*pow( z-1.0,4.0)+4536.0*pow( z-1.0,10.0))

The expression is a bit large and its more than a thousand of them so
CSE probably is not a good idea. However, one can run CSE on a group of
functions to start with and then go trough all in a systematic fashion.
SFC will not do this automatically, but it should still be feasible.
Doing this in the "proper way" could also merit a paper. It should
also be possible to pre-compute the quantities in quadrature points
etc.  

SyFi is sufficiently fast that it is possible to experiment with
different algorithms here. Computing the element functions, computing
one entry in the element matrix etc typically takes less than 1s, even
with this analytical approach. 

> That is a good idea, and it can be merged back in UFC for UFC 2.0.

We would _really_ like to have hp in UFC, but it might be that the best
way to start is a HpUFC. This would alow experimentation. I think this
could be a very good project. 

Kent