From: bangerth Date: Sun, 27 May 2007 02:01:17 +0000 (+0000) Subject: Finally add Yaqi's tutorial program, as step-28. Change references to it. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=01ad676a1b3bf23c82ab3680257e61d4561ccd24;p=dealii-svn.git Finally add Yaqi's tutorial program, as step-28. Change references to it. git-svn-id: https://svn.dealii.org/trunk@14709 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/parameter_handler.h b/deal.II/base/include/base/parameter_handler.h index fe205aecf5..37eb5e4dfb 100644 --- a/deal.II/base/include/base/parameter_handler.h +++ b/deal.II/base/include/base/parameter_handler.h @@ -772,7 +772,7 @@ namespace Patterns * * The ParameterHandler class is discussed in detail in the @ref * step_19 "step-19" example program, and is used in more realistic - * situation in @ref step_22 "step-22". + * situation in @ref step_28 "step-28". * *

Declaring entries

* diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index bad47a009b..1ac18b4633 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -268,7 +268,7 @@ class GridGenerator * * @note For an example of the * use of this function see the - * @ref step_22 "step-22" + * @ref step_28 "step-28" * tutorial program. */ template diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 0c92e0b520..b98c15cd1b 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -849,7 +849,7 @@ namespace internal * This material_id may be set upon construction of a * triangulation (through the CellData data structure), or later * through use of cell iterators. For a typical use of this - * functionality, see the @ref step_22 "step-22" tutorial + * functionality, see the @ref step_28 "step-28" tutorial * program. The functions of the GridGenerator namespace typically * set the material ID of all cells to zero. When reading a * triangulation, the material id must be specified in the input diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index afe9f67b19..26bfa2b015 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -3439,8 +3439,8 @@ class CellAccessor : public TriaObjectAccessor * cell. * * For a typical use of this - * function, see the @ref step_22 - * "step-22" tutorial program. + * function, see the @ref step_28 + * "step-28" tutorial program. */ unsigned char material_id () const; @@ -3449,8 +3449,8 @@ class CellAccessor : public TriaObjectAccessor * cell. * * For a typical use of this - * function, see the @ref step_22 - * "step-22" tutorial program. + * function, see the @ref step_28 + * "step-28" tutorial program. */ void set_material_id (const unsigned char new_material_id) const; diff --git a/deal.II/examples/step-28/Makefile b/deal.II/examples/step-28/Makefile new file mode 100644 index 0000000000..c620808460 --- /dev/null +++ b/deal.II/examples/step-28/Makefile @@ -0,0 +1,154 @@ +# $Id: Makefile,v 1.20 2005/11/23 15:54:25 wolf Exp $ + + +# For the small projects Makefile, you basically need to fill in only +# four fields. +# +# The first is the name of the application. It is assumed that the +# application name is the same as the base file name of the single C++ +# file from which the application is generated. +target = project + +# The second field determines whether you want to run your program in +# debug or optimized mode. The latter is significantly faster, but no +# run-time checking of parameters and internal states is performed, so +# you should set this value to `on' while you develop your program, +# and to `off' when running production computations. +debug-mode = on + + +# As third field, we need to give the path to the top-level deal.II +# directory. You need to adjust this to your needs. Since this path is +# probably the most often needed one in the Makefile internals, it is +# designated by a single-character variable, since that can be +# reference using $D only, i.e. without the parentheses that are +# required for most other parameters, as e.g. in $(target). +D = ../../ + + +# The last field specifies the names of data and other files that +# shall be deleted when calling `make clean'. Object and backup files, +# executables and the like are removed anyway. Here, we give a list of +# files in the various output formats that deal.II supports. +clean-up-files = *gmv *gnuplot *gpl *eps *pov + + + + +# +# +# Usually, you will not need to change something beyond this point. +# +# +# The next statement tell the `make' program where to find the +# deal.II top level directory and to include the file with the global +# settings +include $D/common/Make.global_options + + +# Since the whole project consists of only one file, we need not +# consider difficult dependencies. We only have to declare the +# libraries which we want to link to the object file, and there need +# to be two sets of libraries: one for the debug mode version of the +# application and one for the optimized mode. Here we have selected +# the versions for 2d. Note that the order in which the libraries are +# given here is important and that your applications won't link +# properly if they are given in another order. +# +# You may need to augment the lists of libraries when compiling your +# program for other dimensions, or when using third party libraries +libs.g = $(lib-deal2-2d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-2d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# We now use the variable defined above which switch between debug and +# optimized mode to select the set of libraries to link with. Included +# in the list of libraries is the name of the object file which we +# will produce from the single C++ file. Note that by default we use +# the extension .g.o for object files compiled in debug mode and .o for +# object files in optimized mode (or whatever the local default on your +# system is instead of .o). +ifeq ($(debug-mode),on) + libraries = $(target).g.$(OBJEXT) $(libs.g) +else + libraries = $(target).$(OBJEXT) $(libs.o) +endif + + +# Now comes the first production rule: how to link the single object +# file produced from the single C++ file into the executable. Since +# this is the first rule in the Makefile, it is the one `make' selects +# if you call it without arguments. +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS) + + +# To make running the application somewhat independent of the actual +# program name, we usually declare a rule `run' which simply runs the +# program. You can then run it by typing `make run'. This is also +# useful if you want to call the executable with arguments which do +# not change frequently. You may then want to add them to the +# following rule: +run: $(target) + @echo ============================ Running $< + @./$(target)$(EXEEXT) + + +# As a last rule to the `make' program, we define what to do when +# cleaning up a directory. This usually involves deleting object files +# and other automatically created files such as the executable itself, +# backup files, and data files. Since the latter are not usually quite +# diverse, you needed to declare them at the top of this file. +clean: + -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) + + +# Since we have not yet stated how to make an object file from a C++ +# file, we should do so now. Since the many flags passed to the +# compiler are usually not of much interest, we suppress the actual +# command line using the `at' sign in the first column of the rules +# and write the string indicating what we do instead. +./%.g.$(OBJEXT) : + @echo ==============debug========= $( $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then : else rm $@ ; fi + + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-28/doc/data-q1 b/deal.II/examples/step-28/doc/data-q1 new file mode 100644 index 0000000000..c06f6d7eac --- /dev/null +++ b/deal.II/examples/step-28/doc/data-q1 @@ -0,0 +1,10 @@ + 0 2450 0.907410259645 4.077380000000 + 1 4289 0.907042395776 10.572393000000 + 2 8225 0.906937183449 25.223165000000 + 3 14757 0.906900660844 53.679839000000 + 4 26236 0.906886610239 108.345529000000 + 5 62804 0.906857043552 253.919398000000 + 6 86389 0.906844968762 466.115140000000 + 7 165629 0.906839125266 900.314131000000 + 8 301264 0.906836698748 1776.255968000000 + 9 444166 0.906835969493 3026.378920000000 diff --git a/deal.II/examples/step-28/doc/data-q2 b/deal.II/examples/step-28/doc/data-q2 new file mode 100644 index 0000000000..df6c0e0665 --- /dev/null +++ b/deal.II/examples/step-28/doc/data-q2 @@ -0,0 +1,10 @@ + 0 9522 0.906841960371 17.634319000000 + 1 15428 0.906837901031 33.703876000000 + 2 23566 0.906836075928 65.945975000000 + 3 33930 0.906835500110 120.341705000000 + 4 60358 0.906835001796 246.337551000000 + 5 86798 0.906834858174 403.059725000000 + 6 98538 0.906834824060 545.820023000000 + 7 117254 0.906834787555 734.265375000000 + 8 147302 0.906834761604 1014.816724000000 + 9 194442 0.906834746216 1366.371280000000 diff --git a/deal.II/examples/step-28/doc/data-q3 b/deal.II/examples/step-28/doc/data-q3 new file mode 100644 index 0000000000..4b9d50616c --- /dev/null +++ b/deal.II/examples/step-28/doc/data-q3 @@ -0,0 +1,10 @@ + 0 21218 0.906835491999 91.037161000000 + 1 34067 0.906834908620 131.201055000000 + 2 35069 0.906834859254 187.733460000000 + 3 50476 0.906834817384 279.669484000000 + 4 73375 0.906834783106 408.506898000000 + 5 90529 0.906834744225 590.097292000000 + 6 114417 0.906834730885 852.591387000000 + 7 151059 0.906834726355 1161.735390000000 + 8 215292 0.906834723745 1576.734300000000 + 9 275614 0.906834722268 2160.737518000000 diff --git a/deal.II/examples/step-28/doc/data-q4 b/deal.II/examples/step-28/doc/data-q4 new file mode 100644 index 0000000000..15a4abe16c --- /dev/null +++ b/deal.II/examples/step-28/doc/data-q4 @@ -0,0 +1,10 @@ + 0 37538 0.906834753543 373.753181000000 + 1 64160 0.906834728274 517.920264000000 + 2 88822 0.906834726164 1138.467927000000 + 3 157266 0.906834722755 2460.279981000000 + 4 175774 0.906834722205 4364.633475000000 + 5 227672 0.906834721746 4964.747243999999 + 6 264428 0.906834721654 5461.831674999999 + 7 328106 0.906834721515 6010.156317000000 + 8 395036 0.906834721416 12474.065654000000 + 9 455498 0.906834721380 19863.514286999998 diff --git a/deal.II/examples/step-28/doc/data-q5 b/deal.II/examples/step-28/doc/data-q5 new file mode 100644 index 0000000000..c4aa83ff28 --- /dev/null +++ b/deal.II/examples/step-28/doc/data-q5 @@ -0,0 +1,10 @@ + 0 58482 0.906834731399 1121.873449000000 + 1 97921 0.906834723230 1409.121780000000 + 2 98655 0.906834723009 2989.811479000000 + 3 131934 0.906834722626 3469.306585000000 + 4 188259 0.906834722379 4011.864103000000 + 5 255493 0.906834721986 4578.781919000000 + 6 271824 0.906834721755 10632.535608000000 + 7 307409 0.906834721567 11699.187453000000 + 8 316691 0.906834721509 12528.627359000000 + 9 363748 0.906834721439 13773.087172000000 diff --git a/deal.II/examples/step-28/doc/data-q6 b/deal.II/examples/step-28/doc/data-q6 new file mode 100644 index 0000000000..4c5094d228 --- /dev/null +++ b/deal.II/examples/step-28/doc/data-q6 @@ -0,0 +1,15 @@ + 0 84050 0.906834722842531 3184.047951000000012 + 1 140828 0.906834721670094 3965.101212999999916 + 2 158782 0.906834721605666 8725.235561999999845 + 3 230198 0.906834721530541 16501.528386000001774 + 4 320258 0.906834721476372 18215.300853000000643 + 5 385470 0.906834721418083 21774.706741000001784 + 6 438090 0.906834721389658 25271.607131000000663 + 7 462942 0.906834721380828 27951.595711000001756 + 8 502434 0.906834721371133 33127.243892999998934 + 9 522814 0.906834721367467 36477.579565000000002 + 10 552538 0.906834721365696 39634.395656000000599 + 11 570310 0.906834721364977 43278.281700999999885 + 12 593214 0.906834721364745 47564.244135000000824 + 13 654882 0.906834721363022 53225.507492000004277 + 14 716338 0.906834721360716 60232.646245000003546 diff --git a/deal.II/examples/step-28/doc/gnuplot.1 b/deal.II/examples/step-28/doc/gnuplot.1 new file mode 100644 index 0000000000..42a47f1d78 --- /dev/null +++ b/deal.II/examples/step-28/doc/gnuplot.1 @@ -0,0 +1,26 @@ +set term png +set logsc xy +set da sty lp + +set ylabel "Error" + +set output "error-vs-dofs.png" +set xlabel "Degrees of freedom" + +pl "data-q1" us 2:(abs($3-0.906834721360716)) title "Q1 elements", \ + "data-q2" us 2:(abs($3-0.906834721360716)) title "Q2 elements", \ + "data-q3" us 2:(abs($3-0.906834721360716)) title "Q3 elements", \ + "data-q4" us 2:(abs($3-0.906834721360716)) title "Q4 elements", \ + "data-q5" us 2:(abs($3-0.906834721360716)) title "Q5 elements", \ + "data-q6" us 2:(abs($3-0.906834721360716)) title "Q6 elements" + + +set output "error-vs-time.png" +set xlabel "Time in seconds" + +pl "data-q1" us 4:(abs($3-0.906834721360716)) title "Q1 elements", \ + "data-q2" us 4:(abs($3-0.906834721360716)) title "Q2 elements", \ + "data-q3" us 4:(abs($3-0.906834721360716)) title "Q3 elements", \ + "data-q4" us 4:(abs($3-0.906834721360716)) title "Q4 elements", \ + "data-q5" us 4:(abs($3-0.906834721360716)) title "Q5 elements", \ + "data-q6" us 4:(abs($3-0.906834721360716)) title "Q6 elements" diff --git a/deal.II/examples/step-28/doc/intro.dox b/deal.II/examples/step-28/doc/intro.dox new file mode 100644 index 0000000000..f693a95d3c --- /dev/null +++ b/deal.II/examples/step-28/doc/intro.dox @@ -0,0 +1 @@ +

