From a3eca510120c508444c7c40a070d95c05ed2217b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 24 Dec 2019 07:13:56 -0700 Subject: [PATCH] More edits for step-71, and also address Timo's comments. --- doc/doxygen/references.bib | 27 +++++++ examples/step-71/doc/intro.dox | 36 ++++++++-- examples/step-71/doc/results.dox | 118 ++++++++++++++++++++----------- 3 files changed, 133 insertions(+), 48 deletions(-) diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 0850244146..c7a5a3c543 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -428,6 +428,33 @@ MRREVIEWER = {Jose Luis Gracia}, journal = {Computer Methods in Applied Mechanics and Engineering} } +@article{Brenner2009, + doi = {10.1093/imanum/drn057}, + url = {https://doi.org/10.1093/imanum/drn057}, + year = {2009}, + month = mar, + publisher = {Oxford University Press ({OUP})}, + volume = {30}, + number = {3}, + pages = {777--798}, + author = {S. C. Brenner and T. Gudi and L.-y. Sung}, + title = {An a posteriori error estimator for a quadratic C0-interior penalty method for the biharmonic problem}, + journal = {{IMA} Journal of Numerical Analysis} +} + +@article{Wells2007, + doi = {10.1016/j.cma.2007.03.008}, + url = {https://doi.org/10.1016/j.cma.2007.03.008}, + year = {2007}, + month = jul, + publisher = {Elsevier {BV}}, + volume = {196}, + number = {35-36}, + pages = {3370--3380}, + author = {Garth N. Wells and Nguyen Tien Dung}, + title = {A C0 discontinuous Galerkin formulation for Kirchhoff plates}, + journal = {Computer Methods in Applied Mechanics and Engineering} +} % ------------------------------------ % References used elsewhere diff --git a/examples/step-71/doc/intro.dox b/examples/step-71/doc/intro.dox index 5f90e3d17a..72943194ab 100644 --- a/examples/step-71/doc/intro.dox +++ b/examples/step-71/doc/intro.dox @@ -217,7 +217,7 @@ differentiable) shape functions with an interior penalty formulation. We base this program on the $C^0$ IP method presented by Susanne Brenner and Li-Yeng Sung in the paper "C$^0$ Interior Penalty Method for Linear Fourth Order Boundary Value Problems on polygonal -domains'' (@cite Brenner2005), where the method is +domains'' @cite Brenner2005 , where the method is derived for the biharmonic equation with "clamped" boundary conditions. @@ -359,10 +359,10 @@ discussed below. Ideally, we would like to measure convergence in the "energy norm" $\|D^2(u-u_h)\|$. However, this does not work because, again, the discrete solution $u_h$ does not have two (weak) derivatives. Instead, -one can define a discrete ($C^0$ IP) norm that is "equivalent" to the +one can define a discrete ($C^0$ IP) semi-norm that is "equivalent" to the energy norm, as follows: @f{align*}{ - \|u_h\|_{h}^2 := + |u_h|_{h}^2 := \sum\limits_{K \in \mathbb{T}} \big|u_h\big|_{H^2(K)}^2 + \sum\limits_{e \in \mathbb{F} } @@ -370,10 +370,10 @@ energy norm, as follows: \ \!\bigg[ \!\bigg[ \frac{\partial u_h}{\partial \mathbf n}\bigg]\!\bigg] \right\|_{L^2(e)}^2. @f} -In this norm, the theory in the paper mentioned above yields that we +In this (semi)norm, the theory in the paper mentioned above yields that we can expect @f{align*}{ - \|u-u_h\|_{h}^2 = {\cal O}(h^{p-1}), + |u-u_h|_{h}^2 = {\cal O}(h^{p-1}), @f} much as one would expect given the convergence rates we know are true for the usual discretizations of the Laplace equation. @@ -394,12 +394,34 @@ large. In other words, the only reasonably choices for $p$ are $p\le m+3$ because larger polynomial degrees do not result in higher convergence orders. +For the purposes of this program, we're a bit too lazy to actually +implement this equivalent norm -- though it's not very difficult and +would make for a good exercise. Instead, we'll simply check in the +program what the "broken" $H^2$ seminorm +@f{align*}{ + |u_h|^\circ_{H^2}^2 + := + \sum\limits_{K \in \mathbb{T}} \big|u_h\big|_{H^2(K)}^2 + = + \sum\limits_{K \in \mathbb{T}} \big|D^2 u_h\big|_{L_2}^2 +@f} +yields. The convergence rate in this norm can, from a theoretical +perspective, of course not be worse than the one for +$|\cdot|_h$, but it could be worse. It could also be the case that +we get the optimal convergence rate even though there is a bug in the +program, and that that bug would only show up in sub-optimal rates for +the additional terms present in $|\cdot|_h$. But, one might hope +that if we get the optimal rate in the broken norm and the norms +discussed below, then the program is indeed correct. The results +section will demonstrate that we obtain optimal rates in all norms +shown. + Convergence in the $L_2$-norm The optimal convergence rate in the $L_2$-norm is $\mathcal{O}(h^{p+1})$ -provided $p \ge 3$. More details can be found in Theorem 4.6 of @cite -Engel2002 . +provided $p \ge 3$. More details can be found in Theorem 4.6 of +@cite Engel2002 . The default in the program below is to choose $p=2$. In that case, the theorem does not apply, and indeed one only gets $\mathcal{O}(h^2)$ diff --git a/examples/step-71/doc/results.dox b/examples/step-71/doc/results.dox index 2501d2cdcb..3c6fa3d00b 100644 --- a/examples/step-71/doc/results.dox +++ b/examples/step-71/doc/results.dox @@ -1,34 +1,42 @@

