Error when setting up bounds / Not the right number of controls

[Fin Bauer]

I want to change the example poisson-topology.py slightly by setting up different bounds. Specifically, I want to guarantee that my solution is material at certain points by setting lb=ub=1. I do this, for example, via the following function:

def setup_bounds(a):
    aVec = a.vector()
    lb = np.zeros(len(aVec))
    ub = np.ones(len(aVec))
    lb[::100] = 1
    bounds = zip(lb,ub)
    return bounds

Then, when I run the program with these bounds, i.e. I run

from dolfin import *
from dolfin_adjoint import *
import pyipopt
import numpy as np

# turn off redundant output in parallel
parameters["std_out_all_processes"] = False

V = Constant(0.4)      # volume bound on the control
p = Constant(5)        # power used in the solid isotropic material
                       # with penalisation (SIMP) rule, to encourage the control
                       # solution to attain either 0 or 1
eps = Constant(1.0e-3) # epsilon used in the solid isotropic material
alpha = Constant(1.0e-8) # regularisation coefficient in functional


def k(a):
    """Solid isotropic material with penalisation (SIMP) conductivity
  rule, equation (11)."""
    return eps + (1 - eps) * a**p

n = 50
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, name="SourceTerm") # 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 setup_bounds(a):
    aVec = a.vector()
    lb = np.zeros(len(aVec))
    ub = np.ones(len(aVec))
    lb[::100] = 1
    bounds = zip(lb,ub)
    return bounds

if __name__ == "__main__":
    a = interpolate(V, A, name="Control") # initial guess.
    T = forward(a)                        # solve the forward problem once.


    controls = File("output/control_iterations.pvd")
    a_viz = Function(A, name="ControlVisualisation")
    def eval_cb(j, a):
        a_viz.assign(a)
        controls << a_viz

    J = Functional(f*T*dx + alpha * inner(grad(a), grad(a))*dx)
    m = Control(a)
    Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)

    bounds = setup_bounds(a)

    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):
            self.tmpvec.vector()[:] = m
            integral = self.smass.inner(self.tmpvec.vector())
            if MPI.rank(mpi_comm_world()) == 0:
                print "Current control integral: ", integral
            return [self.V - integral]

        def jacobian(self, m):
            return [-self.smass]

        def output_workspace(self):
            return [0.0]

        def length(self):
            """Return the number of components in the constraint vector (here, one)."""
            return 1

    problem = MinimizationProblem(Jhat, bounds=bounds, constraints=VolumeConstraint(V))
    parameters = {"acceptable_tol": 1.0e-200, "maximum_iterations": 100}

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

I get the error

TypeError: bounds should be of length number of controls of the ReducedFunctional

because Jhat.controls is of length 1 while my number of bounds is of course way larger.

Is there a way to fix this problem, by somehow respecifying the number of controls for example or by passing bounds differently?

I'm using Fenics 1.6 and Dolfin-Adjoint 1.6

[Simon Funke]

The bounds should be a (lower_bound, upper_bound) tuple, where lower_bound and upper_bound are of the same type as the control. In your case:

def setup_bounds(a):                                                                                                                                                                                                 
    lb = interpolate(Expression("x[0]*x[1]"), a.function_space())                                                                                                                                                    
    ub = interpolate(Constant(1), a.function_space())                                                                                                                                                                
    return lb, ub