Subject: [UFL-dev] polynomial order of form

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 13 Feb 2009 11:53:19 +0100

Hello,

Any thoughts on how to implement a better rule for determining the polynomial
order of a form?

as it is now the forms:

a = v*u*dx
a = f*v*u*dx
a = f*f*f*.....*v*u*dx
a = inner(grad(v), grad(u))*dx

all result in the same order (2 for linear basis functions)

Also, what is the strategy for handling the form

a = (inner(grad(v), grad(u)) + v*u)*dx ?

will the order of the form be determined by the highest order, or will it be
possible to split the integrand according to order? The latter could get
complicated for forms like:

a = f*(inner(grad(v), grad(u)) + v*u))*dx

I do notice the # TODO in formdata.py, I'm just curious.

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Fri, 13 Feb 2009 13:06:09 +0100

> all result in the same order (2 for linear basis functions)

What do we want? Add the degrees of all functions and basis functions
that are multiplied? Should be possible to do.

> a = f*(inner(grad(v), grad(u)) + v*u))*dx

If I got it right, lhs and rhs can handle splitting of arbitrary forms by arity,
so it should be possible to do the same for terms with different order.

> I do notice the # TODO in formdata.py, I'm just curious.

Note that I use the convention
FIXME: High priority issue
TODO: Low priority issue

Since this obviously can be postphoned until after we get things
working in general.

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 13 Feb 2009 13:20:07 +0100

> that are multiplied? Should be possible to do.

Yes, I think it will be enough if we do it like that, taking into account
derivatives too.
 
> so it should be possible to do the same for terms with different order.

OK, sounds great. I think this is exactly the kind of functionality I was
looking for.
 
> FIXME: High priority issue

OK, good to know.
 
> Since this obviously can be postphoned until after we get things

Depends how you define 'working in general', I would say getting a correct
integration of a = f*f*f*f*v*dx is part of 'working in general' :)

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Fri, 13 Feb 2009 13:32:53 +0100

> integration of a = f*f*f*f*v*dx is part of 'working in general' :)

'working in general' ~ using ufl and dolfin to solve a few pdes
more complicated than projection/poisson (I had those working
for christmas but sfc is currenty broken).

Integration order should be possible to override by form compiler
options anyway, which should be fine as a first step.

Martin

---

From: Kent Andre <[email-masked]>
Date: Fri, 13 Feb 2009 13:38:30 +0100

> >> that are multiplied? Should be possible to do.

In general I don't think this is what you 'want' to do. 
It is common to use the same quadrature rules for 
a = inner(grad(v), grad(u))*dx
and 
a = f*inner(grad(v), grad(u))*dx

even though f is a function. 

It would however be nice to choose quadrature rule by setting some kind
of meta data. 

Kent

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 13 Feb 2009 12:42:28 +0000

> It would however be nice to choose quadrature rule by setting some kind

We discussed attaching this data to the integration operators a while 
ago and there was general agreement.

Garth

> Kent
> 
> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 13 Feb 2009 14:03:00 +0100

Quoting "Garth N. Wells" <[email-masked]>:

> 
> 
> Kent Andre wrote:
> > On fr., 2009-02-13 at 13:32 +0100, Martin Sandve AlnÃ¦s wrote:
> >>> Quoting Martin Sandve AlnÃ¦s <[email-masked]>:
> >>>

Sure, but I do believe the default should be correct integration.
If it is possible to have different orders it might have an impact on the
performance of the code which will allow the default for the compilers to be
'correct AND fast' code.
Through the meta data the user is given the opportunity of changing the rules,
typically decreasing the order to gain speed by sacrificing accuracy.
 
> ago and there was general agreement.

Yes, and it would not be a problem to override the default values if any
instructions are set.

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Fri, 13 Feb 2009 14:30:00 +0100

> Yes, and it would not be a problem to override the default values if any

Attaching metadata to an integral is implemented, but we
haven't  decided what format the metadata should be in.

This currently works:

metadata = "arbitrary data"
metadata2 = "other arbitrary data"

a = f*dx(0, metadata) + g*dx(0, metadata2) + g*dx(0, metadata)

which will result in a Form with two Integrals with integrands
(f+g) and (g) since the first and last term has the same domain and data.

Martin

---

From: Kent Andre <[email-masked]>
Date: Fri, 13 Feb 2009 14:32:50 +0100

> Sure, but I do believe the default should be correct integration.

As it is now it is correct in the sense that the integration is accurate
up to a certain order and this accuracy is consistent with the order of
the element independently of the number of coefficients or variational
form.

Anyway, I think this is just a decision to be made. One can argue in
both directions. I'm fine with you choice. 

Kent

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Fri, 13 Feb 2009 14:37:56 +0100

Either way, we can have both algorithms available so the user can
either switch between the two or specify the exact order.

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 13 Feb 2009 15:06:35 +0100

> (f+g) and (g) since the first and last term has the same domain and data.

It looks good and I think that we're fully covered if we combine the above with
the algorithms to split the integrand according to order.
Again, I was just curious to know which functionality that comes with UFL but I
think we can safely close this discussion.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Sun, 15 Feb 2009 22:10:00 +0100

> (f+g) and (g) since the first and last term has the same domain and data.

I think there needs to be a convention as part of UFL for how to
specify the integration rule. Otherwise different form compilers will
invent different conventions and interpret the metadata differently.

Here's a suggestions:

  metadata = ("quadrature", degree)
  metadata = ("tensor", degree)

Degree would be the degree of polynomials which will be integrated
exactly by the quadrature rule. For tensor representation we use
quadrature on the reference element so it makes sense there as well.

It might even be useful to add a class IntegrationRule with "enums"
for quadrature and tensor representation. Then an Integral would have

  _integrand
  _measure

and Measure would have

  _domain_type
  _domain_id
  _integration_rule

and possibly metadata.

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 16 Feb 2009 09:57:26 +0100

> invent different conventions and interpret the metadata differently.

Agree.

>  metadata = ("tensor", degree)

Here you're already into FFC territory...

Integration order is an important mathematical
property, the method for code generation is not.
"tensor" is in the latter category.

> for quadrature and tensor representation. Then an Integral would have

In the current implementation _integration_rule == _metadata, but it's
of course not a problem to replace that with two variables or rename it.

I think that probably any metadata other than integration order
will be form compiler specific. We can let Measure have

  _domain_type
  _domain_id
  _integration_order
  _metadata

and allow:

  dx(domain_id)
  dx(domain_id, integration_order)
  dx(domain_id, integration_order, metadata)

where metadata contains compiler hints like "tensor" and
can be form compiler specific.

Then the integration order can be interpreted correctly by
any form compiler, and metadata can be ignored.

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 16 Feb 2009 10:04:51 +0100

> any form compiler, and metadata can be ignored.

This sounds like a very good solution.
If the metadata is irrelevant to a given form compiler it could throw an error,
or maybe a warning is better, rather than just ignoring it?

Kristian

---

From: Anders Logg <[email-masked]>
Date: Mon, 16 Feb 2009 10:10:43 +0100

> "tensor" is in the latter category.

Hints for the degree of the quadrature is also form compiler
specific. There may be other ways to evaluate the integral than
quadrature, for example by random sampling, symbolic integration (like
in SyFi) or using a tensor representation (FFC).

-- 
Anders

-- 
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 16 Feb 2009 09:28:09 +0000

> in SyFi) or using a tensor representation (FFC).

Tensor representation still uses quadrature, it's just that it's 
performed a priori rather than at runtime.

Garth

> 
> 
> ------------------------------------------------------------------------
>

---

From: Anders Logg <[email-masked]>
Date: Mon, 16 Feb 2009 10:32:22 +0100

> performed a priori rather than at runtime.

Yes, but it may potentially be performed without use of quadrature.

I think the suggestion by Martin is ok. I just want to point out that
quadrature is not universally used by all form compilers.

