dolfin team mailing list archive

[Marie Rognes <meg@xxxxxxxxx>]

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

On the other side, if someone wants to supply their own error control,
they probably want to control the entire adaptive algorithms as
well. 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 (!)

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.


--
Marie