ufl team mailing list archive

[Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx>]

On 8 August 2011 15:19, Martin Sandve Alnæs <martinal@xxxxxxxxx> wrote:
> Thanks Kristian.
>
> In the abs case, would this change be sufficient?
>
> @@ -388,7 +388,8 @@
>     def abs(self, o, a):
>         f, fprime = a
>         o = self.reuse_if_possible(o, f)
> -        oprime = conditional(eq(f, 0), 0, Product(sign(f), fprime))
> +        oprime = sign(f)*fprime
>         return (o, oprime)
>
> sign returns 0 if f is 0:
>
> def sign(x):
>    "UFL operator: Take the sign (+1 or -1) of x."
>    return conditional(eq(x, 0), 0, conditional(lt(x, 0), -1, +1))
>
> In the generic conditional case, I'm not sure if this is
> easily solvable in the current UFL implementation.

For the given problem the above change should do the trick.
If other combinations of conditional turns out to crash compilation we
might need to look at other ways of solving the issue as you point
out.
It's difficult to predict so let's wait and see if some one hits a
corner case we can investigate.

Kristian

> Martin
>
> On 6 August 2011 17:00, Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx> wrote:
>> On 4 August 2011 14:56, Chaffra <question166940@xxxxxxxxxxxxxxxxxxxxx> wrote:
>>> New question #166940 on DOLFIN:
>>> https://answers.launchpad.net/dolfin/+question/166940
>>>
>>> Dear all,
>>>
>>> I want to be able to solve poisson's equation with some abs(absolute value) functions in it, but ffc breaks with
>>>
>>> Exception: False value of Condtional should only be one function: {((1, 'k', 3, 3),): 'FE0[ip][k]*C[1]'}
>>
>> I tried the following simplified version of your code as a pure UFL
>> file and compiled it with FFC which results in a different, but
>> related crash (probably because I'm on the dev version of FEniCS):
>>
>> V  = FiniteElement("Lagrange", triangle, 1)
>> u  = Coefficient(V)
>> du = TrialFunction(V)
>> v  = TestFunction(V)
>> L1 = abs(u)*v*dx
>> J = derivative(L1,u,du)
>>
>> Found Argument in Conditional(EQ(Coefficient(FiniteElement('Lagrange',
>> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
>> (), {}), Product(Argument(FiniteElement('Lagrange', Cell('triangle',
>> Space(2)), 1, None), 1),
>> Conditional(EQ(Coefficient(FiniteElement('Lagrange', Cell('triangle',
>> Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
>> Conditional(LT(Coefficient(FiniteElement('Lagrange', Cell('triangle',
>> Space(2)), 1, None), 0), Zero((), (), {})), IntValue(-1, (), (), {}),
>> IntValue(1, (), (), {}))))), this is an invalid expression.
>> Traceback (most recent call last):
>>  File "/home/oelgaard/software/fenics/bin/ffc", line 195, in <module>
>>    sys.exit(main(sys.argv[1:]))
>>  File "/home/oelgaard/software/fenics/bin/ffc", line 176, in main
>>    compile_form(ufd.forms, ufd.object_names, prefix, parameters)
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ffc/compiler.py",
>> line 150, in compile_form
>>    analysis = analyze_forms(forms, object_names, parameters, common_cell)
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ffc/analysis.py",
>> line 64, in analyze_forms
>>    common_cell) for form in forms)
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ffc/analysis.py",
>> line 64, in <genexpr>
>>    common_cell) for form in forms)
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ffc/analysis.py",
>> line 151, in _analyze_form
>>    ffc_assert(len(compute_form_arities(preprocessed_form)) == 1,
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/formtransformations.py",
>> line 332, in compute_form_arities
>>    parts = compute_form_with_arity(form, arity)
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/formtransformations.py",
>> line 317, in compute_form_with_arity
>>    res = transform_integrands(form, _transform)
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/transformations.py",
>> line 799, in transform_integrands
>>    integrand = transform(itg.integrand())
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/formtransformations.py",
>> line 313, in _transform
>>    e, provides = pe.visit(e)
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/transformations.py",
>> line 144, in visit
>>    r = h(o, *map(self.visit, o.operands()))
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/transformations.py",
>> line 148, in visit
>>    r = h(o)
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/algorithms/formtransformations.py",
>> line 68, in expr
>>    error("Found Argument in %s, this is an invalid expression." % repr(x))
>>  File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/ufl/log.py",
>> line 148, in error
>>    raise self._exception_type(self._format_raw(*message))
>> ufl.log.UFLException: Found Argument in
>> Conditional(EQ(Coefficient(FiniteElement('Lagrange', Cell('triangle',
>> Space(2)), 1, None), 0), Zero((), (), {})), Zero((), (), {}),
>> Product(Argument(FiniteElement('Lagrange', Cell('triangle', Space(2)),
>> 1, None), 1), Conditional(EQ(Coefficient(FiniteElement('Lagrange',
>> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})), Zero((),
>> (), {}), Conditional(LT(Coefficient(FiniteElement('Lagrange',
>> Cell('triangle', Space(2)), 1, None), 0), Zero((), (), {})),
>> IntValue(-1, (), (), {}), IntValue(1, (), (), {}))))), this is an
>> invalid expression.
>>
>> The problem is the UFL derivative of abs(u) which can be inspected by:
>>
>> from ufl.algorithms import expand_derivatives
>> V  = FiniteElement("Lagrange", triangle, 1)
>> u  = Coefficient(V)
>> du = TrialFunction(V)
>> v  = TestFunction(V)
>> L1 = abs(u)*v*dx
>> J1 = derivative(L1,u,du)
>> print expand_derivatives(J1)
>>
>> { v_{-2} * (((w_0) == (0)) ? (0) : (v_{-1} * (((w_0) == (0)) ? (0) :
>> (((w_0) < (0)) ? (-1) : (1))))) } * dx0
>>
>> Here, 'v_{-2}' is 'v', 'w_0' is 'u' and 'v_{-1}' is 'du'.
>>
>> This is identical to:
>>
>> c = conditional( eq(u, 0), 0, du*conditional(eq(u, 0), 0,
>> conditional(lt(u,0),-1,1) ) )
>> J = v*c*dx
>> print J
>>
>> Notice that 'du' is present inside the first conditional (the second
>> return value) which triggers the error.
>>
>> I'm not sure that it makes sense to allow 'du' in the condition which
>> is being tested e.g.:
>>
>> c = conditional( eq(du, 0), 0, 1 )
>>
>> but there should be no problems with allowing it in the return values like:
>>
>> c = conditional( eq(u, 0), du, 2*du )
>>
>> However, with the particular form in question, the return value of 'c'
>> will be either 0, du or -du, so J in the case of u == 0 will be a
>> linear form!!!
>>
>> Things would work fine if UFL could move 'du' outside the conditional
>> such that it is equal to:
>>
>> c = du*conditional( eq(u, 0), 0, conditional(eq(u, 0), 0,
>> conditional(lt(u,0),-1,1) ) )
>> J = v*c*dx
>>
>> which compiles fine and looks correct.
>>
>> Kristian
>>
>>> Any idea?
>>>
>>> Here's an example based on demo_poisson.py:
>>>
>>> from dolfin import *
>>>
>>> # Create mesh and define function space
>>> mesh = UnitSquare(32, 32)
>>> V = FunctionSpace(mesh, "Lagrange", 1)
>>>
>>> # Define Dirichlet boundary (x = 0 or x = 1)
>>> def boundary(x):
>>>    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
>>>
>>> # Define boundary condition
>>> u0 = Constant(0.0)
>>> bc = DirichletBC(V, u0, boundary)
>>>
>>> # Define variational problem
>>> u = Function(V)
>>> du = TrialFunction(V)
>>>
>>> v = TestFunction(V)
>>> f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
>>> g = Expression("sin(5*x[0])")
>>> a = inner(grad(u), grad(v))*dx
>>> L = f*v*dx + g*v*ds + abs(u)*v*dx
>>> F = a-L
>>> J = derivative(F,u,du)
>>>
>>> # Compute solution
>>> problem = VariationalProblem(F, J, bc)
>>> u = problem.solve()
>>>
>>> # Save solution in VTK format
>>> file = File("poisson.pvd")
>>> file << u
>>>
>>> # Plot solution
>>> plot(u, interactive=True)
>>>
>>>
>>>
>>> --
>>> You received this question notification because you are a member of
>>> DOLFIN Team, which is an answer contact for DOLFIN.
>>>
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~ufl
>> Post to     : ufl@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~ufl
>> More help   : https://help.launchpad.net/ListHelp
>>
>