From: Wolfgang Bangerth Date: Mon, 13 Oct 2014 23:59:45 +0000 (-0500) Subject: Give it a once over. X-Git-Tag: v8.2.0-rc1~107^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=be478badd044e230ea57f5b9ddebb934e53b627b;p=dealii.git Give it a once over. No functional changes, just updates to documentation and style. --- diff --git a/examples/step-52/doc/intro.dox b/examples/step-52/doc/intro.dox index 2998354dec..ba3cc54946 100644 --- a/examples/step-52/doc/intro.dox +++ b/examples/step-52/doc/intro.dox @@ -2,7 +2,7 @@ This program was contributed by Bruno Turcksin and Damien Lebrun-Grandie. -Note: In order to run this program, deal.II must be configured to use +@note In order to run this program, deal.II must be configured to use the UMFPACK sparse direct solver. Refer to the ReadMe for instructions how to do this. @@ -45,8 +45,12 @@ $[0,b]\times[0,b]$ of the form \phi(x,t) = A\sin(\omega t)(bx-x^2). @f} 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 +any particular time, and all the error will be due to the time +discretization. We do this because it is then easy to observe the order of +convergence of the various time stepping schemes we will consider, without +having to separate spatial and temporal errors. + +We impose the following boundary conditions: homogeneous Dirichlet for $x=0$ and $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: @@ -54,9 +58,15 @@ in fact of the form stated above: 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$. +Because the solution is a sine in time, we know that the exact solution +satisfies $\phi\left(x,\pi\right) = 0$. Therefore, the error at time $t=\pi$ is simply the norm of the numerical -solution and is particularly easily evaluated. +solution, i.e., $\|e(\cdot,t=\pi)\|_{L_2} = \|\phi_h(\cdot,t=\pi)\|_{L_2}$, +and is particularly easily evaluated. In the code, we evaluate the $l_2$ norm +of the vector of nodal values of $\phi_h$ instead of the $L_2$ norm of the +associated spatial function, since the former is simpler to compute; however, +on uniform meshes, the two are just related by a constant and we can +consequently observe the temporal convergence order with either.

Runge-Kutta methods

@@ -85,7 +95,7 @@ spatially discretized equation has the form M\frac{dU}{dt} = f(t,U), @f} where $M$ is the mass matrix and $f(t,U)$ is the spatially discretized version -of the $q(t,u(x,t))$ (where $q$ is typically the place where spatial +of $q(t,u(x,t))$ (where $q$ is typically the place where spatial derivatives appear, but this is not of much concern for the moment given that we only consider time derivatives). In other words, this form fits the general scheme above if we write @@ -113,11 +123,16 @@ Butcher tableaus.) At the time of the writing of this tutorial, the methods implemented in deal.II can be divided in three categories:
    -
  1. explicit Runge-Kutta -
  2. embedded (or adaptive) Runge-Kutta -
  3. implicit Runge-Kutta +
  4. Explicit Runge-Kutta; in order for a method to be explicit, it is +necessary that in the formula above defining $k_i$, $k_i$ does not appear +on the right hand side. In other words, these methods have to satisfy +$a_{ii}=0, i=1,\ldots,s$. +
  5. Embedded (or adaptive) Runge-Kutta; we will discuss their properties below. +
  6. Implicit Runge-Kutta; this class of methods require the solution of a +possibly nonlinear system the stages $k_i$ above, i.e., they have +$a_{ii}\neq 0$ for at least one of the stages $i=1,\ldots,s$.
