JIT-expressions 1) used in automated differentiation 2) using custom functions

Question

I made a modification to the nonlinear Poisson demo in the code below so that the nonlinearity comes from an Expression that depends on the Function u_k, rather than from UFL friendly operators like exp, sin, etc. I have two questions regarding this:

1) Shouldn't it be possible to somehow mark that the expression q below depends on u_k, so that if one tries to use dolfin.derivative(F, u_k, du), it automatically takes this into account when performing automated differentiation. The derivative of the expression should also be provided in that case of course.

2) Using expressions this way is slow, since it uses Python. I would like to do something like

from external_module import q

and then in the expression, use it as follows:

def eval(self, value, x):
    value[0] = q(u_k(x))

but since q is a compiled function (using Cython), I somehow feel like I should use a just-in-time compiled expression for this. Is this possible (and how, since there are not that many examples of this and C++ isn't exactly my strong suit)?

Many thanks

The modified nonlinear Poisson code:

from dolfin import *

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)

V = FunctionSpace(mesh, "CG", 1)

# Define boundary condition
g = Constant(1.0)
bc = DirichletBC(V, g, DirichletBoundary())

# Define variational problem
u_k = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])")

# define F and J
method = 2
if method == 1: # as in the demo-code
    F = inner((1 + u_k**2)*grad(u_k), grad(v))*dx - f*v*dx
    J = derivative(F, u_k, du)

elif method == 2: # using custom expressions
    class CustomNonlinearity(Expression):
        def eval(self, value, x):
            value[0] = 1 + u_k(x)**2
    q = CustomNonlinearity()

    class DCustomNonlinearity(Expression):
        def eval(self, value, x):
            value[0] = 2 * u_k(x)
    Dq = DCustomNonlinearity()

    F = inner(q*grad(u_k), grad(v))*dx - f*v*dx
    J = inner(q*grad(du), grad(v)) * dx + inner(Dq * du * grad(u_k), grad(v)) * dx

# Compute solution
problem = NonlinearVariationalProblem(F, u_k, bc, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

# Plot solution and solution gradient
plot(u_k, title="Solution")
plot(grad(u_k), title="Solution gradient")
interactive()

Answer 1

You need to look at dolfin-adjoint for such tools:
http://dolfin-adjoint-doc.readthedocs.org/en/latest/

Comments

From a glance at the documentation, it seems that it might help me with question 1.

For question 2, I am looking for something like this:
http://fenicsproject.org/qa/3498/use-bessel-functions-in-jit-compiled-expressions

For example:

from dolfin import *
code = '''
#include "my_functions.h"

namespace dolfin {
    class MyFun : public Expression
    {
        public:
            MyFun(): Expression() {};
        void eval(Array<double>& values, const Array<double>& x) const {
            argument = u_k(x[0])
            values[0] = fd(argument);
        }
    };
}'''

f=Expression(code)

For this, I would need to

1. Find a way to tell the compiler it has to look for a function fd in a module my_functions.so

2. Tell the Expression class that it takes another argument, i.e. the function u_k (preferably its underlying C++ object)

It would be nice to know if this is possible.

---

I posted a new question today, focusing on the second part of this question.

---

Sure, it is possible. You may use either the
compile_extension_module function in FEniCS or
Instant directly.
See more about what happens behind the curtain eg here:
http://hplgit.github.io/fenics-mixed/doc/pub/sphinx-basicstrap/._part0003_fenics-mixed.html