← Back to team overview

dolfin team mailing list archive

Re: VariationalProblem interface

 

On Wed, Jun 15, 2011 at 01:50:07PM +0200, Martin Sandve Alnæs wrote:
> Not sure, will have to get back to that.

ok.

> But wether we can use (a==L), (J==-F), (F==0) or just (a, L), (J, -F), (F, 0)
> is basically a purely aestetic difference with no functional difference.
> So the rest of the discussion can continue independently of this.

Yes.

--
Anders


> Martin
>
> 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.
> >
> > class Equation:
> >
> >    def __init__(self, lhs, rhs):
> >        self.lhs = lhs
> >        self.rhs = rhs
> >
> >    def __eq__(self, other):
> >        if isinstance(other, bool):
> >            return self.__nonzero__() == other
> >        else:
> >            return repr(self) == repr(other)
> >
> >    def __nonzero__(self):
> >        return repr(self.lhs) == repr(self.rhs)
> >
> > 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
> >> >
> >


Follow ups

References