From: turcksin Date: Fri, 2 May 2014 20:25:11 +0000 (+0000) Subject: Start a tutorial using the Runge-Kutta methods implemented in TimeStepping.h X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f73a036aa4ba0a762aeb5f60d92b0b1cd2547110;p=dealii-svn.git Start a tutorial using the Runge-Kutta methods implemented in TimeStepping.h git-svn-id: https://svn.dealii.org/trunk@32878 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-52/CMakeLists.txt b/deal.II/examples/step-52/CMakeLists.txt new file mode 100644 index 0000000000..ec29e873c4 --- /dev/null +++ b/deal.II/examples/step-52/CMakeLists.txt @@ -0,0 +1,31 @@ +## +# CMake script for the step-52 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-52") + +# Declare all source files the target consists of: +SET(TARGET_SRC + ${TARGET}.cc + # You can specify additional files here! + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/deal.II/examples/step-52/doc/builds-on b/deal.II/examples/step-52/doc/builds-on new file mode 100644 index 0000000000..48a0f73876 --- /dev/null +++ b/deal.II/examples/step-52/doc/builds-on @@ -0,0 +1 @@ +step-4 diff --git a/deal.II/examples/step-52/doc/intro.dox b/deal.II/examples/step-52/doc/intro.dox new file mode 100644 index 0000000000..4ae88a593d --- /dev/null +++ b/deal.II/examples/step-52/doc/intro.dox @@ -0,0 +1,102 @@ +
+ +This program was contributed by Bruno Turcksin and Damien Lebrun-Grandie. + + +

Introducion

+ +This program shows how to use Runge-Kutta methods to solve a time-dependent +problem. + +

Problem statement

+ +In this example, we solve the energy-integrated 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 +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) +- \Sigma_a(x) \phi(x,t) + S(x,t) +@f} +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: +@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 +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: +@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 +\f$\phi\left(x,\frac{\pi}{\omega}\right) = 0$. Therefore, we can easily +compute the error at this time since it is simply the norm of the solution +found. + +

Runge-Kutta

+ +The Runke-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). +@f} +When using finite elements, the previous equation becomes: +@f{eqnarray*} +M\frac{dy}{dt} = f(t,y), +@f} +where $M$ is the mass matrix. Therefore, we have: +@f{eqnarray*} +\frac{dy}{dt} = M^{-1}f(t,y). +@f} +Runke-Kutta methods can be written as: +@f{eqnarray*} +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) +@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: +
    +
  1. explicit Runge-Kutta +
  2. embedded (or adaptive) Runge-Kutta +
  3. implicit Runge-Kutta +
+ +

Explicit Runge-Kutta

+These methods that include for 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. + +

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

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

Remarks

+To simplify the problem, we solve the domain in 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, we cannot forget to: +
    +
  1. 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. +
  2. update the mass matrix and its inverse. +
diff --git a/deal.II/examples/step-52/doc/kind b/deal.II/examples/step-52/doc/kind new file mode 100644 index 0000000000..86a44aa1ef --- /dev/null +++ b/deal.II/examples/step-52/doc/kind @@ -0,0 +1 @@ +time dependent diff --git a/deal.II/examples/step-52/doc/results.dox b/deal.II/examples/step-52/doc/results.dox new file mode 100644 index 0000000000..4a70bbe4bc --- /dev/null +++ b/deal.II/examples/step-52/doc/results.dox @@ -0,0 +1,30 @@ +

Results

