@f}
augmented by appropriate boundary conditions. Here, $v$ is the velocity of
neutrons (for simplicity we assume it is equal to 1), $D$ is the diffusion coefficient,
-$\Sigma_a$ is the absorption cross section, and $S$ is a source. Because we are only i
-nterested in the time dependence, we assume that $D$ and $\Sigma_a$ are constant. We are looking
-for a solution on a square domain $[0,b]\times[0,b]$ of the form:
+$\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. We are looking for a solution on a square domain $[0,b]\times[0,b]$ of
+the form:
@f{eqnarray*}
\phi(x,t) = A\sin(\omega t)(bx-x^2).
@f}
</ol>
<h4>Explicit Runge-Kutta</h4>
-These methods that include for forward Euler, third order Runge-Kutta, and
+These methods that include 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.
Only embedded explicit methods have been implemented at the time of the writing.
<h4>Implicit Runge-Kutta</h4>
-These methods include backward Euler, implicit midpoint, Crank-Nicolson, and the
+These methods include backward Euler, implicit midpoint, Crank-Nicolson, and a
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)$ or equivalently
+$\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$. These methods are
always stable.
To use the Runge-Kutta methods, we need to be able to evaluate:
@f{eqnarray*}
f = \oint D b_i \frac{\partial b_j}{\partial n} d\boldsymbol{r} - \int D \nabla
-b_i \nabla b_j \phi_j d\boldsymbol{r} -\int \Sigma_a b_i b_j \phi_j \phi_j
+b_i \nabla b_j \phi_j d\boldsymbol{r} -\int \Sigma_a b_i b_j \phi_j
d\boldsymbol{r} + \int b_j S d\boldsymbol{r}
@f}
and $\frac{\partial f}{\partial y}$. Because $f$ is linear in $y$ (or $\phi$ in
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, it is important to note to remember to:
+mesh must be adapted, it is important 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.
const double initial_time,
const double final_time);
- // Driver for the embedded explicit methods. Returns the number of steps
+ // Driver for the embedded explicit methods. This function returns the number of steps
// executed.
unsigned int embedded_explicit_method(TimeStepping::runge_kutta_method method,
const unsigned int n_time_steps,
- // We choose quadratic finite elements so there are no spatial error and we
- // initialize the parameters.
+ // We choose quadratic finite elements and we initialize the parameters.
Diffusion::Diffusion()
:
fe_degree(2),
cell = dof_handler.begin_active(),
endc = dof_handler.end();
- // Compute $-\int D \nabla b \cdot \nabla b - \int \Sigma_a b b $ and the mass matrix
- // $\int b b$.
+ // 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}$.
for (; cell!=endc; ++cell)
{
cell_matrix = 0.;
{
Vector<double> tmp(dof_handler.n_dofs());
tmp = 0.;
- // Compute $tmp=system_matrix y$.
+ // Compute $tmp=system\_matrix\cdot y$.
system_matrix.vmult(tmp,y);
const QGauss<2> quadrature_formula(fe_degree+1);
// 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 = \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My$.
inverse_mass_minus_tau_Jacobian.vmult(result,tmp);
return result;
projects / dealii-svn.git / commitdiff