dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Jun 15, 2011 at 01:45:04PM +0200, Martin Sandve Alnæs wrote:
> On 15 June 2011 11:47, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
> >
> >
> > On 15/06/11 10:33, Martin Sandve Alnęs wrote:
> >> On 15 June 2011 10:09, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
> >>>
> >>>
> >>> On 14/06/11 22:02, Martin Sandve Alnęs wrote:
> >>>> I think one of the main problems with VariationalProblem
> >>>> is the mix of problem and solver in one abstraction.
> >>>>
> >>>
> >>> This point may be worth discussing. It depends on how the 'Problem' ius
> >>> defined. This may or may not include the solution method.
> >>
> >> This is a rather crucial point to decide on before the rest of the
> >> design can be finalized.
> >>
> >
> > I gave his some a while ago, and came to the weak conclusion that the
> > 'Problem' should include the solver, otherwise the Problem is nothing
> > more than a tuple
> >
> >  (linear form, bilinear form, bcs, linear=true)
> >
> > and it seemed like overkill to introduce a class for this.
>
> For my suggested VariationalProblem variant I would argue
> that the unknown function is also a natural part of the problem:
>
> Find u such that // Not 'find f', not 'find c', so 'u' is actually a
> key part of the problem.
> (c * grad u, grad v) = (f, v) in Omega // pde
> u = g on Gamma // bcs
>
> and having a class to gather five variables doesn't seem so strange to me.
> At least not if that class represents a natural abstraction, and can be passed
> around many places (many solvers), since it then will save work and simplify
> interfaces.
>
> Also seems to me that
> mysolver1.solve(problem);
> mysolver2.solve(problem);
> is more flexible w.r.t. solver experimentation.

I agree with all of this, but I'm still attracted by the idea of
adding the solve(a == L) shortcut on top of this.

--
Anders


