crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'master' into KappFiss

Browse source 
Conflicts:
	src/constants.F90
	src/tally.F90
	src/utils/statepoint.py
This commit is contained in:
Paul Romano 2013-01-17 12:39:09 -05:00
commit e0978a80ea
9 changed files with 503 additions and 102 deletions
Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -13,6 +13,9 @@ src/openmc
# emacs backups
*~
# OpenMC statepoints
*.binary
# Documentation builds
docs/build

120
docs/source/usersguide/install.rst
View file

@ -126,11 +126,33 @@ This will build an executable named ``openmc``.
Compiling on Windows
--------------------
To compile OpenMC on a Windows operating system, you will need to first install
Cygwin_, a Linux-like environment for Windows. When configuring Cygwin, make
sure you install both the gfortran compiler as well as git. Once you have
obtained the source code, run the following commands from within the source code
root directory:
Using Cygwin
------------
One option for compiling OpenMC on a Windows operating system is to use Cygwin_,
a Linux-like environment for Windows. You will need to first `install
Cygwin`_. When you are asked to select packages, make sure the following are
selected:
* Devel: gcc4-core
* Devel: gcc4-fortran
* Devel: make
If you plan on obtaining the source code directly using git, select the
following packages:
* Devel: git
* Devel: git-completion (Optional)
* Devel: gitk (Optional)
In order to use the Python scripts provided with OpenMC, you will also need to
install Python. This can be done within Cygwin or directly in Windows. To
install within Cygwin, select the following packages:
* Python: python (Version > 2.7 recommended)
Once you have obtained the source code, run the following commands from within
the source code root directory:
.. code-block:: sh
@ -139,29 +161,93 @@ root directory:
This will build an executable named ``openmc``.
.. _Cygwin: http://www.cygwin.com/
.. _Cygwin: http://cygwin.com/
.. _install Cygwin: http://cygwin.com/setup.exe
Using MinGW
-----------
An alternate option for installing OpenMC on Windows is using MinGW_, which
stands for Minimalist GNU for Windows. An executable for installing the MinGW
distribution is available on SourceForge_. When installing MinGW, make sure the
following components are selected:
* MinGW Compiler Suite: Fortran Compiler
* MSYS Basic System
Once MinGW is installed, copy the OpenMC source distribution to your MinGW home
directory (usually C:\\MinGW\\msys\\1.0\\home\\YourUsername). Once you have
the source code in place, run the following commands from within the MinGW shell
in the root directory of the OpenMC distribution:
.. code-block:: sh
cd src
make
This will build an executable named ``openmc``.
.. _MinGW: http://www.mingw.org
.. _SourceForge: http://sourceforge.net/projects/mingw
---------------------------
Cross-Section Configuration
Cross Section Configuration
---------------------------
In order to run a simulation with OpenMC, you will need cross-section data for
each nuclide in your problem. Since OpenMC uses ACE format cross-sections, you
can use nuclear data distributed with MCNP_ or Serpent_.
In order to run a simulation with OpenMC, you will need cross section data for
each nuclide in your problem. Since OpenMC uses ACE format cross sections, you
can use nuclear data that was processed with NJOY, such as that distributed with
MCNP_ or Serpent_.
Using JEFF Cross Sections from OECD/NEA
---------------------------------------
The NEA_ provides processed ACE data from the JEFF_ nuclear library upon
request. A DVD of the data can be requested here_. To use this data with OpenMC,
the following steps must be taken:
1. Copy and unzip the data on the DVD to a directory on your computer.
2. In the root directory, a file named ``xsdir``, or some variant thereof,
should be present. This file contains a listing of all the cross sections and
is used by MCNP. This file should be converted to a ``cross_sections.xml``
file for use with OpenMC. A Python script is provided in the OpenMC
distribution for this purpose:
.. code-block:: sh
openmc/src/utils/convert_xsdir.py xsdir31 cross_sections.xml
3. In the converted ``cross_sections.xml`` file, change the contents of the
<directory> element to the absolute path of the directory containing the
actual ACE files.
4. Additionally, you may need to change any occurrences of upper-case "ACE"
within the ``cross_sections.xml`` file to lower-case.
5. Either set the :ref:`cross_sections` in a settings.xml file or the
:envvar:`CROSS_SECTIONS` environment variable to the absolute path of the
``cross_sections.xml`` file.
Using Cross Sections from MCNP
------------------------------
To use cross sections distributed with MCNP, change the <directory> element in
the ``cross_sections.xml`` file in the root directory of the OpenMC distribution
to the location of the MCNP cross-sections. Then, either set the
to the location of the MCNP cross sections. Then, either set the
:ref:`cross_sections` in a settings.xml file or the :envvar:`CROSS_SECTIONS`
environment variable to the absolute path of the ``cross_sections.xml`` file.
Similarly, to use cross-sections distributed with Serpent, change the
<directory> element in the ``cross_sections_serpent.xml`` file in the root
directory of the OpenMC distribution to the location of the Serpent
cross-sections. Then, either set the :ref:`cross_sections` in a settings.xml
file or the :envvar:`CROSS_SECTIONS` environment variable to the absolute path
of the ``cross_sections_serpent.xml`` file.
Using Cross Sections from Serpent
---------------------------------
To use cross sections distributed with Serpent, change the <directory> element
in the ``cross_sections_serpent.xml`` file in the root directory of the OpenMC
distribution to the location of the Serpent cross sections. Then, either set the
:ref:`cross_sections` in a settings.xml file or the :envvar:`CROSS_SECTIONS`
environment variable to the absolute path of the ``cross_sections_serpent.xml``
file.
.. _NEA: http://www.oecd-nea.org
.. _JEFF: http://www.oecd-nea.org/dbdata/jeff/
.. _here: http://www.oecd-nea.org/dbdata/pubs/jeff312-cd.html
.. _MCNP: http://mcnp.lanl.gov
.. _Serpent: http://montecarlo.vtt.fi

