non-trivial finite element program that worked right from the start. It is
therefore necessary to find ways to verify whether a computed solution is
correct or not. Usually, this is done by choosing the set-up of a simulation
-in such a way that we know the exact continuous solution and evaluate the difference
-between continuous and computed discrete solution. If this difference
+in such a way that *we know the exact continuous solution* and *evaluate the difference
+between continuous and computed discrete solution*. If this difference
converges to zero with the right order of convergence, this is already a good
indication of correctness, although there may be other sources of error
persisting which have only a small contribution to the total error or are of
higher order. In the context of finite element simulations, this technique
of picking the solution by choosing appropriate right hand sides and
boundary conditions
-is often called the <i>Method of Manufactured Solution</i>.
+is often called the <i>Method of Manufactured Solution</i>. (We will come
+back to how exactly we construct the solution in this method below, after
+discussing the equation we want to solve.)
In this example, we will not go into the theories of systematic software
-verification which is a very complicated problem. Rather we will demonstrate
+verification which is a complicated problem in general. Rather we will demonstrate
the tools which deal.II can offer in this respect. This is basically centered
around the functionality of a single function, VectorTools::integrate_difference().
This function computes the difference between a given continuous function and
operator; the <a
href="https://en.wikipedia.org/wiki/Helmholtz_equation">equation with
the "bad sign"</a> is $-\Delta u - \alpha u$ and results from modeling
-time-harmonic processes. The operator is not positive
+time-harmonic processes. For the equation with the "bad sign", the
+operator $-\Delta-\alpha I$ is not positive
definite if $\alpha>0$ is large, and this leads to all sorts of issues
we need not discuss here. The operator may also not be invertible --
i.e., the equation does not have a unique solution -- if $\alpha$
happens to be one of the eigenvalues of $-\Delta$.)
-Because we want to verify the convergence of our numerical solution $u_h$,
-we want a setup so that we know the exact solution $u$. This is where
-the Method of Manufactured Solutions comes in. To this end, let us
-choose a function
-@f[
- \bar u(x) = \sum_{i=1}^3 \exp\left(-\frac{|x-x_i|^2}{\sigma^2}\right)
-@f]
-where the centers $x_i$ of the exponentials are
- $x_1=(-\frac 12,\frac 12)$,
- $x_2=(-\frac 12,-\frac 12)$, and
- $x_3=(\frac 12,-\frac 12)$,
-and the half width is set to $\sigma=\frac {1}{8}$. The method of manufactured
-solution then says: choose
-@f{align*}
- f &= -\Delta \bar u + \bar u, \\
- g_1 &= \bar u|_{\Gamma_1}, \\
- g_2 &= {\mathbf n}\cdot \nabla\bar u|_{\Gamma_2}.
-@f}
-With this particular choice, we infer that of course the solution of the
-original problem happens to be $u=\bar u$. In other words, by choosing
-the right hand sides of the equation and the boundary conditions in a
-particular way, we have manufactured ourselves a problem to which we
-know the solution. This allows us then to compute the error of our
-numerical solution. In the code below, we represent $\bar u$ by the
-<code>Solution</code> class, and other classes will be used to
-denote $\bar u|_{\Gamma_1}=g_1$ and ${\mathbf n}\cdot \nabla\bar u|_{\Gamma_2}=g_2$.
-
Using the above definitions, we can state the weak formulation of the
equation, which reads: find $u\in H^1_g=\{v\in H^1: v|_{\Gamma_1}=g_1\}$ such
that
following program.
+<h3>The method of manufactured solutions</h3>
+
+
+Because we want to verify the convergence of our numerical solution $u_h$,
+we want a setup so that we know the exact solution $u$. This is where
+the Method of Manufactured Solutions comes in: Let us
+choose a function
+@f[
+ \bar u(\mathbf x) =
+ \sum_{i=1}^3 \exp\left(-\frac{|\mathbf x-\mathbf x_i|^2}{\sigma^2}\right)
+@f]
+where the centers $x_i$ of the exponentials are
+ $\mathbf x_1=(-\frac 12,\frac 12)$,
+ $\mathbf x_2=(-\frac 12,-\frac 12)$, and
+ $\mathbf x_3=(\frac 12,-\frac 12)$,
+and the half width is set to $\sigma=\frac {1}{8}$. The method of manufactured
+solution then says: choose
+@f{align*}
+ f &= -\Delta \bar u + \bar u, \\
+ g_1 &= \bar u|_{\Gamma_1}, \\
+ g_2 &= {\mathbf n}\cdot \nabla\bar u|_{\Gamma_2}.
+@f}
+With this particular choice for $f,g_1,g_2$, the solution of the
+original problem must necessarily be $u=\bar u$. In other words, by choosing
+the right hand sides of the equation and the boundary conditions in a
+particular way, we have manufactured ourselves a problem to which we
+know the solution -- a very useful case given that in all but the very
+simplest cases, PDEs do not have solutions we can just write down.
+This then allows us to compute the error of our
+numerical solution. In the code below, we represent $\bar u$ by the
+<code>Solution</code> class, and other classes will be used to
+denote $\bar u|_{\Gamma_1}=g_1$ and ${\mathbf n}\cdot \nabla\bar u|_{\Gamma_2}=g_2$.
+
+@note In principle, you can choose whatever you want for the function $\bar u$
+ above -- here we have simply chosen a sum of three exponentials. In practice,
+ there are two considerations you want to take into account: (i) The function
+ must be simple enough so that you can compute derivatives of the function
+ with not too much effort, for example in order to determine what
+ $f = -\Delta \bar u + \bar u$ is. Since the derivative of an exponential
+ is relatively straightforward to compute, the choice above satisfies this
+ requirement, whereas a function of the kind
+ $\bar u(\mathbf x) = \atan\left(\|\mathbf x\|^{\|\mathbf x\|}\right)$
+ would have presented greater difficulties.
+ (ii) You *don't* want $\bar u$ be a polynomial of low degree. That is
+ because if you choose the polynomial degree of your finite element
+ sufficiently high, you can *exactly* represent this $\bar u$ with
+ the numerical solution $u_h$, making the error zero regardless of
+ how coarse or fine the mesh is. Verifying that this is so is a useful
+ step, but it will not allow you to verify the correct order of
+ convergence of $\|u-u_h\|$ as a function of the mesh size $h$ in
+ the general case of arbitrary $f$.
+ (iii) The typical finite element error estimates assume sufficiently
+ smooth solutions, i.e., sufficiently smooth domains, right-hand sides
+ $f$ and boundary conditions. As a consequence, you should choose
+ a smooth solution $\bar u$ -- for example, it shouldn't
+ have kinks. (iv) You want a solution whose variations can be resolved
+ on the meshes you consider to test convergence. For example, if you
+ were to choose $\bar u(\mathbf x)=\sin(1000 x_1)\sin(1000 x_2)$,
+ you shouldn't be surprised
+ if you don't observe that the error decreases at the expected rate
+ until your mesh is fine enough to actually resolve the high-frequency
+ oscillations with substantially more than 1,000 mesh cells in each
+ coordinate direction.
+
+ The solution $\bar u$ we choose here satisfies all of these
+ requirements: (i) It is relatively straightforward to
+ differentiate; (ii) it is not a polynomial; (iii) it is smooth;
+ and (iv) it has a length scale of $\sigma=\frac {1}{8}$ which,
+ on the domain $[-1,1]^d$ is relatively straightforward to
+ resolve with 16 or more cells in each coordinate direction.
+
+
+
<h3>A note on good programming practice</h3>
Besides the mathematical topics outlined above, we also want to use this
projects / dealii.git / commitdiff