dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

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

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.

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

Garth

> --
> Anders