--- /dev/null
+\documentclass{article}
+\usepackage{amssymb,amsmath}
+\makeatletter
+\newcommand{\rmnum}[1]{\romannumeral #1}
+\newcommand{\Rmnum}[1]{\expandafter\@slowromancap\romannumeral #1@}
+\makeatother
+\begin{document}
+
+What is new in this example:
+\begin{enumerate}
+\item Solve multigroup neutron diffusion problem with multiple different meshes
+\item Solve an eigenvalue problem
+\item Setting up complicated material properties for nuclear fuel assemblies
+\end{enumerate}
+
+\subsection{Introduction}
+
+In this example, we intend to solve the multigroup diffusion approximation of
+the neutron transport equation. Essentially, the way to view this is as follows: In a
+nuclear reactor, neutrons are speeding around at different energies, get
+absorbed or scattered, or start a new fission
+event. If viewed at long enough length scales, the movement of neutrons can be
+considered a diffusion process.
+
+A mathematical description of this would group neutrons into energy bins, and
+consider the balance equations for the neutron fluxes in each of these
+bins, or energy groups. The scattering, absorption, and fission events would
+then be operators within the diffusion equation describing the neutron
+fluxes. Assume we have energy groups $g=1,\ldots,G$, where by convention we
+assume that the neutrons with the highest energy are in group 1 and those with
+the lowest energy in group $G$. Then the neutron flux of each group satisfies the
+following equations:
+\begin{eqnarray*}
+\frac 1{v_g}\frac{\partial \phi_g(x,t)}{\partial t}
+&=&
+\nabla \cdot(D_g(x) \nabla \phi_g(x,t))
+-
+\Sigma_{r,g}(x)\phi_g(x,t)
+\\
+&& \qquad
++
+\chi_g\sum_{g'=1}^G\nu\Sigma_{f,g'}(x)\phi_{g'}(x,t)
++
+\sum_{g'\ne g}\Sigma_{s,g'\to g}(x)\phi_{g'}(x,t)
++
+s_{\mathrm{ext},g}(x,t)
+\end{eqnarray*}
+augmented by appropriate boundary conditions. Here, $v_g$ is the velocity of
+neutrons within group $g$. In other words, the change in
+time in flux of neutrons in group $g$ is governed by the following
+processes:
+\begin{itemize}
+\item Diffusion $\nabla \cdot(D_g(x) \nabla \phi_g(x,t))$. Here, $D_g$ is the
+ (spatially variable) diffusion coefficient.
+\item Absorption $\Sigma_{r,g}(x)\phi_g(x,t)$ (note the
+ negative sign). The coefficient $\Sigma_{r,g}$ is called the \textit{removal
+ cross section}.
+\item Nuclear fission $\chi_g\sum_{g'=1}^G\nu\Sigma_{f,g'}(x)\phi_{g'}(x,t)$.
+ The production of neutrons of energy $g$ is
+ proportional to the flux of neutrons of energy $g'$ times the
+ probability $\Sigma_{f,g'}$ that neutrons of energy $g'$ cause a fission
+ event times the number $\nu$ of neutrons produced in each fission event
+ times the probability that a neutron produced in this event has energy
+ $g$. $\nu\Sigma_{f,g'}$ is called the \textit{fission cross section} and
+ $\chi_g$ the \textit{fission spectrum}. We will denote the term
+ $\chi_g\nu\Sigma_{f,g'}$ as the \textit{fission distribution cross
+ section} in the program.
+\item Scattering $\sum_{g'\ne g}\Sigma_{s,g'\to g}(x)\phi_{g'}(x,t)$
+ of neutrons of energy $g'$ producing neutrons
+ of energy $g$. $\Sigma_{s,g'\to g}$ is called the \textit{scattering cross
+ section}. The case of elastic, in-group scattering $g'=g$ exists, too, but
+ we subsume this into the removal cross section. The case $g'<g$ is called
+ down-scattering, since a neutron loses energy in such an event. On the
+ other hand, $g'>g$ corresponds to up-scattering: a neutron gains energy in
+ a scattering event from the thermal motion of the atoms surrounding it;
+ up-scattering is therefore only an important process for neutrons with
+ kinetic energies that are already on the same order as the thermal kinetic
+ energy (i.e. in the sub $eV$ range).
+\item An extraneous source $s_{\mathrm{ext},g}$.
+\end{itemize}
+
+For realistic simulations in reactor analysis, one may want to split the
+continuous spectrum of neutron energies into many energy groups, often up to
+100. However, if neutron energy spectra are known well enough for some type of
+reactor (for example Pressurized Water Reactors, PWR), it is possible to obtain
+satisfactory results with only 2 energy groups.
+
+In the program shown in this tutorial program, we provide the structure to
+compute with as many energy groups as desired. However, to keep computing
+times moderate and in order to avoid tabulating hundreds of coefficients, we
+only provide the coefficients for above equations for a two-group simulation,
+i.e. $g=1,2$. We do, however, consider a realistic situation by assuming that
+the coefficients are not constant, but rather depend on the materials that are
+assembled into reactor fuel assemblies in rather complicated ways (see
+below).
+
+
+\subsection{The eigenvalue problem}
+
+If we consider all energy groups at once, we may write above equations in the
+following operator form:
+\begin{equation}
+\frac 1v \frac{\partial \phi}{\partial t}
+=
+-L\phi
++
+F\phi
++
+X\phi
++
+s_{\mathrm{ext}},
+\end{equation}
+where $L,F,X$ are sinking, fission, and scattering operators,
+respectively. $L$ here includes both the diffusion and removal terms. Note
+that $L$ is symmetric, whereas $F$ and $X$ are not.
+
+It is well known that this equation admits a stable solution if all
+eigenvalues of the operator $-L+F+X$ are negative. This can be readily seen by
+multiplying the equation by $\phi$ and integrating over the domain, leading to
+\begin{equation}
+ \frac 1{2v} \frac{\partial}{\partial t} \|\phi\|^2 = ((-L+F+X)\phi,\phi).
+\end{equation}
+Stability means that the solution does not grow, i.e. we want the left hand
+side to be less than zero, which is the case if the eigenvalues of the
+operator on the right are all negative. For obvious reasons, it is
+not very desirable if a nuclear reactor produces neutron fluxes that grow
+exponentially, so eigenvalue analyses are the bread-and-butter of nuclear
+engineers. The main point of the program is therefore to consider the
+eigenvalue problem
+\begin{equation}
+ (L-F-X) \phi = \lambda \phi,
+\end{equation}
+where we want to make sure that all eigenvalues are positive. Note that $L$,
+being the diffusion operator plus the absorption (removal), is positive
+definite; the condition that all eigenvalues are positive therefore means that
+we want to make sure that fission and inter-group scattering are weak enough
+to not shift the spectrum into the negative.
+
+In nuclear engineering, one typically looks at a slightly different
+formulation of the eigenvalue problem. To this end, we do not just multiply
+with $\phi$ and integrate, but rather multiply with $\phi(L-X)^{-1}$. We then
+get the following evolution equation:
+\begin{equation}
+ \frac 1{2v} \frac{\partial}{\partial t} \|\phi\|^2_{(L-X)^{-1}} = ((L-X)^{-1}(-L+F+X)\phi,\phi).
+\end{equation}
+Stability is the guaranteed if the eigenvalues of the following problem are
+all negative:
+\begin{equation}
+ (L-X)^{-1}(-L+F+X)\phi = \lambda_F \phi,
+\end{equation}
+which is equivalent to the eigenvalue problem
+\begin{equation}
+ (L-X)\phi = \frac 1{\lambda_F+1} F \phi.
+\end{equation}
+The typical formulation in nuclear engineering is to write this as
+\begin{equation}
+ (L-X) \phi = \frac 1{k_{\mathrm{eff}}} F \phi,
+\end{equation}
+where $k_{\mathrm{eff}}=\frac 1{\lambda^F+1}$.
+Intuitively, $k_{\mathrm{eff}}$ is something like the multiplication
+factor for neutrons per typical time scale and should be less than or equal to
+one for stable operation of a reactor: if it is less than one, the chain reaction will
+die down, whereas nuclear bombs for example have a $k$-eigenvalue larger than
+one. A stable reactor should have $k_{\mathrm{eff}}=1$.
+
+[For those who wonder how this can be achieved in practice without
+inadvertently getting slightly larger than one and triggering a nuclear bomb:
+first, fission processes happen on different time scales. While most neutrons
+are releases very quickly after a fission event, a small number of neutrons
+are only released by daughter nuclei after several further decays, up to 10-60
+seconds after the fission was initiated. If one is therefore slightly beyond
+$k_{\mathrm{eff}}=1$, one therefore has many seconds to react until all the
+neutrons created in fission re-enter the fission cycle. Nevertheless, control
+rods in nuclear reactors absorbing neutrons -- and therefore reducing
+$k_{\mathrm{eff}}$ -- are designed in such a way that they are all the way in
+the reactor in at most 2 seconds.
+
+One therefore has on the order of 10-60 seconds to regulate the nuclear reaction
+if $k_{\mathrm{eff}}$ should be larger than one for some time, as indicated by
+a growing neutron flux. Regulation can be achieved by continuously monitoring
+the neutron flux, and if necessary increase or reduce neutron flux by moving
+neutron-absorbing control rods a few millimeters into or out of the
+reactor. On a longer scale, the water cooling the reactor contains boron, a
+good neutron absorber. Every few hours, boron concentrations are adjusted by
+adding boron or diluting the coolant.
+
+Finally, some of the absorption and scattering reactions have some
+stability built in; for example, higher neutron fluxes result in locally
+higher temperatures, which lowers the density of water and therefore reduces
+the number of scatterers that are necessary to moderate neutrons from high to
+low energies before they can start fission events themselves.]
+
+In this tutorial program, we solve above $k$-eigenvalue problem for two energy
+groups, and we are looking for the largest multiplication factor
+$k_{\mathrm{eff}}$, which is proportional to the inverse of the minimum
+eigenvalue plus one. To solve the eigenvalue problem, we generally
+use a modified version of the \emph{inverse power method}. The algorithm looks
+like this:
+
+\begin{enumerate}
+\item Initialize $\phi_g$ and $k_{\mathrm{eff}}$ with $\phi_g^{(0)}$
+ and $k_{\mathrm{eff}}^{(0)}$ and let $n=1$.
+
+\item Define the so-called \textit{fission source} by
+ \begin{equation}
+ s_f^{(n-1)}(x)
+ =
+ \frac{1}{k_{\mathrm{eff}}^{(n-1)}}
+ \sum_{g'=1}^G\nu\Sigma_{f,g'}(x)\phi_{g'}^{(n-1)}(x).
+ \end{equation}
+
+\item Solve for all group fluxes $\phi_g,g=1,\ldots,G$ using
+ \begin{equation}
+ -\nabla \cdot D_g\nabla \phi_g^{(n)}
+ +
+ \Sigma_{r,g}\phi_g^{(n)}
+ =
+ \chi_g s_f^{(n-1)}
+ +
+ \sum_{g'< g} \Sigma_{s,g'\to g} \phi_{g'}^{(n)}
+ +
+ \sum_{g'> g}\Sigma_{s,g'\to g}\phi_{g'}^{(n-1)}.
+ \end{equation}
+
+\item Update
+ \begin{equation}
+ k_{\mathrm{eff}}^{(n)}
+ =
+ \sum_{g'=1}^G
+ \int_{\Omega}\nu\Sigma_{f,g'}(x)
+ \phi_{g'}^{(n)}(x)dx.
+ \end{equation}
+
+\item Compare $k_{\mathrm{eff}}^{(n)}$ with $k_{\mathrm{eff}}^{(n-1)}$.
+ If the change greater than a prescribed tolerance then set $n=n+1$ repeat
+ the iteration starting at step 2, otherwise end the iteration.
+\end{enumerate}
+
+Note that in this scheme, we do not solve group fluxes exactly in each power
+iteration, but rather consider previously compute $\phi_{g'}^{(n)}$ only for
+down-scattering events $g'<g$. Up-scattering is only treated by using old
+iterators $\phi_{g'}^{(n-1)}$, in essence assuming that the scattering
+operator is triangular. This is physically motivated since up-scattering does
+not play a too important role in neutron scattering. In addition, practices
+shows that the inverse power iteration is stable even using this
+simplification.
+
+Note also that one can use lots of extrapolation techniques to accelerate the
+power iteration laid out above. However, none of these are implemented in this
+example.
+
+
+\subsection{Meshes and mesh refinement}
+
+One may wonder whether it is appropriate to solve for the solutions of the
+individual energy group equations on the same meshes. The question boils down
+to this: will $\phi_g$ and $\phi_{g'}$ have similar smoothness properties? If
+this is the case, then it is appropriate to use the same mesh for the two; a
+typical application could be chemical combustion, where typically the
+concentrations of all or most chemical species change rapidly within the flame
+front. As it turns out, and as will be apparent by looking at the
+graphs shown in the results section of this tutorial program, this isn't the
+case here, however: since the diffusion coefficient is different for different
+energy groups, fast neutrons (in bins with a small group number $g$) have a very
+smooth flux function, whereas slow neutrons (in bins with a large group
+number) are much more affected by the local material properties and have a
+correspondingly rough solution if the coefficient are rough as in the case we
+compute here. Consequently, we will want to use different meshes to compute
+each energy group.
+
+This has two implications that we will have to consider: First, we need to
+find a way to refine the meshes individually. Second, assembling the source
+terms for the inverse power iteration, where we have to integrate solution
+$\phi_{g'}^{(n)}$ defined on mesh $g'$ against the shape functions defined on
+mesh $g$, becomes a much more complicated task.
+
+
+\subsubsection{Mesh refinement}
+
+We use the usual paradigm: solve on a given mesh, then evaluate an error
+indicator for each cell of each mesh we have. Because it is so convenient, we
+again use the \emph{a posteriori} error estimator by Kelly, Gago, Zienkiewicz
+and Babuska which approximates the error per cell by integrating the jump of
+the gradient of the solution along the faces of each cell. Using this, we
+obtain indicators
+\begin{equation}
+\eta_{g,K}, \qquad g=1,2,\ldots,G,\qquad K\in{\mathbb T}_g,
+\end{equation}
+where ${\mathbb T}_g$ is the triangulation used in the solution of
+$\phi_g$. The question is what to do with this. For one, it is clear that
+refining only those cells with the highest error indicators might lead to bad
+results. To understand this, it is important to realize that $\eta_{g,K}$
+scales with the second derivative of $\phi_g$. In other words, if we have two
+energy groups $g=1,2$ whose solutions are equally smooth but where one is
+larger by a factor of 10,000, for example, then only the cells of that mesh
+will be refined, whereas the mesh for the solution of small magnitude will
+remain coarse. This is probably not what one wants, since we can consider both
+components of the solution equally important.
+
+In essence, we would therefore have to scale $\eta_{g,K}$ by an importance
+factor $z_g$ that says how important it is to resolve $\phi_g$ to any given
+accuracy. Such important factors can be computed using duality techniques
+(see, for example, the step-14 tutorial program, and the
+reference to the book by Bangerth and Rannacher cited there). We won't go
+there, however, and simply assume that all energy groups are equally
+important, and will therefore normalize the error indicators $\eta_{g,K}$ for
+group $g$ by the maximum of the solution $\phi_g$. We then refine the cells
+whose errors satisfy
+\begin{equation}
+ \frac{\eta_{g,K}}{\|\phi_g\|_\infty}
+ >
+ \alpha_1
+ \displaystyle{\max_{\substack{1\le g\le G\\K\in {\mathbb T}_g}}
+ \frac{\eta_{g,K}}{\|\phi_g\|_\infty}}
+\end{equation}
+and coarsen the cells where
+\begin{equation}
+ \frac{\eta_{g,K}}{\|\phi_g\|_\infty}
+ <
+ \alpha_2
+ \displaystyle{\max_{\substack{1\le g\le G\\K\in {\mathbb T}_g}}
+ \frac{\eta_{g,K}}{\|\phi_g\|_\infty}}.
+\end{equation}
+We chose $\alpha_1=0.3$ and $\alpha_2=0.01$ in the code. Note that this will,
+of course, lead to different meshes for the different energy groups.
+
+The strategy above essentially means the following: If for energy group $g$
+there are many cells $K\in {\mathbb T}_g$ on which the error is large, for
+example because the solution is globally very rough, then many cells will be
+above the threshold. On the other hand, if there are a few cells with large
+and many with small errors, for example because the solution is overall rather
+smooth except at a few places, then only the few cells with large errors will
+be refined. Consequently, the strategy allows for meshes that track the global
+smoothness properties of the corresponding solutions rather well.
+
+
+\subsubsection{Assembling terms on different meshes}
+
+As pointed out above, the multigroup refinement strategy results in
+different meshes for the different solutions $\phi_g$. So what's the problem?
+In essence it goes like this: in step 3 of the eigenvalue iteration, we have
+form the weak form for the equation to compute $\phi_g^{(n)}$ as usual by
+multiplication with test functions $\varphi_g^i$ defined on the mesh for
+energy group $g$; in the process, we have to
+compute the right hand side vector that contains terms of the following form:
+\begin{equation}
+ F_i = \int_\Omega f(x) \varphi_g^i(x) \phi_{g'}(x) \ dx,
+\end{equation}
+where $f(x)$ is one of the coefficient functions $\Sigma_{s,g'\to g}$ or
+$\nu\chi_g\Sigma_{f,g'}$ used in the right hand side
+of eigenvalue equation. The difficulty now is that $\phi_{g'}$ is defined on
+the mesh for energy group $g'$, i.e. it can be expanded as
+$\phi_{g'}(x)=\sum_j\phi_{g'}^j \varphi_{g'}^j(x)$, with basis functions
+$\varphi_{g'}^j(x)$ defined on mesh $g'$. The contribution to the right hand
+side can therefore be written as
+\begin{equation}
+ F_i = \sum_j \left\{\int_\Omega f(x) \varphi_g^i(x) \varphi_{g'}^j(x)
+ \ dx \right\} \phi_{g'}^j ,
+\end{equation}
+On the other hand, the test functions $\varphi_g^i(x)$ are defined on mesh
+$g$. This means that we can't just split the integral $\Omega$ into integrals
+over the cells of either mesh $g$ or $g'$, since the respectively other basis
+functions may not be defined on these cells.
+
+The solution to this problem lies in the fact that both the meshes for $g$ and
+$g'$ are derived by adaptive refinement from a common coarse mesh. We can
+therefore always find a set of cells, which we denote by ${\mathbb T}_g \cap
+{\mathbb T}_{g'}$, that satisfy the following conditions:
+\begin{itemize}
+\item the union of the cells covers the entire domain, and
+\item a cell $K \in {\mathbb T}_g \cap {\mathbb T}_{g'}$ is active on at least
+ one of the two meshes.
+\end{itemize}
+A way to construct this set is to take each cell of coarse mesh and do the
+following steps: (i) if the cell is active on either ${\mathbb T}_g$ or
+${\mathbb T}_{g'}$, then add this cell to the set; (ii) otherwise, i.e. if
+this cell has children on both meshes, then do step (i) for each of the
+children of this cell. In fact, deal.II has a function
+\texttt{GridTools::get\_finest\_common\_cells} that computes exactly this set
+of cells that are active on at least one of two meshes.
+
+With this, we can write above integral as follows:
+\begin{equation}
+ F_i
+ =
+ \sum_{K \in {\mathbb T}_g \cap {\mathbb T}_{g'}}
+ \sum_j \left\{\int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x)
+ \ dx \right\} \phi_{g'}^j.
+\end{equation}
+ In the code, we
+compute the right hand side in the function
+\texttt{NeutronDiffusionProblem::assemble\_rhs}, where (among other things) we
+loop over the set of common most refined cells, calling the function
+\texttt{NeutronDiffusionProblem::assemble\_common\_cell} on each pair of
+these cells.
+
+By construction, there are now three cases to be considered:
+\begin{itemize}
+\item[(i)] The cell $K$ is active on both meshes, i.e. both the basis
+ functions $\varphi_g^i$ as well as $\varphi_{g'}^j$ are defined on $K$.
+\item[(ii)] The cell $K$ is active on mesh $g$, but not $g'$, i.e. the
+ $\varphi_g^i$ are defined on $K$, whereas the $\varphi_{g'}^j$ are defined
+ on children of $K$.
+\item[(iii)] The cell $K$ is active on mesh $g'$, but not $g$, with opposite
+ conclusions than in (ii).
+\end{itemize}
+
+To compute the right hand side above, we then need to have different code for
+these three cases, as follows:
+\begin{itemize}
+\item[(i)] If the cell $K$ is active on both meshes, then we can directly
+ evaluate the integral. In fact, we don't even have to bother with the basis
+ functions $\varphi_{g'}$, since all we need is the values of $\phi_{g'}$ at
+ the quadrature points. We can do this using the
+ \texttt{FEValues::get\_function\_values} function. This is done directly in
+ the \texttt{NeutronDiffusionProblem::assemble\_common\_cell} function.
+
+\item[(ii)] If the cell $K$ is active on mesh $g$, but not $g'$, then the
+ basis functions $\varphi_{g'}^j$ are only defined either on the children
+ $K_c,0\le c<2^{\texttt{dim}}$, or on children of these children if cell $K$
+ is refined more than once more on mesh $g'$.
+
+ Let us assume for a second that $K$ is only once more refined on mesh $g'$
+ than on mesh $g$. Using the fact that we use embedded finite element spaces
+ where each basis function on one mesh can be written as a linear combination
+ of basis functions on the next refined mesh, we can expand the restriction
+ of $\phi_g^i$ to child cell $K_c$ into the basis functions defined on that
+ child cell (i.e. on cells on which the basis functions $\varphi_{g'}^l$ are
+ defined):
+ \begin{equation}
+ \phi_g^i|_{K_c} = B_c^{il} \varphi_{g'}^l|_{K_c}.
+ \end{equation}
+ Here, and in the following, summation over indices appearing twice is
+ implied. The matrix $B_c$ is the matrix that interpolated data from a cell
+ to its $c$-th child.
+
+ Then we can write the contribution of cell $K$ to the right hand side
+ component $F_i$ as
+ \begin{eqnarray*}
+ F_i|_K
+ &=&
+ \left\{ \int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x)
+ \ dx \right\} \phi_{g'}^j
+ \\
+ &=&
+ \left\{
+ \sum_{0\le c<2^{\texttt{dim}}}
+ B_c^{il} \int_{K_c} f(x) \varphi_{g'}^l(x) \varphi_{g'}^j(x)
+ \ dx \right\} \phi_{g'}^j.
+ \end{eqnarray*}
+ In matrix notation, this can be written as
+ \begin{eqnarray*}
+ F_i|_K
+ =
+ \sum_{0\le c<2^{\texttt{dim}}}
+ F_i|_{K_c},
+ \qquad
+ \qquad
+ F_i|_{K_c} = B_c^{il} M_{K_c}^{lj} \phi_{g'}^j
+ = (B_c M_{K_c})^{il} \phi_{g'}^j,
+ \end{eqnarray*}
+ where $M_{K_c}^{lj}=\int_{K_c} f(x) \varphi_{g'}^l(x) \varphi_{g'}^j(x)$ is
+ the weighted mass matrix on child $c$ of cell $K$.
+
+ The next question is what happens if a child $K_c$ of $K$ is not
+ active. Then, we have to apply the process recursively, i.e. we have to
+ interpolate the basis functions $\varphi_g^i$ onto child $K_c$ of $K$, then
+ onto child $K_{cc'}$ of that cell, onto child $K_{cc'c''}$ of that one, etc,
+ until we find an active cell. We then have to sum up all the contributions
+ from all the children, grandchildren, etc, of cell $K$, with contributions
+ of the form
+ \begin{equation}
+ F_i|_{K_{cc'}} = (B_cB_{c'} M_{K_{cc'}})^{ij} \phi_{g'}^j,
+ \end{equation}
+ or
+ \begin{equation}
+ F_i|_{K_{cc'c''}} = (B_c B_{c'} B_{c''}M_{K_{cc'c''}})^{ij}
+ \phi_{g'}^j,
+ \end{equation}
+ etc. We do this process recursively, i.e. if we sit on cell $K$ and see that
+ it has children on grid $g'$, then we call a function
+ \texttt{assemble\_case\_2} with an identity matrix; the function will
+ multiply it's argument from the left with the prolongation matrix; if the
+ cell has further children, it will call itself with this new matrix,
+ otherwise it will perform the integration.
+
+\item[(iii)] The last case is where $K$ is active on mesh $g'$ but not mesh
+ $g$. In that case, we have to express basis function $\varphi_{g'}^j$ in
+ terms of the basis functions defined on the children of cell $K$, rather
+ than $\varphi_g^i$ as before. This of course works in exactly the same
+ way. If the children of $K$ are active on mesh $g$, then
+ leading to the expression
+ \begin{eqnarray*}
+ F_i|_K
+ &=&
+ \left\{ \int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x)
+ \ dx \right\} \phi_{g'}^j
+ \\
+ &=&
+ \left\{
+ \sum_{0\le c<2^{\texttt{dim}}}
+ \int_{K_c} f(x) \varphi_{g'}^i(x) B_c^{jl} \varphi_{g'}^l(x)
+ \ dx \right\} \phi_{g'}^j.
+ \end{eqnarray*}
+ In matrix notation, this expression now reads as
+ \begin{eqnarray*}
+ F_i|_K
+ =
+ \sum_{0\le c<2^{\texttt{dim}}}
+ F_i|_{K_c},
+ \qquad
+ \qquad
+ F_i|_{K_c} = M_{K_c}^{il} B_c^{jl} \phi_{g'}^j
+ =
+ (M_{K_c} B_c^T)^{ij} \phi_{g'}^j,
+ \end{eqnarray*}
+ and correspondingly for cases where cell $K$ is refined more than once on
+ mesh $g$:
+ \begin{equation}
+ F_i|_{K_{cc'}} = (M_{K_{cc'}} B_{c'}^T B_c^T)^{ij} \phi_{g'}^j,
+ \end{equation}
+ or
+ \begin{equation}
+ F_i|_{K_{cc'c''}} = (M_{K_{cc'c''}} B_{c''}^T B_{c'}^T B_c^T)^{ij}
+ \phi_{g'}^j,
+ \end{equation}
+ etc. In other words, the process works in exactly the same way as before,
+ except that we have to take the transpose of the prolongation matrices and
+ need to multiply it to the mass matrix from the other side.
+\end{itemize}
+
+
+The expressions for cases (ii) and (iii) can be understood as repeatedly
+interpolating either the left or right basis functions in the scalar product
+$(f \varphi_g^i, \varphi_{g'}^j)_K$ onto child cells, and then finally
+forming the inner product (the mass matrix) on the final cell. To make the
+symmetry in these cases more obvious, we can write them like this: for case
+(ii), we have
+\begin{equation}
+ F_i|_{K_{cc'\cdots c^{(k)}}}
+ = [B_c B_{c'} \cdots B_{c^{(k)}} M_{K_{cc'\cdots c^{(k)}}}]^{ij}
+ \phi_{g'}^j,
+\end{equation}
+whereas for case (iii) we get
+\begin{equation}
+ F_i|_{K_{cc'\cdots c^{(k)}}}
+ = [(B_c B_{c'} \cdots B_{c^{(k)}} M_{K_{cc'\cdots c^{(k)}}})^T]^{ij}
+ \phi_{g'}^j,
+\end{equation}
+
+
+
+\subsection{Description of the test case}
+
+A nuclear reactor core is composed of different types of assemblies. An
+assembly is essentially the smallest unit that can be moved in and out of a
+reactor, and is usually rectangular or square. However, assemblies are not
+fixed units, as they are assembled from a complex lattice of different fuel
+rods, control rods, and instrumentation elements that are held in place
+relative to each other by spacers that are permanently attached to the rods.
+To make things more complicated, there are different kinds of assemblies that
+are used at the same time in a reactor, where assemblies differ in the type
+and arrangement of rods they are made up of.
+
+Obviously, the arrangement of assemblies as well as the arrangement of rods
+inside them affect the distribution of neutron fluxes in the reactor (a fact
+that will be obvious by looking at the solution shown below in the results
+sections of this program). Fuel rods, for example, differ from each other in
+the enrichment of U-235 or Pu-239. Control rods, on the other hand, have zero
+fission, but nonzero scattering and absorption cross sections.
+
+This whole arrangement would make the description or spatially dependent
+material parameters very complicated. It will not become much simpler, but we
+will make one approximation: we merge the volume inhabited by each cylindrical
+rod and the surrounding water into volumes of quadratic cross section into
+so-called ``pin cells'' for which homogenized material data are obtained with
+nuclear database and knowledge of neutron spectrum. The homogenization makes
+all material data piecewise constant on the solution domain for a reactor with
+fresh fuel. Spatially dependent material parameters are then looked up for the
+quadratic assembly in which a point is located, and then for the quadratic pin
+cell within this assembly.
+
+In this tutorial program, we simulate a quarter of a reactor consisting of $4
+\times 4$ assemblies. We use symmetry (Neumann) boundary conditions to reduce
+the problem to one quarter of the domain, and consequently only simulate a
+$2\times 2$ set of assemblies. Two of them will be UO${}_2$ fuel, the other
+two of them MOX fuel. Each of these assemblies consists of $17\times 17$ rods
+of different compositions. In total, we therefore create a $34\times 34$
+lattice of rods. To make things simpler later on, we reflect this fact by
+creating a coarse mesh of $34\times 34$ cells (even though the domain is a
+square, for which we would usually use a single cell). In deal.II, each cell
+has a \texttt{material\_id} which one may use to associated each cell with a
+particular number identifying the material from which this cell's volume is
+made of; we will use this material ID to identify which of the 8 different
+kinds of rods that are used in this testcase make up a particular cell. Note
+that upon mesh refinement, the children of a cell inherit the material ID,
+making it simple to track the material even after mesh refinement.
+
+The arrangement of the rods will be clearly visible in the images shown in
+the results section. The cross sections for materials and for both energy
+groups are taken from a OECD/NEA benchmark problem. The detailed configuration
+and material data is given in the code.
+
+
+\subsection{What the program does (and how it does that)}
+
+As a coarse overview of what exactly the program does, here is the basic
+layout: starting on a coarse mesh that is the same for each energy group, we
+compute inverse eigenvalue iterations to compute the $k$-eigenvalue on a given
+set of meshes. We stop these iterations when the change in the eigenvalue
+drops below a certain tolerance, and then write out the meshes and solutions
+for each energy group for inspection by a graphics program. Because the meshes
+for the solutions are different, we have to generate a separate output file
+for each energy group, rather than being able to add all energy group
+solutions into the same file.
+
+After this, we evaluate the error indicators as explained in one of the sections
+above for each of the meshes, and refine and coarsen the cells of each mesh
+independently. Since the eigenvalue iterations are fairly expensive, we don't
+want to start all over on the new mesh; rather, we use the SolutionTransfer
+class to interpolate the solution on the previous mesh to the next one upon
+mesh refinement. A simple experiment will convince you that this is a lot
+cheaper than if we omitted this step. After doing so, we resume our eigenvalue
+iterations on the next set of meshes.
+
+The program is controlled by a parameter file, using the ParameterHandler
+class already mentioned in the step-19 example program. We will show a
+parameter file in the results section of this section. For the moment suffice
+it to say that it controls the polynomial degree of the finite elements used,
+the number of energy groups (even though all that is presently implemented are
+the coefficients for a 2-group problem), the tolerance where to stop the
+inverse eigenvalue iteration, and the number of refinement cycles we will do.
+
+\end{document}
--- /dev/null
+<a name="Results"></a> <h1>Results</h1>
+
+
+The output of this program consist of the console output, a file
+named ``convergence_table'' to record main results of mesh iteration, the eps
+files including the grids, and the solutions given in gnuplot format.
+
+When we set Polynomial_Order to 2, we got following console output:
+@code
+Cycle 0:
+ Numbers of active cells: 1156 1156
+ Numbers of degrees of freedom: 4761 4761
+Iter number:1 k_eff=319.375676634307 flux ratio=6.836246075631 max_thermal=1.433899030144
+Iter number:2 k_eff=0.834072546055 flux ratio=5.204601882141 max_thermal=0.004630925876
+Iter number:3 k_eff=0.862826188043 flux ratio=4.645051765984 max_thermal=0.005380396338
+Iter number:4 k_eff=0.877887920967 flux ratio=4.318030683875 max_thermal=0.006005512201
+Iter number:5 k_eff=0.887161559547 flux ratio=4.256596788174 max_thermal=0.006639443035
+Iter number:6 k_eff=0.893254525197 flux ratio=4.296498905676 max_thermal=0.007161016401
+Iter number:7 k_eff=0.897386466621 flux ratio=4.323736110066 max_thermal=0.007541125053
+Iter number:8 k_eff=0.900235644733 flux ratio=4.342491852394 max_thermal=0.007813654241
+Iter number:9 k_eff=0.902217719823 flux ratio=4.355367629620 max_thermal=0.008007335384
+Iter number:10 k_eff=0.903602785157 flux ratio=4.364212965582 max_thermal=0.008144201718
+Iter number:11 k_eff=0.904572678811 flux ratio=4.370302672219 max_thermal=0.008240563668
+Iter number:12 k_eff=0.905252379018 flux ratio=4.374506018233 max_thermal=0.008308245052
+Iter number:13 k_eff=0.905728767660 flux ratio=4.377414535866 max_thermal=0.008355707257
+Iter number:14 k_eff=0.906062594755 flux ratio=4.379431495993 max_thermal=0.008388956163
+Iter number:15 k_eff=0.906296449777 flux ratio=4.380832749068 max_thermal=0.008412232394
+Iter number:16 k_eff=0.906460217413 flux ratio=4.381807689696 max_thermal=0.008428519956
+Iter number:17 k_eff=0.906574868174 flux ratio=4.382486812297 max_thermal=0.008439913928
+Iter number:18 k_eff=0.906655112098 flux ratio=4.382960306180 max_thermal=0.008447883119
+Iter number:19 k_eff=0.906711262646 flux ratio=4.383290664907 max_thermal=0.008453456292
+Iter number:20 k_eff=0.906750547146 flux ratio=4.383521281424 max_thermal=0.008457353551
+Iter number:21 k_eff=0.906778027984 flux ratio=4.383682335660 max_thermal=0.008460078740
+Iter number:22 k_eff=0.906797249757 flux ratio=4.383794844804 max_thermal=0.008461984300
+Iter number:23 k_eff=0.906810693563 flux ratio=4.383873459628 max_thermal=0.008463316724
+Iter number:24 k_eff=0.906820095658 flux ratio=4.383928400683 max_thermal=0.008464248386
+Iter number:25 k_eff=0.906826670830 flux ratio=4.383966802041 max_thermal=0.008464899825
+Iter number:26 k_eff=0.906831268888 flux ratio=4.383993645532 max_thermal=0.008465355326
+Iter number:27 k_eff=0.906834484255 flux ratio=4.384012411177 max_thermal=0.008465673822
+Iter number:28 k_eff=0.906836732678 flux ratio=4.384025530521 max_thermal=0.008465896521
+Iter number:29 k_eff=0.906838304919 flux ratio=4.384034702833 max_thermal=0.008466052239
+Iter number:30 k_eff=0.906839404318 flux ratio=4.384041115801 max_thermal=0.008466161120
+Iter number:31 k_eff=0.906840173074 flux ratio=4.384045599636 max_thermal=0.008466237253
+Iter number:32 k_eff=0.906840710623 flux ratio=4.384048734710 max_thermal=0.008466290487
+Iter number:33 k_eff=0.906841086501 flux ratio=4.384050926767 max_thermal=0.008466327710
+Iter number:34 k_eff=0.906841349329 flux ratio=4.384052459477 max_thermal=0.008466353737
+Iter number:35 k_eff=0.906841533109 flux ratio=4.384053531173 max_thermal=0.008466371936
+Iter number:36 k_eff=0.906841661615 flux ratio=4.384054280525 max_thermal=0.008466384662
+Iter number:37 k_eff=0.906841751471 flux ratio=4.384054804489 max_thermal=0.008466393560
+Iter number:38 k_eff=0.906841814301 flux ratio=4.384055170858 max_thermal=0.008466399781
+Iter number:39 k_eff=0.906841858234 flux ratio=4.384055427034 max_thermal=0.008466404132
+Iter number:40 k_eff=0.906841888954 flux ratio=4.384055606159 max_thermal=0.008466407174
+Iter number:41 k_eff=0.906841910434 flux ratio=4.384055731409 max_thermal=0.008466409301
+Iter number:42 k_eff=0.906841925454 flux ratio=4.384055818987 max_thermal=0.008466410788
+Iter number:43 k_eff=0.906841935956 flux ratio=4.384055880225 max_thermal=0.008466411828
+Iter number:44 k_eff=0.906841943300 flux ratio=4.384055923044 max_thermal=0.008466412555
+Iter number:45 k_eff=0.906841948435 flux ratio=4.384055952984 max_thermal=0.008466413064
+Iter number:46 k_eff=0.906841952025 flux ratio=4.384055973920 max_thermal=0.008466413419
+Iter number:47 k_eff=0.906841954536 flux ratio=4.384055988559 max_thermal=0.008466413668
+Iter number:48 k_eff=0.906841956291 flux ratio=4.384055998794 max_thermal=0.008466413842
+Iter number:49 k_eff=0.906841957518 flux ratio=4.384056005952 max_thermal=0.008466413963
+Iter number:50 k_eff=0.906841958377 flux ratio=4.384056010956 max_thermal=0.008466414048
+Iter number:51 k_eff=0.906841958977 flux ratio=4.384056014456 max_thermal=0.008466414108
+Iter number:52 k_eff=0.906841959397 flux ratio=4.384056016902 max_thermal=0.008466414149
+Iter number:53 k_eff=0.906841959690 flux ratio=4.384056018613 max_thermal=0.008466414178
+Iter number:54 k_eff=0.906841959895 flux ratio=4.384056019810 max_thermal=0.008466414199
+Iter number:55 k_eff=0.906841960039 flux ratio=4.384056020646 max_thermal=0.008466414213
+Iter number:56 k_eff=0.906841960139 flux ratio=4.384056021231 max_thermal=0.008466414223
+Iter number:57 k_eff=0.906841960209 flux ratio=4.384056021640 max_thermal=0.008466414230
+Iter number:58 k_eff=0.906841960258 flux ratio=4.384056021926 max_thermal=0.008466414235
+Iter number:59 k_eff=0.906841960292 flux ratio=4.384056022126 max_thermal=0.008466414238
+Iter number:60 k_eff=0.906841960316 flux ratio=4.384056022266 max_thermal=0.008466414240
+Iter number:61 k_eff=0.906841960333 flux ratio=4.384056022364 max_thermal=0.008466414242
+Iter number:62 k_eff=0.906841960345 flux ratio=4.384056022432 max_thermal=0.008466414243
+Iter number:63 k_eff=0.906841960353 flux ratio=4.384056022480 max_thermal=0.008466414244
+Iter number:64 k_eff=0.906841960359 flux ratio=4.384056022513 max_thermal=0.008466414245
+Iter number:65 k_eff=0.906841960363 flux ratio=4.384056022537 max_thermal=0.008466414245
+Iter number:66 k_eff=0.906841960366 flux ratio=4.384056022553 max_thermal=0.008466414245
+Iter number:67 k_eff=0.906841960368 flux ratio=4.384056022564 max_thermal=0.008466414246
+Iter number:68 k_eff=0.906841960369 flux ratio=4.384056022572 max_thermal=0.008466414246
+Iter number:69 k_eff=0.906841960370 flux ratio=4.384056022578 max_thermal=0.008466414246
+Iter number:70 k_eff=0.906841960371 flux ratio=4.384056022582 max_thermal=0.008466414246
+Cycle 1:
+ Numbers of active cells: 1156 2380
+ Numbers of degrees of freedom: 4761 10667
+Iter number:1 k_eff=0.906838267472 flux ratio=4.385474405124 max_thermal=0.008463675976
+Iter number:2 k_eff=0.906837892433 flux ratio=4.385486158840 max_thermal=0.008463675386
+Iter number:3 k_eff=0.906837848258 flux ratio=4.385487761080 max_thermal=0.008463681343
+Iter number:4 k_eff=0.906837849549 flux ratio=4.385488316617 max_thermal=0.008463685560
+Iter number:5 k_eff=0.906837859133 flux ratio=4.385488608157 max_thermal=0.008463688398
+Iter number:6 k_eff=0.906837869078 flux ratio=4.385488780872 max_thermal=0.008463690308
+Iter number:7 k_eff=0.906837877437 flux ratio=4.385488887750 max_thermal=0.008463691599
+Iter number:8 k_eff=0.906837883930 flux ratio=4.385488955514 max_thermal=0.008463692477
+Iter number:9 k_eff=0.906837888778 flux ratio=4.385488999282 max_thermal=0.008463693077
+Iter number:10 k_eff=0.906837892317 flux ratio=4.385489027995 max_thermal=0.008463693490
+Iter number:11 k_eff=0.906837894865 flux ratio=4.385489047078 max_thermal=0.008463693775
+Iter number:12 k_eff=0.906837896682 flux ratio=4.385489059900 max_thermal=0.008463693972
+Iter number:13 k_eff=0.906837897972 flux ratio=4.385489068591 max_thermal=0.008463694108
+Iter number:14 k_eff=0.906837898882 flux ratio=4.385489074523 max_thermal=0.008463694203
+Iter number:15 k_eff=0.906837899524 flux ratio=4.385489078594 max_thermal=0.008463694269
+Iter number:16 k_eff=0.906837899975 flux ratio=4.385489081400 max_thermal=0.008463694315
+Iter number:17 k_eff=0.906837900292 flux ratio=4.385489083340 max_thermal=0.008463694347
+Iter number:18 k_eff=0.906837900514 flux ratio=4.385489084686 max_thermal=0.008463694369
+Iter number:19 k_eff=0.906837900670 flux ratio=4.385489085620 max_thermal=0.008463694385
+Iter number:20 k_eff=0.906837900779 flux ratio=4.385489086270 max_thermal=0.008463694396
+Iter number:21 k_eff=0.906837900855 flux ratio=4.385489086723 max_thermal=0.008463694404
+Iter number:22 k_eff=0.906837900909 flux ratio=4.385489087039 max_thermal=0.008463694409
+Iter number:23 k_eff=0.906837900946 flux ratio=4.385489087259 max_thermal=0.008463694413
+Iter number:24 k_eff=0.906837900972 flux ratio=4.385489087413 max_thermal=0.008463694415
+Iter number:25 k_eff=0.906837900990 flux ratio=4.385489087521 max_thermal=0.008463694417
+Iter number:26 k_eff=0.906837901003 flux ratio=4.385489087596 max_thermal=0.008463694418
+Iter number:27 k_eff=0.906837901012 flux ratio=4.385489087648 max_thermal=0.008463694419
+Iter number:28 k_eff=0.906837901018 flux ratio=4.385489087685 max_thermal=0.008463694420
+Iter number:29 k_eff=0.906837901023 flux ratio=4.385489087710 max_thermal=0.008463694420
+Iter number:30 k_eff=0.906837901026 flux ratio=4.385489087728 max_thermal=0.008463694421
+Iter number:31 k_eff=0.906837901028 flux ratio=4.385489087741 max_thermal=0.008463694421
+Iter number:32 k_eff=0.906837901030 flux ratio=4.385489087749 max_thermal=0.008463694421
+Iter number:33 k_eff=0.906837901031 flux ratio=4.385489087755 max_thermal=0.008463694421
+Iter number:34 k_eff=0.906837901031 flux ratio=4.385489087760 max_thermal=0.008463694421
+Cycle 2:
+ Numbers of active cells: 1156 4219
+ Numbers of degrees of freedom: 4761 18805
+Iter number:1 k_eff=0.906836032131 flux ratio=4.385463219198 max_thermal=0.008463744346
+Iter number:2 k_eff=0.906835885276 flux ratio=4.385464355771 max_thermal=0.008463756043
+Iter number:3 k_eff=0.906835925468 flux ratio=4.385465086805 max_thermal=0.008463765603
+Iter number:4 k_eff=0.906835957285 flux ratio=4.385465633530 max_thermal=0.008463771996
+Iter number:5 k_eff=0.906835986095 flux ratio=4.385465985078 max_thermal=0.008463776366
+Iter number:6 k_eff=0.906836010004 flux ratio=4.385466206119 max_thermal=0.008463779390
+Iter number:7 k_eff=0.906836028489 flux ratio=4.385466348192 max_thermal=0.008463781494
+Iter number:8 k_eff=0.906836042190 flux ratio=4.385466442031 max_thermal=0.008463782961
+Iter number:9 k_eff=0.906836052101 flux ratio=4.385466505347 max_thermal=0.008463783987
+Iter number:10 k_eff=0.906836059171 flux ratio=4.385466548687 max_thermal=0.008463784704
+Iter number:11 k_eff=0.906836064172 flux ratio=4.385466578622 max_thermal=0.008463785205
+Iter number:12 k_eff=0.906836067692 flux ratio=4.385466599413 max_thermal=0.008463785556
+Iter number:13 k_eff=0.906836070164 flux ratio=4.385466613899 max_thermal=0.008463785801
+Iter number:14 k_eff=0.906836071896 flux ratio=4.385466624011 max_thermal=0.008463785972
+Iter number:15 k_eff=0.906836073108 flux ratio=4.385466631076 max_thermal=0.008463786092
+Iter number:16 k_eff=0.906836073957 flux ratio=4.385466636016 max_thermal=0.008463786176
+Iter number:17 k_eff=0.906836074550 flux ratio=4.385466639471 max_thermal=0.008463786235
+Iter number:18 k_eff=0.906836074965 flux ratio=4.385466641887 max_thermal=0.008463786276
+Iter number:19 k_eff=0.906836075255 flux ratio=4.385466643577 max_thermal=0.008463786305
+Iter number:20 k_eff=0.906836075458 flux ratio=4.385466644759 max_thermal=0.008463786325
+Iter number:21 k_eff=0.906836075600 flux ratio=4.385466645585 max_thermal=0.008463786339
+Iter number:22 k_eff=0.906836075699 flux ratio=4.385466646163 max_thermal=0.008463786349
+Iter number:23 k_eff=0.906836075769 flux ratio=4.385466646568 max_thermal=0.008463786356
+Iter number:24 k_eff=0.906836075817 flux ratio=4.385466646850 max_thermal=0.008463786360
+Iter number:25 k_eff=0.906836075851 flux ratio=4.385466647048 max_thermal=0.008463786364
+Iter number:26 k_eff=0.906836075875 flux ratio=4.385466647186 max_thermal=0.008463786366
+Iter number:27 k_eff=0.906836075891 flux ratio=4.385466647283 max_thermal=0.008463786368
+Iter number:28 k_eff=0.906836075903 flux ratio=4.385466647351 max_thermal=0.008463786369
+Iter number:29 k_eff=0.906836075911 flux ratio=4.385466647398 max_thermal=0.008463786370
+Iter number:30 k_eff=0.906836075917 flux ratio=4.385466647431 max_thermal=0.008463786370
+Iter number:31 k_eff=0.906836075921 flux ratio=4.385466647454 max_thermal=0.008463786371
+Iter number:32 k_eff=0.906836075924 flux ratio=4.385466647470 max_thermal=0.008463786371
+Iter number:33 k_eff=0.906836075926 flux ratio=4.385466647482 max_thermal=0.008463786371
+Iter number:34 k_eff=0.906836075927 flux ratio=4.385466647489 max_thermal=0.008463786371
+Iter number:35 k_eff=0.906836075928 flux ratio=4.385466647495 max_thermal=0.008463786371
+Iter number:36 k_eff=0.906836075928 flux ratio=4.385466647499 max_thermal=0.008463786371
+Cycle 3:
+ Numbers of active cells: 1507 6133
+ Numbers of degrees of freedom: 6629 27301
+Iter number:1 k_eff=0.906835269231 flux ratio=4.385403797661 max_thermal=0.008463904814
+Iter number:2 k_eff=0.906835336040 flux ratio=4.385404137890 max_thermal=0.008463908915
+Iter number:3 k_eff=0.906835394509 flux ratio=4.385404254946 max_thermal=0.008463912078
+Iter number:4 k_eff=0.906835430587 flux ratio=4.385404329325 max_thermal=0.008463914341
+Iter number:5 k_eff=0.906835453664 flux ratio=4.385404388953 max_thermal=0.008463915964
+Iter number:6 k_eff=0.906835468672 flux ratio=4.385404437567 max_thermal=0.008463917123
+Iter number:7 k_eff=0.906835478615 flux ratio=4.385404475817 max_thermal=0.008463917947
+Iter number:8 k_eff=0.906835485309 flux ratio=4.385404504893 max_thermal=0.008463918530
+Iter number:9 k_eff=0.906835489870 flux ratio=4.385404526447 max_thermal=0.008463918942
+Iter number:10 k_eff=0.906835493002 flux ratio=4.385404542148 max_thermal=0.008463919232
+Iter number:11 k_eff=0.906835495166 flux ratio=4.385404553451 max_thermal=0.008463919436
+Iter number:12 k_eff=0.906835496666 flux ratio=4.385404561520 max_thermal=0.008463919579
+Iter number:13 k_eff=0.906835497709 flux ratio=4.385404567248 max_thermal=0.008463919680
+Iter number:14 k_eff=0.906835498435 flux ratio=4.385404571298 max_thermal=0.008463919750
+Iter number:15 k_eff=0.906835498941 flux ratio=4.385404574153 max_thermal=0.008463919800
+Iter number:16 k_eff=0.906835499294 flux ratio=4.385404576162 max_thermal=0.008463919834
+Iter number:17 k_eff=0.906835499541 flux ratio=4.385404577572 max_thermal=0.008463919858
+Iter number:18 k_eff=0.906835499713 flux ratio=4.385404578562 max_thermal=0.008463919875
+Iter number:19 k_eff=0.906835499833 flux ratio=4.385404579256 max_thermal=0.008463919887
+Iter number:20 k_eff=0.906835499917 flux ratio=4.385404579742 max_thermal=0.008463919895
+Iter number:21 k_eff=0.906835499976 flux ratio=4.385404580083 max_thermal=0.008463919901
+Iter number:22 k_eff=0.906835500017 flux ratio=4.385404580321 max_thermal=0.008463919905
+Iter number:23 k_eff=0.906835500046 flux ratio=4.385404580488 max_thermal=0.008463919908
+Iter number:24 k_eff=0.906835500066 flux ratio=4.385404580604 max_thermal=0.008463919910
+Iter number:25 k_eff=0.906835500080 flux ratio=4.385404580686 max_thermal=0.008463919911
+Iter number:26 k_eff=0.906835500090 flux ratio=4.385404580743 max_thermal=0.008463919912
+Iter number:27 k_eff=0.906835500097 flux ratio=4.385404580783 max_thermal=0.008463919913
+Iter number:28 k_eff=0.906835500101 flux ratio=4.385404580811 max_thermal=0.008463919914
+Iter number:29 k_eff=0.906835500105 flux ratio=4.385404580831 max_thermal=0.008463919914
+Iter number:30 k_eff=0.906835500107 flux ratio=4.385404580844 max_thermal=0.008463919914
+Iter number:31 k_eff=0.906835500109 flux ratio=4.385404580854 max_thermal=0.008463919914
+Iter number:32 k_eff=0.906835500110 flux ratio=4.385404580860 max_thermal=0.008463919914
+Iter number:33 k_eff=0.906835500111 flux ratio=4.385404580865 max_thermal=0.008463919914
+Cycle 4:
+ Numbers of active cells: 2734 10732
+ Numbers of degrees of freedom: 12263 48095
+Iter number:1 k_eff=0.906834846364 flux ratio=4.385381150927 max_thermal=0.008463963132
+Iter number:2 k_eff=0.906834885938 flux ratio=4.385381403919 max_thermal=0.008463966541
+Iter number:3 k_eff=0.906834926531 flux ratio=4.385381460228 max_thermal=0.008463969391
+Iter number:4 k_eff=0.906834951140 flux ratio=4.385381536948 max_thermal=0.008463971436
+Iter number:5 k_eff=0.906834967134 flux ratio=4.385381607000 max_thermal=0.008463972881
+Iter number:6 k_eff=0.906834977868 flux ratio=4.385381661804 max_thermal=0.008463973895
+Iter number:7 k_eff=0.906834985190 flux ratio=4.385381702101 max_thermal=0.008463974604
+Iter number:8 k_eff=0.906834990235 flux ratio=4.385381730930 max_thermal=0.008463975099
+Iter number:9 k_eff=0.906834993731 flux ratio=4.385381751290 max_thermal=0.008463975445
+Iter number:10 k_eff=0.906834996164 flux ratio=4.385381765579 max_thermal=0.008463975687
+Iter number:11 k_eff=0.906834997860 flux ratio=4.385381775578 max_thermal=0.008463975856
+Iter number:12 k_eff=0.906834999045 flux ratio=4.385381782564 max_thermal=0.008463975973
+Iter number:13 k_eff=0.906834999873 flux ratio=4.385381787442 max_thermal=0.008463976056
+Iter number:14 k_eff=0.906835000452 flux ratio=4.385381790849 max_thermal=0.008463976113
+Iter number:15 k_eff=0.906835000856 flux ratio=4.385381793228 max_thermal=0.008463976153
+Iter number:16 k_eff=0.906835001140 flux ratio=4.385381794889 max_thermal=0.008463976182
+Iter number:17 k_eff=0.906835001338 flux ratio=4.385381796049 max_thermal=0.008463976201
+Iter number:18 k_eff=0.906835001476 flux ratio=4.385381796860 max_thermal=0.008463976215
+Iter number:19 k_eff=0.906835001573 flux ratio=4.385381797427 max_thermal=0.008463976224
+Iter number:20 k_eff=0.906835001641 flux ratio=4.385381797823 max_thermal=0.008463976231
+Iter number:21 k_eff=0.906835001688 flux ratio=4.385381798099 max_thermal=0.008463976236
+Iter number:22 k_eff=0.906835001721 flux ratio=4.385381798293 max_thermal=0.008463976239
+Iter number:23 k_eff=0.906835001744 flux ratio=4.385381798428 max_thermal=0.008463976241
+Iter number:24 k_eff=0.906835001760 flux ratio=4.385381798523 max_thermal=0.008463976243
+Iter number:25 k_eff=0.906835001772 flux ratio=4.385381798589 max_thermal=0.008463976244
+Iter number:26 k_eff=0.906835001780 flux ratio=4.385381798635 max_thermal=0.008463976245
+Iter number:27 k_eff=0.906835001785 flux ratio=4.385381798667 max_thermal=0.008463976246
+Iter number:28 k_eff=0.906835001789 flux ratio=4.385381798690 max_thermal=0.008463976246
+Iter number:29 k_eff=0.906835001792 flux ratio=4.385381798706 max_thermal=0.008463976246
+Iter number:30 k_eff=0.906835001794 flux ratio=4.385381798717 max_thermal=0.008463976246
+Iter number:31 k_eff=0.906835001795 flux ratio=4.385381798724 max_thermal=0.008463976246
+Iter number:32 k_eff=0.906835001796 flux ratio=4.385381798730 max_thermal=0.008463976247
+Iter number:33 k_eff=0.906835001796 flux ratio=4.385381798734 max_thermal=0.008463976247
+Cycle 5:
+ Numbers of active cells: 3928 15598
+ Numbers of degrees of freedom: 17501 69297
+Iter number:1 k_eff=0.906834756419 flux ratio=4.384853631027 max_thermal=0.008464995625
+Iter number:2 k_eff=0.906834805316 flux ratio=4.384853670539 max_thermal=0.008464997340
+Iter number:3 k_eff=0.906834826349 flux ratio=4.384853696977 max_thermal=0.008464998471
+Iter number:4 k_eff=0.906834837235 flux ratio=4.384853725424 max_thermal=0.008464999257
+Iter number:5 k_eff=0.906834843989 flux ratio=4.384853750658 max_thermal=0.008464999812
+Iter number:6 k_eff=0.906834848438 flux ratio=4.384853770680 max_thermal=0.008465000204
+Iter number:7 k_eff=0.906834851444 flux ratio=4.384853785717 max_thermal=0.008465000479
+Iter number:8 k_eff=0.906834853502 flux ratio=4.384853796687 max_thermal=0.008465000673
+Iter number:9 k_eff=0.906834854923 flux ratio=4.384853804560 max_thermal=0.008465000809
+Iter number:10 k_eff=0.906834855908 flux ratio=4.384853810156 max_thermal=0.008465000904
+Iter number:11 k_eff=0.906834856593 flux ratio=4.384853814111 max_thermal=0.008465000971
+Iter number:12 k_eff=0.906834857071 flux ratio=4.384853816895 max_thermal=0.008465001018
+Iter number:13 k_eff=0.906834857404 flux ratio=4.384853818850 max_thermal=0.008465001050
+Iter number:14 k_eff=0.906834857636 flux ratio=4.384853820222 max_thermal=0.008465001073
+Iter number:15 k_eff=0.906834857799 flux ratio=4.384853821183 max_thermal=0.008465001089
+Iter number:16 k_eff=0.906834857912 flux ratio=4.384853821856 max_thermal=0.008465001100
+Iter number:17 k_eff=0.906834857991 flux ratio=4.384853822327 max_thermal=0.008465001108
+Iter number:18 k_eff=0.906834858047 flux ratio=4.384853822657 max_thermal=0.008465001114
+Iter number:19 k_eff=0.906834858085 flux ratio=4.384853822888 max_thermal=0.008465001117
+Iter number:20 k_eff=0.906834858113 flux ratio=4.384853823049 max_thermal=0.008465001120
+Iter number:21 k_eff=0.906834858132 flux ratio=4.384853823162 max_thermal=0.008465001122
+Iter number:22 k_eff=0.906834858145 flux ratio=4.384853823241 max_thermal=0.008465001123
+Iter number:23 k_eff=0.906834858154 flux ratio=4.384853823296 max_thermal=0.008465001124
+Iter number:24 k_eff=0.906834858161 flux ratio=4.384853823335 max_thermal=0.008465001125
+Iter number:25 k_eff=0.906834858165 flux ratio=4.384853823362 max_thermal=0.008465001125
+Iter number:26 k_eff=0.906834858168 flux ratio=4.384853823380 max_thermal=0.008465001126
+Iter number:27 k_eff=0.906834858170 flux ratio=4.384853823394 max_thermal=0.008465001126
+Iter number:28 k_eff=0.906834858172 flux ratio=4.384853823403 max_thermal=0.008465001126
+Iter number:29 k_eff=0.906834858173 flux ratio=4.384853823409 max_thermal=0.008465001126
+Iter number:30 k_eff=0.906834858174 flux ratio=4.384853823414 max_thermal=0.008465001126
+Cycle 6:
+ Numbers of active cells: 4486 17755
+ Numbers of degrees of freedom: 19933 78605
+Iter number:1 k_eff=0.906834797396 flux ratio=4.384850615384 max_thermal=0.008465007174
+Iter number:2 k_eff=0.906834814661 flux ratio=4.384850631692 max_thermal=0.008465007547
+Iter number:3 k_eff=0.906834818206 flux ratio=4.384850637654 max_thermal=0.008465007740
+Iter number:4 k_eff=0.906834820171 flux ratio=4.384850642364 max_thermal=0.008465007876
+Iter number:5 k_eff=0.906834821415 flux ratio=4.384850646407 max_thermal=0.008465007974
+Iter number:6 k_eff=0.906834822244 flux ratio=4.384850649675 max_thermal=0.008465008044
+Iter number:7 k_eff=0.906834822807 flux ratio=4.384850652188 max_thermal=0.008465008093
+Iter number:8 k_eff=0.906834823192 flux ratio=4.384850654060 max_thermal=0.008465008128
+Iter number:9 k_eff=0.906834823457 flux ratio=4.384850655427 max_thermal=0.008465008153
+Iter number:10 k_eff=0.906834823641 flux ratio=4.384850656411 max_thermal=0.008465008170
+Iter number:11 k_eff=0.906834823769 flux ratio=4.384850657114 max_thermal=0.008465008182
+Iter number:12 k_eff=0.906834823858 flux ratio=4.384850657613 max_thermal=0.008465008191
+Iter number:13 k_eff=0.906834823919 flux ratio=4.384850657966 max_thermal=0.008465008197
+Iter number:14 k_eff=0.906834823962 flux ratio=4.384850658214 max_thermal=0.008465008201
+Iter number:15 k_eff=0.906834823992 flux ratio=4.384850658389 max_thermal=0.008465008204
+Iter number:16 k_eff=0.906834824013 flux ratio=4.384850658512 max_thermal=0.008465008206
+Iter number:17 k_eff=0.906834824028 flux ratio=4.384850658598 max_thermal=0.008465008207
+Iter number:18 k_eff=0.906834824038 flux ratio=4.384850658659 max_thermal=0.008465008208
+Iter number:19 k_eff=0.906834824046 flux ratio=4.384850658701 max_thermal=0.008465008209
+Iter number:20 k_eff=0.906834824051 flux ratio=4.384850658731 max_thermal=0.008465008209
+Iter number:21 k_eff=0.906834824054 flux ratio=4.384850658752 max_thermal=0.008465008210
+Iter number:22 k_eff=0.906834824057 flux ratio=4.384850658766 max_thermal=0.008465008210
+Iter number:23 k_eff=0.906834824058 flux ratio=4.384850658776 max_thermal=0.008465008210
+Iter number:24 k_eff=0.906834824059 flux ratio=4.384850658783 max_thermal=0.008465008210
+Iter number:25 k_eff=0.906834824060 flux ratio=4.384850658788 max_thermal=0.008465008210
+Cycle 7:
+ Numbers of active cells: 5434 21370
+ Numbers of degrees of freedom: 23979 93275
+Iter number:1 k_eff=0.906834695333 flux ratio=4.384848325238 max_thermal=0.008465012709
+Iter number:2 k_eff=0.906834781972 flux ratio=4.384848347736 max_thermal=0.008465013711
+Iter number:3 k_eff=0.906834782962 flux ratio=4.384848354071 max_thermal=0.008465013889
+Iter number:4 k_eff=0.906834784295 flux ratio=4.384848360506 max_thermal=0.008465014025
+Iter number:5 k_eff=0.906834785267 flux ratio=4.384848365788 max_thermal=0.008465014122
+Iter number:6 k_eff=0.906834785953 flux ratio=4.384848369732 max_thermal=0.008465014190
+Iter number:7 k_eff=0.906834786435 flux ratio=4.384848372568 max_thermal=0.008465014237
+Iter number:8 k_eff=0.906834786773 flux ratio=4.384848374574 max_thermal=0.008465014271
+Iter number:9 k_eff=0.906834787009 flux ratio=4.384848375984 max_thermal=0.008465014294
+Iter number:10 k_eff=0.906834787174 flux ratio=4.384848376972 max_thermal=0.008465014310
+Iter number:11 k_eff=0.906834787289 flux ratio=4.384848377663 max_thermal=0.008465014322
+Iter number:12 k_eff=0.906834787370 flux ratio=4.384848378146 max_thermal=0.008465014330
+Iter number:13 k_eff=0.906834787426 flux ratio=4.384848378484 max_thermal=0.008465014335
+Iter number:14 k_eff=0.906834787466 flux ratio=4.384848378720 max_thermal=0.008465014339
+Iter number:15 k_eff=0.906834787493 flux ratio=4.384848378885 max_thermal=0.008465014342
+Iter number:16 k_eff=0.906834787512 flux ratio=4.384848379000 max_thermal=0.008465014344
+Iter number:17 k_eff=0.906834787526 flux ratio=4.384848379081 max_thermal=0.008465014345
+Iter number:18 k_eff=0.906834787535 flux ratio=4.384848379137 max_thermal=0.008465014346
+Iter number:19 k_eff=0.906834787542 flux ratio=4.384848379176 max_thermal=0.008465014347
+Iter number:20 k_eff=0.906834787547 flux ratio=4.384848379204 max_thermal=0.008465014347
+Iter number:21 k_eff=0.906834787550 flux ratio=4.384848379223 max_thermal=0.008465014348
+Iter number:22 k_eff=0.906834787552 flux ratio=4.384848379237 max_thermal=0.008465014348
+Iter number:23 k_eff=0.906834787554 flux ratio=4.384848379246 max_thermal=0.008465014348
+Iter number:24 k_eff=0.906834787555 flux ratio=4.384848379253 max_thermal=0.008465014348
+Iter number:25 k_eff=0.906834787556 flux ratio=4.384848379257 max_thermal=0.008465014348
+Cycle 8:
+ Numbers of active cells: 6856 27001
+ Numbers of degrees of freedom: 30285 117017
+Iter number:1 k_eff=0.906834743244 flux ratio=4.384846479257 max_thermal=0.008465017253
+Iter number:2 k_eff=0.906834753823 flux ratio=4.384846519298 max_thermal=0.008465017557
+Iter number:3 k_eff=0.906834756845 flux ratio=4.384846524675 max_thermal=0.008465017731
+Iter number:4 k_eff=0.906834758500 flux ratio=4.384846529782 max_thermal=0.008465017851
+Iter number:5 k_eff=0.906834759516 flux ratio=4.384846533887 max_thermal=0.008465017934
+Iter number:6 k_eff=0.906834760177 flux ratio=4.384846537011 max_thermal=0.008465017993
+Iter number:7 k_eff=0.906834760619 flux ratio=4.384846539312 max_thermal=0.008465018034
+Iter number:8 k_eff=0.906834760921 flux ratio=4.384846540972 max_thermal=0.008465018062
+Iter number:9 k_eff=0.906834761129 flux ratio=4.384846542155 max_thermal=0.008465018083
+Iter number:10 k_eff=0.906834761274 flux ratio=4.384846542991 max_thermal=0.008465018097
+Iter number:11 k_eff=0.906834761374 flux ratio=4.384846543580 max_thermal=0.008465018106
+Iter number:12 k_eff=0.906834761444 flux ratio=4.384846543993 max_thermal=0.008465018113
+Iter number:13 k_eff=0.906834761493 flux ratio=4.384846544282 max_thermal=0.008465018118
+Iter number:14 k_eff=0.906834761527 flux ratio=4.384846544485 max_thermal=0.008465018121
+Iter number:15 k_eff=0.906834761551 flux ratio=4.384846544627 max_thermal=0.008465018124
+Iter number:16 k_eff=0.906834761567 flux ratio=4.384846544726 max_thermal=0.008465018125
+Iter number:17 k_eff=0.906834761579 flux ratio=4.384846544795 max_thermal=0.008465018127
+Iter number:18 k_eff=0.906834761587 flux ratio=4.384846544844 max_thermal=0.008465018127
+Iter number:19 k_eff=0.906834761593 flux ratio=4.384846544878 max_thermal=0.008465018128
+Iter number:20 k_eff=0.906834761597 flux ratio=4.384846544901 max_thermal=0.008465018128
+Iter number:21 k_eff=0.906834761599 flux ratio=4.384846544918 max_thermal=0.008465018129
+Iter number:22 k_eff=0.906834761601 flux ratio=4.384846544930 max_thermal=0.008465018129
+Iter number:23 k_eff=0.906834761603 flux ratio=4.384846544938 max_thermal=0.008465018129
+Iter number:24 k_eff=0.906834761604 flux ratio=4.384846544944 max_thermal=0.008465018129
+Iter number:25 k_eff=0.906834761604 flux ratio=4.384846544947 max_thermal=0.008465018129
+Cycle 9:
+ Numbers of active cells: 9166 35416
+ Numbers of degrees of freedom: 40087 154355
+Iter number:1 k_eff=0.906834746805 flux ratio=4.384846149270 max_thermal=0.008465018800
+Iter number:2 k_eff=0.906834739570 flux ratio=4.384846071221 max_thermal=0.008465018869
+Iter number:3 k_eff=0.906834742439 flux ratio=4.384846070361 max_thermal=0.008465018998
+Iter number:4 k_eff=0.906834743795 flux ratio=4.384846072627 max_thermal=0.008465019087
+Iter number:5 k_eff=0.906834744601 flux ratio=4.384846075162 max_thermal=0.008465019149
+Iter number:6 k_eff=0.906834745118 flux ratio=4.384846077315 max_thermal=0.008465019194
+Iter number:7 k_eff=0.906834745462 flux ratio=4.384846078976 max_thermal=0.008465019225
+Iter number:8 k_eff=0.906834745695 flux ratio=4.384846080201 max_thermal=0.008465019246
+Iter number:9 k_eff=0.906834745854 flux ratio=4.384846081084 max_thermal=0.008465019262
+Iter number:10 k_eff=0.906834745964 flux ratio=4.384846081712 max_thermal=0.008465019272
+Iter number:11 k_eff=0.906834746041 flux ratio=4.384846082157 max_thermal=0.008465019280
+Iter number:12 k_eff=0.906834746094 flux ratio=4.384846082469 max_thermal=0.008465019285
+Iter number:13 k_eff=0.906834746131 flux ratio=4.384846082689 max_thermal=0.008465019289
+Iter number:14 k_eff=0.906834746157 flux ratio=4.384846082842 max_thermal=0.008465019291
+Iter number:15 k_eff=0.906834746175 flux ratio=4.384846082950 max_thermal=0.008465019293
+Iter number:16 k_eff=0.906834746188 flux ratio=4.384846083025 max_thermal=0.008465019294
+Iter number:17 k_eff=0.906834746197 flux ratio=4.384846083078 max_thermal=0.008465019295
+Iter number:18 k_eff=0.906834746203 flux ratio=4.384846083115 max_thermal=0.008465019296
+Iter number:19 k_eff=0.906834746207 flux ratio=4.384846083141 max_thermal=0.008465019296
+Iter number:20 k_eff=0.906834746210 flux ratio=4.384846083159 max_thermal=0.008465019296
+Iter number:21 k_eff=0.906834746213 flux ratio=4.384846083171 max_thermal=0.008465019297
+Iter number:22 k_eff=0.906834746214 flux ratio=4.384846083180 max_thermal=0.008465019297
+Iter number:23 k_eff=0.906834746215 flux ratio=4.384846083186 max_thermal=0.008465019297
+Iter number:24 k_eff=0.906834746216 flux ratio=4.384846083191 max_thermal=0.008465019297
+Cycle 10:
+ Numbers of active cells: 10594 41197
+ Numbers of degrees of freedom: 45467 179469
+Iter number:1 k_eff=0.906780038935 flux ratio=4.384864081404 max_thermal=0.008464912703
+Iter number:2 k_eff=0.906833758951 flux ratio=4.384869346165 max_thermal=0.008465298273
+Iter number:3 k_eff=0.906834914205 flux ratio=4.384864425562 max_thermal=0.008465215350
+Iter number:4 k_eff=0.906835274876 flux ratio=4.384858287964 max_thermal=0.008465154357
+Iter number:5 k_eff=0.906835297124 flux ratio=4.384853804292 max_thermal=0.008465112059
+Iter number:6 k_eff=0.906835209298 flux ratio=4.384850953669 max_thermal=0.008465083180
+Iter number:7 k_eff=0.906835102972 flux ratio=4.384849195436 max_thermal=0.008465063473
+Iter number:8 k_eff=0.906835009160 flux ratio=4.384848100961 max_thermal=0.008465049973
+Iter number:9 k_eff=0.906834935068 flux ratio=4.384847404790 max_thermal=0.008465040684
+Iter number:10 k_eff=0.906834879514 flux ratio=4.384846951907 max_thermal=0.008465034267
+Iter number:11 k_eff=0.906834838999 flux ratio=4.384846651620 max_thermal=0.008465029821
+Iter number:12 k_eff=0.906834809916 flux ratio=4.384846449558 max_thermal=0.008465026733
+Iter number:13 k_eff=0.906834789237 flux ratio=4.384846312108 max_thermal=0.008465024585
+Iter number:14 k_eff=0.906834774618 flux ratio=4.384846217876 max_thermal=0.008465023089
+Iter number:15 k_eff=0.906834764321 flux ratio=4.384846152915 max_thermal=0.008465022045
+Iter number:16 k_eff=0.906834757086 flux ratio=4.384846107956 max_thermal=0.008465021317
+Iter number:17 k_eff=0.906834752010 flux ratio=4.384846076754 max_thermal=0.008465020809
+Iter number:18 k_eff=0.906834748452 flux ratio=4.384846055055 max_thermal=0.008465020454
+Iter number:19 k_eff=0.906834745960 flux ratio=4.384846039944 max_thermal=0.008465020207
+Iter number:20 k_eff=0.906834744216 flux ratio=4.384846029409 max_thermal=0.008465020033
+Iter number:21 k_eff=0.906834742995 flux ratio=4.384846022060 max_thermal=0.008465019912
+Iter number:22 k_eff=0.906834742141 flux ratio=4.384846016929 max_thermal=0.008465019828
+Iter number:23 k_eff=0.906834741544 flux ratio=4.384846013347 max_thermal=0.008465019768
+Iter number:24 k_eff=0.906834741126 flux ratio=4.384846010844 max_thermal=0.008465019727
+Iter number:25 k_eff=0.906834740834 flux ratio=4.384846009095 max_thermal=0.008465019698
+Iter number:26 k_eff=0.906834740629 flux ratio=4.384846007873 max_thermal=0.008465019678
+Iter number:27 k_eff=0.906834740486 flux ratio=4.384846007019 max_thermal=0.008465019664
+Iter number:28 k_eff=0.906834740386 flux ratio=4.384846006422 max_thermal=0.008465019654
+Iter number:29 k_eff=0.906834740316 flux ratio=4.384846006004 max_thermal=0.008465019647
+Iter number:30 k_eff=0.906834740267 flux ratio=4.384846005712 max_thermal=0.008465019642
+Iter number:31 k_eff=0.906834740233 flux ratio=4.384846005508 max_thermal=0.008465019639
+Iter number:32 k_eff=0.906834740209 flux ratio=4.384846005366 max_thermal=0.008465019637
+Iter number:33 k_eff=0.906834740193 flux ratio=4.384846005266 max_thermal=0.008465019635
+Iter number:34 k_eff=0.906834740181 flux ratio=4.384846005196 max_thermal=0.008465019634
+Iter number:35 k_eff=0.906834740173 flux ratio=4.384846005148 max_thermal=0.008465019633
+Iter number:36 k_eff=0.906834740167 flux ratio=4.384846005113 max_thermal=0.008465019632
+Iter number:37 k_eff=0.906834740163 flux ratio=4.384846005090 max_thermal=0.008465019632
+Iter number:38 k_eff=0.906834740160 flux ratio=4.384846005073 max_thermal=0.008465019632
+Iter number:39 k_eff=0.906834740158 flux ratio=4.384846005061 max_thermal=0.008465019632
+Iter number:40 k_eff=0.906834740157 flux ratio=4.384846005053 max_thermal=0.008465019631
+Iter number:41 k_eff=0.906834740156 flux ratio=4.384846005047 max_thermal=0.008465019631
+Iter number:42 k_eff=0.906834740155 flux ratio=4.384846005044 max_thermal=0.008465019631
+Cycle 11:
+ Numbers of active cells: 11749 47074
+ Numbers of degrees of freedom: 50261 204523
+Iter number:1 k_eff=0.906805395149 flux ratio=4.384872231023 max_thermal=0.008464861813
+Iter number:2 k_eff=0.906833353627 flux ratio=4.384863725577 max_thermal=0.008465049652
+Iter number:3 k_eff=0.906834008279 flux ratio=4.384854144970 max_thermal=0.008465027770
+Iter number:4 k_eff=0.906834374542 flux ratio=4.384849413955 max_thermal=0.008465020786
+Iter number:5 k_eff=0.906834559132 flux ratio=4.384847412067 max_thermal=0.008465018746
+Iter number:6 k_eff=0.906834648393 flux ratio=4.384846595090 max_thermal=0.008465018343
+Iter number:7 k_eff=0.906834691266 flux ratio=4.384846266355 max_thermal=0.008465018450
+Iter number:8 k_eff=0.906834712193 flux ratio=4.384846137013 max_thermal=0.008465018681
+Iter number:9 k_eff=0.906834722750 flux ratio=4.384846088708 max_thermal=0.008465018907
+Iter number:10 k_eff=0.906834728325 flux ratio=4.384846072840 max_thermal=0.008465019094
+Iter number:11 k_eff=0.906834731426 flux ratio=4.384846069459 max_thermal=0.008465019237
+Iter number:12 k_eff=0.906834733245 flux ratio=4.384846070461 max_thermal=0.008465019343
+Iter number:13 k_eff=0.906834734362 flux ratio=4.384846072648 max_thermal=0.008465019420
+Iter number:14 k_eff=0.906834735076 flux ratio=4.384846074842 max_thermal=0.008465019475
+Iter number:15 k_eff=0.906834735545 flux ratio=4.384846076677 max_thermal=0.008465019514
+Iter number:16 k_eff=0.906834735859 flux ratio=4.384846078098 max_thermal=0.008465019542
+Iter number:17 k_eff=0.906834736073 flux ratio=4.384846079156 max_thermal=0.008465019561
+Iter number:18 k_eff=0.906834736220 flux ratio=4.384846079925 max_thermal=0.008465019575
+Iter number:19 k_eff=0.906834736321 flux ratio=4.384846080478 max_thermal=0.008465019585
+Iter number:20 k_eff=0.906834736391 flux ratio=4.384846080872 max_thermal=0.008465019591
+Iter number:21 k_eff=0.906834736440 flux ratio=4.384846081150 max_thermal=0.008465019596
+Iter number:22 k_eff=0.906834736474 flux ratio=4.384846081347 max_thermal=0.008465019599
+Iter number:23 k_eff=0.906834736498 flux ratio=4.384846081485 max_thermal=0.008465019602
+Iter number:24 k_eff=0.906834736515 flux ratio=4.384846081582 max_thermal=0.008465019603
+Iter number:25 k_eff=0.906834736526 flux ratio=4.384846081650 max_thermal=0.008465019604
+Iter number:26 k_eff=0.906834736534 flux ratio=4.384846081698 max_thermal=0.008465019605
+Iter number:27 k_eff=0.906834736540 flux ratio=4.384846081731 max_thermal=0.008465019606
+Iter number:28 k_eff=0.906834736544 flux ratio=4.384846081755 max_thermal=0.008465019606
+Iter number:29 k_eff=0.906834736547 flux ratio=4.384846081771 max_thermal=0.008465019606
+Iter number:30 k_eff=0.906834736549 flux ratio=4.384846081783 max_thermal=0.008465019607
+Iter number:31 k_eff=0.906834736550 flux ratio=4.384846081791 max_thermal=0.008465019607
+Iter number:32 k_eff=0.906834736551 flux ratio=4.384846081796 max_thermal=0.008465019607
+Iter number:33 k_eff=0.906834736552 flux ratio=4.384846081800 max_thermal=0.008465019607
+@endcode
+
+We see that power iteration does converge faster after cycle 0 due to the initialization
+with solution from last mesh iteration.
+The contents of ``convergence_table'' are,
+@code
+0 4761 4761 0.906841960371 4.384056022582
+1 4761 10667 0.906837901031 4.385489087760
+2 4761 18805 0.906836075928 4.385466647499
+3 6629 27301 0.906835500111 4.385404580865
+4 12263 48095 0.906835001796 4.385381798734
+5 17501 69297 0.906834858174 4.384853823414
+6 19933 78605 0.906834824060 4.384850658788
+7 23979 93275 0.906834787556 4.384848379257
+8 30285 117017 0.906834761604 4.384846544947
+9 40087 154355 0.906834746216 4.384846083191
+10 45467 179469 0.906834740155 4.384846005044
+11 50261 204523 0.906834736552 4.384846081800
+@endcode
+The meanings of coloms are: number of mesh iteration, numbers of degrees of
+ freedom of fast energy group, numbers of DoFs of thermal group, converged
+k-effective and the ratio between maximum of fast flux and maximum of thermal one.
+
+The grids of fast and thermal energy groups at mesh iteration #9 are shown
+in following figure.
+
+@image html step-28.grid-0.9.order2.png
+@image html step-28.grid-1.9.order2.png
+
+We see that the grid of thermal group is much finner than the one of fast group.
+The solutions on these grids are, (Note: flux are normalized with total fission
+source equal to 1)
+
+@image html step-28.solution-0.9.order2.png
+@image html step-28.solution-1.9.order2.png
+
+Then we plot the convergence data with polynomial order being equal to 1,2 and 3.
+
+@image html step-28.convergence.png
+
+The estimated ``exact'' k-effective = 0.906834721253 which is simply from last
+mesh iteration of polynomial order 3 minus 2e-10. We see that h-adaptive calculations
+deliver an algebraic convergence. And the higher polynomial order is, the faster mesh
+iteration converges. In our problem, we need smaller number of DoFs to achieve same
+accuracy with higher polynoimal order.
--- /dev/null
+//TODO: remove direct (non-eigenvalue) problem possibility (remove
+// ExtraneousSource and some stuff in run(), see isour)
+
+/* $Id: step-6.cc,v 1.43 2006/02/09 02:54:11 wolf Exp $ */
+/* Version: $Name: $ */
+/* */
+/* Copyright (C) 2006 by the deal.II authors and Yaqi Wang */
+/* */
+/* This file is subject to QPL and may not be distributed */
+/* without copyright and license information. Please refer */
+/* to the file deal.II/doc/license.html for the text and */
+/* further information on this license. */
+
+ // @sect3{Include files}
+
+ // We start with a bunch of include
+ // files that have already been
+ // explained in previous tutorial
+ // programs:
+#include <base/timer.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/thread_management.h>
+#include <base/parameter_handler.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_refinement.h>
+#include <grid/grid_out.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+#include <iostream>
+
+#include <base/utilities.h>
+
+ // We use the next include file to
+ // access block vectors which provide
+ // us a convenient way to manage
+ // solution and right hand side
+ // vectors of all energy groups:
+#include <lac/block_vector.h>
+
+ // This include file is for
+ // transferring solutions from one
+ // mesh to another different mesh. We
+ // use it when we are initializing
+ // solutions after each mesh
+ // iteration:
+#include <numerics/solution_transfer.h>
+
+ // When integrating functions defined
+ // on one mesh against shape
+ // functions defined on a different
+ // mesh, we need a function @p
+ // get_finest_common_cells (as
+ // discussed in the introduction)
+ // which is defined in the following
+ // header file:
+#include <grid/grid_tools.h>
+
+ // Here are two more C++ standard
+ // headers that we use to define list
+ // data types as well as to fine-tune
+ // the output we generate:
+#include <list>
+#include <iomanip>
+
+ // The last step is as in all
+ // previous programs:
+using namespace dealii;
+
+
+ // @sect3{Material data}
+
+ // First up, we need to define a
+ // class that provides material data
+ // (including diffusion coefficients,
+ // removal cross sections, scattering
+ // cross sections, fission cross
+ // sections and fission spectra) to
+ // the main class.
+ //
+ // The parameter to the constructor
+ // determines for how many energy
+ // groups we set up the relevant
+ // tables. At present, this program
+ // only includes data for 2 energy
+ // groups, but a more sophisticated
+ // program may be able to initialize
+ // the data structures for more
+ // groups as well, depending on how
+ // many energy groups are selected in
+ // the parameter file.
+ //
+ // For each of the different
+ // coefficient types, there is one
+ // function that returns the value of
+ // this coefficient for a particular
+ // energy group (or combination of
+ // energy groups, as for the
+ // distribution cross section
+ // $\chi_g\nu\Sigma_{f,g'}$ or
+ // scattering cross section
+ // $\Sigma_{s,g'\to g}$). In addition
+ // to the energy group or groups,
+ // these coefficients depend on the
+ // type of fuel or control rod, as
+ // explained in the introduction. The
+ // functions therefore take an
+ // additional parameter, @p
+ // material_id, that identifies the
+ // particular kind of rod. Within
+ // this program, we use
+ // <code>n_materials=8</code>
+ // different kinds of rods.
+ //
+ // Except for the scattering cross
+ // section, each of the coefficients
+ // therefore can be represented as an
+ // entry in a two-dimensional array
+ // of floating point values indexed
+ // by the energy group number as well
+ // as the material ID. The Table
+ // class template is the ideal way to
+ // store such data. Finally, the
+ // scattering coefficient depends on
+ // both two energy group indices and
+ // therefore needs to be stored in a
+ // three-dimensional array, for which
+ // we again use the Table class,
+ // where this time the first template
+ // argument (denoting the
+ // dimensionality of the array) of
+ // course needs to be three:
+class MaterialData
+{
+ public:
+ MaterialData (const unsigned int n_groups);
+
+ double get_diffusion_coefficient (const unsigned int group,
+ const unsigned int material_id) const;
+ double get_removal_XS (const unsigned int group,
+ const unsigned int material_id) const;
+ double get_fission_XS (const unsigned int group,
+ const unsigned int material_id) const;
+ double get_fission_dist_XS (const unsigned int group_1,
+ const unsigned int group_2,
+ const unsigned int material_id) const;
+ double get_scattering_XS (const unsigned int group_1,
+ const unsigned int group_2,
+ const unsigned int material_id) const;
+ double get_fission_spectrum (const unsigned int group,
+ const unsigned int material_id) const;
+
+ private:
+ const unsigned int n_groups;
+ const unsigned int n_materials;
+
+ Table<2,double> diffusion;
+ Table<2,double> sigma_r;
+ Table<2,double> nu_sigma_f;
+ Table<3,double> sigma_s;
+ Table<2,double> chi;
+};
+
+ // The constructor of the class is
+ // used to initialize all the
+ // material data arrays. It takes the
+ // number of energy groups as an
+ // argument (an throws an error if
+ // that value is not equal to two,
+ // since at presently only data for
+ // two energy groups is implemented;
+ // however, using this, the function
+ // remains flexible and extendible
+ // into the future). In the member
+ // initialization part at the
+ // beginning, it also resizes the
+ // arrays to their correct sizes.
+ //
+ // At present, material data is
+ // stored for 8 different types of
+ // material. This, as well, may
+ // easily be extended in the future.
+MaterialData::MaterialData (const unsigned int n_groups)
+ :
+ n_groups (n_groups),
+ n_materials (8),
+ diffusion (n_materials, n_groups),
+ sigma_r (n_materials, n_groups),
+ nu_sigma_f (n_materials, n_groups),
+ sigma_s (n_materials, n_groups, n_groups),
+ chi (n_materials, n_groups)
+{
+ switch (n_groups)
+ {
+ case 2:
+ {
+ for (unsigned int m=0; m<n_materials; ++m)
+ {
+ diffusion[m][0] = 1.2;
+ diffusion[m][1] = 0.4;
+ chi[m][0] = 1.0;
+ chi[m][1] = 0.0;
+ sigma_r[m][0] = 0.03;
+ for (unsigned int group_1=0; group_1<n_groups; ++group_1)
+ for (unsigned int group_2=0; group_2<n_groups; ++ group_2)
+ sigma_s[m][group_1][group_2] = 0.0;
+ }
+
+
+ diffusion[5][1] = 0.2;
+
+ sigma_r[4][0] = 0.026;
+ sigma_r[5][0] = 0.051;
+ sigma_r[6][0] = 0.026;
+ sigma_r[7][0] = 0.050;
+
+ sigma_r[0][1] = 0.100;
+ sigma_r[1][1] = 0.200;
+ sigma_r[2][1] = 0.250;
+ sigma_r[3][1] = 0.300;
+ sigma_r[4][1] = 0.020;
+ sigma_r[5][1] = 0.040;
+ sigma_r[6][1] = 0.020;
+ sigma_r[7][1] = 0.800;
+
+ nu_sigma_f[0][0] = 0.0050;
+ nu_sigma_f[1][0] = 0.0075;
+ nu_sigma_f[2][0] = 0.0075;
+ nu_sigma_f[3][0] = 0.0075;
+ nu_sigma_f[4][0] = 0.000;
+ nu_sigma_f[5][0] = 0.000;
+ nu_sigma_f[6][0] = 1e-7;
+ nu_sigma_f[7][0] = 0.00;
+
+ nu_sigma_f[0][1] = 0.125;
+ nu_sigma_f[1][1] = 0.300;
+ nu_sigma_f[2][1] = 0.375;
+ nu_sigma_f[3][1] = 0.450;
+ nu_sigma_f[4][1] = 0.000;
+ nu_sigma_f[5][1] = 0.000;
+ nu_sigma_f[6][1] = 3e-6;
+ nu_sigma_f[7][1] = 0.00;
+
+ sigma_s[0][0][1] = 0.020;
+ sigma_s[1][0][1] = 0.015;
+ sigma_s[2][0][1] = 0.015;
+ sigma_s[3][0][1] = 0.015;
+ sigma_s[4][0][1] = 0.025;
+ sigma_s[5][0][1] = 0.050;
+ sigma_s[6][0][1] = 0.025;
+ sigma_s[7][0][1] = 0.010;
+
+ break;
+ }
+
+
+ default:
+ Assert (false,
+ ExcMessage ("Presently, only data for 2 groups is implemented"));
+ }
+}
+
+
+ // Next are the functions that return
+ // the coefficient values for given
+ // materials and energy groups. All
+ // they do is to make sure that the
+ // given arguments are within the
+ // allowed ranges, and then look the
+ // respective value up in the
+ // corresponding tables:
+double
+MaterialData::get_diffusion_coefficient (const unsigned int group,
+ const unsigned int material_id) const
+{
+ Assert (group < n_groups,
+ ExcIndexRange (group, 0, n_groups));
+ Assert (material_id < n_materials,
+ ExcIndexRange (material_id, 0, n_materials));
+
+ return diffusion[material_id][group];
+}
+
+
+
+double
+MaterialData::get_removal_XS (const unsigned int group,
+ const unsigned int material_id) const
+{
+ Assert (group < n_groups,
+ ExcIndexRange (group, 0, n_groups));
+ Assert (material_id < n_materials,
+ ExcIndexRange (material_id, 0, n_materials));
+
+ return sigma_r[material_id][group];
+}
+
+
+double
+MaterialData::get_fission_XS (const unsigned int group,
+ const unsigned int material_id) const
+{
+ Assert (group < n_groups,
+ ExcIndexRange (group, 0, n_groups));
+ Assert (material_id < n_materials,
+ ExcIndexRange (material_id, 0, n_materials));
+
+ return nu_sigma_f[material_id][group];
+}
+
+
+
+double
+MaterialData::get_scattering_XS (const unsigned int group_1,
+ const unsigned int group_2,
+ const unsigned int material_id) const
+{
+ Assert (group_1 < n_groups,
+ ExcIndexRange (group_1, 0, n_groups));
+ Assert (group_2 < n_groups,
+ ExcIndexRange (group_2, 0, n_groups));
+ Assert (material_id < n_materials,
+ ExcIndexRange (material_id, 0, n_materials));
+
+ return sigma_s[material_id][group_1][group_2];
+}
+
+
+
+double
+MaterialData::get_fission_spectrum (const unsigned int group,
+ const unsigned int material_id) const
+{
+ Assert (group < n_groups,
+ ExcIndexRange (group, 0, n_groups));
+ Assert (material_id < n_materials,
+ ExcIndexRange (material_id, 0, n_materials));
+
+ return chi[material_id][group];
+}
+
+
+ // The function computing the fission
+ // distribution cross section is
+ // slightly different, since it
+ // computes its value as the product
+ // of two other coefficients. We
+ // don't need to check arguments
+ // here, since this already happens
+ // when we call the two other
+ // functions involved, even though it
+ // would probably not hurt either:
+double
+MaterialData::get_fission_dist_XS (const unsigned int group_1,
+ const unsigned int group_2,
+ const unsigned int material_id) const
+{
+ return (get_fission_spectrum(group_1, material_id) *
+ get_fission_XS(group_2, material_id));
+}
+
+
+
+ // @sect3{The <code>EnergyGroup</code> class}
+
+ // The first interesting class is the
+ // one that contains everything that
+ // is specific to a single energy
+ // group. To group things that belong
+ // together into individual objects,
+ // we declare a structure that holds
+ // the Triangulation and DoFHandler
+ // objects for the mesh used for a
+ // single energy group, and a number
+ // of other objects and member
+ // functions that we will discuss in
+ // the following sections.
+ //
+ // The main reason for this class is
+ // as follows: for both the forward
+ // problem (with a specified right
+ // hand side) as well as for the
+ // eigenvalue problem, one typically
+ // solves a sequence of problems for
+ // a single energy group each, rather
+ // than the fully coupled
+ // problem. This becomes
+ // understandable once one realizes
+ // that the system matrix for a
+ // single energy group is symmetric
+ // and positive definite (it is
+ // simply a diffusion operator),
+ // whereas the matrix for the fully
+ // coupled problem is generally
+ // nonsymmetric and not definite. It
+ // is also very large and quite full
+ // if more than a few energy groups
+ // are involved.
+ //
+ // Let us first look at the equation
+ // to solve in the case of an
+ // external right hand side (for the time
+ // independent case):
+ // @f{eqnarray*}
+ // -\nabla \cdot(D_g(x) \nabla \phi_g(x))
+ // +
+ // \Sigma_{r,g}(x)\phi_g(x)
+ // =
+ // \chi_g\sum_{g'=1}^G\nu\Sigma_{f,g'}(x)\phi_{g'}(x)
+ // +
+ // \sum_{g'\ne g}\Sigma_{s,g'\to g}(x)\phi_{g'}(x)
+ // +
+ // s_{\mathrm{ext},g}(x)
+ // @f}
+ //
+ // We would typically solve this
+ // equation by moving all the terms
+ // on the right hand side with $g'=g$
+ // to the left hand side, and solving
+ // for $\phi_g$. Of course, we don't
+ // know $\phi_{g'}$ yet, since the
+ // equations for those variables
+ // include right hand side terms
+ // involving $\phi_g$. What one
+ // typically does in such situations
+ // is to iterate: compute
+ // @f{eqnarray*}
+ // -\nabla \cdot(D_g(x) \nabla \phi^{(n)}_g(x))
+ // &+&
+ // \Sigma_{r,g}(x)\phi^{(n)}_g(x)
+ // \\ &=&
+ // \chi_g\sum_{g'=1}^{g-1}\nu\Sigma_{f,g'}(x)\phi^{(n)}_{g'}(x)
+ // +
+ // \chi_g\sum_{g'=g}^G\nu\Sigma_{f,g'}(x)\phi^{(n-1)}_{g'}(x)
+ // +
+ // \sum_{g'\ne g, g'<g}\Sigma_{s,g'\to g}(x)\phi^{(n)}_{g'}(x)
+ // +
+ // \sum_{g'\ne g, g'>g}\Sigma_{s,g'\to g}(x)\phi^{(n-1)}_{g'}(x)
+ // +
+ // s_{\mathrm{ext},g}(x)
+ // @f}
+ //
+ // In other words, we solve the
+ // equation one by one, using values
+ // for $\phi_{g'}$ from the previous
+ // iteration $n-1$ if $g'\ge g$ and
+ // already computed values for
+ // $\phi_{g'}$ from the present
+ // iteration if $g'<g$.
+ //
+ // When computing the eigenvalue, we
+ // do a very similar iteration,
+ // except that we have no external
+ // right hand side and that the
+ // solution is scaled after each
+ // iteration as explained in the
+ // introduction.
+ //
+ // In either case, these two cases
+ // can be treated jointly if all we
+ // do is to equip the following class
+ // with these abilities: (i) form the
+ // left hand side matrix, (ii) form
+ // the in-group right hand side
+ // contribution, i.e. involving the
+ // extraneous source, and (iii) form
+ // that contribution to the right
+ // hand side that stems from group
+ // $g'$. This class does exactly
+ // these tasks (as well as some
+ // book-keeping, such as mesh
+ // refinement, setting up matrices
+ // and vectors, etc). On the other
+ // hand, the class itself has no idea
+ // how many energy groups there are,
+ // and in particular how they
+ // interact, i.e. the decision of how
+ // the outer iteration looks (and
+ // consequently whether we solve an
+ // eigenvalue or a direct problem) is
+ // left to the
+ // NeutronDiffusionProblem class
+ // further down below in this
+ // program.
+ //
+ // So let us go through the class and
+ // its interface:
+template <int dim>
+class EnergyGroup
+{
+ public:
+
+ // @sect5{Public member functions}
+ //
+ // The class has a good number of
+ // public member functions, since
+ // its the way it operates is
+ // controlled from the outside,
+ // and therefore all functions
+ // that do something significant
+ // need to be called from another
+ // class. Let's start off with
+ // book-keeping: the class
+ // obviously needs to know which
+ // energy group it represents,
+ // which material data to use,
+ // and from what coarse grid to
+ // start. The constructor takes
+ // this information and
+ // initializes the relevant
+ // member variables with that
+ // (see below).
+ //
+ // Then we also need functions
+ // that set up the linear system,
+ // i.e. correctly size the matrix
+ // and its sparsity pattern, etc,
+ // given a finite element object
+ // to use. The
+ // <code>setup_linear_system</code>
+ // function does that. Finally,
+ // for this initial block, there
+ // are two functions that return
+ // the number of active cells and
+ // degrees of freedom used in
+ // this object -- using this, we
+ // can make the triangulation and
+ // DoF handler member variables
+ // private, and do not have to
+ // grant external use to it,
+ // enhancing encapsulation:
+ EnergyGroup (const unsigned int group,
+ const MaterialData &material_data,
+ const Triangulation<dim> &coarse_grid,
+ const FiniteElement<dim> &fe);
+
+ void setup_linear_system ();
+
+ unsigned int n_active_cells () const;
+ unsigned int n_dofs () const;
+
+ // Then there are functions that
+ // assemble the linear system for
+ // each iteration and the present
+ // energy group. Note that the
+ // matrix is independent of the
+ // iteration number, so only has
+ // to be computed once for each
+ // refinement cycle. The
+ // situation is a bit more
+ // involved for the right hand
+ // side that has to be updated in
+ // each inverse power iteration,
+ // and that is further
+ // complicated by the fact that
+ // computing it may involve
+ // several different meshes as
+ // explained in the
+ // introduction. To make things
+ // more flexible with regard to
+ // solving the forward or the
+ // eigenvalue problem, we split
+ // the computation of the right
+ // hand side into a function that
+ // assembles the extraneous
+ // source and in-group
+ // contributions (which we will
+ // call with a zero function as
+ // source terms for the
+ // eigenvalue problem) and one
+ // that computes contributions to
+ // the right hand side from
+ // another energy group:
+ void assemble_system_matrix ();
+ void assemble_ingroup_rhs (const Function<dim> &extraneous_source);
+ void assemble_cross_group_rhs (const EnergyGroup<dim> &g_prime);
+
+ // Next we need a set of
+ // functions that actually
+ // compute the solution of a
+ // linear system, and do
+ // something with it (such as
+ // computing the fission source
+ // contribution mentioned in the
+ // introduction, writing
+ // graphical information to an
+ // output file, computing error
+ // indicators, or actually
+ // refining the grid based on
+ // these criteria and thresholds
+ // for refinement and
+ // coarsening). All these
+ // functions will later be called
+ // from the driver class
+ // <code>NeutronDiffusionProblem</code>,
+ // or any other class you may
+ // want to implement to solve a
+ // problem involving the neutron
+ // flux equations:
+ void solve ();
+
+ double get_fission_source () const;
+
+ void output_results (const unsigned int cycle) const;
+
+ void estimate_errors (Vector<float> &error_indicators) const;
+
+ void refine_grid (const Vector<float> &error_indicators,
+ const double refine_threshold,
+ const double coarsen_threshold);
+
+ // @sect5{Public data members}
+ //
+ // As is good practice in object
+ // oriented programming, we hide
+ // most data members by making
+ // them private. However, we have
+ // to grant the class that drives
+ // the process access to the
+ // solution vector as well as the
+ // solution of the previous
+ // iteration, since in the power
+ // iteration, the solution vector
+ // is scaled in every iteration
+ // by the present guess of the
+ // eigenvalue we are looking for:
+ public:
+
+ Vector<double> solution;
+ Vector<double> solution_old;
+
+
+ // @sect5{Private data members}
+ //
+ // The rest of the data members
+ // are private. Compared to all
+ // the previous tutorial
+ // programs, the only new data
+ // members are an integer storing
+ // which energy group this object
+ // represents, and a reference to
+ // the material data object that
+ // this object's constructor gets
+ // passed from the driver
+ // class. Likewise, the
+ // constructor gets a reference
+ // to the finite element object
+ // we are to use.
+ //
+ // Finally, we have to apply
+ // boundary values to the linear
+ // system in each iteration,
+ // i.e. quite frequently. Rather
+ // than interpolating them every
+ // time, we interpolate them once
+ // on each new mesh and then
+ // store them along with all the
+ // other data of this class:
+ private:
+
+ const unsigned int group;
+ const MaterialData &material_data;
+
+ Triangulation<dim> triangulation;
+ const FiniteElement<dim> &fe;
+ DoFHandler<dim> dof_handler;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> system_rhs;
+
+ std::map<unsigned int,double> boundary_values;
+ ConstraintMatrix hanging_node_constraints;
+
+
+ // @sect5{Private member functionss}
+ //
+ // There is one private member
+ // function in this class. It
+ // recursively walks over cells
+ // of two meshes to compute the
+ // cross-group right hand side
+ // terms. The algorithm for this
+ // is explained in the
+ // introduction to this
+ // program. The arguments to this
+ // function are a reference to an
+ // object representing the energy
+ // group against which we want to
+ // integrate a right hand side
+ // term, an iterator to a cell of
+ // the mesh used for the present
+ // energy group, an iterator to a
+ // corresponding cell on the
+ // other mesh, and the matrix
+ // that interpolates the degrees
+ // of freedom from the coarser of
+ // the two cells to the finer
+ // one:
+ private:
+
+ void
+ assemble_cross_group_rhs_recursive (const EnergyGroup<dim> &g_prime,
+ const typename DoFHandler<dim>::cell_iterator &cell_g,
+ const typename DoFHandler<dim>::cell_iterator &cell_g_prime,
+ const FullMatrix<double> prolongation_matrix);
+};
+
+
+ // @sect4{Implementation of the <code>EnergyGroup</code> class}
+
+ // The first few functions of this
+ // class are mostly
+ // self-explanatory. The constructor
+ // only sets a few data members and
+ // creates a copy of the given
+ // triangulation as the base for the
+ // triangulation used for this energy
+ // group. The next two functions
+ // simply return data from private
+ // data members, thereby enabling us
+ // to make these data members
+ // private.
+template <int dim>
+EnergyGroup<dim>::EnergyGroup (const unsigned int group,
+ const MaterialData &material_data,
+ const Triangulation<dim> &coarse_grid,
+ const FiniteElement<dim> &fe)
+ :
+ group (group),
+ material_data (material_data),
+ fe (fe),
+ dof_handler (triangulation)
+{
+ triangulation.copy_triangulation (coarse_grid);
+ dof_handler.distribute_dofs (fe);
+}
+
+
+
+template <int dim>
+unsigned int
+EnergyGroup<dim>::n_active_cells () const
+{
+ return triangulation.n_active_cells ();
+}
+
+
+
+template <int dim>
+unsigned int
+EnergyGroup<dim>::n_dofs () const
+{
+ return dof_handler.n_dofs ();
+}
+
+
+
+ // @sect5{<code>EnergyGroup::setup_linear_system</code>}
+ //
+ // The first "real" function is the
+ // one that sets up the mesh,
+ // matrices, etc, on the new mesh or
+ // after mesh refinement. We use this
+ // function to initialize sparse
+ // system matrices, and the right
+ // hand side vector. If the solution
+ // vector has never been set before
+ // (as indicated by a zero size), we
+ // also initialize it and set it to a
+ // default value. We don't do that if
+ // it already has a non-zero size
+ // (i.e. this function is called
+ // after mesh refinement) since in
+ // that case we want to preserve the
+ // solution across mesh refinement
+ // (something we do in the
+ // <code>EnergyGroup::refine_grid</code>
+ // function).
+template <int dim>
+void
+EnergyGroup<dim>::setup_linear_system ()
+{
+ const unsigned int n_dofs = dof_handler.n_dofs();
+
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+
+ system_matrix.clear ();
+
+ sparsity_pattern.reinit (n_dofs, n_dofs,
+ dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+ hanging_node_constraints.condense (sparsity_pattern);
+ sparsity_pattern.compress ();
+
+ system_matrix.reinit (sparsity_pattern);
+
+ system_rhs.reinit (n_dofs);
+
+ if (solution.size() == 0)
+ {
+ solution.reinit (n_dofs);
+ solution_old.reinit(n_dofs);
+ solution_old = 1.0;
+ solution = solution_old;
+ }
+
+
+ // At the end of this function, we
+ // update the list of boundary
+ // nodes and their values, by first
+ // clearing this list and the
+ // re-interpolating boundary values
+ // (remember that this function is
+ // called after first setting up
+ // the mesh, and each time after
+ // mesh refinement).
+ //
+ // To understand the code, it is
+ // necessary to realize that we
+ // create the mesh using the
+ // <code>GridGenerator::subdivided_hyper_rectangle</code>
+ // function (in
+ // <code>NeutronDiffusionProblem::initialize_problem</code>)
+ // where we set the last parameter
+ // to <code>true</code>. This means that
+ // boundaries of the domain are
+ // "colored", i.e. the four (or
+ // six, in 3d) sides of the domain
+ // are assigned different boundary
+ // indicators. As it turns out, the
+ // bottom boundary gets indicator
+ // zero, the top one boundary
+ // indicator one, and left and
+ // right boundaries get indicators
+ // two and three, respectively.
+ //
+ // In this program, we simulate
+ // only one, namely the top right,
+ // quarter of a reactor. That is,
+ // we want to interpolate boundary
+ // conditions only on the top and
+ // right boundaries, while do
+ // nothing on the bottom and left
+ // boundaries (i.e. impose natural,
+ // no-flux Neumann boundary
+ // conditions). This is most easily
+ // generalized to arbitrary
+ // dimension by saying that we want
+ // to interpolate on those
+ // boundaries with indicators 1, 3,
+ // ..., which we do in the
+ // following loop (note that calls
+ // to
+ // <code>VectorTools::interpolate_boundary_values</code>
+ // are additive, i.e. they do not
+ // first clear the boundary value
+ // map):
+ boundary_values.clear();
+
+ for (unsigned int i=0; i<dim; ++i)
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 2*i+1,
+ ZeroFunction<dim>(),
+ boundary_values);
+}
+
+
+
+ // @sect5{<code>EnergyGroup::assemble_system_matrix</code>}
+ //
+ // Next we need functions assembling
+ // the system matrix and right hand
+ // sides. Assembling the matrix is
+ // straightforward given the
+ // equations outlined in the
+ // introduction as well as what we've
+ // seen in previous example
+ // programs. Note the use of
+ // <code>cell->material_id()</code> to get at
+ // the kind of material from which a
+ // cell is made up of. Note also how
+ // we set the order of the quadrature
+ // formula so that it is always
+ // appropriate for the finite element
+ // in use.
+ //
+ // Finally, note that since we only
+ // assemble the system matrix here,
+ // we can't yet eliminate boundary
+ // values (we need the right hand
+ // side vector for this). We defer
+ // this to the <code>EnergyGroup::solve</code>
+ // function, at which point all the
+ // information is available.
+template <int dim>
+void
+EnergyGroup<dim>::assemble_system_matrix ()
+{
+ const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+ FEValues<dim> 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.n_quadrature_points;
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ for (; cell!=endc; ++cell)
+ {
+ cell_matrix = 0;
+
+ fe_values.reinit (cell);
+
+ const double diffusion_coefficient
+ = material_data.get_diffusion_coefficient (group, cell->material_id());
+ const double removal_XS
+ = material_data.get_removal_XS (group,cell->material_id());
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ 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)
+ +
+ removal_XS *
+ 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);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+ }
+
+ hanging_node_constraints.condense (system_matrix);
+}
+
+
+
+ // @sect5{<code>EnergyGroup::assemble_ingroup_rhs</code>}
+ //
+ // As explained in the documentation
+ // of the <code>EnergyGroup</code> class, we
+ // split assembling the right hand
+ // side into two parts: the ingroup
+ // and the cross-group
+ // couplings. First, we need a
+ // function to assemble the right
+ // hand side of one specific group
+ // here, i.e. including an extraneous
+ // source (that we will set to zero
+ // for the eigenvalue problem) as
+ // well as the ingroup fission
+ // contributions. (In-group
+ // scattering has already been
+ // accounted for with the definition
+ // of removal cross section.) The
+ // function's workings are pretty
+ // standard as far as assembling
+ // right hand sides go, and therefore
+ // does not require more comments
+ // except that we mention that the
+ // right hand side vector is set to
+ // zero at the beginning of the
+ // function -- something we are not
+ // going to do for the cross-group
+ // terms that simply add to the right
+ // hand side vector.
+template <int dim>
+void EnergyGroup<dim>::assemble_ingroup_rhs (const Function<dim> &extraneous_source)
+{
+ system_rhs.reinit (dof_handler.n_dofs());
+
+ const QGauss<dim> quadrature_formula (fe.degree + 1);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_q_points |
+ update_JxW_values);
+
+ Vector<double> cell_rhs (dofs_per_cell);
+ std::vector<double> extraneous_source_values (n_q_points);
+ std::vector<double> solution_old_values (n_q_points);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ for (; cell!=endc; ++cell)
+ {
+ cell_rhs = 0;
+
+ fe_values.reinit (cell);
+
+ const double fission_dist_XS
+ = material_data.get_fission_dist_XS (group, group, cell->material_id());
+
+ extraneous_source.value_list (fe_values.get_quadrature_points(),
+ extraneous_source_values);
+
+ fe_values.get_function_values (solution_old, solution_old_values);
+
+ cell->get_dof_indices (local_dof_indices);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_rhs(i) += ((extraneous_source_values[q_point]
+ +
+ fission_dist_XS *
+ solution_old_values[q_point]) *
+ fe_values.shape_value(i,q_point) *
+ fe_values.JxW(q_point));
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+}
+
+
+
+ // @sect5{<code>EnergyGroup::assemble_cross_group_rhs</code>}
+ //
+ // The more interesting function for
+ // assembling the right hand side
+ // vector for the equation of a
+ // single energy group is the one
+ // that couples energy group $g$ and
+ // $g'$. As explained in the
+ // introduction, we first have to
+ // find the set of cells common to
+ // the meshes of the two energy
+ // groups. First we call
+ // <code>get_finest_common_cells</code> to
+ // obtain this list of pairs of
+ // common cells from both
+ // meshes. Both cells in a pair may
+ // not be active but at least one of
+ // them is. We then hand each of
+ // these cell pairs off to a function
+ // tha computes the right hand side
+ // terms recursively.
+ //
+ // Note that ingroup coupling is
+ // handled already before, so we exit
+ // the function early if $g=g'$.
+template <int dim>
+void EnergyGroup<dim>::assemble_cross_group_rhs (const EnergyGroup<dim> &g_prime)
+{
+ if (group == g_prime.group)
+ return;
+
+ const std::list<std::pair<typename DoFHandler<dim>::cell_iterator,
+ typename DoFHandler<dim>::cell_iterator> >
+ cell_list
+ = GridTools::get_finest_common_cells (dof_handler,
+ g_prime.dof_handler);
+
+ typename std::list<std::pair<typename DoFHandler<dim>::cell_iterator,
+ typename DoFHandler<dim>::cell_iterator> >
+ ::const_iterator
+ cell_iter = cell_list.begin();
+
+ for (; cell_iter!=cell_list.end(); ++cell_iter)
+ {
+ FullMatrix<double> unit_matrix (fe.dofs_per_cell);
+ for (unsigned int i=0; i<unit_matrix.m(); ++i)
+ unit_matrix(i,i) = 1;
+ assemble_cross_group_rhs_recursive (g_prime,
+ cell_iter->first,
+ cell_iter->second,
+ unit_matrix);
+ }
+}
+
+
+
+ // @sect5{<code>EnergyGroup::assemble_cross_group_rhs_recursive</code>}
+ //
+ // This is finally the function that
+ // handles assembling right hand side
+ // terms on potentially different
+ // meshes recursively, using the
+ // algorithm described in the
+ // introduction. The function takes a
+ // reference to the object
+ // representing energy group $g'$, as
+ // well as iterators to corresponding
+ // cells in the meshes for energy
+ // groups $g$ and $g'$. At first,
+ // i.e. when this function is called
+ // from the one above, these two
+ // cells will be matching cells on
+ // two meshes; however, one of the
+ // two may be further refined, and we
+ // will call the function recursively
+ // with one of the two iterators
+ // replaced by one of the children of
+ // the original cell.
+ //
+ // The last argument is the matrix
+ // product matrix $B_{c^{(k)}}^T
+ // \cdots B_{c'}^T B_c^T$ from the
+ // introduction that interpolates
+ // from the coarser of the two cells
+ // to the finer one. If the two cells
+ // match, then this is the identity
+ // matrix -- exactly what we pass to
+ // this function initially.
+ //
+ // The function has to consider two
+ // cases: that both of the two cells
+ // are not further refined, i.e. have
+ // no children, in which case we can
+ // finally assemble the right hand
+ // side contributions of this pair of
+ // cells; and that one of the two
+ // cells is further refined, in which
+ // case we have to keep recursing by
+ // looping over the children of the
+ // one cell that is not active. These
+ // two cases will be discussed below:
+template <int dim>
+void
+EnergyGroup<dim>::
+assemble_cross_group_rhs_recursive (const EnergyGroup<dim> &g_prime,
+ const typename DoFHandler<dim>::cell_iterator &cell_g,
+ const typename DoFHandler<dim>::cell_iterator &cell_g_prime,
+ const FullMatrix<double> prolongation_matrix)
+{
+ // The first case is that both
+ // cells are no further refined. In
+ // that case, we can assemble the
+ // relevant terms (see the
+ // introduction). This involves
+ // assembling the mass matrix on
+ // the finer of the two cells (in
+ // fact there are two mass matrices
+ // with different coefficients, one
+ // for the fission distribution
+ // cross section
+ // $\chi_g\nu\Sigma_{f,g'}$ and one
+ // for the scattering cross section
+ // $\Sigma_{s,g'\to g}$). This is
+ // straight forward, but note how
+ // we determine which of the two
+ // cells is ther finer one by
+ // looking at the refinement level
+ // of the two cells:
+ if (!cell_g->has_children() && !cell_g_prime->has_children())
+ {
+ const QGauss<dim> quadrature_formula (fe.degree+1);
+ const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_JxW_values);
+
+ if (cell_g->level() > cell_g_prime->level())
+ fe_values.reinit (cell_g);
+ else
+ fe_values.reinit (cell_g_prime);
+
+ const double fission_dist_XS
+ = material_data.get_fission_dist_XS (group, g_prime.group,
+ cell_g_prime->material_id());
+
+ const double scattering_XS
+ = material_data.get_scattering_XS (g_prime.group, group,
+ cell_g_prime->material_id());
+
+ FullMatrix<double> local_mass_matrix_f (fe.dofs_per_cell,
+ fe.dofs_per_cell);
+ FullMatrix<double> local_mass_matrix_g (fe.dofs_per_cell,
+ fe.dofs_per_cell);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ {
+ local_mass_matrix_f(i,j) += (fission_dist_XS *
+ fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point) *
+ fe_values.JxW(q_point));
+ local_mass_matrix_g(i,j) += (scattering_XS *
+ fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point) *
+ fe_values.JxW(q_point));
+ }
+
+ // Now we have all the
+ // interpolation (prolongation)
+ // matrices as well as local
+ // mass matrices, so we only
+ // have to form the product
+ // @f[
+ // F_i|_{K_{cc'\cdots
+ // c^{(k)}}} = [B_c B_{c'}
+ // \cdots B_{c^{(k)}}
+ // M_{K_{cc'\cdots
+ // c^{(k)}}}]^{ij}
+ // \phi_{g'}^j,
+ // @f]
+ // or
+ // @f[
+ // F_i|_{K_{cc'\cdots
+ // c^{(k)}}} = [(B_c B_{c'}
+ // \cdots B_{c^{(k)}}
+ // M_{K_{cc'\cdots
+ // c^{(k)}}})^T]^{ij}
+ // \phi_{g'}^j,
+ // @f]
+ // depending on which of the two
+ // cells is the finer. We do this
+ // using either the matrix-vector
+ // product provided by the <code>vmult</code>
+ // function, or the product with the
+ // transpose matrix using <code>Tvmult</code>.
+ // After doing so, we transfer the
+ // result into the global right hand
+ // side vector of energy group $g$.
+ Vector<double> g_prime_new_values (fe.dofs_per_cell);
+ Vector<double> g_prime_old_values (fe.dofs_per_cell);
+ cell_g_prime->get_dof_values (g_prime.solution_old, g_prime_old_values);
+ cell_g_prime->get_dof_values (g_prime.solution, g_prime_new_values);
+
+ Vector<double> cell_rhs (fe.dofs_per_cell);
+ Vector<double> tmp (fe.dofs_per_cell);
+
+ if (cell_g->level() > cell_g_prime->level())
+ {
+ prolongation_matrix.vmult (tmp, g_prime_old_values);
+ local_mass_matrix_f.vmult (cell_rhs, tmp);
+
+ prolongation_matrix.vmult (tmp, g_prime_new_values);
+ local_mass_matrix_g.vmult_add (cell_rhs, tmp);
+ }
+ else
+ {
+ local_mass_matrix_f.vmult (tmp, g_prime_old_values);
+ prolongation_matrix.Tvmult (cell_rhs, tmp);
+
+ local_mass_matrix_g.vmult (tmp, g_prime_new_values);
+ prolongation_matrix.Tvmult_add (cell_rhs, tmp);
+ }
+
+ std::vector<unsigned int> local_dof_indices (fe.dofs_per_cell);
+ cell_g->get_dof_indices (local_dof_indices);
+
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+
+ // The alternative is that one of
+ // the two cells is further
+ // refined. In that case, we have
+ // to loop over all the children,
+ // multiply the existing
+ // interpolation (prolongation)
+ // product of matrices from the
+ // left with the interpolation from
+ // the present cell to its child
+ // (using the matrix-matrix
+ // multiplication function
+ // <code>mmult</code>), and then hand the
+ // result off to this very same
+ // function again, but with the
+ // cell that has children replaced
+ // by one of its children:
+ else
+ for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell;++child)
+ {
+ FullMatrix<double> new_matrix (fe.dofs_per_cell, fe.dofs_per_cell);
+ fe.get_prolongation_matrix(child).mmult (new_matrix,
+ prolongation_matrix);
+
+ if (cell_g->has_children())
+ assemble_cross_group_rhs_recursive (g_prime,
+ cell_g->child(child), cell_g_prime,
+ new_matrix);
+ else
+ assemble_cross_group_rhs_recursive (g_prime,
+ cell_g, cell_g_prime->child(child),
+ new_matrix);
+ }
+}
+
+
+ // @sect5{<code>EnergyGroup::get_fission_source</code>}
+ //
+ // In the (inverse) power iteration,
+ // we use the integrated fission
+ // source to update the
+ // $k$-eigenvalue. Given its
+ // definition, the following function
+ // is essentially self-explanatory:
+template <int dim>
+double EnergyGroup<dim>::get_fission_source () const
+{
+ const QGauss<dim> quadrature_formula (fe.degree + 1);
+ const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_JxW_values);
+
+ std::vector<double> solution_values (n_q_points);
+
+ double fission_source = 0;
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ fe_values.reinit (cell);
+
+ const double fission_XS
+ = material_data.get_fission_XS(group, cell->material_id());
+
+ fe_values.get_function_values (solution, solution_values);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ fission_source += (fission_XS *
+ solution_values[q_point] *
+ fe_values.JxW(q_point));
+ }
+
+ return fission_source;
+}
+
+
+ // @sect5{<code>EnergyGroup::solve</code>}
+ //
+ // Next a function that solves the
+ // linear system assembled
+ // before. Things are pretty much
+ // standard, except that we delayed
+ // applying boundary values until we
+ // get here, since in all the
+ // previous functions we were still
+ // adding up contributions the right
+ // hand side vector.
+template <int dim>
+void
+EnergyGroup<dim>::solve ()
+{
+ hanging_node_constraints.condense (system_rhs);
+ MatrixTools::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+
+ SolverControl solver_control (system_matrix.m(),
+ 1e-12*system_rhs.l2_norm());
+ SolverCG<> cg (solver_control);
+
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.2);
+
+ cg.solve (system_matrix, solution, system_rhs, preconditioner);
+
+ hanging_node_constraints.distribute (solution);
+}
+
+
+
+ // @sect5{<code>EnergyGroup::estimate_errors</code>}
+ //
+ // Mesh refinement is split into two
+ // functions. The first estimates the
+ // error for each cell, normalizes it
+ // by the magnitude of the solution,
+ // and returns it in the vector given
+ // as an argument. The calling
+ // function collects all error
+ // indicators from all energy groups,
+ // and computes thresholds for
+ // refining and coarsening cells.
+template <int dim>
+void EnergyGroup<dim>::estimate_errors (Vector<float> &error_indicators) const
+{
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ QGauss<dim-1> (fe.degree + 1),
+ typename FunctionMap<dim>::type(),
+ solution,
+ error_indicators);
+ error_indicators /= solution.linfty_norm();
+}
+
+
+
+ // @sect5{<code>EnergyGroup::refine_grid</code>}
+ //
+ // The second part is to refine the
+ // grid given the error indicators
+ // compute in the previous function
+ // and error thresholds above which
+ // cells shall be refined or below
+ // which cells shall be
+ // coarsened. Note that we do not use
+ // any of the functions in
+ // <code>GridRefinement</code> here,
+ // but rather set refinement flags
+ // ourselves.
+ //
+ // After setting these flags, we use
+ // the SolutionTransfer class to move
+ // the solution vector from the old
+ // to the new mesh. The procedure
+ // used here is described in detail
+ // in the documentation of that
+ // class:
+template <int dim>
+void EnergyGroup<dim>::refine_grid (const Vector<float> &error_indicators,
+ const double refine_threshold,
+ const double coarsen_threshold)
+{
+ typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active(),
+ endc = triangulation.end();
+
+ for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
+ if (error_indicators(cell_index) > refine_threshold)
+ cell->set_refine_flag ();
+ else if (error_indicators(cell_index) < coarsen_threshold)
+ cell->set_coarsen_flag ();
+
+ SolutionTransfer<dim, double> soltrans(dof_handler);
+
+ triangulation.prepare_coarsening_and_refinement();
+ soltrans.prepare_for_coarsening_and_refinement(solution);
+
+ triangulation.execute_coarsening_and_refinement ();
+ dof_handler.distribute_dofs (fe);
+
+ solution.reinit (dof_handler.n_dofs());
+ soltrans.interpolate(solution_old, solution);
+
+ solution_old.reinit (dof_handler.n_dofs());
+ solution_old = solution;
+}
+
+
+ // @sect5{<code>EnergyGroup::output_results</code>}
+ //
+ // The last function of this class
+ // outputs meshes and solutions after
+ // each mesh iteration. This has been
+ // shown many times before. The only
+ // thing worth pointing out is the
+ // use of the
+ // <code>Utilities::int_to_string</code>
+ // function to convert an integer
+ // into its string
+ // representation. The second
+ // argument of that function denotes
+ // how many digits we shall use -- if
+ // this value was larger than one,
+ // then the number would be padded by
+ // leading zeros.
+template <int dim>
+void
+EnergyGroup<dim>::output_results (const unsigned int cycle) const
+{
+ {
+ const std::string filename = std::string("grid-") +
+ Utilities::int_to_string(group,1) +
+ "." +
+ Utilities::int_to_string(cycle,1) +
+ ".eps";
+ std::ofstream output (filename.c_str());
+
+ GridOut grid_out;
+ grid_out.write_eps (triangulation, output);
+ }
+
+ {
+ const std::string filename = std::string("solution-") +
+ Utilities::int_to_string(group,1) +
+ "." +
+ Utilities::int_to_string(cycle,1) +
+ ".gmv";
+
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+ data_out.build_patches ();
+
+ std::ofstream output (filename.c_str());
+ data_out.write_gmv (output);
+ }
+}
+
+
+
+ // @sect3{The <code>NeutronDiffusionProblem</code> class template}
+
+ // This is the main class of the
+ // program, not because it implements
+ // all the functionality (in fact,
+ // most of it is implemented in the
+ // <code>EnergyGroup</code> class)
+ // but because it contains the
+ // driving algorithm that determines
+ // what to compute and when. It is
+ // mostly as shown in many of the
+ // other tutorial programs in that it
+ // has a public <code>run</code>
+ // function and private functions
+ // doing all the rest. In several
+ // places, we have to do something
+ // for all energy groups, in which
+ // case we will spawn threads for
+ // each group to let these things run
+ // in parallel if deal.II was
+ // configured for multithreading.
+ //
+ // The biggest difference to previous
+ // example programs is that we also
+ // declare a nested class that has
+ // member variables for all the
+ // run-time parameters that can be
+ // passed to the program in an input
+ // file. Right now, these are the
+ // number of energy groups, the
+ // number of refinement cycles, the
+ // polynomial degree of the finite
+ // element to be used, and the
+ // tolerance used to determine when
+ // convergence of the inverse power
+ // iteration has occurred. In
+ // addition, we have a constructor of
+ // this class that sets all these
+ // values to their default values, a
+ // function
+ // <code>declare_parameters</code>
+ // that described to the
+ // ParameterHandler class already
+ // used in @ref step_19 "step-19"
+ // what parameters are accepted in
+ // the input file, and a function
+ // <code>get_parameters</code> that
+ // can extract the values of these
+ // parameters from a ParameterHandler
+ // object.
+template <int dim>
+class NeutronDiffusionProblem
+{
+ public:
+ class Parameters
+ {
+ public:
+ Parameters ();
+
+ static void declare_parameters (ParameterHandler &prm);
+ void get_parameters (ParameterHandler &prm);
+
+ unsigned int n_groups;
+ unsigned int n_refinement_cycles;
+
+ unsigned int fe_degree;
+
+ double convergence_tolerance;
+ };
+
+
+
+ NeutronDiffusionProblem (const Parameters ¶meters);
+ ~NeutronDiffusionProblem ();
+
+ void run ();
+
+ private:
+ // @sect5{Private member functions}
+
+ // There are not that many member
+ // functions in this class since
+ // most of the functionality has
+ // been moved into the
+ // <code>EnergyGroup</code> class
+ // and is simply called from the
+ // <code>run()</code> member
+ // function of this class. The
+ // ones that remain have
+ // self-explanatory names:
+ void initialize_problem();
+
+ void refine_grid ();
+
+ double get_total_fission_source () const;
+
+
+ // @sect5{Private member variables}
+
+ // Next, we have a few member
+ // variables. In particular,
+ // these are (i) a reference to
+ // the parameter object (owned by
+ // the main function of this
+ // program, and passed to the
+ // constructor of this class),
+ // (ii) an object describing the
+ // material parameters for the
+ // number of energy groups
+ // requested in the input file,
+ // and (iii) the finite element
+ // to be used by all energy
+ // groups:
+ const Parameters ¶meters;
+ const MaterialData material_data;
+ FE_Q<dim> fe;
+
+ // Furthermore, we have (iv) the
+ // value of the computed
+ // eigenvalue at the present
+ // iteration. This is, in fact,
+ // the only part of the solution
+ // that is shared between all
+ // energy groups -- all other
+ // parts of the solution, such as
+ // neutron fluxes are particular
+ // to one or the other energy
+ // group, and are therefore
+ // stored in objects that
+ // describe a single energy
+ // group:
+ double k_eff;
+
+ // Finally, (v), we have an array
+ // of pointers to the energy
+ // group objects. The length of
+ // this array is, of course,
+ // equal to the number of energy
+ // groups specified in the
+ // parameter file.
+ std::vector<EnergyGroup<dim>*> energy_groups;
+};
+
+
+ // @sect4{Implementation of the <code>NeutronDiffusionProblem::Parameters</code> class}
+
+ // Before going on to the
+ // implementation of the outer class,
+ // we have to implement the functions
+ // of the parameters structure. This
+ // is pretty straightforward and, in
+ // fact, looks pretty much the same
+ // for all such parameters classes
+ // using the ParameterHandler
+ // capabilities. We will therefore
+ // not comment further on this:
+template <int dim>
+NeutronDiffusionProblem<dim>::Parameters::Parameters ()
+ :
+ n_groups (2),
+ n_refinement_cycles (5),
+ fe_degree (2),
+ convergence_tolerance (1e-12)
+{}
+
+
+
+template <int dim>
+void
+NeutronDiffusionProblem<dim>::Parameters::
+declare_parameters (ParameterHandler &prm)
+{
+ prm.declare_entry ("Number of energy groups", "2",
+ Patterns::Integer (),
+ "The number of energy different groups considered");
+ prm.declare_entry ("Refinement cycles", "5",
+ Patterns::Integer (),
+ "Number of refinement cycles to be performed");
+ prm.declare_entry ("Finite element degree", "2",
+ Patterns::Integer (),
+ "Polynomial degree of the finite element to be used");
+ prm.declare_entry ("Power iteration tolerance", "1e-12",
+ Patterns::Double (),
+ "Inner power iterations are stopped when the change in k_eff falls "
+ "below this tolerance");
+}
+
+
+
+template <int dim>
+void
+NeutronDiffusionProblem<dim>::Parameters::
+get_parameters (ParameterHandler &prm)
+{
+ n_groups = prm.get_integer ("Number of energy groups");
+ n_refinement_cycles = prm.get_integer ("Refinement cycles");
+ fe_degree = prm.get_integer ("Finite element degree");
+ convergence_tolerance = prm.get_double ("Power iteration tolerance");
+}
+
+
+
+
+ // @sect4{Implementation of the <code>NeutronDiffusionProblem</code> class}
+
+ // Now for the
+ // <code>NeutronDiffusionProblem</code>
+ // class. The constructor and
+ // destructor have nothing of much
+ // interest:
+template <int dim>
+NeutronDiffusionProblem<dim>::
+NeutronDiffusionProblem (const Parameters ¶meters)
+ :
+ parameters (parameters),
+ material_data (parameters.n_groups),
+ fe (parameters.fe_degree)
+{}
+
+
+
+template <int dim>
+NeutronDiffusionProblem<dim>::~NeutronDiffusionProblem ()
+{
+ for (unsigned int group=0; group<energy_groups.size(); ++group)
+ delete energy_groups[group];
+
+ energy_groups.resize (0);
+}
+
+ // @sect5{<code>NeutronDiffusionProblem::initialize_problem</code>}
+ //
+ // The first function of interest is
+ // the one that sets up the geometry
+ // of the reactor core. This is
+ // described in more detail in the
+ // introduction.
+ //
+ // The first part of the function
+ // defines geometry data, and then
+ // creates a coarse mesh that has as
+ // many cells as there are fuel rods
+ // (or pin cells, for that matter) in
+ // that part of the reactor core that
+ // we simulate. As mentioned when
+ // interpolating boundary values
+ // above, the last parameter to the
+ // <code>GridGenerator::subdivided_hyper_rectangle</code>
+ // function specifies that sides of
+ // the domain shall have unique
+ // boundary indicators that will
+ // later allow us to determine in a
+ // simple way which of the boundaries
+ // have Neumann and which have
+ // Dirichlet conditions attached to
+ // them.
+template <int dim>
+void NeutronDiffusionProblem<dim>::initialize_problem()
+{
+ const unsigned int rods_per_assembly_x = 17,
+ rods_per_assembly_y = 17;
+ const double pin_pitch_x = 1.26,
+ pin_pitch_y = 1.26;
+ const double assembly_height = 200;
+
+ const unsigned int assemblies_x = 2,
+ assemblies_y = 2,
+ assemblies_z = 1;
+
+ const Point<dim> bottom_left = Point<dim>();
+ const Point<dim> upper_right = (dim == 2
+ ?
+ Point<dim> (assemblies_x*rods_per_assembly_x*pin_pitch_x,
+ assemblies_y*rods_per_assembly_y*pin_pitch_y)
+ :
+ Point<dim> (assemblies_x*rods_per_assembly_x*pin_pitch_x,
+ assemblies_y*rods_per_assembly_y*pin_pitch_y,
+ assemblies_z*assembly_height));
+
+ std::vector<unsigned int> n_subdivisions;
+ n_subdivisions.push_back (assemblies_x*rods_per_assembly_x);
+ if (dim >= 2)
+ n_subdivisions.push_back (assemblies_y*rods_per_assembly_y);
+ if (dim >= 3)
+ n_subdivisions.push_back (assemblies_z);
+
+ Triangulation<dim> coarse_grid;
+ GridGenerator::subdivided_hyper_rectangle (coarse_grid,
+ n_subdivisions,
+ bottom_left,
+ upper_right,
+ true);
+
+
+ // The second part of the function
+ // deals with material numbers of
+ // pin cells of each type of
+ // assembly. Here, we define four
+ // different types of assembly, for
+ // which we describe the
+ // arrangement of fuel rods in the
+ // following tables.
+ //
+ // The assemblies described here
+ // are taken from the benchmark
+ // mentioned in the introduction
+ // and are (in this order):
+ // <ol>
+ // <li>'UX' Assembly: UO2 fuel assembly
+ // with 24 guide tubes and a central
+ // Moveable Fission Chamber
+ // <li>'UA' Assembly: UO2 fuel assembly
+ // with 24 AIC and a central
+ // Moveable Fission Chamber
+ // <li>'PX' Assembly: MOX fuel assembly
+ // with 24 guide tubes and a central
+ // Moveable Fission Chamber
+ // <li>'R' Assembly: a reflector.
+ // </ol>
+ //
+ // Note that the numbers listed
+ // here and taken from the
+ // benchmark description are, in
+ // good old Fortran fashion,
+ // one-based. We will later
+ // subtract one from each number
+ // when assigning materials to
+ // individual cells to convert
+ // things into the C-style
+ // zero-based indexing.
+ const unsigned int n_assemblies=4;
+ const unsigned int
+ assembly_materials[n_assemblies][rods_per_assembly_x][rods_per_assembly_y]
+ = {
+ {
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 5, 1, 1, 5, 1, 1, 7, 1, 1, 5, 1, 1, 5, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }
+ },
+ {
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 8, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 8, 1, 1, 8, 1, 1, 7, 1, 1, 8, 1, 1, 8, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 8, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 8, 1, 1, 8, 1, 1, 8, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+ { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }
+ },
+ {
+ { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 },
+ { 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2 },
+ { 2, 3, 3, 3, 3, 5, 3, 3, 5, 3, 3, 5, 3, 3, 3, 3, 2 },
+ { 2, 3, 3, 5, 3, 4, 4, 4, 4, 4, 4, 4, 3, 5, 3, 3, 2 },
+ { 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2 },
+ { 2, 3, 5, 4, 4, 5, 4, 4, 5, 4, 4, 5, 4, 4, 5, 3, 2 },
+ { 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 },
+ { 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 },
+ { 2, 3, 5, 4, 4, 5, 4, 4, 7, 4, 4, 5, 4, 4, 5, 3, 2 },
+ { 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 },
+ { 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 },
+ { 2, 3, 5, 4, 4, 5, 4, 4, 5, 4, 4, 5, 4, 4, 5, 3, 2 },
+ { 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2 },
+ { 2, 3, 3, 5, 3, 4, 4, 4, 4, 4, 4, 4, 3, 5, 3, 3, 2 },
+ { 2, 3, 3, 3, 3, 5, 3, 3, 5, 3, 3, 5, 3, 3, 3, 3, 2 },
+ { 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2 },
+ { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 }
+ },
+ {
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
+ { 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 }
+ }
+ };
+
+ // After the description of the
+ // materials that make up an
+ // assembly, we have to specify the
+ // arrangement of assemblies within
+ // the core. We use a symmetric
+ // pattern that in fact only uses
+ // the 'UX' and 'PX' assemblies:
+ const unsigned int core[assemblies_x][assemblies_y][assemblies_z]
+ = {{{0}, {2}}, {{2}, {0}}};
+
+ // We are now in a position to
+ // actually set material IDs for
+ // each cell. To this end, we loop
+ // over all cells, look at the
+ // location of the cell's center,
+ // and determine which assembly and
+ // fuel rod this would be in. (We
+ // add a few checks to see that the
+ // locations we compute are within
+ // the bounds of the arrays in
+ // which we have to look up
+ // materials.) At the end of the
+ // loop, we set material
+ // identifiers accordingly:
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = coarse_grid.begin_active();
+ cell!=coarse_grid.end();
+ ++cell)
+ {
+ const Point<dim> cell_center = cell->center();
+
+ const unsigned int tmp_x = int(cell_center[0]/pin_pitch_x);
+ const unsigned int ax = tmp_x/rods_per_assembly_x;
+ const unsigned int cx = tmp_x - ax * rods_per_assembly_x;
+
+ const unsigned tmp_y = int(cell_center[1]/pin_pitch_y);
+ const unsigned int ay = tmp_y/rods_per_assembly_y;
+ const unsigned int cy = tmp_y - ay * rods_per_assembly_y;
+
+ const unsigned int az = (dim == 2
+ ?
+ 0
+ :
+ int (cell_center[dim-1]/assembly_height));
+
+ Assert (ax < assemblies_x, ExcInternalError());
+ Assert (ay < assemblies_y, ExcInternalError());
+ Assert (az < assemblies_z, ExcInternalError());
+
+ Assert (core[ax][ay][az] < n_assemblies, ExcInternalError());
+
+ Assert (cx < rods_per_assembly_x, ExcInternalError());
+ Assert (cy < rods_per_assembly_y, ExcInternalError());
+
+ cell->set_material_id(assembly_materials[core[ax][ay][az]][cx][cy] - 1);
+ }
+
+ // With the coarse mesh so
+ // initialized, we create the
+ // appropriate number of energy
+ // group objects and let them
+ // initialize their individual
+ // meshes with the coarse mesh
+ // generated above:
+ energy_groups.resize (parameters.n_groups);
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ energy_groups[group] = new EnergyGroup<dim> (group, material_data,
+ coarse_grid, fe);
+}
+
+
+ // @sect5{<code>NeutronDiffusionProblem::get_total_fission_source</code>}
+ //
+ // In the eigenvalue computation, we
+ // need to calculate total fission
+ // neutron source after each power
+ // iteration. The total power then is
+ // used to renew k-effective.
+ //
+ // Since the total fission source is
+ // a sum over all the energy groups,
+ // and since each of these sums can
+ // be computed independently, we
+ // actually do this in parallel. One
+ // of the problems is that the
+ // function in the
+ // <code>EnergyGroup</code> class
+ // that computes the fission source
+ // returns a value. If we now simply
+ // spin off a new thread, we have to
+ // later capture the return value of
+ // the function run on that
+ // thread. The way this can be done
+ // is to use the return value of the
+ // Threads::spawn function, which is
+ // of type Threads::Thread@<double@>
+ // if the function spawned returns a
+ // double. We can the later ask this
+ // object for the returned value
+ // (when doing so, the
+ // Threads::Thread@<double@>::return_value
+ // function first waits for the
+ // thread to finish).
+ //
+ // The way this function then works
+ // is to first spawn one thread for
+ // each energy group we work with,
+ // then one-by-one collecting the
+ // returned values of each thread and
+ // return the sum.
+template <int dim>
+double NeutronDiffusionProblem<dim>::get_total_fission_source () const
+{
+ std::vector<Threads::Thread<double> > threads;
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ threads.push_back (Threads::spawn (*energy_groups[group],
+ &EnergyGroup<dim>::get_fission_source) ());
+
+ double fission_source = 0;
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ fission_source += threads[group].return_value ();
+
+ return fission_source;
+}
+
+
+
+
+ // @sect5{<code>NeutronDiffusionProblem::refine_grid</code>}
+ //
+ // The next function lets the
+ // individual energy group objects
+ // refine their meshes. Much of this,
+ // again, is a task that can be done
+ // independently in parallel: first,
+ // let all the energy group objects
+ // calculate their error indicators
+ // in parallel, then compute the
+ // maximum error indicator over all
+ // energy groups and determine
+ // thresholds for refinement and
+ // coarsening of cells, and then ask
+ // all the energy groups to refine
+ // their meshes accordingly, again in
+ // parallel.
+template <int dim>
+void NeutronDiffusionProblem<dim>::refine_grid ()
+{
+ std::vector<unsigned int> n_cells (parameters.n_groups);
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ n_cells[group] = energy_groups[group]->n_active_cells();
+
+ BlockVector<float> group_error_indicators(n_cells);
+
+ {
+ Threads::ThreadGroup<> threads;
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ threads += Threads::spawn (*energy_groups[group], &EnergyGroup<dim>::estimate_errors)
+ (group_error_indicators.block(group));
+ threads.join_all ();
+ }
+
+ const float max_error = group_error_indicators.linfty_norm();
+ const float refine_threshold = 0.3*max_error;
+ const float coarsen_threshold = 0.01*max_error;
+
+ {
+ Threads::ThreadGroup<> threads;
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ threads += Threads::spawn (*energy_groups[group], &EnergyGroup<dim>::refine_grid)
+ (group_error_indicators.block(group),
+ refine_threshold,
+ coarsen_threshold);
+ threads.join_all ();
+ }
+}
+
+
+ // @sect5{<code>NeutronDiffusionProblem::run</code>}
+ //
+ // Finally, this is the function
+ // where the meat is: iterate on a
+ // sequence of meshes, and on each of
+ // them do a power iteration to
+ // compute the eigenvalue.
+ //
+ // Given the description of the
+ // algorithm in the introduction,
+ // there is actually not much to
+ // comment on:
+template <int dim>
+void NeutronDiffusionProblem<dim>::run ()
+{
+ std::cout << std::setprecision (12) << std::fixed;
+
+ double k_eff_old = k_eff;
+
+ Timer timer;
+ timer.start ();
+
+ for (unsigned int cycle=0; cycle<parameters.n_refinement_cycles; ++cycle)
+ {
+ std::cout << "Cycle " << cycle << ':' << std::endl;
+
+ if (cycle == 0)
+ initialize_problem();
+ else
+ {
+ refine_grid ();
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ energy_groups[group]->solution *= k_eff;
+ }
+
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ energy_groups[group]->setup_linear_system ();
+
+ std::cout << " Numbers of active cells: ";
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ std::cout << energy_groups[group]->n_active_cells()
+ << ' ';
+ std::cout << std::endl;
+ std::cout << " Numbers of degrees of freedom: ";
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ std::cout << energy_groups[group]->n_dofs()
+ << ' ';
+ std::cout << std::endl << std::endl;
+
+
+ Threads::ThreadGroup<> threads;
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ threads += Threads::spawn
+ (*energy_groups[group], &EnergyGroup<dim>::assemble_system_matrix)
+ ();
+ threads.join_all ();
+
+ double max_old = 0;
+
+ // indicate this is a eigenvalue problem
+ unsigned int isour = 0;
+ // store relative error between two
+ // successive power iterations
+ double error;
+
+ unsigned int iteration = 1;
+ do
+ {
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ {
+ energy_groups[group]->assemble_ingroup_rhs (ZeroFunction<dim>());
+
+ for (unsigned int bgroup=0; bgroup<parameters.n_groups; ++bgroup)
+ energy_groups[group]->assemble_cross_group_rhs (*energy_groups[bgroup]);
+
+ energy_groups[group]->solve ();
+ }
+
+ if (isour==1)
+ {
+ double max_current = 0;
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ {
+ max_current = std::max (max_current,
+ energy_groups[group]->solution.linfty_norm());
+
+ energy_groups[group]->solution_old = energy_groups[group]->solution;
+ }
+
+ error = fabs(max_current-max_old)/max_current;
+ max_old = max_current;
+ }
+ else
+ {
+ k_eff = get_total_fission_source();
+ error = fabs(k_eff-k_eff_old)/fabs(k_eff);
+ std::cout << " Iteration " << iteration
+ << ": k_eff=" << k_eff
+ << std::endl;
+ k_eff_old=k_eff;
+
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ {
+ energy_groups[group]->solution_old = energy_groups[group]->solution;
+ energy_groups[group]->solution_old /= k_eff;
+ }
+ }
+
+ ++iteration;
+ }
+ while((error > parameters.convergence_tolerance)
+ &&
+ (iteration < 500));
+
+ for (unsigned int group=0; group<parameters.n_groups; ++group)
+ energy_groups[group]->output_results (cycle);
+
+ std::cout << std::endl;
+ std::cout << " Cycle=" << cycle
+ << ", n_dofs=" << energy_groups[0]->n_dofs() + energy_groups[1]->n_dofs()
+ << ", k_eff=" << k_eff
+ << ", time=" << timer()
+ << std::endl;
+
+
+ std::cout << std::endl << std::endl;
+ }
+}
+
+
+
+ // @sect3{The <code>main()</code> function}
+ //
+ // The last thing in the program in
+ // the <code>main()</code>
+ // function. The structure is as in
+ // most other tutorial programs, with
+ // the only exception that we here
+ // handle a parameter file. To this
+ // end, we first look at the command
+ // line arguments passed to this
+ // function: if no input file is
+ // specified on the command line,
+ // then use "project.prm", otherwise
+ // take the filename given as the
+ // first argument on the command
+ // line.
+ //
+ // With this, we create a
+ // ParameterHandler object, let the
+ // <code>NeutronDiffusionProblem::Parameters</code>
+ // class declare all the parameters
+ // it wants to see in the input file
+ // (or, take the default values, if
+ // nothing is listed in the parameter
+ // file), then read the input file,
+ // ask the parameters object to
+ // extract the values, and finally
+ // hand everything off to an object
+ // of type
+ // <code>NeutronDiffusionProblem</code>
+ // for computation of the eigenvalue:
+int main (int argc, char ** argv)
+{
+ const unsigned int dim = 2;
+
+ try
+ {
+ deallog.depth_console (0);
+
+ std::string filename;
+ if (argc < 2)
+ filename = "project.prm";
+ else
+ filename = argv[1];
+
+
+ ParameterHandler parameter_handler;
+
+ NeutronDiffusionProblem<dim>::Parameters parameters;
+ parameters.declare_parameters (parameter_handler);
+
+ parameter_handler.read_input (filename);
+
+ parameters.get_parameters (parameter_handler);
+
+
+ NeutronDiffusionProblem<dim> neutron_diffusion_problem (parameters);
+ neutron_diffusion_problem.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;
+}
+