syfi team mailing list archive

[Jed Brown <jed@xxxxxxxx>]

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