diff --git a/docs/source/publications.rst b/docs/source/publications.rst index 9d3491b8ab..99a00d13f2 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -4,6 +4,10 @@ Publications ============ +- Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "Accelerated sampling + of the free gas resonance elastic scattering kernel," *Ann. Nucl. Energy*, + **69**, 116--124 (2014). ``_ + - Benoit Forget, Sheng Xu, and Kord Smith, "Direct Doppler broadening in Monte Carlo simulations using the multipole representation," *Ann. Nucl. Energy*, **64**, 78--85 (2014). ``_ diff --git a/docs/source/releasenotes/index.rst b/docs/source/releasenotes/index.rst index b40126b243..9799a5cfc7 100644 --- a/docs/source/releasenotes/index.rst +++ b/docs/source/releasenotes/index.rst @@ -10,6 +10,7 @@ bugs fixed, and known issues for each successive release. .. toctree:: :maxdepth: 1 + notes_0.6.1 notes_0.6.0 notes_0.5.4 notes_0.5.3 diff --git a/docs/source/releasenotes/notes_0.6.1.rst b/docs/source/releasenotes/notes_0.6.1.rst new file mode 100644 index 0000000000..b0651a8629 --- /dev/null +++ b/docs/source/releasenotes/notes_0.6.1.rst @@ -0,0 +1,63 @@ +.. _notes_0.6.1: + +============================== +Release Notes for OpenMC 0.6.1 +============================== + +------------------- +System Requirements +------------------- + +There are no special requirements for running the OpenMC code. As of this +release, OpenMC has been tested on a variety of Linux distributions, Mac OS X, +and Microsoft Windows 7. Memory requirements will vary depending on the size of +the problem at hand (mostly on the number of nuclides in the problem). + +------------ +New Features +------------ + +- Coarse mesh finite difference acceleration no longer requires PETSc +- Statepoint file numbering is now zero-padded +- Python scripts now compatible with Python 2 or 3 +- Ability to run particle restarts in fixed source calculations +- Capability to filter box source by fissionable materials +- Nuclide/element names are now case insensitive in input files +- Improved treatment of resonance scattering for heavy nuclides + +--------- +Bug Fixes +--------- + +- 03e890_: Check for energy-dependent multiplicities in ACE files +- 4439de_: Fix distance-to-surface calculation for general plane surface +- 5808ed_: Account for differences in URR band probabilities at different energies +- 2e60c0_: Allow zero atom/weight percents in materials +- 3e0870_: Don't use PWD environment variable when setting path to input files +- dc4776_: Handle probability table resampling correctly +- 01178b_: Fix metastables nuclides in NNDC cross_sections.xml file +- 62ec43_: Don't read tallies.xml when OpenMC is run in plotting mode +- 2a95ef_: Prevent segmentation fault on "current" score without mesh filter + +.. _03e890: https://github.com/mit-crpg/openmc/commit/03e890 +.. _4439de: https://github.com/mit-crpg/openmc/commit/4439de +.. _5808ed: https://github.com/mit-crpg/openmc/commit/5808ed +.. _2e60c0: https://github.com/mit-crpg/openmc/commit/2e60c0 +.. _3e0870: https://github.com/mit-crpg/openmc/commit/3e0870 +.. _dc4776: https://github.com/mit-crpg/openmc/commit/dc4776 +.. _01178b: https://github.com/mit-crpg/openmc/commit/01178b +.. _62ec43: https://github.com/mit-crpg/openmc/commit/62ec43 +.. _2a95ef: https://github.com/mit-crpg/openmc/commit/2a95ef + +------------ +Contributors +------------ + +This release contains new contributions from the following people: + +- `Sterling Harper `_ +- `Bryan Herman `_ +- `Adam Nelson `_ +- `Paul Romano `_ +- `Jon Walsh `_ +- `Will Boyd `_ diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 09caa43b5c..4552860787 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -260,6 +260,65 @@ or sub-elements and can be set to either "false" or "true". *Default*: true +```` Element +---------------------------------- + +The ``resonance_scattering`` element can contain one or more of the following +attributes or sub-elements: + + :scatterer: + An element with attributes/sub-elements called ``nuclide``, ``method``, + ``xs_label``, ``xs_label_0K``, ``E_min``, and ``E_max``. The ``nuclide`` + attribute is the name, as given by the ``name`` attribute within the + ``nuclide`` sub-element of the ``material`` element in ``materials.xml``, + of the nuclide to which a resonance scattering treatment is to be applied. + The ``method`` attribute gives the type of resonance scattering treatment + that is to be applied to the ``nuclide``. Acceptable inputs - none of + which are case-sensitive - for the ``method`` attribute are ``ARES``, + ``CXS``, ``WCM``, and ``DBRC``. Descriptions of each of these methods + are documented here_. The ``xs_label`` attribute gives the label for the + cross section data of the ``nuclide`` at a given temperature. The + ``xs_label_0K`` gives the label for the 0 K cross section data for the + ``nuclide``. The ``E_min`` attribute gives the minimum energy above + which the ``method`` is applied. The ``E_max`` attribute gives the + maximum energy below which the ``method`` is applied. One example would + be as follows: + + .. _here: http://dx.doi.org/10.1016/j.anucene.2014.01.017 + + .. code-block:: xml + + + + U-238 + ARES + 92238.72c + 92238.00c + 5.0e-6 + 40.0e-6 + + + Pu-239 + dbrc + 94239.72c + 94239.00c + 0.01e-6 + 210.0e-6 + + + + .. note:: If the ``resonance_scattering`` element is not given, the free gas, + constant cross section (``cxs``) scattering model, which has + historically been used by Monte Carlo codes to sample target + velocities, is used to treat the target motion of all nuclides. If + ``resonance_scattering`` is present, the ``cxs`` method is applied + below ``E_min`` and the target-at-rest (asymptotic) kernel is used + above ``E_max``. An arbitrary number of ``scatterer`` elements may + be specified, each corresponding to a single nuclide at a single + temperature. + + *Defaults*: None (scatterer), ARES (method), 0.01 eV (E_min), 1.0 keV (E_max) + ```` Element ---------------------- @@ -311,6 +370,12 @@ attributes/sub-elements: parallelepiped and the last three of which specify the upper-right corner. Source sites are sampled uniformly through that parallelepiped. + To filter a "box" spatial distribution by fissionable material, specify + "fission" tag instead of "box". The ``parameters`` should be given as six + real numbers, the first three of which specify the lower-left corner of a + parallelepiped and the last three of which specify the upper-right + corner. Source sites are sampled uniformly through that parallelepiped. + For a "point" spatial distribution, ``parameters`` should be given as three real numbers which specify the (x,y,z) location of an isotropic point source @@ -368,6 +433,13 @@ attributes/sub-elements: *Default*: 0.988 2.249 + :write_initial: + An element specifying whether to write out the initial source bank used at + the beginning of the first batch. The output file is named + "initial_source.binary(h5)" + + *Default*: false + ```` Element ------------------------- @@ -1313,15 +1385,6 @@ Currently, it allows users to accelerate fission source convergence during inactive neutron batches. To run CMFD, the ```` element in ``settings.xml`` should be set to "true". -```` Element --------------------------- - -The ```` element controls the batch where CMFD tallies should be -reset. CMFD tallies should be reset before active batches so they are accumulated -without bias. - - *Default*: 0 - ```` Element ------------------- @@ -1343,7 +1406,25 @@ The ```` element sets one additional CMFD output column. Options are: * "source" - prints the RMS [%] between the OpenMC fission source and CMFD fission source. - *Default*: None + *Default*: balance + +```` Element +------------------------ + +The ```` element controls whether :math:`\widehat{D}` nonlinear +CMFD parameters should be reset to zero before solving CMFD eigenproblem. +It can be turned on with "true" and off with "false". + + *Default*: false + +```` Element +------------------------- + +The ```` element controls whether an effective downscatter cross +section should be used when using 2-group CMFD. It can be turned on with "true" +and off with "false". + + *Default*: false ```` Element ---------------------- @@ -1354,24 +1435,16 @@ It can be turned on with "true" and off with "false". *Default*: false -```` Element ----------------------- +```` Element +-------------------- -The ```` element controls if cmfd tallies should be accumulated -during inactive batches. For some applications, CMFD tallies may not be -needed until the start of active batches. This option can be turned on -with "true" and off with "false" +The ```` element specifies two parameters. The first is +the absolute inner tolerance for Gauss-Seidel iterations when performing CMFD +and the second is the relative inner tolerance for Gauss-Seidel iterations +for CMFD calculations. It is only used in the standalone CMFD power iteration +solver and not when PETSc is active. - *Default*: true - -```` Element ----------------------------- - -The ```` element controls when CMFD tallies are reset during -inactive batches. The integer set here is the interval at which this reset -occurs. The amout of resets is controlled with the ```` element. - - *Defualt*: 9999 + *Default*: 1.e-10 1.e-5 ```` Element ------------------------- @@ -1380,9 +1453,16 @@ The ```` element is used to view the convergence of linear GMRES iterations in PETSc. This option can be turned on with "true" and turned off with "false". - *Default*: false +```` Element +-------------------- + +The ```` element specifies the tolerance on the eigenvalue when performing +CMFD power iteration. + + *Default*: 1.e-8 + ```` Element ------------------ @@ -1451,14 +1531,6 @@ not impact the calculation. *Default*: 1.0 -```` Element -------------------------- - -The ```` element controls the number of CMFD tally resets that -occur during inactive CMFD batches. - - *Default*: 9999 - ```` Element --------------------------- @@ -1471,16 +1543,8 @@ This option can be turned on with "true" and turned off with "false". ------------------------- The ```` element can be turned on with "true" to have an adjoint -calculation be performed on the last batch when CMFD is active. - - *Default*: false - -```` Element --------------------------- - -The ```` element is used to view the convergence of the nonlinear SNES -function in PETSc. This option can be turned on with "true" and turned off with "false". - +calculation be performed on the last batch when CMFD is active. OpenMC should be +compiled with PETSc when using this option. *Default*: false @@ -1493,6 +1557,41 @@ By setting "power", power iteration is used and by setting "jfnk", JFNK is used. *Default*: power +```` Element +-------------------- + +The ```` element specifies an optional Wielandt shift parameter for +accelerating power iterations. It can only be used when PETSc is not active. +It is by default very large so the impact of the shift is effectively zero. + + *Default*: 1e6 + +```` Element +---------------------- + +The ```` element specifies an optional spectral radius that can be set to +accelerate the convergence of Gauss-Seidel iterations during CMFD power iteration +solve. Note this is only used in the standalone CMFD solver and does not affect +the calculation when PETSc is active. + + *Default*: power + +```` Element +------------------ + +The ```` element specifies the tolerance on the fission source when performing +CMFD power iteration. + + *Default*: 1.e-8 + +```` Element +------------------------- + +The ```` element contains a list of batch numbers in which CMFD tallies +should be reset. + + *Default*: None + ```` Element ---------------------------- diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fc6c6630cc..8a9424e79d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -277,7 +277,7 @@ file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_*.py) # Check to see if PETSC is compiled for CMFD tests if (NOT ${PETSC_ENABLED}) - file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd*.py) + file(GLOB_RECURSE CMFD_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/../tests/test_cmfd_jfnk.py) foreach(cmfd_test in ${CMFD_TESTS}) list(REMOVE_ITEM TESTS ${cmfd_test}) endforeach(cmfd_test) diff --git a/src/ace.F90 b/src/ace.F90 index 9c088693de..6a712ccd74 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -7,11 +7,11 @@ module ace use error, only: fatal_error, warning use fission, only: nu_total use global - use list_header, only: ListElemInt, ListInt + use list_header, only: ListInt use material_header, only: Material use output, only: write_message use set_header, only: SetChar - use string, only: to_str + use string, only: to_str, to_lower implicit none @@ -36,6 +36,7 @@ contains integer :: i ! index in materials array integer :: j ! index over nuclides in material integer :: k ! index over S(a,b) tables in material + integer :: n ! index over resonant scatterers integer :: i_listing ! index in xs_listings array integer :: i_nuclide ! index in nuclides integer :: i_sab ! index in sab_tables @@ -67,8 +68,8 @@ contains name = mat % names(j) if (.not. already_read % contains(name)) then - i_listing = xs_listing_dict % get_key(name) - i_nuclide = nuclide_dict % get_key(name) + i_listing = xs_listing_dict % get_key(to_lower(name)) + i_nuclide = nuclide_dict % get_key(to_lower(name)) name = xs_listings(i_listing) % name alias = xs_listings(i_listing) % alias @@ -80,6 +81,30 @@ contains ! array call read_ace_table(i_nuclide, i_listing) + ! 0K resonant scatterer information, if treating resonance scattering + if (treat_res_scat) then + do n = 1, n_res_scatterers_total + if (name == nuclides_0K(n) % name) then + nuclides(i_nuclide) % resonant = .true. + nuclides(i_nuclide) % name_0K = nuclides_0K(n) % name_0K + nuclides(i_nuclide) % name_0K = trim(nuclides(i_nuclide) % & + & name_0K) + nuclides(i_nuclide) % scheme = nuclides_0K(n) % scheme + nuclides(i_nuclide) % scheme = trim(nuclides(i_nuclide) % & + & scheme) + nuclides(i_nuclide) % E_min = nuclides_0K(n) % E_min + nuclides(i_nuclide) % E_max = nuclides_0K(n) % E_max + if (.not. already_read % contains(nuclides(i_nuclide) % & + & name_0K)) then + i_listing = xs_listing_dict % get_key(nuclides(i_nuclide) % & + & name_0K) + call read_ace_table(i_nuclide, i_listing) + end if + exit + end if + end do + end if + ! Add name and alias to dictionary call already_read % add(name) call already_read % add(alias) @@ -91,8 +116,8 @@ contains name = mat % sab_names(k) if (.not. already_read % contains(name)) then - i_listing = xs_listing_dict % get_key(name) - i_sab = sab_dict % get_key(name) + i_listing = xs_listing_dict % get_key(to_lower(name)) + i_sab = sab_dict % get_key(to_lower(name)) ! Read the ACE table into the appropriate entry on the sab_tables ! array @@ -172,6 +197,25 @@ contains ! Avoid some valgrind leak errors call already_read % clear() + ! Loop around material + MATERIAL_LOOP3: do i = 1, n_materials + + ! Get material + mat => materials(i) + + ! Loop around nuclides in material + NUCLIDE_LOOP2: do j = 1, mat % n_nuclides + + ! Check for fission in nuclide + if (nuclides(mat % nuclide(j)) % fissionable) then + mat % fissionable = .true. + exit NUCLIDE_LOOP2 + end if + + end do NUCLIDE_LOOP2 + + end do MATERIAL_LOOP3 + end subroutine read_xs !=============================================================================== @@ -198,6 +242,7 @@ contains real(8) :: awrs(16) ! list of atomic weight ratios (not used) real(8) :: awr ! atomic weight ratio for table logical :: file_exists ! does ACE library exist? + logical :: data_0K ! are we reading 0K data? character(7) :: readable ! is ACE library readable? character(10) :: name ! name of ACE table character(10) :: date_ ! date ACE library was processed @@ -301,19 +346,32 @@ contains select case(listing % type) case (ACE_NEUTRON) + + ! only read in a resonant scatterers info once nuc => nuclides(i_table) - nuc % name = name - nuc % awr = awr - nuc % kT = kT - nuc % zaid = NXS(2) + data_0K = .false. + if (trim(adjustl(name)) == nuc % name_0K) then + data_0K = .true. + else + nuc % name = name + nuc % awr = awr + nuc % kT = kT + nuc % zaid = NXS(2) + end if ! read all blocks - call read_esz(nuc) - call read_nu_data(nuc) - call read_reactions(nuc) - call read_angular_dist(nuc) - call read_energy_dist(nuc) - call read_unr_res(nuc) + call read_esz(nuc, data_0K) + + ! don't read unnecessary 0K data for resonant scatterers + if (data_0K) then + continue + else + call read_nu_data(nuc) + call read_reactions(nuc) + call read_angular_dist(nuc) + call read_energy_dist(nuc) + call read_unr_res(nuc) + end if ! Currently subcritical fixed source calculations are not allowed. Thus, ! if any fissionable material is found in a fixed source calculation, @@ -325,9 +383,12 @@ contains ! for fissionable nuclides, precalculate microscopic nu-fission cross ! sections so that we don't need to call the nu_total function during - ! cross section lookups + ! cross section lookups (except if we're dealing w/ 0K data for resonant + ! scatterers) - if (nuc % fissionable) call generate_nu_fission(nuc) + if (nuc % fissionable .and. .not. data_0K) then + call generate_nu_fission(nuc) + end if case (ACE_THERMAL) sab => sab_tables(i_table) @@ -358,43 +419,82 @@ contains ! total xs, absorption xs, elastic scattering xs, and heating numbers. !=============================================================================== - subroutine read_esz(nuc) + subroutine read_esz(nuc, data_0K) type(Nuclide), pointer :: nuc + logical :: data_0K ! are we reading 0K data? + integer :: NE ! number of energy points for total and elastic cross sections + integer :: i ! index in 0K elastic xs array for this nuclide + + real(8) :: xs_cdf_sum = ZERO ! xs cdf value ! determine number of energy points NE = NXS(3) - nuc % n_grid = NE ! allocate storage for energy grid and cross section arrays - allocate(nuc % energy(NE)) - allocate(nuc % total(NE)) - allocate(nuc % elastic(NE)) - allocate(nuc % fission(NE)) - allocate(nuc % nu_fission(NE)) - allocate(nuc % absorption(NE)) - ! initialize cross sections - nuc % total = ZERO - nuc % elastic = ZERO - nuc % fission = ZERO - nuc % nu_fission = ZERO - nuc % absorption = ZERO + ! read in 0K data if we've already read in non-0K data + if (data_0K) then + nuc % n_grid_0K = NE + allocate(nuc % energy_0K(NE)) + allocate(nuc % elastic_0K(NE)) + allocate(nuc % xs_cdf(NE)) + nuc % elastic_0K = ZERO + nuc % xs_cdf = ZERO + XSS_index = 1 + nuc % energy_0K = get_real(NE) - ! Read data from XSS -- only the energy grid, elastic scattering and heating - ! cross section values are actually read from here. The total and absorption - ! cross sections are reconstructed from the partial reaction data. + ! Skip total and absorption + XSS_index = XSS_index + 2*NE + + ! Continue reading elastic scattering and heating + nuc % elastic_0K = get_real(NE) - XSS_index = 1 - nuc % energy = get_real(NE) + do i = 1, nuc % n_grid_0K - 1 - ! Skip total and absorption - XSS_index = XSS_index + 2*NE + ! Negative cross sections result in a CDF that is not monotonically + ! increasing. Set all negative xs values to ZERO. + if (nuc % elastic_0K(i) < ZERO) nuc % elastic_0K(i) = ZERO - ! Continue reading elastic scattering and heating - nuc % elastic = get_real(NE) + ! build xs cdf + xs_cdf_sum = xs_cdf_sum + (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) & + & + sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO & + & * (nuc % energy_0K(i+1) - nuc % energy_0K(i)) + nuc % xs_cdf(i) = xs_cdf_sum + end do + + else ! read in non-0K data + nuc % n_grid = NE + allocate(nuc % energy(NE)) + allocate(nuc % total(NE)) + allocate(nuc % elastic(NE)) + allocate(nuc % fission(NE)) + allocate(nuc % nu_fission(NE)) + allocate(nuc % absorption(NE)) + + ! initialize cross sections + nuc % total = ZERO + nuc % elastic = ZERO + nuc % fission = ZERO + nuc % nu_fission = ZERO + nuc % absorption = ZERO + + ! Read data from XSS -- only the energy grid, elastic scattering and heating + ! cross section values are actually read from here. The total and absorption + ! cross sections are reconstructed from the partial reaction data. + + XSS_index = 1 + nuc % energy = get_real(NE) + + ! Skip total and absorption + XSS_index = XSS_index + 2*NE + + ! Continue reading elastic scattering and heating + nuc % elastic = get_real(NE) + + end if end subroutine read_esz @@ -1463,16 +1563,11 @@ contains integer :: i ! index in nuclides array integer :: j ! index in nuclides array - type(ListElemInt), pointer :: nuc_list => null() ! pointer to nuclide list do i = 1, n_nuclides_total - allocate(nuclides(i) % nuc_list) - nuc_list => nuclides(i) % nuc_list do j = 1, n_nuclides_total if (nuclides(i) % zaid == nuclides(j) % zaid) then - nuc_list % data = j - allocate(nuc_list % next) - nuc_list => nuc_list % next + call nuclides(i) % nuc_list % append(j) end if end do end do diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 2067030c1e..76492dfe80 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -2,7 +2,7 @@ module ace_header use constants, only: MAX_FILE_LEN use endf_header, only: Tab1 - use list_header, only: ListElemInt + use list_header, only: ListInt implicit none @@ -97,7 +97,7 @@ module ace_header real(8) :: kT ! temperature in MeV (k*T) ! Linked list of indices in nuclides array of instances of this same nuclide - type(ListElemInt), pointer :: nuc_list => null() + type(ListInt) :: nuc_list ! Energy grid information integer :: n_grid ! # of nuclide grid points @@ -112,6 +112,17 @@ module ace_header real(8), allocatable :: absorption(:) ! absorption (MT > 100) real(8), allocatable :: heating(:) ! heating + ! Resonance scattering info + logical :: resonant = .false. ! resonant scatterer? + character(10) :: name_0K = '' ! name of 0K nuclide, e.g. 92235.00c + character(16) :: scheme ! target velocity sampling scheme + integer :: n_grid_0K ! number of 0K energy grid points + real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs + real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section + real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section + real(8) :: E_min ! lower cutoff energy for res scattering + real(8) :: E_max ! upper cutoff energy for res scattering + ! Fission information logical :: fissionable ! nuclide is fissionable? logical :: has_partial_fission ! nuclide has partial fission reactions? @@ -147,6 +158,22 @@ module ace_header procedure :: clear => nuclide_clear ! Deallocates Nuclide end type Nuclide +!=============================================================================== +! NUCLIDE0K temporarily contains all 0K cross section data and other parameters +! needed to treat resonance scattering before transferring them to NUCLIDE +!=============================================================================== + + type Nuclide0K + + character(10) :: nuclide ! name of nuclide, e.g. U-238 + character(16) :: scheme = 'ares' ! target velocity sampling scheme + character(10) :: name ! name of nuclide, e.g. 92235.03c + character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c + real(8) :: E_min = 0.01e-6 ! lower cutoff energy for res scattering + real(8) :: E_max = 1000.0e-6 ! upper cutoff energy for res scattering + + end type Nuclide0K + !=============================================================================== ! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic ! thermal scattering collisions which utilize a continuous secondary energy @@ -342,12 +369,24 @@ module ace_header integer :: i ! Loop counter - if (allocated(this % grid_index)) deallocate(this % grid_index) + if (allocated(this % grid_index)) & + deallocate(this % grid_index) if (allocated(this % energy)) & deallocate(this % energy, this % total, this % elastic, & - this % fission, this % nu_fission, this % absorption) - if (allocated(this % heating)) deallocate(this % heating) + & this % fission, this % nu_fission, this % absorption) + + if (allocated(this % energy_0K)) & + deallocate(this % energy_0K) + + if (allocated(this % elastic_0K)) & + deallocate(this % elastic_0K) + + if (allocated(this % xs_cdf)) & + deallocate(this % xs_cdf) + + if (allocated(this % heating)) & + deallocate(this % heating) if (allocated(this % index_fission)) deallocate(this % index_fission) @@ -377,6 +416,8 @@ module ace_header deallocate(this % reactions) end if + call this % nuc_list % clear() + end subroutine nuclide_clear end module ace_header diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index b5612dcd05..abba2d47d7 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -10,8 +10,6 @@ module cmfd_data private public :: set_up_cmfd, neutron_balance - logical :: dhat_reset = .false. - contains !============================================================================== @@ -103,6 +101,8 @@ contains cmfd % hxyz(2,:,:,:) = m % width(2) ! set y width cmfd % hxyz(3,:,:,:) = m % width(3) ! set z width + cmfd % keff_bal = ZERO + ! Begin loop around tallies TAL: do ital = 1, n_cmfd_tallies @@ -211,6 +211,9 @@ contains ! Bank source cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + & t % results(2,score_index) % sum + cmfd % keff_bal = cmfd % keff_bal + & + t % results(2,score_index) % sum / & + dble(t % n_realizations) end do INGROUP @@ -623,7 +626,9 @@ contains subroutine compute_dhat() use constants, only: CMFD_NOACCEL, ZERO - use global, only: cmfd, cmfd_coremap + use global, only: cmfd, cmfd_coremap, message, dhat_reset + use output, only: write_message + use string, only: to_str integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction @@ -743,7 +748,9 @@ contains cmfd%dhat(l,g,i,j,k) = dhat ! check for dhat reset - if (dhat_reset) cmfd%dhat(l,g,i,j,k) = ZERO + if (dhat_reset) then + cmfd%dhat(l,g,i,j,k) = ZERO + end if end do LEAK @@ -755,6 +762,12 @@ contains end do ZLOOP + ! write that dhats are zero + if (dhat_reset) then + message = 'Dhats reset to zero.' + call write_message(1) + end if + end subroutine compute_dhat !=============================================================================== @@ -763,8 +776,8 @@ contains function get_reflector_albedo(l, g, i, j, k) - use constants, only: ALBEDO_REJECT - use global, only: cmfd, cmfd_hold_weights + use constants, only: ONE + use global, only: cmfd real(8) :: get_reflector_albedo ! reflector albedo integer, intent(in) :: i ! iteration counter for x @@ -786,8 +799,7 @@ contains ! Calculate albedo if ((shift_idx == 1 .and. current(2*l ) < 1.0e-10_8) .or. & (shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then - albedo = ALBEDO_REJECT - cmfd_hold_weights = .true. + albedo = ONE else albedo = (current(2*l-1)/current(2*l))**(shift_idx) end if @@ -797,137 +809,6 @@ contains end function get_reflector_albedo -!=============================================================================== -! FIX_NEUTRON_BALANCE is a method to adjust parameters to have perfect balance -!=============================================================================== -#ifdef DEVELOPMENTAL - subroutine fix_neutron_balance() - - use constants, only: ONE, ZERO, CMFD_NOACCEL - use global, only: cmfd, keff - use, intrinsic :: ISO_FORTRAN_ENV - - integer :: nx ! number of mesh cells in x direction - integer :: ny ! number of mesh cells in y direction - integer :: nz ! number of mesh cells in z direction - integer :: ng ! number of energy groups - integer :: i ! iteration counter for x - integer :: j ! iteration counter for y - integer :: k ! iteration counter for z - integer :: l ! iteration counter for surface - real(8) :: leak1 ! leakage rate in group 1 - real(8) :: leak2 ! leakage rate in group 2 - real(8) :: flux1 ! group 1 volume int flux - real(8) :: flux2 ! group 2 volume int flux - real(8) :: sigt1 ! group 1 total xs - real(8) :: sigt2 ! group 2 total xs - real(8) :: sigs11 ! scattering transfer 1 --> 1 - real(8) :: sigs21 ! scattering transfer 2 --> 1 - real(8) :: sigs12 ! scattering transfer 1 --> 2 - real(8) :: sigs22 ! scattering transfer 2 --> 2 - real(8) :: nsigf11 ! fission transfer 1 --> 1 - real(8) :: nsigf21 ! fission transfer 2 --> 1 - real(8) :: nsigf12 ! fission transfer 1 --> 2 - real(8) :: nsigf22 ! fission transfer 2 --> 2 - real(8) :: siga1 ! group 1 abs xs - real(8) :: siga2 ! group 2 abs xs - real(8) :: sigs12_eff ! effective downscatter xs - - ! Extract spatial and energy indices from object - nx = cmfd % indices(1) - ny = cmfd % indices(2) - nz = cmfd % indices(3) - ng = cmfd % indices(4) - - ! Return if not two groups - if (ng /= 2) return - - ! Begin loop around space and energy groups - ZLOOP: do k = 1, nz - - YLOOP: do j = 1, ny - - XLOOP: do i = 1, nx - - ! Check for active mesh - if (allocated(cmfd%coremap)) then - if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) cycle - end if - - ! Compute leakage in groups 1 and 2 - leak1 = ZERO - leak2 = ZERO - LEAK: do l = 1, 3 - - leak1 = leak1 + ((cmfd % current(4*l,1,i,j,k) - & - cmfd % current(4*l-1,1,i,j,k))) - & - ((cmfd % current(4*l-2,1,i,j,k) - & - cmfd % current(4*l-3,1,i,j,k))) - - leak2 = leak2 + ((cmfd % current(4*l,2,i,j,k) - & - cmfd % current(4*l-1,2,i,j,k))) - & - ((cmfd % current(4*l-2,2,i,j,k) - & - cmfd % current(4*l-3,2,i,j,k))) - - - end do LEAK - - ! Extract cross sections and flux from object - flux1 = cmfd % flux(1,i,j,k) - flux2 = cmfd % flux(2,i,j,k) - sigt1 = cmfd % totalxs(1,i,j,k) - sigt2 = cmfd % totalxs(2,i,j,k) - sigs11 = cmfd % scattxs(1,1,i,j,k) - sigs21 = cmfd % scattxs(2,1,i,j,k) - sigs12 = cmfd % scattxs(1,2,i,j,k) - sigs22 = cmfd % scattxs(2,2,i,j,k) - nsigf11 = cmfd % nfissxs(1,1,i,j,k) - nsigf21 = cmfd % nfissxs(2,1,i,j,k) - nsigf12 = cmfd % nfissxs(1,2,i,j,k) - nsigf22 = cmfd % nfissxs(2,2,i,j,k) - - ! Check for no fission into group 2 - if (.not.(nsigf12 < 1e-6_8 .and. nsigf22 < 1e-6_8)) then - write(OUTPUT_UNIT,'(A,1PE11.4,1X,1PE11.4)') 'Fission in G=2', & - nsigf12,nsigf22 - end if - - ! Compute absorption xs - siga1 = sigt1 - sigs11 - sigs12 - siga2 = sigt2 - sigs22 - sigs21 - - ! Compute effective downscatter xs - sigs12_eff = (ONE/keff*nsigf11*flux1 - leak1 - siga1*flux1 & - - ONE/keff*nsigf21/siga2*leak2 ) / ( flux1*(ONE & - - ONE/keff*nsigf21/siga2)) - - ! Redefine flux 2 - flux2 = (sigs12_eff*flux1 - leak2)/siga2 - cmfd % flux(2,i,j,k) = flux2 - - ! Recompute total cross sections (use effective and no upscattering) - sigt1 = siga1 + sigs11 + sigs12_eff - sigt2 = siga2 + sigs22 - - ! Record total xs - cmfd % totalxs(1,i,j,k) = sigt1 - cmfd % totalxs(2,i,j,k) = sigt2 - - ! Record effective downscatter xs - cmfd % scattxs(1,2,i,j,k) = sigs12_eff - - ! Zero out upscatter cross section - cmfd % scattxs(2,1,i,j,k) = ZERO - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - end subroutine fix_neutron_balance -#endif - !=============================================================================== ! COMPUTE_EFFECTIVE_DOWNSCATTER changes downscatter rate for zero upscatter !=============================================================================== diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 6234d9fcf9..d14a3af243 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -22,6 +22,7 @@ contains use cmfd_data, only: set_up_cmfd use cmfd_power_solver, only: cmfd_power_execute use cmfd_jfnk_solver, only: cmfd_jfnk_execute + use cmfd_solver, only: cmfd_solver_execute use error, only: warning, fatal_error ! CMFD single processor on master @@ -37,6 +38,7 @@ contains call process_cmfd_options() ! Call solver +#ifdef PETSC if (trim(cmfd_solver_type) == 'power') then call cmfd_power_execute() elseif (trim(cmfd_solver_type) == 'jfnk') then @@ -45,6 +47,9 @@ contains message = 'solver type became invalid after input processing' call fatal_error() end if +#else + call cmfd_solver_execute() +#endif ! Save k-effective cmfd % k_cmfd(current_batch) = cmfd % keff @@ -64,7 +69,7 @@ contains call calc_fission_source() ! calculate weight factors - if (cmfd_feedback) call cmfd_reweight(.true.) + call cmfd_reweight(.true.) ! stop cmfd timer if (master) call time_cmfd % stop() @@ -77,36 +82,23 @@ contains subroutine cmfd_init_batch() - use global, only: cmfd_begin, cmfd_on, cmfd_tally_on, & - cmfd_inact_flush, cmfd_act_flush, cmfd_run, & - current_batch, cmfd_hold_weights + use global, only: cmfd_begin, cmfd_on, & + cmfd_reset, cmfd_run, & + current_batch ! Check to activate CMFD diffusion and possible feedback ! this guarantees that when cmfd begins at least one batch of tallies are ! accumulated if (cmfd_run .and. cmfd_begin == current_batch) then cmfd_on = .true. - cmfd_tally_on = .true. end if ! If this is a restart run and we are just replaying batches leave if (restart_run .and. current_batch <= restart_batch) return - ! Check to flush cmfd tallies for active batches, no more inactive flush - if (cmfd_run .and. cmfd_act_flush == current_batch) then + ! Check to reset tallies + if (cmfd_run .and. cmfd_reset % contains(current_batch)) then call cmfd_tally_reset() - cmfd_tally_on = .true. - cmfd_inact_flush(2) = -1 - end if - - ! Check to flush cmfd tallies during inactive batches (>= on number of - ! flushes important as the code will flush on the first batch which we - ! dont want to count) - if (cmfd_run .and. mod(current_batch,cmfd_inact_flush(1)) & - == 0 .and. cmfd_inact_flush(2) > 0 .and. cmfd_begin < current_batch) then - cmfd_hold_weights = .true. - call cmfd_tally_reset() - cmfd_inact_flush(2) = cmfd_inact_flush(2) - 1 end if end subroutine cmfd_init_batch @@ -139,6 +131,7 @@ contains use constants, only: CMFD_NOACCEL, ZERO, TWO use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch + use string, only: to_str #ifdef MPI use global, only: mpi_err @@ -267,6 +260,7 @@ contains use mesh_header, only: StructuredMesh use mesh, only: count_bank_sites, get_mesh_indices use search, only: binary_search + use string, only: to_str #ifdef MPI use global, only: mpi_err @@ -287,7 +281,6 @@ contains logical :: in_mesh ! source site is inside mesh type(StructuredMesh), pointer :: m ! point to mesh - real(8), allocatable :: egrid(:) ! energy grid ! Associate pointer m => meshes(n_user_meshes + 1) @@ -308,19 +301,16 @@ contains cmfd % weightfactors = ONE end if - ! Allocate energy grid and reverse cmfd energy grid - if (.not. allocated(egrid)) allocate(egrid(ng + 1)) - egrid = (/(cmfd % egrid(ng - i + 2), i = 1, ng + 1)/) - ! Compute new weight factors if (new_weights) then - ! Zero out weights - cmfd%weightfactors = ZERO + ! Set weight factors to a default 1.0 + cmfd%weightfactors = ONE - ! Count bank sites in mesh - call count_bank_sites(m, source_bank, cmfd%sourcecounts, egrid, & + ! Count bank sites in mesh and reverse due to egrid structure + call count_bank_sites(m, source_bank, cmfd%sourcecounts, cmfd % egrid, & sites_outside=outside, size_bank=work) + cmfd % sourcecounts = cmfd%sourcecounts(ng:1:-1,:,:,:) ! Check for sites outside of the mesh if (master .and. outside) then @@ -336,12 +326,14 @@ contains end where end if + if (.not. cmfd_feedback) return + ! Broadcast weight factors to all procs #ifdef MPI call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, & MPI_COMM_WORLD, mpi_err) #endif - end if + end if ! begin loop over source bank do i = 1, int(work,4) @@ -378,9 +370,6 @@ contains end do - ! Deallocate all - if (allocated(egrid)) deallocate(egrid) - end subroutine cmfd_reweight !=============================================================================== diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index 73bc0e0163..c9bc0f32af 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -83,6 +83,9 @@ module cmfd_header ! List of CMFD k real(8), allocatable :: k_cmfd(:) + ! Balance keff + real(8) :: keff_bal + end type cmfd_type contains diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 7e6ca18e53..67e8f0d0f3 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -62,16 +62,20 @@ contains use error, only: fatal_error, warning use global use output, only: write_message - use string, only: lower_case + use string, only: to_lower use xml_interface use, intrinsic :: ISO_FORTRAN_ENV + integer :: i integer :: ng + integer :: n_params integer, allocatable :: iarray(:) + integer, allocatable :: int_array(:) logical :: file_exists ! does cmfd.xml exist? logical :: found character(MAX_LINE_LEN) :: filename character(MAX_LINE_LEN) :: temp_str + real(8) :: gs_tol(2) type(Node), pointer :: doc => null() type(Node), pointer :: node_mesh => null() @@ -151,91 +155,121 @@ contains ! Set feedback logical if (check_for_node(doc, "feedback")) then call get_node_value(doc, "feedback", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_feedback = .true. + cmfd_feedback = .true. end if ! Set downscatter logical if (check_for_node(doc, "downscatter")) then call get_node_value(doc, "downscatter", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_downscatter = .true. + cmfd_downscatter = .true. + end if + + ! Reset dhat parameters + if (check_for_node(doc, "dhat_reset")) then + call get_node_value(doc, "dhat_reset", temp_str) + temp_str = to_lower(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & + dhat_reset = .true. end if ! Set the solver type if (check_for_node(doc, "solver")) & - call get_node_value(doc, "solver", cmfd_solver_type) + call get_node_value(doc, "solver", cmfd_solver_type) ! Set monitoring if (check_for_node(doc, "snes_monitor")) then call get_node_value(doc, "snes_monitor", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_snes_monitor = .true. + cmfd_snes_monitor = .true. end if if (check_for_node(doc, "ksp_monitor")) then call get_node_value(doc, "ksp_monitor", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_ksp_monitor = .true. + cmfd_ksp_monitor = .true. end if if (check_for_node(doc, "power_monitor")) then call get_node_value(doc, "power_monitor", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_power_monitor = .true. + cmfd_power_monitor = .true. end if ! Output logicals if (check_for_node(doc, "write_matrices")) then - call get_node_value(doc, "write_matices", temp_str) - call lower_case(temp_str) + call get_node_value(doc, "write_matrices", temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_write_matrices = .true. + cmfd_write_matrices = .true. end if ! Run an adjoint calc if (check_for_node(doc, "run_adjoint")) then call get_node_value(doc, "run_adjoint", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & +#ifndef PETSC + message = 'Must use PETSc when running adjoint option.' + call fatal_error() +#endif cmfd_run_adjoint = .true. end if ! Batch to begin cmfd if (check_for_node(doc, "begin")) & - call get_node_value(doc, "begin", cmfd_begin) + call get_node_value(doc, "begin", cmfd_begin) - ! Tally during inactive batches - if (check_for_node(doc, "inactive")) then - call get_node_value(doc, "inactive", temp_str) - call lower_case(temp_str) - if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') & - cmfd_tally_on = .false. + ! Check for cmfd tally resets + if (check_for_node(doc, "tally_reset")) then + n_cmfd_resets = get_arraysize_integer(doc, "tally_reset") + else + n_cmfd_resets = 0 + end if + if (n_cmfd_resets > 0) then + allocate(int_array(n_cmfd_resets)) + call get_node_array(doc, "tally_reset", int_array) + do i = 1, n_cmfd_resets + call cmfd_reset % add(int_array(i)) + end do + deallocate(int_array) end if - - ! Inactive batch flush window - if (check_for_node(doc, "inactive_flush")) & - call get_node_value(doc, "inactive_flush", cmfd_inact_flush(1)) - if (check_for_node(doc, "num_flushes")) & - call get_node_value(doc, "num_flushes", cmfd_inact_flush(2)) - - ! Last flush before active batches - if (check_for_node(doc, "active_flush")) & - call get_node_value(doc, "active_flush", cmfd_act_flush) ! Get display if (check_for_node(doc, "display")) & - call get_node_value(doc, "display", cmfd_display) + call get_node_value(doc, "display", cmfd_display) if (trim(cmfd_display) == 'dominance' .and. & - trim(cmfd_solver_type) /= 'power') then + trim(cmfd_solver_type) /= 'power') then message = 'Dominance Ratio only aviable with power iteration solver' call warning() cmfd_display = '' end if + ! Read in spectral radius estimate and tolerances + if (check_for_node(doc, "spectral")) & + call get_node_value(doc, "spectral", cmfd_spectral) + if (check_for_node(doc, "shift")) & + call get_node_value(doc, "shift", cmfd_shift) + if (check_for_node(doc, "ktol")) & + call get_node_value(doc, "ktol", cmfd_ktol) + if (check_for_node(doc, "stol")) & + call get_node_value(doc, "stol", cmfd_stol) + if (check_for_node(doc, "gauss_seidel_tolerance")) then + n_params = get_arraysize_double(doc, "gauss_seidel_tolerance") + if (n_params /= 2) then + message = 'Gauss Seidel tolerance is not 2 parameters & + &(absolute, relative).' + call fatal_error() + end if + call get_node_array(doc, "gauss_seidel_tolerance", gs_tol) + cmfd_atoli = gs_tol(1) + cmfd_rtoli = gs_tol(2) + end if + ! Create tally objects call create_cmfd_tally(doc) @@ -405,9 +439,9 @@ contains ! Set reset property if (check_for_node(doc, "reset")) then call get_node_value(doc, "reset", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - t % reset = .true. + t % reset = .true. end if ! Set up mesh filter diff --git a/src/cmfd_power_solver.F90 b/src/cmfd_power_solver.F90 index 65b60bc891..f80c122b71 100644 --- a/src/cmfd_power_solver.F90 +++ b/src/cmfd_power_solver.F90 @@ -242,7 +242,7 @@ contains subroutine convergence(iter) - use constants, only: ONE, TINY_BIT + use constants, only: ONE, ZERO use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV @@ -255,7 +255,7 @@ contains kerr = abs(k_o - k_n)/k_n ! Calculate max error in source - where (s_n % val > TINY_BIT) + where (s_n % val > ZERO) serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2 end where serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val)) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 new file mode 100644 index 0000000000..b3e8fb8a88 --- /dev/null +++ b/src/cmfd_solver.F90 @@ -0,0 +1,819 @@ +module cmfd_solver + +! This module contains routines to execute the power iteration solver + + use cmfd_loss_operator, only: init_loss_matrix, build_loss_matrix + use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix + use matrix_header, only: Matrix + use vector_header, only: Vector + + implicit none + private + public :: cmfd_solver_execute + + real(8) :: k_n ! new k-eigenvalue + real(8) :: k_o ! old k-eigenvalue + real(8) :: k_s ! shift of eigenvalue + real(8) :: k_ln ! new shifted eigenvalue + real(8) :: k_lo ! old shifted eigenvalue + real(8) :: norm_n ! current norm of source vector + real(8) :: norm_o ! old norm of source vector + real(8) :: kerr ! error in keff + real(8) :: serr ! error in source + real(8) :: ktol ! tolerance on keff + real(8) :: stol ! tolerance on source + logical :: adjoint_calc ! run an adjoint calculation + type(Matrix) :: loss ! cmfd loss matrix + type(Matrix) :: prod ! cmfd prod matrix + type(Vector) :: phi_n ! new flux vector + type(Vector) :: phi_o ! old flux vector + type(Vector) :: s_n ! new source vector + type(Vector) :: s_o ! old flux vector + type(Vector) :: serr_v ! error in source + + ! CMFD linear solver interface + procedure(linsolve), pointer :: cmfd_linsolver => null() + abstract interface + subroutine linsolve(A, b, x, tol, i) + import :: Matrix + import :: Vector + type(Matrix), intent(inout) :: A + type(Vector), intent(inout) :: b + type(Vector), intent(inout) :: x + real(8), intent(in) :: tol + integer, intent(out) :: i + end subroutine linsolve + end interface + +contains + +!=============================================================================== +! CMFD_SOLVER_EXECUTE sets up and runs power iteration solver for CMFD +!=============================================================================== + + subroutine cmfd_solver_execute(adjoint) + + use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve + + logical, optional, intent(in) :: adjoint ! adjoint calc + + logical :: physical_adjoint = .false. + + ! Check for adjoint execution + adjoint_calc = .false. + if (present(adjoint)) adjoint_calc = adjoint + + ! Check for physical adjoint + if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') & + physical_adjoint = .true. + + ! Start timer for build + call time_cmfdbuild % start() + + ! Initialize matrices and vectors + call init_data(physical_adjoint) + + ! Check for mathematical adjoint calculation + if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') & + call compute_adjoint() + + ! Stop timer for build + call time_cmfdbuild % stop() + + ! Begin power iteration + call time_cmfdsolve % start() + call execute_power_iter() + call time_cmfdsolve % stop() + + ! Extract results + call extract_results() + + ! Deallocate data + call finalize() + + end subroutine cmfd_solver_execute + +!=============================================================================== +! INIT_DATA allocates matrices and vectors for CMFD solution +!=============================================================================== + + subroutine init_data(adjoint) + + use constants, only: ONE, ZERO + use error, only: fatal_error + use global, only: cmfd, cmfd_shift, keff, cmfd_ktol, cmfd_stol, & + cmfd_write_matrices + + logical, intent(in) :: adjoint + + integer :: n ! problem size + real(8) :: guess ! initial guess + real(8) :: dw ! eigenvalue shift + + ! Set up matrices + call init_loss_matrix(loss) + call init_prod_matrix(prod) + + ! Get problem size + n = loss % n + + ! Set up flux vectors + call phi_n % create(n) + call phi_o % create(n) + + ! Set up source vectors + call s_n % create(n) + call s_o % create(n) + call serr_v % create(n) + + ! Set initial guess + guess = ONE + phi_n % val = guess + phi_o % val = guess + k_n = keff + k_o = k_n + dw = cmfd_shift + k_s = k_o + dw + k_ln = ONE/(ONE/k_n - ONE/k_s) + k_lo = k_ln + + ! Fill in loss matrix + call build_loss_matrix(loss, adjoint=adjoint) + + ! Fill in production matrix + call build_prod_matrix(prod, adjoint=adjoint) + + ! Finalize setup of CSR matrices + call loss % assemble() + call prod % assemble() + if (cmfd_write_matrices) then + call loss % write('loss.dat') + call prod % write('prod.dat') + end if + + ! Set norms to 0 + norm_n = ZERO + norm_o = ZERO + + ! Set up solver + select case(cmfd % indices(4)) + case(1) + cmfd_linsolver => cmfd_linsolver_1g + case(2) + cmfd_linsolver => cmfd_linsolver_2g + case default + cmfd_linsolver => cmfd_linsolver_ng + end select + + ! Set tolerances + ktol = cmfd_ktol + stol = cmfd_stol + + end subroutine init_data + +!=============================================================================== +! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem +!=============================================================================== + + subroutine compute_adjoint() + + use error, only: fatal_error +#ifdef PETSC + use global, only: cmfd_write_matrices +#else + use global, only: message +#endif + +#ifdef PETSC + ! Transpose matrices + call loss % transpose() + call prod % transpose() + + ! Write out matrix in binary file (debugging) + if (cmfd_write_matrices) then + call loss % write_petsc_binary('adj_lossmat.bin') + call prod % write_petsc_binary('adj_prodmat.bin') + end if +#else + message = 'Adjoint calculations only allowed with PETSc' + call fatal_error() +#endif + + end subroutine compute_adjoint + +!=============================================================================== +! EXECUTE_POWER_ITER is the main power iteration routine +! for the cmfd calculation +!=============================================================================== + + subroutine execute_power_iter() + + use constants, only: ONE + use error, only: fatal_error + use global, only: cmfd_atoli, cmfd_rtoli, message + + integer :: i ! iteration counter + integer :: innerits ! # of inner iterations + integer :: totalits ! total number of inners + logical :: iconv ! did the problem converged + real(8) :: atoli ! absolute minimum tolerance + real(8) :: rtoli ! relative tolerance based on source conv + real(8) :: toli ! the current tolerance of inners + + ! Reset convergence flag + iconv = .false. + + ! Set up tolerances + atoli = cmfd_atoli + rtoli = cmfd_rtoli + toli = rtoli*100._8 + + ! Perform shift + call wielandt_shift() + totalits = 0 + + ! Begin power iteration + do i = 1, 10000 + + ! Check if reached iteration 10000 + if (i == 10000) then + message = 'Reached maximum iterations in CMFD power iteration solver.' + call fatal_error() + end if + + ! Compute source vector + call prod % vector_multiply(phi_o, s_o) + + ! Normalize source vector + s_o % val = s_o % val / k_lo + + ! Compute new flux vector + call cmfd_linsolver(loss, s_o, phi_n, toli, innerits) + + ! Compute new source vector + call prod % vector_multiply(phi_n, s_n) + + ! Compute new shifted eigenvalue + k_ln = sum(s_n % val) / sum(s_o % val) + + ! Compute new eigenvalue + k_n = ONE/(ONE/k_ln + ONE/k_s) + + ! Renormalize the old source + s_o % val = s_o % val * k_lo + + ! Check convergence + call convergence(i, innerits, iconv) + totalits = totalits + innerits + + ! Break loop if converged + if (iconv) exit + + ! Record old values + phi_o % val = phi_n % val + k_o = k_n + k_lo = k_ln + norm_o = norm_n + + ! Get new tolerance for inners + toli = max(atoli, rtoli*serr) + + end do + + end subroutine execute_power_iter + +!=============================================================================== +! WIELANDT SHIFT +!=============================================================================== + + subroutine wielandt_shift() + + use constants, only: ONE + + integer :: irow ! row counter + integer :: icol ! col counter + integer :: jcol ! current col index in prod matrix + + ! perform subtraction + jcol = 1 + ROWS: do irow = 1, loss % n + COLS: do icol = loss % get_row(irow), loss % get_row(irow + 1) - 1 + if (loss % get_col(icol) == prod % get_col(jcol) .and. & + jcol < prod % get_row(irow + 1)) then + loss % val(icol) = loss % val(icol) - ONE/k_s*prod % val(jcol) + jcol = jcol + 1 + end if + end do COLS + end do ROWS + + end subroutine wielandt_shift + +!=============================================================================== +! CONVERGENCE checks the convergence of the CMFD problem +!=============================================================================== + + subroutine convergence(iter, innerits, iconv) + + use constants, only: ONE, ZERO + use global, only: cmfd_power_monitor, master + use, intrinsic :: ISO_FORTRAN_ENV + + integer, intent(in) :: iter ! outer iteration number + integer, intent(in) :: innerits ! inner iteration nubmer + logical, intent(out) :: iconv ! convergence logical + + ! Reset convergence flag + iconv = .false. + + ! Calculate error in keff + kerr = abs(k_o - k_n)/k_n + + ! Calculate max error in source + where (s_n % val > ZERO) + serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2 + end where + serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val)) + + ! Check for convergence + if(kerr < ktol .and. serr < stol) iconv = .true. + + ! Save the L2 norm of the source + norm_n = serr + + ! Print out to user + if (cmfd_power_monitor .and. master) then + write(OUTPUT_UNIT,FMT='(I0,":",T10,"k-eff: ",F0.8,T30,"k-error: ", & + &1PE12.5,T55, "src-error: ",1PE12.5,T80,I0)') iter, k_n, kerr, & + serr, innerits + end if + + end subroutine convergence + +!=============================================================================== +! CMFD_LINSOLVER_1g solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver_1g(A, b, x, tol, its) + + use constants, only: ONE, ZERO + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message + + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations + + integer :: g ! group index + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: igs ! Gauss-Seidel iteration counter + integer :: irb ! Red/Black iteration switch + integer :: irow ! row iteration + integer :: icol ! iteration counter over columns + integer :: didx ! index for diagonal component + logical :: found ! did we find col + real(8) :: tmp1 ! temporary sum g1 + real(8) :: x1 ! new g1 value of x + real(8) :: err ! error in convergence of solution + real(8) :: w ! overrelaxation parameter + type(Vector) :: tmpx ! temporary solution vector + + ! Set overrelaxation parameter + w = ONE + + ! Dimensions + ng = 1 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = A % n + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + + ! Copy over x vector + call tmpx % copy(x) + + ! Perform red/black gs iterations + REDBLACK: do irb = 0,1 + + ! Begin loop around matrix rows + ROWS: do irow = 1, n + + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + ! Filter out black cells (even) + if (mod(i+j+k,2) == irb) cycle + + ! Get the index of the diagonals for both rows + call A % search_indices(irow, irow, didx, found) + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + do icol = A % get_row(irow), didx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = didx + 1, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + + ! Solve for new x + x1 = (b % val(irow) - tmp1)/A % val(didx) + + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + + end do ROWS + + end do REDBLACK + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs + if (err < tol) exit + + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver_1g + +!=============================================================================== +! CMFD_LINSOLVER_2G solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver_2g(A, b, x, tol, its) + + use constants, only: ONE, ZERO + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message + + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations + + integer :: g ! group index + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: d1idx ! index of row "1" diagonal + integer :: d2idx ! index of row "2" diagonal + integer :: igs ! Gauss-Seidel iteration counter + integer :: irb ! Red/Black iteration switch + integer :: irow ! row iteration + integer :: icol ! iteration counter over columns + logical :: found ! did we find col + real(8) :: m11 ! block diagonal component 1,1 + real(8) :: m12 ! block diagonal component 1,2 + real(8) :: m21 ! block diagonal component 2,1 + real(8) :: m22 ! block diagonal component 2,2 + real(8) :: dm ! determinant of block diagonal + real(8) :: d11 ! inverse component 1,1 + real(8) :: d12 ! inverse component 1,2 + real(8) :: d21 ! inverse component 2,1 + real(8) :: d22 ! inverse component 2,2 + real(8) :: tmp1 ! temporary sum g1 + real(8) :: tmp2 ! temporary sum g2 + real(8) :: x1 ! new g1 value of x + real(8) :: x2 ! new g2 value of x + real(8) :: err ! error in convergence of solution + real(8) :: w ! overrelaxation parameter + type(Vector) :: tmpx ! temporary solution vector + + ! Set tolerance and overrelaxation parameter + w = ONE + + ! Dimensions + ng = 2 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = A % n + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + + ! Copy over x vector + call tmpx % copy(x) + + ! Perform red/black gs iterations + REDBLACK: do irb = 0,1 + + ! Begin loop around matrix rows + ROWS: do irow = 1, n, 2 + + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + ! Filter out black cells (even) + if (mod(i+j+k,2) == irb) cycle + + ! Get the index of the diagonals for both rows + call A % search_indices(irow, irow, d1idx, found) + call A % search_indices(irow + 1, irow + 1, d2idx, found) + + ! Get block diagonal + m11 = A % val(d1idx) ! group 1 diagonal + m12 = A % val(d1idx + 1) ! group 1 right of diagonal (sorted by col) + m21 = A % val(d2idx - 1) ! group 2 left of diagonal (sorted by col) + m22 = A % val(d2idx) ! group 2 diagonal + + ! Analytically invert the diagonal + dm = m11*m22 - m12*m21 + d11 = m22/dm + d12 = -m12/dm + d21 = -m21/dm + d22 = m11/dm + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + tmp2 = ZERO + do icol = A % get_row(irow), d1idx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = A % get_row(irow + 1), d2idx - 2 + tmp2 = tmp2 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = d1idx + 2, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = d2idx + 1, A % get_row(irow + 2) - 1 + tmp2 = tmp2 + A % val(icol)*x % val(A % get_col(icol)) + end do + + ! Adjust with RHS vector + tmp1 = b % val(irow) - tmp1 + tmp2 = b % val(irow + 1) - tmp2 + + ! Solve for new x + x1 = d11*tmp1 + d12*tmp2 + x2 = d21*tmp1 + d22*tmp2 + + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + x % val(irow + 1) = (ONE - w)*x % val(irow + 1) + w*x2 + + end do ROWS + + end do REDBLACK + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs + if (err < tol) exit + + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver_2g + +!=============================================================================== +! CMFD_LINSOLVER_ng solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver_ng(A, b, x, tol, its) + + use constants, only: ONE, ZERO + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message + + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations + + integer :: g ! group index + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: igs ! Gauss-Seidel iteration counter + integer :: irow ! row iteration + integer :: icol ! iteration counter over columns + integer :: didx ! index for diagonal component + logical :: found ! did we find col + real(8) :: tmp1 ! temporary sum g1 + real(8) :: x1 ! new g1 value of x + real(8) :: err ! error in convergence of solution + real(8) :: w ! overrelaxation parameter + type(Vector) :: tmpx ! temporary solution vector + + ! Set overrelaxation parameter + w = ONE + + ! Dimensions + ng = 1 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = A % n + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + + ! Copy over x vector + call tmpx % copy(x) + + ! Begin loop around matrix rows + ROWS: do irow = 1, n + + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + ! Get the index of the diagonals for both rows + call A % search_indices(irow, irow, didx, found) + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + do icol = A % get_row(irow), didx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + do icol = didx + 1, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + + ! Solve for new x + x1 = (b % val(irow) - tmp1)/A % val(didx) + + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + + end do ROWS + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs + + if (err < tol) exit + + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver_ng + +!=============================================================================== +! EXTRACT_RESULTS takes results and puts them in CMFD global data object +!=============================================================================== + + subroutine extract_results() + + use global, only: cmfd, cmfd_write_matrices, current_batch + + character(len=25) :: filename ! name of file to write data + integer :: n ! problem size + + ! Get problem size + n = loss % n + + ! Allocate in cmfd object if not already allocated + if (adjoint_calc) then + if (.not. allocated(cmfd%adj_phi)) allocate(cmfd%adj_phi(n)) + else + if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n)) + end if + + ! Save values + if (adjoint_calc) then + cmfd % adj_phi = phi_n % val + else + cmfd % phi = phi_n % val + end if + + ! Save eigenvalue + if(adjoint_calc) then + cmfd%adj_keff = k_n + else + cmfd%keff = k_n + end if + + ! Normalize phi to 1 + if (adjoint_calc) then + cmfd%adj_phi = cmfd%adj_phi/sqrt(sum(cmfd%adj_phi*cmfd%adj_phi)) + else + cmfd%phi = cmfd%phi/sqrt(sum(cmfd%phi*cmfd%phi)) + end if + + ! Save dominance ratio + cmfd % dom(current_batch) = norm_n/norm_o + + ! Write out results + if (cmfd_write_matrices) then + if (adjoint_calc) then + filename = 'adj_fluxvec.bin' + else + filename = 'fluxvec.bin' + end if +#ifdef PETSC + call phi_n % write_petsc_binary(filename) +#endif + end if + + end subroutine extract_results + +!=============================================================================== +! MATRIX_TO_INDICES converts a matrix index to spatial and group indicies +!=============================================================================== + + subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + use global, only: cmfd, cmfd_coremap + + integer, intent(out) :: i ! iteration counter for x + integer, intent(out) :: j ! iteration counter for y + integer, intent(out) :: k ! iteration counter for z + integer, intent(out) :: g ! iteration counter for groups + integer, intent(in) :: irow ! iteration counter over row (0 reference) + integer, intent(in) :: nx ! maximum number of x cells + integer, intent(in) :: ny ! maximum number of y cells + integer, intent(in) :: nz ! maximum number of z cells + integer, intent(in) :: ng ! maximum number of groups + + ! Check for core map + if (cmfd_coremap) then + + ! Get indices from indexmap + g = mod(irow-1, ng) + 1 + i = cmfd % indexmap((irow-1)/ng+1,1) + j = cmfd % indexmap((irow-1)/ng+1,2) + k = cmfd % indexmap((irow-1)/ng+1,3) + + else + + ! Compute indices + g = mod(irow-1, ng) + 1 + i = mod(irow-1, ng*nx)/ng + 1 + j = mod(irow-1, ng*nx*ny)/(ng*nx)+ 1 + k = mod(irow-1, ng*nx*ny*nz)/(ng*nx*ny) + 1 + + end if + + end subroutine matrix_to_indices + +!=============================================================================== +! FINALIZE frees all memory associated with power iteration +!=============================================================================== + + subroutine finalize() + + ! Destroy all objects + call loss % destroy() + call prod % destroy() + call phi_n % destroy() + call phi_o % destroy() + call s_n % destroy() + call s_o % destroy() + call serr_v % destroy + + end subroutine finalize + +end module cmfd_solver diff --git a/src/constants.F90 b/src/constants.F90 index 0d52e8e811..ae9c77287a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -335,8 +335,9 @@ module constants ! Source spatial distribution types integer, parameter :: & - SRC_SPACE_BOX = 1, & ! Source in a rectangular prism - SRC_SPACE_POINT = 2 ! Source at a single point + SRC_SPACE_BOX = 1, & ! Source in a rectangular prism + SRC_SPACE_POINT = 2, & ! Source at a single point + SRC_SPACE_FISSION = 3 ! Source in prism filtered by fissionable mats ! Source angular distribution types integer, parameter :: & @@ -392,9 +393,6 @@ module constants ! constant to represent a zero flux "albedo" real(8), parameter :: ZERO_FLUX = 999.0_8 - ! constant to represent albedo rejection - real(8), parameter :: ALBEDO_REJECT = 999.0_8 - ! constant for writing out no residual real(8), parameter :: CMFD_NORES = 99999.0_8 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 18a23a6b0f..fddc79dff0 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -342,21 +342,22 @@ contains integer, intent(in) :: i_nuclide ! index into nuclides array real(8), intent(in) :: E ! energy - integer :: i_energy ! index for energy - integer :: i_low ! band index at lower bounding energy - integer :: i_up ! band index at upper bounding energy - real(8) :: f ! interpolation factor - real(8) :: r ! pseudo-random number - real(8) :: elastic ! elastic cross section - real(8) :: capture ! (n,gamma) cross section - real(8) :: fission ! fission cross section - real(8) :: inelastic ! inelastic cross section - logical :: same_nuc ! do we know the xs for this nuclide at this energy? + integer :: i ! loop index + integer :: i_energy ! index for energy + integer :: i_low ! band index at lower bounding energy + integer :: i_up ! band index at upper bounding energy + integer :: same_nuc_idx ! index of same nuclide + real(8) :: f ! interpolation factor + real(8) :: r ! pseudo-random number + real(8) :: elastic ! elastic cross section + real(8) :: capture ! (n,gamma) cross section + real(8) :: fission ! fission cross section + real(8) :: inelastic ! inelastic cross section + logical :: same_nuc ! do we know the xs for this nuclide at this energy? type(UrrData), pointer, save :: urr => null() type(Nuclide), pointer, save :: nuc => null() type(Reaction), pointer, save :: rxn => null() - type(ListElemInt), pointer :: nuc_list => null() -!$omp threadprivate(urr, nuc, rxn, nuc_list) +!$omp threadprivate(urr, nuc, rxn) micro_xs(i_nuclide) % use_ptable = .true. @@ -381,18 +382,16 @@ contains ! this energy but a different temperature, use the original random number to ! preserve correlation of temperature in probability tables same_nuc = .false. - nuc_list => nuc % nuc_list - do - if (E /= ZERO .and. E == micro_xs(nuc_list % data) % last_E) then + do i = 1, nuc % nuc_list % size() + if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % get_item(i)) % last_E) then same_nuc = .true. + same_nuc_idx = i exit end if - nuc_list => nuc_list % next - if (.not. associated(nuc_list % next)) exit end do if (same_nuc) then - r = micro_xs(nuc_list % data) % last_prn + r = micro_xs(nuc % nuc_list % get_item(same_nuc_idx)) % last_prn else r = prn() micro_xs(i_nuclide) % last_prn = r @@ -513,4 +512,41 @@ contains end subroutine find_energy_index +!=============================================================================== +! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section +! for a given nuclide at the trial relative energy used in resonance scattering +!=============================================================================== + + function elastic_xs_0K(E, nuc) result(xs_out) + + type(Nuclide), pointer :: nuc ! target nuclide at temperature + integer :: i_grid ! index on nuclide energy grid + real(8) :: f ! interp factor on nuclide energy grid + real(8), intent(inout) :: E ! trial energy + real(8) :: xs_out ! 0K xs at trial energy + + ! Determine index on nuclide energy grid + if (E < nuc % energy_0K(1)) then + i_grid = 1 + elseif (E > nuc % energy_0K(nuc % n_grid_0K)) then + i_grid = nuc % n_grid_0K - 1 + else + i_grid = binary_search(nuc % energy_0K, nuc % n_grid_0K, E) + end if + + ! check for rare case where two energy points are the same + if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then + i_grid = i_grid + 1 + end if + + ! calculate interpolation factor + f = (E - nuc % energy_0K(i_grid)) & + & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid)) + + ! Calculate microscopic nuclide elastic cross section + xs_out = (ONE - f) * nuc % elastic_0K(i_grid) & + & + f * nuc % elastic_0K(i_grid + 1) + + end function elastic_xs_0K + end module cross_section diff --git a/src/global.F90 b/src/global.F90 index ff27413eac..0660bf9e52 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,7 +1,7 @@ module global use ace_header, only: Nuclide, SAlphaBeta, xsListing, NuclideMicroXS, & - MaterialMacroXS + MaterialMacroXS, Nuclide0K use bank_header, only: Bank use cmfd_header use constants @@ -299,6 +299,9 @@ module global ! Particle restart run logical :: particle_restart_run = .false. + ! Write out initial source + logical :: write_initial_source = .false. + ! ============================================================================ ! CMFD VARIABLES @@ -329,9 +332,6 @@ module global integer :: n_cmfd_meshes = 1 ! # of structured meshes integer :: n_cmfd_tallies = 3 ! # of user-defined tallies - ! Flag to hold cmfd weight adjustment factors - logical :: cmfd_hold_weights = .false. - ! Eigenvalue solver type character(len=10) :: cmfd_solver_type = 'power' @@ -344,11 +344,9 @@ module global ! Batch to begin cmfd integer :: cmfd_begin = 1 - ! When and how long to flush cmfd tallies during inactive batches - integer :: cmfd_inact_flush(2) = (/9999,1/) - - ! Batch to last flush before active batches - integer :: cmfd_act_flush = 0 + ! Tally reset list + integer :: n_cmfd_resets + type(SetInt) :: cmfd_reset ! Compute effective downscatter cross section logical :: cmfd_downscatter = .false. @@ -366,11 +364,18 @@ module global ! CMFD run logicals logical :: cmfd_on = .false. - logical :: cmfd_tally_on = .true. ! CMFD display info character(len=25) :: cmfd_display = 'balance' + ! Estimate of spectral radius of CMFD matrices and tolerances + real(8) :: cmfd_spectral = ZERO + real(8) :: cmfd_shift = 1.e6 + real(8) :: cmfd_ktol = 1.e-8_8 + real(8) :: cmfd_stol = 1.e-8_8 + real(8) :: cmfd_atoli = 1.e-10_8 + real(8) :: cmfd_rtoli = 1.e-5_8 + ! Information about state points to be written integer :: n_state_points = 0 type(SetInt) :: statepoint_batch @@ -384,6 +389,13 @@ module global logical :: output_xs = .false. logical :: output_tallies = .true. + ! ============================================================================ + ! RESONANCE SCATTERING VARIABLES + + logical :: treat_res_scat = .false. ! is resonance scattering treated? + integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers + type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info + !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, & !$omp& trace, thread_id, current_work, matching_bins) @@ -417,6 +429,11 @@ contains end do deallocate(nuclides) end if + + if (allocated(nuclides_0K)) then + deallocate(nuclides_0K) + end if + if (allocated(sab_tables)) deallocate(sab_tables) if (allocated(xs_listings)) deallocate(xs_listings) if (allocated(micro_xs)) deallocate(micro_xs) diff --git a/src/initialize.F90 b/src/initialize.F90 index 9401224f89..718177fb1e 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -12,6 +12,7 @@ module initialize use global use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict, read_plots_xml + use material_header, only: Material use output, only: title, header, write_summary, print_version, & print_usage, write_xs_summary, print_plot, & write_message diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 7e657271e7..6a9d5e90db 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -11,7 +11,7 @@ module input_xml use output, only: write_message use plot_header use random_lcg, only: prn - use string, only: lower_case, to_str, str_to_int, str_to_real, & + use string, only: to_lower, to_str, str_to_int, str_to_real, & starts_with, ends_with use tally_header, only: TallyObject, TallyFilter use tally_initialize, only: add_tallies @@ -61,16 +61,19 @@ contains character(MAX_FILE_LEN) :: env_variable character(MAX_WORD_LEN) :: type character(MAX_LINE_LEN) :: filename - type(Node), pointer :: doc => null() - type(Node), pointer :: node_mode => null() - type(Node), pointer :: node_source => null() - type(Node), pointer :: node_dist => null() - type(Node), pointer :: node_cutoff => null() - type(Node), pointer :: node_entropy => null() - type(Node), pointer :: node_ufs => null() - type(Node), pointer :: node_sp => null() - type(Node), pointer :: node_output => null() - type(Node), pointer :: node_verb => null() + type(Node), pointer :: doc => null() + type(Node), pointer :: node_mode => null() + type(Node), pointer :: node_source => null() + type(Node), pointer :: node_dist => null() + type(Node), pointer :: node_cutoff => null() + type(Node), pointer :: node_entropy => null() + type(Node), pointer :: node_ufs => null() + type(Node), pointer :: node_sp => null() + type(Node), pointer :: node_output => null() + type(Node), pointer :: node_verb => null() + type(Node), pointer :: node_res_scat => null() + type(Node), pointer :: node_scatterer => null() + type(NodeList), pointer :: node_scat_list => null() ! Display output message message = "Reading settings XML file..." @@ -262,6 +265,14 @@ contains call fatal_error() end if + ! Check if we want to write out source + if (check_for_node(node_source, "write_initial")) then + call get_node_value(node_source, "write_initial", temp_str) + temp_str = to_lower(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & + write_initial_source = .true. + end if + ! Check for external source file if (check_for_node(node_source, "file")) then ! Copy path of source file @@ -287,11 +298,13 @@ contains type = '' if (check_for_node(node_dist, "type")) & call get_node_value(node_dist, "type", type) - call lower_case(type) - select case (trim(type)) + select case (to_lower(type)) case ('box') external_source % type_space = SRC_SPACE_BOX coeffs_reqd = 6 + case ('fission') + external_source % type_space = SRC_SPACE_FISSION + coeffs_reqd = 6 case ('point') external_source % type_space = SRC_SPACE_POINT coeffs_reqd = 3 @@ -337,8 +350,7 @@ contains type = '' if (check_for_node(node_dist, "type")) & call get_node_value(node_dist, "type", type) - call lower_case(type) - select case (trim(type)) + select case (to_lower(type)) case ('isotropic') external_source % type_angle = SRC_ANGLE_ISOTROPIC coeffs_reqd = 0 @@ -389,8 +401,7 @@ contains type = '' if (check_for_node(node_dist, "type")) & call get_node_value(node_dist, "type", type) - call lower_case(type) - select case (trim(type)) + select case (to_lower(type)) case ('monoenergetic') external_source % type_energy = SRC_ENERGY_MONO coeffs_reqd = 1 @@ -440,7 +451,7 @@ contains ! Survival biasing if (check_for_node(doc, "survival_biasing")) then call get_node_value(doc, "survival_biasing", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & survival_biasing = .true. end if @@ -448,7 +459,7 @@ contains ! Probability tables if (check_for_node(doc, "ptables")) then call get_node_value(doc, "ptables", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') & urr_ptables_on = .false. end if @@ -679,19 +690,19 @@ contains ! Check if the user has specified to write binary source file if (check_for_node(node_sp, "separate")) then call get_node_value(node_sp, "separate", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') source_separate = .true. end if if (check_for_node(node_sp, "write")) then call get_node_value(node_sp, "write", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'false' .or. & trim(temp_str) == '0') source_write = .false. end if if (check_for_node(node_sp, "overwrite_latest")) then call get_node_value(node_sp, "overwrite_latest", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') then source_latest = .true. @@ -726,7 +737,7 @@ contains ! batch if (check_for_node(doc, "no_reduce")) then call get_node_value(doc, "no_reduce", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & reduce_tallies = .false. end if @@ -735,7 +746,7 @@ contains ! uncertainties rather than standard deviations if (check_for_node(doc, "confidence_intervals")) then call get_node_value(doc, "confidence_intervals", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') confidence_intervals = .true. end if @@ -749,7 +760,7 @@ contains ! Check for summary option if (check_for_node(node_output, "summary")) then call get_node_value(node_output, "summary", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') output_summary = .true. end if @@ -757,7 +768,7 @@ contains ! Check for cross sections option if (check_for_node(node_output, "cross_sections")) then call get_node_value(node_output, "cross_sections", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. & trim(temp_str) == '1') output_xs = .true. end if @@ -765,7 +776,7 @@ contains ! Check for ASCII tallies output option if (check_for_node(node_output, "tallies")) then call get_node_value(node_output, "tallies", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'false' .or. & trim(temp_str) == '0') output_tallies = .false. end if @@ -774,23 +785,97 @@ contains ! Check for cmfd run if (check_for_node(doc, "run_cmfd")) then call get_node_value(doc, "run_cmfd", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then cmfd_run = .true. -#ifndef PETSC - if (master) then - message = 'CMFD is not available, compile OpenMC with PETSc' - call fatal_error() - end if -#endif + end if + end if + + ! Resonance scattering parameters + if (check_for_node(doc, "resonance_scattering")) then + call get_node_ptr(doc, "resonance_scattering", node_res_scat) + call get_node_list(node_res_scat, "scatterer", node_scat_list) + + ! check that a nuclide is specified + if (get_list_size(node_scat_list) >= 1) then + treat_res_scat = .true. + n_res_scatterers_total = get_list_size(node_scat_list) + + ! store 0K info for resonant scatterers + allocate(nuclides_0K(n_res_scatterers_total)) + do i = 1, n_res_scatterers_total + call get_list_item(node_scat_list, i, node_scatterer) + + ! check to make sure a nuclide is specified + if (.not. check_for_node(node_scatterer, "nuclide")) then + message = "No nuclide specified for scatterer " // trim(to_str(i)) & + // " in settings.xml file!" + call fatal_error() + end if + call get_node_value(node_scatterer, "nuclide", & + nuclides_0K(i) % nuclide) + + if (check_for_node(node_scatterer, "method")) then + call get_node_value(node_scatterer, "method", & + nuclides_0K(i) % scheme) + end if + + ! check to make sure xs name for which method is applied is given + if (.not. check_for_node(node_scatterer, "xs_label")) then + message = "Must specify the temperature dependent name of " // '' & + //"scatterer " // trim(to_str(i)) // " given in cross_sections.xml" + call fatal_error() + end if + call get_node_value(node_scatterer, "xs_label", & + nuclides_0K(i) % name) + + ! check to make sure 0K xs name for which method is applied is given + if (.not. check_for_node(node_scatterer, "xs_label_0K")) then + message = "Must specify the 0K name of " // '' & + //"scatterer "// trim(to_str(i)) // " given in cross_sections.xml" + call fatal_error() + end if + call get_node_value(node_scatterer, "xs_label_0K", & + nuclides_0K(i) % name_0K) + + if (check_for_node(node_scatterer, "E_min")) then + call get_node_value(node_scatterer, "E_min", & + nuclides_0K(i) % E_min) + end if + + ! check that E_min is non-negative + if (nuclides_0K(i) % E_min < ZERO) then + message = "Lower resonance scattering energy bound is negative" + call fatal_error() + end if + + if (check_for_node(node_scatterer, "E_max")) then + call get_node_value(node_scatterer, "E_max", & + nuclides_0K(i) % E_max) + end if + + ! check that E_max is not less than E_min + if (nuclides_0K(i) % E_max < nuclides_0K(i) % E_min) then + message = "Lower resonance scattering energy bound exceeds upper" + call fatal_error() + end if + + nuclides_0K(i) % nuclide = trim(nuclides_0K(i) % nuclide) + nuclides_0K(i) % scheme = to_lower(trim(nuclides_0K(i) % scheme)) + nuclides_0K(i) % name = trim(nuclides_0K(i) % name) + nuclides_0K(i) % name_0K = trim(nuclides_0K(i) % name_0K) + end do + else + message = "No resonant scatterers are specified within the " // "" & + // "resonance_scattering element in settings.xml" + call fatal_error() end if end if ! Natural element expansion option if (check_for_node(doc, "natural_elements")) then call get_node_value(doc, "natural_elements", temp_str) - call lower_case(temp_str) - select case (temp_str) + select case (to_lower(temp_str)) case ('endf/b-vii.0') default_expand = ENDF_BVII0 case ('endf/b-vii.1') @@ -924,8 +1009,7 @@ contains word = '' if (check_for_node(node_cell, "material")) & call get_node_value(node_cell, "material", word) - call lower_case(word) - select case(word) + select case(to_lower(word)) case ('void') c % material = MATERIAL_VOID @@ -1094,8 +1178,7 @@ contains word = '' if (check_for_node(node_surf, "type")) & call get_node_value(node_surf, "type", word) - call lower_case(word) - select case(trim(word)) + select case(to_lower(word)) case ('x-plane') s % type = SURF_PX coeffs_reqd = 1 @@ -1156,8 +1239,7 @@ contains word = '' if (check_for_node(node_surf, "boundary")) & call get_node_value(node_surf, "boundary", word) - call lower_case(word) - select case (trim(word)) + select case (to_lower(word)) case ('transmission', 'transmit', '') s % bc = BC_TRANSMIT case ('vacuum') @@ -1648,8 +1730,7 @@ contains end if ! Adjust material density based on specified units - call lower_case(units) - select case(trim(units)) + select case(to_lower(units)) case ('g/cc', 'g/cm3') mat % density = -val case ('kg/m3') @@ -1804,7 +1885,7 @@ contains ALL_NUCLIDES: do j = 1, mat % n_nuclides ! Check that this nuclide is listed in the cross_sections.xml file name = trim(list_names % get_item(j)) - if (.not. xs_listing_dict % has_key(name)) then + if (.not. xs_listing_dict % has_key(to_lower(name))) then message = "Could not find nuclide " // trim(name) // & " in cross_sections.xml file!" call fatal_error() @@ -1819,20 +1900,20 @@ contains end if ! Find xs_listing and set the name/alias according to the listing - index_list = xs_listing_dict % get_key(name) + index_list = xs_listing_dict % get_key(to_lower(name)) name = xs_listings(index_list) % name alias = xs_listings(index_list) % alias ! If this nuclide hasn't been encountered yet, we need to add its name ! and alias to the nuclide_dict - if (.not. nuclide_dict % has_key(name)) then + if (.not. nuclide_dict % has_key(to_lower(name))) then index_nuclide = index_nuclide + 1 mat % nuclide(j) = index_nuclide - call nuclide_dict % add_key(name, index_nuclide) - call nuclide_dict % add_key(alias, index_nuclide) + call nuclide_dict % add_key(to_lower(name), index_nuclide) + call nuclide_dict % add_key(to_lower(alias), index_nuclide) else - mat % nuclide(j) = nuclide_dict % get_key(name) + mat % nuclide(j) = nuclide_dict % get_key(to_lower(name)) end if ! Copy name and atom/weight percent @@ -1891,7 +1972,7 @@ contains mat % sab_names(j) = name ! Check that this nuclide is listed in the cross_sections.xml file - if (.not. xs_listing_dict % has_key(name)) then + if (.not. xs_listing_dict % has_key(to_lower(name))) then message = "Could not find S(a,b) table " // trim(name) // & " in cross_sections.xml file!" call fatal_error() @@ -1899,17 +1980,17 @@ contains ! Find index in xs_listing and set the name and alias according to the ! listing - index_list = xs_listing_dict % get_key(name) + index_list = xs_listing_dict % get_key(to_lower(name)) name = xs_listings(index_list) % name ! If this S(a,b) table hasn't been encountered yet, we need to add its ! name and alias to the sab_dict - if (.not. sab_dict % has_key(name)) then + if (.not. sab_dict % has_key(to_lower(name))) then index_sab = index_sab + 1 mat % i_sab_tables(j) = index_sab - call sab_dict % add_key(name, index_sab) + call sab_dict % add_key(to_lower(name), index_sab) else - mat % i_sab_tables(j) = sab_dict % get_key(name) + mat % i_sab_tables(j) = sab_dict % get_key(to_lower(name)) end if end do end if @@ -2020,7 +2101,7 @@ contains ! Check for setting if (check_for_node(doc, "assume_separate")) then call get_node_value(doc, "assume_separate", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & assume_separate = .true. end if @@ -2053,8 +2134,7 @@ contains temp_str = '' if (check_for_node(node_mesh, "type")) & call get_node_value(node_mesh, "type", temp_str) - call lower_case(temp_str) - select case (trim(temp_str)) + select case (to_lower(temp_str)) case ('rect', 'rectangle', 'rectangular') m % type = LATTICE_RECT case ('hex', 'hexagon', 'hexagonal') @@ -2234,7 +2314,7 @@ contains temp_str = '' if (check_for_node(node_filt, "type")) & call get_node_value(node_filt, "type", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) ! Determine number of bins if (check_for_node(node_filt, "bins")) then @@ -2462,14 +2542,14 @@ contains end if ! Check to make sure nuclide specified is in problem - if (.not. nuclide_dict % has_key(word)) then + if (.not. nuclide_dict % has_key(to_lower(word))) then message = "The nuclide " // trim(word) // " from tally " // & trim(to_str(t % id)) // " is not present in any material." call fatal_error() end if ! Set bin to index in nuclides array - t % nuclide_bins(j) = nuclide_dict % get_key(word) + t % nuclide_bins(j) = nuclide_dict % get_key(to_lower(word)) end do ! Set number of nuclide bins @@ -2500,7 +2580,7 @@ contains ! (i.e., scatter-p#, flux-y#) n_new = 0 do j = 1, n_words - call lower_case(sarray(j)) + sarray(j) = to_lower(sarray(j)) ! Find if scores(j) is of the form 'moment-p' or 'moment-y' present in ! MOMENT_STRS(:) ! If so, check the order, store if OK, then reset the number to 'n' @@ -2766,6 +2846,12 @@ contains ! Get index of mesh filter k = t % find_filter(FILTER_MESH) + ! Check to make sure mesh filter was specified + if (k == 0) then + message = "Cannot tally surface current without a mesh filter." + call fatal_error() + end if + ! Get pointer to mesh i_mesh = t % filters(k) % int_bins(1) m => meshes(i_mesh) @@ -2942,7 +3028,7 @@ contains temp_str = 'slice' if (check_for_node(node_plot, "type")) & call get_node_value(node_plot, "type", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("slice") pl % type = PLOT_TYPE_SLICE @@ -3009,7 +3095,7 @@ contains temp_str = 'xy' if (check_for_node(node_plot, "basis")) & call get_node_value(node_plot, "basis", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("xy") pl % basis = PLOT_BASIS_XY @@ -3056,7 +3142,7 @@ contains temp_str = "cell" if (check_for_node(node_plot, "color")) & call get_node_value(node_plot, "color", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("cell") @@ -3379,12 +3465,21 @@ contains end if ! create dictionary entry for both name and alias - call xs_listing_dict % add_key(listing % name, i) + call xs_listing_dict % add_key(to_lower(listing % name), i) if (check_for_node(node_ace, "alias")) then - call xs_listing_dict % add_key(listing % alias, i) + call xs_listing_dict % add_key(to_lower(listing % alias), i) end if end do + ! Check that 0K nuclides are listed in the cross_sections.xml file + do i = 1, n_res_scatterers_total + if (.not. xs_listing_dict % has_key(trim(nuclides_0K(i) % name_0K))) then + message = "Could not find nuclide " // trim(nuclides_0K(i) % name_0K) // & + " in cross_sections.xml file!" + call fatal_error() + end if + end do + ! Close cross sections XML file call close_xmldoc(doc) @@ -3410,9 +3505,8 @@ contains character(2) :: element_name element_name = name(1:2) - call lower_case(element_name) - select case (element_name) + select case (to_lower(element_name)) case ('h') call list_names % append('1001.' // xs) call list_density % append(density * 0.999885_8) diff --git a/src/material_header.F90 b/src/material_header.F90 index cbfd917498..fca735fd2a 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -21,6 +21,10 @@ module material_header ! Temporary names read during initialization character(12), allocatable :: names(:) ! isotope names character(12), allocatable :: sab_names(:) ! name of S(a,b) table + + ! Does this material contain fissionable nuclides? + logical :: fissionable = .false. + end type Material end module material_header diff --git a/src/matrix_header.F90 b/src/matrix_header.F90 index 24307b10bb..3e3d4ad0db 100644 --- a/src/matrix_header.F90 +++ b/src/matrix_header.F90 @@ -20,19 +20,22 @@ module matrix_header # endif logical :: petsc_active contains - procedure :: create => matrix_create - procedure :: destroy => matrix_destroy - procedure :: add_value => matrix_add_value - procedure :: new_row => matrix_new_row - procedure :: assemble => matrix_assemble - procedure :: get_row => matrix_get_row - procedure :: get_col => matrix_get_col + procedure :: create => matrix_create + procedure :: destroy => matrix_destroy + procedure :: add_value => matrix_add_value + procedure :: new_row => matrix_new_row + procedure :: assemble => matrix_assemble + procedure :: get_row => matrix_get_row + procedure :: get_col => matrix_get_col procedure :: vector_multiply => matrix_vector_multiply -#ifdef PETSC - procedure :: transpose => matrix_transpose + procedure :: search_indices => matrix_search_indices + procedure :: write => matrix_write + procedure :: copy => matrix_copy +# ifdef PETSC procedure :: setup_petsc => matrix_setup_petsc procedure :: write_petsc_binary => matrix_write_petsc_binary -#endif + procedure :: transpose => matrix_transpose +# endif end type matrix #ifdef PETSC @@ -359,4 +362,87 @@ contains end subroutine matrix_vector_multiply +!=============================================================================== +! MATRIX_SEARCH_INDICES searches for an index in column corresponding to a row +!=============================================================================== + + subroutine matrix_search_indices(self, row, col, idx, found) + + class(Matrix), intent(inout) :: self + integer, intent(in) :: row + integer, intent(in) :: col + integer, intent(out) :: idx + logical, intent(out) :: found + + integer :: j + + found = .false. + + COLS: do j = self % get_row(row), self % get_row(row + 1) - 1 + + if (self % get_col(j) == col) then + idx = j + found = .true. + exit + end if + + end do COLS + + end subroutine matrix_search_indices + +!=============================================================================== +! MATRIX_WRITE writes a matrix to file +!=============================================================================== + + subroutine matrix_write(self, filename) + + character(*), intent(in) :: filename + class(Matrix), intent(inout) :: self + + integer :: unit_ + integer :: i + integer :: j + + open(newunit=unit_, file=filename) + + do i = 1, self % n + do j = self % get_row(i), self % get_row(i + 1) - 1 + write(unit_,*) i, self % get_col(j), self % val(j) + end do + end do + + close(unit_) + + end subroutine matrix_write + +!=============================================================================== +! MATRIX_COPY copies a matrix +!=============================================================================== + + subroutine matrix_copy(self, mattocopy) + + class(Matrix), intent(inout) :: self + type(Matrix), intent(in) :: mattocopy + + ! Set n and nnz + self % n_count = mattocopy % n_count + self % nz_count = mattocopy % nz_count + self % n = mattocopy % n + self % nnz = mattocopy % nnz + + ! Allocate vectors + if (.not.allocated(self % row)) allocate(self % row(self % n + 1)) + if (.not.allocated(self % col)) allocate(self % col(self % nnz)) + if (.not.allocated(self % val)) allocate(self % val(self % nnz)) + + ! Set PETSc active to false + self % petsc_active = .false. + + ! Copy over data + self % row = mattocopy % row + self % col = mattocopy % col + self % val = mattocopy % val + + end subroutine matrix_copy + end module matrix_header diff --git a/src/output.F90 b/src/output.F90 index c0695a7f3e..e6e1c3d487 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -14,7 +14,7 @@ module output use mesh, only: mesh_indices_to_bin, bin_to_mesh_indices use particle_header, only: LocalCoord, Particle use plot_header - use string, only: upper_case, to_str + use string, only: to_upper, to_str use tally_header, only: TallyObject implicit none @@ -131,8 +131,7 @@ contains if (mod(len_trim(msg),2) == 0) m = m + 1 ! convert line to upper case - line = msg - call upper_case(line) + line = to_upper(msg) ! print header based on level select case (header_level) diff --git a/src/physics.F90 b/src/physics.F90 index 47da60f2b9..04cc66afa2 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -2,6 +2,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants + use cross_section, only: elastic_xs_0K use endf, only: reaction_name use error, only: fatal_error, warning use fission, only: nu_total, nu_delayed @@ -349,7 +350,7 @@ contains ! Perform collision physics for elastic scattering call elastic_scatter(i_nuclide, rxn, & - p % E, p % coord0 % uvw, p % mu) + p % E, p % coord0 % uvw, p % mu, p % wgt) end if p % event_MT = ELASTIC @@ -408,13 +409,14 @@ contains ! target. !=============================================================================== - subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab) + subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt) integer, intent(in) :: i_nuclide type(Reaction), pointer :: rxn real(8), intent(inout) :: E real(8), intent(inout) :: uvw(3) real(8), intent(out) :: mu_lab + real(8), intent(inout) :: wgt real(8) :: awr ! atomic weight ratio of target real(8) :: mu_cm ! cosine of polar angle in center-of-mass @@ -437,7 +439,8 @@ contains ! Sample velocity of target nucleus if (.not. micro_xs(i_nuclide) % use_ptable) then - call sample_target_velocity(nuc, v_t, E, uvw) + call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, & + & micro_xs(i_nuclide) % elastic) else v_t = ZERO end if @@ -729,19 +732,270 @@ contains end subroutine sab_scatter !=============================================================================== -! SAMPLE_TARGET_VELOCITY samples the target velocity based on the free gas -! scattering formulation used by most Monte Carlo codes. Excellent documentation -! for this method can be found in FRA-TM-123. +! SAMPLE_TARGET_VELOCITY samples the target velocity. The constant cross section +! free gas model is the default method. Methods for correctly accounting +! for the energy dependence of cross sections in treating resonance elastic +! scattering such as the DBRC, WCM, and a new, accelerated scheme are also +! implemented here. !=============================================================================== - subroutine sample_target_velocity(nuc, v_target, E, uvw) + subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff) - type(Nuclide), pointer :: nuc + type(Nuclide), pointer :: nuc ! target nuclide at temperature T + + real(8), intent(out) :: v_target(3) ! target velocity + real(8), intent(in) :: v_neut(3) ! neutron velocity + real(8), intent(in) :: E ! particle energy + real(8), intent(in) :: uvw(3) ! direction cosines + real(8), intent(inout) :: wgt ! particle weight + + real(8) :: awr ! target/neutron mass ratio + real(8) :: kT ! equilibrium temperature of target in MeV + real(8) :: E_rel ! trial relative energy + real(8) :: xs_0K ! 0K xs at E_rel + real(8) :: xs_eff ! effective elastic xs at temperature T + real(8) :: wcf ! weight correction factor + real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1) + real(8) :: E_low ! lowest practical relative energy + real(8) :: E_up ! highest practical relative energy + real(8) :: E_mode ! most probable Maxwellian energy + real(8) :: E_t_max ! highest practical target energy + real(8) :: E_t ! trial target energy + real(8) :: xs_max ! max 0K xs over practical relative energies + real(8) :: xs_low ! 0K xs at lowest practical relative energy + real(8) :: xs_up ! 0K xs at highest practical relative energy + real(8) :: m ! slope for interpolation + real(8) :: R_dbrc ! DBRC rejection criterion + real(8) :: R_speed ! target speed rejection criterion + real(8) :: cdf_low ! xs cdf at lowest practical relative energy + real(8) :: cdf_up ! xs cdf at highest practical relative energy + real(8) :: cdf_rel ! trial xs cdf value + real(8) :: p_mode ! probability at most probable energy + real(8) :: p_t ! probability at trial target energy + real(8) :: mu ! cosine between neutron and target velocities + + integer :: i_E_low ! 0K index to lowest practical relative energy + integer :: i_E_up ! 0K index to highest practical relative energy + integer :: i_E_rel ! index to trial relative energy + + logical :: reject ! resample if true + + character(80) :: sampling_scheme ! method of target velocity sampling + + kT = nuc % kT + awr = nuc % awr + + ! check if nuclide is a resonant scatterer + if (nuc % resonant) then + + ! sampling scheme to use + sampling_scheme = nuc % scheme + + ! upper resonance scattering energy bound (target is at rest above this E) + if (E > nuc % E_max) then + v_target = ZERO + return + + ! lower resonance scattering energy bound (should be no resonances below) + else if (E < nuc % E_min) then + sampling_scheme = 'cxs' + end if + + ! otherwise, use free gas model + else + if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then + v_target = ZERO + return + else + sampling_scheme = 'cxs' + end if + end if + + ! use appropriate target velocity sampling method + select case (sampling_scheme) + + case ('cxs') + + ! sample target velocity with the constant cross section (cxs) approx. + call sample_cxs_target_velocity(nuc, v_target, E, uvw) + + case ('wcm') + + ! sample target velocity with the constant cross section (cxs) approx. + call sample_cxs_target_velocity(nuc, v_target, E, uvw) + + ! adjust weight as prescribed by the weight correction method (wcm) + E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) + xs_0K = elastic_xs_0K(E_rel, nuc) + wcf = xs_0K / xs_eff + wgt = wcf * wgt + + case ('dbrc') + E_red = sqrt((awr * E) / kT) + E_low = (((E_red - 4.0_8)**2) * kT) / awr + E_up = (((E_red + 4.0_8)**2) * kT) / awr + + ! find lower and upper energy bound indices + ! lower index + if (E_low < nuc % energy_0K(1)) then + i_E_low = 1 + elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_low = nuc % n_grid_0K - 1 + else + i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low) + end if + + ! upper index + if (E_up < nuc % energy_0K(1)) then + i_E_up = 1 + elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_up = nuc % n_grid_0K - 1 + else + i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up) + end if + + ! interpolate xs since we're not exactly at the energy indices + xs_low = nuc % elastic_0K(i_E_low) + m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low)) + xs_up = nuc % elastic_0K(i_E_up) + m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & + & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) + xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) + + ! get max 0K xs value over range of practical relative energies + xs_max = max(xs_low, & + & maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up) + + reject = .true. + + ! sample target velocities until one is accepted by the DBRC + do + + ! sample target velocity with the constant cross section (cxs) approx. + call sample_cxs_target_velocity(nuc, v_target, E, uvw) + + ! perform Doppler broadening rejection correction (dbrc) + E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) + xs_0K = elastic_xs_0K(E_rel, nuc) + R_dbrc = xs_0K / xs_max + if (prn() < R_dbrc) reject = .false. + if (.not. reject) exit + end do + + case ('ares') + E_red = sqrt((awr * E) / kT) + E_low = (((E_red - 4.0_8)**2) * kT) / awr + E_up = (((E_red + 4.0_8)**2) * kT) / awr + + ! find lower and upper energy bound indices + ! lower index + if (E_low < nuc % energy_0K(1)) then + i_E_low = 1 + elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_low = nuc % n_grid_0K - 1 + else + i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low) + end if + + ! upper index + if (E_up < nuc % energy_0K(1)) then + i_E_up = 1 + elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_up = nuc % n_grid_0K - 1 + else + i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up) + end if + + ! interpolate xs CDF since we're not exactly at the energy indices + ! cdf value at lower bound attainable energy + if (i_E_low > 1) then + m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + cdf_low = nuc % xs_cdf(i_E_low - 1) & + & + m * (E_low - nuc % energy_0K(i_E_low)) + else + m = nuc % xs_cdf(i_E_low) & + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + cdf_low = m * (E_low - nuc % energy_0K(i_E_low)) + if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO + end if + + ! cdf value at upper bound attainable energy + m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & + & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) + cdf_up = nuc % xs_cdf(i_E_up - 1) & + & + m * (E_up - nuc % energy_0K(i_E_up)) + + ! values used to sample the Maxwellian + E_mode = kT + p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) & + & * exp(-E_mode / kT) + E_t_max = 16.0_8 * E_mode + + reject = .true. + + do + + ! perform Maxwellian rejection sampling + E_t = E_t_max * prn()**2 + p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) & + & * exp(-E_t / kT) + R_speed = p_t / p_mode + + if (prn() < R_speed) then + + ! sample a relative energy using the xs cdf + cdf_rel = cdf_low + prn() * (cdf_up - cdf_low) + i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), & + & i_E_up - i_E_low + 2, cdf_rel) + E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1) + m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) & + & - nuc % xs_cdf(i_E_low + i_E_rel - 2)) & + & / (nuc % energy_0K(i_E_low + i_E_rel) & + & - nuc % energy_0K(i_E_low + i_E_rel - 1)) + E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m + + ! perform rejection sampling on cosine between + ! neutron and target velocities + mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t)) + + if (abs(mu) < ONE) then + + ! set and accept target velocity + E_t = E_t / awr + v_target = sqrt(E_t) * rotate_angle(uvw, mu) + reject = .false. + end if + end if + + if (.not. reject) exit + end do + + case default + message = "Not a recognized resonance scattering treatment!" + call fatal_error() + end select + + end subroutine sample_target_velocity + +!=============================================================================== +! SAMPLE_CXS_TARGET_VELOCITY samples a target velocity based on the free gas +! scattering formulation, used by most Monte Carlo codes, in which cross section +! is assumed to be constant in energy. Excellent documentation for this method +! can be found in FRA-TM-123. +!=============================================================================== + + subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw) + + type(Nuclide), pointer :: nuc ! target nuclide at temperature real(8), intent(out) :: v_target(3) real(8), intent(in) :: E real(8), intent(in) :: uvw(3) real(8) :: kT ! equilibrium temperature of target in MeV + real(8) :: awr ! target/neutron mass ratio real(8) :: alpha ! probability of sampling f2 over f1 real(8) :: mu ! cosine of angle between neutron and target vel real(8) :: r1, r2 ! pseudo-random numbers @@ -752,18 +1006,10 @@ contains real(8) :: beta_vt_sq ! (beta * speed of target)^2 real(8) :: vt ! speed of target - ! Determine equilibrium temperature in MeV kT = nuc % kT + awr = nuc % awr - ! Check if energy is above threshold - if (E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then - v_target = ZERO - return - end if - - ! calculate beta - beta_vn = sqrt(nuc%awr * E / kT) - + beta_vn = sqrt(awr * E / kT) alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) do @@ -795,20 +1041,20 @@ contains ! Determine rejection probability accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) + /(beta_vn + beta_vt) ! Perform rejection sampling on vt and mu if (prn() < accept_prob) exit end do - ! determine speed of target nucleus - vt = sqrt(beta_vt_sq*kT/nuc % awr) + ! Determine speed of target nucleus + vt = sqrt(beta_vt_sq*kT/awr) - ! determine velocity vector of target nucleus based on neutron's velocity + ! Determine velocity vector of target nucleus based on neutron's velocity ! and the sampled angle between them v_target = vt * rotate_angle(uvw, mu) - end subroutine sample_target_velocity + end subroutine sample_cxs_target_velocity !=============================================================================== ! CREATE_FISSION_SITES determines the average total, prompt, and delayed diff --git a/src/relaxng/cmfd.rnc b/src/relaxng/cmfd.rnc index cbcb7709ae..bbb7fc9e6b 100644 --- a/src/relaxng/cmfd.rnc +++ b/src/relaxng/cmfd.rnc @@ -22,15 +22,9 @@ element cmfd { element feedback { xsd:boolean }? & - element n_cmfd_procs { xsd:int }? & - - element reset { xsd:boolean }? & - - element balance { xsd:boolean }? & - element downscatter { xsd:boolean }? & - element run_2grp { xsd:boolean }? & + element dhat_reset { xsd:boolean }? & element solver { xsd:string }? & @@ -40,8 +34,6 @@ element cmfd { element power_monitor { xsd:boolean }? & - element write_balance { xsd:boolean }? & - element write_matrices { xsd:boolean }? & element run_adjoint { xsd:boolean }? & @@ -50,9 +42,18 @@ element cmfd { element begin { xsd:int }? & - element inactive { xsd:boolean }? & + element tally_reset { list { xsd:int+ } }? & - element active_flush { xsd:int }? & + element display { xsd:string }? & + + element spectral { xsd:double }? & + + element shift { xsd:double }? & + + element ktol { xsd: double }? & + + element stol { xsd: double }? & + + element gauss_seidel_tolerance { list { xsd:double+ } }? - element keff_tol { xsd:double }? } diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index ea3ecabf54..707ecbbc80 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -85,7 +85,8 @@ element settings { attribute interplation { xsd:string { maxLength = "10" } })? & (element parameters { list { xsd:double+ } } | attribute parameters { list { xsd:double+ } })? - }? + }? & + (element write_initial { xsd:boolean } | attribute write_initial { xsd:boolean })? }? & element state_point { @@ -129,5 +130,22 @@ element settings { attribute lower_left { list { xsd:double+ } }) & (element upper_right { list { xsd:double+ } } | attribute upper_right { list { xsd:double+ } }) + }? & + + element resonance_scattering { + element scatterer { + (element nuclide { xsd:string { maxLength = "12" } } | + attribute nuclide { xsd:string { maxLength = "12" } }) & + (element method { xsd:string { maxLength = "16" } } | + attribute method { xsd:string { maxLength = "16" } }) & + (element xs_label { xsd:string { maxLength = "12" } } | + attribute xs_label { xsd:string { maxLength = "12" } }) & + (element xs_label_0K { xsd:string { maxLength = "12" } } | + attribute xs_label_0K { xsd:string { maxLength = "12" } }) & + (element E_min { xsd:double } | + attribute E_min { xsd:double }) & + (element E_max { xsd:double } | + attribute E_max { xsd:double })? + }* }? } diff --git a/src/source.F90 b/src/source.F90 index 9dcfd01445..9575020e6d 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -27,6 +27,7 @@ contains subroutine initialize_source() + character(MAX_FILE_LEN) :: filename integer(8) :: i ! loop index over bank sites integer(8) :: id ! particle id integer(4) :: itmp ! temporary integer @@ -76,6 +77,20 @@ contains end do end if + ! Write out initial source + if (write_initial_source) then + message = 'Writing out initial source guess...' + call write_message(1) +#ifdef HDF5 + filename = trim(path_output) // 'initial_source.h5' +#else + filename = trim(path_output) // 'initial_source.binary' +#endif + call sp % file_create(filename, serial = .false.) + call sp % write_source_bank() + call sp % file_close() + end if + end subroutine initialize_source !=============================================================================== @@ -135,6 +150,42 @@ contains end do call p % clear() + case (SRC_SPACE_FISSION) + ! Repeat sampling source location until a good site has been found + found = .false. + do while (.not.found) + ! Set particle defaults + call p % initialize() + + ! Coordinates sampled uniformly over a box + p_min = external_source % params_space(1:3) + p_max = external_source % params_space(4:6) + r = (/ (prn(), i = 1,3) /) + site % xyz = p_min + r*(p_max - p_min) + + ! Fill p with needed data + p % coord0 % xyz = site % xyz + p % coord0 % uvw = [ ONE, ZERO, ZERO ] + + ! Now search to see if location exists in geometry + call find_cell(p, found) + if (.not. found) then + num_resamples = num_resamples + 1 + if (num_resamples == MAX_EXTSRC_RESAMPLES) then + message = "Maximum number of external source spatial resamples & + &reached!" + call fatal_error() + end if + cycle + end if + if (p % material == MATERIAL_VOID) then + found = .false. + cycle + end if + if (.not. materials(p % material) % fissionable) found = .false. + end do + call p % clear() + case (SRC_SPACE_POINT) ! Point source site % xyz = external_source % params_space diff --git a/src/string.F90 b/src/string.F90 index f78b94d07f..d4f8873a0a 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -152,37 +152,47 @@ contains ! LOWER_CASE converts a string to all lower case characters !=============================================================================== - elemental subroutine lower_case(word) + elemental function to_lower(word) result(word_lower) - character(*), intent(inout) :: word + character(*), intent(in) :: word + character(len=len(word)) :: word_lower integer :: i integer :: ic do i = 1, len(word) ic = ichar(word(i:i)) - if (ic >= 65 .and. ic <= 90) word(i:i) = char(ic+32) + if (ic >= 65 .and. ic <= 90) then + word_lower(i:i) = char(ic+32) + else + word_lower(i:i) = word(i:i) + end if end do - end subroutine lower_case + end function to_lower !=============================================================================== ! UPPER_CASE converts a string to all upper case characters !=============================================================================== - elemental subroutine upper_case(word) + elemental function to_upper(word) result(word_upper) - character(*), intent(inout) :: word + character(*), intent(in) :: word + character(len=len(word)) :: word_upper integer :: i integer :: ic do i = 1, len(word) ic = ichar(word(i:i)) - if (ic >= 97 .and. ic <= 122) word(i:i) = char(ic-32) + if (ic >= 97 .and. ic <= 122) then + word_upper(i:i) = char(ic-32) + else + word_upper(i:i) = word(i:i) + end if end do - end subroutine upper_case + end function to_upper !=============================================================================== ! ZERO_PADDED returns a string of the input integer padded with zeros to the @@ -317,7 +327,7 @@ end function zero_padded ! the loop automatically exits when n_digits = 10. n_digits = n_digits + 1 end do - + end function count_digits !=============================================================================== @@ -349,21 +359,21 @@ end function zero_padded end function int8_to_str !=============================================================================== -! STR_TO_INT converts a string to an integer. +! STR_TO_INT converts a string to an integer. !=============================================================================== function str_to_int(str) result(num) character(*), intent(in) :: str integer(8) :: num - + character(5) :: fmt integer :: w integer :: ioError ! Determine width of string w = len_trim(str) - + ! Create format specifier for reading string write(UNIT=fmt, FMT='("(I",I2,")")') w @@ -404,7 +414,7 @@ end function zero_padded integer :: decimal ! number of places after decimal integer :: width ! total field width - real(8) :: num2 ! absolute value of number + real(8) :: num2 ! absolute value of number character(9) :: fmt ! format specifier for writing number ! set default field width diff --git a/src/vector_header.F90 b/src/vector_header.F90 index e050981c76..05cc7e3cde 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -21,10 +21,11 @@ module vector_header procedure :: create => vector_create procedure :: destroy => vector_destroy procedure :: add_value => vector_add_value -#ifdef PETSC + procedure :: copy => vector_copy +# ifdef PETSC procedure :: setup_petsc => vector_setup_petsc procedure :: write_petsc_binary => vector_write_petsc_binary -#endif +# endif end type Vector #ifdef PETSC @@ -127,4 +128,28 @@ contains end subroutine vector_write_petsc_binary #endif +!=============================================================================== +! VECTOR_COPY allocates a separate vector and copies +!=============================================================================== + + subroutine vector_copy(self, vectocopy) + + class(Vector), target, intent(inout) :: self + type(Vector), intent(in) :: vectocopy + + ! Preallocate vector + if (.not.allocated(self % data)) allocate(self % data(vectocopy % n)) + self % val => self % data(1:vectocopy % n) + + ! Set n + self % n = vectocopy % n + + ! Copy values + self % val = vectocopy % val + + ! Petsc is default not active + self % petsc_active = .false. + + end subroutine vector_copy + end module vector_header diff --git a/tests/test_cmfd_feed/cmfd.xml b/tests/test_cmfd_feed/cmfd.xml index 51def74088..228b260a46 100644 --- a/tests/test_cmfd_feed/cmfd.xml +++ b/tests/test_cmfd_feed/cmfd.xml @@ -12,5 +12,5 @@ dominance power true - + 1.e-15 1.e-20 diff --git a/tests/test_cmfd_nofeed/cmfd.xml b/tests/test_cmfd_nofeed/cmfd.xml index c36a396043..eb6c2c721a 100644 --- a/tests/test_cmfd_nofeed/cmfd.xml +++ b/tests/test_cmfd_nofeed/cmfd.xml @@ -12,5 +12,6 @@ dominance power false + 1.e-15 1.e-20 diff --git a/tests/test_resonance_scattering/geometry.xml b/tests/test_resonance_scattering/geometry.xml new file mode 100644 index 0000000000..612e46132e --- /dev/null +++ b/tests/test_resonance_scattering/geometry.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/tests/test_resonance_scattering/materials.xml b/tests/test_resonance_scattering/materials.xml new file mode 100644 index 0000000000..52a8c04be2 --- /dev/null +++ b/tests/test_resonance_scattering/materials.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/test_resonance_scattering/results.py b/tests/test_resonance_scattering/results.py new file mode 100644 index 0000000000..8ff10971cd --- /dev/null +++ b/tests/test_resonance_scattering/results.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python + +import sys + +# import statepoint +sys.path.append('../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat new file mode 100644 index 0000000000..62b0ac2239 --- /dev/null +++ b/tests/test_resonance_scattering/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +6.452021E-02 1.738736E-03 diff --git a/tests/test_resonance_scattering/settings.xml b/tests/test_resonance_scattering/settings.xml new file mode 100644 index 0000000000..7ce4f23ac7 --- /dev/null +++ b/tests/test_resonance_scattering/settings.xml @@ -0,0 +1,27 @@ + + + + + + U-238 + cxs + 92238.71c + 92238.71c + 5.0e-6 + 40.0e-6 + + + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py new file mode 100644 index 0000000000..6fdbf87459 --- /dev/null +++ b/tests/test_resonance_scattering/test_resonance_scattering.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + try: + test_run() + test_created_statepoint() + test_results() + finally: + teardown()