Hi all!
I am trying to perform an example of topology optimization of stokes flow for a diffuser. I am using the code below.
I am getting the following message:
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.461e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.122e-01 (tol = 1.000e-10) r (rel) = 7.684e-03 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 1.954e-01 (tol = 1.000e-10) r (rel) = 1.338e-02 (tol = 1.000e-09)....
Then, after 50 iterations, it stops.
Can anybody help me?
Thanks!
C. Okubo
from dolfin import *
from dolfin_adjoint import *
import pyipopt
## Geometry
L = 1.0
h = 1.0
## Optimisation Parameters
mu = Constant(1.0) # Viscosity
alphaunderbar = 2.5 * mu / (100**2) # Parameter for \alpha
alphabar = 2.5 * mu / (0.01**2) # Parameter for \alpha
q = Constant(0.01) # q for difficulty/discrete-valuedness
V = Constant(1.0/2) * L # 1/2 of the domain
def alpha(rho): return alphabar + (alphaunderbar - alphabar) * rho * (1 + q) / (rho + q)
## Mesh, Control and Solution Spaces
N = 50
mesh = RectangleMesh(Point(0.0, -h/2.0), Point(L, h/2.0), N, N)
A = FunctionSpace(mesh, "CG", 1) # Control
U = VectorFunctionSpace(mesh, "CG", 2) # Velocity
P = FunctionSpace(mesh, "CG", 1) # Pressure
W = MixedFunctionSpace([U, P]) # Mixed Taylor-Hood function
## Boundary conditions
tol = 1E-10
umax = 1.0
inflow = Expression(("umax*(1-(x[1]*x[1])/(h*h))","0.0"), umax = umax, h = h/2.0)
outflow = Expression(("umax*(1-(x[1]*x[1])/(h*h))","0.0"), umax = 3*umax, h = h/6.0)
L_out = L-tol
h_out = h/6.0
def Inflow(x, on_boundary): return on_boundary and x[0] < tol
def Outflow(x, on_boundary): return on_boundary and x[0] > L_out and -h_out < x[1] < h_out
bc1 = DirichletBC(W.sub(0), inflow, Inflow)
bc2 = DirichletBC(W.sub(0), outflow, Outflow)
## Forward problem solution
def forward(rho):
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
bc = [bc1, bc2]
F = (alpha(rho) * inner(u, v) * dx + inner(grad(u), grad(v)) * dx + inner(grad(p), v) * dx + inner(div(u), q) * dx)
solve(F == 0, w, bcs=bc)
return w
## MAIN
if __name__ == "__main__":
rho = interpolate(Constant(float(V)/L), A)
w = forward(rho)
(u, p) = split(w)
J = Functional(0.5 * inner(alpha(rho) * u, u) * dx + mu * inner(grad(u), grad(u)) * dx)
m = Control(rho) # Control
Jhat = ReducedFunctional(J, m) # Reduced Functional
lb = 0.0 # Inferior
ub = 1.0 # Superior
# Volume constraints
class VolumeConstraint(InequalityConstraint):
"""A class that enforces the volume constraint g(a) = V - a*dx >= 0."""
def __init__(self, V):
self.V = float(V)
self.smass = assemble(TestFunction(A) * Constant(1) * dx)
self.tmpvec = Function(A)
def function(self, m):
print "Evaluting constraint residual"
self.tmpvec.vector()[:] = m
integral = self.smass.inner(self.tmpvec.vector())
print "Current control integral: ", integral
return [self.V - integral]
def jacobian(self, m):
print "Computing constraint Jacobian"
return [-self.smass]
def output_workspace(self):
return [0.0]
# Solve the optimisation problem with q = 0.01
problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V))
parameters = {'maximum_iterations': 20}
solver = IPOPTSolver(problem, parameters=parameters)
rho_opt = solver.solve()
File("stokes_diffuser/control_solution_guess.xdmf") << rho_opt
q.assign(0.1)
rho.assign(rho_opt)
adj_reset()
File("intermediate-guess-%s.xdmf" % N) << rho
w = forward(rho)
(u, p) = split(w)
# Define the reduced functionals
rho_viz = Function(A)
J = Functional(0.5 * inner(alpha(rho) * u, u) * dx + mu * inner(grad(u), grad(u)) * dx)
m = Control(rho)
Jhat = ReducedFunctional(J, m)
problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V))
parameters = {'maximum_iterations': 100}
solver = IPOPTSolver(problem, parameters=parameters)
rho_opt = solver.solve()
File("stokes_diffuser/control_solution_final.xdmf") << rho_opt
