JIT dying when using project or interpolate inside forward problem

Topology optimization example is working fine, but then I try to modify it by using project or interpolate I'm getting JIT errors during optimization phase. The following code provides minimal not working example (it is based to topology optimization code from the examples): It is also, along with error responses, attached to this message.

from __future__ import print_function
from dolfin import *
from dolfin_adjoint import *
import numpy as np
import pickle
from collections import defaultdict
import sys
import pyipopt

parameters["std_out_all_processes"] = False

V = Constant(0.4)      
p = Constant(5)        

eps = Constant(1.0e-3) 
alpha = Constant(1.0e-8) 

def k(a):
    return eps + (1 - eps) * a**p

n = 250
mesh  = UnitSquareMesh(n, n)
mesh2 = UnitSquareMesh(2*n, 2*n)

A = FunctionSpace(mesh, "CG", 1)
A2 = FunctionSpace(mesh2, "CG", 1)

class WestNorth(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] == 0.0 or x[1] == 1.0) and on_boundary

bc = [DirichletBC(A2, 0.0, WestNorth())]
f = interpolate(Constant(1.0e-2), A2, name="SourceTerm")

def forward(a):
    T = Function(A2, name="Temperature")
    v = TestFunction(A2)
    #aa = interpolate(a, A2) 
    aa = project(a, A2) 
    F = inner(grad(v), k(aa)*grad(T))*dx - f*v*dx
    solve(F == 0, T, bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7,
                                                              "maximum_iterations": 20}})
    return T

if __name__ == "__main__":
    a = interpolate(V, A, name="Control")
    T = forward(a)                  

    J = Functional(f*T*dx)
    m = Control(a)
    Jhat = ReducedFunctional(J, m)

    lb = 0.0
    ub = 1.0

    class VolumeConstraint(InequalityConstraint):
        def __init__(self, V):
            self.V  = float(V)
            self.smass  = assemble(TestFunction(A) * Constant(1) * dx)
            self.tmpvec = Function(A)

        def function(self, m):
            reduced_functional_numpy.set_local(self.tmpvec, m)
            integral = self.smass.inner(self.tmpvec.vector())
            if MPI.rank(mpi_comm_world()) == 0:
                print("Current control integral: ", integral)
            return [self.V - integral]

        def jacobian(self, m):
            return [-self.smass]

        def output_workspace(self):
            return [0.0]

        def length(self):
            return 1

    problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V))
    parameters = {"acceptable_tol": 1.0e-3, "maximum_iterations": 100, "print_level":0}
    solver = IPOPTSolver(problem, parameters=parameters)
    a_opt = solver.solve()
    File("output/a_opt.pvd") << a_opt

Comments (1)