+ +The output of this program consist of the console output and solutions given in +vtu format. + +The console output is: +@code +Forward Euler error: 1.00883 +Third order Runge-Kutta error: 0.000227982 +Fourth order Runge-Kutta error: 1.90541e-06 +Backward Euler error: 1.03428 +Implicit Midpoint error: 0.00862702 +Crank-Nicolson error: 0.00862675 +SDIRK error: 0.0042349 +Heun-Euler error: 0.0073012 +Number of steps done: 284 +Bogacki-Shampine error: 0.000207511 +Number of steps done: 200 +Dopri error: 4.01775e-09 +Number of steps done: 200 +Fehlberg error: 9.89504e-09 +Number of steps done: 200 +Cash-Karp error: 2.5579e-10 +Number of steps done: 200 +@endcode + +Like expected the high-order methods give a more accurate solutions. We see that +the Heun-Euler method adapted the number of time steps in order to satisfy the +tolerance. The others embedded methods did not need to change the number of time +steps. diff --git a/deal.II/examples/step-52/doc/tooltip b/deal.II/examples/step-52/doc/tooltip new file mode 100644 index 0000000000..d35ea621a8 --- /dev/null +++ b/deal.II/examples/step-52/doc/tooltip @@ -0,0 +1 @@ +Time-dependent diffusion equation. Neutron transport diff --git a/deal.II/examples/step-52/step-52.cc b/deal.II/examples/step-52/step-52.cc new file mode 100644 index 0000000000..1031e49cdd --- /dev/null +++ b/deal.II/examples/step-52/step-52.cc @@ -0,0 +1,629 @@ +/* --------------------------------------------------------------------- + * $Id: step-52.cc 30526 2013-08-29 20:06:27Z felix.gruber $ + * + * Copyright (C) 2014 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * + * Authors: Damien Lebrun-Grandie, Bruno Turcksin, 2014 + */ + +// @sect3{Include files} + +// The first task as usal is to include the functionality of these well-known +// deal.II library files and some C++ header files. +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +// This is the only include file that is new: It includes all the Runge-Kutta +// methods. +#include + + +// The next step is like in all previous tutorial programs: We put everything +// into a namespace of its own and then import the deal.II classes and functions +// into it. +namespace Step52 +{ + using namespace dealii; + + // @sect3{Diffusion} + + // Now, here comes the declaration of the main class. + class Diffusion + { + public: + Diffusion(); + + void run(); + + private: + // Create the sparsity_pattern and initialize system_matrix. + void setup_system(); + + // Assemble the part of the matrix of the system that does not depend on + // the time. + void assemble_system(); + + // Compute the intensity of the source at the given point. + double get_source(double time,const Point<2> &point) const; + + // Evaluate the diffusion equation \f$M^{-1}(f(t,y))\f$ + Vector evaluate_diffusion(const double time, const Vector &y) const; + + // Evaluate \f$\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 \f$ + Vector id_minus_tau_J_inverse(const double time, const double tau, + const Vector &y); + + // Output the results as vtu + void output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const; + + // 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 + void implicit_method(TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time); + + // Driver for the embedded explicit methods. Returns the number of steps + // executed. + 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 fe_degree; + + double diffusion_coefficient; + double absorption_xs; + + Triangulation<2> triangulation; + + FE_Q<2> fe; + + DoFHandler<2> dof_handler; + + ConstraintMatrix constraint_matrix; + + SparsityPattern sparsity_pattern; + + SparseMatrix system_matrix; + SparseMatrix mass_matrix; + SparseMatrix mass_minus_tau_Jacobian; + + SparseDirectUMFPACK inverse_mass_matrix; + + Vector solution; + }; + + + + // We choose quadratic finite elements so that there are no spatial error. + Diffusion::Diffusion() + : + fe_degree(2), + diffusion_coefficient(1./30.), + absorption_xs(1.), + fe(fe_degree), + dof_handler(triangulation) + {} + + + + void Diffusion::setup_system() + { + dof_handler.distribute_dofs(fe); + + // Create the constraint matrix. + VectorTools::interpolate_boundary_values(dof_handler,1,ZeroFunction<2>(),constraint_matrix); + constraint_matrix.close(); + + // Create the sparsity_pattern. + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler,c_sparsity,constraint_matrix); + sparsity_pattern.copy_from(c_sparsity); + + system_matrix.reinit(sparsity_pattern); + mass_matrix.reinit(sparsity_pattern); + mass_minus_tau_Jacobian.reinit(sparsity_pattern); + solution.reinit(dof_handler.n_dofs()); + } + + + + void Diffusion::assemble_system() + { + system_matrix = 0.; + mass_matrix = 0.; + + const QGauss<2> quadrature_formula(fe_degree+1); + + FEValues<2> fe_values(fe, quadrature_formula, + update_values | update_gradients | update_JxW_values); + + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler<2>::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + // Compute \f$-\int D \nabla b \cdot \nabla b - \int \Sigma_a b b\f$ and \f$\int b b\f$ + for (; cell!=endc; ++cell) + { + cell_matrix = 0.; + cell_mass_matrix = 0.; + + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + + constraint_matrix.distribute_local_to_global(cell_matrix,local_dof_indices,system_matrix); + constraint_matrix.distribute_local_to_global(cell_mass_matrix,local_dof_indices,mass_matrix); + } + + // Compute the inverse of the mass matrix. + inverse_mass_matrix.initialize(mass_matrix); + } + + + + double Diffusion::get_source(double time,const Point<2> &point) const + { + const double pi = 3.14159265358979323846; + const double intensity = 10.; + const double frequency = 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; + } + + + + Vector Diffusion::evaluate_diffusion(const double time, const Vector &y) const + { + Vector tmp(dof_handler.n_dofs()); + tmp = 0.; + // 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, + update_values | update_quadrature_points | update_JxW_values); + + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + Vector cell_source(dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler<2>::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + cell_source = 0.; + + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + + // Add the source term to the tmp vector. + constraint_matrix.distribute_local_to_global(cell_source,local_dof_indices,tmp); + } + + + Vector value(dof_handler.n_dofs()); + inverse_mass_matrix.vmult(value,tmp); + + return value; + } + + + + 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); + inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian); + mass_matrix.vmult(tmp,y); + inverse_mass_minus_tau_Jacobian.vmult(result,tmp); + + return result; + } + + + + void Diffusion::output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const + { + std::string method_name; + + switch (method) + { + case TimeStepping::FORWARD_EULER : + { + method_name = "forward_euler"; + break; + } + case TimeStepping::RK_THIRD_ORDER : + { + method_name = "rk3"; + break; + } + case TimeStepping::RK_CLASSIC_FOURTH_ORDER : + { + method_name = "rk4"; + break; + } + case TimeStepping::BACKWARD_EULER : + { + method_name = "backward_euler"; + break; + } + case TimeStepping::IMPLICIT_MIDPOINT : + { + method_name = "implicit_midpoint"; + break; + } + case TimeStepping::SDIRK_TWO_STAGES : + { + method_name = "sdirk"; + break; + } + case TimeStepping::HEUN_EULER : + { + method_name = "heun_euler"; + break; + } + case TimeStepping::BOGACKI_SHAMPINE : + { + method_name = "bocacki_shampine"; + break; + } + case TimeStepping::DOPRI : + { + method_name = "dopri"; + break; + } + case TimeStepping::FEHLBERG : + { + method_name = "fehlberg"; + break; + } + case TimeStepping::CASH_KARP : + { + method_name = "cash_karp"; + break; + } + default : + { + break; + } + } + + DataOut<2> data_out; + + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "flux"); + + data_out.build_patches(); + + const std::string filename = "solution-" + method_name + "-" + + Utilities::int_to_string (time_step, 3) + + ".vtu"; + std::ofstream output(filename.c_str()); + data_out.write_vtu(output); + } + + + + void Diffusion::explicit_method(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; + solution = 0.; + + TimeStepping::ExplicitRungeKutta > explicit_runge_kutta(method); + output_results(0,method); + for (unsigned int i=0; i (n_time_steps); + double time = initial_time; + solution = 0.; + + TimeStepping::ImplicitRungeKutta > implicit_runge_kutta(method); + output_results(0,method); + for (unsigned int i=0; i (n_time_steps); + double time = initial_time; + // Factor multiplying the current time step when the error is below the + // threshold. + const double coarsen_param = 1.2; + // Factor multiplying the current time step when the error is above the + // threshold. + const double refine_param = 0.8; + // Smallest time step acceptable. + const double min_delta = 1e-8; + // Largest time step acceptable. + const double max_delta = 10*time_step; + // Threshold above which the time step is refined. + const double refine_tol = 1e-1; + // Threshold below which the time step is coarsen. + 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); + unsigned int n_steps=0; + while (timefinal_time) + time_step = final_time-time; + + // 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), + time,time_step,solution); + + // We output the results every 10 time steps. + if ((n_steps+1)%10==0) + output_results(n_steps+1,method); + + // Update the time step + time_step = embedded_explicit_runge_kutta.get_status().delta_t_guess; + ++n_steps; + } + + return n_steps; + } + + + + 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. + 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(); + + for (; cell!=endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) + { + if ((cell->face(f)->center()[0]==0.) || (cell->face(f)->center()[0]==5.)) + cell->face(f)->set_boundary_indicator(1); + else + cell->face(f)->set_boundary_indicator(0); + } + + setup_system(); + + assemble_system(); + + unsigned int n_steps = 0; + const unsigned int n_time_steps = 200; + const double initial_time = 0.; + const double final_time = 10.; + + // Use forward Euler + explicit_method(TimeStepping::FORWARD_EULER,n_time_steps,initial_time,final_time); + std::cout<<"Forward Euler error: "<main function is similar to previous examples as +// well, and need not be commented on. +int main () +{ + try + { + Step52::Diffusion diffusion; + diffusion.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; + + return 0; +}