]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Start to look over step-52. Adjust step-26 accordingly.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 22 Jun 2014 11:27:56 +0000 (11:27 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 22 Jun 2014 11:27:56 +0000 (11:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@33071 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-26/doc/intro.dox
deal.II/examples/step-26/doc/results.dox
deal.II/examples/step-52/doc/intro.dox

index 7aa1ab6697eb958ba931de5236d3fe28a149cea5..11e9fe0e42a7869983469578915d8de0a99f1fb4 100644 (file)
@@ -47,7 +47,8 @@ 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.
+below, but make it so that playing with this parameter remains simple. (If you
+are interested in playing with higher order methods, take a look at step-52.)
 
 Given this time discretization, space discretization happens as it always
 does, by multiplying with test functions, integrating by parts, and then
index 798cc8f298e2906a1c15440177df8ecbac71b106..58a105902d162b838b95c3da086cec7bb9769e2a 100644 (file)
@@ -83,6 +83,15 @@ 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.
 
 
+<h4>Better time stepping methods</h4>
+
+We here use one of the simpler time stepping methods, namely the second order
+in time Crank-Nicolson method. However, more accurate methods such as
+Runge-Kutta methods are available and should be used as they do not represent
+much additional effort. It is not difficult to implement this for the current
+program, but a more systematic treatment is also given in step-52.
+
+
 <h4>Better refinement criteria</h4>
 
 If you look at the meshes in the movie above, it is clear that they are not
index 926c81cdb42e644c4ce97201bd0c858b2f03c6e4..d552ba90a8432f56d216625a79fd53b00fb18cb6 100644 (file)
@@ -12,46 +12,63 @@ problem.
 
 In this example, we solve the one-group time-dependent diffusion
 approximation of the neutron transport equation (see step-28 for the
-time-independent multigroup diffusion). We assume that the medium is not
+time-independent multigroup diffusion). This is a model for how neutrons move
+around highly scattering media, and consequently it is a variant of the
+time-dependent diffusion equation -- which is just a different name for the
+heat equation discussed in step-26, plus some extra terms.
+We assume that the medium is not
 fissible and therefore, the neutron flux satisfies the following equation:
 @f{eqnarray*}
-\frac{1}{v}\frac{\partial \phi(x,t)}{\partial t} = \nabla D(x) \nabla \phi(x,t)
+\frac{1}{v}\frac{\partial \phi(x,t)}{\partial t} = \nabla \cdot D(x) \nabla \phi(x,t)
 - \Sigma_a(x) \phi(x,t) + S(x,t)
 @f}
 augmented by appropriate boundary conditions. Here, $v$ is the velocity of
-neutrons (for simplicity we assume it is equal to 1), $D$ is the diffusion coefficient, 
+neutrons (for simplicity we assume it is equal to 1 which can be achieved by
+simply scaling the time variable), $D$ is the diffusion coefficient, 
 $\Sigma_a$ is the absorption cross section, and $S$ is a source. Because we are
 only interested in the time dependence, we assume that $D$ and $\Sigma_a$ are
-constant. We are looking for a solution on a square domain $[0,b]\times[0,b]$ of
-the form:
+constant. 
+
+Since this program only intends to demonstrate how to use advanced time
+stepping algorithms, we will only look for the solutions of relatively simple
+problems. Specifically, we are looking for a solution on a square domain
+$[0,b]\times[0,b]$ of the form
 @f{eqnarray*}
 \phi(x,t) = A\sin(\omega t)(bx-x^2).
 @f}
-By using quadratic finite elements, all the error will be due to the time discretization. We
+By using quadratic finite elements, we can represent this function exactly at
+any particular time, and all the error will be due to the time discretization. We
 impose the following boundary conditions: homogeneous Dirichlet for $x=0$ and
-$x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. The source is
-given by:
+$x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. We choose the
+source term so that the corresponding solution is
+in fact of the form stated above:
 @f{eqnarray*}
 S=A\left(\frac{1}{v}\omega \cos(\omega t)(bx -x^2) + \sin(\omega t)
 \left(\Sigma_a (bx-x^2)+2D\right) \right).
 @f}
 Because the solution is a sine, we know that $\phi\left(x,\pi\right) = 0$.
-Therefore, the error at this time is simply the norm of the numerical solution.
+Therefore, the error at time $t=\pi$ is simply the norm of the numerical
+solution and is particularly easily evaluated.
+
 
 <h3>Runge-Kutta</h3>
 
 The Runge-Kutta methods implemented in deal.II assume that the equation to be
 solved can be written as:
 @f{eqnarray*}
-\frac{dy}{dt} = f(t,y).
+\frac{dy}{dt} = g(t,y).
 @f}
-When using finite elements, the previous equation becomes:
+When using finite elements, discretized time derivatives always result in the
+presence of a mass matrix. If the solution vector $y(t)$ is given by the vector
+of nodal coefficients $U(t)$ the previous equation, then spatially discretized
+equations always have the form
 @f{eqnarray*}
-M\frac{dy}{dt} = f(t,y),
+M\frac{dU}{dt} = f(t,U),
 @f}
-where $M$ is the mass matrix. Therefore, we have:
+where $M$ is the mass matrix. In other words, this fits the scheme above if we
+write 
 @f{eqnarray*}
-\frac{dy}{dt} = M^{-1}f(t,y).
+\frac{dy}{dt} = g(t,y) = M^{-1}f(t,y).
 @f}
 Runke-Kutta methods can be written as:
 @f{eqnarray*}
