Skip to content

Split forms on mixed function spaces - #194

Merged
kynan merged 27 commits into
masterfrom
split-forms
Apr 3, 2014
Merged

kynan merged 27 commits into
masterfrom
split-forms

Conversation

@kynan

@kynan kynan commented Mar 13, 2014

Copy link
Copy Markdown
Contributor

When passing a form to the ffc_interface, the form's arguments are
split and a form per block in the mixed space is returned.

Comment thread python/firedrake/ffc_interface.py Outdated

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we really want to unconditionally silence FFC?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is what we used to do in PyOP2, but maybe not.

Comment thread python/firedrake/ffc_interface.py Outdated

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

list of kernels?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In fact, it's more complicated than that right. It's a list of block indices and kernels, yes?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True, that docstring needs updating.

@wence-

wence- commented Mar 17, 2014

Copy link
Copy Markdown
Contributor

This doesn't work for problems cast in nonlinear form:

from firedrake import *

mesh = UnitSquareMesh(2, 2)
V1 = FunctionSpace(mesh, "BDM", 1, name="V")
V2 = FunctionSpace(mesh, "DG", 0, name="P")
W = V1 * V2

# Define variational problem
lmbda = 1
x = Function(W)
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
f = Function(V2)

f.interpolate(Expression("(1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)"))
F = (p*q - q*div(u) + lmbda*inner(v, u) + div(v)*p) * dx - f*q*dx


a = lhs(F)
L = rhs(F)
F = action(a, x) - L

solve(F == 0, x)

Results in:

   solve(F == 0, x)
  File "/data/lmitche1/src/firedrake/python/firedrake/solving.py", line 756, in solve
    _solve_varproblem(*args, **kwargs)
  File "/data/lmitche1/src/firedrake/python/firedrake/solving.py", line 798, in _solve_varproblem
    nullspace=nullspace)
  File "/data/lmitche1/src/firedrake/python/firedrake/solving.py", line 111, in __init__
    self._jac_tensor = assemble(self._problem.J_ufl, bcs=self._problem.bcs)
  File "/data/lmitche1/src/firedrake/python/firedrake/solving.py", line 365, in assemble
    return _assemble(f, tensor=tensor, bcs=_extract_bcs(bcs))
  File "/data/lmitche1/src/firedrake/python/firedrake/solving.py", line 384, in _assemble
    kernels = compile_form(f, "form")
  File "/data/lmitche1/src/firedrake/python/firedrake/ffc_interface.py", line 150, in compile_form
    for forms in FormSplitter().split(form):
  File "/data/lmitche1/src/firedrake/python/firedrake/ffc_interface.py", line 73, in split
    integrand = self.visit(it.integrand())
  File "/data/lmitche1/src/fenics/ufl/ufl/algorithms/transformer.py", line 101, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/data/lmitche1/src/fenics/ufl/ufl/algorithms/transformer.py", line 101, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/data/lmitche1/src/fenics/ufl/ufl/algorithms/transformer.py", line 101, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/data/lmitche1/src/fenics/ufl/ufl/algorithms/transformer.py", line 105, in visit
    r = h(o)
  File "/data/lmitche1/src/firedrake/python/firedrake/ffc_interface.py", line 86, in argument
    if self._idx[o.count()] == i:
KeyError: 0

@kynan

kynan commented Mar 17, 2014

Copy link
Copy Markdown
Contributor Author

@wence- I have fixed that case and added the test you were running. The form splitter wasn't correctly dealing with the case of renumbered arguments.

@wence-

wence- commented Mar 18, 2014

Copy link
Copy Markdown
Contributor

Thanks, here's another one:

from firedrake import *

mesh = UnitSquareMesh(2, 2)
V1 = FunctionSpace(mesh, "BDM", 1, name="V")
V2 = FunctionSpace(mesh, "DG", 0, name="P")
W = V1 * V2

# Define variational problem
lmbda = 1
x = Function(W)
u, p = split(x)
v, q = TestFunctions(W)

F = (inner(u, v) + v[1]*p)*dx

