diff --git a/docs/source/io_formats/materials.rst b/docs/source/io_formats/materials.rst index fc57280f8..5f0aa3a5e 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 crystalline versus amorphous). ``fraction`` + defaults to unity. *Default*: None diff --git a/openmc/material.py b/openmc/material.py index 762332584..03ef7fbc4 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% amorphous 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,9 @@ 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]) + if sab[1] != 1.0: + subelement.set("fraction", str(sab[1])) return element diff --git a/src/cross_section.F90 b/src/cross_section.F90 index ce7d97133..68955f2df 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 @@ -60,21 +61,25 @@ contains ! Add contribution from each nuclide in material do i = 1, mat % n_nuclides - ! ======================================================================== + ! ====================================================================== ! 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 + ! 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 @@ -84,7 +89,7 @@ contains end if end if - ! ======================================================================== + ! ====================================================================== ! CALCULATE MICROSCOPIC CROSS SECTION ! Determine microscopic cross sections for this nuclide @@ -92,13 +97,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) - 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) + .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 - ! ======================================================================== + ! ====================================================================== ! ADD TO MACROSCOPIC CROSS SECTION ! Copy atom density of nuclide in material @@ -133,13 +139,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 +251,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) % 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) & @@ -271,20 +281,21 @@ 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 + micro_xs(i_nuclide) % sab_frac = ZERO ! 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 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 @@ -295,7 +306,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 @@ -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,16 +413,20 @@ contains end if end associate + ! Store the S(a,b) cross sections. + 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 - & - micro_xs(i_nuclide) % elastic + inelastic + elastic - micro_xs(i_nuclide) % elastic = inelastic + elastic + micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total & + + 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 - ! Store S(a,b) elastic cross section for sampling later - micro_xs(i_nuclide) % elastic_sab = 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/input_xml.F90 b/src/input_xml.F90 index 3e1370843..aa41dc4b5 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) & @@ -5159,9 +5167,11 @@ 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 + 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 @@ -5179,6 +5189,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 @@ -5192,17 +5203,34 @@ 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) 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 @@ -5215,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 @@ -5223,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 @@ -5232,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/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..5c5131abc 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -102,33 +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 - real(8) :: absorption ! microscopic absorption xs - real(8) :: fission ! microscopic fission xs - real(8) :: nu_fission ! microscopic production xs + ! 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 9482772ef..3d557b44f 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 - micro_xs(i_nuclide) % thermal 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 = micro_xs(i_nuclide) % elastic + 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 @@ -559,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) % elastic) 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/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 @@ - - - - - - - - + + + + + + + + + + + + + + + + + + + + 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..02e6813f0 --- /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..2a99303f6 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.447580E-01 1.806149E-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()