dolfin team mailing list archive

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

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