diff --git a/LICENSE b/LICENSE index e69de29bb2..8c94262bd2 100644 --- a/LICENSE +++ b/LICENSE @@ -0,0 +1,18 @@ +Copyright (c) 2011-2012 Massachusetts Institute of Technology + +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to +use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software is furnished to do so, +subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER +IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/docs/source/index.rst b/docs/source/index.rst index 828f5335d1..18db6e8628 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -25,6 +25,7 @@ Contents :maxdepth: 1 install + releasenotes/index methods/index usersguide/index developers diff --git a/docs/source/releasenotes/index.rst b/docs/source/releasenotes/index.rst new file mode 100644 index 0000000000..8f6c633421 --- /dev/null +++ b/docs/source/releasenotes/index.rst @@ -0,0 +1,15 @@ +.. _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.4.2 + notes_0.4.1 + notes_0.4.0 diff --git a/docs/source/releasenotes/notes_0.4.0.rst b/docs/source/releasenotes/notes_0.4.0.rst new file mode 100644 index 0000000000..3e0a444e35 --- /dev/null +++ b/docs/source/releasenotes/notes_0.4.0.rst @@ -0,0 +1,36 @@ +.. _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. diff --git a/docs/source/releasenotes/notes_0.4.1.rst b/docs/source/releasenotes/notes_0.4.1.rst new file mode 100644 index 0000000000..2492f23a98 --- /dev/null +++ b/docs/source/releasenotes/notes_0.4.1.rst @@ -0,0 +1,55 @@ +.. _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 diff --git a/docs/source/releasenotes/notes_0.4.2.rst b/docs/source/releasenotes/notes_0.4.2.rst new file mode 100644 index 0000000000..8d4c1f2d2a --- /dev/null +++ b/docs/source/releasenotes/notes_0.4.2.rst @@ -0,0 +1,37 @@ +.. _notes_0.4.2: + +============================== +Release Notes for OpenMC 0.4.2 +============================== + +.. note:: + These release notes are for an upcoming release of OpenMC and are still + subject to change. + +------------------- +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 versions 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 +------------ + +- Ability to specify 'point' source type. +- Ability to change random number seed. +- User's can now specify units='sum' on a 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 +--------- + +- `2a5b9c`_: Fixed regression in plotting. + +.. _2a5b9c: https://github.com/mit-crpg/openmc/commit/2a5b9c diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 89431ba114..3e0aec34e7 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -307,16 +307,23 @@ settings.xml file. The ``criticality`` element indicates that a criticality calculation should be performed. It has the following attributes/sub-elements: - :cycles: - The number of total fission source iterations. + :batches: + The total number of batches, where each batch corresponds to multiple + fission source iterations. Batching is done to eliminate correlation between + realizations of random variables. *Default*: None + :generations_per_batch: + The number of total fission source iterations per batch. + + *Default*: 1 + :inactive: - The number of inactive fission source iterations. In general, the starting - cycles in a criticality calculation can not be used to contribute to tallies - since the fission source distribution and eigenvalue are generally not - converged immediately + The number of inactive batches. In general, the starting cycles in a + criticality calculation can not be used to contribute to tallies since the + fission source distribution and eigenvalue are generally not converged + immediately. *Default*: None @@ -353,6 +360,18 @@ default. This element has the following attributes/sub-elements: *Default*: 1.0 +``energy_grid`` Element +----------------------- + +The ``energy_grid`` element determines the treatment of the energy grid during a +simulation. Setting this element to "nuclide" will cause OpenMC to use a +nuclide's energy grid when determining what points to interpolate between for +determining cross sections (i.e. non-unionized energy grid). To use a unionized +energy grid, set this element to "union". Note that the unionized energy grid +treatment is slightly different than that employed in Serpent. + + *Default*: union + ``entropy`` Element ------------------- @@ -385,6 +404,14 @@ sub-elements and can be set to either "off" or "on". *Default*: on +``seed`` Element +---------------- + +The ``seed`` element is used to set the seed used for the linear congruential +pseudo-random number generator. + + *Default*: 1 + ``source`` Element ------------------ @@ -573,40 +600,115 @@ tallies. This element should be followed by "yes" or "no" Geometry Plotting Specification -- plot.xml ------------------------------------------- -A rudimentary plotting capability is available in OpenMC by specifying a -plot.xml file and subsequently running with the command-line flag ``-plot``. The -root element of the plot.xml is simply ```` and four sub-elements can be -defined to configure the plotting range and resolution. +A basic 2D plotting capability is available in OpenMC by creating a +plots.xml file and subsequently running with the command-line flag ``-plot``. The +root element of the plot.xml is simply ```` and any number output +figures can be defined with ```` sub-elements. -``origin`` Element +``plot`` Element ------------------ -The ``origin`` element has no attributes/sub-elements and indicates the -Cartesian coordinates of the center of the plot. +Each plot must contain a combination of the following attributes or sub-elements: - *Default*: None + :id: + The unique ``id`` of the plot. -``width`` Element ------------------ + *Default*: None - Required entry -The ``width`` element has no attributes/sub-elements and indicates the width of -the plot in each of the basis directions. + :filename: + Filename for the output plot file. - *Default*: None + *Default*: "plot" -``basis`` Element ------------------ + :color: + Keyword for plot coloring. This can only be either ``cell`` or ``mat``, + which colors regions by cells and materials, respectively. -The ``basis`` element has no attributes/sub-elements and indicates the specified -basis for plotting. The only option option currently accepted is "xy". + *Default*: ``cell`` - *Default*: xy + :origin: + Specifies the XYZ coordinate of the center of the plot. Should be 3 floats + separated by spaces. -``pixel`` Element ------------------ + *Default*: None - Required entry -The ``pixel`` element has no attributes/sub-elements and indicates the distance -between horizontal rays sent through the geometry to record surface crossings. A -smaller ``pixel`` will result in a higher-resolution plot. + :width: + Specifies the width of the plot along each of the basis directions. + Should be 2 or 3 floats separated by spaces for 2D plots and 3D plots, + respectively. - *Default*: 0.01 + *Default*: None - Required entry + + :type: + Keyword for type of plot to be produced. Currently only ``slice`` plots are + implemented, which create 2D pixel maps saved in the PPM file format. PPM + files can be displayed in most viewers (e.g. the default Gnome viewer, + IrfanView, etc.). + + .. note:: Since the PPM format is saved without any kind of compression, + the resulting file sizes can be quite large. Saving the image in + the PNG format can often times reduce the file size by orders of + magnitude without any loss of image quality. + + *Default*: "slice" + +``plot`` elements of ``type`` ``slice`` also contain the following attributes or +sub-elements: + + :basis: + Keyword specifying the plane of the plot for ``slice`` type plots. Can be + one of: ``xy``, ``xz``, ``yz``. + + *Default*: ``xy`` + + :pixels: + Specifies the number of pixes to be used along each of the basis directions for + ``slice`` plots. Should be 2 integers separated by spaces. + + .. warning:: The ``pixels`` input determines the output file size. For the PPM + format, 10 million pixels will result in a file just under 30 MB in + size. + + .. warning:: If the aspect ratio defined in ``pixels`` does not match the aspect + ratio defined in ``width`` the plot may appear stretched or squeezed. + + .. warning:: Geometry features along a basis direction smaller than ``width``/``pixels`` + along that basis direction may not appear in the plot. + + *Default*: None - Required entry for ``slice`` plots + + :background: + Specifies the RGB color of the regions where no OpenMC cell can be found. Should + be 3 integers deparated by spaces. + + *Default*: 0 0 0 (white) + + :col_spec: + Any number of this optional tag may be included in each ``plot`` element, which can + override the default random colors for cells or materials. Each ``col_spec`` + element must contain ``id`` and ``rgb`` sub-elements. + + :id: + Specifies the cell or material unique id for the color specification. + + :rgb: + Specifies the custom color for the cell or material. Should be 3 intergers separated + by spaces. + + *Default*: None + + :mask: + The special ``mask`` sub-element allows for the selective plotting of *only* + user-specified cells or materials. Only one ``mask`` element is allowed per ``plot`` + element, and it must contain as atributes or sub-elements a background masking color and + a list of cells or materials to plot: + + :components: + List of unique ``id`` numbers of the cells or materials to plot. Should be any number + of integers separated by spaces. + + :background: + Color to apply to all cells or materials not in the ``components`` list of cells or + materials to plot. This overrides any ``col_spec`` color specifications. + + *Default*: None diff --git a/examples/basic/materials.xml b/examples/basic/materials.xml index 88e20c69ef..e50a423112 100644 --- a/examples/basic/materials.xml +++ b/examples/basic/materials.xml @@ -1,17 +1,18 @@ - + 70c + - + - - - + + + diff --git a/examples/basic/settings.xml b/examples/basic/settings.xml index a9a83dfe54..73dd3cfbe0 100644 --- a/examples/basic/settings.xml +++ b/examples/basic/settings.xml @@ -1,14 +1,9 @@ - - - /home/paulromano/openmc/cross_sections_serpent.xml - - - 15 + 15 5 10000 diff --git a/examples/basic/tallies.xml b/examples/basic/tallies.xml new file mode 100644 index 0000000000..a434deae11 --- /dev/null +++ b/examples/basic/tallies.xml @@ -0,0 +1,28 @@ + + + + + + 100 + + total scatter nu-scatter absorption fission nu-fission + + + + + 100 + 0 20.0 + + total scatter nu-scatter absorption fission nu-fission + + + + + 100 + 0 20.0 + 0 20.0 + + scatter nu-scatter nu-fission + + + diff --git a/examples/lattice/nested/geometry.xml b/examples/lattice/nested/geometry.xml index 4ce62280a2..f223cf3bb4 100644 --- a/examples/lattice/nested/geometry.xml +++ b/examples/lattice/nested/geometry.xml @@ -14,7 +14,7 @@ rectangular 2 2 - -1.0 -1.0 + -1.0 -1.0 1.0 1.0 1 2 @@ -26,7 +26,7 @@ rectangular 2 2 - -2.0 -2.0 + -2.0 -2.0 2.0 2.0 5 5 diff --git a/examples/lattice/nested/plot.xml b/examples/lattice/nested/plot.xml deleted file mode 100644 index d6a0ca3cdf..0000000000 --- a/examples/lattice/nested/plot.xml +++ /dev/null @@ -1,6 +0,0 @@ - - - 0. 0. 0. - 3.999 3.999 - 0.01 - diff --git a/examples/lattice/nested/plots.xml b/examples/lattice/nested/plots.xml new file mode 100644 index 0000000000..a8c6fbb7c5 --- /dev/null +++ b/examples/lattice/nested/plots.xml @@ -0,0 +1,10 @@ + + + + + 0. 0. 0. + 4.0 4.0 + 400 400 + + + diff --git a/examples/lattice/nested/settings.xml b/examples/lattice/nested/settings.xml index aa0b98ee4a..2c440b3049 100644 --- a/examples/lattice/nested/settings.xml +++ b/examples/lattice/nested/settings.xml @@ -3,7 +3,7 @@ - 20 + 20 10 10000 diff --git a/examples/lattice/nested/tallies.xml b/examples/lattice/nested/tallies.xml index a389dc905c..271e81fe79 100644 --- a/examples/lattice/nested/tallies.xml +++ b/examples/lattice/nested/tallies.xml @@ -4,7 +4,7 @@ rectangular 4 4 - -2.0 -2.0 + -2.0 -2.0 1.0 1.0 @@ -13,7 +13,7 @@ 1 - total + total - \ No newline at end of file + diff --git a/examples/lattice/simple/geometry.xml b/examples/lattice/simple/geometry.xml index 4a77229392..b5bad90b9c 100644 --- a/examples/lattice/simple/geometry.xml +++ b/examples/lattice/simple/geometry.xml @@ -12,7 +12,7 @@ rectangular 4 4 - -2.0 -2.0 + -2.0 -2.0 1.0 1.0 1 2 1 2 diff --git a/examples/lattice/simple/plot.xml b/examples/lattice/simple/plot.xml deleted file mode 100644 index d6a0ca3cdf..0000000000 --- a/examples/lattice/simple/plot.xml +++ /dev/null @@ -1,6 +0,0 @@ - - - 0. 0. 0. - 3.999 3.999 - 0.01 - diff --git a/examples/lattice/simple/plots.xml b/examples/lattice/simple/plots.xml new file mode 100644 index 0000000000..a8c6fbb7c5 --- /dev/null +++ b/examples/lattice/simple/plots.xml @@ -0,0 +1,10 @@ + + + + + 0. 0. 0. + 4.0 4.0 + 400 400 + + + diff --git a/examples/lattice/simple/settings.xml b/examples/lattice/simple/settings.xml index aa0b98ee4a..2c440b3049 100644 --- a/examples/lattice/simple/settings.xml +++ b/examples/lattice/simple/settings.xml @@ -3,7 +3,7 @@ - 20 + 20 10 10000 diff --git a/examples/lattice/simple/tallies.xml b/examples/lattice/simple/tallies.xml index a389dc905c..271e81fe79 100644 --- a/examples/lattice/simple/tallies.xml +++ b/examples/lattice/simple/tallies.xml @@ -4,7 +4,7 @@ rectangular 4 4 - -2.0 -2.0 + -2.0 -2.0 1.0 1.0 @@ -13,7 +13,7 @@ 1 - total + total - \ No newline at end of file + diff --git a/examples/reflective/materials.xml b/examples/reflective/materials.xml index eb55455f7a..3ae574ff3e 100644 --- a/examples/reflective/materials.xml +++ b/examples/reflective/materials.xml @@ -1,10 +1,11 @@ - + 70c + - + diff --git a/examples/reflective/settings.xml b/examples/reflective/settings.xml index 3dd2f6e7d3..a3eefb21e1 100644 --- a/examples/reflective/settings.xml +++ b/examples/reflective/settings.xml @@ -1,14 +1,9 @@ - - - /home/paulromano/openmc/cross_sections_serpent.xml - - - 500 + 500 10 10000 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 59a0307533..51a6f6d081 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -86,6 +86,7 @@ energy_grid.o: output.o finalize.o: global.o finalize.o: hdf5_interface.o finalize.o: output.o +finalize.o: source.o finalize.o: tally.o finalize.o: timing.o @@ -145,6 +146,7 @@ initialize.o: random_lcg.o initialize.o: source.o initialize.o: string.o initialize.o: tally.o +initialize.o: tally_header.o initialize.o: timing.o input_xml.o: cmfd_input.o @@ -194,12 +196,13 @@ main.o: intercycle.o main.o: output.o main.o: physics.o main.o: plot.o -main.o: random_lcg.o main.o: source.o main.o: string.o main.o: tally.o main.o: timing.o +mesh.o: constants.o +mesh.o: global.o mesh.o: mesh_header.o mesh.o: particle_header.o @@ -226,6 +229,7 @@ physics.o: geometry.o physics.o: geometry_header.o physics.o: global.o physics.o: interpolation.o +physics.o: material_header.o physics.o: output.o physics.o: particle_header.o physics.o: random_lcg.o @@ -233,6 +237,8 @@ physics.o: search.o physics.o: string.o physics.o: tally.o +plot_header.o: constants.o + plot.o: constants.o plot.o: error.o plot.o: geometry.o diff --git a/src/Makefile b/src/Makefile index bb25277980..a5f899057b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -23,6 +23,12 @@ OPTIMIZE = no USE_MPI = no USE_HDF5 = yes +#=============================================================================== +# Add git SHA-1 hash +#=============================================================================== + +GIT_SHA1 = $(shell git log -1 | head -n 1 | awk '{print $$2}') + #=============================================================================== # GNU Fortran compiler options #=============================================================================== @@ -233,7 +239,7 @@ neat: .PHONY: all xml-fortran clean neat distclean %.o: %.F90 - $(F90) $(F90FLAGS) -Ixml-fortran -Ixml-fortran/templates -c $< + $(F90) $(F90FLAGS) -DGIT_SHA1="\"$(GIT_SHA1)\"" -Ixml-fortran -Ixml-fortran/templates -c $< #=============================================================================== # Dependencies diff --git a/src/ace.F90 b/src/ace.F90 index 0406c7dbca..af256c9e21 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -22,8 +22,6 @@ module ace real(8), allocatable :: XSS(:) ! Cross section data integer :: XSS_index ! current index in XSS data - type(DictionaryCI), pointer :: already_read => null() - private :: NXS private :: JXS private :: XSS @@ -47,6 +45,7 @@ contains type(Material), pointer :: mat => null() type(Nuclide), pointer :: nuc => null() type(SAB_Table), pointer :: sab => null() + type(DictionaryCI), pointer :: already_read => null() ! ========================================================================== ! COUNT NUMBER OF TABLES AND CREATE DICTIONARIES @@ -154,6 +153,8 @@ contains name = xs_listings(index_list) % name alias = xs_listings(index_list) % alias + ! Read the ACE table into the appropriate entry on the nuclides + ! array call read_ace_table(index_nuclides, index_list) ! Print out information on table to cross_sections.out file @@ -172,6 +173,8 @@ contains index_list = dict_get_key(xs_listing_dict, name) index_sab = dict_get_key(sab_dict, name) + ! Read the ACE table into the appropriate entry on the sab_tables + ! array call read_ace_table(index_sab, index_list) call dict_add_key(already_read, name, 0) diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 0b42aa155a..233ff611d4 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -9,6 +9,11 @@ module bank_header !=============================================================================== type Bank + ! The 'sequence' attribute is used here to ensure that the data listed + ! appears in the given order. This is important for MPI purposes when bank + ! sites are sent from one processor to another. + sequence + integer(8) :: id ! Unique ID real(8) :: xyz(3) ! location of bank particle real(8) :: uvw(3) ! diretional cosines diff --git a/src/constants.F90 b/src/constants.F90 index 89490eb183..90828d3019 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -7,7 +7,7 @@ module constants integer, parameter :: VERSION_MAJOR = 0 integer, parameter :: VERSION_MINOR = 4 - integer, parameter :: VERSION_RELEASE = 0 + integer, parameter :: VERSION_RELEASE = 1 ! ============================================================================ ! ADJUSTABLE PARAMETERS @@ -243,7 +243,7 @@ module constants EVENT_FISSION = 3 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 15 + integer, parameter :: N_SCORE_TYPES = 16 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -252,14 +252,15 @@ module constants SCORE_SCATTER_1 = -5, & ! first scattering moment SCORE_SCATTER_2 = -6, & ! second scattering moment SCORE_SCATTER_3 = -7, & ! third scattering moment - SCORE_N_1N = -8, & ! (n,1n) rate - SCORE_N_2N = -9, & ! (n,2n) rate - SCORE_N_3N = -10, & ! (n,3n) rate - SCORE_N_4N = -11, & ! (n,4n) rate - SCORE_ABSORPTION = -12, & ! absorption rate - SCORE_FISSION = -13, & ! fission rate - SCORE_NU_FISSION = -14, & ! neutron production rate - SCORE_CURRENT = -15 ! partial current + SCORE_DIFFUSION = -8, & ! diffusion coefficient + SCORE_N_1N = -9, & ! (n,1n) rate + SCORE_N_2N = -10, & ! (n,2n) rate + SCORE_N_3N = -11, & ! (n,3n) rate + SCORE_N_4N = -12, & ! (n,4n) rate + SCORE_ABSORPTION = -13, & ! absorption rate + SCORE_FISSION = -14, & ! fission rate + SCORE_NU_FISSION = -15, & ! neutron production rate + SCORE_CURRENT = -16 ! partial current ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 @@ -312,11 +313,17 @@ module constants integer, parameter :: ERROR_INT = -huge(0) real(8), parameter :: ERROR_REAL = -huge(0.0_8) * 0.917826354_8 + ! Energy grid methods + integer, parameter :: & + GRID_NUCLIDE = 1, & ! non-unionized energy grid (MCNP) + GRID_UNION = 2, & ! union grid with pointers + GRID_LETHARGY = 3 ! lethargy mapping (MC21) + ! Source types integer, parameter :: & SRC_BOX = 1, & ! Source in a rectangular prism - SRC_CELL = 2, & ! Source in a cell - SRC_SURFACE = 3 ! Source on a surface + SRC_POINT = 2, & ! Source at a single point + SRC_FILE = 3 ! Source from a file integer, parameter :: & PROB_SOURCE = 1, & ! External source problem @@ -327,5 +334,6 @@ module constants integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file + integer, parameter :: UNIT_SOURCE = 15 ! unit # for writing source file end module constants diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 966ff1cd94..aa155eefe7 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -37,7 +37,7 @@ contains mat => materials(p % material) ! Find energy index on unionized grid - call find_energy_index() + if (grid_method == GRID_UNION) call find_energy_index() ! Check if there's an S(a,b) table for this material if (mat % has_sab_table) then @@ -107,11 +107,28 @@ contains ! Set pointer to nuclide nuc => nuclides(index_nuclide) - ! TODO: If not using unionized energy grid, we need to find the index on the - ! nuclide energy grid using lethargy mapping or whatever other technique + ! Determine index on nuclide energy grid + select case (grid_method) + case (GRID_UNION) + ! If we're using the unionized grid with pointers, finding the index on + ! the nuclide energy grid is as simple as looking up the pointer - ! get index on nuclide energy grid - IE = nuc % grid_index(p % IE) + IE = nuc % grid_index(p % IE) + + case (GRID_NUCLIDE) + ! If we're not using the unionized grid, we have to do a binary search on + ! the nuclide energy grid in order to determine which points to + ! interpolate between + + if (p % E < nuc % energy(1)) then + IE = 1 + elseif (p % E > nuc % energy(nuc % n_grid)) then + IE = nuc % n_grid - 1 + else + IE = binary_search(nuc % energy, nuc % n_grid, p % E) + end if + + end select ! check for rare case where two energy points are the same if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1 @@ -171,7 +188,7 @@ contains end if ! Set last evaluated energy - micro_xs(index_nuclide) % last_E = p % E + if (index_sab == 0) micro_xs(index_nuclide) % last_E = p % E end subroutine calculate_nuclide_xs diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 3fb1d940d2..4828338a2e 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -37,7 +37,7 @@ contains ! Create dictionary for keeping track of cross sections already added call dict_create(already_added) - ! Loop over all files + ! Loop over all materials do i = 1, n_materials mat => materials(i) diff --git a/src/error.F90 b/src/error.F90 index 4f4c53d78a..59cc6d80e0 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -76,9 +76,10 @@ contains write(ERROR_UNIT,*) ! Write information on current cycle and particle - if (current_cycle > 0) then - write(ERROR_UNIT,'(1X,A,I11) ') 'Cycle: ', current_cycle - write(ERROR_UNIT,'(1X,A,I11)') 'Particle: ', p % id + if (current_batch > 0) then + write(ERROR_UNIT,'(1X,A,I12) ') 'Batch: ', current_batch + write(ERROR_UNIT,'(1X,A,I12) ') 'Generation:', current_gen + write(ERROR_UNIT,'(1X,A,I12)') 'Particle: ', p % id write(ERROR_UNIT,*) end if diff --git a/src/finalize.F90 b/src/finalize.F90 index c4dcaf0623..e239c6db9f 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -2,6 +2,7 @@ module finalize use global use output, only: print_runtime + use source, only: write_source_binary use tally, only: write_tallies, tally_statistics use timing, only: timer_start, timer_stop @@ -23,10 +24,13 @@ contains ! Start finalization timer call timer_start(time_finalize) - ! Calculate statistics for tallies and write to tallies.out if (.not. plotting) then + ! Calculate statistics for tallies and write to tallies.out call tally_statistics() if (master) call write_tallies() + + ! Write out binary source + if (write_source) call write_source_binary() end if ! stop timers and show timing statistics diff --git a/src/geometry.F90 b/src/geometry.F90 index 3457719556..0abda64c13 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -136,6 +136,12 @@ contains ! Set cell on this level p % coord % cell = index_cell + ! Show cell information on trace + if (verbosity >= 10 .or. trace) then + message = " Entering cell " // trim(to_str(c % id)) + call write_message() + end if + if (c % type == CELL_NORMAL) then ! ================================================================= ! AT LOWEST UNIVERSE, TERMINATE SEARCH diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index aa1fcca5f6..fbcf0353a9 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -64,7 +64,6 @@ module geometry_header & surfaces(:) ! List of surfaces bounding cell -- note that ! parentheses, union, etc operators will be listed ! here too - integer, allocatable :: rgb(:) ! plotting color end type Cell ! array index of universe 0 diff --git a/src/global.F90 b/src/global.F90 index 5ed20ccd12..0068021c64 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -83,7 +83,8 @@ module global type(DictionaryCI), pointer :: xs_listing_dict => null() ! Unionized energy grid - integer :: n_grid ! number of points on unionized grid + integer :: grid_method ! how to treat the energy grid + integer :: n_grid ! number of points on unionized grid real(8), allocatable :: e_grid(:) ! energies on unionized grid ! Unreoslved resonance probablity tables @@ -126,11 +127,13 @@ module global ! ============================================================================ ! CRITICALITY SIMULATION VARIABLES - integer(8) :: n_particles = 10000 ! # of particles per cycle - integer :: n_cycles = 500 ! # of cycles - integer :: n_inactive = 50 ! # of inactive cycles - integer :: n_active ! # of active cycles - integer :: current_cycle = 0 ! current cycle + integer(8) :: n_particles ! # of particles per generation + integer :: n_batches ! # of batches + integer :: n_inactive ! # of inactive batches + integer :: n_active ! # of active batches + integer :: gen_per_batch = 1 ! # of generations per batch + integer :: current_batch = 0 ! current batch + integer :: current_gen = 0 ! current generation ! Random Number seed integer(8) :: seed = 1_8 @@ -147,16 +150,19 @@ module global integer(8) :: work ! number of particles per processor integer(8) :: maxwork ! maximum number of particles per processor - ! cycle keff + ! single-genreation keff real(8) :: keff = ONE real(8) :: keff_std ! Shannon entropy logical :: entropy_on = .false. - real(8) :: entropy ! value of shannon entropy - real(8), allocatable :: entropy_p(:) ! fraction of source sites in each cell + real(8) :: entropy ! value of shannon entropy + real(8), allocatable :: entropy_p(:,:,:) ! % of source sites in each cell type(StructuredMesh), pointer :: entropy_mesh + ! Write source at end of simulation + logical :: write_source = .false. + ! ============================================================================ ! PARALLEL PROCESSING VARIABLES @@ -212,6 +218,9 @@ module global ! Message used in message/warning/fatal_error character(MAX_LINE_LEN) :: message + ! Random number seed + integer(8) :: seed = 1_8 + ! Problem type integer :: problem_type = PROB_CRITICALITY @@ -221,7 +230,8 @@ module global ! Trace for single particle logical :: trace - integer :: trace_cycle + integer :: trace_batch + integer :: trace_gen integer(8) :: trace_particle ! ============================================================================ diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 4c5091f05f..56e921822e 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -101,19 +101,22 @@ contains call hdf5_make_double(hdf5_output_file, "n_particles", real(n_particles,8)) ! Use H5LT interface to write n_cycles, n_inactive, and n_active - call hdf5_make_integer(hdf5_output_file, "n_cycles", n_cycles) + call hdf5_make_integer(hdf5_output_file, "n_batches", n_batches) call hdf5_make_integer(hdf5_output_file, "n_inactive", n_inactive) call hdf5_make_integer(hdf5_output_file, "n_active", n_active) + call hdf5_make_integer(hdf5_output_file, "gen_per_batch", gen_per_batch) ! Add description of each variable call h5ltset_attribute_string_f(hdf5_output_file, "n_particles", & "description", "Number of particles per cycle", hdf5_err) - call h5ltset_attribute_string_f(hdf5_output_file, "n_cycles", & - "description", "Total number of cycles", hdf5_err) + call h5ltset_attribute_string_f(hdf5_output_file, "n_batches", & + "description", "Total number of batches", hdf5_err) call h5ltset_attribute_string_f(hdf5_output_file, "n_inactive", & "description", "Number of inactive cycles", hdf5_err) call h5ltset_attribute_string_f(hdf5_output_file, "n_active", & "description", "Number of active cycles", hdf5_err) + call h5ltset_attribute_string_f(hdf5_output_file, "gen_per_batch", & + "description", "Number of generations per batch", hdf5_err) end if call hdf5_write_geometry() @@ -812,7 +815,7 @@ contains call h5ltset_attribute_string_f(timing_group, "time_transport", & "description", "Time in transport only (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_intercycle", & - "description", "Total time between cycles (s)", hdf5_err) + "description", "Total time between generations (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_tallies", & "description", "Time between cycles accumulating tallies (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_sample", & @@ -820,16 +823,16 @@ contains call h5ltset_attribute_string_f(timing_group, "time_sendrecv", & "description", "Time between cycles SEND/RECVing source sites (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_inactive", & - "description", "Total time in inactive cycles (s)", hdf5_err) + "description", "Total time in inactive batches (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_active", & - "description", "Total time in active cycles (s)", hdf5_err) + "description", "Total time in active batches (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_finalize", & "description", "Total time for finalization (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_total", & "description", "Total time elapsed (s)", hdf5_err) ! Write calculation rate - total_particles = n_particles * n_cycles + total_particles = n_particles * n_batches * gen_per_batch speed = real(total_particles) / (time_inactive % elapsed + & time_active % elapsed) call hdf5_make_double(timing_group, "neutrons_per_second", speed) diff --git a/src/initialize.F90 b/src/initialize.F90 index db2ca4b585..4c007f860d 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -3,23 +3,24 @@ module initialize use ace, only: read_xs use bank_header, only: Bank use constants - use datatypes, only: dict_create, dict_add_key, dict_get_key, & - dict_has_key, dict_keys + use datatypes, only: dict_create, dict_get_key, dict_has_key, & + dict_keys use datatypes_header, only: ListKeyValueII, DictionaryII use energy_grid, only: unionized_grid use error, only: fatal_error use geometry, only: neighbor_lists - use geometry_header, only: Cell, Surface, Universe, Lattice, BASE_UNIVERSE + use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE use global use input_xml, only: read_input_xml, read_cross_sections_xml, & - cells_in_univ_dict + cells_in_univ_dict, read_plots_xml use output, only: title, header, print_summary, print_geometry, & print_plot, create_summary_file, & create_xs_summary_file use random_lcg, only: initialize_prng use source, only: initialize_source use string, only: to_str, starts_with, ends_with, lower_case - use tally, only: create_tally_map, TallyObject + use tally, only: create_tally_map + use tally_header, only: TallyObject use timing, only: timer_start, timer_stop #ifdef MPI @@ -33,8 +34,6 @@ module initialize implicit none - type(DictionaryII), pointer :: build_dict => null() - contains !=============================================================================== @@ -81,6 +80,10 @@ contains ! Initialize random number generator call initialize_prng() + ! Read plots.xml if it exists -- this has to be done separate from the other + ! XML files because we need the PRNG to be initialized first + if (plotting) call read_plots_xml() + ! Set up universe structures call prepare_universes() @@ -106,9 +109,11 @@ contains call timer_stop(time_read_xs) ! Construct unionized energy grid from cross-sections - call timer_start(time_unionize) - call unionized_grid() - call timer_stop(time_unionize) + if (grid_method == GRID_UNION) then + call timer_start(time_unionize) + call unionized_grid() + call timer_stop(time_unionize) + end if ! Create tally map call create_tally_map() @@ -282,9 +287,6 @@ contains ! Create special dictionary used in input_xml call dict_create(cells_in_univ_dict) - - ! Create special dictionary for building universes - call dict_create(build_dict) end subroutine create_dictionaries diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c7adcae12f..146c04ca5b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2,8 +2,7 @@ module input_xml use cmfd_input, only: read_cmfd_xml use constants - use datatypes, only: dict_create, dict_add_key, dict_has_key, & - dict_get_key + use datatypes, only: dict_add_key, dict_has_key, dict_get_key use error, only: fatal_error, warning use geometry_header, only: Cell, Surface, Lattice use global @@ -34,7 +33,6 @@ contains call read_materials_xml() call read_tallies_xml() call read_cmfd_xml() - if (plotting) call read_plots_xml() end subroutine read_input_xml @@ -69,6 +67,10 @@ contains ! Initialize XML scalar variables cross_sections_ = "" verbosity_ = 0 + energy_grid_ = "union" + seed_ = 0_8 + write_source_ = "" + source_ % type = "" ! Parse settings.xml file call read_xml_file_settings_t(filename) @@ -93,7 +95,7 @@ contains end if ! Criticality information - if (criticality % cycles > 0) then + if (criticality % batches > 0) then problem_type = PROB_CRITICALITY ! Check number of particles @@ -104,41 +106,71 @@ contains n_particles = str_to_int(criticality % particles) ! Copy cycle information - n_cycles = criticality % cycles - n_inactive = criticality % inactive - n_active = n_cycles - n_inactive + n_batches = criticality % batches + n_inactive = criticality % inactive + n_active = n_batches - n_inactive + gen_per_batch = criticality % generations_per_batch + else + message = "Need to specify number of batches with tag." + call fatal_error() end if +<<<<<<< HEAD:src/input_xml.F90 ! Random number seed if (dble(criticality % seed) > 0.0_8) seed = dble(criticality % seed) +======= + ! Copy random number seed if specified + if (seed_ > 0) seed = seed_ + + ! Energy grid methods + select case (energy_grid_) + case ('nuclide') + grid_method = GRID_NUCLIDE + case ('union') + grid_method = GRID_UNION + case ('lethargy') + message = "Lethargy mapped energy grid not yet supported." + call fatal_error() + case default + message = "Unknown energy grid method: " // energy_grid_ + call fatal_error() + end select +>>>>>>> master:src/input_xml.F90 ! Verbosity if (verbosity_ > 0) verbosity = verbosity_ - if (associated(source_ % coeffs)) then + if (len(source_ % type) > 0) then ! Determine external source type type = source_ % type call lower_case(type) - select case (trim(type)) + select case (type) case ('box') external_source % type = SRC_BOX coeffs_reqd = 6 + case ('point') + external_source % type = SRC_POINT + coeffs_reqd = 3 + case ('file') + external_source % type = SRC_FILE case default message = "Invalid source type: " // trim(source_ % type) call fatal_error() end select ! Coefficients for external surface - n = size(source_ % coeffs) - if (n < coeffs_reqd) then - message = "Not enough coefficients specified for external source." - call fatal_error() - elseif (n > coeffs_reqd) then - message = "Too many coefficients specified for external source." - call fatal_error() - else - allocate(external_source % values(n)) - external_source % values = source_ % coeffs + if (type /= 'file') then + n = size(source_ % coeffs) + if (n < coeffs_reqd) then + message = "Not enough coefficients specified for external source." + call fatal_error() + elseif (n > coeffs_reqd) then + message = "Too many coefficients specified for external source." + call fatal_error() + else + allocate(external_source % values(n)) + external_source % values = source_ % coeffs + end if end if end if @@ -156,8 +188,9 @@ contains ! Particle trace if (associated(trace_)) then - trace_cycle = trace_(1) - trace_particle = trace_(2) + trace_batch = trace_(1) + trace_gen = trace_(2) + trace_particle = trace_(3) end if ! Shannon Entropy mesh @@ -210,6 +243,9 @@ contains entropy_on = .true. end if + ! Check if the user has specified to write binary source file + if (trim(write_source_) == 'on') write_source = .true. + end subroutine read_settings_xml !=============================================================================== @@ -266,14 +302,6 @@ contains c % material = cell_(i) % material c % fill = cell_(i) % fill - ! Set plot color - if (plotting) then - allocate(c % rgb(3)) - c % rgb(1) = prn()*255 - c % rgb(2) = prn()*255 - c % rgb(3) = prn()*255 - end if - ! Check to make sure that either material or fill was specified if (c % material == 0 .and. c % fill == 0) then message = "Neither material nor fill was specified for cell " // & @@ -505,6 +533,7 @@ contains integer :: n ! number of nuclides real(8) :: val ! value entered for density logical :: file_exists ! does materials.xml exist? + logical :: sum_density ! density is taken to be sum of nuclide densities character(3) :: default_xs ! default xs identifier (e.g. 70c) character(12) :: name ! name of nuclide character(MAX_WORD_LEN) :: units ! units on density @@ -544,40 +573,43 @@ contains ! Copy material id m % id = material_(i) % id - ! Set plot color - if (plotting) then - allocate(m % rgb(3)) - m % rgb(1) = prn()*255 - m % rgb(2) = prn()*255 - m % rgb(3) = prn()*255 - end if - - ! Copy density -- the default value for the units is given in the - ! material_t.xml file and doesn't need to be specified here, hence case - ! default results in an error. - val = material_(i) % density % value - if (val <= ZERO) then - message = "Need to specify a positive density on material " // & - trim(to_str(m % id)) // "." - call fatal_error() - end if - + ! Copy value and units + val = material_(i) % density % value units = material_(i) % density % units - call lower_case(units) - select case(trim(units)) - case ('g/cc', 'g/cm3') - m % density = -val - case ('kg/m3') - m % density = -0.001 * val - case ('atom/b-cm') - m % density = val - case ('atom/cm3', 'atom/cc') - m % density = 1.0e-24 * val - case default - message = "Unkwown units '" // trim(material_(i) % density % units) & - // "' specified on material " // trim(to_str(m % id)) - call fatal_error() - end select + + if (units == 'sum') then + ! If the user gave the units as 'sum', then the total density of the + ! material is taken to be the sum of the atom fractions listed on the + ! nuclides + + sum_density = .true. + + else + ! Check for erroneous density + sum_density = .false. + if (val <= ZERO) then + message = "Need to specify a positive density on material " // & + trim(to_str(m % id)) // "." + call fatal_error() + end if + + ! Adjust material density based on specified units + call lower_case(units) + select case(trim(units)) + case ('g/cc', 'g/cm3') + m % density = -val + case ('kg/m3') + m % density = -0.001 * val + case ('atom/b-cm') + m % density = val + case ('atom/cm3', 'atom/cc') + m % density = 1.0e-24 * val + case default + message = "Unkwown units '" // trim(material_(i) % density % units) & + // "' specified on material " // trim(to_str(m % id)) + call fatal_error() + end select + end if ! Check to ensure material has at least one nuclide if (.not. associated(material_(i) % nuclides)) then @@ -652,6 +684,9 @@ contains end if end do + ! Determine density if it is a sum value + if (sum_density) m % density = sum(m % atom_percent) + ! Add material to dictionary call dict_add_key(material_dict, m % id, i) @@ -986,6 +1021,11 @@ contains case ('scatter-3') t % score_bins(j) % scalar = SCORE_SCATTER_3 + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + case ('diffusion') + t % score_bins(j) % scalar = SCORE_DIFFUSION + ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('n1n') @@ -1116,10 +1156,11 @@ contains use xml_data_plots_t - integer i + integer i, j + integer n_cols, col_id logical :: file_exists ! does plots.xml file exist? character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml - type(Plot), pointer :: pl => null() + type(Plot), pointer :: pl => null() ! Check if plots.xml exists filename = trim(path_input) // "plots.xml" @@ -1145,26 +1186,33 @@ contains ! Copy data into plots pl % id = plot_(i) % id - pl % aspect = plot_(i) % aspect + ! Set output file path + pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & + "_" // trim(plot_(i) % filename) // ".ppm" + + ! Copy plot pixel size if (size(plot_(i) % pixels) == 2) then pl % pixels = plot_(i) % pixels else - message = " must be length 2 in plot " // to_str(i) + message = " must be length 2 in plot " // to_str(pl % id) call fatal_error() end if - select case (plot_(i) % color) - case ("cell") - pl % color = PLOT_COLOR_CELLS - case ("mat") - pl % color = PLOT_COLOR_MATS - case default - message = "Unsupported plot color '" // plot_(i) % color & - // "' in plot " // trim(to_str(i)) - call fatal_error() - end select + ! Copy plot background color + if (associated(plot_(i) % background)) then + if (size(plot_(i) % background) == 3) then + pl % not_found % rgb = plot_(i) % background + else + message = "Bad background RGB " & + // "in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + else + pl % not_found % rgb = (/ 255, 255, 255 /) + end if + ! Copy plot type select case (plot_(i) % type) case ("slice") pl % type = PLOT_TYPE_SLICE @@ -1172,10 +1220,11 @@ contains ! pl % type = PLOT_TYPE_POINTS case default message = "Unsupported plot type '" // plot_(i) % type & - // "' in plot " // trim(to_str(i)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select + ! Copy plot basis select case (plot_(i) % basis) case ("xy") pl % basis = PLOT_BASIS_XY @@ -1185,7 +1234,7 @@ contains pl % basis = PLOT_BASIS_YZ case default message = "Unsupported plot basis '" // plot_(i) % basis & - // "' in plot " // trim(to_str(i)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -1194,7 +1243,7 @@ contains pl % origin = plot_(i) % origin else message = "Origin must be length 3 " & - // "in plot " // trim(to_str(i)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -1206,10 +1255,88 @@ contains pl % width(2) = plot_(i) % width(2) else message = "Bad plot width " & - // "in plot " // trim(to_str(i)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if + ! Copy plot color type and initialize all colors randomly + select case (plot_(i) % color) + case ("cell") + + pl % color_by = PLOT_COLOR_CELLS + allocate(pl % colors(n_cells)) + do j = 1, n_cells + pl % colors(j) % rgb(1) = prn()*255 + pl % colors(j) % rgb(2) = prn()*255 + pl % colors(j) % rgb(3) = prn()*255 + end do + + case ("mat", "material") + + pl % color_by = PLOT_COLOR_MATS + allocate(pl % colors(n_materials)) + do j = 1, n_materials + pl % colors(j) % rgb(1) = prn()*255 + pl % colors(j) % rgb(2) = prn()*255 + pl % colors(j) % rgb(3) = prn()*255 + end do + + case default + message = "Unsupported plot color type '" // plot_(i) % color & + // "' in plot " // trim(to_str(pl % id)) + call fatal_error() + end select + + ! Copy user specified colors + if (associated(plot_(i) % col_spec_)) then + n_cols = size(plot_(i) % col_spec_) + do j = 1, n_cols + if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then + message = "Bad RGB " & + // "in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + col_id = plot_(i) % col_spec_(j) % id + + if (pl % color_by == PLOT_COLOR_CELLS) then + + if (dict_has_key(cell_dict, col_id)) then + pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb + else + message = "Could not find cell " // trim(to_str(col_id)) // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + else if (pl % color_by == PLOT_COLOR_MATS) then + + if (dict_has_key(material_dict, col_id)) then + pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb + else + message = "Could not find material " // trim(to_str(col_id)) // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + end if + end do + end if + + ! Alter colors based on mask information + if (associated(plot_(i) % mask_)) then + if (size(plot_(i) % mask_) > 1) then + message = "Mutliple masks" // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + else if (.not. size(plot_(i) % mask_) == 0) then + do j=1,size(pl % colors) + if (.not. any(j .eq. plot_(i) % mask_(1) % components)) then + pl % colors(j) % rgb = plot_(i) % mask_(1) % background + end if + end do + end if + end if end do diff --git a/src/intercycle.F90 b/src/intercycle.F90 index adfeaa0288..1961b84a69 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -4,13 +4,13 @@ module intercycle use error, only: fatal_error, warning use global - use mesh, only: get_mesh_bin + use mesh, only: count_fission_sites use mesh_header, only: StructuredMesh use output, only: write_message use random_lcg, only: prn, set_particle_seed, prn_skip use search, only: binary_search use string, only: to_str - use tally, only: accumulate_cycle_estimate + use tally, only: accumulate_batch_estimate use tally_header, only: TallyObject use timing, only: timer_start, timer_stop @@ -22,8 +22,8 @@ contains !=============================================================================== ! SYNCHRONIZE_BANK samples source sites from the fission sites that were -! accumulated during the cycle. This routine is what allows this Monte Carlo to -! scale to large numbers of processors where other codes cannot. +! accumulated during the generation. This routine is what allows this Monte +! Carlo to scale to large numbers of processors where other codes cannot. !=============================================================================== subroutine synchronize_bank() @@ -73,7 +73,7 @@ contains total = n_bank #endif - ! If there are not that many particles per cycle, it's possible that no + ! If there are not that many particles per generation, it's possible that no ! fission sites were created at all on a single processor. Rather than add ! extra logic to treat this circumstance, we really want to ensure the user ! runs enough particles to avoid this in the first place. @@ -83,12 +83,12 @@ contains call fatal_error() end if - ! Make sure all processors start at the same point for random sampling by - ! using the current_cycle as the starting seed. Then skip ahead in the - ! sequence using the starting index in the 'global' fission bank for each - ! processor Skip ahead however many random numbers are needed + ! Make sure all processors start at the same point for random sampling. Then + ! skip ahead in the sequence using the starting index in the 'global' + ! fission bank for each processor. - call set_particle_seed(int(current_cycle,8)) + call set_particle_seed(int((current_batch - 1)*gen_per_batch + & + current_gen,8)) call prn_skip(start) ! Determine how many fission sites we need to sample from the source bank @@ -257,20 +257,22 @@ contains ! Send we initiated a series of asynchronous ISENDs and IRECVs, now we have ! to ensure that the data has actually been communicated before moving on to - ! the next cycle + ! the next generation call MPI_WAITALL(n_request, request, MPI_STATUSES_IGNORE, mpi_err) - ! Deallocate space for bank_position on the last cycle - if (current_cycle == n_cycles) deallocate(bank_position) + ! Deallocate space for bank_position on the very last generation + if (current_batch == n_batches .and. current_gen == gen_per_batch) & + deallocate(bank_position) #else source_bank = temp_sites(1:n_particles) #endif call timer_stop(time_ic_sendrecv) - ! Deallocate space for the temporary source bank on the last cycle - if (current_cycle == n_cycles) deallocate(temp_sites) + ! Deallocate space for the temporary source bank on the last generation + if (current_batch == n_batches .and. current_gen == gen_per_batch) & + deallocate(temp_sites) end subroutine synchronize_bank @@ -281,12 +283,10 @@ contains subroutine shannon_entropy() - integer :: i ! index for bank sites + integer :: i, j, k ! index for bank sites integer :: n ! # of boxes in each dimension - integer :: bin ! index in entropy_p - integer(8) :: total_bank ! total # of fission bank sites - integer, save :: n_box ! total # of boxes on mesh - logical :: outside_box ! were there sites outside entropy box? + real(8) :: total ! total weight of fission bank sites + logical :: sites_outside ! were there sites outside entropy box? type(StructuredMesh), pointer :: m => null() ! Get pointer to entropy mesh @@ -310,96 +310,61 @@ contains m % dimension = n end if - ! Determine total number of mesh boxes - n_box = product(m % dimension) - ! allocate and determine width allocate(m % width(3)) m % width = (m % upper_right - m % lower_left) / m % dimension ! allocate p - allocate(entropy_p(n_box)) + allocate(entropy_p(m % dimension(1), m % dimension(2), & + m % dimension(3))) end if - ! initialize p - entropy_p = ZERO - outside_box = .false. - - ! loop over fission sites and count how many are in each mesh box - FISSION_SITES: do i = 1, int(n_bank,4) - ! determine scoring bin for entropy mesh - call get_mesh_bin(m, fission_bank(i) % xyz, bin) - - ! if outside mesh, skip particle - if (bin == NO_BIN_FOUND) then - outside_box = .true. - cycle - end if - - ! add to appropriate mesh box - entropy_p(bin) = entropy_p(bin) + 1 - end do FISSION_SITES + ! count number of fission sites over mesh + call count_fission_sites(m, entropy_p, total, sites_outside) ! display warning message if there were sites outside entropy box - if (outside_box) then + if (sites_outside) then message = "Fission source site(s) outside of entropy box." call warning() end if -#ifdef MPI - ! collect values from all processors - if (master) then - call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) - else - call MPI_REDUCE(entropy_p, entropy_p, n_box, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) - end if - - ! determine total number of bank sites - call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & - MPI_COMM_WORLD, mpi_err) -#else - total_bank = n_bank -#endif - ! sum values to obtain shannon entropy if (master) then - ! Normalize to number of bank sites - entropy_p = entropy_p / total_bank + ! Normalize to total weight of bank sites + entropy_p = entropy_p / total entropy = 0 - do i = 1, n_box - if (entropy_p(i) > 0) then - entropy = entropy - entropy_p(i) * log(entropy_p(i))/log(2.0) - end if + do i = 1, m % dimension(1) + do j = 1, m % dimension(2) + do k = 1, m % dimension(3) + if (entropy_p(i,j,k) > ZERO) then + entropy = entropy - entropy_p(i,j,k) * & + log(entropy_p(i,j,k))/log(2.0) + end if + end do + end do end do end if end subroutine shannon_entropy !=============================================================================== -! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the -! mean and standard deviation of the mean for active cycles and displays them +! CALCULATE_KEFF calculates the single batch estimate of keff as well as the +! mean and standard deviation of the mean for active batches and displays them !=============================================================================== subroutine calculate_keff() - integer :: n ! active cycle number - real(8) :: k_cycle ! single cycle estimate of keff + integer :: n ! active batch number + real(8) :: k_batch ! single batch estimate of keff #ifdef MPI real(8) :: global_temp(N_GLOBAL_TALLIES) #endif - message = "Calculate cycle keff..." + message = "Calculate batch keff..." call write_message(8) - ! Since the creation of bank sites was originally weighted by the last - ! cycle keff, we need to multiply by that keff to get the current cycle's - ! value #ifdef MPI - global_tallies(K_ANALOG) % value = n_bank * keff - ! Copy global tallies into array to be reduced global_temp = global_tallies(:) % value @@ -416,59 +381,65 @@ contains ! Reset value on other processors global_tallies(:) % value = ZERO end if -#else - global_tallies(K_ANALOG) % value = n_bank * keff #endif ! Collect statistics and print output if (master) then - k_cycle = global_tallies(K_ANALOG) % value/n_particles + k_batch = global_tallies(K_ANALOG) % value/(n_particles*gen_per_batch) - if (current_cycle > n_inactive) then - ! Active cycle number - n = current_cycle - n_inactive + if (tallies_on) then + ! Active batch number + n = current_batch - n_inactive - ! Accumulate single cycle realizations of k - call accumulate_cycle_estimate(global_tallies) + ! Accumulate single batch realizations of k + call accumulate_batch_estimate(global_tallies) - ! Determine mean and standard deviation of mean + ! Determine sample mean of keff keff = global_tallies(K_ANALOG) % sum/n - keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - keff*keff)/n) - ! Display output for this cycle - if (current_cycle > n_inactive + 1) then + ! Display output for this batch + if (current_batch > n_inactive + 1) then + ! Determine standard deviation of the sample mean of keff + keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - & + keff*keff)/(n - 1)) + if (entropy_on) then if (cmfd_on) then - write(UNIT=OUTPUT_UNIT, FMT=104) current_cycle, k_cycle, & + write(UNIT=OUTPUT_UNIT, FMT=104) current_batch, k_batch, & entropy, keff, keff_std, cmfd%keff, cmfd%entropy else - write(UNIT=OUTPUT_UNIT, FMT=103) current_cycle, k_cycle, & + write(UNIT=OUTPUT_UNIT, FMT=103) current_batch, k_batch, & entropy, keff, keff_std end if else if (cmfd_on) then - write(UNIT=OUTPUT_UNIT, FMT=105) current_cycle, k_cycle, & + write(UNIT=OUTPUT_UNIT, FMT=105) current_batch, k_batch, & keff, keff_std, cmfd%keff else - write(UNIT=OUTPUT_UNIT, FMT=101) current_cycle, k_cycle, & + write(UNIT=OUTPUT_UNIT, FMT=101) current_batch, k_batch, & keff, keff_std end if end if else if (entropy_on) then - write(UNIT=OUTPUT_UNIT, FMT=102) current_cycle, k_cycle, entropy + write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy else - write(UNIT=OUTPUT_UNIT, FMT=100) current_cycle, k_cycle + write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch end if end if else - ! Display output for inactive cycle + ! Display output for inactive batch if (entropy_on) then - write(UNIT=OUTPUT_UNIT, FMT=102) current_cycle, k_cycle, entropy + write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy else - write(UNIT=OUTPUT_UNIT, FMT=100) current_cycle, k_cycle + write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch end if - keff = k_cycle + + ! Set keff + keff = k_batch + + ! Reset tally values + global_tallies(:) % value = ZERO end if end if diff --git a/src/main.F90 b/src/main.F90 index 5f75d57f92..58e911d6f5 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -9,10 +9,9 @@ program main use output, only: write_message, header use plotter, only: run_plot use physics, only: transport - use random_lcg, only: set_particle_seed use source, only: get_source_particle use string, only: to_str - use tally, only: synchronize_tallies + use tally, only: synchronize_tallies, add_to_score use timing, only: timer_start, timer_stop #ifdef MPI @@ -55,25 +54,26 @@ contains ! Display column titles if (entropy_on) then +<<<<<<< HEAD:src/main.F90 if (cmfd_on) then - message = " Cycle k(cycle) Entropy Average k CMFD k CMFD Ent." + message = " Cycle k(batch) Entropy Average k CMFD k CMFD Ent." call write_message(1) message = " ===== ======== ======= =================== ========== =========" call write_message(1) else - message = " Cycle k(cycle) Entropy Average k" + message = " Cycle k(batch) Entropy Average k" call write_message(1) message = " ===== ======== ======= ===================" call write_message(1) end if else if (cmfd_on) then - message = " Cycle k(cycle) Average k CMFD k" + message = " Cycle k(batch) Average k CMFD k" call write_message(1) message = " ===== ======== =================== ==========" call write_message(1) else - message = " Cycle k(cycle) Average k" + message = " Cycle k(batch) Average k" call write_message(1) message = " ===== ======== ===================" call write_message(1) @@ -81,39 +81,56 @@ contains end if ! ========================================================================== - ! LOOP OVER CYCLES - CYCLE_LOOP: do current_cycle = 1, n_cycles + ! LOOP OVER BATCHES + BATCH_LOOP: do current_batch = 1, n_batches - message = "Simulating cycle " // trim(to_str(current_cycle)) // "..." + message = "Simulating batch " // trim(to_str(current_batch)) // "..." call write_message(8) - ! Set all tallies to zero - n_bank = 0 - ! ======================================================================= - ! LOOP OVER HISTORIES + ! LOOP OVER GENERATIONS + GENERATION_LOOP: do current_gen = 1, gen_per_batch - ! Start timer for transport - call timer_start(time_transport) + ! Set all tallies to zero + n_bank = 0 - HISTORY_LOOP: do i = 1, work + ! ==================================================================== + ! LOOP OVER HISTORIES - ! grab source particle from bank - call get_source_particle(i) + ! Start timer for transport + call timer_start(time_transport) - ! transport particle - call transport() + HISTORY_LOOP: do i = 1, work - end do HISTORY_LOOP + ! grab source particle from bank + call get_source_particle(i) - ! Accumulate time for transport - call timer_stop(time_transport) + ! transport particle + call transport() - ! ======================================================================= - ! WRAP UP FISSION BANK AND COMPUTE TALLIES, KEFF, ETC + end do HISTORY_LOOP - ! Start timer for inter-cycle synchronization - call timer_start(time_intercycle) + ! Accumulate time for transport + call timer_stop(time_transport) + + ! ==================================================================== + ! WRAP UP FISSION BANK AND COMPUTE TALLIES, KEFF, ETC + + ! Start timer for inter-cycle synchronization + call timer_start(time_intercycle) + + ! Distribute fission bank across processors evenly + call synchronize_bank() + + ! Add to analog estimate of keff -- since the creation of bank sites + ! was originally weighted by the last cycle keff, we need to multiply + ! by that keff to get the current cycle's value + call add_to_score(global_tallies(K_ANALOG), n_bank * keff) + + ! Stop timer for inter-cycle synchronization + call timer_stop(time_intercycle) + + end do GENERATION_LOOP ! Collect tallies if (tallies_on) then @@ -125,30 +142,22 @@ contains ! Calculate shannon entropy if (entropy_on) call shannon_entropy() - ! Distribute fission bank across processors evenly - call synchronize_bank() - ! Collect results and statistics call calculate_keff() - ! print cycle information - ! run cmfd - if (current_cycle > n_inactive .and. cmfd_on) then + if (current_batch > n_inactive .and. cmfd_on) then call execute_cmfd() end if ! Turn tallies on once inactive cycles are complete - if (current_cycle == n_inactive) then + if (current_batch == n_inactive) then tallies_on = .true. call timer_stop(time_inactive) call timer_start(time_active) end if - ! Stop timer for inter-cycle synchronization - call timer_stop(time_intercycle) - - end do CYCLE_LOOP + end do BATCH_LOOP call timer_stop(time_active) diff --git a/src/material_header.F90 b/src/material_header.F90 index d8e9ca4ea2..916578923b 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -21,8 +21,6 @@ module material_header character(12) :: sab_name ! name of S(a,b) table integer :: sab_table = 0 ! index in sab_tables integer :: sab_nuclide = 0 ! index of nuclide which has S(a,b) table - - integer, allocatable :: rgb(:) ! plotting color end type Material end module material_header diff --git a/src/mesh.F90 b/src/mesh.F90 index 22c2f28c86..37847d58ac 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1,9 +1,14 @@ module mesh use constants + use global use mesh_header use particle_header, only: Particle +#ifdef MPI + use mpi +#endif + implicit none contains @@ -137,4 +142,174 @@ contains end subroutine bin_to_mesh_indices +!=============================================================================== +! COUNT_FISSION_SITES determines the number of fission bank sites in each cell +! of a given mesh. This can be used for a variety of purposes (Shannon entropy, +! CMFD, uniform fission source weighting) +!=============================================================================== + + subroutine count_fission_sites(m, count, total, sites_outside) + + type(StructuredMesh), pointer :: m ! mesh to count sites + real(8), intent(out) :: count(:,:,:) ! weight of sites in each cell + real(8), intent(out) :: total ! total weight of sites + logical, optional :: sites_outside ! were there sites outside mesh? + + integer :: i ! loop index for local fission sites + integer :: ijk(3) ! indices on mesh + real(8) :: weight ! accumulated weight of sites + logical :: in_mesh ! was single site outside mesh? + logical :: outside ! was any site outside mesh? +#ifdef MPI + integer :: n ! total size of count variable +#endif + + ! initialize variables + count = ZERO + weight = ZERO + outside = .false. + + ! loop over fission sites and count how many are in each mesh box + FISSION_SITES: do i = 1, int(n_bank,4) + ! determine scoring bin for entropy mesh + call get_mesh_indices(m, fission_bank(i) % xyz, ijk, in_mesh) + + ! if outside mesh, skip particle + if (.not. in_mesh) then + outside = .true. + cycle + end if + + ! add weight + weight = weight + ONE + + ! add to appropriate mesh box + ! TODO: if tracking weight through bank, add weight instead + count(ijk(1),ijk(2),ijk(3)) = count(ijk(1),ijk(2),ijk(3)) + 1 + end do FISSION_SITES + +#ifdef MPI + ! determine total number of mesh cells + n = size(count,1) * size(count,2) * size(count,3) + + ! collect values from all processors + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, count, n, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) + else + call MPI_REDUCE(count, count, n, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) + end if + + ! Check if there were sites outside the mesh for any processor + if (present(sites_outside)) then + call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, & + MPI_COMM_WORLD, mpi_err) + end if + + ! determine total weight of bank sites + call MPI_REDUCE(weight, total, 1, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) +#else + total = weight + sites_outside = outside +#endif + + end subroutine count_fission_sites + +!=============================================================================== +! MESH_INTERSECT determines if a line between xyz0 and xyz1 intersects the outer +! boundary of the given mesh. This is important for determining whether a track +! will score to a mesh tally. +!=============================================================================== + + function mesh_intersects(m, xyz0, xyz1) result(intersects) + + type(StructuredMesh), pointer :: m + real(8), intent(in) :: xyz0(3) + real(8), intent(in) :: xyz1(3) + logical :: intersects + + real(8) :: x0, y0, z0 ! track start point + real(8) :: x1, y1, z1 ! track end point + real(8) :: xi, yi, zi ! track intersection point with mesh + real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh + real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh + + ! Copy coordinates of starting point + x0 = xyz0(1) + y0 = xyz0(2) + z0 = xyz0(3) + + ! Copy coordinates of ending point + x1 = xyz1(1) + y1 = xyz1(2) + z1 = xyz1(3) + + ! Copy coordinates of mesh lower_left + xm0 = m % lower_left(1) + ym0 = m % lower_left(2) + zm0 = m % lower_left(3) + + ! Copy coordinates of mesh upper_right + xm1 = m % upper_right(1) + ym1 = m % upper_right(2) + zm1 = m % upper_right(3) + + ! Check if line intersects bottom surface -- calculate the intersection + ! point (y,z) + yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then + intersects = .true. + return + end if + + ! Check if line intersects left surface -- calculate the intersection point + ! (x,z) + xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then + intersects = .true. + return + end if + + ! Check if line intersects front surface -- calculate the intersection point + ! (x,y) + xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0) + if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then + intersects = .true. + return + end if + + ! Check if line intersects top surface -- calculate the intersection + ! point (y,z) + yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then + intersects = .true. + return + end if + + ! Check if line intersects right surface -- calculate the intersection point + ! (x,z) + xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then + intersects = .true. + return + end if + + ! Check if line intersects back surface -- calculate the intersection point + ! (x,y) + xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0) + if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then + intersects = .true. + return + end if + + end function mesh_intersects + end module mesh diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index b70bbd45b4..8088972ce0 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -8,13 +8,13 @@ module mesh_header !=============================================================================== type StructuredMesh - integer :: id - integer :: type - integer :: n_dimension - integer, allocatable :: dimension(:) - real(8), allocatable :: lower_left(:) - real(8), allocatable :: upper_right(:) - real(8), allocatable :: width(:) + integer :: id ! user-specified id + integer :: type ! rectangular, hexagonal + integer :: n_dimension ! rank of mesh + integer, allocatable :: dimension(:) ! number of cells in each direction + real(8), allocatable :: lower_left(:) ! lower-left corner of mesh + real(8), allocatable :: upper_right(:) ! upper-right corner of mesh + real(8), allocatable :: width(:) ! width of each mesh cell end type StructuredMesh end module mesh_header diff --git a/src/output.F90 b/src/output.F90 index dc3a9f4fbb..aa9b414a25 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -29,8 +29,8 @@ contains subroutine title() - character(10) :: today_date - character(8) :: today_time + character(10) :: date_ + character(10) :: time_ write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & ' .d88888b. 888b d888 .d8888b.', & @@ -50,20 +50,28 @@ contains ' Developed At: Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE +#ifdef GIT_SHA1 + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1 +#endif ! Write the date and time - call get_today(today_date, today_time) + call date_and_time(DATE=date_, TIME=time_) + date_ = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) + time_ = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A)') & - trim(today_date), trim(today_time) + trim(date_), trim(time_) ! Write information to summary file call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) write(UNIT=UNIT_SUMMARY, FMT=*) & - "Copyright: 2011 Massachusetts Institute of Technology" + "Copyright: 2011-2012 Massachusetts Institute of Technology" write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE +#ifdef GIT_SHA1 + write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1 +#endif write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A,1X,A)') & - trim(today_date), trim(today_time) + trim(date_), trim(time_) ! Write information on number of processors #ifdef MPI @@ -86,8 +94,8 @@ contains integer, optional :: unit ! unit to write to integer, optional :: level ! specified header level - integer :: n - integer :: m + integer :: n ! number of = signs on left + integer :: m ! number of = signs on right integer :: unit_ ! unit to write to integer :: header_level ! actual header level character(MAX_LINE_LEN) :: line @@ -137,14 +145,15 @@ contains subroutine write_message(level) - integer, optional :: level + integer, optional :: level ! verbosity level - integer :: n_lines - integer :: i + integer :: n_lines ! number of lines needed + integer :: i ! index for lines ! Only allow master to print to screen if (.not. master .and. present(level)) return + ! TODO: Take care of line wrapping so words don't get cut off if (.not. present(level) .or. level <= verbosity) then n_lines = (len_trim(message)-1)/79 + 1 do i = 1, n_lines @@ -154,55 +163,13 @@ contains end subroutine write_message -!=============================================================================== -! GET_TODAY determines the date and time at which the program began execution -! and returns it in a readable format -!=============================================================================== - - subroutine get_today(today_date, today_time) - - character(10), intent(out) :: today_date - character(8), intent(out) :: today_time - - integer :: val(8) - character(8) :: date_ - character(10) :: time_ - character(5) :: zone - - call date_and_time(date_, time_, zone, val) - ! val(1) = year (YYYY) - ! val(2) = month (MM) - ! val(3) = day (DD) - ! val(4) = timezone - ! val(5) = hours (HH) - ! val(6) = minutes (MM) - ! val(7) = seconds (SS) - ! val(8) = milliseconds - - if (val(2) < 10) then - if (val(3) < 10) then - today_date = date_(6:6) // "/" // date_(8:8) // "/" // date_(1:4) - else - today_date = date_(6:6) // "/" // date_(7:8) // "/" // date_(1:4) - end if - else - if (val(3) < 10) then - today_date = date_(5:6) // "/" // date_(8:8) // "/" // date_(1:4) - else - today_date = date_(5:6) // "/" // date_(7:8) // "/" // date_(1:4) - end if - end if - today_time = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) - - end subroutine get_today - !=============================================================================== ! PRINT_PARTICLE displays the attributes of a particle !=============================================================================== subroutine print_particle() - integer :: i + integer :: i ! index for coordinate levels type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() @@ -263,6 +230,7 @@ contains write(ou,*) ' Surface = ' // to_str(sign(s % id, p % surface)) end if + ! Display weight, energy, grid index, and interpolation factor write(ou,*) ' Weight = ' // to_str(p % wgt) write(ou,*) ' Energy = ' // to_str(p % E) write(ou,*) ' IE = ' // to_str(p % IE) @@ -298,27 +266,35 @@ contains subroutine print_cell(c, unit) type(Cell), pointer :: c - integer, optional :: unit + integer, optional :: unit ! specified unit to write to - integer :: temp - integer :: i - integer :: unit_ + integer :: index_cell ! index in cells array + integer :: i ! loop index for surfaces + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Universe), pointer :: u => null() type(Lattice), pointer :: l => null() type(Material), pointer :: m => null() + ! Set unit to stdout if not already set if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id for cell write(unit_,*) 'Cell ' // to_str(c % id) - temp = dict_get_key(cell_dict, c % id) - write(unit_,*) ' Array Index = ' // to_str(temp) + + ! Find index in cells array and write + index_cell = dict_get_key(cell_dict, c % id) + write(unit_,*) ' Array Index = ' // to_str(index_cell) + + ! Write what universe this cell is in u => universes(c % universe) write(unit_,*) ' Universe = ' // to_str(u % id) + + ! Write information on fill for cell select case (c % type) case (CELL_NORMAL) write(unit_,*) ' Fill = NONE' @@ -329,12 +305,16 @@ contains l => lattices(c % fill) write(unit_,*) ' Fill = Lattice ' // to_str(l % id) end select + + ! Write information on material if (c % material == 0) then write(unit_,*) ' Material = NONE' else m => materials(c % material) write(unit_,*) ' Material = ' // to_str(m % id) end if + + ! Write surface specification string = "" do i = 1, c % n_surfaces select case (c % surfaces(i)) @@ -364,24 +344,31 @@ contains type(Universe), pointer :: univ integer, optional :: unit - integer :: i - integer :: unit_ + integer :: i ! loop index for cells in this universe + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() - type(Universe), pointer :: base_u + type(Universe), pointer :: base_u => null() + ! Set default unit to stdout if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Get a pointer to the base universe base_u => universes(BASE_UNIVERSE) + ! Write user-specified id for this universe write(unit_,*) 'Universe ' // to_str(univ % id) + + ! If this is the base universe, indicate so if (associated(univ, base_u)) then write(unit_,*) ' Base Universe' end if + + ! Write list of cells in this universe string = "" do i = 1, univ % n_cells c => cells(univ % cells(i)) @@ -401,14 +388,16 @@ contains type(Lattice), pointer :: lat integer, optional :: unit - integer :: unit_ + integer :: unit_ ! unit to write to + ! set default unit if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write information about lattice write(unit_,*) 'Lattice ' // to_str(lat % id) write(unit_,*) ' n_x = ' // to_str(lat % n_x) write(unit_,*) ' n_y = ' // to_str(lat % n_y) @@ -427,19 +416,23 @@ contains subroutine print_surface(surf, unit) type(Surface), pointer :: surf - integer, optional :: unit + integer, optional :: unit ! specified unit to write to - integer :: i - integer :: unit_ + integer :: i ! loop index for coefficients + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string + ! set default unit if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id of surface write(unit_,*) 'Surface ' // to_str(surf % id) + + ! Write type of surface select case (surf % type) case (SURF_PX) string = "X Plane" @@ -466,12 +459,14 @@ contains end select write(unit_,*) ' Type = ' // trim(string) + ! Write coefficients for this surface string = "" do i = 1, size(surf % coeffs) string = trim(string) // ' ' // to_str(surf % coeffs(i), 4) end do write(unit_,*) ' Coefficients = ' // trim(string) + ! Write neighboring cells on positive side of this surface string = "" if (allocated(surf % neighbor_pos)) then do i = 1, size(surf % neighbor_pos) @@ -480,6 +475,7 @@ contains end if write(unit_,*) ' Positive Neighbors = ' // trim(string) + ! Write neighboring cells on negative side of this surface string = "" if (allocated(surf % neighbor_neg)) then do i = 1, size(surf % neighbor_neg) @@ -487,6 +483,8 @@ contains end do end if write(unit_,*) ' Negative Neighbors =' // trim(string) + + ! Write boundary condition for this surface select case (surf % bc) case (BC_TRANSMIT) write(unit_,*) ' Boundary Condition = Transmission' @@ -510,12 +508,13 @@ contains type(Material), pointer :: mat integer, optional :: unit - integer :: i - integer :: unit_ - real(8) :: density + integer :: i ! loop index for nuclides + integer :: unit_ ! unit to write to + real(8) :: density ! density in atom/b-cm character(MAX_LINE_LEN) :: string type(Nuclide), pointer :: nuc => null() + ! set default unit to stdout if not specified if (present(unit)) then unit_ = unit else @@ -556,9 +555,9 @@ contains type(TallyObject), pointer :: t integer, optional :: unit - integer :: i - integer :: id - integer :: unit_ + integer :: i ! index for filter or score bins + integer :: id ! user-specified id + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() @@ -566,14 +565,17 @@ contains type(Material), pointer :: m => null() type(StructuredMesh), pointer :: sm => null() + ! set default unit to stdout if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id of tally write(unit_,*) 'Tally ' // to_str(t % id) + ! Write any cells bins if present if (t % n_filter_bins(FILTER_CELL) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_CELL) @@ -584,6 +586,7 @@ contains write(unit_, *) ' Cell Bins:' // trim(string) end if + ! Write any surface bins if present if (t % n_filter_bins(FILTER_SURFACE) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_SURFACE) @@ -594,6 +597,7 @@ contains write(unit_, *) ' Surface Bins:' // trim(string) end if + ! Write any universe bins if present if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_UNIVERSE) @@ -604,6 +608,7 @@ contains write(unit_, *) ' Material Bins:' // trim(string) end if + ! Write any material bins if present if (t % n_filter_bins(FILTER_MATERIAL) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_MATERIAL) @@ -614,6 +619,7 @@ contains write(unit_, *) ' Material Bins:' // trim(string) end if + ! Write any mesh bins if present if (t % n_filter_bins(FILTER_MESH) > 0) then string = "" id = t % mesh @@ -625,6 +631,7 @@ contains write(unit_, *) ' Mesh Bins:' // trim(string) end if + ! Write any birth region bins if present if (t % n_filter_bins(FILTER_CELLBORN) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_CELLBORN) @@ -635,6 +642,7 @@ contains write(unit_, *) ' Birth Region Bins:' // trim(string) end if + ! Write any incoming energy bins if present if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_ENERGYIN) + 1 @@ -644,6 +652,7 @@ contains write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if + ! Write any outgoing energy bins if present if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_ENERGYOUT) + 1 @@ -653,6 +662,7 @@ contains write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if + ! Write any score bins if present if (t % n_score_bins > 0) then string = "" do i = 1, t % n_score_bins @@ -671,7 +681,7 @@ contains string = trim(string) // ' nu-fission' end select end do - write(unit_,*) ' Macro Reactions:' // trim(string) + write(unit_,*) ' Scores:' // trim(string) end if write(unit_,*) @@ -684,7 +694,7 @@ contains subroutine print_geometry() - integer :: i + integer :: i ! loop index for various arrays type(Surface), pointer :: s => null() type(Cell), pointer :: c => null() type(Universe), pointer :: u => null() @@ -732,12 +742,12 @@ contains type(Nuclide), pointer :: nuc integer, optional :: unit - integer :: i - integer :: unit_ - integer :: size_total - integer :: size_xs - integer :: size_angle - integer :: size_energy + integer :: i ! loop index over nuclides + integer :: unit_ ! unit to write to + integer :: size_total ! memory used by nuclide (bytes) + integer :: size_xs ! memory used for cross-sections (bytes) + integer :: size_angle ! memory used for angle distributions (bytes) + integer :: size_energy ! memory used for energy distributions (bytes) type(Reaction), pointer :: rxn => null() type(UrrData), pointer :: urr => null() @@ -821,7 +831,7 @@ contains subroutine print_summary() - integer :: i + integer :: i ! loop index character(15) :: string type(Material), pointer :: m => null() type(TallyObject), pointer :: t => null() @@ -830,8 +840,9 @@ contains call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY) if (problem_type == PROB_CRITICALITY) then write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality' - write(UNIT_SUMMARY,101) 'Number of Cycles:', n_cycles - write(UNIT_SUMMARY,101) 'Number of Inactive Cycles:', n_inactive + write(UNIT_SUMMARY,101) 'Number of Batches:', n_batches + write(UNIT_SUMMARY,101) 'Number of Inactive Batches:', n_inactive + write(UNIT_SUMMARY,101) 'Generations per Batch:', gen_per_batch elseif (problem_type == PROB_SOURCE) then write(UNIT_SUMMARY,100) 'Problem type:', 'External Source' end if @@ -892,41 +903,32 @@ contains subroutine print_plot() - integer i - type(Plot), pointer :: pl => null() + integer :: i ! loop index for plots + type(Plot), pointer :: pl => null() ! Display header for plotting call header("PLOTTING SUMMARY") - do i=1,n_plots + do i = 1, n_plots pl => plots(i) - ! Print plot id + ! Write plot id write(ou,100) "Plot ID:", trim(to_str(pl % id)) - ! Print plotting origin + ! Write plotting origin write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // & " " // trim(to_str(pl % origin(2))) // " " // & trim(to_str(pl % origin(3))) - ! Print plotting width + ! Write plotting width if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Width:", trim(to_str(pl % width(1))) // & " " // trim(to_str(pl % width(2))) - write(ou,100) "Coloring:", trim(to_str(pl % color)) + write(ou,100) "Coloring:", trim(to_str(pl % color_by)) write(ou,100) "Basis:", trim(to_str(pl % basis)) write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // & trim(to_str(pl % pixels(2))) - - else if (pl % type == PLOT_TYPE_POINTS) then - - write(ou,100) "Width:", trim(to_str(pl % width(1))) // & - " " // trim(to_str(pl % width(2))) // " " & - // trim(to_str(pl % width(3))) - write(ou,100) "Coloring:", trim(to_str(pl % color)) - write(ou,100) "Ray Spacing:", trim(to_str(pl % aspect)) - end if write(ou,*) @@ -945,8 +947,8 @@ contains subroutine print_runtime() - integer(8) :: total_particles - real(8) :: speed + integer(8) :: total_particles ! total # of particles simulated + real(8) :: speed ! # of neutrons/second character(15) :: string ! display header block @@ -959,10 +961,16 @@ contains write(ou,100) "Total time in simulation", time_inactive % elapsed + & time_active % elapsed write(ou,100) " Time in transport only", time_transport % elapsed +<<<<<<< HEAD:src/output.F90 write(ou,100) " Time in inactive cycles", time_inactive % elapsed write(ou,100) " Time in active cycles", time_active % elapsed write(ou,100) " Time between cycles", time_intercycle % elapsed if(cmfd_on) write(ou,100) "Total CMFD time", time_cmfd % elapsed +======= + write(ou,100) " Time in inactive batches", time_inactive % elapsed + write(ou,100) " Time in active batches", time_active % elapsed + write(ou,100) " Time between generations", time_intercycle % elapsed +>>>>>>> master:src/output.F90 write(ou,100) " Accumulating tallies", time_ic_tallies % elapsed write(ou,100) " Sampling source sites", time_ic_sample % elapsed write(ou,100) " SEND/RECV source sites", time_ic_sendrecv % elapsed @@ -970,7 +978,7 @@ contains write(ou,100) "Total time elapsed", time_total % elapsed ! display calculate rate and final keff - total_particles = n_particles * n_cycles + total_particles = n_particles * n_batches * gen_per_batch speed = real(total_particles) / (time_inactive % elapsed + & time_active % elapsed) string = to_str(speed) @@ -1023,7 +1031,8 @@ contains end subroutine create_summary_file !=============================================================================== -! CREATE_XS_SUMMARY_FILE +! CREATE_XS_SUMMARY_FILE creates an output file to write information about the +! cross section tables used in the simulation. !=============================================================================== subroutine create_xs_summary_file() diff --git a/src/particle_header.F90 b/src/particle_header.F90 index a9cfa12c24..4330dfb882 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -52,7 +52,8 @@ module particle_header real(8) :: last_E ! last energy ! What event last took place - integer :: event + integer :: event ! scatter, absorption, fission + integer :: event_nuclide ! index in nuclides array ! Post-collision physical data integer :: n_bank ! number of fission sites banked diff --git a/src/physics.F90 b/src/physics.F90 index 8ed65c8c7f..ca9d389ec1 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -3,7 +3,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants use cross_section, only: calculate_xs - use endf, only: reaction_name, is_fission, is_scatter + use endf, only: reaction_name use error, only: fatal_error, warning use fission, only: nu_total, nu_prompt, nu_delayed use geometry, only: find_cell, distance_to_boundary, cross_surface, & @@ -11,6 +11,7 @@ module physics use geometry_header, only: Universe, BASE_UNIVERSE use global use interpolation, only: interpolate_tab1 + use material_header, only: Material use output, only: write_message use particle_header, only: LocalCoord use random_lcg, only: prn @@ -39,6 +40,11 @@ contains logical :: found_cell ! found cell which particle is in? type(LocalCoord), pointer :: coord => null() + if (verbosity >= 9 .or. trace) then + message = "Simulating Particle " // trim(to_str(p % id)) + call write_message() + end if + if (p % coord % cell == NONE) then call find_cell(found_cell) ! Particle couldn't be located @@ -51,17 +57,6 @@ contains p % cell_born = p % coord % cell end if - if (verbosity >= 9 .or. trace) then - message = "Simulating Particle " // trim(to_str(p % id)) - call write_message() - end if - - if (verbosity >= 10 .or. trace) then - message = " Born in cell " // trim(to_str(& - cells(p % coord % cell) % id)) - call write_message() - end if - ! Initialize number of events to zero n_event = 0 @@ -133,7 +128,7 @@ contains p % surface = NONE call collision() - ! Save coordinates at collision for tallying purposes + ! Save coordinates for tallying purposes p % last_xyz = p % coord0 % xyz ! Set last material to none since cross sections will need to be @@ -188,9 +183,18 @@ contains ! Display information about collision if (verbosity >= 10 .or. trace) then -!!$ message = " " // trim(reaction_name(MT)) // ". Energy = " // & -!!$ trim(to_str(p % E * 1e6_8)) // " eV." -!!$ call write_message() + select case (p % event) + case (EVENT_SCATTER) + message = " Scattered off of" + case (EVENT_FISSION) + message = " Fissioned in" + case (EVENT_ABSORB) + message = " Absorbed in" + end select + message = trim(message) // " " // trim(nuclides(& + p % event_nuclide) % name) // ". Energy = " // & + trim(to_str(p % E * 1e6_8)) // " eV." + call write_message() end if ! check for very low energy @@ -268,6 +272,9 @@ contains IE = micro_xs(index_nuclide) % index_grid f = micro_xs(index_nuclide) % interp_factor + ! Save which nuclide particle had collision with + p % event_nuclide = index_nuclide + ! ========================================================================== ! DISAPPEARANCE REACTIONS (ANALOG) OR IMPLICIT CAPTURE (SURVIVAL BIASING) @@ -2055,77 +2062,4 @@ contains end function watt_spectrum -!=============================================================================== -! WIGNER samples a Wigner distribution of energy level spacings. Note that this -! scheme is C50 in the Monte Carlo Sampler from Los Alamos (LA-9721-MS). -!=============================================================================== - - function wigner(D_avg) result (D) - - real(8), intent(in) :: D_avg ! average level spacing - real(8) :: D ! sampled level spacing - - real(8) :: c - - c = -4.*D_avg*D_avg/PI * log(prn()) - D = sqrt(c) - - end function wigner - -!=============================================================================== -! CHI_SQUARED samples a chi-squared distribution with n degrees of freedom. The -! distribution of resonance widths in the unresolved region is given by a -! chi-squared distribution. For the special case of n=1, this is a Porter-Thomas -! distribution. For cases with n odd, rule C64 is used whereas for cases with n -! even, rule C45 is used. -!=============================================================================== - - function chi_squared(n, G_avg) result(G) - - integer, intent(in) :: n ! number of degrees of freedom - real(8), intent(in), optional :: G_avg ! average resonance width - - integer :: i ! loop index - real(8) :: G ! sampled random variable (or resonance width) - real(8) :: x, y, c ! dummy variables - real(8) :: r1, r2 ! psuedorandom numbers - - select case (mod(n,2)) - case (0) - ! Even number of degrees of freedom can be sampled via rule C45. We can - ! sample x as -2/n*log(product(r_i, i = 1 to n/2)) - x = ONE - do i = 1, n/2 - x = x * prn() - end do - x = -2./n * log(x) - - case (1) - ! Odd number of degrees of freedom can be sampled via rule C64. We can - ! sample x as -2/n*(log(r)*cos^2(pi/2*r) + log(product(r_i, i = 1 to - ! floor(n/2))) - - ! Note that we take advantage of integer division on n/2 - y = ONE - do i = 1, n/2 - y = y * prn() - end do - - r1 = prn() - r2 = prn() - c = cos(PI/2.*r2) - x = -2./n * (log(y) + log(r1)*c*c) - end select - - ! If sampling a chi-squared distribution for a resonance width and the - ! average resonance width has been given, return the sampled resonance - ! width. - if (present(G_avg)) then - G = x * G_avg - else - G = x - end if - - end function chi_squared - end module physics diff --git a/src/plot.F90 b/src/plot.F90 index cd2aec0caa..7582130454 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -2,12 +2,11 @@ module plotter use constants use error, only: fatal_error - use geometry, only: find_cell, distance_to_boundary, cross_surface, & - cross_lattice, cell_contains + use geometry, only: find_cell use geometry_header, only: Universe, BASE_UNIVERSE use global use output, only: write_message - use particle_header, only: LocalCoord, deallocate_coord + use particle_header, only: deallocate_coord use plot_header use ppmlib, only: Image, init_image, allocate_image, & deallocate_image, set_pixel @@ -18,17 +17,16 @@ module plotter contains - !=============================================================================== -! RUN_PLOT +! RUN_PLOT controls the logic for making one or many plots !=============================================================================== subroutine run_plot() - integer :: i - type(Plot), pointer :: pl => null() + integer :: i ! loop index for plots + type(Plot), pointer :: pl => null() - do i=1,n_plots + do i = 1, n_plots pl => plots(i) ! Display output message @@ -36,31 +34,34 @@ contains call write_message(5) if (pl % type == PLOT_TYPE_SLICE) then - ! create 2d image call create_ppm(pl) - end if end do end subroutine run_plot !=============================================================================== -! create_ppm +! CREATE_PPM creates an image based on user input from a plots.xml +! specification in the portable pixmap format (PPM) !=============================================================================== + subroutine create_ppm(pl) - type(Plot), pointer :: pl + type(Plot), pointer :: pl - type(Image) :: img - integer :: in_i, out_i - integer :: x, y - integer :: r, g, b - real(8) :: in_pixel, out_pixel - real(8) :: xyz(3) - logical :: found_cell - type(Cell), pointer :: c => null() + integer :: in_i + integer :: out_i + integer :: x, y ! pixel location + integer :: r, g, b ! colors (red, green, blue) from 0-255 + real(8) :: in_pixel + real(8) :: out_pixel + real(8) :: xyz(3) + logical :: found_cell + type(Image) :: img + type(Cell), pointer :: c => null() + ! Initialize and allocate space for image call init_image(img) call allocate_image(img, pl % pixels(1), pl % pixels(2)) @@ -70,7 +71,7 @@ contains in_pixel = pl % width(1)/dble(pl % pixels(1)) out_pixel = pl % width(2)/dble(pl % pixels(2)) xyz(1) = pl % origin(1) - pl % width(1) / 2.0 - xyz(2) = pl % origin(2) - pl % width(2) / 2.0 + xyz(2) = pl % origin(2) + pl % width(2) / 2.0 xyz(3) = pl % origin(3) else if (pl % basis == PLOT_BASIS_XZ) then in_i = 1 @@ -79,7 +80,7 @@ contains out_pixel = pl % width(2)/dble(pl % pixels(2)) xyz(1) = pl % origin(1) - pl % width(1) / 2.0 xyz(2) = pl % origin(2) - xyz(3) = pl % origin(3) - pl % width(2) / 2.0 + xyz(3) = pl % origin(3) + pl % width(2) / 2.0 else if (pl % basis == PLOT_BASIS_YZ) then in_i = 2 out_i = 3 @@ -87,7 +88,7 @@ contains out_pixel = pl % width(2)/dble(pl % pixels(2)) xyz(1) = pl % origin(1) xyz(2) = pl % origin(2) - pl % width(1) / 2.0 - xyz(3) = pl % origin(3) - pl % width(2) / 2.0 + xyz(3) = pl % origin(3) + pl % width(2) / 2.0 end if ! allocate and initialize particle @@ -97,8 +98,8 @@ contains p % coord % uvw = (/ 1, 0, 0 /) p % coord % universe = BASE_UNIVERSE - do y=1, img%height - do x=1, img%width + do y = 1, img % height + do x = 1, img % width call deallocate_coord(p % coord0 % next) p % coord => p % coord0 @@ -106,65 +107,74 @@ contains call find_cell(found_cell) if (.not. found_cell) then - r = 255 - g = 255 - b = 255 + ! If no cell, revert to default color + r = pl % not_found % rgb(1) + g = pl % not_found % rgb(2) + b = pl % not_found % rgb(3) else - c => cells(p % coord % cell) - if (pl % color == PLOT_COLOR_MATS) then - r = materials(c % material) % rgb(1) - g = materials(c % material) % rgb(2) - b = materials(c % material) % rgb(3) - else if (pl % color == PLOT_COLOR_CELLS) then - r = c % rgb(1) - g = c % rgb(2) - b = c % rgb(3) - else - r = 0 - g = 0 - b = 0 - end if + if (pl % color_by == PLOT_COLOR_MATS) then + ! Assign color based on material + c => cells(p % coord % cell) + r = pl % colors(c % material) % rgb(1) + g = pl % colors(c % material) % rgb(2) + b = pl % colors(c % material) % rgb(3) + else if (pl % color_by == PLOT_COLOR_CELLS) then + ! Assign color based on cell + r = pl % colors(p % coord % cell) % rgb(1) + g = pl % colors(p % coord % cell) % rgb(2) + b = pl % colors(p % coord % cell) % rgb(3) + else + r = 0 + g = 0 + b = 0 + end if end if + ! Create a pixel at (x,y) with color (r,g,b) call set_pixel(img, x, y, r, g, b) + ! Advance pixel in first direction p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel end do + ! Advance pixel in second direction p % coord0 % xyz(in_i) = xyz(in_i) - p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) + out_pixel + p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel end do + ! Write out the ppm to a file call output_ppm(pl,img) + ! Free up space call deallocate_image(img) end subroutine create_ppm - !=============================================================================== -! output_ppm +! OUTPUT_PPM writes out a previously generated image to a PPM file !=============================================================================== - subroutine output_ppm(pl,img) - type(Plot), pointer :: pl - type(Image), intent(in) :: img + subroutine output_ppm(pl, img) - integer :: i, j - character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file + type(Plot), pointer :: pl + type(Image), intent(in) :: img - path_plot = trim(path_input) // "slice" // trim(to_str(pl % id)) // ".ppm" - open(UNIT=UNIT_PLOT, FILE=path_plot) + integer :: i ! loop index for height + integer :: j ! loop index for width + ! Open PPM file for writing + open(UNIT=UNIT_PLOT, FILE=pl % path_plot) + + ! Write header write(UNIT_PLOT, '(A2)') 'P6' write(UNIT_PLOT, '(I0,'' '',I0)') img%width, img%height write(UNIT_PLOT, '(A)') '255' - do j=1, img%height - do i=1, img%width + ! Write color for each pixel + do j = 1, img % height + do i = 1, img % width write(UNIT_PLOT, '(3A1)', advance='no') achar(img%red(i,j)), & - achar(img%green(i,j)), & - achar(img%blue(i,j)) + achar(img%green(i,j)), achar(img%blue(i,j)) end do end do @@ -173,5 +183,4 @@ contains end subroutine output_ppm - end module plotter diff --git a/src/plot_header.F90 b/src/plot_header.F90 index bfe6f9b071..6268ffadf9 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -1,20 +1,32 @@ module plot_header + use constants + implicit none !=============================================================================== -! PLOT hold plot information +! ObjectColor holds color information for plotted objects +!=============================================================================== + + type ObjectColor + integer :: rgb(3) + end type ObjectColor + +!=============================================================================== +! PLOT holds plot information !=============================================================================== type Plot integer :: id ! Unique ID + character(MAX_LINE_LEN) :: path_plot ! path for plot file integer :: type ! Type - integer :: color ! quantity to color regions by + integer :: color_by ! quantity to color regions by real(8) :: origin(3) ! xyz center of plot location - real(8) :: aspect ! spacing between rays in raytracer real(8) :: width(3) ! xyz widths of plot integer :: basis ! direction of plot slice integer :: pixels(2) ! pixel width/height of plot slice + type(ObjectColor) :: not_found ! color for positions where no cell found + type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats end type Plot integer :: PLOT_TYPE_SLICE = 1 diff --git a/src/ppmlib.F90 b/src/ppmlib.F90 index 2c9687617a..d67268412e 100644 --- a/src/ppmlib.F90 +++ b/src/ppmlib.F90 @@ -13,6 +13,10 @@ module ppmlib contains +!=============================================================================== +! INIT_IMAGE initializes the Image derived type +!=============================================================================== + subroutine init_image(img) type(Image), intent(out) :: img @@ -26,29 +30,46 @@ contains end subroutine init_image +!=============================================================================== +! ALLOCATE_IMAGE sets the width and height of an image and allocates color +! arrays +!=============================================================================== + subroutine allocate_image(img, w, h) type(Image), intent(inout) :: img - integer, intent(in) :: w, h + integer, intent(in) :: w ! width of image + integer, intent(in) :: h ! height of image + ! allocate red, green, and blue array allocate(img % red(w, h)) allocate(img % green(w, h)) allocate(img % blue(w, h)) + ! set width and height img % width = w img % height = h end subroutine allocate_image + +!=============================================================================== +! DEALLOCATE_IMAGE +!=============================================================================== subroutine deallocate_image(img) type(Image) :: img - if ( associated(img % red) ) deallocate(img % red) - if ( associated(img % green) ) deallocate(img % green) - if ( associated(img % blue) ) deallocate(img % blue) + if (associated(img % red)) deallocate(img % red) + if (associated(img % green)) deallocate(img % green) + if (associated(img % blue)) deallocate(img % blue) + end subroutine deallocate_image +!=============================================================================== +! INSIDE_IMAGE determines whether a point (x,y) is inside the image +!=============================================================================== + function inside_image(img, x, y) result(inside) @@ -58,13 +79,15 @@ contains inside = .false. - if ( ( x < img % width) .and. & - ( y < img % height ) .and. & - ( x >= 0 ) .and. & - ( y >= 0 ) ) inside = .true. + if ((x < img % width) .and. (y < img % height) .and. & + (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 +!=============================================================================== function valid_image(img) result(valid) @@ -73,29 +96,32 @@ contains valid = .false. - if ( img % width == 0 ) return - if ( img % height == 0 ) return - if ( .not. associated(img % red) .or. & + if (img % width == 0) return + if (img % height == 0) return + if (.not. associated(img % red) .or. & .not. associated(img % green) .or. & - .not. associated(img % blue) ) return + .not. associated(img % blue)) return valid = .true. end function valid_image - +!=============================================================================== +! SET_PIXEL sets the +!=============================================================================== + subroutine set_pixel(img, x, y, r, g, b) type(Image), intent(inout) :: img - integer, intent(in) :: x, y - integer, intent(in) :: r, g, b + integer, intent(in) :: x, y ! coordinates + integer, intent(in) :: r, g, b ! red, green, and blue - if ( inside_image(img, x, y) .and. valid_image(img)) then - img % red(x+1,y+1) = mod(abs(r), 256) - img % green(x+1, y+1) = mod(abs(g), 256) - img % blue(x+1, y+1) = mod(abs(b), 256) + if (inside_image(img, x, y) .and. valid_image(img)) then + img % red(x+1,y+1) = mod(abs(r), 256) + img % green(x+1, y+1) = mod(abs(g), 256) + img % blue(x+1, y+1) = mod(abs(b), 256) end if - end subroutine set_pixel + end subroutine set_pixel end module ppmlib diff --git a/src/search.F90 b/src/search.F90 index 3cdb51bfcb..5968113fec 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -1,6 +1,5 @@ module search - use constants, only: ONE, MAX_LINE_LEN use error, only: fatal_error use global, only: message diff --git a/src/source.F90 b/src/source.F90 index 1b7ec51f86..cd31dc2ec2 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,7 +1,7 @@ module source use bank_header, only: Bank - use constants, only: ONE + use constants use error, only: fatal_error use geometry_header, only: BASE_UNIVERSE use global @@ -11,6 +11,10 @@ module source use random_lcg, only: prn, set_particle_seed use string, only: to_str +#ifdef MPI + use mpi +#endif + implicit none contains @@ -43,6 +47,13 @@ contains ! Determine maximum amount of particles to simulate on each processor maxwork = ceiling(real(n_particles)/n_procs,8) + ! ID's of first and last source particles + bank_first = rank*maxwork + 1 + bank_last = min((rank+1)*maxwork, n_particles) + + ! number of particles for this processor + work = bank_last - bank_first + 1 + ! Allocate source bank allocate(source_bank(maxwork), STAT=alloc_err) if (alloc_err /= 0) then @@ -69,16 +80,16 @@ contains call fatal_error() end if + ! Read the source from a binary file instead of sampling from some assumed + ! source distribution + if (external_source % type == SRC_FILE) then + call read_source_binary() + return + end if + ! Initialize first cycle source bank do i = 0, n_procs - 1 if (rank == i) then - ! ID's of first and last source particles - bank_first = rank*maxwork + 1 - bank_last = min((rank+1)*maxwork, n_particles) - - ! number of particles for this processor - work = bank_last - bank_first + 1 - do j = 1, work id = bank_first + j - 1 source_bank(j) % id = id @@ -93,6 +104,8 @@ contains p_max = external_source % values(4:6) r = (/ (prn(), k = 1,3) /) source_bank(j) % xyz = p_min + r*(p_max - p_min) + case (SRC_POINT) + source_bank(j) % xyz = external_source % values end select ! sample angle @@ -145,12 +158,13 @@ contains p % id = bank_first + index_source - 1 ! set random number seed - particle_seed = (current_cycle - 1)*n_particles + p % id + particle_seed = ((current_batch - 1)*gen_per_batch + & + current_gen - 1)*n_particles + p % id call set_particle_seed(particle_seed) ! set particle trace trace = .false. - if (current_cycle == trace_cycle .and. & + if (current_batch == trace_batch .and. current_gen == trace_gen .and. & p % id == trace_particle) trace = .true. end subroutine get_source_particle @@ -186,4 +200,128 @@ contains end subroutine initialize_particle +!=============================================================================== +! WRITE_SOURCE writes out the final source distribution to a binary file that +! can be used as a starting source in a new simulation +!=============================================================================== + + subroutine write_source_binary() + +#ifdef MPI + integer :: fh ! file handle + integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file) + + ! ========================================================================== + ! PARALLEL I/O USING MPI-2 ROUTINES + + ! Open binary source file for reading + call MPI_FILE_OPEN(MPI_COMM_WORLD, 'source.binary', MPI_MODE_CREATE + & + MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpi_err) + + if (master) then + offset = 0 + call MPI_FILE_WRITE_AT(fh, offset, n_particles, 1, MPI_INTEGER8, & + MPI_STATUS_IGNORE, mpi_err) + end if + + ! Set proper offset for source data on this processor + offset = 8*(1 + rank*maxwork*8) + + ! Write all source sites + call MPI_FILE_WRITE_AT(fh, offset, source_bank(1), work, MPI_BANK, & + MPI_STATUS_IGNORE, mpi_err) + + ! Close binary source file + call MPI_FILE_CLOSE(fh, mpi_err) + +#else + ! ========================================================================== + ! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES + + ! Open binary source file for writing + open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='replace', & + ACCESS='stream') + + ! Write the number of particles + write(UNIT=UNIT_SOURCE) n_particles + + ! Write information from the source bank + write(UNIT=UNIT_SOURCE) source_bank(1:work) + + ! Close binary source file + close(UNIT=UNIT_SOURCE) +#endif + + end subroutine write_source_binary + +!=============================================================================== +! READ_SOURCE reads a source distribution from a source.binary file and +! initializes the source bank +!=============================================================================== + + subroutine read_source_binary() + + integer(8) :: n_sites ! number of sites in binary file +#ifdef MPI + integer :: fh ! file handle + integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file) +#endif + +#ifdef MPI + ! ========================================================================== + ! PARALLEL I/O USING MPI-2 ROUTINES + + ! Open binary source file for reading + call MPI_FILE_OPEN(MPI_COMM_WORLD, 'source.binary', MPI_MODE_RDONLY, & + MPI_INFO_NULL, fh, mpi_err) + + ! Read number of source sites in file + offset = 0 + call MPI_FILE_READ_AT(fh, offset, n_sites, 1, MPI_INTEGER8, & + MPI_STATUS_IGNORE, mpi_err) + + ! Check that number of source sites matches + if (n_sites /= n_particles) then + message = "No support yet for source files of different size than & + &specified number of particles per generation." + call fatal_error() + end if + + ! Set proper offset for source data on this processor + offset = 8*(1 + rank*maxwork*8) + + ! Read all source sites + call MPI_FILE_READ_AT(fh, offset, source_bank(1), work, MPI_BANK, & + MPI_STATUS_IGNORE, mpi_err) + + ! Close binary source file + call MPI_FILE_CLOSE(fh, mpi_err) + +#else + ! ========================================================================== + ! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES + + ! Open binary source file for reading + open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='old', & + ACCESS='stream') + + ! Read number of source sites in file + read(UNIT=UNIT_SOURCE) n_sites + + ! Check that number of source sites matches + if (n_sites /= n_particles) then + message = "No support yet for source files of different size than & + &specified number of particles per generation." + call fatal_error() + end if + + ! Read position, angle, and energy + read(UNIT=UNIT_SOURCE) source_bank(1:work) + + ! Close binary source file + close(UNIT=UNIT_SOURCE) +#endif + + end subroutine read_source_binary + end module source diff --git a/src/source_header.F90 b/src/source_header.F90 index 077e6338d8..638a960244 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -8,7 +8,7 @@ module source_header !=============================================================================== type ExtSource - integer :: type ! type of source, e.g. 'box' or 'cell' + integer :: type ! type of source, e.g. 'box' or 'point' real(8), allocatable :: values(:) ! values for particular source type end type ExtSource diff --git a/src/tally.F90 b/src/tally.F90 index 663cb25b85..932dee2bd3 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3,7 +3,8 @@ module tally use constants use error, only: fatal_error use global - use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices + use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, & + mesh_indices_to_bin, mesh_intersects use mesh_header, only: StructuredMesh use output, only: header use search, only: binary_search @@ -345,6 +346,20 @@ contains score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu) + case (SCORE_DIFFUSION) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Temporarily store the scattering cross section + score = material_xs % total - material_xs % absorption + + ! Since this only gets tallied at every scattering event, the flux + ! estimator is 1/Sigma_s. Therefore, the diffusion coefficient + ! times flux is 1/(3*Sigma_s*(Sigma_t - mu*Sigma_s)). + + score = last_wgt / (3.0_8 * score * (material_xs % total - & + mu * score)) + case (SCORE_N_1N) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle @@ -533,6 +548,14 @@ contains do i = 1, n_tracklength_tallies t => tallies(tracklength_tallies(i)) + ! Check if this tally has a mesh filter -- if so, we treat it separately + ! since multiple bins can be scored to with a single track + + if (t % n_filter_bins(FILTER_MESH) > 0) then + call score_tl_on_mesh(tracklength_tallies(i), distance) + cycle + end if + ! ======================================================================= ! DETERMINE SCORING BIN COMBINATION @@ -592,6 +615,239 @@ contains end subroutine score_tracklength_tally +!=============================================================================== +! SCORE_TL_ON_MESH calculate fluxes and reaction rates based on the track-length +! estimate of the flux specifically for tallies that have mesh filters. For +! these tallies, it is possible to score to multiple mesh cells for each track. +!=============================================================================== + + subroutine score_tl_on_mesh(index_tally, d_track) + + integer, intent(in) :: index_tally + real(8), intent(in) :: d_track + + integer :: i ! loop index for filter/score bins + integer :: j ! loop index for direction + integer :: k ! loop index for mesh cell crossings + integer :: ijk0(3) ! indices of starting coordinates + integer :: ijk1(3) ! indices of ending coordinates + integer :: ijk_cross(3) ! indices of mesh cell crossed + integer :: n_cross ! number of surface crossings + integer :: bins(N_FILTER_TYPES) ! scoring bin combination + integer :: score_index ! single index for single bin + integer :: score_bin ! scoring bin, e.g. SCORE_FLUX + real(8) :: flux ! tracklength estimate of flux + real(8) :: score ! actual score (e.g., flux*xs) + real(8) :: uvw(3) ! cosine of angle of particle + real(8) :: xyz0(3) ! starting/intermediate coordinates + real(8) :: xyz1(3) ! ending coordinates of particle + real(8) :: xyz_cross(3) ! coordinates of next boundary + real(8) :: d(3) ! distance to each bounding surface + real(8) :: distance ! distance traveled in mesh cell + logical :: found_bin ! was a scoring bin found? + logical :: start_in_mesh ! starting coordinates inside mesh? + logical :: end_in_mesh ! ending coordinates inside mesh? + type(TallyObject), pointer :: t => null() + type(StructuredMesh), pointer :: m => null() + + found_bin = .true. + bins = 1 + t => tallies(index_tally) + + ! ========================================================================== + ! CHECK IF THIS TRACK INTERSECTS THE MESH + + ! Copy starting and ending location of particle + xyz0 = p % coord0 % xyz - (d_track - TINY_BIT) * p % coord0 % uvw + xyz1 = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw + + ! Determine indices for starting and ending location + m => meshes(t % mesh) + call get_mesh_indices(m, xyz0, ijk0, start_in_mesh) + call get_mesh_indices(m, xyz1, ijk1, end_in_mesh) + + ! Check if start or end is in mesh -- if not, check if track still + ! intersects with mesh + if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then + if (.not. mesh_intersects(m, xyz0, xyz1)) return + end if + + ! Reset starting and ending location + xyz0 = p % coord0 % xyz - d_track * p % coord0 % uvw + xyz1 = p % coord0 % xyz + + ! ========================================================================= + ! CHECK FOR SCORING COMBINATION FOR FILTERS OTHER THAN MESH + + FILTER_LOOP: do i = 1, t % n_filters + + select case (t % filters(i)) + case (FILTER_UNIVERSE) + ! determine next universe bin + ! TODO: Account for multiple universes when performing this filter + bins(FILTER_UNIVERSE) = get_next_bin(FILTER_UNIVERSE, & + p % coord % universe, index_tally) + if (bins(FILTER_UNIVERSE) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_MATERIAL) + bins(FILTER_MATERIAL) = get_next_bin(FILTER_MATERIAL, & + p % material, index_tally) + if (bins(FILTER_MATERIAL) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_CELL) + ! determine next cell bin + ! TODO: Account for cells in multiple levels when performing this filter + bins(FILTER_CELL) = get_next_bin(FILTER_CELL, & + p % coord % cell, index_tally) + if (bins(FILTER_CELL) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_CELLBORN) + ! determine next cellborn bin + bins(FILTER_CELLBORN) = get_next_bin(FILTER_CELLBORN, & + p % cell_born, index_tally) + if (bins(FILTER_CELLBORN) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_SURFACE) + ! determine next surface bin + bins(FILTER_SURFACE) = get_next_bin(FILTER_SURFACE, & + p % surface, index_tally) + if (bins(FILTER_SURFACE) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_ENERGYIN) + ! determine incoming energy bin + k = t % n_filter_bins(FILTER_ENERGYIN) + ! check if energy of the particle is within energy bins + if (p % E < t % energy_in(1) .or. & + p % E > t % energy_in(k + 1)) then + found_bin = .false. + return + end if + + ! search to find incoming energy bin + bins(FILTER_ENERGYIN) = binary_search(t % energy_in, k + 1, p % E) + + end select + end do FILTER_LOOP + + ! ========================================================================== + ! DETERMINE WHICH MESH CELLS TO SCORE TO + + ! Calculate number of surface crossings + n_cross = sum(abs(ijk1 - ijk0)) + 1 + + ! Copy particle's direction + uvw = p % coord0 % uvw + + ! Bounding coordinates + do j = 1, 3 + if (uvw(j) > 0) then + xyz_cross(j) = m % lower_left(j) + ijk0(j) * m % width(j) + else + xyz_cross(j) = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) + end if + end do + + MESH_LOOP: do k = 1, n_cross + found_bin = .false. + + ! Calculate distance to each bounding surface. We need to treat special + ! case where the cosine of the angle is zero since this would result in a + ! divide-by-zero. + + if (k == n_cross) xyz_cross = xyz1 + + do j = 1, 3 + if (uvw(j) == 0) then + d(j) = INFINITY + else + d(j) = (xyz_cross(j) - xyz0(j))/uvw(j) + end if + end do + + ! Determine the closest bounding surface of the mesh cell by calculating + ! the minimum distance + + j = minloc(d, 1) + distance = d(j) + + ! Now use the minimum distance and diretion of the particle to determine + ! which surface was crossed + + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then + ijk_cross = ijk0 + found_bin = .true. + end if + + ! Increment indices and determine new crossing point + if (uvw(j) > 0) then + ijk0(j) = ijk0(j) + 1 + xyz_cross(j) = xyz_cross(j) + m % width(j) + else + ijk0(j) = ijk0(j) - 1 + xyz_cross(j) = xyz_cross(j) - m % width(j) + end if + + ! ======================================================================= + ! SCORE TO THIS MESH CELL + + if (found_bin) then + ! Calculate track-length estimate of flux + flux = p % wgt * distance + + ! Determine mesh bin + bins(FILTER_MESH) = mesh_indices_to_bin(m, ijk_cross) + + ! Determining scoring index + score_index = sum((bins - 1) * t % stride) + 1 + + ! Determine score for each bin + SCORE_LOOP: do i = 1, t % n_score_bins + ! determine what type of score bin + score_bin = t % score_bins(i) % scalar + + ! Determine cross section + select case(score_bin) + case (SCORE_FLUX) + score = flux + case (SCORE_TOTAL) + score = material_xs % total * flux + case (SCORE_SCATTER) + score = (material_xs % total - material_xs % absorption) * flux + case (SCORE_ABSORPTION) + score = material_xs % absorption * flux + case (SCORE_FISSION) + score = material_xs % fission * flux + case (SCORE_NU_FISSION) + score = material_xs % nu_fission * flux + end select + + ! Add score to tally + call add_to_score(t % scores(score_index, i), score) + end do SCORE_LOOP + end if + + ! Calculate new coordinates + xyz0 = xyz0 + distance * uvw + + end do MESH_LOOP + + end subroutine score_tl_on_mesh + !=============================================================================== ! GET_SCORING_BINS determines a combination of filter bins that should be scored ! for a tally based on the particle's current attributes. @@ -606,6 +862,7 @@ contains integer :: i ! loop index for filters integer :: n ! number of bins for single filter integer :: mesh_bin ! index for mesh bin + real(8) :: E ! particle energy type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() @@ -677,15 +934,22 @@ contains case (FILTER_ENERGYIN) ! determine incoming energy bin n = t % n_filter_bins(FILTER_ENERGYIN) + + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + E = p % E + else + E = p % last_E + end if + ! check if energy of the particle is within energy bins - if (p % last_E < t % energy_in(1) .or. & - p % last_E > t % energy_in(n + 1)) then + if (E < t % energy_in(1) .or. E > t % energy_in(n + 1)) then found_bin = .false. return end if ! search to find incoming energy bin - bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) + bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, E) case (FILTER_ENERGYOUT) ! determine outgoing energy bin @@ -1048,28 +1312,32 @@ contains end subroutine add_to_score !=============================================================================== -! ACCUMULATE_CYCLE_ESTIMATE +! ACCUMULATE_BATCH_ESTIMATE accumulates scores from many histories (or many +! generations) into a single realization of a random variable. !=============================================================================== - elemental subroutine accumulate_cycle_estimate(score) + elemental subroutine accumulate_batch_estimate(score) type(TallyScore), intent(inout) :: score + real(8) :: val + ! Add the sum and square of the sum of contributions from each cycle ! within a cycle to the variables sum and sum_sq. This will later allow us ! to calculate a variance on the tallies - score % sum = score % sum + score % value/n_particles - score % sum_sq = score % sum_sq + (score % value/n_particles)**2 + val = score % value/(n_particles*gen_per_batch) + score % sum = score % sum + val + score % sum_sq = score % sum_sq + val*val - ! Reset the single cycle estimate + ! Reset the single batch estimate score % value = ZERO - end subroutine accumulate_cycle_estimate + end subroutine accumulate_batch_estimate !=============================================================================== ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history -! within the cycle to a new random variable +! within the batch to a new random variable !=============================================================================== subroutine synchronize_tallies() @@ -1086,7 +1354,7 @@ contains t => tallies(i) ! Loop over all filter and scoring bins - call accumulate_cycle_estimate(t % scores) + call accumulate_batch_estimate(t % scores) end do end subroutine synchronize_tallies @@ -1183,6 +1451,7 @@ contains score_name(abs(SCORE_SCATTER_1)) = "First Scattering Moment" score_name(abs(SCORE_SCATTER_2)) = "Second Scattering Moment" score_name(abs(SCORE_SCATTER_3)) = "Third Scattering Moment" + score_name(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient" score_name(abs(SCORE_N_1N)) = "(n,1n) Rate" score_name(abs(SCORE_N_2N)) = "(n,2n) Rate" score_name(abs(SCORE_N_3N)) = "(n,3n) Rate" @@ -1522,16 +1791,21 @@ contains end function get_label !=============================================================================== -! CALCULATE_STATISTICS +! CALCULATE_STATISTICS determines the sample mean and the standard deviation of +! the mean for a TallyScore. !=============================================================================== elemental subroutine calculate_statistics(score) type(TallyScore), intent(inout) :: score - ! Calculate mean and standard deviation of the mean + ! Calculate sample mean and standard deviation of the mean -- note that we + ! have used Bessel's correction so that the estimator of the variance of the + ! sample mean is unbiased. + score % sum = score % sum/n_active - score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2)/n_active) + score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / & + (n_active - 1)) end subroutine calculate_statistics diff --git a/src/timing.F90 b/src/timing.F90 index f6e91c0583..483057125d 100644 --- a/src/timing.F90 +++ b/src/timing.F90 @@ -13,7 +13,7 @@ module timing type Timer logical :: running = .false. ! is timer running? integer :: start_counts = 0 ! counts when started - real(8) :: elapsed = 0. ! total time elapsed in seconds + real(8) :: elapsed = ZERO ! total time elapsed in seconds !!$ contains !!$ procedure :: start => timer_start !!$ procedure :: get_value => timer_get_value @@ -24,7 +24,7 @@ module timing contains !=============================================================================== -! TIMER_START +! TIMER_START starts running a timer and measures the current time !=============================================================================== subroutine timer_start(self) @@ -38,7 +38,7 @@ contains end subroutine timer_start !=============================================================================== -! TIMER_GET_VALUE +! TIMER_GET_VALUE returns the current value of the timer !=============================================================================== function timer_get_value(self) result(elapsed) @@ -61,7 +61,7 @@ contains end function timer_get_value !=============================================================================== -! TIMER_STOP +! TIMER_STOP stops the timer and sets the elapsed time !=============================================================================== subroutine timer_stop(self) @@ -78,7 +78,7 @@ contains end subroutine timer_stop !=============================================================================== -! TIMER_RESET +! TIMER_RESET resets a timer to have a zero value !=============================================================================== subroutine timer_reset(self) diff --git a/src/utils/read_source.py b/src/utils/read_source.py new file mode 100755 index 0000000000..a7cb801c7f --- /dev/null +++ b/src/utils/read_source.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python + +import struct + +class SourceFile(object): + def __init__(self, filename): + # Open source file for reading + self.f = open(filename, 'r') + + # Read number of source sites + self.n_sites = struct.unpack('q', self.f.read(8))[0] + + # Create list to store source sites + self.sites = [] + + for i in range(self.n_sites): + # Read position, angle, and energy + uid = struct.unpack('q', self.f.read(8)) + xyz = struct.unpack('3d', self.f.read(24)) + uvw = struct.unpack('3d', self.f.read(24)) + E = struct.unpack('d', self.f.read(8)) + + # Create source site and append to list + self.sites.append(SourceSite(xyz, uvw, E)) + +class SourceSite(object): + def __init__(self, xyz, uvw, E): + self.xyz = xyz + self.uvw = uvw + self.E = E + + def __repr__(self): + return "".format(self.xyz, self.E) diff --git a/src/xml-fortran/templates/plots_t.xml b/src/xml-fortran/templates/plots_t.xml index af8f2c6e27..b3b86fc44e 100644 --- a/src/xml-fortran/templates/plots_t.xml +++ b/src/xml-fortran/templates/plots_t.xml @@ -3,15 +3,28 @@ + + + + + + + + + + + - + + + diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 6061343067..961397d443 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -4,10 +4,11 @@ - + + @@ -27,13 +28,16 @@ - + + + + + + - - - + - +