dolfin team mailing list archive

[Marie Rognes <meg@xxxxxxxxx>]

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

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.


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

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

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)

(2) Introducing an additional sub-class Unknown of Function, thus
identifying the unknown
    (Looks a bit weird, does not increase consistency with linear problems.)


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


--
Marie

>> As before, compile with:
>>
>>   ffc -e -l dolfin Ex2.ufl
>>
>> And the main program can read:
>>
>> main2.cpp:
>>
>>   #include "Ex2.h"
>>
>>   ...
>>
>>   Ex1::FunctionSpace V(mesh)
>>   Ex1::LinearForm F(V, V)
>>
>>   Function u(V);
>>   F.u = u;
>>
>>   VariationalProblem(F, ...);
>>   Ex1::GoalFunctional M;
>>   pde.solve(u, tol, M);
>>
>>
>> Marie has not given much thought to how this should be wrapped "back"
>> to python, but is sure that the python-wrapper-guru(s) will find out.
>>     
> :)
>
> Johan
>
>   
>> --
>> Marie
>>
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~dolfin
>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help   : https://help.launchpad.net/ListHelp
>>