From: bangerth Date: Tue, 19 Mar 2013 12:08:31 +0000 (+0000) Subject: Write one section. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=592957e5574e56ffbd5f66e177b521609f3f843f;p=dealii-svn.git Write one section. git-svn-id: https://svn.dealii.org/trunk@28941 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-26/doc/intro.dox b/deal.II/examples/step-26/doc/intro.dox index c59cfdf609..51dcc23237 100644 --- a/deal.II/examples/step-26/doc/intro.dox +++ b/deal.II/examples/step-26/doc/intro.dox @@ -1,2 +1,115 @@

Introduction

+ + +

Verifying whether the code is correct

+ +There are a number of things one can typically get wrong when implementing a +finite element code. In particular, for time dependent problems, the following +are common sources of bugs: +- The time integration, for example by getting the coefficients in front of + the terms involving the current and previous time steps wrong (e.g., mixing + up a factor $\theta$ for $1-\theta$). +- Handling the right hand side, for example forgetting a factor of $k_n$ or + $\theta$. +- Mishandling the boundary values, again for example forgetting a factor of + $k_n$ or $\theta$, or forgetting to apply nonzero boundary values not only + to the right hand side but also to the system matrix. + +A less common problem is getting the initial conditions wrong because one can +typically see that it is wrong by just outputting the first time step. In any +case, in order to verify the correctness of the code, it is helpful to have a +testing protocol that allows us to verify each of these components +separately. This means: +- Testing the code with nonzero initial conditions but zero right hand side + and boundary values and verifying that the time evolution is correct. +- Then testing with zero initial conditions and boundary values but nonzero + right hand side and again ensuring correctness. +- Finally, testing with zero initial conditions and right hand side but + nonzero boundary values. + +This sounds complicated, but fortunately, for linear partial differential +equations without coefficients (or constant coefficients) like the one here, +there is a fairly standard protocol that rests on the following observation: +if you choose as your domain a square $[0,1]^2$ (or, with slight +modifications, a rectangle), then the exact solution can be written as +@f{align*} + u(x,y,t) = a(t) \sin(n_x \pi x) \sin(n_y \pi y) +@f} +(with integer constants $n_x,n_y$) +if only the initial condition, right hand side and boundary values are all +of the form $\sin(n_x \pi x) \sin(n_y \pi y)$ as well. This is due to the fact +that the function $\sin(n_x \pi x) \sin(n_y \pi y)$ is an eigenfunction of the +Laplace operator and allows us to compute things like the time factor $a(t)$ +analytically and, consequently, compare with what we get numerically. + +As an example, let us consider the situation where we have +$u_0(x,y)=\sin(n_x \pi x) \sin(n_x \pi y)$ and +$f(x,y,t)=0$. With the claim (ansatz) of the form for +$u(x,y,t)$ above, we get that +@f{align*} + \left(\frac{\partial}{\partial t} -\Delta\right) + u(x,y,t) + &= + \left(\frac{\partial}{\partial t} -\Delta\right) + a(t) \sin(n_x \pi x) \sin(n_y \pi y) + \\ + &= + \left(a'(t) + (n_x^2+n_y^2)\pi^2 a(t) \right) \sin(n_x \pi x) \sin(n_y \pi y). +@f} +For this to be equal to $f(x,y,t)=0$, we need that +@f{align*} + a'(t) + (n_x^2+n_y^2)\pi^2 a(t) = 0 +@f} +and due to the initial conditions, $a(0)=1$. This differential equation can be +integrated to yield +@f{align*} + a(t) = - e^{-(n_x^2+n_y^2)\pi^2 t}. +@f} +In other words, if the initial condition is a product of sines, then the +solution has exactly the same shape of a product of sines that decays to zero +with a known time dependence. This is something that is easy to test if you +have a sufficiently fine mesh and sufficiently small time step. + +What is typically going to happen if you get the time integration scheme wrong +(e.g., by having the wrong factors of $\theta$ or $k$ in front of the various +terms) is that you don't get the right temporal behavior of the +solution. Double check the various factors until you get the right +behavior. You may also want to verify that the temporal decay rate (as +determined, for example, by plotting the value of the solution at a fixed +point) does not double or halve each time you double or halve the time step or +mesh size. You know that it's not the handling of the +boundary conditions or right hand side because these were both zero. + +If you have so verified that the time integrator is correct, take the +situation where the right hand side is nonzero but the initial conditions are +zero: $u_0(x,y)=0$ and +$f(x,y,t)=\sin(n_x \pi x) \sin(n_x \pi y)$. Again, +@f{align*} + \left(\frac{\partial}{\partial t} -\Delta\right) + u(x,y,t) + &= + \left(\frac{\partial}{\partial t} -\Delta\right) + a(t) \sin(n_x \pi x) \sin(n_y \pi y) + \\ + &= + \left(a'(t) + (n_x^2+n_y^2)\pi^2 a(t) \right) \sin(n_x \pi x) \sin(n_y \pi y), +@f} +and for this to be equal to $f(x,y,t)$, we need that +@f{align*} + a'(t) + (n_x^2+n_y^2)\pi^2 a(t) = 1 +@f} +and due to the initial conditions, $a(0)=0$. Integrating this equation in time +yields +@f{align*} + a(t) = \frac{1}{(n_x^2+n_y^2)\pi^2} \left[ 1 - e^{-(n_x^2+n_y^2)\pi^2 t} \right]. +@f} + +Again, if you have the wrong factors of $\theta$ or $k$ in front of the right +hand side terms you will either not get the right temporal behavior of the +solution, or it will converge to a maximum value other than +$\frac{1}{(n_x^2+n_y^2)\pi^2}$. + +Once we have verified that the time integration and right hand side handling +are correct using this scheme, we can go on to verifying that we have the +boundary values correct, using a very similar approach.