Hi Lawrence,

Ah yes I meant to say pass in SNESVISetVariableBounds if that's what you're patch handles. I will try this out.

And thanks Andrew for your help. I ended up porting the undocumented SUPG advection-diffusion FEniCS example to Firedrake.

One more question:

4 I also want to work with a tensor diffusion coefficient. In FEniCS I could normally do something like

D = as_matrix([[1.0],[0,1]])

but i specifically want tensorial dispersion from velocity. E.g.,

D = alpha_t*norm(velocity)*Identity + (alpha_L-alpha_t)*outer(v,v)/norm(v)

Where alpha_t and alpha_L are the transverse and longitudinal dispersivity respectively (ignoring molecular diffusion for the moment). How would do I express this D matrix?

Thanks,

Justin

On Tue, Sep 1, 2015 at 1:39 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:

> 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)

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake