dolfin team mailing list archive

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

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

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

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

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

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.

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.

Garth


> --
> Anders