dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Wednesday October 20 2010 14:27:55 Marie Rognes wrote:
> On 20. okt. 2010 23:14, Johan Hake wrote:
> > On Wednesday October 20 2010 13:24:17 Marie Rognes wrote:
> >> On 20. okt. 2010 21:51, Johan Hake wrote:
> >>>> Progress has been made towards the "real" implementation (aka in the
> >>>> c++ core) of automated goal-oriented error control and adaptivity (cf
> >>>> email thread "Fwd: [Branch ~dolfin-core/dolfin/main] Rev 4734:
> >>>> Introduce functionality for automated error control and adaptivity"
> >>>> for background info.)
> >>>> 
> >>>> Summary:
> >>>> ----------------
> >>>> 
> >>>> Marie suggests the following interface:
> >>>> 
> >>>> Same as for standard VariationalProblems, but with
> >>>> 
> >>>>      pde.solve(u, tol, M)
> >>>> 
> >>>> where "tol" should be a (positive) real number and M a Form,
> >>>> representing the quantity of interest/goal.
> >>>> 
> >>>> Main question:
> >>>> -----------------------
> >>>> 
> >>>> Could the rest of you live happily with the above suggestion?
> >>>> 
> >>>> Longer explanation:
> >>>> -------------------------------
> >>>> 
> >>>> For solving a variational problem adaptively such that the error in a
> >>>> certain goal is controlled, the following algorithm is more or less
> >>>> 
> >>>> standard:
> >>>>   # Start with initial mesh T_0
> >>>>   
> >>>>   for i in (0, max_iterations):
> >>>>     (1)  find solution to variational problem on mesh T_i
> >>>>     (2)  estimate global error of solution
> >>>>     
> >>>>     (2b) stop if error is below tolerance
> >>>>     
> >>>>     (3)  estimate local errors (indicators)
> >>>>     (4)  refine mesh T_i -> T_i+1 based on indicators
> >>>> 
> >>>> Steps (2) and (3) depend on the variational problem, and requires the
> >>>> solving of additional variational problems (which can be derived from
> >>>> the original variational problem), and the evaluation of specific
> >>>> additional forms. In other words, these steps require a certain amount
> >>>> of expert knowledge. However, we can automate this by generating the
> >>>> required code.
> >>>> 
> >>>> Here is an example of how this could work:
> >>>> 
> >>>> Ex1.ufl:
> >>>>   V = FiniteElement("CG", "triangle", 1)
> >>>>   f = Coefficient(V)
> >>>>   u = TrialFunction(V)
> >>>>   v = TestFunction(V)
> >>>>   a = dot(grad(u), grad(v))*dx
> >>>>   L = f*v*dx
> >>>>   M = u*dx
> >>>> 
> >>>> Compile with:
> >>>>   ffc -e -l dolfin Ex1.ufl
> >>>> 
> >>>> (where -e == --with-error-control))
> >>>> 
> >>>> FFC will then generate "DOLFIN code" for the forms and finite elements
> >>>> involved in a, L and M, additional forms/finite element spaces, and a
> >>>> (sub-)class
> >>>> 
> >>>>     Ex1::ErrorControl : dolfin::ErrorControl
> >>>> 
> >>>> where dolfin::ErrorControl supplies the methods
> >>>> 
> >>>>     estimate_error(u)
> >>>>     compute_error_indicators(u)
> >>>> 
> >>>> The functionality of this class will take care of steps 2 and 3.
> >>>> 
> >>>> The main file can then look as follows:
> >>>> 
> >>>> main1.cpp:
> >>>>  #include "Ex1.h"
> >>>>  
> >>>>   ...
> >>>>   
> >>>>   Ex1::FunctionSpace V(mesh)
> >>>>   Ex1::BilinearForm a(V, V)
> >>>>   Ex1::LinearForm L(V);
> >>>>   
> >>>>   VariationalProblem(a, L, ...);
> >>>>   
> >>>>   Function u(V);
> >>>>   Ex1::GoalFunctional M(V);
> >>>>   pde.solve(u, tol, M);
> >>>> 
> >>>> Marie thinks this looks pretty clean and knows that it is doable
> >>>> (there is a (more or less) working prototype in the launchpad branches
> >>>> dolfin-error-control/ffc-error-control). (The error control
> >>>> class/object can attached to the GoalFunctional and hence available
> >>>> for the VariationalProblem.)
> >>>> 
> >>>> Remark: If we want to increase explicitness/flexibility we could do
> >>>> 
> >>>> something like instead:
> >>>>   Ex1::ErrorControl ec(...);
> >>>>   pde.solve(u, tol, M, ec);
> >>> 
> >>> What kind of guy would ErrorControl be?
> >> 
> >> The good guy?  ;)
> > 
> > Obviously!
> > 
> >> ErrorControl (or maybe ErrorController or ErrorEstimator) should take
> >> care of estimating the error of an approximate solution u (to be used as
> >> a stopping criterion for mesh refinement) and computing error indicators
> >> (to be used for actual mesh refinement).
> >> 
> >> Take a look at
> >> 
> >>     dolfin/adaptivity/ErrorControl.h
> >> 
> >> from
> >> 
> >>     bzr branch lp:~dolfin-core/dolfin/error-control
> >> 
> >> With the current plan, each generated sub-class will define the
> >> required, problem-specific forms.
> > 
> > Which a user need to do if she want to implement her own scheme?
> 
> I think yes is the answer here (but I'm not sure if I fully understand
> what you ask).
> 
> If a user want to use the brilliant setup devised by us, FFC/DOLFIN will
> take care of everything. If not, the user is on his own, more or less.

