dolfin team mailing list archive

["Garth N. Wells" <g.n.wells@xxxxxxxxxx>]

On Mon, 2006-02-27 at 08:32 -0600, Anders Logg wrote:
> Garth, I think the nonlinear solver and the nonlinear PDE class look
> very good now.
> 
> I just have one remaining question. It doesn't seem like F() and J()
> are used by the Newton solver, only form(). Is that right?
> 

They're not used yet. They will be used later when the Newton solver can
be tuned. For example, a modified Newton method might call J() only at
the first iteration and F() at every iteration. The extra functions
could be removed until I get around to adding functionality to
NewtonSolver.

> Maybe one can put something in NonlinearProblem::form() that calls the
> user-supplied F() and J() functions if form() has not been
> implemented? The NewtonSolver would thus call form() and which in turn
> calls F() and J(), unless form() has been overloaded by the user
> (which would be more efficient).
> 

Sure.  I'll take a look at it.

Any ideas about NonlinearPDE having NewtonSolver as a private variable?
It's possible that it could be removed once NonlinearPDE has more
sophisticated solve options, like doing Newton-Raphson itself rather
than single Newton solves as it does now.

Garth   

> /Anders
> 
> On Fri, Feb 24, 2006 at 05:49:09PM +0100, Garth N. Wells wrote:
> > On Fri, 2006-02-24 at 09:47 -0600, Anders Logg wrote:
> > > I haven't looked much at the design of the nonlinear solver before,
> > > but I've looked at it some now. I would suggest some minor changes.
> > > Here are some questions/suggestions:
> > > 
> > > 1. The main class should be NewtonSolver and it should solve
> > >    a nonlinear system F(x) = 0. The class should be in src/kernel/nls.
> > > 
> > This has been cleaned up
> > 
> > > 2. The main function should be
> > > 
> > >    // Solve F(x) = 0
> > >    uint NewtonSolver::solve(NonlinearFunction& F, Vector& x);
> > > 
> > It is now.
> > 
> > > 3. The class NewtonSolver should not know anything about PDE, bilinear
> > >    forms, variational problems etc, it just solves the discrete
> > >    system F(x) = 0.
> > > 
> > It doesn't now.
> > 
> > > 4. All classes in src/kernel/nls only handle the system F(x) = 0,
> > >    no PDEs etc. Looks like we only need NewtonSolver and
> > >    NonlinerFunction right now.
> > > 
> > Yes. NonlinearFunction has changed to NonlinearProblem.
> > 
> > > 5. Maybe NonlinearFunction::form() is not needed? NonlinearFunction
> > >    should only need F() and J().
> > > 
> > 
> > I would prefer to keep NonlinearFunction::form(). It will generally be
> > more efficient to compute the Jacobian and the RHS together,
> > particularly when complex nonlinear terms for both the RHS and the
> > Jacobian are computed together in a separate procedure (such as in
> > plasticity, for example).
> > 
> > > 6. Why does NewtonSolver inherit from KrylovSolver?
> > > 
> > 
> > This was to allow details of the linear solver to be set. NewtonSolver
> > doesn't inherit from KrylovSolver any more. NewtonSolver now has a
> > private member KrylovSolver, so that it's paramters can be set, and it
> > is not recreated at each solve.
> > 
> > 
> > > 7. Move all PDE stuff to src/kernel/pde. It's enough to have form
> > >    stuff in src/kernel/form.
> > > 
> > Sure. I was planning to do this.
> > 
> > > 8. NonlinearPDE knows about forms etc and it uses NewtonSolver to
> > >    solve the nonlinear PDE. I don't think newton_solver needs to
> > >    be a member variable. It's probably enough to create a NewtonSolver
> > >    in the NonlinearPDE::solve() function, same as a GMRES object is
> > >    created in LinearPDE::solve().
> > > 
> > 
> > The reason I didn't do this is that the NewtonSolver typically re-used
> > many times 
> > for a given PDE (see src/demo/pde/nonlinear/nonlinear-possion). The
> > Jacobian matrix and RHS vector are private varaibles of NewtonSolver,
> > and newton_solver is a member variable of NonlinearPDE, so that a matrix
> > for the Jacobian, and the vectors for the RHS and the incremental
> > solution are not re-created at each step of a nonlinear solve. 
> > 
> > 
> > 
> > 
> > > 9. NonlinearPDE needs to create an object of type NonlinearFunction
> > >    that it can pass to NewtonSolver::solve(). This can be done by
> > >    creating a sub class of NonlinearFunction that represents the
> > >    nonlinear function for a given BilinearForm and LinearForm. The sub
> > >    class can be called for example NonlinearPDEFunction and be put in
> > >    a separate file in src/kernel/pde.
> > > 
> > >    The constructor of NonlinearPDEFunction takes a NonlinearPDE as
> > >    an argument, so in NonlinearPDE::solve(), we would do something
> > >    like
> > > 
> > >    NonlinearPDEFunction F(this);
> > >    NewtonSolver newton;
> > >    newton.solve(F, x);
> > > 
> > > So all in all, we need the following classes:
> > > 
> > >    NewtonSolver         - solves F(x) = 0
> > >    NonlinearFunction    - represents the function F
> > > 
> > >    NonlinearPDE         - represents a nonlinear PDE
> > >    NonlinearPDEFunction - represents F for a nonlinear PDE
> > > 
> > 
> > OK
> > 
> > Garth
> > 
> > > /Anders
> > > 
> > > 
> > > On Fri, Feb 24, 2006 at 02:29:36PM +0100, Garth N. Wells wrote:
> > > > On Fri, 2006-02-24 at 14:06 +0100, Johan Jansson wrote:
> > > > > On Fri, Feb 24, 2006 at 01:45:33PM +0100, Garth N. Wells wrote:
> > > > > > Hi,
> > > > > > 
> > > > > > I'm working on integrating NonlinearPDE into PDE, and it's coming along
> > > > > > well. I do have a question about the correntness/style of some code.
> > > > > > I've sketched out the code below. The class NonlinearPDE has a member
> > > > > > object of type NewtonSolver (called newton_solver). When asking
> > > > > > newton_solver to solve the problem, I pass nonlinear_pde to it (because
> > > > > > nonlinear_pde knows how to form it's Jacobian and RHS vector). Is this a
> > > > > > problem? I could tidy it up by making a another class - is there a
> > > > > > problem with embedded classes (I recall something to so with swig)?
> > > > > > 
> > > > > > Garth
> > > > > > 
> > > > > > class NonlinearPDE : public GenericPDE
> > > > > > {
> > > > > > public:
> > > > > >     .
> > > > > >     .
> > > > > >     void solve(NonlinearPDE
> > > > > > 
> > > > > >     // Form Jacobian and RHS (x is current solution)
> > > > > >     void form(Matrix& A, Vector& b, const Vector& x)
> > > > > > 
> > > > > >     void solve(Function &u) 
> > > > > >     {     
> > > > > >       newton_solver.solve(*this, u);
> > > > > >     }
> > > > > > 
> > > > > > private:
> > > > > > 
> > > > > >   NewtonSolver newton_solver;
> > > > > >  
> > > > > > }
> > > > > > 
> > > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > > Hi!
> > > > > 
> > > > > What you are doing is perfectly all right. NonlinearPDE should also
> > > > > inherit from NonlinearProblem if I understand correctly:
> > > > > 
> > > > > class NonlinearPDE : public GenericPDE, public NonlinearProblem
> > > > > 
> > > > > 
> > > > 
> > > > Exactly. I'm using the name NonlinearFunction at the moment, but
> > > > NonlinearProblem is a better name. An even better name would not refer
> > > > to "nonlinear". Any ideas for names for a base class for objects which
> > > > simply know how to form a matrix and a vector for a given problem?
> > > > 
> > > > Garth 
> > > > 
> > > > 
> > > > > 
> > > > > SWIG has problems parsing nested class definitions:
> > > > > 
> > > > > class A
> > > > > {
> > > > >   ...
> > > > > 
> > > > >   class B
> > > > >   {
> > > > >     ...
> > > > >   }
> > > > > 
> > > > >   ...
> > > > > }
> > > > > 
> > > > > This is usually not a problem though, SWIG cannot see the interface of
> > > > > B, but it can see the interface of A, and that's usually enough. We
> > > > > shouldn't worry about this, it's a deficiency of SWIG which will
> > > > > eventually be fixed.
> > > > > 
> > > > >   Johan
> > > > 
> > > > 
> > >