On Fri, 2 Oct 2009, Derek Gaston wrote:

> Just a quick observation that this could continue to get extended... for 
> instance "we know what a diffusion operator looks like... we can assemble 
> that for you!".... which starts to look a lot like the system I've grown here 
> at INL....

You didn't number this, but let's call it 3:

3. The ability to add additional operators to a System should be
generalized.  Then "add a spatially dependent forcing function" would
just be a special case of that.

> 1.  Speed.  As soon as you do this you're going to introduce another set of 
> element loops, FE reinits, quadrature point location recalculations, and 
> possibly variable calculations (for coupled manufactured solutions?).

This isn't a problem for us on verification runs, but it interacts with
issue 3, in that a general system that adds redundant loops would have
an overhead that for some codes would be unacceptable for production
runs.

C++ makes it possible to do generalized compile-time dial-an-operator
without redundant loops, recalculations, or even much heavy parsing in
the inner loop, but the necessary template tricks may not look pretty.
I'd want to experiment with creating such an API quite a bit before
considering adding anything like that to the library, and it would be
difficult to apply to existing code.

I think the fundamental problem may be that we don't have a clear
separation between the System (which handles allocation of variables,
vectors, matrices) and the Physics (which would express the equations
that are evaluated to find residual and jacobian components).
FEMSystem gives a little more abstraction to physics in libMesh
applications, but it doesn't separate this intertwining.  If people
could create two Physics objects that operated on the same System's
data, then we would never consider adding a new method like
add_forcing_function().  Just call the constructor 
SumPhysics(Physics &A, Physics &B), where B is a special case that
assembles your forcing function.

> 2.  Speaking of coupling... how are you going to handle that (but maybe 
> people don't do coupled manufactured solutions... I don't know).

We're planning on doing things fully coupled, so (other than making
sure that f(x) can be called by variable name) there are no worries
there.  For people doing loosely coupled solves, it would be necessary
for them to override EquationSystems::solve() to call each system in
the right order and to loop until convergence.  That's the sort of
code everyone tends to just shove in main() now, but EquationSystems
is designed to be subclassable for just that reason.  (also for use
with UniformErrorEstimator).

How to do the API for transient problems properly (for the adjoint
methods, too) is another similar question.  Let's call that "4." Still
thinking about that one.
---
Roy

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay 
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Libmesh-devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-devel

Reply via email to