On Tue, 1 Sep 2015 at 09:09 Justin Chang <jychang48@gmail.com> wrote:

Hi Lawrence,

Ah yes I meant to say pass in SNESVISetVariableBounds <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/SNES/SNESVISetVariableBounds.html> 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.

If you feel like documenting it, we'd love to have it on the Firedrake website ;).

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