crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merged with develop

Browse source
This commit is contained in:
Will Boyd 2015-09-05 17:58:17 -04:00
commit f5374974e6
189 changed files with 14626 additions and 9318 deletions
Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -60,3 +60,6 @@ data/nndc
# PyCharm project configuration files
.idea
.idea/*
# IPython notebook checkpoints
.ipynb_checkpoints

22
CMakeLists.txt
View file

@ -287,6 +287,22 @@ foreach(test ${TESTS})
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> -p ${TEST_PATH})
elseif(${test} MATCHES "test_filter_distribcell")
# Add each case for distribcell tests
add_test(NAME ${TEST_NAME}_case-1
WORKING_DIRECTORY ${TEST_PATH}/case-1
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-1)
add_test(NAME ${TEST_NAME}_case-2
WORKING_DIRECTORY ${TEST_PATH}/case-2
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-2)
add_test(NAME ${TEST_NAME}_case-3
WORKING_DIRECTORY ${TEST_PATH}/case-3
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-3)
add_test(NAME ${TEST_NAME}_case-4
WORKING_DIRECTORY ${TEST_PATH}/case-4
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-4)
# If a restart test is encounted, need to run with -r and restart file(s)
elseif(${test} MATCHES "restart")
@ -299,9 +315,9 @@ foreach(test ${TESTS})
elseif(${test} MATCHES "test_sourcepoint_restart")
set(RESTART_FILE statepoint.07.h5 source.07.h5)
elseif(${test} MATCHES "test_particle_restart_eigval")
set(RESTART_FILE particle_12_616.h5)
set(RESTART_FILE particle_9_555.h5)
elseif(${test} MATCHES "test_particle_restart_fixed")
set(RESTART_FILE particle_7_6144.h5)
set(RESTART_FILE particle_7_928.h5)
else(${test} MATCHES "test_statepoint_restart")
message(FATAL_ERROR "Restart test ${test} not recognized")
endif(${test} MATCHES "test_statepoint_restart")
@ -314,7 +330,7 @@ foreach(test ${TESTS})
elseif(${test} MATCHES "test_sourcepoint_restart")
set(RESTART_FILE statepoint.07.binary source.07.binary)
elseif(${test} MATCHES "test_particle_restart_eigval")
set(RESTART_FILE particle_12_616.binary)
set(RESTART_FILE particle_9_555.binary)
elseif(${test} MATCHES "test_particle_restart_fixed")
set(RESTART_FILE particle_7_6144.binary)
else(${test} MATCHES "test_statepoint_restart")

7
docs/source/conf.py
View file

@ -28,7 +28,8 @@ extensions = ['sphinx.ext.autodoc',
'sphinx.ext.napoleon',
'sphinx.ext.pngmath',
'sphinxcontrib.tikz',
'sphinx_numfig']
'sphinx_numfig',
'notebook_sphinxext']
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
@ -51,9 +52,9 @@ copyright = u'2011-2015, Massachusetts Institute of Technology'
# built documents.
#
# The short X.Y version.
version = "0.6"
version = "0.7"
# The full version, including alpha/beta/rc tags.
release = "0.6.2"
release = "0.7.0"
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.

2
docs/source/developers.rst
View file

@ -11,7 +11,7 @@ Active development of the OpenMC Monte Carlo code is currently led by:
* `Nick Horelik <https://github.com/nhorelik>`_
* `Adam Nelson <https://github.com/nelsonag>`_
* `Jon Walsh <https://github.com/walshjon>`_
* `Sterling Harper <https://github.com/walshjon>`_
* `Sterling Harper <https://github.com/smharper>`_
* `Will Boyd <https://github.com/wbinventor>`_
* `Benoit Forget <http://web.mit.edu/nse/people/faculty/forget.html>`_
* `Kord Smith <http://web.mit.edu/nse/people/faculty/smith.html>`_

2
docs/source/devguide/styleguide.rst
View file

@ -153,6 +153,8 @@ Avoid extraneous whitespace in the following situations:
Yes: if (variable == 2) then
No: if ( variable==2 ) then
Do not leave trailing whitespace at the end of a line.
------
Python
------

2
docs/source/index.rst
View file

@ -28,7 +28,7 @@ free to send a message to the User's Group `mailing list`_.
:maxdepth: 1
quickinstall
releasenotes/index
releasenotes
methods/index
usersguide/index
devguide/index

10
docs/source/methods/physics.rst
View file

@ -187,11 +187,11 @@ secondary photons from nuclear de-excitation are tracked in OpenMC.
------------------------
These types of reactions are just treated as inelastic scattering and as such
are subject to the same procedure as described in
:ref:`inelastic-scatter`. Rather than tracking multiple secondary neutrons, the
weight of the outgoing neutron is multiplied by the number of secondary
neutrons, e.g. for :math:`(n,2n)`, only one outgoing neutron is tracked but its
weight is doubled.
are subject to the same procedure as described in :ref:`inelastic-scatter`. For
reactions with integral multiplicity, e.g., :math:`(n,2n)`, an appropriate
number of secondary neutrons are created. For reactions that have a multiplicity
given as a function of the incoming neutron energy (which occasionally occurs
for MT=5), the weight of the outgoing neutron is multiplied by the multiplcity.
.. _fission:

2467
docs/source/pythonapi/examples/pandas-dataframes.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

11
docs/source/pythonapi/examples/pandas-dataframes.rst Normal file
View file

@ -0,0 +1,11 @@
=================
Pandas Dataframes
=================
.. only:: html
.. notebook:: pandas-dataframes.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

1657
docs/source/pythonapi/examples/tally-arithmetic.ipynb Normal file
View file

File diff suppressed because it is too large Load diff

11
docs/source/pythonapi/examples/tally-arithmetic.rst Normal file
View file

@ -0,0 +1,11 @@
================
Tally Arithmetic
================
.. only:: html
.. notebook:: tally-arithmetic.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

53
docs/source/pythonapi/index.rst
View file

@ -4,29 +4,68 @@
Python API
==========
--------
Contents
--------
OpenMC includes a rich Python API that enables programmatic pre- and
post-processing. The easiest way to begin using the API is to take a look at the
example Jupyter_ notebooks provided. However, this assumes that you are already
familiar with Python and common third-party packages such as NumPy_. If you have
never programmed in Python before, there are many good tutorials available
online. We recommend going through the modules from Codecademy_ and/or the
`Scipy lectures`_. The full API documentation serves to provide more information
on a given module or class.
**Handling nuclear data:**
.. toctree::
:maxdepth: 1
ace
**Creating input files:**
.. toctree::
:maxdepth: 1
cmfd
element
executor
filter
geometry
material
mesh
nuclide
opencg_compatible
particle_restart
plots
settings
statepoint
summary
surface
tallies
trigger
universe
**Running OpenMC:**
.. toctree::
:maxdepth: 1
executor
**Post-processing:**
.. toctree::
:maxdepth: 1
particle_restart
statepoint
summary
tallies
**Example Jupyter Notebooks:**
.. toctree::
:maxdepth: 1
examples/pandas-dataframes
examples/tally-arithmetic
.. _Jupyter: https://jupyter.org/
.. _NumPy: http://www.numpy.org/
.. _Codecademy: https://www.codecademy.com/tracks/python
.. _Scipy lectures: https://scipy-lectures.github.io/

67
docs/source/releasenotes.rst Normal file
View file

@ -0,0 +1,67 @@
.. _releasenotes:
==============================
Release Notes for OpenMC 0.7.0
==============================
-------------------
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
------------
- Complete Python API
- Python 3 compatability for all scripts
- All scripts consistently named openmc-* and installed together
- New 'distribcell' tally filter for repeated cells
- Ability to specify outer lattice universe
- XML input validation utility (openmc-validate-xml)
- Support for hexagonal lattices
- Material union energy grid method
- Tally triggers
- Remove dependence on PETSc
- Significant OpenMP performance improvements
- Support for Fortran 2008 MPI interface
- Use of Travis CI for continuous integration
- Simplifications and improvements to test suite
---------
Bug Fixes
---------
- b5f712_: Fix bug in spherical harmonics tallies
- e6675b_: Ensure all constants are double precision
- 04e2c1_: Fix potential bug in sample_nuclide routine
- 6121d9_: Fix bugs related to particle track files
- 2f0e89_: Fixes for nuclide specification in tallies
.. _b5f712: https://github.com/mit-crpg/openmc/commit/b5f712
.. _e6675b: https://github.com/mit-crpg/openmc/commit/e6675b
.. _04e2c1: https://github.com/mit-crpg/openmc/commit/04e2c1
.. _6121d9: https://github.com/mit-crpg/openmc/commit/6121d9
.. _2f0e89: https://github.com/mit-crpg/openmc/commit/2f0e89
------------
Contributors
------------
This release contains new contributions from the following people:
- `Will Boyd <wbinventor@gmail.com>`_
- `Matt Ellis <mellis13@mit.edu>`_
- `Sterling Harper <sterlingmharper@mit.edu>`_
- `Bryan Herman <bherman@mit.edu>`_
- `Nicholas Horelik <nicholas.horelik@gmail.com>`_
- `Colin Josey <cjosey@mit.edu>`_
- `William Lyu <PaleNeutron@users.noreply.github.com>`_
- `Adam Nelson <nelsonag@umich.edu>`_
- `Paul Romano <paul.k.romano@gmail.com>`_
- `Anthony Scopatz <scopatz@gmail.com>`_
- `Jon Walsh <walshjon@mit.edu>`_

25
docs/source/releasenotes/index.rst
View file

@ -1,25 +0,0 @@
.. _releasenotes:
=============
Release Notes
=============
The release notes for OpenMC give a list of system requirements, new features,
bugs fixed, and known issues for each successive release.
.. toctree::
:maxdepth: 1
notes_0.6.2
notes_0.6.1
notes_0.6.0
notes_0.5.4
notes_0.5.3
notes_0.5.2
notes_0.5.1
notes_0.5.0
notes_0.4.4
notes_0.4.3
notes_0.4.2
notes_0.4.1
notes_0.4.0

36
docs/source/releasenotes/notes_0.4.0.rst
View file

@ -1,36 +0,0 @@
.. _notes_0.4.0:
==============================
Release Notes for OpenMC 0.4.0
==============================
-------------------
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 as well as
Mac OS X. However, it has not been tested yet on any releases of Microsoft
Windows. Memory requirements will vary depending on the size of the problem at
hand (mostly on the number of nuclides in the problem).
------------
New Features
------------
- The probability table method for treatment of energy self-shielding in the
unresolved resonance range has been implemented and is now turned on by
default.
- Calculation of Shannon entropy for assessing convergence of the fission source
distribution.
- Ability to compile with the PGI Fortran compiler.
- Ability to run on IBM BlueGene/P machines.
- Completely rewrote how nested universes are handled. Geometry is now much more
robust.
---------
Bug Fixes
---------
- Many geometry errors have been fixed. The Monte Carlo performance benchmark
can now be successfully run in OpenMC.

55
docs/source/releasenotes/notes_0.4.1.rst
View file

@ -1,55 +0,0 @@
.. _notes_0.4.1:
==============================
Release Notes for OpenMC 0.4.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 as well as
Mac OS X. However, it has not been tested yet on any releases of Microsoft
Windows. Memory requirements will vary depending on the size of the problem at
hand (mostly on the number of nuclides in the problem).
------------
New Features
------------
- A batching method has been implemented so that statistics can be calculated
based on multiple generations instead of a single generation. This can help to
overcome problems with underpredicted variance in problems where there is
correlation between successive fission source iterations.
- Users now have the option to select a non-unionized energy grid for problems
with many nuclides where the use of a unionized grid is not feasible.
- Improved plotting capability (Nick Horelik). The plotting input is now in
``plots.xml`` instead of ``plot.xml``.
- Added multiple estimators for k-effective and added a global tally for
leakage.
- Moved cross section-related output into cross_sections.out.
- Improved timing capabilities.
- Can now use more than 2**31 - 1 particles per generation.
- Improved fission bank synchronization method. This also necessitated changing
the source bank to be of type Bank rather than of type Particle.
- Added HDF5 output (not complete yet).
- Major changes to tally implementation.
---------
Bug Fixes
---------
- `b206a8`_: Fixed subtle error in the sampling of energy distributions.
- `800742`_: Fixed error in sampling of angle and rotating angles.
- `a07c08`_: Fixed bug in linear-linear interpolation during sampling energy.
- `a75283`_: Fixed energy and energyout tally filters to support many bins.
- `95cfac`_: Fixed error in cell neighbor searches.
- `83a803`_: Fixed bug related to probability tables.
.. _b206a8: https://github.com/mit-crpg/openmc/commit/b206a8
.. _800742: https://github.com/mit-crpg/openmc/commit/800742
.. _a07c08: https://github.com/mit-crpg/openmc/commit/a07c08
.. _a75283: https://github.com/mit-crpg/openmc/commit/a75283
.. _95cfac: https://github.com/mit-crpg/openmc/commit/95cfac
.. _83a803: https://github.com/mit-crpg/openmc/commit/83a803

56
docs/source/releasenotes/notes_0.4.2.rst
View file

@ -1,56 +0,0 @@
.. _notes_0.4.2:
==============================
Release Notes for OpenMC 0.4.2
==============================
-------------------
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
------------
- Ability to specify void materials.
- Option to not reduce tallies across processors at end of each batch.
- Uniform fission site method for reducing variance on local tallies.
- Reading/writing binary source files.
- Added more messages for <trace> or high verbosity.
- Estimator for diffusion coefficient.
- Ability to specify 'point' source type.
- Ability to change random number seed.
- Users can now specify units='sum' on a <density> tag. This tells the code that
the total material density is the sum of the atom fractions listed for each
nuclide on the material.
---------
Bug Fixes
---------
- a27f8f_: Fixed runtime error bug when using Intel compiler with DEBUG on.
- afe121_: Fixed minor bug in fission bank algorithms.
- e0968e_: Force re-evaluation of cross-sections when each particle is born.
- 298db8_: Fixed bug in surface currents when using energy-in filter.
- 2f3bbe_: Fixed subtle bug in S(a,b) cross section calculation.
- 671f30_: Fixed surface currents on mesh not encompassing geometry.
- b2c40e_: Fixed bug in incoming energy filter for track-length tallies.
- 5524fd_: Mesh filter now works with track-length tallies.
- d050c7_: Added Bessel's correction to make estimate of variance unbiased.
- 2a5b9c_: Fixed regression in plotting.
.. _a27f8f: https://github.com/mit-crpg/openmc/commit/a27f8f
.. _afe121: https://github.com/mit-crpg/openmc/commit/afe121
.. _e0968e: https://github.com/mit-crpg/openmc/commit/e0968e
.. _298db8: https://github.com/mit-crpg/openmc/commit/298db8
.. _2f3bbe: https://github.com/mit-crpg/openmc/commit/2f3bbe
.. _671f30: https://github.com/mit-crpg/openmc/commit/671f30
.. _b2c40e: https://github.com/mit-crpg/openmc/commit/b2c40e
.. _5524fd: https://github.com/mit-crpg/openmc/commit/5524fd
.. _d050c7: https://github.com/mit-crpg/openmc/commit/d050c7
.. _2a5b9c: https://github.com/mit-crpg/openmc/commit/2a5b9c

53
docs/source/releasenotes/notes_0.4.3.rst
View file

@ -1,53 +0,0 @@
.. _notes_0.4.3:
==============================
Release Notes for OpenMC 0.4.3
==============================
-------------------
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
------------
- Option to report confidence intervals for tally results.
- Rotation and translation for filled cells.
- Ability to explicitly specify <estimator> for tallies.
- Ability to store state points and use them to restart runs.
- Fixed source calculations (no subcritical multiplication however).
- Expanded options for external source distribution.
- Ability to tally reaction rates for individual nuclides within a material.
- Reduced memory usage by removing redundant storage or some cross-sections.
- 3bd35b_: Log-log interpolation for URR probability tables.
- Support to specify labels on tallies (nelsonag_).
---------
Bug Fixes
---------
- 33f29a_: Handle negative values in probability table.
- 1c472d_: Fixed survival biasing with probability tables.
- 3c6e80_: Fixed writing tallies with no filters.
- 460ef1_: Invalid results for duplicate tallies.
- 0069d5_: Fixed bug with 0 inactive batches.
- 7af2cf_: Fixed bug in score_analog_tallies.
- 85a60e_: Pick closest angular distribution for law 61.
- 3212f5_: Fixed issue with blank line at beginning of XML files.
.. _nelsonag: https://github.com/nelsonag
.. _33f29a: https://github.com/mit-crpg/openmc/commit/33f29a
.. _1c472d: https://github.com/mit-crpg/openmc/commit/1c472d
.. _3c6e80: https://github.com/mit-crpg/openmc/commit/3c6e80
.. _3bd35b: https://github.com/mit-crpg/openmc/commit/3bd35b
.. _0069d5: https://github.com/mit-crpg/openmc/commit/0069d5
.. _7af2cf: https://github.com/mit-crpg/openmc/commit/7af2cf
.. _460ef1: https://github.com/mit-crpg/openmc/commit/460ef1
.. _85a60e: https://github.com/mit-crpg/openmc/commit/85a60e
.. _3212f5: https://github.com/mit-crpg/openmc/commit/3212f5

45
docs/source/releasenotes/notes_0.4.4.rst
View file

@ -1,45 +0,0 @@
.. _notes_0.4.4:
==============================
Release Notes for OpenMC 0.4.4
==============================
-------------------
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
------------
- Ability to write state points when using <no_reduce>.
- Real-time XML validation in GNU Emacs with RELAX NG schemata.
- Writing state points every n batches with <state_point interval="..." />
- Suppress creation of summary.out and cross_sections.out by default with option
to turn them on with <output> tag in settings.xml file.
- Ability to create HDF5 state points.
- Binary source file is now part of state point file by default.
- Enhanced state point usage and added state point Python scripts.
- Turning confidence intervals on affects k-effective.
- Option to specify <upper_right> for tally meshes.
---------
Bug Fixes
---------
- 4654ee_: Fixed plotting with void cells.
- 7ee461_: Fixed bug with multi-line input using type='word'.
- 792eb3_: Fixed degrees of freedom for confidence intervals.
- 7fd617_: Fixed bug with restart runs in parallel.
- dc4a8f_: Fixed bug with fixed source restart runs.
.. _4654ee: https://github.com/mit-crpg/openmc/commit/4654ee
.. _7ee461: https://github.com/mit-crpg/openmc/commit/7ee461
.. _792eb3: https://github.com/mit-crpg/openmc/commit/792eb3
.. _7fd617: https://github.com/mit-crpg/openmc/commit/7fd617
.. _dc4a8f: https://github.com/mit-crpg/openmc/commit/dc4a8f

52
docs/source/releasenotes/notes_0.5.0.rst
View file

@ -1,52 +0,0 @@
.. _notes_0.5.0:
==============================
Release Notes for OpenMC 0.5.0
==============================
-------------------
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
------------
- All user input options that formerly accepted "off" or "on" should now be
"false" or "true" (the proper XML schema datatype).
- The <criticality> element is deprecated and was replaced with <eigenvalue>.
- Added 'events' score that returns number of events that scored to a tally.
- Restructured tally filter implementation and user input.
- Source convergence acceleration via CMFD (implemented with PETSc).
- Ability to read source files in parallel when number of particles is greater
than that number of source sites.
- Cone surface types.
---------
Bug Fixes
---------
- 737b90_: Coincident surfaces from separate universes / particle traveling
tangent to a surface.
- a819b4_: Output of surface neighbors in summary.out file.
- b11696_: Reading long attribute lists in XML input.
- 2bd46a_: Search for tallying nuclides when no default_xs specified.
- 7a1f08_: Fix word wrapping when writing messages.
- c0e3ec_: Prevent underflow when compiling with MPI=yes and DEBUG=yes.
- 6f8d9d_: Set default tally labels.
- 6a3a5e_: Fix problem with corner-crossing in lattices.
.. _737b90: https://github.com/mit-crpg/openmc/commit/737b90
.. _a819b4: https://github.com/mit-crpg/openmc/commit/a819b4
.. _b11696: https://github.com/mit-crpg/openmc/commit/b11696
.. _2bd46a: https://github.com/mit-crpg/openmc/commit/2bd46a
.. _7a1f08: https://github.com/mit-crpg/openmc/commit/7a1f08
.. _c0e3ec: https://github.com/mit-crpg/openmc/commit/c0e3ec
.. _6f8d9d: https://github.com/mit-crpg/openmc/commit/6f8d9d
.. _6a3a5e: https://github.com/mit-crpg/openmc/commit/6a3a5e

45
docs/source/releasenotes/notes_0.5.1.rst
View file

@ -1,45 +0,0 @@
.. _notes_0.5.1:
==============================
Release Notes for OpenMC 0.5.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
------------
- Absorption and combined estimators for k-effective.
- Natural elements can now be specified in materials using <element> rather than
<nuclide>.
- Support for multiple S(a,b) tables in a single material (e.g. BeO).
- Test suite using Python nosetests.
- Proper install capability with 'make install'.
- Lattices can now be 2 or 3 dimensions.
- New scatter-PN score type.
- New kappa-fission score type.
- Ability to tally any reaction by specifying MT.
---------
Bug Fixes
---------
- 94103e_: Two checks for outgoing energy filters.
- e77059_: Fix reaction name for MT=849.
- b0fe88_: Fix distance to surface for cones.
- 63bfd2_: Fix tracklength tallies with cell filter and universes.
- 88daf7_: Fix analog tallies with survival biasing.
.. _94103e: https://github.com/mit-crpg/openmc/commit/94103e
.. _e77059: https://github.com/mit-crpg/openmc/commit/e77059
.. _b0fe88: https://github.com/mit-crpg/openmc/commit/b0fe88
.. _63bfd2: https://github.com/mit-crpg/openmc/commit/63bfd2
.. _88daf7: https://github.com/mit-crpg/openmc/commit/88daf7

57
docs/source/releasenotes/notes_0.5.2.rst
View file

@ -1,57 +0,0 @@
.. _notes_0.5.2:
==============================
Release Notes for OpenMC 0.5.2
==============================
-------------------
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
------------
- Python script for mesh tally plotting
- Isotopic abundances based on IUPAC 2009 when using <element>
- Particle restart files for debugging
- Code will abort after certain number of lost particles (defaults to 10)
- Region outside lattice can be filled with material (void by default)
- 3D voxel plots
- Full HDF5/PHDF5 support (including support in statepoint.py)
- Cell overlap checking with -g command line flag (or when plotting)
---------
Bug Fixes
---------
- 7632f3_: Fixed bug in statepoint.py for multiple generations per batch.
- f85ac4_: Fix infinite loop bug in error module.
- 49c36b_: Don't convert surface ids if surface filter is for current tallies.
- 5ccc78_: Fix bug in reassignment of bins for mesh filter.
- b1f52f_: Fixed bug in plot color specification.
- eae7e5_: Fixed many memory leaks.
- 10c1cc_: Minor CMFD fixes.
- afdb50_: Add compatibility for XML comments without whitespace.
- a3c593_: Fixed bug in use of free gas scattering for H-1.
- 3a66e3_: Fixed bug in 2D mesh tally implementation.
- ab0793_: Corrected PETSC_NULL references to their correct types.
- 182ebd_: Use analog estimator with energyout filter.
.. _7632f3: https://github.com/mit-crpg/openmc/commit/7632f3
.. _f85ac4: https://github.com/mit-crpg/openmc/commit/f85ac4
.. _49c36b: https://github.com/mit-crpg/openmc/commit/49c36b
.. _5ccc78: https://github.com/mit-crpg/openmc/commit/5ccc78
.. _b1f52f: https://github.com/mit-crpg/openmc/commit/b1f52f
.. _eae7e5: https://github.com/mit-crpg/openmc/commit/eae7e5
.. _10c1cc: https://github.com/mit-crpg/openmc/commit/10c1cc
.. _afdb50: https://github.com/mit-crpg/openmc/commit/afdb50
.. _a3c593: https://github.com/mit-crpg/openmc/commit/a3c593
.. _3a66e3: https://github.com/mit-crpg/openmc/commit/3a66e3
.. _ab0793: https://github.com/mit-crpg/openmc/commit/ab0793
.. _182ebd: https://github.com/mit-crpg/openmc/commit/182ebd

49
docs/source/releasenotes/notes_0.5.3.rst
View file

@ -1,49 +0,0 @@
.. _notes_0.5.3:
==============================
Release Notes for OpenMC 0.5.3
==============================
-------------------
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
------------
- Output interface enhanced to allow multiple files handles to be opened
- Particle restart file linked to output interface
- Particle restarts and state point restarts are both identified with the -r
command line flag.
- Particle instance no longer global, passed to all physics routines
- Physics routines refactored to rely less on global memory, more arguments
passed in
- CMFD routines refactored and now can compute dominance ratio on the fly
- PETSc 3.4.2 or higher must be used and compiled with fortran datatype support
- Memory leaks fixed except for ones from xml-fortran package
- Test suite enhanced to test output with different compiler options
- Description of OpenMC development workflow added
- OpenMP shared-memory parallelism added
- Special run mode --tallies removed.
---------
Bug Fixes
---------
- 2b1e8a_: Normalize direction vector after reflecting particle.
- 5853d2_: Set blank default for cross section listing alias.
- e178c7_: Fix infinite loop with words greater than 80 characters in write_message.
- c18a6e_: Check for valid secondary mode on S(a,b) tables.
- 82c456_: Fix bug where last process could have zero particles.
.. _2b1e8a: https://github.com/mit-crpg/openmc/commit/2b1e8a
.. _5853d2: https://github.com/mit-crpg/openmc/commit/5853d2
.. _e178c7: https://github.com/mit-crpg/openmc/commit/e178c7
.. _c18a6e: https://github.com/mit-crpg/openmc/commit/c18a6e
.. _82c456: https://github.com/mit-crpg/openmc/commit/82c456

64
docs/source/releasenotes/notes_0.5.4.rst
View file

@ -1,64 +0,0 @@
.. _notes_0.5.4:
==============================
Release Notes for OpenMC 0.5.4
==============================
-------------------
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
------------
- Source sites outside geometry are resampled
- XML-Fortran backend replaced by FoX XML
- Ability to write particle track files
- Handle lost particles more gracefully (via particle track files)
- Multiple random number generator streams
- Mesh tally plotting utility converted to use Tkinter rather than PyQt
- Script added to download ACE data from NNDC
- Mixed ASCII/binary cross_sections.xml now allowed
- Expanded options for writing source bank
- Re-enabled ability to use source file as starting source
- S(a,b) recalculation avoided when same nuclide and S(a,b) table are accessed
---------
Bug Fixes
---------
- 32c03c_: Check for valid data in cross_sections.xml
- c71ef5_: Fix bug in statepoint.py
- 8884fb_: Check for all ZAIDs for S(a,b) tables
- b38af0_: Fix XML reading on multiple levels of input
- d28750_: Fix bug in convert_xsdir.py
- cf567c_: ENDF/B-VI data checked for compatibility
- 6b9461_: Fix p_valid sampling inside of sample_energy
.. _32c03c: https://github.com/mit-crpg/openmc/commit/32c03c
.. _c71ef5: https://github.com/mit-crpg/openmc/commit/c71ef5
.. _8884fb: https://github.com/mit-crpg/openmc/commit/8884fb
.. _b38af0: https://github.com/mit-crpg/openmc/commit/b38af0
.. _d28750: https://github.com/mit-crpg/openmc/commit/d28750
.. _cf567c: https://github.com/mit-crpg/openmc/commit/cf567c
.. _6b9461: https://github.com/mit-crpg/openmc/commit/6b9461
------------
Contributors
------------
This release contains new contributions from the following people:
- `Sterling Harper <smharper@mit.edu>`_
- `Bryan Herman <bherman@mit.edu>`_
- `Nick Horelik <nhorelik@mit.edu>`_
- `Adam Nelson <nelsonag@umich.edu>`_
- `Paul Romano <paul.k.romano@gmail.com>`_
- `Tuomas Viitanen <tuomas.viitanen@vtt.fi>`_
- `Jon Walsh <walshjon@mit.edu>`_

59
docs/source/releasenotes/notes_0.6.0.rst
View file

@ -1,59 +0,0 @@
.. _notes_0.6.0:
==============================
Release Notes for OpenMC 0.6.0
==============================
-------------------
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
------------
- Legendre and spherical harmonic expansion tally scores
- CMake is now default build system
- Regression test suite based on CTests and NNDC cross sections
- FoX is now a git submodule
- Support for older cross sections (e.g. MCNP 66c)
- Progress bar for plots
- Expanded support for natural elements via <natural_elements> in settings.xml
---------
Bug Fixes
---------
- 41f7ca_: Fixed erroneous results from survival biasing
- 038736_: Fix tallies over void materials
- 46f9e8_: Check for negative values in probability tables
- d1ca35_: Fixed sampling of angular distribution
- 0291c0_: Fixed indexing error in plotting
- d7a7d0_: Fix bug with <element> specifying xs attribute
- 85b3cb_: Fix out-of-bounds error with OpenMP threading
.. _41f7ca: https://github.com/mit-crpg/openmc/commit/41f7ca
.. _038736: https://github.com/mit-crpg/openmc/commit/038736
.. _46f9e8: https://github.com/mit-crpg/openmc/commit/46f9e8
.. _d1ca35: https://github.com/mit-crpg/openmc/commit/d1ca35
.. _0291c0: https://github.com/mit-crpg/openmc/commit/0291c0
.. _d7a7d0: https://github.com/mit-crpg/openmc/commit/d7a7d0
.. _85b3cb: https://github.com/mit-crpg/openmc/commit/85b3cb
------------
Contributors
------------
This release contains new contributions from the following people:
- `Sterling Harper <smharper@mit.edu>`_
- `Bryan Herman <bherman@mit.edu>`_
- `Nick Horelik <nhorelik@mit.edu>`_
- `Adam Nelson <nelsonag@umich.edu>`_
- `Paul Romano <paul.k.romano@gmail.com>`_
- `Jon Walsh <walshjon@mit.edu>`_

65
docs/source/releasenotes/notes_0.6.1.rst
View file

@ -1,65 +0,0 @@
.. _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 (CMFD) 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
- 93e482_: Check for negative values in probability tables
.. _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
.. _93e482: https://github.com/mit-crpg/openmc/commit/93e482
------------
Contributors
------------
This release contains new contributions from the following people:
- `Sterling Harper <smharper@mit.edu>`_
- `Bryan Herman <bherman@mit.edu>`_
- `Adam Nelson <nelsonag@umich.edu>`_
- `Paul Romano <paul.k.romano@gmail.com>`_
- `Jon Walsh <walshjon@mit.edu>`_
- `Will Boyd <wbinventor@gmail.com>`_

58
docs/source/releasenotes/notes_0.6.2.rst
View file

@ -1,58 +0,0 @@
.. _notes_0.6.2:
==============================
Release Notes for OpenMC 0.6.2
==============================
-------------------
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
------------
- Meshline plotting capability
- Support for plotting cells/materials on middle universe levels
- Ability to model cells with no surfaces
- Compatibility with PETSc 3.5
- Compatability with OpenMPI 1.7/1.8
- Improved overall performance via logarithmic-mapped energy grid search
- Improved multi-threaded performance with atomic operations
- Support for fixed source problems with fissionable materials
---------
Bug Fixes
---------
- 26fb93_: Fix problem with -t, --track command-line flag
- 2f07c0_: Improved evaporation spectrum algorithm
- e6abb9_: Fix segfault when tallying in a void material
- 291b45_: Handle metastable nuclides in NNDC data and multiplicities in MT=5 data
.. _26fb93: https://github.com/mit-crpg/openmc/commit/26fb93
.. _2f07c0: https://github.com/mit-crpg/openmc/commit/2f07c0
.. _e6abb9: https://github.com/mit-crpg/openmc/commit/e6abb9
.. _291b45: https://github.com/mit-crpg/openmc/commit/291b45
------------
Contributors
------------
This release contains new contributions from the following people:
- `Will Boyd <wbinventor@gmail.com>`_
- `Matt Ellis <mellis13@mit.edu>`_
- `Sterling Harper <smharper@mit.edu>`_
- `Bryan Herman <bherman@mit.edu>`_
- `Nicholas Horelik <nicholas.horelik@gmail.com>`_
- `Anton Leontiev <bunder@t-25.ru>`_
- `Adam Nelson <nelsonag@umich.edu>`_
- `Paul Romano <paul.k.romano@gmail.com>`_
- `Jon Walsh <walshjon@mit.edu>`_
- `John Xia <john.danger.xia@gmail.com>`_

68
docs/sphinxext/LICENSE Normal file
View file

@ -0,0 +1,68 @@
The file notebook_sphinxext.py was derived from code in PyNE and yt.
PyNE has the following license:
-------------------------------------------------------------------------------
Copyright 2011-2015, the PyNE Development Team. All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are
permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this list of
conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice, this list
of conditions and the following disclaimer in the documentation and/or other materials
provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE PYNE DEVELOPMENT TEAM ``AS IS'' AND ANY EXPRESS OR IMPLIED
WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
The views and conclusions contained in the software and documentation are those of the
authors and should not be interpreted as representing official policies, either expressed
or implied, of the stakeholders of the PyNE project or the employers of PyNE developers.
-------------------------------------------------------------------------------
yt has the following license:
-------------------------------------------------------------------------------
yt is licensed under the terms of the Modified BSD License (also known as New
or Revised BSD), as follows:
Copyright (c) 2013-, yt Development Team
Copyright (c) 2006-2013, Matthew Turk <matthewturk@gmail.com>
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this
list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
Neither the name of the yt Development Team nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-------------------------------------------------------------------------------

120
docs/sphinxext/notebook_sphinxext.py Normal file
View file

@ -0,0 +1,120 @@
import sys
import os.path
import re
import time
from docutils import io, nodes, statemachine, utils
try:
from docutils.utils.error_reporting import ErrorString # the new way
except ImportError:
from docutils.error_reporting import ErrorString # the old way
from docutils.parsers.rst import Directive, convert_directive_function
from docutils.parsers.rst import directives, roles, states
from docutils.parsers.rst.roles import set_classes
from docutils.transforms import misc
try:
from IPython.nbconver.exporters import html
except ImportError:
from IPython.nbconvert import html
class Notebook(Directive):
"""Use nbconvert to insert a notebook into the environment.
This is based on the Raw directive in docutils
"""
required_arguments = 1
optional_arguments = 0
final_argument_whitespace = True
option_spec = {}
has_content = False
def run(self):
# check if raw html is supported
if not self.state.document.settings.raw_enabled:
raise self.warning('"%s" directive disabled.' % self.name)
# set up encoding
attributes = {'format': 'html'}
encoding = self.options.get(
'encoding', self.state.document.settings.input_encoding)
e_handler = self.state.document.settings.input_encoding_error_handler
# get path to notebook
source_dir = os.path.dirname(
os.path.abspath(self.state.document.current_source))
nb_path = os.path.normpath(os.path.join(source_dir,
self.arguments[0]))
nb_path = utils.relative_path(None, nb_path)
# convert notebook to html
exporter = html.HTMLExporter(template_file='full')
output, resources = exporter.from_filename(nb_path)
header = output.split('<head>', 1)[1].split('</head>',1)[0]
body = output.split('<body>', 1)[1].split('</body>',1)[0]
# add HTML5 scoped attribute to header style tags
header = header.replace('<style', '<style scoped="scoped"')
header = header.replace('body {\n overflow: visible;\n padding: 8px;\n}\n',
'')
header = header.replace("code,pre{", "code{")
# Filter out styles that conflict with the sphinx theme.
filter_strings = [
'navbar',
'body{',
'alert{',
'uneditable-input{',
'collapse{',
]
filter_strings.extend(['h%s{' % (i+1) for i in range(6)])
line_begin = [
'pre{',
'p{margin'
]
filterfunc = lambda x: not any([s in x for s in filter_strings])
header_lines = filter(filterfunc, header.split('\n'))
filterfunc = lambda x: not any([x.startswith(s) for s in line_begin])
header_lines = filter(filterfunc, header_lines)
header = '\n'.join(header_lines)
# concatenate raw html lines
lines = ['<div class="ipynotebook">']
lines.append(header)
lines.append(body)
lines.append('</div>')
text = '\n'.join(lines)
# add dependency
self.state.document.settings.record_dependencies.add(nb_path)
attributes['source'] = nb_path
# create notebook node
nb_node = notebook('', text, **attributes)
(nb_node.source, nb_node.line) = \
self.state_machine.get_source_and_line(self.lineno)
return [nb_node]
class notebook(nodes.raw):
pass
def visit_notebook_node(self, node):
self.visit_raw(node)
def depart_notebook_node(self, node):
self.depart_raw(node)
def setup(app):
app.add_node(notebook,
html=(visit_notebook_node, depart_notebook_node))
app.add_directive('notebook', Notebook)

150
openmc/checkvalue.py
View file

@ -1,3 +1,36 @@
from collections import Iterable
from numbers import Integral, Real
import numpy as np
def _isinstance(value, expected_type):
"""A Numpy-aware replacement for isinstance
This function will be obsolete when Numpy v. >= 1.9 is established.
"""
# Declare numpy numeric types.
np_ints = (np.int_, np.intc, np.intp, np.int8, np.int16, np.int32, np.int64,
np.uint8, np.uint16, np.uint32, np.uint64)
np_floats = (np.float_, np.float16, np.float32, np.float64)
# Include numpy integers, if necessary.
if type(expected_type) is tuple:
if Integral in expected_type:
expected_type = expected_type + np_ints
elif expected_type is Integral:
expected_type = (Integral, ) + np_ints
# Include numpy floats, if necessary.
if type(expected_type) is tuple:
if Real in expected_type:
expected_type = expected_type + np_floats
elif expected_type is Real:
expected_type = (Real, ) + np_floats
# Now, make the instance check.
return isinstance(value, expected_type)
def check_type(name, value, expected_type, expected_iter_type=None):
"""Ensure that an object is of an expected type. Optionally, if the object is
iterable, check that each element is of a particular type.
@ -16,20 +49,97 @@ def check_type(name, value, expected_type, expected_iter_type=None):
"""
if not isinstance(value, expected_type):
msg = 'Unable to set {0} to {1} which is not of type {2}'.format(
if not _isinstance(value, expected_type):
msg = 'Unable to set "{0}" to "{1}" which is not of type "{2}"'.format(
name, value, expected_type.__name__)
raise ValueError(msg)
if expected_iter_type:
for item in value:
if not isinstance(item, expected_iter_type):
msg = 'Unable to set {0} to {1} since each item must be ' \
'of type {2}'.format(name, value,
if not _isinstance(item, expected_iter_type):
msg = 'Unable to set "{0}" to "{1}" since each item must be ' \
'of type "{2}"'.format(name, value,
expected_iter_type.__name__)
raise ValueError(msg)
def check_iterable_type(name, value, expected_type, min_depth=1, max_depth=1):
"""Ensure that an object is an iterable containing an expected type.
Parameters
----------
name : str
Description of value being checked
value : Iterable
Iterable, possibly of other iterables, that should ultimately contain
the expected type
expected_type : type
type that the iterable should contain
min_depth : int
The minimum number of layers of nested iterables there should be before
reaching the ultimately contained items
max_depth : int
The maximum number of layers of nested iterables there should be before
reaching the ultimately contained items
"""
# Initialize the tree at the very first item.
tree = [value]
index = [0]
# Traverse the tree.
while index[0] != len(tree[0]):
# If we are done with this level of the tree, go to the next branch on
# the level above this one.
if index[-1] == len(tree[-1]):
del index[-1]
del tree[-1]
index[-1] += 1
continue
# Get a string representation of the current index in case we raise an
# exception.
form = '[' + '{:d}, '*(len(index)-1) + '{:d}]'
ind_str = form.format(*index)
# What is the current item we are looking at?
current_item = tree[-1][index[-1]]
# If this item is of the expected type, then we've reached the bottom
# level of this branch.
if _isinstance(current_item, expected_type):
# Is this deep enough?
if len(tree) < min_depth:
msg = 'Error setting "{0}": The item at {1} does not meet the '\
'minimum depth of {2}'.format(name, ind_str, min_depth)
raise ValueError(msg)
# This item is okay. Move on to the next item.
index[-1] += 1
# If this item is not of the expected type, then it's either an error or
# another level of the tree that we need to pursue deeper.
else:
if isinstance(current_item, Iterable):
# The tree goes deeper here, let's explore it.
tree.append(current_item)
index.append(0)
# But first, have we exceeded the max depth?
if len(tree) > max_depth:
msg = 'Error setting {0}: Found an iterable at {1}, items '\
'in that iterable excceed the maximum depth of {2}' \
.format(name, ind_str, max_depth)
raise ValueError(msg)
else:
# This item is completely unexpected.
msg = "Error setting {0}: Items must be of type '{1}', but " \
"item at {2} is of type '{3}'"\
.format(name, expected_type.__name__, ind_str,
type(current_item).__name__)
raise ValueError(msg)
def check_length(name, value, length_min, length_max=None):
"""Ensure that a sized object has length within a given range.
@ -49,16 +159,16 @@ def check_length(name, value, length_min, length_max=None):
if length_max is None:
if len(value) != length_min:
msg = 'Unable to set {0} to {1} since it must be of ' \
'length {2}'.format(name, value, length_min)
msg = 'Unable to set "{0}" to "{1}" since it must be of ' \
'length "{2}"'.format(name, value, length_min)
raise ValueError(msg)
elif not length_min <= len(value) <= length_max:
if length_min == length_max:
msg = 'Unable to set {0} to {1} since it must be of ' \
'length {2}'.format(name, value, length_min)
msg = 'Unable to set "{0}" to "{1}" since it must be of ' \
'length "{2}"'.format(name, value, length_min)
else:
msg = 'Unable to set {0} to {1} since it must have length ' \
'between {2} and {3}'.format(name, value, length_min,
msg = 'Unable to set "{0}" to "{1}" since it must have length ' \
'between "{2}" and "{3}"'.format(name, value, length_min,
length_max)
raise ValueError(msg)
@ -78,7 +188,7 @@ def check_value(name, value, accepted_values):
"""
if value not in accepted_values:
msg = 'Unable to set {0} to {1} since it is not in {2}'.format(
msg = 'Unable to set "{0}" to "{1}" since it is not in "{2}"'.format(
name, value, accepted_values)
raise ValueError(msg)
@ -100,13 +210,13 @@ def check_less_than(name, value, maximum, equality=False):
if equality:
if value > maximum:
msg = 'Unable to set {0} to {1} since it is greater than ' \
'{2}'.format(name, value, maximum)
msg = 'Unable to set "{0}" to "{1}" since it is greater than ' \
'"{2}"'.format(name, value, maximum)
raise ValueError(msg)
else:
if value >= maximum:
msg = 'Unable to set {0} to {1} since it is greater than ' \
'or equal to {2}'.format(name, value, maximum)
msg = 'Unable to set "{0}" to "{1}" since it is greater than ' \
'or equal to "{2}"'.format(name, value, maximum)
raise ValueError(msg)
def check_greater_than(name, value, minimum, equality=False):
@ -127,11 +237,11 @@ def check_greater_than(name, value, minimum, equality=False):
if equality:
if value < minimum:
msg = 'Unable to set {0} to {1} since it is less than ' \
'{2}'.format(name, value, minimum)
msg = 'Unable to set "{0}" to "{1}" since it is less than ' \
'"{2}"'.format(name, value, minimum)
raise ValueError(msg)
else:
if value <= minimum:
msg = 'Unable to set {0} to {1} since it is less than ' \
'or equal to {2}'.format(name, value, minimum)
msg = 'Unable to set "{0}" to "{1}" since it is less than ' \
'or equal to "{2}"'.format(name, value, minimum)
raise ValueError(msg)

8
openmc/clean_xml.py
View file

@ -82,11 +82,11 @@ def clean_xml_indentation(element, level=0):
if not element.tail or not element.tail.strip():
element.tail = i
for element in element:
clean_xml_indentation(element, level+1)
for sub_element in element:
clean_xml_indentation(sub_element, level+1)
if not element.tail or not element.tail.strip():
element.tail = i
if not sub_element.tail or not sub_element.tail.strip():
sub_element.tail = i
else:
if level and (not element.tail or not element.tail.strip()):

340
openmc/cross.py Normal file
View file

@ -0,0 +1,340 @@
from openmc import Filter, Nuclide
class CrossScore(object):
"""A special-purpose tally score used to encapsulate all combinations of two
tally's scores as a outer product for tally arithmetic.
Parameters
----------
left_score : str or CrossScore
The left score in the outer product
right_score : str or CrossScore
The right score in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
combine two tally's scores with this CrossNuclide
Attributes
----------
left_score : str or CrossScore
The left score in the outer product
right_score : str or CrossScore
The right score in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
combine two tally's scores with this CrossNuclide
"""
def __init__(self, left_score=None, right_score=None, binary_op=None):
self._left_score = None
self._right_score = None
self._binary_op = None
if left_score is not None:
self.left_score = left_score
if right_score is not None:
self.right_score = right_score
if binary_op is not None:
self.binary_op = binary_op
@property
def left_score(self):
return self._left_score
@property
def right_score(self):
return self._right_score
@property
def binary_op(self):
return self._binary_op
@left_score.setter
def left_score(self, left_score):
self._left_score = left_score
@right_score.setter
def right_score(self, right_score):
self._right_score = right_score
@binary_op.setter
def binary_op(self, binary_op):
self._binary_op = binary_op
def __eq__(self, other):
return str(other) == str(self)
def __repr__(self):
string = '({0} {1} {2})'.format(self.left_score,
self.binary_op, self.right_score)
return string
class CrossNuclide(object):
"""A special-purpose nuclide used to encapsulate all combinations of two
tally's nuclides as a outer product for tally arithmetic.
Parameters
----------
left_nuclide : Nuclide or CrossNuclide
The left nuclide in the outer product
right_nuclide : Nuclide or CrossNuclide
The right nuclide in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
combine two tally's nuclides with this CrossNuclide
Attributes
----------
left_nuclide : Nuclide or CrossNuclide
The left nuclide in the outer product
right_nuclide : Nuclide or CrossNuclide
The right nuclide in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
combine two tally's nuclides with this CrossNuclide
"""
def __init__(self, left_nuclide=None, right_nuclide=None, binary_op=None):
self._left_nuclide = None
self._right_nuclide = None
self._binary_op = None
if left_nuclide is not None:
self.left_nuclide = left_nuclide
if right_nuclide is not None:
self.right_nuclide = right_nuclide
if binary_op is not None:
self.binary_op = binary_op
@property
def left_nuclide(self):
return self._left_nuclide
@property
def right_nuclide(self):
return self._right_nuclide
@property
def binary_op(self):
return self._binary_op
@left_nuclide.setter
def left_nuclide(self, left_nuclide):
self._left_nuclide = left_nuclide
@right_nuclide.setter
def right_nuclide(self, right_nuclide):
self._right_nuclide = right_nuclide
@binary_op.setter
def binary_op(self, binary_op):
self._binary_op = binary_op
def __eq__(self, other):
return str(other) == str(self)
def __repr__(self):
string = ''
# If the Summary was linked, the left nuclide is a Nuclide object
if isinstance(self.left_nuclide, Nuclide):
string += '(' + self.left_nuclide.name
# If the Summary was not linked, the left nuclide is the ZAID
else:
string += '(' + str(self.left_nuclide)
string += ' ' + self.binary_op + ' '
# If the Summary was linked, the right nuclide is a Nuclide object
if isinstance(self.right_nuclide, Nuclide):
string += self.right_nuclide.name + ')'
# If the Summary was not linked, the right nuclide is the ZAID
else:
string += str(self.right_nuclide) + ')'
return string
class CrossFilter(object):
"""A special-purpose filter used to encapsulate all combinations of two
tally's filter bins as a outer product for tally arithmetic.
Parameters
----------
left_filter : Filter or CrossFilter
The left filter in the outer product
right_filter : Filter or CrossFilter
The right filter in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
combine two tally's filter bins with this CrossFilter
Attributes
----------
left_filter : Filter or CrossFilter
The left filter in the outer product
right_filter : Filter or CrossFilter
The right filter in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
combine two tally's filter bins with this CrossFilter
"""
def __init__(self, left_filter=None, right_filter=None, binary_op=None):
left_type = left_filter.type
right_type = right_filter.type
self.type = '({0} {1} {2})'.format(left_type, binary_op, right_type)
self._bins = {}
self._bins['left'] = left_filter.bins
self._bins['right'] = right_filter.bins
self._num_bins = left_filter.num_bins * right_filter.num_bins
self._left_filter = None
self._right_filter = None
self._binary_op = None
if left_filter is not None:
self.left_filter = left_filter
if right_filter is not None:
self.right_filter = right_filter
if binary_op is not None:
self.binary_op = binary_op
def __hash__(self):
return hash((self.type, self.bins))
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._left_filter = self.left_filter
clone._right_filter = self.right_filter
clone._type = self.type
clone._bins = self.bins
clone._num_bins = self.num_bins
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
@property
def left_filter(self):
return self._left_filter
@property
def right_filter(self):
return self._right_filter
@property
def binary_op(self):
return self._binary_op
@property
def type(self):
return self._type
@property
def bins(self):
return (self._bins['left'], self._bins['right'])
@property
def num_bins(self):
return self._num_bins
@property
def stride(self):
return self.left_filter.stride * self.right_filter.stride
@type.setter
def type(self, filter_type):
self._type = filter_type
@left_filter.setter
def left_filter(self, left_filter):
self._left_filter = left_filter
@right_filter.setter
def right_filter(self, right_filter):
self._right_filter = right_filter
@binary_op.setter
def binary_op(self, binary_op):
self._binary_op = binary_op
def __eq__(self, other):
return str(other) == str(self)
def split_filters(self):
split_filters = []
# If left Filter is not a CrossFilter, simply append to list
if isinstance(self.left_filter, Filter):
split_filters.append(self.left_filter)
# Recursively descend CrossFilter tree to collect all Filters
else:
split_filters.extend(self.left_filter.split_filters())
# If right Filter is not a CrossFilter, simply append to list
if isinstance(self.right_filter, Filter):
split_filters.append(self.right_filter)
# Recursively descend CrossFilter tree to collect all Filters
else:
split_filters.extend(self.right_filter.split_filters())
return split_filters
def get_bin_index(self, filter_bin):
"""Returns the index in the CrossFilter for some bin.
Parameters
----------
filter_bin : 2-tuple
A 2-tuple where each value corresponds to the bin of interest
in the left and right filter, respectively. A bin is the integer
ID for 'material', 'surface', 'cell', 'cellborn', and 'universe'
Filters. The bin is an integer for the cell instance ID for
'distribcell' Filters. The bin is a 2-tuple of floats for 'energy'
and 'energyout' filters corresponding to the energy boundaries of
the bin of interest. The bin is a (x,y,z) 3-tuple for 'mesh'
filters corresponding to the mesh cell of interest.
Returns
-------
filter_index : int
The index in the Tally data array for this filter bin.
"""
left_index = self.left_filter.get_bin_index(filter_bin[0])
right_index = self.right_filter.get_bin_index(filter_bin[0])
filter_index = left_index * self.right_filter.num_bins + right_index
return filter_index
def __repr__(self):
string = 'CrossFilter\n'
filter_type = '({0} {1} {2})'.format(self.left_filter.type,
self.binary_op,
self.right_filter.type)
filter_bins = '({0} {1} {2})'.format(self.left_filter.bins,
self.binary_op,
self.right_filter.bins)
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', filter_type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', filter_bins)
return string

2
openmc/executor.py
View file

@ -50,7 +50,7 @@ class Executor(object):
check_type("Executor's working directory", working_directory,
basestring)
if not os.path.isdir(working_directory):
msg = 'Unable to set Executor\'s working directory to {0} ' \
msg = 'Unable to set Executor\'s working directory to "{0}" ' \
'which does not exist'.format(working_directory)
raise ValueError(msg)

78
openmc/filter.py
View file

@ -6,7 +6,8 @@ import numpy as np
from openmc import Mesh
from openmc.constants import *
from openmc.checkvalue import check_type
from openmc.checkvalue import check_type, check_iterable_type, \
check_greater_than
class Filter(object):
"""A filter used to constrain a tally to a specific criterion, e.g. only tally
@ -42,15 +43,15 @@ class Filter(object):
def __eq__(self, filter2):
# Check type
if self._type != filter2._type:
if self.type != filter2.type:
return False
# Check number of bins
elif len(self._bins) != len(filter2._bins):
elif len(self.bins) != len(filter2.bins):
return False
# Check bin edges
elif not np.allclose(self._bins, filter2._bins):
elif not np.allclose(self.bins, filter2.bins):
return False
else:
@ -65,12 +66,12 @@ class Filter(object):
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._type = self._type
clone._bins = copy.deepcopy(self._bins, memo)
clone._num_bins = self._num_bins
clone._mesh = copy.deepcopy(self._mesh, memo)
clone._offset = self._offset
clone._stride = self._stride
clone._type = self.type
clone._bins = copy.deepcopy(self.bins, memo)
clone._num_bins = self.num_bins
clone._mesh = copy.deepcopy(self.mesh, memo)
clone._offset = self.offset
clone._stride = self.stride
memo[id(self)] = clone
@ -132,36 +133,30 @@ class Filter(object):
else:
bins = list(bins)
if self._type in ['cell', 'cellborn', 'surface', 'material',
if self.type in ['cell', 'cellborn', 'surface', 'material',
'universe', 'distribcell']:
check_iterable_type('filter bins', bins, Integral)
for edge in bins:
if not isinstance(edge, Integral):
msg = 'Unable to add bin "{0}" to a {1} Filter since ' \
'it is not an integer'.format(edge, self._type)
raise ValueError(msg)
elif edge < 0:
msg = 'Unable to add bin "{0}" to a {1} Filter since ' \
'it is negative'.format(edge, self._type)
raise ValueError(msg)
check_greater_than('filter bin', edge, 0, equality=True)
elif self._type in ['energy', 'energyout']:
for edge in bins:
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a {1} Filter ' \
msg = 'Unable to add bin edge "{0}" to a "{1}" Filter ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, self._type)
'value'.format(edge, self.type)
raise ValueError(msg)
elif edge < 0.:
msg = 'Unable to add bin edge "{0}" to a {1} Filter ' \
'since it is a negative value'.format(edge, self._type)
msg = 'Unable to add bin edge "{0}" to a "{1}" Filter ' \
'since it is a negative value'.format(edge, self.type)
raise ValueError(msg)
# Check that bin edges are monotonically increasing
for index in range(len(bins)):
if index > 0 and bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a {1} Filter ' \
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, self._type)
'increasing'.format(bins, self.type)
raise ValueError(msg)
# mesh filters
@ -180,17 +175,13 @@ class Filter(object):
raise ValueError(msg)
# If all error checks passed, add bin edges
self._bins = bins
self._bins = np.array(bins)
# FIXME
@num_bins.setter
def num_bins(self, num_bins):
if not isinstance(num_bins, Integral) or num_bins < 0:
msg = 'Unable to set the number of bins "{0}" for a {1} Filter ' \
'since it is not a positive ' \
'integer'.format(num_bins, self._type)
raise ValueError(msg)
check_type('filter num_bins', num_bins, Integral)
check_greater_than('filter num_bins', num_bins, 0, equality=True)
self._num_bins = num_bins
@mesh.setter
@ -199,7 +190,7 @@ class Filter(object):
self._mesh = mesh
self.type = 'mesh'
self.bins = self._mesh._id
self.bins = self.mesh.id
@offset.setter
def offset(self, offset):
@ -210,8 +201,8 @@ class Filter(object):
def stride(self, stride):
check_type('filter stride', stride, Integral)
if stride < 0:
msg = 'Unable to set stride "{0}" for a {1} Filter since it is a ' \
'negative value'.format(stride, self._type)
msg = 'Unable to set stride "{0}" for a "{1}" Filter since it ' \
'is a negative value'.format(stride, self.type)
raise ValueError(msg)
self._stride = stride
@ -269,14 +260,15 @@ class Filter(object):
"""
if not self.can_merge(filter):
msg = 'Unable to merge {0} with {1} filters'.format(self._type, filter._type)
msg = 'Unable to merge "{0}" with "{1}" ' \
'filters'.format(self.type, filter.type)
raise ValueError(msg)
# Create deep copy of filter to return as merged filter
merged_filter = copy.deepcopy(self)
# Merge unique filter bins
merged_bins = list(set(self._bins + filter._bins))
merged_bins = list(set(self.bins + filter.bins))
merged_filter.bins = merged_bins
merged_filter.num_bins = len(merged_bins)
@ -322,17 +314,17 @@ class Filter(object):
# Use lower energy bound to find index for energy Filters
elif self.type in ['energy', 'energyout']:
val = self.bins.index(filter_bin[0])
val = np.where(self.bins == filter_bin[0])[0][0]
filter_index = val
# Filter bins for distribcell are the "IDs" of each unique placement
# of the Cell in the Geometry (integers starting at 0)
elif self._type == 'distribcell':
elif self.type == 'distribcell':
filter_index = filter_bin
# Use ID for all other Filters (e.g., material, cell, etc.)
else:
val = self.bins.index(filter_bin)
val = np.where(self.bins == filter_bin)[0][0]
filter_index = val
except ValueError:
@ -344,7 +336,7 @@ class Filter(object):
def __repr__(self):
string = 'Filter\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self._bins)
string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offset)
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins)
string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self.offset)
return string

4
openmc/geometry.py
View file

@ -34,8 +34,8 @@ class Geometry(object):
def root_universe(self, root_universe):
check_type('root universe', root_universe, openmc.Universe)
if root_universe._id != 0:
msg = 'Unable to add root Universe {0} to Geometry since ' \
'it has ID={1} instead of ' \
msg = 'Unable to add root Universe "{0}" to Geometry since ' \
'it has ID="{1}" instead of ' \
'ID=0'.format(root_universe, root_universe._id)
raise ValueError(msg)

56
openmc/material.py
View file

@ -26,7 +26,7 @@ def reset_auto_material_id():
# Units for density supported by OpenMC
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'at/b-cm', 'at/cm3', 'sum']
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum']
# Constant for density when not needed
NO_DENSITY = 99999.
@ -122,7 +122,7 @@ class Material(object):
else:
check_type('material ID', material_id, Integral)
if material_id in MATERIAL_IDS:
msg = 'Unable to set Material ID to {0} since a Material with ' \
msg = 'Unable to set Material ID to "{0}" since a Material with ' \
'this ID was already initialized'.format(material_id)
raise ValueError(msg)
check_greater_than('material ID', material_id, 0)
@ -132,7 +132,7 @@ class Material(object):
@name.setter
def name(self, name):
check_type('name for Material ID={0}'.format(self._id),
check_type('name for Material ID="{0}"'.format(self._id),
name, basestring)
self._name = name
@ -149,12 +149,12 @@ class Material(object):
"""
check_type('the density for Material ID={0}'.format(self._id),
check_type('the density for Material ID="{0}"'.format(self._id),
density, Real)
check_value('density units', units, DENSITY_UNITS)
if density == NO_DENSITY and units is not 'sum':
msg = 'Unable to set the density Material ID={0} ' \
msg = 'Unable to set the density Material ID="{0}" ' \
'because a density must be set when not using ' \
'sum unit'.format(self._id)
raise ValueError(msg)
@ -169,8 +169,8 @@ class Material(object):
'version of openmc')
if not isinstance(filename, basestring) and filename is not None:
msg = 'Unable to add OTF material file to Material ID={0} with a ' \
'non-string name {1}'.format(self._id, filename)
msg = 'Unable to add OTF material file to Material ID="{0}" with a ' \
'non-string name "{1}"'.format(self._id, filename)
raise ValueError(msg)
self._distrib_otf_file = filename
@ -198,18 +198,18 @@ class Material(object):
"""
if not isinstance(nuclide, (openmc.Nuclide, str)):
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'non-Nuclide value {1}'.format(self._id, nuclide)
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'non-Nuclide value "{1}"'.format(self._id, nuclide)
raise ValueError(msg)
elif not isinstance(percent, Real):
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'non-floating point value {1}'.format(self._id, percent)
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'non-floating point value "{1}"'.format(self._id, percent)
raise ValueError(msg)
elif percent_type not in ['ao', 'wo', 'at/g-cm']:
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'percent type {1}'.format(self._id, percent_type)
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'percent type "{1}"'.format(self._id, percent_type)
raise ValueError(msg)
if isinstance(nuclide, openmc.Nuclide):
@ -232,7 +232,7 @@ class Material(object):
"""
if not isinstance(nuclide, openmc.Nuclide):
msg = 'Unable to remove a Nuclide {0} in Material ID={1} ' \
msg = 'Unable to remove a Nuclide "{0}" in Material ID="{1}" ' \
'since it is not a Nuclide'.format(self._id, nuclide)
raise ValueError(msg)
@ -255,18 +255,18 @@ class Material(object):
"""
if not isinstance(element, openmc.Element):
msg = 'Unable to add an Element to Material ID={0} with a ' \
'non-Element value {1}'.format(self._id, element)
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'non-Element value "{1}"'.format(self._id, element)
raise ValueError(msg)
if not isinstance(percent, Real):
msg = 'Unable to add an Element to Material ID={0} with a ' \
'non-floating point value {1}'.format(self._id, percent)
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'non-floating point value "{1}"'.format(self._id, percent)
raise ValueError(msg)
if percent_type not in ['ao', 'wo']:
msg = 'Unable to add an Element to Material ID={0} with a ' \
'percent type {1}'.format(self._id, percent_type)
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'percent type "{1}"'.format(self._id, percent_type)
raise ValueError(msg)
# Copy this Element to separate it from same Element in other Materials
@ -301,13 +301,13 @@ class Material(object):
"""
if not isinstance(name, basestring):
msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
'non-string table name {1}'.format(self._id, name)
msg = 'Unable to add an S(a,b) table to Material ID="{0}" with a ' \
'non-string table name "{1}"'.format(self._id, name)
raise ValueError(msg)
if not isinstance(xs, basestring):
msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
'non-string cross-section identifier {1}'.format(self._id, xs)
msg = 'Unable to add an S(a,b) table to Material ID="{0}" with a ' \
'non-string cross-section identifier "{1}"'.format(self._id, xs)
raise ValueError(msg)
self._sab.append((name, xs))
@ -336,7 +336,7 @@ class Material(object):
return nuclides
def _repr__(self):
def __repr__(self):
string = 'Material\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
@ -529,7 +529,7 @@ class MaterialsFile(object):
"""
if not isinstance(material, Material):
msg = 'Unable to add a non-Material {0} to the ' \
msg = 'Unable to add a non-Material "{0}" to the ' \
'MaterialsFile'.format(material)
raise ValueError(msg)
@ -546,7 +546,7 @@ class MaterialsFile(object):
"""
if not isinstance(materials, Iterable):
msg = 'Unable to create OpenMC materials.xml file from {0} which ' \
msg = 'Unable to create OpenMC materials.xml file from "{0}" which ' \
'is not iterable'.format(materials)
raise ValueError(msg)
@ -564,7 +564,7 @@ class MaterialsFile(object):
"""
if not isinstance(material, Material):
msg = 'Unable to remove a non-Material {0} from the ' \
msg = 'Unable to remove a non-Material "{0}" from the ' \
'MaterialsFile'.format(material)
raise ValueError(msg)

6
openmc/mesh.py
View file

@ -148,14 +148,14 @@ class Mesh(object):
@name.setter
def name(self, name):
check_type('name for mesh ID={0}'.format(self._id), name, basestring)
check_type('name for mesh ID="{0}"'.format(self._id), name, basestring)
self._name = name
@type.setter
def type(self, meshtype):
check_type('type for mesh ID={0}'.format(self._id),
check_type('type for mesh ID="{0}"'.format(self._id),
meshtype, basestring)
check_value('type for mesh ID={0}'.format(self._id),
check_value('type for mesh ID="{0}"'.format(self._id),
meshtype, ['rectangular', 'hexagonal'])
self._type = meshtype

40
openmc/opencg_compatible.py
View file

@ -78,7 +78,7 @@ def get_opencg_material(openmc_material):
"""
if not isinstance(openmc_material, openmc.Material):
msg = 'Unable to create an OpenCG Material from {0} ' \
msg = 'Unable to create an OpenCG Material from "{0}" ' \
'which is not an OpenMC Material'.format(openmc_material)
raise ValueError(msg)
@ -118,7 +118,7 @@ def get_openmc_material(opencg_material):
"""
if not isinstance(opencg_material, opencg.Material):
msg = 'Unable to create an OpenMC Material from {0} ' \
msg = 'Unable to create an OpenMC Material from "{0}" ' \
'which is not an OpenCG Material'.format(opencg_material)
raise ValueError(msg)
@ -165,7 +165,7 @@ def is_opencg_surface_compatible(opencg_surface):
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to check if OpenCG Surface is compatible' \
'since {0} is not a Surface'.format(opencg_surface)
'since "{0}" is not a Surface'.format(opencg_surface)
raise ValueError(msg)
if opencg_surface._type in ['x-squareprism',
@ -191,7 +191,7 @@ def get_opencg_surface(openmc_surface):
"""
if not isinstance(openmc_surface, openmc.Surface):
msg = 'Unable to create an OpenCG Surface from {0} ' \
msg = 'Unable to create an OpenCG Surface from "{0}" ' \
'which is not an OpenMC Surface'.format(openmc_surface)
raise ValueError(msg)
@ -277,7 +277,7 @@ def get_openmc_surface(opencg_surface):
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from {0} which ' \
msg = 'Unable to create an OpenMC Surface from "{0}" which ' \
'is not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
@ -335,7 +335,7 @@ def get_openmc_surface(opencg_surface):
else:
msg = 'Unable to create an OpenMC Surface from an OpenCG ' \
'Surface of type {0} since it is not a compatible ' \
'Surface of type "{0}" since it is not a compatible ' \
'Surface type in OpenMC'.format(opencg_surface._type)
raise ValueError(msg)
@ -368,7 +368,7 @@ def get_compatible_opencg_surfaces(opencg_surface):
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from {0} which ' \
msg = 'Unable to create an OpenMC Surface from "{0}" which ' \
'is not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
@ -421,7 +421,7 @@ def get_compatible_opencg_surfaces(opencg_surface):
else:
msg = 'Unable to create a compatible OpenMC Surface an OpenCG ' \
'Surface of type {0} since it already a compatible ' \
'Surface of type "{0}" since it already a compatible ' \
'Surface type in OpenMC'.format(opencg_surface._type)
raise ValueError(msg)
@ -450,7 +450,7 @@ def get_opencg_cell(openmc_cell):
"""
if not isinstance(openmc_cell, openmc.Cell):
msg = 'Unable to create an OpenCG Cell from {0} which ' \
msg = 'Unable to create an OpenCG Cell from "{0}" which ' \
'is not an OpenMC Cell'.format(openmc_cell)
raise ValueError(msg)
@ -518,17 +518,17 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
"""
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create compatible OpenMC Cell from {0} which ' \
msg = 'Unable to create compatible OpenMC Cell from "{0}" which ' \
'is not an OpenCG Cell'.format(opencg_cell)
raise ValueError(msg)
elif not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create compatible OpenMC Cell since {0} is ' \
msg = 'Unable to create compatible OpenMC Cell since "{0}" is ' \
'not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
elif halfspace not in [-1, +1]:
msg = 'Unable to create compatible Cell since {0}' \
msg = 'Unable to create compatible Cell since "{0}"' \
'is not a +/-1 halfspace'.format(halfspace)
raise ValueError(msg)
@ -626,7 +626,7 @@ def make_opencg_cells_compatible(opencg_universe):
"""
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to make compatible OpenCG Cells for {0} which ' \
msg = 'Unable to make compatible OpenCG Cells for "{0}" which ' \
'is not an OpenCG Universe'.format(opencg_universe)
raise ValueError(msg)
@ -685,7 +685,7 @@ def get_openmc_cell(opencg_cell):
"""
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create an OpenMC Cell from {0} which ' \
msg = 'Unable to create an OpenMC Cell from "{0}" which ' \
'is not an OpenCG Cell'.format(opencg_cell)
raise ValueError(msg)
@ -749,7 +749,7 @@ def get_opencg_universe(openmc_universe):
"""
if not isinstance(openmc_universe, openmc.Universe):
msg = 'Unable to create an OpenCG Universe from {0} which ' \
msg = 'Unable to create an OpenCG Universe from "{0}" which ' \
'is not an OpenMC Universe'.format(openmc_universe)
raise ValueError(msg)
@ -796,7 +796,7 @@ def get_openmc_universe(opencg_universe):
"""
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to create an OpenMC Universe from {0} which ' \
msg = 'Unable to create an OpenMC Universe from "{0}" which ' \
'is not an OpenCG Universe'.format(opencg_universe)
raise ValueError(msg)
@ -846,7 +846,7 @@ def get_opencg_lattice(openmc_lattice):
"""
if not isinstance(openmc_lattice, openmc.Lattice):
msg = 'Unable to create an OpenCG Lattice from {0} which ' \
msg = 'Unable to create an OpenCG Lattice from "{0}" which ' \
'is not an OpenMC Lattice'.format(openmc_lattice)
raise ValueError(msg)
@ -926,7 +926,7 @@ def get_openmc_lattice(opencg_lattice):
"""
if not isinstance(opencg_lattice, opencg.Lattice):
msg = 'Unable to create an OpenMC Lattice from {0} which ' \
msg = 'Unable to create an OpenMC Lattice from "{0}" which ' \
'is not an OpenCG Lattice'.format(opencg_lattice)
raise ValueError(msg)
@ -997,7 +997,7 @@ def get_opencg_geometry(openmc_geometry):
"""
if not isinstance(openmc_geometry, openmc.Geometry):
msg = 'Unable to get OpenCG geometry from {0} which is ' \
msg = 'Unable to get OpenCG geometry from "{0}" which is ' \
'not an OpenMC Geometry object'.format(openmc_geometry)
raise ValueError(msg)
@ -1037,7 +1037,7 @@ def get_openmc_geometry(opencg_geometry):
"""
if not isinstance(opencg_geometry, opencg.Geometry):
msg = 'Unable to get OpenMC geometry from {0} which is ' \
msg = 'Unable to get OpenMC geometry from "{0}" which is ' \
'not an OpenCG Geometry object'.format(opencg_geometry)
raise ValueError(msg)

12
openmc/plots.py
View file

@ -210,18 +210,18 @@ class Plot(object):
for key in col_spec:
if key < 0:
msg = 'Unable to create Plot ID={0} with col_spec ID {1} ' \
msg = 'Unable to create Plot ID="{0}" with col_spec ID "{1}" ' \
'which is less than 0'.format(self._id, key)
raise ValueError(msg)
elif not isinstance(col_spec[key], Iterable):
msg = 'Unable to create Plot ID={0} with col_spec RGB values' \
' {1} which is not iterable'.format(self._id, col_spec[key])
msg = 'Unable to create Plot ID="{0}" with col_spec RGB values' \
' "{1}" which is not iterable'.format(self._id, col_spec[key])
raise ValueError(msg)
elif len(col_spec[key]) != 3:
msg = 'Unable to create Plot ID={0} with col_spec RGB ' \
'values of length {1} since 3 values must be ' \
msg = 'Unable to create Plot ID="{0}" with col_spec RGB ' \
'values of length "{1}" since 3 values must be ' \
'input'.format(self._id, len(col_spec[key]))
raise ValueError(msg)
@ -332,7 +332,7 @@ class PlotsFile(object):
"""
if not isinstance(plot, Plot):
msg = 'Unable to add a non-Plot {0} to the PlotsFile'.format(plot)
msg = 'Unable to add a non-Plot "{0}" to the PlotsFile'.format(plot)
raise ValueError(msg)
self._plots.append(plot)

22
openmc/settings.py
View file

@ -426,28 +426,28 @@ class SettingsFile(object):
@keff_trigger.setter
def keff_trigger(self, keff_trigger):
if not isinstance(keff_trigger, dict):
msg = 'Unable to set a trigger on keff from {0} which ' \
msg = 'Unable to set a trigger on keff from "{0}" which ' \
'is not a Python dictionary'.format(keff_trigger)
raise ValueError(msg)
elif 'type' not in keff_trigger:
msg = 'Unable to set a trigger on keff from {0} which ' \
msg = 'Unable to set a trigger on keff from "{0}" which ' \
'does not have a "type" key'.format(keff_trigger)
raise ValueError(msg)
elif keff_trigger['type'] not in ['variance', 'std_dev', 'rel_err']:
msg = 'Unable to set a trigger on keff with ' \
'type {0}'.format(keff_trigger['type'])
'type "{0}"'.format(keff_trigger['type'])
raise ValueError(msg)
elif 'threshold' not in keff_trigger:
msg = 'Unable to set a trigger on keff from {0} which ' \
msg = 'Unable to set a trigger on keff from "{0}" which ' \
'does not have a "threshold" key'.format(keff_trigger)
raise ValueError(msg)
elif not isinstance(keff_trigger['threshold'], Real):
msg = 'Unable to set a trigger on keff with ' \
'threshold {0}'.format(keff_trigger['threshold'])
'threshold "{0}"'.format(keff_trigger['threshold'])
raise ValueError(msg)
self._keff_trigger = keff_trigger
@ -574,20 +574,20 @@ class SettingsFile(object):
@output.setter
def output(self, output):
if not isinstance(output, dict):
msg = 'Unable to set output to {0} which is not a Python ' \
msg = 'Unable to set output to "{0}" which is not a Python ' \
'dictionary of string keys and boolean values'.format(output)
raise ValueError(msg)
for element in output:
keys = ['summary', 'cross_sections', 'tallies', 'distribmats']
if element not in keys:
msg = 'Unable to set output to {0} which is unsupported by ' \
msg = 'Unable to set output to "{0}" which is unsupported by ' \
'OpenMC'.format(element)
raise ValueError(msg)
if not isinstance(output[element], (bool, np.bool)):
msg = 'Unable to set output for {0} to a non-boolean ' \
'value {1}'.format(element, output[element])
msg = 'Unable to set output for "{0}" to a non-boolean ' \
'value "{1}"'.format(element, output[element])
raise ValueError(msg)
self._output = output
@ -753,7 +753,7 @@ class SettingsFile(object):
def track(self, track):
check_type('track', track, Iterable, Integral)
if len(track) % 3 != 0:
msg = 'Unable to set the track to {0} since its length is ' \
msg = 'Unable to set the track to "{0}" since its length is ' \
'not a multiple of 3'.format(track)
raise ValueError(msg)
for t in zip(track[::3], track[1::3], track[2::3]):
@ -832,7 +832,7 @@ class SettingsFile(object):
len_nodemap = np.prod(self._dd_mesh_dimension)
if len(nodemap) < len_nodemap or len(nodemap) > len_nodemap:
msg = 'Unable to set DD nodemap with length {0} which ' \
msg = 'Unable to set DD nodemap with length "{0}" which ' \
'does not have the same dimensionality as the domain ' \
'mesh'.format(len(nodemap))
raise ValueError(msg)

4
openmc/statepoint.py
View file

@ -372,7 +372,7 @@ class StatePoint(object):
path='{0}{1}/n_bins'.format(subbase, j))[0]
if n_bins <= 0:
msg = 'Unable to create Filter {0} for Tally ID={2} ' \
msg = 'Unable to create Filter "{0}" for Tally ID="{1}" ' \
'since no bins were specified'.format(j, tally_key)
raise ValueError(msg)
@ -748,7 +748,7 @@ class StatePoint(object):
"""
if not isinstance(summary, openmc.summary.Summary):
msg = 'Unable to link statepoint with {0} which ' \
msg = 'Unable to link statepoint with "{0}" which ' \
'is not a Summary object'.format(summary)
raise ValueError(msg)

72
openmc/summary.py
View file

@ -411,16 +411,70 @@ class Summary(object):
if outer != -22:
lattice.outer = self.universes[outer]
# Build array of Universe pointers for the Lattice
universes = \
np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe)
# Build array of Universe pointers for the Lattice. Note that
# we need to convert between the HDF5's square array of
# (x, alpha, z) to the Python API's format of a ragged nested
# list of (z, ring, theta).
universes = []
for z in range(lattice.num_axial):
# Add a list for this axial level.
universes.append([])
x = lattice.num_rings - 1
a = 2*lattice.num_rings - 2
for r in range(lattice.num_rings - 1, 0, -1):
# Add a list for this ring.
universes[-1].append([])
for i in range(universe_ids.shape[0]):
for j in range(universe_ids.shape[1]):
for k in range(universe_ids.shape[2]):
if universe_ids[i, j, k] != -1:
universes[i, j, k] = self.get_universe_by_id(
universe_ids[i, j, k])
# Climb down the top-right.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
x += 1
a -= 1
# Climb down the right.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
a -= 1
# Climb down the bottom-right.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
x -= 1
# Climb up the bottom-left.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
x -= 1
a += 1
# Climb up the left.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
a += 1
# Climb up the top-left.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
x += 1
# Move down to the next ring.
a -= 1
# Convert the ids into Universe objects.
universes[-1][-1] = [self.get_universe_by_id(u_id)
for u_id in universes[-1][-1]]
# Handle the degenerate center ring separately.
u_id = universe_ids[z, a, x]
universes[-1].append([self.get_universe_by_id(u_id)])
# Add the universes to the lattice.
if len(pitch) == 2:
# Lattice is 3D
lattice.universes = universes
else:
# Lattice is 2D; extract the only axial level
lattice.universes = universes[0]
if offset_size > 0:
lattice.offsets = offsets

950
openmc/tallies.py
View file

File diff suppressed because it is too large Load diff

12
openmc/temp.py Normal file
View file

@ -0,0 +1,12 @@
from checkvalue import *
from checkvalue import _isinstance
import numpy as np
zs = np.zeros((2,))
print _isinstance(zs[0], Integral)
print _isinstance(zs[0], Real)
print _isinstance(zs[0], (Integral, Real))
print check_iterable_type('thing', zs, (Real, Integral))

120
openmc/universe.py
View file

@ -7,7 +7,7 @@ import sys
import numpy as np
import openmc
from openmc.checkvalue import check_type, check_length, check_greater_than
import openmc.checkvalue as cv
if sys.version_info[0] >= 3:
basestring = str
@ -111,13 +111,13 @@ class Cell(object):
self._id = AUTO_CELL_ID
AUTO_CELL_ID += 1
else:
check_type('cell ID', cell_id, Integral)
check_greater_than('cell ID', cell_id, 0)
cv.check_type('cell ID', cell_id, Integral)
cv.check_greater_than('cell ID', cell_id, 0)
self._id = cell_id
@name.setter
def name(self, name):
check_type('cell name', name, basestring)
cv.check_type('cell name', name, basestring)
self._name = name
@fill.setter
@ -126,8 +126,8 @@ class Cell(object):
if fill.strip().lower() == 'void':
self._type = 'void'
else:
msg = 'Unable to set Cell ID={0} to use a non-Material or ' \
'Universe fill {1}'.format(self._id, fill)
msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \
'Universe fill "{1}"'.format(self._id, fill)
raise ValueError(msg)
elif isinstance(fill, openmc.Material):
@ -140,27 +140,27 @@ class Cell(object):
self._type = 'lattice'
else:
msg = 'Unable to set Cell ID={0} to use a non-Material or ' \
'Universe fill {1}'.format(self._id, fill)
msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \
'Universe fill "{1}"'.format(self._id, fill)
raise ValueError(msg)
self._fill = fill
@rotation.setter
def rotation(self, rotation):
check_type('cell rotation', rotation, Iterable, Real)
check_length('cell rotation', rotation, 3)
cv.check_type('cell rotation', rotation, Iterable, Real)
cv.check_length('cell rotation', rotation, 3)
self._rotation = rotation
@translation.setter
def translation(self, translation):
check_type('cell translation', translation, Iterable, Real)
check_length('cell translation', translation, 3)
cv.check_type('cell translation', translation, Iterable, Real)
cv.check_length('cell translation', translation, 3)
self._translation = translation
@offsets.setter
def offsets(self, offsets):
check_type('cell offsets', offsets, Iterable)
cv.check_type('cell offsets', offsets, Iterable)
self._offsets = offsets
def add_surface(self, surface, halfspace):
@ -177,13 +177,13 @@ class Cell(object):
"""
if not isinstance(surface, openmc.Surface):
msg = 'Unable to add Surface {0} to Cell ID={1} since it is ' \
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" since it is ' \
'not a Surface object'.format(surface, self._id)
raise ValueError(msg)
if halfspace not in [-1, +1]:
msg = 'Unable to add Surface {0} to Cell ID={1} with halfspace ' \
'{2} since it is not +/-1'.format(surface, self._id, halfspace)
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" with halfspace ' \
'"{2}" since it is not +/-1'.format(surface, self._id, halfspace)
raise ValueError(msg)
# If the Cell does not already contain the Surface, add it
@ -201,7 +201,7 @@ class Cell(object):
"""
if not isinstance(surface, openmc.Surface):
msg = 'Unable to remove Surface {0} from Cell ID={1} since it is ' \
msg = 'Unable to remove Surface "{0}" from Cell ID="{1}" since it is ' \
'not a Surface object'.format(surface, self._id)
raise ValueError(msg)
@ -432,13 +432,13 @@ class Universe(object):
self._id = AUTO_UNIVERSE_ID
AUTO_UNIVERSE_ID += 1
else:
check_type('universe ID', universe_id, Integral)
check_greater_than('universe ID', universe_id, 0, True)
cv.check_type('universe ID', universe_id, Integral)
cv.check_greater_than('universe ID', universe_id, 0, True)
self._id = universe_id
@name.setter
def name(self, name):
check_type('universe name', name, basestring)
cv.check_type('universe name', name, basestring)
self._name = name
def add_cell(self, cell):
@ -452,7 +452,7 @@ class Universe(object):
"""
if not isinstance(cell, Cell):
msg = 'Unable to add a Cell to Universe ID={0} since {1} is not ' \
msg = 'Unable to add a Cell to Universe ID="{0}" since "{1}" is not ' \
'a Cell'.format(self._id, cell)
raise ValueError(msg)
@ -472,7 +472,7 @@ class Universe(object):
"""
if not isinstance(cells, Iterable):
msg = 'Unable to add Cells to Universe ID={0} since {1} is not ' \
msg = 'Unable to add Cells to Universe ID="{0}" since "{1}" is not ' \
'iterable'.format(self._id, cells)
raise ValueError(msg)
@ -490,7 +490,7 @@ class Universe(object):
"""
if not isinstance(cell, Cell):
msg = 'Unable to remove a Cell from Universe ID={0} since {1} is ' \
msg = 'Unable to remove a Cell from Universe ID="{0}" since "{1}" is ' \
'not a Cell'.format(self._id, cell)
raise ValueError(msg)
@ -671,24 +671,25 @@ class Lattice(object):
self._id = AUTO_UNIVERSE_ID
AUTO_UNIVERSE_ID += 1
else:
check_type('lattice ID', lattice_id, Integral)
check_greater_than('lattice ID', lattice_id, 0)
cv.check_type('lattice ID', lattice_id, Integral)
cv.check_greater_than('lattice ID', lattice_id, 0)
self._id = lattice_id
@name.setter
def name(self, name):
check_type('lattice name', name, basestring)
cv.check_type('lattice name', name, basestring)
self._name = name
@outer.setter
def outer(self, outer):
check_type('outer universe', outer, Universe)
cv.check_type('outer universe', outer, Universe)
self._outer = outer
@universes.setter
def universes(self, universes):
check_type('lattice universes', universes, Iterable)
self._universes = np.asarray(universes, dtype=Universe)
cv.check_iterable_type('lattice universes', universes, Universe,
min_depth=2, max_depth=3)
self._universes = universes
def get_unique_universes(self):
"""Determine all unique universes in the lattice
@ -701,13 +702,22 @@ class Lattice(object):
"""
unique_universes = np.unique(self._universes.ravel())
universes = {}
univs = dict()
for k in range(len(self._universes)):
for j in range(len(self._universes[k])):
if isinstance(self._universes[k][j], Universe):
u = self._universes[k][j]
univs[u._id] = u
else:
for i in range(len(self._universes[k][j])):
u = self._universes[k][j][i]
assert isinstance(u, Universe)
univs[u._id] = u
for universe in unique_universes:
universes[universe._id] = universe
if self.outer is not None:
univs[self.outer._id] = self.outer
return universes
return univs
def get_all_nuclides(self):
"""Return all nuclides contained in the lattice
@ -825,29 +835,29 @@ class RectLattice(Lattice):
@dimension.setter
def dimension(self, dimension):
check_type('lattice dimension', dimension, Iterable, Integral)
check_length('lattice dimension', dimension, 2, 3)
cv.check_type('lattice dimension', dimension, Iterable, Integral)
cv.check_length('lattice dimension', dimension, 2, 3)
for dim in dimension:
check_greater_than('lattice dimension', dim, 0)
cv.check_greater_than('lattice dimension', dim, 0)
self._dimension = dimension
@lower_left.setter
def lower_left(self, lower_left):
check_type('lattice lower left corner', lower_left, Iterable, Real)
check_length('lattice lower left corner', lower_left, 2, 3)
cv.check_type('lattice lower left corner', lower_left, Iterable, Real)
cv.check_length('lattice lower left corner', lower_left, 2, 3)
self._lower_left = lower_left
@offsets.setter
def offsets(self, offsets):
check_type('lattice offsets', offsets, Iterable)
cv.check_type('lattice offsets', offsets, Iterable)
self._offsets = offsets
@Lattice.pitch.setter
def pitch(self, pitch):
check_type('lattice pitch', pitch, Iterable, Real)
check_length('lattice pitch', pitch, 2, 3)
cv.check_type('lattice pitch', pitch, Iterable, Real)
cv.check_length('lattice pitch', pitch, 2, 3)
for dim in pitch:
check_greater_than('lattice pitch', dim, 0.0)
cv.check_greater_than('lattice pitch', dim, 0.0)
self._pitch = pitch
def get_offset(self, path, filter_offset):
@ -1041,28 +1051,28 @@ class HexLattice(Lattice):
@num_rings.setter
def num_rings(self, num_rings):
check_type('number of rings', num_rings, Integral)
check_greater_than('number of rings', num_rings, 0)
cv.check_type('number of rings', num_rings, Integral)
cv.check_greater_than('number of rings', num_rings, 0)
self._num_rings = num_rings
@num_axial.setter
def num_axial(self, num_axial):
check_type('number of axial', num_axial, Integral)
check_greater_than('number of axial', num_axial, 0)
cv.check_type('number of axial', num_axial, Integral)
cv.check_greater_than('number of axial', num_axial, 0)
self._num_axial = num_axial
@center.setter
def center(self, center):
check_type('lattice center', center, Iterable, Real)
check_length('lattice center', center, 2, 3)
cv.check_type('lattice center', center, Iterable, Real)
cv.check_length('lattice center', center, 2, 3)
self._center = center
@Lattice.pitch.setter
def pitch(self, pitch):
check_type('lattice pitch', pitch, Iterable, Real)
check_length('lattice pitch', pitch, 1, 2)
cv.check_type('lattice pitch', pitch, Iterable, Real)
cv.check_length('lattice pitch', pitch, 1, 2)
for dim in pitch:
check_greater_than('lattice pitch', dim, 0)
cv.check_greater_than('lattice pitch', dim, 0)
self._pitch = pitch
@Lattice.universes.setter
@ -1091,14 +1101,14 @@ class HexLattice(Lattice):
# Set the number of axial positions.
if n_dims == 3:
self.num_axial = self._universes.shape[0]
self.num_axial = len(self._universes)
else:
self._num_axial = None
# Set the number of rings and make sure this number is consistent for
# all axial positions.
if n_dims == 3:
self.num_rings = len(self._universes[0])
self.num_rings = len(self._universes)
for rings in self._universes:
if len(rings) != self._num_rings:
msg = 'HexLattice ID={0:d} has an inconsistent number of ' \
@ -1106,7 +1116,7 @@ class HexLattice(Lattice):
raise ValueError(msg)
else:
self.num_rings = self._universes.shape[0]
self.num_rings = len(self._universes)
# Make sure there are the correct number of elements in each ring.
if n_dims == 3:

18
scripts/openmc-plot-mesh-tally
View file

@ -133,7 +133,11 @@ class MeshPlotter(tk.Frame):
self.mesh = selectedTally.filters_by_name['mesh'].mesh
# Get mesh dimensions
self.nx, self.ny, self.nz = self.mesh.dimension
if len(self.mesh.dimension) == 2:
self.nx, self.ny = self.mesh.dimension
self.nz = 1
else:
self.nx, self.ny, self.nz = self.mesh.dimension
# Repopulate comboboxes baesd on current basis selection
text = self.basisBox.get()
@ -210,7 +214,7 @@ class MeshPlotter(tk.Frame):
mesh_filter = f
continue
index = self.filterBoxes[f.type].current()
spec_list.append((f, index))
spec_list.append((f.type, (index,)))
text = self.basisBox.get()
if text == 'xy':
@ -229,11 +233,11 @@ class MeshPlotter(tk.Frame):
else:
meshtuple = (i + 1, axial_level, j + 1)
filters, filter_bins = zip(*spec_list + [
(mesh_filter, meshtuple)])
mean = selectedTally.get_value(
self.scoreBox.get(), filters, filter_bins)
stdev = selectedTally.get_value(
self.scoreBox.get(), filters, filter_bins,
(mesh_filter.type, (meshtuple,))])
mean = selectedTally.get_values(
[self.scoreBox.get()], filters, filter_bins)
stdev = selectedTally.get_values(
[self.scoreBox.get()], filters, filter_bins,
value='std_dev')
if mbvalue == 'Mean':
matrix[i, j] = mean

2
scripts/openmc-statepoint-3d
View file

@ -219,7 +219,7 @@ def main(file_, o):
for y in range(1,ny+1):
for z in range(1,nz+1):
filterspec[0][1] = (x,y,z)
val = sp.get_value(tally.id-1, filterspec, sid)[o.valerr]
val = sp.get_values(tally.id-1, filterspec, sid)[o.valerr]
if o.vtk:
# vtk cells go z, y, x, so we store it now and enter it later
i = (z-1)*nx*ny + (y-1)*nx + x-1

53
scripts/openmc-track-to-vtk
View file

@ -1,4 +1,4 @@
#!/usr/bin/env python2
#!/usr/bin/env python
"""Convert binary particle track to VTK poly data.
Usage information can be obtained by running 'track.py --help':
@ -64,26 +64,43 @@ def main():
for fname in args.input:
# Write coordinate values to points array.
if fname.endswith('.binary'):
track = open(fname, 'rb').read()
coords = [struct.unpack("ddd", track[24*i : 24*(i+1)])
for i in range(len(track)/24)]
n_points = len(coords)
for triplet in coords:
points.InsertNextPoint(triplet)
track = open(fname, 'rb')
# Determine number of particles and tracks/particle
n_particles = struct.unpack('i', track.read(4))[0]
n_coords = struct.unpack('i'*n_particles, track.read(4*n_particles))
coords = []
for i in range(n_particles):
# Read coordinates for each particle
coords.append([struct.unpack('ddd', track.read(24))
for j in range(n_coords[i])])
# Add coordinates to points data
for triplet in coords[i]:
points.InsertNextPoint(triplet)
else:
coords = h5py.File(fname).get('coordinates')
n_points = coords.shape[0]
for i in range(n_points):
points.InsertNextPoint(coords[i, :])
track = h5py.File(fname)
n_particles = track['n_particles'].value[0]
n_coords = track['n_coords']
coords = []
for i in range(n_particles):
coords.append(track['coordinates_' + str(i + 1)].value)
for j in range(n_coords[i]):
points.InsertNextPoint(coords[i][j,:])
# Create VTK line and assign points to line.
line = vtk.vtkPolyLine()
line.GetPointIds().SetNumberOfIds(n_points)
for i in range(n_points):
line.GetPointIds().SetId(i, point_offset+i)
for i in range(n_particles):
# Create VTK line and assign points to line.
line = vtk.vtkPolyLine()
line.GetPointIds().SetNumberOfIds(n_coords[i])
for j in range(n_coords[i]):
line.GetPointIds().SetId(j, point_offset + j)
# Add line to cell array
cells.InsertNextCell(line)
point_offset += n_coords[i]
cells.InsertNextCell(line)
point_offset += n_points
data = vtk.vtkPolyData()
data.SetPoints(points)
data.SetLines(cells)

2
setup.py
View file

@ -10,7 +10,7 @@ except ImportError:
have_setuptools = False
kwargs = {'name': 'openmc',
'version': '0.6.2',
'version': '0.7.0',
'packages': ['openmc'],
'scripts': glob.glob('scripts/openmc-*'),

27
src/ace.F90
View file

@ -88,16 +88,16 @@ contains
nuclides(i_nuclide) % resonant = .true.
nuclides(i_nuclide) % name_0K = nuclides_0K(n) % name_0K
nuclides(i_nuclide) % name_0K = trim(nuclides(i_nuclide) % &
& name_0K)
& name_0K)
nuclides(i_nuclide) % scheme = nuclides_0K(n) % scheme
nuclides(i_nuclide) % scheme = trim(nuclides(i_nuclide) % &
& scheme)
& scheme)
nuclides(i_nuclide) % E_min = nuclides_0K(n) % E_min
nuclides(i_nuclide) % E_max = nuclides_0K(n) % E_max
if (.not. already_read % contains(nuclides(i_nuclide) % &
& name_0K)) then
& name_0K)) then
i_listing = xs_listing_dict % get_key(nuclides(i_nuclide) % &
& name_0K)
& name_0K)
call read_ace_table(i_nuclide, i_listing)
end if
exit
@ -446,9 +446,10 @@ contains
if (nuc % elastic_0K(i) < ZERO) nuc % elastic_0K(i) = ZERO
! build xs cdf
xs_cdf_sum = xs_cdf_sum + (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) &
& + sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO &
& * (nuc % energy_0K(i+1) - nuc % energy_0K(i))
xs_cdf_sum = xs_cdf_sum &
+ (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) &
+ sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO &
* (nuc % energy_0K(i+1) - nuc % energy_0K(i))
nuc % xs_cdf(i) = xs_cdf_sum
end do
@ -752,31 +753,31 @@ contains
! Set flag and allocate space for Tab1 to store yield
rxn % multiplicity_with_E = .true.
allocate(rxn % multiplicity_E)
XSS_index = JXS(11) + rxn % multiplicity - 101
NR = nint(XSS(XSS_index))
rxn % multiplicity_E % n_regions = NR
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(rxn % multiplicity_E % nbt(NR))
allocate(rxn % multiplicity_E % int(NR))
end if
! read ENDF interpolation parameters
XSS_index = XSS_index + 1
if (NR > 0) then
rxn % multiplicity_E % nbt = get_int(NR)
rxn % multiplicity_E % int = get_int(NR)
end if
! allocate space for yield data
XSS_index = XSS_index + 2*NR
NE = nint(XSS(XSS_index))
rxn % multiplicity_E % n_pairs = NE
allocate(rxn % multiplicity_E % x(NE))
allocate(rxn % multiplicity_E % y(NE))
! read yield data
XSS_index = XSS_index + 1
rxn % multiplicity_E % x = get_real(NE)
@ -1480,7 +1481,7 @@ contains
table % inelastic_data(i) % e_out_pdf(j) = XSS(XSS_index + 2)
table % inelastic_data(i) % e_out_cdf(j) = XSS(XSS_index + 3)
table % inelastic_data(i) % mu(:, j) = &
XSS(XSS_index + 4: XSS_index + 4 + NMU - 1)
XSS(XSS_index + 4: XSS_index + 4 + NMU - 1)
XSS_index = XSS_index + 4 + NMU - 1
end do
end do

30
src/cmfd_data.F90
View file

@ -104,23 +104,23 @@ contains
cmfd % keff_bal = ZERO
! Begin loop around tallies
TAL: do ital = 1, n_cmfd_tallies
! Begin loop around tallies
TAL: do ital = 1, n_cmfd_tallies
! Associate tallies and mesh
t => cmfd_tallies(ital)
i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1)
m => meshes(i_mesh)
! Associate tallies and mesh
t => cmfd_tallies(ital)
i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1)
m => meshes(i_mesh)
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
i_filter_eout = t % find_filter(FILTER_ENERGYOUT)
i_filter_surf = t % find_filter(FILTER_SURFACE)
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
i_filter_eout = t % find_filter(FILTER_ENERGYOUT)
i_filter_surf = t % find_filter(FILTER_SURFACE)
! Begin loop around space
ZLOOP: do k = 1,nz
! Begin loop around space
ZLOOP: do k = 1,nz
YLOOP: do j = 1,ny
YLOOP: do j = 1,ny
XLOOP: do i = 1,nx
@ -600,7 +600,7 @@ contains
dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + &
cell_hxyz(xyz_idx)*neig_dc)
end if
end if
end if
@ -797,7 +797,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
(shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then
albedo = ONE
else
albedo = (current(2*l-1)/current(2*l))**(shift_idx)

2
src/cmfd_execute.F90
View file

@ -121,7 +121,7 @@ contains
! Allocate cmfd source if not already allocated and allocate buffer
if (.not. allocated(cmfd % cmfd_src)) &
allocate(cmfd % cmfd_src(ng,nx,ny,nz))
allocate(cmfd % cmfd_src(ng,nx,ny,nz))
! Reset cmfd source to 0
cmfd % cmfd_src = ZERO

10
src/cmfd_header.F90
View file

@ -1,4 +1,4 @@
module cmfd_header
module cmfd_header
use constants, only: CMFD_NOACCEL, ZERO, ONE
@ -17,7 +17,7 @@ module cmfd_header
! Core overlay map
integer, allocatable :: coremap(:,:,:)
integer, allocatable :: indexmap(:,:)
integer :: mat_dim = CMFD_NOACCEL
integer :: mat_dim = CMFD_NOACCEL
! Energy grid
real(8), allocatable :: egrid(:)
@ -51,7 +51,7 @@ module cmfd_header
! Source sites in each mesh box
real(8), allocatable :: sourcecounts(:,:,:,:)
! Weight adjustment factors
! Weight adjustment factors
real(8), allocatable :: weightfactors(:,:,:,:)
! Eigenvector/eigenvalue from cmfd run
@ -118,7 +118,7 @@ contains
if (.not. allocated(this % nfissxs)) allocate(this % nfissxs(ng,ng,nx,ny,nz))
if (.not. allocated(this % diffcof)) allocate(this % diffcof(ng,nx,ny,nz))
! Allocate dtilde and dhat
! Allocate dtilde and dhat
if (.not. allocated(this % dtilde)) allocate(this % dtilde(6,ng,nx,ny,nz))
if (.not. allocated(this % dhat)) allocate(this % dhat(6,ng,nx,ny,nz))
@ -167,7 +167,7 @@ contains
end subroutine allocate_cmfd
!===============================================================================
! DEALLOCATE_CMFD frees all memory of cmfd type
! DEALLOCATE_CMFD frees all memory of cmfd type
!===============================================================================
subroutine deallocate_cmfd(this)

18
src/cmfd_input.F90
View file

@ -120,7 +120,7 @@ contains
allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), &
cmfd % indices(3)))
if (get_arraysize_integer(node_mesh, "map") /= &
product(cmfd % indices(1:3))) then
product(cmfd % indices(1:3))) then
call fatal_error('CMFD coremap not to correct dimensions')
end if
allocate(iarray(get_arraysize_integer(node_mesh, "map")))
@ -156,7 +156,7 @@ contains
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.
dhat_reset = .true.
end if
! Set monitoring
@ -180,7 +180,7 @@ contains
call get_node_value(doc, "run_adjoint", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_run_adjoint = .true.
cmfd_run_adjoint = .true.
end if
! Batch to begin cmfd
@ -289,7 +289,7 @@ contains
! Determine number of dimensions for mesh
n = get_arraysize_integer(node_mesh, "dimension")
if (n /= 2 .and. n /= 3) then
call fatal_error("Mesh must be two or three dimensions.")
call fatal_error("Mesh must be two or three dimensions.")
end if
m % n_dimension = n
@ -318,14 +318,14 @@ contains
! Make sure both upper-right or width were specified
if (check_for_node(node_mesh, "upper_right") .and. &
check_for_node(node_mesh, "width")) then
check_for_node(node_mesh, "width")) then
call fatal_error("Cannot specify both <upper_right> and <width> on a &
&tally mesh.")
end if
! Make sure either upper-right or width was specified
if (.not.check_for_node(node_mesh, "upper_right") .and. &
.not.check_for_node(node_mesh, "width")) then
.not.check_for_node(node_mesh, "width")) then
call fatal_error("Must specify either <upper_right> and <width> on a &
&tally mesh.")
end if
@ -333,7 +333,7 @@ contains
if (check_for_node(node_mesh, "width")) then
! Check to ensure width has same dimensions
if (get_arraysize_double(node_mesh, "width") /= &
get_arraysize_double(node_mesh, "lower_left")) then
get_arraysize_double(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <width> must be the same as the &
&number of entries on <lower_left>.")
end if
@ -351,7 +351,7 @@ contains
elseif (check_for_node(node_mesh, "upper_right")) then
! Check to ensure width has same dimensions
if (get_arraysize_double(node_mesh, "upper_right") /= &
get_arraysize_double(node_mesh, "lower_left")) then
get_arraysize_double(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <upper_right> must be the same &
&as the number of entries on <lower_left>.")
end if
@ -548,7 +548,7 @@ contains
! Deallocate filter bins
deallocate(filters(1) % int_bins)
if (check_for_node(node_mesh, "energy")) &
deallocate(filters(2) % real_bins)
deallocate(filters(2) % real_bins)
end do

20
src/cmfd_loss_operator.F90
View file

@ -54,7 +54,7 @@ contains
n_e = 4*(nx - 2) + 4*(ny - 2) + 4*(nz - 2) ! define # of edges
n_s = 2*(nx - 2)*(ny - 2) + 2*(nx - 2)*(nz - 2) &
+ 2*(ny - 2)*(nz - 2) ! define # of sides
n_i = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors
n_i = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors
nz_c = ng*n_c*(4 + ng - 1) ! define # nonzero corners
nz_e = ng*n_e*(5 + ng - 1) ! define # nonzero edges
nz_s = ng*n_s*(6 + ng - 1) ! define # nonzero sides
@ -131,7 +131,7 @@ contains
! Check for global boundary
if (bound(l) /= nxyz(xyz_idx,dir_idx)) then
! Check for coremap
! Check for coremap
if (cmfd_coremap) then
! Check for neighbor that is non-acceleartred
@ -139,11 +139,11 @@ contains
CMFD_NOACCEL) then
! Get neighbor matrix index
call indices_to_matrix(g,neig_idx(1), neig_idx(2), &
call indices_to_matrix(g,neig_idx(1), neig_idx(2), &
neig_idx(3), neig_mat_idx, ng, nx, ny)
! Record nonzero
nnz = nnz + 1
nnz = nnz + 1
end if
@ -218,11 +218,11 @@ contains
real(8) :: jn ! direction dependent leakage coeff to neig
real(8) :: jo(6) ! leakage coeff in front of cell flux
real(8) :: jnet ! net leakage from jo
real(8) :: val ! temporary variable before saving to
real(8) :: val ! temporary variable before saving to
! Check for adjoint
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
if (present(adjoint)) adjoint_calc = adjoint
! Get maximum number of cells in each direction
nx = cmfd%indices(1)
@ -257,11 +257,11 @@ contains
dhat = ZERO
end if
! Create boundary vector
! Create boundary vector
bound = (/i,i,j,j,k,k/)
! Begin loop over leakages
! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z
! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z
LEAK: do l = 1,6
! Define (x,y,z) and (-,+) indices
@ -350,7 +350,7 @@ contains
scattxshg = cmfd%scattxs(h, g, i, j, k)
end if
! Negate the scattering xs
! Negate the scattering xs
val = -scattxshg
! Record value in matrix
@ -358,7 +358,7 @@ contains
end do SCATTR
end do ROWS
end do ROWS
! CSR requires n+1 row
call loss_matrix % new_row()

4
src/cmfd_prod_operator.F90
View file

@ -97,7 +97,7 @@ contains
end if
! loop around all other groups
! loop around all other groups
NFISS: do h = 1, ng
@ -122,7 +122,7 @@ contains
end do NFISS
end do ROWS
end do ROWS
! CSR requires n+1 row
call prod_matrix % new_row()

6
src/cmfd_solver.F90
View file

@ -65,7 +65,7 @@ contains
! Check for physical adjoint
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') &
physical_adjoint = .true.
physical_adjoint = .true.
! Start timer for build
call time_cmfdbuild % start()
@ -75,7 +75,7 @@ contains
! Check for mathematical adjoint calculation
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') &
call compute_adjoint()
call compute_adjoint()
! Stop timer for build
call time_cmfdbuild % stop()
@ -286,7 +286,7 @@ contains
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
jcol < prod % get_row(irow + 1)) then
loss % val(icol) = loss % val(icol) - ONE/k_s*prod % val(jcol)
jcol = jcol + 1
end if

7
src/constants.F90
View file

@ -7,8 +7,8 @@ module constants
! OpenMC major, minor, and release numbers
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 6
integer, parameter :: VERSION_RELEASE = 2
integer, parameter :: VERSION_MINOR = 7
integer, parameter :: VERSION_RELEASE = 0
! Revision numbers for binary files
integer, parameter :: REVISION_STATEPOINT = 13
@ -44,6 +44,9 @@ module constants
integer, parameter :: MAX_EVENTS = 10000
integer, parameter :: MAX_SAMPLE = 100000
! Maximum number of secondary particles created
integer, parameter :: MAX_SECONDARY = 1000
! Maximum number of words in a single line, length of line, and length of
! single word
integer, parameter :: MAX_WORDS = 500

10
src/cross_section.F90
View file

@ -446,7 +446,7 @@ contains
! Calculate elastic cross section/factor
elastic = ZERO
if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
i_up)))
@ -455,7 +455,7 @@ contains
! Calculate fission cross section/factor
fission = ZERO
if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
i_up)))
@ -464,7 +464,7 @@ contains
! Calculate capture cross section/factor
capture = ZERO
if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
i_up)))
@ -571,11 +571,11 @@ contains
! calculate interpolation factor
f = (E - nuc % energy_0K(i_grid)) &
& / (nuc % energy_0K(i_grid + 1) - 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)
& + f * nuc % elastic_0K(i_grid + 1)
end function elastic_xs_0K

22
src/dict_header.F90
View file

@ -56,7 +56,7 @@ module dict_header
!===============================================================================
! DICT* is a dictionary of (key,value) pairs with convenience methods as
! type-bound procedures. DictCharInt has character(*) keys and integer values,
! and DictIntInt has integer keys and values.
! and DictIntInt has integer keys and values.
!===============================================================================
type, public :: DictCharInt
@ -105,7 +105,7 @@ contains
if (associated(elem)) then
elem % value = value
else
! Get hash
! Get hash
hash = dict_hash_key_ci(key)
! Create new element
@ -135,7 +135,7 @@ contains
if (associated(elem)) then
elem % value = value
else
! Get hash
! Get hash
hash = dict_hash_key_ii(key)
! Create new element
@ -156,13 +156,13 @@ contains
!===============================================================================
function dict_get_elem_ci(this, key) result(elem)
class(DictCharInt) :: this
character(*), intent(in) :: key
type(ElemKeyValueCI), pointer :: elem
integer :: hash
! Check for dictionary not being allocated
if (.not. associated(this % table)) then
allocate(this % table(HASH_SIZE))
@ -178,13 +178,13 @@ contains
end function dict_get_elem_ci
function dict_get_elem_ii(this, key) result(elem)
class(DictIntInt) :: this
integer, intent(in) :: key
type(ElemKeyValueII), pointer :: elem
integer :: hash
! Check for dictionary not being allocated
if (.not. associated(this % table)) then
allocate(this % table(HASH_SIZE))
@ -279,7 +279,7 @@ contains
character(*), intent(in) :: key
integer :: val
integer :: i
val = 0
@ -298,7 +298,7 @@ contains
integer, intent(in) :: key
integer :: val
val = 0
! Added the absolute val on val-1 since the sum in the do loop is
@ -420,7 +420,7 @@ contains
integer :: i
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
if (associated(this % table)) then
do i = 1, size(this % table)
current => this % table(i) % list
@ -440,5 +440,5 @@ contains
end if
end subroutine dict_clear_ii
end module dict_header

205
src/doppler.F90
View file

@ -52,62 +52,115 @@ contains
! Loop over incoming neutron energies
ENERGY_NEUTRON: do i = 1, n
sigma = ZERO
y = x(i)
y_sq = y*y
y_inv = ONE / y
y_inv_sq = y_inv / y
sigma = ZERO
y = x(i)
y_sq = y*y
y_inv = ONE / y
y_inv_sq = y_inv / y
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4
k = i
a = ZERO
call calculate_F(F_a, a)
k = i
a = ZERO
call calculate_F(F_a, a)
do while (a >= -4.0 .and. k > 1)
! Move to next point
F_b = F_a
k = k - 1
a = x(k) - y
do while (a >= -4.0 .and. k > 1)
! Move to next point
F_b = F_a
k = k - 1
a = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k+1) - xs(k)) / (x(k+1)**2 - x(k)**2)
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k+1) - xs(k)) / (x(k+1)**2 - x(k)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! =======================================================================
! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE
! =======================================================================
! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE
if (k == 1 .and. a >= -4.0) then
! Since x = 0, this implies that a = -y
F_b = F_a
a = -y
if (k == 1 .and. a >= -4.0) then
! Since x = 0, this implies that a = -y
F_b = F_a
a = -y
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma + xs(k)*x(k)*(y_inv_sq*H(1) + y_inv*H(0))
end if
! Add contribution to broadened cross section
sigma = sigma + xs(k)*x(k)*(y_inv_sq*H(1) + y_inv*H(0))
end if
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4
k = i
b = ZERO
call calculate_F(F_b, b)
k = i
b = ZERO
call calculate_F(F_b, b)
do while (b <= 4.0 .and. k < n)
do while (b <= 4.0 .and. k < n)
! Move to next point
F_a = F_b
k = k + 1
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_b, b)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! =======================================================================
! EXTEND CROSS SECTION TO INFINITY ASSUMING CONSTANT SHAPE
if (k == n .and. b <= 4.0) then
! Calculate F function at last energy point
a = x(k) - y
call calculate_F(F_a, a)
! Add contribution to broadened cross section
sigma = sigma + xs(k) * (y_inv_sq*F_a(2) + 2.0*y_inv*F_a(1) + F_a(0))
end if
! =======================================================================
! EVALUATE SECOND TERM FROM x(k) + y = 0 to +4
if (y <= 4.0) then
! Swap signs on y
y = -y
y_inv = -y_inv
k = 1
! Calculate a and b based on 0 and x(1)
a = -y
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
call calculate_F(F_b, b)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma - xs(k) * x(k) * (y_inv_sq*H(1) + y_inv*H(0))
! Now progress forward doing the remainder of the second term
do while (b <= 4.0)
! Move to next point
F_a = F_b
k = k + 1
@ -119,69 +172,17 @@ contains
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) &
+ y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
sigma = sigma - Ak*(xs(k) - slope*x(k)**2) - slope*Bk
end do
end if
! =======================================================================
! EXTEND CROSS SECTION TO INFINITY ASSUMING CONSTANT SHAPE
if (k == n .and. b <= 4.0) then
! Calculate F function at last energy point
a = x(k) - y
call calculate_F(F_a, a)
! Add contribution to broadened cross section
sigma = sigma + xs(k) * (y_inv_sq*F_a(2) + 2.0*y_inv*F_a(1) + F_a(0))
end if
! =======================================================================
! EVALUATE SECOND TERM FROM x(k) + y = 0 to +4
if (y <= 4.0) then
! Swap signs on y
y = -y
y_inv = -y_inv
k = 1
! Calculate a and b based on 0 and x(1)
a = -y
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
call calculate_F(F_b, b)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma - xs(k) * x(k) * (y_inv_sq*H(1) + y_inv*H(0))
! Now progress forward doing the remainder of the second term
do while (b <= 4.0)
! Move to next point
F_a = F_b
k = k + 1
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_b, b)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma - Ak*(xs(k) - slope*x(k)**2) - slope*Bk
end do
end if
! Set broadened cross section
sigmaNew(i) = sigma
! Set broadened cross section
sigmaNew(i) = sigma
end do ENERGY_NEUTRON

298
src/eigenvalue.F90
View file

@ -4,273 +4,24 @@ module eigenvalue
use message_passing
#endif
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use error, only: fatal_error, warning
use global
use math, only: t_percentile
use mesh, only: count_bank_sites
use mesh_header, only: StructuredMesh
use output, only: write_message, header, print_columns, &
print_batch_keff, print_generation
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use source, only: get_source_particle
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use trigger, only: check_triggers
use tracking, only: transport
implicit none
private
public :: run_eigenvalue
real(8) :: keff_generation ! Single-generation k on each
! processor
real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
real(8) :: keff_generation ! Single-generation k on each processor
real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
contains
!===============================================================================
! RUN_EIGENVALUE encompasses all the main logic where iterations are performed
! over the batches, generations, and histories in a k-eigenvalue calculation.
!===============================================================================
subroutine run_eigenvalue()
type(Particle) :: p
integer(8) :: i_work
if (master) call header("K EIGENVALUE SIMULATION", level=1)
! Display column titles
if(master) call print_columns()
! Turn on inactive timer
call time_inactive % start()
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_max_batches
call initialize_batch()
! Handle restart runs
if (restart_run .and. current_batch <= restart_batch) then
call replay_batch_history()
cycle BATCH_LOOP
end if
! =======================================================================
! LOOP OVER GENERATIONS
GENERATION_LOOP: do current_gen = 1, gen_per_batch
call initialize_generation()
! Start timer for transport
call time_transport % start()
! ====================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
! grab source particle from bank
call get_source_particle(p, current_work)
! transport particle
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
call finalize_generation()
end do GENERATION_LOOP
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
end do BATCH_LOOP
call time_active % stop()
! ==========================================================================
! END OF RUN WRAPUP
if (master) call header("SIMULATION FINISHED", level=1)
! Clear particle
call p % clear()
end subroutine run_eigenvalue
!===============================================================================
! INITIALIZE_BATCH
!===============================================================================
subroutine initialize_batch()
call write_message("Simulating batch " // trim(to_str(current_batch)) &
&// "...", 8)
! Reset total starting particle weight used for normalizing tallies
total_weight = ZERO
if (current_batch == n_inactive + 1) then
! Switch from inactive batch timer to active batch timer
call time_inactive % stop()
call time_active % start()
! Enable active batches (and tallies_on if it hasn't been enabled)
active_batches = .true.
tallies_on = .true.
! Add user tallies to active tallies list
!$omp parallel
call setup_active_usertallies()
!$omp end parallel
end if
! check CMFD initialize batch
if (cmfd_run) call cmfd_init_batch()
end subroutine initialize_batch
!===============================================================================
! INITIALIZE_GENERATION
!===============================================================================
subroutine initialize_generation()
! set overall generation number
overall_gen = gen_per_batch*(current_batch - 1) + current_gen
! Reset number of fission bank sites
n_bank = 0
! Count source sites if using uniform fission source weighting
if (ufs) call count_source_for_ufs()
! Store current value of tracklength k
keff_generation = global_tallies(K_TRACKLENGTH) % value
end subroutine initialize_generation
!===============================================================================
! FINALIZE_GENERATION
!===============================================================================
subroutine finalize_generation()
! Update global tallies with the omp private accumulation variables
!$omp parallel
!$omp critical
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + global_tally_collision
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + global_tally_leakage
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + global_tally_absorption
!$omp end critical
! reset private tallies
global_tally_tracklength = 0
global_tally_collision = 0
global_tally_leakage = 0
global_tally_absorption = 0
!$omp end parallel
#ifdef _OPENMP
! Join the fission bank from each thread into one global fission bank
call join_bank_from_threads()
#endif
! Distribute fission bank across processors evenly
call time_bank % start()
call synchronize_bank()
call time_bank % stop()
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Collect results and statistics
call calculate_generation_keff()
call calculate_average_keff()
! Write generation output
if (master .and. current_gen /= gen_per_batch) call print_generation()
end subroutine finalize_generation
!===============================================================================
! FINALIZE_BATCH handles synchronization and accumulation of tallies,
! calculation of Shannon entropy, getting single-batch estimate of keff, and
! turning on tallies when appropriate
!===============================================================================
subroutine finalize_batch()
! Collect tallies
call time_tallies % start()
call synchronize_tallies()
call time_tallies % stop()
! Reset global tally results
if (.not. active_batches) then
call reset_result(global_tallies)
n_realizations = 0
end if
! Perform CMFD calculation if on
if (cmfd_on) call execute_cmfd()
! Display output
if (master) call print_batch_keff()
! Calculate combined estimate of k-effective
if (master) call calculate_combined_keff()
! Check_triggers
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call write_state_point()
end if
! Write out source point if it's been specified for this batch
if ((sourcepoint_batch % contains(current_batch) .or. source_latest) .and. &
source_write) then
call write_source_point()
end if
if (master .and. current_batch == n_max_batches) then
! Make sure combined estimate of k-effective is calculated at the last
! batch in case no state point is written
call calculate_combined_keff()
end if
end subroutine finalize_batch
!===============================================================================
! SYNCHRONIZE_BANK samples source sites from the fission sites that were
! accumulated during the generation. This routine is what allows this Monte
@ -830,40 +581,11 @@ contains
end subroutine count_source_for_ufs
!===============================================================================
! REPLAY_BATCH_HISTORY displays keff and entropy for each generation within a
! batch using data read from a state point file
!===============================================================================
subroutine replay_batch_history
! Write message at beginning
if (current_batch == 1) then
call write_message("Replaying history from state point...", 1)
end if
do current_gen = 1, gen_per_batch
overall_gen = overall_gen + 1
call calculate_average_keff()
! print out batch keff
if (current_gen < gen_per_batch) then
if (master) call print_generation()
else
if (master) call print_batch_keff()
end if
end do
! Write message at end
if (current_batch == restart_batch) then
call write_message("Resuming simulation...", 1)
end if
end subroutine replay_batch_history
#ifdef _OPENMP
!===============================================================================
! JOIN_BANK_FROM_THREADS
! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission
! bank that can be sampled. Note that this operation is necessarily sequential
! to preserve the order of the bank when using varying numbers of threads.
!===============================================================================
subroutine join_bank_from_threads()
@ -880,8 +602,8 @@ contains
!$omp do ordered schedule(static)
do i = 1, n_threads
!$omp ordered
master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank)
total = total + n_bank
master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank)
total = total + n_bank
!$omp end ordered
end do
!$omp end do
@ -891,10 +613,10 @@ contains
! Now copy the shared fission bank sites back to the master thread's copy.
if (thread_id == 0) then
n_bank = total
fission_bank(1:n_bank) = master_fission_bank(1:n_bank)
n_bank = total
fission_bank(1:n_bank) = master_fission_bank(1:n_bank)
else
n_bank = 0
n_bank = 0
end if
!$omp end parallel

16
src/endf_header.F90
View file

@ -3,7 +3,7 @@ module endf_header
implicit none
!===============================================================================
! TAB1 represents a one-dimensional interpolable function
! TAB1 represents a one-dimensional interpolable function
!===============================================================================
type Tab1
@ -13,28 +13,28 @@ module endf_header
integer :: n_pairs ! # of pairs of (x,y) values
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
! Type-Bound procedures
contains
procedure :: clear => tab1_clear ! deallocates a Tab1 Object.
end type Tab1
contains
!===============================================================================
! TAB1_CLEAR deallocates the items in Tab1
!===============================================================================
subroutine tab1_clear(this)
class(Tab1), intent(inout) :: this ! The Tab1 to clear
if (allocated(this % nbt)) &
deallocate(this % nbt, this % int)
if (allocated(this % x)) &
deallocate(this % x, this % y)
end subroutine tab1_clear
end module endf_header

165
src/fixed_source.F90
View file

@ -1,165 +0,0 @@
module fixed_source
#ifdef MPI
use message_passing
#endif
use constants, only: ZERO, MAX_LINE_LEN
use global
use output, only: write_message, header
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use source, only: sample_external_source, copy_source_attributes
use state_point, only: write_state_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies
use trigger, only: check_triggers
use tracking, only: transport
implicit none
contains
subroutine run_fixedsource()
integer(8) :: i ! index over histories in single cycle
type(Particle) :: p
if (master) call header("FIXED SOURCE TRANSPORT SIMULATION", level=1)
! Allocate particle and dummy source site
!$omp parallel
allocate(source_site)
!$omp end parallel
! Turn timer and tallies on
tallies_on = .true.
!$omp parallel
call setup_active_usertallies()
!$omp end parallel
call time_active % start()
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_max_batches
! In a restart run, skip any batches that have already been simulated
if (restart_run .and. current_batch <= restart_batch) then
if (current_batch > n_inactive) n_realizations = n_realizations + 1
cycle BATCH_LOOP
end if
call initialize_batch()
! Start timer for transport
call time_transport % start()
! =======================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p)
PARTICLE_LOOP: do i = 1, work
! Set unique particle ID
p % id = (current_batch - 1)*n_particles + work_index(rank) + i
! set particle trace
trace = .false.
if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
work_index(rank) + i == trace_particle) trace = .true.
! set random number seed
call set_particle_seed(p % id)
! grab source particle from bank
call sample_source_particle(p)
! transport particle
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
end do BATCH_LOOP
call time_active % stop()
! ==========================================================================
! END OF RUN WRAPUP
if (master) call header("SIMULATION FINISHED", level=1)
end subroutine run_fixedsource
!===============================================================================
! INITIALIZE_BATCH
!===============================================================================
subroutine initialize_batch()
call write_message("Simulating batch " // trim(to_str(current_batch)) &
&// "...", 1)
! Reset total starting particle weight used for normalizing tallies
total_weight = ZERO
end subroutine initialize_batch
!===============================================================================
! FINALIZE_BATCH
!===============================================================================
subroutine finalize_batch()
! Collect and accumulate tallies
call time_tallies % start()
call synchronize_tallies()
call time_tallies % stop()
! Check_triggers
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call write_state_point()
end if
end subroutine finalize_batch
!===============================================================================
! SAMPLE_SOURCE_PARTICLE
!===============================================================================
subroutine sample_source_particle(p)
type(Particle), intent(inout) :: p
! Set particle
call p % initialize()
! Sample the external source distribution
call sample_external_source(source_site)
! Copy source attributes to the particle
call copy_source_attributes(p, source_site)
! Determine whether to create track file
if (write_all_tracks) p % write_track = .true.
end subroutine sample_source_particle
end module fixed_source

69
src/geometry_header.F90
View file

@ -9,13 +9,13 @@ module geometry_header
!===============================================================================
type Universe
integer :: id ! Unique ID
integer :: type ! Type
integer :: n_cells ! # of cells within
integer, allocatable :: cells(:) ! List of cells within
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
integer :: id ! Unique ID
integer :: type ! Type
integer :: n_cells ! # of cells within
integer, allocatable :: cells(:) ! List of cells within
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
end type Universe
!===============================================================================
@ -119,14 +119,14 @@ module geometry_header
!===============================================================================
type Surface
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of surface
real(8), allocatable :: coeffs(:) ! Definition of surface
integer, allocatable :: &
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of surface
real(8), allocatable :: coeffs(:) ! Definition of surface
integer, allocatable :: &
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
end type Surface
!===============================================================================
@ -134,24 +134,29 @@ module geometry_header
!===============================================================================
type Cell
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of cell (normal, universe, lattice)
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
integer :: instances ! number of instances of this cell in the geom
integer :: material ! Material within cell (0 for universe)
integer :: n_surfaces ! Number of surfaces within
integer, allocatable :: offset (:) ! Distribcell offset for tally counter
integer, allocatable :: &
& surfaces(:) ! List of surfaces bounding cell -- note that
! parentheses, union, etc operators will be listed
! here too
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of cell (normal, universe,
! lattice)
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
integer :: instances ! number of instances of this cell in
! the geom
integer :: material ! Material within cell (0 for
! universe)
integer :: n_surfaces ! Number of surfaces within
integer, allocatable :: offset (:) ! Distribcell offset for tally
! counter
integer, allocatable :: &
& surfaces(:) ! List of surfaces bounding cell
! -- note that parentheses, union,
! etc operators will be listed here
! too
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)
real(8), allocatable :: rotation(:)
real(8), allocatable :: rotation_matrix(:,:)
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)
real(8), allocatable :: rotation(:)
real(8), allocatable :: rotation_matrix(:,:)
end type Cell
! array index of universe 0

10
src/global.F90
View file

@ -283,10 +283,6 @@ module global
! Mode to run in (fixed source, eigenvalue, plotting, etc)
integer :: run_mode = NONE
! Fixed source particle bank
type(Bank), pointer :: source_site => null()
!$omp threadprivate(source_site)
! Restart run
logical :: restart_run = .false.
integer :: restart_batch
@ -529,9 +525,9 @@ contains
! Deallocate entropy mesh
if (associated(entropy_mesh)) then
if (allocated(entropy_mesh % lower_left)) &
deallocate(entropy_mesh % lower_left)
deallocate(entropy_mesh % lower_left)
if (allocated(entropy_mesh % upper_right)) &
deallocate(entropy_mesh % upper_right)
deallocate(entropy_mesh % upper_right)
if (allocated(entropy_mesh % width)) deallocate(entropy_mesh % width)
deallocate(entropy_mesh)
end if
@ -541,7 +537,7 @@ contains
if (associated(ufs_mesh)) then
if (allocated(ufs_mesh % lower_left)) deallocate(ufs_mesh % lower_left)
if (allocated(ufs_mesh % upper_right)) &
deallocate(ufs_mesh % upper_right)
deallocate(ufs_mesh % upper_right)
if (allocated(ufs_mesh % width)) deallocate(ufs_mesh % width)
deallocate(ufs_mesh)
end if

2
src/hdf5_interface.F90
View file

@ -7,7 +7,7 @@ module hdf5_interface
use, intrinsic :: ISO_C_BINDING
#ifdef MPI
use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL
use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL
#endif
implicit none

28
src/hdf5_summary.F90
View file

@ -186,7 +186,7 @@ contains
call su % write_data("lattice", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call su % write_data(lattices(c % fill) % obj % id, "lattice", &
group="geometry/cells/cell " // trim(to_str(c % id)))
group="geometry/cells/cell " // trim(to_str(c % id)))
end select
! Write list of bounding surfaces
@ -373,7 +373,7 @@ contains
! Write lattice universes.
allocate(lattice_universes(lat % n_cells(1), lat % n_cells(2), &
&lat % n_cells(3)))
&lat % n_cells(3)))
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
@ -394,23 +394,23 @@ contains
! Write number of lattice cells.
call su % write_data(lat % n_rings, "n_rings", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call su % write_data(lat % n_rings, "n_axial", &
call su % write_data(lat % n_axial, "n_axial", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
! Write lattice center, pitch and outer universe.
if (lat % is_3d) then
call su % write_data(lat % center, "center", length=3, &
call su % write_data(lat % center, "center", length=3, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
else
call su % write_data(lat % center, "center", length=2, &
call su % write_data(lat % center, "center", length=2, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
if (lat % is_3d) then
call su % write_data(lat % pitch, "pitch", length=3, &
call su % write_data(lat % pitch, "pitch", length=2, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
else
call su % write_data(lat % pitch, "pitch", length=2, &
call su % write_data(lat % pitch, "pitch", length=1, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
@ -429,7 +429,7 @@ contains
! Write lattice universes.
allocate(lattice_universes(2*lat % n_rings - 1, 2*lat % n_rings - 1, &
&lat % n_axial))
&lat % n_axial))
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
@ -447,7 +447,7 @@ contains
end do
end do
call su % write_data(lattice_universes, "universes", &
&length=(/lat % n_axial, 2*lat % n_rings-1, 2*lat % n_rings-1/), &
&length=(/2*lat % n_rings-1, 2*lat % n_rings-1, lat % n_axial/), &
&group="geometry/lattices/lattice " // trim(to_str(lat % id)))
deallocate(lattice_universes)
end select
@ -647,12 +647,12 @@ contains
// "/filter " // trim(to_str(j)))
case(FILTER_ENERGYIN)
call su % write_data("energy", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_ENERGYOUT)
call su % write_data("energyout", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
end select
end do FILTER_LOOP

45
src/initialize.F90
View file

@ -20,7 +20,6 @@ module initialize
write_message
use output_interface
use random_lcg, only: initialize_prng
use source, only: initialize_source
use state_point, only: load_state_point
use string, only: to_str, str_to_int, starts_with, ends_with
use tally_header, only: TallyObject, TallyResult, TallyFilter
@ -141,13 +140,9 @@ contains
! Determine how much work each processor should do
call calculate_work()
! Allocate banks and create source particles -- for a fixed source
! calculation, the external source distribution is sampled during the
! run, not at initialization
if (run_mode == MODE_EIGENVALUE) then
call allocate_banks()
if (.not. restart_run) call initialize_source()
end if
! Allocate source bank, and for eigenvalue simulations also allocate the
! fission bank
call allocate_banks()
! If this is a restart run, load the state point data and binary source
! file
@ -904,31 +899,33 @@ contains
call fatal_error("Failed to allocate source bank.")
end if
if (run_mode == MODE_EIGENVALUE) then
#ifdef _OPENMP
! If OpenMP is being used, each thread needs its own private fission
! bank. Since the private fission banks need to be combined at the end of a
! generation, there is also a 'master_fission_bank' that is used to collect
! the sites from each thread.
! If OpenMP is being used, each thread needs its own private fission
! bank. Since the private fission banks need to be combined at the end of
! a generation, there is also a 'master_fission_bank' that is used to
! collect the sites from each thread.
n_threads = omp_get_max_threads()
n_threads = omp_get_max_threads()
!$omp parallel
thread_id = omp_get_thread_num()
thread_id = omp_get_thread_num()
if (thread_id == 0) then
allocate(fission_bank(3*work))
else
allocate(fission_bank(3*work/n_threads))
end if
if (thread_id == 0) then
allocate(fission_bank(3*work))
else
allocate(fission_bank(3*work/n_threads))
end if
!$omp end parallel
allocate(master_fission_bank(3*work), STAT=alloc_err)
allocate(master_fission_bank(3*work), STAT=alloc_err)
#else
allocate(fission_bank(3*work), STAT=alloc_err)
allocate(fission_bank(3*work), STAT=alloc_err)
#endif
! Check for allocation errors
if (alloc_err /= 0) then
call fatal_error("Failed to allocate fission bank.")
! Check for allocation errors
if (alloc_err /= 0) then
call fatal_error("Failed to allocate fission bank.")
end if
end if
end subroutine allocate_banks

238
src/input_xml.F90
View file

@ -472,7 +472,7 @@ contains
! Check for type of energy distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('monoenergetic')
external_source % type_energy = SRC_ENERGY_MONO
@ -798,7 +798,7 @@ contains
if (.not. source_separate) then
do i = 1, n_source_points
if (.not. statepoint_batch % contains(sourcepoint_batch % &
get_item(i))) then
get_item(i))) then
call fatal_error('Sourcepoint batches are not a subset&
& of statepoint batches.')
end if
@ -811,7 +811,7 @@ contains
call get_node_value(doc, "no_reduce", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
reduce_tallies = .false.
reduce_tallies = .false.
end if
! Check if the user has specified to use confidence intervals for
@ -884,11 +884,11 @@ contains
&// trim(to_str(i)) // " in settings.xml file!")
end if
call get_node_value(node_scatterer, "nuclide", &
nuclides_0K(i) % nuclide)
nuclides_0K(i) % nuclide)
if (check_for_node(node_scatterer, "method")) then
call get_node_value(node_scatterer, "method", &
nuclides_0K(i) % scheme)
nuclides_0K(i) % scheme)
end if
! check to make sure xs name for which method is applied is given
@ -898,7 +898,7 @@ contains
&// " given in cross_sections.xml")
end if
call get_node_value(node_scatterer, "xs_label", &
nuclides_0K(i) % name)
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
@ -906,11 +906,11 @@ contains
&// trim(to_str(i)) // " given in cross_sections.xml")
end if
call get_node_value(node_scatterer, "xs_label_0K", &
nuclides_0K(i) % name_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)
nuclides_0K(i) % E_min)
end if
! check that E_min is non-negative
@ -921,7 +921,7 @@ contains
if (check_for_node(node_scatterer, "E_max")) then
call get_node_value(node_scatterer, "E_max", &
nuclides_0K(i) % E_max)
nuclides_0K(i) % E_max)
end if
! check that E_max is not less than E_min
@ -1081,7 +1081,7 @@ contains
! Read material
word = ''
if (check_for_node(node_cell, "material")) &
call get_node_value(node_cell, "material", word)
call get_node_value(node_cell, "material", word)
select case(to_lower(word))
case ('void')
c % material = MATERIAL_VOID
@ -1248,7 +1248,7 @@ contains
! Copy and interpret surface type
word = ''
if (check_for_node(node_surf, "type")) &
call get_node_value(node_surf, "type", word)
call get_node_value(node_surf, "type", word)
select case(to_lower(word))
case ('x-plane')
s % type = SURF_PX
@ -1306,7 +1306,7 @@ contains
! Boundary conditions
word = ''
if (check_for_node(node_surf, "boundary")) &
call get_node_value(node_surf, "boundary", word)
call get_node_value(node_surf, "boundary", word)
select case (to_lower(word))
case ('transmission', 'transmit', '')
s % bc = BC_TRANSMIT
@ -1447,7 +1447,7 @@ contains
do k = 0, n_y - 1
do j = 1, n_x
lat % universes(j, n_y - k, m) = &
&temp_int_array(j + n_x*k + n_x*n_y*(m-1))
&temp_int_array(j + n_x*k + n_x*n_y*(m-1))
end do
end do
end do
@ -1714,7 +1714,7 @@ contains
! Copy default cross section if present
if (check_for_node(doc, "default_xs")) &
call get_node_value(doc, "default_xs", default_xs)
call get_node_value(doc, "default_xs", default_xs)
! Get pointer to list of XML <material>
call get_node_list(doc, "material", node_mat_list)
@ -1860,7 +1860,7 @@ contains
! store full name
call get_node_value(node_nuc, "name", temp_str)
if (check_for_node(node_nuc, "xs")) &
call get_node_value(node_nuc, "xs", name)
call get_node_value(node_nuc, "xs", name)
name = trim(temp_str) // "." // trim(name)
! save name and density to list
@ -1869,7 +1869,7 @@ contains
! Check if no atom/weight percents were specified or if both atom and
! weight percents were specified
if (.not.check_for_node(node_nuc, "ao") .and. &
.not.check_for_node(node_nuc, "wo")) then
.not.check_for_node(node_nuc, "wo")) then
call fatal_error("No atom or weight percent specified for nuclide " &
&// trim(name))
elseif (check_for_node(node_nuc, "ao") .and. &
@ -1919,7 +1919,7 @@ contains
! Check if no atom/weight percents were specified or if both atom and
! weight percents were specified
if (.not.check_for_node(node_ele, "ao") .and. &
.not.check_for_node(node_ele, "wo")) then
.not.check_for_node(node_ele, "wo")) then
call fatal_error("No atom or weight percent specified for element " &
&// trim(name))
elseif (check_for_node(node_ele, "ao") .and. &
@ -2036,7 +2036,7 @@ contains
! Determine name of S(a,b) table
if (.not.check_for_node(node_sab, "name") .or. &
.not.check_for_node(node_sab, "xs")) then
.not.check_for_node(node_sab, "xs")) then
call fatal_error("Need to specify <name> and <xs> for S(a,b) &
&table.")
end if
@ -2183,7 +2183,7 @@ contains
call get_node_value(doc, "assume_separate", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
assume_separate = .true.
assume_separate = .true.
end if
! ==========================================================================
@ -2211,7 +2211,7 @@ contains
! Read mesh type
temp_str = ''
if (check_for_node(node_mesh, "type")) &
call get_node_value(node_mesh, "type", temp_str)
call get_node_value(node_mesh, "type", temp_str)
select case (to_lower(temp_str))
case ('rect', 'rectangle', 'rectangular')
m % type = LATTICE_RECT
@ -2253,14 +2253,14 @@ contains
! Make sure both upper-right or width were specified
if (check_for_node(node_mesh, "upper_right") .and. &
check_for_node(node_mesh, "width")) then
check_for_node(node_mesh, "width")) then
call fatal_error("Cannot specify both <upper_right> and <width> on a &
&tally mesh.")
end if
! Make sure either upper-right or width was specified
if (.not.check_for_node(node_mesh, "upper_right") .and. &
.not.check_for_node(node_mesh, "width")) then
.not.check_for_node(node_mesh, "width")) then
call fatal_error("Must specify either <upper_right> and <width> on a &
&tally mesh.")
end if
@ -2268,7 +2268,7 @@ contains
if (check_for_node(node_mesh, "width")) then
! Check to ensure width has same dimensions
if (get_arraysize_double(node_mesh, "width") /= &
get_arraysize_double(node_mesh, "lower_left")) then
get_arraysize_double(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <width> must be the same as &
&the number of entries on <lower_left>.")
end if
@ -2286,7 +2286,7 @@ contains
elseif (check_for_node(node_mesh, "upper_right")) then
! Check to ensure width has same dimensions
if (get_arraysize_double(node_mesh, "upper_right") /= &
get_arraysize_double(node_mesh, "lower_left")) then
get_arraysize_double(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <upper_right> must be the &
&same as the number of entries on <lower_left>.")
end if
@ -2349,7 +2349,7 @@ contains
! Copy tally name
if (check_for_node(node_tal, "name")) &
call get_node_value(node_tal, "name", t % name)
call get_node_value(node_tal, "name", t % name)
! =======================================================================
! READ DATA FOR FILTERS
@ -2371,13 +2371,13 @@ contains
! Convert filter type to lower case
temp_str = ''
if (check_for_node(node_filt, "type")) &
call get_node_value(node_filt, "type", temp_str)
call get_node_value(node_filt, "type", temp_str)
temp_str = to_lower(temp_str)
! Determine number of bins
if (check_for_node(node_filt, "bins")) then
if (trim(temp_str) == 'energy' .or. &
trim(temp_str) == 'energyout') then
trim(temp_str) == 'energyout') then
n_words = get_arraysize_double(node_filt, "bins")
else
n_words = get_arraysize_integer(node_filt, "bins")
@ -2644,7 +2644,7 @@ contains
if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
@ -2687,7 +2687,7 @@ contains
if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
@ -2710,7 +2710,7 @@ contains
if (starts_with(score_name,trim(MOMENT_N_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
@ -3032,7 +3032,7 @@ contains
! Get the trigger type - "variance", "std_dev" or "rel_err"
if (check_for_node(node_trigger, "type")) then
call get_node_value(node_trigger, "type", temp_str)
temp_str = to_lower(temp_str)
temp_str = to_lower(temp_str)
else
call fatal_error("Must specify trigger type for tally " // &
trim(to_str(t % id)) // " in tally XML file.")
@ -3131,7 +3131,7 @@ contains
! Increment the overall trigger index
trig_ind = trig_ind + 1
end if
end if
end do SCORE_LOOP
! Deallocate the list of tally scores used to create triggers
@ -3249,7 +3249,7 @@ contains
! Copy plot type
temp_str = 'slice'
if (check_for_node(node_plot, "type")) &
call get_node_value(node_plot, "type", temp_str)
call get_node_value(node_plot, "type", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("slice")
@ -3264,7 +3264,7 @@ contains
! Set output file path
filename = trim(to_str(pl % id)) // "_plot"
if (check_for_node(node_plot, "filename")) &
call get_node_value(node_plot, "filename", filename)
call get_node_value(node_plot, "filename", filename)
select case (pl % type)
case (PLOT_TYPE_SLICE)
pl % path_plot = trim(path_input) // trim(filename) // ".ppm"
@ -3309,7 +3309,7 @@ contains
if (pl % type == PLOT_TYPE_SLICE) then
temp_str = 'xy'
if (check_for_node(node_plot, "basis")) &
call get_node_value(node_plot, "basis", temp_str)
call get_node_value(node_plot, "basis", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("xy")
@ -3364,7 +3364,7 @@ contains
! Copy plot color type and initialize all colors randomly
temp_str = "cell"
if (check_for_node(node_plot, "color")) &
call get_node_value(node_plot, "color", temp_str)
call get_node_value(node_plot, "color", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("cell")
@ -3478,7 +3478,7 @@ contains
! Ensure that there is a linewidth for this meshlines specification
if (check_for_node(node_meshlines, "linewidth")) then
call get_node_value(node_meshlines, "linewidth", &
pl % meshlines_width)
pl % meshlines_width)
else
call fatal_error("Must specify a linewidth for meshlines &
&specification in plot " // trim(to_str(pl % id)))
@ -3494,7 +3494,7 @@ contains
end if
call get_node_array(node_meshlines, "color", &
pl % meshlines_color % rgb)
pl % meshlines_color % rgb)
else
pl % meshlines_color % rgb = (/ 0, 0, 0 /)
@ -3520,8 +3520,8 @@ contains
end if
i_mesh = cmfd_tallies(1) % &
filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % &
int_bins(1)
filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % &
int_bins(1)
pl % meshlines_mesh => meshes(i_mesh)
case ('entropy')
@ -3680,9 +3680,9 @@ contains
! Check if cross_sections.xml exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find cross_sections.xml file
call fatal_error("Cross sections XML file '" &
&// trim(path_cross_sections) // "' does not exist!")
! Could not find cross_sections.xml file
call fatal_error("Cross sections XML file '" &
&// trim(path_cross_sections) // "' does not exist!")
end if
call write_message("Reading cross sections XML file...", 5)
@ -3691,28 +3691,28 @@ contains
call open_xmldoc(doc, path_cross_sections)
if (check_for_node(doc, "directory")) then
! Copy directory information if present
call get_node_value(doc, "directory", directory)
! Copy directory information if present
call get_node_value(doc, "directory", directory)
else
! If no directory is listed in cross_sections.xml, by default select the
! directory in which the cross_sections.xml file resides
i = index(path_cross_sections, "/", BACK=.true.)
directory = path_cross_sections(1:i)
! If no directory is listed in cross_sections.xml, by default select the
! directory in which the cross_sections.xml file resides
i = index(path_cross_sections, "/", BACK=.true.)
directory = path_cross_sections(1:i)
end if
! determine whether binary/ascii
temp_str = ''
if (check_for_node(doc, "filetype")) &
call get_node_value(doc, "filetype", temp_str)
call get_node_value(doc, "filetype", temp_str)
if (trim(temp_str) == 'ascii') then
filetype = ASCII
filetype = ASCII
elseif (trim(temp_str) == 'binary') then
filetype = BINARY
filetype = BINARY
elseif (len_trim(temp_str) == 0) then
filetype = ASCII
filetype = ASCII
else
call fatal_error("Unknown filetype in cross_sections.xml: " &
&// trim(temp_str))
call fatal_error("Unknown filetype in cross_sections.xml: " &
&// trim(temp_str))
end if
! copy default record length and entries for binary files
@ -3727,83 +3727,83 @@ contains
! Allocate xs_listings array
if (n_listings == 0) then
call fatal_error("No ACE table listings present in cross_sections.xml &
&file!")
call fatal_error("No ACE table listings present in cross_sections.xml &
&file!")
else
allocate(xs_listings(n_listings))
allocate(xs_listings(n_listings))
end if
do i = 1, n_listings
listing => xs_listings(i)
listing => xs_listings(i)
! Get pointer to ace table XML node
call get_list_item(node_ace_list, i, node_ace)
! Get pointer to ace table XML node
call get_list_item(node_ace_list, i, node_ace)
! copy a number of attributes
call get_node_value(node_ace, "name", listing % name)
if (check_for_node(node_ace, "alias")) &
call get_node_value(node_ace, "alias", listing % alias)
call get_node_value(node_ace, "zaid", listing % zaid)
call get_node_value(node_ace, "awr", listing % awr)
if (check_for_node(node_ace, "temperature")) &
call get_node_value(node_ace, "temperature", listing % kT)
call get_node_value(node_ace, "location", listing % location)
! copy a number of attributes
call get_node_value(node_ace, "name", listing % name)
if (check_for_node(node_ace, "alias")) &
call get_node_value(node_ace, "alias", listing % alias)
call get_node_value(node_ace, "zaid", listing % zaid)
call get_node_value(node_ace, "awr", listing % awr)
if (check_for_node(node_ace, "temperature")) &
call get_node_value(node_ace, "temperature", listing % kT)
call get_node_value(node_ace, "location", listing % location)
! determine type of cross section
if (ends_with(listing % name, 'c')) then
listing % type = ACE_NEUTRON
elseif (ends_with(listing % name, 't')) then
listing % type = ACE_THERMAL
end if
! determine type of cross section
if (ends_with(listing % name, 'c')) then
listing % type = ACE_NEUTRON
elseif (ends_with(listing % name, 't')) then
listing % type = ACE_THERMAL
end if
! set filetype, record length, and number of entries
if (check_for_node(node_ace, "filetype")) then
temp_str = ''
call get_node_value(node_ace, "filetype", temp_str)
if (temp_str == 'ascii') then
listing % filetype = ASCII
else if (temp_str == 'binary') then
listing % filetype = BINARY
end if
else
listing % filetype = filetype
end if
! set filetype, record length, and number of entries
if (check_for_node(node_ace, "filetype")) then
temp_str = ''
call get_node_value(node_ace, "filetype", temp_str)
if (temp_str == 'ascii') then
listing % filetype = ASCII
else if (temp_str == 'binary') then
listing % filetype = BINARY
end if
else
listing % filetype = filetype
end if
! Set record length and entries for binary files
if (filetype == BINARY) then
listing % recl = recl
listing % entries = entries
end if
! Set record length and entries for binary files
if (filetype == BINARY) then
listing % recl = recl
listing % entries = entries
end if
! determine metastable state
if (.not.check_for_node(node_ace, "metastable")) then
listing % metastable = .false.
else
listing % metastable = .true.
end if
! determine metastable state
if (.not.check_for_node(node_ace, "metastable")) then
listing % metastable = .false.
else
listing % metastable = .true.
end if
! determine path of cross section table
if (check_for_node(node_ace, "path")) then
call get_node_value(node_ace, "path", temp_str)
else
call fatal_error("Path missing for isotope " // listing % name)
end if
! determine path of cross section table
if (check_for_node(node_ace, "path")) then
call get_node_value(node_ace, "path", temp_str)
else
call fatal_error("Path missing for isotope " // listing % name)
end if
if (starts_with(temp_str, '/')) then
listing % path = trim(temp_str)
else
if (ends_with(directory,'/')) then
listing % path = trim(directory) // trim(temp_str)
else
listing % path = trim(directory) // '/' // trim(temp_str)
end if
end if
if (starts_with(temp_str, '/')) then
listing % path = trim(temp_str)
else
if (ends_with(directory,'/')) then
listing % path = trim(directory) // trim(temp_str)
else
listing % path = trim(directory) // '/' // trim(temp_str)
end if
end if
! create dictionary entry for both name and alias
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(to_lower(listing % alias), i)
end if
! create dictionary entry for both name and alias
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(to_lower(listing % alias), i)
end if
end do
! Check that 0K nuclides are listed in the cross_sections.xml file

18
src/list_header.F90
View file

@ -46,7 +46,7 @@ module list_header
private
integer :: count = 0 ! Number of elements in list
! Used in get_item for fast sequential lookups
! Used in get_item for fast sequential lookups
integer :: last_index = huge(0)
type(ListElemInt), pointer :: last_elem => null()
@ -68,7 +68,7 @@ module list_header
private
integer :: count = 0 ! Number of elements in list
! Used in get_item for fast sequential lookups
! Used in get_item for fast sequential lookups
integer :: last_index = huge(0)
type(ListElemReal), pointer :: last_elem => null()
@ -90,7 +90,7 @@ module list_header
private
integer :: count = 0 ! Number of elements in list
! Used in get_item for fast sequential lookups
! Used in get_item for fast sequential lookups
integer :: last_index = huge(0)
type(ListElemChar), pointer :: last_elem => null()
@ -199,7 +199,7 @@ contains
type(ListElemInt), pointer :: current => null()
type(ListElemInt), pointer :: next => null()
if (this % count > 0) then
current => this % head
do while (associated(current))
@ -225,7 +225,7 @@ contains
type(ListElemReal), pointer :: current => null()
type(ListElemReal), pointer :: next => null()
if (this % count > 0) then
current => this % head
do while (associated(current))
@ -251,7 +251,7 @@ contains
type(ListElemChar), pointer :: current => null()
type(ListElemChar), pointer :: next => null()
if (this % count > 0) then
current => this % head
do while (associated(current))
@ -519,7 +519,7 @@ contains
! Allocate new element
allocate(new_elem)
new_elem % data = data
! Put it before the i-th element
new_elem % prev => elem % prev
new_elem % next => elem
@ -574,7 +574,7 @@ contains
! Allocate new element
allocate(new_elem)
new_elem % data = data
! Put it before the i-th element
new_elem % prev => elem % prev
new_elem % next => elem
@ -629,7 +629,7 @@ contains
! Allocate new element
allocate(new_elem)
new_elem % data = data
! Put it before the i-th element
new_elem % prev => elem % prev
new_elem % next => elem

9
src/main.F90
View file

@ -1,13 +1,12 @@
program main
use constants
use eigenvalue, only: run_eigenvalue
use finalize, only: finalize_run
use fixed_source, only: run_fixedsource
use global
use initialize, only: initialize_run
use particle_restart, only: run_particle_restart
use plot, only: run_plot
use simulation, only: run_simulation
implicit none
@ -16,10 +15,8 @@ program main
! start problem based on mode
select case (run_mode)
case (MODE_FIXEDSOURCE)
call run_fixedsource()
case (MODE_EIGENVALUE)
call run_eigenvalue()
case (MODE_FIXEDSOURCE, MODE_EIGENVALUE)
call run_simulation()
case (MODE_PLOTTING)
call run_plot()
case (MODE_PARTICLE)

241
src/math.F90
View file

@ -143,20 +143,20 @@ contains
pnx = 7.875_8 * (x ** 5) - 8.75_8 * x * x * x + 1.875 * x
case(6)
pnx = 14.4375_8 * (x ** 6) - 19.6875_8 * (x ** 4) + &
6.5625_8 * x * x - 0.3125_8
6.5625_8 * x * x - 0.3125_8
case(7)
pnx = 26.8125_8 * (x ** 7) - 43.3125_8 * (x ** 5) + &
19.6875_8 * x * x * x - 2.1875_8 * x
19.6875_8 * x * x * x - 2.1875_8 * x
case(8)
pnx = 50.2734375_8 * (x ** 8) - 93.84375_8 * (x ** 6) + &
54.140625 * (x ** 4) - 9.84375_8 * x * x + 0.2734375_8
54.140625 * (x ** 4) - 9.84375_8 * x * x + 0.2734375_8
case(9)
pnx = 94.9609375_8 * (x ** 9) - 201.09375_8 * (x ** 7) + &
140.765625_8 * (x ** 5) - 36.09375_8 * x * x * x + 2.4609375_8 * x
140.765625_8 * (x ** 5) - 36.09375_8 * x * x * x + 2.4609375_8 * x
case(10)
pnx = 180.42578125_8 * (x ** 10) - 427.32421875_8 * (x ** 8) + &
351.9140625_8 * (x ** 6) - 117.3046875_8 * (x ** 4) + &
13.53515625_8 * x * x - 0.24609375_8
351.9140625_8 * (x ** 6) - 117.3046875_8 * (x ** 4) + &
13.53515625_8 * x * x - 0.24609375_8
case default
pnx = ONE ! correct for case(0), incorrect for the rest
end select
@ -215,12 +215,12 @@ contains
rn(2) = 1.93649167310371_8 * w*(w2m1) * sin(TWO*phi)
! l = 3, m = -1
rn(3) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * &
sin(phi)
sin(phi)
! l = 3, m = 0
rn(4) = 2.5_8 * w**3 - 1.5_8 * w
! l = 3, m = 1
rn(5) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * &
cos(phi)
cos(phi)
! l = 3, m = 2
rn(6) = 1.93649167310371_8 * w*(w2m1) * cos(TWO*phi)
! l = 3, m = 3
@ -232,18 +232,18 @@ contains
rn(2) = 2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * sin(THREE* phi)
! l = 4, m = -2
rn(3) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * &
sin(TWO*phi)
sin(TWO*phi)
! l = 4, m = -1
rn(4) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w) * &
sin(phi)
rn(4) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)&
* sin(phi)
! l = 4, m = 0
rn(5) = 4.375_8 * w**4 - 3.75_8 * w**2 + 0.375_8
! l = 4, m = 1
rn(6) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w) * &
cos(phi)
rn(6) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)&
* cos(phi)
! l = 4, m = 2
rn(7) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * &
cos(TWO*phi)
cos(TWO*phi)
! l = 4, m = 3
rn(8) = 2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * cos(THREE* phi)
! l = 4, m = 4
@ -255,24 +255,26 @@ contains
rn(2) = 2.21852991866236_8 * w*(w2m1)**2 * sin(4.0_8*phi)
! l = 5, m = -3
rn(3) = 0.00996023841111995_8 * (w2m1)**(THREE/TWO)* &
((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi)
((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi)
! l = 5, m = -2
rn(4) = 0.0487950036474267_8 * (w2m1)*((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * &
sin(TWO*phi)
rn(4) = 0.0487950036474267_8 * (w2m1) &
* ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * sin(TWO*phi)
! l = 5, m = -1
rn(5) = 0.258198889747161_8*sqrt(w2m1)* &
((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) * sin(phi)
((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) &
* sin(phi)
! l = 5, m = 0
rn(6) = 7.875_8 * w**5 - 8.75_8 * w**3 + 1.875_8 * w
! l = 5, m = 1
rn(7) = 0.258198889747161_8*sqrt(w2m1)* &
((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) * cos(phi)
((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) &
* cos(phi)
! l = 5, m = 2
rn(8) = 0.0487950036474267_8 * (w2m1)* &
((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi)
((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi)
! l = 5, m = 3
rn(9) = 0.00996023841111995_8 * (w2m1)**(THREE/TWO)* &
((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi)
((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi)
! l = 5, m = 4
rn(10) = 2.21852991866236_8 * w*(w2m1)**2 * cos(4.0_8*phi)
! l = 5, m = 5
@ -284,30 +286,34 @@ contains
rn(2) = 2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * sin(5.0_8*phi)
! l = 6, m = -4
rn(3) = 0.00104990131391452_8 * (w2m1)**2 * &
((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi)
((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi)
! l = 6, m = -3
rn(4) = 0.00575054632785295_8 * (w2m1)**(THREE/TWO) * &
((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi)
((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi)
! l = 6, m = -2
rn(5) = 0.0345032779671177_8 * (w2m1) * &
((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * sin(TWO*phi)
((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) &
* sin(TWO*phi)
! l = 6, m = -1
rn(6) = 0.218217890235992_8*sqrt(w2m1) * &
((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) * sin(phi)
((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) &
* sin(phi)
! l = 6, m = 0
rn(7) = 14.4375_8 * w**6 - 19.6875_8 * w**4 + 6.5625_8 * w**2 - 0.3125_8
! l = 6, m = 1
rn(8) = 0.218217890235992_8*sqrt(w2m1) * &
((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) * cos(phi)
((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) &
* cos(phi)
! l = 6, m = 2
rn(9) = 0.0345032779671177_8 * (w2m1) * &
((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * cos(TWO*phi)
((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) &
* cos(TWO*phi)
! l = 6, m = 3
rn(10) = 0.00575054632785295_8 * (w2m1)**(THREE/TWO) * &
((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi)
((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi)
! l = 6, m = 4
rn(11) = 0.00104990131391452_8 * (w2m1)**2 * &
((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi)
((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi)
! l = 6, m = 5
rn(12) = 2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * cos(5.0_8*phi)
! l = 6, m = 6
@ -319,42 +325,43 @@ contains
rn(2) = 2.42182459624969_8 * w*(w2m1)**3 * sin(6.0_8*phi)
! l = 7, m = -5
rn(3) = 9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* &
((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi)
((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi)
! l = 7, m = -4
rn(4) = 0.000548293079133141_8 * (w2m1)**2* &
((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi)
((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi)
! l = 7, m = -3
rn(5) = 0.00363696483726654_8 * (w2m1)**(THREE/TWO)* &
((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
sin(THREE*phi)
((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
sin(THREE*phi)
! l = 7, m = -2
rn(6) = 0.025717224993682_8 * (w2m1)* &
((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
sin(TWO*phi)
((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
sin(TWO*phi)
! l = 7, m = -1
rn(7) = 0.188982236504614_8*sqrt(w2m1)* &
((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
(945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi)
((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
(945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi)
! l = 7, m = 0
rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 * w
rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 &
* w
! l = 7, m = 1
rn(9) = 0.188982236504614_8*sqrt(w2m1)* &
((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
(945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi)
((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
(945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi)
! l = 7, m = 2
rn(10) = 0.025717224993682_8 * (w2m1)* &
((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
cos(TWO*phi)
((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
cos(TWO*phi)
! l = 7, m = 3
rn(11) = 0.00363696483726654_8 * (w2m1)**(THREE/TWO)* &
((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
cos(THREE*phi)
((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
cos(THREE*phi)
! l = 7, m = 4
rn(12) = 0.000548293079133141_8 * (w2m1)**2 * &
((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi)
((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi)
! l = 7, m = 5
rn(13) = 9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* &
((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi)
((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi)
! l = 7, m = 6
rn(14) = 2.42182459624969_8 * w*(w2m1)**3 * cos(6.0_8*phi)
! l = 7, m = 7
@ -366,48 +373,50 @@ contains
rn(2) = 2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * sin(7.0_8*phi)
! l = 8, m = -6
rn(3) = 6.77369783729086d-6*(w2m1)**3* &
((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi)
((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi)
! l = 8, m = -5
rn(4) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)* &
((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi)
((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi)
! l = 8, m = -4
rn(5) = 0.000316557156832328_8 * (w2m1)**2* &
((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * sin(4.0_8*phi)
((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 &
+ 10395.0_8/8.0_8) * sin(4.0_8*phi)
! l = 8, m = -3
rn(6) = 0.00245204119306875_8 * (w2m1)**(THREE/TWO)* &
((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + (10395.0_8/8.0_8)*w) * sin(THREE*phi)
((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 &
+ (10395.0_8/8.0_8)*w) * sin(THREE*phi)
! l = 8, m = -2
rn(7) = 0.0199204768222399_8 * (w2m1)* &
((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + &
(10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi)
((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + &
(10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi)
! l = 8, m = -1
rn(8) = 0.166666666666667_8*sqrt(w2m1)* &
((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
(3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi)
((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
(3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi)
! l = 8, m = 0
rn(9) = 50.2734375_8 * w**8 - 93.84375_8 * w**6 + 54.140625_8 * w**4 -&
9.84375_8 * w**2 + 0.2734375_8
9.84375_8 * w**2 + 0.2734375_8
! l = 8, m = 1
rn(10) = 0.166666666666667_8*sqrt(w2m1)* &
((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
(3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi)
((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
(3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi)
! l = 8, m = 2
rn(11) = 0.0199204768222399_8 * (w2m1)*((45045.0_8/16.0_8)*w**6- &
45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - &
315.0_8/16.0_8) * cos(TWO*phi)
45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - &
315.0_8/16.0_8) * cos(TWO*phi)
! l = 8, m = 3
rn(12) = 0.00245204119306875_8 * (w2m1)**(THREE/TWO)* &
((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + &
(10395.0_8/8.0_8)*w) * cos(THREE*phi)
((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + &
(10395.0_8/8.0_8)*w) * cos(THREE*phi)
! l = 8, m = 4
rn(13) = 0.000316557156832328_8 * (w2m1)**2*((675675.0_8/8.0_8)*w**4 - &
135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi)
135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi)
! l = 8, m = 5
rn(14) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 - &
135135.0_8/TWO*w) * cos(5.0_8*phi)
rn(14) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 -&
135135.0_8/TWO*w) * cos(5.0_8*phi)
! l = 8, m = 6
rn(15) = 6.77369783729086d-6*(w2m1)**3*((2027025.0_8/TWO)*w**2 - &
135135.0_8/TWO) * cos(6.0_8*phi)
135135.0_8/TWO) * cos(6.0_8*phi)
! l = 8, m = 7
rn(16) = 2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * cos(7.0_8*phi)
! l = 8, m = 8
@ -419,54 +428,56 @@ contains
rn(2) = 2.58397773170915_8 * w*(w2m1)**4 * sin(8.0_8*phi)
! l = 9, m = -7
rn(3) = 4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* &
((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi)
((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi)
! l = 9, m = -6
rn(4) = 3.02928976464514d-6*(w2m1)**3* &
((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi)
((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi)
! l = 9, m = -5
rn(5) = 2.34647776186144d-5*(w2m1)**(5.0_8/TWO)* &
((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + &
135135.0_8/8.0_8) * sin(5.0_8*phi)
((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + &
135135.0_8/8.0_8) * sin(5.0_8*phi)
! l = 9, m = -4
rn(6) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - &
675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi)
675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi)
! l = 9, m = -3
rn(7) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)* &
((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + &
(135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi)
((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + &
(135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi)
! l = 9, m = -2
rn(8) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7- &
135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/16.0_8 * w)* &
sin(TWO*phi)
135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 &
- 3465.0_8/16.0_8 * w) * sin(TWO*phi)
! l = 9, m = -1
rn(9) = 0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - &
45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 * w**2 + 315.0_8/128.0_8) * sin(phi)
45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 &
* w**2 + 315.0_8/128.0_8) * sin(phi)
! l = 9, m = 0
rn(10) = 94.9609375_8 * w**9 - 201.09375_8 * w**7 + 140.765625_8 * w**5- &
36.09375_8 * w**3 + 2.4609375_8 * w
36.09375_8 * w**3 + 2.4609375_8 * w
! l = 9, m = 1
rn(11) = 0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - &
45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 * w**2 + 315.0_8/128.0_8) * cos(phi)
45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 &
* w**2 + 315.0_8/128.0_8) * cos(phi)
! l = 9, m = 2
rn(12) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7 - &
135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/ 16.0_8 * w) * &
cos(TWO*phi)
135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 &
- 3465.0_8/ 16.0_8 * w) * cos(TWO*phi)
! l = 9, m = 3
rn(13) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)*w**6 - &
675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8)* &
cos(THREE*phi)
rn(13) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)&
*w**6 - 675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 &
- 3465.0_8/16.0_8)* cos(THREE*phi)
! l = 9, m = 4
rn(14) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - &
675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi)
675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi)
! l = 9, m = 5
rn(15) = 2.34647776186144d-5*(w2m1)**(5.0_8/TWO)*((11486475.0_8/8.0_8)* &
w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi)
w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi)
! l = 9, m = 6
rn(16) = 3.02928976464514d-6*(w2m1)**3*((11486475.0_8/TWO)*w**3 - &
2027025.0_8/TWO*w) * cos(6.0_8*phi)
2027025.0_8/TWO*w) * cos(6.0_8*phi)
! l = 9, m = 7
rn(17) = 4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* &
((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi)
((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi)
! l = 9, m = 8
rn(18) = 2.58397773170915_8 * w*(w2m1)**4 * cos(8.0_8*phi)
! l = 9, m = 9
@ -478,65 +489,65 @@ contains
rn(2) = 2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * sin(9.0_8*phi)
! l = 10, m = -8
rn(3) = 2.49953651452314d-8*(w2m1)**4*((654729075.0_8/TWO)*w**2 - &
34459425.0_8/TWO) * sin(8.0_8*phi)
34459425.0_8/TWO) * sin(8.0_8*phi)
! l = 10, m = -7
rn(4) = 1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* &
((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi)
((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi)
! l = 10, m = -6
rn(5) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - &
34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi)
34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi)
! l = 10, m = -5
rn(6) = 1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* &
((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
(2027025.0_8/8.0_8)*w) * sin(5.0_8*phi)
((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
(2027025.0_8/8.0_8)*w) * sin(5.0_8*phi)
! l = 10, m = -4
rn(7) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - &
11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
45045.0_8/16.0_8) * sin(4.0_8*phi)
11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
45045.0_8/16.0_8) * sin(4.0_8*phi)
! l = 10, m = -3
rn(8) = 0.00127230170115096_8 * (w2m1)**(THREE/TWO)* &
((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
(675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi)
((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
(675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi)
! l = 10, m = -2
rn(9) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - &
765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - &
45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi)
765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - &
45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi)
! l = 10, m = -1
rn(10) = 0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - &
109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - &
15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi)
109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - &
15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi)
! l = 10, m = 0
rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 * w**6 - &
117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8
rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 &
* w**6 - 117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8
! l = 10, m = 1
rn(12) = 0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - &
109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ &
32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi)
109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ &
32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi)
! l = 10, m = 2
rn(13) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - &
765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -&
45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi)
765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -&
45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi)
! l = 10, m = 3
rn(14) = 0.00127230170115096_8 * (w2m1)**(THREE/TWO)* &
((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
(675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi)
((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
(675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi)
! l = 10, m = 4
rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - &
11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
45045.0_8/16.0_8) * cos(4.0_8*phi)
rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 -&
11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
45045.0_8/16.0_8) * cos(4.0_8*phi)
! l = 10, m = 5
rn(16) = 1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* &
((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
(2027025.0_8/8.0_8)*w) * cos(5.0_8*phi)
((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
(2027025.0_8/8.0_8)*w) * cos(5.0_8*phi)
! l = 10, m = 6
rn(17) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - &
34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi)
34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi)
! l = 10, m = 7
rn(18) = 1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* &
((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi)
((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi)
! l = 10, m = 8
rn(19) = 2.49953651452314d-8*(w2m1)**4* &
((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi)
((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi)
! l = 10, m = 9
rn(20) = 2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * cos(9.0_8*phi)
! l = 10, m = 10

26
src/matrix_header.F90
View file

@ -11,19 +11,19 @@ module matrix_header
integer, allocatable :: row(:) ! csr row vector
integer, allocatable :: col(:) ! column vector
real(8), allocatable :: val(:) ! matrix value vector
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 :: vector_multiply => matrix_vector_multiply
procedure :: search_indices => matrix_search_indices
procedure :: write => matrix_write
procedure :: copy => matrix_copy
procedure :: transpose => matrix_transpose
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 :: vector_multiply => matrix_vector_multiply
procedure :: search_indices => matrix_search_indices
procedure :: write => matrix_write
procedure :: copy => matrix_copy
procedure :: transpose => matrix_transpose
end type matrix
contains

28
src/mpiio_interface.F90
View file

@ -129,13 +129,13 @@ contains
integer, intent(inout) :: buffer ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_integer
@ -341,13 +341,13 @@ contains
real(8), intent(inout) :: buffer ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_double

65
src/output.F90
View file

@ -7,7 +7,7 @@ module output
use endf, only: reaction_name
use error, only: fatal_error, warning
use geometry_header, only: Cell, Universe, Surface, Lattice, RectLattice, &
&HexLattice, BASE_UNIVERSE
HexLattice, BASE_UNIVERSE
use global
use math, only: t_percentile
use mesh_header, only: StructuredMesh
@ -1372,7 +1372,7 @@ contains
if (cmfd_run) then
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
if (cmfd_display /= '') &
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
end if
write(UNIT=ou, FMT=*)
@ -1432,20 +1432,20 @@ contains
! write out cmfd keff if it is active and other display info
if (cmfd_on) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % k_cmfd(current_batch)
cmfd % k_cmfd(current_batch)
select case(trim(cmfd_display))
case('entropy')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % entropy(current_batch)
cmfd % entropy(current_batch)
case('balance')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % balance(current_batch)
cmfd % balance(current_batch)
case('source')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % src_cmp(current_batch)
cmfd % src_cmp(current_batch)
case('dominance')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % dom(current_batch)
cmfd % dom(current_batch)
end select
end if
@ -1547,11 +1547,15 @@ contains
write(ou,100) "Total time in simulation", time_inactive % elapsed + &
time_active % elapsed
write(ou,100) " Time in transport only", time_transport % elapsed
write(ou,100) " Time in inactive batches", time_inactive % elapsed
if (run_mode == MODE_EIGENVALUE) then
write(ou,100) " Time in inactive batches", time_inactive % elapsed
end if
write(ou,100) " Time in active batches", time_active % elapsed
write(ou,100) " Time synchronizing fission bank", time_bank % elapsed
write(ou,100) " Sampling source sites", time_bank_sample % elapsed
write(ou,100) " SEND/RECV source sites", time_bank_sendrecv % elapsed
if (run_mode == MODE_EIGENVALUE) then
write(ou,100) " Time synchronizing fission bank", time_bank % elapsed
write(ou,100) " Sampling source sites", time_bank_sample % elapsed
write(ou,100) " SEND/RECV source sites", time_bank_sendrecv % elapsed
end if
write(ou,100) " Time accumulating tallies", time_tallies % elapsed
if (cmfd_run) write(ou,100) " Time in CMFD", time_cmfd % elapsed
if (cmfd_run) write(ou,100) " Building matrices", &
@ -1567,7 +1571,7 @@ contains
speed_inactive = real(n_particles * (n_inactive - restart_batch) * &
gen_per_batch) / time_inactive % elapsed
speed_active = real(n_particles * n_active * gen_per_batch) / &
time_active % elapsed
time_active % elapsed
else
speed_inactive = ZERO
speed_active = real(n_particles * (n_batches - restart_batch) * &
@ -1884,21 +1888,22 @@ contains
select case(t % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
score_name = 'P' // trim(to_str(t % moment_order(k))) // " " // &
score_names(abs(t % score_bins(k)))
score_names(abs(t % score_bins(k)))
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
do n_order = 0, t % moment_order(k)
score_index = score_index + 1
score_name = 'P' // trim(to_str(n_order)) // " " //&
score_names(abs(t % score_bins(k)))
score_names(abs(t % score_bins(k)))
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) &
% sum_sq))
end do
k = k + t % moment_order(k)
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
@ -1908,11 +1913,13 @@ contains
do nm_order = -n_order, n_order
score_index = score_index + 1
score_name = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order)) // " " // score_names(abs(t % score_bins(k)))
trim(to_str(nm_order)) // " " &
// score_names(abs(t % score_bins(k)))
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index)&
% sum_sq))
end do
end do
k = k + (t % moment_order(k) + 1)**2 - 1
@ -1923,9 +1930,9 @@ contains
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), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
end select
end do
indent = indent - 2
@ -2334,8 +2341,8 @@ contains
! Loop over lattice coordinates
do k = 1, n_x
do l = 1, n_y
do m = 1, n_z
do l = 1, n_y
do m = 1, n_z
if (final >= lat % offset(map, k, l, m) + offset) then
if (k == n_x .and. l == n_y .and. m == n_z) then

28
src/output_interface.F90
View file

@ -31,7 +31,7 @@ module output_interface
#endif
#endif
logical :: serial ! Serial I/O when using MPI/PHDF5
contains
contains
generic, public :: write_data => write_double, &
write_double_1Darray, &
write_double_2Darray, &
@ -111,7 +111,7 @@ contains
self % serial = serial
else
self % serial = .true.
end if
end if
#ifdef HDF5
# ifdef MPI
@ -153,7 +153,7 @@ contains
self % serial = serial
else
self % serial = .true.
end if
end if
#ifdef HDF5
# ifdef MPI
@ -201,18 +201,18 @@ contains
#ifdef HDF5
# ifdef MPI
call hdf5_file_close(self % hdf5_fh)
call hdf5_file_close(self % hdf5_fh)
# else
call hdf5_file_close(self % hdf5_fh)
call hdf5_file_close(self % hdf5_fh)
# endif
#elif MPI
if (self % serial) then
close(UNIT=self % unit_fh)
else
call mpi_close_file(self % mpi_fh)
end if
if (self % serial) then
close(UNIT=self % unit_fh)
else
call mpi_close_file(self % mpi_fh)
end if
#else
close(UNIT=self % unit_fh)
close(UNIT=self % unit_fh)
#endif
end subroutine file_close
@ -475,7 +475,7 @@ contains
call hdf5_read_double_1Darray(self % hdf5_grp, name_, buffer, length)
else
call hdf5_read_double_1Darray_parallel(self % hdf5_grp, name_, buffer, &
length, collect_)
length, collect_)
end if
# else
call hdf5_read_double_1Darray(self % hdf5_grp, name_, buffer, length)
@ -663,8 +663,8 @@ contains
if (self % serial) then
call hdf5_write_double_3Darray(self % hdf5_grp, name_, buffer, length)
else
call hdf5_write_double_3Darray_parallel(self % hdf5_grp, name_, buffer, length, &
collect_)
call hdf5_write_double_3Darray_parallel(self % hdf5_grp, name_, buffer, &
length, collect_)
end if
# else
call hdf5_write_double_3Darray(self % hdf5_grp, name_, buffer, length)

66
src/particle_header.F90
View file

@ -1,6 +1,8 @@
module particle_header
use constants, only: NEUTRON, ONE, NONE, ZERO
use bank_header, only: Bank
use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY
use error, only: fatal_error
use geometry_header, only: BASE_UNIVERSE
implicit none
@ -79,9 +81,15 @@ module particle_header
! Track output
logical :: write_track = .false.
! Secondary particles created
integer :: n_secondary = 0
type(Bank) :: secondary_bank(MAX_SECONDARY)
contains
procedure :: initialize => initialize_particle
procedure :: clear => clear_particle
procedure :: initialize_from_source
procedure :: create_secondary
end type Particle
contains
@ -122,7 +130,7 @@ contains
end subroutine initialize_particle
!===============================================================================
! CLEAR_PARTICLE
! CLEAR_PARTICLE resets all coordinate levels for the particle
!===============================================================================
subroutine clear_particle(this)
@ -138,7 +146,7 @@ contains
end subroutine clear_particle
!===============================================================================
! RESET_COORD
! RESET_COORD clears data from a single coordinate level
!===============================================================================
elemental subroutine reset_coord(this)
@ -154,4 +162,56 @@ contains
end subroutine reset_coord
!===============================================================================
! INITIALIZE_FROM_SOURCE initializes a particle from data stored in a source
! site. The source site may have been produced from an external source, from
! fission, or simply as a secondary particle.
!===============================================================================
subroutine initialize_from_source(this, src)
class(Particle), intent(inout) :: this
type(Bank), intent(in) :: src
! set defaults
call this % initialize()
! copy attributes from source bank site
this % wgt = src % wgt
this % last_wgt = src % wgt
this % coord(1) % xyz = src % xyz
this % coord(1) % uvw = src % uvw
this % last_xyz = src % xyz
this % last_uvw = src % uvw
this % E = src % E
this % last_E = src % E
end subroutine initialize_from_source
!===============================================================================
! CREATE_SECONDARY stores the current phase space attributes of the particle in
! the secondary bank and increments the number of sites in the secondary bank.
!===============================================================================
subroutine create_secondary(this, uvw, type)
class(Particle), intent(inout) :: this
real(8), intent(in) :: uvw(3)
integer, intent(in) :: type
integer :: n
! Check to make sure that the hard-limit on secondary particles is not
! exceeded.
if (this % n_secondary == MAX_SECONDARY) then
call fatal_error("Too many secondary particles created.")
end if
n = this % n_secondary + 1
this % secondary_bank(n) % wgt = this % wgt
this % secondary_bank(n) % xyz(:) = this % coord(1) % xyz
this % secondary_bank(n) % uvw(:) = uvw
this % secondary_bank(n) % E = this % E
this % n_secondary = n
end subroutine create_secondary
end module particle_header

7
src/particle_restart_write.F90
View file

@ -43,12 +43,7 @@ contains
call pr % file_create(filename)
! Get information about source particle
select case (run_mode)
case (MODE_EIGENVALUE)
src => source_bank(current_work)
case (MODE_FIXEDSOURCE)
src => source_site
end select
src => source_bank(current_work)
! Write data to file
call pr % write_data(FILETYPE_PARTICLE_RESTART, 'filetype')

72
src/physics.F90
View file

@ -387,8 +387,7 @@ contains
end do
! Perform collision physics for inelastic scattering
call inelastic_scatter(nuc, rxn, p % E, p % coord(1) % uvw, &
p % mu, p % wgt, materials(p % material) % p0(i_nuc_mat))
call inelastic_scatter(nuc, rxn, p, materials(p % material) % p0(i_nuc_mat))
p % event_MT = rxn % MT
end if
@ -439,7 +438,7 @@ contains
! Sample velocity of target nucleus
if (.not. micro_xs(i_nuclide) % use_ptable) then
call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, &
& micro_xs(i_nuclide) % elastic)
& micro_xs(i_nuclide) % elastic)
else
v_t = ZERO
end if
@ -605,7 +604,7 @@ contains
! accompanying PDF and CDF is utilized)
if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. &
(sab % secondary_mode == SAB_SECONDARY_SKEWED)) then
(sab % secondary_mode == SAB_SECONDARY_SKEWED)) then
if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then
! All bins equally likely
@ -877,16 +876,16 @@ contains
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
xs_up = nuc % elastic_0K(i_E_up)
m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
! get max 0K xs value over range of practical relative energies
xs_max = max(xs_low, &
& maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up)
& maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up)
reject = .true.
@ -932,26 +931,26 @@ contains
! cdf value at lower bound attainable energy
if (i_E_low > 1) then
m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = nuc % xs_cdf(i_E_low - 1) &
& + m * (E_low - nuc % energy_0K(i_E_low))
& + m * (E_low - nuc % energy_0K(i_E_low))
else
m = nuc % xs_cdf(i_E_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = m * (E_low - nuc % energy_0K(i_E_low))
if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO
end if
! cdf value at upper bound attainable energy
m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
cdf_up = nuc % xs_cdf(i_E_up - 1) &
& + m * (E_up - nuc % energy_0K(i_E_up))
& + m * (E_up - nuc % energy_0K(i_E_up))
! values used to sample the Maxwellian
E_mode = kT
p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_mode / kT)
& * exp(-E_mode / kT)
E_t_max = 16.0_8 * E_mode
reject = .true.
@ -961,7 +960,7 @@ contains
! perform Maxwellian rejection sampling
E_t = E_t_max * prn()**2
p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_t / kT)
& * exp(-E_t / kT)
R_speed = p_t / p_mode
if (prn() < R_speed) then
@ -969,12 +968,12 @@ contains
! sample a relative energy using the xs cdf
cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
& i_E_up - i_E_low + 2, cdf_rel)
& i_E_up - i_E_low + 2, cdf_rel)
E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
& - nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
& / (nuc % energy_0K(i_E_low + i_E_rel) &
& - nuc % energy_0K(i_E_low + i_E_rel - 1))
& - nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
& / (nuc % energy_0K(i_E_low + i_E_rel) &
& - nuc % energy_0K(i_E_low + i_E_rel - 1))
E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
! perform rejection sampling on cosine between
@ -1060,7 +1059,7 @@ contains
! Determine rejection probability
accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) &
/(beta_vn + beta_vt)
/(beta_vn + beta_vt)
! Perform rejection sampling on vt and mu
if (prn() < accept_prob) exit
@ -1302,18 +1301,18 @@ contains
! than fission), i.e. level scattering, (n,np), (n,na), etc.
!===============================================================================
subroutine inelastic_scatter(nuc, rxn, E, uvw, mu, wgt, iso_lab)
subroutine inelastic_scatter(nuc, rxn, p, iso_lab)
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
real(8), intent(inout) :: E ! energy in lab (incoming/outgoing)
real(8), intent(inout) :: uvw(3) ! directional cosines
real(8), intent(out) :: mu ! cosine of scattering angle in lab
real(8), intent(inout) :: wgt ! particle weight
type(Particle), intent(inout) :: p
logical, intent(in) :: iso_lab
integer :: i ! loop index
integer :: law ! secondary energy distribution law
real(8) :: mu ! cosine of scattering angle in lab
real(8) :: A ! atomic weight ratio of nuclide
real(8) :: E ! energy in lab (incoming/outgoing)
real(8) :: E_in ! incoming energy
real(8) :: E_cm ! outgoing energy in center-of-mass
real(8) :: Q ! Q-value of reaction
@ -1322,8 +1321,8 @@ contains
real(8) :: phi ! azimuthal angle
! copy energy, direction of neutron
E_in = E
uvw_in(:) = uvw(:)
E_in = p % E
uvw_in(:) = p % coord(1) % uvw(:)
! determine A and Q
A = nuc % awr
@ -1360,22 +1359,29 @@ contains
! compute cosine of scattering angle in LAB frame by taking dot product of
! neutron's pre- and post-collision unit vectors
if (iso_lab) then
uvw(1) = TWO * prn() - ONE
p % coord(1) % uvw(1) = TWO * prn() - ONE
phi = TWO * PI * prn()
uvw(2) = cos(phi) * sqrt(ONE - uvw(1)*uvw(1))
uvw(3) = sin(phi) * sqrt(ONE - uvw(1)*uvw(1))
mu = dot_product(uvw_in, uvw)
p % coord(1) % uvw(2) = cos(phi) * sqrt(ONE - uvw_in(1) * uvw_in(1))
p % coord(1) % uvw(3) = sin(phi) * sqrt(ONE - uvw_in(1) * uvw_in(1))
mu = dot_product(uvw_in, p % coord(1) % uvw)
else
uvw = rotate_angle(uvw_in, mu)
! Set outgoing energy and scattering angle
p % E = E
p % mu = mu
! change direction of particle
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
end if
! change weight of particle based on yield
if (rxn % multiplicity_with_E) then
yield = interpolate_tab1(rxn % multiplicity_E, E_in)
p % wgt = yield * p % wgt
else
yield = rxn % multiplicity
do i = 1, rxn % multiplicity - 1
call p % create_secondary(p % coord(1) % uvw, NEUTRON)
end do
end if
wgt = yield * wgt
end subroutine inelastic_scatter

2
src/plot.F90
View file

@ -250,7 +250,7 @@ contains
do j = ijk_ll(inner), ijk_ur(inner)
! check if we're in the mesh for this ijk
if (i > 0 .and. i <= m % dimension(outer) .and. &
j > 0 .and. j <= m % dimension(inner)) then
j > 0 .and. j <= m % dimension(inner)) then
! get xyz's of lower left and upper right of this mesh cell
xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1)

4
src/plot_header.F90
View file

@ -24,7 +24,7 @@ module plot_header
integer :: color_by ! quantity to color regions by
real(8) :: origin(3) ! xyz center of plot location
real(8) :: width(3) ! xyz widths of plot
integer :: basis ! direction of plot slice
integer :: basis ! direction of plot slice
integer :: pixels(3) ! pixel width/height of plot slice
integer :: meshlines_width ! pixel width of meshlines
integer :: level ! universe depth to plot the cells of
@ -37,7 +37,7 @@ module plot_header
! Plot type
integer, parameter :: PLOT_TYPE_SLICE = 1
integer, parameter :: PLOT_TYPE_VOXEL = 2
! Plot level
integer, parameter :: PLOT_LEVEL_LOWEST = -1

12
src/ppmlib.F90
View file

@ -1,7 +1,7 @@
module ppmlib
implicit none
!===============================================================================
! Image holds RGB information for output PPM image
!===============================================================================
@ -12,7 +12,7 @@ module ppmlib
end type Image
contains
!===============================================================================
! INIT_IMAGE initializes the Image derived type
!===============================================================================
@ -55,7 +55,7 @@ contains
!===============================================================================
! DEALLOCATE_IMAGE
!===============================================================================
subroutine deallocate_image(img)
type(Image) :: img
@ -70,7 +70,7 @@ contains
! INSIDE_IMAGE determines whether a point (x,y) is inside the image
!===============================================================================
function inside_image(img, x, y) result(inside)
type(Image), intent(in) :: img
@ -83,7 +83,7 @@ contains
(x >= 0) .and. (y >= 0)) inside = .true.
end function inside_image
!===============================================================================
! VALID_IMAGE checks whether the image has a width and height and if its color
! arrays are allocated
@ -123,5 +123,5 @@ contains
end if
end subroutine set_pixel
end module ppmlib

20
src/progress_header.F90
View file

@ -22,12 +22,12 @@ module progress_header
private
character(len=72) :: bar="???% | " // &
" |"
contains
procedure :: set_value => bar_set_value
contains
procedure :: set_value => bar_set_value
end type ProgressBar
contains
!===============================================================================
! IS_TERMINAL checks if output is currently being output to a terminal. Relies
! on a POSIX implementation of isatty, and defaults to false if that is not
@ -46,7 +46,7 @@ contains
#endif
end function is_terminal
!===============================================================================
! BAR_SET_VALUE prints the progress bar without advancing. The value is
! specified as percent completion, from 0 to 100. If the value is ever set to
@ -57,7 +57,7 @@ contains
class(ProgressBar), intent(inout) :: self
real(8), intent(in) :: val
integer :: i
if (.not. is_terminal()) return
@ -84,19 +84,19 @@ contains
write(OUTPUT_UNIT, '(A1,A1,A72)', ADVANCE='no') '+', char(13), self % bar
flush(OUTPUT_UNIT)
if (val >= 100.) then
! make new line
write(OUTPUT_UNIT, "(A)") ""
flush(OUTPUT_UNIT)
! reset the bar in case we want to use this instance again
self % bar = "???% | " // &
" |"
end if
end subroutine bar_set_value
end module progress_header

2
src/random_lcg.F90
View file

@ -93,7 +93,7 @@ contains
end do
end subroutine set_particle_seed
!===============================================================================
! PRN_SKIP advances the random number seed 'n' times from the current seed
!===============================================================================

6
src/search.F90
View file

@ -62,7 +62,7 @@ contains
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
&search.")
end if
end do
@ -114,7 +114,7 @@ contains
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
&search.")
end if
end do
@ -166,7 +166,7 @@ contains
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
&search.")
end if
end do

371
src/simulation.F90 Normal file
View file

@ -0,0 +1,371 @@
module simulation
#ifdef MPI
use mpi
#endif
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use eigenvalue, only: count_source_for_ufs, calculate_average_keff, &
calculate_combined_keff, calculate_generation_keff, &
shannon_entropy, synchronize_bank, keff_generation
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
#endif
use global
use output, only: write_message, header, print_columns, &
print_batch_keff, print_generation
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use source, only: initialize_source
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use trigger, only: check_triggers
use tracking, only: transport
implicit none
private
public :: run_simulation
contains
!===============================================================================
! RUN_SIMULATION encompasses all the main logic where iterations are performed
! over the batches, generations, and histories in a fixed source or k-eigenvalue
! calculation.
!===============================================================================
subroutine run_simulation()
type(Particle) :: p
integer(8) :: i_work
if (.not. restart_run) call initialize_source()
! Display header
if (master) then
if (run_mode == MODE_FIXEDSOURCE) then
call header("FIXED SOURCE TRANSPORT SIMULATION", level=1)
elseif (run_mode == MODE_EIGENVALUE) then
call header("K EIGENVALUE SIMULATION", level=1)
call print_columns()
end if
end if
! Turn on inactive timer
call time_inactive % start()
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_max_batches
call initialize_batch()
! Handle restart runs
if (restart_run .and. current_batch <= restart_batch) then
call replay_batch_history()
cycle BATCH_LOOP
end if
! =======================================================================
! LOOP OVER GENERATIONS
GENERATION_LOOP: do current_gen = 1, gen_per_batch
call initialize_generation()
! Start timer for transport
call time_transport % start()
! ====================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
! grab source particle from bank
call initialize_history(p, current_work)
! transport particle
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
call finalize_generation()
end do GENERATION_LOOP
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
end do BATCH_LOOP
call time_active % stop()
! ==========================================================================
! END OF RUN WRAPUP
if (master) call header("SIMULATION FINISHED", level=1)
! Clear particle
call p % clear()
end subroutine run_simulation
!===============================================================================
! INITIALIZE_HISTORY
!===============================================================================
subroutine initialize_history(p, index_source)
type(Particle), intent(inout) :: p
integer(8), intent(in) :: index_source
integer(8) :: particle_seed ! unique index for particle
integer :: i
! set defaults
call p % initialize_from_source(source_bank(index_source))
! set identifier for particle
p % id = work_index(rank) + index_source
! set random number seed
particle_seed = (overall_gen - 1)*n_particles + p % id
call set_particle_seed(particle_seed)
! set particle trace
trace = .false.
if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
p % id == trace_particle) trace = .true.
! Set particle track.
p % write_track = .false.
if (write_all_tracks) then
p % write_track = .true.
else if (allocated(track_identifiers)) then
do i=1, size(track_identifiers(1,:))
if (current_batch == track_identifiers(1,i) .and. &
&current_gen == track_identifiers(2,i) .and. &
&p % id == track_identifiers(3,i)) then
p % write_track = .true.
exit
end if
end do
end if
end subroutine initialize_history
!===============================================================================
! INITIALIZE_BATCH
!===============================================================================
subroutine initialize_batch()
if (run_mode == MODE_FIXEDSOURCE) then
call write_message("Simulating batch " // trim(to_str(current_batch)) &
// "...", 1)
end if
! Reset total starting particle weight used for normalizing tallies
total_weight = ZERO
if (current_batch == n_inactive + 1) then
! Switch from inactive batch timer to active batch timer
call time_inactive % stop()
call time_active % start()
! Enable active batches (and tallies_on if it hasn't been enabled)
active_batches = .true.
tallies_on = .true.
! Add user tallies to active tallies list
!$omp parallel
call setup_active_usertallies()
!$omp end parallel
end if
! check CMFD initialize batch
if (run_mode == MODE_EIGENVALUE) then
if (cmfd_run) call cmfd_init_batch()
end if
end subroutine initialize_batch
!===============================================================================
! INITIALIZE_GENERATION
!===============================================================================
subroutine initialize_generation()
! set overall generation number
overall_gen = gen_per_batch*(current_batch - 1) + current_gen
if (run_mode == MODE_EIGENVALUE) then
! Reset number of fission bank sites
n_bank = 0
! Count source sites if using uniform fission source weighting
if (ufs) call count_source_for_ufs()
! Store current value of tracklength k
keff_generation = global_tallies(K_TRACKLENGTH) % value
end if
end subroutine initialize_generation
!===============================================================================
! FINALIZE_GENERATION
!===============================================================================
subroutine finalize_generation()
! Update global tallies with the omp private accumulation variables
!$omp parallel
!$omp critical
if (run_mode == MODE_EIGENVALUE) then
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + global_tally_collision
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + global_tally_absorption
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength
end if
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + global_tally_leakage
!$omp end critical
! reset private tallies
if (run_mode == MODE_EIGENVALUE) then
global_tally_collision = 0
global_tally_absorption = 0
global_tally_tracklength = 0
end if
global_tally_leakage = 0
!$omp end parallel
if (run_mode == MODE_EIGENVALUE) then
#ifdef _OPENMP
! Join the fission bank from each thread into one global fission bank
call join_bank_from_threads()
#endif
! Distribute fission bank across processors evenly
call time_bank % start()
call synchronize_bank()
call time_bank % stop()
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Collect results and statistics
call calculate_generation_keff()
call calculate_average_keff()
! Write generation output
if (master .and. current_gen /= gen_per_batch) call print_generation()
end if
end subroutine finalize_generation
!===============================================================================
! FINALIZE_BATCH handles synchronization and accumulation of tallies,
! calculation of Shannon entropy, getting single-batch estimate of keff, and
! turning on tallies when appropriate
!===============================================================================
subroutine finalize_batch()
! Collect tallies
call time_tallies % start()
call synchronize_tallies()
call time_tallies % stop()
! Reset global tally results
if (.not. active_batches) then
call reset_result(global_tallies)
n_realizations = 0
end if
if (run_mode == MODE_EIGENVALUE) then
! Perform CMFD calculation if on
if (cmfd_on) call execute_cmfd()
! Display output
if (master) call print_batch_keff()
! Calculate combined estimate of k-effective
if (master) call calculate_combined_keff()
end if
! Check_triggers
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call write_state_point()
end if
! Write out source point if it's been specified for this batch
if ((sourcepoint_batch % contains(current_batch) .or. source_latest) .and. &
source_write) then
call write_source_point()
end if
if (master .and. current_batch == n_max_batches .and. &
run_mode == MODE_EIGENVALUE) then
! Make sure combined estimate of k-effective is calculated at the last
! batch in case no state point is written
call calculate_combined_keff()
end if
end subroutine finalize_batch
!===============================================================================
! REPLAY_BATCH_HISTORY displays keff and entropy for each generation within a
! batch using data read from a state point file
!===============================================================================
subroutine replay_batch_history
! Write message at beginning
if (current_batch == 1) then
call write_message("Replaying history from state point...", 1)
end if
if (run_mode == MODE_EIGENVALUE) then
do current_gen = 1, gen_per_batch
overall_gen = overall_gen + 1
call calculate_average_keff()
! print out batch keff
if (current_gen < gen_per_batch) then
if (master) call print_generation()
else
if (master) call print_batch_keff()
end if
end do
end if
! Write message at end
if (current_batch == restart_batch) then
call write_message("Resuming simulation...", 1)
end if
end subroutine replay_batch_history
end module simulation

75
src/source.F90
View file

@ -78,7 +78,7 @@ contains
! Write out initial source
if (write_initial_source) then
call write_message('Writing out initial source guess...', 1)
call write_message('Writing out initial source...', 1)
#ifdef HDF5
filename = trim(path_output) // 'initial_source.h5'
#else
@ -92,7 +92,8 @@ contains
end subroutine initialize_source
!===============================================================================
! SAMPLE_EXTERNAL_SOURCE
! SAMPLE_EXTERNAL_SOURCE samples the user-specified external source and stores
! the position, angle, and energy in a Bank type.
!===============================================================================
subroutine sample_external_source(site)
@ -245,74 +246,4 @@ contains
end subroutine sample_external_source
!===============================================================================
! GET_SOURCE_PARTICLE returns the next source particle
!===============================================================================
subroutine get_source_particle(p, index_source)
type(Particle), intent(inout) :: p
integer(8), intent(in) :: index_source
integer(8) :: particle_seed ! unique index for particle
integer :: i
type(Bank), pointer :: src
! set defaults
call p % initialize()
! Copy attributes from source to particle
src => source_bank(index_source)
call copy_source_attributes(p, src)
! set identifier for particle
p % id = work_index(rank) + index_source
! set random number seed
particle_seed = (overall_gen - 1)*n_particles + p % id
call set_particle_seed(particle_seed)
! set particle trace
trace = .false.
if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
p % id == trace_particle) trace = .true.
! Set particle track.
p % write_track = .false.
if (write_all_tracks) then
p % write_track = .true.
else if (allocated(track_identifiers)) then
do i=1, size(track_identifiers(1,:))
if (current_batch == track_identifiers(1,i) .and. &
&current_gen == track_identifiers(2,i) .and. &
&p % id == track_identifiers(3,i)) then
p % write_track = .true.
exit
end if
end do
end if
end subroutine get_source_particle
!===============================================================================
! COPY_SOURCE_ATTRIBUTES
!===============================================================================
subroutine copy_source_attributes(p, src)
type(Particle), intent(inout) :: p
type(Bank), pointer :: src
! copy attributes from source bank site
p % wgt = src % wgt
p % last_wgt = src % wgt
p % coord(1) % xyz = src % xyz
p % coord(1) % uvw = src % uvw
p % last_xyz = src % xyz
p % last_uvw = src % uvw
p % E = src % E
p % last_E = src % E
end subroutine copy_source_attributes
end module source

28
src/state_point.F90
View file

@ -52,7 +52,7 @@ contains
! Set filename for state point
filename = trim(path_output) // 'statepoint.' // &
& zero_padded(current_batch, count_digits(n_max_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
! Append appropriate extension
#ifdef HDF5
@ -256,7 +256,7 @@ contains
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j))
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT) then
tally % filters(j) % type == FILTER_ENERGYOUT) then
call sp % write_data(tally % filters(j) % real_bins, "bins", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j), &
@ -309,8 +309,8 @@ contains
do n_order = 0, tally % moment_order(k)
moment_name = 'P' // trim(to_str(n_order))
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
end do
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
@ -318,7 +318,7 @@ contains
do n_order = 0, tally % moment_order(k)
do nm_order = -n_order, n_order
moment_name = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order))
trim(to_str(nm_order))
call sp % write_data(moment_name, "order" // &
trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
@ -412,7 +412,7 @@ contains
! Set filename
filename = trim(path_output) // 'source.' // &
& zero_padded(current_batch, count_digits(n_max_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
#ifdef HDF5
filename = trim(filename) // '.h5'
@ -434,7 +434,7 @@ contains
! Set filename for state point
filename = trim(path_output) // 'statepoint.' // &
& zero_padded(current_batch, count_digits(n_max_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
#ifdef HDF5
filename = trim(filename) // '.h5'
#else
@ -607,12 +607,12 @@ contains
tally % results(:,:) % sum_sq = tally_temp(2,:,:)
end if
! Put in temporary tally result
allocate(tallyresult_temp(m,n))
tallyresult_temp(:,:) % sum = tally_temp(1,:,:)
tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:)
! Put in temporary tally result
allocate(tallyresult_temp(m,n))
tallyresult_temp(:,:) % sum = tally_temp(1,:,:)
tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:)
! Write reduced tally results to file
! Write reduced tally results to file
call sp % write_tally_result(tally % results, "results", &
group="tallies/tally " // trim(to_str(tally % id)), n1=m, n2=n)
@ -687,7 +687,7 @@ contains
call sp % read_data(int_array(2), "version_minor")
call sp % read_data(int_array(3), "version_release")
if (int_array(1) /= VERSION_MAJOR .or. int_array(2) /= VERSION_MINOR &
.or. int_array(3) /= VERSION_RELEASE) then
.or. int_array(3) /= VERSION_RELEASE) then
if (master) call warning("State point file was created with a different &
&version of OpenMC.")
end if
@ -843,7 +843,7 @@ contains
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j))
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT) then
tally % filters(j) % type == FILTER_ENERGYOUT) then
call sp % read_data(tally % filters(j) % real_bins, "bins", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j), &

4
src/string.F90
View file

@ -7,7 +7,7 @@ module string
implicit none
interface to_str
module procedure int4_to_str, int8_to_str, real_to_str
module procedure int4_to_str, int8_to_str, real_to_str
end interface
contains
@ -50,7 +50,7 @@ contains
n = n + 1
if (i_end - i_start + 1 > len(words(n))) then
if (master) call warning("The word '" // string(i_start:i_end) &
&// "' is longer than the space allocated for it.")
&// "' is longer than the space allocated for it.")
end if
words(n) = string(i_start:i_end)
! reset indices

6
src/tally_header.F90
View file

@ -124,7 +124,7 @@ module tally_header
! Number of realizations of tally random variables
integer :: n_realizations = 0
! Tally precision triggers
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
@ -197,10 +197,10 @@ module tally_header
this % reset = .false.
this % n_realizations = 0
if (allocated(this % triggers)) &
deallocate (this % triggers)
this % n_triggers = 0
end subroutine tallyobject_clear

12
src/timer_header.F90
View file

@ -15,11 +15,11 @@ module timer_header
logical :: running = .false. ! is timer running?
integer :: start_counts = 0 ! counts when started
real(8), public :: elapsed = ZERO ! total time elapsed in seconds
contains
procedure :: start => timer_start
procedure :: get_value => timer_get_value
procedure :: stop => timer_stop
procedure :: reset => timer_reset
contains
procedure :: start => timer_start
procedure :: get_value => timer_get_value
procedure :: stop => timer_stop
procedure :: reset => timer_reset
end type Timer
contains
@ -83,7 +83,7 @@ contains
!===============================================================================
subroutine timer_reset(self)
class(Timer), intent(inout) :: self
self % running = .false.

80
src/track_output.F90
View file

@ -12,18 +12,22 @@ module track_output
implicit none
integer, private :: n_tracks ! total number of tracks
real(8), private, allocatable :: coords(:,:) ! track coordinates
!$omp threadprivate(n_tracks, coords)
type, private :: TrackCoordinates
real(8), allocatable :: coords(:,:)
end type TrackCoordinates
type(TrackCoordinates), private, allocatable :: tracks(:)
!$omp threadprivate(tracks)
contains
!===============================================================================
! INITIALIZE_PARTICLE_TRACK
! INITIALIZE_PARTICLE_TRACK allocates the array to store particle track
! information
!===============================================================================
subroutine initialize_particle_track()
n_tracks = 0
allocate(tracks(1))
end subroutine initialize_particle_track
!===============================================================================
@ -32,23 +36,50 @@ contains
subroutine write_particle_track(p)
type(Particle), intent(in) :: p
real(8), allocatable :: new_coords(:, :)
integer :: i
integer :: n_tracks
! Add another column to coords
n_tracks = n_tracks + 1
if (allocated(coords)) then
allocate(new_coords(3, n_tracks))
new_coords(:, 1:n_tracks-1) = coords
call move_alloc(FROM=new_coords, TO=coords)
i = size(tracks)
if (allocated(tracks(i) % coords)) then
n_tracks = size(tracks(i) % coords, 2)
allocate(new_coords(3, n_tracks + 1))
new_coords(:, 1:n_tracks) = tracks(i) % coords
call move_alloc(FROM=new_coords, TO=tracks(i) % coords)
else
allocate(coords(3,1))
n_tracks = 0
allocate(tracks(i) % coords(3, 1))
end if
! Write current coordinates into the newest column.
coords(:, n_tracks) = p % coord(1) % xyz
n_tracks = n_tracks + 1
tracks(i) % coords(:, n_tracks) = p % coord(1) % xyz
end subroutine write_particle_track
!===============================================================================
! ADD_PARTICLE_TRACK creates a new entry in the track coordinates for a
! secondary particle
!===============================================================================
subroutine add_particle_track()
type(TrackCoordinates), allocatable :: new_tracks(:)
integer :: i
! Determine current number of particle tracks
i = size(tracks)
! Create array one larger than current
allocate(new_tracks(i + 1))
! Copy memory and move allocation
new_tracks(1:i) = tracks(i)
call move_alloc(FROM=new_tracks, TO=tracks)
end subroutine add_particle_track
!===============================================================================
! FINALIZE_PARTICLE_TRACK writes the particle track array to disk.
!===============================================================================
@ -60,6 +91,10 @@ contains
character(MAX_FILE_LEN) :: fname
type(BinaryOutput) :: binout
integer :: i
integer, allocatable :: n_coords(:)
integer :: n_particle_tracks
#ifdef HDF5
fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) &
// '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) &
@ -69,13 +104,26 @@ contains
// '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) &
// '.binary'
#endif
! Determine total number of particles and number of coordinates for each
n_particle_tracks = size(tracks)
allocate(n_coords(n_particle_tracks))
do i = 1, n_particle_tracks
n_coords(i) = size(tracks(i) % coords, 2)
end do
!$omp critical (FinalizeParticleTrack)
call binout % file_create(fname)
length = [3, n_tracks]
call binout % write_data(coords, 'coordinates', length=length)
call binout % write_data(n_particle_tracks, 'n_particles')
call binout % write_data(n_coords, 'n_coords', length=n_particle_tracks)
do i = 1, n_particle_tracks
length(:) = [3, n_coords(i)]
call binout % write_data(tracks(i) % coords, 'coordinates_' // &
trim(to_str(i)), length=length)
end do
call binout % file_close()
!$omp end critical (FinalizeParticleTrack)
deallocate(coords)
deallocate(tracks)
end subroutine finalize_particle_track
end module track_output

44
src/tracking.F90
View file

@ -15,7 +15,7 @@ module tracking
use tally, only: score_analog_tally, score_tracklength_tally, &
score_surface_current
use track_output, only: initialize_particle_track, write_particle_track, &
finalize_particle_track
add_particle_track, finalize_particle_track
implicit none
@ -45,20 +45,6 @@ contains
call write_message("Simulating Particle " // trim(to_str(p % id)))
end if
! If the cell hasn't been determined based on the particle's location,
! initiate a search for the current cell
if (p % coord(p % n_coord) % cell == NONE) then
call find_cell(p, found_cell)
! Particle couldn't be located
if (.not. found_cell) then
call fatal_error("Could not locate particle " // trim(to_str(p % id)))
end if
! set birth cell attribute
p % cell_born = p % coord(p % n_coord) % cell
end if
! Initialize number of events to zero
n_event = 0
@ -74,7 +60,19 @@ contains
call initialize_particle_track()
endif
do while (p % alive)
EVENT_LOOP: do
! If the cell hasn't been determined based on the particle's location,
! initiate a search for the current cell. This generally happens at the
! beginning of the history and again for any secondary particles
if (p % coord(p % n_coord) % cell == NONE) then
call find_cell(p, found_cell)
if (.not. found_cell) then
call fatal_error("Could not locate particle " // trim(to_str(p % id)))
end if
! set birth cell attribute
if (p % cell_born == NONE) p % cell_born = p % coord(p % n_coord) % cell
end if
! Write particle track.
if (p % write_track) call write_particle_track(p)
@ -197,7 +195,19 @@ contains
p % alive = .false.
end if
end do
! Check for secondary particles if this particle is dead
if (.not. p % alive) then
if (p % n_secondary > 0) then
call p % initialize_from_source(p % secondary_bank(p % n_secondary))
p % n_secondary = p % n_secondary - 1
! Enter new particle in particle track file
if (p % write_track) call add_particle_track()
else
exit EVENT_LOOP
end if
end if
end do EVENT_LOOP
! Finish particle track output.
if (p % write_track) then

6
src/trigger.F90
View file

@ -43,7 +43,7 @@ contains
! When trigger threshold is reached, write information
if (satisfy_triggers) then
call write_message("Triggers satisfied for batch " // &
trim(to_str(current_batch)))
trim(to_str(current_batch)))
! When trigger is not reached write convergence info for user
elseif (name == "eigenvalue") then
@ -83,7 +83,7 @@ contains
! and finds the maximum uncertainty/threshold ratio for all triggers
!===============================================================================
subroutine check_tally_triggers(max_ratio, tally_id, name)
subroutine check_tally_triggers(max_ratio, tally_id, name)
! Variables to reflect distance to trigger convergence criteria
real(8), intent(inout) :: max_ratio ! max uncertainty/thresh ratio
@ -299,7 +299,7 @@ contains
! precision trigger(s).
!===============================================================================
subroutine compute_tally_current(t, trigger)
subroutine compute_tally_current(t, trigger)
integer :: i ! mesh index for x
integer :: j ! mesh index for y

6
src/trigger_header.F90
View file

@ -14,16 +14,16 @@ module trigger_header
character(len=52) :: score_name ! the name of the score
integer :: score_index ! the index of the score
real(8) :: variance=0.0 ! temp variance container
real(8) :: std_dev =0.0 ! temp std. dev. container
real(8) :: std_dev =0.0 ! temp std. dev. container
real(8) :: rel_err =0.0 ! temp rel. err. container
end type TriggerObject
!===============================================================================
! KTRIGGER describes a user-specified precision trigger for k-effective
!===============================================================================
type KTrigger
integer :: trigger_type = 0
real(8) :: threshold = 0
real(8) :: threshold = 0
end type KTrigger
end module trigger_header

12
src/vector_header.F90
View file

@ -9,12 +9,12 @@ module vector_header
integer :: n ! number of rows/cols in matrix
real(8), allocatable :: data(:) ! where vector data is stored
real(8), pointer :: val(:) ! pointer to vector data
contains
procedure :: create => vector_create
procedure :: destroy => vector_destroy
procedure :: add_value => vector_add_value
procedure :: copy => vector_copy
! TODO: procedure :: write => vector_write
contains
procedure :: create => vector_create
procedure :: destroy => vector_destroy
procedure :: add_value => vector_add_value
procedure :: copy => vector_copy
! TODO: procedure :: write => vector_write
end type Vector
contains

50
src/xml_interface.F90
View file

@ -92,7 +92,7 @@ contains
! Check if node exists
if (associated(temp_ptr)) return
! Check for a sub-element
! Check for a sub-element
elem_list => getChildrenByTagName(ptr, trim(node_name))
! Get the length of the list
@ -122,7 +122,7 @@ contains
! Set found to false
found_ = .false.
! Check for a sub-element
! Check for a sub-element
elem_list => getChildrenByTagName(in_ptr, trim(node_name))
! Get the length of the list
@ -147,7 +147,7 @@ contains
type(Node), pointer, intent(in) :: in_ptr
type(NodeList), pointer, intent(out) :: out_ptr
! Check for a sub-element
! Check for a sub-element
out_ptr => getChildrenByTagName(in_ptr, trim(node_name))
end subroutine get_node_list
@ -176,7 +176,7 @@ contains
type(NodeList), pointer, intent(in) :: in_ptr
type(Node), pointer, intent(out) :: out_ptr
! Check for a sub-element
! Check for a sub-element
out_ptr => item(in_ptr, idx - 1)
end subroutine get_list_item
@ -199,7 +199,7 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
@ -210,7 +210,7 @@ contains
else
call extractDataContent(temp_ptr, result_int)
end if
end subroutine get_node_value_integer
!===============================================================================
@ -231,18 +231,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_long)
else
call extractDataContent(temp_ptr, result_long)
end if
end subroutine get_node_value_long
!===============================================================================
@ -263,18 +263,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_double)
else
call extractDataContent(temp_ptr, result_double)
end if
end subroutine get_node_value_double
!===============================================================================
@ -295,18 +295,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_int)
else
call extractDataContent(temp_ptr, result_int)
end if
end subroutine get_node_array_integer
!===============================================================================
@ -327,18 +327,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_double)
else
call extractDataContent(temp_ptr, result_double)
end if
end subroutine get_node_array_double
!===============================================================================
@ -359,18 +359,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_string)
else
call extractDataContent(temp_ptr, result_string)
end if
end subroutine get_node_array_string
!===============================================================================
@ -391,18 +391,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " // &
getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_str)
else
call extractDataContent(temp_ptr, result_str)
end if
end subroutine get_node_value_string
!===============================================================================
@ -513,7 +513,7 @@ contains
! Check if node exists
if (associated(out_ptr)) return
! Check for a sub-element
! Check for a sub-element
elem_list => getChildrenByTagName(in_ptr, trim(node_name))
! Get the length of the list

232
tests/check_source.py Executable file
View file

@ -0,0 +1,232 @@
#!/usr/bin/env python
from __future__ import print_function
import glob
from string import whitespace
from sys import exit
from textwrap import TextWrapper
################################################################################
# Fortran reading/parsing backend.
################################################################################
class LineOfCode(object):
"""Contains and provides basic info about a line of Fortran 90 code."""
def __init__(self, number, content, is_continuation=False,
string_terminator=None):
# Initialize member variables.
self.number = number
self.content = content
self.is_continuation = is_continuation
self.is_continued = False
assert (string_terminator is None or
string_terminator == "'" or
string_terminator == '"')
self.initial_string_terminator = string_terminator
self.final_string_terminator = None
self.is_comment_only = False
# Parse the string.
self.parse()
def __repr__(self):
rep = 'LineOfCode: line number = {0:d}\n'.format(self.number)
rep += '\tis_continuation = ' + str(self.is_continuation) + '\n'
rep += '\tis_continued = ' + str(self.is_continued) + '\n'
rep += ('\tinitial_string_terminator = '
+ str(self.initial_string_terminator) + '\n')
rep += ('\tfinal_string_terminator = '
+ str(self.final_string_terminator) + '\n')
rep += '\tis_comment_only = ' + str(self.is_comment_only) + '\n'
rep += '\tContent:\n' + self.content
return rep
def __str__(self):
return repr(self)
def parse(self):
# Initialize the string variables.
if self.initial_string_terminator == "'":
in_a_single_quote = True
in_a_double_quote = False
elif self.initial_string_terminator == '"':
in_a_single_quote = False
in_a_double_quote = True
else:
in_a_single_quote = False
in_a_double_quote = False
# Initialize other variables.
in_a_comment = False
ends_with_ampersand = False
has_real_content = False
# Parse through the line.
for char in self.content:
# Check for the start or end of strings.
if char == "'" and in_a_single_quote:
in_a_single_quote = False
elif char == "'" and not in_a_double_quote:
in_a_single_quote = True
elif char == '"' and in_a_double_quote:
in_a_double_quote = False
elif char == '"' and not in_a_single_quote:
in_a_double_quote = True
# Check for the start of a comment.
if (char == "!" and not in_a_single_quote
and not in_a_double_quote):
in_a_comment = True
break # Don't care about anything in a comment.
# Check for a continuation ampersand.
if char == "&":
ends_with_ampersand = True
elif all([char != w for w in whitespace]):
ends_with_ampersand = False
# Check to see if there is any real content (non-whitespace, not in
# a comment)
if (not in_a_comment) and all([char != w for w in whitespace]):
has_real_content = True
# Make sure nothing went terribly wrong.
assert not (in_a_single_quote and in_a_double_quote)
assert not (in_a_single_quote and in_a_comment)
assert not (in_a_double_quote and in_a_comment)
# Is this line continued onto the next?
if ends_with_ampersand:
self.is_continued = True
# Is a multiline string continued on the next line?
if in_a_single_quote:
assert self.is_continued
self.final_string_terminator = "'"
if in_a_double_quote:
assert self.is_continued
self.final_string_terminator = '"'
# Is this line a comment-only line?
if (not has_real_content) and in_a_comment:
self.is_comment_only = True
def get_indent(self):
"""Return the number of indentation spaces prefixing this line."""
return len(self.content) - len(self.content.lstrip(' '))
def contains_whitespace_only(self):
"""Return True/False if all characters in the line are whitespace."""
if len(self.content.strip('\n')) == 0: return False
is_ws = [ any([char == ws for ws in whitespace])
for char in self.content.strip('\n') ]
return all(is_ws)
def has_trailing_whitespace(self):
"""Return True/False if this line ends with a whitespace character."""
stripped = self.content.strip('\n')
if len(stripped) == 0: return False
return any([stripped[-1] == ws for ws in whitespace])
def contains_tab(self):
"""Return True/False if a tab character appears in the line."""
return '\t' in self.content
def read_lines_of_code(fname):
line_num = 0
cont = False
str_term = None
with open(fname) as fin:
for line in fin:
loc = LineOfCode(line_num, line, is_continuation=cont,
string_terminator=str_term)
cont = loc.is_continued
str_term = loc.final_string_terminator
line_num += 1
yield loc
################################################################################
# Error checking.
################################################################################
def print_error(fname, line_number, err_msg):
header = "Error in file {0}, line {1}:".format(fname, line_number + 1)
tw = TextWrapper(width=80, initial_indent=' '*4, subsequent_indent=' '*4)
body = '\n'.join(tw.wrap(err_msg))
print(header + '\n' + body + '\n')
def check_source(fname):
"""Make sure the given Fortran source file meets OpenMC standards.
If errors are found, messages will be printed to the screen describing the
error. The function will return True if no errors were found or False
otherwise.
"""
good_code = True
base_indent = 0
for loc in read_lines_of_code(fname):
# Check for extra whitespace errors.
if loc.contains_whitespace_only():
good_code = False
print_error(fname, loc.number, "Line contains whitespace "
"characters but no content. Please remove whitespace.")
elif loc.has_trailing_whitespace():
good_code = False
print_error(fname, loc.number, "Line has trailing whitespace after"
" the content. Please remove trailing whitespace.")
if loc.contains_tab():
good_code = False
print_error(fname, loc.number, "Line contains a tab character. "
"Please replace with single whitespace characters.")
# Check indentation.
current_indent = loc.get_indent()
if ((not loc.is_continuation) and (not loc.is_comment_only) and
(not loc.contains_whitespace_only()) and current_indent % 2 != 0):
good_code = False
print_error(fname, loc.number, "Line is indented an odd number of "
"spaces. All non-continuation lines should be indented an "
"even number of spaces.")
if loc.is_continuation and current_indent < base_indent + 5:
good_code = False
print_error(fname, loc.number, "Continuation lines must be "
"indented by at least 5 spaces, but this line is indented {0}"
" spaces.".format(current_indent - base_indent))
# Set base indentation for next lines.
if not loc.is_continuation:
base_indent = current_indent
return good_code
################################################################################
# Main loop.
################################################################################
if __name__ == '__main__':
# Get a list of the F90 source files.
f90_files = glob.glob('../src/*.F90')
# Make sure we found the source files.
assert len(f90_files) != 0, 'No .F90 source files found.'
# Make sure all the F90 source files meet our standards.
good_code = [check_source(fname) for fname in f90_files]
if not all(good_code):
print("ERROR: The Fortran source code does not meet OpenMC's standards")
exit(-1)
else:
print("SUCCESS! Looks like the Fortran source meets our standards")
exit(0)

2
tests/test_basic/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
3.011726E-01 1.841644E-03
3.021779E-01 3.813358E-03

Some files were not shown because too many files have changed in this diff Show more