<h3>Problem statement</h3>
-In this example, we solve the energy-integrated time-dependent diffusion
+In this example, we solve the one-group 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:
- \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 <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. The domain is square
-$[0,b]\times[0,b]$ and we are looking for a solution of the form:
+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:
@f{eqnarray*}
\phi(x,t) = A\sin(\omega t)(bx-x^2).
@f}
-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
+By using quadratic finite elements, all the error will be due to the time discretization. We
+impose the following boundary conditions: homogeneous Dirichlet for $x=0$ and
$x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. The source is
given by:
@f{eqnarray*}
<h4>Implicit Runge-Kutta</h4>
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)$ or equivalently
-$\left(M - \Deltat \frac{\partial f}{\partial y}\right)^{-1} M$. These methods are
+$\left(I-\Delta t M^{-1} \frac{\partial f}{\partial y}\right)$ or equivalently
+$\left(M - \Delta t \frac{\partial f}{\partial y}\right)^{-1} M$. These methods are
always stable.
+<h3>Weak form</h3>
+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
+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
+this case) $\frac{\partial f}{\partial y} y = f$.
+
<h3>Remarks</h3>
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
public:
Diffusion();
- // This function is the driver that will call the other ones.
+ // This function is the driver that will call the other functions.
void run();
private:
// the time.
void assemble_system();
- // Compute the intensity of the source at the given point.
+ // Compute the intensity of the source at a given time for a given point.
double get_source(double time,const Point<2> &point) const;
// Evaluate the diffusion equation $M^{-1}(f(t,y))$ at a given time and
{
Vector<double> tmp(dof_handler.n_dofs());
tmp = 0.;
- // Compute system_matrix*y
+ // Compute $tmp=system_matrix y$.
system_matrix.vmult(tmp,y);
const QGauss<2> quadrature_formula(fe_degree+1);
// @sect5{<code>Diffusion::output_results}
+ //
+ // We output the solution in vtu files.
void Diffusion::output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const
{
std::string method_name;
// @sect5{<code>Diffusion::run</code>}
void Diffusion::run()
{
- // Create the grid (a square [0,5]x[0,5]) and refine the mesh four times.
- // The final gird has 16 by 16 cells, for a total of 256.
+ // 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);
projects / dealii-svn.git / commitdiff