dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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?

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

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

-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/