> Martin
>
>
> >>>> Find u in V such that
> >>>>   F = 0 in Omega // pde
> >>>>   u = g on Gamma // bcs
> >>>>
> >>>> Here's my suggested variation of VariationalProblem:
> >>>>
> >>>> # Generic form:
> >>>> problem = VariationalProblem(unknown, lhsform, rhsform, bcs)
> >>>> solve(problem)
> >>>>
> >>>> # Examples and variations follow, all consistent with the above:
> >>>>
> >>>> # Typical linear system
> >>>> p = VariationalProblem(u, a, L, bcs)
> >>>> # Note that ordering and sign is consistent with a linear system
> >>>> p = VariationalProblem(u, J, -F, bcs)
> >>>> # Possibly(!) nonlinear system F(u) = 0
> >>>> p = VariationalProblem(u, F, 0, bcs)
> >>>>
> >>>> # Note the 0 which maybe makes this less ambiguous(?)
> >>>>
> >>>
> >>> This brings us almost full circle to where we are now.
> >>
> >> I see the closeness to the current design as
> >> an advantage, requiring fairly small changes
> >> to existing code.
> >>
> >>> If we prefix
> >>> 'Linear' and 'Nonlinear' we use a descriptive word rather than just '0'.
> >>
> >> Not at all. The 0 does not mean nonlinear, the problems defined
> >> above are neither specified as linear nor nonlinear. This will allow
> >> passing the same problem to a linear solver or a nonlinear solver,
> >> like you requested (passing a linear problem to a nonlinear solver
> >> for testing). This is what can be achieved by separating problem
> >> and solver.
> >>
> >>> Having a free function 'solve' removes the potential for more powerful
> >>> solvers since it removes the possibility for data structure and solver
> >>> re-use and fine tuning.
> >>
> >> Not at all. The existence of a function solve does not remove anything
> >> from anywhere.
> >>
> >> Since a problem is independent of solution method, it can be passed
> >> to a user defined solver. You could even tune your solver and apply it
> >> to multiple problems, or define one problem and pass it to multiple solvers
> >> for comparing.
> >>
> >> Also, if you read the pseudocode for solve, you'll see I attempted
> >> to allow a userdefined solver to be passed, while reusing the fancy
> >> (Python-only) features such as differentiation.
> >>
> >>
> >>>> # The fully automated variant of solve:
> >>>> solve(p) # Autodetect linear or nonlinear and differentiate if necessary
> >>>>
> >>>
> >>> In Python. It is a nice feature that we should add, but it doesn't
> >>> impact on the C++ side.
> >>
> >> True. But not necessary if you specify that
> >> you want to use a linear or nonlinear solver.
> >> And the F=0 formulation will be useful for
> >> nonlinear methods that doesn't need J.
> >>
> >>>> # With optional behaviour changes:
> >>>> solve(p, nonlinear=True) # Force nonlinear solver
> >>>> solve(p, nonlinear=False) # Force linear solver
> >>>> solve(p, u=u2) # place solution in u2 instead of p.u, useful for linear solvers
> >>>> solve(p, solver=mysolver) # Use nondefault solver (need solver API)
> >>>>
> >>>
> >>> I think that there is a consensus to not use "linear/nonlinear =
> >>> true/false" arguments. We had this for a long time and then removed it.
> >>
> >> I'm fine with any other way to specify the choice, as long as it is
> >> not the order of arguments :)
> >>
> >> mysolver.solve(problem) # apply arbitrary solver object
> >> linear_solve(problem) # apply default linear solver
> >> nonlinear_solve(problem) # apply default nonlinear solver
> >> solve(problem) # automatic detection in python, safe default in C++ if
> >> we cannot automate that
> >>
> >> Martin
> >>
> >>> Garth
> >>>
> >>>> # Pseudocode for solve, showing it's feasible:
> >>>> def solve(p, nonlinear=None, solver="auto"):
> >>>>     # Allow alternative solvers in some way
> >>>>     if solver != "auto":
> >>>>         check properties of solver and adjust what is necessary below
> >>>>
> >>>>     # Differentiate if necessary (skip if we use a solver that doesn't need it):
> >>>>     if isinstance(p.rhs, int) and p.rhs == 0:
> >>>>         lhs, rhs = derivative(p.lhs, p.u), -p.lhs
> >>>>     else:
> >>>>         lhs, rhs = p.lhs, p.rhs
> >>>>
> >>>>     # Detect nonlinear or linear equation unless forced
> >>>>     if nonlinear is None: # Set to True or False to force this
> >>>>         # I have this function in my private sandbox
> >>>>         nonlinear = ufl.algorithms.depends_on(lhs, p.u)
> >>>>
> >>>>     # Solve with linear or nonlinear default solver
> >>>>     if solver == "auto":
> >>>>         if nonlinear:
> >>>>             return solve_nonlinear(p.u, lhs, rhs, p.bcs)
> >>>>         else:
> >>>>             return solve_linear(p.u, lhs, rhs, p.bcs)
> >>>>     else:
> >>>>         return solver.solve(...)
> >>>>
> >>>>
> >>>> Good properties of this solution:
> >>>> - Makes the VariationalProblem independent of the solver method
> >>>> - Allows explicit or automatic choice of nonlinear/linear solver
> >>>> - Consistent ordering of forms given to VariationalProblem
> >>>> - Binds u to the problem, which is in accordance with (*) and allows
> >>>> differentiation
> >>>> There may be something I've forgotten of course, so speak up.
> >>>>
> >>>> The binding of u to the problem is not really necessary in
> >>>> the case of a linear problem, but makes things consistent.
> >>>> With solve(p, u=u2) we still allow reuse of the VariationalProblem.
> >>>>
> >>>> Martin
> >>>>
> >>>> On 14 June 2011 22:09, Anders Logg <logg@xxxxxxxxx> wrote:
> >>>>> That's too bad... I will think of an alternate route.
> >>>>>
> >>>>>
> >>>>>
> >>>>> 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
> >>>>>>>
> >>>>>
> >>>>
> >>>> _______________________________________________
> >>>> Mailing list: https://launchpad.net/~dolfin
> >>>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> >>>> Unsubscribe : https://launchpad.net/~dolfin
> >>>> More help   : https://help.launchpad.net/ListHelp
> >>>
> >>>
> >>> _______________________________________________
> >>> Mailing list: https://launchpad.net/~dolfin
> >>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> >>> Unsubscribe : https://launchpad.net/~dolfin
> >>> More help   : https://help.launchpad.net/ListHelp
> >>>
> >
> >
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp