dolfin team mailing list archive

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

On 13/06/11 14:54, Anders Logg wrote:
> On Mon, Jun 13, 2011 at 02:42:57PM +0100, Garth N. Wells wrote:
>>
>>
>> On 13/06/11 14:10, Anders Logg wrote:
>>> On Mon, Jun 13, 2011 at 01:09:20PM +0100, Garth N. Wells wrote:
>>>>
>>>>
>>>> On 13/06/11 12:51, Marie E. Rognes wrote:
>>>>> On 06/13/2011 10:03 AM, Garth N. Wells wrote:
>>>>>>
>>>>>>
>>>>>> On 12/06/11 22:54, Anders Logg wrote:
>>>>>>> On Wed, Jun 08, 2011 at 11:47:10PM +0200, Anders Logg wrote:
>>>>>>>> I'm on a cell phone and can't engage in any discussions but I just
>>>>>>>> want to throw something in. Marie and I discussed at some point
>>>>>>>> another option which is to write everything as F = 0 and let UFL
>>>>>>>> compute J even for linear problems. J would be a optional
>>>>>>>> variable. Then we could have a common interface and also a common
>>>>>>>> solver.
>>>>>>>
>>>>>>
>>>>>> I think that this would be a nice option, but I don't think that it can
>>>>>> work without a more elaborate interface than just
>>>>>>
>>>>>>    pde = VariationalProblem(F, bcs)
>>>>>>
>>>>>> because UFL cannot know which coefficient in a form is unknown, and
>>>>>> therefore which function to compute the linearisation with respect to.
>>>>>> Moreover,
>>>>>>
>>>>>>   - UFL cannot compute the linearisation for all forms
>>>>>>   - In some cases it's not desirable to let UFL do the linearisastion
>>>>>>
>>>>>
>>>>> Some specific responses: First, the moment solve is called (with a
>>>>> coefficient), which function to compute the linearization with respect
>>>>> to, is known.
>>>>
>>>> Which means that the simple call
>>>>
>>>>   u = pde.solve()
>>>>
>>>> will not be possible.
>>>
>>> Why do you call it "pde" and not "problem"?
>>
>> The name has nothing to do with the discussion.
> 
> This discussion is about just that: which names to use.
> 

Yes, the class name, but not what I dream up for an object:

  xxx = FooVariationalFoo(....)

  u = xxx.solve()

>> I could call it xxx.
>>> Once upon a time it was
>>> called that, so yet another option would be
>>>
>>>   LinearPDE
>>>   NonlinearPDE
>>>
>>> which is short enough.
>>>
>>> What speaks against this, and likely one of the reasons we switched to
>>> VariationProblem is that the above looks weird when used to define a
>>> projection which does not involve any derivatives.
>>>
>>> And here's yet another option:
>>>
>>>   LinearProblem
>>>   NonlinearProblem
>>>
>>
>> This looks ok, but it's a secondary issue for if we agree that we should
>> have separate classes for linear and nonlinear equations (or have we
>> agreed this?).
> 
> Not yet. I think we should settle first one whether we want
> 
> 1. a common class for (non)linear problems
> 2. two classes
> 

Agree. I strongly support 2.

What I oppose vehemently is distinguishing between linear and nonlinear
by the order of the lhs/rhs arguments. As pointed out by Martin, the line

   VariationalProblem(C, D)

says nothing to the reader about the type of problem, and it has led to
confusion.

The code

   LinearVariationalProblem(C, D)

or

   NonlinearVariationalProblem(C, D)

(class names subject to change) are clear, and the code can test for the
arity of C and D and throw an error if the order is mixed up.