:)

> > Is this guy
> > generated for you by UFL and FFC, or is the one in
> > 
> >   dolfin/adaptivity/ErrorControl .h
> > 
> > always used in your prefered example.
> 
> FFC will generate (and the ffc/error-control branch does) a sub-class of
> this guy adapted to a given VariationalProblem F (or a/L) and a given
> goal M. (Same design as for Forms and FunctionSpaces.)

Ok, but in your prefered example you just instantiate the GoalFunctional. 
Where is the subclassed ErrorControl guy, or is the GoalFunctional him in 
disguise?

> >>>> On the other side, if someone wants to supply their own error control,
> >>>> they probably want to control the entire adaptive algorithms as
> >>>> well.
> >>> 
> >>> Can you give a schematical example of how a user would do that for both
> >>> of these suggestions?
> >> 
> >> With the Marie-preferred suggestion, the code is as in main1.cpp.
> >> 
> >> With the other (non-preferred suggestion), I say that if the user wants
> >> to do adaptivity in their own way, I suggest the user does adaptivity in
> >> their own way, namely, implements their own adaptive loop, their own
> >> error estimates, marking strategy etc. I would claim that the error
> >> estimation stage is the more difficult, so if the user knows how to do
> >> that, the remainder should be pretty easy. In essence, these are the
> >> steps in
> >> 
> >>      dolfin/adaptivity/AdaptiveSolver::solve(...)
> >> 
> >> (in the above mentioned branch)
> > 
> > Agree! It just took some time to grasp the concepts, I am probably not
> > there yet ;)
> > 
> >>>> So, I suggest not. Other well-founded opinions are of course
> >>>> welcome ;)
> >>>> 
> >>>> For nonlinear problems, however, a certain amount of differentiations
> >>>> of forms will be required. We therefore need to know what the unknown
> >>>> is, so that we can take the derivative of the forms with respect to
> >>>> it. Anders and Marie have discussed various alternatives quite a
> >>>> bit. In order not to "change everything", we suggest the following,
> >>>> namely, that the user is required to state which coefficient the
> >>>> unknown is:
> >>>> 
> >>>> Ex2.ufl:
> >>>>   V = FiniteElement("CG", "triangle", 1)
> >>>>   f = Coefficient(V)
> >>>>   u = Coefficient(V)
> >>>>   v = TestFunction(V)
> >>>>   F = u*dot(grad(u), grad(v))*dx - f*v*dx
> >>>>   M = u*dx
> >>>>   
> >>>>   # Need to know what the unknown is
> >>>>   unknown = u (!)
> >>> 
> >>> It is probably what we need. Is M the prefered name for an error
> >>> functional? It makes me think of a Matrix with a particular name...
> >> 
> >> Good point. We can easily adjust that.
> >> 
> >>> I also first thought of making unknown a property of F, like:
> >>>   F.unknown = u
> >>> 
> >>> and thought this might be usefull in the python interface too.
> >> 
> >> That is another possibility. If so, u would also need to be attached to
> >> M. I would prefer to keep the labelling as terse as possible.
> > 
> > Sure. For now just one error functional can be defined in a form file.
> > The one with name M, right? Would there be a need to be able to define
> > several functionals in the same file? If so should there be some sort of
> > form attribute for this purpose?
> > 
> >> Anders and I have been through (but rejected) a couple of other ideas,
> >> including
> >> 
> >> (1) Making TrialFunction a special sub-class of Function instead of of
> >> Argument
> >> 
> >>     (Very sensible, consistentifies linear and nonlinear, but way
> >>     radical)
> > 
> > What would a functional of u then be?
> > 
> >   u = TrialFunction(V)
> >   f = assemble(u*dx) ?
> 
> Maybe.

Which means?

> >> (2) Introducing an additional sub-class Unknown of Function, thus
> >> identifying the unknown
> >> 
> >>     (Looks a bit weird, does not increase consistency with linear
> >> 
> >> problems.)
> > 
> > I think this could work. The point with a linear problem is that the
> > problem can be described without the Unknown (=> Linear). Which means
> > that a linear problem wont have an Unknown.
> 
> Yes, but linear problems have a TrialFunction which plays the same role.

But isn't a TrialFunction only meaningfull in a bilinear form? 

> >>> But this is what you use the 'u' argument in pde.solve(u) to indicate?
> >> 
> >> Yes and no. Without ever having looked at the nonlinear solver in
> >> DOLFIN, I assume that the  u argument in pde.solve(u) is (primarily)
> >> needed for the updating in the Newton loop.  The same need is there for
> >> the adaptive solve.
> > 
> > Sure, but you also use it to figure our the Unknown in the pure python
> > code? How could you otherwise do it?
> 
> In the current dolfin Python adaptivity module? We give "u" as an
> argument to AdaptiveVariationalProblem and do all the code generation on
> the fly there.

Ok.

Johan

> --
> Marie