From: bangerth Date: Fri, 8 Sep 2006 19:00:52 +0000 (+0000) Subject: More text, some fixes X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0418264b4375a1a46b87ac84b471091a93bc5b9e;p=dealii-svn.git More text, some fixes git-svn-id: https://svn.dealii.org/trunk@13858 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-23/doc/intro.dox b/deal.II/examples/step-23/doc/intro.dox index 3bce82cb2a..11072ca942 100644 --- a/deal.II/examples/step-23/doc/intro.dox +++ b/deal.II/examples/step-23/doc/intro.dox @@ -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. +

Energy conservation

-

Who are Courant, Friedrichs, and Levy?

- -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 + + + \frac 12 \left. +@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 -

How the program works

+

Who are Courant, Friedrichs, and Levy?

-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 — 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.

The testcase

-... +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

Results

+ + +@image html step-23.movie.gif "Animation of the solution of step-23." diff --git a/deal.II/examples/step-23/step-23.cc b/deal.II/examples/step-23/step-23.cc index db6a1ef6b2..ae2ce932ab 100644 --- a/deal.II/examples/step-23/step-23.cc +++ b/deal.II/examples/step-23/step-23.cc @@ -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::value (const Point &/*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 -class BoundaryValues : public Function +class BoundaryValuesU : public Function { public: - BoundaryValues () : Function() {}; + BoundaryValuesU () : Function() {}; virtual double value (const Point &p, const unsigned int component = 0) const; @@ -281,8 +281,21 @@ class BoundaryValues : public Function template -double BoundaryValues::value (const Point &p, - const unsigned int component) const +class BoundaryValuesV : public Function +{ + public: + BoundaryValuesV () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + + +template +double BoundaryValuesU::value (const Point &p, + const unsigned int component) const { Assert (component == 0, ExcInternalError()); @@ -297,6 +310,24 @@ double BoundaryValues::value (const Point &p, +template +double BoundaryValuesV::value (const Point &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 WaveEquation class} @@ -312,7 +343,10 @@ double BoundaryValues::value (const Point &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 WaveEquation::WaveEquation () : fe (1), @@ -341,13 +375,13 @@ void WaveEquation::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::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 - // system_rhs 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 system_rhs + // 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::run () // function for all times at any given // spatial location. Vector tmp (solution_u.size()); + Vector forcing_terms (solution_u.size()); for (timestep_number=1, time=time_step; time<=5; @@ -599,12 +634,17 @@ void WaveEquation::run () rhs_function.set_time (time); VectorTools::create_right_hand_side (dof_handler, QGauss(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(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::run () // boundary values as we usually // do. The result is then handed off to // the solve_u() function: - BoundaryValues boundary_values_function; - boundary_values_function.set_time (time); + { + BoundaryValuesU boundary_values_u_function; + boundary_values_u_function.set_time (time); - std::map 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 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::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(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(2), - rhs_function, tmp); - system_rhs.add ((1-theta) * time_step, tmp); + system_rhs += forcing_terms; + { + BoundaryValuesV boundary_values_v_function; + boundary_values_v_function.set_time (time); + + std::map 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$ and + // $\left$ 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; }