dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jun 14, 2011 at 09:25:17AM +0100, Garth N. Wells wrote:
>
>
> On 14/06/11 08:53, Anders Logg wrote:
> >
> > 14 jun 2011 kl. 09:18 skrev "Garth N. Wells" <gnw20@xxxxxxxxx>:
> >
> >>
> >>
> >> On 14/06/11 08:03, Marie E. Rognes wrote:
> >>> On 06/13/2011 11:16 PM, Anders Logg wrote:
> >>>>>> But while we are heading in that direction, how about
> >>>>>> abolishing the *Problem class(es) altogether, and just use
> >>>>>> LinearVariationalSolver and
> >>>>>> NonlinearVariationalSolver/NewtonSolver taking as input (a,
> >>>>>> L,
> >>>>> bc)
> >>>>>> and (F, dF, bcs), respectively.
> >>>> This will be in line with an old blueprint. We noted some time
> >>>> ago that problems/solvers are designed differently for linear
> >>>> systems Ax = b than for variational problems a(u, v) = L(v).
> >>>> For linear systems, we have solvers while for variational
> >>>> problems we have both problem and solver classes.
> >>>>
> >>>>>> I mean, the main difference lies in how to solve the
> >>>>>> problems, right?
> >>>> It looks like the only property a VariationalProblem has in
> >>>> addition to (forms, bc) + solver parameters is the parameter
> >>>> symmetric=true/false.
> >>>>
> >>>> If we go this route, we could mimic the design of the linear
> >>>> algebra solvers and provide two different options, one that
> >>>> offers more control, solver = KrylovSolver() + solver.solve(),
> >>>> and one quick option that just calls solve:
> >>>>
> >>>> 1. complex option
> >>>>
> >>>> solver = LinearVariationalSolver() # which arguments to
> >>>> constructor? solver.parameters["foo"] = ... u = solver.solve()
> >>>>
> >>
> >> I favour this option, but I think that the name
> >> 'LinearVariationalSolver' is misleading since it's not a
> >> 'variational solver', but solves variational problems, nor should
> >> it be confused with a LinearSolver that solves Ax = f.
> >> LinearVariationalProblem is a better name. For total control, we
> >> could have a LinearVariationalProblem constructor that accepts a
> >> GenericLinearSolver as an argument (as the NewtonSolver does).
> >>
> >>>
> >>> For the eigensolvers, all arguments go in the call to solve.
> >>>
> >>>> 2. simple option
> >>>>
> >>>> u = solve(a, L, bc)
> >>>>
> >>>
> >>
> >> I think that saving one line of code and making the code less
> >> explicit isn't worthwhile. I can foresee users trying to solve
> >> nonlinear problems with this.
> >
> > With the syntax suggested below it would be easy to check for errors.
> >
> >
> >>
> >>> Just for linears?
> >>>
> >>>> 3. very tempting option (simple to implement in both C++ and
> >>>> Python)
> >>>>
> >>>> u = solve(a == L, bc)    # linear u = solve(F == 0, J, bc) #
> >>>> nonlinear
> >>>>
> >>>
> >>
> >> I don't like this on the same grounds that I don't like the
> >> present design. Also, I don't follow the above syntax
> >
> > I'm not surprised you don't like it. But don't understand why. It's
> > very clear which is linear and which is nonlinear. And it would be
> > easy to check for errors. And it would just be a thin layer on top of
> > the very explicit linear/nonlinear solver classes. And it would
> > follow the exact same design as for la with solver classes plus a
> > quick access solve function.
> >
>
> Is not clear to me - possibly because, as I wrote above, I don't
> understand the syntax. What does the '==' mean?

Here's how I see it:

1. Linear problems

  solve(a == L, bc)

  solve the linear variational problem a = L subject to bc

2. Nonlinear problems

  solve(F == 0, bc)

  solve the nonlinear variational problem F = 0 subject to bc

It would be easy to in the first case check that the first operand (a)
is a bilinear form and the second (L) is a linear form.

And it would be easy to check in the second case that the first
operand (F) is a linear form and the second is an integer that must be
zero.

In both cases one can print an informative error message and catch any
pitfalls.

The nonlinear case would in C++ accept an additional argument J for
the Jacobian (and in Python an optional additional argument):

  solve(F == 0, J, bc);

The comparison operator == would for a == L return an object of class
LinearVariationalProblem and in the second case
NonlinearVariationalProblem. These two would just be simple classes
holding shared pointers to the forms. Then we can overload solve() to
take either of the two and pass the call on to either
LinearVariationalSolver or NonlinearVariationalSolver.

I'm starting to think this would be an ideal solution. It's compact,
fairly intuitive, and it's possible to catch errors.

The only problem I see is overloading operator== in Python if that
has implications for UFL that Martin objects to... :-)

--
Anders