From 1c5016167a5712b8dc3c6b0acf246839e12ef116 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 21 Jun 2017 11:24:37 -0400 Subject: [PATCH 01/12] Allow partly S(a,b), partly amorphus scattering --- openmc/material.py | 16 ++++++++-- src/cross_section.F90 | 71 ++++++++++++++++++++++++----------------- src/input_xml.F90 | 18 +++++++++-- src/material_header.F90 | 3 +- src/nuclide_header.F90 | 3 +- src/physics.F90 | 60 +++++++++++++++++++--------------- src/tally.F90 | 3 +- 7 files changed, 111 insertions(+), 63 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 762332584..2c633a5ac 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -618,13 +618,18 @@ class Material(IDManagerMixin): if element == elm[0]: self._elements.remove(elm) - def add_s_alpha_beta(self, name): + def add_s_alpha_beta(self, name, fraction=1.0): r"""Add an :math:`S(\alpha,\beta)` table to the material Parameters ---------- name : str Name of the :math:`S(\alpha,\beta)` table + fraction : float + The fraction of relevant nuclei that are affected by the + :math:`S(\alpha,\beta)` table. For example, if the material is a + block of carbon that is 60% graphite and 40% amorphus then add a + graphite :math:`S(\alpha,\beta)` table with fraction=0.6. """ @@ -638,13 +643,17 @@ class Material(IDManagerMixin): 'non-string table name "{}"'.format(self._id, name) raise ValueError(msg) + cv.check_type('S(a,b) fraction', fraction, Real) + cv.check_greater_than('S(a,b) fraction', fraction, 0.0, True) + cv.check_less_than('S(a,b) fraction', fraction, 1.0, True) + new_name = openmc.data.get_thermal_name(name) if new_name != name: msg = 'OpenMC S(a,b) tables follow the GND naming convention. ' \ 'Table "{}" is being renamed as "{}".'.format(name, new_name) warnings.warn(msg) - self._sab.append(new_name) + self._sab.append((new_name, fraction)) def make_isotropic_in_lab(self): for nuclide, percent, percent_type in self._nuclides: @@ -972,7 +981,8 @@ class Material(IDManagerMixin): if len(self._sab) > 0: for sab in self._sab: subelement = ET.SubElement(element, "sab") - subelement.set("name", sab) + subelement.set("name", sab[0]) + subelement.set("fraction", sab[1]) return element diff --git a/src/cross_section.F90 b/src/cross_section.F90 index ce7d97133..1306fb669 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -36,6 +36,7 @@ contains integer :: i_grid ! index into logarithmic mapping array or material ! union grid real(8) :: atom_density ! atom density of a nuclide + real(8) :: sab_frac ! fraction of atoms affected by S(a,b) logical :: check_sab ! should we check for S(a,b) table? ! Set all material macroscopic cross sections to zero @@ -71,6 +72,7 @@ contains if (i == mat % i_sab_nuclides(j)) then ! Get index in sab_tables i_sab = mat % i_sab_tables(j) + sab_frac = mat % sab_fracs(j) ! If particle energy is greater than the highest energy for the S(a,b) ! table, don't use the S(a,b) table @@ -84,7 +86,7 @@ contains end if end if - ! ======================================================================== + ! ====================================================================== ! CALCULATE MICROSCOPIC CROSS SECTION ! Determine microscopic cross sections for this nuclide @@ -93,12 +95,14 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E & .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, & + p % sqrtkT, sab_frac) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, & + p % sqrtkT, sab_frac) end if - ! ======================================================================== + ! ====================================================================== ! ADD TO MACROSCOPIC CROSS SECTION ! Copy atom density of nuclide in material @@ -110,7 +114,8 @@ contains ! Add contributions to material macroscopic scattering cross section material_xs % elastic = material_xs % elastic + & - atom_density * micro_xs(i_nuclide) % elastic + atom_density * (micro_xs(i_nuclide) % elastic & + + micro_xs(i_nuclide) % scatter_sab) ! Add contributions to material macroscopic absorption cross section material_xs % absorption = material_xs % absorption + & @@ -133,13 +138,15 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT) - 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 + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT, & + sab_frac) + 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_log_union ! index into logarithmic mapping array or ! material union energy grid - real(8), intent(in) :: sqrtkT ! Square root of kT, material dependent + real(8), intent(in) :: sqrtkT ! square root of kT, material dependent + real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) logical :: use_mp ! true if XS can be calculated with windowed multipole integer :: i_temp ! index for temperature @@ -243,8 +250,10 @@ contains micro_xs(i_nuclide) % interp_factor = f ! Initialize nuclide cross-sections to zero - micro_xs(i_nuclide) % fission = ZERO - micro_xs(i_nuclide) % nu_fission = ZERO + micro_xs(i_nuclide) % fission = ZERO + micro_xs(i_nuclide) % nu_fission = ZERO + micro_xs(i_nuclide) % scatter_sab = ZERO + micro_xs(i_nuclide) % elastic_sab = ZERO ! Calculate microscopic nuclide total cross section micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) & @@ -277,14 +286,15 @@ contains ! Initialize URR probability table treatment to false micro_xs(i_nuclide) % use_ptable = .false. - ! If there is S(a,b) data for this nuclide, we need to do a few - ! things. Since the total cross section was based on non-S(a,b) data, we - ! need to correct it by subtracting the non-S(a,b) elastic cross section and - ! then add back in the calculated S(a,b) elastic+inelastic cross section. + ! If there is S(a,b) data for this nuclide, we need to set the sab_scatter + ! and sab_elastic cross sections and correct the total and elastic cross + ! sections. - if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT) + if (i_sab > 0) then + call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT, sab_frac) + end if - ! if the particle is in the unresolved resonance range and there are + ! If the particle is in the unresolved resonance range and there are ! probability tables, we need to determine cross sections from the table if (urr_ptables_on .and. nuc % urr_present .and. .not. use_mp) then @@ -303,16 +313,16 @@ contains !=============================================================================== ! CALCULATE_SAB_XS determines the elastic and inelastic scattering -! cross-sections in the thermal energy range. These cross sections replace -! whatever data were taken from the normal Nuclide table. +! cross-sections in the thermal energy range. These cross sections replace a +! fraction of whatever data were taken from the normal Nuclide table. !=============================================================================== - subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT) - + subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT, sab_frac) 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 real(8), intent(in) :: sqrtkT ! temperature + real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) integer :: i_grid ! index on S(a,b) energy grid integer :: i_temp ! temperature index @@ -341,7 +351,8 @@ contains ! Randomly sample between temperature i and i+1 f = (kT - sab_tables(i_sab) % kTs(i_temp)) / & - (sab_tables(i_sab) % kTs(i_temp + 1) - sab_tables(i_sab) % kTs(i_temp)) + (sab_tables(i_sab) % kTs(i_temp + 1) & + - sab_tables(i_sab) % kTs(i_temp)) if (f > prn()) i_temp = i_temp + 1 end if @@ -402,13 +413,15 @@ contains end if end associate - ! Correct total and elastic cross sections - micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % elastic + inelastic + elastic - micro_xs(i_nuclide) % elastic = inelastic + elastic + ! Store the S(a,b) cross sections. + micro_xs(i_nuclide) % scatter_sab = sab_frac * (elastic + inelastic) + micro_xs(i_nuclide) % elastic_sab = sab_frac * elastic - ! Store S(a,b) elastic cross section for sampling later - micro_xs(i_nuclide) % elastic_sab = elastic + ! Correct total and elastic cross sections + micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total & + + sab_frac * (elastic + inelastic - micro_xs(i_nuclide) % elastic) + micro_xs(i_nuclide) % elastic = & + (ONE - sab_frac) * micro_xs(i_nuclide) % elastic ! Save temperature index micro_xs(i_nuclide) % index_temp_sab = i_temp diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 5ff4d119e..2e7219a31 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2565,6 +2565,7 @@ contains ! table is indeed applied to multiple nuclides. allocate(mat % sab_names(n_sab)) allocate(mat % i_sab_tables(n_sab)) + allocate(mat % sab_fracs(n_sab)) do j = 1, n_sab ! Get pointer to S(a,b) table @@ -2578,6 +2579,13 @@ contains name = trim(name) mat % sab_names(j) = name + ! Read the fraction of nuclei affected by this S(a,b) table + if (check_for_node(node_sab, "fraction")) then + call get_node_value(node_sab, "fraction", mat % sab_fracs(j)) + else + mat % sab_fracs(j) = ONE + end if + ! Check that this nuclide is listed in the cross_sections.xml file if (.not. library_dict % has_key(to_lower(name))) then call fatal_error("Could not find S(a,b) table " // trim(name) & @@ -5151,8 +5159,9 @@ contains integer :: temp_nuclide ! temporary value for sorting integer :: temp_table ! temporary value for sorting logical :: found - type(VectorInt) :: i_sab_tables - type(VectorInt) :: i_sab_nuclides + type(VectorInt) :: i_sab_tables + type(VectorInt) :: i_sab_nuclides + type(VectorReal) :: sab_fracs do i = 1, size(materials) ! Skip materials with no S(a,b) tables @@ -5170,6 +5179,7 @@ contains if (any(sab % nuclides == nuclides(mat % nuclide(j)) % name)) then call i_sab_tables % push_back(mat % i_sab_tables(k)) call i_sab_nuclides % push_back(j) + call sab_fracs % push_back(mat % sab_fracs(k)) found = .true. end if end do FIND_NUCLIDE @@ -5185,15 +5195,19 @@ contains ! Update i_sab_tables and i_sab_nuclides deallocate(mat % i_sab_tables) + deallocate(mat % sab_fracs) m = i_sab_tables % size() allocate(mat % i_sab_tables(m)) allocate(mat % i_sab_nuclides(m)) + allocate(mat % sab_fracs(m)) mat % i_sab_tables(:) = i_sab_tables % data(1:m) mat % i_sab_nuclides(:) = i_sab_nuclides % data(1:m) + mat % sab_fracs = sab_fracs % data(1:m) ! Clear entries in vectors for next material call i_sab_tables % clear() call i_sab_nuclides % clear() + call sab_fracs % clear() ! If there are multiple S(a,b) tables, we need to make sure that the ! entries in i_sab_nuclides are sorted or else they won't be applied diff --git a/src/material_header.F90 b/src/material_header.F90 index fdc36548f..17fe8d248 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -22,10 +22,11 @@ module material_header ! Unionized energy grid information integer, allocatable :: nuclide_grid_index(:,:) ! nuclide e_grid pointers - ! S(a,b) data references + ! S(a,b) data integer :: n_sab = 0 ! number of S(a,b) tables integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide integer, allocatable :: i_sab_tables(:) ! index in sab_tables + real(8), allocatable :: sab_fracs(:) ! how often to use S(a,b) ! Temporary names read during initialization character(20), allocatable :: names(:) ! isotope names diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index ea6eab328..41b1008a8 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -112,10 +112,11 @@ module nuclide_header real(8) :: last_E = ZERO ! last evaluated energy real(8) :: interp_factor ! interpolation factor on nuc. energy grid real(8) :: total ! microscropic total xs - real(8) :: elastic ! microscopic elastic scattering xs + real(8) :: elastic ! microscopic elastic scattering xs (non S(a,b)) real(8) :: absorption ! microscopic absorption xs real(8) :: fission ! microscopic fission xs real(8) :: nu_fission ! microscopic production xs + real(8) :: scatter_sab ! microscopic scattering xs due to S(a,b) ! Information for S(a,b) use integer :: index_sab ! index in sab_tables (zero means no table) diff --git a/src/physics.F90 b/src/physics.F90 index 9482772ef..15b53fd7f 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -324,6 +324,7 @@ contains real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering real(8) :: phi ! azimuthal angle for iso-in-lab scattering real(8) :: kT ! temperature in eV + logical :: sampled ! whether or not a reaction type has been sampled type(Nuclide), pointer :: nuc ! copy incoming direction @@ -337,36 +338,43 @@ contains ! For tallying purposes, this routine might be called directly. In that ! case, we need to sample a reaction via the cutoff variable - prob = ZERO cutoff = prn() * (micro_xs(i_nuclide) % total - & micro_xs(i_nuclide) % absorption) + sampled = .false. - prob = prob + micro_xs(i_nuclide) % elastic + prob = micro_xs(i_nuclide) % elastic if (prob > cutoff) then ! ======================================================================= - ! ELASTIC SCATTERING - - if (micro_xs(i_nuclide) % index_sab /= NONE) then - ! S(a,b) scattering - call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, & - p % E, p % coord(1) % uvw, p % mu) + ! NON-S(A,B) ELASTIC SCATTERING + ! Determine temperature + if (nuc % mp_present) then + kT = p % sqrtkT**2 else - ! Determine temperature - if (nuc % mp_present) then - kT = p % sqrtkT**2 - else - kT = nuc % kTs(micro_xs(i_nuclide) % index_temp) - end if - - ! Perform collision physics for elastic scattering - call elastic_scatter(i_nuclide, nuc % reactions(1), kT, & - p % E, p % coord(1) % uvw, p % mu, p % wgt) + kT = nuc % kTs(micro_xs(i_nuclide) % index_temp) end if - p % event_MT = ELASTIC + ! Perform collision physics for elastic scattering + call elastic_scatter(i_nuclide, nuc % reactions(1), kT, p % E, & + p % coord(1) % uvw, p % mu, p % wgt) - else + p % event_MT = ELASTIC + sampled = .true. + end if + + prob = prob + micro_xs(i_nuclide) % scatter_sab + if (prob > cutoff .and. .not. sampled) then + ! ======================================================================= + ! S(A,B) SCATTERING + + call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, p % E, & + p % coord(1) % uvw, p % mu) + + p % event_MT = ELASTIC + sampled = .true. + end if + + if (.not. sampled) then ! ======================================================================= ! INELASTIC SCATTERING @@ -388,7 +396,7 @@ contains if (rx % MT == N_FISSION .or. rx % MT == N_F .or. rx % MT == N_NF & .or. rx % MT == N_2NF .or. rx % MT == N_3NF) cycle - ! some materials have gas production cross sections with MT > 200 that + ! 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 (rx % MT >= 200 .or. rx % MT == N_LEVEL) cycle @@ -406,24 +414,24 @@ contains ! Perform collision physics for inelastic scattering call inelastic_scatter(nuc, nuc%reactions(i), p) - p % event_MT = nuc%reactions(i)%MT + p % event_MT = nuc % reactions(i) % MT end if ! Set event component p % event = EVENT_SCATTER - ! sample new outgoing angle for isotropic in lab scattering + ! Sample new outgoing angle for isotropic in lab scattering if (materials(p % material) % p0(i_nuc_mat)) then - ! sample isotropic-in-lab outgoing direction + ! Sample isotropic-in-lab outgoing direction uvw_new(1) = TWO * prn() - ONE phi = TWO * PI * prn() uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) p % mu = dot_product(uvw_old, uvw_new) - ! change direction of particle + ! Change direction of particle p % coord(1) % uvw = uvw_new end if @@ -560,7 +568,7 @@ contains ! Determine whether inelastic or elastic scattering will occur if (prn() < micro_xs(i_nuclide) % elastic_sab / & - micro_xs(i_nuclide) % elastic) then + micro_xs(i_nuclide) % scatter_sab) then ! elastic scattering ! Get index and interpolation factor for elastic grid diff --git a/src/tally.F90 b/src/tally.F90 index 0ee041cbe..d827d9bac 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1051,7 +1051,8 @@ contains else if (i_nuclide > 0) then - score = micro_xs(i_nuclide) % elastic * atom_density * flux + score = (micro_xs(i_nuclide) % elastic & + + micro_xs(i_nuclide) % scatter_sab) * atom_density * flux else score = material_xs % elastic * flux end if From 182e8107d32fb5f044b00ba279db025bad7a8e1c Mon Sep 17 00:00:00 2001 From: tjlaboss Date: Tue, 27 Jun 2017 08:54:43 -0600 Subject: [PATCH 02/12] Float -> string for ElementTree --- openmc/material.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index 2c633a5ac..6ac153c86 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -982,7 +982,7 @@ class Material(IDManagerMixin): for sab in self._sab: subelement = ET.SubElement(element, "sab") subelement.set("name", sab[0]) - subelement.set("fraction", sab[1]) + subelement.set("fraction", str(sab[1])) return element From 9a7f9e0dc05de969b13afe833c6f85f418ca0ee7 Mon Sep 17 00:00:00 2001 From: tjlaboss Date: Thu, 29 Jun 2017 09:02:02 -0600 Subject: [PATCH 03/12] Updated test suite to use S(a,b) fraction inputs_true.dat for many of the tests now contain the new syntax --- tests/test_asymmetric_lattice/inputs_true.dat | 18 +++++++++--------- tests/test_diff_tally/inputs_true.dat | 18 +++++++++--------- tests/test_filter_energyfun/inputs_true.dat | 18 +++++++++--------- tests/test_filter_mesh/inputs_true.dat | 18 +++++++++--------- tests/test_iso_in_lab/inputs_true.dat | 18 +++++++++--------- .../test_mgxs_library_ce_to_mg/inputs_true.dat | 2 +- .../test_mgxs_library_condense/inputs_true.dat | 2 +- .../inputs_true.dat | 2 +- tests/test_mgxs_library_hdf5/inputs_true.dat | 2 +- tests/test_mgxs_library_mesh/inputs_true.dat | 18 +++++++++--------- .../inputs_true.dat | 2 +- .../test_mgxs_library_nuclides/inputs_true.dat | 2 +- tests/test_multipole/inputs_true.dat | 2 +- tests/test_periodic/inputs_true.dat | 2 +- tests/test_tallies/inputs_true.dat | 18 +++++++++--------- tests/test_tally_aggregation/inputs_true.dat | 18 +++++++++--------- tests/test_tally_arithmetic/inputs_true.dat | 18 +++++++++--------- tests/test_tally_slice_merge/inputs_true.dat | 18 +++++++++--------- tests/test_triso/inputs_true.dat | 8 ++++---- tests/test_volume_calc/inputs_true.dat | 2 +- 20 files changed, 103 insertions(+), 103 deletions(-) diff --git a/tests/test_asymmetric_lattice/inputs_true.dat b/tests/test_asymmetric_lattice/inputs_true.dat index 672f6bee1..50803105d 100644 --- a/tests/test_asymmetric_lattice/inputs_true.dat +++ b/tests/test_asymmetric_lattice/inputs_true.dat @@ -78,7 +78,7 @@ - + @@ -86,7 +86,7 @@ - + @@ -114,7 +114,7 @@ - + @@ -129,7 +129,7 @@ - + @@ -144,7 +144,7 @@ - + @@ -159,7 +159,7 @@ - + @@ -174,7 +174,7 @@ - + @@ -187,7 +187,7 @@ - + @@ -200,7 +200,7 @@ - + diff --git a/tests/test_diff_tally/inputs_true.dat b/tests/test_diff_tally/inputs_true.dat index ff134bd5c..a2c5638cf 100644 --- a/tests/test_diff_tally/inputs_true.dat +++ b/tests/test_diff_tally/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_filter_energyfun/inputs_true.dat b/tests/test_filter_energyfun/inputs_true.dat index d857451cb..d6f7c1f79 100644 --- a/tests/test_filter_energyfun/inputs_true.dat +++ b/tests/test_filter_energyfun/inputs_true.dat @@ -171,7 +171,7 @@ - + @@ -179,7 +179,7 @@ - + @@ -207,7 +207,7 @@ - + @@ -222,7 +222,7 @@ - + @@ -237,7 +237,7 @@ - + @@ -252,7 +252,7 @@ - + @@ -267,7 +267,7 @@ - + @@ -280,7 +280,7 @@ - + @@ -293,7 +293,7 @@ - + diff --git a/tests/test_filter_mesh/inputs_true.dat b/tests/test_filter_mesh/inputs_true.dat index 6d14f9e7e..153f4dcd9 100644 --- a/tests/test_filter_mesh/inputs_true.dat +++ b/tests/test_filter_mesh/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_iso_in_lab/inputs_true.dat b/tests/test_iso_in_lab/inputs_true.dat index 098056f5b..2311d8588 100644 --- a/tests/test_iso_in_lab/inputs_true.dat +++ b/tests/test_iso_in_lab/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_mgxs_library_ce_to_mg/inputs_true.dat b/tests/test_mgxs_library_ce_to_mg/inputs_true.dat index 8e8cde281..675ef50f1 100644 --- a/tests/test_mgxs_library_ce_to_mg/inputs_true.dat +++ b/tests/test_mgxs_library_ce_to_mg/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_mgxs_library_condense/inputs_true.dat b/tests/test_mgxs_library_condense/inputs_true.dat index d2f28d0fe..38a40c1cc 100644 --- a/tests/test_mgxs_library_condense/inputs_true.dat +++ b/tests/test_mgxs_library_condense/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_mgxs_library_distribcell/inputs_true.dat b/tests/test_mgxs_library_distribcell/inputs_true.dat index d7a4a186a..8b6522bdd 100644 --- a/tests/test_mgxs_library_distribcell/inputs_true.dat +++ b/tests/test_mgxs_library_distribcell/inputs_true.dat @@ -60,7 +60,7 @@ - + diff --git a/tests/test_mgxs_library_hdf5/inputs_true.dat b/tests/test_mgxs_library_hdf5/inputs_true.dat index d2f28d0fe..38a40c1cc 100644 --- a/tests/test_mgxs_library_hdf5/inputs_true.dat +++ b/tests/test_mgxs_library_hdf5/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_mgxs_library_mesh/inputs_true.dat b/tests/test_mgxs_library_mesh/inputs_true.dat index 4756b27bf..8bfcc64b6 100644 --- a/tests/test_mgxs_library_mesh/inputs_true.dat +++ b/tests/test_mgxs_library_mesh/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_mgxs_library_no_nuclides/inputs_true.dat b/tests/test_mgxs_library_no_nuclides/inputs_true.dat index d2f28d0fe..38a40c1cc 100644 --- a/tests/test_mgxs_library_no_nuclides/inputs_true.dat +++ b/tests/test_mgxs_library_no_nuclides/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_mgxs_library_nuclides/inputs_true.dat b/tests/test_mgxs_library_nuclides/inputs_true.dat index 826eb5b62..03e3b3ff6 100644 --- a/tests/test_mgxs_library_nuclides/inputs_true.dat +++ b/tests/test_mgxs_library_nuclides/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 0d1fe99bd..4b349ecba 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -25,7 +25,7 @@ - + diff --git a/tests/test_periodic/inputs_true.dat b/tests/test_periodic/inputs_true.dat index 61958701b..26e700004 100644 --- a/tests/test_periodic/inputs_true.dat +++ b/tests/test_periodic/inputs_true.dat @@ -16,7 +16,7 @@ - + diff --git a/tests/test_tallies/inputs_true.dat b/tests/test_tallies/inputs_true.dat index a85491e94..664e6e983 100644 --- a/tests/test_tallies/inputs_true.dat +++ b/tests/test_tallies/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_tally_aggregation/inputs_true.dat b/tests/test_tally_aggregation/inputs_true.dat index 7a8bf4613..6ef6bdb4d 100644 --- a/tests/test_tally_aggregation/inputs_true.dat +++ b/tests/test_tally_aggregation/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_tally_arithmetic/inputs_true.dat b/tests/test_tally_arithmetic/inputs_true.dat index 743ca3c58..e9c156287 100644 --- a/tests/test_tally_arithmetic/inputs_true.dat +++ b/tests/test_tally_arithmetic/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_tally_slice_merge/inputs_true.dat b/tests/test_tally_slice_merge/inputs_true.dat index b1c089c96..2c90e475a 100644 --- a/tests/test_tally_slice_merge/inputs_true.dat +++ b/tests/test_tally_slice_merge/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat index fdbc1cb5f..862c2db9a 100644 --- a/tests/test_triso/inputs_true.dat +++ b/tests/test_triso/inputs_true.dat @@ -403,12 +403,12 @@ - + - + @@ -420,12 +420,12 @@ - + - + diff --git a/tests/test_volume_calc/inputs_true.dat b/tests/test_volume_calc/inputs_true.dat index 28f1cbd9f..35c866028 100644 --- a/tests/test_volume_calc/inputs_true.dat +++ b/tests/test_volume_calc/inputs_true.dat @@ -16,7 +16,7 @@ - + From b3ee0652068025a2902447c480beb7a29c2d976a Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 29 Jun 2017 13:57:25 -0400 Subject: [PATCH 04/12] Add fraction to S(a,b) test; move to PyAPI --- tests/test_salphabeta/geometry.xml | 15 ----- tests/test_salphabeta/inputs_true.dat | 60 ++++++++++++++++++ tests/test_salphabeta/materials.xml | 41 ------------- tests/test_salphabeta/results_true.dat | 2 +- tests/test_salphabeta/settings.xml | 15 ----- tests/test_salphabeta/test_salphabeta.py | 77 +++++++++++++++++++++++- 6 files changed, 136 insertions(+), 74 deletions(-) delete mode 100644 tests/test_salphabeta/geometry.xml create mode 100644 tests/test_salphabeta/inputs_true.dat delete mode 100644 tests/test_salphabeta/materials.xml delete mode 100644 tests/test_salphabeta/settings.xml diff --git a/tests/test_salphabeta/geometry.xml b/tests/test_salphabeta/geometry.xml deleted file mode 100644 index 2b978b914..000000000 --- a/tests/test_salphabeta/geometry.xml +++ /dev/null @@ -1,15 +0,0 @@ - - - - - - - - - - - - - - - diff --git a/tests/test_salphabeta/inputs_true.dat b/tests/test_salphabeta/inputs_true.dat new file mode 100644 index 000000000..c42b90e8c --- /dev/null +++ b/tests/test_salphabeta/inputs_true.dat @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 5 + 0 + + + -4 -4 -4 4 4 4 + + + diff --git a/tests/test_salphabeta/materials.xml b/tests/test_salphabeta/materials.xml deleted file mode 100644 index bfe0a6224..000000000 --- a/tests/test_salphabeta/materials.xml +++ /dev/null @@ -1,41 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/test_salphabeta/results_true.dat b/tests/test_salphabeta/results_true.dat index bdfb98990..e1121bce5 100644 --- a/tests/test_salphabeta/results_true.dat +++ b/tests/test_salphabeta/results_true.dat @@ -1,2 +1,2 @@ k-combined: -8.544160E-01 1.274133E-02 +8.403447E-01 2.461538E-02 diff --git a/tests/test_salphabeta/settings.xml b/tests/test_salphabeta/settings.xml deleted file mode 100644 index 70b4e802f..000000000 --- a/tests/test_salphabeta/settings.xml +++ /dev/null @@ -1,15 +0,0 @@ - - - - eigenvalue - 10 - 5 - 1000 - - - - -4 -4 -4 4 4 4 - - - - diff --git a/tests/test_salphabeta/test_salphabeta.py b/tests/test_salphabeta/test_salphabeta.py index b04fcc6eb..600c70332 100644 --- a/tests/test_salphabeta/test_salphabeta.py +++ b/tests/test_salphabeta/test_salphabeta.py @@ -3,9 +3,82 @@ import os import sys sys.path.insert(0, os.pardir) -from testing_harness import TestHarness + +from testing_harness import PyAPITestHarness +import openmc +import openmc.model + + +def make_model(): + model = openmc.model.Model() + + # Materials + m1 = openmc.Material() + m1.set_density('g/cc', 4.5) + m1.add_nuclide('U235', 1.0) + m1.add_nuclide('H1', 1.0) + m1.add_s_alpha_beta('c_H_in_H2O', fraction=0.5) + + m2 = openmc.Material() + m2.set_density('g/cc', 4.5) + m2.add_nuclide('U235', 1.0) + m2.add_nuclide('C0', 1.0) + m2.add_s_alpha_beta('c_Graphite') + + m3 = openmc.Material() + m3.set_density('g/cc', 4.5) + m3.add_nuclide('U235', 1.0) + m3.add_nuclide('Be9', 1.0) + m3.add_nuclide('O16', 1.0) + m3.add_s_alpha_beta('c_Be_in_BeO') + m3.add_s_alpha_beta('c_O_in_BeO') + + m4 = openmc.Material() + m4.set_density('g/cm3', 5.90168) + m4.add_nuclide('H1', 0.3) + m4.add_nuclide('Zr90', 0.15) + m4.add_nuclide('Zr91', 0.1) + m4.add_nuclide('Zr92', 0.1) + m4.add_nuclide('Zr94', 0.05) + m4.add_nuclide('Zr96', 0.05) + m4.add_nuclide('U235', 0.1) + m4.add_nuclide('U238', 0.15) + m4.add_s_alpha_beta('c_Zr_in_ZrH') + m4.add_s_alpha_beta('c_H_in_ZrH') + + model.materials += [m1, m2, m3, m4] + + # Geometry + x0 = openmc.XPlane(x0=-10, boundary_type='vacuum') + x1 = openmc.XPlane(x0=-5) + x2 = openmc.XPlane(x0=0) + x3 = openmc.XPlane(x0=5) + x4 = openmc.XPlane(x0=10, boundary_type='vacuum') + + root_univ = openmc.Universe() + + surfs = (x0, x1, x2, x3, x4) + mats = (m1, m2, m3, m4) + cells = [] + for i in range(4): + cell = openmc.Cell() + cell.region = +surfs[i] & -surfs[i+1] + cell.fill = mats[i] + root_univ.add_cell(cell) + + model.geometry.root_universe = root_univ + + # Settings + model.settings.batches = 5 + model.settings.inactive = 0 + model.settings.particles = 1000 + model.settings.source = openmc.Source(space=openmc.stats.Box( + [-4, -4, -4], [4, 4, 4])) + + return model if __name__ == '__main__': - harness = TestHarness('statepoint.10.h5') + model = make_model() + harness = PyAPITestHarness('statepoint.5.h5', model) harness.main() From 7ba35310c8da200d896d418933d30237b3c1d414 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 29 Jun 2017 15:49:18 -0400 Subject: [PATCH 05/12] Update docs, RelaxNG for partial S(a,b) --- docs/source/io_formats/materials.rst | 6 +++++- src/relaxng/materials.rnc | 3 ++- src/relaxng/materials.rng | 28 ++++++++++++++++++++-------- 3 files changed, 27 insertions(+), 10 deletions(-) diff --git a/docs/source/io_formats/materials.rst b/docs/source/io_formats/materials.rst index fc57280f8..286342a2b 100644 --- a/docs/source/io_formats/materials.rst +++ b/docs/source/io_formats/materials.rst @@ -107,9 +107,13 @@ Each ``material`` element can have the following attributes or sub-elements: multi-group :ref:`energy_mode`. :sab: - Associates an S(a,b) table with the material. This element has one + Associates an S(a,b) table with the material. This element has an attribute/sub-element called ``name``. The ``name`` attribute is the name of the S(a,b) table that should be associated with the material. + There is also an optional ``fraction`` element which indicates what fraction + of the relevant nuclides will be affected by the S(a,b) table (e.g. which + fraction of a material is crystaline versus amorphus). ``fraction`` + defaults to unity. *Default*: None diff --git a/src/relaxng/materials.rnc b/src/relaxng/materials.rnc index 4ad3ea0f0..d55656536 100644 --- a/src/relaxng/materials.rnc +++ b/src/relaxng/materials.rnc @@ -28,7 +28,8 @@ element materials { }* & element sab { - (element name { xsd:string } | attribute name { xsd:string }) + (element name { xsd:string } | attribute name { xsd:string }) & + (element fraction { xsd:double } | attribute fraction { xsd:double })? }* }+ & diff --git a/src/relaxng/materials.rng b/src/relaxng/materials.rng index 0303d361e..0fa81c11b 100644 --- a/src/relaxng/materials.rng +++ b/src/relaxng/materials.rng @@ -119,14 +119,26 @@ - - - - - - - - + + + + + + + + + + + + + + + + + + + + From 8dcd8dd4f0d20c1e8a1729b0fa699ecfa24ad52e Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 7 Jul 2017 15:55:19 -0400 Subject: [PATCH 06/12] Fix spelling errors --- docs/source/io_formats/materials.rst | 2 +- openmc/material.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/io_formats/materials.rst b/docs/source/io_formats/materials.rst index 286342a2b..5f0aa3a5e 100644 --- a/docs/source/io_formats/materials.rst +++ b/docs/source/io_formats/materials.rst @@ -112,7 +112,7 @@ Each ``material`` element can have the following attributes or sub-elements: is the name of the S(a,b) table that should be associated with the material. There is also an optional ``fraction`` element which indicates what fraction of the relevant nuclides will be affected by the S(a,b) table (e.g. which - fraction of a material is crystaline versus amorphus). ``fraction`` + fraction of a material is crystalline versus amorphous). ``fraction`` defaults to unity. *Default*: None diff --git a/openmc/material.py b/openmc/material.py index 6ac153c86..24cd3343a 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -628,7 +628,7 @@ class Material(IDManagerMixin): fraction : float The fraction of relevant nuclei that are affected by the :math:`S(\alpha,\beta)` table. For example, if the material is a - block of carbon that is 60% graphite and 40% amorphus then add a + block of carbon that is 60% graphite and 40% amorphous then add a graphite :math:`S(\alpha,\beta)` table with fraction=0.6. """ From e2fcd68d1f4dffaa739fe3f6bda98d6913146e2b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 10 Jul 2017 14:43:18 -0400 Subject: [PATCH 07/12] Fatal error for nuclides in multiple S(a,b) tables --- src/input_xml.F90 | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2e7219a31..ef742ed45 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -5193,6 +5193,19 @@ contains end if end do ASSIGN_SAB + ! Make sure each nuclide only appears in one table. + do j = 1, i_sab_nuclides % size() + do k = j+1, i_sab_nuclides % size() + if (i_sab_nuclides % data(j) == i_sab_nuclides % data(k)) then + call fatal_error(trim( & + nuclides(mat % nuclide(i_sab_nuclides % data(j))) % name) & + // " in material " // trim(to_str(mat % id)) // " was found & + &in multiple S(a,b) tables. Each nuclide can only appear in & + &one S(a,b) table per material.") + end if + end do + end do + ! Update i_sab_tables and i_sab_nuclides deallocate(mat % i_sab_tables) deallocate(mat % sab_fracs) From 27d1a1b56251e1bca366ad3191ba584a3df11244 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 12 Jul 2017 11:40:21 -0400 Subject: [PATCH 08/12] Fix sorting of S(a,b) fractions --- src/input_xml.F90 | 6 +++++- tests/test_salphabeta/results_true.dat | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3aae18f19..aa41dc4b5 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -5167,6 +5167,7 @@ contains integer :: m ! position for sorting integer :: temp_nuclide ! temporary value for sorting integer :: temp_table ! temporary value for sorting + real(8) :: temp_frac ! temporary value for sorting logical :: found type(VectorInt) :: i_sab_tables type(VectorInt) :: i_sab_nuclides @@ -5224,7 +5225,7 @@ contains allocate(mat % sab_fracs(m)) mat % i_sab_tables(:) = i_sab_tables % data(1:m) mat % i_sab_nuclides(:) = i_sab_nuclides % data(1:m) - mat % sab_fracs = sab_fracs % data(1:m) + mat % sab_fracs(:) = sab_fracs % data(1:m) ! Clear entries in vectors for next material call i_sab_tables % clear() @@ -5242,6 +5243,7 @@ contains m = k temp_nuclide = mat % i_sab_nuclides(k) temp_table = mat % i_sab_tables(k) + temp_frac = mat % i_sab_tables(k) MOVE_OVER: do ! Check if insertion value is greater than (m-1)th value @@ -5250,6 +5252,7 @@ contains ! Move values over until hitting one that's not larger mat % i_sab_nuclides(m) = mat % i_sab_nuclides(m-1) mat % i_sab_tables(m) = mat % i_sab_tables(m-1) + mat % sab_fracs(m) = mat % sab_fracs(m-1) m = m - 1 ! Exit if we've reached the beginning of the list @@ -5259,6 +5262,7 @@ contains ! Put the original value into its new position mat % i_sab_nuclides(m) = temp_nuclide mat % i_sab_tables(m) = temp_table + mat % sab_fracs(m) = temp_frac end do SORT_SAB end if diff --git a/tests/test_salphabeta/results_true.dat b/tests/test_salphabeta/results_true.dat index e1121bce5..2a99303f6 100644 --- a/tests/test_salphabeta/results_true.dat +++ b/tests/test_salphabeta/results_true.dat @@ -1,2 +1,2 @@ k-combined: -8.403447E-01 2.461538E-02 +8.447580E-01 1.806149E-02 From 391cb442ef22984af7459ab60e485a313f06693b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 12 Jul 2017 11:53:50 -0400 Subject: [PATCH 09/12] Only put S(a,b) fraction in XML if not unity --- openmc/material.py | 3 ++- tests/test_asymmetric_lattice/inputs_true.dat | 18 +++++++++--------- tests/test_diff_tally/inputs_true.dat | 18 +++++++++--------- tests/test_filter_energyfun/inputs_true.dat | 18 +++++++++--------- tests/test_filter_mesh/inputs_true.dat | 18 +++++++++--------- tests/test_iso_in_lab/inputs_true.dat | 18 +++++++++--------- .../test_mgxs_library_ce_to_mg/inputs_true.dat | 2 +- .../test_mgxs_library_condense/inputs_true.dat | 2 +- .../inputs_true.dat | 2 +- tests/test_mgxs_library_hdf5/inputs_true.dat | 2 +- tests/test_mgxs_library_mesh/inputs_true.dat | 18 +++++++++--------- .../inputs_true.dat | 2 +- .../test_mgxs_library_nuclides/inputs_true.dat | 2 +- tests/test_multipole/inputs_true.dat | 2 +- tests/test_periodic/inputs_true.dat | 2 +- tests/test_salphabeta/inputs_true.dat | 10 +++++----- tests/test_tallies/inputs_true.dat | 18 +++++++++--------- tests/test_tally_aggregation/inputs_true.dat | 18 +++++++++--------- tests/test_tally_arithmetic/inputs_true.dat | 18 +++++++++--------- tests/test_tally_slice_merge/inputs_true.dat | 18 +++++++++--------- tests/test_triso/inputs_true.dat | 8 ++++---- tests/test_volume_calc/inputs_true.dat | 2 +- 22 files changed, 110 insertions(+), 109 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 24cd3343a..03ef7fbc4 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -982,7 +982,8 @@ class Material(IDManagerMixin): for sab in self._sab: subelement = ET.SubElement(element, "sab") subelement.set("name", sab[0]) - subelement.set("fraction", str(sab[1])) + if sab[1] != 1.0: + subelement.set("fraction", str(sab[1])) return element diff --git a/tests/test_asymmetric_lattice/inputs_true.dat b/tests/test_asymmetric_lattice/inputs_true.dat index 50803105d..672f6bee1 100644 --- a/tests/test_asymmetric_lattice/inputs_true.dat +++ b/tests/test_asymmetric_lattice/inputs_true.dat @@ -78,7 +78,7 @@ - + @@ -86,7 +86,7 @@ - + @@ -114,7 +114,7 @@ - + @@ -129,7 +129,7 @@ - + @@ -144,7 +144,7 @@ - + @@ -159,7 +159,7 @@ - + @@ -174,7 +174,7 @@ - + @@ -187,7 +187,7 @@ - + @@ -200,7 +200,7 @@ - + diff --git a/tests/test_diff_tally/inputs_true.dat b/tests/test_diff_tally/inputs_true.dat index a2c5638cf..ff134bd5c 100644 --- a/tests/test_diff_tally/inputs_true.dat +++ b/tests/test_diff_tally/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_filter_energyfun/inputs_true.dat b/tests/test_filter_energyfun/inputs_true.dat index d6f7c1f79..d857451cb 100644 --- a/tests/test_filter_energyfun/inputs_true.dat +++ b/tests/test_filter_energyfun/inputs_true.dat @@ -171,7 +171,7 @@ - + @@ -179,7 +179,7 @@ - + @@ -207,7 +207,7 @@ - + @@ -222,7 +222,7 @@ - + @@ -237,7 +237,7 @@ - + @@ -252,7 +252,7 @@ - + @@ -267,7 +267,7 @@ - + @@ -280,7 +280,7 @@ - + @@ -293,7 +293,7 @@ - + diff --git a/tests/test_filter_mesh/inputs_true.dat b/tests/test_filter_mesh/inputs_true.dat index 153f4dcd9..6d14f9e7e 100644 --- a/tests/test_filter_mesh/inputs_true.dat +++ b/tests/test_filter_mesh/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_iso_in_lab/inputs_true.dat b/tests/test_iso_in_lab/inputs_true.dat index 2311d8588..098056f5b 100644 --- a/tests/test_iso_in_lab/inputs_true.dat +++ b/tests/test_iso_in_lab/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_mgxs_library_ce_to_mg/inputs_true.dat b/tests/test_mgxs_library_ce_to_mg/inputs_true.dat index 675ef50f1..8e8cde281 100644 --- a/tests/test_mgxs_library_ce_to_mg/inputs_true.dat +++ b/tests/test_mgxs_library_ce_to_mg/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_mgxs_library_condense/inputs_true.dat b/tests/test_mgxs_library_condense/inputs_true.dat index 38a40c1cc..d2f28d0fe 100644 --- a/tests/test_mgxs_library_condense/inputs_true.dat +++ b/tests/test_mgxs_library_condense/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_mgxs_library_distribcell/inputs_true.dat b/tests/test_mgxs_library_distribcell/inputs_true.dat index 8b6522bdd..d7a4a186a 100644 --- a/tests/test_mgxs_library_distribcell/inputs_true.dat +++ b/tests/test_mgxs_library_distribcell/inputs_true.dat @@ -60,7 +60,7 @@ - + diff --git a/tests/test_mgxs_library_hdf5/inputs_true.dat b/tests/test_mgxs_library_hdf5/inputs_true.dat index 38a40c1cc..d2f28d0fe 100644 --- a/tests/test_mgxs_library_hdf5/inputs_true.dat +++ b/tests/test_mgxs_library_hdf5/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_mgxs_library_mesh/inputs_true.dat b/tests/test_mgxs_library_mesh/inputs_true.dat index 8bfcc64b6..4756b27bf 100644 --- a/tests/test_mgxs_library_mesh/inputs_true.dat +++ b/tests/test_mgxs_library_mesh/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_mgxs_library_no_nuclides/inputs_true.dat b/tests/test_mgxs_library_no_nuclides/inputs_true.dat index 38a40c1cc..d2f28d0fe 100644 --- a/tests/test_mgxs_library_no_nuclides/inputs_true.dat +++ b/tests/test_mgxs_library_no_nuclides/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_mgxs_library_nuclides/inputs_true.dat b/tests/test_mgxs_library_nuclides/inputs_true.dat index 03e3b3ff6..826eb5b62 100644 --- a/tests/test_mgxs_library_nuclides/inputs_true.dat +++ b/tests/test_mgxs_library_nuclides/inputs_true.dat @@ -33,7 +33,7 @@ - + diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 4b349ecba..0d1fe99bd 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -25,7 +25,7 @@ - + diff --git a/tests/test_periodic/inputs_true.dat b/tests/test_periodic/inputs_true.dat index 26e700004..61958701b 100644 --- a/tests/test_periodic/inputs_true.dat +++ b/tests/test_periodic/inputs_true.dat @@ -16,7 +16,7 @@ - + diff --git a/tests/test_salphabeta/inputs_true.dat b/tests/test_salphabeta/inputs_true.dat index c42b90e8c..02e6813f0 100644 --- a/tests/test_salphabeta/inputs_true.dat +++ b/tests/test_salphabeta/inputs_true.dat @@ -22,15 +22,15 @@ - + - - + + @@ -42,8 +42,8 @@ - - + + diff --git a/tests/test_tallies/inputs_true.dat b/tests/test_tallies/inputs_true.dat index 664e6e983..a85491e94 100644 --- a/tests/test_tallies/inputs_true.dat +++ b/tests/test_tallies/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_tally_aggregation/inputs_true.dat b/tests/test_tally_aggregation/inputs_true.dat index 6ef6bdb4d..7a8bf4613 100644 --- a/tests/test_tally_aggregation/inputs_true.dat +++ b/tests/test_tally_aggregation/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_tally_arithmetic/inputs_true.dat b/tests/test_tally_arithmetic/inputs_true.dat index e9c156287..743ca3c58 100644 --- a/tests/test_tally_arithmetic/inputs_true.dat +++ b/tests/test_tally_arithmetic/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_tally_slice_merge/inputs_true.dat b/tests/test_tally_slice_merge/inputs_true.dat index 2c90e475a..b1c089c96 100644 --- a/tests/test_tally_slice_merge/inputs_true.dat +++ b/tests/test_tally_slice_merge/inputs_true.dat @@ -170,7 +170,7 @@ - + @@ -178,7 +178,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -221,7 +221,7 @@ - + @@ -236,7 +236,7 @@ - + @@ -251,7 +251,7 @@ - + @@ -266,7 +266,7 @@ - + @@ -279,7 +279,7 @@ - + @@ -292,7 +292,7 @@ - + diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat index 862c2db9a..fdbc1cb5f 100644 --- a/tests/test_triso/inputs_true.dat +++ b/tests/test_triso/inputs_true.dat @@ -403,12 +403,12 @@ - + - + @@ -420,12 +420,12 @@ - + - + diff --git a/tests/test_volume_calc/inputs_true.dat b/tests/test_volume_calc/inputs_true.dat index 35c866028..28f1cbd9f 100644 --- a/tests/test_volume_calc/inputs_true.dat +++ b/tests/test_volume_calc/inputs_true.dat @@ -16,7 +16,7 @@ - + From 53f38806d8bec0ac6cd21b18991ba216c453d6e3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 12 Jul 2017 14:39:24 -0400 Subject: [PATCH 10/12] Restructure S(a,b) treatment in NuclideMicroXS Also make sure cross sections are re-evaluated if sab_frac changes --- src/cross_section.F90 | 52 +++++++++++++++++++++--------------------- src/nuclide_header.F90 | 45 ++++++++++++++++++------------------ src/physics.F90 | 8 +++---- src/tally.F90 | 3 +-- 4 files changed, 54 insertions(+), 54 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1306fb669..1222624d4 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -61,22 +61,24 @@ contains ! Add contribution from each nuclide in material do i = 1, mat % n_nuclides - ! ======================================================================== + ! ====================================================================== ! CHECK FOR S(A,B) TABLE i_sab = 0 - ! Check if this nuclide matches one of the S(a,b) tables specified -- this - ! relies on i_sab_nuclides being in sorted order + ! Check if this nuclide matches one of the S(a,b) tables specified. + ! This relies on i_sab_nuclides being in sorted order if (check_sab) then if (i == mat % i_sab_nuclides(j)) then ! Get index in sab_tables i_sab = mat % i_sab_tables(j) sab_frac = mat % sab_fracs(j) - ! If particle energy is greater than the highest energy for the S(a,b) - ! table, don't use the S(a,b) table - if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) i_sab = 0 + ! If particle energy is greater than the highest energy for the + ! S(a,b) table, then don't use the S(a,b) table + if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then + i_sab = 0 + end if ! Increment position in i_sab_nuclides j = j + 1 @@ -94,10 +96,9 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E & - .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, & - p % sqrtkT, sab_frac) - else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then + .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT & + .or. i_sab /= micro_xs(i_nuclide) % index_sab & + .or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, & p % sqrtkT, sab_frac) end if @@ -114,8 +115,7 @@ contains ! Add contributions to material macroscopic scattering cross section material_xs % elastic = material_xs % elastic + & - atom_density * (micro_xs(i_nuclide) % elastic & - + micro_xs(i_nuclide) % scatter_sab) + atom_density * micro_xs(i_nuclide) % elastic ! Add contributions to material macroscopic absorption cross section material_xs % absorption = material_xs % absorption + & @@ -250,10 +250,10 @@ contains micro_xs(i_nuclide) % interp_factor = f ! Initialize nuclide cross-sections to zero - micro_xs(i_nuclide) % fission = ZERO - micro_xs(i_nuclide) % nu_fission = ZERO - micro_xs(i_nuclide) % scatter_sab = ZERO - micro_xs(i_nuclide) % elastic_sab = ZERO + micro_xs(i_nuclide) % fission = ZERO + micro_xs(i_nuclide) % nu_fission = ZERO + micro_xs(i_nuclide) % thermal = ZERO + micro_xs(i_nuclide) % thermal_elastic = ZERO ! Calculate microscopic nuclide total cross section micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) & @@ -280,11 +280,10 @@ contains end if ! Initialize sab treatment to false - micro_xs(i_nuclide) % index_sab = NONE - micro_xs(i_nuclide) % elastic_sab = ZERO + micro_xs(i_nuclide) % index_sab = NONE ! Initialize URR probability table treatment to false - micro_xs(i_nuclide) % use_ptable = .false. + micro_xs(i_nuclide) % use_ptable = .false. ! If there is S(a,b) data for this nuclide, we need to set the sab_scatter ! and sab_elastic cross sections and correct the total and elastic cross @@ -305,7 +304,6 @@ contains end if micro_xs(i_nuclide) % last_E = E - micro_xs(i_nuclide) % last_index_sab = i_sab micro_xs(i_nuclide) % last_sqrtkT = sqrtkT end associate @@ -414,17 +412,19 @@ contains end associate ! Store the S(a,b) cross sections. - micro_xs(i_nuclide) % scatter_sab = sab_frac * (elastic + inelastic) - micro_xs(i_nuclide) % elastic_sab = sab_frac * elastic + micro_xs(i_nuclide) % thermal = sab_frac * (elastic + inelastic) + micro_xs(i_nuclide) % thermal_elastic = sab_frac * elastic ! Correct total and elastic cross sections micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total & - + sab_frac * (elastic + inelastic - micro_xs(i_nuclide) % elastic) - micro_xs(i_nuclide) % elastic = & - (ONE - sab_frac) * micro_xs(i_nuclide) % elastic + + micro_xs(i_nuclide) % thermal & + - sab_frac * micro_xs(i_nuclide) % elastic + micro_xs(i_nuclide) % elastic = micro_xs(i_nuclide) % thermal & + + (ONE - sab_frac) * micro_xs(i_nuclide) % elastic - ! Save temperature index + ! Save temperature index and thermal fraction micro_xs(i_nuclide) % index_temp_sab = i_temp + micro_xs(i_nuclide) % sab_frac = sab_frac end subroutine calculate_sab_xs diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 41b1008a8..5c5131abc 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -102,34 +102,35 @@ module nuclide_header end type Nuclide !=============================================================================== -! NUCLIDEMICROXS contains cached microscopic cross sections for a -! particular nuclide at the current energy +! NUCLIDEMICROXS contains cached microscopic cross sections for a particular +! nuclide at the current energy !=============================================================================== type NuclideMicroXS - integer :: index_grid ! index on nuclide energy grid - integer :: index_temp ! temperature index for nuclide - real(8) :: last_E = ZERO ! last evaluated energy - real(8) :: interp_factor ! interpolation factor on nuc. energy grid - real(8) :: total ! microscropic total xs - real(8) :: elastic ! microscopic elastic scattering xs (non S(a,b)) - real(8) :: absorption ! microscopic absorption xs - real(8) :: fission ! microscopic fission xs - real(8) :: nu_fission ! microscopic production xs - real(8) :: scatter_sab ! microscopic scattering xs due to S(a,b) + ! Microscopic cross sections in barns + real(8) :: total + real(8) :: elastic ! If sab_frac is not 1 or 0, then this value is + ! averaged over bound and non-bound nuclei + real(8) :: absorption + real(8) :: fission + real(8) :: nu_fission + real(8) :: thermal ! Bound thermal elastic & inelastic scattering + real(8) :: thermal_elastic ! Bound thermal elastic scattering - ! Information for S(a,b) use - integer :: index_sab ! index in sab_tables (zero means no table) - integer :: last_index_sab = 0 ! index in sab_tables last used by this nuclide - integer :: index_temp_sab ! temperature index for sab_tables - real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table + ! Indicies and factors needed to compute cross sections from the data tables + integer :: index_grid ! Index on nuclide energy grid + integer :: index_temp ! Temperature index for nuclide + real(8) :: interp_factor ! Interpolation factor on nuc. energy grid + integer :: index_sab = NONE ! Index in sab_tables + integer :: index_temp_sab ! Temperature index for sab_tables + real(8) :: sab_frac ! Fraction of atoms affected by S(a,b) + logical :: use_ptable ! In URR range with probability tables? - ! Information for URR probability table use - logical :: use_ptable ! in URR range with probability tables? - - ! Information for Doppler broadening + ! Energy and temperature last used to evaluate these cross sections. If + ! these values have changed, then the cross sections must be re-evaluated. + real(8) :: last_E = ZERO ! Last evaluated energy real(8) :: last_sqrtkT = ZERO ! Last temperature in sqrt(Boltzmann - ! constant * temperature (eV)) + ! constant * temperature (eV)) end type NuclideMicroXS !=============================================================================== diff --git a/src/physics.F90 b/src/physics.F90 index 15b53fd7f..3d557b44f 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -342,7 +342,7 @@ contains micro_xs(i_nuclide) % absorption) sampled = .false. - prob = micro_xs(i_nuclide) % elastic + prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal if (prob > cutoff) then ! ======================================================================= ! NON-S(A,B) ELASTIC SCATTERING @@ -362,7 +362,7 @@ contains sampled = .true. end if - prob = prob + micro_xs(i_nuclide) % scatter_sab + prob = micro_xs(i_nuclide) % elastic if (prob > cutoff .and. .not. sampled) then ! ======================================================================= ! S(A,B) SCATTERING @@ -567,8 +567,8 @@ contains associate (sab => sab_tables(i_sab) % data(i_temp)) ! Determine whether inelastic or elastic scattering will occur - if (prn() < micro_xs(i_nuclide) % elastic_sab / & - micro_xs(i_nuclide) % scatter_sab) then + if (prn() < micro_xs(i_nuclide) % thermal_elastic / & + micro_xs(i_nuclide) % thermal) then ! elastic scattering ! Get index and interpolation factor for elastic grid diff --git a/src/tally.F90 b/src/tally.F90 index d827d9bac..0ee041cbe 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1051,8 +1051,7 @@ contains else if (i_nuclide > 0) then - score = (micro_xs(i_nuclide) % elastic & - + micro_xs(i_nuclide) % scatter_sab) * atom_density * flux + score = micro_xs(i_nuclide) % elastic * atom_density * flux else score = material_xs % elastic * flux end if From 4ce4235099e41867d25dd739efe242887eb38b17 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 13 Jul 2017 12:16:46 -0400 Subject: [PATCH 11/12] Make sure micro_xs % sab_frac is defined Otherwise, cross sections can be unnecessarily re-evaluated --- src/cross_section.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1222624d4..0b3c1fdb0 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -281,6 +281,7 @@ contains ! Initialize sab treatment to false micro_xs(i_nuclide) % index_sab = NONE + micro_xs(i_nuclide) % sab_frac = ZERO ! Initialize URR probability table treatment to false micro_xs(i_nuclide) % use_ptable = .false. From 0ced86c595947c64b72386fc41d2fe53938a21e1 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 13 Jul 2017 12:31:06 -0400 Subject: [PATCH 12/12] Make sure sab_frac is defined --- src/cross_section.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 0b3c1fdb0..68955f2df 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -65,6 +65,7 @@ contains ! CHECK FOR S(A,B) TABLE i_sab = 0 + sab_frac = ZERO ! Check if this nuclide matches one of the S(a,b) tables specified. ! This relies on i_sab_nuclides being in sorted order