On 31 Aug 2015, at 22:31, Justin Chang <jychang48@gmail.com> wrote:
Hi all,
2) I want to play around with PETSc's Variational Inequality solver. In my own PETSc code all I had to do was provide a Jacobian matrix, a residual vector, pass in the options -snes_type vinewtonssls (or it might be vinewotnrsls), and pass in SNESVINEWTONSSLS . From my last discussion with some PETSc folks it seems this VI feature is relatively now, so I was wondering if that last item can be invoked through the version of petsc4py you guys are using?
If it's just as simple as passing -snes_type vinewtonssls, then I think you're good to go. But presumably you also need to provide the SNES with the variable bounds. I think the DOLFIN interface for this (which we should mimic) is to do: solver = NLVSolver(...) solver.solve(lower_bound, upper_bound) where you provide Functions in the solution space defining your variable bounds (and can omit one or the other if you don't have a bound in that direction). The SNES solver needs to be informed of these, so if they're provided one must do: solver.snes.setVariableBounds(lb, ub) Here's a sketch patch (untested): Use it like: solver = ... lb = Function(V) lb.assign(-10) # Ensure solution is bounded below by -10. solver.solve(lb=lb) diff --git a/firedrake/variational_solver.py b/firedrake/variational_solver.py index ec47ce7..342c3e9 100644 --- a/firedrake/variational_solver.py +++ b/firedrake/variational_solver.py @@ -155,7 +155,18 @@ class NonlinearVariationalSolver(object): @timed_function("SNES solver execution") @profile - def solve(self): + def solve(self, lb=None, ub=None): + """Solve the variational problem. + + You can provide optional lower and upper bounds on the + solution vector if using PETSc's variational inequality + solvers. + + :arg lb: An optional lower bound on the solution vector. If + not provided, use -inf. + :arg ub: An optional upper bound on the solution vector. If + not provided, use +inf. + """ dm = self.snes.getDM() dm.setAppCtx(weakref.proxy(self._ctx)) dm.setCreateMatrix(self._ctx.create_matrix) @@ -168,6 +179,20 @@ class NonlinearVariationalSolver(object): # solve, ensure these are passed through to the snes. solving_utils.update_parameters(self, self.snes) + if lb is not None or ub is not None: + if ub: + with ub.dat.vec_ro as v: + ub_vec = v + else: + ub_vec = dm.createGlobalVector() + ub_vec.set(PETSc.INFINITY) + if lb: + with lb.dat.vec_ro as v: + lb_vec = v + else: + lb_vec = dm.createGlobalVector() + lb_vec.set(PETSc.NINFINITY) + self.snes.setVariableBounds(lb_vec, ub_vec) with self._problem.u.dat.vec as v: self.snes.solve(None, v)