Optimization only on subdomain

[Luís F]

I'm using IPOPT and the MinimizationProblem function to solve a topology optimization problem. Thus, I would like to know if it is possible to solve the optimization only in a subdomain.

For example, in the stokes topology optimization, optimizing only half of the domain:

p0 = Point()
p1 = Point(1.0,1.0, 0.0)
mesh = RectangleMesh(1.0,1.0, 10, 10,"right")

U = VectorFunctionSpace(mesh, "CG", 2) # velocity function space    -U
P = FunctionSpace(mesh, "CG", 1)       # pressure function space    -P
W = MixedFunctionSpace([U, P])         # mixed Taylor-Hood function space

class OmegaSub(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] > 0.5

osub = OmegaSub()

# create sumdomains
domains = CellFunction('size_t', mesh)
domains.set_all(0)
osub.mark(domains, 1)

# create submeshes
submesh = SubMesh(mesh, domains, 1)

# Define new measures associated with the interior domains
dx = Measure('dx', domain=mesh, subdomain_data=domains)

ASub = FunctionSpace(submesh, "CG", 1)

rhosub = interpolate(Constant(1.0), ASub)

Then,

J = Functional (...)

m = Control( rhosub )
Jhat = ReducedFunctional(J, m)

problem = MinimizationProblem(Jhat, bounds=(0.0,1.0), constraints=VolumeConstraint(V))
parameters = {"maximum_iterations": 100}
solver = IPOPTSolver(problem, parameters=parameters)
rho_opt = solver.solve()

Thank you,

Luís

[Simon Funke]

Sure, that should work just fine.

[Luís F]