Results

-We run the program with a right hand side that will produce the -solution $u = \sin(\pi x) \sin(\pi y)$ and with clamped -boundary conditions in the domain $\Omega = (0,1)^2$. +We run the program with right hand side and boundary values as +discussed in the introduction. These will produce the +solution $u = \sin(\pi x) \sin(\pi y)$ on the domain $\Omega = (0,1)^2$. We test this setup using $Q_2$, $Q_3$, and $Q_4$ elements, which one can change via the `fe_degree` variable in the `main()` function. With mesh -refinement, the $L_2$ convergence rates, $H_1$-seminorm convergence -and $H_2$-seminorm convergence of $u$ -should then be around 2, 2, 1 for $Q_2$; 4, 3, 2 for +refinement, the $L_2$ convergence rates, $H^1$-seminorm rate, +and $H^2$-seminorm convergence of $u$ +should then be around 2, 2, 1 for $Q_2$ (with the $L_2$ norm +sub-optimal as discussed in the introduction); 4, 3, 2 for $Q_3$; and 5, 4, 3 for $Q_4$, respectively. -From the papers by Brenner et al., it is not immediately clear what -the penalty parameter $\gamma$ should be. Educated guesses, comparing +From the literature, it is not immediately clear what +the penalty parameter $\gamma$ should be. For example, +@cite Brenner2009 state that it needs to be larger than one, and +choose $\gamma=5$. The FEniCS/Dolphin tutorial chooses it as +$\gamma=8$. @cite Wells2007 uses a value for $\gamma$ larger than the +number of edges belonging to an element for Kirchhoff plates (see +their Section 4.2). This suggests that maybe +$\gamma = 1$, $2$, are too small; on the other hand, a value +$p(p+1)$ would be reasonable, +where $p$ is the degree of polynomials. The last of these choices is +the one one would expect to work by comparing to the discontinuous Galerkin formulations for the Laplace equation, -suggest that $\gamma = 1$, $2$, and $p(p+1)$ would all be reasonable, -where $p$ is the degree of polynomials. This is easy to change -in the code from its current default. Below we show results for -all of these. +and it will turn out to also work here. +But we should check what value of $\gamma$ is right, and we will do so +below; changing $\gamma$ is easy in the two `face_worker` and +`boundary_worker` functions defined in `assemble_system()`.

Test results on Q2 with γ = p(p+1)

-

Convergence table

- We run the code with differently refined meshes and get the following convergence rates. - + @@ -44,17 +52,19 @@ and get the following convergence rates.
Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H_1}$ Conv. rates $|u-u_h|_{H_2}$ Conv. rates Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H^1}$ Conv. rates $|u-u_h|_{H^2}$ Conv. rates
2 1.539e-02 8.528e-02 1.602
We can see that the $L_2$ convergence rates are around 2, -$H_1$-seminorm convergence rates are around 2, -and $H_2$-seminorm convergence rates are around 1. +$H^1$-seminorm convergence rates are around 2, +and $H^2$-seminorm convergence rates are around 1. The latter two +match the theoretically expected rates; for the former, we have no +theorem but are not surprised that it is sub-optimal given the remark +in the introduction. +

Test results on Q3 with γ = p(p+1)

-

Convergence table

- - + @@ -70,18 +80,16 @@ and $H_2$-seminorm convergence rates are around 1.
Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H_1}$ Conv. rates $|u-u_h|_{H_2}$ Conv. rates Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H^1}$ Conv. rates $|u-u_h|_{H^2}$ Conv. rates
2 2.187e-04 4.46269e-03 1.638e-01
We can see that the $L_2$ convergence rates are around 4, -$H_1$-seminorm convergence rates are around 3, -and $H_2$-seminorm convergence rates are around 2. +$H^1$-seminorm convergence rates are around 3, +and $H^2$-seminorm convergence rates are around 2. This, of course, matches our theoretical expectations. -

