From caef6830ccfa2cf79dcaf116c1207da1d9473714 Mon Sep 17 00:00:00 2001 From: turcksin Date: Tue, 20 May 2014 19:18:52 +0000 Subject: [PATCH] Improve documentation on step-52. git-svn-id: https://svn.dealii.org/trunk@32948 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-52/doc/intro.dox | 26 ++++---- deal.II/examples/step-52/step-52.cc | 90 +++++++++++++++++--------- 2 files changed, 73 insertions(+), 43 deletions(-) diff --git a/deal.II/examples/step-52/doc/intro.dox b/deal.II/examples/step-52/doc/intro.dox index 11d19584c0..e2706e98eb 100644 --- a/deal.II/examples/step-52/doc/intro.dox +++ b/deal.II/examples/step-52/doc/intro.dox @@ -22,12 +22,13 @@ augmented by appropriate boundary conditions. Here, $v$ is the velocity of neutrons, $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. In this -example, we are only interested in the error in time and thus, we are looking -for a solution of the form: +example, we are only interested in the error in time. The domain is square +$[0,b]\times[0,b] and we are looking for a solution of the form: @f{eqnarray*} \phi(x,t) = A\sin(\omega t)(bx-x^2). @f} -By using quadratic finite elements, we will not have any spatial error. We +By using quadratic finite elements, we will not have any spatial error and all +the error will come from the time discretization. We impose the following boundary conditions: homogeneous Dirichlet fo $x=0$ and $x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. The source is given by: @@ -42,7 +43,7 @@ found.

Runge-Kutta

