dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jun 14, 2011 at 12:50:07PM -0700, Johan Hake wrote:
> On Tuesday June 14 2011 12:03:31 Marie E. Rognes wrote:
> > On 06/14/2011 08:35 PM, Garth N. Wells wrote:
> > > On 14/06/11 19:24, Anders Logg wrote:
> > >> On Tue, Jun 14, 2011 at 10:19:20AM -0700, Johan Hake wrote:
> > >>> On Tuesday June 14 2011 03:33:59 Anders Logg wrote:
> > >>>> 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... :-)
> > >>>
> > >>> Wow, you really like magical syntaxes ;)
> > >>
> > >> Yes, a pretty syntax has been a priority for me ever since we
> > >> started. I think it is worth a lot.
> > >
> > > Magic and pretty are not the same thing.
> > >
> > >>> The problem with this syntax is that who on earth would expect a
> > >>> VariationalProblem to be the result of an == operator...
> > >>
> > >> I don't think that's an issue. Figuring out how to solve variational
> > >> problems is not something one picks up by reading the Programmer's
> > >> Reference. It's something that will be stated on the first page of any
> > >> FEniCS tutorial or user manual.
> > >>
> > >> I think the solve(a == L) is the one missing piece to make the form
> > >> language complete. We have all the nice syntax for expressing forms in
> > >> a declarative way, but then it ends with
> > >>
> > >> problem = VariationalProblem(a, L)
> > >> problem.solve()
> > >>
> > >> which I think looks ugly. It's not as extreme as this example taken
> > >> from cppunit, but it follows the same "create object, call method on
> > >>
> > >> object" paradigm which I think is ugly:
> > >>    TestResult result;
> > >>    TestResultCollector collected_results;
> > >>    result.addListener(&collected_results);
> > >>    TestRunner runner;
> > >>    runner.addTest(CppUnit::TestFactoryRegistry::getRegistry().makeTest()
> > >>    ); runner.run(result);
> > >>    CompilerOutputter outputter(&collected_results, std::cerr);
> > >>    outputter.write ();
> > >>>
> > >>> I see the distinction between FEniCS developers who have programming
> > >>> versus math in mind when designing the api ;)
> > >>
> > >> It's always been one of the top priorities in our work on FEniCS to
> > >> build an API with the highest possible level of mathematical
> > >> expressiveness to the API. That sometimes leads to challenges, like
> > >> needing to develop a special form language, form compilers, JIT
> > >> compilation, the Expression class etc, but that's the sort of thing
> > >> we're pretty good at and one of the main selling points of FEniCS.
>
> We all agree on that, and you have been the top proponent for this. But FEniCS
> would not have been so pretty and usefull without resistance.

I feel the resistance has grown recently. :-)

> > > This is an exaggeration to me. The code
> > >
> > >    problem = [Linear]VariationalProblem(a, L)
> > >    u = problem.solve()
> > >
> > > is compact and explicit. It's a stretch to call it ugly.
>
> Agree.
>
> > >>> Also __eq__ is already used in ufl.Form to compare two forms.
> > >>
> > >> I think it would be worth replacing the use of form0 == form1 by
> > >> repr(form0) == repr(form1) in UFL to be able to use __eq__ for this:
> > >>
> > >> class Equation:
> > >>    def __init__(self, lhs, rhs):
> > >>        self.lhs = lhs
> > >>        self.rhs = rhs
> > >>
> > >> class Form:
> > >>    def __eq__(self, other):
> > >>        return Equation(self, other)
> > >>
> > >> I understand there are other priorities, and others don't care as much
> > >> as I do about how fancy we can make the DOLFIN Python and C++ interface,
> > >> but I think this would make a nice final touch to the interface.
> > >
> > > I don't see value in it. In fact the opposite - it introduces complexity
> > > and a degree of ambiguity.
> >
> > Evidently, we all see things differently.
>
> Agree.
>
> > I fully support Anders in that
> > mathematical expressiveness is one of the key features of FEniCS, and I
> > think that without pushing these types of boundaries with regard to the
> > language, it will end up as yet another finite element library.
>
> Agree in general, but disagree for this particular case.
>
> > Could we compromise on having the two versions, one explicit (based on
> > LinearVariational[Problem|Solver] or something of the kind) and one
> > terse (based on solve(x == y)) ?
>
> Yes, I have seen this suggestion coming up :)
>
> I am really not that vigourus about this, as long as one can create intuitive
> objects which can be used and passed around, not just solved as in a script. I
> would like to here Martin's opinon to the UFL changes.

That's important. If Martin says no it's a showstopper, unfortunately.

> It really is one line versus two lines for a solve. Possible one line using
> linear_solve or nonlinear_solve, instead of just one solve.

No, it's about two lines versus either one or two lines depending on
your choice, just as for linear algebra. I very much like that it's
possible to do

  solve(A, x, b)

without needing to first create the solver object. Users who feel they
need to control lots of parameters for the solver, reuse
factorizations etc can choose to create the solver object. Not
everyone cares, and then it's good that one can just call solve(),
like backslash in MATLAB.

> Would introducing a namespace be a solution?
>
>   nonlinear::solve(...);
>   linear::solve(...);
>
> then:
>
>   using dolfin::[non]linear;
>   solve(...);

I think my suggestion is more elegant. :-)

--
Anders