Subject: [FFC-dev] [HG FFC] Use new UFL function estimate_max_polynomial_degree to estimate quadrature

From: [email-masked] (FFC)
Date: Sun, 19 Apr 2009 18:23:03 +0200 (CEST)

One or more new changesets pushed to the primary ffc repository.
A short summary of the last three changesets is included below.

changeset:   1497:38268dcd85de2bfcae1094ed96ea584293320491
tag:         tip
user:        Anders Logg <[email-masked]>
date:        Sun Apr 19 18:23:00 2009 +0200
files:       ffc/compiler/uflcompiler.py
description:
Use new UFL function estimate_max_polynomial_degree to estimate quadrature
degree. Seems to give wrong results for quadrature representation.

changeset:   1496:2eaeb80e1d1fdf45208005bd1cc23016ca872919
user:        "Garth N. Wells <[email-masked]>"
date:        Sun Apr 19 14:26:22 2009 +0800
files:       ffc/compiler/uflcompiler.py
description:
Fix silly error.

changeset:   1495:4710241c0954645d8150ec6fcd1eb3bff4f6d333
user:        "Garth N. Wells <[email-masked]>"
date:        Sun Apr 19 13:47:10 2009 +0800
files:       ffc/compiler/uflcompiler.py
description:
Small fix.

----------------------------------------------------------------------
For more details, visit http://www.fenics.org/hg/ffc

---

From: Anders Logg <[email-masked]>
Date: Sun, 19 Apr 2009 18:26:15 +0200

I tried to use the new UFL function estimate_max_polynomial_degree
but the results seem to be incorrect for quadrature representation.
(Or the references are wrong which I doubt.)

Kristian, can you take a look at this?

Try running the test in simple_verify_tensors. It reports the
following errors:

FunctionOperatorsBilinearForm_cell_integral_0:      *** (diff = 0.010098)
MassBilinearForm_cell_integral_0:                   *** (diff = 1e-09)
NeumannProblemLinearForm_cell_integral_0:           *** (diff = 0.0029685)
NeumannProblemLinearForm_exterior_facet_integral_0: *** (diff = 0.00292)
OptimizationBilinearForm_cell_integral_0:           *** (diff = 2.19169)
PoissonDGBilinearForm_interior_facet_integral_0:    *** (diff = 0.3504)
PoissonSystemLinearForm_cell_integral_0:            *** (diff = 0.0029685)
QuadratureElementBilinearForm_cell_integral_0:      *** (diff = 0.249892)
QuadratureElementLinearForm_cell_integral_0:        *** (diff = 0.10106)
StokesLinearForm_cell_integral_0:                   *** (diff = 0.00140453)

-- 
Anders

On Sun, Apr 19, 2009 at 06:23:03PM +0200, FFC wrote:
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 20 Apr 2009 10:59:02 +0200

> NeumannProblemLinearForm_cell_integral_0:           *** (diff = 0.0029685)

This difference is because:
element = VectorElement("Lagrange", triangle, 1)

v = TestFunction(element)
f = Function(element)
L = inner(v, f)*dx

results in:
Input form: { sum_{i_0} ((v_0)[i_0]) * ((w_0)[i_0])  } *
dx0<{'ffc_representation': 'tensor', 'quadrature_order': 1}>

where the quadrature order should be 2.

The tensor representation computes the degree elsewhere (which is 2) as far as I
can tell from using the -d1 option.

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 20 Apr 2009 11:05:12 +0200

Try applying this to the form after expand_compounds etc has been applied.

In [20]: estimate_max_polynomial_degree(L.form_data().form)
Out[20]: 2

In [21]: estimate_max_polynomial_degree(L.form_data().original_form)
Out[21]: 1

This is because rules have not been defined for all types.

In SFC I don't look at the original form at all,
using form.form_data().form from the start.

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 20 Apr 2009 11:25:45 +0200

> using form.form_data().form from the start.

OK, thanks.

Kristian

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 20 Apr 2009 11:26:40 +0200

> StokesLinearForm_cell_integral_0:                   *** (diff = 0.00140453)

This should be fixed now such that we only get the expected differences.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Apr 2009 13:06:52 +0200

> Try applying this to the form after expand_compounds etc has been applied.

If expand_compounds is needed by an algorithm, then wouldn't it be
better if that algorithm called expand_compounds?

In general I find it the use of form, form.form_data() and form_data.form 
a bit confusing.

I would like the program flow to be to call form.form_data() once and
not having to call it in multiple places. Even validate_form is
calling form.form_data().

-- 
Anders

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

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 20 Apr 2009 13:21:07 +0200

Yes, I agree that the data flow can be confusing,
and I've intended to discuss it for a while.
But there's so much to do and so little time, we
can't do any major restructuring right now anyway...

In SFC, I do this as the first thing:

  fd = input_form.form_data() # not recomputed if done before
  form = fd.form # has applied all common expansions and renumberings