>>>>> Second, the Jacobian could definitely be provided as an
>>>>> optional argument for those cases where it is not feasible/desirable
>>>>> for UFL to do the linearization. See below for longer, more general
>>>>> rant.
>>>>>
>>>>>>> Any further thoughts on this? It's important we get this right.
>>>>>>>
>>>>>>> I can live with
>>>>>>>
>>>>>>>    LinearVariationalProblem
>>>>>>>    NonlinearVariationalProblem
>>>>>>>
>>>>>>> if everyone else thinks that's the best option,
>>>>>
>>>>> I don't think it is the best option.
>>>>>
>>>>>>> but don't really like
>>>>>>> it since it's very long.
>>>>>
>>>>> My reasons do not really relate to the length of the name, but rather
>>>>> that I do not see why we should make the difference between linear and
>>>>> nonlinear variational problems greater in the _highest level
>>>>> interface_ of DOLFIN. I think this just makes things more cumbersome
>>>>> (a) in terms of debugging and (b) in terms of explaining how to solve
>>>>> variational problems to FEniCS beginners.
>>>>>
>>>>
>>>> My experience is the complete opposite of this.
>>>
>>> I agree with Marie here. It's a bit weird when we do
>>>
>>>   u = TrialFunction(V)
>>>   ...
>>>   u = Function(V)
>>>   problem.solve(u) # or u = problem.solve()
>>>
>>
>> I agree that this is not elegant, but I think that it's a secondary
>> issue (but still a something to improve).
> 
> I don't think this is an issue we can fix in small increments. We need
> to find out what kind of interface we want, then implement it and then
> not touch it ever again.
>

The above is part of the Python interface. If the C++ interface is
resolved, we can turn our attention to prettifying Python issues, like
the above.

>>> It's almost as confusing as the V=V thing we had for a while.
>>>
>>>>> My typical pedagogical use case looks something like this (in Python):
>>>>> Say you have explained to someone how to solve stationary Stokes with
>>>>> Taylor-Hood, and then want to extend to Navier-Stokes. Then ideally
>>>>> all that should be necessary is to add the nonlinear term to the
>>>>> form(s). But no. In addition, you at least have to replace the
>>>>> TrialFunction with a Function (which typically leads to a longer
>>>>> discussion on what a TrialFunction represents and why that is used for
>>>>> linear problems when the same ansatz is made for the two in most text
>>>>> books), and then you have to call solve with the same Function instead
>>>>> of just returning it (which is also on my top 3 list of pitfalls). I
>>>>> don't see how adding
>>>>> LinearVariationalProblem/NonlinearVariationalProblem to the list makes
>>>>> it better.  The same (but the other way) goes when wanting to debug
>>>>> Navier-Stokes by reducing to Stokes.
>>>>>
>>>>> In essence, I've started treating all variational problems as
>>>>> nonlinear.
>>>>>
>>>>
>>>> This might be fine for toy problems, but not for real problems. There is
>>>> overhead in a nonlinear solve to check for convergence, which is not
>>>> required if the eqn is linear.
>>>
>>> I don't see why this should be restricted to toy problems. Defining
>>> the Jacobian is always an option and then it would be just like today.
>>>
>>
>> Because if I'm solving a very large linear problem, I'm going to be
>> annoyed if the code insists on unnecessarily assembling the residual.
> 
> Aren't all "real world problems" nonlinear? ;-)
> 

Fortunately no.

>>> Perhaps it would be possible for UFL to check that a problem is linear
>>> (by differentiating twice and getting zero) and then we could remove
>>> the overhead of the extra residual evaluation.
>>>
>>
>> Too much magic over simplicity for my liking. We also need to bear in
>> mind that we're talking about the C++ interface too. I think a better
>> starting point is the C++ interface. The Python interface will follow
>> naturally, and sugar can be added later.
> 
> I agree with that. Say for a moment that we choose to implement the
> common (non)linear interface. Then what should the interface look
> like in C++?
> 
> 1. Poisson::BilinearForm a(V, V);
>    Poisson::LinearForm L(V);
> 
> 2. Poisson::Residual F(V);
>    Poisson::Jacobian J(V, V);
> 
> 3. Poisson::VariationalProblem F(V, V);
> 
>> What if the equation is linear, but I want to solve it as a nonlinear
>> equation? I do this all the time as a first step for testing when
>> testing a nonlinear solver? I know we can always 'add another option' to
>> handle various case, but it all add complexity.
> 
> Isn't that an argument for having a common nonlinear solver for all
> problems?
> 

No. It's an argument that we shouldn't try to make one class do
everything, and that the code should be explicit.

>> Another factor is that there is more to solving nonlinear equations than
>> vanilla versions of Newton's method. It would nice in the future to have
>> a NonlinearFoo that does more than what we now, and without increasing
>> the code complexity for a straightforward linear eqn solve.
> 
> Perhaps, but one could also argue that is yet another argument in
> favor of a single nonlinear solver. There is some code duplication
> already in the two linear/nonlinear solvers, setting options etc.
>

Ditto what I wrote above. Trying to make one class do too much leads to
complications.

Garth

> --
> Anders