Add Primitive operator

[Chaffra Affouda]

How hard would it be to add a primitive operator to UFL. I am thinking something like

#new integration.py file
#!python

@ufl_type(is_abstract=True, num_ops=1, inherit_indices_from_operand=0, is_scalar=True)
class Primitive(Operator):
    __slots__ = ("_name",)
    def __init__(self,f, **kwargs):
        #dolfin.cpp.Expression.__init__(self)
        self._name = 'Primitive'
        Operator.__init__(self,operands=(f,))


    def __new__(cls, f):
        if isinstance(f, (ScalarValue, Zero)):
            return FloatValue(0.0)
        return Operator.__new__(cls)

    def reconstruct(self, f):
        "Return a new object of the same type with new operands."
        return self._ufl_class_(f)

    def __str__(self):
        return "%s(%s)" % (self._name, self.ufl_operands[0])

    def __repr__(self):
        return "%s(%r)" % (self._name, self.ufl_operands[0])

you can then create the appropriate handle in class GenericDerivativeRuleset

def primitive(self, o, df):
        return o.ufl_operands[0]*df

I can then create a primitive special function in dolfin like

class PrimitiveFunction(Primitive, dolfin.cpp.Expression):
    #__slots__ = ("_name","_f","_constant","_mesh","_subdomains","_F")
    def __init__(self,f, **kwargs):
        dolfin.cpp.Expression.__init__(self)
        self._f = f
        mesh = f.function_space().mesh()
        self._mesh = mesh
        subdomains = CellFunction('size_t',mesh)
        self._subdomains = subdomains
        dx = Measure('dx',domain=mesh)[subdomains]
        self._F = Form(f*dx(1))

        Primitive.__init__(self,f)



    def evaluate(self, x, mapping, component, index_values, derivatives=()):
        return self.ufl_evaluate(x, component, derivatives)

    def eval(self, values, x, component=0):

        inside_function = lambda y: y[component]<=x[component]
        domain = AutoSubDomain(inside_function=inside_function)
        self._subdomains.set_all(0)
        domain.mark(self._subdomains,1)
        F = self._F
        values[0] = assemble(F)

    def eval_cell(self, values, x, cell):
        self.eval(values,x)

    def ufl_evaluate(self,x, component=(0,), derivatives=()):
        values = np.zeros(1,dtype='d')
        self.eval(values,x, component=component[0])
        return values[0]

Where I am stuck is how do I translate the Primitive operator evaluation in FFC. I created a new function in QuadratureTransformerBase but I don't know what to add here for the JIT compiler to return the integral of the operand at a given point?? The problem is that the evaluation of the primitive relies heavily on dolfin.

def primitive(self,o, *operands):
        print("\n\nVisiting Primitive: " + repr(o))
        #What else?

[Martin Sandve Alnæs]

You don't need ufl_evaluate. That's for symbolic evaluation and is mainly used for testing and debugging.

The jit compiler returns code that is evaluated inside the UFC integral code that FFC generates. That code (currently) only has access to the local degrees of freedom for each ufl.Coefficient/dolfin.Function in the form.

We've been thinking about adding a callback mechanism so the form can call your primitive function (or anything that can be evaluated really) for each quadrature point, however it's not high on anybodys todo list at the moment. It's a rather crosscutting feature that requires changes to all FEniCS components, so it must be carefully designed to fit with other work in progress.

If I understand you correctly here, you want to define an operator F(f, y), where f is a Function, and express derivative(F, f) in UFL? I think it you can do that within the current limitations. I'm posting the solution to the QA forum instead because others may be interested.

[Martin Sandve Alnæs]

The QA forum is really slow, I'm pasting the code here instead.

@chaffra do you have a good example of what do use this for? I'd like to refine this into a dolfin demo.

