dolfin team mailing list archive

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

On 14/06/11 20:53, Anders Logg wrote:
> On Tue, Jun 14, 2011 at 09:03:31PM +0200, 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.
> 
> That's true, but some magic is usually required to make pretty.
> 
> Being able to write dot(grad(u), grad(v))*dx is also a bit magic.
> The step from there to solve(a == L) is short.
> 
>>>>> 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.
> 
> Yes, of course it's a stretch. It's not very ugly, but enough to
> bother me.
> 
>>>>> 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.
> 
> Complexity yes (but not much, it would require say around 50-100
> additional lines of code that I will gladly contribute), but I don't
> think it's ambiguous. We could perform very rigorous and helpful
> checks on the input arguments.
> 
>> Evidently, we all see things differently. 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.
>>
>> 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)) ?
> 
> That's what I'm suggesting. The solve(x == y) would just rely on the
> more "explicit" version and do
> 
>   <lots of checks>
>   LinearVariationalSolver solver(x, y, ...);
>   solver.solve(u);
> 
> So in essence what I'm asking for is please let me add that tiny layer
> on top of what we already have + remove the Problem classes (replaced
> by the Solver classes).
>

It seems that we (almost?) all agree to add to the C++ side:

  LinearVariationalProblem

  NonlinearVariationalProblem

I think that 'FooVariationalProblem' is probably a better name than
'FooVariationalSolver'. Could 'variational' (minimisation) and 'solver'
imply some abstract linear algebra problem? Most accurate is probably
'FooVariationalProblemSolver', but it's a bit long.

On the Python side, I don't mind others trying something fancy, but I'll
reserve my judgement ;).

Garth


> --
> Anders
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp