some problems concerning step-14

[曾元圆]

Hi, nowadays I'm learning about the goal oriented error estimator and I read the tutorial of step-14 (http://www.dealii.org/developer/doxygen/deal.II/step_14.html.) But I'm confused about some problems and hope you can help me with these:

1. When deriving the error with respect to the functional, why must we change J(e)=a(e,z) to J(e)=a(e,z-z_h) ?

I know the dual solution z must be approximated in a richer space than the primal solution, otherwise J(e) will be 0. But why not just solve the dual problem in a richer space without subtracting its interpolation to the primal space? I didn't see the necessity to introduce z_h into the formula.

2. Why must we change J(e)=∑(f+△ uh,z-zh)-(∂uh,z-zh) to J(e)=∑(f+△ uh,z-zh)-(1/2[∂uh],z-zh)? Is it just an implementation consideration for saving computational effort? Can this kind of rewriting be generally adopted in other kind of problems(e.g in advection problem where the face integrals in J(e) relies on the upstream information)?

Thank you very much!

[曾元圆]

3. In the member function: WeightedResidual<dim>::refine_grid(), you firstly defined a Vector error_indicators whose dimension is n_active_cells of the triangulation, then you call the function "estimate_error(error_indicators)" to compute the error_indicators. But in the function "estimate_error(Vector<float>& error_indicators)", I noticed that you reinitialize the error_indicators by calling "error_indicators.reinit(DualSolver<dim>::dof_handler.get_triangulation().n_active_cells());", so why you do this?   Since the PrimalSolver and DualSolver are both derived from the base class and they share the same triangulation, then what is the purpose of using DualSolver<dim>::dof_handler.get_triangulation to reinit error_indicator?

[Wolfgang Bangerth]

On 03/02/2018 03:36 AM, 曾元圆 wrote:
> Hi, nowadays I'm learning about the goal oriented error estimator and I
> read the tutorial of step-14
> (http://www.dealii.org/developer/doxygen/deal.II/step_14.html.) But I'm
> confused about some problems and hope you can help me with these: 1. When
> deriving the error with respect to the functional, why must we change
> J(e)=a(e,z) to J(e)=a(e,z-z_h) ? I know the dual solution z must be
> approximated in a richer space than the primal solution, otherwise J(e)
> will be 0. But why not just solve the dual problem in a richer space
> without subtracting its interpolation to the primal space? I didn't see the
> necessity to introduce z_h into the formula.

You are correct: the values you get from both of these formulas are exactly
the same. So it is not *necessary* to introduce z_h if you are interested in
computing the *error*.

We introduce the interpolant because z-z_h is a quantity that is only large
where the dual solution is rough, where z may be large also in other places.
Doing so ensures that the error estimator is *localized*: It is large exactly
where the primal and dual solutions are rough, i.e., where we expect the error
to be caused. As mentioned above, if you sum the contributions of all cells,
you will get the same value whether you introduce z_h or not, but the
contribution of each cell is going to be different. If you want the
contributions of each cell to serve as a good mesh refinement criterion, then
it turns out that you need to introduce z_h.

> 2. Why must we change J(e)=∑(f+△ uh,z-zh)-(∂uh,z-zh) to J(e)=∑(f+△
> uh,z-zh)-(1/2[∂uh],z-zh)? Is it just an implementation consideration for
> saving computational effort?

For the same kind of reason. For a smooth (primal) solution, the term
(∂uh,z-zh)
may be large because the normal derivative of the primal solution may simply
be what it is -- think of, for example, a linear exact solution u that can be
perfectly approximated by u_h. So if you leave this term as is, this would
suggest that the error on this cell is large. But that's wrong -- the error is
actually quite small because you can approximate linear functions well.

On the other hand, if you do the rewrite (which again leaves the *sum* over
all cells the same, but change the contributions of each cell), then
(1/2[∂uh],z-zh)
is going to be small because while the normal derivative of u_h may be large,
the *jump* of the normal derivative is small if the solution is linear or
nearly so.

Another way of seeing this is to think of both of the terms in J(e) as
residual times dual weight.
The residuals here are f+△ uh and 1/2[∂uh]. You want to define these residuals
in such a way that they are zero for the exact solution. That is true for
these two residuals: f+△ u is zero because 'u' satisfies the equation, and
1/2[∂u] is zero because the solutions of the Laplace equation have continuous
gradients (if f is smooth enough).

On the other hand, the term ∂u is not zero, even for the exact solution.

> Can this kind of rewriting be generally adopted in other kind of
> problems(e.g in advection problem where the face integrals in J(e) relies
> on the upstream information)?

Yes.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Wolfgang Bangerth]

Yes, this seems duplicative and unnecessary. The reinit call in
estimate_error() in line 2356 could be removed.

Do you want to try to submit a patch for this and become one of the authors of
deal.II? Or would you prefer me doing this? (We always appreciate newcomers'
patches, and this case would be no different -- and I think this would
actually make for a really nice first contribution!)

Best
Wolfgang

[曾元圆]

Thank you so much Bangerth!

