Moving degree estimation into compute_form_data

[Lawrence Mitchell]

When the symbolic processing pipeline is turned on, the estimated polynomial degree of forms swiftly goes through the roof. Especially with piola-mapped elements on manifolds where the form ends up with terms like (J^T J)^-1. I think it would be good to add a quadrature-selection step to compute_form_data that attaches estimated degrees of the non-transformed form after derivative expansion has been applied, but before pullbacks.

How best to do this? At the moment I do:

    # Estimate polynomial degree of integrands now, before applying
    # any pullbacks and geometric lowering.  Otherwise quad degrees
    # blow up horrifically.
    def gather_and_estimate_degrees(form):
        from collections import defaultdict
        from ufl.utils.sorting import canonicalize_metadata
        from operator import add
        itgs = defaultdict(list)
        for itg in form.integrals():
            md = hash(canonicalize_metadata(itg.metadata()))
            itgs[(itg.integral_type(),
                  itg.ufl_domain(),
                  itg.subdomain_id(),
                  md)].append(itg)
        new_itgs = []
        for itg_list in itgs.values():
            itg = itg_list[0].reconstruct(reduce(add, [i.integrand() for i in itg_list]))
            md = itg.metadata().copy()
            degree = md.get("quadrature_degree", "auto")
            if degree == "auto":
                degree = estimate_total_polynomial_degree(itg.integrand())
            md["quadrature_degree"] = degree
            new_itgs.append(itg.reconstruct(metadata=md))
        return Form(new_itgs)

    form = gather_and_estimate_degrees(form)

But this is somewhat different from what used to happen. In particular, the form compiler pipelines appear to look at the built integral data which does some merging/splitting of integrands. Can one pull that step forward before the pullback application and then add degree estimation then?

Thoughts?

[Martin Sandve Alnæs]

I think it needs to happen between some of the current passes in compute form data, in particular after the first apply derivatives call but probably before geometry lowering.

Also, what you do here is almost but not quite consistent with the way integrals are combined in ufl, in particular the "everywhere" and "otherwise" subdomain ids need special treatment.

[Lawrence Mitchell]

Yes, I currently do this after apply_derivatives but before geometry lowering.

At this point, I think there are no "otherwise" subdomain ids. I looked at using rearrange_integrals_by_single_subdomains but I realised I can't just use that out of the box, because it expects never to see an input integral containing an "otherwise" subdomain. If I understand things correctly, I think I want to effectively do the rearrangement that happens in the guts of build_integral_data. But that function is only expecting to be called once. How about this:

Split out the integral rearrangement from the building of integral data in build_integral_data. So that I convert the input form into a "canonical" grouped by integrals/subdomain ids type. Now I estimate the polynomial degree for these integrals (this could happen optionally in compute_form_data). Now, after the rest of the pipeline I now build the integral data, but this is now simplified, because the rearrangement has already happened.

[Martin Sandve Alnæs]

Yes, that sounds good.

To make sure you understand why it can only be called once, consider:

M = f*dx(0) +  f**2*dx(1) + g*dx

# is equivalent to
M == f*dx(0) +  f**2*dx(1) + g*dx("everywhere")

# is transformed to
M == (f+g)*dx(0) + (f**2 + g)*dx(1) + g*dx("otherwise")
# where "otherwise" means subdomain not in (0,1)

In ufc and the dolfin assembler "otherwise" is called the default integral, i.e. it's used for any subdomain where no specific integral is defined in the form. Adding integrals to the form at this stage would make this definition invalid, and instead of just "otherwise" we would need to represent "domain - subdomains 0 and 1".

[Lawrence Mitchell]

OK, I'll have a go over the next few days.

[Lawrence Mitchell]

See https://bitbucket.org/fenics-project/ufl/pull-requests/40/split-form-grouping-out-of/diff

[Lawrence Mitchell]

I think this is now also resolved, see https://bitbucket.org/fenics-project/ufl/src/8014fbf900e1b8025cc91e48f1ecdb8da8554c53/ufl/algorithms/compute_form_data.py?at=master&fileviewer=file-view-default#compute_form_data.py-248