dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

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.

Garth

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