-Many well known time stepping schemes that one does not typically associated +Many well known time stepping schemes that one does not typically associate with the names Runge or Kutta can in fact be written in a way so that they, too, can be expressed in these categories. They oftentimes represent the lowest-order members of these families. @@ -159,13 +174,13 @@ form $\alpha y = f(t,y)$ for $y$ in each (sub-)timestep. Internally, this is done using a Newton-type method and, consequently, they require that the user provide functions that can 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$. +$\left(I-\tau M^{-1} \frac{\partial f}{\partial y}\right)^{-1}$ or equivalently +$\left(M - \tau \frac{\partial f}{\partial y}\right)^{-1} M$. The particular form of this operator results from the fact that 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 + \left(M - \tau \frac{\partial f}{\partial y}\right) \Delta y = -M h(t,y) @f} for some (given) $h(t,y)$. Implicit methods are @@ -174,12 +189,17 @@ course affect the accuracy of the solution, even if the numerical solution remains stable and bounded. Methods in this class include backward Euler, implicit midpoint, -Crank-Nicolson, and a two stage SDIRK. +Crank-Nicolson, and the two stage SDIRK method (short for "singly diagonally +implicit Runge-Kutta", a term coined to indicate that the diagonal elements +$a_{ii}$ defining the time stepping method are all equal; this property +allows for the Newton matrix $I-\tau M^{-1}\frac{\partial f}{\partial y}$ to +be re-used between stages because $\tau$ is the same every time).

Spatially discrete formulation

-By expanding the solution as always using shape functions $\psi_j$ and writing +By expanding the solution of our model problem +as always using shape functions $\psi_j$ and writing @f{eqnarray*} \phi_h(x,t) = \sum_j U_j(t) \psi_j(x), @f} @@ -197,13 +217,16 @@ where @f} See also step-24 and step-26 to understand how we arrive here. %Boundary terms are not necessary due to the chosen boundary conditions for -the current problem. To use the Runge-Kutta methods, we can then recast this +the current problem. To use the Runge-Kutta methods, we recast this as follows: @f{eqnarray*} -f(y) = -{\cal D}y - {\cal A}y + {\cal S} +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. +with its derivative, +@f{eqnarray*} +\frac{\partial f}{\partial y} = -{\cal D} - {\cal A}. +@f}

Notes on the testcase

