augmented by appropriate boundary conditions. Here, $v$ is the velocity of
neutrons, $D$ is the diffusion coefficient, $\Sigma_a$ is the <i>absorption
cross section</i>, 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. The domain is square
+time dependence, we assume that $D$ and $\Sigma_a$ are constant. 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 and all
+By using quadratic finite elements, there 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
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,\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.
+Because the solution is a sine, we know that $\phi\left(x,\pi\right) = 0$.
+Therefore, the error at this time is simply the norm of the numerical solution.
<h3>Runge-Kutta</h3>
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:
+mesh must be adapted, it is important to note to remember to:
<ol>
<li> 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.
// 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<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);
// Output the results as vtu files.
void output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const;
// @sect5{<code>Diffusion:evaluate_diffusion</code>}
//
- // Evaluate the diffusion weak form give a time t and a vector y.
+ // Evaluate the weak form of the diffusion equation at a given time t and for a given vector y.
Vector<double> Diffusion::evaluate_diffusion(const double time, const Vector<double> &y) const
{
Vector<double> tmp(dof_handler.n_dofs());
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 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$
+ // Compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp$.
inverse_mass_minus_tau_Jacobian.vmult(result,tmp);
return result;
}
- // sect5{<code>Diffusion::explicit_method</code>}
+ // @sect5{<code>Diffusion::explicit_method</code>}
void Diffusion::explicit_method(TimeStepping::runge_kutta_method method,
const unsigned int n_time_steps,
const double initial_time,
- // sect5{<code>Diffusion::implicit_method</code>}
+ // @sect5{<code>Diffusion::implicit_method</code>}
void Diffusion::implicit_method(TimeStepping::runge_kutta_method method,
const unsigned int n_time_steps,
const double initial_time,
- // sect5{<code>Diffusion::embedded_explicit_method</code>}
+ // @sect5{<code>Diffusion::embedded_explicit_method</code>}
unsigned int Diffusion::embedded_explicit_method(TimeStepping::runge_kutta_method method,
const unsigned int n_time_steps,
const double initial_time,
unsigned int n_steps=0;
while (time<final_time)
{
- // Choose the last time step to reach final_time.
+ // Choose the last time step to exactly reach the final time.
if (time+time_step>final_time)
time_step = final_time-time;
if ((n_steps+1)%10==0)
output_results(n_steps+1,method);
- // Update the time step
+ // Update the time step.
time_step = embedded_explicit_runge_kutta.get_status().delta_t_guess;
++n_steps;
}
- // sect5{<code>Diffusion::run</code>}
+ // @sect5{<code>Diffusion::run</code>}
void Diffusion::run()
{
// Create the grid (a square [0,5]x[0,5]) and refine the mesh four times.
// Use implicit midpoint.
implicit_method(TimeStepping::IMPLICIT_MIDPOINT,n_time_steps,initial_time,final_time);
std::cout<<"Implicit Midpoint error: "<<solution.l2_norm()<<std::endl;
- // Use Crank-NICOLSON.
+ // Use Crank-Nicolson.
implicit_method(TimeStepping::CRANK_NICOLSON,n_time_steps,initial_time,final_time);
std::cout<<"Crank-Nicolson error: "<<solution.l2_norm()<<std::endl;
// Use two stages SDIRK.
projects / dealii-svn.git / commitdiff