-The Runke-Kutta methods implemented in deal.II assume that the equation to be +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). @@ -61,11 +62,11 @@ y_{n+1} = y_n + \sum_{i=1}^s b_i k_i @f} where @f{eqnarray*} -k_i = h M^{-1} f(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j) +k_i = h M^{-1} f\left(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j\right) @f} with $a_{ij}$, $b_i$, and $c_i$ are known coefficient and $h$ is the time step -used. The methods currently implemented in deal.II can be divided in three -categories: +used. 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 @@ -80,17 +81,18 @@ methods become unstable when the time step chosen is too large.

    Embedded Runge-Kutta

    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 it can be -coarsen. Only embedded explicit methods have been implemented so far. +estimate the error and decide if the time step needs to be refined or coarsen. +Only embedded explicit methods have been implemented at the time of the writing.

    Implicit Runge-Kutta

    These methods include backward Euler, implicit midpoint, Crank-Nicolson, and the 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) = \left(M - \Delta -t \frac{\partial f}{\partial y}\right)^{-1} M$. These methods are always stable. +$\left(I-\Delta t M^{-1} \frac{\partial f}{\partial Y}\right)$ or equivalently +$\left(M - \Deltat \frac{\partial f}{\partial y}\right)^{-1} M$. These methods are +always stable.

    Remarks

    -To simplify the problem, we solve the domain in two dimensional and the mesh is +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 diff --git a/deal.II/examples/step-52/step-52.cc b/deal.II/examples/step-52/step-52.cc index f13222ed70..f49dfb75d1 100644 --- a/deal.II/examples/step-52/step-52.cc +++ b/deal.II/examples/step-52/step-52.cc @@ -20,7 +20,7 @@ // @sect3{Include files} -// The first task as usal is to include the functionality of these well-known +// The first task as usual is to include the functionality of these well-known // deal.II library files and some C++ header files. #include #include @@ -70,6 +70,7 @@ namespace Step52 public: Diffusion(); + // This function is the driver that will call the other ones. void run(); private: @@ -83,24 +84,26 @@ namespace Step52 // Compute the intensity of the source at the given point. double get_source(double time,const Point<2> &point) const; - // Evaluate the diffusion equation $M^{-1}(f(t,y))$ + // Evaluate the diffusion equation $M^{-1}(f(t,y))$ at a given time and + // for a given y. Vector evaluate_diffusion(const double time, const Vector &y) const; - // Evaluate $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1} = - // \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M $ + // Evaluate $\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. Vector id_minus_tau_J_inverse(const double time, const double tau, const Vector &y); - // Output the results as vtu + // Output the results as vtu files. void output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const; - // Driver for the explicit methods + // Driver for the explicit methods. void explicit_method(TimeStepping::runge_kutta_method method, const unsigned int n_time_steps, const double initial_time, const double final_time); - // Driver for the implicit methods + // Driver for the implicit methods. void implicit_method(TimeStepping::runge_kutta_method method, const unsigned int n_time_steps, const double initial_time, @@ -114,6 +117,7 @@ namespace Step52 const double final_time); + // The next parameters are self-explanatory. unsigned int fe_degree; double diffusion_coefficient; @@ -140,7 +144,8 @@ namespace Step52 - // We choose quadratic finite elements so that there are no spatial error. + // We choose quadratic finite elements so there are no spatial error and we + // initialize the parameters. Diffusion::Diffusion() : fe_degree(2), @@ -152,6 +157,7 @@ namespace Step52 + // @sect5{Diffusion::setup_system} void Diffusion::setup_system() { dof_handler.distribute_dofs(fe); @@ -165,6 +171,7 @@ namespace Step52 DoFTools::make_sparsity_pattern(dof_handler,c_sparsity,constraint_matrix); sparsity_pattern.copy_from(c_sparsity); + // Initialize the matrices and the solution vector. system_matrix.reinit(sparsity_pattern); mass_matrix.reinit(sparsity_pattern); mass_minus_tau_Jacobian.reinit(sparsity_pattern); @@ -173,6 +180,7 @@ namespace Step52 + // @sect5{Diffusion::assemble_system} void Diffusion::assemble_system() { system_matrix = 0.; @@ -196,7 +204,8 @@ namespace Step52 cell = dof_handler.begin_active(), endc = dof_handler.end(); - // Compute $-\int D \nabla b \cdot \nabla b - \int \Sigma_a b b $ and $\int b b $ + // Compute $-\int D \nabla b \cdot \nabla b - \int \Sigma_a b b $ and the mass matrix + // $\int b b$. for (; cell!=endc; ++cell) { cell_matrix = 0.; @@ -229,6 +238,9 @@ namespace Step52 + // @sect5{Diffusion::get_source} + // + // Compute the source for a given time and a given point. double Diffusion::get_source(double time,const Point<2> &point) const { const double pi = 3.14159265358979323846; @@ -246,6 +258,9 @@ namespace Step52 + // @sect5{Diffusion:evaluate_diffusion} + // + // Evaluate the diffusion weak form give a time t and a vector y. Vector Diffusion::evaluate_diffusion(const double time, const Vector &y) const { Vector tmp(dof_handler.n_dofs()); @@ -253,8 +268,6 @@ namespace Step52 // Compute system_matrix*y system_matrix.vmult(tmp,y); - - // Compute the source term const QGauss<2> quadrature_formula(fe_degree+1); FEValues<2> fe_values(fe, quadrature_formula, @@ -300,7 +313,9 @@ namespace Step52 } - + // @sect5{Diffusion::id_minus_tau_J_inverse} + // + // We compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$. Vector Diffusion::id_minus_tau_J_inverse(const double time, const double tau, const Vector &y) { @@ -308,10 +323,17 @@ namespace Step52 Vector result(y); SparseDirectUMFPACK inverse_mass_minus_tau_Jacobian; + // Compute $M-\tau \frac{\partial f}{\partial y}$. mass_minus_tau_Jacobian.copy_from(mass_matrix); mass_minus_tau_Jacobian.add(-tau,system_matrix); + + // Inverse the matrix to get $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1}$ inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian); + + // Compute $tmp=My$. mass_matrix.vmult(tmp,y); + + // Compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp$ inverse_mass_minus_tau_Jacobian.vmult(result,tmp); return result; @@ -319,6 +341,7 @@ namespace Step52 + // @sect5{Diffusion::output_results} void Diffusion::output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const { std::string method_name; @@ -401,7 +424,7 @@ namespace Step52 } - + // sect5{Diffusion::explicit_method} void Diffusion::explicit_method(TimeStepping::runge_kutta_method method, const unsigned int n_time_steps, const double initial_time, @@ -415,7 +438,7 @@ namespace Step52 output_results(0,method); for (unsigned int i=0; iDiffusion::implicit_method} void Diffusion::implicit_method(TimeStepping::runge_kutta_method method, const unsigned int n_time_steps, const double initial_time, @@ -442,7 +466,7 @@ namespace Step52 output_results(0,method); for (unsigned int i=0; iDiffusion::embedded_explicit_method} unsigned int Diffusion::embedded_explicit_method(TimeStepping::runge_kutta_method method, const unsigned int n_time_steps, const double initial_time, @@ -491,7 +516,7 @@ namespace Step52 if (time+time_step>final_time) time_step = final_time-time; - // Because we use a member function, we need to bind this to the + // Because we use a member function, we need to bind $this$ to the // function. time = embedded_explicit_runge_kutta.evolve_one_time_step( std_cxx1x::bind(&Diffusion::evaluate_diffusion,this,std_cxx1x::_1,std_cxx1x::_2), @@ -511,14 +536,15 @@ namespace Step52 + // sect5{Diffusion::run} void Diffusion::run() { // Create the grid (a square [0,5]x[0,5]) and refine the mesh four times. - // The final gird has 16 times 16 cells, for a total of 256. + // The final gird 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 + // 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(); @@ -542,56 +568,58 @@ namespace Step52 const double initial_time = 0.; const double final_time = 10.; - // Use forward Euler + // Use forward Euler. explicit_method(TimeStepping::FORWARD_EULER,n_time_steps,initial_time,final_time); std::cout<<"Forward Euler error: "<main() function} +// // The following main function is similar to previous examples as // well, and need not be commented on. int main () -- 2.39.5