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. > A problem is usually formulated: > > Find u in V such that > a = L in Omega // pde > u = g on Gamma // bcs > > 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. If we prefix 'Linear' and 'Nonlinear' we use a descriptive word rather than just '0'. 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. > # 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. > # 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. 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 <[email protected]> wrote: >> That's too bad... I will think of an alternate route. >> >> -- >> Anders >> >> >> 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 <[email protected]> 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"<[email protected]>: >>>>>>>>>>>> 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 : [email protected] >>>> Unsubscribe : https://launchpad.net/~dolfin >>>> More help : https://help.launchpad.net/ListHelp >>>> >> > > _______________________________________________ > Mailing list: https://launchpad.net/~dolfin > Post to : [email protected] > Unsubscribe : https://launchpad.net/~dolfin > More help : https://help.launchpad.net/ListHelp _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