Hi Simon, I've tried the subdomain optimization but I'm getting this error in the Volume Constraint:

    problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V))
  File "/media/luis/Novo volume/FEniCS/Juan/stokes_topology_ex.py", line 132, in __init__
    self.smass = assemble(TestFunction(ASub) * Constant(1) * dx)
  File "/usr/lib/python2.7/dist-packages/dolfin_adjoint/assembly.py", line 19, in assemble
    output = backend_assemble(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 178, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 58, in _create_dolfin_form
    return Form(form, form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/form.py", line 82, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 64, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 128, in jit
    return form_compiler.jit(form, parameters=p)
  File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 74, in jit
    return jit_form(ufl_object, parameters)
  File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 103, in jit_form
    module_name = "ffc_form_" + jit_object.signature()
  File "/usr/lib/python2.7/dist-packages/ffc/jitobject.py", line 78, in signature
    form_signature = self.form.signature()
  File "/usr/lib/python2.7/dist-packages/ufl/form.py", line 206, in signature
    self._compute_signature()
  File "/usr/lib/python2.7/dist-packages/ufl/form.py", line 384, in _compute_signature
    self._signature = compute_form_signature(self, self._compute_renumbering())
  File "/usr/lib/python2.7/dist-packages/ufl/algorithms/signature.py", line 176, in compute_form_signature
    terminal_hashdata = compute_terminal_hashdata(integrands, renumbering)
  File "/usr/lib/python2.7/dist-packages/ufl/algorithms/signature.py", line 74, in compute_terminal_hashdata
    data = expr.signature_data(renumbering)
  File "/usr/lib/python2.7/dist-packages/ufl/argument.py", line 100, in signature_data
    edata = self.element().signature_data(renumbering)
  File "/usr/lib/python2.7/dist-packages/ufl/finiteelement/finiteelement.py", line 107, in signature_data
    ("no domain" if self._domain is None else self._domain.signature_data(renumbering)))
  File "/usr/lib/python2.7/dist-packages/ufl/domain.py", line 191, in signature_data
    count = renumbering[self]
KeyError: Domain(Cell('triangle', 2), label='dolfin_mesh_with_id_14', data='<data with id 14>')

The modified stokes_topology example code is:

# -*- coding: utf-8 -*-

from dolfin import *
from dolfin_adjoint import *

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 value that controls difficulty/discrete-valuedness of solution

def alpha(rho):
    """Inverse permeability as a function of rho, equation (40)"""
    return alphabar + (alphaunderbar - alphabar) * rho * (1 + q) / (rho + q)


N = 20
delta = 1.5  # The aspect ratio of the domain, 1 high and \delta wide
V = Constant(1.0/3) * delta/2  # want the fluid to occupy 1/3 of the domain

mesh = RectangleMesh(Point(0.0, 0.0,0.0), Point(delta, 1.0,0.0), N, N)
A = FunctionSpace(mesh, "CG", 1)        # control function space
U = VectorFunctionSpace(mesh, "CG", 2)  # velocity function space
P = FunctionSpace(mesh, "CG", 1)        # pressure function space
W = MixedFunctionSpace([U, P])          # mixed Taylor-Hood function space


### SUBMESH ----------------------------------------------------------------###

class OmegaSub(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 1.0/2

osub = OmegaSub()

# create sumdomains
domains = CellFunction('size_t', mesh)
domains.set_all(0)
osub.mark(domains, 1)

# create submeshes
submesh = SubMesh(mesh, domains, 1)
plot(domains,key='subdomain')


dx = Measure('dx', domain=mesh, subdomain_data=domains)

ASub = FunctionSpace(submesh, "CG", 1)
rho_sub = interpolate(Constant(1.0),ASub) #Define RHO only on submesh

### SUBMESH END-------------------------------------------------------------###


class InflowOutflow(Expression):
    def eval(self, values, x):
        values[1] = 0.0
        values[0] = 0.0
        l = 1.0/6.0
        gbar = 1.0

        if x[0] == 0.0 or x[0] == delta:
            if (1.0/4 - l/2) < x[1] < (1.0/4 + l/2):
                t = x[1] - 1.0/4
                values[0] = gbar*(1 - (2*t/l)**2)
            if (3.0/4 - l/2) < x[1] < (3.0/4 + l/2):
                t = x[1] - 3.0/4
                values[0] = gbar*(1 - (2*t/l)**2)

    def value_shape(self):
        return (2,)

class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary 

walls = Wall()
facet_marker = FacetFunction("uint", mesh)
facet_marker.set_all(0)
walls.mark(facet_marker,1) 

# Define new measures associated with the exterior boundaries
ds = Measure("ds")[facet_marker]



def forward(rho_sub):
    """Solve the forward problem for a given fluid distribution rho(x)."""
    w = Function(W)
    (u, p) = split(w)
    (v, q) = TestFunctions(W)

    F = (alpha(rho_sub) * inner(u, v) * dx(1) + inner(grad(u), grad(v)) * dx +
         inner(grad(p), v) * dx  + inner(div(u), q) * dx)
    bc = DirichletBC(W.sub(0), InflowOutflow(),facet_marker,1)
    solve(F == 0, w, bcs=bc)

    return w

if __name__ == "__main__":
    rho_sub = interpolate(Constant(float(V)/delta), ASub, name="Control")
    w   = forward(rho_sub)
    (u, p) = split(w)


    controls = File("output/control_iterations_guess.pvd")
    allctrls = File("output/allcontrols.pvd")
    rho_viz = Function(ASub, name="ControlVisualisation")
    def eval_cb(j, rho):
        rho_viz.assign(rho)
        controls << rho_viz
        allctrls << rho_viz

    #Changed rho to rho sub and dx to dx(1)
    J = Functional(0.5 * inner(alpha(rho_sub) * u, u) * dx(1) + mu * inner(grad(u), grad(u)) * dx(1))
    m = Control(rho_sub)
    Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)

    # Bound constraints
    lb = 0.0
    ub = 1.0

    # 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(ASub) * Constant(1) * dx)
            self.tmpvec = Function(ASub)

        def function(self, m):
            print "Evaluting constraint residual"
            self.tmpvec.vector()[:] = m

            # Compute the integral of the control over the domain
            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("output/control_solution_guess.xdmf") << rho_opt

    q.assign(0.1)
    rho_sub.assign(rho_opt)
    adj_reset()

    File("intermediate-guess-%s.xdmf" % N) << rho_sub

    w = forward(rho_sub)
    (u, p) = split(w)

    # Define the reduced functionals
    controls = File("output/control_iterations_final.pvd")
    rho_viz = Function(ASub, name="ControlVisualisation")
    def eval_cb(j, rho):
        rho_viz.assign(rho)
        controls << rho_viz
        allctrls << rho_viz

    #Changed rho to rho sub and dx to dx(1)
    J = Functional(0.5 * inner(alpha(rho_sub) * u, u) * dx(1) + mu * inner(grad(u), grad(u)) * dx(1))
    m = Control(rho_sub)
    Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)

    problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V))
    parameters = {'maximum_iterations': 50}

    solver = IPOPTSolver(problem, parameters=parameters)
    rho_opt = solver.solve()

    File("output/control_solution_final.xdmf") << rho_opt

[Simon Funke]

An easy fix is to define the control Function on the entire domain, i.e. ASub = FunctionSpace(mesh, "CG", 1).

Your equations only use the control function in the upper part of the domain (thanks to your measure), and consequently the functional gradient will always be zero in the lower domain.

[Francisco Oliveira]

Hi Simon,

Its very usefull this method to select the region to optimize. I tried to change the original code as above, but still have some mistakes , could you to check it if possible. Thanks in advanced, Francisco

ERROR

---

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/usr/lib/python2.7/dist-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 540, in runfile execfile(filename, namespace) File "/Sample/stokes-topology_SubDomain.py", line 113, in <module> w = forward(rho_sub) File "/Sample/stokes-topology_SubDomain.py", line 107, in forward solve(F == 0, w, bcs=bc) File "/usr/lib/python2.7/dist-packages/dolfin_adjoint/solving.py", line 357, in solve ret = backend.solve(args, kwargs) File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 297, in solve _solve_varproblem(args, **kwargs) File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 345, in _solve_varproblem solver.solve() RuntimeError:

