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/usersguide/input.rst b/docs/source/usersguide/input.rst index cc2d5e7597..1a8d1f2ac7 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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 033ce2a020..3f336767cd 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -287,7 +287,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..6d4f3c740c 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -7,7 +7,7 @@ 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 @@ -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 @@ -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) @@ -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_input.F90 b/src/cmfd_input.F90 index 32625034c3..f006c67f42 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -155,7 +155,7 @@ contains call get_node_value(doc, "feedback", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_feedback = .true. + cmfd_feedback = .true. end if ! Set downscatter logical @@ -163,7 +163,7 @@ contains call get_node_value(doc, "downscatter", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_downscatter = .true. + cmfd_downscatter = .true. end if ! Reset dhat parameters @@ -184,26 +184,26 @@ contains ! 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) 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) 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) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_power_monitor = .true. + cmfd_power_monitor = .true. end if ! Output logicals @@ -211,7 +211,7 @@ contains call get_node_value(doc, "write_matrices", temp_str) call lower_case(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 @@ -228,7 +228,7 @@ contains ! 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) ! Check for cmfd tally resets if (check_for_node(doc, "tally_reset")) then @@ -247,9 +247,9 @@ contains ! 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 = '' @@ -257,17 +257,17 @@ contains ! Read in spectral radius estimate and tolerances if (check_for_node(doc, "spectral")) & - call get_node_value(doc, "spectral", cmfd_spectral) + call get_node_value(doc, "spectral", cmfd_spectral) if (check_for_node(doc, "shift")) & - call get_node_value(doc, "shift", cmfd_shift) + call get_node_value(doc, "shift", cmfd_shift) if (check_for_node(doc, "ktol")) & - call get_node_value(doc, "ktol", cmfd_ktol) + call get_node_value(doc, "ktol", cmfd_ktol) if (check_for_node(doc, "stol")) & - call get_node_value(doc, "stol", cmfd_stol) + call get_node_value(doc, "stol", cmfd_stol) if (check_for_node(doc, "atoli")) & - call get_node_value(doc, "atoli", cmfd_atoli) + call get_node_value(doc, "atoli", cmfd_atoli) if (check_for_node(doc, "rtoli")) & - call get_node_value(doc, "rtoli", cmfd_rtoli) + call get_node_value(doc, "rtoli", cmfd_rtoli) ! Create tally objects call create_cmfd_tally(doc) @@ -440,7 +440,7 @@ contains call get_node_value(doc, "reset", temp_str) call lower_case(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_solver.F90 b/src/cmfd_solver.F90 index 45bc918a14..38e3436a24 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -11,7 +11,6 @@ module cmfd_solver private public :: cmfd_solver_execute - logical :: iconv ! did the problem converged real(8) :: k_n ! new k-eigenvalue real(8) :: k_o ! old k-eigenvalue real(8) :: k_s ! shift of eigenvalue @@ -38,11 +37,11 @@ module cmfd_solver subroutine linsolve(A, b, x, tol, i) import :: Matrix import :: Vector - type(Matrix) :: A - type(Vector) :: b - type(Vector) :: x - real(8) :: tol - integer :: i + 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 @@ -56,7 +55,7 @@ contains use global, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve - logical, optional :: adjoint ! adjoint calc + logical, optional, intent(in) :: adjoint ! adjoint calc logical :: physical_adjoint = .false. @@ -105,7 +104,7 @@ contains use global, only: cmfd, cmfd_shift, keff, cmfd_ktol, cmfd_stol, & cmfd_write_matrices - logical :: adjoint + logical, intent(in) :: adjoint integer :: n ! problem size real(8) :: guess ! initial guess @@ -144,7 +143,7 @@ contains ! Fill in production matrix call build_prod_matrix(prod, adjoint=adjoint) - ! Setup petsc for everything + ! Finalize setup of CSR matrices call loss % assemble() call prod % assemble() if (cmfd_write_matrices) then @@ -210,11 +209,13 @@ contains subroutine execute_power_iter() use constants, only: ONE - use global, only: cmfd_atoli, cmfd_rtoli + 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 @@ -234,6 +235,12 @@ contains ! 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) @@ -256,7 +263,7 @@ contains s_o % val = s_o % val * k_lo ! Check convergence - call convergence(i, innerits) + call convergence(i, innerits, iconv) totalits = totalits + innerits ! Break loop if converged @@ -305,14 +312,15 @@ contains ! CONVERGENCE checks the convergence of the CMFD problem !=============================================================================== - subroutine convergence(iter, innerits) + subroutine convergence(iter, innerits, iconv) use constants, only: ONE, ZERO use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV - integer :: iter ! outer iteration number - integer :: innerits ! inner iteration nubmer + 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. @@ -350,11 +358,11 @@ contains use constants, only: ONE, ZERO use global, only: cmfd, cmfd_spectral - type(Matrix) :: A ! coefficient matrix - type(Vector) :: b ! right hand side vector - type(Vector) :: x ! unknown vector - real(8) :: tol ! tolerance on final error - integer :: its ! number of inner iterations + 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 @@ -430,7 +438,6 @@ contains ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) its = igs -!print *, igs, err, w if (err < tol) exit ! Calculation new overrelaxation parameter @@ -451,11 +458,11 @@ contains use constants, only: ONE, ZERO use global, only: cmfd, cmfd_spectral - type(Matrix) :: A ! coefficient matrix - type(Vector) :: b ! right hand side vector - type(Vector) :: x ! unknown vector - real(8) :: tol ! tolerance on final error - integer :: its ! number of inner iterations + 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 @@ -570,8 +577,6 @@ contains ! Check convergence err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) its = igs -!print *, igs, err, w -!read * if (err < tol) exit ! Calculation new overrelaxation parameter @@ -744,15 +749,15 @@ contains use global, only: cmfd, cmfd_coremap - integer :: i ! iteration counter for x - integer :: j ! iteration counter for y - integer :: k ! iteration counter for z - integer :: g ! iteration counter for groups - integer :: irow ! iteration counter over row (0 reference) - integer :: nx ! maximum number of x cells - integer :: ny ! maximum number of y cells - integer :: nz ! maximum number of z cells - integer :: ng ! maximum number of groups + 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 diff --git a/src/constants.F90 b/src/constants.F90 index f7d419eaac..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 :: & 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 e2509ce1a7..4bb025ff0c 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 @@ -392,6 +392,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) @@ -425,6 +432,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 3e67ae8662..7d4adad493 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -11,6 +11,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 842c104b36..09e6721c91 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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..." @@ -300,6 +303,9 @@ contains 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 @@ -788,6 +794,88 @@ contains 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 = trim(nuclides_0K(i) % scheme) + call lower_case(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) @@ -3189,6 +3277,15 @@ contains 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) 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 3f5aa204f7..3e3d4ad0db 100644 --- a/src/matrix_header.F90 +++ b/src/matrix_header.F90 @@ -368,11 +368,11 @@ contains subroutine matrix_search_indices(self, row, col, idx, found) - class(Matrix) :: self - integer :: row - integer :: col - integer :: idx - logical :: found + class(Matrix), intent(inout) :: self + integer, intent(in) :: row + integer, intent(in) :: col + integer, intent(out) :: idx + logical, intent(out) :: found integer :: j @@ -396,8 +396,8 @@ contains subroutine matrix_write(self, filename) - character(*) :: filename - class(Matrix) :: self + character(*), intent(in) :: filename + class(Matrix), intent(inout) :: self integer :: unit_ integer :: i @@ -421,8 +421,8 @@ contains subroutine matrix_copy(self, mattocopy) - class(Matrix) :: self - type(Matrix) :: mattocopy + class(Matrix), intent(inout) :: self + type(Matrix), intent(in) :: mattocopy ! Set n and nnz self % n_count = mattocopy % n_count 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/settings.rnc b/src/relaxng/settings.rnc index ea3ecabf54..dd4b91dc7d 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -129,5 +129,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 090ea93d1e..f6de6042b1 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -146,6 +146,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/vector_header.F90 b/src/vector_header.F90 index 4d8cb2cf5c..05cc7e3cde 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -28,7 +28,9 @@ module vector_header # endif end type Vector +#ifdef PETSC integer :: petsc_err ! petsc error code +#endif contains @@ -132,8 +134,8 @@ contains subroutine vector_copy(self, vectocopy) - class(Vector), target :: self - type(Vector) :: vectocopy + class(Vector), target, intent(inout) :: self + type(Vector), intent(in) :: vectocopy ! Preallocate vector if (.not.allocated(self % data)) allocate(self % data(vectocopy % n)) 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()