Lawrence,

Thanks for all of your help. Here’s what I ended up with: 
——————
def ObjGrad_D(tao, petsc_x, petsc_g):
  with b_D_nobc.dat.vec as b_D_vec:
    A_D_nobc.M.handle.mult(petsc_x,petsc_g)
    xtHx = petsc_x.dot(petsc_g)
    xtf = petsc_x.dot(b_D_vec)
    petsc_g.axpy(-1.0,b_D_vec)
    arr = petsc_g.array
    for bc in bc_D:
      nodes = bc.nodes[bc.nodes < bc.function_space().dof_dset.size]
      if isinstance(bc.function_arg, Function):
        with bc.function_arg.dat.vec_ro as bc_vec:
          arr[nodes] = bc_vec.array_r[nodes]
      else:
        arr[nodes] = float(bc.function_arg)
    return 0.5*xtHx - xtf 
———————
However it seems something is still off about the boundary conditions I am applying. I have attached three different scenarios for my two solutions psi_A and psi_B. The boundary conditions are:

psi_A: set to 1 at x == 0.0 and 0.0 <= y <= 0.5, 0 everywhere else
psi_B: set to 1 at x == 0.0 and 0.5 < y <= 1.0, 0 everywhere else

- no_opt_psiA/B.png illustrates solving the standard diffusion without any optimization. As far as I know, these solutions are “correct” in terms of the boundary conditions.

- opt_wo_bounds_psiA/B.png illustrates solving the standard diffusion with TAO’s optimization but without any variable bounds. I would expect this to have the same solution as the no_opt cases but it seems to me the boundary conditions are messed up; face ID’s 3 and 4 (i.e., top and bottom) should have homogeneous bcs but it seems they are not registering.

- opt_w_bounds_psiA/B.png now includes the variable bounds set by these lines: 

lb = Function(W)
ub = Function(W)
lb.assign(0)
ub.assign(1)
tao = PETSc.TAO().create(MPI.COMM_WORLD)
with lb.dat.vec as lb_vec, ub.dat.vec as ub_vec:
  tao.setVariableBounds(lb_vec,ub_vec)

The boundary conditions are also screwed up at the top and bottom faces. In this case, psi_B (corresponding to W.sub(1)) is even more screwed up - the non-homogeneous part of psi_B isn’t even registering (hence concentrations of < 1).

Replacing the RHS of the arr[nodes] assignments with 0 didn’t change the solution either. Any idea what’s going on here?