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