UFL

dot does not check the argument shape dimensions

Bug #1155210 reported by Andrew McRae
6
This bug affects 1 person
Affects Status Importance Assigned to Milestone
UFL
Fix Released
Critical
Martin Sandve Alnæs

Bug Description

When attempting to call as_vector() with too few arguments (in this case, 1 instead of 2), the program runs if we only solving for a single field, but breaks when trying to solve for multiple functions simultaneously. Code below.

=====================================================
from dolfin import *
mesh = UnitSquareMesh(8,8)
S = FunctionSpace(mesh, "BDM", 1)
V = FunctionSpace(mesh, "DG", 0)
W = S*V

u0 = interpolate(Constant((-1.0,1.0)),S)
h0 = interpolate(Constant("1.0"),V)

u = TrialFunction(S)
v = TestFunction(S)
a = dot(as_vector([-u[1]]),v)*dx # works
#a = dot(as_vector([-u[1],0.0]),v)*dx # also works
L = dot(u0,v)*dx
u = Function(S)
solve(a==L,u)
print "First part compiles"
plot(u)
interactive()

u,h = TrialFunctions(W)
v,phi = TestFunctions(W)
a = dot(as_vector([-u[1]]),v)*dx + h*phi*dx # fails
#a = dot(as_vector([-u[1],0.0]),v)*dx + h*phi*dx # works
L = dot(u0,v)*dx + h0*phi*dx
w = Function(W)
solve(a==L,w)
u,h = w.split(deepcopy=True)
print "Second part compiles"
plot(u)
interactive()
=====================================================

atm112@ae-amcrae:~/Dropbox/FEniCS/scratch$ python asvector2.py
Solving linear variational problem.
First part compiles
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "asvector2.py", line 27, in <module>
    solve(a==L,w)
  File "/home/atm112/local/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/fem/solving.py", line 268, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/atm112/local/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/fem/solving.py", line 292, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/home/atm112/local/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/fem/solving.py", line 80, in __init__
    a = Form(a, form_compiler_parameters=form_compiler_parameters)
  File "/home/atm112/local/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/fem/form.py", line 56, in __init__
    common_cell)
  File "/home/atm112/local/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/compilemodules/jit.py", line 66, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/atm112/local/fenics/src/dolfin/current/local/lib/python2.7/site-packages/dolfin/compilemodules/jit.py", line 154, in jit
    return jit_compile(form, parameters=p, common_cell=common_cell)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/jitcompiler.py", line 77, in jit
    return jit_form(ufl_object, parameters, common_cell)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/jitcompiler.py", line 197, in jit_form
    parameters=parameters)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/compiler.py", line 156, in compile_form
    ir = compute_ir(analysis, parameters)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/representation.py", line 96, in compute_ir
    for (i, fd) in enumerate(form_datas)]
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/representation.py", line 216, in _compute_integral_ir
    parameters)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/quadrature/quadraturerepresentation.py", line 87, in compute_integral_ir
    itg_data.domain_type, form_data.cell)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/quadrature/quadraturerepresentation.py", line 145, in _transform_integrals_by_type
    terms = _transform_integrals(transformer, integrals_dict, domain_type)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/quadrature/quadraturerepresentation.py", line 363, in _transform_integrals
    terms = transformer.generate_terms(integral.integrand(), domain_type)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/quadrature/quadraturetransformerbase.py", line 739, in generate_terms
    terms = self.visit(integrand)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ufl/algorithms/transformer.py", line 101, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ufl/algorithms/transformer.py", line 105, in visit
    r = h(o)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/quadrature/quadraturetransformerbase.py", line 549, in index_sum
    ops.append(self.visit(summand))
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ufl/algorithms/transformer.py", line 101, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ufl/algorithms/transformer.py", line 105, in visit
    r = h(o)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/quadrature/quadraturetransformerbase.py", line 512, in indexed
    code = self.visit(indexed_expr)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ufl/algorithms/transformer.py", line 105, in visit
    r = h(o)
  File "/home/atm112/local/fenics/lib/python2.7/site-packages/ffc/quadrature/quadraturetransformerbase.py", line 714, in list_tensor
    op = o.operands()[c0]
IndexError: tuple index out of range

Revision history for this message
Andrew McRae (andymc) wrote :

Aside: the output seems to be crazy, at least compared to the related program where we move the rotation onto the RHS. Not sure if bug or genuine...

=============================
from dolfin import *
mesh = UnitSquareMesh(8,8)
S = FunctionSpace(mesh, "BDM", 1)
V = FunctionSpace(mesh, "DG", 0)
W = S*V

u0 = interpolate(Constant((-1.0,1.0)),S)
h0 = interpolate(Constant("1.0"),V)

u = TrialFunction(S)
v = TestFunction(S)
a = dot(u,v)*dx
#L = dot(as_vector([-u0[1]]),v)*dx # works
L = dot(as_vector([-u0[1],0.0]),v)*dx # also works
u = Function(S)
solve(a==L,u)
print "First part compiles"
plot(u)
interactive()

u,h = TrialFunctions(W)
v,phi = TestFunctions(W)
a = dot(u,v)*dx + h*phi*dx
#L = dot(as_vector([-u0[1]]),v)*dx + h0*phi*dx # fails
L = dot(as_vector([-u0[1],0.0]),v)*dx + h0*phi*dx # works
w = Function(W)
solve(a==L,w)
u,h = w.split(deepcopy=True)
print "Second part compiles"
plot(u)
interactive()

Revision history for this message
Martin Sandve Alnæs (martinal) wrote :

The bug is in dot, not checking for wrong argument dimensions. Fix on the way.

Changed in ufl:
importance: Undecided → Critical
no longer affects: ffc
Changed in ufl:
status: New → In Progress
assignee: nobody → Martin Sandve Alnæs (martinal)
milestone: none → 1.2.0
summary: - inconsistent behaviour of as_vector when too few arguments provided
+ dot does not check the argument shape dimensions
Changed in ufl:
status: In Progress → Fix Committed
Johannes Ring (johannr)
Changed in ufl:
status: Fix Committed → Fix Released
To post a comment you must log in.
This report contains Public information  
Everyone can see this information.

Other bug subscribers

Remote bug watches

Bug watches keep track of this bug in other bug trackers.