<i>This program was contributed by Bruno Turcksin and Damien Lebrun-Grandie.</i>
-<b>Note:</b> 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 <a
href="../../readme.html#umfpack">ReadMe</a> for instructions how to do this.
\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:
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.
<h3>Runge-Kutta methods</h3>
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
At the time of the writing of this tutorial, the methods implemented in
deal.II can be divided in three categories:
<ol>
-<li> explicit Runge-Kutta
-<li> embedded (or adaptive) Runge-Kutta
-<li> implicit Runge-Kutta
+<li> 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$.
+<li> Embedded (or adaptive) Runge-Kutta; we will discuss their properties below.
+<li> 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$.
</ol>
-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.
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
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).
<h3>Spatially discrete formulation</h3>
-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}
@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}
<h3>Notes on the testcase</h3>
<ol>
<li> 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.
<li> You will need to update the mass matrix and its inverse every time the
{
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 <code>evaluate_diffusion</code> and
- // <code>id_minus_tau_J_inverse</code>. <code>evaluate_diffusion</code> evaluates the
- // diffusion equation, $M^{-1}(f(t,y))$, at a given time, for a given $\tau$
- // and a given $y$. <code>id_minus_tau_J_inverse</code> 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 <code>Diffusion</code> 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
+ // <code>evaluate_diffusion()</code> and
+ // <code>id_minus_tau_J_inverse()</code>. <code>evaluate_diffusion()</code>
+ // evaluates the diffusion equation, $M^{-1}(f(t,y))$, at a given time and a
+ // given $y$. <code>id_minus_tau_J_inverse()</code> 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
{
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<double> evaluate_diffusion(const double time, const Vector<double> &y) const;
+ Vector<double> evaluate_diffusion (const double time,
+ const Vector<double> &y) const;
- Vector<double> id_minus_tau_J_inverse(const double time,
- const double tau,
- const Vector<double> &y);
+ Vector<double> id_minus_tau_J_inverse (const double time,
+ const double tau,
+ const Vector<double> &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;
:
fe_degree(2),
diffusion_coefficient(1./30.),
- absorption_xs(1.),
+ absorption_cross_section(1.),
fe(fe_degree),
dof_handler(triangulation)
{}
- // @sect5{<code>Diffusion::setup_system</code>}
+ // @sect4{<code>Diffusion::setup_system</code>}
// 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);
- // @sect5{<code>Diffusion::assemble_system</code>}
- // 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{<code>Diffusion::assemble_system</code>} 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 <code>inverse_mass_matrix</code> variable will then store the inverse
+ // of the mass matrix so that $M^{-1}$ can be applied to a vector using the
+ // <code>vmult()</code> 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.;
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
- cell_matrix(i,j) += ((-diffusion_coefficient * fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) - absorption_xs *
- fe_values.shape_value(i,q_point) * fe_values.shape_value(j,q_point)) *
+ cell_matrix(i,j) += ((-diffusion_coefficient *
+ fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point)
+ - absorption_cross_section *
+ fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point)) *
fe_values.JxW(q_point));
- cell_mass_matrix(i,j) += fe_values.shape_value(i,q_point) *
- fe_values.shape_value(j,q_point) *
- fe_values.JxW(q_point);
+ cell_mass_matrix(i,j) += fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point) *
+ fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
- // @sect5{<code>Diffusion::get_source</code>}
+ // @sect4{<code>Diffusion::get_source</code>}
//
- // 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{<code>Diffusion:evaluate_diffusion</code>}
+ // @sect4{<code>Diffusion:evaluate_diffusion</code>}
//
- // Now, the weak form of the diffusion equation is evaluated at a given
- // time $t$ and for a given vector $y$.
- Vector<double> Diffusion::evaluate_diffusion(const double time, const Vector<double> &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
+ // <code>system_matrix</code>) 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<double> Diffusion::evaluate_diffusion (const double time,
+ const Vector<double> &y) const
{
Vector<double> tmp(dof_handler.n_dofs());
tmp = 0.;
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);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
- double source = get_source(time,fe_values.quadrature_point(q_point)) ;
+ const double source = get_source(time,
+ fe_values.quadrature_point(q_point));
for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_source(i) += source * fe_values.shape_value(i,q_point) *
+ cell_source(i) += source *
+ fe_values.shape_value(i,q_point) *
fe_values.JxW(q_point);
}
cell->get_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<double> value(dof_handler.n_dofs());
inverse_mass_matrix.vmult(value,tmp);
}
- // @sect5{<code>Diffusion::id_minus_tau_J_inverse</code>}
+ // @sect4{<code>Diffusion::id_minus_tau_J_inverse</code>}
//
// 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<double> Diffusion::id_minus_tau_J_inverse(const double time, const double tau,
- const Vector<double> &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<double> Diffusion::id_minus_tau_J_inverse (const double time,
+ const double tau,
+ const Vector<double> &y)
{
- Vector<double> tmp(dof_handler.n_dofs());
- Vector<double> 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<double> tmp(dof_handler.n_dofs());
+ mass_matrix.vmult(tmp, y);
+ Vector<double> result(y);
inverse_mass_minus_tau_Jacobian.vmult(result,tmp);
return result;
- // @sect5{<code>Diffusion::output_results</code>}
+ // @sect4{<code>Diffusion::output_results</code>}
//
- // 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;
}
}
- // @sect5{<code>Diffusion::explicit_method</code>}
- // This function is the driver for all the explicit method. It calls
- // <code>evolve_one_time_step</code> which performs one time step.
- // <code>evolve_one_time_step</code> needs to evaluate $M^{-1}(f(t,y))$,
- // i.e, it needs <code>evaluate_diffusion</code>. Because
- // <code>evaluate_diffusion</code> 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{<code>Diffusion::explicit_method</code>}
+ //
+ // This function is the driver for all the explicit methods. It calls
+ // <code>evolve_one_time_step</code> which performs one time step. For
+ // explicit methods, <code>evolve_one_time_step</code> needs to evaluate
+ // $M^{-1}(f(t,y))$, i.e, it needs <code>evaluate_diffusion</code>. Because
+ // <code>evaluate_diffusion</code> is a member function, it needs to be bound
+ // to <code>this</code>. 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<double> (n_time_steps);
double time = initial_time;
for (unsigned int i=0; i<n_time_steps; ++i)
{
time = explicit_runge_kutta.evolve_one_time_step(
- std_cxx11::bind(&Diffusion::evaluate_diffusion,this,std_cxx11::_1,std_cxx11::_2),
+ std_cxx11::bind(&Diffusion::evaluate_diffusion,
+ this,
+ std_cxx11::_1,
+ std_cxx11::_2),
time,time_step,solution);
if ((i+1)%10==0)
- // @sect5{<code>Diffusion::implicit_method</code>}
+ // @sect4{<code>Diffusion::implicit_method</code>}
// This function is equivalent to <code>explicit_method</code> 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<double> (n_time_steps);
double time = initial_time;
for (unsigned int i=0; i<n_time_steps; ++i)
{
time = implicit_runge_kutta.evolve_one_time_step(
- std_cxx11::bind(&Diffusion::evaluate_diffusion,this,std_cxx11::_1,std_cxx11::_2),
- std_cxx11::bind(&Diffusion::id_minus_tau_J_inverse,this,std_cxx11::_1,std_cxx11::_2,
+ std_cxx11::bind(&Diffusion::evaluate_diffusion,
+ this,
+ std_cxx11::_1,
+ std_cxx11::_2),
+ std_cxx11::bind(&Diffusion::id_minus_tau_J_inverse,
+ this,
+ std_cxx11::_1,
+ std_cxx11::_2,
std_cxx11::_3),
time,time_step,solution);
- // @sect5{<code>Diffusion::embedded_explicit_method</code>}
+ // @sect4{<code>Diffusion::embedded_explicit_method</code>}
// 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
// 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.
// </code>delta_t_guess</code> 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<double> (n_time_steps);
double time = initial_time;
const double coarsen_tol = 1e-5;
solution = 0.;
- TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> > embedded_explicit_runge_kutta(method,
- coarsen_param,refine_param,min_delta,max_delta,refine_tol,coarsen_tol);
- output_results(0,method);
+ TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> >
+ 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 (time<final_time)
{
- // The last time step is chosend such that the final time is exactly
- // reached.
if (time+time_step>final_time)
time_step = final_time-time;
- // @sect5{<code>Diffusion::run</code>}
- void Diffusion::run()
+ // @sect4{<code>Diffusion::run</code>}
+ //
+ // 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();
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,
projects / dealii.git / commitdiff