"""Demo of a global operator with FEniCS.

Given

  f = f(x)
  I(f) = f
  \Omega(y) = circle around y

Sets up the functional

  F = F(f, y) = \int_{\Omega(y)} I(f) dx

and computes its derivative

  dF/df = \int_{\Omega(y)} dI/df dx

with automatic differentiation of the integrand I(f).
Next, sets up the functional

  M = \int_\Omega F^2 dx

and computes the Gateaux derivative

  L_1 = d/ds M(f + s*v)

with automatic differentiation by specifying dF/df manually.

Finally compares L_1 with the manually derived

  L_2 = \int_\Omega 2 F dF/df v dx

to show that they are the same.

Try changing f(x), I(f), \Omega(y), M(F) to see how the automatic differentiation works.
"""


from dolfin import *


class FunctionalExpression(Expression):
    """Computes a functional for each eval(x) after calling an update function.

    This allows the update function to update mesh, subdomains, and functions
    that the functional depends on.

    In the example below, the integration subdomain is updated depending on x.
    """

    def init(self, functional, update_functional):
        self._update_functional = update_functional
        self._functional = Form(functional)

    def eval(self, values, x):
        self._update_functional(x)
        values[0] = assemble(self._functional)


# Setup discretization. Note that the number of dofs is N = O(n^d),
# and the time to assemble the form involving the global operator
# is on the order O(N^2) = O(n^(2d)).
n = 16
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
f = Function(V)

f.interpolate(Expression("x[0] < .5 ? 1.0: 0.0")) # TRYME


# Define form for F(f)
# NB! Note the dependency on the cell_markers object here,
# which is updated in update_functional
cell_markers = CellFunction("size_t", mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)

F_form = f*dx(1) # TRYME


def update_functional(x):
    """Update the global cell_markers object to mark a circle around x with 1.

    Note that this global object is referenced by the form for F and used when re-assembling F.
    """
    global cell_markers

    inside = lambda y: sum((x-y)**2) < 0.5**2 # TRYME

    cell_markers.set_all(0)
    domain = AutoSubDomain(inside_function=inside)
    domain.mark(cell_markers, 1)


# F = F(y, f) = f(x)*dx(1) where dx(1) depends on y and is defined by mark_subdomain
F = FunctionalExpression(element=f.element())
F.init(F_form, update_functional)

# Compute partial derivative of F w.r.t. f.
# (NB! Note the difference between diff(f, f) == 1 and derivative(f, f, v) == v)
dFdf_form = diff(F_form, f)
dFdf = FunctionalExpression(element=f.element())
dFdf.init(dFdf_form, update_functional)


# Define a functional
M = F**2*dx # TRYME

# Apply semi-automatic differentiation, by specifying the
# partial derivative of F w.r.t. f as another function dFdf:
L1 = derivative(M, f, v, coefficient_derivatives={F: dFdf})

# Manually derived derivative of F**2*dx
L2 = 2*F*dFdf*v*dx # TRYME: If you change M, this will of course be wrong


# Assemble and compare
print "Assembling b1 (this may take a while)"
b1 = assemble(L1)
print "Assembling b2 (this may take a while)"
b2 = assemble(L2)
bdiff = Vector(b2)
bdiff.axpy(-1.0, b1)
print "|b2-b1|_rel:", norm(bdiff) / max(norm(b1), norm(b2))


# Interpolate Integral objects into functions for postprocessing
print "Interpolating F (this may take a while)"
F = interpolate(F, V)
print "Interpolating dFdf (this may take a while)"
dFdf = interpolate(dFdf, V)

# Write and plot functions
File("f.xdmf") << f
File("dFdf.xdmf") << dFdf
File("F.xdmf") << F
plot(f, title="f")
plot(dFdf, mesh=mesh, title="dFdf")
plot(F, mesh=mesh, title="F")
interactive()

[Martin Sandve Alnæs]

I'm closing this issue, opened a dolfin demo issue instead

https://bitbucket.org/fenics-project/dolfin/issue/413/add-a-demo-of-global-operators-with-semi

[Chaffra Affouda]

Yes the example I was going for is detailed on the QA website

http://fenicsproject.org/qa/5920/general-strategy-for-solving-energy-minimization-problem

I think what you wrote applies to what I am trying to acheive. I will give it a try and post my findings on the qa. Thanks a lot!