syfi team mailing list archive

[Jed Brown <jed@xxxxxxxx>]

On Fri 2008-05-02 10:20, kent-and@xxxxxxxxx wrote:
> > 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.
> 
> 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.

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.

> >
> > 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
> 
> 
> A starting point could be a high-order matrix-free system
> preconditioned with a low order algebraic multigrid cycle.

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

> The hard part in the implementation might be the mapping of the high-order
> 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