dolfin team mailing list archive

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

I think one of the main problems with VariationalProblem
is the mix of problem and solver in one abstraction.

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(?)

# The fully automated variant of solve:
solve(p) # Autodetect linear or nonlinear and differentiate if necessary

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

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