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