]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
More text, some fixes
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 8 Sep 2006 19:00:52 +0000 (19:00 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 8 Sep 2006 19:00:52 +0000 (19:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@13858 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-23/doc/intro.dox
deal.II/examples/step-23/doc/results.dox
deal.II/examples/step-23/step-23.cc

index 3bce82cb2a2346bf049bcbbf20150c337a144f7e..11072ca942f4a46f64ee8cd6051f50fcdbaa3625 100644 (file)
@@ -153,11 +153,11 @@ then reformulate the original wave equation as follows:
 The advantage of this formulation is that it now only contains first
 time derivatives for both variables, for which it is simple to write
 down time stepping schemes. Note that we do not have boundary
-conditions for $v$, a fact that is surprising at first as one would
-think that we could enforce $v=\frac{\partial g}{\partial t}$ on the
-boundary; however, it turns out that this leads to really bad
-approximations to the solution, as simple numerical experiments can
-show
+conditions for $v$ at first. However, we could enforce $v=\frac{\partial
+g}{\partial t}$ on the boundary. It turns out in numerical examples that this
+as actually necessary: without doing so the solution doesn't look particularly
+wrong, but the Crank-Nicolson scheme does not conserve energy if one doesn't
+enforce these boundary conditions
 
 With this formulation, let us introduce the following time
 discretization where a superscript $n$ indicates the number of a time
@@ -204,8 +204,9 @@ rearranging terms. We then get
 \f{eqnarray*}
   \left[ 1-k^2\theta^2\Delta \right] u^n &=&
         \left[ 1+k^2\theta(1-\theta)\Delta\right] u^{n-1} + k v^{n-1} 
-        - k^2\theta\left[\theta f^n + (1-\theta) f^{n-1}\right],\\
-   v^n &=& v^{n-1} + k\Delta\left[ \theta u^n + (1-\theta) u^{n-1}\right].
+        + k^2\theta\left[\theta f^n + (1-\theta) f^{n-1}\right],\\
+   v^n &=& v^{n-1} + k\Delta\left[ \theta u^n + (1-\theta) u^{n-1}\right]
+   + k\left[\theta f^n + (1-\theta) f^{n-1}\right].
 \f}
 In this form, we see that if we are given the solution
 $u^{n-1},v^{n-1}$ of the previous timestep, that we can the solve for
@@ -227,10 +228,10 @@ the entire domain, and integrate by parts where necessary. This leads
 to
 \f{eqnarray*}
   (u^n,\varphi) + k^2\theta^2(\nabla u^n,\nabla \varphi) &=&
-  (u^{n-1},\varphi) - k^2\theta^2(\nabla u^{n-1},\nabla \varphi)
+  (u^{n-1},\varphi) - k^2\theta(1-\theta)(\nabla u^{n-1},\nabla \varphi)
   +
   k(v^{n-1},\varphi)
-  - k^2\theta
+  + k^2\theta
   \left[
   \theta (f^n,\varphi) + (1-\theta) (f^{n-1},\varphi)
   \right],
@@ -240,7 +241,11 @@ to
    (v^{n-1},\varphi)
     - 
     k\left[ \theta (\nabla u^n,\nabla\varphi) + 
-    (1-\theta) (\nabla u^{n-1},\nabla \varphi)\right].
+    (1-\theta) (\nabla u^{n-1},\nabla \varphi)\right]
+  + k
+  \left[
+  \theta (f^n,\varphi) + (1-\theta) (f^{n-1},\varphi)
+  \right].
 \f}
 
 It is then customary to approximate $u^n(x) \approx u^n_h(x) = \sum_i
@@ -260,11 +265,11 @@ If we plug these expansions into above equations and test with the
 test functions from the present mesh, we get the following linear
 system: 
 \f{eqnarray*}
-  M^nU^n + k^2\theta^2 A^nU^n &=&
-  M^{n,n-1}U^{n-1} - k^2\theta^2 A^{n,n-1}U^{n-1}
+  (M^n + k^2\theta^2 A^n)U^n &=&
+  M^{n,n-1}U^{n-1} - k^2\theta(1-\theta) A^{n,n-1}U^{n-1}
   +
   kM^{n,n-1}V^{n-1}
