I decided to radically edit my answer based on some of the comments.
I haven't used TAO. From perusing the documentation, it seems like the only way that TAO can handle constrained optimization problems (excluding the special case of only box constraints) is to convert the problem into a variational inequality using the Karush-Kuhn-Tucker (KKT) conditions, which are necessary under constraint qualification (the type that I usually see is the Slater point condition), and sufficient under convexity of the objective and constraints (more generally, Type 1 invexity). If $f$ is nonconvex, the variational inequality formulation using the KKT conditions is NOT equivalent to the original optimization problem, so strictly speaking, if you want a global optimum for the optimization problem, you should not express it as a variational inequality. It would be difficult to find a global optimum anyway for PDE-constrained optimization (see below), so maybe ignoring this detail is fine. Given what Wolfgang has said, I'd be skeptical of using TAO; I'm already skeptical because it doesn't implement methods for solving nonlinear programs (NLPs) as NLPs, rather than variational inequalities.
I am not an expert on PDE-constrained optimization; co-workers and associates of mine work on ODE-constrained optimization problems. I know that for intrusive formulations, Larry Biegler (and others) will use collocation methods to discretize the PDE and make it a very large, sparse NLP, and then he will solve it using interior point methods. To really solve the problem to global optimality, you would also need to generate convex relaxations, but as far as I know, this approach is not taken because PDE-constrained optimization problems lead to such large NLPs that it would be difficult to solve them to global optimality. I mention these details only because problem formulation heavily influences choice of solver package. For nonintrusive formulations, repeated PDE solves yield gradient information for optimization algorithms.
Some people who study ODE-constrained optimization problems use a similar approach of discretizing the problem using collocation and a numerical method, and then relaxing the resulting NLP to yield a convex formulation used in a global optimization algorithm. An alternate approach to ODE-constrained optimization is to relax the problem, and then discretize the ODE, which is the approach taken in my lab. It could be possible to relax certain classes of PDE-constrained optimization problems, but I don't know of any extant work done on that problem. (It was a potential project in my lab at one point.)
Ultimately, what matters is not the differentiability of the original PDE, but the differentiability of the discretization with respect to the decision variables.
If the discretized problem is twice-differentiable with respect to the decision variables, the following packages will calculate a local solution:
- IPOPT is an open source interior point solver developed by Andreas Wächter at IBM. It's a very high quality code. As an interior-point solver, it is better for objective functions with large, sparse Jacobian matrices, and would be useful for PDE-constrained optimization
- SNOPT is a commercial sequential quadratic programming solver that is another high quality code. It is better for objective functions with small, dense Jacobian matrices, so I would not expect it to be useful for PDE-constrained optimization, but you could give it a try.
- NLopt is a small, open source code written by Steven Johnson at MIT that contains basic implementations of a number of nonlinear optimization algorithms. All of the algorithms should be adequate for solving bound-constrained problems.
fmincon in Matlab implements a number of algorithms (including interior point and sequential quadratic programming) for nonlinear optimization
- GAMS and AMPL are both commercial modeling languages used to formulate optimization problems, and contain interfaces to a large number of nonlinear programming solvers. I know that GAMS has a trial version that can be used for smaller problems, and problem instances can also be submitted to the NEOS server for solution as well.
However, it's possible that the discretization is only once differentiable with respect to the decision variables, in which case, you should use projected steepest descent or some other first-order optimization method when calculating a local solution. Since a lot of study focuses on problems where second-order methods can be used (and when you can use them, their superior convergence properties make them a better choice), I couldn't find many implementations of steepest descent that weren't solutions to homework problems. The GNU Scientific Library has an implementation, but it's only for demonstration purposes. You would probably need to code up your own implementation.
If the problem is only continuous with respect to the decision variables, then you could use direct methods to solve it locally. There is an excellent survey on direct methods by Kolda, Lewis, and Torczon. The most widely known of these methods is the Nelder-Mead simplex algorithm. It is not guaranteed to converge to a local minimum in multiple dimensions, but it has found considerable practical use anyway.
Choice of package really depends on the language you want to use to solve the problem, if solving the bound-constrained problem is only part of an algorithm you want to implement (or if it's the only step in your algorithm, in which case modeling languages become more viable for production code), the type and size of the problem, and if you need any parallelism.