Topology optimization example (http://www.dolfin-adjoint.org/en/latest/documentation/poisson-topology/poisson-topology.html#poisson-topology-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):
In case I'm using interpolate inside forward problem it says
libadjoint.exceptions.LibadjointErrorHashFailed: Hash failed:
Control:0:0:Forward has no entry.
and in case of project it says
return ("Mesh", renumbering[self], self._ufl_coordinate_element)
KeyError: Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1),
dim=2), 4)
I've tired to ask on dolphin-adjoint issues page, but got no answers, so I'm asking here. Would be glad to receive any help.
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