-  - k^2\theta
+  + k^2\theta
   \left[
   \theta F^n + (1-\theta) F^{n-1}
   \right],
@@ -274,7 +279,11 @@ system:
    M^{n,n-1}V^{n-1}
     - 
     k\left[ \theta A^n U^n +
-    (1-\theta) A^{n,n-1} U^{n-1}\right],
+    (1-\theta) A^{n,n-1} U^{n-1}\right]
+   + k
+  \left[
+  \theta F^n + (1-\theta) F^{n-1}
+  \right],
 \f}
 where
 @f{eqnarray*}
@@ -304,23 +313,86 @@ integrals that contain shape functions defined on two meshes. This is a
 somewhat messy process that we omit here, but that is treated in some
 detail in @ref step_22 "step-22".
 
-As a final note, we remember that we have posed boundary values for
-$u$, but not $v$. In practice that means that we have to apply
-boundary values for the first equation above (i.e. the one for $U^n$),
-whereas we do not have to do that for the second one.
 
+<h3>Energy conservation</h3>
 
-<h3>Who are Courant, Friedrichs, and Levy?<h3>
-
-CFL condition
+One way to compare the quality of a time stepping scheme is to see whether the
+numerical approximation preserves conservation properties of the continuous
+equation. For the wave equation, the natural quantity to look at is the
+energy. By multiplying the wave equation by $u_t$, integrating over $\Omega$,
+and integrating by parts where necessary, we find that
+@f[
+       \frac{\partial}{\partial t}
+       \left[\frac 12 \int_\Omega \left(\frac{\partial u}{\partial
+       t}\right)^2 + (\nabla u)^2 \; dx\right]
+       =
+       \int_\Omega f \; dx
+       +
+       \int_{\partial\Omega} n\cdot\nabla u
+       \frac{\partial g}{\partial t} \; dx.
+@f]
+By consequence, in absence of body forces and constant boundary values, we get
+that
+@f[
+       E(t) = \frac 12 \int_\Omega \left(\frac{\partial u}{\partial
+       t}\right)^2 + (\nabla u)^2 \; dx
+@f]
+is a conserved quantity. We will compute this quantity after each time
+step. It is straightforward to see that if we replace $u$ by its finite
+element approximation, and $\frac{\partial u}{\partial t}$ by the finite
+element approximation of the velocity $v$, then
+@f[
+       E(t_n) = \frac 12 \left<V^n, M^n V^n\right>
+       +
+       \frac 12 \left<U^n, A^n U^n\right>.
+@f]
+As we will see in the results section, the Crank-Nicolson scheme does indeed
+conserve the energy, whereas neither the forward nor the backward Euler scheme
+do. 
 
-Energy conservation
 
-<h3>How the program works</h3>
+<h3>Who are Courant, Friedrichs, and Levy?</h3>
 
-Given the above formulation, ....
+One of the reasons why the wave equation is nasty to solve numerically is that
+explicit time discretizations are only stable if the time step is small
+enough. In particular, it is coupled to the spatial mesh width $h$. For the
+lowest order discretization we use here, the relationship reads
+@f[
+       k\le \frac ch
+@f]
+where $c$ is the wave speed, which in our formulation of the wave equation has
+been normalized to one. Consequently, unless we use the implicit schemes with
+$\theta>0$, our solutions will not be numerically stable if we violate this
+restriction. Implicit schemes do not have this restriction for stability, but
+they become inaccurate if the time step is too large.
+
+This condition was first recognized by Courant, Friedrichs, and Levy &mdash; in
+1943, long before computers become available for numerical computations! In
+the program, we will refine a square $[-1,1]$ seven times uniformly, giving a
+mesh size of $h=\frac 1{64}$, which is what we set the time step to. The fact
+that we set the time step and mesh size individually in two different places
+is error prone: it is too easy to refine the mesh once more but forget to also
+adjust the time step. @ref step_24 "step-24" shows a better way how to keep
+these things in synch.
 
 
 <h3>The testcase</h3>
 
