As some of you may know, I am now working at NASA's Johnson Space Center doing CFD for hypersonic aerodynamic applications.  The network setup is a little backwards around here, so it will be a couple of days before I can check in/out at work.  For some reason they restrict where I ssh.

Anyway, I checked in a lot of changes last week when I got back from North Carolina.  Now, most of these were trivial (changing Copyright 2002-2003 to 2002-2004) and don't warrant discussing, but the new solver stuff is a little more interesting.

The basic idea is this:

1.) Each simulation will have one EquationSystems class.  This can contain multiple Systems of various types. 
    (It might be worth making the EquationSystems a singleton?)

2.) Calling EquationSystems::solve() simply calls solve on each of the systems contained in it.  This goes for init(), reinit(), etc...  Therefore the EquationSystems interface is quite similar to that of each individual system.

3.) The new interface maintains backward compatibility.  However, we can now do some more interesting stuff, too.  In the past, all systems got a matrix, which is super-stupid and wasteful for field data or post-processed quantities.  Now the idea is that something like this will work:

EquationSystems es(mesh);

es.add_system<ImplicitSystem>                    ("S1"); // Create a standard implicit system. solve() solves Ku=f
es.add_system<Nonlinear<ImplicitSystem> > ("S2"); // Create a nonlinear implicit system. solve() solves K(u)u=f(u)
es.add_system<ExplicitSystem>                   ("S3"); // An explicit system, no matrix.

4.) Now, the solvers can operate on an individual system *or* the entire EquationSystems. Above "Nonlinear<>" is a solver acting on an individual system.  But consider this example:

Transient<> solver (es);


Here "Transient<>" is a solver acting on an entire EquationSystems object. In this case, the solver will execute a timestep loop and call es.solve() at each time step.  Now, suppose instead I want a transient loop with adaptivity at each time step?

No problem:

Transient<Adaptive<> > solver(es);


Oh, what about a transient loop with adaptivity at each timestep, and a nonlinear loop inside each adaptive loop? Well,

Transient<Adaptive<Nonlinear<> > > solver(es);


This is where it is all going.  The interface needs a little more work...  The user should be able to add pre and post-processing functions that get called at specified increments to do things like write the solution at each time step, etc...  I just wanted to check this in so that we can get more people interested in it. 

Steffen:  I think that in this framework we can easily add an EigenSolver, and then have:

Eigen<> solver (es);