>From that point I use fd and form, never touch input_form
again, and never call form_data() again. Do that in FFC, and
you should be safe.

Most algorithms work on expanded expressions,
but that could be documented better. Some assume
already expanded expressions, and some call
expand functions, and this is not consistent.
Don't hesitate to ask or report if something specific is unclear.

Martin

On Mon, Apr 20, 2009 at 1:06 PM, Anders Logg <[email-masked]> wrote:
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAknsV0wACgkQTuwUCDsYZdEYsACfeuTvMVgE/UAyn0Hmv3ve6kn8
> HeEAoIuHs4p/1Zdj+HguzFaerxZbYi/9
> =PhO4
> -----END PGP SIGNATURE-----
>
>

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Apr 2009 13:27:41 +0200

> again, and never call form_data() again. Do that in FFC, and

We do something similar, but it gets a bit ugly since we want to
modify the integral metadata (set correct representation and
quadrature order if not already set by a user). So then we first call
form_data() only to be able to estimate the quadrature order. Then we
need to call form_data() again on the modified form after having
modified the integral metadata.

  # Validate form
  validate_form(form)

  # Set metadata for all integrals (requires form_data().form)
  new_form = Form([])
  for integral in form.form_data().form.integrals():
      new_form += Form([_set_integral_metadata(integral, options)])
  form = new_form

  # Extract form data
  form_data = form.form_data()
  form = form_data.form

> Don't hesitate to ask or report if something specific is unclear.

Is it possible to either set a flag that can be checked by algorithms
requiring expanded expressions, or that these algorithms call
expand_compounds?

It's confusing if something is required by an algorithm and it doesn't
report an error if that requirement is not met.

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Apr 2009 13:28:50 +0200

> This should be fixed now such that we only get the expected differences.

Garth, can you check that the example you had works now?

Then we can go ahead and remove the old stuff.

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 20 Apr 2009 13:39:06 +0200

> modified the integral metadata.

This calling of form_data on the second form will result in
a new computation, which is costly.

If the metadata is the only thing you modify, can't
you store the modified metadata outside the form?

> report an error if that requirement is not met.

That's a solution I have considered, but not found the time to go into detail.
Such flags must be possibly be propagated many places, the implementation
could be very error prone if not done carefully.

Martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 20 Apr 2009 13:46:19 +0200

> > > NeumannProblemLinearForm_cell_integral_0:           *** (diff =

Once this has been confirmed and *.form removed, we should generate new
references using the *.ufl forms such that the last differences disappear.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Apr 2009 13:55:39 +0200

> references using the *.ufl forms such that the last differences disappear.

Yes!

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Apr 2009 14:02:21 +0200

> you store the modified metadata outside the form?

Yes, good idea. I changed this so it's now stored as a dictionary in
form_data named metadata. It can be accessed by

  form_data.metadata[integral]

FFC attaches a number of extra data to form_data after it has been
created by UFL. Some is FFC specific, such as form_data.ffc_elements
(FFC wrappers for FIAT elements corresponding to the UFL elements) but
some others may be worthwile to include in ufl.FormData.

Take look at analyze_form in ffc/compiler/uflcompiler.py if some of
these can be placed in ufl.FormData.

> could be very error prone if not done carefully.

Yes, let's check back on this later, but it's very important to put in
some checks, especially since the functions sometimes return something
that looks valid but is incorrect (like wrong integrand degree).

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 20 Apr 2009 14:10:12 +0200

> that looks valid but is incorrect (like wrong integrand degree).

We can remove the default behaviour of DegreeEstimator
and implement rules for all operators that should be supported.

Martin

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 20 Apr 2009 18:38:47 +0100

> Garth, can you check that the example you had works now?

The debugging output that I was looking at has changed (a day or two ago 
'dx(0, .... degree) was printed to the screen). If the current debug 
line 'q = XX' is the polynomial order of the form, then for the tests 
that I made it is working correctlt.

> Then we can go ahead and remove the old stuff.

Very good. Remember that the option '-O' should be used in combination 
with the quadrature representation (which is the default representation).

Garth

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

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Apr 2009 19:44:15 +0200

> that I made it is working correctlt.

Yes, that's the polynomial degree. I left it in there while debugging.
The debugging output will be cleaned up as well.

> > Then we can go ahead and remove the old stuff.

I inperpret this as everything set for starting to break compatibility
with the old .form format.

> with the quadrature representation (which is the default representation).

What do you mean? Should optimization be enabled by default?

-- 
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 20 Apr 2009 18:46:35 +0100

> What do you mean? Should optimization be enabled by default?

Probably, but we need to be able to turn it off (or control certain 
types of optimisations seperately) since it can be slow for complicated 
forms.

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 20 Apr 2009 19:49:07 +0200

> types of optimisations seperately) since it can be slow for complicated 

ok. If it is potentially very slow then it sounds like it should not
be enabled by default. But we can enable it when compiling the demos
in DOLFIN.

-- 
Anders