diff --git a/docs/source/_images/loss.png b/docs/source/_images/loss.png new file mode 100644 index 0000000000..37d4169a18 Binary files /dev/null and b/docs/source/_images/loss.png differ diff --git a/docs/source/_images/prod.png b/docs/source/_images/prod.png new file mode 100644 index 0000000000..a4c2887d8f Binary files /dev/null and b/docs/source/_images/prod.png differ diff --git a/docs/source/conf.py b/docs/source/conf.py index 1943432165..16be3567d0 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -23,7 +23,7 @@ sys.path.insert(0, os.path.abspath('../sphinxext')) # Add any Sphinx extension module names here, as strings. They can be extensions # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. -extensions = ['sphinx.ext.pngmath'] +extensions = ['sphinx.ext.pngmath', 'sphinxcontrib.tikz'] # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -188,7 +188,14 @@ latex_documents = [ u'Massachusetts Institute of Technology', 'manual'), ] -latex_elements = {'preamble': '\\usepackage{enumitem}\\setlistdepth{9}'} +latex_elements = { +'preamble': ''' +\usepackage{enumitem} +\setlistdepth{9} +\usepackage{tikz} +\usetikzlibrary{shapes,snakes,shadows,arrows,calc,decorations.markings,patterns,fit,matrix,spy} +''' +} # The name of an image file (relative to this directory) to place at the top of # the title page. diff --git a/docs/source/methods/cmfd.rst b/docs/source/methods/cmfd.rst new file mode 100644 index 0000000000..3681df3b75 --- /dev/null +++ b/docs/source/methods/cmfd.rst @@ -0,0 +1,561 @@ +.. _methods_cmfd: + +================================================================ +Nonlinear Diffusion Acceleration - Coarse Mesh Finite Difference +================================================================ + +This page section discusses how nonlinear diffusion acceleration (NDA) using +coarse mesh finite difference (CMFD) is implemented into OpenMC. Before we get +into the theory, general notation for this section is discussed. + +-------- +Notation +-------- + +Before deriving NDA relationships, notation is explained. If a parameter has a +:math:`\overline{\cdot}`, it is surface area-averaged and if it has a +:math:`\overline{\overline\cdot}`, it is volume-averaged. When describing a +specific cell in the geometry, indices :math:`(i,j,k)` are used which correspond +to directions :math:`(x,y,z)`. In most cases, the same operation is performed in +all three directions. To compactly write this, an arbitrary direction set +:math:`(u,v,w)` that corresponds to cell indices :math:`(l,m,n)` is used. Note +that :math:`u` and :math:`l` do not have to correspond to :math:`x` and +:math:`i`. However, if :math:`u` and :math:`l` correspond to :math:`y` and +:math:`j`, :math:`v` and :math:`w` correspond to :math:`x` and :math:`z` +directions. An example of this is shown in the following expression: + +.. math:: + :label: not1 + + \sum\limits_{u\in(x,y,z)}\left\langle\overline{J}^{u,g}_{l+1/2,m,n} + \Delta_m^v\Delta_n^w\right\rangle + +Here, :math:`u` takes on each direction one at a time. The parameter :math:`J` +is surface area-averaged over the transverse indices :math:`m` and :math:`n` +located at :math:`l+1/2`. Usually, spatial indices are listed as subscripts and +the direction as a superscript. Energy group indices represented by :math:`g` +and :math:`h` are also listed as superscripts here. The group :math:`g` is the +group of interest and, if present, :math:`h` is all groups. Finally, any +parameter surrounded by :math:`\left\langle\cdot\right\rangle` represents a +tally quantity that can be edited from a Monte Carlo (MC) solution. + +------ +Theory +------ + +NDA is a diffusion model that has equivalent physics to a transport model. There +are many different methods that can be classified as NDA. The CMFD method is a +type of NDA that represents second order multigroup diffusion equations on a +coarse spatial mesh. Whether a transport model or diffusion model is used to +represent the distribution of neutrons, these models must satisfy the *neutron +balance equation*. This balance is represented by the following formula for a +specific energy group :math:`g` in cell :math:`(l,m,n)`: + +.. math:: + :label: eq_neut_bal + + \sum\limits_{u\in(x,y,z)}\left(\left\langle\overline{J}^{u,g}_{l+1/2,m,n} + \Delta_m^v\Delta_n^w\right\rangle - + \left\langle\overline{J}^{u,g}_{l-1/2,m,n} + \Delta_m^v\Delta_n^w\right\rangle\right) + + + \left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g + \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle + = \\ + \sum\limits_{h=1}^G\left\langle + \overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow + g}\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w + \right\rangle + + + \frac{1}{k_{eff}}\sum\limits_{h=1}^G + \left\langle\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow + g}\overline{\overline\phi}_{l,m,n}^h + \Delta_l^u\Delta_m^v\Delta_n^w\right\rangle. + +In eq. :eq:`eq_neut_bal` the parameters are defined as: + +* :math:`\left\langle\overline{J}^{u,g}_{l\pm + 1/2,m,n}\Delta_m^v\Delta_n^w\right\rangle` --- surface area-integrated net + current over surface :math:`(l\pm 1/2,m,n)` with surface normal in direction + :math:`u` in energy group :math:`g`. By dividing this quantity by the transverse + area, :math:`\Delta_m^v\Delta_n^w`, the surface area-averaged net current can + be computed. +* :math:`\left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g + \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` + --- volume-integrated total reaction rate over energy group :math:`g`. +* :math:`\left\langle\overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow + g} + \overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` + --- volume-integrated scattering production rate of neutrons that begin with + energy in group :math:`h` and exit reaction in group :math:`g`. This reaction + rate also includes the energy transfer of reactions (except fission) that + produce multiple neutrons such as (n, 2n); hence, the need for :math:`\nu_s` + to represent neutron multiplicity. +* :math:`k_{eff}` --- core multiplication factor. +* :math:`\left\langle\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow + g}\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` + --- volume-integrated fission production rate of neutrons from fissions in + group :math:`h` that exit in group :math:`g`. + +Each quantity in :math:`\left\langle\cdot\right\rangle` represents a scalar value that +is obtained from an MC tally. A good verification step when using an MC code is +to make sure that tallies satisfy this balance equation within statistics. No +NDA acceleration can be performed if the balance equation is not satisfied. + +There are three major steps to consider when performing NDA: (1) calculation of +macroscopic cross sections and nonlinear parameters, (2) solving an eigenvalue +problem with a system of linear equations, and (3) modifying MC source +distribution to align with the NDA solution on a chosen mesh. This process is +illustrated as a flow chart below. After a batch of neutrons +is simulated, NDA can take place. Each of the steps described above is described +in detail in the following sections. + +.. tikz:: Flow chart of NDA process. Note "XS" is used for cross section and + "DC" is used for diffusion coefficient. + :libs: shapes, snakes, shadows, arrows, calc, decorations.markings, patterns, fit, matrix, spy + :include: cmfd_tikz/cmfd_flow.tikz + +Calculation of Macroscopic Cross Sections +----------------------------------------- + +A diffusion model needs macroscopic cross sections and diffusion coefficients to +solve for multigroup fluxes. Cross sections are derived by conserving reaction +rates predicted by MC tallies. From Eq. :eq:`eq_neut_bal`, total, scattering +production and fission production macroscopic cross sections are needed. They are +defined from MC tallies as follows: + +.. math:: + :label: xs1 + + \overline{\overline\Sigma}_{t_{l,m,n}}^g \equiv + \frac{\left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g + \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle} + {\left\langle\overline{\overline\phi}_{l,m,n}^g + \Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}, + +.. math:: + :label: xs2 + + \overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow g} \equiv + \frac{\left\langle\overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow + g}\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle} + {\left\langle\overline{\overline\phi}_{l,m,n}^h + \Delta_l^u\Delta_m^v\Delta_n^w\right\rangle} + +and + +.. math:: + :label: xs3 + + \overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow g} \equiv + \frac{\left\langle\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow + g}\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle} + {\left\langle\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}. + +In order to fully conserve neutron balance, leakage rates also need to be +preserved. In standard diffusion theory, leakage rates are represented by +diffusion coefficients. Unfortunately, it is not easy in MC to calculate a +single diffusion coefficient for a cell that describes leakage out of each +surface. Luckily, it does not matter what definition of diffusion coefficient is +used because nonlinear equivalence parameters will correct for this +inconsistency. However, depending on the diffusion coefficient definition +chosen, different convergence properties of NDA equations are observed. +Here, we introduce a diffusion coefficient that is derived for a coarse energy +transport reaction rate. This definition can easily be constructed from +MC tallies provided that angular moments of scattering reaction rates can +be obtained. The diffusion coefficient is defined as follows: + +.. math:: + :label: eq_transD + + \overline{\overline D}_{l,m,n}^g = + \frac{\left\langle\overline{\overline\phi}_{l,m,n}^g + \Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}{3 + \left\langle\overline{\overline\Sigma}_{tr_{l,m,n}}^g + \overline{\overline\phi}_{l,m,n}^g + \Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}, + +where + +.. math:: + :label: xs4 + + \left\langle\overline{\overline\Sigma}_{tr_{l,m,n}}^g + \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle + = + \left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g + \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle + \\ - + \left\langle\overline{\overline{\nu_s\Sigma}}_{s1_{l,m,n}}^g + \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle. + +Note that the transport reaction rate is calculated from the total reaction rate +reduced by the :math:`P_1` scattering production reaction rate. Equation :eq:`eq_transD` +does not represent the best definition of diffusion coefficients from MC; +however, it is very simple and usually fits into MC tally frameworks +easily. Different methods to calculate more accurate diffusion coefficients can +found in [Herman]_. + +CMFD Equations +-------------- + +The first part of this section is devoted to discussing second-order finite +volume discretization of multigroup diffusion equations. This will be followed +up by the formulation of CMFD equations that are used in this NDA +scheme. When performing second-order finite volume discretization of the +diffusion equation, we need information that relates current to flux. In this +numerical scheme, each cell is coupled only to its direct neighbors. Therefore, +only two types of coupling exist: (1) cell-to-cell coupling and (2) +cell-to-boundary coupling. The derivation of this procedure is referred to as +finite difference diffusion equations and can be found in literature such +as [Hebert]_. These current/flux relationships are as follows: + +* cell-to-cell coupling + +.. math:: + :label: eq_cell_cell + + \overline{J}^{u,g}_{l\pm1/2,m,n} = -\frac{2\overline{\overline + D}_{l\pm1,m,n}^g\overline{\overline + D}_{l,m,n}^g}{\overline{\overline D}_{l\pm1,m,n}^g\Delta_l^u + + \overline{\overline + D}_{l,m,n}^g\Delta_{l\pm1}^u} + \left(\pm\overline{\overline{\phi}}_{l\pm1,m,n}^g\mp + \overline{\overline{\phi}}_{l,m,n}^g\right), + +* cell-to-boundary coupling + +.. math:: + :label: eq_cell_bound + + \overline{J}^{u,g}_{l\pm1/2,m,n} = \pm\frac{2\overline{\overline + D}_{l,m,n}^g\left(1 - + \beta_{l\pm1/2,m,n}^{u,g}\right)}{4\overline{\overline + D}_{l,m,n}^g\left(1 + \beta_{l\pm1/2,m,n}^{u,g}\right) + \left(1 - + \beta_{l\pm1/2,m,n}^{u,g}\right)\Delta_l^u}\overline{\overline{\phi}}_{l,m,n}^{g}. + +In Eqs. :eq:`eq_cell_cell` and :eq:`eq_cell_bound`, the :math:`\pm` refers to +left (:math:`-x`) or right (:math:`+x`) surface in the :math:`x` direction, +back (:math:`-y`) or front (:math:`+y`) surface in the :math:`y` direction and +bottom (:math:`-z`) or top (:math:`+z`) surface in the :math:`z` direction. For +cell-to-boundary coupling, a general albedo, :math:`\beta_{l\pm1/2,m,n}^{u,g}`, +is used. The albedo is defined as the ratio of incoming (:math:`-` superscript) +to outgoing (:math:`+` superscript) partial current on any surface represented +as + +.. math:: + :label: eq_albedo + + \beta_{l\pm1/2,m,n}^{u,g} = + \frac{\overline{J}^{u,g-}_{l\pm1/2,m,n}}{\overline{J}^{u,g+}_{l\pm1/2,m,n}}. + +Common boundary conditions are: vacuum (:math:`\beta=0`), reflective +(:math:`\beta=1`) and zero flux (:math:`\beta=-1`). Both eq. :eq:`eq_cell_cell` +and eq. :eq:`eq_cell_bound` can be written in this generic form, + +.. math:: + :label: eq_dtilde + + \overline{J}^{u,g}_{l\pm1/2,m,n} = \widetilde{D}_{l,m,n}^{u,g} \left(\dots\right). + +The parameter :math:`\widetilde{D}_{l,m,n}^{u,g}` represents the linear +coupling term between current and flux. These current relationships can be +sustituted into eq. :eq:`eq_neut_bal` to produce a linear system of multigroup +diffusion equations for each spatial cell and energy group. However, a solution +to these equations is not consistent with a higher order transport solution +unless equivalence factors are present. This is because both the diffusion +approximation, governed by Fick's Law, and spatial trunction error will produce +differences. Therefore, a nonlinear parameter, +:math:`\widehat{D}_{l,m,n}^{u,g}`, is added to eqs. :eq:`eq_cell_cell` and +:eq:`eq_cell_bound`. These equations are, respectively, + +.. math:: + :label: eq_dhat_cell + + \overline{J}^{u,g}_{l\pm1/2,m,n} = -\widetilde{D}_{l,m,n}^{u,g} + \left(\pm\overline{\overline{\phi}}_{l\pm1,m,n}^g\mp + \overline{\overline{\phi}}_{l,m,n}^g\right) + \widehat{D}_{l,m,n}^{u,g} + \left(\overline{\overline{\phi}}_{l\pm1,m,n}^g + + \overline{\overline{\phi}}_{l,m,n}^g\right) + +and + +.. math:: + :label: eq_dhat_bound + + \overline{J}^{u,g}_{l\pm1/2,m,n} = \pm\widetilde{D}_{l,m,n}^{u,g} + \overline{\overline{\phi}}_{l,m,n}^{g} + \widehat{D}_{l,m,n}^{u,g} + \overline{\overline{\phi}}_{l,m,n}^{g}. + +The only unknown in each of these equations is the equivalence parameter. The +current, linear coupling term and flux can either be obtained or derived from +MC tallies. Thus, it is called nonlinear because it is dependent on the flux +which is updated on the next iteration. + +Equations :eq:`eq_dhat_cell` and :eq:`eq_dhat_bound` can be substituted into +eq. :eq:`eq_neut_bal` to create a linear system of equations that is consistent +with transport physics. One example of this equation is written for an +interior cell, + +.. math:: + :label: eq_cmfd_sys + + \sum_{u\in + x,y,x}\frac{1}{\Delta_l^u}\left[\left(-\tilde{D}_{l-1/2,m,n}^{u,g} - + \hat{D}_{l-1/2,m,n}^{u,g}\right)\overline{\overline{\phi}}_{l-1,m,n}^g\right. + + \left(\tilde{D}_{l-1/2,m,n}^{u,g} + + \tilde{D}_{l+1/2,m,n}^{u,g} - \hat{D}_{l-1/2,m,n}^{u,g} + + \hat{D}_{l+1/2,m,n}^{u,g}\right)\overline{\overline{\phi}}_{l,m,n}^g + \\ + + \left. \left(-\tilde{D}_{l+1/2,m,n}^{u,g} + + \hat{D}_{l+1/2,m,n}^{u,g}\right)\overline{\overline{\phi}}_{l+1,m,n}^g + \right] + + \overline{\overline\Sigma}_{t_{l,m,n}}^g\overline{\overline{\phi}}_{l,m,n}^g + - \sum\limits_{h=1}^G\overline{\overline{\nu_s\Sigma}}^{h\rightarrow + g}_{s_{l,m,n}}\overline{\overline{\phi}}_{l,m,n}^h = + \frac{1}{k}\sum\limits_{h=1}^G\overline{\overline{\nu_f\Sigma}}^{h\rightarrow + g}_{f_{l,m,n}}\overline{\overline{\phi}}_{l,m,n}^h. + +It should be noted that before substitution, eq. :eq:`eq_neut_bal` was divided +by the volume of the cell, :math:`\Delta_l^u\Delta_m^v\Delta_n^w`. Equation +:eq:`eq_cmfd_sys` can be represented in operator form as + +.. math:: + :label: eq_CMFDopers + + \mathbb{M}\mathbf{\Phi} = \frac{1}{k}\mathbb{F}\mathbf{\Phi}, + +where :math:`\mathbb{M}` is the neutron loss matrix operator, +:math:`\mathbb{F}` is the neutron production matrix operator, +:math:`\mathbf{\Phi}` is the multigroup flux vector and :math:`k` is the +eigenvalue. This generalized eigenvalue problem is solved to obtain fundamental +mode multigroup fluxes and eigenvalue. In order to produce consistent results +with transport theory from these equations, the neutron balance equation must +have been satisfied by MC tallies. The desire is that CMFD equations will +produce a more accurate source than MC after each fission source generation. + +CMFD Feedback +------------- + +Now that a more accurate representation of the expected source distribution is +estimated from CMFD, it needs to be communicated back to MC. The first step +in this process is to generate a probability mass function that provides +information about how probable it is for a neutron to be born in a given cell +and energy group. This is represented as + +.. math:: + :label: eq_cmfd_psrc + + p_{l,m,n}^g = + \frac{\sum_{h=1}^{G}\overline{\overline{\nu_f\Sigma}}^{h\rightarrow + g}_{f_{l,m,n}}\overline{\overline{\phi}}_{l,m,n}^h\Delta_l^u\Delta_m^v + \Delta_n^w}{\sum_n\sum_m\sum_l\sum_{h=1}^{G}\overline{ + \overline{\nu_f\Sigma}}^{h\rightarrow + g}_{f_{l,m,n}}\overline{\overline{\phi}}_{l,m,n}^h\Delta_l^u\Delta_m^v + \Delta_n^w}. + +This equation can be multiplied by the number of source neutrons to obtain an +estimate of the expected number of neutrons to be born in a given cell and +energy group. This distribution can be compared to the MC source distribution +to generate weight adjusted factors defined as + +.. math:: + :label: eq_waf + + f_{l,m,n}^g = \frac{Np_{l,m,n}^g}{\sum\limits_s w_s};\quad s\in + \left(g,l,m,n\right). + +The MC source distribution is represented on the same coarse mesh as +CMFD by summing all neutrons' weights, :math:`w_s`, in a given cell and +energy group. MC source weights can then be modified by this weight +adjustment factor so that it matches the CMFD solution on the coarse +mesh, + +.. math:: + :label: src_mod + + w^\prime_s = w_s\times f_{l,m,n}^g;\quad s\in \left(g,l,m,n\right). + +It should be noted that heterogeneous information about local coordinates and +energy remain constant throughout this modification process. + +------------------------ +Implementation in OpenMC +------------------------ + +The section describes how CMFD was implemented in OpenMC. Before the simulation +begins, a user sets up a CMFD input file that contains the following basic +information: + +* CMFD mesh (space and energy), +* boundary conditions at edge of mesh (albedos), +* acceleration region (subset of mesh, optional), +* fission source generation (FSG)/batch that CMFD should begin, and +* whether CMFD feedback should be applied. + +It should be noted that for more difficult simulations (e.g., light water +reactors), there are other options available to users such as tally resetting +parameters, effective down-scatter usage, tally estimator, etc. For more +information please see :ref:`usersguide_cmfd`. + +Of the options described above, the optional acceleration subset region is an +uncommon feature. Because OpenMC only has a structured Cartesian mesh, mesh +cells may overlay regions that don't contain fissionable material and may be so +far from the core that the neutron flux is very low. If these regions were +included in the CMFD solution, bad estimates of diffusion parameters may result +and affect CMFD feedback. To deal with this, a user can carve out an active +acceleration region from their structured Cartesian mesh. This is illustrated +in diagram below. When placing a CMFD mesh over a geometry, the boundary +conditions must be known at the global edges of the mesh. If the geometry is +complex like the one below, one may have to cover the whole geometry including +the reactor pressure vessel because we know that there is a zero incoming +current boundary condition at the outer edge of the pressure vessel. This is +not viable in practice because neutrons in simulations may not reach mesh cells +that are near the pressure vessel. To circumvent this, one can shrink the mesh +to cover just the core region as shown in the diagram. However, one must still +estimate the boundary conditions at the global boundaries, but at these +locations, they are not readily known. In OpenMC, one can carve out the active +core region from the entire structured Cartesian mesh. This is shown in the +diagram below by the darkened region over the core. The albedo boundary +conditions at the active core/reflector boundary can be tallied indirectly +during the MC simulation with incoming and outgoing partial currents. This +allows the user to not have to worry about neutrons producing adequate tallies +in mesh cells far away from the core. + +.. tikz:: Diagram of CMFD acceleration mesh + :libs: shapes, snakes, shadows, arrows, calc, decorations.markings, patterns, fit, matrix, spy + :include: cmfd_tikz/meshfig.tikz + +During an MC simulation, CMFD tallies are accumulated. The basic tallies needed +are listed in Table :ref:`tab_tally`. Each tally is performed on a spatial and +energy mesh basis. The surface area-integrated net current is tallied on every +surface of the mesh. OpenMC tally objects are created by the CMFD code +internally, and cross sections are calculated at each CMFD feedback iteration. +The first CMFD iteration, controlled by the user, occurs just after tallies are +communicated to the master processor. Once tallies are collapsed, cross +sections, diffusion coefficients and equivalence parameters are calculated. This +is performed only on the acceleration region if that option has been activated +by the user. Once all diffusion parameters are calculated, CMFD matrices are +formed where energy groups are the inner most iteration index. In OpenMC, +compressed row storage sparse matrices are used due to the sparsity of CMFD +operators. An example of this sparsity is shown for the 3-D BEAVRS model in +figures :ref:`fig_loss` and :ref:`fig_prod` [BEAVRS]_. These matrices represent +an assembly radial mesh, 24 cell mesh in the axial direction and two energy +groups. The loss matrix is 99.92% sparse and the production matrix is 99.99% +sparse. Although the loss matrix looks like it is tridiagonal, it is really a +seven banded matrix with a block diagonal matrix for scattering. The production +matrix is a :math:`2\times 2` block diagonal; however, zeros are present because +no fission neutrons appear with energies in the thermal group. + +.. _tab_tally: + +.. table:: OpenMC CMFD tally list + + +--------------------------------------------------------------------------------------------+----------------+---------------------------+ + +--------------------------------------------------------------------------------------------+----------------+---------------------------+ + | tally | score | filter | + +============================================================================================+================+===========================+ + | \ :math:`\left\langle\overline{\overline\phi}_{l,m,n}^g | flux | mesh, energy | + | \Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | | + +--------------------------------------------------------------------------------------------+----------------+---------------------------+ + | \ :math:`\left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g | total | mesh, energy | + | \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | | + +--------------------------------------------------------------------------------------------+----------------+---------------------------+ + | \ :math:`\left\langle\overline{\overline{\nu_s\Sigma}}_{s1_{l,m,n}}^g | nu-scatter-1 | mesh, energy | + | \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | | + +--------------------------------------------------------------------------------------------+----------------+---------------------------+ + | \ :math:`\left\langle\overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow g} | nu-scatter | mesh, energy, energyout | + | \overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | | + +--------------------------------------------------------------------------------------------+----------------+---------------------------+ + | \ :math:`\left\langle\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow g} | nu-fission | mesh, energy, energyout | + | \overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | | + +--------------------------------------------------------------------------------------------+----------------+---------------------------+ + | \ :math:`\left\langle\overline{J}^{u,g}_{l\pm 1/2,m,n}\Delta_m^v\Delta_n^w\right\rangle` | current | mesh, energy | + +--------------------------------------------------------------------------------------------+----------------+---------------------------+ + +.. _fig_loss: + +.. figure:: ../_images/loss.png + :scale: 50 + + Sparsity of Neutron Loss Operator + +.. _fig_prod: + +.. figure:: ../_images/prod.png + :scale: 50 + + Sparsity of Neutron Production Operator + +To solve the eigenvalue problem with these matrices, different source iteration +and linear solvers can be used. The most common source iteration solver used is +standard power iteration as described in [Gill]_. To accelerate these source +iterations, a Wielandt shift scheme can be used as discussed in [Park]_. PETSc +solvers were first implemented to perform the linear solution in parallel that +occurs once per source iteration. When using PETSc, different types of parallel +linear solvers and preconditioners can be used. By default, OpenMC uses an +incomplete LU preconditioner and a GMRES Krylov solver. After some initial +studies of parallelization with PETSc, it was observed that because CMFD +matrices are very sparse, solution times do not scale well. An additional +Gauss-Seidel linear solver with Chebyshev acceleration was added that is +similar to the one used for CMFD in CASMO [Rhodes]_ and [Smith]_. This solver +was implemented with a custom section for two energy groups. Because energy +group is the inner most index, a block diagonal is formed when using more than +one group. For two groups, it is easy to invert this diagonal analytically +inside the Gauss-Seidel iterative solver. For more than two groups, this +analytic inversion can still be performed, but with more computational effort. +A standard Gauss-Seidel solver is used for more than two groups. + +Besides a power iteration, a Jacobian-free Newton-Krylov method was also +implemented to obtain eigenvalue and multigroup fluxes as described in [Gill]_ +and [Knoll]_. This method is not the primary one used, but has gotten recent +attention due to its coupling advantages to other physics such as thermal +hydraulics. Once multigroup fluxes are obtained, a normalized fission source is +calculated in the code using eq. :eq:`eq_cmfd_psrc` directly. + +The next step in the process is to compute weight adjustment factors. These are +calculated by taking the ratio of the expected number of neutrons from the CMFD +source distribution to the current number of neutrons in each mesh. It is +straightforward to compute the CMFD number of neutrons because it is the +product between the total starting initial weight of neutrons and the CMFD +normalized fission source distribution. To compute the number of neutrons from +the current MC source, OpenMC sums the statistical +weights of neutrons from the source bank on a given spatial and energy mesh. +Once weight adjustment factors were calculated, each neutron's statistical +weight in the source bank was modified according to its location and energy. +Examples of CMFD simulations using OpenMC can be found in [Herman_Thesis]_. + +---------- +References +---------- + +.. [BEAVRS] Nick Horelik, Bryan Herman. *Benchmark for Evaluation And Verification of Reactor + Simulations*. Massachusetts Institute of Technology, http://crpg.mit.edu/pub/beavrs + , 2013. + +.. [Gill] Daniel F. Gill. *Newton-Krylov methods for the solution of the k-eigenvalue problem in + multigroup neutronics calculations*. Ph.D. thesis, Pennsylvania State University, 2010. + +.. [Hebert] Alain Hebert. *Applied reactor physics*. Presses Internationales Polytechnique, + Montreal, 2009. + +.. [Herman] Bryan R. Herman, Benoit Forget, Kord Smith, and Brian N. Aviles. Improved + diffusion coefficients generated from Monte Carlo codes. In *Proceedings of M&C + 2013*, Sun Valley, ID, USA, May 5 - 9, 2013. + +.. [Herman_Thesis] Bryan R. Herman. *Monte Carlo and Thermal Hydraulic Coupling using + Low-Order Nonlinear Diffusion Acceleration*. Sc.D. thesis, + Massachusetts Institute of Technology, 2014. + +.. [Knoll] D.A. Knoll, H. Park, and C. Newman. *Acceleration of k-eigenvalue/criticality + calculations using the Jacobian-free Newton-Krylov method*. Nuclear Science and + Engineering, 167:133–140, 2011. + +.. [Park] H. Park, D.A. Knoll, and C.K. Newman. *Nonlinear acceleration of transport + criticality problems*. Nuclear Science and Engineering, 172:52–65, 2012. + +.. [Rhodes] Joel Rhodes and Malte Edenius. *CASMO-4 --- A Fuel Assembly Burnup Program. + User’s Manual*. Studsvik of America, ssp-09/443-u rev 0, proprietary edition, 2001. + +.. [Smith] Kord S Smith and Joel D Rhodes III. *Full-core, 2-D, LWR core calculations with + CASMO-4E*. In Proceedings of PHYSOR 2002, Seoul, Korea, October 7 - 10, 2002. diff --git a/docs/source/methods/cmfd_tikz/cmfd_flow.tikz b/docs/source/methods/cmfd_tikz/cmfd_flow.tikz new file mode 100644 index 0000000000..0611918052 --- /dev/null +++ b/docs/source/methods/cmfd_tikz/cmfd_flow.tikz @@ -0,0 +1,19 @@ +\begin{tikzpicture} + \matrix[every node/.style={draw, thick, minimum width=3cm, minimum height=1cm, align=center}, column sep=2cm, row sep=1cm] (m) { +\node[draw, fill=red!40] (start) {Batch $i$ \\ tally NDA}; & \\ + \node[draw, diamond, aspect=2, fill=green!40] (cmfd) {Run NDA?}; & \node[draw, fill=red!40] (end) {Batch $i + 1$ \\ tally NDA}; \\ +\node[draw, fill=blue!40] (xs) {Calculate XS \& DC}; & \node[draw, fill=blue!40] (modify) {Modify MC Source}; \\ +\node[draw, fill=blue!40] (nonlinear) {Calculate Equivalence}; & \node[draw, fill=blue!40] (eqs) {Solve NDA eqs.};\\ +}; + +\begin{scope}[every path/.style={->,very thick,draw}] + \draw (start.south) -- (cmfd.north); + \draw (cmfd.east) -- node[above] {no} (end.west); + \draw (cmfd.south) -- node[right] {yes} (xs.north); + \draw (xs.south) -- (nonlinear.north); + \draw (nonlinear.east) -- (eqs.west); + \draw (eqs.north) -- (modify.south); + \draw (modify.north) -- (end.south); + \end{scope} + +\end{tikzpicture} \ No newline at end of file diff --git a/docs/source/methods/cmfd_tikz/meshfig.tikz b/docs/source/methods/cmfd_tikz/meshfig.tikz new file mode 100644 index 0000000000..1e5bbfa476 --- /dev/null +++ b/docs/source/methods/cmfd_tikz/meshfig.tikz @@ -0,0 +1,628 @@ + + % these dimensions are determined in arrow_dimms.ods + + \def\scale{1.0} + + \def\latWidth{0.2808363589*\scale} + + \def\RPVOR{3*\scale} + \def\rectW{0.75*\scale} + \def\RPVIR{2.8694005485*\scale} + \def\BarrelIR{2.4547472901*\scale} + \def\BarrelOR{2.5293848766*\scale} + \def\ShieldOR{2.6040224631*\scale} + + \def\bafCIRx{0.9829272561*\scale} + \def\bafCIRy{2.1062726917*\scale} + \def\bafCORx{1.0119529842*\scale} + \def\bafCORy{2.1352984197*\scale} + \def\bafMIRx{1.8254363328*\scale} + \def\bafMIRy{1.5445999739*\scale} + \def\bafMORx{1.8544620609*\scale} + \def\bafMORy{1.573625702*\scale} + + \tikzset{Assembly/.style={ + inner sep=0pt, + text width=\latWidth in, + minimum size=\latWidth in, + draw=black, + align=center + } + } + + \def\tkzRPV{(0,0) circle (\RPVIR) (0,0) circle (\RPVOR)} + \def\tkzBarrel{(0,0) circle (\BarrelIR) (0,0) circle (\BarrelOR)} + \def\tkzShields{(0,0) circle (\BarrelOR) (0,0) circle (\ShieldOR)} + + \def\tkzBaffCOR{(-\bafCORx, -\bafCORy) rectangle (\bafCORx, \bafCORy)} + \def\tkzBaffCIR{(-\bafCIRx, -\bafCIRy) rectangle (\bafCIRx, \bafCIRy)} + \def\tkzBaffMOR{(-\bafMORx, -\bafMORy) rectangle (\bafMORx, \bafMORy)} + \def\tkzBaffMIR{(-\bafMIRx, -\bafMIRy) rectangle (\bafMIRx, \bafMIRy) } + \def\tkzBaffleC{ \tkzBaffCIR \tkzBaffCOR } + \def\tkzBaffleM{ \tkzBaffMIR \tkzBaffMOR } + + \def\tkzBaffCClip{\tkzBaffCIR (-\RPVOR, -\RPVOR) rectangle (\RPVOR, \RPVOR)} + \def\tkzBaffMClip{\tkzBaffMIR (-\RPVOR, -\RPVOR) rectangle (\RPVOR, \RPVOR)} + + \def\highenr{blue!50} + \def\midenr{yellow!50} + \def\lowenr{red!50} + \def\lightgray{black!25} + \def\darkgray{black!80} + + \begin{tikzpicture}[x=1in,y=1in, xshift=3in] +\scalebox{0.6}{ + % draw RPV, barrel, and shield panels + + \path[fill=black,even odd rule] \tkzRPV; + \path[fill=black,even odd rule] \tkzBarrel; + \begin{scope} + \clip[rotate around={45:(0,0)}] (-\RPVOR, -\rectW) rectangle (\RPVOR, \rectW) (-\rectW, \RPVOR) rectangle (\rectW, -\RPVOR); + \path[fill=black,even odd rule] \tkzShields; + \end{scope} + + + % draw assembly row/column headers + + \draw[red, thick] ($(-7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {R} -- ($(-7*\latWidth,4*\latWidth)$); + \draw[red, thick] ($(-6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {P} -- ($(-6*\latWidth,6*\latWidth)$); + \draw[red, thick] ($(-5*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {N} -- ($(-5*\latWidth,7*\latWidth)$); + \draw[red, thick] ($(-4*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {M} -- ($(-4*\latWidth,7*\latWidth)$); + \draw[red, thick] ($(-3*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {L} -- ($(-3*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(-2*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {K} -- ($(-2*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(-1*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {J} -- ($(-1*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(-0*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {H} -- ($(-0*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(1*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {G} -- ($(1*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(2*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {F} -- ($(2*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(3*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {E} -- ($(3*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(4*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {D} -- ($(4*\latWidth,7*\latWidth)$); + \draw[red, thick] ($(5*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {C} -- ($(5*\latWidth,7*\latWidth)$); + \draw[red, thick] ($(6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {B} -- ($(6*\latWidth,6*\latWidth)$); + \draw[red, thick] ($(7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {A} -- ($(7*\latWidth,4*\latWidth)$); + + \begin{scope}[rotate=90] + \draw[red, thick] ($(-7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {15} -- ($(-7*\latWidth,4*\latWidth)$); + \draw[red, thick] ($(-6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {14} -- ($(-6*\latWidth,6*\latWidth)$); + \draw[red, thick] ($(-5*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {13} -- ($(-5*\latWidth,7*\latWidth)$); + \draw[red, thick] ($(-4*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {12} -- ($(-4*\latWidth,7*\latWidth)$); + \draw[red, thick] ($(-3*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {11} -- ($(-3*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(-2*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {10} -- ($(-2*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(-1*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {9} -- ($(-1*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(-0*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {8} -- ($(-0*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(1*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {7} -- ($(1*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(2*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {6} -- ($(2*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(3*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {5} -- ($(3*\latWidth,8*\latWidth)$); + \draw[red, thick] ($(4*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {4} -- ($(4*\latWidth,7*\latWidth)$); + \draw[red, thick] ($(5*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {3} -- ($(5*\latWidth,7*\latWidth)$); + \draw[red, thick] ($(6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {2} -- ($(6*\latWidth,6*\latWidth)$); + \draw[red, thick] ($(7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {1} -- ($(7*\latWidth,4*\latWidth)$); + \end{scope} + + % draw fuel assembly nodes + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-5*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-4*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-3*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-2*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-1*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-0*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 1*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 2*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 3*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 4*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 5*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,8*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-5*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-4*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-3*\latWidth,7*\latWidth)$) {}; % L1 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-2*\latWidth,7*\latWidth)$) {6}; % K1 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-1*\latWidth,7*\latWidth)$) {}; % J1 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-0*\latWidth,7*\latWidth)$) {6}; % H1 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 1*\latWidth,7*\latWidth)$) {}; % G1 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 2*\latWidth,7*\latWidth)$) {6}; % F1 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 3*\latWidth,7*\latWidth)$) {}; % E1 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 4*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 5*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,7*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-5*\latWidth,6*\latWidth)$) {}; % N2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-4*\latWidth,6*\latWidth)$) {}; % M2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-3*\latWidth,6*\latWidth)$) {16}; % L2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-2*\latWidth,6*\latWidth)$) {}; % K2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-1*\latWidth,6*\latWidth)$) {20}; % J2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-0*\latWidth,6*\latWidth)$) {}; % H2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 1*\latWidth,6*\latWidth)$) {20}; % G2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 2*\latWidth,6*\latWidth)$) {}; % F2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 3*\latWidth,6*\latWidth)$) {16}; % E2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 4*\latWidth,6*\latWidth)$) {}; % D2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 5*\latWidth,6*\latWidth)$) {}; % C2 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,6*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-6*\latWidth,5*\latWidth)$) {}; % P3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-5*\latWidth,5*\latWidth)$) {15}; % N3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-4*\latWidth,5*\latWidth)$) {16}; % M3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-3*\latWidth,5*\latWidth)$) {}; % L3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-2*\latWidth,5*\latWidth)$) {16}; % K3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-1*\latWidth,5*\latWidth)$) {}; % J3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-0*\latWidth,5*\latWidth)$) {16}; % H3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 1*\latWidth,5*\latWidth)$) {}; % G3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 2*\latWidth,5*\latWidth)$) {16}; % F3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 3*\latWidth,5*\latWidth)$) {}; % E3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 4*\latWidth,5*\latWidth)$) {16}; % D3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 5*\latWidth,5*\latWidth)$) {15}; % C3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 6*\latWidth,5*\latWidth)$) {}; % B3 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,5*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,5*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-6*\latWidth,4*\latWidth)$) {}; % P4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-5*\latWidth,4*\latWidth)$) {16}; % N4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-4*\latWidth,4*\latWidth)$) {}; % M4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-3*\latWidth,4*\latWidth)$) {16}; % L4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-2*\latWidth,4*\latWidth)$) {}; % K4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-1*\latWidth,4*\latWidth)$) {12}; % J4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-0*\latWidth,4*\latWidth)$) {}; % H4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 1*\latWidth,4*\latWidth)$) {12}; % G4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 2*\latWidth,4*\latWidth)$) {}; % F4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 3*\latWidth,4*\latWidth)$) {16}; % E4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 4*\latWidth,4*\latWidth)$) {}; % D4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 5*\latWidth,4*\latWidth)$) {16}; % C4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 6*\latWidth,4*\latWidth)$) {}; % B4 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,4*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,4*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-7*\latWidth,3*\latWidth)$) {}; % R5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-7*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-6*\latWidth,3*\latWidth)$) {16}; % P5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-5*\latWidth,3*\latWidth)$) {}; % N5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-4*\latWidth,3*\latWidth)$) {16}; % M5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-3*\latWidth,3*\latWidth)$) {}; % L5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-2*\latWidth,3*\latWidth)$) {12}; % K5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-1*\latWidth,3*\latWidth)$) {}; % J5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-0*\latWidth,3*\latWidth)$) {12}; % H5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 1*\latWidth,3*\latWidth)$) {}; % G5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 2*\latWidth,3*\latWidth)$) {12}; % F5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 3*\latWidth,3*\latWidth)$) {}; % E5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 4*\latWidth,3*\latWidth)$) {16}; % D5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 5*\latWidth,3*\latWidth)$) {}; % C5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 6*\latWidth,3*\latWidth)$) {16}; % B5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 7*\latWidth,3*\latWidth)$) {}; % A5 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 7*\latWidth,3*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,3*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-7*\latWidth,2*\latWidth)$) {6}; % R6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-7*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-6*\latWidth,2*\latWidth)$) {}; % P6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-5*\latWidth,2*\latWidth)$) {16}; % N6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-4*\latWidth,2*\latWidth)$) {}; % M6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-3*\latWidth,2*\latWidth)$) {12}; % L6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-2*\latWidth,2*\latWidth)$) {}; % K6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-1*\latWidth,2*\latWidth)$) {12}; % J6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-0*\latWidth,2*\latWidth)$) {}; % H6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 1*\latWidth,2*\latWidth)$) {12}; % G6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 2*\latWidth,2*\latWidth)$) {}; % F6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 3*\latWidth,2*\latWidth)$) {12}; % E6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 4*\latWidth,2*\latWidth)$) {}; % D6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 5*\latWidth,2*\latWidth)$) {16}; % C6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 6*\latWidth,2*\latWidth)$) {}; % B6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 7*\latWidth,2*\latWidth)$) {6}; % A6 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 7*\latWidth,2*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,2*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-7*\latWidth,1*\latWidth)$) {}; % R7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-7*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-6*\latWidth,1*\latWidth)$) {20}; % P7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-5*\latWidth,1*\latWidth)$) {}; % N7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-4*\latWidth,1*\latWidth)$) {12}; % M7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-3*\latWidth,1*\latWidth)$) {}; % L7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-2*\latWidth,1*\latWidth)$) {12}; % K7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-1*\latWidth,1*\latWidth)$) {}; % J7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-0*\latWidth,1*\latWidth)$) {16}; % H7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 1*\latWidth,1*\latWidth)$) {}; % G7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 2*\latWidth,1*\latWidth)$) {12}; % F7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 3*\latWidth,1*\latWidth)$) {}; % E7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 4*\latWidth,1*\latWidth)$) {12}; % D7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 5*\latWidth,1*\latWidth)$) {}; % C7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 6*\latWidth,1*\latWidth)$) {20}; % B7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 7*\latWidth,1*\latWidth)$) {}; % A7 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 7*\latWidth,1*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,1*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-7*\latWidth,0*\latWidth)$) {6}; % R8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-7*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-6*\latWidth,0*\latWidth)$) {}; % P8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-5*\latWidth,0*\latWidth)$) {16}; % N8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-4*\latWidth,0*\latWidth)$) {}; % M8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-3*\latWidth,0*\latWidth)$) {12}; % L8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-2*\latWidth,0*\latWidth)$) {}; % K8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-1*\latWidth,0*\latWidth)$) {16}; % J8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-0*\latWidth,0*\latWidth)$) {}; % H8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 1*\latWidth,0*\latWidth)$) {16}; % G8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 2*\latWidth,0*\latWidth)$) {}; % F8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 3*\latWidth,0*\latWidth)$) {12}; % E8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 4*\latWidth,0*\latWidth)$) {}; % D8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 5*\latWidth,0*\latWidth)$) {16}; % C8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 6*\latWidth,0*\latWidth)$) {}; % B8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 7*\latWidth,0*\latWidth)$) {6}; % A8 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 7*\latWidth,0*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,0*\latWidth)$) {}; + + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-7*\latWidth,-1*\latWidth)$) {}; % R9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-7*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-6*\latWidth,-1*\latWidth)$) {20}; % P9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-5*\latWidth,-1*\latWidth)$) {}; % N9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-4*\latWidth,-1*\latWidth)$) {12}; % M9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-3*\latWidth,-1*\latWidth)$) {}; % L9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-2*\latWidth,-1*\latWidth)$) {12}; % K9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-1*\latWidth,-1*\latWidth)$) {}; % J9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-0*\latWidth,-1*\latWidth)$) {16}; % H9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 1*\latWidth,-1*\latWidth)$) {}; % G9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 2*\latWidth,-1*\latWidth)$) {12}; % F9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 3*\latWidth,-1*\latWidth)$) {}; % E9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 4*\latWidth,-1*\latWidth)$) {12}; % D9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 5*\latWidth,-1*\latWidth)$) {}; % C9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 6*\latWidth,-1*\latWidth)$) {20}; % B9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 7*\latWidth,-1*\latWidth)$) {}; % A9 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 7*\latWidth,-1*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-1*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-7*\latWidth,-2*\latWidth)$) {6}; % R10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-7*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-6*\latWidth,-2*\latWidth)$) {}; % P10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-5*\latWidth,-2*\latWidth)$) {16}; % N10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-4*\latWidth,-2*\latWidth)$) {}; % M10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-3*\latWidth,-2*\latWidth)$) {12}; % L10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-2*\latWidth,-2*\latWidth)$) {}; % K10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-1*\latWidth,-2*\latWidth)$) {12}; % J10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-0*\latWidth,-2*\latWidth)$) {}; % H10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 1*\latWidth,-2*\latWidth)$) {12}; % G10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 2*\latWidth,-2*\latWidth)$) {}; % F10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 3*\latWidth,-2*\latWidth)$) {12}; % E10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 4*\latWidth,-2*\latWidth)$) {}; % D10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 5*\latWidth,-2*\latWidth)$) {16}; % C10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 6*\latWidth,-2*\latWidth)$) {}; % B10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 7*\latWidth,-2*\latWidth)$) {6}; % A10 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 7*\latWidth,-2*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-2*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-7*\latWidth,-3*\latWidth)$) {}; % R11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-7*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-6*\latWidth,-3*\latWidth)$) {16}; % P11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-5*\latWidth,-3*\latWidth)$) {}; % N11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-4*\latWidth,-3*\latWidth)$) {16}; % M11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-3*\latWidth,-3*\latWidth)$) {}; % L11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-2*\latWidth,-3*\latWidth)$) {12}; % K11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-1*\latWidth,-3*\latWidth)$) {}; % J11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-0*\latWidth,-3*\latWidth)$) {12}; % H11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 1*\latWidth,-3*\latWidth)$) {}; % G11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 2*\latWidth,-3*\latWidth)$) {12}; % F11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 3*\latWidth,-3*\latWidth)$) {}; % E11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 4*\latWidth,-3*\latWidth)$) {16}; % D11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 5*\latWidth,-3*\latWidth)$) {}; % C11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 6*\latWidth,-3*\latWidth)$) {16}; % B11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 7*\latWidth,-3*\latWidth)$) {}; % A11 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 7*\latWidth,-3*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-3*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-6*\latWidth,-4*\latWidth)$) {}; % P12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-5*\latWidth,-4*\latWidth)$) {16}; % N12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-4*\latWidth,-4*\latWidth)$) {}; % M12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-3*\latWidth,-4*\latWidth)$) {16}; % L12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-2*\latWidth,-4*\latWidth)$) {}; % K12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-1*\latWidth,-4*\latWidth)$) {12}; % J12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-0*\latWidth,-4*\latWidth)$) {}; % H12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 1*\latWidth,-4*\latWidth)$) {12}; % G12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 2*\latWidth,-4*\latWidth)$) {}; % F12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 3*\latWidth,-4*\latWidth)$) {16}; % E12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 4*\latWidth,-4*\latWidth)$) {}; % D12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 5*\latWidth,-4*\latWidth)$) {16}; % C12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 6*\latWidth,-4*\latWidth)$) {}; % B12 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,-4*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-4*\latWidth)$) {}; + + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-6*\latWidth,-5*\latWidth)$) {}; % P13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-6*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-5*\latWidth,-5*\latWidth)$) {15}; % N13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-4*\latWidth,-5*\latWidth)$) {16}; % M13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-3*\latWidth,-5*\latWidth)$) {}; % L13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-2*\latWidth,-5*\latWidth)$) {16}; % K13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-1*\latWidth,-5*\latWidth)$) {}; % J13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($(-0*\latWidth,-5*\latWidth)$) {16}; % H13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 1*\latWidth,-5*\latWidth)$) {}; % G13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 2*\latWidth,-5*\latWidth)$) {16}; % F13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 3*\latWidth,-5*\latWidth)$) {}; % E13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\midenr] at ($( 4*\latWidth,-5*\latWidth)$) {16}; % D13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 5*\latWidth,-5*\latWidth)$) {15}; % C13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 6*\latWidth,-5*\latWidth)$) {}; % B13 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 6*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,-5*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-5*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-5*\latWidth,-6*\latWidth)$) {}; % N14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-5*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-4*\latWidth,-6*\latWidth)$) {}; % M14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-4*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-3*\latWidth,-6*\latWidth)$) {16}; % L14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-2*\latWidth,-6*\latWidth)$) {}; % K14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-1*\latWidth,-6*\latWidth)$) {20}; % J14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($(-0*\latWidth,-6*\latWidth)$) {}; % H14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 1*\latWidth,-6*\latWidth)$) {20}; % G14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\lowenr] at ($( 2*\latWidth,-6*\latWidth)$) {}; % F14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 3*\latWidth,-6*\latWidth)$) {16}; % E14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 4*\latWidth,-6*\latWidth)$) {}; % D14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 4*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 5*\latWidth,-6*\latWidth)$) {}; % C14 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 5*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,-6*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-6*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-5*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-4*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-3*\latWidth,-7*\latWidth)$) {}; % L15 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-3*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-2*\latWidth,-7*\latWidth)$) {6}; % K15 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-2*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-1*\latWidth,-7*\latWidth)$) {}; % J15 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-1*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($(-0*\latWidth,-7*\latWidth)$) {6}; % H15 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($(-0*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 1*\latWidth,-7*\latWidth)$) {}; % G15 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 1*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 2*\latWidth,-7*\latWidth)$) {6}; % F15 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 2*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\highenr] at ($( 3*\latWidth,-7*\latWidth)$) {}; % E15 + \node [Assembly, fill=\darkgray, opacity=0.7] at ($( 3*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 4*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 5*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,-7*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-7*\latWidth)$) {}; + + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-5*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-4*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-3*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-2*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-1*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-0*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 1*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 2*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 3*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 4*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 5*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,-8*\latWidth)$) {}; + \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-8*\latWidth)$) {}; + + % draw baffle north/south + + \begin{scope}[even odd rule] + \clip[rotate=90] \tkzBaffMClip; + \path[fill=black] \tkzBaffleC; + \end{scope} + \begin{scope}[even odd rule] + \clip \tkzBaffCClip; + \clip \tkzBaffMClip; + \path[fill=black, rotate=90] \tkzBaffleM; + \end{scope} + + % draw baffle east/west + + \begin{scope}[rotate=90] + \begin{scope}[even odd rule] + \clip[rotate=90] \tkzBaffMClip; + \path[fill=black] \tkzBaffleC; + \end{scope} + \begin{scope}[even odd rule] + \clip \tkzBaffCClip; + \clip \tkzBaffMClip; + \path[fill=black, rotate=90] \tkzBaffleM; + \end{scope} + \end{scope}} + \end{tikzpicture} diff --git a/docs/source/methods/index.rst b/docs/source/methods/index.rst index 0ff4fb1988..1df4f324a3 100644 --- a/docs/source/methods/index.rst +++ b/docs/source/methods/index.rst @@ -16,3 +16,4 @@ Theory and Methodology tallies eigenvalue parallelization + cmfd diff --git a/docs/source/releasenotes/index.rst b/docs/source/releasenotes/index.rst index b40126b243..9799a5cfc7 100644 --- a/docs/source/releasenotes/index.rst +++ b/docs/source/releasenotes/index.rst @@ -10,6 +10,7 @@ bugs fixed, and known issues for each successive release. .. toctree:: :maxdepth: 1 + notes_0.6.1 notes_0.6.0 notes_0.5.4 notes_0.5.3 diff --git a/docs/source/releasenotes/notes_0.6.1.rst b/docs/source/releasenotes/notes_0.6.1.rst new file mode 100644 index 0000000000..b0651a8629 --- /dev/null +++ b/docs/source/releasenotes/notes_0.6.1.rst @@ -0,0 +1,63 @@ +.. _notes_0.6.1: + +============================== +Release Notes for OpenMC 0.6.1 +============================== + +------------------- +System Requirements +------------------- + +There are no special requirements for running the OpenMC code. As of this +release, OpenMC has been tested on a variety of Linux distributions, Mac OS X, +and Microsoft Windows 7. Memory requirements will vary depending on the size of +the problem at hand (mostly on the number of nuclides in the problem). + +------------ +New Features +------------ + +- Coarse mesh finite difference acceleration no longer requires PETSc +- Statepoint file numbering is now zero-padded +- Python scripts now compatible with Python 2 or 3 +- Ability to run particle restarts in fixed source calculations +- Capability to filter box source by fissionable materials +- Nuclide/element names are now case insensitive in input files +- Improved treatment of resonance scattering for heavy nuclides + +--------- +Bug Fixes +--------- + +- 03e890_: Check for energy-dependent multiplicities in ACE files +- 4439de_: Fix distance-to-surface calculation for general plane surface +- 5808ed_: Account for differences in URR band probabilities at different energies +- 2e60c0_: Allow zero atom/weight percents in materials +- 3e0870_: Don't use PWD environment variable when setting path to input files +- dc4776_: Handle probability table resampling correctly +- 01178b_: Fix metastables nuclides in NNDC cross_sections.xml file +- 62ec43_: Don't read tallies.xml when OpenMC is run in plotting mode +- 2a95ef_: Prevent segmentation fault on "current" score without mesh filter + +.. _03e890: https://github.com/mit-crpg/openmc/commit/03e890 +.. _4439de: https://github.com/mit-crpg/openmc/commit/4439de +.. _5808ed: https://github.com/mit-crpg/openmc/commit/5808ed +.. _2e60c0: https://github.com/mit-crpg/openmc/commit/2e60c0 +.. _3e0870: https://github.com/mit-crpg/openmc/commit/3e0870 +.. _dc4776: https://github.com/mit-crpg/openmc/commit/dc4776 +.. _01178b: https://github.com/mit-crpg/openmc/commit/01178b +.. _62ec43: https://github.com/mit-crpg/openmc/commit/62ec43 +.. _2a95ef: https://github.com/mit-crpg/openmc/commit/2a95ef + +------------ +Contributors +------------ + +This release contains new contributions from the following people: + +- `Sterling Harper `_ +- `Bryan Herman `_ +- `Adam Nelson `_ +- `Paul Romano `_ +- `Jon Walsh `_ +- `Will Boyd `_ diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index c97ab8ac1e..09aea601bb 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -261,7 +261,7 @@ or sub-elements and can be set to either "false" or "true". *Default*: true ```` Element ----------------------- +---------------------------------- The ``resonance_scattering`` element can contain one or more of the following attributes or sub-elements: @@ -269,7 +269,7 @@ attributes or sub-elements: :scatterer: An element with attributes/sub-elements called ``nuclide``, ``method``, ``xs_label``, ``xs_label_0K``, ``E_min``, and ``E_max``. The ``nuclide`` - attribute is the name, as given by the ``name`` attribute within the + attribute is the name, as given by the ``name`` attribute within the ``nuclide`` sub-element of the ``material`` element in ``materials.xml``, of the nuclide to which a resonance scattering treatment is to be applied. The ``method`` attribute gives the type of resonance scattering treatment @@ -433,6 +433,13 @@ attributes/sub-elements: *Default*: 0.988 2.249 + :write_initial: + An element specifying whether to write out the initial source bank used at + the beginning of the first batch. The output file is named + "initial_source.binary(h5)" + + *Default*: false + ```` Element ------------------------- @@ -1323,6 +1330,8 @@ attributes or sub-elements. These are not used in "voxel" plots: *Default*: None +.. _usersguide_cmfd: + ------------------------------ CMFD Specification -- cmfd.xml ------------------------------ @@ -1332,15 +1341,6 @@ Currently, it allows users to accelerate fission source convergence during inactive neutron batches. To run CMFD, the ```` element in ``settings.xml`` should be set to "true". -```` Element --------------------------- - -The ```` element controls the batch where CMFD tallies should be -reset. CMFD tallies should be reset before active batches so they are accumulated -without bias. - - *Default*: 0 - ```` Element ------------------- @@ -1362,7 +1362,25 @@ The ```` element sets one additional CMFD output column. Options are: * "source" - prints the RMS [%] between the OpenMC fission source and CMFD fission source. - *Default*: None + *Default*: balance + +```` Element +------------------------ + +The ```` element controls whether :math:`\widehat{D}` nonlinear +CMFD parameters should be reset to zero before solving CMFD eigenproblem. +It can be turned on with "true" and off with "false". + + *Default*: false + +```` Element +------------------------- + +The ```` element controls whether an effective downscatter cross +section should be used when using 2-group CMFD. It can be turned on with "true" +and off with "false". + + *Default*: false ```` Element ---------------------- @@ -1373,24 +1391,16 @@ It can be turned on with "true" and off with "false". *Default*: false -```` Element ----------------------- +```` Element +------------------------------------ -The ```` element controls if cmfd tallies should be accumulated -during inactive batches. For some applications, CMFD tallies may not be -needed until the start of active batches. This option can be turned on -with "true" and off with "false" +The ```` element specifies two parameters. The first is +the absolute inner tolerance for Gauss-Seidel iterations when performing CMFD +and the second is the relative inner tolerance for Gauss-Seidel iterations +for CMFD calculations. It is only used in the standalone CMFD power iteration +solver and not when PETSc is active. - *Default*: true - -```` Element ----------------------------- - -The ```` element controls when CMFD tallies are reset during -inactive batches. The integer set here is the interval at which this reset -occurs. The amout of resets is controlled with the ```` element. - - *Defualt*: 9999 + *Default*: 1.e-10 1.e-5 ```` Element ------------------------- @@ -1399,9 +1409,16 @@ The ```` element is used to view the convergence of linear GMRES iterations in PETSc. This option can be turned on with "true" and turned off with "false". - *Default*: false +```` Element +-------------------- + +The ```` element specifies the tolerance on the eigenvalue when performing +CMFD power iteration. + + *Default*: 1.e-8 + ```` Element ------------------ @@ -1470,14 +1487,6 @@ not impact the calculation. *Default*: 1.0 -```` Element -------------------------- - -The ```` element controls the number of CMFD tally resets that -occur during inactive CMFD batches. - - *Default*: 9999 - ```` Element --------------------------- @@ -1490,16 +1499,8 @@ This option can be turned on with "true" and turned off with "false". ------------------------- The ```` element can be turned on with "true" to have an adjoint -calculation be performed on the last batch when CMFD is active. - - *Default*: false - -```` Element --------------------------- - -The ```` element is used to view the convergence of the nonlinear SNES -function in PETSc. This option can be turned on with "true" and turned off with "false". - +calculation be performed on the last batch when CMFD is active. OpenMC should be +compiled with PETSc when using this option. *Default*: false @@ -1512,6 +1513,41 @@ By setting "power", power iteration is used and by setting "jfnk", JFNK is used. *Default*: power +```` Element +-------------------- + +The ```` element specifies an optional Wielandt shift parameter for +accelerating power iterations. It can only be used when PETSc is not active. +It is by default very large so the impact of the shift is effectively zero. + + *Default*: 1e6 + +```` Element +---------------------- + +The ```` element specifies an optional spectral radius that can be set to +accelerate the convergence of Gauss-Seidel iterations during CMFD power iteration +solve. Note this is only used in the standalone CMFD solver and does not affect +the calculation when PETSc is active. + + *Default*: power + +```` Element +------------------ + +The ```` element specifies the tolerance on the fission source when performing +CMFD power iteration. + + *Default*: 1.e-8 + +```` Element +------------------------- + +The ```` element contains a list of batch numbers in which CMFD tallies +should be reset. + + *Default*: None + ```` Element ---------------------------- diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fc6c6630cc..8a9424e79d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -277,7 +277,7 @@ file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_*.py) # Check to see if PETSC is compiled for CMFD tests if (NOT ${PETSC_ENABLED}) - file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd*.py) + file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd_jfnk.py) foreach(cmfd_test in ${CMFD_TESTS}) list(REMOVE_ITEM TESTS ${cmfd_test}) endforeach(cmfd_test) diff --git a/src/ace.F90 b/src/ace.F90 index fb42180da0..6a712ccd74 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -7,11 +7,11 @@ module ace use error, only: fatal_error, warning use fission, only: nu_total use global - use list_header, only: ListElemInt, ListInt + use list_header, only: ListInt use material_header, only: Material use output, only: write_message use set_header, only: SetChar - use string, only: to_str + use string, only: to_str, to_lower implicit none @@ -68,8 +68,8 @@ contains name = mat % names(j) if (.not. already_read % contains(name)) then - i_listing = xs_listing_dict % get_key(name) - i_nuclide = nuclide_dict % get_key(name) + i_listing = xs_listing_dict % get_key(to_lower(name)) + i_nuclide = nuclide_dict % get_key(to_lower(name)) name = xs_listings(i_listing) % name alias = xs_listings(i_listing) % alias @@ -116,8 +116,8 @@ contains name = mat % sab_names(k) if (.not. already_read % contains(name)) then - i_listing = xs_listing_dict % get_key(name) - i_sab = sab_dict % get_key(name) + i_listing = xs_listing_dict % get_key(to_lower(name)) + i_sab = sab_dict % get_key(to_lower(name)) ! Read the ACE table into the appropriate entry on the sab_tables ! array @@ -1563,16 +1563,11 @@ contains integer :: i ! index in nuclides array integer :: j ! index in nuclides array - type(ListElemInt), pointer :: nuc_list => null() ! pointer to nuclide list do i = 1, n_nuclides_total - allocate(nuclides(i) % nuc_list) - nuc_list => nuclides(i) % nuc_list do j = 1, n_nuclides_total if (nuclides(i) % zaid == nuclides(j) % zaid) then - nuc_list % data = j - allocate(nuc_list % next) - nuc_list => nuc_list % next + call nuclides(i) % nuc_list % append(j) end if end do end do diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 6d339891c9..76492dfe80 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -2,7 +2,7 @@ module ace_header use constants, only: MAX_FILE_LEN use endf_header, only: Tab1 - use list_header, only: ListElemInt + use list_header, only: ListInt implicit none @@ -97,7 +97,7 @@ module ace_header real(8) :: kT ! temperature in MeV (k*T) ! Linked list of indices in nuclides array of instances of this same nuclide - type(ListElemInt), pointer :: nuc_list => null() + type(ListInt) :: nuc_list ! Energy grid information integer :: n_grid ! # of nuclide grid points @@ -114,7 +114,7 @@ module ace_header ! Resonance scattering info logical :: resonant = .false. ! resonant scatterer? - character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c + character(10) :: name_0K = '' ! name of 0K nuclide, e.g. 92235.00c character(16) :: scheme ! target velocity sampling scheme integer :: n_grid_0K ! number of 0K energy grid points real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs @@ -416,6 +416,8 @@ module ace_header deallocate(this % reactions) end if + call this % nuc_list % clear() + end subroutine nuclide_clear end module ace_header diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index b5612dcd05..abba2d47d7 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -10,8 +10,6 @@ module cmfd_data private public :: set_up_cmfd, neutron_balance - logical :: dhat_reset = .false. - contains !============================================================================== @@ -103,6 +101,8 @@ contains cmfd % hxyz(2,:,:,:) = m % width(2) ! set y width cmfd % hxyz(3,:,:,:) = m % width(3) ! set z width + cmfd % keff_bal = ZERO + ! Begin loop around tallies TAL: do ital = 1, n_cmfd_tallies @@ -211,6 +211,9 @@ contains ! Bank source cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + & t % results(2,score_index) % sum + cmfd % keff_bal = cmfd % keff_bal + & + t % results(2,score_index) % sum / & + dble(t % n_realizations) end do INGROUP @@ -623,7 +626,9 @@ contains subroutine compute_dhat() use constants, only: CMFD_NOACCEL, ZERO - use global, only: cmfd, cmfd_coremap + use global, only: cmfd, cmfd_coremap, message, dhat_reset + use output, only: write_message + use string, only: to_str integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction @@ -743,7 +748,9 @@ contains cmfd%dhat(l,g,i,j,k) = dhat ! check for dhat reset - if (dhat_reset) cmfd%dhat(l,g,i,j,k) = ZERO + if (dhat_reset) then + cmfd%dhat(l,g,i,j,k) = ZERO + end if end do LEAK @@ -755,6 +762,12 @@ contains end do ZLOOP + ! write that dhats are zero + if (dhat_reset) then + message = 'Dhats reset to zero.' + call write_message(1) + end if + end subroutine compute_dhat !=============================================================================== @@ -763,8 +776,8 @@ contains function get_reflector_albedo(l, g, i, j, k) - use constants, only: ALBEDO_REJECT - use global, only: cmfd, cmfd_hold_weights + use constants, only: ONE + use global, only: cmfd real(8) :: get_reflector_albedo ! reflector albedo integer, intent(in) :: i ! iteration counter for x @@ -786,8 +799,7 @@ contains ! Calculate albedo if ((shift_idx == 1 .and. current(2*l ) < 1.0e-10_8) .or. & (shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then - albedo = ALBEDO_REJECT - cmfd_hold_weights = .true. + albedo = ONE else albedo = (current(2*l-1)/current(2*l))**(shift_idx) end if @@ -797,137 +809,6 @@ contains end function get_reflector_albedo -!=============================================================================== -! FIX_NEUTRON_BALANCE is a method to adjust parameters to have perfect balance -!=============================================================================== -#ifdef DEVELOPMENTAL - subroutine fix_neutron_balance() - - use constants, only: ONE, ZERO, CMFD_NOACCEL - use global, only: cmfd, keff - use, intrinsic :: ISO_FORTRAN_ENV - - integer :: nx ! number of mesh cells in x direction - integer :: ny ! number of mesh cells in y direction - integer :: nz ! number of mesh cells in z direction - integer :: ng ! number of energy groups - integer :: i ! iteration counter for x - integer :: j ! iteration counter for y - integer :: k ! iteration counter for z - integer :: l ! iteration counter for surface - real(8) :: leak1 ! leakage rate in group 1 - real(8) :: leak2 ! leakage rate in group 2 - real(8) :: flux1 ! group 1 volume int flux - real(8) :: flux2 ! group 2 volume int flux - real(8) :: sigt1 ! group 1 total xs - real(8) :: sigt2 ! group 2 total xs - real(8) :: sigs11 ! scattering transfer 1 --> 1 - real(8) :: sigs21 ! scattering transfer 2 --> 1 - real(8) :: sigs12 ! scattering transfer 1 --> 2 - real(8) :: sigs22 ! scattering transfer 2 --> 2 - real(8) :: nsigf11 ! fission transfer 1 --> 1 - real(8) :: nsigf21 ! fission transfer 2 --> 1 - real(8) :: nsigf12 ! fission transfer 1 --> 2 - real(8) :: nsigf22 ! fission transfer 2 --> 2 - real(8) :: siga1 ! group 1 abs xs - real(8) :: siga2 ! group 2 abs xs - real(8) :: sigs12_eff ! effective downscatter xs - - ! Extract spatial and energy indices from object - nx = cmfd % indices(1) - ny = cmfd % indices(2) - nz = cmfd % indices(3) - ng = cmfd % indices(4) - - ! Return if not two groups - if (ng /= 2) return - - ! Begin loop around space and energy groups - ZLOOP: do k = 1, nz - - YLOOP: do j = 1, ny - - XLOOP: do i = 1, nx - - ! Check for active mesh - if (allocated(cmfd%coremap)) then - if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) cycle - end if - - ! Compute leakage in groups 1 and 2 - leak1 = ZERO - leak2 = ZERO - LEAK: do l = 1, 3 - - leak1 = leak1 + ((cmfd % current(4*l,1,i,j,k) - & - cmfd % current(4*l-1,1,i,j,k))) - & - ((cmfd % current(4*l-2,1,i,j,k) - & - cmfd % current(4*l-3,1,i,j,k))) - - leak2 = leak2 + ((cmfd % current(4*l,2,i,j,k) - & - cmfd % current(4*l-1,2,i,j,k))) - & - ((cmfd % current(4*l-2,2,i,j,k) - & - cmfd % current(4*l-3,2,i,j,k))) - - - end do LEAK - - ! Extract cross sections and flux from object - flux1 = cmfd % flux(1,i,j,k) - flux2 = cmfd % flux(2,i,j,k) - sigt1 = cmfd % totalxs(1,i,j,k) - sigt2 = cmfd % totalxs(2,i,j,k) - sigs11 = cmfd % scattxs(1,1,i,j,k) - sigs21 = cmfd % scattxs(2,1,i,j,k) - sigs12 = cmfd % scattxs(1,2,i,j,k) - sigs22 = cmfd % scattxs(2,2,i,j,k) - nsigf11 = cmfd % nfissxs(1,1,i,j,k) - nsigf21 = cmfd % nfissxs(2,1,i,j,k) - nsigf12 = cmfd % nfissxs(1,2,i,j,k) - nsigf22 = cmfd % nfissxs(2,2,i,j,k) - - ! Check for no fission into group 2 - if (.not.(nsigf12 < 1e-6_8 .and. nsigf22 < 1e-6_8)) then - write(OUTPUT_UNIT,'(A,1PE11.4,1X,1PE11.4)') 'Fission in G=2', & - nsigf12,nsigf22 - end if - - ! Compute absorption xs - siga1 = sigt1 - sigs11 - sigs12 - siga2 = sigt2 - sigs22 - sigs21 - - ! Compute effective downscatter xs - sigs12_eff = (ONE/keff*nsigf11*flux1 - leak1 - siga1*flux1 & - - ONE/keff*nsigf21/siga2*leak2 ) / ( flux1*(ONE & - - ONE/keff*nsigf21/siga2)) - - ! Redefine flux 2 - flux2 = (sigs12_eff*flux1 - leak2)/siga2 - cmfd % flux(2,i,j,k) = flux2 - - ! Recompute total cross sections (use effective and no upscattering) - sigt1 = siga1 + sigs11 + sigs12_eff - sigt2 = siga2 + sigs22 - - ! Record total xs - cmfd % totalxs(1,i,j,k) = sigt1 - cmfd % totalxs(2,i,j,k) = sigt2 - - ! Record effective downscatter xs - cmfd % scattxs(1,2,i,j,k) = sigs12_eff - - ! Zero out upscatter cross section - cmfd % scattxs(2,1,i,j,k) = ZERO - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - end subroutine fix_neutron_balance -#endif - !=============================================================================== ! COMPUTE_EFFECTIVE_DOWNSCATTER changes downscatter rate for zero upscatter !=============================================================================== diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 6234d9fcf9..d14a3af243 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -22,6 +22,7 @@ contains use cmfd_data, only: set_up_cmfd use cmfd_power_solver, only: cmfd_power_execute use cmfd_jfnk_solver, only: cmfd_jfnk_execute + use cmfd_solver, only: cmfd_solver_execute use error, only: warning, fatal_error ! CMFD single processor on master @@ -37,6 +38,7 @@ contains call process_cmfd_options() ! Call solver +#ifdef PETSC if (trim(cmfd_solver_type) == 'power') then call cmfd_power_execute() elseif (trim(cmfd_solver_type) == 'jfnk') then @@ -45,6 +47,9 @@ contains message = 'solver type became invalid after input processing' call fatal_error() end if +#else + call cmfd_solver_execute() +#endif ! Save k-effective cmfd % k_cmfd(current_batch) = cmfd % keff @@ -64,7 +69,7 @@ contains call calc_fission_source() ! calculate weight factors - if (cmfd_feedback) call cmfd_reweight(.true.) + call cmfd_reweight(.true.) ! stop cmfd timer if (master) call time_cmfd % stop() @@ -77,36 +82,23 @@ contains subroutine cmfd_init_batch() - use global, only: cmfd_begin, cmfd_on, cmfd_tally_on, & - cmfd_inact_flush, cmfd_act_flush, cmfd_run, & - current_batch, cmfd_hold_weights + use global, only: cmfd_begin, cmfd_on, & + cmfd_reset, cmfd_run, & + current_batch ! Check to activate CMFD diffusion and possible feedback ! this guarantees that when cmfd begins at least one batch of tallies are ! accumulated if (cmfd_run .and. cmfd_begin == current_batch) then cmfd_on = .true. - cmfd_tally_on = .true. end if ! If this is a restart run and we are just replaying batches leave if (restart_run .and. current_batch <= restart_batch) return - ! Check to flush cmfd tallies for active batches, no more inactive flush - if (cmfd_run .and. cmfd_act_flush == current_batch) then + ! Check to reset tallies + if (cmfd_run .and. cmfd_reset % contains(current_batch)) then call cmfd_tally_reset() - cmfd_tally_on = .true. - cmfd_inact_flush(2) = -1 - end if - - ! Check to flush cmfd tallies during inactive batches (>= on number of - ! flushes important as the code will flush on the first batch which we - ! dont want to count) - if (cmfd_run .and. mod(current_batch,cmfd_inact_flush(1)) & - == 0 .and. cmfd_inact_flush(2) > 0 .and. cmfd_begin < current_batch) then - cmfd_hold_weights = .true. - call cmfd_tally_reset() - cmfd_inact_flush(2) = cmfd_inact_flush(2) - 1 end if end subroutine cmfd_init_batch @@ -139,6 +131,7 @@ contains use constants, only: CMFD_NOACCEL, ZERO, TWO use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch + use string, only: to_str #ifdef MPI use global, only: mpi_err @@ -267,6 +260,7 @@ contains use mesh_header, only: StructuredMesh use mesh, only: count_bank_sites, get_mesh_indices use search, only: binary_search + use string, only: to_str #ifdef MPI use global, only: mpi_err @@ -287,7 +281,6 @@ contains logical :: in_mesh ! source site is inside mesh type(StructuredMesh), pointer :: m ! point to mesh - real(8), allocatable :: egrid(:) ! energy grid ! Associate pointer m => meshes(n_user_meshes + 1) @@ -308,19 +301,16 @@ contains cmfd % weightfactors = ONE end if - ! Allocate energy grid and reverse cmfd energy grid - if (.not. allocated(egrid)) allocate(egrid(ng + 1)) - egrid = (/(cmfd % egrid(ng - i + 2), i = 1, ng + 1)/) - ! Compute new weight factors if (new_weights) then - ! Zero out weights - cmfd%weightfactors = ZERO + ! Set weight factors to a default 1.0 + cmfd%weightfactors = ONE - ! Count bank sites in mesh - call count_bank_sites(m, source_bank, cmfd%sourcecounts, egrid, & + ! Count bank sites in mesh and reverse due to egrid structure + call count_bank_sites(m, source_bank, cmfd%sourcecounts, cmfd % egrid, & sites_outside=outside, size_bank=work) + cmfd % sourcecounts = cmfd%sourcecounts(ng:1:-1,:,:,:) ! Check for sites outside of the mesh if (master .and. outside) then @@ -336,12 +326,14 @@ contains end where end if + if (.not. cmfd_feedback) return + ! Broadcast weight factors to all procs #ifdef MPI call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, & MPI_COMM_WORLD, mpi_err) #endif - end if + end if ! begin loop over source bank do i = 1, int(work,4) @@ -378,9 +370,6 @@ contains end do - ! Deallocate all - if (allocated(egrid)) deallocate(egrid) - end subroutine cmfd_reweight !=============================================================================== diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index 73bc0e0163..c9bc0f32af 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -83,6 +83,9 @@ module cmfd_header ! List of CMFD k real(8), allocatable :: k_cmfd(:) + ! Balance keff + real(8) :: keff_bal + end type cmfd_type contains diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 7e6ca18e53..67e8f0d0f3 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -62,16 +62,20 @@ contains use error, only: fatal_error, warning use global use output, only: write_message - use string, only: lower_case + use string, only: to_lower use xml_interface use, intrinsic :: ISO_FORTRAN_ENV + integer :: i integer :: ng + integer :: n_params integer, allocatable :: iarray(:) + integer, allocatable :: int_array(:) logical :: file_exists ! does cmfd.xml exist? logical :: found character(MAX_LINE_LEN) :: filename character(MAX_LINE_LEN) :: temp_str + real(8) :: gs_tol(2) type(Node), pointer :: doc => null() type(Node), pointer :: node_mesh => null() @@ -151,91 +155,121 @@ contains ! Set feedback logical if (check_for_node(doc, "feedback")) then call get_node_value(doc, "feedback", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_feedback = .true. + cmfd_feedback = .true. end if ! Set downscatter logical if (check_for_node(doc, "downscatter")) then call get_node_value(doc, "downscatter", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_downscatter = .true. + cmfd_downscatter = .true. + end if + + ! Reset dhat parameters + if (check_for_node(doc, "dhat_reset")) then + call get_node_value(doc, "dhat_reset", temp_str) + temp_str = to_lower(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & + dhat_reset = .true. end if ! Set the solver type if (check_for_node(doc, "solver")) & - call get_node_value(doc, "solver", cmfd_solver_type) + call get_node_value(doc, "solver", cmfd_solver_type) ! Set monitoring if (check_for_node(doc, "snes_monitor")) then call get_node_value(doc, "snes_monitor", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_snes_monitor = .true. + cmfd_snes_monitor = .true. end if if (check_for_node(doc, "ksp_monitor")) then call get_node_value(doc, "ksp_monitor", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_ksp_monitor = .true. + cmfd_ksp_monitor = .true. end if if (check_for_node(doc, "power_monitor")) then call get_node_value(doc, "power_monitor", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_power_monitor = .true. + cmfd_power_monitor = .true. end if ! Output logicals if (check_for_node(doc, "write_matrices")) then - call get_node_value(doc, "write_matices", temp_str) - call lower_case(temp_str) + call get_node_value(doc, "write_matrices", temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_write_matrices = .true. + cmfd_write_matrices = .true. end if ! Run an adjoint calc if (check_for_node(doc, "run_adjoint")) then call get_node_value(doc, "run_adjoint", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & +#ifndef PETSC + message = 'Must use PETSc when running adjoint option.' + call fatal_error() +#endif cmfd_run_adjoint = .true. end if ! Batch to begin cmfd if (check_for_node(doc, "begin")) & - call get_node_value(doc, "begin", cmfd_begin) + call get_node_value(doc, "begin", cmfd_begin) - ! Tally during inactive batches - if (check_for_node(doc, "inactive")) then - call get_node_value(doc, "inactive", temp_str) - call lower_case(temp_str) - if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') & - cmfd_tally_on = .false. + ! Check for cmfd tally resets + if (check_for_node(doc, "tally_reset")) then + n_cmfd_resets = get_arraysize_integer(doc, "tally_reset") + else + n_cmfd_resets = 0 + end if + if (n_cmfd_resets > 0) then + allocate(int_array(n_cmfd_resets)) + call get_node_array(doc, "tally_reset", int_array) + do i = 1, n_cmfd_resets + call cmfd_reset % add(int_array(i)) + end do + deallocate(int_array) end if - - ! Inactive batch flush window - if (check_for_node(doc, "inactive_flush")) & - call get_node_value(doc, "inactive_flush", cmfd_inact_flush(1)) - if (check_for_node(doc, "num_flushes")) & - call get_node_value(doc, "num_flushes", cmfd_inact_flush(2)) - - ! Last flush before active batches - if (check_for_node(doc, "active_flush")) & - call get_node_value(doc, "active_flush", cmfd_act_flush) ! Get display if (check_for_node(doc, "display")) & - call get_node_value(doc, "display", cmfd_display) + call get_node_value(doc, "display", cmfd_display) if (trim(cmfd_display) == 'dominance' .and. & - trim(cmfd_solver_type) /= 'power') then + trim(cmfd_solver_type) /= 'power') then message = 'Dominance Ratio only aviable with power iteration solver' call warning() cmfd_display = '' end if + ! Read in spectral radius estimate and tolerances + if (check_for_node(doc, "spectral")) & + call get_node_value(doc, "spectral", cmfd_spectral) + if (check_for_node(doc, "shift")) & + call get_node_value(doc, "shift", cmfd_shift) + if (check_for_node(doc, "ktol")) & + call get_node_value(doc, "ktol", cmfd_ktol) + if (check_for_node(doc, "stol")) & + call get_node_value(doc, "stol", cmfd_stol) + if (check_for_node(doc, "gauss_seidel_tolerance")) then + n_params = get_arraysize_double(doc, "gauss_seidel_tolerance") + if (n_params /= 2) then + message = 'Gauss Seidel tolerance is not 2 parameters & + &(absolute, relative).' + call fatal_error() + end if + call get_node_array(doc, "gauss_seidel_tolerance", gs_tol) + cmfd_atoli = gs_tol(1) + cmfd_rtoli = gs_tol(2) + end if + ! Create tally objects call create_cmfd_tally(doc) @@ -405,9 +439,9 @@ contains ! Set reset property if (check_for_node(doc, "reset")) then call get_node_value(doc, "reset", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - t % reset = .true. + t % reset = .true. end if ! Set up mesh filter diff --git a/src/cmfd_power_solver.F90 b/src/cmfd_power_solver.F90 index 65b60bc891..f80c122b71 100644 --- a/src/cmfd_power_solver.F90 +++ b/src/cmfd_power_solver.F90 @@ -242,7 +242,7 @@ contains subroutine convergence(iter) - use constants, only: ONE, TINY_BIT + use constants, only: ONE, ZERO use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV @@ -255,7 +255,7 @@ contains kerr = abs(k_o - k_n)/k_n ! Calculate max error in source - where (s_n % val > TINY_BIT) + where (s_n % val > ZERO) serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2 end where serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val)) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 new file mode 100644 index 0000000000..b3e8fb8a88 --- /dev/null +++ b/src/cmfd_solver.F90 @@ -0,0 +1,819 @@ +module cmfd_solver + +! This module contains routines to execute the power iteration solver + + use cmfd_loss_operator, only: init_loss_matrix, build_loss_matrix + use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix + use matrix_header, only: Matrix + use vector_header, only: Vector + + implicit none + private + public :: cmfd_solver_execute + + real(8) :: k_n ! new k-eigenvalue + real(8) :: k_o ! old k-eigenvalue + real(8) :: k_s ! shift of eigenvalue + real(8) :: k_ln ! new shifted eigenvalue + real(8) :: k_lo ! old shifted eigenvalue + real(8) :: norm_n ! current norm of source vector + real(8) :: norm_o ! old norm of source vector + real(8) :: kerr ! error in keff + real(8) :: serr ! error in source + real(8) :: ktol ! tolerance on keff + real(8) :: stol ! tolerance on source + logical :: adjoint_calc ! run an adjoint calculation + type(Matrix) :: loss ! cmfd loss matrix + type(Matrix) :: prod ! cmfd prod matrix + type(Vector) :: phi_n ! new flux vector + type(Vector) :: phi_o ! old flux vector + type(Vector) :: s_n ! new source vector + type(Vector) :: s_o ! old flux vector + type(Vector) :: serr_v ! error in source + + ! CMFD linear solver interface + procedure(linsolve), pointer :: cmfd_linsolver => null() + abstract interface + subroutine linsolve(A, b, x, tol, i) + import :: Matrix + import :: Vector + type(Matrix), intent(inout) :: A + type(Vector), intent(inout) :: b + type(Vector), intent(inout) :: x + real(8), intent(in) :: tol + integer, intent(out) :: i + end subroutine linsolve + end interface + +contains + +!=============================================================================== +! CMFD_SOLVER_EXECUTE sets up and runs power iteration solver for CMFD +!=============================================================================== + + subroutine cmfd_solver_execute(adjoint) + + use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve + + logical, optional, intent(in) :: adjoint ! adjoint calc + + logical :: physical_adjoint = .false. + + ! Check for adjoint execution + adjoint_calc = .false. + if (present(adjoint)) adjoint_calc = adjoint + + ! Check for physical adjoint + if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') & + physical_adjoint = .true. + + ! Start timer for build + call time_cmfdbuild % start() + + ! Initialize matrices and vectors + call init_data(physical_adjoint) + + ! Check for mathematical adjoint calculation + if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') & + call compute_adjoint() + + ! Stop timer for build + call time_cmfdbuild % stop() + + ! Begin power iteration + call time_cmfdsolve % start() + call execute_power_iter() + call time_cmfdsolve % stop() + + ! Extract results + call extract_results() + + ! Deallocate data + call finalize() + + end subroutine cmfd_solver_execute + +!=============================================================================== +! INIT_DATA allocates matrices and vectors for CMFD solution +!=============================================================================== + + subroutine init_data(adjoint) + + use constants, only: ONE, ZERO + use error, only: fatal_error + use global, only: cmfd, cmfd_shift, keff, cmfd_ktol, cmfd_stol, & + cmfd_write_matrices + + logical, intent(in) :: adjoint + + integer :: n ! problem size + real(8) :: guess ! initial guess + real(8) :: dw ! eigenvalue shift + + ! Set up matrices + call init_loss_matrix(loss) + call init_prod_matrix(prod) + + ! Get problem size + n = loss % n + + ! Set up flux vectors + call phi_n % create(n) + call phi_o % create(n) + + ! Set up source vectors + call s_n % create(n) + call s_o % create(n) + call serr_v % create(n) + + ! Set initial guess + guess = ONE + phi_n % val = guess + phi_o % val = guess + k_n = keff + k_o = k_n + dw = cmfd_shift + k_s = k_o + dw + k_ln = ONE/(ONE/k_n - ONE/k_s) + k_lo = k_ln + + ! Fill in loss matrix + call build_loss_matrix(loss, adjoint=adjoint) + + ! Fill in production matrix + call build_prod_matrix(prod, adjoint=adjoint) + + ! Finalize setup of CSR matrices + call loss % assemble() + call prod % assemble() + if (cmfd_write_matrices) then + call loss % write('loss.dat') + call prod % write('prod.dat') + end if + + ! Set norms to 0 + norm_n = ZERO + norm_o = ZERO + + ! Set up solver + select case(cmfd % indices(4)) + case(1) + cmfd_linsolver => cmfd_linsolver_1g + case(2) + cmfd_linsolver => cmfd_linsolver_2g + case default + cmfd_linsolver => cmfd_linsolver_ng + end select + + ! Set tolerances + ktol = cmfd_ktol + stol = cmfd_stol + + end subroutine init_data + +!=============================================================================== +! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem +!=============================================================================== + + subroutine compute_adjoint() + + use error, only: fatal_error +#ifdef PETSC + use global, only: cmfd_write_matrices +#else + use global, only: message +#endif + +#ifdef PETSC + ! Transpose matrices + call loss % transpose() + call prod % transpose() + + ! Write out matrix in binary file (debugging) + if (cmfd_write_matrices) then + call loss % write_petsc_binary('adj_lossmat.bin') + call prod % write_petsc_binary('adj_prodmat.bin') + end if +#else + message = 'Adjoint calculations only allowed with PETSc' + call fatal_error() +#endif + + end subroutine compute_adjoint + +!=============================================================================== +! EXECUTE_POWER_ITER is the main power iteration routine +! for the cmfd calculation +!=============================================================================== + + subroutine execute_power_iter() + + use constants, only: ONE + use error, only: fatal_error + use global, only: cmfd_atoli, cmfd_rtoli, message + + integer :: i ! iteration counter + integer :: innerits ! # of inner iterations + integer :: totalits ! total number of inners + logical :: iconv ! did the problem converged + real(8) :: atoli ! absolute minimum tolerance + real(8) :: rtoli ! relative tolerance based on source conv + real(8) :: toli ! the current tolerance of inners + + ! Reset convergence flag + iconv = .false. + + ! Set up tolerances + atoli = cmfd_atoli + rtoli = cmfd_rtoli + toli = rtoli*100._8 + + ! Perform shift + call wielandt_shift() + totalits = 0 + + ! Begin power iteration + do i = 1, 10000 + + ! Check if reached iteration 10000 + if (i == 10000) then + message = 'Reached maximum iterations in CMFD power iteration solver.' + call fatal_error() + end if + + ! Compute source vector + call prod % vector_multiply(phi_o, s_o) + + ! Normalize source vector + s_o % val = s_o % val / k_lo + + ! Compute new flux vector + call cmfd_linsolver(loss, s_o, phi_n, toli, innerits) + + ! Compute new source vector + call prod % vector_multiply(phi_n, s_n) + + ! Compute new shifted eigenvalue + k_ln = sum(s_n % val) / sum(s_o % val) + + ! Compute new eigenvalue + k_n = ONE/(ONE/k_ln + ONE/k_s) + + ! Renormalize the old source + s_o % val = s_o % val * k_lo + + ! Check convergence + call convergence(i, innerits, iconv) + totalits = totalits + innerits + + ! Break loop if converged + if (iconv) exit + + ! Record old values + phi_o % val = phi_n % val + k_o = k_n + k_lo = k_ln + norm_o = norm_n + + ! Get new tolerance for inners + toli = max(atoli, rtoli*serr) + + end do + + end subroutine execute_power_iter + +!=============================================================================== +! WIELANDT SHIFT +!=============================================================================== + + subroutine wielandt_shift() + + use constants, only: ONE + + integer :: irow ! row counter + integer :: icol ! col counter + integer :: jcol ! current col index in prod matrix + + ! perform subtraction + jcol = 1 + ROWS: do irow = 1, loss % n + COLS: do icol = loss % get_row(irow), loss % get_row(irow + 1) - 1 + if (loss % get_col(icol) == prod % get_col(jcol) .and. & + jcol < prod % get_row(irow + 1)) then + loss % val(icol) = loss % val(icol) - ONE/k_s*prod % val(jcol) + jcol = jcol + 1 + end if + end do COLS + end do ROWS + + end subroutine wielandt_shift + +!=============================================================================== +! CONVERGENCE checks the convergence of the CMFD problem +!=============================================================================== + + subroutine convergence(iter, innerits, iconv) + + use constants, only: ONE, ZERO + use global, only: cmfd_power_monitor, master + use, intrinsic :: ISO_FORTRAN_ENV + + integer, intent(in) :: iter ! outer iteration number + integer, intent(in) :: innerits ! inner iteration nubmer + logical, intent(out) :: iconv ! convergence logical + + ! Reset convergence flag + iconv = .false. + + ! Calculate error in keff + kerr = abs(k_o - k_n)/k_n + + ! Calculate max error in source + where (s_n % val > ZERO) + serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2 + end where + serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val)) + + ! Check for convergence + if(kerr < ktol .and. serr < stol) iconv = .true. + + ! Save the L2 norm of the source + norm_n = serr + + ! Print out to user + if (cmfd_power_monitor .and. master) then + write(OUTPUT_UNIT,FMT='(I0,":",T10,"k-eff: ",F0.8,T30,"k-error: ", & + &1PE12.5,T55, "src-error: ",1PE12.5,T80,I0)') iter, k_n, kerr, & + serr, innerits + end if + + end subroutine convergence + +!=============================================================================== +! CMFD_LINSOLVER_1g solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver_1g(A, b, x, tol, its) + + use constants, only: ONE, ZERO + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message + + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations + + integer :: g ! group index + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: igs ! Gauss-Seidel iteration counter + integer :: irb ! Red/Black iteration switch + integer :: irow ! row iteration + integer :: icol ! iteration counter over columns + integer :: didx ! index for diagonal component + logical :: found ! did we find col + real(8) :: tmp1 ! temporary sum g1 + real(8) :: x1 ! new g1 value of x + real(8) :: err ! error in convergence of solution + real(8) :: w ! overrelaxation parameter + type(Vector) :: tmpx ! temporary solution vector + + ! Set overrelaxation parameter + w = ONE + + ! Dimensions + ng = 1 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = A % n + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + + ! Copy over x vector + call tmpx % copy(x) + + ! Perform red/black gs iterations + REDBLACK: do irb = 0,1 + + ! Begin loop around matrix rows + ROWS: do irow = 1, n + + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + ! Filter out black cells (even) + if (mod(i+j+k,2) == irb) cycle + + ! Get the index of the diagonals for both rows + call A % search_indices(irow, irow, didx, found) + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + do icol = A % get_row(irow), didx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = didx + 1, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + + ! Solve for new x + x1 = (b % val(irow) - tmp1)/A % val(didx) + + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + + end do ROWS + + end do REDBLACK + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs + if (err < tol) exit + + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver_1g + +!=============================================================================== +! CMFD_LINSOLVER_2G solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver_2g(A, b, x, tol, its) + + use constants, only: ONE, ZERO + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message + + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations + + integer :: g ! group index + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: d1idx ! index of row "1" diagonal + integer :: d2idx ! index of row "2" diagonal + integer :: igs ! Gauss-Seidel iteration counter + integer :: irb ! Red/Black iteration switch + integer :: irow ! row iteration + integer :: icol ! iteration counter over columns + logical :: found ! did we find col + real(8) :: m11 ! block diagonal component 1,1 + real(8) :: m12 ! block diagonal component 1,2 + real(8) :: m21 ! block diagonal component 2,1 + real(8) :: m22 ! block diagonal component 2,2 + real(8) :: dm ! determinant of block diagonal + real(8) :: d11 ! inverse component 1,1 + real(8) :: d12 ! inverse component 1,2 + real(8) :: d21 ! inverse component 2,1 + real(8) :: d22 ! inverse component 2,2 + real(8) :: tmp1 ! temporary sum g1 + real(8) :: tmp2 ! temporary sum g2 + real(8) :: x1 ! new g1 value of x + real(8) :: x2 ! new g2 value of x + real(8) :: err ! error in convergence of solution + real(8) :: w ! overrelaxation parameter + type(Vector) :: tmpx ! temporary solution vector + + ! Set tolerance and overrelaxation parameter + w = ONE + + ! Dimensions + ng = 2 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = A % n + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + + ! Copy over x vector + call tmpx % copy(x) + + ! Perform red/black gs iterations + REDBLACK: do irb = 0,1 + + ! Begin loop around matrix rows + ROWS: do irow = 1, n, 2 + + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + ! Filter out black cells (even) + if (mod(i+j+k,2) == irb) cycle + + ! Get the index of the diagonals for both rows + call A % search_indices(irow, irow, d1idx, found) + call A % search_indices(irow + 1, irow + 1, d2idx, found) + + ! Get block diagonal + m11 = A % val(d1idx) ! group 1 diagonal + m12 = A % val(d1idx + 1) ! group 1 right of diagonal (sorted by col) + m21 = A % val(d2idx - 1) ! group 2 left of diagonal (sorted by col) + m22 = A % val(d2idx) ! group 2 diagonal + + ! Analytically invert the diagonal + dm = m11*m22 - m12*m21 + d11 = m22/dm + d12 = -m12/dm + d21 = -m21/dm + d22 = m11/dm + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + tmp2 = ZERO + do icol = A % get_row(irow), d1idx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = A % get_row(irow + 1), d2idx - 2 + tmp2 = tmp2 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = d1idx + 2, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = d2idx + 1, A % get_row(irow + 2) - 1 + tmp2 = tmp2 + A % val(icol)*x % val(A % get_col(icol)) + end do + + ! Adjust with RHS vector + tmp1 = b % val(irow) - tmp1 + tmp2 = b % val(irow + 1) - tmp2 + + ! Solve for new x + x1 = d11*tmp1 + d12*tmp2 + x2 = d21*tmp1 + d22*tmp2 + + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + x % val(irow + 1) = (ONE - w)*x % val(irow + 1) + w*x2 + + end do ROWS + + end do REDBLACK + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs + if (err < tol) exit + + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver_2g + +!=============================================================================== +! CMFD_LINSOLVER_ng solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver_ng(A, b, x, tol, its) + + use constants, only: ONE, ZERO + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message + + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations + + integer :: g ! group index + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: igs ! Gauss-Seidel iteration counter + integer :: irow ! row iteration + integer :: icol ! iteration counter over columns + integer :: didx ! index for diagonal component + logical :: found ! did we find col + real(8) :: tmp1 ! temporary sum g1 + real(8) :: x1 ! new g1 value of x + real(8) :: err ! error in convergence of solution + real(8) :: w ! overrelaxation parameter + type(Vector) :: tmpx ! temporary solution vector + + ! Set overrelaxation parameter + w = ONE + + ! Dimensions + ng = 1 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = A % n + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + + ! Copy over x vector + call tmpx % copy(x) + + ! Begin loop around matrix rows + ROWS: do irow = 1, n + + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + ! Get the index of the diagonals for both rows + call A % search_indices(irow, irow, didx, found) + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + do icol = A % get_row(irow), didx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = didx + 1, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + + ! Solve for new x + x1 = (b % val(irow) - tmp1)/A % val(didx) + + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + + end do ROWS + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs + + if (err < tol) exit + + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver_ng + +!=============================================================================== +! EXTRACT_RESULTS takes results and puts them in CMFD global data object +!=============================================================================== + + subroutine extract_results() + + use global, only: cmfd, cmfd_write_matrices, current_batch + + character(len=25) :: filename ! name of file to write data + integer :: n ! problem size + + ! Get problem size + n = loss % n + + ! Allocate in cmfd object if not already allocated + if (adjoint_calc) then + if (.not. allocated(cmfd%adj_phi)) allocate(cmfd%adj_phi(n)) + else + if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n)) + end if + + ! Save values + if (adjoint_calc) then + cmfd % adj_phi = phi_n % val + else + cmfd % phi = phi_n % val + end if + + ! Save eigenvalue + if(adjoint_calc) then + cmfd%adj_keff = k_n + else + cmfd%keff = k_n + end if + + ! Normalize phi to 1 + if (adjoint_calc) then + cmfd%adj_phi = cmfd%adj_phi/sqrt(sum(cmfd%adj_phi*cmfd%adj_phi)) + else + cmfd%phi = cmfd%phi/sqrt(sum(cmfd%phi*cmfd%phi)) + end if + + ! Save dominance ratio + cmfd % dom(current_batch) = norm_n/norm_o + + ! Write out results + if (cmfd_write_matrices) then + if (adjoint_calc) then + filename = 'adj_fluxvec.bin' + else + filename = 'fluxvec.bin' + end if +#ifdef PETSC + call phi_n % write_petsc_binary(filename) +#endif + end if + + end subroutine extract_results + +!=============================================================================== +! MATRIX_TO_INDICES converts a matrix index to spatial and group indicies +!=============================================================================== + + subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + use global, only: cmfd, cmfd_coremap + + integer, intent(out) :: i ! iteration counter for x + integer, intent(out) :: j ! iteration counter for y + integer, intent(out) :: k ! iteration counter for z + integer, intent(out) :: g ! iteration counter for groups + integer, intent(in) :: irow ! iteration counter over row (0 reference) + integer, intent(in) :: nx ! maximum number of x cells + integer, intent(in) :: ny ! maximum number of y cells + integer, intent(in) :: nz ! maximum number of z cells + integer, intent(in) :: ng ! maximum number of groups + + ! Check for core map + if (cmfd_coremap) then + + ! Get indices from indexmap + g = mod(irow-1, ng) + 1 + i = cmfd % indexmap((irow-1)/ng+1,1) + j = cmfd % indexmap((irow-1)/ng+1,2) + k = cmfd % indexmap((irow-1)/ng+1,3) + + else + + ! Compute indices + g = mod(irow-1, ng) + 1 + i = mod(irow-1, ng*nx)/ng + 1 + j = mod(irow-1, ng*nx*ny)/(ng*nx)+ 1 + k = mod(irow-1, ng*nx*ny*nz)/(ng*nx*ny) + 1 + + end if + + end subroutine matrix_to_indices + +!=============================================================================== +! FINALIZE frees all memory associated with power iteration +!=============================================================================== + + subroutine finalize() + + ! Destroy all objects + call loss % destroy() + call prod % destroy() + call phi_n % destroy() + call phi_o % destroy() + call s_n % destroy() + call s_o % destroy() + call serr_v % destroy + + end subroutine finalize + +end module cmfd_solver diff --git a/src/constants.F90 b/src/constants.F90 index 7d0657b695..ae9c77287a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -393,9 +393,6 @@ module constants ! constant to represent a zero flux "albedo" real(8), parameter :: ZERO_FLUX = 999.0_8 - ! constant to represent albedo rejection - real(8), parameter :: ALBEDO_REJECT = 999.0_8 - ! constant for writing out no residual real(8), parameter :: CMFD_NORES = 99999.0_8 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 0620cf11f8..8e4c19cccd 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -113,7 +113,7 @@ contains atom_density * micro_xs(i_nuclide) % elastic ! Add contributions to material macroscopic absorption cross section - material_xs % absorption = material_xs % absorption + & + material_xs % absorption = material_xs % absorption + & atom_density * micro_xs(i_nuclide) % absorption ! Add contributions to material macroscopic fission cross section @@ -123,7 +123,7 @@ contains ! Add contributions to material macroscopic nu-fission cross section material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission - + ! Add contributions to material macroscopic energy release from fission material_xs % kappa_fission = material_xs % kappa_fission + & atom_density * micro_xs(i_nuclide) % kappa_fission @@ -214,7 +214,7 @@ contains ! Calculate microscopic nuclide nu-fission cross section micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( & i_grid) + f * nuc % nu_fission(i_grid+1) - + ! Calculate microscopic nuclide kappa-fission cross section ! The ENDF standard (ENDF-102) states that MT 18 stores ! the fission energy as the Q_value (fission(1)) @@ -276,7 +276,7 @@ contains f = ZERO else i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) - f = (E - sab%inelastic_e_in(i_grid)) / & + f = (E - sab%inelastic_e_in(i_grid)) / & (sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid)) end if @@ -342,21 +342,22 @@ contains integer, intent(in) :: i_nuclide ! index into nuclides array real(8), intent(in) :: E ! energy - integer :: i_energy ! index for energy - integer :: i_low ! band index at lower bounding energy - integer :: i_up ! band index at upper bounding energy - real(8) :: f ! interpolation factor - real(8) :: r ! pseudo-random number - real(8) :: elastic ! elastic cross section - real(8) :: capture ! (n,gamma) cross section - real(8) :: fission ! fission cross section - real(8) :: inelastic ! inelastic cross section - logical :: same_nuc ! do we know the xs for this nuclide at this energy? + integer :: i ! loop index + integer :: i_energy ! index for energy + integer :: i_low ! band index at lower bounding energy + integer :: i_up ! band index at upper bounding energy + integer :: same_nuc_idx ! index of same nuclide + real(8) :: f ! interpolation factor + real(8) :: r ! pseudo-random number + real(8) :: elastic ! elastic cross section + real(8) :: capture ! (n,gamma) cross section + real(8) :: fission ! fission cross section + real(8) :: inelastic ! inelastic cross section + logical :: same_nuc ! do we know the xs for this nuclide at this energy? type(UrrData), pointer, save :: urr => null() type(Nuclide), pointer, save :: nuc => null() type(Reaction), pointer, save :: rxn => null() - type(ListElemInt), pointer :: nuc_list => null() -!$omp threadprivate(urr, nuc, rxn, nuc_list) +!$omp threadprivate(urr, nuc, rxn) micro_xs(i_nuclide) % use_ptable = .true. @@ -381,18 +382,16 @@ contains ! this energy but a different temperature, use the original random number to ! preserve correlation of temperature in probability tables same_nuc = .false. - nuc_list => nuc % nuc_list - do - if (E /= ZERO .and. E == micro_xs(nuc_list % data) % last_E) then + do i = 1, nuc % nuc_list % size() + if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % get_item(i)) % last_E) then same_nuc = .true. + same_nuc_idx = i exit end if - nuc_list => nuc_list % next - if (.not. associated(nuc_list % next)) exit end do if (same_nuc) then - r = micro_xs(nuc_list % data) % last_prn + r = micro_xs(nuc % nuc_list % get_item(same_nuc_idx)) % last_prn else r = prn() micro_xs(i_nuclide) % last_prn = r @@ -476,6 +475,11 @@ contains fission = fission * micro_xs(i_nuclide) % fission end if + ! Check for negative values + if (elastic < ZERO) elastic = ZERO + if (fission < ZERO) fission = ZERO + if (capture < ZERO) capture = ZERO + ! Set elastic, absorption, fission, and total cross sections. Note that the ! total cross section is calculated as sum of partials rather than using the ! table-provided value @@ -539,11 +543,11 @@ contains if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then i_grid = i_grid + 1 end if - + ! calculate interpolation factor f = (E - nuc % energy_0K(i_grid)) & & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid)) - + ! Calculate microscopic nuclide elastic cross section xs_out = (ONE - f) * nuc % elastic_0K(i_grid) & & + f * nuc % elastic_0K(i_grid + 1) diff --git a/src/global.F90 b/src/global.F90 index b136c2c166..b5f9c4a01e 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -299,6 +299,9 @@ module global ! Particle restart run logical :: particle_restart_run = .false. + ! Write out initial source + logical :: write_initial_source = .false. + ! ============================================================================ ! CMFD VARIABLES @@ -329,9 +332,6 @@ module global integer :: n_cmfd_meshes = 1 ! # of structured meshes integer :: n_cmfd_tallies = 3 ! # of user-defined tallies - ! Flag to hold cmfd weight adjustment factors - logical :: cmfd_hold_weights = .false. - ! Eigenvalue solver type character(len=10) :: cmfd_solver_type = 'power' @@ -344,11 +344,9 @@ module global ! Batch to begin cmfd integer :: cmfd_begin = 1 - ! When and how long to flush cmfd tallies during inactive batches - integer :: cmfd_inact_flush(2) = (/9999,1/) - - ! Batch to last flush before active batches - integer :: cmfd_act_flush = 0 + ! Tally reset list + integer :: n_cmfd_resets + type(SetInt) :: cmfd_reset ! Compute effective downscatter cross section logical :: cmfd_downscatter = .false. @@ -366,11 +364,18 @@ module global ! CMFD run logicals logical :: cmfd_on = .false. - logical :: cmfd_tally_on = .true. ! CMFD display info character(len=25) :: cmfd_display = 'balance' + ! Estimate of spectral radius of CMFD matrices and tolerances + real(8) :: cmfd_spectral = ZERO + real(8) :: cmfd_shift = 1.e6 + real(8) :: cmfd_ktol = 1.e-8_8 + real(8) :: cmfd_stol = 1.e-8_8 + real(8) :: cmfd_atoli = 1.e-10_8 + real(8) :: cmfd_rtoli = 1.e-5_8 + ! Information about state points to be written integer :: n_state_points = 0 type(SetInt) :: statepoint_batch diff --git a/src/input_xml.F90 b/src/input_xml.F90 index f852d2243d..a55fc363f4 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -11,7 +11,7 @@ module input_xml use output, only: write_message use plot_header use random_lcg, only: prn - use string, only: lower_case, to_str, str_to_int, str_to_real, & + use string, only: to_lower, to_str, str_to_int, str_to_real, & starts_with, ends_with use tally_header, only: TallyObject, TallyFilter use tally_initialize, only: add_tallies @@ -265,6 +265,14 @@ contains call fatal_error() end if + ! Check if we want to write out source + if (check_for_node(node_source, "write_initial")) then + call get_node_value(node_source, "write_initial", temp_str) + temp_str = to_lower(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & + write_initial_source = .true. + end if + ! Check for external source file if (check_for_node(node_source, "file")) then ! Copy path of source file @@ -290,8 +298,7 @@ contains type = '' if (check_for_node(node_dist, "type")) & call get_node_value(node_dist, "type", type) - call lower_case(type) - select case (trim(type)) + select case (to_lower(type)) case ('box') external_source % type_space = SRC_SPACE_BOX coeffs_reqd = 6 @@ -343,8 +350,7 @@ contains type = '' if (check_for_node(node_dist, "type")) & call get_node_value(node_dist, "type", type) - call lower_case(type) - select case (trim(type)) + select case (to_lower(type)) case ('isotropic') external_source % type_angle = SRC_ANGLE_ISOTROPIC coeffs_reqd = 0 @@ -395,8 +401,7 @@ contains type = '' if (check_for_node(node_dist, "type")) & call get_node_value(node_dist, "type", type) - call lower_case(type) - select case (trim(type)) + select case (to_lower(type)) case ('monoenergetic') external_source % type_energy = SRC_ENERGY_MONO coeffs_reqd = 1 @@ -446,7 +451,7 @@ contains ! Survival biasing if (check_for_node(doc, "survival_biasing")) then call get_node_value(doc, "survival_biasing", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & survival_biasing = .true. end if @@ -454,7 +459,7 @@ contains ! Probability tables if (check_for_node(doc, "ptables")) then call get_node_value(doc, "ptables", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') & urr_ptables_on = .false. end if @@ -691,19 +696,19 @@ contains ! Check if the user has specified to write binary source file if (check_for_node(node_sp, "separate")) then call get_node_value(node_sp, "separate", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') source_separate = .true. end if if (check_for_node(node_sp, "write")) then call get_node_value(node_sp, "write", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'false' .or. & trim(temp_str) == '0') source_write = .false. end if if (check_for_node(node_sp, "overwrite_latest")) then call get_node_value(node_sp, "overwrite_latest", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') then source_latest = .true. @@ -738,7 +743,7 @@ contains ! batch if (check_for_node(doc, "no_reduce")) then call get_node_value(doc, "no_reduce", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & reduce_tallies = .false. end if @@ -747,7 +752,7 @@ contains ! uncertainties rather than standard deviations if (check_for_node(doc, "confidence_intervals")) then call get_node_value(doc, "confidence_intervals", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') confidence_intervals = .true. end if @@ -761,7 +766,7 @@ contains ! Check for summary option if (check_for_node(node_output, "summary")) then call get_node_value(node_output, "summary", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') output_summary = .true. end if @@ -769,7 +774,7 @@ contains ! Check for cross sections option if (check_for_node(node_output, "cross_sections")) then call get_node_value(node_output, "cross_sections", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') output_xs = .true. end if @@ -777,7 +782,7 @@ contains ! Check for ASCII tallies output option if (check_for_node(node_output, "tallies")) then call get_node_value(node_output, "tallies", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'false' .or. & trim(temp_str) == '0') output_tallies = .false. end if @@ -786,15 +791,9 @@ contains ! Check for cmfd run if (check_for_node(doc, "run_cmfd")) then call get_node_value(doc, "run_cmfd", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then cmfd_run = .true. -#ifndef PETSC - if (master) then - message = 'CMFD is not available, compile OpenMC with PETSc' - call fatal_error() - end if -#endif end if end if @@ -802,7 +801,7 @@ contains if (check_for_node(doc, "resonance_scattering")) then call get_node_ptr(doc, "resonance_scattering", node_res_scat) call get_node_list(node_res_scat, "scatterer", node_scat_list) - + ! check that a nuclide is specified if (get_list_size(node_scat_list) >= 1) then treat_res_scat = .true. @@ -812,7 +811,7 @@ contains allocate(nuclides_0K(n_res_scatterers_total)) do i = 1, n_res_scatterers_total call get_list_item(node_scat_list, i, node_scatterer) - + ! check to make sure a nuclide is specified if (.not. check_for_node(node_scatterer, "nuclide")) then message = "No nuclide specified for scatterer " // trim(to_str(i)) & @@ -821,12 +820,12 @@ contains end if call get_node_value(node_scatterer, "nuclide", & nuclides_0K(i) % nuclide) - + if (check_for_node(node_scatterer, "method")) then call get_node_value(node_scatterer, "method", & nuclides_0K(i) % scheme) end if - + ! check to make sure xs name for which method is applied is given if (.not. check_for_node(node_scatterer, "xs_label")) then message = "Must specify the temperature dependent name of " // '' & @@ -835,7 +834,7 @@ contains end if call get_node_value(node_scatterer, "xs_label", & nuclides_0K(i) % name) - + ! check to make sure 0K xs name for which method is applied is given if (.not. check_for_node(node_scatterer, "xs_label_0K")) then message = "Must specify the 0K name of " // '' & @@ -844,7 +843,7 @@ contains end if call get_node_value(node_scatterer, "xs_label_0K", & nuclides_0K(i) % name_0K) - + if (check_for_node(node_scatterer, "E_min")) then call get_node_value(node_scatterer, "E_min", & nuclides_0K(i) % E_min) @@ -860,7 +859,7 @@ contains call get_node_value(node_scatterer, "E_max", & nuclides_0K(i) % E_max) end if - + ! check that E_max is not less than E_min if (nuclides_0K(i) % E_max < nuclides_0K(i) % E_min) then message = "Lower resonance scattering energy bound exceeds upper" @@ -868,8 +867,7 @@ contains end if nuclides_0K(i) % nuclide = trim(nuclides_0K(i) % nuclide) - nuclides_0K(i) % scheme = trim(nuclides_0K(i) % scheme) - call lower_case(nuclides_0K(i) % scheme) + nuclides_0K(i) % scheme = to_lower(trim(nuclides_0K(i) % scheme)) nuclides_0K(i) % name = trim(nuclides_0K(i) % name) nuclides_0K(i) % name_0K = trim(nuclides_0K(i) % name_0K) end do @@ -883,8 +881,7 @@ contains ! Natural element expansion option if (check_for_node(doc, "natural_elements")) then call get_node_value(doc, "natural_elements", temp_str) - call lower_case(temp_str) - select case (temp_str) + select case (to_lower(temp_str)) case ('endf/b-vii.0') default_expand = ENDF_BVII0 case ('endf/b-vii.1') @@ -1017,8 +1014,7 @@ contains word = '' if (check_for_node(node_cell, "material")) & call get_node_value(node_cell, "material", word) - call lower_case(word) - select case(word) + select case(to_lower(word)) case ('void') c % material = MATERIAL_VOID @@ -1187,8 +1183,7 @@ contains word = '' if (check_for_node(node_surf, "type")) & call get_node_value(node_surf, "type", word) - call lower_case(word) - select case(trim(word)) + select case(to_lower(word)) case ('x-plane') s % type = SURF_PX coeffs_reqd = 1 @@ -1249,8 +1244,7 @@ contains word = '' if (check_for_node(node_surf, "boundary")) & call get_node_value(node_surf, "boundary", word) - call lower_case(word) - select case (trim(word)) + select case (to_lower(word)) case ('transmission', 'transmit', '') s % bc = BC_TRANSMIT case ('vacuum') @@ -1312,8 +1306,7 @@ contains word = '' if (check_for_node(node_lat, "type")) & call get_node_value(node_lat, "type", word) - call lower_case(word) - select case (trim(word)) + select case (to_lower(word)) case ('rect', 'rectangle', 'rectangular') lat % type = LATTICE_RECT case ('hex', 'hexagon', 'hexagonal') @@ -1544,8 +1537,7 @@ contains end if ! Adjust material density based on specified units - call lower_case(units) - select case(trim(units)) + select case(to_lower(units)) case ('g/cc', 'g/cm3') mat % density = -val case ('kg/m3') @@ -1700,7 +1692,7 @@ contains ALL_NUCLIDES: do j = 1, mat % n_nuclides ! Check that this nuclide is listed in the cross_sections.xml file name = trim(list_names % get_item(j)) - if (.not. xs_listing_dict % has_key(name)) then + if (.not. xs_listing_dict % has_key(to_lower(name))) then message = "Could not find nuclide " // trim(name) // & " in cross_sections.xml file!" call fatal_error() @@ -1715,20 +1707,20 @@ contains end if ! Find xs_listing and set the name/alias according to the listing - index_list = xs_listing_dict % get_key(name) + index_list = xs_listing_dict % get_key(to_lower(name)) name = xs_listings(index_list) % name alias = xs_listings(index_list) % alias ! If this nuclide hasn't been encountered yet, we need to add its name ! and alias to the nuclide_dict - if (.not. nuclide_dict % has_key(name)) then + if (.not. nuclide_dict % has_key(to_lower(name))) then index_nuclide = index_nuclide + 1 mat % nuclide(j) = index_nuclide - call nuclide_dict % add_key(name, index_nuclide) - call nuclide_dict % add_key(alias, index_nuclide) + call nuclide_dict % add_key(to_lower(name), index_nuclide) + call nuclide_dict % add_key(to_lower(alias), index_nuclide) else - mat % nuclide(j) = nuclide_dict % get_key(name) + mat % nuclide(j) = nuclide_dict % get_key(to_lower(name)) end if ! Copy name and atom/weight percent @@ -1787,7 +1779,7 @@ contains mat % sab_names(j) = name ! Check that this nuclide is listed in the cross_sections.xml file - if (.not. xs_listing_dict % has_key(name)) then + if (.not. xs_listing_dict % has_key(to_lower(name))) then message = "Could not find S(a,b) table " // trim(name) // & " in cross_sections.xml file!" call fatal_error() @@ -1795,17 +1787,17 @@ contains ! Find index in xs_listing and set the name and alias according to the ! listing - index_list = xs_listing_dict % get_key(name) + index_list = xs_listing_dict % get_key(to_lower(name)) name = xs_listings(index_list) % name ! If this S(a,b) table hasn't been encountered yet, we need to add its ! name and alias to the sab_dict - if (.not. sab_dict % has_key(name)) then + if (.not. sab_dict % has_key(to_lower(name))) then index_sab = index_sab + 1 mat % i_sab_tables(j) = index_sab - call sab_dict % add_key(name, index_sab) + call sab_dict % add_key(to_lower(name), index_sab) else - mat % i_sab_tables(j) = sab_dict % get_key(name) + mat % i_sab_tables(j) = sab_dict % get_key(to_lower(name)) end if end do end if @@ -1916,7 +1908,7 @@ contains ! Check for setting if (check_for_node(doc, "assume_separate")) then call get_node_value(doc, "assume_separate", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & assume_separate = .true. end if @@ -1949,8 +1941,7 @@ contains temp_str = '' if (check_for_node(node_mesh, "type")) & call get_node_value(node_mesh, "type", temp_str) - call lower_case(temp_str) - select case (trim(temp_str)) + select case (to_lower(temp_str)) case ('rect', 'rectangle', 'rectangular') m % type = LATTICE_RECT case ('hex', 'hexagon', 'hexagonal') @@ -2133,7 +2124,7 @@ contains temp_str = '' if (check_for_node(node_filt, "type")) & call get_node_value(node_filt, "type", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) ! Determine number of bins if (check_for_node(node_filt, "bins")) then @@ -2361,14 +2352,14 @@ contains end if ! Check to make sure nuclide specified is in problem - if (.not. nuclide_dict % has_key(word)) then + if (.not. nuclide_dict % has_key(to_lower(word))) then message = "The nuclide " // trim(word) // " from tally " // & trim(to_str(t % id)) // " is not present in any material." call fatal_error() end if ! Set bin to index in nuclides array - t % nuclide_bins(j) = nuclide_dict % get_key(word) + t % nuclide_bins(j) = nuclide_dict % get_key(to_lower(word)) end do ! Set number of nuclide bins @@ -2399,7 +2390,7 @@ contains ! (i.e., scatter-p#, flux-y#) n_new = 0 do j = 1, n_words - call lower_case(sarray(j)) + sarray(j) = to_lower(sarray(j)) ! Find if scores(j) is of the form 'moment-p' or 'moment-y' present in ! MOMENT_STRS(:) ! If so, check the order, store if OK, then reset the number to 'n' @@ -2665,6 +2656,12 @@ contains ! Get index of mesh filter k = t % find_filter(FILTER_MESH) + ! Check to make sure mesh filter was specified + if (k == 0) then + message = "Cannot tally surface current without a mesh filter." + call fatal_error() + end if + ! Get pointer to mesh i_mesh = t % filters(k) % int_bins(1) m => meshes(i_mesh) @@ -2846,7 +2843,7 @@ contains temp_str = 'slice' if (check_for_node(node_plot, "type")) & call get_node_value(node_plot, "type", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("slice") pl % type = PLOT_TYPE_SLICE @@ -2913,7 +2910,7 @@ contains temp_str = 'xy' if (check_for_node(node_plot, "basis")) & call get_node_value(node_plot, "basis", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("xy") pl % basis = PLOT_BASIS_XY @@ -2960,7 +2957,7 @@ contains temp_str = "cell" if (check_for_node(node_plot, "color")) & call get_node_value(node_plot, "color", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("cell") @@ -3415,9 +3412,9 @@ contains end if ! create dictionary entry for both name and alias - call xs_listing_dict % add_key(listing % name, i) + call xs_listing_dict % add_key(to_lower(listing % name), i) if (check_for_node(node_ace, "alias")) then - call xs_listing_dict % add_key(listing % alias, i) + call xs_listing_dict % add_key(to_lower(listing % alias), i) end if end do @@ -3455,9 +3452,8 @@ contains character(2) :: element_name element_name = name(1:2) - call lower_case(element_name) - select case (element_name) + select case (to_lower(element_name)) case ('h') call list_names % append('1001.' // xs) call list_density % append(density * 0.999885_8) diff --git a/src/matrix_header.F90 b/src/matrix_header.F90 index 24307b10bb..3e3d4ad0db 100644 --- a/src/matrix_header.F90 +++ b/src/matrix_header.F90 @@ -20,19 +20,22 @@ module matrix_header # endif logical :: petsc_active contains - procedure :: create => matrix_create - procedure :: destroy => matrix_destroy - procedure :: add_value => matrix_add_value - procedure :: new_row => matrix_new_row - procedure :: assemble => matrix_assemble - procedure :: get_row => matrix_get_row - procedure :: get_col => matrix_get_col + procedure :: create => matrix_create + procedure :: destroy => matrix_destroy + procedure :: add_value => matrix_add_value + procedure :: new_row => matrix_new_row + procedure :: assemble => matrix_assemble + procedure :: get_row => matrix_get_row + procedure :: get_col => matrix_get_col procedure :: vector_multiply => matrix_vector_multiply -#ifdef PETSC - procedure :: transpose => matrix_transpose + procedure :: search_indices => matrix_search_indices + procedure :: write => matrix_write + procedure :: copy => matrix_copy +# ifdef PETSC procedure :: setup_petsc => matrix_setup_petsc procedure :: write_petsc_binary => matrix_write_petsc_binary -#endif + procedure :: transpose => matrix_transpose +# endif end type matrix #ifdef PETSC @@ -359,4 +362,87 @@ contains end subroutine matrix_vector_multiply +!=============================================================================== +! MATRIX_SEARCH_INDICES searches for an index in column corresponding to a row +!=============================================================================== + + subroutine matrix_search_indices(self, row, col, idx, found) + + class(Matrix), intent(inout) :: self + integer, intent(in) :: row + integer, intent(in) :: col + integer, intent(out) :: idx + logical, intent(out) :: found + + integer :: j + + found = .false. + + COLS: do j = self % get_row(row), self % get_row(row + 1) - 1 + + if (self % get_col(j) == col) then + idx = j + found = .true. + exit + end if + + end do COLS + + end subroutine matrix_search_indices + +!=============================================================================== +! MATRIX_WRITE writes a matrix to file +!=============================================================================== + + subroutine matrix_write(self, filename) + + character(*), intent(in) :: filename + class(Matrix), intent(inout) :: self + + integer :: unit_ + integer :: i + integer :: j + + open(newunit=unit_, file=filename) + + do i = 1, self % n + do j = self % get_row(i), self % get_row(i + 1) - 1 + write(unit_,*) i, self % get_col(j), self % val(j) + end do + end do + + close(unit_) + + end subroutine matrix_write + +!=============================================================================== +! MATRIX_COPY copies a matrix +!=============================================================================== + + subroutine matrix_copy(self, mattocopy) + + class(Matrix), intent(inout) :: self + type(Matrix), intent(in) :: mattocopy + + ! Set n and nnz + self % n_count = mattocopy % n_count + self % nz_count = mattocopy % nz_count + self % n = mattocopy % n + self % nnz = mattocopy % nnz + + ! Allocate vectors + if (.not.allocated(self % row)) allocate(self % row(self % n + 1)) + if (.not.allocated(self % col)) allocate(self % col(self % nnz)) + if (.not.allocated(self % val)) allocate(self % val(self % nnz)) + + ! Set PETSc active to false + self % petsc_active = .false. + + ! Copy over data + self % row = mattocopy % row + self % col = mattocopy % col + self % val = mattocopy % val + + end subroutine matrix_copy + end module matrix_header diff --git a/src/output.F90 b/src/output.F90 index 834c28e188..1493b82943 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -13,7 +13,7 @@ module output use mesh, only: mesh_indices_to_bin, bin_to_mesh_indices use particle_header, only: LocalCoord, Particle use plot_header - use string, only: upper_case, to_str + use string, only: to_upper, to_str use tally_header, only: TallyObject implicit none @@ -130,8 +130,7 @@ contains if (mod(len_trim(msg),2) == 0) m = m + 1 ! convert line to upper case - line = msg - call upper_case(line) + line = to_upper(msg) ! print header based on level select case (header_level) diff --git a/src/relaxng/cmfd.rnc b/src/relaxng/cmfd.rnc index cbcb7709ae..bbb7fc9e6b 100644 --- a/src/relaxng/cmfd.rnc +++ b/src/relaxng/cmfd.rnc @@ -22,15 +22,9 @@ element cmfd { element feedback { xsd:boolean }? & - element n_cmfd_procs { xsd:int }? & - - element reset { xsd:boolean }? & - - element balance { xsd:boolean }? & - element downscatter { xsd:boolean }? & - element run_2grp { xsd:boolean }? & + element dhat_reset { xsd:boolean }? & element solver { xsd:string }? & @@ -40,8 +34,6 @@ element cmfd { element power_monitor { xsd:boolean }? & - element write_balance { xsd:boolean }? & - element write_matrices { xsd:boolean }? & element run_adjoint { xsd:boolean }? & @@ -50,9 +42,18 @@ element cmfd { element begin { xsd:int }? & - element inactive { xsd:boolean }? & + element tally_reset { list { xsd:int+ } }? & - element active_flush { xsd:int }? & + element display { xsd:string }? & + + element spectral { xsd:double }? & + + element shift { xsd:double }? & + + element ktol { xsd: double }? & + + element stol { xsd: double }? & + + element gauss_seidel_tolerance { list { xsd:double+ } }? - element keff_tol { xsd:double }? } diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index dd4b91dc7d..707ecbbc80 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -85,7 +85,8 @@ element settings { attribute interplation { xsd:string { maxLength = "10" } })? & (element parameters { list { xsd:double+ } } | attribute parameters { list { xsd:double+ } })? - }? + }? & + (element write_initial { xsd:boolean } | attribute write_initial { xsd:boolean })? }? & element state_point { diff --git a/src/source.F90 b/src/source.F90 index 5f928436b9..9575020e6d 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -27,6 +27,7 @@ contains subroutine initialize_source() + character(MAX_FILE_LEN) :: filename integer(8) :: i ! loop index over bank sites integer(8) :: id ! particle id integer(4) :: itmp ! temporary integer @@ -76,6 +77,20 @@ contains end do end if + ! Write out initial source + if (write_initial_source) then + message = 'Writing out initial source guess...' + call write_message(1) +#ifdef HDF5 + filename = trim(path_output) // 'initial_source.h5' +#else + filename = trim(path_output) // 'initial_source.binary' +#endif + call sp % file_create(filename, serial = .false.) + call sp % write_source_bank() + call sp % file_close() + end if + end subroutine initialize_source !=============================================================================== diff --git a/src/string.F90 b/src/string.F90 index f78b94d07f..d4f8873a0a 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -152,37 +152,47 @@ contains ! LOWER_CASE converts a string to all lower case characters !=============================================================================== - elemental subroutine lower_case(word) + elemental function to_lower(word) result(word_lower) - character(*), intent(inout) :: word + character(*), intent(in) :: word + character(len=len(word)) :: word_lower integer :: i integer :: ic do i = 1, len(word) ic = ichar(word(i:i)) - if (ic >= 65 .and. ic <= 90) word(i:i) = char(ic+32) + if (ic >= 65 .and. ic <= 90) then + word_lower(i:i) = char(ic+32) + else + word_lower(i:i) = word(i:i) + end if end do - end subroutine lower_case + end function to_lower !=============================================================================== ! UPPER_CASE converts a string to all upper case characters !=============================================================================== - elemental subroutine upper_case(word) + elemental function to_upper(word) result(word_upper) - character(*), intent(inout) :: word + character(*), intent(in) :: word + character(len=len(word)) :: word_upper integer :: i integer :: ic do i = 1, len(word) ic = ichar(word(i:i)) - if (ic >= 97 .and. ic <= 122) word(i:i) = char(ic-32) + if (ic >= 97 .and. ic <= 122) then + word_upper(i:i) = char(ic-32) + else + word_upper(i:i) = word(i:i) + end if end do - end subroutine upper_case + end function to_upper !=============================================================================== ! ZERO_PADDED returns a string of the input integer padded with zeros to the @@ -317,7 +327,7 @@ end function zero_padded ! the loop automatically exits when n_digits = 10. n_digits = n_digits + 1 end do - + end function count_digits !=============================================================================== @@ -349,21 +359,21 @@ end function zero_padded end function int8_to_str !=============================================================================== -! STR_TO_INT converts a string to an integer. +! STR_TO_INT converts a string to an integer. !=============================================================================== function str_to_int(str) result(num) character(*), intent(in) :: str integer(8) :: num - + character(5) :: fmt integer :: w integer :: ioError ! Determine width of string w = len_trim(str) - + ! Create format specifier for reading string write(UNIT=fmt, FMT='("(I",I2,")")') w @@ -404,7 +414,7 @@ end function zero_padded integer :: decimal ! number of places after decimal integer :: width ! total field width - real(8) :: num2 ! absolute value of number + real(8) :: num2 ! absolute value of number character(9) :: fmt ! format specifier for writing number ! set default field width diff --git a/src/vector_header.F90 b/src/vector_header.F90 index e050981c76..05cc7e3cde 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -21,10 +21,11 @@ module vector_header procedure :: create => vector_create procedure :: destroy => vector_destroy procedure :: add_value => vector_add_value -#ifdef PETSC + procedure :: copy => vector_copy +# ifdef PETSC procedure :: setup_petsc => vector_setup_petsc procedure :: write_petsc_binary => vector_write_petsc_binary -#endif +# endif end type Vector #ifdef PETSC @@ -127,4 +128,28 @@ contains end subroutine vector_write_petsc_binary #endif +!=============================================================================== +! VECTOR_COPY allocates a separate vector and copies +!=============================================================================== + + subroutine vector_copy(self, vectocopy) + + class(Vector), target, intent(inout) :: self + type(Vector), intent(in) :: vectocopy + + ! Preallocate vector + if (.not.allocated(self % data)) allocate(self % data(vectocopy % n)) + self % val => self % data(1:vectocopy % n) + + ! Set n + self % n = vectocopy % n + + ! Copy values + self % val = vectocopy % val + + ! Petsc is default not active + self % petsc_active = .false. + + end subroutine vector_copy + end module vector_header diff --git a/tests/test_cmfd_feed/cmfd.xml b/tests/test_cmfd_feed/cmfd.xml index 51def74088..228b260a46 100644 --- a/tests/test_cmfd_feed/cmfd.xml +++ b/tests/test_cmfd_feed/cmfd.xml @@ -12,5 +12,5 @@ dominance power true - + 1.e-15 1.e-20 diff --git a/tests/test_cmfd_nofeed/cmfd.xml b/tests/test_cmfd_nofeed/cmfd.xml index c36a396043..eb6c2c721a 100644 --- a/tests/test_cmfd_nofeed/cmfd.xml +++ b/tests/test_cmfd_nofeed/cmfd.xml @@ -12,5 +12,6 @@ dominance power false + 1.e-15 1.e-20