-- 
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 16 Feb 2009 10:34:44 +0100

> performed a priori rather than at runtime.

I guess the degree is just the highest polynomial order that the user wishes
integrated exactly.
If it is possible for the form compiler to comply with this wish it may do so,
if not it will ignore the instruction.
On second thought, I think we do need 'quadrature', 'tensor' to be allowed in
the metadata, otherwise the user can't specify different representations on
different subdomains.

Kristian

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 16 Feb 2009 10:09:55 +0000

> the metadata, otherwise the user can't specify different representations on

Perhaps this information should belong to the form rather than the integral?

     a_tensor = . . .
     a_quad   = . . .
     a = a_tensor + a_quad

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 16 Feb 2009 11:27:25 +0100

> Perhaps this information should belong to the form rather than the integral?

I think it belongs in the integral. All terms need an integral
(measure) so it's practical to attach it to the integral.

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 16 Feb 2009 11:36:11 +0100

>> the metadata, otherwise the user can't specify different representations on

That's why I used the term "integration order" and not "quadrature order".
Basically, by this the user will ask the form compiler to use an integration
method that guarantees a certain order.

I'm fine with lumping this in with other metadata, but I think at least how to
specify the integration order should be standardized. I don't know how many
other options we can standardize, but we can define a format for it:

metadata = {
    # ... The meaning of some parameters is defined by UFL:
    "integration_order": None,
    # ... Other packages use their name to wrap specific options:
    "ffc": ffc_options,
    "sfc": sfc_options,
}

So you can do e.g.
    ffc_options = dict(representation="quadrature")
    sfc_options = dict(integration_method="quadrature")
    metadata = dict(integration_order=3, ffc=ffc_options, sfc=sfc_options)
    a = v*u*dx(0, metadata) + f*v*dx(0)

At any point we can add new standardized options.
I don't see any value in letting representation/integration_method
be a standard parameter, since the values it can take
will be different for each compiler anyway.

> Perhaps this information should belong to the form rather than the integral?

As Kristian said, sometimes different integrals or integral terms has
different optimal compiler options, so attaching options to integrals
is a good feature to have.

Form-global compiler options doesn't have to be attached to UFL objects,
isn't it better to pass them to the form compiler jit function directly?

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 16 Feb 2009 11:48:02 +0100

> method that guarantees a certain order.

I totally agree.
 
> will be different for each compiler anyway.

This looks like a very good solution, any form compiler just needs to get the
relevant options dictionary and disregard all others.
The options supported in these dictionaries should then be described in the
relevant compiler manual. Using this approach UFL doesn't have to know anything
about compiler options, apart from which compilers define these options?

Agree.

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 16 Feb 2009 11:56:32 +0100

> about compiler options, apart from which compilers define these options?

UFL doesn't need to know anything about compiler options,
not even which compilers exist.

Form compilers should just do this for a given integral:

m = integral.metadata()

# Get UFL-defined options
integration_order = m.get("integration_order")
...

# Get my own options
ffc_options = m.get("ffc", {})

Or iterate over what's given:

for k,v in m:
    if k == "ffc": ...
    elif k == "integration_order": ...
    else: warning("Unrecognized option %s" % k)

Martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 16 Feb 2009 12:02:20 +0100

> relevant options dictionary and disregard all others.

I also think this looks good.

-- 
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 16 Feb 2009 12:11:12 +0100

Quoting Martin Sandve Alnæs <[email-masked]>:

> >> >> If it is possible for the form compiler to comply with this wish it may
> do
> >> so,
> >> >> On second thought, I think we do need 'quadrature', 'tensor' to be
> allowed
> >> in
> >> >> the metadata, otherwise the user can't specify different
> representations
> >> on
> >> >> different subdomains.
> >>
> >> I'm fine with lumping this in with other metadata, but I think at least
> how
> >> to

True, I misread the
'dict(integration_order=3, ffc=ffc_options, sfc=sfc_options)'
as a call with named arguments.

Kristian