Hello,
I would like to perform some topology optimization using the IPOPT optimizer and the automatic differentiation functionalities provided by dolfin-adjoint.
However, I would like the optimizer to account for a filtered/modified version of the gradient of my functional and not directly the one resulting from the operations saved on the Tape. This is referred to as “sensitivity-filtering” in the topology optimization field.
Do you have suggestions or advice on how to pass my custom gradient to the solver? And how can it take this filtering operation into account for every iteration of the optimization?
Here is a MWE based on the tutorial for topology optimization of heat conduction problems. I have indicated the default gradient and the desired one.
from __future__ import print_function
from dolfin import *
from dolfin_adjoint import *
from pyadjoint.drivers import compute_gradient
import matplotlib.pyplot as plt
import numpy as np
#import sklearn.metrics.pairwise as sp
try:
from pyadjoint import ipopt
except ImportError:
info_red("""This example depends on IPOPT and pyipopt. \
When compiling IPOPT, make sure to link against HSL, as it \
is a necessity for practical problems.""")
raise
parameters["std_out_all_processes"] = True
p = Constant(5) # power used in the solid isotropic material with penalisation (SIMP) rule
eps = Constant(1.0e-3) # epsilon used in the solid isotropic material
rmin = 0.05 # filter radius
alpha = Constant(1.0e-8) # regularisation coefficient in functional
V_star = Constant(0.4) # volume bound on the control
def k(a):
"""Solid isotropic material with penalisation (SIMP) conductivity
rule, equation (11)."""
return eps + (1 - eps) * a**p
n = 100
mesh = UnitSquareMesh(n, n)
A = FunctionSpace(mesh, "CG", 1) # function space for control
P = FunctionSpace(mesh, "CG", 1) # function space for solution
class WestNorth(SubDomain):
"""The top and left boundary of the unitsquare, used to enforce the Dirichlet boundary condition."""
def inside(self, x, on_boundary):
return (x[0] == 0.0 or x[1] == 1.0) and on_boundary
bc = [DirichletBC(P, 0.0, WestNorth())]
f = interpolate(Constant(1.0e-2), P) # the volume source term for the PDE
def forward(a):
"""Solve the forward problem for a given material distribution a(x)."""
T = Function(P, name="Temperature")
v = TestFunction(P)
F = inner(grad(v), k(a)*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
def J(a,T):
"""Functional of the optimization problem, with `a` control var and `T` temperature."""
return assemble(f*T*dx + alpha * inner(grad(a), grad(a))*dx)
def filter_sensitivity(s_ini,radius_f,FS):
"""PDE Helmholtz-type filter for sensitivity"""
v_f = TrialFunction(FS)
w_f = TestFunction(FS)
sys_left = (radius_f**2)*inner(grad(v_f), grad(w_f))*dx +v_f*w_f*dx
sys_right = s_ini*w_f*dx
sys_mat, sys_vec = assemble_system(sys_left,sys_right)
s_f = Function(FS)
solve(sys_mat, s_f.vector(), sys_vec, annotate=True)
return s_f
if __name__ == "__main__":
a = interpolate(V_star, A)
T = forward(a)
Jfunc = J(a,T)
m = Control(a)
Jhat = ReducedFunctional(Jfunc, m)
dJda = compute_gradient(Jfunc,Control(a)) # default gradient used to solve the minimization problem "behind the scenes"
dJda_f = filter_sensitivity(dJda,rmin,a.function_space()) # desired gradient (to be performed every iteration)
volume_constraint = UFLInequalityConstraint((V_star - a)*dx, m)
problem = MinimizationProblem(Jhat, bounds=(0.0, 1.0), constraints=volume_constraint)
solver = IPOPTSolver(problem, parameters={"acceptable_tol": 1.0e-3, "maximum_iterations": 100})
a_opt = solver.solve()
# Visualization
Jfunc_opt = J(a_opt,forward(a_opt))
dJda_opt = compute_gradient(Jfunc_opt,Control(a_opt))
dJda_opt_f = filter_sensitivity(dJda_opt,rmin,a_opt.function_space())
pl, ax = plt.subplots(); fig = plt.gcf(); fig.set_size_inches(12, 4)
plt.subplot(1,3,1); p = plot(a_opt, title='Topology', mode='color');cbar = plt.colorbar(p)
plt.subplot(1,3,2); p = plot(dJda_opt, title='Sensitivity', mode='color'); cbar = plt.colorbar(p)
plt.subplot(1,3,3); p = plot(dJda_opt_f, title='Sensitivity filtered', mode='color'); cbar = plt.colorbar(p)
Thank you for your time.