Now I understand why we need to rewrite the error formula on a cell as residual times dual weight. But I'm still a little confused with the reason why we must introduce z_h. 

Just as you mentioned, if we introduce z_h, then z-z_h is a quantity that is only large where the dual solution is rough. But why do we need to care about the accuracy of z here? I think the only  thing we need to care about is the value of z on that cell, because z is a quantity that represents how important the residual on that cell is.

  

My understanding is: now the dual_weight z-z_h does not only represent how important the residual on a certain cell is, but also tells us some information about how good the dual solution on that cell is. But another problem is, does z-z_h still has the same tendency as z? If not, how z-z_h can represent the importance of a certan cell as z can?

 

I'm not sure if my understanding is correct. I tried to run the code using only z as dual_weights, and I found the result almost the same as that using z-z_h.

Finally, I am certainly glad to submit patches to deal.II and make my own contribution. But I didn't fork deal.ii on my github account yet, and this is relatively a small issue, so I will be glad if you can do it for the moment. 

在 2018年3月4日星期日 UTC+8上午1:42:27，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 03/06/2018 08:40 AM, 曾元圆 wrote:
>
> Now I understand why we need to rewrite the error formula on a cell as
> residual times dual weight. But I'm still a little confused with the reason
> why we must introduce z_h.
> Just as you mentioned, if we introduce z_h, then z-z_h is a quantity that is
> only large where the dual solution is rough. But why do we need to care about
> the accuracy of z here? I think the only thing we need to care about is the
> value of z on that cell, because z is a quantity that represents how important
> the residual on that cell is.

No. z tells you how important the *locally generated error is for the global
error functional*. (That is because z is the Green's function associated with
your error functional.) But you don't have the local error. All you have is
the local residual.

> My understanding is: now the dual_weight z-z_h does not only represent how
> important the residual on a certain cell is, but also tells us some
> information about how good the dual solution on that cell is. But another
> problem is, does z-z_h still has the same tendency as z?

Almost. Think of it as z-phi_h where you can choose phi_h as you want. For
example, on each cell you can think of choosing phi_h so that it cancels the
constant and linear term of the Taylor expansion of z. Then z-z_h would
contain the quadratic and higher order Taylor terms, i.e. something like
z''*(x-x0)^2 where x0 can be chosen as a point on the cell.

> If not, how z-z_h can
> represent the importance of a certan cell as z can?
> I'm not sure if my understanding is correct. I tried to run the code using
> only z as dual_weights, and I found the result almost the same as that using
> z-z_h.

Nice idea to try this out. Do you get "almost the same" overall error
estimate, or "almost the same" mesh?

I think all of these are good questions to ask. Although I have worked on this
for a long time, I can not actually give you a particularly good answer for
all of this. I am sure others who are more versed in the theory of errors,
residuals, etc could tell you the precise reason for why it is in fact
necessary to subtract z_h. The best I can say is that that's the way I've
always seen it done, and while I have a vague idea why that is so (see above),
I can't say that I can describe it well enough to explain it.

> Finally, I am certainly glad to submit patches to deal.II and make my own
> contribution. But I didn't fork deal.ii on my github account yet, and this is
> relatively a small issue, so I will be glad if you can do it for the moment.

OK, I will take care of this then.

[Thomas Wick]

In addition what Wolfgang said:

1. It would be indeed interesting to
see whether neglecting z_h really yields
the same error, the same effectivity indices,
and/or the same mesh.

2. Going back to your initial questions:
Inserting z_h is the key when classical
a posteriori bounds in terms of the mesh size h
are of interest.

The final goal is usually to obtain an error estimate
in terms of the mesh size h in order to
quantify the order of convergence of your scheme.

For this reason you need to insert z_h such that

|| z - z_h||

to apply interpolation estimates that give you
some h^{a} on the right hand side:

|| z-z_h || = O(h^a)

with the order a.


3. In practice you have indeed different choices
how to evaluate J(u) - J(u_h).

Also some people do not integrate back
into the strong form and work
with a weak form of the error estimator, which
has the advantage that no second-order operators
and partial integration needs to be applied.

From this point of view, it may be that neglecting
z_h could work in practice. But as said above,
a careful study for some model problems
would be useful.


Best Thomas W.

--
++--------------------------------------------++
Prof. Dr. Thomas Wick
Institut für Angewandte Mathematik (IfAM)
Leibniz Universität Hannover
Welfengarten 1
30167 Hannover, Germany

Tel.: +49 511 762 3360
Email: thoma...@ifam.uni-hannover.de
www: http://www.ifam.uni-hannover.de/wick
www: http://www.cmap.polytechnique.fr/~wick/
++--------------------------------------------++
--

[Wolfgang Bangerth]

On 03/06/2018 08:40 AM, 曾元圆 wrote:
>

> Finally, I am certainly glad to submit patches to deal.II and make my own
> contribution. But I didn't fork deal.ii on my github account yet, and this is
> relatively a small issue, so I will be glad if you can do it for the moment.

See here:
https://github.com/dealii/dealii/pull/6029
In any case, we do appreciate contributions by everyone, so if you find
something like this in the future, please do feel encouraged to submit a patch
yourself!