-...
+Although the program has all the hooks to deal with nonzero initial and
+boundary conditions and body forces, we take a simple case where the domain is
+a square $[-1,1]^2$ and
+@f{eqnarray*}
+       f &=& 0,
+       \\
+       u_0 &=& 0,
+       \\
+       u_1 &=& 0,
+       \\
+       g &=& \left\{{\sin (4\pi t) \atop 0} 
+       \qquad {\textrm{for}\ t\le \frac 12, x=-1, -\frac 13<y<\frac 13 
+       \atop \textrm{otherwise}}
+       \right.
+@f}
+This corresponds to a membrane initially at rest and clamped all around, where
+someone is waving a part of the clamped boundary once up and down, thereby
+shooting a wave into the domain. 
index e67fccc1c1403f0a550ec7d6f60ab4918aaf24b5..f349855cdd9e25c4a2276bc79435c392ac482627 100644 (file)
@@ -1,2 +1,5 @@
 <a name="Results"></a>
 <h1>Results</h1>
+
+
+@image html step-23.movie.gif "Animation of the solution of step-23."
index db6a1ef6b22b4cf58037d6d8028f31b0dbb86ea1..ae2ce932ab237b5d5e7399f16ede884e146108e3 100644 (file)
@@ -180,17 +180,16 @@ class WaveEquation
 
                                 // @sect3{Equation data}
 
-                                // Before we go on filling in the
-                                // details of the main class, let us
-                                // define the equation data
-                                // corresponding to the problem,
-                                // i.e. initial and boundary values
-                                // as well as a right hand side
-                                // class. We do so using classes
-                                // derived from the Function class
-                                // template that has been used many
-                                // times before, so the following
-                                // should not be a surprise.
+                                // Before we go on filling in the details of
+                                // the main class, let us define the equation
+                                // data corresponding to the problem,
+                                // i.e. initial and boundary values for both
+                                // the solution $u$ as well as its time
+                                // derivative $v$, as well as a right hand
+                                // side class. We do so using classes derived
+                                // from the Function class template that has
+                                // been used many times before, so the
+                                // following should not be a surprise.
                                 //
                                 // Let's start with initial values
                                 // and choose zero for both the value
