On Wed, Jun 15, 2011 at 01:45:04PM +0200, Martin Sandve Alnæs wrote: > On 15 June 2011 11:47, Garth N. Wells <gn...@cam.ac.uk> wrote: > > > > > > On 15/06/11 10:33, Martin Sandve Alnęs wrote: > >> On 15 June 2011 10:09, Garth N. Wells <gn...@cam.ac.uk> 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 <l...@simula.no> 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 <l...@simula.no> 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"<gn...@cam.ac.uk>: > >>>>>>>>>>>>>>> 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@lists.launchpad.net > >>>>>>> Unsubscribe : https://launchpad.net/~dolfin > >>>>>>> More help : https://help.launchpad.net/ListHelp > >>>>>>> > >>>>> > >>>> > >>>> _______________________________________________ > >>>> Mailing list: https://launchpad.net/~dolfin > >>>> Post to : dolfin@lists.launchpad.net > >>>> Unsubscribe : https://launchpad.net/~dolfin > >>>> More help : https://help.launchpad.net/ListHelp > >>> > >>> > >>> _______________________________________________ > >>> Mailing list: https://launchpad.net/~dolfin > >>> Post to : dolfin@lists.launchpad.net > >>> Unsubscribe : https://launchpad.net/~dolfin > >>> More help : https://help.launchpad.net/ListHelp > >>> > > > > > > _______________________________________________ > Mailing list: https://launchpad.net/~dolfin > Post to : dolfin@lists.launchpad.net > Unsubscribe : https://launchpad.net/~dolfin > More help : https://help.launchpad.net/ListHelp _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@lists.launchpad.net Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
