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. 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. >>>>> >>>>> -- >>>>> 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 <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