@@ -264,14 +263,15 @@ double RightHandSide<dim>::value (const Point<dim>  &/*p*/,
 
 
 
-                                // Finally, we have boundary
-                                // values. They are as described in
-                                // the introduction:
+                                // Finally, we have boundary values for $u$
+                                // and $v$. They are as described in the
+                                // introduction, one being the time
+                                // derivative of the other:
 template <int dim>
-class BoundaryValues : public Function<dim> 
+class BoundaryValuesU : public Function<dim> 
 {
   public:
-    BoundaryValues () : Function<dim>() {};
+    BoundaryValuesU () : Function<dim>() {};
     
     virtual double value (const Point<dim>   &p,
                          const unsigned int  component = 0) const;
@@ -281,8 +281,21 @@ class BoundaryValues : public Function<dim>
 
 
 template <int dim>
-double BoundaryValues<dim>::value (const Point<dim> &p,
-                                  const unsigned int component) const 
+class BoundaryValuesV : public Function<dim> 
+{
+  public:
+    BoundaryValuesV () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+
+
+template <int dim>
+double BoundaryValuesU<dim>::value (const Point<dim> &p,
+                                   const unsigned int component) const 
 {
   Assert (component == 0, ExcInternalError());
 
@@ -297,6 +310,24 @@ double BoundaryValues<dim>::value (const Point<dim> &p,
 
 
 
+template <int dim>
+double BoundaryValuesV<dim>::value (const Point<dim> &p,
+                                   const unsigned int component) const 
+{
+  Assert (component == 0, ExcInternalError());
+
+  if ((this->get_time() <= 0.5) &&
+      (p[0] < 0) &&
+      (p[1] < 1./3) &&
+      (p[1] > -1./3))
+    return (std::cos (this->get_time() * 4 * deal_II_numbers::PI) *
+           4 * deal_II_numbers::PI);
+  else
+    return 0;
+}
+
+
+
 
                                 // @sect3{Implementation of the <code>WaveEquation</code> class}
 
@@ -312,7 +343,10 @@ double BoundaryValues<dim>::value (const Point<dim> &p,
                                 // functions that solve linear systems, or
                                 // that generate output).
                                 //
-                                // Let's start with the constructor:
+                                // Let's start with the constructor (for an
+                                // explanation of the choice of time step,
+                                // see the section on Courant, Friedrichs,
+                                // and Levy in the introduction):
 template <int dim>
 WaveEquation<dim>::WaveEquation () :
                 fe (1),
@@ -341,13 +375,13 @@ void WaveEquation<dim>::setup_system ()
            << std::endl
            << "Total number of cells: "
            << triangulation.n_cells()
-           << std::endl
            << std::endl;
 
   dof_handler.distribute_dofs (fe);
 
-  std::cout << "   Number of degrees of freedom: "
+  std::cout << "Number of degrees of freedom: "
            << dof_handler.n_dofs()
+           << std::endl
            << std::endl;
 
   sparsity_pattern.reinit (dof_handler.n_dofs(),
@@ -543,21 +577,21 @@ void WaveEquation<dim>::run ()
                                   // ($T=5$ in this case). In each time step,
                                   // we first have to solve for $U^n$, using
                                   // the equation $(M^n + k^2\theta^2 A^n)U^n
-                                  // = M^{n,n-1}U^{n-1} - k^2\theta^2
-                                  // A^{n,n-1}U^{n-1} + kM^{n,n-1}V^{n-1} -
-                                  // k^2\theta \left[ \theta F^n + (1-\theta)
+                                  // =$ $(M^{n,n-1} - k^2\theta(1-\theta)
+                                  // A^{n,n-1})U^{n-1} + kM^{n,n-1}V^{n-1} +$
+                                  // $k\theta \left[k \theta F^n + k(1-\theta)
                                   // F^{n-1} \right]$. Note that we use the
                                   // same mesh for all time steps, so that
                                   // $M^n=M^{n,n-1}=M$ and
                                   // $A^n=A^{n,n-1}=A$. What we therefore
                                   // have to do first is to add up $MU^{n-1}
-                                  // - k^2\theta^2 AU^{n-1} + kMV^{n-1}$ and
-                                  // put the result into the
-                                  // <code>system_rhs</code> vector. (For
-                                  // these additions, we need a temporary
-                                  // vector that we declare before the loop
-                                  // to avoid repeated memory allocations in
-                                  // each time step.)
+                                  // - k^2\theta(1-\theta) AU^{n-1} + kMV^{n-1}$ and
+                                  // the forcing terms, and put the result
+                                  // into the <code>system_rhs</code>
+                                  // vector. (For these additions, we need a
+                                  // temporary vector that we declare before
+                                  // the loop to avoid repeated memory
+                                  // allocations in each time step.)
                                   //
                                   // The one thing to realize here is how we
                                   // communicate the time variable to the
@@ -578,6 +612,7 @@ void WaveEquation<dim>::run ()
                                   // function for all times at any given
                                   // spatial location.
   Vector<double> tmp (solution_u.size());
+  Vector<double> forcing_terms (solution_u.size());
   
   for (timestep_number=1, time=time_step;
        time<=5;
@@ -599,12 +634,17 @@ void WaveEquation<dim>::run ()
       rhs_function.set_time (time);
       VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
                                           rhs_function, tmp);
-      system_rhs.add (theta * theta * time_step * time_step, tmp);
+      forcing_terms = tmp;
+      forcing_terms *= theta * time_step;
+      
+      system_rhs.add (theta * time_step, forcing_terms);
 
       rhs_function.set_time (time-time_step);
       VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
                                           rhs_function, tmp);
-      system_rhs.add (theta * (1-theta) * time_step * time_step, tmp);
+
+      forcing_terms.add ((1-theta) * time_step, tmp);
+      system_rhs.add (theta * time_step, forcing_terms);
 
                                       // After so constructing the right hand
                                       // side vector of the first equation,
@@ -618,29 +658,33 @@ void WaveEquation<dim>::run ()
                                       // boundary values as we usually
                                       // do. The result is then handed off to
                                       // the solve_u() function:
-      BoundaryValues<dim> boundary_values_function;
-      boundary_values_function.set_time (time);
+      {
+       BoundaryValuesU<dim> boundary_values_u_function;
+       boundary_values_u_function.set_time (time);
       
-      std::map<unsigned int,double> boundary_values;
-      VectorTools::interpolate_boundary_values (dof_handler,
-                                               0,
-                                               boundary_values_function,
-                                               boundary_values);
-      MatrixTools::apply_boundary_values (boundary_values,
-                                         system_matrix,
-                                         solution_u,
-                                         system_rhs);
+       std::map<unsigned int,double> boundary_values;
+       VectorTools::interpolate_boundary_values (dof_handler,
+                                                 0,
+                                                 boundary_values_u_function,
+                                                 boundary_values);
+       MatrixTools::apply_boundary_values (boundary_values,
+                                           system_matrix,
+                                           solution_u,
+                                           system_rhs);
+      }
       solve_u ();
 
 
-                                      // The second step -- solving for $V^n$
-                                      // works similarly, except that this
-                                      // time the matrix on the left is the
-                                      // mass matrix, the right hand side is
-                                      // $MV^{n-1} - k\left[ \theta A U^n +
-                                      // (1-\theta) AU^{n-1}\right]$, and
-                                      // there are no boundary values to be
-                                      // applied.
+                                      // The second step, i.e. solving for
+                                      // $V^n$, works similarly, except that
+                                      // this time the matrix on the left is
+                                      // the mass matrix, the right hand side
+                                      // is $MV^{n-1} - k\left[ \theta A U^n
+                                      // + (1-\theta) AU^{n-1}\right]$ plus
+                                      // forcing terms. Boundary values are
+                                      // applied in the same way as before,
+                                      // except that now we have to use the
+                                      // BoundaryValuesV class:
       laplace_matrix.vmult (system_rhs, solution_u);
       system_rhs *= -theta * time_step;
 
@@ -650,27 +694,46 @@ void WaveEquation<dim>::run ()
       laplace_matrix.vmult (tmp, old_solution_u);
       system_rhs.add (-time_step * (1-theta), tmp);
 
-      rhs_function.set_time (time);
-      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
-                                          rhs_function, tmp);
-      system_rhs.add (theta * time_step, tmp);
-
-      rhs_function.set_time (time-time_step);
-      VectorTools::create_right_hand_side (dof_handler, QGauss<dim>(2),
-                                          rhs_function, tmp);
-      system_rhs.add ((1-theta) * time_step, tmp);
+      system_rhs += forcing_terms;
 
+      {
+       BoundaryValuesV<dim> boundary_values_v_function;
+       boundary_values_v_function.set_time (time);
+      
+       std::map<unsigned int,double> boundary_values;
+       VectorTools::interpolate_boundary_values (dof_handler,
+                                                 0,
+                                                 boundary_values_v_function,
+                                                 boundary_values);
+       MatrixTools::apply_boundary_values (boundary_values,
+                                           mass_matrix,
+                                           solution_v,
+                                           system_rhs);
+      }
       solve_v ();
 
                                       // Finally, after both solution
                                       // components have been computed, we
-                                      // output the result, and go on to the
-                                      // next time step after shifting the
-                                      // present solution into the vectors
-                                      // that hold the solution at the
-                                      // previous time step:
+                                      // output the result, compute the
+                                      // energy in the solution, and go on to
+                                      // the next time step after shifting
+                                      // the present solution into the
+                                      // vectors that hold the solution at
+                                      // the previous time step. Note the
+                                      // function
+                                      // SparseMatrix::matrix_norm_square
+                                      // that can compute
+                                      // $\left<V^n,MV^n\right>$ and
+                                      // $\left<U^n,AU^n\right>$ in one step,
+                                      // saving us the expense of a temporary
+                                      // vector and several lines of code:
       output_results ();
 
+      std::cout << "  Total energy: "
+               << (mass_matrix.matrix_norm_square (solution_v) +
+                   laplace_matrix.matrix_norm_square (solution_u)) / 2
+               << std::endl;
+      
       old_solution_u = solution_u;
       old_solution_v = solution_v;
     }

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.