From 199d4af9c5e7a274faf36363c9e9e10076376bcf Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 8 Nov 2015 17:24:16 -0500 Subject: [PATCH] Add diff tally for keff---density or nuc. density --- src/constants.F90 | 10 +- src/initialize.F90 | 1 + src/input_xml.F90 | 60 ++++++- src/material_header.F90 | 1 + src/output.F90 | 1 + src/tally.F90 | 110 ++++++++++++ src/tally_header.F90 | 18 ++ src/tracking.F90 | 9 +- tests/test_diff_density/build_xml.py | 150 ++++++++++++++++ tests/test_diff_density/geometry.xml | 16 ++ tests/test_diff_density/materials.xml | 45 +++++ tests/test_diff_density/results_true.dat | 5 + tests/test_diff_density/settings.xml | 18 ++ tests/test_diff_density/tallies.xml | 7 + tests/test_diff_density/test_diff_density.py | 11 ++ tests/test_diff_nuclide_density/build_xml.py | 164 ++++++++++++++++++ tests/test_diff_nuclide_density/geometry.xml | 16 ++ tests/test_diff_nuclide_density/materials.xml | 50 ++++++ .../results_true.dat | 5 + tests/test_diff_nuclide_density/settings.xml | 18 ++ tests/test_diff_nuclide_density/tallies.xml | 10 ++ .../test_diff_nuclide_density.py | 11 ++ 22 files changed, 732 insertions(+), 4 deletions(-) create mode 100644 tests/test_diff_density/build_xml.py create mode 100644 tests/test_diff_density/geometry.xml create mode 100644 tests/test_diff_density/materials.xml create mode 100644 tests/test_diff_density/results_true.dat create mode 100644 tests/test_diff_density/settings.xml create mode 100644 tests/test_diff_density/tallies.xml create mode 100755 tests/test_diff_density/test_diff_density.py create mode 100644 tests/test_diff_nuclide_density/build_xml.py create mode 100644 tests/test_diff_nuclide_density/geometry.xml create mode 100644 tests/test_diff_nuclide_density/materials.xml create mode 100644 tests/test_diff_nuclide_density/results_true.dat create mode 100644 tests/test_diff_nuclide_density/settings.xml create mode 100644 tests/test_diff_nuclide_density/tallies.xml create mode 100755 tests/test_diff_nuclide_density/test_diff_nuclide_density.py diff --git a/src/constants.F90 b/src/constants.F90 index 375c517e76..8e8e2db76c 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -259,7 +259,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 22 + integer, parameter :: N_SCORE_TYPES = 23 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -282,7 +282,8 @@ module constants SCORE_NU_SCATTER_YN = -19, & ! angular flux-weighted nu-scattering moment (0:N) SCORE_EVENTS = -20, & ! number of events SCORE_DELAYED_NU_FISSION = -21, & ! delayed neutron production rate - SCORE_INVERSE_VELOCITY = -22 ! flux-weighted inverse velocity + SCORE_INVERSE_VELOCITY = -22, & ! flux-weighted inverse velocity + SCORE_KEFF = -23 ! multiplication factor ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 @@ -348,6 +349,11 @@ module constants K_TRACKLENGTH = 3, & LEAKAGE = 4 + ! Differential tally dependent variables + integer, parameter :: & + DIFF_DENSITY = 1, & + DIFF_NUCLIDE_DENSITY = 2 + ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS diff --git a/src/initialize.F90 b/src/initialize.F90 index 9d63eb4e93..6130be0e25 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -831,6 +831,7 @@ contains ! Change density in g/cm^3 to atom/b-cm. Since all values are now in atom ! percent, the sum needs to be re-evaluated as 1/sum(x*awr) if (.not. density_in_atom) then + mat % density_gpcc = -mat % density sum_percent = ZERO do j = 1, mat%n_nuclides index_list = xs_listing_dict%get_key(mat%names(j)) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 58e4013355..052748d139 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2197,6 +2197,7 @@ contains type(Node), pointer :: node_tal => null() type(Node), pointer :: node_filt => null() type(Node), pointer :: node_trigger=>null() + type(Node), pointer :: node_deriv => null() type(NodeList), pointer :: node_mesh_list => null() type(NodeList), pointer :: node_tal_list => null() type(NodeList), pointer :: node_filt_list => null() @@ -2407,7 +2408,7 @@ contains t % estimator = ESTIMATOR_TRACKLENGTH - ! Copy material id + ! Copy tally id if (check_for_node(node_tal, "id")) then call get_node_value(node_tal, "id", t % id) else @@ -2843,6 +2844,52 @@ contains end do end if + ! ======================================================================= + ! READ DATA FOR DERIVATIVES + + if (check_for_node(node_tal, "derivative")) then + call get_node_ptr(node_tal, "derivative", node_deriv) + allocate(t % deriv) + + temp_str = "" + call get_node_value(node_deriv, "variable", temp_str) + temp_str = to_lower(temp_str) + + select case(temp_str) + + case("density") + t % deriv % dep_var = DIFF_DENSITY + call get_node_value(node_deriv, "material", t % deriv % diff_material) + + case("nuclide_density") + t % deriv % dep_var = DIFF_NUCLIDE_DENSITY + call get_node_value(node_deriv, "material", t % deriv % diff_material) + + call get_node_value(node_deriv, "nuclide", word) + word = trim(to_lower(word)) + pair_list => nuclide_dict % keys() + do while (associated(pair_list)) + if (starts_with(pair_list % key, word)) then + word = pair_list % key(1:150) + exit + end if + + ! Advance to next + pair_list => pair_list % next + end do + + ! Check if no nuclide was found + if (.not. associated(pair_list)) then + call fatal_error("Could not find the nuclide " & + &// trim(word) // " specified in tally " & + &// trim(to_str(t % id)) // " in any material.") + end if + deallocate(pair_list) + + t % deriv % diff_nuclide = nuclide_dict % get_key(word) + end select + end if + ! ======================================================================= ! READ DATA FOR SCORES @@ -2991,6 +3038,11 @@ contains call fatal_error("Cannot tally flux with an outgoing energy & &filter.") end if + + if (allocated(t % deriv)) then + t % estimator = ESTIMATOR_COLLISION + end if + case ('flux-yn') ! Prohibit user from tallying flux for an individual nuclide if (.not. (t % n_nuclide_bins == 1 .and. & @@ -3188,6 +3240,12 @@ contains case ('events') t % score_bins(j) = SCORE_EVENTS + case ('keff') + t % score_bins(j) = SCORE_KEFF + if (allocated(t % deriv)) then + t % estimator = ESTIMATOR_COLLISION + end if + case ('elastic', '(n,elastic)') t % score_bins(j) = ELASTIC case ('(n,2nd)') diff --git a/src/material_header.F90 b/src/material_header.F90 index c2f4b51bc0..24e459cbb5 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -13,6 +13,7 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm + real(8) :: density_gpcc ! total density in g/cm^3 ! Energy grid information integer :: n_grid ! # of union material grid points diff --git a/src/output.F90 b/src/output.F90 index 55cd5b2a5b..0f667c121b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -986,6 +986,7 @@ contains score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment" score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate" score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" + score_names(abs(SCORE_KEFF)) = "k_eff, multiplaction factor" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/tally.F90 b/src/tally.F90 index fcd6dfdd9b..7fce90cce5 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -654,6 +654,13 @@ contains end if end if + case (SCORE_KEFF) + if (t % estimator == ESTIMATOR_COLLISION) then + score = p % last_wgt * material_xs % nu_fission / material_xs % total + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + score = flux * material_xs % nu_fission + end if + case default if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need @@ -732,6 +739,35 @@ contains end select + !######################################################################### + ! Add derivative information on score for differential tallies. + + if (allocated(t % deriv) .and. t % estimator == ESTIMATOR_COLLISION) then + select case (score_bin) + + case (SCORE_KEFF) + select case (t % deriv % dep_var) + + case (DIFF_DENSITY) + score = score * t % deriv % accumulator + + case (DIFF_NUCLIDE_DENSITY) + mat => materials(p % material) + if (mat % id == t % deriv % diff_material .and. & + micro_xs(t % deriv % diff_nuclide) % nu_fission /= ZERO) then + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == t % deriv % diff_nuclide) exit + end do + score = score * (t % deriv % accumulator + ONE & + / mat % atom_density(l)) + else + score = score * t % deriv % accumulator + end if + + end select + end select + end if + !######################################################################### ! Expand score if necessary and add to tally results. @@ -2210,6 +2246,80 @@ contains end subroutine score_surface_current +!=============================================================================== +! SCORE_DIFF_ACCUMULATORS +!=============================================================================== + + subroutine score_diff_accumulators(p, distance, collision) + type(particle), intent(in) :: p + real(8), intent(in) :: distance + logical, intent(in) :: collision + + integer :: i, j + type(TallyObject), pointer :: t + type(Material), pointer :: mat + + do i = 1, n_user_tallies + t => user_tallies(i) + + if (allocated(t % deriv)) then + select case (t % deriv % dep_var) + + case (DIFF_DENSITY) + if (p % material == MATERIAL_VOID) cycle + mat => materials(p % material) + if (mat % id == t % deriv % diff_material) then + if (collision) then + t % deriv % accumulator = t % deriv % accumulator & + + ONE / mat % density_gpcc & + - distance * material_xs % total / mat % density_gpcc + else + t % deriv % accumulator = t % deriv % accumulator & + - distance * material_xs % total / mat % density_gpcc + end if + end if + + case (DIFF_NUCLIDE_DENSITY) + if (p % material == MATERIAL_VOID) cycle + mat => materials(p % material) + if (mat % id == t % deriv % diff_material) then + do j = 1, mat % n_nuclides + if (mat % nuclide(j) == t % deriv % diff_nuclide) exit + end do + if (mat % nuclide(j) /= t % deriv % diff_nuclide) then + call fatal_error("Couldn't find the right nuclide.") + end if + t % deriv % accumulator = t % deriv % accumulator & + - distance * micro_xs(t % deriv % diff_nuclide) % total + if (collision) then + if (p % event == EVENT_SCATTER) then + t % deriv % accumulator = t % deriv % accumulator & + + micro_xs(t % deriv % diff_nuclide) % elastic & + / material_xs % elastic + else if (p % event == EVENT_ABSORB) then + t % deriv % accumulator = t % deriv % accumulator & + + micro_xs(t % deriv % diff_nuclide) % absorption & + / material_xs % absorption + end if + end if + end if + end select + end if + end do + end subroutine + + subroutine clear_diff_accumulators() + integer :: i + type(TallyObject), pointer :: t + + do i = 1, n_user_tallies + t => user_tallies(i) + if (allocated(t % deriv)) then + t % deriv % accumulator = ZERO + end if + end do + end subroutine + !=============================================================================== ! GET_NEXT_BIN determines the next scoring bin for a particular filter variable !=============================================================================== diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 01dcd9bdb5..cfbe0d5376 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -64,6 +64,18 @@ module tally_header procedure :: clear => tallyfilter_clear ! Deallocates TallyFilter end type TallyFilter + +!=============================================================================== +! TALLYDERIVATIVE +!=============================================================================== + + type TallyDerivative + real(8) :: accumulator + integer :: dep_var + integer :: diff_material + integer :: diff_nuclide + end type TallyDerivative + !=============================================================================== ! TALLYOBJECT describes a user-specified tally. The region of phase space to ! tally in is given by the TallyFilters and the results are stored in a @@ -130,6 +142,9 @@ module tally_header integer :: n_triggers = 0 ! # of triggers type(TriggerObject), allocatable :: triggers(:) ! Array of triggers + ! Derivative for differentially tallies + type(TallyDerivative), allocatable :: deriv + ! Type-Bound procedures contains procedure :: clear => tallyobject_clear ! Deallocates TallyObject @@ -204,6 +219,9 @@ module tally_header this % n_triggers = 0 + if (allocated(this % deriv)) & + deallocate (this % deriv) + end subroutine tallyobject_clear end module tally_header diff --git a/src/tracking.F90 b/src/tracking.F90 index 2e4503e139..9502746da7 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -13,7 +13,8 @@ module tracking use random_lcg, only: prn use string, only: to_str use tally, only: score_analog_tally, score_tracklength_tally, & - score_collision_tally, score_surface_current + score_collision_tally, score_surface_current, & + score_diff_accumulators, clear_diff_accumulators use track_output, only: initialize_particle_track, write_particle_track, & add_particle_track, finalize_particle_track @@ -60,6 +61,9 @@ contains call initialize_particle_track() endif + ! Every particle starts with no accumulated flux derivative. + call clear_diff_accumulators() + EVENT_LOOP: do ! If the cell hasn't been determined based on the particle's location, ! initiate a search for the current cell. This generally happens at the @@ -114,6 +118,9 @@ contains distance * material_xs % nu_fission end if + ! Score flux derivative accumulators for differential tallies. + call score_diff_accumulators(p, distance, d_boundary > d_collision) + if (d_collision > d_boundary) then ! ==================================================================== ! PARTICLE CROSSES SURFACE diff --git a/tests/test_diff_density/build_xml.py b/tests/test_diff_density/build_xml.py new file mode 100644 index 0000000000..d5e50c7e82 --- /dev/null +++ b/tests/test_diff_density/build_xml.py @@ -0,0 +1,150 @@ +import openmc + + +def make_mats(mod_density=0.7420582): + water = openmc.Material(name='water') + water.set_density('g/cm3', mod_density) + water.add_nuclide(openmc.Nuclide('H-1'), 2.0) + water.add_nuclide(openmc.Nuclide('O-16'), 1.0) + + fuel = openmc.Material(name='fuel 2.4%') + fuel.set_density('g/cm3', 10.29769) + fuel.add_nuclide(openmc.Nuclide('U-234'), 5.7987e-06) + fuel.add_nuclide(openmc.Nuclide('U-235'), 7.2175e-04) + fuel.add_nuclide(openmc.Nuclide('U-238'), 2.2253e-02) + fuel.add_nuclide(openmc.Nuclide('O-16'), 4.5750e-02) + fuel.add_nuclide(openmc.Nuclide('O-17'), 9.4222e-05) + + zirc4 = openmc.Material(name='zircaloy 4') + zirc4.set_density('g/cm3', 6.55) + zirc4.add_nuclide(openmc.Nuclide('O-16'), 3.0743e-04) + zirc4.add_nuclide(openmc.Nuclide('O-17'), 1.1711e-07) + # O-18 omitted + zirc4.add_nuclide(openmc.Nuclide('Cr-50'), 3.2962e-06) + zirc4.add_nuclide(openmc.Nuclide('Cr-52'), 6.3564e-05) + zirc4.add_nuclide(openmc.Nuclide('Cr-53'), 7.2076e-06) + zirc4.add_nuclide(openmc.Nuclide('Cr-54'), 1.7941e-06) + zirc4.add_nuclide(openmc.Nuclide('Fe-54'), 8.6699e-06) + zirc4.add_nuclide(openmc.Nuclide('Fe-56'), 1.3610e-04) + zirc4.add_nuclide(openmc.Nuclide('Fe-57'), 3.1431e-06) + zirc4.add_nuclide(openmc.Nuclide('Fe-58'), 4.1829e-07) + zirc4.add_nuclide(openmc.Nuclide('Zr-90'), 2.1827e-02) + zirc4.add_nuclide(openmc.Nuclide('Zr-91'), 4.7600e-03) + zirc4.add_nuclide(openmc.Nuclide('Zr-92'), 7.2758e-03) + zirc4.add_nuclide(openmc.Nuclide('Zr-94'), 7.3734e-03) + zirc4.add_nuclide(openmc.Nuclide('Zr-96'), 1.1879e-03) + zirc4.add_nuclide(openmc.Nuclide('Sn-112'), 4.6735e-06) + zirc4.add_nuclide(openmc.Nuclide('Sn-114'), 3.1799e-06) + zirc4.add_nuclide(openmc.Nuclide('Sn-115'), 1.6381e-06) + zirc4.add_nuclide(openmc.Nuclide('Sn-116'), 7.0055e-05) + zirc4.add_nuclide(openmc.Nuclide('Sn-117'), 3.7003e-05) + zirc4.add_nuclide(openmc.Nuclide('Sn-118'), 1.1669e-04) + zirc4.add_nuclide(openmc.Nuclide('Sn-119'), 4.1387e-05) + zirc4.add_nuclide(openmc.Nuclide('Sn-120'), 1.5697e-04) + zirc4.add_nuclide(openmc.Nuclide('Sn-122'), 2.2308e-05) + zirc4.add_nuclide(openmc.Nuclide('Sn-124'), 2.7897e-05) + + materials_file = openmc.MaterialsFile() + materials_file.default_xs = '71c' + materials_file.add_materials([water, fuel, zirc4]) + + return (materials_file, + {'mod':water, + 'fuel':fuel, + 'clad':zirc4}) + + +def make_geom(mats): + # Instantiate surfaces. + rfo = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.39218, name='fuel outer') + rci = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.40005, name='clad inner') + rco = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720, name='clad outer') + + x0 = openmc.XPlane(x0=-0.62992) + x1 = openmc.XPlane(x0=0.62992) + y0 = openmc.YPlane(y0=-0.62992) + y1 = openmc.YPlane(y0=0.62992) + z0 = openmc.ZPlane(z0=-1.0) + z1 = openmc.ZPlane(z0=1.0) + x0.boundary_type = 'reflective' + x1.boundary_type = 'reflective' + y0.boundary_type = 'reflective' + y1.boundary_type = 'reflective' + z0.boundary_type = 'reflective' + z1.boundary_type = 'reflective' + + # Instantiate cells. + fuel_c = openmc.Cell(name='fuel') + fuel_c.add_surface(rfo, halfspace=-1) + fuel_c.add_surface(z0, halfspace=+1) + fuel_c.add_surface(z1, halfspace=-1) + + gap = openmc.Cell(name='gap') + gap.add_surface(rci, halfspace=-1) + gap.add_surface(rfo, halfspace=+1) + gap.add_surface(z0, halfspace=+1) + gap.add_surface(z1, halfspace=-1) + + clad_c = openmc.Cell(name='clad') + clad_c.add_surface(rco, halfspace=-1) + clad_c.add_surface(rci, halfspace=+1) + clad_c.add_surface(z0, halfspace=+1) + clad_c.add_surface(z1, halfspace=-1) + + mod_c = openmc.Cell(name='moderator') + mod_c.add_surface(rco, halfspace=+1) + mod_c.add_surface(x0, halfspace=+1) + mod_c.add_surface(x1, halfspace=-1) + mod_c.add_surface(y0, halfspace=+1) + mod_c.add_surface(y1, halfspace=-1) + mod_c.add_surface(z0, halfspace=+1) + mod_c.add_surface(z1, halfspace=-1) + + # Add materials to cells. + fuel_c.fill = mats['fuel'] + gap.fill = 'void' + clad_c.fill = mats['clad'] + mod_c.fill = mats['mod'] + + # Instantiate universes. + u0 = openmc.Universe(universe_id=0) + u0.add_cells([fuel_c, gap, clad_c, mod_c]) + + # Write the XML file. + geometry = openmc.Geometry() + geometry.root_universe = u0 + geometry_file = openmc.GeometryFile() + geometry_file.geometry = geometry + + return geometry_file + + +def make_settings(**kwargs): + if 'batches' in kwargs: + batches = kwargs['batches'] + + settings_file = openmc.SettingsFile() + settings_file.batches = kwargs.setdefault('batches', 100) + settings_file.inactive = kwargs.setdefault('inactive', 10) + settings_file.particles = kwargs.setdefault('particles', 1000) + settings_file.set_source_space('box', (-0.6, -0.6, -0.9, + 0.6, 0.6, 0.9)) + settings_file.entropy_lower_left = (-0.7, -0.7, -1.1) + settings_file.entropy_upper_right = (0.7, 0.7, 1.1) + settings_file.entropy_dimension = (10, 10, 10) + return settings_file + + +def make_inputs(mod_density=0.7420582, **kwargs): + (mat_file, mats) = make_mats(mod_density) + mat_file.export_to_xml() + + geo_file = make_geom(mats) + geo_file.export_to_xml() + + sets_file = make_settings(**kwargs) + sets_file.export_to_xml() + + +if __name__ == '__main__': + make_inputs(inactive=5, batches=10, particles=100) diff --git a/tests/test_diff_density/geometry.xml b/tests/test_diff_density/geometry.xml new file mode 100644 index 0000000000..e7e5b097ba --- /dev/null +++ b/tests/test_diff_density/geometry.xml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/tests/test_diff_density/materials.xml b/tests/test_diff_density/materials.xml new file mode 100644 index 0000000000..955ea92e7d --- /dev/null +++ b/tests/test_diff_density/materials.xml @@ -0,0 +1,45 @@ + + + 71c + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_diff_density/results_true.dat b/tests/test_diff_density/results_true.dat new file mode 100644 index 0000000000..5e515c4613 --- /dev/null +++ b/tests/test_diff_density/results_true.dat @@ -0,0 +1,5 @@ +k-combined: +1.502605E+00 2.273838E-02 +tally 1: +5.280104E+00 +1.379258E+01 diff --git a/tests/test_diff_density/settings.xml b/tests/test_diff_density/settings.xml new file mode 100644 index 0000000000..a41ba9725b --- /dev/null +++ b/tests/test_diff_density/settings.xml @@ -0,0 +1,18 @@ + + + + 100 + 10 + 5 + + + + -0.6 -0.6 -0.9 0.6 0.6 0.9 + + + + 10 10 10 + -0.7 -0.7 -1.1 + 0.7 0.7 1.1 + + diff --git a/tests/test_diff_density/tallies.xml b/tests/test_diff_density/tallies.xml new file mode 100644 index 0000000000..214f33ab37 --- /dev/null +++ b/tests/test_diff_density/tallies.xml @@ -0,0 +1,7 @@ + + + + keff + + + diff --git a/tests/test_diff_density/test_diff_density.py b/tests/test_diff_density/test_diff_density.py new file mode 100755 index 0000000000..ed6addec45 --- /dev/null +++ b/tests/test_diff_density/test_diff_density.py @@ -0,0 +1,11 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness + + +if __name__ == '__main__': + harness = TestHarness('statepoint.10.*', True) + harness.main() diff --git a/tests/test_diff_nuclide_density/build_xml.py b/tests/test_diff_nuclide_density/build_xml.py new file mode 100644 index 0000000000..f77a337cb1 --- /dev/null +++ b/tests/test_diff_nuclide_density/build_xml.py @@ -0,0 +1,164 @@ +import openmc + + +def make_mats(**kwargs): + water_dict = { + 'B-10': 8.0042e-06, + 'B-11': 3.2218e-05, + 'H-1': 4.9457e-02, + 'H-2': 7.4196e-06, + 'O-16': 2.4672e-02, + 'O-17': 9.3982e-06 + 5.0701e-05} + water = openmc.Material(name='water') + water.set_density('sum') + for nuclide in water_dict: + water.add_nuclide(nuclide, water_dict[nuclide]) + + U235_dens = kwargs.setdefault('U235_dens', 5.5814e-04) + fuel_dict = { + 'B-10': 8.0042e-06, + 'U-234': 4.4842e-06, + 'U-235': U235_dens, + 'U-238': 2.2407e-02, + 'O-16': 4.5828e-02, + 'O-17': 1.7457e-05 + 9.4176e-05} + fuel = openmc.Material(name='fuel 2.4%') + fuel.set_density('sum') + for nuclide in fuel_dict: + fuel.add_nuclide(nuclide, fuel_dict[nuclide]) + + zirc4_dict = { + 'O-16': 3.0743e-04, + 'O-17': 1.1711e-07 + 6.3176e-07, + 'Cr-50': 3.2962e-06, + 'Cr-52': 6.3564e-05, + 'Cr-53': 7.2076e-06, + 'Cr-54': 1.7941e-06, + 'Fe-54': 8.6699e-06, + 'Fe-56': 1.3610e-04, + 'Fe-57': 3.1431e-06, + 'Fe-58': 4.1829e-07, + 'Zr-90': 2.1827e-02, + 'Zr-91': 4.7600e-03, + 'Zr-92': 7.2758e-03, + 'Zr-94': 7.3734e-03, + 'Zr-96': 1.1879e-03, + 'Sn-112': 4.6735e-06, + 'Sn-114': 3.1799e-06, + 'Sn-115': 1.6381e-06, + 'Sn-116': 7.0055e-05, + 'Sn-117': 3.7003e-05, + 'Sn-118': 1.1669e-04, + 'Sn-119': 4.1387e-05, + 'Sn-120': 1.5697e-04, + 'Sn-122': 2.2308e-05, + 'Sn-124': 2.7897e-05} + zirc4 = openmc.Material(name='zircaloy 4') + zirc4.set_density('sum') + for nuclide in zirc4_dict: + zirc4.add_nuclide(nuclide, zirc4_dict[nuclide]) + + materials_file = openmc.MaterialsFile() + materials_file.default_xs = '71c' + materials_file.add_materials([water, fuel, zirc4]) + + return (materials_file, + {'mod':water, + 'fuel':fuel, + 'clad':zirc4}) + + +def make_geom(mats): + # Instantiate surfaces. + rfo = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.39218, name='fuel outer') + rci = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.40005, name='clad inner') + rco = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720, name='clad outer') + + x0 = openmc.XPlane(x0=-0.62992) + x1 = openmc.XPlane(x0=0.62992) + y0 = openmc.YPlane(y0=-0.62992) + y1 = openmc.YPlane(y0=0.62992) + z0 = openmc.ZPlane(z0=-1.0) + z1 = openmc.ZPlane(z0=1.0) + x0.boundary_type = 'reflective' + x1.boundary_type = 'reflective' + y0.boundary_type = 'reflective' + y1.boundary_type = 'reflective' + z0.boundary_type = 'reflective' + z1.boundary_type = 'reflective' + + # Instantiate cells. + fuel_c = openmc.Cell(name='fuel') + fuel_c.add_surface(rfo, halfspace=-1) + fuel_c.add_surface(z0, halfspace=+1) + fuel_c.add_surface(z1, halfspace=-1) + + gap = openmc.Cell(name='gap') + gap.add_surface(rci, halfspace=-1) + gap.add_surface(rfo, halfspace=+1) + gap.add_surface(z0, halfspace=+1) + gap.add_surface(z1, halfspace=-1) + + clad_c = openmc.Cell(name='clad') + clad_c.add_surface(rco, halfspace=-1) + clad_c.add_surface(rci, halfspace=+1) + clad_c.add_surface(z0, halfspace=+1) + clad_c.add_surface(z1, halfspace=-1) + + mod_c = openmc.Cell(name='moderator') + mod_c.add_surface(rco, halfspace=+1) + mod_c.add_surface(x0, halfspace=+1) + mod_c.add_surface(x1, halfspace=-1) + mod_c.add_surface(y0, halfspace=+1) + mod_c.add_surface(y1, halfspace=-1) + mod_c.add_surface(z0, halfspace=+1) + mod_c.add_surface(z1, halfspace=-1) + + # Add materials to cells. + fuel_c.fill = mats['fuel'] + gap.fill = 'void' + clad_c.fill = mats['clad'] + mod_c.fill = mats['mod'] + + # Instantiate universes. + u0 = openmc.Universe(universe_id=0) + u0.add_cells([fuel_c, gap, clad_c, mod_c]) + + # Write the XML file. + geometry = openmc.Geometry() + geometry.root_universe = u0 + geometry_file = openmc.GeometryFile() + geometry_file.geometry = geometry + + return geometry_file + + +def make_settings(**kwargs): + if 'batches' in kwargs: + batches = kwargs['batches'] + + settings_file = openmc.SettingsFile() + settings_file.batches = kwargs.setdefault('batches', 1610) + settings_file.inactive = kwargs.setdefault('inactive', 10) + settings_file.particles = kwargs.setdefault('particles', 1000) + settings_file.set_source_space('box', (-0.6, -0.6, -0.9, + 0.6, 0.6, 0.9)) + settings_file.entropy_lower_left = (-0.7, -0.7, -1.1) + settings_file.entropy_upper_right = (0.7, 0.7, 1.1) + settings_file.entropy_dimension = (10, 10, 10) + return settings_file + + +def make_inputs(**kwargs): + (mat_file, mats) = make_mats(**kwargs) + mat_file.export_to_xml() + + geo_file = make_geom(mats) + geo_file.export_to_xml() + + sets_file = make_settings(**kwargs) + sets_file.export_to_xml() + + +if __name__ == '__main__': + make_inputs(inactive=5, batches=10, particles=100) diff --git a/tests/test_diff_nuclide_density/geometry.xml b/tests/test_diff_nuclide_density/geometry.xml new file mode 100644 index 0000000000..e7e5b097ba --- /dev/null +++ b/tests/test_diff_nuclide_density/geometry.xml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/tests/test_diff_nuclide_density/materials.xml b/tests/test_diff_nuclide_density/materials.xml new file mode 100644 index 0000000000..02b55f4523 --- /dev/null +++ b/tests/test_diff_nuclide_density/materials.xml @@ -0,0 +1,50 @@ + + + 71c + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_diff_nuclide_density/results_true.dat b/tests/test_diff_nuclide_density/results_true.dat new file mode 100644 index 0000000000..68b611ae15 --- /dev/null +++ b/tests/test_diff_nuclide_density/results_true.dat @@ -0,0 +1,5 @@ +k-combined: +1.055198E+00 5.785829E-02 +tally 1: +3.081784E+03 +1.905482E+06 diff --git a/tests/test_diff_nuclide_density/settings.xml b/tests/test_diff_nuclide_density/settings.xml new file mode 100644 index 0000000000..a41ba9725b --- /dev/null +++ b/tests/test_diff_nuclide_density/settings.xml @@ -0,0 +1,18 @@ + + + + 100 + 10 + 5 + + + + -0.6 -0.6 -0.9 0.6 0.6 0.9 + + + + 10 10 10 + -0.7 -0.7 -1.1 + 0.7 0.7 1.1 + + diff --git a/tests/test_diff_nuclide_density/tallies.xml b/tests/test_diff_nuclide_density/tallies.xml new file mode 100644 index 0000000000..321b2c9b0e --- /dev/null +++ b/tests/test_diff_nuclide_density/tallies.xml @@ -0,0 +1,10 @@ + + + + keff + + + + diff --git a/tests/test_diff_nuclide_density/test_diff_nuclide_density.py b/tests/test_diff_nuclide_density/test_diff_nuclide_density.py new file mode 100755 index 0000000000..ed6addec45 --- /dev/null +++ b/tests/test_diff_nuclide_density/test_diff_nuclide_density.py @@ -0,0 +1,11 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness + + +if __name__ == '__main__': + harness = TestHarness('statepoint.10.*', True) + harness.main()