Introduction

diff --git a/deal.II/examples/step-28/doc/intro.tex b/deal.II/examples/step-28/doc/intro.tex new file mode 100644 index 0000000000..4516336213 --- /dev/null +++ b/deal.II/examples/step-28/doc/intro.tex @@ -0,0 +1,634 @@ +\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$ 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' + \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} diff --git a/deal.II/examples/step-28/doc/results.dox b/deal.II/examples/step-28/doc/results.dox new file mode 100644 index 0000000000..f7a468d988 --- /dev/null +++ b/deal.II/examples/step-28/doc/results.dox @@ -0,0 +1,500 @@ +

Results

+ + +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. diff --git a/deal.II/examples/step-28/doc/step-28.convergence.png b/deal.II/examples/step-28/doc/step-28.convergence.png new file mode 100755 index 0000000000..65619bf273 Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.convergence.png differ diff --git a/deal.II/examples/step-28/doc/step-28.error-vs-dofs.png b/deal.II/examples/step-28/doc/step-28.error-vs-dofs.png new file mode 100644 index 0000000000..0240f378d0 Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.error-vs-dofs.png differ diff --git a/deal.II/examples/step-28/doc/step-28.error-vs-time.png b/deal.II/examples/step-28/doc/step-28.error-vs-time.png new file mode 100644 index 0000000000..f3f8924cd4 Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.error-vs-time.png differ diff --git a/deal.II/examples/step-28/doc/step-28.grid-0.9.order2.png b/deal.II/examples/step-28/doc/step-28.grid-0.9.order2.png new file mode 100755 index 0000000000..439a916ca1 Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.grid-0.9.order2.png differ diff --git a/deal.II/examples/step-28/doc/step-28.grid-1.9.order2.png b/deal.II/examples/step-28/doc/step-28.grid-1.9.order2.png new file mode 100755 index 0000000000..72129326dd Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.grid-1.9.order2.png differ diff --git a/deal.II/examples/step-28/doc/step-28.solution-0.9.order2.png b/deal.II/examples/step-28/doc/step-28.solution-0.9.order2.png new file mode 100755 index 0000000000..a01e2381da Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.solution-0.9.order2.png differ diff --git a/deal.II/examples/step-28/doc/step-28.solution-1.9.order2.png b/deal.II/examples/step-28/doc/step-28.solution-1.9.order2.png new file mode 100755 index 0000000000..aa5431851e Binary files /dev/null and b/deal.II/examples/step-28/doc/step-28.solution-1.9.order2.png differ diff --git a/deal.II/examples/step-28/step-28.cc b/deal.II/examples/step-28/step-28.cc new file mode 100644 index 0000000000..304abc2483 --- /dev/null +++ b/deal.II/examples/step-28/step-28.cc @@ -0,0 +1,2347 @@ +//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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + + // 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 + + // 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 + + // 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 + + // 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 +#include + + // 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 + // n_materials=8 + // 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; mEnergyGroup 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-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' +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 + // setup_linear_system + // 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 &coarse_grid, + const FiniteElement &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 &extraneous_source); + void assemble_cross_group_rhs (const EnergyGroup &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 + // NeutronDiffusionProblem, + // 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 &error_indicators) const; + + void refine_grid (const Vector &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 solution; + Vector 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 triangulation; + const FiniteElement &fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector system_rhs; + + std::map 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 &g_prime, + const typename DoFHandler::cell_iterator &cell_g, + const typename DoFHandler::cell_iterator &cell_g_prime, + const FullMatrix prolongation_matrix); +}; + + + // @sect4{Implementation of the EnergyGroup 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 +EnergyGroup::EnergyGroup (const unsigned int group, + const MaterialData &material_data, + const Triangulation &coarse_grid, + const FiniteElement &fe) + : + group (group), + material_data (material_data), + fe (fe), + dof_handler (triangulation) +{ + triangulation.copy_triangulation (coarse_grid); + dof_handler.distribute_dofs (fe); +} + + + +template +unsigned int +EnergyGroup::n_active_cells () const +{ + return triangulation.n_active_cells (); +} + + + +template +unsigned int +EnergyGroup::n_dofs () const +{ + return dof_handler.n_dofs (); +} + + + + // @sect5{EnergyGroup::setup_linear_system} + // + // 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 + // EnergyGroup::refine_grid + // function). +template +void +EnergyGroup::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 + // GridGenerator::subdivided_hyper_rectangle + // function (in + // NeutronDiffusionProblem::initialize_problem) + // where we set the last parameter + // to true. 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 + // VectorTools::interpolate_boundary_values + // are additive, i.e. they do not + // first clear the boundary value + // map): + boundary_values.clear(); + + for (unsigned int i=0; i(), + boundary_values); +} + + + + // @sect5{EnergyGroup::assemble_system_matrix} + // + // 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 + // cell->material_id() 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 EnergyGroup::solve + // function, at which point all the + // information is available. +template +void +EnergyGroup::assemble_system_matrix () +{ + const QGauss quadrature_formula(fe.degree + 1); + + FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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_pointget_dof_indices (local_dof_indices); + + for (unsigned int i=0; iEnergyGroup::assemble_ingroup_rhs} + // + // As explained in the documentation + // of the EnergyGroup 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 +void EnergyGroup::assemble_ingroup_rhs (const Function &extraneous_source) +{ + system_rhs.reinit (dof_handler.n_dofs()); + + const QGauss 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 fe_values (fe, quadrature_formula, + update_values | update_q_points | + update_JxW_values); + + Vector cell_rhs (dofs_per_cell); + std::vector extraneous_source_values (n_q_points); + std::vector solution_old_values (n_q_points); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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_pointEnergyGroup::assemble_cross_group_rhs} + // + // 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 + // get_finest_common_cells 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 +void EnergyGroup::assemble_cross_group_rhs (const EnergyGroup &g_prime) +{ + if (group == g_prime.group) + return; + + const std::list::cell_iterator, + typename DoFHandler::cell_iterator> > + cell_list + = GridTools::get_finest_common_cells (dof_handler, + g_prime.dof_handler); + + typename std::list::cell_iterator, + typename DoFHandler::cell_iterator> > + ::const_iterator + cell_iter = cell_list.begin(); + + for (; cell_iter!=cell_list.end(); ++cell_iter) + { + FullMatrix unit_matrix (fe.dofs_per_cell); + for (unsigned int i=0; ifirst, + cell_iter->second, + unit_matrix); + } +} + + + + // @sect5{EnergyGroup::assemble_cross_group_rhs_recursive} + // + // 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 +void +EnergyGroup:: +assemble_cross_group_rhs_recursive (const EnergyGroup &g_prime, + const typename DoFHandler::cell_iterator &cell_g, + const typename DoFHandler::cell_iterator &cell_g_prime, + const FullMatrix 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 quadrature_formula (fe.degree+1); + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + FEValues 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 local_mass_matrix_f (fe.dofs_per_cell, + fe.dofs_per_cell); + FullMatrix local_mass_matrix_g (fe.dofs_per_cell, + fe.dofs_per_cell); + + for (unsigned int q_point=0; q_pointvmult + // function, or the product with the + // transpose matrix using Tvmult. + // After doing so, we transfer the + // result into the global right hand + // side vector of energy group $g$. + Vector g_prime_new_values (fe.dofs_per_cell); + Vector 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 cell_rhs (fe.dofs_per_cell); + Vector 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 local_dof_indices (fe.dofs_per_cell); + cell_g->get_dof_indices (local_dof_indices); + + for (unsigned int i=0; immult), 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::children_per_cell;++child) + { + FullMatrix 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{EnergyGroup::get_fission_source} + // + // 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 +double EnergyGroup::get_fission_source () const +{ + const QGauss quadrature_formula (fe.degree + 1); + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + FEValues fe_values (fe, quadrature_formula, + update_values | update_JxW_values); + + std::vector solution_values (n_q_points); + + double fission_source = 0; + + typename DoFHandler::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_pointEnergyGroup::solve} + // + // 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 +void +EnergyGroup::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{EnergyGroup::estimate_errors} + // + // 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 +void EnergyGroup::estimate_errors (Vector &error_indicators) const +{ + KellyErrorEstimator::estimate (dof_handler, + QGauss (fe.degree + 1), + typename FunctionMap::type(), + solution, + error_indicators); + error_indicators /= solution.linfty_norm(); +} + + + + // @sect5{EnergyGroup::refine_grid} + // + // 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 + // GridRefinement 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 +void EnergyGroup::refine_grid (const Vector &error_indicators, + const double refine_threshold, + const double coarsen_threshold) +{ + typename Triangulation::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 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{EnergyGroup::output_results} + // + // 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 + // Utilities::int_to_string + // 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 +void +EnergyGroup::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 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 NeutronDiffusionProblem 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 + // EnergyGroup 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 run + // 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 + // declare_parameters + // 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 + // get_parameters that + // can extract the values of these + // parameters from a ParameterHandler + // object. +template +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 + // EnergyGroup class + // and is simply called from the + // run() 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 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*> energy_groups; +}; + + + // @sect4{Implementation of the NeutronDiffusionProblem::Parameters 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 +NeutronDiffusionProblem::Parameters::Parameters () + : + n_groups (2), + n_refinement_cycles (5), + fe_degree (2), + convergence_tolerance (1e-12) +{} + + + +template +void +NeutronDiffusionProblem::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 +void +NeutronDiffusionProblem::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 NeutronDiffusionProblem class} + + // Now for the + // NeutronDiffusionProblem + // class. The constructor and + // destructor have nothing of much + // interest: +template +NeutronDiffusionProblem:: +NeutronDiffusionProblem (const Parameters ¶meters) + : + parameters (parameters), + material_data (parameters.n_groups), + fe (parameters.fe_degree) +{} + + + +template +NeutronDiffusionProblem::~NeutronDiffusionProblem () +{ + for (unsigned int group=0; groupNeutronDiffusionProblem::initialize_problem} + // + // 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 + // GridGenerator::subdivided_hyper_rectangle + // 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 +void NeutronDiffusionProblem::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 bottom_left = Point(); + const Point upper_right = (dim == 2 + ? + Point (assemblies_x*rods_per_assembly_x*pin_pitch_x, + assemblies_y*rods_per_assembly_y*pin_pitch_y) + : + Point (assemblies_x*rods_per_assembly_x*pin_pitch_x, + assemblies_y*rods_per_assembly_y*pin_pitch_y, + assemblies_z*assembly_height)); + + std::vector 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 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): + //
    + //
  1. 'UX' Assembly: UO2 fuel assembly + // with 24 guide tubes and a central + // Moveable Fission Chamber + //
  2. 'UA' Assembly: UO2 fuel assembly + // with 24 AIC and a central + // Moveable Fission Chamber + //
  3. 'PX' Assembly: MOX fuel assembly + // with 24 guide tubes and a central + // Moveable Fission Chamber + //
  4. 'R' Assembly: a reflector. + //
