dolfin team mailing list archive

[Martin Sandve Alnæs <martinal@xxxxxxxxx>]

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.


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