From bb02e301020145443b321373899e4b7a1a91d359 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 6 Oct 2015 16:46:07 +0700 Subject: [PATCH 01/19] Fix various nitpicky things Make procedures pure where possible. Pass argments as intent(in) instead of pointer unless pointer semantics actually needed (never?). Don't initialize local pointers. --- src/ace.F90 | 66 +++++++++++++------------------- src/cross_section.F90 | 12 +++--- src/endf.F90 | 2 +- src/fission.F90 | 19 ++++------ src/interpolation.F90 | 32 +++++++--------- src/math.F90 | 6 +-- src/mesh.F90 | 38 ++++++++----------- src/output.F90 | 53 ++++++++++++-------------- src/physics.F90 | 26 +++++-------- src/search.F90 | 27 +++++++------ src/source.F90 | 3 +- src/string.F90 | 82 +++++++++++++++++++--------------------- src/tally.F90 | 25 ++++++------ src/tally_initialize.F90 | 4 +- src/timer_header.F90 | 12 +----- src/trigger_header.F90 | 10 ++--- src/xml_interface.F90 | 2 +- 17 files changed, 183 insertions(+), 236 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 22fcc6b3d..622d6d7b7 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -45,9 +45,9 @@ contains integer :: temp_table ! temporary value for sorting character(12) :: name ! name of isotope, e.g. 92235.03c character(12) :: alias ! alias of nuclide, e.g. U-235.03c - type(Material), pointer :: mat => null() - type(Nuclide), pointer :: nuc => null() - type(SAlphaBeta), pointer :: sab => null() + type(Material), pointer :: mat + type(Nuclide), pointer :: nuc + type(SAlphaBeta), pointer :: sab type(SetChar) :: already_read ! allocate arrays for ACE table storage and cross section cache @@ -234,7 +234,6 @@ contains !=============================================================================== subroutine read_ace_table(i_table, i_listing) - integer, intent(in) :: i_table ! index in nuclides/sab_tables integer, intent(in) :: i_listing ! index in xs_listings @@ -258,9 +257,9 @@ contains character(10) :: mat ! material identifier character(70) :: comment ! comment for ACE table character(MAX_FILE_LEN) :: filename ! path to ACE cross section library - type(Nuclide), pointer :: nuc => null() - type(SAlphaBeta), pointer :: sab => null() - type(XsListing), pointer :: listing => null() + type(Nuclide), pointer :: nuc + type(SAlphaBeta), pointer :: sab + type(XsListing), pointer :: listing ! determine path, record length, and location of table listing => xs_listings(i_listing) @@ -406,8 +405,6 @@ contains end select deallocate(XSS) - if(associated(nuc)) nullify(nuc) - if(associated(sab)) nullify(sab) end subroutine read_ace_table @@ -417,10 +414,8 @@ contains !=============================================================================== subroutine read_esz(nuc, data_0K) - - type(Nuclide), pointer :: nuc - - logical :: data_0K ! are we reading 0K data? + type(Nuclide), intent(inout) :: nuc + logical, intent(in) :: 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 @@ -507,8 +502,7 @@ contains !=============================================================================== subroutine read_nu_data(nuc) - - type(Nuclide), pointer :: nuc + type(Nuclide), intent(inout) :: nuc integer :: i ! loop index integer :: JXS2 ! location for fission nu data @@ -524,7 +518,7 @@ contains integer :: LOCC ! location of energy distributions for given MT integer :: lc ! locator integer :: length ! length of data to allocate - type(DistEnergy), pointer :: edist => null() + type(DistEnergy), pointer :: edist JXS2 = JXS(2) JXS24 = JXS(24) @@ -707,8 +701,7 @@ contains !=============================================================================== subroutine read_reactions(nuc) - - type(Nuclide), pointer :: nuc + type(Nuclide), intent(inout) :: nuc integer :: i ! loop indices integer :: i_fission ! index in nuc % index_fission @@ -722,7 +715,7 @@ contains integer :: IE ! reaction's starting index on energy grid integer :: NE ! number of energies integer :: NR ! number of interpolation regions - type(Reaction), pointer :: rxn => null() + type(Reaction), pointer :: rxn type(ListInt) :: MTs LMT = JXS(3) @@ -890,8 +883,7 @@ contains !=============================================================================== subroutine read_angular_dist(nuc) - - type(Nuclide), pointer :: nuc + type(Nuclide), intent(inout) :: nuc integer :: JXS8 ! location of angular distribution locators integer :: JXS9 ! location of angular distributions @@ -902,7 +894,7 @@ contains integer :: i ! index in reactions array integer :: j ! index over incoming energies integer :: length ! length of data array to allocate - type(Reaction), pointer :: rxn => null() + type(Reaction), pointer :: rxn JXS8 = JXS(8) JXS9 = JXS(9) @@ -985,13 +977,12 @@ contains !=============================================================================== subroutine read_energy_dist(nuc) - - type(Nuclide), pointer :: nuc + type(Nuclide), intent(inout) :: nuc integer :: LED ! location of energy distribution locators integer :: LOCC ! location of energy distributions for given MT integer :: i ! loop index - type(Reaction), pointer :: rxn => null() + type(Reaction), pointer :: rxn LED = JXS(10) @@ -1019,10 +1010,9 @@ contains !=============================================================================== recursive subroutine get_energy_dist(edist, loc_law, delayed_n) - - type(DistEnergy), pointer :: edist ! energy distribution - integer, intent(in) :: loc_law ! locator for data - logical, optional :: delayed_n ! is this for delayed neutrons? + type(DistEnergy), intent(inout) :: edist ! energy distribution + integer, intent(in) :: loc_law ! locator for data + logical, intent(in), optional :: delayed_n ! is this for delayed neutrons? integer :: LDIS ! location of all energy distributions integer :: LNW ! location of next energy distribution if multiple @@ -1102,7 +1092,6 @@ contains !=============================================================================== function length_energy_dist(lc, law, LOCC, lid) result(length) - integer, intent(in) :: lc ! location in XSS array integer, intent(in) :: law ! energy distribution law integer, intent(in) :: LOCC ! location of energy distribution @@ -1146,7 +1135,7 @@ contains NR = int(XSS(lc + 1)) NE = int(XSS(lc + 2 + 2*NR)) allocate(L(NE)) - L = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) + L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) ! Continue with finding data length length = length + 2 + 2*NR + 2*NE @@ -1204,7 +1193,7 @@ contains NR = int(XSS(lc + 1)) NE = int(XSS(lc + 2 + 2*NR)) allocate(L(NE)) - L = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) + L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) ! Continue with finding data length length = length + 2 + 2*NR + 2*NE @@ -1234,7 +1223,7 @@ contains NR = int(XSS(lc + 1)) NE = int(XSS(lc + 2 + 2*NR)) allocate(L(NE)) - L = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) + L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) ! Continue with finding data length length = length + 2 + 2*NR + 2*NE @@ -1285,7 +1274,7 @@ contains ! in a way inconsistent with the current form of the ACE Format Guide ! (MCNP5 Manual, Vol 3) allocate(L(NE)) - L = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) + L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1)) ! Don't currently do anything with L deallocate(L) ! Continue with finding data length @@ -1301,8 +1290,7 @@ contains !=============================================================================== subroutine read_unr_res(nuc) - - type(Nuclide), pointer :: nuc + type(Nuclide), intent(inout) :: nuc integer :: JXS23 ! location of URR data integer :: lc ! locator @@ -1390,8 +1378,7 @@ contains !=============================================================================== subroutine generate_nu_fission(nuc) - - type(Nuclide), pointer :: nuc + type(Nuclide), intent(inout) :: nuc integer :: i ! index on nuclide energy grid real(8) :: E ! energy @@ -1417,8 +1404,7 @@ contains !=============================================================================== subroutine read_thermal_data(table) - - type(SAlphaBeta), pointer :: table + type(SAlphaBeta), intent(inout) :: table integer :: i ! index for incoming energies integer :: j ! index for outgoing energies diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1c56d961e..287077c46 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -551,13 +551,13 @@ contains ! for a given nuclide at the trial relative energy used in resonance scattering !=============================================================================== - function elastic_xs_0K(E, nuc) result(xs_out) + pure function elastic_xs_0K(E, nuc) result(xs_out) + real(8), intent(inout) :: E ! trial energy + type(Nuclide), intent(in) :: nuc ! target nuclide at temperature + real(8) :: xs_out ! 0K xs at trial energy - 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 + integer :: i_grid ! index on nuclide energy grid + real(8) :: f ! interp factor on nuclide energy grid ! Determine index on nuclide energy grid if (E < nuc % energy_0K(1)) then diff --git a/src/endf.F90 b/src/endf.F90 index fb85262b2..ba324722c 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -11,7 +11,7 @@ contains ! REACTION_NAME gives the name of the reaction for a given MT value !=============================================================================== - function reaction_name(MT) result(string) + pure function reaction_name(MT) result(string) integer, intent(in) :: MT character(20) :: string diff --git a/src/fission.F90 b/src/fission.F90 index 7b2911997..3188138f2 100644 --- a/src/fission.F90 +++ b/src/fission.F90 @@ -15,9 +15,8 @@ contains ! given nuclide and incoming neutron energy !=============================================================================== - function nu_total(nuc, E) result(nu) - - type(Nuclide), pointer :: nuc ! nuclide from which to find nu + pure function nu_total(nuc, E) result(nu) + type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu real(8), intent(in) :: E ! energy of incoming neutron real(8) :: nu ! number of total neutrons emitted per fission @@ -26,7 +25,7 @@ contains real(8) :: c ! polynomial coefficient if (nuc % nu_t_type == NU_NONE) then - call fatal_error("No neutron emission data for table: " // nuc % name) + nu = ERROR_REAL elseif (nuc % nu_t_type == NU_POLYNOMIAL) then ! determine number of coefficients NC = int(nuc % nu_t_data(1)) @@ -49,11 +48,10 @@ contains ! for a given nuclide and incoming neutron energy !=============================================================================== - function nu_prompt(nuc, E) result(nu) - - type(Nuclide), pointer :: nuc ! nuclide from which to find nu - real(8), intent(in) :: E ! energy of incoming neutron - real(8) :: nu ! number of prompt neutrons emitted per fission + pure function nu_prompt(nuc, E) result(nu) + type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu + real(8), intent(in) :: E ! energy of incoming neutron + real(8) :: nu ! number of prompt neutrons emitted per fission integer :: i ! loop index integer :: NC ! number of polynomial coefficients @@ -87,8 +85,7 @@ contains ! for a given nuclide and incoming neutron energy !=============================================================================== - function nu_delayed(nuc, E) result(nu) - + pure function nu_delayed(nuc, E) result(nu) type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu real(8), intent(in) :: E ! energy of incoming neutron real(8) :: nu ! number of delayed neutrons emitted per fission diff --git a/src/interpolation.F90 b/src/interpolation.F90 index 5c44ed7c3..9e28bc086 100644 --- a/src/interpolation.F90 +++ b/src/interpolation.F90 @@ -21,7 +21,7 @@ contains ! tabulated x's and y's. !=============================================================================== - function interpolate_tab1_array(data, x, loc_start) result(y) + pure function interpolate_tab1_array(data, x, loc_start) result(y) real(8), intent(in) :: data(:) ! array of data real(8), intent(in) :: x ! x value to find y at @@ -106,18 +106,16 @@ contains select case (interp) case (LINEAR_LINEAR) r = (x - x0)/(x1 - x0) - y = (1 - r)*y0 + r*y1 + y = y0 + r*(y1 - y0) case (LINEAR_LOG) - r = (log(x) - log(x0))/(log(x1) - log(x0)) - y = (1 - r)*y0 + r*y1 + r = log(x/x0)/log(x1/x0) + y = y0 + r*(y1 - y0) case (LOG_LINEAR) r = (x - x0)/(x1 - x0) - y = exp((1-r)*log(y0) + r*log(y1)) + y = y0*exp(r*log(y1/y0)) case (LOG_LOG) - r = (log(x) - log(x0))/(log(x1) - log(x0)) - y = exp((1-r)*log(y0) + r*log(y1)) - case default - call fatal_error("Unsupported interpolation scheme: " // to_str(interp)) + r = log(x/x0)/log(x1/x0) + y = y0*exp(r*log(y1/y0)) end select end function interpolate_tab1_array @@ -129,7 +127,7 @@ contains ! tabulated x's and y's. !=============================================================================== - function interpolate_tab1_object(obj, x) result(y) + pure function interpolate_tab1_object(obj, x) result(y) type(Tab1), intent(in) :: obj ! ENDF Tab1 interpolable function real(8), intent(in) :: x ! x value to find y at @@ -191,18 +189,16 @@ contains select case (interp) case (LINEAR_LINEAR) r = (x - x0)/(x1 - x0) - y = (1 - r)*y0 + r*y1 + y = y0 + r*(y1 - y0) case (LINEAR_LOG) - r = (log(x) - log(x0))/(log(x1) - log(x0)) - y = (1 - r)*y0 + r*y1 + r = log(x/x0)/log(x1/x0) + y = y0 + r*(y1 - y0) case (LOG_LINEAR) r = (x - x0)/(x1 - x0) - y = exp((1-r)*log(y0) + r*log(y1)) + y = y0*exp(r*log(y1/y0)) case (LOG_LOG) - r = (log(x) - log(x0))/(log(x1) - log(x0)) - y = exp((1-r)*log(y0) + r*log(y1)) - case default - call fatal_error("Unsupported interpolation scheme: " // to_str(interp)) + r = log(x/x0)/log(x1/x0) + y = y0*exp(r*log(y1/y0)) end select end function interpolate_tab1_object diff --git a/src/math.F90 b/src/math.F90 index 6b96baa0d..15aa672e1 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -12,7 +12,7 @@ contains ! distribution with a specified probability level !=============================================================================== - function normal_percentile(p) result(z) + elemental function normal_percentile(p) result(z) real(8), intent(in) :: p ! probability level real(8) :: z ! corresponding z-value @@ -71,7 +71,7 @@ contains ! specified probability level and number of degrees of freedom !=============================================================================== - function t_percentile(p, df) result(t) + elemental function t_percentile(p, df) result(t) real(8), intent(in) :: p ! probability level integer, intent(in) :: df ! degrees of freedom @@ -123,7 +123,7 @@ contains ! the return value will be 1.0. !=============================================================================== - pure function calc_pn(n,x) result(pnx) + elemental function calc_pn(n,x) result(pnx) integer, intent(in) :: n ! Legendre order requested real(8), intent(in) :: x ! Independent variable the Legendre is to be diff --git a/src/mesh.F90 b/src/mesh.F90 index 3d0235d18..905772eb0 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -18,9 +18,8 @@ contains ! GET_MESH_BIN determines the tally bin for a particle in a structured mesh !=============================================================================== - subroutine get_mesh_bin(m, xyz, bin) - - type(RegularMesh), pointer :: m ! mesh pointer + pure subroutine get_mesh_bin(m, xyz, bin) + type(RegularMesh), intent(in) :: m ! mesh pointer real(8), intent(in) :: xyz(:) ! coordinates integer, intent(out) :: bin ! tally bin @@ -71,9 +70,8 @@ contains ! GET_MESH_INDICES determines the indices of a particle in a structured mesh !=============================================================================== - subroutine get_mesh_indices(m, xyz, ijk, in_mesh) - - type(RegularMesh), pointer :: m + pure subroutine get_mesh_indices(m, xyz, ijk, in_mesh) + type(RegularMesh), intent(in) :: m real(8), intent(in) :: xyz(:) ! coordinates to check integer, intent(out) :: ijk(:) ! indices in mesh logical, intent(out) :: in_mesh ! were given coords in mesh? @@ -96,11 +94,10 @@ contains ! use in a TallyObject results array !=============================================================================== - function mesh_indices_to_bin(m, ijk, surface_current) result(bin) - - type(RegularMesh), pointer :: m + pure function mesh_indices_to_bin(m, ijk, surface_current) result(bin) + type(RegularMesh), intent(in) :: m integer, intent(in) :: ijk(:) - logical, optional :: surface_current + logical, intent(in), optional :: surface_current integer :: bin integer :: n_y ! number of mesh cells in y direction @@ -130,9 +127,8 @@ contains ! (i,j) or (i,j,k) indices !=============================================================================== - subroutine bin_to_mesh_indices(m, bin, ijk) - - type(RegularMesh), pointer :: m + pure subroutine bin_to_mesh_indices(m, bin, ijk) + type(RegularMesh), intent(in) :: m integer, intent(in) :: bin integer, intent(out) :: ijk(:) @@ -167,9 +163,9 @@ contains type(Bank), intent(in) :: bank_array(:) ! fission or source bank real(8), intent(out) :: cnt(:,:,:,:) ! weight of sites in each ! cell and energy group - real(8), optional :: energies(:) ! energy grid to search - integer(8), optional :: size_bank ! # of bank sites (on each proc) - logical, optional :: sites_outside ! were there sites outside mesh? + real(8), intent(in), optional :: energies(:) ! energy grid to search + integer(8), intent(in), optional :: size_bank ! # of bank sites (on each proc) + logical, intent(inout), optional :: sites_outside ! were there sites outside mesh? integer :: i ! loop index for local fission sites integer :: n_sites ! size of bank array @@ -262,9 +258,8 @@ contains ! track will score to a mesh tally. !=============================================================================== - function mesh_intersects_2d(m, xyz0, xyz1) result(intersects) - - type(RegularMesh), pointer :: m + pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects) + type(RegularMesh), intent(in) :: m real(8), intent(in) :: xyz0(2) real(8), intent(in) :: xyz1(2) logical :: intersects @@ -328,9 +323,8 @@ contains end function mesh_intersects_2d - function mesh_intersects_3d(m, xyz0, xyz1) result(intersects) - - type(RegularMesh), pointer :: m + pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects) + type(RegularMesh), intent(in) :: m real(8), intent(in) :: xyz0(3) real(8), intent(in) :: xyz1(3) logical :: intersects diff --git a/src/output.F90 b/src/output.F90 index 55cd5b2a5..2aa90985a 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -100,10 +100,9 @@ contains !=============================================================================== subroutine header(msg, unit, level) - character(*), intent(in) :: msg ! header message - integer, optional :: unit ! unit to write to - integer, optional :: level ! specified header level + integer, intent(in), optional :: unit ! unit to write to + integer, intent(in), optional :: level ! specified header level integer :: n ! number of = signs on left integer :: m ! number of = signs on right @@ -195,9 +194,8 @@ contains !=============================================================================== subroutine write_message(message, level) - - character(*) :: message - integer, optional :: level ! verbosity level + character(*), intent(in) :: message + integer, intent(in), optional :: level ! verbosity level integer :: i_start ! starting position integer :: i_end ! ending position @@ -250,7 +248,6 @@ contains !=============================================================================== subroutine print_particle(p) - type(Particle), intent(in) :: p integer :: i ! index for coordinate levels @@ -320,9 +317,8 @@ contains !=============================================================================== subroutine print_nuclide(nuc, unit) - - type(Nuclide), pointer :: nuc - integer, optional :: unit + type(Nuclide), intent(in) :: nuc + integer, intent(in), optional :: unit integer :: i ! loop index over nuclides integer :: unit_ ! unit to write to @@ -334,8 +330,8 @@ contains integer :: size_energy ! memory used for a energy distributions (bytes) integer :: size_urr ! memory used for probability tables (bytes) character(11) :: law ! secondary energy distribution law - type(Reaction), pointer :: rxn => null() - type(UrrData), pointer :: urr => null() + type(Reaction), pointer :: rxn + type(UrrData), pointer :: urr ! set default unit for writing information if (present(unit)) then @@ -438,9 +434,8 @@ contains !=============================================================================== subroutine print_sab_table(sab, unit) - - type(SAlphaBeta), pointer :: sab - integer, optional :: unit + type(SAlphaBeta), intent(in) :: sab + integer, intent(in), optional :: unit integer :: size_sab ! memory used by S(a,b) table integer :: unit_ ! unit to write to @@ -526,8 +521,8 @@ contains integer :: i ! loop index integer :: unit_xs ! cross_sections.out file unit character(MAX_FILE_LEN) :: path ! path of summary file - type(Nuclide), pointer :: nuc => null() - type(SAlphaBeta), pointer :: sab => null() + type(Nuclide), pointer :: nuc + type(SAlphaBeta), pointer :: sab ! Create filename for log file path = trim(path_output) // "cross_sections.out" @@ -681,7 +676,7 @@ contains subroutine print_plot() integer :: i ! loop index for plots - type(ObjectPlot), pointer :: pl => null() + type(ObjectPlot), pointer :: pl ! Display header for plotting call header("PLOTTING SUMMARY") @@ -1183,7 +1178,7 @@ contains !=============================================================================== subroutine write_surface_current(t, unit_tally) - type(TallyObject), pointer :: t + type(TallyObject), intent(in) :: t integer, intent(in) :: unit_tally integer :: i ! mesh index for x @@ -1356,9 +1351,9 @@ contains function get_label(t, i_filter) result(label) - type(TallyObject), pointer :: t ! tally object - integer, intent(in) :: i_filter ! index in filters array - character(100) :: label ! user-specified identifier + type(TallyObject), intent(in) :: t ! tally object + integer, intent(in) :: i_filter ! index in filters array + character(100) :: label ! user-specified identifier integer :: i ! index in cells/surfaces/etc array integer :: bin @@ -1422,12 +1417,12 @@ contains recursive subroutine find_offset(map, goal, univ, final, offset, path) - integer, intent(in) :: map ! Index in maps vector - integer, intent(in) :: goal ! The target cell ID - type(Universe), pointer, intent(in) :: univ ! Universe to begin search - integer, intent(in) :: final ! Target offset - integer, intent(inout) :: offset ! Current offset - character(100) :: path ! Path to offset + integer, intent(in) :: map ! Index in maps vector + integer, intent(in) :: goal ! The target cell ID + type(Universe), intent(in) :: univ ! Universe to begin search + integer, intent(in) :: final ! Target offset + integer, intent(inout) :: offset ! Current offset + character(*), intent(inout) :: path ! Path to offset integer :: i, j ! Index over cells integer :: k, l, m ! Indices in lattice @@ -1439,7 +1434,7 @@ contains integer :: temp_offset ! Looped sum of offsets logical :: this_cell = .false. ! Advance in this cell? logical :: later_cell = .false. ! Fill cells after this one? - type(Cell), pointer:: c ! Pointer to current cell + type(Cell), pointer :: c ! Pointer to current cell type(Universe), pointer :: next_univ ! Next universe to loop through class(Lattice), pointer :: lat ! Pointer to current lattice diff --git a/src/physics.F90 b/src/physics.F90 index 9f8fb535b..c06919758 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -420,9 +420,8 @@ contains !=============================================================================== subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt) - integer, intent(in) :: i_nuclide - type(Reaction), pointer :: rxn + type(Reaction), intent(in) :: rxn real(8), intent(inout) :: E real(8), intent(inout) :: uvw(3) real(8), intent(out) :: mu_lab @@ -759,9 +758,7 @@ contains !=============================================================================== subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff) - - type(Nuclide), pointer :: nuc ! target nuclide at temperature T - + type(Nuclide), intent(in) :: 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 @@ -1006,8 +1003,7 @@ contains !=============================================================================== subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw) - - type(Nuclide), pointer :: nuc ! target nuclide at temperature + type(Nuclide), intent(in) :: nuc ! target nuclide at temperature real(8), intent(out) :: v_target(3) real(8), intent(in) :: E real(8), intent(in) :: uvw(3) @@ -1080,7 +1076,6 @@ contains !=============================================================================== subroutine create_fission_sites(p, i_nuclide, i_reaction) - type(Particle), intent(inout) :: p integer, intent(in) :: i_nuclide integer, intent(in) :: i_reaction @@ -1197,8 +1192,8 @@ contains function sample_fission_energy(nuc, rxn, p) result(E_out) - type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn + type(Nuclide), intent(in) :: nuc + type(Reaction), intent(in) :: rxn type(Particle), intent(inout) :: p ! Particle causing fission real(8) :: E_out ! outgoing energy of fission neutron @@ -1323,8 +1318,8 @@ contains !=============================================================================== subroutine inelastic_scatter(nuc, rxn, p) - type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn + type(Nuclide), intent(in) :: nuc + type(Reaction), intent(in) :: rxn type(Particle), intent(inout) :: p integer :: i ! loop index @@ -1409,8 +1404,7 @@ contains !=============================================================================== function sample_angle(rxn, E) result(mu) - - type(Reaction), pointer :: rxn ! reaction + type(Reaction), intent(in) :: rxn ! reaction real(8), intent(in) :: E ! incoming energy real(8) :: xi ! random number on [0,1) @@ -1536,7 +1530,6 @@ contains !=============================================================================== function rotate_angle(uvw0, mu) result(uvw) - real(8), intent(in) :: uvw0(3) ! directional cosine real(8), intent(in) :: mu ! cosine of angle in lab or CM real(8) :: uvw(3) ! rotated directional cosine @@ -1585,8 +1578,7 @@ contains !=============================================================================== recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q) - - type(DistEnergy), pointer :: edist + type(DistEnergy), intent(in) :: edist real(8), intent(in) :: E_in ! incoming energy of neutron real(8), intent(out) :: E_out ! outgoing energy real(8), intent(inout), optional :: mu_out ! outgoing cosine of angle diff --git a/src/search.F90 b/src/search.F90 index d38dfb986..0c345471c 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -18,7 +18,7 @@ contains ! value lies in the array. This is used extensively for energy grid searching !=============================================================================== - function binary_search_real(array, n, val) result(array_index) + pure function binary_search_real(array, n, val) result(array_index) integer, intent(in) :: n real(8), intent(in) :: array(n) @@ -33,7 +33,8 @@ contains R = n if (val < array(L) .or. val > array(R)) then - call fatal_error("Value outside of array during binary search") + array_index = -1 + return end if n_iteration = 0 @@ -49,8 +50,8 @@ contains ! check for large number of iterations n_iteration = n_iteration + 1 if (n_iteration == MAX_ITERATION) then - call fatal_error("Reached maximum number of iterations on binary & - &search.") + array_index = -2 + return end if end do @@ -58,7 +59,7 @@ contains end function binary_search_real - function binary_search_int4(array, n, val) result(array_index) + pure function binary_search_int4(array, n, val) result(array_index) integer, intent(in) :: n integer, intent(in) :: array(n) @@ -73,7 +74,8 @@ contains R = n if (val < array(L) .or. val > array(R)) then - call fatal_error("Value outside of array during binary search") + array_index = -1 + return end if n_iteration = 0 @@ -89,8 +91,8 @@ contains ! check for large number of iterations n_iteration = n_iteration + 1 if (n_iteration == MAX_ITERATION) then - call fatal_error("Reached maximum number of iterations on binary & - &search.") + array_index = -2 + return end if end do @@ -98,7 +100,7 @@ contains end function binary_search_int4 - function binary_search_int8(array, n, val) result(array_index) + pure function binary_search_int8(array, n, val) result(array_index) integer, intent(in) :: n integer(8), intent(in) :: array(n) @@ -113,7 +115,8 @@ contains R = n if (val < array(L) .or. val > array(R)) then - call fatal_error("Value outside of array during binary search") + array_index = -1 + return end if n_iteration = 0 @@ -129,8 +132,8 @@ contains ! check for large number of iterations n_iteration = n_iteration + 1 if (n_iteration == MAX_ITERATION) then - call fatal_error("Reached maximum number of iterations on binary & - &search.") + array_index = -2 + return end if end do diff --git a/src/source.F90 b/src/source.F90 index 6226517f3..a16eb245d 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -96,8 +96,7 @@ contains !=============================================================================== subroutine sample_external_source(site) - - type(Bank), pointer :: site ! source site + type(Bank), intent(inout) :: site ! source site integer :: i ! dummy loop index real(8) :: r(3) ! sampled coordinates diff --git a/src/string.F90 b/src/string.F90 index 9ca630e29..ce130a212 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -25,7 +25,6 @@ contains !=============================================================================== subroutine split_string(string, words, n) - character(*), intent(in) :: string character(*), intent(out) :: words(MAX_WORDS) integer, intent(out) :: n @@ -166,7 +165,7 @@ contains ! string = concatenated string !=============================================================================== - function concatenate(words, n_words) result(string) + pure function concatenate(words, n_words) result(string) integer, intent(in) :: n_words character(*), intent(in) :: words(n_words) @@ -186,8 +185,7 @@ contains ! TO_LOWER converts a string to all lower case characters !=============================================================================== - function to_lower(word) result(word_lower) - + pure function to_lower(word) result(word_lower) character(*), intent(in) :: word character(len=len(word)) :: word_lower @@ -209,8 +207,7 @@ contains ! TO_UPPER converts a string to all upper case characters !=============================================================================== - function to_upper(word) result(word_upper) - + pure function to_upper(word) result(word_upper) character(*), intent(in) :: word character(len=len(word)) :: word_upper @@ -234,39 +231,38 @@ contains ! integers. !=============================================================================== -function zero_padded(num, n_digits) result(str) - integer, intent(in) :: num - integer, intent(in) :: n_digits - character(11) :: str + function zero_padded(num, n_digits) result(str) + integer, intent(in) :: num + integer, intent(in) :: n_digits + character(11) :: str - character(8) :: zp_form + character(8) :: zp_form - ! Make sure n_digits is reasonable. 10 digits is the maximum needed for the - ! largest integer(4). - if (n_digits > 10) then - call fatal_error('zero_padded called with an unreasonably large & - &n_digits (>10)') - end if + ! Make sure n_digits is reasonable. 10 digits is the maximum needed for the + ! largest integer(4). + if (n_digits > 10) then + call fatal_error('zero_padded called with an unreasonably large & + &n_digits (>10)') + end if - ! Write a format string of the form '(In.m)' where n is the max width and - ! m is the min width. If a sign is present, then n must be one greater - ! than m. - if (num < 0) then - write(zp_form, '("(I", I0, ".", I0, ")")') n_digits+1, n_digits - else - write(zp_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits - end if + ! Write a format string of the form '(In.m)' where n is the max width and + ! m is the min width. If a sign is present, then n must be one greater + ! than m. + if (num < 0) then + write(zp_form, '("(I", I0, ".", I0, ")")') n_digits+1, n_digits + else + write(zp_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits + end if - ! Format the number. - write(str, zp_form) num -end function zero_padded + ! Format the number. + write(str, zp_form) num + end function zero_padded !=============================================================================== ! IS_NUMBER determines whether a string of characters is all 0-9 characters !=============================================================================== - function is_number(word) result(number) - + pure function is_number(word) result(number) character(*), intent(in) :: word logical :: number @@ -286,10 +282,9 @@ end function zero_padded ! sequence of characters !=============================================================================== - logical function starts_with(str, seq) - - character(*) :: str ! string to check - character(*) :: seq ! sequence of characters + pure logical function starts_with(str, seq) + character(*), intent(in) :: str ! string to check + character(*), intent(in) :: seq ! sequence of characters integer :: i integer :: i_start @@ -321,10 +316,9 @@ end function zero_padded ! of characters !=============================================================================== - logical function ends_with(str, seq) - - character(*) :: str ! string to check - character(*) :: seq ! sequence of characters + pure logical function ends_with(str, seq) + character(*), intent(in) :: str ! string to check + character(*), intent(in) :: seq ! sequence of characters integer :: i_start integer :: str_len @@ -350,7 +344,7 @@ end function zero_padded ! integer. !=============================================================================== - function count_digits(num) result(n_digits) + pure function count_digits(num) result(n_digits) integer, intent(in) :: num integer :: n_digits @@ -368,7 +362,7 @@ end function zero_padded ! INT4_TO_STR converts an integer(4) to a string. !=============================================================================== - function int4_to_str(num) result(str) + pure function int4_to_str(num) result(str) integer, intent(in) :: num character(11) :: str @@ -382,7 +376,7 @@ end function zero_padded ! INT8_TO_STR converts an integer(8) to a string. !=============================================================================== - function int8_to_str(num) result(str) + pure function int8_to_str(num) result(str) integer(8), intent(in) :: num character(21) :: str @@ -396,7 +390,7 @@ end function zero_padded ! STR_TO_INT converts a string to an integer. !=============================================================================== - function str_to_int(str) result(num) + pure function str_to_int(str) result(num) character(*), intent(in) :: str integer(8) :: num @@ -421,7 +415,7 @@ end function zero_padded ! STR_TO_REAL converts an arbitrary string to a real(8) !=============================================================================== - function str_to_real(string) result(num) + pure function str_to_real(string) result(num) character(*), intent(in) :: string real(8) :: num @@ -440,7 +434,7 @@ end function zero_padded ! are used. !=============================================================================== - function real_to_str(num, sig_digits) result(string) + pure function real_to_str(num, sig_digits) result(string) real(8), intent(in) :: num ! number to convert integer, optional, intent(in) :: sig_digits ! # of significant digits diff --git a/src/tally.F90 b/src/tally.F90 index da4e82800..56fff7e62 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -37,13 +37,13 @@ contains subroutine score_general(p, t, start_index, filter_index, i_nuclide, & atom_density, flux) - type(Particle), intent(in) :: p - type(TallyObject), pointer, intent(inout) :: t - integer, intent(in) :: start_index - integer, intent(in) :: i_nuclide - integer, intent(in) :: filter_index ! for % results - real(8), intent(in) :: flux ! flux estimate - real(8), intent(in) :: atom_density ! atom/b-cm + type(Particle), intent(in) :: p + type(TallyObject), intent(inout) :: t + integer, intent(in) :: start_index + integer, intent(in) :: i_nuclide + integer, intent(in) :: filter_index ! for % results + real(8), intent(in) :: flux ! flux estimate + real(8), intent(in) :: atom_density ! atom/b-cm integer :: i ! loop index for scoring bins integer :: l ! loop index for nuclides in material @@ -1018,9 +1018,8 @@ contains !=============================================================================== subroutine score_fission_eout(p, t, i_score) - type(Particle), intent(in) :: p - type(TallyObject), pointer :: t + type(TallyObject), intent(inout) :: t integer, intent(in) :: i_score ! index for score integer :: i ! index of outgoing energy filter @@ -1348,9 +1347,9 @@ contains logical :: end_in_mesh ! ending coordinates inside mesh? real(8) :: theta real(8) :: phi - type(TallyObject), pointer :: t + type(TallyObject), pointer :: t type(RegularMesh), pointer :: m - type(Material), pointer :: mat + type(Material), pointer :: mat t => tallies(i_tally) matching_bins(1:t%n_filters) = 1 @@ -1724,7 +1723,7 @@ contains integer :: offset ! offset for distribcell real(8) :: E ! particle energy real(8) :: theta, phi ! Polar and Azimuthal Angles, respectively - type(TallyObject), pointer :: t + type(TallyObject), pointer :: t type(RegularMesh), pointer :: m found_bin = .true. @@ -1948,7 +1947,7 @@ contains logical :: x_same ! same starting/ending x index (i) logical :: y_same ! same starting/ending y index (j) logical :: z_same ! same starting/ending z index (k) - type(TallyObject), pointer :: t + type(TallyObject), pointer :: t type(RegularMesh), pointer :: m TALLY_LOOP: do i = 1, active_current_tallies % size() diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index b7dc5ed98..73aa18c60 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -38,7 +38,7 @@ contains integer :: j ! loop index for filters integer :: n ! temporary stride integer :: max_n_filters = 0 ! maximum number of filters - type(TallyObject), pointer :: t => null() + type(TallyObject), pointer :: t TALLY_LOOP: do i = 1, n_tallies ! Get pointer to tally @@ -88,7 +88,7 @@ contains integer :: k ! loop index for bins integer :: bin ! filter bin entries integer :: type ! type of tally filter - type(TallyObject), pointer :: t => null() + type(TallyObject), pointer :: t ! allocate tally map array -- note that we don't need a tally map for the ! energy_in and energy_out filters diff --git a/src/timer_header.F90 b/src/timer_header.F90 index 0bf1b7aef..6b0580f42 100644 --- a/src/timer_header.F90 +++ b/src/timer_header.F90 @@ -29,13 +29,11 @@ contains !=============================================================================== subroutine timer_start(self) - class(Timer), intent(inout) :: self ! Turn timer on and measure starting time self % running = .true. call system_clock(self % start_counts) - end subroutine timer_start !=============================================================================== @@ -43,7 +41,6 @@ contains !=============================================================================== function timer_get_value(self) result(elapsed) - class(Timer), intent(in) :: self ! the timer real(8) :: elapsed ! total elapsed time @@ -58,7 +55,6 @@ contains else elapsed = self % elapsed end if - end function timer_get_value !=============================================================================== @@ -66,30 +62,26 @@ contains !=============================================================================== subroutine timer_stop(self) - class(Timer), intent(inout) :: self ! Check to make sure timer was running if (.not. self % running) return ! Stop timer and add time - self % elapsed = timer_get_value(self) + self % elapsed = self % get_value() self % running = .false. - end subroutine timer_stop !=============================================================================== ! TIMER_RESET resets a timer to have a zero value !=============================================================================== - subroutine timer_reset(self) - + pure subroutine timer_reset(self) class(Timer), intent(inout) :: self self % running = .false. self % start_counts = 0 self % elapsed = ZERO - end subroutine timer_reset end module timer_header diff --git a/src/trigger_header.F90 b/src/trigger_header.F90 index e137829bd..96421314c 100644 --- a/src/trigger_header.F90 +++ b/src/trigger_header.F90 @@ -1,6 +1,6 @@ module trigger_header - use constants, only: NONE, N_FILTER_TYPES + use constants, only: NONE, N_FILTER_TYPES, ZERO implicit none @@ -13,9 +13,9 @@ module trigger_header real(8) :: threshold ! a convergence threshold character(len=52) :: score_name ! the name of the score integer :: score_index ! the index of the score - real(8) :: variance=0.0 ! temp variance container - real(8) :: std_dev =0.0 ! temp std. dev. container - real(8) :: rel_err =0.0 ! temp rel. err. container + real(8) :: variance = ZERO ! temp variance container + real(8) :: std_dev = ZERO ! temp std. dev. container + real(8) :: rel_err = ZERO ! temp rel. err. container end type TriggerObject !=============================================================================== @@ -23,7 +23,7 @@ module trigger_header !=============================================================================== type KTrigger integer :: trigger_type = 0 - real(8) :: threshold = 0 + real(8) :: threshold = ZERO end type KTrigger end module trigger_header diff --git a/src/xml_interface.F90 b/src/xml_interface.F90 index 4ae05ee28..8f0370b15 100644 --- a/src/xml_interface.F90 +++ b/src/xml_interface.F90 @@ -117,7 +117,7 @@ contains type(Node), pointer, intent(out) :: out_ptr logical :: found_ - type(NodeList), pointer :: elem_list => null() + type(NodeList), pointer :: elem_list ! Set found to false found_ = .false. From 3386191764293ac2176d0afd56d9683b58bd3d41 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 6 Oct 2015 21:18:11 +0700 Subject: [PATCH 02/19] Get rid of union_grid_index module variable. --- src/cross_section.F90 | 54 +++++++++++++++++++------------------------ 1 file changed, 24 insertions(+), 30 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 287077c46..9e084e170 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -13,10 +13,6 @@ module cross_section use search, only: binary_search implicit none - save - - integer :: union_grid_index -!$omp threadprivate(union_grid_index) contains @@ -33,7 +29,8 @@ contains integer :: i_nuclide ! index into nuclides array integer :: i_sab ! index into sab_tables array integer :: j ! index in mat % i_sab_nuclides - integer :: u ! index into logarithmic mapping array + integer :: i_grid ! index into logarithmic mapping array or material + ! union grid real(8) :: atom_density ! atom density of a nuclide logical :: check_sab ! should we check for S(a,b) table? type(Material), pointer :: mat ! current material @@ -52,11 +49,10 @@ contains mat => materials(p % material) ! Find energy index on energy grid - u = 0 if (grid_method == GRID_MAT_UNION) then - call find_energy_index(p % E, p % material) + i_grid = find_energy_index(mat, p % E) else if (grid_method == GRID_LOGARITHM) then - u = int(log(p % E/energy_min_neutron)/log_spacing) + i_grid = int(log(p % E/energy_min_neutron)/log_spacing) end if ! Determine if this material has S(a,b) tables @@ -99,9 +95,9 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid) end if ! ======================================================================== @@ -142,18 +138,19 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, u) - + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, i_log_union) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array + real(8), intent(in) :: E ! energy integer, intent(in) :: i_mat ! index into materials array integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material - integer, intent(in) :: u ! index into logarithmic mapping array + integer, intent(in) :: i_log_union ! index into logarithmic mapping array or + ! material union energy grid + integer :: i_grid ! index on nuclide energy grid integer :: i_low ! lower logarithmic mapping index integer :: i_high ! upper logarithmic mapping index - real(8), intent(in) :: E ! energy - real(8) :: f ! interp factor on nuclide energy grid + real(8) :: f ! interp factor on nuclide energy grid type(Nuclide), pointer :: nuc type(Material), pointer :: mat @@ -165,7 +162,7 @@ contains select case (grid_method) case (GRID_MAT_UNION) - i_grid = mat % nuclide_grid_index(i_nuc_mat, union_grid_index) + i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union) case (GRID_LOGARITHM) ! Determine the energy grid index using a logarithmic mapping to reduce @@ -178,8 +175,8 @@ contains else ! Determine bounding indices based on which equal log-spaced interval ! the energy is in - i_low = nuc % grid_index(u) - i_high = nuc % grid_index(u + 1) + 1 + i_low = nuc % grid_index(i_log_union) + i_high = nuc % grid_index(i_log_union + 1) + 1 ! Perform binary search over reduced range i_grid = binary_search(nuc % energy(i_low:i_high), & @@ -526,25 +523,22 @@ contains ! energy !=============================================================================== - subroutine find_energy_index(E, i_mat) - - real(8), intent(in) :: E ! energy of particle - integer, intent(in) :: i_mat ! material index - type(Material), pointer :: mat ! pointer to current material - - mat => materials(i_mat) + pure function find_energy_index(mat, E) result(i) + type(Material), intent(in) :: mat ! pointer to current material + real(8), intent(in) :: E ! energy of particle + integer :: i ! energy grid index ! if the energy is outside of energy grid range, set to first or last ! index. Otherwise, do a binary search through the union energy grid. if (E <= mat % e_grid(1)) then - union_grid_index = 1 + i = 1 elseif (E > mat % e_grid(mat % n_grid)) then - union_grid_index = mat % n_grid - 1 + i = mat % n_grid - 1 else - union_grid_index = binary_search(mat % e_grid, mat % n_grid, E) + i = binary_search(mat % e_grid, mat % n_grid, E) end if - end subroutine find_energy_index + end function find_energy_index !=============================================================================== ! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section @@ -552,7 +546,7 @@ contains !=============================================================================== pure function elastic_xs_0K(E, nuc) result(xs_out) - real(8), intent(inout) :: E ! trial energy + real(8), intent(in) :: E ! trial energy type(Nuclide), intent(in) :: nuc ! target nuclide at temperature real(8) :: xs_out ! 0K xs at trial energy From e6eb35bae8e95a2a29208749a37f500fa4453c89 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Oct 2015 16:25:29 +0700 Subject: [PATCH 03/19] Create dictionary on Nuclide type mapping MT -> index in reactions This dictionary is used to speed up tallying when a user requests a tally of a specific reaction, e.g. (n,gamma). --- src/ace.F90 | 1 + src/ace_header.F90 | 10 +++++--- src/tally.F90 | 62 +++++++++++++++++++--------------------------- 3 files changed, 34 insertions(+), 39 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 622d6d7b7..c0c4395cf 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -816,6 +816,7 @@ contains ! Create set of MT values do i = 1, size(nuc % reactions) call MTs % append(nuc % reactions(i) % MT) + call nuc%reaction_index%add_key(nuc%reactions(i)%MT, i) end do ! Create total, absorption, and fission cross sections diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 467887c19..1221ee81d 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -1,8 +1,9 @@ module ace_header - use constants, only: MAX_FILE_LEN, ZERO - use endf_header, only: Tab1 - use list_header, only: ListInt + use constants, only: MAX_FILE_LEN, ZERO + use dict_header, only: DictIntInt + use endf_header, only: Tab1 + use list_header, only: ListInt implicit none @@ -154,6 +155,8 @@ module ace_header ! Reactions integer :: n_reaction ! # of reactions type(Reaction), pointer :: reactions(:) => null() + type(DictIntInt) :: reaction_index ! map MT values to index in reactions + ! array; used at tally-time ! Type-Bound procedures contains @@ -417,6 +420,7 @@ module ace_header end if call this % nuc_list % clear() + call this % reaction_index % clear() end subroutine nuclide_clear diff --git a/src/tally.F90 b/src/tally.F90 index 56fff7e62..2f178fe46 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -691,26 +691,20 @@ contains score = ZERO if (i_nuclide > 0) then - ! TODO: The following search for the matching reaction could - ! be replaced by adding a dictionary on each Nuclide instance - ! of the form {MT: i_reaction, ...} - REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction - ! Get pointer to reaction + if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) rxn => nuclides(i_nuclide) % reactions(m) - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - if (i_energy >= rxn % threshold) then - score = ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - exit REACTION_LOOP + + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor + if (i_energy >= rxn % threshold) then + score = ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density * flux end if - end do REACTION_LOOP + end if else ! Get pointer to current material @@ -718,28 +712,24 @@ contains do l = 1, mat % n_nuclides ! Get atom density atom_density_ = mat % atom_density(l) + ! Get index in nuclides array i_nuc = mat % nuclide(l) - ! TODO: The following search for the matching reaction could - ! be replaced by adding a dictionary on each Nuclide - ! instance of the form {MT: i_reaction, ...} - do m = 1, nuclides(i_nuc) % n_reaction - ! Get pointer to reaction + + if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then + m = nuclides(i_nuc)%reaction_index%get_key(score_bin) rxn => nuclides(i_nuc) % reactions(m) - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuc) % index_grid - f = micro_xs(i_nuc) % interp_factor - if (i_energy >= rxn % threshold) then - score = score + ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density_ * flux - end if - exit + + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuc) % index_grid + f = micro_xs(i_nuc) % interp_factor + if (i_energy >= rxn % threshold) then + score = score + ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density_ * flux end if - end do + end if end do end if From 8d36d577e0391d68a361b71fae951d0b239c60d3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Oct 2015 17:07:41 +0700 Subject: [PATCH 04/19] Make Nuclide%nuc_list a type(VectorInt) instead of type(ListInt) --- src/ace.F90 | 2 +- src/ace_header.F90 | 5 ++--- src/cross_section.F90 | 4 ++-- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index c0c4395cf..dbed66757 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -1587,7 +1587,7 @@ contains do i = 1, n_nuclides_total do j = 1, n_nuclides_total if (nuclides(i) % zaid == nuclides(j) % zaid) then - call nuclides(i) % nuc_list % append(j) + call nuclides(i) % nuc_list % push_back(j) end if end do end do diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 1221ee81d..c5e9e1479 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -3,7 +3,7 @@ module ace_header use constants, only: MAX_FILE_LEN, ZERO use dict_header, only: DictIntInt use endf_header, only: Tab1 - use list_header, only: ListInt + use stl_vector, only: VectorInt implicit none @@ -100,7 +100,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(ListInt) :: nuc_list + type(VectorInt) :: nuc_list ! Energy grid information integer :: n_grid ! # of nuclide grid points @@ -419,7 +419,6 @@ module ace_header deallocate(this % reactions) end if - call this % nuc_list % clear() call this % reaction_index % clear() end subroutine nuclide_clear diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 9e084e170..b9c76b503 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -405,7 +405,7 @@ contains ! preserve correlation of temperature in probability tables same_nuc = .false. do i = 1, nuc % nuc_list % size() - if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % get_item(i)) % last_E) then + if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % data(i)) % last_E) then same_nuc = .true. same_nuc_idx = i exit @@ -413,7 +413,7 @@ contains end do if (same_nuc) then - r = micro_xs(nuc % nuc_list % get_item(same_nuc_idx)) % last_prn + r = micro_xs(nuc % nuc_list % data(same_nuc_idx)) % last_prn else r = prn() micro_xs(i_nuclide) % last_prn = r From f3b1de3b0d4c73beb0259a80451adf67903dc235 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Oct 2015 20:51:41 +0700 Subject: [PATCH 05/19] Get rid of some unnecessary deallocates --- src/ace_header.F90 | 25 ------------------------- 1 file changed, 25 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index c5e9e1479..950a998a3 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -375,31 +375,6 @@ module ace_header integer :: i ! Loop counter - if (allocated(this % energy)) & - deallocate(this % energy, this % total, this % elastic, & - & 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) - - if (allocated(this % nu_t_data)) deallocate(this % nu_t_data) - if (allocated(this % nu_p_data)) deallocate(this % nu_p_data) - if (allocated(this % nu_d_data)) deallocate(this % nu_d_data) - - if (allocated(this % nu_d_precursor_data)) & - deallocate(this % nu_d_precursor_data) - if (associated(this % nu_d_edist)) then do i = 1, size(this % nu_d_edist) call this % nu_d_edist(i) % clear() From ceea23d01aa76645ee2d0bce4c4aebb66ba06fc9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Oct 2015 21:24:14 +0700 Subject: [PATCH 06/19] Remove a bunch of clear methods. Make multiplicity_E allocatable. --- src/ace_header.F90 | 44 +------------------------------------------- 1 file changed, 1 insertion(+), 43 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 950a998a3..a10c6fe87 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -18,10 +18,6 @@ module ace_header integer, allocatable :: type(:) ! type of distribution integer, allocatable :: location(:) ! location of each table real(8), allocatable :: data(:) ! angular distribution data - - ! Type-Bound procedures - contains - procedure :: clear => distangle_clear ! Deallocates DistAngle end type DistAngle !=============================================================================== @@ -52,7 +48,7 @@ module ace_header integer :: MT ! ENDF MT value real(8) :: Q_value ! Reaction Q value integer :: multiplicity ! Number of secondary particles released - type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield + type(Tab1), allocatable :: multiplicity_E ! Energy-dependent neutron yield integer :: threshold ! Energy grid index of threshold logical :: scatter_in_cm ! scattering system in center-of-mass? logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity @@ -80,10 +76,6 @@ module ace_header logical :: multiply_smooth ! multiply by smooth cross section? real(8), allocatable :: energy(:) ! incident energies real(8), allocatable :: prob(:,:,:) ! actual probabibility tables - - ! Type-Bound procedures - contains - procedure :: clear => urrdata_clear ! Deallocates UrrData end type UrrData !=============================================================================== @@ -169,14 +161,12 @@ module ace_header !=============================================================================== 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_8 ! lower cutoff energy for res scattering real(8) :: E_max = 1000.0e-6_8 ! upper cutoff energy for res scattering - end type Nuclide0K !=============================================================================== @@ -296,19 +286,6 @@ module ace_header contains -!=============================================================================== -! DISTANGLE_CLEAR resets and deallocates data in Reaction. -!=============================================================================== - - subroutine distangle_clear(this) - - class(DistAngle), intent(inout) :: this ! The DistAngle object to clear - - if (allocated(this % energy)) & - deallocate(this % energy, this % type, this % location, this % data) - - end subroutine distangle_clear - !=============================================================================== ! DISTENERGY_CLEAR resets and deallocates data in DistEnergy. !=============================================================================== @@ -339,32 +316,13 @@ module ace_header class(Reaction), intent(inout) :: this ! The Reaction object to clear - if (allocated(this % sigma)) deallocate(this % sigma) - - if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E) - if (associated(this % edist)) then call this % edist % clear() deallocate(this % edist) end if - call this % adist % clear() - end subroutine reaction_clear -!=============================================================================== -! URRDATA_CLEAR resets and deallocates data in Reaction. -!=============================================================================== - - subroutine urrdata_clear(this) - - class(UrrData), intent(inout) :: this ! The UrrData object to clear - - if (allocated(this % energy)) & - deallocate(this % energy, this % prob) - - end subroutine urrdata_clear - !=============================================================================== ! NUCLIDE_CLEAR resets and deallocates data in Nuclide. !=============================================================================== From 5751ce4a6688b646d56fc5018457ba9d2cc111a5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Oct 2015 21:30:02 +0700 Subject: [PATCH 07/19] Get rid of tally clearing routines --- src/global.F90 | 9 +---- src/tally_header.F90 | 79 -------------------------------------------- 2 files changed, 1 insertion(+), 87 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index a4c80daa7..1d8099630 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -463,14 +463,7 @@ contains ! Deallocate tally-related arrays if (allocated(global_tallies)) deallocate(global_tallies) if (allocated(meshes)) deallocate(meshes) - if (allocated(tallies)) then - ! First call the clear routines - do i = 1, size(tallies) - call tallies(i) % clear() - end do - ! Now deallocate the tally array - deallocate(tallies) - end if + if (allocated(tallies)) deallocate(tallies) if (allocated(matching_bins)) deallocate(matching_bins) if (allocated(tally_maps)) deallocate(tally_maps) diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 01dcd9bdb..18b921952 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -58,10 +58,6 @@ module tally_header integer :: offset = 0 ! Only used for distribcell filters integer, allocatable :: int_bins(:) real(8), allocatable :: real_bins(:) ! Only used for energy filters - - ! Type-Bound procedures - contains - procedure :: clear => tallyfilter_clear ! Deallocates TallyFilter end type TallyFilter !=============================================================================== @@ -129,81 +125,6 @@ module tally_header ! Tally precision triggers integer :: n_triggers = 0 ! # of triggers type(TriggerObject), allocatable :: triggers(:) ! Array of triggers - - ! Type-Bound procedures - contains - procedure :: clear => tallyobject_clear ! Deallocates TallyObject end type TallyObject - contains - -!=============================================================================== -! TALLYFILTER_CLEAR deallocates a TallyFilter element and sets it to its as -! initialized state. -!=============================================================================== - - subroutine tallyfilter_clear(this) - class(TallyFilter), intent(inout) :: this ! The TallyFilter to be cleared - - this % type = NONE - this % n_bins = 0 - if (allocated(this % int_bins)) & - deallocate(this % int_bins) - if (allocated(this % real_bins)) & - deallocate(this % real_bins) - - end subroutine tallyfilter_clear - -!=============================================================================== -! TALLYOBJECT_CLEAR deallocates a TallyObject element and sets it to its as -! initialized state. -!=============================================================================== - - subroutine tallyobject_clear(this) - class(TallyObject), intent(inout) :: this ! The TallyObject to be cleared - - integer :: i ! Loop Index - - ! This routine will go through each item in TallyObject and set the value - ! to its default, as-initialized values, including deallocations. - this % name = "" - - if (allocated(this % filters)) then - do i = 1, size(this % filters) - call this % filters(i) % clear() - end do - deallocate(this % filters) - end if - - if (allocated(this % stride)) & - deallocate(this % stride) - - this % find_filter = 0 - - this % n_nuclide_bins = 0 - if (allocated(this % nuclide_bins)) & - deallocate(this % nuclide_bins) - this % all_nuclides = .false. - - this % n_score_bins = 0 - if (allocated(this % score_bins)) & - deallocate(this % score_bins) - if (allocated(this % moment_order)) & - deallocate(this % moment_order) - this % n_user_score_bins = 0 - - if (allocated(this % results)) & - deallocate(this % results) - - this % reset = .false. - - this % n_realizations = 0 - - if (allocated(this % triggers)) & - deallocate (this % triggers) - - this % n_triggers = 0 - - end subroutine tallyobject_clear - end module tally_header From 9b87c1b77dfb1d5fc55d8c0de6d14941c0ae6c30 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Oct 2015 21:45:36 +0700 Subject: [PATCH 08/19] Rewrite neighbor_lists using type(VectorInt) --- src/geometry.F90 | 88 +++++++++++++++--------------------------------- 1 file changed, 27 insertions(+), 61 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index e922479b1..75165158b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -9,6 +9,7 @@ module geometry use particle_header, only: LocalCoord, Particle use particle_restart_write, only: write_particle_restart use surface_header + use stl_vector, only: VectorInt use string, only: to_str use tally, only: score_surface_current @@ -871,81 +872,46 @@ contains subroutine neighbor_lists() - integer :: i ! index in cells/surfaces array - integer :: j ! index of surface in cell - integer :: i_surface ! index in count arrays - integer, allocatable :: count_positive(:) ! # of cells on positive side - integer, allocatable :: count_negative(:) ! # of cells on negative side - logical :: positive ! positive side specified in surface list - type(Cell), pointer :: c + integer :: i ! index in cells/surfaces array + integer :: j ! index in region specification + integer :: k ! surface half-space spec + type(VectorInt), allocatable :: neighbor_pos(:) + type(VectorInt), allocatable :: neighbor_neg(:) call write_message("Building neighboring cells lists for each surface...", & - &4) + 4) - allocate(count_positive(n_surfaces)) - allocate(count_negative(n_surfaces)) - count_positive = 0 - count_negative = 0 + allocate(neighbor_pos(n_surfaces)) + allocate(neighbor_neg(n_surfaces)) do i = 1, n_cells - c => cells(i) + do j = 1, size(cells(i)%region) + ! Get token from region specification and skip any tokens that + ! correspond to operators rather than regions + k = cells(i)%region(j) + if (abs(k) >= OP_UNION) cycle - ! loop over each region specification - do j = 1, size(c%region) - i_surface = c % region(j) - positive = (i_surface > 0) - - ! Skip any tokens that correspond to operators rather than regions - i_surface = abs(i_surface) - if (i_surface >= OP_UNION) cycle - - if (positive) then - count_positive(i_surface) = count_positive(i_surface) + 1 + ! Add this cell ID to neighbor list for k-th surface + if (k > 0) then + call neighbor_pos(abs(k))%push_back(i) else - count_negative(i_surface) = count_negative(i_surface) + 1 + call neighbor_neg(abs(k))%push_back(i) end if end do end do - ! allocate neighbor lists for each surface do i = 1, n_surfaces - if (count_positive(i) > 0) then - allocate(surfaces(i)%obj%neighbor_pos(count_positive(i))) - end if - if (count_negative(i) > 0) then - allocate(surfaces(i)%obj%neighbor_neg(count_negative(i))) - end if + ! Copy positive neighbors to Surface instance + j = neighbor_pos(i)%size() + allocate(surfaces(i)%obj%neighbor_pos(j)) + surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:j) + + ! Copy negative neighbors to Surface instance + j = neighbor_neg(i)%size() + allocate(surfaces(i)%obj%neighbor_neg(j)) + surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:j) end do - count_positive = 0 - count_negative = 0 - - ! loop over all cells - do i = 1, n_cells - c => cells(i) - - ! loop through the region specification - do j = 1, size(c%region) - i_surface = c % region(j) - positive = (i_surface > 0) - - ! Skip any tokens that correspond to operators rather than regions - i_surface = abs(i_surface) - if (i_surface >= OP_UNION) cycle - - if (positive) then - count_positive(i_surface) = count_positive(i_surface) + 1 - surfaces(i_surface)%obj%neighbor_pos(count_positive(i_surface)) = i - else - count_negative(i_surface) = count_negative(i_surface) + 1 - surfaces(i_surface)%obj%neighbor_neg(count_negative(i_surface)) = i - end if - end do - end do - - deallocate(count_positive) - deallocate(count_negative) - end subroutine neighbor_lists !=============================================================================== From a3199df9d96d338d245458fed947e41a016ff26d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 7 Oct 2015 22:21:02 +0700 Subject: [PATCH 09/19] Make Nuclide%reactions allocatable by using associate constructs --- src/ace.F90 | 341 +++++++++++++++++++++--------------------- src/ace_header.F90 | 5 +- src/cross_section.F90 | 14 +- src/output.F90 | 47 +++--- src/physics.F90 | 60 ++++---- src/tally.F90 | 197 ++++++++++++------------ 6 files changed, 324 insertions(+), 340 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index dbed66757..4cd53a251 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -715,7 +715,6 @@ contains integer :: IE ! reaction's starting index on energy grid integer :: NE ! number of energies integer :: NR ! number of interpolation regions - type(Reaction), pointer :: rxn type(ListInt) :: MTs LMT = JXS(3) @@ -733,14 +732,15 @@ contains ! Store elastic scattering cross-section on reaction one -- note that the ! sigma array is not allocated or stored for elastic scattering since it is ! already stored in nuc % elastic - rxn => nuc % reactions(1) - rxn % MT = 2 - rxn % Q_value = ZERO - rxn % multiplicity = 1 - rxn % threshold = 1 - rxn % scatter_in_cm = .true. - rxn % has_angle_dist = .false. - rxn % has_energy_dist = .false. + associate (rxn => nuc % reactions(1)) + rxn % MT = 2 + rxn % Q_value = ZERO + rxn % multiplicity = 1 + rxn % threshold = 1 + rxn % scatter_in_cm = .true. + rxn % has_angle_dist = .false. + rxn % has_energy_dist = .false. + end associate ! Add contribution of elastic scattering to total cross section nuc % total = nuc % total + nuc % elastic @@ -753,64 +753,64 @@ contains i_fission = 0 do i = 1, NMT - rxn => nuc % reactions(i+1) + associate (rxn => nuc % reactions(i+1)) + ! set defaults + rxn % has_angle_dist = .false. + rxn % has_energy_dist = .false. - ! set defaults - rxn % has_angle_dist = .false. - rxn % has_energy_dist = .false. + ! read MT number, Q-value, and neutrons produced + rxn % MT = int(XSS(LMT + i - 1)) + rxn % Q_value = XSS(JXS4 + i - 1) + rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1))) + rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0) - ! read MT number, Q-value, and neutrons produced - rxn % MT = int(XSS(LMT + i - 1)) - rxn % Q_value = XSS(JXS4 + i - 1) - rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1))) - rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0) + ! Read energy-dependent multiplicities + if (rxn % multiplicity > 100) then + ! Set flag and allocate space for Tab1 to store yield + rxn % multiplicity_with_E = .true. + allocate(rxn % multiplicity_E) - ! Read energy-dependent multiplicities - if (rxn % multiplicity > 100) then - ! Set flag and allocate space for Tab1 to store yield - rxn % multiplicity_with_E = .true. - allocate(rxn % multiplicity_E) + XSS_index = JXS(11) + rxn % multiplicity - 101 + NR = nint(XSS(XSS_index)) + rxn % multiplicity_E % n_regions = NR - XSS_index = JXS(11) + rxn % multiplicity - 101 - NR = nint(XSS(XSS_index)) - rxn % multiplicity_E % n_regions = NR + ! allocate space for ENDF interpolation parameters + if (NR > 0) then + allocate(rxn % multiplicity_E % nbt(NR)) + allocate(rxn % multiplicity_E % int(NR)) + end if - ! allocate space for ENDF interpolation parameters - if (NR > 0) then - allocate(rxn % multiplicity_E % nbt(NR)) - allocate(rxn % multiplicity_E % int(NR)) + ! read ENDF interpolation parameters + XSS_index = XSS_index + 1 + if (NR > 0) then + rxn % multiplicity_E % nbt = get_int(NR) + rxn % multiplicity_E % int = get_int(NR) + end if + + ! allocate space for yield data + XSS_index = XSS_index + 2*NR + NE = nint(XSS(XSS_index)) + rxn % multiplicity_E % n_pairs = NE + allocate(rxn % multiplicity_E % x(NE)) + allocate(rxn % multiplicity_E % y(NE)) + + ! read yield data + XSS_index = XSS_index + 1 + rxn % multiplicity_E % x = get_real(NE) + rxn % multiplicity_E % y = get_real(NE) end if - ! read ENDF interpolation parameters - XSS_index = XSS_index + 1 - if (NR > 0) then - rxn % multiplicity_E % nbt = get_int(NR) - rxn % multiplicity_E % int = get_int(NR) - end if + ! read starting energy index + LOCA = int(XSS(LXS + i - 1)) + IE = int(XSS(JXS7 + LOCA - 1)) + rxn % threshold = IE - ! allocate space for yield data - XSS_index = XSS_index + 2*NR - NE = nint(XSS(XSS_index)) - rxn % multiplicity_E % n_pairs = NE - allocate(rxn % multiplicity_E % x(NE)) - allocate(rxn % multiplicity_E % y(NE)) - - ! read yield data - XSS_index = XSS_index + 1 - rxn % multiplicity_E % x = get_real(NE) - rxn % multiplicity_E % y = get_real(NE) - end if - - ! read starting energy index - LOCA = int(XSS(LXS + i - 1)) - IE = int(XSS(JXS7 + LOCA - 1)) - rxn % threshold = IE - - ! read number of energies cross section values - NE = int(XSS(JXS7 + LOCA)) - allocate(rxn % sigma(NE)) - XSS_index = JXS7 + LOCA + 1 - rxn % sigma = get_real(NE) + ! read number of energies cross section values + NE = int(XSS(JXS7 + LOCA)) + allocate(rxn % sigma(NE)) + XSS_index = JXS7 + LOCA + 1 + rxn % sigma = get_real(NE) + end associate end do ! Create set of MT values @@ -821,56 +821,57 @@ contains ! Create total, absorption, and fission cross sections do i = 2, size(nuc % reactions) - rxn => nuc % reactions(i) - IE = rxn % threshold - NE = size(rxn % sigma) + associate (rxn => nuc % reactions(i)) + IE = rxn % threshold + NE = size(rxn % sigma) - ! Skip total inelastic level scattering, gas production cross sections - ! (MT=200+), etc. - if (rxn % MT == N_LEVEL) cycle - if (rxn % MT > N_5N2P .and. rxn % MT < N_P0) cycle + ! Skip total inelastic level scattering, gas production cross sections + ! (MT=200+), etc. + if (rxn % MT == N_LEVEL) cycle + if (rxn % MT > N_5N2P .and. rxn % MT < N_P0) cycle - ! Skip level cross sections if total is available - if (rxn % MT >= N_P0 .and. rxn % MT <= N_PC .and. MTs % contains(N_P)) cycle - if (rxn % MT >= N_D0 .and. rxn % MT <= N_DC .and. MTs % contains(N_D)) cycle - if (rxn % MT >= N_T0 .and. rxn % MT <= N_TC .and. MTs % contains(N_T)) cycle - if (rxn % MT >= N_3HE0 .and. rxn % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle - if (rxn % MT >= N_A0 .and. rxn % MT <= N_AC .and. MTs % contains(N_A)) cycle - if (rxn % MT >= N_2N0 .and. rxn % MT <= N_2NC .and. MTs % contains(N_2N)) cycle + ! Skip level cross sections if total is available + if (rxn % MT >= N_P0 .and. rxn % MT <= N_PC .and. MTs % contains(N_P)) cycle + if (rxn % MT >= N_D0 .and. rxn % MT <= N_DC .and. MTs % contains(N_D)) cycle + if (rxn % MT >= N_T0 .and. rxn % MT <= N_TC .and. MTs % contains(N_T)) cycle + if (rxn % MT >= N_3HE0 .and. rxn % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle + if (rxn % MT >= N_A0 .and. rxn % MT <= N_AC .and. MTs % contains(N_A)) cycle + if (rxn % MT >= N_2N0 .and. rxn % MT <= N_2NC .and. MTs % contains(N_2N)) cycle - ! Add contribution to total cross section - nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma + ! Add contribution to total cross section + nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma - ! Add contribution to absorption cross section - if (is_disappearance(rxn % MT)) then - nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma - end if + ! Add contribution to absorption cross section + if (is_disappearance(rxn % MT)) then + nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma + end if - ! Information about fission reactions - if (rxn % MT == N_FISSION) then - allocate(nuc % index_fission(1)) - elseif (rxn % MT == N_F) then - allocate(nuc % index_fission(PARTIAL_FISSION_MAX)) - nuc % has_partial_fission = .true. - end if + ! Information about fission reactions + if (rxn % MT == N_FISSION) then + allocate(nuc % index_fission(1)) + elseif (rxn % MT == N_F) then + allocate(nuc % index_fission(PARTIAL_FISSION_MAX)) + nuc % has_partial_fission = .true. + end if - ! Add contribution to fission cross section - if (is_fission(rxn % MT)) then - nuc % fissionable = .true. - nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma + ! Add contribution to fission cross section + if (is_fission(rxn % MT)) then + nuc % fissionable = .true. + nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma - ! Also need to add fission cross sections to absorption - nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma + ! Also need to add fission cross sections to absorption + nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma - ! If total fission reaction is present, there's no need to store the - ! reaction cross-section since it was copied to nuc % fission - if (rxn % MT == N_FISSION) deallocate(rxn % sigma) + ! If total fission reaction is present, there's no need to store the + ! reaction cross-section since it was copied to nuc % fission + if (rxn % MT == N_FISSION) deallocate(rxn % sigma) - ! Keep track of this reaction for easy searching later - i_fission = i_fission + 1 - nuc % index_fission(i_fission) = i - nuc % n_fission = nuc % n_fission + 1 - end if + ! Keep track of this reaction for easy searching later + i_fission = i_fission + 1 + nuc % index_fission(i_fission) = i + nuc % n_fission = nuc % n_fission + 1 + end if + end associate end do ! Clear MTs set @@ -895,7 +896,6 @@ contains integer :: i ! index in reactions array integer :: j ! index over incoming energies integer :: length ! length of data array to allocate - type(Reaction), pointer :: rxn JXS8 = JXS(8) JXS9 = JXS(9) @@ -903,71 +903,72 @@ contains ! loop over all reactions with secondary neutrons -- NXS(5) does not include ! elastic scattering do i = 1, NXS(5) + 1 - rxn => nuc%reactions(i) + associate (rxn => nuc%reactions(i)) - ! find location of angular distribution - LOCB = int(XSS(JXS8 + i - 1)) - if (LOCB == -1) then - ! Angular distribution data are specified through LAWi = 44 in the DLW - ! block - cycle - elseif (LOCB == 0) then - ! No angular distribution data are given for this reaction, isotropic - ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0) - cycle - end if - rxn % has_angle_dist = .true. - - ! allocate space for incoming energies and locations - NE = int(XSS(JXS9 + LOCB - 1)) - rxn % adist % n_energy = NE - allocate(rxn % adist % energy(NE)) - allocate(rxn % adist % type(NE)) - allocate(rxn % adist % location(NE)) - - ! read incoming energy grid and location of nucs - XSS_index = JXS9 + LOCB - rxn % adist % energy = get_real(NE) - rxn % adist % location = get_int(NE) - - ! determine dize of data block - length = 0 - do j = 1, NE - LC = rxn % adist % location(j) - if (LC == 0) then - ! isotropic - rxn % adist % type(j) = ANGLE_ISOTROPIC - elseif (LC > 0) then - ! 32 equiprobable bins - rxn % adist % type(j) = ANGLE_32_EQUI - length = length + 33 - elseif (LC < 0) then - ! tabular distribution - rxn % adist % type(j) = ANGLE_TABULAR - NP = int(XSS(JXS9 + abs(LC))) - length = length + 2 + 3*NP + ! find location of angular distribution + LOCB = int(XSS(JXS8 + i - 1)) + if (LOCB == -1) then + ! Angular distribution data are specified through LAWi = 44 in the DLW + ! block + cycle + elseif (LOCB == 0) then + ! No angular distribution data are given for this reaction, isotropic + ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0) + cycle end if - end do + rxn % has_angle_dist = .true. - ! allocate angular distribution data and read - allocate(rxn % adist % data(length)) + ! allocate space for incoming energies and locations + NE = int(XSS(JXS9 + LOCB - 1)) + rxn % adist % n_energy = NE + allocate(rxn % adist % energy(NE)) + allocate(rxn % adist % type(NE)) + allocate(rxn % adist % location(NE)) - ! read angular distribution -- currently this does not actually parse the - ! angular distribution tables for each incoming energy, that must be done - ! on-the-fly - XSS_index = JXS9 + LOCB + 2 * NE - rxn % adist % data = get_real(length) + ! read incoming energy grid and location of nucs + XSS_index = JXS9 + LOCB + rxn % adist % energy = get_real(NE) + rxn % adist % location = get_int(NE) - ! change location pointers since they are currently relative to JXS(9) - LC = LOCB + 2 * NE + 1 - do j = 1, NE - ! For consistency, leave location as 0 if type is isotropic. - ! This is not necessary for current correctness, but can avoid - ! future issues - if (rxn % adist % location(j) /= 0) then - rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC - end if - end do + ! determine dize of data block + length = 0 + do j = 1, NE + LC = rxn % adist % location(j) + if (LC == 0) then + ! isotropic + rxn % adist % type(j) = ANGLE_ISOTROPIC + elseif (LC > 0) then + ! 32 equiprobable bins + rxn % adist % type(j) = ANGLE_32_EQUI + length = length + 33 + elseif (LC < 0) then + ! tabular distribution + rxn % adist % type(j) = ANGLE_TABULAR + NP = int(XSS(JXS9 + abs(LC))) + length = length + 2 + 3*NP + end if + end do + + ! allocate angular distribution data and read + allocate(rxn % adist % data(length)) + + ! read angular distribution -- currently this does not actually parse the + ! angular distribution tables for each incoming energy, that must be done + ! on-the-fly + XSS_index = JXS9 + LOCB + 2 * NE + rxn % adist % data = get_real(length) + + ! change location pointers since they are currently relative to JXS(9) + LC = LOCB + 2 * NE + 1 + do j = 1, NE + ! For consistency, leave location as 0 if type is isotropic. + ! This is not necessary for current correctness, but can avoid + ! future issues + if (rxn % adist % location(j) /= 0) then + rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC + end if + end do + end associate end do end subroutine read_angular_dist @@ -983,23 +984,23 @@ contains integer :: LED ! location of energy distribution locators integer :: LOCC ! location of energy distributions for given MT integer :: i ! loop index - type(Reaction), pointer :: rxn LED = JXS(10) ! Loop over all reactions do i = 1, NXS(5) - rxn => nuc % reactions(i+1) ! skip over elastic scattering - rxn % has_energy_dist = .true. + associate (rxn => nuc % reactions(i+1)) ! skip over elastic scattering + rxn % has_energy_dist = .true. - ! find location of energy distribution data - LOCC = int(XSS(LED + i - 1)) + ! find location of energy distribution data + LOCC = int(XSS(LED + i - 1)) - ! allocate energy distribution - allocate(rxn % edist) + ! allocate energy distribution + allocate(rxn % edist) - ! read data for energy distribution - call get_energy_dist(rxn % edist, LOCC) + ! read data for energy distribution + call get_energy_dist(rxn % edist, LOCC) + end associate end do end subroutine read_energy_dist diff --git a/src/ace_header.F90 b/src/ace_header.F90 index a10c6fe87..9eadca594 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -146,7 +146,7 @@ module ace_header ! Reactions integer :: n_reaction ! # of reactions - type(Reaction), pointer :: reactions(:) => null() + type(Reaction), allocatable :: reactions(:) type(DictIntInt) :: reaction_index ! map MT values to index in reactions ! array; used at tally-time @@ -341,11 +341,10 @@ module ace_header end if if (associated(this % urr_data)) then - call this % urr_data % clear() deallocate(this % urr_data) end if - if (associated(this % reactions)) then + if (allocated(this % reactions)) then do i = 1, size(this % reactions) call this % reactions(i) % clear() end do diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b9c76b503..03bcefca8 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -379,7 +379,6 @@ contains logical :: same_nuc ! do we know the xs for this nuclide at this energy? type(UrrData), pointer :: urr type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn micro_xs(i_nuclide) % use_ptable = .true. @@ -475,18 +474,17 @@ contains ! Determine treatment of inelastic scattering inelastic = ZERO if (urr % inelastic_flag > 0) then - ! Get pointer to inelastic scattering reaction - rxn => nuc % reactions(nuc % urr_inelastic) - ! Get index on energy grid and interpolation factor i_energy = micro_xs(i_nuclide) % index_grid f = micro_xs(i_nuclide) % interp_factor ! Determine inelastic scattering cross section - if (i_energy >= rxn % threshold) then - inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + & - f * rxn % sigma(i_energy - rxn%threshold + 2) - end if + associate (rxn => nuc % reactions(nuc % urr_inelastic)) + if (i_energy >= rxn % threshold) then + inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + & + f * rxn % sigma(i_energy - rxn%threshold + 2) + end if + end associate end if ! Multiply by smooth cross-section if needed diff --git a/src/output.F90 b/src/output.F90 index 2aa90985a..19ae6268b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -330,7 +330,6 @@ contains integer :: size_energy ! memory used for a energy distributions (bytes) integer :: size_urr ! memory used for probability tables (bytes) character(11) :: law ! secondary energy distribution law - type(Reaction), pointer :: rxn type(UrrData), pointer :: urr ! set default unit for writing information @@ -359,32 +358,32 @@ contains ! Information on each reaction write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)' do i = 1, nuc % n_reaction - rxn => nuc % reactions(i) + associate (rxn => nuc % reactions(i)) + ! Determine size of angle distribution + if (rxn % has_angle_dist) then + size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 + else + size_angle = 0 + end if - ! Determine size of angle distribution - if (rxn % has_angle_dist) then - size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 - else - size_angle = 0 - end if + ! Determine size of energy distribution and law + if (rxn % has_energy_dist) then + size_energy = size(rxn % edist % data) * 8 + law = to_str(rxn % edist % law) + else + size_energy = 0 + law = 'None' + end if - ! Determine size of energy distribution and law - if (rxn % has_energy_dist) then - size_energy = size(rxn % edist % data) * 8 - law = to_str(rxn % edist % law) - else - size_energy = 0 - law = 'None' - end if + write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') & + reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & + law(1:4), rxn % threshold, size_angle, size_energy - write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') & - reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & - law(1:4), rxn % threshold, size_angle, size_energy - - ! Accumulate data size - size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8 - size_angle_total = size_angle_total + size_angle - size_energy_total = size_energy_total + size_energy + ! Accumulate data size + size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8 + size_angle_total = size_angle_total + size_angle + size_energy_total = size_energy_total + size_energy + end associate end do ! Add memory required for summary reactions (total, absorption, fission, diff --git a/src/physics.F90 b/src/physics.F90 index c06919758..a01a3a30c 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -185,7 +185,6 @@ contains !=============================================================================== subroutine sample_fission(i_nuclide, i_reaction) - integer, intent(in) :: i_nuclide ! index in nuclides array integer, intent(out) :: i_reaction ! index in nuc % reactions array @@ -195,7 +194,6 @@ contains real(8) :: prob real(8) :: cutoff type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn ! Get pointer to nuclide nuc => nuclides(i_nuclide) @@ -220,14 +218,15 @@ contains FISSION_REACTION_LOOP: do i = 1, nuc % n_fission i_reaction = nuc % index_fission(i) - rxn => nuc % reactions(i_reaction) - ! if energy is below threshold for this reaction, skip it - if (i_grid < rxn % threshold) cycle + associate (rxn => nuc % reactions(i_reaction)) + ! if energy is below threshold for this reaction, skip it + if (i_grid < rxn % threshold) cycle - ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & - + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + ! add to cumulative probability + prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + end associate ! Create fission bank sites if fission occurs if (prob > cutoff) exit FISSION_REACTION_LOOP @@ -312,11 +311,10 @@ contains real(8) :: f real(8) :: prob real(8) :: cutoff - type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering real(8) :: phi ! azimuthal angle for iso-in-lab scattering + type(Nuclide), pointer :: nuc ! copy incoming direction uvw_old(:) = p % coord(1) % uvw @@ -343,11 +341,8 @@ contains p % E, p % coord(1) % uvw, p % mu) else - ! get pointer to elastic scattering reaction - rxn => nuc % reactions(1) - ! Perform collision physics for elastic scattering - call elastic_scatter(i_nuclide, rxn, & + call elastic_scatter(i_nuclide, nuc % reactions(1), & p % E, p % coord(1) % uvw, p % mu, p % wgt) end if @@ -370,28 +365,28 @@ contains &// trim(nuc % name)) end if - rxn => nuc % reactions(i) + associate (rxn => nuc % reactions(i)) + ! Skip fission reactions + if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF & + .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle - ! Skip fission reactions - if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF & - .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle + ! some materials have gas production cross sections with MT > 200 that + ! are duplicates. Also MT=4 is total level inelastic scattering which + ! should be skipped + if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle - ! some materials have gas production cross sections with MT > 200 that - ! are duplicates. Also MT=4 is total level inelastic scattering which - ! should be skipped - if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle + ! if energy is below threshold for this reaction, skip it + if (i_grid < rxn % threshold) cycle - ! if energy is below threshold for this reaction, skip it - if (i_grid < rxn % threshold) cycle - - ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & - + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + ! add to cumulative probability + prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + end associate end do ! Perform collision physics for inelastic scattering - call inelastic_scatter(nuc, rxn, p) - p % event_MT = rxn % MT + call inelastic_scatter(nuc, nuc%reactions(i), p) + p % event_MT = nuc%reactions(i)%MT end if @@ -1090,11 +1085,9 @@ contains real(8) :: weight ! weight adjustment for ufs method logical :: in_mesh ! source site in ufs mesh? type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn ! Get pointers nuc => nuclides(i_nuclide) - rxn => nuc % reactions(i_reaction) ! TODO: Heat generation from fission @@ -1165,7 +1158,8 @@ contains ! Sample secondary energy distribution for fission reaction and set energy ! in fission bank - fission_bank(i) % E = sample_fission_energy(nuc, rxn, p) + fission_bank(i) % E = sample_fission_energy(nuc, nuc%reactions(& + i_reaction), p) ! Set the delayed group of the neutron fission_bank(i) % delayed_group = p % delayed_group diff --git a/src/tally.F90 b/src/tally.F90 index 2f178fe46..4790a9900 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -66,9 +66,6 @@ contains real(8) :: macro_scatt ! material macro scatt xs real(8) :: uvw(3) ! particle direction real(8) :: E ! particle energy - type(Material), pointer :: mat - type(Reaction), pointer :: rxn - type(Nuclide), pointer :: nuc i = 0 SCORE_LOOP: do q = 1, t % n_user_score_bins @@ -220,24 +217,20 @@ contains ! of one. score = p % last_wgt else - do m = 1, nuclides(p % event_nuclide) % n_reaction - ! Check if this is the desired MT - if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then - ! Found the reaction, set our pointer and move on with life - rxn => nuclides(p % event_nuclide) % reactions(m) - exit - end if - end do + m = nuclides(p%event_nuclide)%reaction_index% & + get_key(p % event_MT) ! Get multiplicity and apply to score - if (rxn % multiplicity_with_E) then - ! Then the multiplicity was already incorporated in to p % wgt - ! per the scattering routine, - score = p % wgt - else - ! Grab the multiplicity from the rxn - score = p % last_wgt * rxn % multiplicity - end if + associate (rxn => nuclides(p%event_nuclide)%reactions(m)) + if (rxn % multiplicity_with_E) then + ! Then the multiplicity was already incorporated in to p % wgt + ! per the scattering routine, + score = p % wgt + else + ! Grab the multiplicity from the rxn + score = p % last_wgt * rxn % multiplicity + end if + end associate end if @@ -257,24 +250,20 @@ contains ! of one. score = p % last_wgt else - do m = 1, nuclides(p % event_nuclide) % n_reaction - ! Check if this is the desired MT - if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then - ! Found the reaction, set our pointer and move on with life - rxn => nuclides(p % event_nuclide) % reactions(m) - exit - end if - end do + m = nuclides(p%event_nuclide)%reaction_index% & + get_key(p % event_MT) ! Get multiplicity and apply to score - if (rxn % multiplicity_with_E) then - ! Then the multiplicity was already incorporated in to p % wgt - ! per the scattering routine, - score = p % wgt - else - ! Grab the multiplicity from the rxn - score = p % last_wgt * rxn % multiplicity - end if + associate (rxn => nuclides(p%event_nuclide)%reactions(m)) + if (rxn % multiplicity_with_E) then + ! Then the multiplicity was already incorporated in to p % wgt + ! per the scattering routine, + score = p % wgt + else + ! Grab the multiplicity from the rxn + score = p % last_wgt * rxn % multiplicity + end if + end associate end if @@ -294,24 +283,20 @@ contains ! of one. score = p % last_wgt else - do m = 1, nuclides(p % event_nuclide) % n_reaction - ! Check if this is the desired MT - if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then - ! Found the reaction, set our pointer and move on with life - rxn => nuclides(p % event_nuclide) % reactions(m) - exit - end if - end do + m = nuclides(p%event_nuclide)%reaction_index% & + get_key(p % event_MT) ! Get multiplicity and apply to score - if (rxn % multiplicity_with_E) then - ! Then the multiplicity was already incorporated in to p % wgt - ! per the scattering routine, - score = p % wgt - else - ! Grab the multiplicity from the rxn - score = p % last_wgt * rxn % multiplicity - end if + associate (rxn => nuclides(p%event_nuclide)%reactions(m)) + if (rxn % multiplicity_with_E) then + ! Then the multiplicity was already incorporated in to p % wgt + ! per the scattering routine, + score = p % wgt + else + ! Grab the multiplicity from the rxn + score = p % last_wgt * rxn % multiplicity + end if + end associate end if @@ -466,9 +451,6 @@ contains ! delayed-nu-fission if (micro_xs(p % event_nuclide) % absorption > ZERO) then - ! Get the event nuclide - nuc => nuclides(p % event_nuclide) - ! Check if the delayed group filter is present if (dg_filter > 0) then @@ -480,11 +462,11 @@ contains d = t % filters(dg_filter) % int_bins(d_bin) ! Compute the yield for this delayed group - yield = yield_delayed(nuc, E, d) + yield = yield_delayed(nuclides(p % event_nuclide), E, d) ! Compute the score and tally to bin score = p % absorb_wgt * yield * micro_xs(p % event_nuclide) & - % fission * nu_delayed(nuc, E) / & + % fission * nu_delayed(nuclides(p % event_nuclide), E) / & micro_xs(p % event_nuclide) % absorption call score_fission_delayed_dg(t, d_bin, score, score_index) end do @@ -494,7 +476,7 @@ contains ! by multiplying the absorbed weight by the fraction of the ! delayed-nu-fission xs to the absorption xs score = p % absorb_wgt * micro_xs(p % event_nuclide) & - % fission * nu_delayed(nuc, E) / & + % fission * nu_delayed(nuclides(p % event_nuclide), E) / & micro_xs(p % event_nuclide) % absorption end if end if @@ -535,9 +517,6 @@ contains ! Check if tally is on a single nuclide if (i_nuclide > 0) then - ! Get the nuclide of interest - nuc => nuclides(i_nuclide) - ! Check if the delayed group filter is present if (dg_filter > 0) then @@ -548,11 +527,19 @@ contains d = t % filters(dg_filter) % int_bins(d_bin) ! Compute the yield for this delayed group +<<<<<<< HEAD yield = yield_delayed(nuc, E, d) ! Compute the score and tally to bin score = micro_xs(i_nuclide) % fission * yield & * nu_delayed(nuc, E) * atom_density * flux +======= + yield = yield_delayed(nuclides(i_nuclide), p % E, d) + + ! Compute the score and tally to bin + score = micro_xs(i_nuclide) % fission * yield & + * nu_delayed(nuclides(i_nuclide), p % E) * atom_density * flux +>>>>>>> Make Nuclide%reactions allocatable by using associate constructs call score_fission_delayed_dg(t, d_bin, score, score_index) end do cycle SCORE_LOOP @@ -560,27 +547,29 @@ contains ! If the delayed group filter is not present, compute the score ! by multiplying the delayed-nu-fission macro xs by the flux +<<<<<<< HEAD score = micro_xs(i_nuclide) % fission * nu_delayed(nuc, E)& * atom_density * flux +======= + score = micro_xs(i_nuclide) % fission * & + nu_delayed(nuclides(i_nuclide), p % E) * atom_density * flux +>>>>>>> Make Nuclide%reactions allocatable by using associate constructs end if ! Tally is on total nuclides else - ! Get pointer to current material - mat => materials(p % material) - ! Check if the delayed group filter is present if (dg_filter > 0) then ! Loop over all nuclides in the current material - do l = 1, mat % n_nuclides + do l = 1, materials(p % material) % n_nuclides ! Get atom density - atom_density_ = mat % atom_density(l) + atom_density_ = materials(p % material) % atom_density(l) ! Get index in nuclides array - i_nuc = mat % nuclide(l) + i_nuc = materials(p % material) % nuclide(l) ! Loop over all delayed group bins and tally to them individually do d_bin = 1, t % filters(dg_filter) % n_bins @@ -588,15 +577,20 @@ contains ! Get the delayed group for this bin d = t % filters(dg_filter) % int_bins(d_bin) - ! Get the current nuclide - nuc => nuclides(i_nuc) - ! Get the yield for the desired nuclide and delayed group +<<<<<<< HEAD yield = yield_delayed(nuc, E, d) ! Compute the score and tally to bin score = micro_xs(i_nuc) % fission * yield & * nu_delayed(nuc, E) * atom_density_ * flux +======= + yield = yield_delayed(nuclides(i_nuc), p % E, d) + + ! Compute the score and tally to bin + score = micro_xs(i_nuc) % fission * yield & + * nu_delayed(nuclides(i_nuc), p % E) * atom_density_ * flux +>>>>>>> Make Nuclide%reactions allocatable by using associate constructs call score_fission_delayed_dg(t, d_bin, score, score_index) end do end do @@ -606,13 +600,13 @@ contains score = ZERO ! Loop over all nuclides in the current material - do l = 1, mat % n_nuclides + do l = 1, materials(p % material) % n_nuclides ! Get atom density - atom_density_ = mat % atom_density(l) + atom_density_ = materials(p % material) % atom_density(l) ! Get index in nuclides array - i_nuc = mat % nuclide(l) + i_nuc = materials(p % material) % nuclide(l) ! Accumulate the contribution from each nuclide score = score + micro_xs(i_nuc) % fission & @@ -693,42 +687,41 @@ contains if (i_nuclide > 0) then if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) - rxn => nuclides(i_nuclide) % reactions(m) - - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - if (i_energy >= rxn % threshold) then - score = ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - end if - - else - ! Get pointer to current material - mat => materials(p % material) - do l = 1, mat % n_nuclides - ! Get atom density - atom_density_ = mat % atom_density(l) - - ! Get index in nuclides array - i_nuc = mat % nuclide(l) - - if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then - m = nuclides(i_nuc)%reaction_index%get_key(score_bin) - rxn => nuclides(i_nuc) % reactions(m) + associate (rxn => nuclides(i_nuclide) % reactions(m)) ! Retrieve index on nuclide energy grid and interpolation ! factor - i_energy = micro_xs(i_nuc) % index_grid - f = micro_xs(i_nuc) % interp_factor + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor if (i_energy >= rxn % threshold) then - score = score + ((ONE - f) * rxn % sigma(i_energy - & + score = ((ONE - f) * rxn % sigma(i_energy - & rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density_ * flux + rxn%threshold + 2)) * atom_density * flux end if + end associate + end if + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then + m = nuclides(i_nuc)%reaction_index%get_key(score_bin) + associate (rxn => nuclides(i_nuc) % reactions(m)) + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuc) % index_grid + f = micro_xs(i_nuc) % interp_factor + if (i_energy >= rxn % threshold) then + score = score + ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density_ * flux + end if + end associate end if end do end if From 441fd4f00dfb3cd6f79abc0ad2887b04dd5dbfd8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 30 Oct 2015 15:47:01 -0500 Subject: [PATCH 10/19] Don't pre-compute kappa-fission cross sections. This also fixes a bug in the kappa-fission score. Before, kappa-fission was computed as Q*fission, but this was done before URR cross sections were determined. Thus, if fission changed in calculate_urr_xs, this wasn't reflected in the kappa-fission score. Now, since it is all done at tally-time, there is no inconsistency. --- src/ace_header.F90 | 2 - src/cross_section.F90 | 13 --- src/tally.F90 | 86 ++++++++++--------- tests/test_many_scores/results_true.dat | 4 +- .../test_score_kappafission/results_true.dat | 20 ++--- 5 files changed, 59 insertions(+), 66 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 9eadca594..6c27747d8 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -258,7 +258,6 @@ module ace_header real(8) :: absorption ! microscopic absorption xs real(8) :: fission ! microscopic fission xs real(8) :: nu_fission ! microscopic production xs - real(8) :: kappa_fission ! microscopic energy-released from fission ! Information for S(a,b) use integer :: index_sab ! index in sab_tables (zero means no table) @@ -281,7 +280,6 @@ module ace_header real(8) :: absorption ! macroscopic absorption xs real(8) :: fission ! macroscopic fission xs real(8) :: nu_fission ! macroscopic production xs - real(8) :: kappa_fission ! macroscopic energy-released from fission end type MaterialMacroXS contains diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 03bcefca8..969a397f0 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -41,7 +41,6 @@ contains material_xs % absorption = ZERO material_xs % fission = ZERO material_xs % nu_fission = ZERO - material_xs % kappa_fission = ZERO ! Exit subroutine if material is void if (p % material == MATERIAL_VOID) return @@ -125,10 +124,6 @@ contains ! Add contributions to material macroscopic nu-fission cross section material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission - - ! Add contributions to material macroscopic energy release from fission - material_xs % kappa_fission = material_xs % kappa_fission + & - atom_density * micro_xs(i_nuclide) % kappa_fission end do end subroutine calculate_xs @@ -216,7 +211,6 @@ contains ! Initialize nuclide cross-sections to zero micro_xs(i_nuclide) % fission = ZERO micro_xs(i_nuclide) % nu_fission = ZERO - micro_xs(i_nuclide) % kappa_fission = ZERO ! Calculate microscopic nuclide total cross section micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) & @@ -238,13 +232,6 @@ contains ! Calculate microscopic nuclide nu-fission cross section micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( & i_grid) + f * nuc % nu_fission(i_grid+1) - - ! Calculate microscopic nuclide kappa-fission cross section - ! The ENDF standard (ENDF-102) states that MT 18 stores - ! the fission energy as the Q_value (fission(1)) - micro_xs(i_nuclide) % kappa_fission = & - nuc % reactions(nuc % index_fission(1)) % Q_value * & - micro_xs(i_nuclide) % fission end if ! If there is S(a,b) data for this nuclide, we need to do a few diff --git a/src/tally.F90 b/src/tally.F90 index 4790a9900..c4ca2028b 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -527,19 +527,11 @@ contains d = t % filters(dg_filter) % int_bins(d_bin) ! Compute the yield for this delayed group -<<<<<<< HEAD - yield = yield_delayed(nuc, E, d) + yield = yield_delayed(nuclides(i_nuclide), E, d) ! Compute the score and tally to bin score = micro_xs(i_nuclide) % fission * yield & - * nu_delayed(nuc, E) * atom_density * flux -======= - yield = yield_delayed(nuclides(i_nuclide), p % E, d) - - ! Compute the score and tally to bin - score = micro_xs(i_nuclide) % fission * yield & - * nu_delayed(nuclides(i_nuclide), p % E) * atom_density * flux ->>>>>>> Make Nuclide%reactions allocatable by using associate constructs + * nu_delayed(nuclides(i_nuclide), E) * atom_density * flux call score_fission_delayed_dg(t, d_bin, score, score_index) end do cycle SCORE_LOOP @@ -547,13 +539,8 @@ contains ! If the delayed group filter is not present, compute the score ! by multiplying the delayed-nu-fission macro xs by the flux -<<<<<<< HEAD - score = micro_xs(i_nuclide) % fission * nu_delayed(nuc, E)& - * atom_density * flux -======= score = micro_xs(i_nuclide) % fission * & - nu_delayed(nuclides(i_nuclide), p % E) * atom_density * flux ->>>>>>> Make Nuclide%reactions allocatable by using associate constructs + nu_delayed(nuclides(i_nuclide), E) * atom_density * flux end if ! Tally is on total nuclides @@ -578,19 +565,11 @@ contains d = t % filters(dg_filter) % int_bins(d_bin) ! Get the yield for the desired nuclide and delayed group -<<<<<<< HEAD - yield = yield_delayed(nuc, E, d) + yield = yield_delayed(nuclides(i_nuc), E, d) ! Compute the score and tally to bin score = micro_xs(i_nuc) % fission * yield & - * nu_delayed(nuc, E) * atom_density_ * flux -======= - yield = yield_delayed(nuclides(i_nuc), p % E, d) - - ! Compute the score and tally to bin - score = micro_xs(i_nuc) % fission * yield & - * nu_delayed(nuclides(i_nuc), p % E) * atom_density_ * flux ->>>>>>> Make Nuclide%reactions allocatable by using associate constructs + * nu_delayed(nuclides(i_nuc), E) * atom_density_ * flux call score_fission_delayed_dg(t, d_bin, score, score_index) end do end do @@ -618,38 +597,67 @@ contains case (SCORE_KAPPA_FISSION) + ! Determine kappa-fission cross section on the fly. The ENDF standard + ! (ENDF-102) states that MT 18 stores the fission energy as the Q_value + ! (fission(1)) + + score = ZERO + if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing) then ! No fission events occur if survival biasing is on -- need to ! calculate fraction of absorptions that would have resulted in ! fission scale by kappa-fission - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * & - micro_xs(p % event_nuclide) % kappa_fission / & - micro_xs(p % event_nuclide) % absorption - else - score = ZERO - end if + associate (nuc => nuclides(p % event_nuclide)) + if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + nuc % fissionable) then + score = p % absorb_wgt * & + nuc%reactions(nuc%index_fission(1))%Q_value * & + micro_xs(p % event_nuclide) % fission / & + micro_xs(p % event_nuclide) % absorption + end if + end associate else ! Skip any non-absorption events if (p % event == EVENT_SCATTER) cycle SCORE_LOOP ! All fission events will contribute, so again we can use ! particle's weight entering the collision as the estimate for ! the fission energy production rate - score = p % last_wgt * & - micro_xs(p % event_nuclide) % kappa_fission / & - micro_xs(p % event_nuclide) % absorption + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % fissionable) then + score = p % last_wgt * & + nuc%reactions(nuc%index_fission(1))%Q_value * & + micro_xs(p % event_nuclide) % fission / & + micro_xs(p % event_nuclide) % absorption + end if + end associate end if else if (i_nuclide > 0) then - score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux + associate (nuc => nuclides(i_nuclide)) + if (nuc % fissionable) then + score = nuc%reactions(nuc%index_fission(1))%Q_value * & + micro_xs(i_nuclide)%fission * atom_density * flux + end if + end associate else - score = material_xs % kappa_fission * flux + do l = 1, materials(p%material)%n_nuclides + ! Determine atom density and index of nuclide + atom_density_ = materials(p%material)%atom_density(l) + i_nuc = materials(p%material)%nuclide(l) + + ! If nuclide is fissionable, accumulate kappa fission + associate(nuc => nuclides(i_nuc)) + if (nuc % fissionable) then + score = score + nuc%reactions(nuc%index_fission(1))%Q_value * & + micro_xs(i_nuc)%fission * atom_density_ * flux + end if + end associate + end do end if end if - case (SCORE_EVENTS) ! Simply count number of scoring events score = ONE diff --git a/tests/test_many_scores/results_true.dat b/tests/test_many_scores/results_true.dat index 9309d0964..0d5dd6e32 100644 --- a/tests/test_many_scores/results_true.dat +++ b/tests/test_many_scores/results_true.dat @@ -33,8 +33,8 @@ tally 1: 7.620560E-01 1.816851E+00 1.102658E+00 -1.338067E+02 -5.986137E+03 +1.337996E+02 +5.985519E+03 2.247257E+01 1.683779E+02 1.512960E-01 diff --git a/tests/test_score_kappafission/results_true.dat b/tests/test_score_kappafission/results_true.dat index 75d37cd87..976eefa36 100644 --- a/tests/test_score_kappafission/results_true.dat +++ b/tests/test_score_kappafission/results_true.dat @@ -1,17 +1,17 @@ k-combined: 9.903196E-01 4.279617E-02 tally 1: -2.266048E+02 -1.049833E+04 +2.266169E+02 +1.049923E+04 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.366139E+02 -3.859561E+03 +1.366590E+02 +3.861651E+03 tally 2: -2.402814E+02 -1.174003E+04 +2.403775E+02 +1.175130E+04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -19,11 +19,11 @@ tally 2: 1.270420E+02 3.297537E+03 tally 3: -2.217075E+02 -1.003168E+04 +2.217588E+02 +1.003581E+04 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.375693E+02 -3.872389E+03 +1.376303E+02 +3.875598E+03 From f17a1f436a9e320910034b2699dabcb0287245e3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 11 Nov 2015 11:42:59 -0600 Subject: [PATCH 11/19] Damn you gfortran 4.6. Comment out a perfectly-legitimate deallocate. --- src/ace_header.F90 | 7 ------- src/endf_header.F90 | 22 ---------------------- src/global.F90 | 7 ++++++- 3 files changed, 6 insertions(+), 30 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 6c27747d8..985371ff1 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -292,12 +292,6 @@ module ace_header class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear - ! Clear p_valid - call this % p_valid % clear() - - if (allocated(this % data)) & - deallocate(this % data) - if (associated(this % next)) then ! recursively clear this item call this % next % clear() @@ -346,7 +340,6 @@ module ace_header do i = 1, size(this % reactions) call this % reactions(i) % clear() end do - deallocate(this % reactions) end if call this % reaction_index % clear() diff --git a/src/endf_header.F90 b/src/endf_header.F90 index 54af0f738..af62231a5 100644 --- a/src/endf_header.F90 +++ b/src/endf_header.F90 @@ -13,28 +13,6 @@ module endf_header integer :: n_pairs ! # of pairs of (x,y) values real(8), allocatable :: x(:) ! values of abscissa real(8), allocatable :: y(:) ! values of ordinate - - ! Type-Bound procedures - contains - procedure :: clear => tab1_clear ! deallocates a Tab1 Object. end type Tab1 - contains - -!=============================================================================== -! TAB1_CLEAR deallocates the items in Tab1 -!=============================================================================== - - subroutine tab1_clear(this) - - class(Tab1), intent(inout) :: this ! The Tab1 to clear - - if (allocated(this % nbt)) & - deallocate(this % nbt, this % int) - - if (allocated(this % x)) & - deallocate(this % x, this % y) - - end subroutine tab1_clear - end module endf_header diff --git a/src/global.F90 b/src/global.F90 index 1d8099630..88ace73b6 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -436,7 +436,12 @@ contains do i = 1, size(nuclides) call nuclides(i) % clear() end do - deallocate(nuclides) + + ! WARNING: The following statement should work but doesn't under gfortran + ! 4.6 because of a bug. Technically, commenting this out leaves a memory + ! leak. + + ! deallocate(nuclides) end if if (allocated(nuclides_0K)) then From 7661c8902f1abc5f4f9a649f3f9808dfd0306ab7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Nov 2015 06:57:56 -0600 Subject: [PATCH 12/19] Purify yield_delayed and fix typo. --- src/ace.F90 | 2 +- src/fission.F90 | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 4cd53a251..c6d9b0221 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -913,7 +913,7 @@ contains cycle elseif (LOCB == 0) then ! No angular distribution data are given for this reaction, isotropic - ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0) + ! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0) cycle end if rxn % has_angle_dist = .true. diff --git a/src/fission.F90 b/src/fission.F90 index 3188138f2..4c7613db3 100644 --- a/src/fission.F90 +++ b/src/fission.F90 @@ -108,8 +108,7 @@ contains ! a given nuclide and incoming neutron energy in a given delayed group. !=============================================================================== - function yield_delayed(nuc, E, g) result(yield) - + pure function yield_delayed(nuc, E, g) result(yield) type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu real(8), intent(in) :: E ! energy of incoming neutron real(8) :: yield ! delayed neutron precursor yield From 0927f2530065f91ebf951d1f44f5596fc14460c6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Nov 2015 07:57:06 -0600 Subject: [PATCH 13/19] Fix bug in reimplementation of neighbor_lists. --- src/geometry.F90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 75165158b..3f8b43d4c 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -903,13 +903,17 @@ contains do i = 1, n_surfaces ! Copy positive neighbors to Surface instance j = neighbor_pos(i)%size() - allocate(surfaces(i)%obj%neighbor_pos(j)) - surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:j) + if (j > 0) then + allocate(surfaces(i)%obj%neighbor_pos(j)) + surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:j) + end if ! Copy negative neighbors to Surface instance j = neighbor_neg(i)%size() - allocate(surfaces(i)%obj%neighbor_neg(j)) - surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:j) + if (j > 0) then + allocate(surfaces(i)%obj%neighbor_neg(j)) + surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:j) + end if end do end subroutine neighbor_lists From 4d8015d11d5ad828f5e6ee1bab24bfe806b269af Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 19 Nov 2015 17:41:01 -0500 Subject: [PATCH 14/19] Implemented initial fixed source mode in Python API --- openmc/settings.py | 76 +++++++++++++++++++++++++--------------------- 1 file changed, 41 insertions(+), 35 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 519b5c7cf..8cb438a33 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -20,6 +20,8 @@ class SettingsFile(object): Attributes ---------- + run_mode : {'eigenvalue' or 'fixed source'} + The type of calculation to perform (default is 'eigenvalue') batches : int Number of batches to simulate generations_per_batch : int @@ -122,7 +124,9 @@ class SettingsFile(object): """ def __init__(self): - # Eigenvalue subelement + + # Run mode subelement (default is 'eigenvalue') + self._run_mode = 'eigenvalue' self._batches = None self._generations_per_batch = None self._inactive = None @@ -196,9 +200,13 @@ class SettingsFile(object): self._dd_count_interactions = False self._settings_file = ET.Element("settings") - self._eigenvalue_subelement = None + self._run_mode_subelement = None self._source_element = None + @property + def run_mode(self): + return self._run_mode + @property def batches(self): return self._batches @@ -399,6 +407,14 @@ class SettingsFile(object): def dd_count_interactions(self): return self._dd_count_interactions + @run_mode.setter + def run_mode(self, run_mode): + if not 'run_mode' in ['eigenvalue', 'fixed source']: + msg = 'Unable to set run mode to "{0}". Only "eigenvalue" ' \ + 'and "fixed source" are supported."'.format(run_mode) + raise ValueError(msg) + self._run_mode = run_mode + @batches.setter def batches(self, batches): check_type('batches', batches, Integral) @@ -861,57 +877,47 @@ class SettingsFile(object): self._dd_count_interactions = interactions - def _create_eigenvalue_subelement(self): - self._create_particles_subelement() - self._create_batches_subelement() - self._create_inactive_subelement() - self._create_generations_per_batch_subelement() - self._create_keff_trigger_subelement() + def _create_run_mode_subelement(self): + + if self.run_mode == 'eigenvalue': + self._run_mode_subelement = \ + ET.SubElement(self._settings_file, "eigenvalue") + self._create_batches_subelement() + self._create_generations_per_batch_subelement() + self._create_inactive_subelement() + self._create_particles_subelement() + self._create_keff_trigger_subelement() + else: + if self._run_mode_subelement is None: + self._run_mode_subelement = \ + ET.SubElement(self._settings_file, "fixed_source") + self._create_batches_subelement() + self._create_particles_subelement() def _create_batches_subelement(self): if self._batches is not None: - if self._eigenvalue_subelement is None: - self._eigenvalue_subelement = ET.SubElement(self._settings_file, - "eigenvalue") - - element = ET.SubElement(self._eigenvalue_subelement, "batches") + element = ET.SubElement(self._run_mode_subelement, "batches") element.text = str(self._batches) def _create_generations_per_batch_subelement(self): if self._generations_per_batch is not None: - if self._eigenvalue_subelement is None: - self._eigenvalue_subelement = ET.SubElement(self._settings_file, - "eigenvalue") - - element = ET.SubElement(self._eigenvalue_subelement, + element = ET.SubElement(self._run_mode_subelement, "generations_per_batch") element.text = str(self._generations_per_batch) def _create_inactive_subelement(self): if self._inactive is not None: - if self._eigenvalue_subelement is None: - self._eigenvalue_subelement = ET.SubElement(self._settings_file, - "eigenvalue") - - element = ET.SubElement(self._eigenvalue_subelement, "inactive") + element = ET.SubElement(self._run_mode_subelement, "inactive") element.text = str(self._inactive) def _create_particles_subelement(self): if self._particles is not None: - if self._eigenvalue_subelement is None: - self._eigenvalue_subelement = ET.SubElement(self._settings_file, - "eigenvalue") - - element = ET.SubElement(self._eigenvalue_subelement, "particles") + element = ET.SubElement(self._run_mode_subelement, "particles") element.text = str(self._particles) def _create_keff_trigger_subelement(self): if self._keff_trigger is not None: - if self._eigenvalue_subelement is None: - self._eigenvalue_subelement = ET.SubElement(self._settings_file, - "eigenvalue") - - element = ET.SubElement(self._eigenvalue_subelement, "keff_trigger") + element = ET.SubElement(self._run_mode_subelement, "keff_trigger") for key in self._keff_trigger: subelement = ET.SubElement(element, key) @@ -1182,10 +1188,10 @@ class SettingsFile(object): self._settings_file.clear() self._source_subelement = None self._trigger_subelement = None - self._eigenvalue_subelement = None + self._run_mode_subelement = None self._source_element = None - self._create_eigenvalue_subelement() + self._create_run_mode_subelement() self._create_source_subelement() self._create_output_subelement() self._create_statepoint_subelement() From fca18c7de454dc66deb33450ccd46c4eb0afbb38 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 19 Nov 2015 17:47:29 -0500 Subject: [PATCH 15/19] Reverted to original ordering for Python API Settings XML attribute creation --- openmc/settings.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 8cb438a33..742442ba9 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -882,17 +882,17 @@ class SettingsFile(object): if self.run_mode == 'eigenvalue': self._run_mode_subelement = \ ET.SubElement(self._settings_file, "eigenvalue") - self._create_batches_subelement() - self._create_generations_per_batch_subelement() - self._create_inactive_subelement() self._create_particles_subelement() + self._create_batches_subelement() + self._create_inactive_subelement() + self._create_generations_per_batch_subelement() self._create_keff_trigger_subelement() else: if self._run_mode_subelement is None: self._run_mode_subelement = \ ET.SubElement(self._settings_file, "fixed_source") - self._create_batches_subelement() self._create_particles_subelement() + self._create_batches_subelement() def _create_batches_subelement(self): if self._batches is not None: From f61408fb21c5c5400f2fe08d29667f94e5249f48 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 19 Nov 2015 20:58:30 -0500 Subject: [PATCH 16/19] Changed Python API SettingsFile not run_mode in to more Pythonic run_mode not in --- openmc/settings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/settings.py b/openmc/settings.py index 742442ba9..9eb54b9eb 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -409,7 +409,7 @@ class SettingsFile(object): @run_mode.setter def run_mode(self, run_mode): - if not 'run_mode' in ['eigenvalue', 'fixed source']: + if 'run_mode' not in ['eigenvalue', 'fixed source']: msg = 'Unable to set run mode to "{0}". Only "eigenvalue" ' \ 'and "fixed source" are supported."'.format(run_mode) raise ValueError(msg) From 088e8652440b78add7b23ef8c7f3cbfa978f96d9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 20 Nov 2015 19:45:55 -0600 Subject: [PATCH 17/19] Respond to comments by @wbinventor and @nelsonag on #500 --- src/ace.F90 | 18 +++++++++--------- src/cross_section.F90 | 10 +++++----- src/fission.F90 | 10 +++++----- src/geometry.F90 | 17 +++++++++-------- src/output.F90 | 9 ++++----- src/tally.F90 | 24 ++++++++++++------------ 6 files changed, 44 insertions(+), 44 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index c6d9b0221..e97338b55 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -733,13 +733,13 @@ contains ! sigma array is not allocated or stored for elastic scattering since it is ! already stored in nuc % elastic associate (rxn => nuc % reactions(1)) - rxn % MT = 2 - rxn % Q_value = ZERO - rxn % multiplicity = 1 - rxn % threshold = 1 - rxn % scatter_in_cm = .true. - rxn % has_angle_dist = .false. - rxn % has_energy_dist = .false. + rxn%MT = 2 + rxn%Q_value = ZERO + rxn%multiplicity = 1 + rxn%threshold = 1 + rxn%scatter_in_cm = .true. + rxn%has_angle_dist = .false. + rxn%has_energy_dist = .false. end associate ! Add contribution of elastic scattering to total cross section @@ -1013,8 +1013,8 @@ contains recursive subroutine get_energy_dist(edist, loc_law, delayed_n) type(DistEnergy), intent(inout) :: edist ! energy distribution - integer, intent(in) :: loc_law ! locator for data - logical, intent(in), optional :: delayed_n ! is this for delayed neutrons? + integer, intent(in) :: loc_law ! locator for data + logical, intent(in), optional :: delayed_n ! is this for delayed neutrons? integer :: LDIS ! location of all energy distributions integer :: LNW ! location of next energy distribution if multiple diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 969a397f0..f874eb2a7 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -510,8 +510,8 @@ contains pure function find_energy_index(mat, E) result(i) type(Material), intent(in) :: mat ! pointer to current material - real(8), intent(in) :: E ! energy of particle - integer :: i ! energy grid index + real(8), intent(in) :: E ! energy of particle + integer :: i ! energy grid index ! if the energy is outside of energy grid range, set to first or last ! index. Otherwise, do a binary search through the union energy grid. @@ -531,9 +531,9 @@ contains !=============================================================================== pure function elastic_xs_0K(E, nuc) result(xs_out) - real(8), intent(in) :: E ! trial energy - type(Nuclide), intent(in) :: nuc ! target nuclide at temperature - real(8) :: xs_out ! 0K xs at trial energy + real(8), intent(in) :: E ! trial energy + type(Nuclide), intent(in) :: nuc ! target nuclide at temperature + real(8) :: xs_out ! 0K xs at trial energy integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid diff --git a/src/fission.F90 b/src/fission.F90 index 4c7613db3..a005d9f5b 100644 --- a/src/fission.F90 +++ b/src/fission.F90 @@ -17,8 +17,8 @@ contains pure function nu_total(nuc, E) result(nu) type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu - real(8), intent(in) :: E ! energy of incoming neutron - real(8) :: nu ! number of total neutrons emitted per fission + real(8), intent(in) :: E ! energy of incoming neutron + real(8) :: nu ! number of total neutrons emitted per fission integer :: i ! loop index integer :: NC ! number of polynomial coefficients @@ -50,8 +50,8 @@ contains pure function nu_prompt(nuc, E) result(nu) type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu - real(8), intent(in) :: E ! energy of incoming neutron - real(8) :: nu ! number of prompt neutrons emitted per fission + real(8), intent(in) :: E ! energy of incoming neutron + real(8) :: nu ! number of prompt neutrons emitted per fission integer :: i ! loop index integer :: NC ! number of polynomial coefficients @@ -87,7 +87,7 @@ contains pure function nu_delayed(nuc, E) result(nu) type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu - real(8), intent(in) :: E ! energy of incoming neutron + real(8), intent(in) :: E ! energy of incoming neutron real(8) :: nu ! number of delayed neutrons emitted per fission if (nuc % nu_d_type == NU_NONE) then diff --git a/src/geometry.F90 b/src/geometry.F90 index 3f8b43d4c..22f5c3a15 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -875,6 +875,7 @@ contains integer :: i ! index in cells/surfaces array integer :: j ! index in region specification integer :: k ! surface half-space spec + integer :: n ! size of vector type(VectorInt), allocatable :: neighbor_pos(:) type(VectorInt), allocatable :: neighbor_neg(:) @@ -902,17 +903,17 @@ contains do i = 1, n_surfaces ! Copy positive neighbors to Surface instance - j = neighbor_pos(i)%size() - if (j > 0) then - allocate(surfaces(i)%obj%neighbor_pos(j)) - surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:j) + n = neighbor_pos(i)%size() + if (n > 0) then + allocate(surfaces(i)%obj%neighbor_pos(n)) + surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:n) end if ! Copy negative neighbors to Surface instance - j = neighbor_neg(i)%size() - if (j > 0) then - allocate(surfaces(i)%obj%neighbor_neg(j)) - surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:j) + n = neighbor_neg(i)%size() + if (n > 0) then + allocate(surfaces(i)%obj%neighbor_neg(n)) + surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:n) end if end do diff --git a/src/output.F90 b/src/output.F90 index 19ae6268b..c31b20b70 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -194,8 +194,8 @@ contains !=============================================================================== subroutine write_message(message, level) - character(*), intent(in) :: message - integer, intent(in), optional :: level ! verbosity level + character(*), intent(in) :: message ! message to write + integer, intent(in), optional :: level ! verbosity level integer :: i_start ! starting position integer :: i_end ! ending position @@ -1349,10 +1349,9 @@ contains !=============================================================================== function get_label(t, i_filter) result(label) - type(TallyObject), intent(in) :: t ! tally object - integer, intent(in) :: i_filter ! index in filters array - character(100) :: label ! user-specified identifier + integer, intent(in) :: i_filter ! index in filters array + character(100) :: label ! user-specified identifier integer :: i ! index in cells/surfaces/etc array integer :: bin diff --git a/src/tally.F90 b/src/tally.F90 index c4ca2028b..09e30a242 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -608,13 +608,13 @@ contains ! No fission events occur if survival biasing is on -- need to ! calculate fraction of absorptions that would have resulted in ! fission scale by kappa-fission - associate (nuc => nuclides(p % event_nuclide)) - if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & - nuc % fissionable) then - score = p % absorb_wgt * & + associate (nuc => nuclides(p%event_nuclide)) + if (micro_xs(p%event_nuclide)%absorption > ZERO .and. & + nuc%fissionable) then + score = p%absorb_wgt * & nuc%reactions(nuc%index_fission(1))%Q_value * & - micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption + micro_xs(p%event_nuclide)%fission / & + micro_xs(p%event_nuclide)%absorption end if end associate else @@ -623,12 +623,12 @@ contains ! All fission events will contribute, so again we can use ! particle's weight entering the collision as the estimate for ! the fission energy production rate - associate (nuc => nuclides(p % event_nuclide)) - if (nuc % fissionable) then - score = p % last_wgt * & + associate (nuc => nuclides(p%event_nuclide)) + if (nuc%fissionable) then + score = p%last_wgt * & nuc%reactions(nuc%index_fission(1))%Q_value * & - micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption + micro_xs(p%event_nuclide)%fission / & + micro_xs(p%event_nuclide)%absorption end if end associate end if @@ -636,7 +636,7 @@ contains else if (i_nuclide > 0) then associate (nuc => nuclides(i_nuclide)) - if (nuc % fissionable) then + if (nuc%fissionable) then score = nuc%reactions(nuc%index_fission(1))%Q_value * & micro_xs(i_nuclide)%fission * atom_density * flux end if From db550a05370bc81189134f9cc2fb1cab1596cc3d Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Tue, 24 Nov 2015 21:08:54 -0500 Subject: [PATCH 18/19] Fixed bug casting OpenCG rotations to integers is now double --- openmc/opencg_compatible.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index 93c0e5fae..6430b2424 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -726,7 +726,7 @@ def get_openmc_cell(opencg_cell): openmc_cell.fill = get_openmc_material(fill) if opencg_cell.rotation: - rotation = np.asarray(opencg_cell.rotation, dtype=np.int) + rotation = np.asarray(opencg_cell.rotation, dtype=np.float64) openmc_cell.rotation = rotation if opencg_cell.translation: From a11950f94ae890f8c78b2be736ff006d71e19fc4 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Tue, 24 Nov 2015 22:01:12 -0500 Subject: [PATCH 19/19] Now over-riding openmc/opencg geometries in MGXS Library when loading from StatePoint --- openmc/mgxs/library.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 87c6665b2..fa49d24e6 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -361,6 +361,8 @@ class Library(object): raise ValueError(msg) self._sp_filename = statepoint._f.filename + self._openmc_geometry = statepoint.summary.openmc_geometry + self._opencg_geometry = None # Load tallies for each MGXS for each domain and mgxs type for domain in self.domains: