# Dijitso JIT compilation failed on assembly of subdomain problem

Hi,

I am quite new to Fenics and dolfin-adjoint, and I’ve been trying to solve a 3D- heat transfer topology optimization problem.
When defining an assembly for minimization, I get an error:

`DijitsoError: Dijitso JIT compilation failed, see '/..../box/jitfailure-ffc_form_26b5aad55790b45b60fe7a7a80201fa4da1014b6' for details`

bellow is a minimum working example:

``````from fenics import *
import numpy as np
from mpi4py import MPI

rank = MPI.COMM_WORLD.rank

P1 = Point(0, 0, 0)
P2 = Point(10, 10, 10)
mesh = BoxMesh(P1, P2, 10, 10, 10)

V = Constant(0.4)  # volume bound on the control
p = Constant(5)  # power used in the solid isotropic material
eps = Constant(1.0e-3)  # epsilon used in the solid isotropic material
alpha = Constant(1.0e-8)  # regularisation coefficient in functional
A = FunctionSpace(mesh, "CG", 1)  # function space for control
P = FunctionSpace(mesh, "CG", 1)  # function space for solution

# SIMP rule
def k(rho):
return eps + (1 - eps) * rho **p

# Mark subdomains (baseplate, heatplate, and design space in between)
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
subdomains.set_all(0)
class BasePlate(SubDomain):
def inside(self, x, on_boundary):
return x[2] < 2

class HeaterPlate(SubDomain):
def inside(self, x, on_boundary):
return x[2] > 8

BasePlate().mark(subdomains, 1)  # baseplate mark
HeaterPlate().mark(subdomains, 2)  # heater mark
#File('subdomains.pvd') << subdomains

submesh = SubMesh(mesh, subdomains, 0)

# Define a measure for optimization subdomain
dx = Measure('dx', domain=mesh, subdomain_data= subdomains)

Asub = FunctionSpace(submesh, "CG", 1)
rhosub = interpolate(Constant(1.0), Asub)

# Boundary condition at baseplate
class BPBottom(SubDomain):
def inside(self, x, on_boundary):
return x[2] < DOLFIN_EPS and on_boundary
bc = [DirichletBC(P, 20, BPBottom())]
f = interpolate(Constant(1.0e-2), P)  # the volume source term for the PDE

def forward(rho,rho_sub):
"""Solve the forward problem for a given material distribution rho(x)."""
T = Function(P, name="Temperature")
v = TestFunction(P)

F = (
- f * v * dx(2)                                # heater source term
)
solve(F == 0, T, bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7,
"maximum_iterations": 20}})

return T

rho_sub = interpolate(Constant(0.4), Asub)
rho = interpolate(V,A)
T = forward(rho,rho_sub)
file = File("solution.pvd")
file << T

## error happens here!!

J = assemble(f * T * dx(2) + alpha * inner(grad(rho_sub), grad(rho_sub)) * dx(0))
m = Control(rhosub)

Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)
``````

I’ve tried including all subdomains in the equation but with no luck.
Perhaps someone could also suggest a better way to define control on a subdomain?

Installed on Ubuntu 20.04.2 LTS with:

inside a conda enviroment.

I can reproduce this issue with dolfin and dolfin-adjoint v 2019.1.0, but not with the development version. Please upgrade dolfin to the latest development versions.

Thank you for the quick reply.

I still get the same error with the version installed the development version with:

`````` pip3 install git+https://github.com/dolfin-adjoint/pyadjoint.git@master
``````

After a fresh installation of Fenics on a new environment.
Isn’t that the development version of dolfin-adjoint?

On a slight tangent. Does dolfin-adjoint work with dolfinx?

It is not the development version of dolfin-adjont that is required, but the development version of dolfin (i.e. version 2019.2.0.dev0).

Dolfin-adjoint does not currently support dolfin-x, as the API has changed drastically. A long term goal is to extend dolfin-adjoint to dolfin-x, but the API has to stabilize (a stable release has to be made before such work can start).

I tried the development version and it works there. Thanks for the tip.

Hi @dokken

I am still having issues while assembling subdomains. More specifically, in assembling a TestFunction of a submesh.

Here’s a MWE, where the VolumeConstraint class fails at __init__():

``````from fenics import *
import numpy as np
from mpi4py import MPI

rank = MPI.COMM_WORLD.rank
parameters["std_out_all_processes"] = False

class Rho(UserExpression):
def __init__(self, subdomains, **kwargs):
super(Rho, self).__init__()
self.subdomains = subdomains

def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] == 3:
values[0] = 0.4
else:
values[0] = 1

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

class BasePlate(SubDomain):
def inside(self, x, on_boundary):
return x[2] < 2

class HeaterPlate(SubDomain):
def inside(self, x, on_boundary):
return x[2] > 8

# Boundary condition at baseplate
class BPBottom(SubDomain):
def inside(self, x, on_boundary):
return x[2] < DOLFIN_EPS and on_boundary

def forward(rho, rho_sub):
"""Solve the forward problem for a given material distribution rho(x)."""
T = Function(P, name="Temperature")
v = TestFunction(P)
bc = [DirichletBC(P, 20, BPBottom())]
f = interpolate(Constant(1.0e-2), P)  # the volume source term for the PDE

F = (
- f * v * dx(2)  # heater source term
)
solve(F == 0, T, bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7,
"maximum_iterations": 20}})
return T

def eval_cb(j, a):
a_viz.assign(a)
controls << a_viz

class VolumeConstraint(InequalityConstraint):
"""A class that enforces the volume constraint g(a) = V - a*dx >= 0."""
def __init__(self, V):
self.V = V
# problem happens here !!
self.smass = assemble(TestFunction(ASub) * Constant(1) * dx(3))
self.tmpvec = Function(ASub)
def function(self, m):
set_local(self.tmpvec, m)
integral = self.smass.inner(self.tmpvec.vector())

if MPI.COMM_WORLD.rank == 0:
print("Current control integral: ", integral)
return [self.V - assemble(rho_sub*dx(3))]

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

def output_workspace(self):
return [0.0]

def length(self):
"""Return the number of components in the constraint vector (here, 3)."""
return 1

P1 = Point(0, 0, 0)
P2 = Point(10, 10, 10)
mesh = BoxMesh(P1, P2, 10, 10, 10)

p = Constant(5)  # power used in the solid isotropic material
eps = Constant(1.0e-3)  # epsilon used in the solid isotropic material
alpha = Constant(1.0e-8)  # regularisation coefficient in functional
A = FunctionSpace(mesh, "CG", 1)  # function space for control
P = FunctionSpace(mesh, "CG", 1)  # function space for solution

subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
subdomains.set_all(3)            # design domain 3
BasePlate().mark(subdomains, 1)  # baseplate  1
HeaterPlate().mark(subdomains, 2)  # heater  2
# File('subdomains.pvd') << subdomains

submesh = SubMesh(mesh, subdomains, 3)
ASub = FunctionSpace(submesh, "CG", 1)
rho_sub = interpolate(Rho(subdomains), ASub)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
rho = interpolate(Rho(subdomains), A)

# rho_file = File('rho.pvd')
# rho_file << rho

bc = [DirichletBC(P, 20, BPBottom())]
f = interpolate(Constant(1.0e-2), P)  # the volume source term for the PDE
T = forward(rho, rho_sub)
file = File("solution.pvd")
file << T

#
controls = File("output/control_iterations.pvd")
a_viz = Function(A, name="ControlVisualisation")

J = assemble(f * T * dx(2)
)
m = Control(rho_sub)
Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)

lb = 0.0
ub = 1.0
V3 = assemble(rho_sub*dx(3))  # volume bound
parameters = {"acceptable_tol": 1.0e-3, "maximum_iterations": 100}
problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V3))
solver = IPOPTSolver(problem, parameters=parameters)
a_opt = solver.solve()

#
# File("output/final_solution.pvd") << a_opt
# xdmf_filename = XDMFFile(MPI.comm_world, "output/final_solution.xdmf")
# xdmf_filename.write(a_opt)

``````

Is there a better way to apply control over a single subdomain?

Thanks again for your help. I really appreciate the work you’re all doing here.

You are mixing the integration measures of the `SubMesh` and your mesh. Changing:

to

``````        self.smass = assemble(TestFunction(ASub) * Constant(1) * Measure("dx"))
``````

which will progress the solver further, but will result in another error, as you are mixing sub-meshes and meshes all over the place. For me to be of any more help, you need to supply a minimal example with a full description of what you would like to do. The code above is far from minimal and is lacking a mathematical description.

Hi,

Sorry for the trivial mistakes, I am still quite new to fenics.

The problem I am trying to solve is a 3D heat-transfer topology optimization problem based on poisson-topology.

The design domain sits between a volume heat-source on top and a volume heat sink at the bottom. Here marker 1 is the heat sink, marker 2 is the volume heat source, and marker 3 is the design domain to be optimized:

The objective would be either to minimize the thermal compliance of the design domain or to minimize the temperature gradient in the volume source. Subject to a volume constraint over the design domain.

I am trying to only optimize the design domain, without changing the design of the heater plate or the heat sink
To do that, I tried using control over a subdomain defined with a SubMesh, but this wasn’t as straight forward as I thought it would be.

I also tried setting 3 volume constraints on each subdomain, where the non-design domains have volume bounds that are equal to their initial values, but those were violated during the optimization

Here’s a minimum working example without any submeshes and fixed measures:

``````from fenics import *
from mpi4py import MPI

rank = MPI.COMM_WORLD.rank
parameters["std_out_all_processes"] = False

class Rho(UserExpression):
def __init__(self, subdomains, **kwargs):
super(Rho, self).__init__()
self.subdomains = subdomains

def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] == 3:
values[0] = 0.4
else:
values[0] = 1

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

class HeatSink(SubDomain):
def inside(self, x, on_boundary):
return x[2] < 2

class HeaterPlate(SubDomain):
def inside(self, x, on_boundary):
return x[2] > 8

def forward(rho):
"""Solve the forward problem for a given material distribution rho(x)."""
T = Function(P, name="Temperature")
v = TestFunction(P)
F = (
- f * v * dx(2)  # heater source term
)
solve(F == 0, T, bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7,
"maximum_iterations": 20}})
return T

def eval_cb(j, a):
a_viz.assign(a)
controls << a_viz

class VolumeConstraint(InequalityConstraint):
"""A class that enforces the volume constraint g(a) = V - a*dx >= 0."""
def __init__(self, V):
self.V = V
self.smass = assemble(TestFunction(A) * Constant(1)* dx(3))
self.tmpvec = Function(A)
def function(self, m):
set_local(self.tmpvec, m)
integral = self.smass.inner(self.tmpvec.vector())

if MPI.COMM_WORLD.rank == 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 the number of components in the constraint vector (here, 3)."""
return 1

P1 = Point(0, 0, 0)
P2 = Point(10, 10, 10)
mesh = BoxMesh(P1, P2, 10, 10, 10)

p = Constant(5)  # power used in the solid isotropic material
eps = Constant(1.0e-3)  # epsilon used in the solid isotropic material
alpha = Constant(1.0e-8)  # regularisation coefficient in functional
A = FunctionSpace(mesh, "CG", 1)  # function space for control
P = FunctionSpace(mesh, "CG", 1)  # function space for solution

subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
subdomains.set_all(3)            # design domain 3
HeatSink().mark(subdomains, 1)  # baseplate  1
HeaterPlate().mark(subdomains, 2)  # heater  2
# File('subdomains.pvd') << subdomains

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
rho = interpolate(Rho(subdomains), A)

# rho_file = File('rho.pvd')
# rho_file << rho

bc = [DirichletBC(P, 20, HeatSink())]
f = interpolate(Constant(1.0e-2), P)  # the volume source term for the PDE
T = forward(rho)
file = File("solution.pvd")
file << T

controls = File("output/control_iterations.pvd")
a_viz = Function(A, name="ControlVisualisation")

J = assemble(f * T * dx(2)
)
m = Control(rho)
Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)

lb = 0.0
ub = 1.0
V3 = assemble(rho*dx(3))  # volume bound
parameters = {"acceptable_tol": 1.0e-3, "maximum_iterations": 100}
problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V3))
solver = IPOPTSolver(problem, parameters=parameters)
a_opt = solver.solve()
``````

Is there any way where I can define control over a subdomain?
Thanks!

It seems Submesh does not work, so I don’t think there is a way to define the control on a subdomain. Instead you can define your design variable for the whole domain but only have it contribute on the subdomain.

You achieve by the same approach as in your previous MWE by decomposing the integral into 3 parts:

``````         inner(grad(v), k(rho) * grad(T)) *dx(1)
``````

but this time `rho_sub` lives in `A` and not `ASub`.
The gradient for `rho_sub` should be zero in the interior of the other subdomains and there should be no change there during optimization.

1 Like

Thank you for taking the time to look into it.