@@ -219,7 +242,7 @@ important to remember the following issues:
  1. You will need to project the solution to the new mesh when the mesh is changed. Of course, the mesh - used should be the same at the beginning and at the end of the time step, + used should be the same from the beginning to the end of each time step, a question that arises because Runge-Kutta methods use multiple evaluations of the equations within each time step.
  2. You will need to update the mass matrix and its inverse every time the diff --git a/examples/step-52/doc/results.dox b/examples/step-52/doc/results.dox index 72a5cb0518..4c2d9ea0b0 100644 --- a/examples/step-52/doc/results.dox +++ b/examples/step-52/doc/results.dox @@ -1,10 +1,10 @@

    Results

    -The point of this program is less to show partcular results, but instead to -show how it is done and that we have already demonstrated using the code -above. Consequently, the output the program yields is relatively sparse and -consists only of the console output and the solutions given in -VTU format for visualization. +The point of this program is less to show particular results, but instead to +show how it is done. This we have already demonstrated simply by discussing +the code above. Consequently, the output the program yields is relatively +sparse and consists only of the console output and the solutions given in VTU +format for visualization. The console output contains both errors and, for some of the methods, the number of steps they performed: @@ -20,7 +20,7 @@ Implicit Midpoint: error=0.00862702 Crank-Nicolson: error=0.00862675 SDIRK: error=0.0042349 -Embedded explicit methods: +Embedded %explicit methods: Heun-Euler: error=0.0073012 steps performed=284 Bogacki-Shampine: error=0.000207511 diff --git a/examples/step-52/doc/tooltip b/examples/step-52/doc/tooltip index d35ea621a8..7754b44e6b 100644 --- a/examples/step-52/doc/tooltip +++ b/examples/step-52/doc/tooltip @@ -1 +1 @@ -Time-dependent diffusion equation. Neutron transport +Time-dependent diffusion equation. Time stepping methods. diff --git a/examples/step-52/step-52.cc b/examples/step-52/step-52.cc index 2c5cd4ed82..51c7d9e302 100644 --- a/examples/step-52/step-52.cc +++ b/examples/step-52/step-52.cc @@ -61,17 +61,18 @@ namespace Step52 { using namespace dealii; - // @sect3{Diffusion} - - // The next piece is the declaration of the main class. Most of the functions in - // this class are not new and have been explained in previous tutorials. The - // only interesting functions are evaluate_diffusion and - // id_minus_tau_J_inverse. evaluate_diffusion evaluates the - // diffusion equation, $M^{-1}(f(t,y))$, at a given time, for a given $\tau$ - // and a given $y$. id_minus_tau_J_inverse evaluates - // $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$ or - // equivalently $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$ at - // a given time, for a given $\tau$ and $y$. This function is needed when an + // @sect3{The Diffusion class} + + // The next piece is the declaration of the main class. Most of the + // functions in this class are not new and have been explained in previous + // tutorials. The only interesting functions are + // evaluate_diffusion() and + // id_minus_tau_J_inverse(). evaluate_diffusion() + // evaluates the diffusion equation, $M^{-1}(f(t,y))$, at a given time and a + // given $y$. id_minus_tau_J_inverse() evaluates $\left(I-\tau + // M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$ or equivalently + // $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$ at a given + // time, for a given $\tau$ and $y$. This function is needed when an // implicit method is used. class Diffusion { @@ -85,40 +86,45 @@ namespace Step52 void assemble_system(); - double get_source(double time,const Point<2> &point) const; + double get_source (const double time, + const Point<2> &point) const; - Vector evaluate_diffusion(const double time, const Vector &y) const; + Vector evaluate_diffusion (const double time, + const Vector &y) const; - Vector id_minus_tau_J_inverse(const double time, - const double tau, - const Vector &y); + Vector id_minus_tau_J_inverse (const double time, + const double tau, + const Vector &y); - void output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const; + void output_results (const unsigned int time_step, + TimeStepping::runge_kutta_method method) const; - // The next three functions are the driver for the explicit methods, the + // The next three functions are the drivers for the explicit methods, the // implicit methods, and the embedded explicit methods respectively. The // driver function for embedded explicit methods returns the number of - // steps executed. - void explicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); + // steps executed given that it only takes the number of time steps passed + // as an argument as a hint, but internally computed the optimal time step + // itself. + void explicit_method (const TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time); - void implicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); + void implicit_method (const TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time); - unsigned int embedded_explicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); + unsigned int embedded_explicit_method (const TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time); unsigned int fe_degree; double diffusion_coefficient; - double absorption_xs; + double absorption_cross_section; Triangulation<2> triangulation; @@ -146,17 +152,17 @@ namespace Step52 : fe_degree(2), diffusion_coefficient(1./30.), - absorption_xs(1.), + absorption_cross_section(1.), fe(fe_degree), dof_handler(triangulation) {} - // @sect5{Diffusion::setup_system} + // @sect4{Diffusion::setup_system} // Now, we create the constraint matrix and the sparsity pattern. Then, we // initialize the matrices and the solution vector. - void Diffusion::setup_system() + void Diffusion::setup_system () { dof_handler.distribute_dofs(fe); @@ -175,12 +181,18 @@ namespace Step52 - // @sect5{Diffusion::assemble_system} - // In this function, we compute - // $-\int D \nabla b_i \cdot \nabla b_j d\boldsymbol{r} - \int \Sigma_a b_i b_j d\boldsymbol{r}$ - // and the mass matrix $\int b_i b_j d\boldsymbol{r}$. The mass matrix is then - // inverted using a direct solver. - void Diffusion::assemble_system() + // @sect4{Diffusion::assemble_system} In this function, we + // compute $-\int D \nabla b_i \cdot \nabla b_j d\boldsymbol{r} - \int + // \Sigma_a b_i b_j d\boldsymbol{r}$ and the mass matrix $\int b_i b_j + // d\boldsymbol{r}$. The mass matrix is then inverted using a direct solver; + // the inverse_mass_matrix variable will then store the inverse + // of the mass matrix so that $M^{-1}$ can be applied to a vector using the + // vmult() function of that object. (Internally, UMFPACK does + // not really store the inverse of the matrix, but its LU factors; applying + // the inverse matrix is then equivalent to doing one forward and one + // backward solves with these two factors, which has the same complexity as + // applying an explicit inverse of the matrix). + void Diffusion::assemble_system () { system_matrix = 0.; mass_matrix = 0.; @@ -214,13 +226,16 @@ namespace Step52 for (unsigned int i=0; iget_dof_indices(local_dof_indices); @@ -234,32 +249,41 @@ namespace Step52 - // @sect5{Diffusion::get_source} + // @sect4{Diffusion::get_source} // - // In this function, the source for a given time and a given point is - // computed. - double Diffusion::get_source(double time,const Point<2> &point) const + // In this function, the source term of the equation for a given time and a + // given point is computed. + double Diffusion::get_source (const double time, + const Point<2> &point) const { - const double pi = 3.14159265358979323846; const double intensity = 10.; - const double frequency = pi/10.; + const double frequency = numbers::PI/10.; const double b = 5.; const double x = point(0); - double source = 0.; - source = intensity*(frequency*std::cos(frequency*time)*(b*x-x*x) + std::sin(frequency*time) * - (absorption_xs*(b*x-x*x)+2.*diffusion_coefficient)); - - return source; + return intensity* (frequency*std::cos(frequency*time)*(b*x-x*x) + + + std::sin(frequency*time) * (absorption_cross_section*(b*x-x*x) + + + 2.*diffusion_coefficient)); } - // @sect5{Diffusion:evaluate_diffusion} + // @sect4{Diffusion:evaluate_diffusion} // - // Now, the weak form of the diffusion equation is evaluated at a given - // time $t$ and for a given vector $y$. - Vector Diffusion::evaluate_diffusion(const double time, const Vector &y) const + // Next, we evaluate the weak form of the diffusion equation at a given time + // $t$ and for a given vector $y$. In other words, as outlined in the + // introduction, we evaluate $M^{-1}(-{\cal D}y - {\cal A}y + {\cal + // S})$. For this, we have to apply the matrix $-{\cal D} - {\cal A}$ + // (previously computed and stored in the variable + // system_matrix) to $y$ and then add the source term which we + // integrate as we usually do. (Integrating up the solution could be done + // using VectorTools::create_right_hand_side() if you wanted to save a few + // lines of code, or wanted to take advantage of doing the integration in + // parallel.) The result is then multiplied by $M^{-1}$. + Vector Diffusion::evaluate_diffusion (const double time, + const Vector &y) const { Vector tmp(dof_handler.n_dofs()); tmp = 0.; @@ -267,7 +291,8 @@ namespace Step52 const QGauss<2> quadrature_formula(fe_degree+1); - FEValues<2> fe_values(fe, quadrature_formula, + FEValues<2> fe_values(fe, + quadrature_formula, update_values | update_quadrature_points | update_JxW_values); @@ -290,18 +315,21 @@ namespace Step52 for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); - constraint_matrix.distribute_local_to_global(cell_source,local_dof_indices,tmp); + constraint_matrix.distribute_local_to_global(cell_source, + local_dof_indices, + tmp); } - Vector value(dof_handler.n_dofs()); inverse_mass_matrix.vmult(value,tmp); @@ -309,28 +337,30 @@ namespace Step52 } - // @sect5{Diffusion::id_minus_tau_J_inverse} + // @sect4{Diffusion::id_minus_tau_J_inverse} // // We compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$. This // is done in several steps: // - compute $M-\tau \frac{\partial f}{\partial y}$ - // - inverse the matrix to get $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1}$ + // - invert the matrix to get $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1}$ // - compute $tmp=My$ - // - compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp = \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My$. - Vector Diffusion::id_minus_tau_J_inverse(const double time, const double tau, - const Vector &y) + // - compute $z=\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp = \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My$ + // - return z. + Vector Diffusion::id_minus_tau_J_inverse (const double time, + const double tau, + const Vector &y) { - Vector tmp(dof_handler.n_dofs()); - Vector result(y); SparseDirectUMFPACK inverse_mass_minus_tau_Jacobian; mass_minus_tau_Jacobian.copy_from(mass_matrix); - mass_minus_tau_Jacobian.add(-tau,system_matrix); + mass_minus_tau_Jacobian.add(-tau, system_matrix); inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian); - mass_matrix.vmult(tmp,y); + Vector tmp(dof_handler.n_dofs()); + mass_matrix.vmult(tmp, y); + Vector result(y); inverse_mass_minus_tau_Jacobian.vmult(result,tmp); return result; @@ -338,71 +368,75 @@ namespace Step52 - // @sect5{Diffusion::output_results} + // @sect4{Diffusion::output_results} // - // We output the solution in vtu files. - void Diffusion::output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const + // The following function then outputs the solution in vtu files indexed by + // the number of the time step and the name of the time stepping method. Of + // course, the (exact) result should really be the same for all time + // stepping method, but the output here at least allows us to compare them. + void Diffusion::output_results (const unsigned int time_step, + TimeStepping::runge_kutta_method method) const { std::string method_name; switch (method) { - case TimeStepping::FORWARD_EULER : + case TimeStepping::FORWARD_EULER: { method_name = "forward_euler"; break; } - case TimeStepping::RK_THIRD_ORDER : + case TimeStepping::RK_THIRD_ORDER: { method_name = "rk3"; break; } - case TimeStepping::RK_CLASSIC_FOURTH_ORDER : + case TimeStepping::RK_CLASSIC_FOURTH_ORDER: { method_name = "rk4"; break; } - case TimeStepping::BACKWARD_EULER : + case TimeStepping::BACKWARD_EULER: { method_name = "backward_euler"; break; } - case TimeStepping::IMPLICIT_MIDPOINT : + case TimeStepping::IMPLICIT_MIDPOINT: { method_name = "implicit_midpoint"; break; } - case TimeStepping::SDIRK_TWO_STAGES : + case TimeStepping::SDIRK_TWO_STAGES: { method_name = "sdirk"; break; } - case TimeStepping::HEUN_EULER : + case TimeStepping::HEUN_EULER: { method_name = "heun_euler"; break; } - case TimeStepping::BOGACKI_SHAMPINE : + case TimeStepping::BOGACKI_SHAMPINE: { method_name = "bocacki_shampine"; break; } - case TimeStepping::DOPRI : + case TimeStepping::DOPRI: { method_name = "dopri"; break; } - case TimeStepping::FEHLBERG : + case TimeStepping::FEHLBERG: { method_name = "fehlberg"; break; } - case TimeStepping::CASH_KARP : + case TimeStepping::CASH_KARP: { method_name = "cash_karp"; break; } - default : + default: { break; } @@ -423,17 +457,18 @@ namespace Step52 } - // @sect5{Diffusion::explicit_method} - // This function is the driver for all the explicit method. It calls - // evolve_one_time_step which performs one time step. - // evolve_one_time_step needs to evaluate $M^{-1}(f(t,y))$, - // i.e, it needs evaluate_diffusion. Because - // evaluate_diffusion is a member function, it needs to be bind to - // $this$. Finally, the solution is output every 10 time steps. - void Diffusion::explicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time) + // @sect4{Diffusion::explicit_method} + // + // This function is the driver for all the explicit methods. It calls + // evolve_one_time_step which performs one time step. For + // explicit methods, evolve_one_time_step needs to evaluate + // $M^{-1}(f(t,y))$, i.e, it needs evaluate_diffusion. Because + // evaluate_diffusion is a member function, it needs to be bound + // to this. Finally, the solution is output every 10 time steps. + void Diffusion::explicit_method (const TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time) { const double time_step = (final_time-initial_time)/static_cast (n_time_steps); double time = initial_time; @@ -444,7 +479,10 @@ namespace Step52 for (unsigned int i=0; iDiffusion::implicit_method} + // @sect4{Diffusion::implicit_method} // This function is equivalent to explicit_method but for implicit // methods. When using implicit methods, we need to evaluate $M^{-1}(f(t,y))$ - // and $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$. - void Diffusion::implicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time) + // and $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$ + // for which we use the two member functions previously introduced. + void Diffusion::implicit_method (const TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time) { const double time_step = (final_time-initial_time)/static_cast (n_time_steps); double time = initial_time; @@ -472,8 +511,14 @@ namespace Step52 for (unsigned int i=0; iDiffusion::embedded_explicit_method} + // @sect4{Diffusion::embedded_explicit_method} // This function is the driver for the embedded explict methods. It requires // more parameters: // - coarsen_param: factor multiplying the current time step when the error @@ -499,10 +544,10 @@ namespace Step52 // is too large, the time step will be reduced. If the error is below the // threshold, a larger time step will be tried for the next time step. // delta_t_guess is the guessed time step produced by the embedded method. - unsigned int Diffusion::embedded_explicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time) + unsigned int Diffusion::embedded_explicit_method(const TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time) { double time_step = (final_time-initial_time)/static_cast (n_time_steps); double time = initial_time; @@ -514,14 +559,21 @@ namespace Step52 const double coarsen_tol = 1e-5; solution = 0.; - TimeStepping::EmbeddedExplicitRungeKutta > embedded_explicit_runge_kutta(method, - coarsen_param,refine_param,min_delta,max_delta,refine_tol,coarsen_tol); - output_results(0,method); + TimeStepping::EmbeddedExplicitRungeKutta > + embedded_explicit_runge_kutta(method, + coarsen_param, + refine_param, + min_delta, + max_delta, + refine_tol, + coarsen_tol); + output_results (0, method); + + // Now for the time loop. The last time step is chosen such that the final + // time is exactly reached. unsigned int n_steps=0; while (timefinal_time) time_step = final_time-time; @@ -541,15 +593,17 @@ namespace Step52 - // @sect5{Diffusion::run} - void Diffusion::run() + // @sect4{Diffusion::run} + // + // The following is the main function of the program. At the top, we create + // the grid (a [0,5]x[0,5] square) and refine it four times to get a mesh + // that has 16 by 16 cells, for a total of 256. We then set the boundary + // indicator to 1 for those parts of the boundary where $x=0$ and $x=5$. + void Diffusion::run () { - // We create the grid (a [0,5]x[0,5] square) and refine the mesh four times. - // The final grid has 16 by 16 cells, for a total of 256. GridGenerator::hyper_cube(triangulation, 0., 5.); triangulation.refine_global(4); - // Set the boundary indicator for x=0 and x=5 to 1. typename Triangulation<2>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); @@ -564,16 +618,23 @@ namespace Step52 cell->face(f)->set_boundary_indicator(0); } + // Next, we set up the linear systems and fill them with content so that + // they can be used throughout the time stepping process: setup_system(); assemble_system(); - unsigned int n_steps = 0; + // Finally, we solve the diffusion problem using several of the + // Runge-Kutta methods implemented in namespace TimeStepping, each time + // outputting the error at the end time. (As explained in the + // introduction, since the exact solution is zero at the final time, the + // error equals the numerical solution and can be computed by just taking + // the $l_2$ norm of the solution vector.) + unsigned int n_steps = 0; const unsigned int n_time_steps = 200; - const double initial_time = 0.; - const double final_time = 10.; + const double initial_time = 0.; + const double final_time = 10.; - // Next, we solve the diffusion problem using different Runge-Kutta methods. std::cout << "Explicit methods:" << std::endl; explicit_method (TimeStepping::FORWARD_EULER, n_time_steps,