---

from dolfin import * from dolfin_adjoint import *

try: import pyipopt 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

turn off redundant output in parallel

parameters["std_out_all_processes"] = False

mu = Constant(1.0) # viscosity alphaunderbar = 2.5 * mu / (1002) # parameter for \alpha alphabar = 2.5 * mu / (0.012) # parameter for \alpha q = Constant(0.01) # q value that controls difficulty/discrete-valuedness of solution

def alpha(rho): """Inverse permeability as a function of rho, equation (40)""" return alphabar + (alphaunderbar - alphabar) * rho * (1 + q) / (rho + q)

N = 20 delta = 1.5 # The aspect ratio of the domain, 1 high and \delta wide V = Constant(1.0/3) * delta # want the fluid to occupy 1/3 of the domain

mesh = RectangleMesh(Point(0.0, 0.0,0.0), Point(delta, 1.0,0.0), N, N) A = FunctionSpace(mesh, "CG", 1) # control function space U = VectorFunctionSpace(mesh, "CG", 2) # velocity function space P = FunctionSpace(mesh, "CG", 1) # pressure function space W = MixedFunctionSpace([U, P]) # mixed Taylor-Hood function space

SUBMESH ----------------------------------------------------------------

class OmegaSub(SubDomain): def inside(self, x, on_boundary): return x[1] >= 1.0/2

osub = OmegaSub()

create sumdomains

domains = CellFunction('size_t', mesh) domains.set_all(0) osub.mark(domains, 1)

create submeshes

submesh = SubMesh(mesh, domains, 1) plot(domains,key='subdomain')

dx = Measure('dx', domain=mesh, subdomain_data=domains)

ASub = FunctionSpace(submesh, "CG", 1) rho_sub = interpolate(Constant(1.0),ASub) #Define RHO only on submesh

SUBMESH END-------------------------------------------------------------

Define the boundary condition on velocity

class InflowOutflow(Expression): def eval(self, values, x): values[1] = 0.0 values[0] = 0.0 l = 1.0/6.0 gbar = 1.0

if x[0] == 0.0 or x[0] == delta:
  if (1.0/4 - l/2) < x[1] < (1.0/4 + l/2):
    t = x[1] - 1.0/4
    values[0] = gbar*(1 - (2*t/l)**2)
  if (3.0/4 - l/2) < x[1] < (3.0/4 + l/2):
    t = x[1] - 3.0/4
    values[0] = gbar*(1 - (2*t/l)**2)

def value_shape(self): return (2,)

class Wall(SubDomain): def inside(self, x, on_boundary): return on_boundary

walls = Wall() facet_marker = FacetFunction("uint", mesh) facet_marker.set_all(0) walls.mark(facet_marker,1)

Define new measures associated with the exterior boundaries

ds = Measure("ds")[facet_marker]

def forward(rho): """Solve the forward problem for a given fluid distribution rho(x).""" w = Function(W) (u, p) = split(w) (v, q) = TestFunctions(W)

F = (alpha(rho) * inner(u, v) * dx + inner(grad(u), grad(v)) * dx + inner(grad(p), v) * dx + inner(div(u), q) * dx) bc = DirichletBC(W.sub(0), InflowOutflow(),facet_marker,1) solve(F == 0, w, bcs=bc)

return w

if name == "main": rho_sub = interpolate(Constant(float(V)/delta), ASub, name="Control") w = forward(rho_sub) (u, p) = split(w)

controls = File("output/control_iterations_guess.pvd") allctrls = File("output/allcontrols.pvd") rho_viz = Function(ASub, name="ControlVisualisation") def eval_cb(j, rho_sub): rho_viz.assign(rho_sub) controls << rho_viz allctrls << rho_viz

J = Functional(0.5 * inner(alpha(rho_sub) * u, u) * dx + mu * inner(grad(u), grad(u)) * dx) m = Control(rho_sub) Jhat = ReducedFunctional(J, m)

lb = 0.0 ub = 1.0

# 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(ASub) * Constant(1) * dx)
  self.tmpvec = Function(ASub)

def function(self, m):
  print "Evaluting constraint residual"
  self.tmpvec.vector()[:] = m

  # Compute the integral of the control over the domain
  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]

problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V)) parameters = {'maximum_iterations': 20}

solver = IPOPTSolver(problem, parameters=parameters) rho_opt = solver.solve()

File("output/control_solution_guess.xdmf") << rho_opt

q.assign(0.1) rho_sub.assign(rho_opt) adj_reset()

File("intermediate-guess-%s.xdmf" % N) << rho_sub

w = forward(rho_sub) (u, p) = split(w)

# Define the reduced functionals controls = File("output/control_iterations_final.pvd") rho_viz = Function(ASub, name="ControlVisualisation") def eval_cb(j, rho_sub): rho_viz.assign(rho_sub) controls << rho_viz allctrls << rho_viz

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("output/control_solution_final.xdmf") << rho_opt