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

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

Quoting Anders Logg <[email-masked]>:

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

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

> 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:
> > 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
> > _______________________________________________
> > FFC-dev mailing list
> > [email-masked]
> > http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

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