dolfin team mailing list archive

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

On 15 June 2011 11:47, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
>
>
> 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.

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