@@ -62,7 +79,12 @@ where
 k_i = h M^{-1} f\left(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j\right)
 @f}
 where $a_{ij}$, $b_i$, and $c_i$ are known coefficients and $h$ is the time step
-used. At the time of the writing of this tutorial, the methods implemented in
+used. Different time stepping methods of the Runge-Kutta class differ in the
+number of stages $s$ and the values they use for the coefficients $a_{ij}$,
+$b_i$, and $c_i$ but are otherwise easy to implement since one can look up
+tabulated values for these coefficients.
+
+At the time of the writing of this tutorial, the methods implemented in
 deal.II can be divided in three categories:
 <ol>
 <li> explicit Runge-Kutta
@@ -70,42 +92,71 @@ deal.II can be divided in three categories:
 <li> implicit Runge-Kutta
 </ol> 
 
-<h4>Explicit Runge-Kutta</h4>
-These methods that include forward Euler, third order Runge-Kutta, and
-fourth order Runge-Kutta, require a function to evaluate $M^{-1}f(t,y)$. These
-methods become unstable when the time step chosen is too large.
+<h4>Explicit Runge-Kutta</h4> 
+These methods that include forward Euler, third order Runge-Kutta, and fourth
+order Runge-Kutta, only require a function to evaluate $M^{-1}f(t,y)$ but not
+(as implicit methods) to solve an equation for $y$ that involves
+$f(t,y)$. These methods become unstable when the time step chosen is too
+large.
 
 <h4>Embedded Runge-Kutta</h4>
 These methods include Heun-Euler, Bogacki-Shampine, Dormand-Prince (ode45 in
 Matlab), Fehlberg, and Cash-Karp. These methods use a low order method to
-estimate the error and decide if the time step needs to be refined or coarsen. 
+estimate the error and decide if the time step needs to be refined or coarsened
 Only embedded explicit methods have been implemented at the time of the writing.
 
 <h4>Implicit Runge-Kutta</h4>
 These methods include backward Euler, implicit midpoint, Crank-Nicolson, and a
 two stages SDIRK. These methods require to evaluate $M^{-1}f(t,y)$ and
 $\left(I-\Delta t M^{-1} \frac{\partial f}{\partial y}\right)^{-1}$ or equivalently 
-$\left(M - \Delta t \frac{\partial f}{\partial y}\right)^{-1} M$. These methods are 
-always stable.
+$\left(M - \Delta t \frac{\partial f}{\partial y}\right)^{-1} M$. This is necessary in order to solve for the solution of (possibly nonlinear)
+systems in every time step, where each Newton step requires the solution of an
+equation of the form
+@f{align*}
+  \left(M - \Delta t \frac{\partial f}{\partial y}\right) \Delta y
+  = -M h(t,y)
+@f}
+for some (given) $h(t,y)$. These methods are 
+always stable
 
-<h3>Weak form</h3>
-To use the Runge-Kutta methods, we need to be able to evaluate:
+
+<h3>Spatially discrete formulation</h3>
+
+By expanding the solution as always using shape functions $\psi_j$ and writing
+@f{eqnarray*}
+\phi_h(x,t) = \sum_j U_j(t) \psi_j(x),
+@f}
+we can then get the spatially discretized version of the diffusion equation as
+@f{eqnarray*}
+  M \frac{dU(t)}{dt}
+  = -{\cal D} U(t) - {\cal A} U(t) + {\cal S}(t)
+@f}
+where
 @f{eqnarray*}
-f = \oint D b_i \frac{\partial b_j}{\partial n} d\boldsymbol{r} - \int D \nabla
-b_i \nabla b_j \phi_j d\boldsymbol{r} -\int \Sigma_a b_i b_j \phi_j 
-d\boldsymbol{r} + \int b_j S d\boldsymbol{r}
+  M_{ij}  &= (\psi_i,\psi_j), \\
+  {\cal D}_{ij}  &= (D\nabla\psi_i,\nabla\psi_j)_\Omega, \\
+  {\cal A}_{ij}  &= (\Sigma_a\psi_i,\Sigma_a\psi_j)_\Omega, \\
+  {\cal S}_{i}(t)  &= (\psi_i,S(x,t))_\Omega.
 @f}
-and $\frac{\partial f}{\partial y}$. Because $f$ is linear in $y$ (or $\phi$ in
-this case) $\frac{\partial f}{\partial y} y = f$.
+Boundary terms are not necessary due to the chosen boundary conditions.
+To use the Runge-Kutta methods, we can then recast this as follows:
+@f{eqnarray*}
+f(y) = -{\cal D)y - {\cal A}y + {\cal S}
+@f}
+In the code, we will need to be able to evaluate this function $f(U)$ along
+with its derivative. However, in view of the linearity of $f$, we will be able
+to use $\frac{\partial f}{\partial y} y = f(y)$.
+
 
 <h3>Remarks</h3>
 To simplify the problem, the domain is two dimensional and the mesh is
 uniform (there is no need to adapt the mesh since we use quadratic finite
 elements and the exact solution is quadratic). Going from a two dimensional
 domain to a three dimensional domain is not very challenging. However if the
-mesh must be adapted, it is important to remember to:
+mesh must be adapted, it is important to remember to do the following:
 <ol>
-<li> project the solution to the new mesh when the mesh is changed. The mesh
-used should be the same at the beginning and at the end of the time step.
+<li> Project the solution to the new mesh when the mesh is changed. The mesh
+     used should be the same at the beginning and at the end of the time step.
 <li> update the mass matrix and its inverse.
 </ol>
+The techniques to do all of this are available in step-26.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.