2
src/DEPENDENCIES
View file

@ -336,6 +336,7 @@ string.o: constants.o
string.o: error.o
string.o: global.o
tally.o: ace_header.o
tally.o: constants.o
tally.o: datatypes_header.o
tally.o: error.o
@ -344,6 +345,7 @@ tally.o: math.o
tally.o: mesh.o
tally.o: mesh_header.o
tally.o: output.o
tally.o: particle_header.o
tally.o: search.o
tally.o: string.o
tally.o: tally_header.o

1
src/OBJECTS
View file

@ -13,6 +13,7 @@ cmfd_output.o \
cmfd_power_solver.o \
cmfd_prod_operator.o \
cmfd_snes_solver.o \
constants.o \
cross_section.o \
datatypes.o \
datatypes_header.o \

17
src/constants.F90
View file

@ -252,7 +252,7 @@ module constants
EVENT_FISSION = 3
! Tally score type
integer, parameter :: N_SCORE_TYPES = 18
integer, parameter :: N_SCORE_TYPES = 15
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
@ -263,15 +263,12 @@ module constants
SCORE_TRANSPORT = -7, & ! transport reaction rate
SCORE_DIFFUSION = -8, & ! diffusion coefficient
SCORE_N_1N = -9, & ! (n,1n) rate
SCORE_N_2N = -10, & ! (n,2n) rate
SCORE_N_3N = -11, & ! (n,3n) rate
SCORE_N_4N = -12, & ! (n,4n) rate
SCORE_ABSORPTION = -13, & ! absorption rate
SCORE_FISSION = -14, & ! fission rate
SCORE_NU_FISSION = -15, & ! neutron production rate
SCORE_K_FISSION = -16, & ! fission energy production rate
SCORE_CURRENT = -17, & ! partial current
SCORE_EVENTS = -18 ! number of events
SCORE_ABSORPTION = -10, & ! absorption rate
SCORE_FISSION = -11, & ! fission rate
SCORE_NU_FISSION = -12, & ! neutron production rate
SCORE_K_FISSION = -13, & ! fission energy production rate
SCORE_CURRENT = -14, & ! partial current
SCORE_EVENTS = -15 ! number of events
! Maximum scattering order supported
integer, parameter :: SCATT_ORDER_MAX = 10

