Re: [Libmesh-users] New Solvers & update...

[<lu...@gi...>]

Ben

Your proposed new interface sounds extremely cool. 

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

Do you want to restrict the equation systems in EquationSystems to one
single mesh. For multi-physics applications it would be useful to be
able to have different mesh refinements for different systems
(e.g. heat transfer and fluid flow). One can even think of problems,
where one domain (e.g. fluid flow) should be a subdomain of the other
(e.g. heat transfer).

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

Is the user still able do decide whether to call es.solve(), or to
call solve() for each systems individually?

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

Here I would like to be able to say either es(mesh) for a shared mesh,
or 

es.system("S1").add_mesh(mesh1)
es.system("S2").add_mesh(mesh2)

Which such an approach, some quantities would have to be mapped from
mesh to mesh. OTOH one would be able to define different (overlapping)
problem domains and refine the each mesh based on different error
estimators.

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

cool.

 > 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);
 > 
 > solver.solve();

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

Where is the adaptivity happening? On the global mesh object? What
happens if two solvers are declared adaptive?

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

This is a very good point. I think that the calculation of
intermediate state variables and the definition of post-variables is
the same. It would be nice to obtain these quantitites remapped on the
current (refined) grid.

Thanks for sharing your thoughts 

Martin