dolfin team mailing list archive

[Marie Rognes <meg@xxxxxxxxx>]

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


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

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

--
Marie