dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

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

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

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

Johan