Test results on Q4 with γ = p(p+1)

- -

Convergence table

+

Test results on Q4 with γ = p(p+1)

- + @@ -97,21 +105,23 @@ This, of course, matches our theoretical expectations.
Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H_1}$ Conv. rates $|u-u_h|_{H_2}$ Conv. rates Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H^1}$ Conv. rates $|u-u_h|_{H^2}$ Conv. rates
2 8.34446e-06 0.000239323 0.0109785
We can see that the $L_2$ norm convergence rates are around 5, -$H_1$-seminorm convergence rates are around 4, -and $H_2$-seminorm convergence rates are around 3. +$H^1$-seminorm convergence rates are around 4, +and $H^2$-seminorm convergence rates are around 3. On the finest mesh, the $L_2$ norm convergence rate is much smaller than our theoretical expectations because the linear solver becomes the limiting factor due -to round-off. But the $L_2$ error is pretty small in that case. +to round-off. Of course the $L_2$ error is also very small already in +that case. -

Test results on Q2 with γ = 1

-

Convergence table

+

Test results on Q2 with γ = 1

+For comparison with the results above, let us now also consider the +case where we simply choose $\gamma=1$: - + @@ -126,20 +136,21 @@ to round-off. But the $L_2$ error is pretty small in that case.
Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H_1}$ Conv. rates $|u-u_h|_{H_2}$ Conv. rates Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H^1}$ Conv. rates $|u-u_h|_{H^2}$ Conv. rates
2 4.86048e-02 3.30386e-01 4.34917 5 8.41058e-04 1.98 3.57705e-02 1.00 3.95719 0.001
-Although $L_2$ norm and $H_1$-seminorm convergence rates of $u$ -follow the theoretical expectations, $H_2$-seminorm does not converge. +Although $L_2$ norm and $H^1$-seminorm convergence rates of $u$ +follow the theoretical expectations, the $H^2$-seminorm does not seem to converge. Comparing results from $\gamma = 1$ and $\gamma = p(p+1)$, it is clear that $\gamma = p(p+1)$ is a better penalty.

Test results on Q2 with γ = 2

-

Convergence table

- +Since $\gamma=1$ is clearly too small, one might conjecture that +$\gamma=2$ might actually work better. Here is what one obtains in +that case: - + @@ -154,10 +165,12 @@ $\gamma = p(p+1)$ is a better penalty.
Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H_1}$ Conv. rates $|u-u_h|_{H_2}$ Conv. rates Number of refinements $\|u-u_h^\circ\|_{L_2}$ Conv. rates $|u-u_h|_{H^1}$ Conv. rates $|u-u_h|_{H^2}$ Conv. rates
2 5.482e-03 7.652e-02 1.756e-00 5 2.822e-04 1.68 1.605e-03 1.90 2.041e-01 1.03
-In this case, all convergence rates of $u$ -follow the theoretical expectations, but, -compared to the results from $\gamma = p(p+1)$, -we don't see good convergence of the $L_2$ errors. +In this case, the convergence rates more or less follow the +theoretical expectations, but, compared to the results from $\gamma = +p(p+1)$, are more variable. That suggests that the penalty parameter +is already too small even for the value $p=2$ chosen here, and one can +readily check that it is indeed too small when using higher polynomial +degrees.

Conclusions for the choice of the penalty parameter

@@ -165,3 +178,26 @@ we don't see good convergence of the $L_2$ errors. The conclusions for which of the "reasonable" choices one should use for the penalty parameter is that $\gamma=p(p+1)$ yields the expected results. It is, consequently, what the code uses as currently written. + + +

Possibilities for extensions

+ +There are a number of obvious extensions to this program that would +make sense: + +- The program uses a square domain and a uniform mesh. Real problems + don't come this way, and one should verify convergence also on + domains with other shapes and, in particular, curved boundaries. One + may also be interested in resolving areas of less regularity by + using adaptive mesh refinement. + +- From a more theoretical perspective, the convergence results above + only used the "broken" $H^2$ seminorm $|\cdot|^\circ_{H^2}$ instead + of the "equivalent" norm $|\cdot|_h$. This is good enough to + convince ourselves that the program isn't fundamentally + broken. However, it might be interesting to measure the error in the + actual norm for which we have theoretical results. Implementing this + addition should not be overly difficult using, for example, the + FEInterfaceValues class combined with MeshWorker::mesh_loop() in the + same spirit as we used for the assembly of the linear system. + -- 2.39.5