assemble(derivative(F, x))
   assemble(derivative(F, x))
  File "/data/lmitche1/src/firedrake/python/firedrake/solving.py", line 364, in assemble
    return _assemble(f, tensor=tensor, bcs=_extract_bcs(bcs))
  File "/data/lmitche1/src/firedrake/python/firedrake/solving.py", line 383, in _assemble
    kernels = compile_form(f, "form")
  File "/data/lmitche1/src/firedrake/python/firedrake/ffc_interface.py", line 166, in compile_form
    kernel, = FFCKernel(form, name + str(i) + str(j)).kernels
  File "/data/lmitche1/src/PyOP2/pyop2/caching.py", line 69, in __new__
    obj.__init__(*args, **kwargs)
  File "/data/lmitche1/src/firedrake/python/firedrake/ffc_interface.py", line 132, in __init__
    ffc_tree = ffc_compile_form(form, prefix=name, parameters=ffc_parameters)
  File "/data/lmitche1/src/fenics/ffc/ffc/compiler.py", line 172, in compile_form
    pyop2_ir = [generate_pyop2_ir(ir, prefix, parameters) for ir in oir[2]]
  File "/data/lmitche1/src/fenics/ffc/ffc/quadrature/quadraturepyop2ir.py", line 54, in generate_pyop2_ir
    body_ir = _tabulate_tensor(ir, parameters)
  File "/data/lmitche1/src/fenics/ffc/ffc/quadrature/quadraturepyop2ir.py", line 148, in _tabulate_tensor
    opt_par, parameters)
  File "/data/lmitche1/src/fenics/ffc/ffc/quadrature/quadraturepyop2ir.py", line 406, in _generate_element_tensor
    nest_ir, ops = _generate_integral_ir(points, terms, sets, optimise_parameters, parameters)
  File "/data/lmitche1/src/fenics/ffc/ffc/quadrature/quadraturepyop2ir.py", line 599, in _generate_integral_ir
    if len(loop) == 0:
UnboundLocalError: local variable 'loop' referenced before assignment

@wence-

wence- commented Mar 18, 2014

Copy link
Copy Markdown
Contributor

This is simplified from Colin's Skamarock-Klemp model.

@kynan

kynan commented Mar 18, 2014

Copy link
Copy Markdown
Contributor Author

@wence- You need this UFL branch (which still waits on review): https://bitbucket.org/mapdes/ufl/pull-request/8/simplification-rules-required-for-form/

What is happening is that UFL doesn't detect that the first block of this form is effectively 0. This form is then fed to FFC, which does some more analysis and ends up with no terms in its IR. @FabioLuporini's PyOP2 IR currently chokes on this with this rather unhelpful error.

I've added some simplification rules to make UFL detect common cases of "empty" blocks so that we're not even passing those on to FFC.

@kynan

kynan commented Mar 21, 2014

Copy link
Copy Markdown
Contributor Author

Is this now good to go?

@dorugeber

Copy link
Copy Markdown
Contributor

Why did this pass tests? It's not rebased on master, so it doesn't have Lawrence's updated facet tests in. The old test should be failing, now that the FFC fix is in.

@kynan

kynan commented Mar 21, 2014

Copy link
Copy Markdown
Contributor Author

@dorugeber Travis tests the result of the pull request i.e. this branch merged with master.

@dorugeber

Copy link
Copy Markdown
Contributor

Thanks for clarifying that.

I assume this branch is fine, but Doru and I had some strange behaviour when we rebased our facet branch on this. (Rebasing on this was our solution to the integral ordering mismatch, and we found an unexpected extra argument was being passed into the kernel)

It shouldn't block this PR from going ahead, assuming everything works at the moment, but we'll probably be asking you for advice in the next couple of days.

kynan added 15 commits April 2, 2014 11:55
Does not work with disk caching for now since Coefficients can't be
serialised.
Adapt to changes in compile_form, which now returns a four-tuple of
index, domain type, coefficients and kernel.

We now need to pull apart the tensor and test and trial functions if
it is a mixed type and extract the component for the block the kernel
contributes to and filter the boundary conditions to only include
those defined on the current block.
FFCKernel goes back to handling only a single kernel on a single
integral.
Instead of making one pass over the form and getting back a list of
subforms on the blocks of the mixed space, make as many passes as
there are blocks in the mixed space and only enabling one component of
each Argument at a time.

When splitting an Argument, create a UFL ListTensor where only the
component of the MixedFunctionSpace selected by the current block
index is included and all other components are set to Zero. These
Zeros are subsequently eliminated as they are used in expressions and
not propagated up the form tree.
@kynan

kynan commented Apr 3, 2014

Copy link
Copy Markdown
Contributor Author

Can this finally land? The Travis failure is unrelated as you can tell and it passes on buildbot.

kynan added a commit that referenced this pull request Apr 3, 2014
Split forms on mixed function spaces
@kynan
kynan merged commit 70749be into master Apr 3, 2014
@kynan
kynan deleted the split-forms branch April 3, 2014 16:18
connorjward pushed a commit that referenced this pull request Dec 4, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants