Separating optimization domain from the rest

Legacy DOLFIN
1

Hello,
In the context of my thesis I am looking to perform Density Based topological optimisations for fluids using fenics 2019 and ipopt. I would like to be able to restrict the optimisation of my design variable alpha only to a specific part of the computational domain. But the objective function is calculated on another part of the domain. I know how to divide my domain into sub-domains, for example I’ve created a “dx” that allows me to integrate the objective where I want, but I can’t get ipopt to modify the design variable in a specific domain. Is it possible to do this?

Thank you and have a nice day !

2

As you haven’t provided any code, this is a fairly general request, that is a bit hard to interpret.

Are you saying that you have an \alpha\in V(\Omega_0), where \Omega_0\subset\Omega, and your functional is J(u(\alpha))=\int_{\Omega_1}f(u(\alpha))\mathrm{d} x where \Omega_0\cap\Omega_1=\emptyset?

If so, what is your issue and how have you coded this up so far?

3

Thank you for your response and your time !
Yes, my problem is almost like this. More precisely, the design variable alpha is defined throughout the rectangular domain of length 1, while my objective function is calculated only on a vertical line at x=0.9 (i.e. close to the output). But I’d like the optimiser to be able to change the value of alpha only upstream (i.e. x < 0.80 for example).
For the moment, I’ve applied an inequality constraint which consists of having an average of alpha = 1 at the output so that the optimiser is obliged to leave fluid there. However, I’ve found quite a few papers that separate the optimisation zone from the objective zone but using fenics and ipopt I’ve found this, it’s not exactly like me but they manage to disable zones of the optimiser: Topology optimization of subsonic compressible flows | Structural and Multidisciplinary Optimization

Here are the main stages of my code, which I think will be of interest to you:

## The Mesh
mesh = RectangleMesh(Point(x_min, y_min), Point(x_max, y_max), N_mesh, 300, diagonal)

## Define sub-domain for BC and Optimizer
class Inlet(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[0] == x_min)
    
class Outlet(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[0] == x_max)

class Walls(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[1] == y_min or x[1] == y_max)
    
class Section_Obj(SubDomain):
        def inside(self, x, on_boundary):
            return (x[0] >= 0.9 * x_max-0.05) and (x[0] <= 0.9 * x_max+0.05)

class UpstreamZone(SubDomain):
        def inside(self, x, on_boundary):
            return (x[0] > 0.9 * x_max-0.1)

    marker_numbers = {'unset' : 0, 'wall' : 1, 'inlet' : 2, 'outlet' : 3, 'section_obj' : 4}
    boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
    boundary_markers.set_all(marker_numbers['unset'])
    Walls().mark(boundary_markers, marker_numbers['wall'])
    Inlet().mark(boundary_markers, marker_numbers['inlet'])
    Outlet().mark(boundary_markers, marker_numbers['outlet'])
    domains = MeshFunction("size_t", mesh, mesh.topology().dim())
    domains2 = MeshFunction("size_t", mesh, mesh.topology().dim())
    domains.set_all(0)
    domains2.set_all(0)
    Section_Obj().mark(domains, 6) # Objective section
    UpstreamZone().mark(domains2,7) #Optimizer section
    dx = Measure('dx', domain=mesh, subdomain_data=domains)
    dx2 = Measure('dx', domain=mesh, subdomain_data=domains2)

## NS weak form
I_ = as_tensor(np.eye(2)); sigma = -p*I_ + mu_*(grad(v) + grad(v).T);
F_PV = (rho_ * inner(v - v_prev1, w_v) / dt * dx
        + rho_ * (inner(grad(v) * v_prev1, w_v)+ inner(grad(v_prev1) * v, w_v)- inner(grad(v_prev1) * v_prev1, w_v)) * dx
        + div(v) * w_p * dx
        + inner(grad(w_v), sigma) * dx
        + inner(k(alpha) * v, w_v) * dx)
dF = derivative(F, u); problem = NonlinearVariationalProblem(F, u, bcs, dF)
        solver = NonlinearVariationalSolver(problem)
        solver.solve()

## Objective function
J = assemble(w_ud*dot(v - u_target, v - u_target) * dx(6) + inner(k(alpha) * v, v)*dx)

##Optimizer
alpha_C = Control(alpha)
Jhat = ReducedFunctional(J, alpha_C, derivative_cb_post = derivative_cb_post)
problem_min = MinimizationProblem(Jhat, bounds = (0.0, 1.0), constraints = [UFLInequalityConstraint((alpha-1)*dx2(7), alpha_C)])
max_iterations = 50 if idx == 0 else 150
 tol = 1e-6
solver_opt = IPOPTSolver(problem_min, parameters = {'maximum_iterations': max_iterations, 'tol':tol, 'print_level' : 5,'linear_solver' : 'ma86' })

# Perform topology optimization
alpha_opt = solver_opt.solve()

Have a nice day !

4

If you want to deactivate alpha in a certain region, you should make the bounds a function of the control space, and fix them in the region of interest.

Here is an example where you fix your control to 0.5 in a specific region:

    ub = Function(A)
    ut = Function(A)
    bc = DirichletBC(A, 0.5, "x[0] <0.2  && x[1] > 0.25 && x[1] < 0.75")
    ub.vector()[:] = 0.0
    ut.vector()[:] = 1.0
    bc.apply(ub.vector())
    bc.apply(ut.vector())

(based on dolfin-adjoint/examples/stokes-topology/stokes-topology.py at 5f020398c59173d879f9dbf3fda36e8acafc4afb · dolfin-adjoint/dolfin-adjoint · GitHub)

6

It works ! thank you !!!