dolfin team mailing list archive

["Marie E. Rognes" <meg@xxxxxxxxx>]

On 21. juni 2011, at 09:58, Anders Logg <logg@xxxxxxxxx> wrote:

> Thanks for looking into this. So it looks like this will work out
> fine. So are there any objections to changing the VariationalProblem
> (which depends on the order of arguments) to the more explicit
> 
>  LinearVariationalSolver
>  NonlinearVariationalSolver
> 
> *solver* classes (instead of problem classes so that it's consistent
> with the way we treat la), and adding on top of this the free function
> solve() as a shortcut:
> 
>  solve(a == L, ...)
>  solve(F == 0, ...)
> 
> If not, I can get started on implementing this.
> 

Very much in favor :) 

Where does u enter in the free functions?

--
Marie


> --
> Anders
> 
> 
> On Tue, Jun 21, 2011 at 09:36:18AM +0200, Martin Sandve Alnæs wrote:
>> On 15 June 2011 13:06, Anders Logg <logg@xxxxxxxxx> wrote:
>>> Would this work? What it does is to define a simple Equation class
>>> that holds a pair of forms (lhs, rhs) but can still be used for
>>> comparisons.
>> 
>> Short story:
>> I'm pretty sure this will work out fine, except that the
>> Equation.__eq__ function must be dropped or changed
>> as explained below.
>> 
>> 
>> The details:
>> 
>>> class Equation:
>>> 
>>>    def __init__(self, lhs, rhs):
>>>        self.lhs = lhs
>>>        self.rhs = rhs
>> 
>> 
>> Allowing the syntax (a == L) == True has severe pitfalls,
>> and I see no use for it. It breaks the Python model.
>> In general, <expression> == True is not good python code,
>> as many boolean expressions are allowed to return non-bool
>> objects. That's why this works out in the first place, so we
>> should respect that :)
>> 
>> Anyway, we only need bool(a == L) to say "if a == L" and use
>> a and L in dicts and sets. So this would have to change:
>> 
>> (class Equation)
>>>    def __eq__(self, other):
>>>        if isinstance(other, bool):
>>>            return self.__nonzero__() == other
>>>        else:
>>>            return repr(self) == repr(other)
>> 
>> It can eventually be replaced with this:
>> 
>>   def __eq__(self, other):
>>       return isinstance(other, Equation) and\
>>           bool(self.lhs == other.lhs) and\
>>           bool(self.rhs == other.rhs)
>> 
>> 
>> As a generic implementation pattern, this should rather
>> use comparison functions from lhs and rhs:
>> 
>>>    def __nonzero__(self):
>>>        return repr(self.lhs) == repr(self.rhs)
>>    def __nonzero__(self):
>>        return self.lhs.eval_equals(self.rhs) # or another name, maybe
>> __cmp__ == 0
>> 
>> But for this particular case, it's probably ok since lhs and rhs are always
>> Form objects, and repr is used in that single __eq__ implementation today.
>> 
>> Note that I have only evaluated this for Form/Equation, my conclusion
>> that I guess this is ok does not extend to any other part of the universe...
>> E.g., if we were to use this trick for Expr as well, it could have severe
>> performance penalties, so that's not something to do in a hurry.
>> 
>> Martin
>> 
>> 
>>> class Form:
>>> 
>>>    def __init__(self, a):
>>>        self.a = a
>>> 
>>>    def __eq__(self, other):
>>>        return Equation(self, other)
>>> 
>>>    def __repr__(self):
>>>        return str(self.a)
>>> 
>>> f1 = Form(1)
>>> f2 = Form(2)
>>> f3 = Form(1)
>>> 
>>> eq = f1 == f2
>>> print eq == False
>>> print eq == True
>>> 
>>> if f1 == f2:
>>>    print "f1 == f2 is True"
>>> 
>>> if f1 == f3:
>>>    print "f1 == f3 is True"
>>> 
>>> 
>>> 
>>> On Tue, Jun 14, 2011 at 10:01:47PM +0200, Martin Sandve Alnæs wrote:
>>>> I agree it is pretty in a way, but here comes the showstopper... Sorry! ;)
>>>> 
>>>> There can be no ufl.Form.__eq__ implementation
>>>> that does not conform to the Python conventions
>>>> of actually comparing objects and returning
>>>> True or False, because that will break usage
>>>> of Form objects in built in Python data structures.
>>>> 
>>>> Martin
>>>> 
>>>> On 14 June 2011 21:53, Anders Logg <logg@xxxxxxxxx> 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).
>>>>> 
>>>>> 
>>>>> _______________________________________________
>>>>> Mailing list: https://launchpad.net/~dolfin
>>>>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>>>> Unsubscribe : https://launchpad.net/~dolfin
>>>>> More help   : https://help.launchpad.net/ListHelp
>>>>> 
>>>