syfi team mailing list archive

[kent-and@xxxxxxxxx]

> 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.

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

Kent