dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Monday June 13 2011 08:09:36 Anders Logg wrote:
> On Mon, Jun 13, 2011 at 05:04:53PM +0200, Kristian Ølgaard wrote:
> > On 13 June 2011 16:37, Anders Logg <logg@xxxxxxxxx> wrote:
> > > On Mon, Jun 13, 2011 at 03:30:22PM +0100, Garth N. Wells wrote:
> > >> 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.
> > > 
> > > I agree with this. It's important that we can check input arguments.
> > > 
> > > What are the votes in favor of either of these two options?
> > > 
> > > 1. One class VariationalProblem for all (non)linear problems.
> > > 
> > > 2. Two classes LinearProblem and NonlinearProblem.

2 gets my vote.

> > I vote for 2 with the arguments of readability/debugging as Martin and
> > Garth have promoted,
> > and like Garth I also don't think that one class should try to do
> > everything.
> 
> ok. Note the new choice of names LinearProblem/NonlinearProblem
> (without Variational).
> 
> This would also let us keep the VariationalProblem class and issue a
> suitable error message that it has been replaced by those other
> classes.

Sounds good!

Johan
 
> --
> Anders
> 
> > Kristian
> > 
> > >> >>>>> 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.
> > > 
> > > _______________________________________________
> > > 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