Program "Killed" when running compute_gradient() on transient problem

dolfin-adjoint
1

Im writting a code to perform topology optimization on a heatsink design using FEniCS 2019.1 on Ubuntu. For this I need to calculate the sensitivities of an objective function that depends on a transient simulation and I do this with the dolphin-adjoint function compute_gradient(). An MWE is given bellow:

# ##################################################### #
#                       LIBRARIES                       #
# ##################################################### #

from fenics import *
from fenics_adjoint import *    # Note: It is important that fenics_adjoint is called after fenics

# ##################################################### #
#                         DATA                          #
# ##################################################### #

# ----- SPACIAL DATA -----
lx = ly = 1    # Dimensions of computational domain

# ----- THERMAL DATA -----
qheat = 2           # Heat rate [W]
Tinitial = 0        # Initial temperature

# ----- DISCRETIZATION DATA -----
nx = ny = 100       # Number of element in x and y direction
degree_phys = 1     # Element order for physical probelem
niter = 100         # Max number of iterations
num_steps = 250     # number of timespteps

# ----- TIME DATA -----
tfin = 20               # Final time
dt = tfin / num_steps  # time step size [s]

# ----- OPTIMIZATION DATA -----
volfrac = 0.3       # Maximum allowable volume fraction used in the optimisation           

# ----- DEFINING DOMAIN AND FUNCTION SPACES -----
# Create mesh and define function space
mesh = RectangleMesh(fenics.Point(-lx/2, -ly/2), fenics.Point(lx/2, ly/2), nx, ny, "crossed")

# Function space
V0 = FunctionSpace(mesh, "DG", 0)                 # Discontinous function space for rho (Values stored in elements)
Vphys = FunctionSpace(mesh, "CG", degree_phys)    # Continuous function space for solution (Values stored in nodes)

# ----- BOUNDARY CONDITIONS -----
# Applying Neumann boundary condition to simulate the heat generation from electronics
class HeatSourceBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol_bc = 1E-14
        return on_boundary and near(x[1], -ly/2, tol_bc) and x[0]>=-lx/4 and x[0]<=lx/4
HS_boundary = HeatSourceBoundary()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)  # store all facets in array
HS_boundary.mark(boundaries, 1)     # Mark all facest that are on the HS_boundary with 1
# Denfine the line intrgrals
dsHS = ds(subdomain_id=1, subdomain_data=boundaries)    

# ##################################################### #
#                     INITIAL DESIGN                    #
# ##################################################### #

# ----- INITIALIZE DENSITIES -----
rho_tilde = Function(V0, name="Filtered Density")       # Filtered design parameter
rho_tilde.interpolate(Constant(volfrac))                # Uniform distribution of material

# ----- DEFINING PROBLEM -----
# Initialize functions
T = TrialFunction(Vphys)
v = TestFunction(Vphys)
T_n = Function(Vphys, name="Previous Temperature") 

# Weak form of transient heat transfer problem
F = rho_tilde*(T-T_n)/dt*v*dx + dot(grad(T),grad(v))*dx - qheat*v*dsHS

# a: unknowns, L: knowns
a, L = lhs(F), rhs(F)

# ##################################################### #
#                      OPTIMISATION                     #
# ##################################################### #

# Initialize the functionspace for the Temperature and sensitivities
T = Function(Vphys, name="Temperature")
df0dx = Function(V0, name="Object function Sensitivity")

# ----- OPTIMZATION LOOP -----
for iter in range(niter):
    # ------ INITILIZE PHYSICS  -----
    T_n.assign(interpolate(Constant(Tinitial), Vphys))
    T.assign(interpolate(Constant(Tinitial), Vphys))
    t = 0
    T_array = [T]
    T_Electronics_array = [assemble(T*dsHS)]    # Save the Temperature of the eletronics

    # ------ SOLVING TRANSIENT PROBLEM  -----
    while (t <= tfin-dt):
        t += float(dt)      # Update time
        solve(a == L, T)    # Solve problem for current time

        T_n.assign(T)           # Update previous tempertature

    # ----- COMPUTING OBJECTIVE FUNCTION AND SENSITIVITIES -----
    f0 = assemble((T-Tinitial)*dx)               # Objective function

    df0dx.assign(compute_gradient(f0, Control(rho_tilde)))      # Sensitivity of objective function wrt rho_tilde

When running the code I get the following message after around 2 hours (on 6GB RAM).

image

By running the code on a computer with more RAM available I can extend this time, leading me to believe that it is a memory issue. Since the issue does not occur when omitting the the compute_gradient() function I believe it to be the culprit.

I have read that dolphin-adjoint uses a tape to keep track on the dependencies in order to compute the sensitivities with compute_gradient(). My theory is that this tape keeps storing dependencies and these decencies carry over optimization iteration leading to a accumulation of data. But I only have a limited knowledge of data structures so I can’t be sure.

Is there anyone that can confirm or deny this suspicion and can help me solve it? Maybe it is it possible to clear this tape after each optimization iteration or tell dolphin-adjoint what to store in the tape?

Thanks in advance.

2

This is not really how you should use dolfin-adjoint.
You should define a Reduced functional prior to the loop, where rho_tilde is the control and evaluate the functional to get the value at the next iteration.
Consider the following minimal Poisson example:


from dolfin import *

from dolfin_adjoint import *
# Create mesh, refined in the center
n = 16
mesh = Mesh(UnitSquareMesh(n, n))

cf = MeshFunction("bool", mesh, mesh.geometry().dim())
subdomain = CompiledSubDomain('std::abs(x[0]-0.5) < 0.25 && std::abs(x[1]-0.5) < 0.25')
subdomain.mark(cf, True)
mesh = refine(mesh, cf)

# Define discrete function spaces and funcions
V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 0)

f = interpolate(Expression("0.0", degree=1, name='Control'), W)
u = Function(V, name='State')
v = TestFunction(V)

# Define and solve the Poisson equation to generate the dolfin-adjoint annotation
F = (inner(grad(u), grad(v)) - f * v) * dx
bc = DirichletBC(V, 0.0, "on_boundary")
solve(F == 0, u, bc)

# Define functional of interest and the reduced functional
x = SpatialCoordinate(mesh)
w = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=3)
d = 1 / (2 * pi ** 2)
d = Expression("d*w", d=d, w=w, degree=3)

alpha = Constant(1e-6)
J = assemble((0.5 * inner(u - d, u - d)) * dx + alpha / 2 * f ** 2 * dx)

control = Control(f)
rf = ReducedFunctional(J, control)

import pyadjoint
with pyadjoint.stop_annotating() as _:
    num_iters = 200
    step = 100
    set_log_level(LogLevel.ERROR)
    rho_i = Function(W)
    rho_i.assign(f)
    for i in range(num_iters):
        Ji = rf(rho_i)
        dJdrhoi = rf.derivative()
        print(f"Iterate {i}: {Ji}")
        rho_i.vector()[:] -= step * dJdrhoi.vector()[:]
3

I have read the example, but I use a MMA function provided by Arjen Deetman
https://github.com/arjendeetman/GCMMA-MMA-Python for the optimization to have a little more control over the process. This requires the sensitivities to be supplied in form of an array as an input, hence the use of compute_gradient().

If possible I would like to keep on using the MMA function by Arjen Deetman.

4

That is the exact same output as

gives you, or dJdrhoi.vector().get_local() to get the corresponding numpy array

5

Oh sorry I might not have read the example you posted completely through, it was exactly what I needed. From the tests performed it seems to works perfectly, thank you.