36
src/input_xml.F90
View file

@ -1221,6 +1221,7 @@ contains
integer :: n_scores ! number of tot scores after adjusting for Pn tally
integer :: n_order ! Scattering order requested
integer :: n_order_pos ! Position of Scattering order in score name string
integer :: MT ! user-specified MT for score
logical :: file_exists ! does tallies.xml file exist?
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
@ -1823,20 +1824,14 @@ contains
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('n2n')
t % score_bins(j) = SCORE_N_2N
t % score_bins(j) = N_2N
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('n3n')
t % score_bins(j) = SCORE_N_3N
t % score_bins(j) = N_3N
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('n4n')
t % score_bins(j) = SCORE_N_4N
t % score_bins(j) = N_4N
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('absorption')
t % score_bins(j) = SCORE_ABSORPTION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
@ -1914,9 +1909,26 @@ contains
t % score_bins(j) = SCORE_EVENTS
case default
message = "Unknown scoring function: " // &
trim(tally_(i) % scores(j))
call fatal_error()
! Assume that user has specified an MT number
MT = str_to_int(score_name)
if (MT /= ERROR_INT) then
! Specified score was an integer
if (MT > 1) then
t % score_bins(j) = MT
else
message = "Invalid MT on <scores>: " // &
trim(tally_(i) % scores(j))
call fatal_error()
end if
else
! Specified score was not an integer
message = "Unknown scoring function: " // &
trim(tally_(i) % scores(j))
call fatal_error()
end if
end select
end do
t % n_score_bins = n_scores

19
src/output.F90
View file