+ // + // 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::active_cell_iterator + cell = coarse_grid.begin_active(); + cell!=coarse_grid.end(); + ++cell) + { + const Point 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 (group, material_data, + coarse_grid, fe); +} + + + // @sect5{NeutronDiffusionProblem::get_total_fission_source} + // + // 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 + // EnergyGroup 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@ + // if the function spawned returns a + // double. We can the later ask this + // object for the returned value + // (when doing so, the + // Threads::Thread@::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 +double NeutronDiffusionProblem::get_total_fission_source () const +{ + std::vector > threads; + for (unsigned int group=0; group::get_fission_source) ()); + + double fission_source = 0; + for (unsigned int group=0; groupNeutronDiffusionProblem::refine_grid} + // + // 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 +void NeutronDiffusionProblem::refine_grid () +{ + std::vector n_cells (parameters.n_groups); + for (unsigned int group=0; groupn_active_cells(); + + BlockVector group_error_indicators(n_cells); + + { + Threads::ThreadGroup<> threads; + for (unsigned int group=0; group::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::refine_grid) + (group_error_indicators.block(group), + refine_threshold, + coarsen_threshold); + threads.join_all (); + } +} + + + // @sect5{NeutronDiffusionProblem::run} + // + // 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 +void NeutronDiffusionProblem::run () +{ + std::cout << std::setprecision (12) << std::fixed; + + double k_eff_old = k_eff; + + Timer timer; + timer.start (); + + for (unsigned int cycle=0; cyclesolution *= k_eff; + } + + for (unsigned int group=0; groupsetup_linear_system (); + + std::cout << " Numbers of active cells: "; + for (unsigned int group=0; groupn_active_cells() + << ' '; + std::cout << std::endl; + std::cout << " Numbers of degrees of freedom: "; + for (unsigned int group=0; groupn_dofs() + << ' '; + std::cout << std::endl << std::endl; + + + Threads::ThreadGroup<> threads; + for (unsigned int group=0; group::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; groupassemble_ingroup_rhs (ZeroFunction()); + + for (unsigned int bgroup=0; bgroupassemble_cross_group_rhs (*energy_groups[bgroup]); + + energy_groups[group]->solve (); + } + + if (isour==1) + { + double max_current = 0; + for (unsigned int group=0; groupsolution.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; groupsolution_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; groupoutput_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 main() function} + // + // The last thing in the program in + // the main() + // 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 + // NeutronDiffusionProblem::Parameters + // 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 + // NeutronDiffusionProblem + // 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::Parameters parameters; + parameters.declare_parameters (parameter_handler); + + parameter_handler.read_input (filename); + + parameters.get_parameters (parameter_handler); + + + NeutronDiffusionProblem 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; +} +