From 6ef6740ebc01f97c4fdc7772ba1b6fe29789e07b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Jul 2016 11:49:25 -0500 Subject: [PATCH] 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