@ -809,12 +809,6 @@ contains
string = trim(string) // ' diffusion'
case (SCORE_N_1N)
string = trim(string) // ' n1n'
case (SCORE_N_2N)
string = trim(string) // ' n2n'
case (SCORE_N_3N)
string = trim(string) // ' n3n'
case (SCORE_N_4N)
string = trim(string) // ' n4n'
case (SCORE_ABSORPTION)
string = trim(string) // ' absorption'
case (SCORE_FISSION)
@ -825,6 +819,8 @@ contains
string = trim(string) // ' k-fission'
case (SCORE_CURRENT)
string = trim(string) // ' current'
case default
string = trim(string) // ' ' // reaction_name(t % score_bins(j))
end select
end do
write(unit_,*) ' Scores:' // trim(string)
@ -1457,13 +1453,10 @@ contains
score_names(abs(SCORE_TRANSPORT)) = "Transport Rate"
score_names(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient"
score_names(abs(SCORE_N_1N)) = "(n,1n) Rate"
score_names(abs(SCORE_N_2N)) = "(n,2n) Rate"
score_names(abs(SCORE_N_3N)) = "(n,3n) Rate"
score_names(abs(SCORE_N_4N)) = "(n,4n) Rate"
score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate"
score_names(abs(SCORE_FISSION)) = "Fission Rate"
score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate"
score_names(abs(SCORE_K_FISSION)) = "Kappa-Fission Rate"
score_names(abs(SCORE_K_FISSION)) = "Kappa-Fission Rate"
score_names(abs(SCORE_EVENTS)) = "Events"
! Create filename for tally output
@ -1626,7 +1619,11 @@ contains
end do
k = k + n_order - 1
else
score_name = score_names(abs(t % score_bins(k)))
if (t % score_bins(k) > 0) then
score_name = reaction_name(t % score_bins(k))
else
score_name = score_names(abs(t % score_bins(k)))
end if
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &

316
src/tally.F90
View file

@ -1,17 +1,20 @@
module tally
use ace_header, only: Reaction
use constants
use datatypes_header, only: ListInt
use error, only: fatal_error
use error, only: fatal_error
use global
use math, only: t_percentile, calc_pn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, &
mesh_indices_to_bin, mesh_intersects
use mesh_header, only: StructuredMesh
use output, only: header
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyResult, TallyMapItem, TallyMapElement
use math, only: t_percentile, calc_pn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, &
get_mesh_indices, mesh_indices_to_bin, &
mesh_intersects
use mesh_header, only: StructuredMesh
use output, only: header
use particle_header, only: LocalCoord
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyResult, TallyMapItem, TallyMapElement
#ifdef MPI
use mpi
@ -245,21 +248,6 @@ contains
! All events that reach this point are (n,1n) reactions
score = last_wgt
case (SCORE_N_2N)
! Skip any event that is not (n,2n)
if (p % event_MT /= N_2N) cycle SCORE_LOOP
score = last_wgt
case (SCORE_N_3N)
! Skip any event that is not (n,3n)
if (p % event_MT /= N_3N) cycle SCORE_LOOP
score = last_wgt
case (SCORE_N_4N)
! Skip any event that is not (n,4n)
if (p % event /= N_4N) cycle SCORE_LOOP
score = last_wgt
case (SCORE_ABSORPTION)
if (survival_biasing) then
! No absorption events actually occur if survival biasing is on --
@ -358,8 +346,12 @@ contains
score = ONE
case default
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
! Any other score is assumed to be a MT number. Thus, we just need
! to check if it matches the MT number of the event
if (p % event_MT /= score_bin) cycle SCORE_LOOP
score = last_wgt
end select
! Add score to tally
@ -455,10 +447,15 @@ contains
integer :: j ! loop index for scoring bins
integer :: k ! loop index for nuclide bins
integer :: l ! loop index for nuclides in material
integer :: m ! loop index for reactions
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array
integer :: i_nuclide ! index in nuclides array (from bins)
integer :: i_nuc ! index in nuclides array (from material)
integer :: i_energy ! index in nuclide energy grid
integer :: score_bin ! scoring type, e.g. SCORE_FLUX
integer :: score_index ! scoring bin index
real(8) :: f ! interpolation factor
real(8) :: flux ! tracklength estimate of flux
real(8) :: score ! actual score (e.g., flux*xs)
real(8) :: atom_density ! atom density of single nuclide in atom/b-cm
@ -466,6 +463,7 @@ contains
type(TallyObject), pointer :: t => null()
type(Material), pointer :: mat => null()
type(ListInt), pointer :: curr_ptr => null()
type(Reaction), pointer :: rxn => null()
! Determine track-length estimate of flux
flux = p % wgt * distance
@ -542,56 +540,171 @@ contains
score_bin = t % score_bins(j)
if (i_nuclide > 0) then
! Determine macroscopic nuclide cross section
! ================================================================
! DETERMINE NUCLIDE CROSS SECTION
select case(score_bin)
case (SCORE_TOTAL)
! Total cross section is pre-calculated
score = micro_xs(i_nuclide) % total * &
atom_density * flux
case (SCORE_SCATTER)
! Scattering cross section is pre-calculated
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
case (SCORE_ABSORPTION)
! Absorption cross section is pre-calculated
score = micro_xs(i_nuclide) % absorption * &
atom_density * flux
case (SCORE_FISSION)
! Fission cross section is pre-calculated
score = micro_xs(i_nuclide) % fission * &
atom_density * flux
case (SCORE_NU_FISSION)
! Nu-fission cross section is pre-calculated
score = micro_xs(i_nuclide) % nu_fission * &
atom_density * flux
case (SCORE_K_FISSION)
score = micro_xs(i_nuclide) % k_fission * &
atom_density * flux
case (SCORE_EVENTS)
! For number of events, just score unity
score = ONE
case default
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
! Any other cross section has to be calculated on-the-fly. For
! cross sections that are used often (e.g. n2n, ngamma, etc. for
! depletion), it might make sense to optimize this section or
! pre-calculate cross sections
if (score_bin > 1) then
! Set default score
score = ZERO
! TODO: The following search for the matching reaction could
! be replaced by adding a dictionary on each Nuclide instance
! of the form {MT: i_reaction, ...}
REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction
! Get pointer to reaction
rxn => nuclides(i_nuclide) % reactions(m)
! Check if this is the desired MT
if (score_bin == rxn % MT) then
! Retrieve index on nuclide energy grid and interpolation
! factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
if (i_energy >= rxn % threshold) then
score = ((ONE - f) * rxn % sigma(i_energy - &
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density * flux
end if
exit REACTION_LOOP
end if
end do REACTION_LOOP
else
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
end if
end select
else
! Determine macroscopic material cross section
! ================================================================
! DETERMINE MATERIAL CROSS SECTION
select case(score_bin)
case (SCORE_FLUX)
! For flux, we need no cross section
score = flux
case (SCORE_TOTAL)
! Total cross section is pre-calculated
score = material_xs % total * flux
case (SCORE_SCATTER)
! Scattering cross section is pre-calculated
score = (material_xs % total - material_xs % absorption) * flux
case (SCORE_ABSORPTION)
! Absorption cross section is pre-calculated
score = material_xs % absorption * flux
case (SCORE_FISSION)
! Fission cross section is pre-calculated
score = material_xs % fission * flux
case (SCORE_NU_FISSION)
! Nu-fission cross section is pre-calculated
score = material_xs % nu_fission * flux
case (SCORE_K_FISSION)
score = material_xs % k_fission * flux
case (SCORE_EVENTS)
! For number of events, just score unity
score = ONE
case default
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
! Any other cross section has to be calculated on-the-fly. This
! is somewhat costly since it requires a loop over each nuclide
! in a material and each reaction in the nuclide
if (score_bin > 1) then
! Set default score
score = ZERO
! Get pointer to current material
mat => materials(p % material)
do l = 1, mat % n_nuclides
! Get atom density
atom_density = mat % atom_density(l)
! Get index in nuclides array
i_nuc = mat % nuclide(l)
! TODO: The following search for the matching reaction could
! be replaced by adding a dictionary on each Nuclide
! instance of the form {MT: i_reaction, ...}
do m = 1, nuclides(i_nuc) % n_reaction
! Get pointer to reaction
rxn => nuclides(i_nuc) % reactions(m)
! Check if this is the desired MT
if (score_bin == rxn % MT) then
! Retrieve index on nuclide energy grid and interpolation
! factor
i_energy = micro_xs(i_nuc) % index_grid
f = micro_xs(i_nuc) % interp_factor
if (i_energy >= rxn % threshold) then
score = score + ((ONE - f) * rxn % sigma(i_energy - &
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density * flux
end if
exit
end if
end do
end do
else
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
end if
end select
end if
@ -639,13 +752,17 @@ contains
integer :: i ! loop index for nuclides in material
integer :: j ! loop index for scoring bin types
integer :: m ! loop index for reactions in nuclide
integer :: i_nuclide ! index in nuclides array
integer :: score_bin ! type of score, e.g. SCORE_FLUX
integer :: score_index ! scoring bin index
integer :: i_energy ! index in nuclide energy grid
real(8) :: f ! interpolation factor
real(8) :: score ! actual scoring tally value
real(8) :: atom_density ! atom density of single nuclide in atom/b-cm
type(TallyObject), pointer :: t => null()
type(Material), pointer :: mat => null()
type(Reaction), pointer :: rxn => null()
! Get pointer to tally
t => tallies(i_tally)
@ -673,22 +790,64 @@ contains
select case(score_bin)
case (SCORE_TOTAL)
score = micro_xs(i_nuclide) % total * atom_density * flux
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * atom_density * flux
case (SCORE_ABSORPTION)
score = micro_xs(i_nuclide) % absorption * atom_density * flux
case (SCORE_FISSION)
score = micro_xs(i_nuclide) % fission * atom_density * flux
case (SCORE_NU_FISSION)
score = micro_xs(i_nuclide) % nu_fission * atom_density * flux
case (SCORE_K_FISSION)
score = micro_xs(i_nuclide) % k_fission * atom_density * flux
case (SCORE_EVENTS)
score = ONE
case default
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
! Any other cross section has to be calculated on-the-fly. For cross
! sections that are used often (e.g. n2n, ngamma, etc. for depletion),
! it might make sense to optimize this section or pre-calculate cross
! sections
if (score_bin > 1) then
! Set default score
score = ZERO
! TODO: The following search for the matching reaction could be
! replaced by adding a dictionary on each Nuclide instance of the
! form {MT: i_reaction, ...}
REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction
! Get pointer to reaction
rxn => nuclides(i_nuclide) % reactions(m)
! Check if this is the desired MT
if (score_bin == rxn % MT) then
! Retrieve index on nuclide energy grid and interpolation factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
if (i_energy >= rxn % threshold) then
score = ((ONE - f) * rxn % sigma(i_energy - &
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density * flux
end if
exit REACTION_LOOP
end if
end do REACTION_LOOP
else
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
end if
end select
! Determine scoring bin index based on what the index of the nuclide
@ -704,7 +863,7 @@ contains
end do NUCLIDE_LOOP
! ==========================================================================
! SCORE ALL INDIVIDUAL NUCLIDE REACTION RATES
! SCORE TOTAL MATERIAL REACTION RATES
! Loop over score types for each bin
MATERIAL_SCORE_LOOP: do j = 1, t % n_score_bins
@ -715,23 +874,78 @@ contains
select case(score_bin)
case (SCORE_FLUX)
score = flux
case (SCORE_TOTAL)
score = material_xs % total * flux
case (SCORE_SCATTER)
score = (material_xs % total - material_xs % absorption) * flux
case (SCORE_ABSORPTION)
score = material_xs % absorption * flux
case (SCORE_FISSION)
score = material_xs % fission * flux
case (SCORE_NU_FISSION)
score = material_xs % nu_fission * flux
case (SCORE_K_FISSION)
score = material_xs % k_fission * flux
case (SCORE_EVENTS)
score = ONE
case default
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
! Any other cross section has to be calculated on-the-fly. This is
! somewhat costly since it requires a loop over each nuclide in a
! material and each reaction in the nuclide
if (score_bin > 1) then
! Set default score
score = ZERO
! Get pointer to current material
mat => materials(p % material)
do i = 1, mat % n_nuclides
! Get atom density
atom_density = mat % atom_density(i)
! Get index in nuclides array
i_nuclide = mat % nuclide(i)
! TODO: The following search for the matching reaction could
! be replaced by adding a dictionary on each Nuclide
! instance of the form {MT: i_reaction, ...}
do m = 1, nuclides(i_nuclide) % n_reaction
! Get pointer to reaction
rxn => nuclides(i_nuclide) % reactions(m)
! Check if this is the desired MT
if (score_bin == rxn % MT) then
! Retrieve index on nuclide energy grid and interpolation
! factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
if (i_energy >= rxn % threshold) then
score = score + ((ONE - f) * rxn % sigma(i_energy - &
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density * flux
end if
exit
end if
end do
end do
else
message = "Invalid score type on tally " // to_str(t % id) // "."
call fatal_error()
end if
end select
! Determine scoring bin index based on what the index of the nuclide
@ -785,6 +999,7 @@ contains
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
type(Material), pointer :: mat => null()
type(LocalCoord), pointer :: coord => null()
t => tallies(i_tally)
t % matching_bins = 1
@ -832,9 +1047,15 @@ contains
case (FILTER_CELL)
! determine next cell bin
! TODO: Account for cells in multiple levels when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
p % coord % cell, i_tally)
coord => p % coord0
do while(associated(coord))
position(FILTER_CELL) = 0
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
coord % cell, i_tally)
if (t % matching_bins(i) /= NO_BIN_FOUND) exit
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin
@ -1063,6 +1284,7 @@ contains
real(8) :: E ! particle energy
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
type(LocalCoord), pointer :: coord => null()
found_bin = .true.
t => tallies(i_tally)
@ -1090,9 +1312,15 @@ contains
case (FILTER_CELL)
! determine next cell bin
! TODO: Account for cells in multiple levels when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
p % coord % cell, i_tally)
coord => p % coord0
do while(associated(coord))
position(FILTER_CELL) = 0
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
coord % cell, i_tally)
if (t % matching_bins(i) /= NO_BIN_FOUND) exit
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin

91
src/utils/statepoint.py
View file

@ -10,12 +10,86 @@ import scipy.stats
filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn',
5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'}
score_types = {-1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter',
-5: 'scatter-n', -6: 'scatter-pn', -7: 'transport',
-8: 'diffusion', -9: 'n1n', -10: 'n2n', -11: 'n3n',
-12: 'n4n', -13: 'absorption', -14: 'fission',
-15: 'nu-fission', -16: 'k-fission', -17: 'current',
-18: 'events'}
score_types = {-1: 'flux',
-2: 'total',
-3: 'scatter',
-4: 'nu-scatter',
-5: 'scatter-n',
-6: 'scatter-pn',
-7: 'transport',
-8: 'diffusion',
-9: 'n1n',
-10: 'absorption',
-11: 'fission',
-12: 'nu-fission',
-13: 'k-fission',
-14: 'current',
-15: 'events',
1: '(n,total)',
2: '(n,elastic)',
4: '(n,level)',
11: '(n,2nd)',
16: '(n,2n)',
17: '(n,3n)',
18: '(n,fission)',
19: '(n,f)',
20: '(n,nf)',
21: '(n,2nf)',
22: '(n,na)',
23: '(n,n3a)',
24: '(n,2na)',
25: '(n,3na)',
28: '(n,np)',
29: '(n,n2a)',
30: '(n,2n2a)',
32: '(n,nd)',
33: '(n,nt)',
34: '(n,nHe-3)',
35: '(n,nd2a)',
36: '(n,nt2a)',
37: '(n,4n)',
38: '(n,3nf)',
41: '(n,2np)',
42: '(n,3np)',
44: '(n,n2p)',
45: '(n,npa)',
91: '(n,nc)',
101: '(n,disappear)',
102: '(n,gamma)',
103: '(n,p)',
104: '(n,d)',
105: '(n,t)',
106: '(n,3He)',
107: '(n,a)',
108: '(n,2a)',
109: '(n,3a)',
111: '(n,2p)',
112: '(n,pa)',
113: '(n,t2a)',
114: '(n,d2a)',
115: '(n,pd)',
116: '(n,pt)',
117: '(n,da)',
201: '(n,Xn)',
202: '(n,Xgamma)',
203: '(n,Xp)',
204: '(n,Xd)',
205: '(n,Xt)',
206: '(n,X3He)',
207: '(n,Xa)',
444: '(damage)',
649: '(n,pc)',
699: '(n,dc)',
749: '(n,tc)',
799: '(n,3Hec)',
849: '(n,tc)'}
score_types.update({MT: '(n,n' + str(MT-50) + ')' for MT in range(51,91)})
score_types.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)})
score_types.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)})
score_types.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)})
score_types.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)})
score_types.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)})
class BinaryFile(object):
def __init__(self, filename):
@ -183,10 +257,11 @@ class StatePoint(BinaryFile):
t.n_nuclides = n_nuclides
t.nuclides = self._get_int(n_nuclides)
# Read score bins
# Read score bins and scattering order
t.n_scores = self._get_int()[0]
t.scores = [score_types[j] for j in self._get_int(t.n_scores)]
t.scatt_order = [score_types[j] for j in self._get_int(t.n_scores)]
t.scatt_order = self._get_int(t.n_scores)
# Read number of user score bins
t.n_user_scores = self._get_int()[0]