From 6ef6740ebc01f97c4fdc7772ba1b6fe29789e07b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Jul 2016 11:49:25 -0500 Subject: [PATCH 01/12] Implement stochastic volume calculation on Fortran side --- docs/source/io_formats/index.rst | 1 + docs/source/io_formats/nuclear_data.rst | 2 +- docs/source/io_formats/volume.rst | 22 ++ docs/source/usersguide/input.rst | 29 ++ src/constants.F90 | 3 +- src/global.F90 | 3 + src/hdf5_interface.F90 | 76 +++++ src/input_xml.F90 | 10 + src/random_lcg.F90 | 1 - src/simulation.F90 | 4 + src/volume_calc.F90 | 389 ++++++++++++++++++++++++ src/volume_header.F90 | 42 +++ 12 files changed, 579 insertions(+), 3 deletions(-) create mode 100644 docs/source/io_formats/volume.rst create mode 100644 src/volume_calc.F90 create mode 100644 src/volume_header.F90 diff --git a/docs/source/io_formats/index.rst b/docs/source/io_formats/index.rst index acab7e8930..1da8872c66 100644 --- a/docs/source/io_formats/index.rst +++ b/docs/source/io_formats/index.rst @@ -30,3 +30,4 @@ Output Files particle_restart track voxel + volume diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index ba6a54eb1e..9d7ff0eb1e 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -1,4 +1,4 @@ -.. _usersguide_nuclear_data: +.. _io_nuclear_data: ======================== Nuclear Data File Format diff --git a/docs/source/io_formats/volume.rst b/docs/source/io_formats/volume.rst new file mode 100644 index 0000000000..10b7a3b732 --- /dev/null +++ b/docs/source/io_formats/volume.rst @@ -0,0 +1,22 @@ +.. _io_volume: + +================== +Volume File Format +================== + +**/** + +:Attributes: - **samples** (*int*) -- Number of samples + - **lower_left** (*double[3]*) -- Lower-left coordinates of + bounding box + - **upper_right** (*double[3]*) -- Upper-right coordinates of + bounding box + +**/cell_/** + +:Datasets: - **volume** (*double[2]*) -- Calculated volume and its uncertainty + in cubic centimeters + - **nuclides** (*char[][]*) -- Names of nuclides identified in the + cell + - **atoms** (*double[][2]*) -- Total number of atoms of each nuclide + and its uncertainty diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index e530097ded..530c45d2fe 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -836,6 +836,35 @@ displayed. This element takes the following attributes: *Default*: 5 +```` Element +------------------------- + +The ```` element indicates that a stochastic volume calculation +should be run at the beginning of the simulation. This element has the following +sub-elements/attributes: + + :cells: + The unique IDs of cells for which the volume should be estimated. + + *Default*: None + + :samples: + The number of samples used to estimate volumes. + + *Default*: None + + :lower_left: + The lower-left Cartesian coordinates of a bounding box that is used to + sample points within. + + *Default*: None + + :upper_right: + The upper-right Cartesian coordinates of a bounding box that is used to + sample points within. + + *Default*: None + -------------------------------------- Geometry Specification -- geometry.xml -------------------------------------- diff --git a/src/constants.F90 b/src/constants.F90 index a22c9ac050..076e1f0d0b 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -380,11 +380,12 @@ module constants ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS - integer, parameter :: N_STREAMS = 4 + integer, parameter :: N_STREAMS = 5 integer, parameter :: STREAM_TRACKING = 1 integer, parameter :: STREAM_TALLIES = 2 integer, parameter :: STREAM_SOURCE = 3 integer, parameter :: STREAM_URR_PTABLE = 4 + integer, parameter :: STREAM_VOLUME = 5 ! ============================================================================ ! MISCELLANEOUS CONSTANTS diff --git a/src/global.F90 b/src/global.F90 index 1b1e5632a4..ba5947f0e7 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -17,6 +17,7 @@ module global use tally_header, only: TallyObject, TallyMap, TallyResult use trigger_header, only: KTrigger use timer_header, only: Timer + use volume_header, only: VolumeCalculation #ifdef MPIF08 use mpi_f08 @@ -35,6 +36,8 @@ module global type(Material), allocatable, target :: materials(:) type(ObjectPlot), allocatable, target :: plots(:) + type(VolumeCalculation), allocatable :: volume_calcs(:) + ! Size of main arrays integer :: n_cells ! # of cells integer :: n_universes ! # of universes diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 5a1e2a2fff..6d1c87d7f5 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -75,8 +75,15 @@ module hdf5_interface module procedure read_attribute_string end interface read_attribute + interface write_attribute + module procedure write_attribute_double + module procedure write_attribute_double_1D + module procedure write_attribute_integer + end interface write_attribute + public :: write_dataset public :: read_dataset + public :: write_attribute public :: read_attribute public :: file_create public :: file_open @@ -2059,6 +2066,25 @@ contains call h5aclose_f(attr_id, hdf5_err) end subroutine read_attribute_double + subroutine write_attribute_double(obj_id, name, buffer) + integer(HID_T), intent(in) :: obj_id + character(*), intent(in) :: name + real(8), intent(in), target :: buffer + + integer :: hdf5_err + integer(HID_T) :: dspace_id + integer(HID_T) :: attr_id + type(C_PTR) :: f_ptr + + call h5screate_f(H5S_SCALAR_F, dspace_id, hdf5_err) + call h5acreate_f(obj_id, trim(name), H5T_NATIVE_DOUBLE, dspace_id, & + attr_id, hdf5_err) + f_ptr = c_loc(buffer) + call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) + call h5aclose_f(attr_id, hdf5_err) + call h5sclose_f(dspace_id, hdf5_err) + end subroutine write_attribute_double + subroutine read_attribute_double_1D(buffer, obj_id, name) real(8), target, allocatable, intent(inout) :: buffer(:) integer(HID_T), intent(in) :: obj_id @@ -2097,6 +2123,37 @@ contains call h5aread_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) end subroutine read_attribute_double_1D_explicit + subroutine write_attribute_double_1D(obj_id, name, buffer) + integer(HID_T), intent(in) :: obj_id + character(*), intent(in) :: name + real(8), target, intent(in) :: buffer(:) + + integer(HSIZE_T) :: dims(1) + + dims(:) = shape(buffer) + call write_attribute_double_1D_explicit(obj_id, dims, name, buffer) + end subroutine write_attribute_double_1D + + subroutine write_attribute_double_1D_explicit(obj_id, dims, name, buffer) + integer(HID_T), intent(in) :: obj_id + integer(HSIZE_T), intent(in) :: dims(1) + character(*), intent(in) :: name + real(8), target, intent(in) :: buffer(dims(1)) + + integer :: hdf5_err + integer(HID_T) :: dspace_id + integer(HID_T) :: attr_id + type(C_PTR) :: f_ptr + + call h5screate_simple_f(1, dims, dspace_id, hdf5_err) + call h5acreate_f(obj_id, trim(name), H5T_NATIVE_DOUBLE, dspace_id, & + attr_id, hdf5_err) + f_ptr = c_loc(buffer) + call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err) + call h5aclose_f(attr_id, hdf5_err) + call h5sclose_f(dspace_id, hdf5_err) + end subroutine write_attribute_double_1D_explicit + subroutine read_attribute_double_2D(buffer, obj_id, name) real(8), target, allocatable, intent(inout) :: buffer(:,:) integer(HID_T), intent(in) :: obj_id @@ -2150,6 +2207,25 @@ contains call h5aclose_f(attr_id, hdf5_err) end subroutine read_attribute_integer + subroutine write_attribute_integer(obj_id, name, buffer) + integer(HID_T), intent(in) :: obj_id + character(*), intent(in) :: name + integer, intent(in), target :: buffer + + integer :: hdf5_err + integer(HID_T) :: dspace_id + integer(HID_T) :: attr_id + type(C_PTR) :: f_ptr + + call h5screate_f(H5S_SCALAR_F, dspace_id, hdf5_err) + call h5acreate_f(obj_id, trim(name), H5T_NATIVE_INTEGER, dspace_id, & + attr_id, hdf5_err) + f_ptr = c_loc(buffer) + call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err) + call h5aclose_f(attr_id, hdf5_err) + call h5sclose_f(dspace_id, hdf5_err) + end subroutine write_attribute_integer + subroutine read_attribute_integer_1D(buffer, obj_id, name) integer, target, allocatable, intent(inout) :: buffer(:) integer(HID_T), intent(in) :: obj_id diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3e5a687124..0c2eca33a5 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -87,8 +87,10 @@ contains type(Node), pointer :: node_scatterer => null() type(Node), pointer :: node_trigger => null() type(Node), pointer :: node_keff_trigger => null() + type(Node), pointer :: node_vol => null() type(NodeList), pointer :: node_scat_list => null() type(NodeList), pointer :: node_source_list => null() + type(NodeList), pointer :: node_vol_list => null() ! Check if settings.xml exists filename = trim(path_input) // "settings.xml" @@ -1119,6 +1121,14 @@ contains end select end if + call get_node_list(doc, "volume_calc", node_vol_list) + n = get_list_size(node_vol_list) + allocate(volume_calcs(n)) + do i = 1, n + call get_list_item(node_vol_list, i, node_vol) + call volume_calcs(i) % from_xml(node_vol) + end do + ! Close settings XML file call close_xmldoc(doc) diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 08f1034ab7..287bbba6d9 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -29,7 +29,6 @@ module random_lcg public :: set_particle_seed public :: advance_prn_seed public :: prn_set_stream - public :: STREAM_TRACKING, STREAM_TALLIES contains diff --git a/src/simulation.F90 b/src/simulation.F90 index 3321dc70fa..f59ce371cc 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -24,6 +24,7 @@ module simulation reset_result use trigger, only: check_triggers use tracking, only: transport + use volume_calc, only: run_volume_calculations implicit none private @@ -42,6 +43,9 @@ contains type(Particle) :: p integer(8) :: i_work + ! Volume calculations + if (size(volume_calcs) > 0) call run_volume_calculations() + if (.not. restart_run) call initialize_source() ! Display header diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 new file mode 100644 index 0000000000..cb36f0b7f3 --- /dev/null +++ b/src/volume_calc.F90 @@ -0,0 +1,389 @@ +module volume_calc + + use hdf5, only: HID_T +#ifdef _OPENMP + use omp_lib +#endif + + use constants + use geometry, only: find_cell + use global + use hdf5_interface, only: file_create, file_close, write_attribute, & + create_group, close_group, write_dataset + use output, only: write_message, header + use message_passing + use particle_header, only: Particle + use random_lcg, only: prn, prn_set_stream, set_particle_seed + use stl_vector, only: VectorInt, VectorReal + use timer_header, only: Timer + use volume_header + + implicit none + private + + public :: run_volume_calculations + +contains + +!=============================================================================== +! RUN_VOLUME_CALCULATIONS runs each of the stochastic volume calculations that +! the user has specified and writes results to HDF5 files +!=============================================================================== + + subroutine run_volume_calculations() + integer :: i, j + integer :: n + real(8), allocatable :: volume(:,:) ! volume mean/stdev in each cell + character(MAX_FILE_LEN) :: filename ! filename for HDF5 file + type(Timer) :: time_volume ! timer for volume calculation + type(VectorInt), allocatable :: nuclide_vec(:) ! indices in nuclides array + type(VectorReal), allocatable :: atoms_vec(:) ! total # of atoms of each nuclide + type(VectorReal), allocatable :: uncertainty_vec(:) ! uncertainty of total # of atoms + + if (master) then + call header("STOCHASTIC VOLUME CALCULATION", level=1) + call time_volume % start() + end if + + do i = 1, size(volume_calcs) + n = size(volume_calcs(i) % cell_id) + allocate(nuclide_vec(n)) + allocate(atoms_vec(n), uncertainty_vec(n)) + allocate(volume(2,n)) + + if (master) then + call write_message("Running volume calculation " // trim(to_str(i)) & + // "...") + end if + + call get_volume(volume_calcs(i), volume, nuclide_vec, atoms_vec, & + uncertainty_vec) + + if (master) then + ! Display cell volumes + do j = 1, size(volume_calcs(i) % cell_id) + call write_message(" Cell " // trim(to_str(volume_calcs(i) % & + cell_id(j))) // ": " // trim(to_str(volume(1,j))) // " +/- " // & + trim(to_str(volume(2,j))) // " cm^3") + end do + call write_message("") + + filename = trim(path_output) // 'volume_' // trim(to_str(i)) // '.h5' + call write_volume(volume_calcs(i), filename, volume, nuclide_vec, & + atoms_vec, uncertainty_vec) + end if + + deallocate(nuclide_vec, atoms_vec, uncertainty_vec, volume) + end do + + ! Show elapsed time + if (master) then + call time_volume % stop() + call write_message("Elapsed time: " // trim(to_str(time_volume % & + get_value())) // " s") + end if + end subroutine run_volume_calculations + +!=============================================================================== +! GET_VOLUME stochastically determines the volume of a set of cells along with +! the average number densities of nuclides within the cell +!=============================================================================== + + subroutine get_volume(this, volume, nuclide_vec, atoms_vec, uncertainty_vec) + type(VolumeCalculation), intent(in) :: this + real(8), intent(out) :: volume(:,:) ! volume mean/stdev in each cell + type(VectorInt), intent(out) :: nuclide_vec(:) ! indices in nuclides array + type(VectorReal), intent(out) :: atoms_vec(:) ! total # of atoms of each nuclide + type(VectorReal), intent(out) :: uncertainty_vec(:) ! uncertainty of total # of atoms + + ! Variables that are private to each thread + integer(8) :: i + integer :: j, k + integer :: i_cell ! index in cell_id array + integer :: i_material ! index in materials array + integer :: level ! local coordinate level + logical :: found_cell + type(VectorInt) :: indices(size(this % cell_id)) ! List of material indices + type(VectorInt) :: hits(size(this % cell_id)) ! Number of hits for each material + type(Particle) :: p + + ! Shared variables + integer :: i_start, i_end ! Starting/ending sample for each process + type(VectorInt) :: master_indices(size(this % cell_id)) + type(VectorInt) :: master_hits(size(this % cell_id)) + + ! Variables used outside of parallel region + integer :: i_nuclide ! index in nuclides array + integer :: total_hits ! total hits for a single cell (summed over materials) + integer :: min_samples ! minimum number of samples per process + integer :: remainder ! leftover samples from uneven divide +#ifdef MPI + integer :: m ! index over materials + integer :: n ! number of materials + integer, allocatable :: data(:) ! array used to send number of hits +#endif + real(8) :: f ! fraction of hits + real(8) :: var_f ! variance of fraction of hits + real(8) :: volume_sample ! total volume of sampled region + real(8) :: atoms(2, size(nuclides)) + + ! Divide work over MPI processes + min_samples = this % samples / n_procs + remainder = mod(this % samples, n_procs) + if (rank < remainder) then + i_start = (min_samples + 1)*rank + i_end = i_start + min_samples + else + i_start = (min_samples + 1)*remainder + (rank - remainder)*min_samples + i_end = i_start + min_samples - 1 + end if + + call p % initialize() + +!$omp parallel private(i, j, k, i_cell, i_material, level, found_cell) & +!$omp& firstprivate(p, indices, hits) + + call prn_set_stream(STREAM_VOLUME) + + ! ========================================================================== + ! SAMPLES LOCATIONS AND COUNT HITS + +!$omp do + SAMPLE_LOOP: do i = i_start, i_end + call set_particle_seed(i) + + p % n_coord = 1 + p % coord(1) % xyz(1) = this % lower_left(1) + prn()*(& + this % upper_right(1) - this % lower_left(1)) + p % coord(1) % xyz(2) = this % lower_left(2) + prn()*(& + this % upper_right(2) - this % lower_left(2)) + p % coord(1) % xyz(3) = this % lower_left(3) + prn()*(& + this % upper_right(3) - this % lower_left(3)) + p % coord(1) % uvw(:) = [HALF, HALF, HALF] + + ! If this location is not in the geometry at all, move on to the next + ! block + call find_cell(p, found_cell) + if (.not. found_cell) cycle + + ! Determine if point is within desired cell + LEVEL_LOOP: do level = 1, p % n_coord + CELL_CHECK_LOOP: do i_cell = 1, size(this % cell_id) + if (cells(p % coord(level) % cell) % id == this % cell_id(i_cell)) then + + ! Determine what material this is + i_material = p % material + + ! Check if we've already had a hit in this material and if so, + ! simply add one + do j = 1, indices(i_cell) % size() + if (indices(i_cell) % data(j) == i_material) then + hits(i_cell) % data(j) = hits(i_cell) % data(j) + 1 + cycle CELL_CHECK_LOOP + end if + end do + + ! If we make it here, that means we haven't yet had a hit in this + ! material. Add an entry to both the indices list and the hits list + call indices(i_cell) % push_back(i_material) + call hits(i_cell) % push_back(1) + end if + end do CELL_CHECK_LOOP + end do LEVEL_LOOP + end do SAMPLE_LOOP + !$omp end do + + ! ========================================================================== + ! REDUCE HITS ONTO MASTER THREAD + + ! At this point, each thread has its own pair of index/hits lists and we now + ! need to reduce them. OpenMP is not nearly smart enough to do this on its + ! own, so we have to manually reduce them. + +#ifdef _OPENMP +!$omp do ordered schedule(static) + THREAD_LOOP: do i = 1, omp_get_num_threads() +!$omp ordered + do i_cell = 1, size(this % cell_id) + INDEX_LOOP: do j = 1, indices(i_cell) % size() + ! Check if this material has been added to the master list and if so, + ! accumulate the number of hits + do k = 1, master_indices(i_cell) % size() + if (indices(i_cell) % data(j) == master_indices(i_cell) % data(k)) then + master_hits(i_cell) % data(k) = & + master_hits(i_cell) % data(k) + hits(i_cell) % data(j) + cycle INDEX_LOOP + end if + end do + + ! If we made it here, this means the material hasn't yet been added to + ! the master list, so add an entry to both the master indices and master + ! hits lists + call master_indices(i_cell) % push_back(indices(i_cell) % data(j)) + call master_hits(i_cell) % push_back(hits(i_cell) % data(j)) + end do INDEX_LOOP + end do +!$omp end ordered + end do THREAD_LOOP +!$omp end do +#else + master_indices = indices + master_hits = hits +#endif + + call prn_set_stream(STREAM_TRACKING) +!$omp end parallel + + ! ========================================================================== + ! REDUCE HITS ONTO MASTER PROCESS + + volume_sample = product(this % upper_right - this % lower_left) + + do i_cell = 1, size(this % cell_id) + atoms(:, :) = ZERO + total_hits = 0 + + if (master) then +#ifdef MPI + do j = 1, n_procs - 1 + call MPI_RECV(n, 1, MPI_INTEGER, j, 0, MPI_COMM_WORLD, & + MPI_STATUS_IGNORE, mpi_err) + + allocate(data(2*n)) + call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, MPI_COMM_WORLD, & + MPI_STATUS_IGNORE, mpi_err) + do k = 0, n - 1 + do m = 1, master_indices(i_cell) % size() + if (data(2*k + 1) == master_indices(i_cell) % data(m)) then + master_hits(i_cell) % data(m) = master_hits(i_cell) % data(m) + & + data(2*k + 2) + end if + end do + end do + deallocate(data) + end do +#endif + + do j = 1, master_indices(i_cell) % size() + total_hits = total_hits + master_hits(i_cell) % data(j) + f = real(master_hits(i_cell) % data(j), 8) / this % samples + var_f = f*(ONE - f) / this % samples + + i_material = master_indices(i_cell) % data(j) + if (i_material == MATERIAL_VOID) cycle + + associate (mat => materials(i_material)) + do k = 1, size(mat % nuclide) + ! Accumulate nuclide density + i_nuclide = mat % nuclide(k) + atoms(1, i_nuclide) = atoms(1, i_nuclide) + & + mat % atom_density(k) * f + atoms(2, i_nuclide) = atoms(2, i_nuclide) + & + mat % atom_density(k)**2 * var_f + end do + end associate + end do + + ! Determine volume + volume(1, i_cell) = real(total_hits, 8) / this % samples * volume_sample + volume(2, i_cell) = sqrt(volume(1, i_cell) * (volume_sample - & + volume(1, i_cell)) / this % samples) + + ! Determine total number of atoms. At this point, we have values in + ! atoms/b-cm. To get to atoms we multiple by 10^24 V. + do j = 1, size(atoms, 2) + atoms(1, j) = 1.0e24_8 * volume_sample * atoms(1, j) + atoms(2, j) = 1.0e24_8 * volume_sample * sqrt(atoms(2, j)) + end do + + ! Convert full arrays to vectors + do j = 1, size(nuclides) + if (atoms(1, j) > ZERO) then + call nuclide_vec(i_cell) % push_back(j) + call atoms_vec(i_cell) % push_back(atoms(1, j)) + call uncertainty_vec(i_cell) % push_back(atoms(2, j)) + end if + end do + + else +#ifdef MPI + n = master_indices(i_cell) % size() + allocate(data(2*n)) + do k = 0, n - 1 + data(2*k + 1) = master_indices(i_cell) % data(k + 1) + data(2*k + 2) = master_hits(i_cell) % data(k + 1) + end do + + call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, MPI_COMM_WORLD, mpi_err) + call MPI_SEND(data, 2*n, MPI_INTEGER, 0, 1, MPI_COMM_WORLD, mpi_err) + deallocate(data) +#endif + end if + end do + + end subroutine get_volume + +!=============================================================================== +! WRITE_VOLUME writes the results of a single stochastic volume calculation to +! an HDF5 file +!=============================================================================== + + subroutine write_volume(this, filename, volume, nuclide_vec, atoms_vec, & + uncertainty_vec) + type(VolumeCalculation), intent(in) :: this + character(*), intent(in) :: filename ! filename for HDF5 file + real(8), intent(in) :: volume(:,:) ! volume mean/stdev in each cell + type(VectorInt), intent(in) :: nuclide_vec(:) ! indices in nuclides array + type(VectorReal), intent(in) :: atoms_vec(:) ! total # of atoms of each nuclide + type(VectorReal), intent(in) :: uncertainty_vec(:) ! uncertainty of total # of atoms + + integer :: i, j + integer :: n + integer(HID_T) :: file_id + integer(HID_T) :: group_id + real(8), allocatable :: atom_data(:,:) ! mean/stdev of total # of atoms for + ! each nuclide + character(MAX_WORD_LEN), allocatable :: nucnames(:) ! names of nuclides + + ! Create HDF5 file + file_id = file_create(filename) + + ! Write basic metadata + call write_attribute(file_id, "samples", this % samples) + call write_attribute(file_id, "lower_left", this % lower_left) + call write_attribute(file_id, "upper_right", this % upper_right) + + do i = 1, size(this % cell_id) + group_id = create_group(file_id, "cell_" // trim(to_str(& + this % cell_id(i)))) + + ! Write volume for cell + call write_dataset(group_id, "volume", volume(:, i)) + + ! Create array of nuclide names from the vector + n = nuclide_vec(i) % size() + if (n > 0) then + allocate(nucnames(n)) + do j = 1, n + nucnames(j) = nuclides(nuclide_vec(i) % data(j)) % name + end do + + ! Create array of total # of atoms with uncertainty for each nuclide + allocate(atom_data(2, n)) + atom_data(1, :) = atoms_vec(i) % data(1:n) + atom_data(2, :) = uncertainty_vec(i) % data(1:n) + + ! Write results + call write_dataset(group_id, "nuclides", nucnames) + call write_dataset(group_id, "atoms", atom_data) + + deallocate(nucnames) + deallocate(atom_data) + end if + + call close_group(group_id) + end do + call file_close(file_id) + end subroutine write_volume + +end module volume_calc diff --git a/src/volume_header.F90 b/src/volume_header.F90 new file mode 100644 index 0000000000..3a34f1231a --- /dev/null +++ b/src/volume_header.F90 @@ -0,0 +1,42 @@ +module volume_header + + use error, only: fatal_error + use xml_interface + + implicit none + + type VolumeCalculation + integer, allocatable :: cell_id(:) + real(8) :: lower_left(3) + real(8) :: upper_right(3) + integer :: samples + contains + procedure :: from_xml => volume_from_xml + end type VolumeCalculation + +contains + + subroutine volume_from_xml(this, node_vol) + class(VolumeCalculation), intent(out) :: this + type(Node), pointer :: node_vol + + integer :: num_cells + + ! Read cell IDs + if (check_for_node(node_vol, "cells")) then + num_cells = get_arraysize_integer(node_vol, "cells") + else + call fatal_error("Must specify at least one cell for a volume calculation") + end if + allocate(this % cell_id(num_cells)) + call get_node_array(node_vol, "cells", this % cell_id) + + ! Read lower-left and upper-right bounding coordinates + call get_node_array(node_vol, "lower_left", this % lower_left) + call get_node_array(node_vol, "upper_right", this % upper_right) + + ! Read number of samples + call get_node_value(node_vol, "samples", this % samples) + end subroutine volume_from_xml + +end module volume_header From cac6ac6f665a3c3d2d72bb1bd87c813064a6b87a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 28 Jul 2016 12:35:15 -0500 Subject: [PATCH 02/12] Add VolumeCalculation object. Update settings.xml RELAX NG schema --- openmc/__init__.py | 5 +- openmc/settings.py | 204 +++++++++++++++++++++------------------ openmc/volume.py | 204 +++++++++++++++++++++++++++++++++++++++ src/relaxng/settings.rnc | 11 +++ src/relaxng/settings.rng | 62 ++++++++++++ 5 files changed, 392 insertions(+), 94 deletions(-) create mode 100644 openmc/volume.py diff --git a/openmc/__init__.py b/openmc/__init__.py index 557e13039f..026ccce114 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -6,6 +6,9 @@ from openmc.nuclide import * from openmc.macroscopic import * from openmc.material import * from openmc.plots import * +from openmc.region import * +from openmc.volume import * +from openmc.source import * from openmc.settings import * from openmc.surface import * from openmc.universe import * @@ -18,8 +21,6 @@ from openmc.cmfd import * from openmc.executor import * from openmc.statepoint import * from openmc.summary import * -from openmc.region import * -from openmc.source import * from openmc.particle_restart import * try: diff --git a/openmc/settings.py b/openmc/settings.py index b9a93bd114..918f97c7da 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1,4 +1,4 @@ -from collections import Iterable +from collections import Iterable, MutableSequence from numbers import Real, Integral import warnings from xml.etree import ElementTree as ET @@ -7,10 +7,8 @@ import sys import numpy as np from openmc.clean_xml import clean_xml_indentation -from openmc.checkvalue import (check_type, check_length, check_value, - check_greater_than, check_less_than) -from openmc import Nuclide -from openmc.source import Source +import openmc.checkvalue as cv +from openmc import Nuclide, VolumeCalculation, Source if sys.version_info[0] >= 3: basestring = str @@ -137,6 +135,8 @@ class Settings(object): resonance cross sections. resonance_scattering : ResonanceScattering or iterable of ResonanceScattering The elastic scattering model to use for resonant isotopes + volume_calculations : VolumeCalculation or iterable of VolumeCalculation + Stochastic volume calculation specifications """ @@ -220,6 +220,8 @@ class Settings(object): self._multipole_active = None self._resonance_scattering = None + self._volume_calculations = cv.CheckedList(VolumeCalculation, + 'volume calculations') @property def run_mode(self): @@ -421,6 +423,10 @@ class Settings(object): def resonance_scattering(self): return self._resonance_scattering + @property + def volume_calculations(self): + return self._volume_calculations + @run_mode.setter def run_mode(self, run_mode): if run_mode not in ['eigenvalue', 'fixed source']: @@ -431,26 +437,26 @@ class Settings(object): @batches.setter def batches(self, batches): - check_type('batches', batches, Integral) - check_greater_than('batches', batches, 0) + cv.check_type('batches', batches, Integral) + cv.check_greater_than('batches', batches, 0) self._batches = batches @generations_per_batch.setter def generations_per_batch(self, generations_per_batch): - check_type('generations per patch', generations_per_batch, Integral) - check_greater_than('generations per batch', generations_per_batch, 0) + cv.check_type('generations per patch', generations_per_batch, Integral) + cv.check_greater_than('generations per batch', generations_per_batch, 0) self._generations_per_batch = generations_per_batch @inactive.setter def inactive(self, inactive): - check_type('inactive batches', inactive, Integral) - check_greater_than('inactive batches', inactive, 0, True) + cv.check_type('inactive batches', inactive, Integral) + cv.check_greater_than('inactive batches', inactive, 0, True) self._inactive = inactive @particles.setter def particles(self, particles): - check_type('particles', particles, Integral) - check_greater_than('particles', particles, 0) + cv.check_type('particles', particles, Integral) + cv.check_greater_than('particles', particles, 0) self._particles = particles @keff_trigger.setter @@ -484,14 +490,14 @@ class Settings(object): @energy_mode.setter def energy_mode(self, energy_mode): - check_value('energy mode', energy_mode, + cv.check_value('energy mode', energy_mode, ['continuous-energy', 'multi-group']) self._energy_mode = energy_mode @max_order.setter def max_order(self, max_order): - check_type('maximum scattering order', max_order, Integral) - check_greater_than('maximum scattering order', max_order, 0, True) + cv.check_type('maximum scattering order', max_order, Integral) + cv.check_greater_than('maximum scattering order', max_order, 0, True) self._max_order = max_order @source.setter @@ -499,7 +505,7 @@ class Settings(object): if isinstance(source, Source): self._source = [source,] else: - check_type('source distribution', source, Iterable, Source) + cv.check_type('source distribution', source, Iterable, Source) self._source = source @output.setter @@ -525,197 +531,197 @@ class Settings(object): @output_path.setter def output_path(self, output_path): - check_type('output path', output_path, basestring) + cv.check_type('output path', output_path, basestring) self._output_path = output_path @verbosity.setter def verbosity(self, verbosity): - check_type('verbosity', verbosity, Integral) - check_greater_than('verbosity', verbosity, 1, True) - check_less_than('verbosity', verbosity, 10, True) + cv.check_type('verbosity', verbosity, Integral) + cv.check_greater_than('verbosity', verbosity, 1, True) + cv.check_less_than('verbosity', verbosity, 10, True) self._verbosity = verbosity @statepoint_batches.setter def statepoint_batches(self, batches): - check_type('statepoint batches', batches, Iterable, Integral) + cv.check_type('statepoint batches', batches, Iterable, Integral) for batch in batches: - check_greater_than('statepoint batch', batch, 0) + cv.check_greater_than('statepoint batch', batch, 0) self._statepoint_batches = batches @statepoint_interval.setter def statepoint_interval(self, interval): - check_type('statepoint interval', interval, Integral) + cv.check_type('statepoint interval', interval, Integral) self._statepoint_interval = interval @sourcepoint_batches.setter def sourcepoint_batches(self, batches): - check_type('sourcepoint batches', batches, Iterable, Integral) + cv.check_type('sourcepoint batches', batches, Iterable, Integral) for batch in batches: - check_greater_than('sourcepoint batch', batch, 0) + cv.check_greater_than('sourcepoint batch', batch, 0) self._sourcepoint_batches = batches @sourcepoint_interval.setter def sourcepoint_interval(self, interval): - check_type('sourcepoint interval', interval, Integral) + cv.check_type('sourcepoint interval', interval, Integral) self._sourcepoint_interval = interval @sourcepoint_separate.setter def sourcepoint_separate(self, source_separate): - check_type('sourcepoint separate', source_separate, bool) + cv.check_type('sourcepoint separate', source_separate, bool) self._sourcepoint_separate = source_separate @sourcepoint_write.setter def sourcepoint_write(self, source_write): - check_type('sourcepoint write', source_write, bool) + cv.check_type('sourcepoint write', source_write, bool) self._sourcepoint_write = source_write @sourcepoint_overwrite.setter def sourcepoint_overwrite(self, source_overwrite): - check_type('sourcepoint overwrite', source_overwrite, bool) + cv.check_type('sourcepoint overwrite', source_overwrite, bool) self._sourcepoint_overwrite = source_overwrite @confidence_intervals.setter def confidence_intervals(self, confidence_intervals): - check_type('confidence interval', confidence_intervals, bool) + cv.check_type('confidence interval', confidence_intervals, bool) self._confidence_intervals = confidence_intervals @cross_sections.setter def cross_sections(self, cross_sections): - check_type('cross sections', cross_sections, basestring) + cv.check_type('cross sections', cross_sections, basestring) self._cross_sections = cross_sections @multipole_library.setter def multipole_library(self, multipole_library): - check_type('cross sections', multipole_library, basestring) + cv.check_type('cross sections', multipole_library, basestring) self._multipole_library = multipole_library @energy_grid.setter def energy_grid(self, energy_grid): - check_value('energy grid', energy_grid, + cv.check_value('energy grid', energy_grid, ['nuclide', 'logarithm', 'material-union']) self._energy_grid = energy_grid @ptables.setter def ptables(self, ptables): - check_type('probability tables', ptables, bool) + cv.check_type('probability tables', ptables, bool) self._ptables = ptables @run_cmfd.setter def run_cmfd(self, run_cmfd): - check_type('run_cmfd', run_cmfd, bool) + cv.check_type('run_cmfd', run_cmfd, bool) self._run_cmfd = run_cmfd @seed.setter def seed(self, seed): - check_type('random number generator seed', seed, Integral) - check_greater_than('random number generator seed', seed, 0) + cv.check_type('random number generator seed', seed, Integral) + cv.check_greater_than('random number generator seed', seed, 0) self._seed = seed @survival_biasing.setter def survival_biasing(self, survival_biasing): - check_type('survival biasing', survival_biasing, bool) + cv.check_type('survival biasing', survival_biasing, bool) self._survival_biasing = survival_biasing @weight.setter def weight(self, weight): - check_type('weight cutoff', weight, Real) - check_greater_than('weight cutoff', weight, 0.0) + cv.check_type('weight cutoff', weight, Real) + cv.check_greater_than('weight cutoff', weight, 0.0) self._weight = weight @weight_avg.setter def weight_avg(self, weight_avg): - check_type('average survival weight', weight_avg, Real) - check_greater_than('average survival weight', weight_avg, 0.0) + cv.check_type('average survival weight', weight_avg, Real) + cv.check_greater_than('average survival weight', weight_avg, 0.0) self._weight_avg = weight_avg @entropy_dimension.setter def entropy_dimension(self, dimension): - check_type('entropy mesh dimension', dimension, Iterable, Integral) - check_length('entropy mesh dimension', dimension, 3) + cv.check_type('entropy mesh dimension', dimension, Iterable, Integral) + cv.check_length('entropy mesh dimension', dimension, 3) self._entropy_dimension = dimension @entropy_lower_left.setter def entropy_lower_left(self, lower_left): - check_type('entropy mesh lower left corner', lower_left, + cv.check_type('entropy mesh lower left corner', lower_left, Iterable, Real) - check_length('entropy mesh lower left corner', lower_left, 3) + cv.check_length('entropy mesh lower left corner', lower_left, 3) self._entropy_lower_left = lower_left @entropy_upper_right.setter def entropy_upper_right(self, upper_right): - check_type('entropy mesh upper right corner', upper_right, + cv.check_type('entropy mesh upper right corner', upper_right, Iterable, Real) - check_length('entropy mesh upper right corner', upper_right, 3) + cv.check_length('entropy mesh upper right corner', upper_right, 3) self._entropy_upper_right = upper_right @trigger_active.setter def trigger_active(self, trigger_active): - check_type('trigger active', trigger_active, bool) + cv.check_type('trigger active', trigger_active, bool) self._trigger_active = trigger_active @trigger_max_batches.setter def trigger_max_batches(self, trigger_max_batches): - check_type('trigger maximum batches', trigger_max_batches, Integral) - check_greater_than('trigger maximum batches', trigger_max_batches, 0) + cv.check_type('trigger maximum batches', trigger_max_batches, Integral) + cv.check_greater_than('trigger maximum batches', trigger_max_batches, 0) self._trigger_max_batches = trigger_max_batches @trigger_batch_interval.setter def trigger_batch_interval(self, trigger_batch_interval): - check_type('trigger batch interval', trigger_batch_interval, Integral) - check_greater_than('trigger batch interval', trigger_batch_interval, 0) + cv.check_type('trigger batch interval', trigger_batch_interval, Integral) + cv.check_greater_than('trigger batch interval', trigger_batch_interval, 0) self._trigger_batch_interval = trigger_batch_interval @no_reduce.setter def no_reduce(self, no_reduce): - check_type('no reduction option', no_reduce, bool) + cv.check_type('no reduction option', no_reduce, bool) self._no_reduce = no_reduce @threads.setter def threads(self, threads): - check_type('number of threads', threads, Integral) - check_greater_than('number of threads', threads, 0) + cv.check_type('number of threads', threads, Integral) + cv.check_greater_than('number of threads', threads, 0) self._threads = threads @trace.setter def trace(self, trace): - check_type('trace', trace, Iterable, Integral) - check_length('trace', trace, 3) - check_greater_than('trace batch', trace[0], 0) - check_greater_than('trace generation', trace[1], 0) - check_greater_than('trace particle', trace[2], 0) + cv.check_type('trace', trace, Iterable, Integral) + cv.check_length('trace', trace, 3) + cv.check_greater_than('trace batch', trace[0], 0) + cv.check_greater_than('trace generation', trace[1], 0) + cv.check_greater_than('trace particle', trace[2], 0) self._trace = trace @track.setter def track(self, track): - check_type('track', track, Iterable, Integral) + cv.check_type('track', track, Iterable, Integral) if len(track) % 3 != 0: msg = 'Unable to set the track to "{0}" since its length is ' \ 'not a multiple of 3'.format(track) raise ValueError(msg) for t in zip(track[::3], track[1::3], track[2::3]): - check_greater_than('track batch', t[0], 0) - check_greater_than('track generation', t[0], 0) - check_greater_than('track particle', t[0], 0) + cv.check_greater_than('track batch', t[0], 0) + cv.check_greater_than('track generation', t[0], 0) + cv.check_greater_than('track particle', t[0], 0) self._track = track @ufs_dimension.setter def ufs_dimension(self, dimension): - check_type('UFS mesh dimension', dimension, Iterable, Integral) - check_length('UFS mesh dimension', dimension, 3) + cv.check_type('UFS mesh dimension', dimension, Iterable, Integral) + cv.check_length('UFS mesh dimension', dimension, 3) for dim in dimension: - check_greater_than('UFS mesh dimension', dim, 1, True) + cv.check_greater_than('UFS mesh dimension', dim, 1, True) self._ufs_dimension = dimension @ufs_lower_left.setter def ufs_lower_left(self, lower_left): - check_type('UFS mesh lower left corner', lower_left, Iterable, Real) - check_length('UFS mesh lower left corner', lower_left, 3) + cv.check_type('UFS mesh lower left corner', lower_left, Iterable, Real) + cv.check_length('UFS mesh lower left corner', lower_left, 3) self._ufs_lower_left = lower_left @ufs_upper_right.setter def ufs_upper_right(self, upper_right): - check_type('UFS mesh upper right corner', upper_right, Iterable, Real) - check_length('UFS mesh upper right corner', upper_right, 3) + cv.check_type('UFS mesh upper right corner', upper_right, Iterable, Real) + cv.check_length('UFS mesh upper right corner', upper_right, 3) self._ufs_upper_right = upper_right @dd_mesh_dimension.setter @@ -724,8 +730,8 @@ class Settings(object): warnings.warn('This feature is not yet implemented in a release ' 'version of openmc') - check_type('DD mesh dimension', dimension, Iterable, Integral) - check_length('DD mesh dimension', dimension, 3) + cv.check_type('DD mesh dimension', dimension, Iterable, Integral) + cv.check_length('DD mesh dimension', dimension, 3) self._dd_mesh_dimension = dimension @@ -735,8 +741,8 @@ class Settings(object): warnings.warn('This feature is not yet implemented in a release ' 'version of openmc') - check_type('DD mesh lower left corner', lower_left, Iterable, Real) - check_length('DD mesh lower left corner', lower_left, 3) + cv.check_type('DD mesh lower left corner', lower_left, Iterable, Real) + cv.check_length('DD mesh lower left corner', lower_left, 3) self._dd_mesh_lower_left = lower_left @@ -746,8 +752,8 @@ class Settings(object): warnings.warn('This feature is not yet implemented in a release ' 'version of openmc') - check_type('DD mesh upper right corner', upper_right, Iterable, Real) - check_length('DD mesh upper right corner', upper_right, 3) + cv.check_type('DD mesh upper right corner', upper_right, Iterable, Real) + cv.check_length('DD mesh upper right corner', upper_right, 3) self._dd_mesh_upper_right = upper_right @@ -757,7 +763,7 @@ class Settings(object): warnings.warn('This feature is not yet implemented in a release ' 'version of openmc') - check_type('DD nodemap', nodemap, Iterable) + cv.check_type('DD nodemap', nodemap, Iterable) nodemap = np.array(nodemap).flatten() @@ -782,7 +788,7 @@ class Settings(object): warnings.warn('This feature is not yet implemented in a release ' 'version of openmc') - check_type('DD allow leakage', allow, bool) + cv.check_type('DD allow leakage', allow, bool) self._dd_allow_leakage = allow @@ -793,25 +799,34 @@ class Settings(object): warnings.warn('This feature is not yet implemented in a release ' 'version of openmc') - check_type('DD count interactions', interactions, bool) + cv.check_type('DD count interactions', interactions, bool) self._dd_count_interactions = interactions @use_windowed_multipole.setter def use_windowed_multipole(self, active): - check_type('use_windowed_multipole', active, bool) + cv.check_type('use_windowed_multipole', active, bool) self._multipole_active = active @resonance_scattering.setter def resonance_scattering(self, res): if isinstance(res, Iterable): - check_type('resonance_scattering', res, Iterable, + cv.check_type('resonance_scattering', res, Iterable, ResonanceScattering) self._resonance_scattering = res else: - check_type('resonance_scattering', res, ResonanceScattering) + cv.check_type('resonance_scattering', res, ResonanceScattering) self._resonance_scattering = [res] + @volume_calculations.setter + def volume_calculations(self, vol_calcs): + name = 'stochastic volume calculations' + if not isinstance(vol_calcs, MutableSequence): + vol_calcs = [vol_calcs] + cv.check_type(name, vol_calcs, MutableSequence) + self._volume_calculations = cv.CheckedList(VolumeCalculation, + name, vol_calcs) + def _create_run_mode_subelement(self): if self.run_mode == 'eigenvalue': @@ -873,6 +888,10 @@ class Settings(object): for source in self.source: self._settings_file.append(source.to_xml()) + def _create_volume_calcs_subelement(self): + for calc in self.volume_calculations: + self._settings_file.append(calc.to_xml()) + def _create_output_subelement(self): if self._output is not None: element = ET.SubElement(self._settings_file, "output") @@ -1154,6 +1173,7 @@ class Settings(object): self._create_dd_subelement() self._create_use_multipole_subelement() self._create_resonance_scattering_element() + self._create_volume_calcs_subelement() # Clean the indentation in the file to be user-readable clean_xml_indentation(self._settings_file) @@ -1216,29 +1236,29 @@ class ResonanceScattering(object): @nuclide.setter def nuclide(self, nuc): - check_type('nuclide', nuc, Nuclide) + cv.check_type('nuclide', nuc, Nuclide) self._nuclide = nuc @nuclide_0K.setter def nuclide_0K(self, nuc): - check_type('nuclide_0K', nuc, Nuclide) + cv.check_type('nuclide_0K', nuc, Nuclide) self._nuclide_0K = nuc @method.setter def method(self, m): - check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM')) + cv.check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM')) self._method = m @E_min.setter def E_min(self, E): - check_type('E_min', E, Real) - check_greater_than('E_min', E, 0, True) + cv.check_type('E_min', E, Real) + cv.check_greater_than('E_min', E, 0, True) self._E_min = E @E_max.setter def E_max(self, E): - check_type('E_max', E, Real) - check_greater_than('E_max', E, 0, True) + cv.check_type('E_max', E, Real) + cv.check_greater_than('E_max', E, 0, True) self._E_max = E def create_xml_subelement(self, xml_element): diff --git a/openmc/volume.py b/openmc/volume.py new file mode 100644 index 0000000000..511af62ed3 --- /dev/null +++ b/openmc/volume.py @@ -0,0 +1,204 @@ +from collections import Iterable, Mapping +from numbers import Real, Integral +from xml.etree import ElementTree as ET + +import numpy as np +import pandas as pd + +from openmc import Cell, Union +import openmc.checkvalue as cv + + +class VolumeCalculation(object): + """Stochastic volume calculation specifications and results. + + Parameters + ---------- + cells : Iterable of Cell + Cells to find volumes of + samples : int + Number of samples used to generate volume estimates + lower_left : Iterable of float + Lower-left coordinates of bounding box used to sample points. If this + argument is not supplied, an attempt is made to automatically determine + a bounding box. + upper_right : Iterable of float + Upper-right coordinates of bounding box used to sample points. If this + argument is not supplied, an attempt is made to automatically determine + a bounding box. + + Attributes + ---------- + cell_ids : Iterable of int + IDs of cells to find volumes of + samples : int + Number of samples used to generate volume estimates + lower_left : Iterable of float + Lower-left coordinates of bounding box used to sample points + upper_right : Iterable of float + Upper-right coordinates of bounding box used to sample points + results : dict + Dictionary whose keys are unique IDs of cells and values are + dictionaries with calculated volumes and total number of atoms for each + nuclide present in the cell. + volumes : dict + Dictionary whose keys are unique IDs of cells and values are the + estimated volumes + atoms_dataframe : pandas.DataFrame + DataFrame showing the estimated number of atoms for each nuclide present + in each cell specified. + + """ + def __init__(self, cells, samples, lower_left=None, + upper_right=None): + self._results = None + + cv.check_type('cells', cells, Iterable, Cell) + self.cell_ids = [c.id for c in cells] + self.samples = samples + + if lower_left is not None: + self.lower_left = lower_left + if upper_right is None: + raise ValueError('Both lower-left and upper-right coordinates ' + 'should be specified') + self.upper_right = upper_right + else: + ll, ur = Union(*[c.region for c in cells]).bounding_box + if np.any(np.isinf(ll)) or np.any(np.isinf(ur)): + raise ValueError('Could not automatically determine bounding box ' + 'for stochastic volume calculation.') + else: + self.lower_left = ll + self.upper_right = ur + + @property + def cell_ids(self): + return self._cell_ids + + @property + def samples(self): + return self._samples + + @property + def lower_left(self): + return self._lower_left + + @property + def upper_right(self): + return self._upper_right + + @property + def results(self): + return self._results + + @property + def volumes(self): + return {uid: results['volume'] for uid, results in self.results.items()} + + @property + def atoms_dataframe(self): + items = [] + columns = ['Cell', 'Nuclide', 'Atoms', 'Uncertainty'] + for cell_id, results in self.results.items(): + for name, atoms in results['atoms']: + items.append((cell_id, name, atoms[0], atoms[1])) + + return pd.DataFrame.from_records(items, columns=columns) + + @cell_ids.setter + def cell_ids(self, cell_ids): + cv.check_type('cell IDs', cell_ids, Iterable, Real) + self._cell_ids = cell_ids + + @samples.setter + def samples(self, samples): + cv.check_type('number of samples', samples, Integral) + cv.check_greater_than('number of samples', samples, 0) + self._samples = samples + + @lower_left.setter + def lower_left(self, lower_left): + name = 'lower-left bounding box coordinates', + cv.check_type(name, lower_left, Iterable, Real) + cv.check_length(name, lower_left, 3) + self._lower_left = lower_left + + @upper_right.setter + def upper_right(self, upper_right): + name = 'upper-right bounding box coordinates' + cv.check_type(name, upper_right, Iterable, Real) + cv.check_length(name, upper_right, 3) + self._upper_right = upper_right + + @results.setter + def results(self, results): + cv.check_type('results', results, Mapping) + self._results = results + + @classmethod + def from_hdf5(cls, filename): + """Load stochastic volume calculation results from HDF5 file. + + Parameters + ---------- + filename : str + Path to volume.h5 file + + Returns + ------- + openmc.VolumeCalculation + Results of the stochastic volume calculation + + """ + import h5py + + with h5py.File(filename, 'r') as f: + samples = f.attrs['samples'] + lower_left = f.attrs['lower_left'] + upper_right = f.attrs['upper_right'] + + results = {} + cell_ids = [] + for obj_name in f: + if obj_name.startswith('cell_'): + cell_id = int(obj_name[5:]) + cell_ids.append(cell_id) + group = f[obj_name] + volume = tuple(group['volume'].value) + nucnames = group['nuclides'].value + atoms = group['atoms'].value + + atom_list = [] + for name_i, atoms_i in zip(nucnames, atoms): + atom_list.append((name_i.decode(), tuple(atoms_i))) + results[cell_id] = {'volume': volume, 'atoms': atom_list} + + # Instantiate some throw-away cells that are used by the constructor to + # assign IDs + cells = [Cell(uid) for uid in cell_ids] + + # Instantiate the class and assign results + vol = cls(cells, samples, lower_left, upper_right) + vol.results = results + return vol + + def to_xml(self): + """Return XML representation of the volume calculation + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing volume calculation data + + """ + element = ET.Element("volume_calc") + cell_elem = ET.SubElement(element, "cells") + cell_elem.text = ' '.join(str(uid) for uid in self.cell_ids) + samples_elem = ET.SubElement(element, "samples") + samples_elem.text = str(self.samples) + ll_elem = ET.SubElement(element, "lower_left") + ll_elem.text = ' '.join(str(x) for x in self.lower_left) + ur_elem = ET.SubElement(element, "upper_right") + ur_elem.text = ' '.join(str(x) for x in self.upper_right) + return element diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 46950c63fb..78a32171b9 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -142,6 +142,17 @@ element settings { element verbosity { xsd:positiveInteger }? & + element volume_calc { + (element cells { list { xsd:positiveInteger+ } } | + attribute cells { list { xsd:positiveInteger+ } }) & + (element samples { xsd:positiveInteger } | + attribute samples { xsd:positiveInteger }) & + (element lower_left { list { xsd:double+ } } | + attribute lower_left { list { xsd:double+ } }) & + (element upper_right { list { xsd:double+ } } | + attribute upper_right { list { xsd:double+ } }) + }+ & + element uniform_fs{ (element dimension { list { xsd:positiveInteger+ } } | attribute dimension { list { xsd:positiveInteger+ } }) & diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 0b50f79699..64fe42239d 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -625,6 +625,68 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 66d772da7997008832711db8ecfe8b58989a8ea2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 28 Jul 2016 13:52:42 -0500 Subject: [PATCH 03/12] Refactor settings to use more CheckedLists. Rename all to_xml methods to_xml_element --- openmc/settings.py | 79 +++++++++++++++++------------------ openmc/source.py | 8 ++-- openmc/stats/multivariate.py | 74 ++++++++++++++++++++++++++------ openmc/stats/univariate.py | 81 ++++++++++++++++++++++++++++++++---- openmc/volume.py | 2 +- 5 files changed, 177 insertions(+), 67 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 918f97c7da..eb56381b07 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -155,7 +155,7 @@ class Settings(object): self._max_order = None # Source subelement - self._source = None + self._source = cv.CheckedList(Source, 'source distributions') self._confidence_intervals = None self._cross_sections = None @@ -216,12 +216,12 @@ class Settings(object): self._settings_file = ET.Element("settings") self._run_mode_subelement = None - self._source_element = None self._multipole_active = None - self._resonance_scattering = None - self._volume_calculations = cv.CheckedList(VolumeCalculation, - 'volume calculations') + self._resonance_scattering = cv.CheckedList( + ResonanceScattering, 'resonance scattering models') + self._volume_calculations = cv.CheckedList( + VolumeCalculation, 'volume calculations') @property def run_mode(self): @@ -502,11 +502,9 @@ class Settings(object): @source.setter def source(self, source): - if isinstance(source, Source): - self._source = [source,] - else: - cv.check_type('source distribution', source, Iterable, Source) - self._source = source + if not isinstance(source, MutableSequence): + source = [source] + self._source = cv.CheckedList(Source, 'source distributions', source) @output.setter def output(self, output): @@ -810,22 +808,17 @@ class Settings(object): @resonance_scattering.setter def resonance_scattering(self, res): - if isinstance(res, Iterable): - cv.check_type('resonance_scattering', res, Iterable, - ResonanceScattering) - self._resonance_scattering = res - else: - cv.check_type('resonance_scattering', res, ResonanceScattering) - self._resonance_scattering = [res] + if not isinstance(res, MutableSequence): + res = [res] + self._resonance_scattering = cv.CheckedList( + ResonanceScattering, 'resonance scattering models', res) @volume_calculations.setter def volume_calculations(self, vol_calcs): - name = 'stochastic volume calculations' if not isinstance(vol_calcs, MutableSequence): vol_calcs = [vol_calcs] - cv.check_type(name, vol_calcs, MutableSequence) - self._volume_calculations = cv.CheckedList(VolumeCalculation, - name, vol_calcs) + self._volume_calculations = cv.CheckedList( + VolumeCalculation, 'stochastic volume calculations', vol_calcs) def _create_run_mode_subelement(self): @@ -884,13 +877,12 @@ class Settings(object): element.text = str(self._max_order) def _create_source_subelement(self): - if self.source is not None: - for source in self.source: - self._settings_file.append(source.to_xml()) + for source in self.source: + self._settings_file.append(source.to_xml_element()) def _create_volume_calcs_subelement(self): for calc in self.volume_calculations: - self._settings_file.append(calc.to_xml()) + self._settings_file.append(calc.to_xml_element()) def _create_output_subelement(self): if self._output is not None: @@ -1121,18 +1113,15 @@ class Settings(object): "use_windowed_multipole") element.text = str(self._multipole_active) - def _create_resonance_scattering_element(self): - if self.resonance_scattering is None: - return - - element = ET.SubElement(self._settings_file, "resonance_scattering") - - for r in self.resonance_scattering: - if r.nuclide.name != r.nuclide_0K.name: - raise ValueError("The nuclide and nuclide_0K attributes of " - "a ResonantScattering object must have " - "identical names.") - r.create_xml_subelement(element) + def _create_resonance_scattering_subelement(self): + if len(self.resonance_scattering) > 0: + elem = ET.SubElement(self._settings_file, 'resonance_scattering') + for r in self.resonance_scattering: + if r.nuclide.name != r.nuclide_0K.name: + raise ValueError("The nuclide and nuclide_0K attributes of " + "a ResonantScattering object must have " + "identical names.") + elem.append(r.to_xml_element()) def export_to_xml(self): """Create a settings.xml file that can be used for a simulation. @@ -1144,7 +1133,6 @@ class Settings(object): self._source_subelement = None self._trigger_subelement = None self._run_mode_subelement = None - self._source_element = None self._create_run_mode_subelement() self._create_source_subelement() @@ -1172,7 +1160,7 @@ class Settings(object): self._create_ufs_subelement() self._create_dd_subelement() self._create_use_multipole_subelement() - self._create_resonance_scattering_element() + self._create_resonance_scattering_subelement() self._create_volume_calcs_subelement() # Clean the indentation in the file to be user-readable @@ -1261,8 +1249,16 @@ class ResonanceScattering(object): cv.check_greater_than('E_max', E, 0, True) self._E_max = E - def create_xml_subelement(self, xml_element): - scatterer = ET.SubElement(xml_element, "scatterer") + def to_xml_element(self): + """Return XML representation of the resonance scattering model + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing resonance scattering model + + """ + scatterer = ET.Element("scatterer") subelement = ET.SubElement(scatterer, 'nuclide') subelement.text = self.nuclide.name if self.method is not None: @@ -1278,3 +1274,4 @@ class ResonanceScattering(object): if self.E_max is not None: subelement = ET.SubElement(scatterer, 'E_max') subelement.text = str(self.E_max) + return scatterer diff --git a/openmc/source.py b/openmc/source.py index 7e8a68accf..ee32cd0f4b 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -103,7 +103,7 @@ class Source(object): cv.check_greater_than('source strength', strength, 0.0, True) self._strength = strength - def to_xml(self): + def to_xml_element(self): """Return XML representation of the source Returns @@ -117,9 +117,9 @@ class Source(object): if self.file is not None: element.set("file", self.file) if self.space is not None: - element.append(self.space.to_xml()) + element.append(self.space.to_xml_element()) if self.angle is not None: - element.append(self.angle.to_xml()) + element.append(self.angle.to_xml_element()) if self.energy is not None: - element.append(self.energy.to_xml('energy')) + element.append(self.energy.to_xml_element('energy')) return element diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index e4eadd7aa4..e49a94ee19 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -50,7 +50,7 @@ class UnitSphere(object): self._reference_uvw = uvw/np.linalg.norm(uvw) @abstractmethod - def to_xml(self): + def to_xml_element(self): return '' @@ -109,13 +109,21 @@ class PolarAzimuthal(UnitSphere): cv.check_type('azimuthal angle', phi, Univariate) self._phi = phi - def to_xml(self): + def to_xml_element(self): + """Return XML representation of the angular distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing angular distribution data + + """ element = ET.Element('angle') element.set("type", "mu-phi") if self.reference_uvw is not None: element.set("reference_uvw", ' '.join(map(str, self.reference_uvw))) - element.append(self.mu.to_xml('mu')) - element.append(self.phi.to_xml('phi')) + element.append(self.mu.to_xml_element('mu')) + element.append(self.phi.to_xml_element('phi')) return element @@ -127,7 +135,15 @@ class Isotropic(UnitSphere): def __init__(self): super(Isotropic, self).__init__() - def to_xml(self): + def to_xml_element(self): + """Return XML representation of the isotropic distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing isotropic distribution data + + """ element = ET.Element('angle') element.set("type", "isotropic") return element @@ -152,7 +168,15 @@ class Monodirectional(UnitSphere): def __init__(self, reference_uvw=[1., 0., 0.]): super(Monodirectional, self).__init__(reference_uvw) - def to_xml(self): + def to_xml_element(self): + """Return XML representation of the monodirectional distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing monodirectional distribution data + + """ element = ET.Element('angle') element.set("type", "monodirectional") if self.reference_uvw is not None: @@ -174,7 +198,7 @@ class Spatial(object): pass @abstractmethod - def to_xml(self): + def to_xml_element(self): return '' @@ -238,12 +262,20 @@ class CartesianIndependent(Spatial): cv.check_type('z coordinate', z, Univariate) self._z = z - def to_xml(self): + def to_xml_element(self): + """Return XML representation of the spatial distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing spatial distribution data + + """ element = ET.Element('space') element.set('type', 'cartesian') - element.append(self.x.to_xml('x')) - element.append(self.y.to_xml('y')) - element.append(self.z.to_xml('z')) + element.append(self.x.to_xml_element('x')) + element.append(self.y.to_xml_element('y')) + element.append(self.z.to_xml_element('z')) return element @@ -308,7 +340,15 @@ class Box(Spatial): cv.check_type('only fissionable', only_fissionable, bool) self._only_fissionable = only_fissionable - def to_xml(self): + def to_xml_element(self): + """Return XML representation of the box distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing box distribution data + + """ element = ET.Element('space') if self.only_fissionable: element.set("type", "fission") @@ -352,7 +392,15 @@ class Point(Spatial): cv.check_length('coordinate', xyz, 3) self._xyz = xyz - def to_xml(self): + def to_xml_element(self): + """Return XML representation of the point distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing point distribution location + + """ element = ET.Element('space') element.set("type", "point") params = ET.SubElement(element, "parameters") diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 24ab98895d..3fdd5e4a29 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -29,7 +29,7 @@ class Univariate(object): pass @abstractmethod - def to_xml(self, element_name): + def to_xml_element(self, element_name): return '' @abstractmethod @@ -92,7 +92,20 @@ class Discrete(Univariate): cv.check_greater_than('discrete probability', pk, 0.0, True) self._p = p - def to_xml(self, element_name): + def to_xml_element(self, element_name): + """Return XML representation of the discrete distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing discrete distribution data + + """ element = ET.Element(element_name) element.set("type", "discrete") @@ -153,7 +166,20 @@ class Uniform(Univariate): t.c = [0., 1.] return t - def to_xml(self, element_name): + def to_xml_element(self, element_name): + """Return XML representation of the uniform distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing uniform distribution data + + """ element = ET.Element(element_name) element.set("type", "uniform") element.set("parameters", '{} {}'.format(self.a, self.b)) @@ -196,7 +222,20 @@ class Maxwell(Univariate): cv.check_greater_than('Maxwell temperature', theta, 0.0) self._theta = theta - def to_xml(self, element_name): + def to_xml_element(self, element_name): + """Return XML representation of the Maxwellian distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing Maxwellian distribution data + + """ element = ET.Element(element_name) element.set("type", "maxwell") element.set("parameters", str(self.theta)) @@ -254,7 +293,20 @@ class Watt(Univariate): cv.check_greater_than('Watt b', b, 0.0) self._b = b - def to_xml(self, element_name): + def to_xml_element(self, element_name): + """Return XML representation of the Watt distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing Watt distribution data + + """ element = ET.Element(element_name) element.set("type", "watt") element.set("parameters", '{} {}'.format(self.a, self.b)) @@ -333,7 +385,20 @@ class Tabular(Univariate): cv.check_value('interpolation', interpolation, _INTERPOLATION_SCHEMES) self._interpolation = interpolation - def to_xml(self, element_name): + def to_xml_element(self, element_name): + """Return XML representation of the tabular distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing tabular distribution data + + """ element = ET.Element(element_name) element.set("type", "tabular") element.set("interpolation", self.interpolation) @@ -386,7 +451,7 @@ class Legendre(Univariate): self._legendre_polynomial = np.polynomial.legendre.Legendre( coefficients) - def to_xml(self, element_name): + def to_xml_element(self, element_name): raise NotImplementedError @@ -440,5 +505,5 @@ class Mixture(Univariate): Iterable, Univariate) self._distribution = distribution - def to_xml(self, element_name): + def to_xml_element(self, element_name): raise NotImplementedError diff --git a/openmc/volume.py b/openmc/volume.py index 511af62ed3..167971bdbd 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -183,7 +183,7 @@ class VolumeCalculation(object): vol.results = results return vol - def to_xml(self): + def to_xml_element(self): """Return XML representation of the volume calculation Returns From 08c9036bb69c781694f65ca6957f83cfe7387d20 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 1 Aug 2016 07:02:16 -0500 Subject: [PATCH 04/12] Ability to link volume calculation results in order to compute microscopic cross sections from mgxs --- openmc/cell.py | 60 +++++++++++++++++++++++++++++++++++++++--- openmc/geometry.py | 31 +++++++++------------- openmc/lattice.py | 17 ++++++------ openmc/material.py | 24 ++++++++++++++++- openmc/mgxs/library.py | 2 +- openmc/mgxs/mgxs.py | 41 +++++++++++++---------------- openmc/statepoint.py | 25 ++++++++++++++++-- openmc/summary.py | 11 ++++++++ openmc/universe.py | 31 ++++++++++++++++------ 9 files changed, 178 insertions(+), 64 deletions(-) diff --git a/openmc/cell.py b/openmc/cell.py index 9055f10e59..23246a0f68 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -85,6 +85,10 @@ class Cell(object): Array of offsets used for distributed cell searches distribcell_index : int Index of this cell in distribcell arrays + volume_information : dict + Estimate of the volume and total number of atoms of each nuclide from a + stochastic volume calculation. This information is set with the + :meth:`Cell.add_volume_information` method. """ @@ -100,6 +104,7 @@ class Cell(object): self._translation = None self._offsets = None self._distribcell_index = None + self._volume_information = None def __contains__(self, point): if self.region is None: @@ -212,6 +217,10 @@ class Cell(object): def distribcell_index(self): return self._distribcell_index + @property + def volume_information(self): + return self._volume_information + @id.setter def id(self, cell_id): if cell_id is None: @@ -351,6 +360,22 @@ class Cell(object): else: self.region = Intersection(self.region, region) + def add_volume_information(self, volume_calc): + """Add volume information to a cell. + + Parameters + ---------- + volume_calc : openmc.VolumeCalculation + Results from a stochastic volume calculation + + """ + for cell_id in volume_calc.results: + if cell_id == self.id: + self._volume_information = volume_calc.results[cell_id] + break + else: + raise ValueError('No volume information found for this cell.') + def get_cell_instance(self, path, distribcell_index): # If the Cell is filled by a Material @@ -368,8 +393,19 @@ class Cell(object): return offset - def get_all_nuclides(self): - """Return all nuclides contained in the cell + def get_nuclides(self): + """Returns all nuclides in the cell + + Returns + ------- + nuclides : list + List of nuclide names + + """ + return self.fill.get_nuclides() if self.fill_type != 'void' else [] + + def get_nuclide_densities(self): + """Return all nuclides contained in the cell and their densities Returns ------- @@ -381,8 +417,24 @@ class Cell(object): nuclides = OrderedDict() - if self.fill_type != 'void': - nuclides.update(self.fill.get_all_nuclides()) + if self.fill_type == 'material': + nuclides.update(self.fill.get_nuclide_densities()) + elif self.fill_type == 'void': + pass + else: + if self.volume_information is not None: + volume = self.volume_information['volume'][0] + for full_name, atoms in self.volume_information['atoms']: + name, xs = full_name.split('.') + nuclide = openmc.Nuclide(name, xs) + density = 1.0e-24 * atoms[0]/volume # density in atoms/b-cm + nuclides[name] = (nuclide, density) + else: + raise RuntimeError( + 'Volume information is needed to calculate microscopic cross ' + 'sections for cell {}. This can be done by running a ' + 'stochastic volume calculation via the ' + 'openmc.VolumeCalculation object'.format(self.id)) return nuclides diff --git a/openmc/geometry.py b/openmc/geometry.py index 14fd48fb1e..8314a9049c 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -50,6 +50,19 @@ class Geometry(object): self._root_universe = root_universe + def add_volume_information(self, volume_calc): + """Add volume information to from a stochastic volume calculation. + + Parameters + ---------- + volume_calc : openmc.VolumeCalculation + Results from a stochastic volume calculation + + """ + for cell in self.get_all_cells(): + if cell.id in volume_calc.results: + cell.add_volume_information(volume_calc) + def export_to_xml(self): """Create a geometry.xml file that can be used for a simulation. @@ -166,24 +179,6 @@ class Geometry(object): universes.sort(key=lambda x: x.id) return universes - def get_all_nuclides(self): - """Return all nuclides assigned to a material in the geometry - - Returns - ------- - list of openmc.Nuclide - Nuclides in the geometry - - """ - - nuclides = OrderedDict() - materials = self.get_all_materials() - - for material in materials: - nuclides.update(material.get_all_nuclides()) - - return nuclides - def get_all_materials(self): """Return all materials assigned to a cell diff --git a/openmc/lattice.py b/openmc/lattice.py index 81144e4d81..b26aab4c67 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -144,25 +144,26 @@ class Lattice(object): return univs - def get_all_nuclides(self): - """Return all nuclides contained in the lattice + def get_nuclides(self): + """Returns all nuclides in the lattice Returns ------- - nuclides : collections.OrderedDict - Dictionary whose keys are nuclide names and values are 2-tuples of - (nuclide, density) + nuclides : list + List of nuclide names """ - nuclides = OrderedDict() + nuclides = [] # Get all unique Universes contained in each of the lattice cells unique_universes = self.get_unique_universes() # Append all Universes containing each cell to the dictionary - for universe_id, universe in unique_universes.items(): - nuclides.update(universe.get_all_nuclides()) + for universe in unique_universes.values(): + for nuclide in universe.get_nuclides(): + if nuclide not in nuclides: + nuclides.append(nuclide) return nuclides diff --git a/openmc/material.py b/openmc/material.py index d7ffd04e94..5603d9d25c 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -492,9 +492,31 @@ class Material(object): for element, percent, percent_type in self._elements: element.scattering = 'iso-in-lab' - def get_all_nuclides(self): + def get_nuclides(self): """Returns all nuclides in the material + Returns + ------- + nuclides : list + List of nuclide names + + """ + + nuclides = [] + + for nuclide, density, density_type in self._nuclides: + nuclides.append(nuclide.name) + + for element, density, density_type in self._elements: + # Expand natural element into isotopes + for isotope, abundance in element.expand(): + nuclides.append(isotope.name) + + return nuclides + + def get_nuclide_densities(self): + """Returns all nuclides in the material and their densities + Returns ------- nuclides : dict diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 84571ce236..ee11d0ef64 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -1008,7 +1008,7 @@ class Library(object): # Create the xsdata object and add it to the mgxs_file for i, domain in enumerate(self.domains): if self.by_nuclide: - nuclides = list(domain.get_all_nuclides().keys()) + nuclides = domain.get_nuclides() else: nuclides = ['total'] for nuclide in nuclides: diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index e1c2e62217..7aa12d0418 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -288,9 +288,7 @@ class MGXS(object): # If this is a by-nuclide cross-section, add nuclides to Tally if self.by_nuclide and score != 'flux': - all_nuclides = self.get_all_nuclides() - for nuclide in all_nuclides: - self._tallies[key].nuclides.append(nuclide) + self._tallies[key].nuclides += self.get_nuclides() else: self._tallies[key].nuclides.append('total') @@ -329,14 +327,14 @@ class MGXS(object): @property def num_nuclides(self): if self.by_nuclide: - return len(self.get_all_nuclides()) + return len(self.get_nuclides()) else: return 1 @property def nuclides(self): if self.by_nuclide: - return self.get_all_nuclides() + return self.get_nuclides() else: return 'sum' @@ -503,7 +501,7 @@ class MGXS(object): mgxs.name = name return mgxs - def get_all_nuclides(self): + def get_nuclides(self): """Get all nuclides in the cross section's spatial domain. Returns @@ -528,8 +526,7 @@ class MGXS(object): # Otherwise, return all nuclides in the spatial domain else: - nuclides = self.domain.get_all_nuclides() - return list(nuclides.keys()) + return self.domain.get_nuclides() def get_nuclide_density(self, nuclide): """Get the atomic number density in units of atoms/b-cm for a nuclide @@ -556,7 +553,7 @@ class MGXS(object): cv.check_type('nuclide', nuclide, basestring) # Get list of all nuclides in the spatial domain - nuclides = self.domain.get_all_nuclides() + nuclides = self.domain.get_nuclide_densities() if nuclide not in nuclides: msg = 'Unable to get density for nuclide "{0}" which is not in ' \ @@ -597,14 +594,14 @@ class MGXS(object): # Sum the atomic number densities for all nuclides if nuclides == 'sum': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() densities = np.zeros(1, dtype=np.float) for nuclide in nuclides: densities[0] += self.get_nuclide_density(nuclide) # Tabulate the atomic number densities for all nuclides elif nuclides == 'all': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() densities = np.zeros(self.num_nuclides, dtype=np.float) for i, nuclide in enumerate(nuclides): densities[i] += self.get_nuclide_density(nuclide) @@ -635,7 +632,7 @@ class MGXS(object): # If computing xs for each nuclide, replace CrossNuclides with originals if self.by_nuclide: self.xs_tally._nuclides = [] - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() for nuclide in nuclides: self.xs_tally.nuclides.append(openmc.Nuclide(nuclide)) @@ -796,7 +793,7 @@ class MGXS(object): # Construct a collection of the nuclides to retrieve from the xs tally if self.by_nuclide: if nuclides == 'all' or nuclides == 'sum' or nuclides == ['sum']: - query_nuclides = self.get_all_nuclides() + query_nuclides = self.get_nuclides() else: query_nuclides = nuclides else: @@ -1164,7 +1161,7 @@ class MGXS(object): # Construct a collection of the nuclides to report if self.by_nuclide: if nuclides == 'all': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() elif nuclides == 'sum': nuclides = ['sum'] else: @@ -1302,7 +1299,7 @@ class MGXS(object): # Construct a collection of the nuclides to report if self.by_nuclide: if nuclides == 'all': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() densities = np.zeros(len(nuclides), dtype=np.float) elif nuclides == 'sum': nuclides = ['sum'] @@ -1484,7 +1481,7 @@ class MGXS(object): if self.by_nuclide and nuclides == 'sum': # Use tally summation to sum across all nuclides - query_nuclides = self.get_all_nuclides() + query_nuclides = self.get_nuclides() xs_tally = self.xs_tally.summation(nuclides=query_nuclides) df = xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) @@ -1790,7 +1787,7 @@ class MatrixMGXS(MGXS): # Construct a collection of the nuclides to retrieve from the xs tally if self.by_nuclide: if nuclides == 'all' or nuclides == 'sum' or nuclides == ['sum']: - query_nuclides = self.get_all_nuclides() + query_nuclides = self.get_nuclides() else: query_nuclides = nuclides else: @@ -1937,7 +1934,7 @@ class MatrixMGXS(MGXS): # Construct a collection of the nuclides to report if self.by_nuclide: if nuclides == 'all': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() if nuclides == 'sum': nuclides = ['sum'] else: @@ -3622,7 +3619,7 @@ class ScatterMatrixXS(MatrixMGXS): # Construct a collection of the nuclides to retrieve from the xs tally if self.by_nuclide: if nuclides == 'all' or nuclides == 'sum' or nuclides == ['sum']: - query_nuclides = self.get_all_nuclides() + query_nuclides = self.get_nuclides() else: query_nuclides = nuclides else: @@ -3789,7 +3786,7 @@ class ScatterMatrixXS(MatrixMGXS): # Construct a collection of the nuclides to report if self.by_nuclide: if nuclides == 'all': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() if nuclides == 'sum': nuclides = ['sum'] else: @@ -4595,7 +4592,7 @@ class Chi(MGXS): nu_fission_out = self.tallies['nu-fission-out'] # Sum out all nuclides - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() nu_fission_in = nu_fission_in.summation(nuclides=nuclides) nu_fission_out = nu_fission_out.summation(nuclides=nuclides) @@ -4615,7 +4612,7 @@ class Chi(MGXS): # Get chi for all nuclides in the domain elif nuclides == 'all': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() xs = self.xs_tally.get_values(filters=filters, filter_bins=filter_bins, nuclides=nuclides, value=value) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 6337746650..58a79fa219 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -2,6 +2,7 @@ import sys import re import os import warnings +import glob import numpy as np @@ -22,8 +23,9 @@ class StatePoint(object): filename : str Path to file to load autolink : bool, optional - Whether to automatically link in metadata from a summary.h5 - file. Defaults to True. + Whether to automatically link in metadata from a summary.h5 file and + stochastic volume calculation results from volume_*.h5 files. Defaults + to True. Attributes ---------- @@ -143,6 +145,12 @@ class StatePoint(object): su = openmc.Summary(path_summary) self.link_with_summary(su) + path_volume = os.path.join(os.path.dirname(filename), 'volume_*.h5') + for path_i in glob.glob(path_volume): + if re.search(r'volume_\d+\.h5', path_i): + vol = openmc.VolumeCalculation.from_hdf5(path_i) + self.add_volume_information(vol) + def close(self): self._f.close() @@ -501,6 +509,19 @@ class StatePoint(object): for tally_id in self.tallies: self.tallies[tally_id].sparse = self.sparse + def add_volume_information(self, volume_calc): + """Add volume information to the geometry within the file + + Parameters + ---------- + volume_calc : openmc.VolumeCalculation + Results from a stochastic volume calculation + + """ + if self.summary is not None: + self.summary.add_volume_information(volume_calc) + + def get_tally(self, scores=[], filters=[], nuclides=[], name=None, id=None, estimator=None, exact_filters=False, exact_nuclides=False, exact_scores=False): diff --git a/openmc/summary.py b/openmc/summary.py index 8c626940e5..6fbace72ea 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -594,6 +594,17 @@ class Summary(object): # Add Tally to the global dictionary of all Tallies self.tallies[tally_id] = tally + def add_volume_information(self, volume_calc): + """Add volume information to the geometry within the summary file + + Parameters + ---------- + volume_calc : openmc.VolumeCalculation + Results from a stochastic volume calculation + + """ + self.openmc_geometry.add_volume_information(volume_calc) + def get_material_by_id(self, material_id): """Return a Material object given the material id diff --git a/openmc/universe.py b/openmc/universe.py index c8e7fcab1e..f2b8641699 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -349,7 +349,27 @@ class Universe(object): # Return the offset computed at all nested Universe levels return offset - def get_all_nuclides(self): + def get_nuclides(self): + """Returns all nuclides in the universe + + Returns + ------- + nuclides : list + List of nuclide names + + """ + + nuclides = [] + + # Append all Nuclides in each Cell in the Universe to the dictionary + for cell in self.cells.values(): + for nuclide in cell.get_nuclides(): + if nuclide not in nuclides: + nuclides.append(nuclide) + + return nuclides + + def get_nuclide_densities(self): """Return all nuclides contained in the universe Returns @@ -360,13 +380,8 @@ class Universe(object): """ - nuclides = OrderedDict() - - # Append all Nuclides in each Cell in the Universe to the dictionary - for cell in self._cells.values(): - nuclides.update(cell.get_all_nuclides()) - - return nuclides + raise NotImplementedError('Determining average nuclide densities over ' + 'an entire universe not yet supported.') def get_all_cells(self): """Return all cells that are contained within the universe From aaefad9e2ea99c8edec4874741c37c128fa9e0e9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 1 Aug 2016 10:25:08 -0500 Subject: [PATCH 05/12] Add volume calculation test --- tests/test_volume_calc/inputs_true.dat | 1 + tests/test_volume_calc/results_true.dat | 13 ++++ tests/test_volume_calc/test_volume_calc.py | 82 ++++++++++++++++++++++ tests/testing_harness.py | 1 + 4 files changed, 97 insertions(+) create mode 100644 tests/test_volume_calc/inputs_true.dat create mode 100644 tests/test_volume_calc/results_true.dat create mode 100644 tests/test_volume_calc/test_volume_calc.py diff --git a/tests/test_volume_calc/inputs_true.dat b/tests/test_volume_calc/inputs_true.dat new file mode 100644 index 0000000000..2a3f25e6d0 --- /dev/null +++ b/tests/test_volume_calc/inputs_true.dat @@ -0,0 +1 @@ +b2de6ac20ca2ca38b00d61adae707d05519f6130eea25ab050220bac005cb6d23fa5812bf876307e92e90a8806f0b371500c5611ab23f80b7c073f118eb31649 \ No newline at end of file diff --git a/tests/test_volume_calc/results_true.dat b/tests/test_volume_calc/results_true.dat new file mode 100644 index 0000000000..e7a39a9312 --- /dev/null +++ b/tests/test_volume_calc/results_true.dat @@ -0,0 +1,13 @@ +k-combined: 4.165451e-02 3.582531e-04 +Cell 1: 31.4693 +/- 0.0721 cm^3 +Cell 2: 2.0933 +/- 0.0310 cm^3 +Cell 3: 2.0486 +/- 0.0307 cm^3 + Cell Nuclide Atoms Uncertainty +0 1 U235.71c 3.481769e+23 7.979991e+20 +1 1 Mo99.71c 3.481769e+22 7.979991e+19 +2 2 H1.71c 1.399770e+23 2.072914e+21 +3 2 O16.71c 6.998852e+22 1.036457e+21 +4 2 B10.71c 6.998852e+18 1.036457e+17 +5 3 H1.71c 1.369920e+23 2.051689e+21 +6 3 O16.71c 6.849599e+22 1.025844e+21 +7 3 B10.71c 6.849599e+18 1.025844e+17 diff --git a/tests/test_volume_calc/test_volume_calc.py b/tests/test_volume_calc/test_volume_calc.py new file mode 100644 index 0000000000..a6d779a25a --- /dev/null +++ b/tests/test_volume_calc/test_volume_calc.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +import openmc + + +class VolumeTest(PyAPITestHarness): + def _build_inputs(self): + # Define materials + water = openmc.Material(1) + water.add_nuclide('H1', 2.0) + water.add_nuclide('O16', 1.0) + water.add_nuclide('B10', 0.0001) + water.add_s_alpha_beta('c_H_in_H2O', '71t') + water.set_density('g/cc', 1.0) + + fuel = openmc.Material(2) + fuel.add_nuclide('U235', 1.0) + fuel.add_nuclide('Mo99', 0.1) + fuel.set_density('g/cc', 4.5) + + materials = openmc.Materials((water, fuel)) + materials.default_xs = '71c' + materials.export_to_xml() + + cyl = openmc.ZCylinder(1, R=1.0, boundary_type='vacuum') + top_sphere = openmc.Sphere(2, z0=5., R=1., boundary_type='vacuum') + top_plane = openmc.ZPlane(3, z0=5.) + bottom_sphere = openmc.Sphere(4, z0=-5., R=1., boundary_type='vacuum') + bottom_plane = openmc.ZPlane(5, z0=-5.) + + # Define geometry + inside_cyl = openmc.Cell(1, fill=fuel, region=-cyl & -top_plane & +bottom_plane) + top_hemisphere = openmc.Cell(2, fill=water, region=-top_sphere & +top_plane) + bottom_hemisphere = openmc.Cell(3, fill=water, region=-bottom_sphere & -top_plane) + root = openmc.Universe(0, cells=(inside_cyl, top_hemisphere, bottom_hemisphere)) + + geometry = openmc.Geometry() + geometry.root_universe = root + geometry.export_to_xml() + + # Set up stochastic volume calculation + vol_calc = openmc.VolumeCalculation( + [inside_cyl, top_hemisphere, bottom_hemisphere], + 100000) + + # Define settings + settings = openmc.Settings() + settings.particles = 1000 + settings.batches = 4 + settings.inactive = 0 + settings.source = openmc.Source(space=openmc.stats.Box( + [-1., -1., -5.], [1., 1., 5.])) + settings.volume_calculations = vol_calc + settings.export_to_xml() + + def _get_results(self): + # Read the statepoint file. + statepoint = os.path.join(os.getcwd(), self._sp_name) + sp = openmc.StatePoint(statepoint) + + # Write out k-combined. + outstr = 'k-combined: {:12.6e} {:12.6e}\n'.format(*sp.k_combined) + + # Read volume calculation results + vol = openmc.VolumeCalculation.from_hdf5( + os.path.join(os.getcwd(), 'volume_1.h5')) + + # Write cell volumes and total # of atoms for each nuclide + for cell_id, results in sorted(vol.results.items()): + outstr += 'Cell {0}: {1[0]:.4f} +/- {1[1]:.4f} cm^3\n'.format( + cell_id, results['volume']) + outstr += str(vol.atoms_dataframe) + '\n' + + return outstr + +if __name__ == '__main__': + harness = VolumeTest('statepoint.4.h5') + harness.main() diff --git a/tests/testing_harness.py b/tests/testing_harness.py index d360184045..19a7ffb06c 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -137,6 +137,7 @@ class TestHarness(object): output.append(os.path.join(os.getcwd(), 'tallies.out')) output.append(os.path.join(os.getcwd(), 'results_test.dat')) output.append(os.path.join(os.getcwd(), 'summary.h5')) + output += glob.glob(os.path.join(os.getcwd(), 'volume_*.h5')) for f in output: if os.path.exists(f): os.remove(f) From 413a58a363ee8fbedeffa652f068f647f0fcb0dc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 1 Aug 2016 11:15:48 -0500 Subject: [PATCH 06/12] Change a few mentions of ACE format cross sections and update manpage --- man/man1/openmc.1 | 16 +++++++++++++--- readme.rst | 2 +- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index e69360a7c1..7f18ab28a2 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -3,7 +3,7 @@ openmc \- Executes the OpenMC Monte Carlo code .SH DESCRIPTION This command is used to execute the OpenMC Monte Carlo code. It is assumed that -a set of XML input files has already been created and that ACE format cross +a set of XML input files has already been created and that HDF5 format cross sections are available. .SH SYNOPSIS \fBopenmc\fR [\fIoptions\fR] [\fIpath\fR] @@ -40,11 +40,21 @@ The behavior of .B openmc is affected by the following environment variables. .TP -.B CROSS_SECTIONS +.B OPENMC_CROSS_SECTIONS Indicates the default path to the cross_sections.xml summary file that is used -to locate ACE format cross section libraries if the user has not specified the +to locate HDF5 format cross section libraries if the user has not specified the tag in .I settings.xml\fP. +.TP +.B OPENMC_MG_CROSS_SECTIONS +Indicates the default path to the mgxs.xml file that contains multi-group cross +section libraries if the user has not specified the tag in +.I settings.xml\fP. +.TP +.B OPENMC_MULTIPOLE_LIBRARY +Indicates the default path to a directory containing windowed multipole data if +the user has not specified the tag in +.I settings.xml\fP. .SH LICENSE Copyright \(co 2011-2016 Massachusetts Institute of Technology. .PP diff --git a/readme.rst b/readme.rst index 03964d9436..90484ad494 100644 --- a/readme.rst +++ b/readme.rst @@ -7,7 +7,7 @@ OpenMC Monte Carlo Particle Transport Code The OpenMC project aims to provide a fully-featured Monte Carlo particle transport code based on modern methods. It is a constructive solid geometry, -continuous-energy transport code that uses ACE format cross sections. The +continuous-energy transport code that uses HDF5 format cross sections. The project started under the Computational Reactor Physics Group at MIT. Complete documentation on the usage of OpenMC is hosted on Read the Docs at From af27542af509b250e2872d48bf2708d7210ecfd5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 1 Aug 2016 21:48:02 -0500 Subject: [PATCH 07/12] Fix a few docstrings --- openmc/cell.py | 2 +- openmc/geometry.py | 2 +- openmc/lattice.py | 2 +- openmc/material.py | 2 +- openmc/universe.py | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/cell.py b/openmc/cell.py index 23246a0f68..c4a0952e91 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -398,7 +398,7 @@ class Cell(object): Returns ------- - nuclides : list + nuclides : list of str List of nuclide names """ diff --git a/openmc/geometry.py b/openmc/geometry.py index 8314a9049c..c827f4b2a7 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -51,7 +51,7 @@ class Geometry(object): self._root_universe = root_universe def add_volume_information(self, volume_calc): - """Add volume information to from a stochastic volume calculation. + """Add volume information from a stochastic volume calculation. Parameters ---------- diff --git a/openmc/lattice.py b/openmc/lattice.py index b26aab4c67..c6a5af4109 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -149,7 +149,7 @@ class Lattice(object): Returns ------- - nuclides : list + nuclides : list of str List of nuclide names """ diff --git a/openmc/material.py b/openmc/material.py index 5603d9d25c..aeebf385c6 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -497,7 +497,7 @@ class Material(object): Returns ------- - nuclides : list + nuclides : list of str List of nuclide names """ diff --git a/openmc/universe.py b/openmc/universe.py index f2b8641699..c8fb89a384 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -354,7 +354,7 @@ class Universe(object): Returns ------- - nuclides : list + nuclides : list of str List of nuclide names """ From 99adb712aba31bac85245940f09019470427a113 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 1 Aug 2016 21:51:31 -0500 Subject: [PATCH 08/12] gfortran 4.6 doesn't like firstprivate objects with allocatable components --- src/stl_vector.F90 | 15 +++++++++++++++ src/volume_calc.F90 | 12 ++++++++++-- 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/stl_vector.F90 b/src/stl_vector.F90 index c7f2246ff8..06f487dc1e 100644 --- a/src/stl_vector.F90 +++ b/src/stl_vector.F90 @@ -114,6 +114,11 @@ contains ! Since integer is trivially destructible, we only need to set size to zero ! and can leave capacity as is this%size_ = 0 + if (allocated(this % data)) then + this%capacity_ = size(this % data) + else + this%capacity_ = 0 + end if end subroutine clear_int subroutine initialize_fill_int(this, n, val) @@ -249,6 +254,11 @@ contains ! Since real is trivially destructible, we only need to set size to zero and ! can leave capacity as is this%size_ = 0 + if (allocated(this % data)) then + this%capacity_ = size(this % data) + else + this%capacity_ = 0 + end if end subroutine clear_real subroutine initialize_fill_real(this, n, val) @@ -384,6 +394,11 @@ contains ! Since char is trivially destructible, we only need to set size to zero and ! can leave capacity as is this%size_ = 0 + if (allocated(this % data)) then + this%capacity_ = size(this % data) + else + this%capacity_ = 0 + end if end subroutine clear_char subroutine initialize_fill_char(this, n, val) diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index cb36f0b7f3..eac96f9ba7 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -140,8 +140,16 @@ contains call p % initialize() -!$omp parallel private(i, j, k, i_cell, i_material, level, found_cell) & -!$omp& firstprivate(p, indices, hits) +!$omp parallel private(i, j, k, i_cell, i_material, level, found_cell, & +!$omp& indices, hits) firstprivate(p) + + ! Reset vectors -- this is really to get around a gfortran 4.6 bug. Ideally, + ! indices and hits would just be firstprivate but 4.6 complains because they + ! have allocatable components... + do i_cell = 1, size(this % cell_id) + call indices(i_cell) % clear() + call hits(i_cell) % clear() + end do call prn_set_stream(STREAM_VOLUME) From a5763ae599eb082daf170d29ab99dbe7fa537ef9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Aug 2016 15:27:46 -0500 Subject: [PATCH 09/12] Refactor volume calculation to allow material/universe volumes as well --- openmc/cell.py | 11 +- openmc/geometry.py | 7 +- openmc/volume.py | 101 ++++++---- src/relaxng/settings.rnc | 6 +- src/relaxng/settings.rng | 16 +- src/volume_calc.F90 | 209 ++++++++++++++------- src/volume_header.F90 | 31 ++- tests/test_volume_calc/inputs_true.dat | 2 +- tests/test_volume_calc/results_true.dat | 24 ++- tests/test_volume_calc/test_volume_calc.py | 31 +-- 10 files changed, 298 insertions(+), 140 deletions(-) diff --git a/openmc/cell.py b/openmc/cell.py index c4a0952e91..fe8b4a52cd 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -369,10 +369,13 @@ class Cell(object): Results from a stochastic volume calculation """ - for cell_id in volume_calc.results: - if cell_id == self.id: - self._volume_information = volume_calc.results[cell_id] - break + if volume_calc.domain_type == 'cell': + for cell_id in volume_calc.results: + if cell_id == self.id: + self._volume_information = volume_calc.results[cell_id] + break + else: + raise ValueError('No volume information found for this cell.') else: raise ValueError('No volume information found for this cell.') diff --git a/openmc/geometry.py b/openmc/geometry.py index c827f4b2a7..b2a9b5e4aa 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -59,9 +59,10 @@ class Geometry(object): Results from a stochastic volume calculation """ - for cell in self.get_all_cells(): - if cell.id in volume_calc.results: - cell.add_volume_information(volume_calc) + if volume_calc.domain_type == 'cell': + for cell in self.get_all_cells(): + if cell.id in volume_calc.results: + cell.add_volume_information(volume_calc) def export_to_xml(self): """Create a geometry.xml file that can be used for a simulation. diff --git a/openmc/volume.py b/openmc/volume.py index 167971bdbd..e9b1dac1d3 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -5,7 +5,7 @@ from xml.etree import ElementTree as ET import numpy as np import pandas as pd -from openmc import Cell, Union +import openmc import openmc.checkvalue as cv @@ -14,8 +14,8 @@ class VolumeCalculation(object): Parameters ---------- - cells : Iterable of Cell - Cells to find volumes of + domains : Iterable of openmc.Cell, openmc.Material, or openmc.Universe + Domains to find volumes of samples : int Number of samples used to generate volume estimates lower_left : Iterable of float @@ -29,8 +29,10 @@ class VolumeCalculation(object): Attributes ---------- - cell_ids : Iterable of int - IDs of cells to find volumes of + ids : Iterable of int + IDs of domains to find volumes of + domain_type : {'cell', 'material', 'universe'} + Type of each domain samples : int Number of samples used to generate volume estimates lower_left : Iterable of float @@ -38,23 +40,31 @@ class VolumeCalculation(object): upper_right : Iterable of float Upper-right coordinates of bounding box used to sample points results : dict - Dictionary whose keys are unique IDs of cells and values are + Dictionary whose keys are unique IDs of domains and values are dictionaries with calculated volumes and total number of atoms for each - nuclide present in the cell. + nuclide present in the domain. volumes : dict - Dictionary whose keys are unique IDs of cells and values are the + Dictionary whose keys are unique IDs of domains and values are the estimated volumes atoms_dataframe : pandas.DataFrame DataFrame showing the estimated number of atoms for each nuclide present - in each cell specified. + in each domain specified. """ - def __init__(self, cells, samples, lower_left=None, + def __init__(self, domains, samples, lower_left=None, upper_right=None): self._results = None - cv.check_type('cells', cells, Iterable, Cell) - self.cell_ids = [c.id for c in cells] + cv.check_type('domains', domains, Iterable, + (openmc.Cell, openmc.Material, openmc.Universe)) + if isinstance(domains[0], openmc.Cell): + self._domain_type = 'cell' + elif isinstance(domains[0], openmc.Material): + self._domain_type = 'material' + elif isinstance(domains[0], openmc.Universe): + self._domain_type = 'universe' + self.ids = [d.id for d in domains] + self.samples = samples if lower_left is not None: @@ -64,17 +74,21 @@ class VolumeCalculation(object): 'should be specified') self.upper_right = upper_right else: - ll, ur = Union(*[c.region for c in cells]).bounding_box - if np.any(np.isinf(ll)) or np.any(np.isinf(ur)): + if self.domain_type == 'cell': + ll, ur = openmc.Union(*[c.region for c in domains]).bounding_box + if np.any(np.isinf(ll)) or np.any(np.isinf(ur)): + raise ValueError('Could not automatically determine bounding ' + 'box for stochastic volume calculation.') + else: + self.lower_left = ll + self.upper_right = ur + else: raise ValueError('Could not automatically determine bounding box ' 'for stochastic volume calculation.') - else: - self.lower_left = ll - self.upper_right = ur @property - def cell_ids(self): - return self._cell_ids + def ids(self): + return self._ids @property def samples(self): @@ -92,6 +106,10 @@ class VolumeCalculation(object): def results(self): return self._results + @property + def domain_type(self): + return self._domain_type + @property def volumes(self): return {uid: results['volume'] for uid, results in self.results.items()} @@ -99,17 +117,18 @@ class VolumeCalculation(object): @property def atoms_dataframe(self): items = [] - columns = ['Cell', 'Nuclide', 'Atoms', 'Uncertainty'] - for cell_id, results in self.results.items(): + columns = [self.domain_type.capitalize(), 'Nuclide', 'Atoms', + 'Uncertainty'] + for uid, results in self.results.items(): for name, atoms in results['atoms']: - items.append((cell_id, name, atoms[0], atoms[1])) + items.append((uid, name, atoms[0], atoms[1])) return pd.DataFrame.from_records(items, columns=columns) - @cell_ids.setter - def cell_ids(self, cell_ids): - cv.check_type('cell IDs', cell_ids, Iterable, Real) - self._cell_ids = cell_ids + @ids.setter + def ids(self, ids): + cv.check_type('domain IDs', ids, Iterable, Real) + self._ids = ids @samples.setter def samples(self, samples): @@ -154,16 +173,17 @@ class VolumeCalculation(object): import h5py with h5py.File(filename, 'r') as f: + domain_type = f.attrs['domain_type'].decode() samples = f.attrs['samples'] lower_left = f.attrs['lower_left'] upper_right = f.attrs['upper_right'] results = {} - cell_ids = [] + ids = [] for obj_name in f: - if obj_name.startswith('cell_'): - cell_id = int(obj_name[5:]) - cell_ids.append(cell_id) + if obj_name.startswith('domain_'): + domain_id = int(obj_name[7:]) + ids.append(domain_id) group = f[obj_name] volume = tuple(group['volume'].value) nucnames = group['nuclides'].value @@ -172,14 +192,19 @@ class VolumeCalculation(object): atom_list = [] for name_i, atoms_i in zip(nucnames, atoms): atom_list.append((name_i.decode(), tuple(atoms_i))) - results[cell_id] = {'volume': volume, 'atoms': atom_list} + results[domain_id] = {'volume': volume, 'atoms': atom_list} - # Instantiate some throw-away cells that are used by the constructor to - # assign IDs - cells = [Cell(uid) for uid in cell_ids] + # Instantiate some throw-away domains that are used by the constructor + # to assign IDs + if domain_type == 'cell': + domains = [openmc.Cell(uid) for uid in ids] + elif domain_type == 'material': + domains = [openmc.Material(uid) for uid in ids] + elif domain_type == 'universe': + domains = [openmc.Universe(uid) for uid in ids] # Instantiate the class and assign results - vol = cls(cells, samples, lower_left, upper_right) + vol = cls(domains, samples, lower_left, upper_right) vol.results = results return vol @@ -193,8 +218,10 @@ class VolumeCalculation(object): """ element = ET.Element("volume_calc") - cell_elem = ET.SubElement(element, "cells") - cell_elem.text = ' '.join(str(uid) for uid in self.cell_ids) + dt_elem = ET.SubElement(element, "domain_type") + dt_elem.text = self.domain_type + id_elem = ET.SubElement(element, "domain_ids") + id_elem.text = ' '.join(str(uid) for uid in self.ids) samples_elem = ET.SubElement(element, "samples") samples_elem.text = str(self.samples) ll_elem = ET.SubElement(element, "lower_left") diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 78a32171b9..e23cbdd6cd 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -143,8 +143,10 @@ element settings { element verbosity { xsd:positiveInteger }? & element volume_calc { - (element cells { list { xsd:positiveInteger+ } } | - attribute cells { list { xsd:positiveInteger+ } }) & + (element domain_type { xsd:string } | + attribute domain_type { xsd:string }) & + (element domain_ids { list { xsd:integer+ } } | + attribute domain_ids { list { xsd:integer+ } }) & (element samples { xsd:positiveInteger } | attribute samples { xsd:positiveInteger }) & (element lower_left { list { xsd:double+ } } | diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 64fe42239d..38acea5985 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -629,17 +629,25 @@ - + + + + + + + + + - + - + - + diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index eac96f9ba7..03e241a66b 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -9,7 +9,7 @@ module volume_calc use geometry, only: find_cell use global use hdf5_interface, only: file_create, file_close, write_attribute, & - create_group, close_group, write_dataset + create_group, close_group, write_dataset, write_attribute_string use output, only: write_message, header use message_passing use particle_header, only: Particle @@ -33,7 +33,8 @@ contains subroutine run_volume_calculations() integer :: i, j integer :: n - real(8), allocatable :: volume(:,:) ! volume mean/stdev in each cell + real(8), allocatable :: volume(:,:) ! volume mean/stdev in each domain + character(10) :: domain_type character(MAX_FILE_LEN) :: filename ! filename for HDF5 file type(Timer) :: time_volume ! timer for volume calculation type(VectorInt), allocatable :: nuclide_vec(:) ! indices in nuclides array @@ -46,7 +47,7 @@ contains end if do i = 1, size(volume_calcs) - n = size(volume_calcs(i) % cell_id) + n = size(volume_calcs(i) % domain_id) allocate(nuclide_vec(n)) allocate(atoms_vec(n), uncertainty_vec(n)) allocate(volume(2,n)) @@ -60,11 +61,20 @@ contains uncertainty_vec) if (master) then - ! Display cell volumes - do j = 1, size(volume_calcs(i) % cell_id) - call write_message(" Cell " // trim(to_str(volume_calcs(i) % & - cell_id(j))) // ": " // trim(to_str(volume(1,j))) // " +/- " // & - trim(to_str(volume(2,j))) // " cm^3") + select case (volume_calcs(i) % domain_type) + case (FILTER_CELL) + domain_type = ' Cell' + case (FILTER_MATERIAL) + domain_type = ' Material' + case (FILTER_UNIVERSE) + domain_type = ' Universe' + end select + + ! Display domain volumes + do j = 1, size(volume_calcs(i) % domain_id) + call write_message(trim(domain_type) // " " // trim(to_str(& + volume_calcs(i) % domain_id(j))) // ": " // trim(to_str(& + volume(1,j))) // " +/- " // trim(to_str(volume(2,j))) // " cm^3") end do call write_message("") @@ -85,13 +95,13 @@ contains end subroutine run_volume_calculations !=============================================================================== -! GET_VOLUME stochastically determines the volume of a set of cells along with -! the average number densities of nuclides within the cell +! GET_VOLUME stochastically determines the volume of a set of domains along with +! the average number densities of nuclides within the domain !=============================================================================== subroutine get_volume(this, volume, nuclide_vec, atoms_vec, uncertainty_vec) type(VolumeCalculation), intent(in) :: this - real(8), intent(out) :: volume(:,:) ! volume mean/stdev in each cell + real(8), intent(out) :: volume(:,:) ! volume mean/stdev in each domain type(VectorInt), intent(out) :: nuclide_vec(:) ! indices in nuclides array type(VectorReal), intent(out) :: atoms_vec(:) ! total # of atoms of each nuclide type(VectorReal), intent(out) :: uncertainty_vec(:) ! uncertainty of total # of atoms @@ -99,22 +109,22 @@ contains ! Variables that are private to each thread integer(8) :: i integer :: j, k - integer :: i_cell ! index in cell_id array + integer :: i_domain ! index in domain_id array integer :: i_material ! index in materials array integer :: level ! local coordinate level logical :: found_cell - type(VectorInt) :: indices(size(this % cell_id)) ! List of material indices - type(VectorInt) :: hits(size(this % cell_id)) ! Number of hits for each material + type(VectorInt) :: indices(size(this % domain_id)) ! List of material indices + type(VectorInt) :: hits(size(this % domain_id)) ! Number of hits for each material type(Particle) :: p ! Shared variables integer :: i_start, i_end ! Starting/ending sample for each process - type(VectorInt) :: master_indices(size(this % cell_id)) - type(VectorInt) :: master_hits(size(this % cell_id)) + type(VectorInt) :: master_indices(size(this % domain_id)) + type(VectorInt) :: master_hits(size(this % domain_id)) ! Variables used outside of parallel region integer :: i_nuclide ! index in nuclides array - integer :: total_hits ! total hits for a single cell (summed over materials) + integer :: total_hits ! total hits for a single domain (summed over materials) integer :: min_samples ! minimum number of samples per process integer :: remainder ! leftover samples from uneven divide #ifdef MPI @@ -140,15 +150,15 @@ contains call p % initialize() -!$omp parallel private(i, j, k, i_cell, i_material, level, found_cell, & +!$omp parallel private(i, j, k, i_domain, i_material, level, found_cell, & !$omp& indices, hits) firstprivate(p) ! Reset vectors -- this is really to get around a gfortran 4.6 bug. Ideally, ! indices and hits would just be firstprivate but 4.6 complains because they ! have allocatable components... - do i_cell = 1, size(this % cell_id) - call indices(i_cell) % clear() - call hits(i_cell) % clear() + do i_domain = 1, size(this % domain_id) + call indices(i_domain) % clear() + call hits(i_domain) % clear() end do call prn_set_stream(STREAM_VOLUME) @@ -174,32 +184,89 @@ contains call find_cell(p, found_cell) if (.not. found_cell) cycle - ! Determine if point is within desired cell - LEVEL_LOOP: do level = 1, p % n_coord - CELL_CHECK_LOOP: do i_cell = 1, size(this % cell_id) - if (cells(p % coord(level) % cell) % id == this % cell_id(i_cell)) then - - ! Determine what material this is - i_material = p % material + if (this % domain_type == FILTER_MATERIAL) then + ! ====================================================================== + ! MATERIAL VOLUME + i_material = p % material + MATERIAL_LOOP: do i_domain = 1, size(this % domain_id) + if (i_material == materials(i_domain) % id) then ! Check if we've already had a hit in this material and if so, ! simply add one - do j = 1, indices(i_cell) % size() - if (indices(i_cell) % data(j) == i_material) then - hits(i_cell) % data(j) = hits(i_cell) % data(j) + 1 - cycle CELL_CHECK_LOOP + do j = 1, indices(i_domain) % size() + if (indices(i_domain) % data(j) == i_material) then + hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1 + cycle MATERIAL_LOOP end if end do ! If we make it here, that means we haven't yet had a hit in this ! material. Add an entry to both the indices list and the hits list - call indices(i_cell) % push_back(i_material) - call hits(i_cell) % push_back(1) + call indices(i_domain) % push_back(i_material) + call hits(i_domain) % push_back(1) end if - end do CELL_CHECK_LOOP - end do LEVEL_LOOP + end do MATERIAL_LOOP + + elseif (this % domain_type == FILTER_CELL) THEN + ! ====================================================================== + ! CELL VOLUME + + do level = 1, p % n_coord + CELL_LOOP: do i_domain = 1, size(this % domain_id) + if (cells(p % coord(level) % cell) % id == this % domain_id(i_domain)) then + + ! Determine what material this is + i_material = p % material + + ! Check if we've already had a hit in this material and if so, + ! simply add one + do j = 1, indices(i_domain) % size() + if (indices(i_domain) % data(j) == i_material) then + hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1 + cycle CELL_LOOP + end if + end do + + ! If we make it here, that means we haven't yet had a hit in this + ! material. Add an entry to both the indices list and the hits list + call indices(i_domain) % push_back(i_material) + call hits(i_domain) % push_back(1) + end if + end do CELL_LOOP + end do + + elseif (this % domain_type == FILTER_UNIVERSE) then + ! ====================================================================== + ! UNIVERSE VOLUME + + do level = 1, p % n_coord + UNIVERSE_LOOP: do i_domain = 1, size(this % domain_id) + if (universes(p % coord(level) % universe) % id == & + this % domain_id(i_domain)) then + + ! Determine what material this is + i_material = p % material + + ! Check if we've already had a hit in this material and if so, + ! simply add one + do j = 1, indices(i_domain) % size() + if (indices(i_domain) % data(j) == i_material) then + hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1 + cycle UNIVERSE_LOOP + end if + end do + + ! If we make it here, that means we haven't yet had a hit in this + ! material. Add an entry to both the indices list and the hits list + call indices(i_domain) % push_back(i_material) + call hits(i_domain) % push_back(1) + end if + end do UNIVERSE_LOOP + end do + + end if end do SAMPLE_LOOP - !$omp end do +!$omp end do ! ========================================================================== ! REDUCE HITS ONTO MASTER THREAD @@ -212,14 +279,14 @@ contains !$omp do ordered schedule(static) THREAD_LOOP: do i = 1, omp_get_num_threads() !$omp ordered - do i_cell = 1, size(this % cell_id) - INDEX_LOOP: do j = 1, indices(i_cell) % size() + do i_domain = 1, size(this % domain_id) + INDEX_LOOP: do j = 1, indices(i_domain) % size() ! Check if this material has been added to the master list and if so, ! accumulate the number of hits - do k = 1, master_indices(i_cell) % size() - if (indices(i_cell) % data(j) == master_indices(i_cell) % data(k)) then - master_hits(i_cell) % data(k) = & - master_hits(i_cell) % data(k) + hits(i_cell) % data(j) + do k = 1, master_indices(i_domain) % size() + if (indices(i_domain) % data(j) == master_indices(i_domain) % data(k)) then + master_hits(i_domain) % data(k) = & + master_hits(i_domain) % data(k) + hits(i_domain) % data(j) cycle INDEX_LOOP end if end do @@ -227,8 +294,8 @@ contains ! If we made it here, this means the material hasn't yet been added to ! the master list, so add an entry to both the master indices and master ! hits lists - call master_indices(i_cell) % push_back(indices(i_cell) % data(j)) - call master_hits(i_cell) % push_back(hits(i_cell) % data(j)) + call master_indices(i_domain) % push_back(indices(i_domain) % data(j)) + call master_hits(i_domain) % push_back(hits(i_domain) % data(j)) end do INDEX_LOOP end do !$omp end ordered @@ -247,7 +314,7 @@ contains volume_sample = product(this % upper_right - this % lower_left) - do i_cell = 1, size(this % cell_id) + do i_domain = 1, size(this % domain_id) atoms(:, :) = ZERO total_hits = 0 @@ -261,9 +328,9 @@ contains call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, MPI_COMM_WORLD, & MPI_STATUS_IGNORE, mpi_err) do k = 0, n - 1 - do m = 1, master_indices(i_cell) % size() - if (data(2*k + 1) == master_indices(i_cell) % data(m)) then - master_hits(i_cell) % data(m) = master_hits(i_cell) % data(m) + & + do m = 1, master_indices(i_domain) % size() + if (data(2*k + 1) == master_indices(i_domain) % data(m)) then + master_hits(i_domain) % data(m) = master_hits(i_domain) % data(m) + & data(2*k + 2) end if end do @@ -272,12 +339,12 @@ contains end do #endif - do j = 1, master_indices(i_cell) % size() - total_hits = total_hits + master_hits(i_cell) % data(j) - f = real(master_hits(i_cell) % data(j), 8) / this % samples + do j = 1, master_indices(i_domain) % size() + total_hits = total_hits + master_hits(i_domain) % data(j) + f = real(master_hits(i_domain) % data(j), 8) / this % samples var_f = f*(ONE - f) / this % samples - i_material = master_indices(i_cell) % data(j) + i_material = master_indices(i_domain) % data(j) if (i_material == MATERIAL_VOID) cycle associate (mat => materials(i_material)) @@ -293,9 +360,9 @@ contains end do ! Determine volume - volume(1, i_cell) = real(total_hits, 8) / this % samples * volume_sample - volume(2, i_cell) = sqrt(volume(1, i_cell) * (volume_sample - & - volume(1, i_cell)) / this % samples) + volume(1, i_domain) = real(total_hits, 8) / this % samples * volume_sample + volume(2, i_domain) = sqrt(volume(1, i_domain) * (volume_sample - & + volume(1, i_domain)) / this % samples) ! Determine total number of atoms. At this point, we have values in ! atoms/b-cm. To get to atoms we multiple by 10^24 V. @@ -307,19 +374,19 @@ contains ! Convert full arrays to vectors do j = 1, size(nuclides) if (atoms(1, j) > ZERO) then - call nuclide_vec(i_cell) % push_back(j) - call atoms_vec(i_cell) % push_back(atoms(1, j)) - call uncertainty_vec(i_cell) % push_back(atoms(2, j)) + call nuclide_vec(i_domain) % push_back(j) + call atoms_vec(i_domain) % push_back(atoms(1, j)) + call uncertainty_vec(i_domain) % push_back(atoms(2, j)) end if end do else #ifdef MPI - n = master_indices(i_cell) % size() + n = master_indices(i_domain) % size() allocate(data(2*n)) do k = 0, n - 1 - data(2*k + 1) = master_indices(i_cell) % data(k + 1) - data(2*k + 2) = master_hits(i_cell) % data(k + 1) + data(2*k + 1) = master_indices(i_domain) % data(k + 1) + data(2*k + 2) = master_hits(i_domain) % data(k + 1) end do call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, MPI_COMM_WORLD, mpi_err) @@ -340,7 +407,7 @@ contains uncertainty_vec) type(VolumeCalculation), intent(in) :: this character(*), intent(in) :: filename ! filename for HDF5 file - real(8), intent(in) :: volume(:,:) ! volume mean/stdev in each cell + real(8), intent(in) :: volume(:,:) ! volume mean/stdev in each domain type(VectorInt), intent(in) :: nuclide_vec(:) ! indices in nuclides array type(VectorReal), intent(in) :: atoms_vec(:) ! total # of atoms of each nuclide type(VectorReal), intent(in) :: uncertainty_vec(:) ! uncertainty of total # of atoms @@ -357,15 +424,23 @@ contains file_id = file_create(filename) ! Write basic metadata + select case (this % domain_type) + case (FILTER_CELL) + call write_attribute_string(file_id, ".", "domain_type", "cell") + case (FILTER_MATERIAL) + call write_attribute_string(file_id, ".", "domain_type", "material") + case (FILTER_UNIVERSE) + call write_attribute_string(file_id, ".", "domain_type", "universe") + end select call write_attribute(file_id, "samples", this % samples) call write_attribute(file_id, "lower_left", this % lower_left) call write_attribute(file_id, "upper_right", this % upper_right) - do i = 1, size(this % cell_id) - group_id = create_group(file_id, "cell_" // trim(to_str(& - this % cell_id(i)))) + do i = 1, size(this % domain_id) + group_id = create_group(file_id, "domain_" // trim(to_str(& + this % domain_id(i)))) - ! Write volume for cell + ! Write volume for domain call write_dataset(group_id, "volume", volume(:, i)) ! Create array of nuclide names from the vector diff --git a/src/volume_header.F90 b/src/volume_header.F90 index 3a34f1231a..c70345f155 100644 --- a/src/volume_header.F90 +++ b/src/volume_header.F90 @@ -1,12 +1,14 @@ module volume_header - use error, only: fatal_error + use constants, only: FILTER_CELL, FILTER_MATERIAL, FILTER_UNIVERSE + use error, only: fatal_error use xml_interface implicit none type VolumeCalculation - integer, allocatable :: cell_id(:) + integer :: domain_type + integer, allocatable :: domain_id(:) real(8) :: lower_left(3) real(8) :: upper_right(3) integer :: samples @@ -20,16 +22,31 @@ contains class(VolumeCalculation), intent(out) :: this type(Node), pointer :: node_vol - integer :: num_cells + integer :: num_domains + character(10) :: temp_str + + ! Check domain type + call get_node_value(node_vol, "domain_type", temp_str) + select case (temp_str) + case ('cell') + this % domain_type = FILTER_CELL + case ('material') + this % domain_type = FILTER_MATERIAL + case ('universe') + this % domain_type = FILTER_UNIVERSE + case default + call fatal_error("Unrecognized domain type for stochastic volume & + &calculation: " // trim(temp_str)) + end select ! Read cell IDs - if (check_for_node(node_vol, "cells")) then - num_cells = get_arraysize_integer(node_vol, "cells") + if (check_for_node(node_vol, "domain_ids")) then + num_domains = get_arraysize_integer(node_vol, "domain_ids") else call fatal_error("Must specify at least one cell for a volume calculation") end if - allocate(this % cell_id(num_cells)) - call get_node_array(node_vol, "cells", this % cell_id) + allocate(this % domain_id(num_domains)) + call get_node_array(node_vol, "domain_ids", this % domain_id) ! Read lower-left and upper-right bounding coordinates call get_node_array(node_vol, "lower_left", this % lower_left) diff --git a/tests/test_volume_calc/inputs_true.dat b/tests/test_volume_calc/inputs_true.dat index 2a3f25e6d0..e146703f31 100644 --- a/tests/test_volume_calc/inputs_true.dat +++ b/tests/test_volume_calc/inputs_true.dat @@ -1 +1 @@ -b2de6ac20ca2ca38b00d61adae707d05519f6130eea25ab050220bac005cb6d23fa5812bf876307e92e90a8806f0b371500c5611ab23f80b7c073f118eb31649 \ No newline at end of file +102569289552d021b6803f404a0c17a9c17a40578fdba43a6ba08b77a731e0368fffa6a8a7abd48555167cb9997c6dba9ec5044c8593b12056957b7e3ec44ed0 \ No newline at end of file diff --git a/tests/test_volume_calc/results_true.dat b/tests/test_volume_calc/results_true.dat index e7a39a9312..da6dfd2afc 100644 --- a/tests/test_volume_calc/results_true.dat +++ b/tests/test_volume_calc/results_true.dat @@ -1,7 +1,8 @@ k-combined: 4.165451e-02 3.582531e-04 -Cell 1: 31.4693 +/- 0.0721 cm^3 -Cell 2: 2.0933 +/- 0.0310 cm^3 -Cell 3: 2.0486 +/- 0.0307 cm^3 +Volume calculation 0 +Domain 1: 31.4693 +/- 0.0721 cm^3 +Domain 2: 2.0933 +/- 0.0310 cm^3 +Domain 3: 2.0486 +/- 0.0307 cm^3 Cell Nuclide Atoms Uncertainty 0 1 U235.71c 3.481769e+23 7.979991e+20 1 1 Mo99.71c 3.481769e+22 7.979991e+19 @@ -11,3 +12,20 @@ Cell 3: 2.0486 +/- 0.0307 cm^3 5 3 H1.71c 1.369920e+23 2.051689e+21 6 3 O16.71c 6.849599e+22 1.025844e+21 7 3 B10.71c 6.849599e+18 1.025844e+17 +Volume calculation 1 +Domain 1: 4.1419 +/- 0.0426 cm^3 +Domain 2: 31.4693 +/- 0.0721 cm^3 + Material Nuclide Atoms Uncertainty +0 1 H1.71c 2.769690e+23 2.850068e+21 +1 1 O16.71c 1.384845e+23 1.425034e+21 +2 1 B10.71c 1.384845e+19 1.425034e+17 +3 2 U235.71c 3.481769e+23 7.979991e+20 +4 2 Mo99.71c 3.481769e+22 7.979991e+19 +Volume calculation 2 +Domain 0: 35.6112 +/- 0.0664 cm^3 + Universe Nuclide Atoms Uncertainty +0 0 H1.71c 2.769690e+23 2.850068e+21 +1 0 O16.71c 1.384845e+23 1.425034e+21 +2 0 B10.71c 1.384845e+19 1.425034e+17 +3 0 U235.71c 3.481769e+23 7.979991e+20 +4 0 Mo99.71c 3.481769e+22 7.979991e+19 diff --git a/tests/test_volume_calc/test_volume_calc.py b/tests/test_volume_calc/test_volume_calc.py index a6d779a25a..e2c3eba797 100644 --- a/tests/test_volume_calc/test_volume_calc.py +++ b/tests/test_volume_calc/test_volume_calc.py @@ -1,6 +1,7 @@ #!/usr/bin/env python import os +import glob import sys sys.path.insert(0, os.pardir) from testing_harness import PyAPITestHarness @@ -43,9 +44,12 @@ class VolumeTest(PyAPITestHarness): geometry.export_to_xml() # Set up stochastic volume calculation - vol_calc = openmc.VolumeCalculation( - [inside_cyl, top_hemisphere, bottom_hemisphere], - 100000) + ll, ur = openmc.Union(*[c.region for c in root.cells.values()]).bounding_box + vol_calcs = [ + openmc.VolumeCalculation(list(root.cells.values()), 100000), + openmc.VolumeCalculation([water, fuel], 100000, ll, ur), + openmc.VolumeCalculation([root], 100000, ll, ur) + ] # Define settings settings = openmc.Settings() @@ -54,7 +58,7 @@ class VolumeTest(PyAPITestHarness): settings.inactive = 0 settings.source = openmc.Source(space=openmc.stats.Box( [-1., -1., -5.], [1., 1., 5.])) - settings.volume_calculations = vol_calc + settings.volume_calculations = vol_calcs settings.export_to_xml() def _get_results(self): @@ -65,15 +69,18 @@ class VolumeTest(PyAPITestHarness): # Write out k-combined. outstr = 'k-combined: {:12.6e} {:12.6e}\n'.format(*sp.k_combined) - # Read volume calculation results - vol = openmc.VolumeCalculation.from_hdf5( - os.path.join(os.getcwd(), 'volume_1.h5')) + for i, filename in enumerate(sorted(glob.glob(os.path.join( + os.getcwd(), 'volume_*.h5')))): + outstr += 'Volume calculation {}\n'.format(i) - # Write cell volumes and total # of atoms for each nuclide - for cell_id, results in sorted(vol.results.items()): - outstr += 'Cell {0}: {1[0]:.4f} +/- {1[1]:.4f} cm^3\n'.format( - cell_id, results['volume']) - outstr += str(vol.atoms_dataframe) + '\n' + # Read volume calculation results + vol = openmc.VolumeCalculation.from_hdf5(filename) + + # Write cell volumes and total # of atoms for each nuclide + for uid, results in sorted(vol.results.items()): + outstr += 'Domain {0}: {1[0]:.4f} +/- {1[1]:.4f} cm^3\n'.format( + uid, results['volume']) + outstr += str(vol.atoms_dataframe) + '\n' return outstr From f85bf914d9013b4109fc377da736a7ce134f3cba Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 9 Aug 2016 09:44:30 -0500 Subject: [PATCH 10/12] Refactor get_volume to not use VectorInt within !$omp parallel --- src/volume_calc.F90 | 163 ++++++++++++++++++++++---------------------- 1 file changed, 83 insertions(+), 80 deletions(-) diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index 03e241a66b..22647c4da9 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -112,9 +112,10 @@ contains integer :: i_domain ! index in domain_id array integer :: i_material ! index in materials array integer :: level ! local coordinate level + integer :: n_mat(size(this % domain_id)) ! Number of materials for each domain + integer, allocatable :: indices(:,:) ! List of material indices for each domain + integer, allocatable :: hits(:,:) ! Number of hits for each material in each domain logical :: found_cell - type(VectorInt) :: indices(size(this % domain_id)) ! List of material indices - type(VectorInt) :: hits(size(this % domain_id)) ! Number of hits for each material type(Particle) :: p ! Shared variables @@ -123,10 +124,10 @@ contains type(VectorInt) :: master_hits(size(this % domain_id)) ! Variables used outside of parallel region - integer :: i_nuclide ! index in nuclides array - integer :: total_hits ! total hits for a single domain (summed over materials) + integer :: i_nuclide ! index in nuclides array + integer :: total_hits ! total hits for a single domain (summed over materials) integer :: min_samples ! minimum number of samples per process - integer :: remainder ! leftover samples from uneven divide + integer :: remainder ! leftover samples from uneven divide #ifdef MPI integer :: m ! index over materials integer :: n ! number of materials @@ -151,15 +152,12 @@ contains call p % initialize() !$omp parallel private(i, j, k, i_domain, i_material, level, found_cell, & -!$omp& indices, hits) firstprivate(p) +!$omp& indices, hits, n_mat) firstprivate(p) - ! Reset vectors -- this is really to get around a gfortran 4.6 bug. Ideally, - ! indices and hits would just be firstprivate but 4.6 complains because they - ! have allocatable components... - do i_domain = 1, size(this % domain_id) - call indices(i_domain) % clear() - call hits(i_domain) % clear() - end do + ! Create space for material indices and number of hits for each + allocate(indices(size(this % domain_id), 8)) + allocate(hits(size(this % domain_id), 8)) + n_mat(:) = 0 call prn_set_stream(STREAM_VOLUME) @@ -185,83 +183,32 @@ contains if (.not. found_cell) cycle if (this % domain_type == FILTER_MATERIAL) then - ! ====================================================================== - ! MATERIAL VOLUME - i_material = p % material - MATERIAL_LOOP: do i_domain = 1, size(this % domain_id) + do i_domain = 1, size(this % domain_id) if (i_material == materials(i_domain) % id) then - ! Check if we've already had a hit in this material and if so, - ! simply add one - do j = 1, indices(i_domain) % size() - if (indices(i_domain) % data(j) == i_material) then - hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1 - cycle MATERIAL_LOOP - end if - end do - - ! If we make it here, that means we haven't yet had a hit in this - ! material. Add an entry to both the indices list and the hits list - call indices(i_domain) % push_back(i_material) - call hits(i_domain) % push_back(1) + call check_hit(i_domain, i_material, indices, hits, n_mat) end if - end do MATERIAL_LOOP + end do elseif (this % domain_type == FILTER_CELL) THEN - ! ====================================================================== - ! CELL VOLUME - do level = 1, p % n_coord - CELL_LOOP: do i_domain = 1, size(this % domain_id) + do i_domain = 1, size(this % domain_id) if (cells(p % coord(level) % cell) % id == this % domain_id(i_domain)) then - - ! Determine what material this is i_material = p % material - - ! Check if we've already had a hit in this material and if so, - ! simply add one - do j = 1, indices(i_domain) % size() - if (indices(i_domain) % data(j) == i_material) then - hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1 - cycle CELL_LOOP - end if - end do - - ! If we make it here, that means we haven't yet had a hit in this - ! material. Add an entry to both the indices list and the hits list - call indices(i_domain) % push_back(i_material) - call hits(i_domain) % push_back(1) + call check_hit(i_domain, i_material, indices, hits, n_mat) end if - end do CELL_LOOP + end do end do elseif (this % domain_type == FILTER_UNIVERSE) then - ! ====================================================================== - ! UNIVERSE VOLUME - do level = 1, p % n_coord - UNIVERSE_LOOP: do i_domain = 1, size(this % domain_id) + do i_domain = 1, size(this % domain_id) if (universes(p % coord(level) % universe) % id == & this % domain_id(i_domain)) then - - ! Determine what material this is i_material = p % material - - ! Check if we've already had a hit in this material and if so, - ! simply add one - do j = 1, indices(i_domain) % size() - if (indices(i_domain) % data(j) == i_material) then - hits(i_domain) % data(j) = hits(i_domain) % data(j) + 1 - cycle UNIVERSE_LOOP - end if - end do - - ! If we make it here, that means we haven't yet had a hit in this - ! material. Add an entry to both the indices list and the hits list - call indices(i_domain) % push_back(i_material) - call hits(i_domain) % push_back(1) + call check_hit(i_domain, i_material, indices, hits, n_mat) end if - end do UNIVERSE_LOOP + end do end do end if @@ -280,13 +227,13 @@ contains THREAD_LOOP: do i = 1, omp_get_num_threads() !$omp ordered do i_domain = 1, size(this % domain_id) - INDEX_LOOP: do j = 1, indices(i_domain) % size() + INDEX_LOOP: do j = 1, n_mat(i_domain) ! Check if this material has been added to the master list and if so, ! accumulate the number of hits do k = 1, master_indices(i_domain) % size() - if (indices(i_domain) % data(j) == master_indices(i_domain) % data(k)) then + if (indices(i_domain, j) == master_indices(i_domain) % data(k)) then master_hits(i_domain) % data(k) = & - master_hits(i_domain) % data(k) + hits(i_domain) % data(j) + master_hits(i_domain) % data(k) + hits(i_domain, j) cycle INDEX_LOOP end if end do @@ -294,16 +241,20 @@ contains ! If we made it here, this means the material hasn't yet been added to ! the master list, so add an entry to both the master indices and master ! hits lists - call master_indices(i_domain) % push_back(indices(i_domain) % data(j)) - call master_hits(i_domain) % push_back(hits(i_domain) % data(j)) + call master_indices(i_domain) % push_back(indices(i_domain, j)) + call master_hits(i_domain) % push_back(hits(i_domain, j)) end do INDEX_LOOP end do !$omp end ordered end do THREAD_LOOP !$omp end do #else - master_indices = indices - master_hits = hits + do i_domain = 1, size(this % domain_id) + do j = 1, n_mat(i_domain) + call master_indices(i_domain) % push_back(indices(i_domain, j)) + call master_hits(i_domain) % push_back(hits(i_domain, j)) + end do + end do #endif call prn_set_stream(STREAM_TRACKING) @@ -396,6 +347,58 @@ contains end if end do + contains + + !=========================================================================== + ! CHECK_HIT is an internal subroutine that checks for whether a material has + ! already been hit for a given domain. If not, it increases the list size by + ! one (taking care of re-allocation if needed). + !=========================================================================== + + subroutine check_hit(i_domain, i_material, indices, hits, n_mat) + integer :: i_domain + integer :: i_material + integer, allocatable :: indices(:,:) + integer, allocatable :: hits(:,:) + integer :: n_mat(:) + + integer, allocatable :: temp(:,:) + logical :: already_hit + integer :: j, k, nm + + ! Check if we've already had a hit in this material and if so, + ! simply add one + already_hit = .false. + nm = n_mat(i_domain) + do j = 1, nm + if (indices(i_domain, j) == i_material) then + hits(i_domain, j) = hits(i_domain, j) + 1 + already_hit = .true. + end if + end do + + if (.not. already_hit) then + ! If we make it here, that means we haven't yet had a hit in this + ! material. First check if the indices and hits arrays are large enough + ! and if not, double them. + if (nm == size(indices, 2)) then + k = 2*size(indices, 2) + allocate(temp(size(this % domain_id), k)) + temp(:, 1:nm) = indices(:, 1:nm) + call move_alloc(FROM=temp, TO=indices) + + allocate(temp(size(this % domain_id), k)) + temp(:, 1:nm) = hits(:, 1:nm) + call move_alloc(FROM=temp, TO=indices) + end if + + ! Add an entry to both the indices list and the hits list + n_mat(i_domain) = n_mat(i_domain) + 1 + indices(i_domain, n_mat(i_domain)) = i_material + hits(i_domain, n_mat(i_domain)) = 1 + end if + end subroutine check_hit + end subroutine get_volume !=============================================================================== From cec3df43d5000f98bbd6a1a49a909876ca73db13 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 11 Aug 2016 10:26:49 -0500 Subject: [PATCH 11/12] If OpenCG cell has no surfaces, don't assign region --- openmc/opencg_compatible.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index 93a257f465..d301cf2b5e 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -689,13 +689,14 @@ def get_openmc_cell(opencg_cell): translation = np.asarray(opencg_cell.translation, dtype=np.float64) openmc_cell.translation = translation - surfaces = [] - operators = [] - for surface, halfspace in opencg_cell.surfaces.values(): - surfaces.append(get_openmc_surface(surface)) - operators.append(operator.neg if halfspace == -1 else operator.pos) - openmc_cell.region = openmc.Intersection( - *[op(s) for op, s in zip(operators, surfaces)]) + if opencg_cell.surfaces: + surfaces = [] + operators = [] + for surface, halfspace in opencg_cell.surfaces.values(): + surfaces.append(get_openmc_surface(surface)) + operators.append(operator.neg if halfspace == -1 else operator.pos) + openmc_cell.region = openmc.Intersection( + *[op(s) for op, s in zip(operators, surfaces)]) # Add the OpenMC Cell to the global collection of all OpenMC Cells OPENMC_CELLS[cell_id] = openmc_cell From 0bfa509b3c330e46bf07d12f3cc629d1c66eaaa5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 11 Aug 2016 11:37:11 -0500 Subject: [PATCH 12/12] Remove extra line per PEP8 --- openmc/statepoint.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 58a79fa219..b20ea25908 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -521,7 +521,6 @@ class StatePoint(object): if self.summary is not None: self.summary.add_volume_information(volume_calc) - def get_tally(self, scores=[], filters=[], nuclides=[], name=None, id=None, estimator=None, exact_filters=False, exact_nuclides=False, exact_scores=False):