From 77712cfa8c068967ab65b5d003254d943f7e3017 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 16 Dec 2013 20:17:46 +0000 Subject: [PATCH] Make some progress with the documentation. git-svn-id: https://svn.dealii.org/trunk@32031 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-26/doc/intro.dox | 321 ++++++++++++++++++++++- deal.II/examples/step-26/doc/results.dox | 99 +++++++ deal.II/examples/step-26/step-26.cc | 65 +++-- 3 files changed, 465 insertions(+), 20 deletions(-) diff --git a/deal.II/examples/step-26/doc/intro.dox b/deal.II/examples/step-26/doc/intro.dox index 51dcc23237..3641c4b904 100644 --- a/deal.II/examples/step-26/doc/intro.dox +++ b/deal.II/examples/step-26/doc/intro.dox @@ -1,8 +1,259 @@

Introduction

+This program implements the heat equation +@f{align*} + \frac{\partial u(\mathbf x, t)}{\partial t} + - + \Delta u(\mathbf x, t) + &= + f(\mathbf x, t), + \qquad\qquad && + \forall \mathbf x \in \Omega, t\in (0,T), + \\ + u(\mathbf x, 0) &= u_0(\mathbf x) && + \forall \mathbf x \in \Omega, \\ + \\ + u(\mathbf x, t) &= g(\mathbf x,t) && + \forall \mathbf x \in \partial\Omega, t \in (0,T). +@f} +In some sense, this equation is simpler than the ones we have discussed in the +preceding programs step-23, step-24, step-25, namely the wave equation. This +is due to the fact that the heat equation smoothes out the solution over time, +and is consequently more forgiving in many regards. For example, when using +implicit time stepping methods, we can actually take large time steps, we have +less trouble with the small disturbances we introduce through adapting the +mesh every few time steps, etc. + +Our goal here will be to solve the equations above using the theta-scheme that +discretizes the equation in time using the following approach, where we would +like $u^n(\mathbf x)$ to approximate $u(\mathbf x, t_n)$ at some time $t_n$: +@f{align*} + \frac{u^n(\mathbf x)-u^{n-1}(\mathbf x)}{k_n} + - + \left[ + (1-\theta)\Delta u^{n-1}(\mathbf x) + + + \theta\Delta u^n(\mathbf x) + \right] + &= + \left[ + (1-\theta)f(\mathbf x, t_{n-1}) + + + \theta f(\mathbf x, t_n) + \right]. +@f} +Here, $k_n=t_n-t_{n-1}$ is the time step size. The theta-scheme generalizes +the explicit Euler ($\theta=0$), implicit Euler ($\theta=1$) and +Crank-Nicolson ($\theta=\frac 12$) time discretizations. Since the latter has +the highest convergence order, we will choose $\theta=\frac 12$ in the program +below, but make it so that playing with this parameter remains simple. + +Given this time discretization, space discretization happens as it always +does, by multiplying with test functions, integrating by parts, and then +restricting everything to a finite dimensional subspace. This yields the +following set of fully discrete equations after multiplying through with +$k_n$: +@f{align*} + M U^n-MU^{n-1} + + + k_n \left[ + (1-\theta)A U^{n-1} + + + \theta A U^n + \right] + &= + k_n + \left[ + (1-\theta)F^{n-1} + + + \theta F^n + \right], +@f} +where $M$ is the mass matrix and $A$ is the stiffness matrix that results from +discretizing the Laplacian. Bringing all known quantities to the right hand +side yields the linear system we have to solve in every step: +@f{align*} + (M + + + k_n \theta A) U^n + &= + MU^{n-1} + - + k_n + (1-\theta)A U^{n-1} + + + k_n + \left[ + (1-\theta)F^{n-1} + + + \theta F^n + \right]. +@f} +The linear system on the left hand side is symmetric and positive definite, so +we should have no trouble solving it with the Conjugate Gradient method. + +We can start the iteration above if we have the set of nodal coefficients +$U^0$ at the initial time. Here, we take the ones we get by interpolating the +initial values $u_0(\mathbf x)$ onto the mesh used for the first time step. We +will also need to choose a time step; we will here just choose it as fixed, +but clearly advanced simulators will want to choose it adaptively. We will +briefly come back to this in the results section +below. + + +

Adapting meshes for time dependent problems

+ +When solving the wave equation and its variants in the previous few programs, +we kept the mesh fixed. Just as for stationary equations, one can make a good +case that this is not the smartest approach and that significant savings can +be had by adapting the mesh. There are, however, significant difficulties +compared to the stationary case. Let us go through them in turn: -

Verifying whether the code is correct

+ + + + +

What could possibly go wrong? 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 @@ -113,3 +364,71 @@ $\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. + + + +

The testcase

+ +Solving the heat equation on a simple domain with a simple right hand side +almost always leads to solutions that are exceedingly boring, since they +become very smooth very quickly and then do not move very much any +more. Rather, we here solve the equation on the L-shaped domain with zero +Dirichlet boundary values and zero initial conditions, but as right hand side +we choose +@f{align*} + f(\mathbf x, t) + = + \left\{ + \begin{array}{ll} + \chi_1(\mathbf x) + & \text{if \(0\le t \le 0.2\tau\) or \(\tau\le t \le 1.2\tau\) or \(2\tau\le t + \le 2.2\tau\), etc} + \\ + \chi_2(\mathbf x) + & \text{if \(0.5\le t \le 0.7\tau\) or \(1.5\tau\le t \le 1.7\tau\) or \(2.5\tau\le t + \le 2.7\tau\), etc} + \\ + 0 + & \text{otherwise} + \end{array} + \right. +@f} +Here, +@f{align*} + \chi_1(\mathbf x) &= + \left\{ + \begin{array}{ll} + 1 + & \text{if \(x>0.5\) and \(y>-0.5\)} + \\ + 0 + & \text{otherwise} + \end{array} + \right. + \\ + \chi_2(\mathbf x) &= + \left\{ + \begin{array}{ll} + 1 + & \text{if \(x>-0.5\) and \(y>0.5\)} + \\ + 0 + & \text{otherwise} + \end{array} + \right. +@f} +In other words, in every period of lenght $\tau$, the right hand side first +flashes on in domain 1, then off completely, then on in domain 2, then off +completely again. This pattern is probably best observed via the little +animation of the solution shown in the results +section. + +If you interpret the heat equation as finding the spatially and temporally +variable temperature distribution of a conducting solid, then the test case +above corresponds to an L-shaped body where we keep the boundary at zero +temperature, and heat alternatingly in two parts of the domain. While heating +is in effect, the temperature rises in these places, after which it diffuses +and diminishes again. The point of these initial conditions is that they +provide us with a solution that has singularities both in time (when sources +switch on and off) as well as time (at the reentrant corner as well as at the +edges and corners of the regions where the source acts). diff --git a/deal.II/examples/step-26/doc/results.dox b/deal.II/examples/step-26/doc/results.dox index f4c6feefb5..6a401f850d 100644 --- a/deal.II/examples/step-26/doc/results.dox +++ b/deal.II/examples/step-26/doc/results.dox @@ -1 +1,100 @@

Results

+ +As in many of the tutorials, the actual output of the program matters less +than how we arrived there. Nonetheless, here it is: +@code +=========================================== +Number of active cells: 48 +Number of degrees of freedom: 65 + +Time step 1 at t=0.002 + 7 CG iterations. + +=========================================== +Number of active cells: 60 +Number of degrees of freedom: 81 + + +Time step 1 at t=0.002 + 7 CG iterations. + +=========================================== +Number of active cells: 105 +Number of degrees of freedom: 136 + + +Time step 1 at t=0.002 + 7 CG iterations. + +[...] + +Time step 249 at t=0.498 + 13 CG iterations. +Time step 250 at t=0.5 + 14 CG iterations. + +=========================================== +Number of active cells: 1803 +Number of degrees of freedom: 2109 +@endcode + +Maybe of more interest is a visualization of the solution and the mesh on wich +it was computed: + + + +The movie shows how the two sources switch on and off and how the mesh reacts +to this. It is quite obvious that the mesh as is is probably not the best we +could come up with. We'll get back to this in the next section. + + + +

Possibilities for extensions

+ +There are at least two areas where one can improve this program significantly: +adaptive time stepping and a better choice of the mesh. + +

Adaptive time stepping

+ +Having chosen an implicit time stepping scheme, we are not bound by any +CFL-like condition on the time step. Furthermore, because the time scales on +which change happens on a given cell in the heat equation are not bound to the +cells diameter (unlike the case with the wave equation, where we had a fixed +speed of information transport that couples the temporal and spatial scales), +we can choose the time step as we please. Or, better, choose it as we deem +necessary for accuracy. + +Looking at the solution, it is clear that the action does not happen uniformly +over time: a lot is changing around the time we switch on a source, things +become less dramatic once a source is on for a little while, and we enter a +long phase of decline when both sources are off. During these times, we could +surely get away with a larger time step than before without sacrificing too +much accuracy. + +The literature has many suggestions on how to choose the time step size +adaptively. Much can be learned, for example, from the way ODE solvers choose +their time steps. One can also be inspired by a posteriori error estimators +that can, ideally, be written in a way that the consist of a temporal and a +spatial contribution to the overall error. If the temporal one is too large, +we should choose a smaller time step. Ideas in this direction can be found, +for example, in the PhD thesis of a former principal developer of deal.II, +Ralf Hartmann, published by the University of Heidelberg, Germany, in 2002. + + +

Better refinement criteria

+ +If you look at the meshes in the movie above, it is clear that they are not +particularly well suited to the task at hand. In fact, they look rather +random. The underlying reason is that we only adapt the mesh once every fifth +time step, and only allow for a single refinement in these cases. Whenever a +source switches on, the solution had been very smooth in this area before and +the mesh was consequently rather coarse. This implies that the next time step +when we refine the mesh, we will get one refinement level more in this area, +and five time steps later another level, etc. But this is not enough: first, +we should refine immediately when a source switches on (after all, in the +current context we at least know what the right hand side is), and we should +allow for more than one refinement level. Of course, all of this can be done +using deal.II, it just requires a bit of algorithmic thinking in how to make +this work! diff --git a/deal.II/examples/step-26/step-26.cc b/deal.II/examples/step-26/step-26.cc index c7dd258c06..28bcb029fa 100644 --- a/deal.II/examples/step-26/step-26.cc +++ b/deal.II/examples/step-26/step-26.cc @@ -19,6 +19,8 @@ */ +// The program starts with the usual include files, all of which you should +// have seen before by now: #include #include #include @@ -51,12 +53,30 @@ #include +// Then the usual placing of all content of this program into a namespace and +// the importation of the deal.II namespace into the one we will work in: namespace Step26 { using namespace dealii; - + // @sect3{The HeatEquation class} + // + // The next piece is the declaration of the main class of this program. It + // follows the well trodden path of previous examples. If you have looked at + // step-6, for example, the only thing worth noting here is that we need to + // build two matrices (the mass and Laplace matrix) and keep the current and + // previous time step's solution. We then also need to store the current + // time, the size of the time step, and the number of the current time + // step. The last of the member variables denotes the theta parameter + // discussed in the introduction that allows us to treat the explicit and + // implicit Euler methods as well as the Crank-Nicolson method and other + // generalizations all in one program. + // + // As far as member functions are concerned, the only possible surprise is + // that the refine_mesh function takes arguments for the + // minimal and maximal mesh refinement level. The purpose of this is + // discussed in the introduction. template class HeatEquation { @@ -95,18 +115,26 @@ namespace Step26 + // @sect3{Equation data} + + // In the following classes and functions, we implement the various pieces + // of data that define this problem (right hand side and boundary values) + // that are used in this program and for which we need function objects. The + // right hand side is chosen as discussed at the end of the + // introduction. For boundary values, we choose zero values, but this is + // easily changed below. template - class RightHandSide: public Function + class RightHandSide : public Function { public: - RightHandSide() + RightHandSide () : Function(), period (0.2) {} - virtual double value(const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; private: const double period; @@ -115,8 +143,8 @@ namespace Step26 template - double RightHandSide::value(const Point &p, - const unsigned int component) const + double RightHandSide::value (const Point &p, + const unsigned int component) const { Assert (component == 0, ExcInternalError()); Assert (dim == 2, ExcNotImplemented()); @@ -145,22 +173,18 @@ namespace Step26 template - class BoundaryValues: public Function + class BoundaryValues : public Function { public: - BoundaryValues() - : - Function() - { - } - - virtual double value(const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; }; + + template - double BoundaryValues::value(const Point &/*p*/, - const unsigned int component) const + double BoundaryValues::value (const Point &/*p*/, + const unsigned int component) const { Assert(component == 0, ExcInternalError()); return 0; @@ -168,8 +192,11 @@ namespace Step26 + // @sect3{The HeatEquation implementation} + // + // The next step then is the implementation of the main class. template - HeatEquation::HeatEquation() + HeatEquation::HeatEquation () : fe(1), dof_handler(triangulation), -- 2.39.5