syfi team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Fri, May 02, 2008 at 12:34:30PM +0200, Martin Sandve Alnæs wrote:
> 2008/5/2 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
> 
> 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.

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