The assertion error tells you that TSFC was expecting things to be in reference space, but some modified terminal was not actually in reference space. This is not so surprising given that pullbacks were disabled, because it is very easy to leave physical space terms in your UFL form. Instead of skipping apply_function_pullbacks altogether, I would just try to write the mass term (which is what you want to pull back manually) in all reference terms, that is: - ReferenceValue(f) instead of f - ReferenceGrad instead of grad - ReferenceDiv instead of div - etc. and then hope that UFL can handle this. "Shipton, Jemma" <j.shipton@imperial.ac.uk> writes:
Thanks Lawrence, I've tried this... I've attached some code for the shallow water equations with hand coded pullbacks in the mass term (other terms don't need them). If I set do_apply_function_pullbacks=False in tsfc/ufl_utils.py then running the code (sometimes) gives the error below. When it doesn't error, the solver just fails to converge. I can't figure out why the behaviour changes.
Any suggestions most welcome! Thanks!
Jemma
Traceback (most recent call last):
File "/Users/jshipton/vanilla_firedrake/src/PyOP2/pyop2/caching.py", line 197, in __new__
return cls._cache_lookup(key)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/tsfc_interface.py", line 53, in _cache_lookup
return cls._cache.get(key) or cls._read_from_disk(key, comm)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/tsfc_interface.py", line 73, in _read_from_disk
raise KeyError("Object with key %s not found" % key)
KeyError: 'Object with key 216f38281a96e8d31c8f18aa7b153641 not found'
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "no_pullbacks.py", line 44, in <module>
solver = LinearVariationalSolver(prob)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/variational_solver.py", line 314, in __init__
super(LinearVariationalSolver, self).__init__(*args, **kwargs)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/variational_solver.py", line 161, in __init__
options_prefix=self.options_prefix)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/solving_utils.py", line 122, in __init__
form_compiler_parameters=self.fcp)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/assemble.py", line 145, in create_assembly_callable
collect_loops=True)
File "</Users/jshipton/vanilla_firedrake/lib/python3.6/site-packages/decorator.py:decorator-gen-49>", line 2, in _assemble
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/utils.py", line 61, in wrapper
return f(*args, **kwargs)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/assemble.py", line 214, in _assemble
kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/tsfc_interface.py", line 205, in compile_form
number_map, interface).kernels
File "/Users/jshipton/vanilla_firedrake/src/PyOP2/pyop2/caching.py", line 199, in __new__
obj = make_obj()
File "/Users/jshipton/vanilla_firedrake/src/PyOP2/pyop2/caching.py", line 189, in make_obj
obj.__init__(*args, **kwargs)
File "/Users/jshipton/vanilla_firedrake/src/firedrake/firedrake/tsfc_interface.py", line 115, in __init__
tree = tsfc_compile_form(form, prefix=name, parameters=parameters, interface=interface)
File "/Users/jshipton/vanilla_firedrake/src/tsfc/tsfc/driver.py", line 61, in compile_form
kernel = compile_integral(integral_data, fd, prefix, parameters, interface=interface)
File "/Users/jshipton/vanilla_firedrake/src/tsfc/tsfc/driver.py", line 153, in compile_integral
integrand = ufl_utils.split_coefficients(integrand, builder.coefficient_split)
File "/Users/jshipton/vanilla_firedrake/src/tsfc/tsfc/ufl_utils.py", line 211, in split_coefficients
return map_expr_dag(splitter, expression)
File "/Users/jshipton/vanilla_firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File "/Users/jshipton/vanilla_firedrake/src/ufl/ufl/corealg/map_dag.py", line 84, in map_expr_dags
r = handlers[v._ufl_typecode_](v)
File "/Users/jshipton/vanilla_firedrake/src/tsfc/tsfc/ufl_utils.py", line 147, in _modified_terminal
return self.modified_terminal(o)
File "/Users/jshipton/vanilla_firedrake/src/tsfc/tsfc/ufl_utils.py", line 181, in modified_terminal
assert mt.reference_value
AssertionError
I also need to stop pullbacks from being applied, how do I do that?
What we actually need is to be able to replace J uhat/det J with J P(uhat)/det J where P is the projection onto a vector in local coordinates, probably it is easiest to just hand write the whole expression as it should appear after pullback.
I you write everything after pullback, you might be able to get away with editing tsfc/ufl_utils.py to set do_apply_function_pullbacks=False. I'm hope that will work.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake from firedrake import * from ufl.geometry import Jacobian, JacobianDeterminant
# Set up mesh = PeriodicSquareMesh(120, 120, 2, quadrilateral=True) dt = 0.0001 outfile = File("out.pvd")
V1_elt = FiniteElement("RTCF", quadrilateral, 2) V = FunctionSpace(mesh, V1_elt) Q = FunctionSpace(mesh, "DG", 1) W = MixedFunctionSpace((V, Q))
# Initialise U0 = Function(W) u0, h0 = U0.split() f = Constant(100.) x = SpatialCoordinate(mesh) a = Constant(0.1) delta_x = Constant(2/120.) h0.interpolate(cos(pi*x[1]/delta_x)*cos(pi*x[0]/delta_x)*exp(-((x[0]-1)**2 + (x[1]-1)**2)/a**2)) w, phi = TestFunctions(W) u, h = TrialFunctions(W) U1 = Function(W) U1.assign(U0) u1, h1 = U1.split() outfile.write(u1, h1)
J = Jacobian(mesh) detJ = JacobianDeterminant(mesh) detJ_inv = (1./detJ)
aeqn = ( inner(dot(J, w), dot(J, u))*detJ_inv*dx + 0.5*dt*f*inner(w, perp(u))*dx + phi*h*detJ*dx - 0.5*dt*h*div(w)*dx + 0.5*dt*phi*div(u)*dx )
Leqn = ( inner(dot(J, w), dot(J, u0))*detJ_inv*dx - 0.5*dt*f*inner(w, perp(u0))*dx + phi*h0*detJ*dx + 0.5*dt*div(w)*h0*dx - 0.5*dt*phi*div(u0)*dx )
prob = LinearVariationalProblem(aeqn, Leqn, U1) solver = LinearVariationalSolver(prob) solver.solve() u1, h1 = U1.split() outfile.write(u1, h1) _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake