Dijitso JIT compilation failed on assembly of subdomain problem

dolfin-adjoint
1

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 *
from fenics_adjoint 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 = (
        inner(grad(v), k(rho_sub) * grad(T)) * dx(0)   # design domain
        +inner(grad(v), k(rho) * grad(T)) * dx(1)      # baseplate
        +inner(grad(v), k(rho) * grad(T)) * dx(2)      # heater
        - 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:

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

inside a conda enviroment.

Thanks in advance!

2

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.

3

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?

4

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).

5

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

6

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 *
from fenics_adjoint 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 = (
         inner(grad(v), k(rho) * grad(T)) *dx(1)
        +inner(grad(v), k(rho) * grad(T)) * dx(2)
        +inner(grad(v), k(rho_sub) * grad(T)) * dx(3)
         - 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):
        from pyadjoint.reduced_functional_numpy import set_local
        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)
             +alpha * inner(grad(rho_sub), grad(rho_sub)) * dx(3)
             )
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.

7

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.

8

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:

Screenshot from 2021-02-22 17-12-56

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 fenics_adjoint 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 = (
         inner(grad(v), k(rho) * grad(T)) * dx
         - 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):
        from pyadjoint.reduced_functional_numpy import set_local
        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)
             +alpha * inner(grad(rho), grad(rho)) * dx(3)
             )
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!

9

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)
        +inner(grad(v), k(rho) * grad(T)) * dx(2)
        +inner(grad(v), k(rho_sub) * grad(T)) * dx(3)

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
10

Thank you for taking the time to look into it.

It works with your method.

11

May I ask what commond did you use to install the devlopment version of dolfin, thank you in advance!

1 Like