From a7fdb386ad7f214918f32be64c14a6cfaa07614c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 24 Dec 2018 15:15:36 -0500 Subject: [PATCH] Remove some dead F90 code --- CMakeLists.txt | 2 - include/openmc/constants.h | 3 + include/openmc/output.h | 13 +- src/cell.cpp | 14 +- src/geometry_header.F90 | 53 +-- src/input_xml.F90 | 43 +-- src/lattice.cpp | 10 - src/list_header.F90 | 491 -------------------------- src/matrix_header.F90 | 390 -------------------- src/output.F90 | 76 +--- src/output.cpp | 103 +++++- src/set_header.F90 | 66 ---- src/tallies/tally_filter.F90 | 1 - src/tallies/tally_filter_header.F90 | 3 + src/tallies/tally_filter_legendre.F90 | 21 -- 15 files changed, 120 insertions(+), 1169 deletions(-) delete mode 100644 src/matrix_header.F90 delete mode 100644 src/tallies/tally_filter_legendre.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 9b9f5a3c95..e29b5429d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -323,7 +323,6 @@ add_library(libopenmc SHARED src/list_header.F90 src/material_header.F90 src/math.F90 - src/matrix_header.F90 src/mesh_header.F90 src/message_passing.F90 src/mgxs_data.F90 @@ -365,7 +364,6 @@ add_library(libopenmc SHARED src/tallies/tally_filter_delayedgroup.F90 src/tallies/tally_filter_distribcell.F90 src/tallies/tally_filter_energy.F90 - src/tallies/tally_filter_legendre.F90 src/tallies/tally_filter_mesh.F90 src/tallies/tally_filter_meshsurface.F90 src/tallies/tally_filter_particle.F90 diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 7697843510..5b0fd8ec4c 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -419,8 +419,11 @@ constexpr int DIFF_DENSITY {1}; constexpr int DIFF_NUCLIDE_DENSITY {2}; constexpr int DIFF_TEMPERATURE {3}; +// Miscellaneous constexpr int C_NONE {-1}; constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE +constexpr int ERROR_INT {-2147483647}; // TODO: use when F90 + // interop is gone // Interpolation rules enum class Interpolation { diff --git a/include/openmc/output.h b/include/openmc/output.h index 6d6d6ba0f4..67bb2a5469 100644 --- a/include/openmc/output.h +++ b/include/openmc/output.h @@ -6,6 +6,8 @@ #include +#include "openmc/particle.h" + namespace openmc { //============================================================================== @@ -17,9 +19,8 @@ namespace openmc { void header(const char* msg, int level); - //============================================================================== -//! Retrieve a time stamp +//! Retrieve a time stamp. //! //! \return current time stamp (format: "yyyy-mm-dd hh:mm:ss") //============================================================================== @@ -27,7 +28,13 @@ void header(const char* msg, int level); std::string time_stamp(); //============================================================================== -//! Display plot information +//! Display the attributes of a particle. +//============================================================================== + +extern "C" void print_particle(Particle* p); + +//============================================================================== +//! Display plot information. //============================================================================== extern "C" void print_plot(); diff --git a/src/cell.cpp b/src/cell.cpp index e8b025adf9..1dec11f5dd 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -804,22 +804,20 @@ extern "C" { } } - int cell_type(Cell* c) {return c->type_;} - -#ifdef DAGMC - + #ifdef DAGMC int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed ) { - moab::EntityHandle surf = surf_xed->dagmc_ptr_->entity_by_id(2,surf_xed->id_); - moab::EntityHandle vol = cur_cell->dagmc_ptr_->entity_by_id(3,cur_cell->id_); + moab::EntityHandle surf = + surf_xed->dagmc_ptr_->entity_by_id(2,surf_xed->id_); + moab::EntityHandle vol = + cur_cell->dagmc_ptr_->entity_by_id(3,cur_cell->id_); moab::EntityHandle new_vol; cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol); return cur_cell->dagmc_ptr_->index_by_handle(new_vol); } - -#endif + #endif int32_t cell_universe(Cell* c) {return c->universe_;} diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index df28193829..953360c0f3 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -3,11 +3,9 @@ module geometry_header use, intrinsic :: ISO_C_BINDING use algorithm, only: find - use constants, only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, & - MATERIAL_VOID - use dict_header, only: DictCharInt, DictIntInt - use hdf5_interface, only: HID_T - use material_header, only: Material, materials, material_dict, n_materials + use constants, only: K_BOLTZMANN, MATERIAL_VOID + use dict_header, only: DictIntInt + use material_header, only: Material, materials use nuclide_header use sab_header use stl_vector, only: VectorReal @@ -40,12 +38,6 @@ module geometry_header integer(C_INT32_T), intent(in), value :: id end subroutine cell_set_id_c - function cell_type_c(cell_ptr) bind(C, name='cell_type') result(type) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT) :: type - end function cell_type_c - function cell_universe_c(cell_ptr) bind(C, name='cell_universe') & result(universe) import C_PTR, C_INT32_T @@ -89,34 +81,12 @@ module geometry_header real(C_DOUBLE) :: sqrtkT end function cell_sqrtkT_c - function lattice_pointer(lat_ind) bind(C) result(ptr) - import C_PTR, C_INT32_T - integer(C_INT32_T), intent(in), value :: lat_ind - type(C_PTR) :: ptr - end function lattice_pointer - - function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: lat_ptr - integer(C_INT32_T) :: id - end function lattice_id_c - subroutine extend_cells_c(n) bind(C) import C_INT32_t integer(C_INT32_T), intent(in), value :: n end subroutine extend_cells_c end interface -!=============================================================================== -! LATTICE abstract type for ordered array of universes. -!=============================================================================== - - type :: Lattice - type(C_PTR) :: ptr - contains - procedure :: id => lattice_id - end type Lattice - !=============================================================================== ! CELL defines a closed volume by its bounding surfaces !=============================================================================== @@ -132,7 +102,6 @@ module geometry_header procedure :: id => cell_id procedure :: set_id => cell_set_id - procedure :: type => cell_type procedure :: universe => cell_universe procedure :: fill => cell_fill procedure :: material_size => cell_material_size @@ -149,20 +118,12 @@ module geometry_header integer(C_INT32_T), bind(C) :: n_universes ! # of universes type(Cell), allocatable, target :: cells(:) - type(Lattice), allocatable, target :: lattices(:) ! Dictionaries which map user IDs to indices in the global arrays type(DictIntInt) :: cell_dict - type(DictIntInt) :: lattice_dict contains - function lattice_id(this) result(id) - class(Lattice), intent(in) :: this - integer(C_INT32_T) :: id - id = lattice_id_c(this % ptr) - end function lattice_id - !=============================================================================== function cell_id(this) result(id) @@ -177,12 +138,6 @@ contains call cell_set_id_c(this % ptr, id) end subroutine cell_set_id - function cell_type(this) result(type) - class(Cell), intent(in) :: this - integer(C_INT) :: type - type = cell_type_c(this % ptr) - end function cell_type - function cell_universe(this) result(universe) class(Cell), intent(in) :: this integer(C_INT32_T) :: universe @@ -301,10 +256,8 @@ contains n_universes = 0 if (allocated(cells)) deallocate(cells) - if (allocated(lattices)) deallocate(lattices) call cell_dict % clear() - call lattice_dict % clear() end subroutine free_memory_geometry diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 67c418509e..190ce52ab2 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -13,7 +13,6 @@ module input_xml use dagmc_header #endif use hdf5_interface - use list_header, only: ListChar, ListInt, ListReal use material_header use mesh_header use message_passing @@ -25,7 +24,7 @@ module input_xml use photon_header use random_lcg, only: prn use surface_header - use set_header, only: SetChar, SetInt + use set_header, only: SetChar use settings use stl_vector, only: VectorInt, VectorReal, VectorChar use string, only: to_lower, to_str, str_to_int, str_to_real, & @@ -291,8 +290,7 @@ contains subroutine read_geometry_xml() - integer :: i - integer :: n, n_rlats, n_hlats + integer :: i, n integer :: univ_id integer :: n_cells_in_univ real(8) :: phi, theta, psi @@ -300,14 +298,10 @@ contains logical :: boundary_exists character(MAX_LINE_LEN) :: filename type(Cell), pointer :: c - class(Lattice), pointer :: lat type(XMLDocument) :: doc type(XMLNode) :: root type(XMLNode) :: node_cell - type(XMLNode) :: node_lat type(XMLNode), allocatable :: node_cell_list(:) - type(XMLNode), allocatable :: node_rlat_list(:) - type(XMLNode), allocatable :: node_hlat_list(:) type(VectorInt) :: univ_ids ! List of all universe IDs type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each ! universe contains @@ -450,39 +444,6 @@ contains call read_lattices(root % ptr) - ! Get pointer to list of XML - call get_node_list(root, "lattice", node_rlat_list) - call get_node_list(root, "hex_lattice", node_hlat_list) - - ! Allocate lattices array - n_rlats = size(node_rlat_list) - n_hlats = size(node_hlat_list) - allocate(lattices(n_rlats + n_hlats)) - - RECT_LATTICES: do i = 1, n_rlats - lat => lattices(i) - lat % ptr = lattice_pointer(i - 1) - - ! Get pointer to i-th lattice - node_lat = node_rlat_list(i) - - ! Add lattice to dictionary - call lattice_dict % set(lat % id(), i) - - end do RECT_LATTICES - - HEX_LATTICES: do i = 1, n_hlats - lat => lattices(n_rlats + i) - lat % ptr = lattice_pointer(n_rlats + i - 1) - - ! Get pointer to i-th lattice - node_lat = node_hlat_list(i) - - ! Add lattice to dictionary - call lattice_dict % set(lat % id(), n_rlats + i) - - end do HEX_LATTICES - ! ========================================================================== ! SETUP UNIVERSES diff --git a/src/lattice.cpp b/src/lattice.cpp index a8966f8988..d424404390 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -888,14 +888,4 @@ read_lattices(pugi::xml_node *node) } } -//============================================================================== -// Fortran compatibility functions -//============================================================================== - -extern "C" { - Lattice* lattice_pointer(int lat_ind) {return model::lattices[lat_ind];} - - int32_t lattice_id(Lattice *lat) {return lat->id_;} -} - } // namespace openmc diff --git a/src/list_header.F90 b/src/list_header.F90 index 8020b96b17..383607834a 100644 --- a/src/list_header.F90 +++ b/src/list_header.F90 @@ -16,18 +16,6 @@ module list_header ! LISTELEM* types hold one piece of data and a pointer to the next piece of data !=============================================================================== - type :: ListElemInt - integer :: data - type(ListElemInt), pointer :: next => null() - type(ListElemInt), pointer :: prev => null() - end type ListElemInt - - type :: ListElemReal - real(8) :: data - type(ListElemReal), pointer :: next => null() - type(ListElemReal), pointer :: prev => null() - end type ListElemReal - type :: ListElemChar character(MAX_WORD_LEN) :: data type(ListElemChar), pointer :: next => null() @@ -41,50 +29,6 @@ module list_header ! code. For the time being, a separate derived type exists for each datatype. !=============================================================================== - type, public :: ListInt - private - integer :: count = 0 ! Number of elements in list - - ! Used in get_item for fast sequential lookups - integer :: last_index = huge(0) - type(ListElemInt), pointer :: last_elem => null() - - ! Pointers to beginning and end of list - type(ListElemInt), public, pointer :: head => null() - type(ListElemInt), public, pointer :: tail => null() - contains - procedure :: append => list_append_int ! Add item to end of list - procedure :: clear => list_clear_int ! Remove all items - procedure :: contains => list_contains_int ! Does list contain? - procedure :: get_item => list_get_item_int ! Get i-th item in list - procedure :: index => list_index_int ! Determine index of given item - procedure :: insert => list_insert_int ! Insert item in i-th position - procedure :: remove => list_remove_int ! Remove specified item - procedure :: size => list_size_int ! Size of list - end type ListInt - - type, public :: ListReal - private - integer :: count = 0 ! Number of elements in list - - ! Used in get_item for fast sequential lookups - integer :: last_index = huge(0) - type(ListElemReal), pointer :: last_elem => null() - - ! Pointers to beginning and end of list - type(ListElemReal), public, pointer :: head => null() - type(ListElemReal), public, pointer :: tail => null() - contains - procedure :: append => list_append_real ! Add item to end of list - procedure :: clear => list_clear_real ! Remove all items - procedure :: contains => list_contains_real ! Does list contain? - procedure :: get_item => list_get_item_real ! Get i-th item in list - procedure :: index => list_index_real ! Determine index of given item - procedure :: insert => list_insert_real ! Insert item in i-th position - procedure :: remove => list_remove_real ! Remove specified item - procedure :: size => list_size_real ! Size of list - end type ListReal - type, public :: ListChar private integer :: count = 0 ! Number of elements in list @@ -114,56 +58,6 @@ contains ! becomes the first item. !=============================================================================== - subroutine list_append_int(this, data) - class(ListInt) :: this - integer :: data - - type(ListElemInt), pointer :: elem - - ! Create element and set dat - allocate(elem) - elem % data = data - - if (.not. associated(this % head)) then - ! If list is empty, set head and tail to new element - this % head => elem - this % tail => elem - else - ! Otherwise append element at end of list - this % tail % next => elem - elem % prev => this % tail - this % tail => this % tail % next - end if - - this % count = this % count + 1 - - end subroutine list_append_int - - subroutine list_append_real(this, data) - class(ListReal) :: this - real(8) :: data - - type(ListElemReal), pointer :: elem - - ! Create element and set dat - allocate(elem) - elem % data = data - - if (.not. associated(this % head)) then - ! If list is empty, set head and tail to new element - this % head => elem - this % tail => elem - else - ! Otherwise append element at end of list - this % tail % next => elem - elem % prev => this % tail - this % tail => this % tail % next - end if - - this % count = this % count + 1 - - end subroutine list_append_real - subroutine list_append_char(this, data) class(ListChar) :: this character(*) :: data @@ -193,58 +87,6 @@ contains ! LIST_CLEAR removes all elements from the list !=============================================================================== - subroutine list_clear_int(this) - class(ListInt) :: this - - type(ListElemInt), pointer :: current => null() - type(ListElemInt), pointer :: next => null() - - if (this % count > 0) then - current => this % head - do while (associated(current)) - ! Set pointer to next element - next => current % next - - ! Deallocate memory for current element - deallocate(current) - - ! Move to next element - current => next - end do - - nullify(this % head) - nullify(this % tail) - this % count = 0 - end if - - end subroutine list_clear_int - - subroutine list_clear_real(this) - class(ListReal) :: this - - type(ListElemReal), pointer :: current => null() - type(ListElemReal), pointer :: next => null() - - if (this % count > 0) then - current => this % head - do while (associated(current)) - ! Set pointer to next element - next => current % next - - ! Deallocate memory for current element - deallocate(current) - - ! Move to next element - current => next - end do - - nullify(this % head) - nullify(this % tail) - this % count = 0 - end if - - end subroutine list_clear_real - subroutine list_clear_char(this) class(ListChar) :: this @@ -276,24 +118,6 @@ contains ! relies on the index method, it is O(n). !=============================================================================== - function list_contains_int(this, data) result(in_list) - class(ListInt) :: this - integer :: data - logical :: in_list - - in_list = (this % index(data) > 0) - - end function list_contains_int - - function list_contains_real(this, data) result(in_list) - class(ListReal) :: this - real(8) :: data - logical :: in_list - - in_list = (this % index(data) > 0) - - end function list_contains_real - function list_contains_char(this, data) result(in_list) class(ListChar) :: this character(*) :: data @@ -308,72 +132,6 @@ contains ! is out of bounds, an error code is returned. ! =============================================================================== - function list_get_item_int(this, i_list) result(data) - class(ListInt) :: this - integer :: i_list - integer :: data - - integer :: last_index - - if (i_list < 1 .or. i_list > this % count) then - ! Check for index out of bounds - data = ERROR_INT - elseif (i_list == 1) then - data = this % head % data - this % last_index = 1 - this % last_elem => this % head - elseif (i_list == this % count) then - data = this % tail % data - this % last_index = this % count - this % last_elem => this % tail - else - if (i_list < this % last_index) then - this % last_index = 1 - this % last_elem => this % head - end if - - do last_index = this % last_index + 1, i_list - this % last_elem => this % last_elem % next - this % last_index = last_index - end do - data = this % last_elem % data - end if - - end function list_get_item_int - - function list_get_item_real(this, i_list) result(data) - class(ListReal) :: this - integer :: i_list - real(8) :: data - - integer :: last_index - - if (i_list < 1 .or. i_list > this % count) then - ! Check for index out of bounds - data = ERROR_REAL - elseif (i_list == 1) then - data = this % head % data - this % last_index = 1 - this % last_elem => this % head - elseif (i_list == this % count) then - data = this % tail % data - this % last_index = this % count - this % last_elem => this % tail - else - if (i_list < this % last_index) then - this % last_index = 1 - this % last_elem => this % head - end if - - do last_index = this % last_index + 1, i_list - this % last_elem => this % last_elem % next - this % last_index = last_index - end do - data = this % last_elem % data - end if - - end function list_get_item_real - function list_get_item_char(this, i_list) result(data) class(ListChar) :: this integer :: i_list @@ -412,48 +170,6 @@ contains ! 'data' is not present in the list, the return value is -1. !=============================================================================== - function list_index_int(this, data) result(i_list) - - class(ListInt) :: this - integer :: data - integer :: i_list - - type(ListElemInt), pointer :: elem - - i_list = 0 - elem => this % head - do while (associated(elem)) - i_list = i_list + 1 - if (data == elem % data) exit - elem => elem % next - end do - - ! Check if we reached the end of the list - if (.not. associated(elem)) i_list = -1 - - end function list_index_int - - function list_index_real(this, data) result(i_list) - - class(ListReal) :: this - real(8) :: data - integer :: i_list - - type(ListElemReal), pointer :: elem - - i_list = 0 - elem => this % head - do while (associated(elem)) - i_list = i_list + 1 - if (data == elem % data) exit - elem => elem % next - end do - - ! Check if we reached the end of the list - if (.not. associated(elem)) i_list = -1 - - end function list_index_real - function list_index_char(this, data) result(i_list) class(ListChar) :: this @@ -480,117 +196,6 @@ contains ! exceeds the size of the list, the data is appends at the end of the list. !=============================================================================== - subroutine list_insert_int(this, i_list, data) - - class(ListInt) :: this - integer :: i_list - integer :: data - - integer :: i - type(ListElemInt), pointer :: elem => null() - type(ListElemInt), pointer :: new_elem => null() - - if (i_list > this % count) then - ! Check whether specified index is greater than number of elements -- if - ! so, just append it to the end of the list - call this % append(data) - - else if (i_list == 1) then - ! Check for new head element - allocate(new_elem) - new_elem % data = data - new_elem % next => this % head - this % head => new_elem - this % count = this % count + 1 - - else - ! Default case with new element somewhere in middle of list - if (i_list >= this % last_index) then - i = this % last_index - elem => this % last_elem - else - i = 0 - elem => this % head - end if - do while (associated(elem)) - i = i + 1 - if (i == i_list - 1) then - ! Allocate new element - allocate(new_elem) - new_elem % data = data - - ! Put it before the i-th element - new_elem % prev => elem % prev - new_elem % next => elem - new_elem % prev % next => new_elem - new_elem % next % prev => new_elem - this % count = this % count + 1 - this % last_index = i_list - this % last_elem => new_elem - exit - end if - i = i + 1 - elem => elem % next - end do - end if - - end subroutine list_insert_int - - subroutine list_insert_real(this, i_list, data) - - class(ListReal) :: this - integer :: i_list - real(8) :: data - - integer :: i - type(ListElemReal), pointer :: elem => null() - type(ListElemReal), pointer :: new_elem => null() - - if (i_list > this % count) then - ! Check whether specified index is greater than number of elements -- if - ! so, just append it to the end of the list - call this % append(data) - - else if (i_list == 1) then - ! Check for new head element - allocate(new_elem) - new_elem % data = data - new_elem % next => this % head - this % head => new_elem - this % count = this % count + 1 - - else - ! Default case with new element somewhere in middle of list - if (i_list >= this % last_index) then - i = this % last_index - elem => this % last_elem - else - i = 0 - elem => this % head - end if - do while (associated(elem)) - if (i == i_list) then - ! Allocate new element - allocate(new_elem) - new_elem % data = data - - ! Put it before the i-th element - new_elem % prev => elem % prev - new_elem % next => elem - new_elem % prev % next => new_elem - new_elem % next % prev => new_elem - this % count = this % count + 1 - this % last_index = i_list - this % last_elem => new_elem - exit - end if - i = i + 1 - elem => elem % next - end do - end if - - end subroutine list_insert_real - subroutine list_insert_char(this, i_list, data) class(ListChar) :: this @@ -651,84 +256,6 @@ contains ! is not in the list, no action is taken. !=============================================================================== - subroutine list_remove_int(this, data) - - class(ListInt) :: this - integer :: data - - type(ListElemInt), pointer :: elem => null() - - elem => this % head - do while (associated(elem)) - ! Check for matching data - if (elem % data == data) then - - ! Determine whether the current element is the head, tail, or a middle - ! element - if (associated(elem, this % head)) then - this % head => elem % next - if (associated(elem, this % tail)) nullify(this % tail) - if (associated(this % head)) nullify(this % head % prev) - deallocate(elem) - else if (associated(elem, this % tail)) then - this % tail => elem % prev - deallocate(this % tail % next) - else - elem % prev % next => elem % next - elem % next % prev => elem % prev - deallocate(elem) - end if - - ! Decrease count and exit - this % count = this % count - 1 - exit - end if - - ! Advance pointers - elem => elem % next - end do - - end subroutine list_remove_int - - subroutine list_remove_real(this, data) - - class(ListReal) :: this - real(8) :: data - - type(ListElemReal), pointer :: elem => null() - - elem => this % head - do while (associated(elem)) - ! Check for matching data - if (elem % data == data) then - - ! Determine whether the current element is the head, tail, or a middle - ! element - if (associated(elem, this % head)) then - this % head => elem % next - if (associated(elem, this % tail)) nullify(this % tail) - if (associated(this % head)) nullify(this % head % prev) - deallocate(elem) - else if (associated(elem, this % tail)) then - this % tail => elem % prev - deallocate(this % tail % next) - else - elem % prev % next => elem % next - elem % next % prev => elem % prev - deallocate(elem) - end if - - ! Decrease count and exit - this % count = this % count - 1 - exit - end if - - ! Advance pointers - elem => elem % next - end do - - end subroutine list_remove_real - subroutine list_remove_char(this, data) class(ListChar) :: this @@ -772,24 +299,6 @@ contains ! LIST_SIZE returns the number of elements in the list !=============================================================================== - function list_size_int(this) result(size) - - class(ListInt) :: this - integer :: size - - size = this % count - - end function list_size_int - - function list_size_real(this) result(size) - - class(ListReal) :: this - integer :: size - - size = this % count - - end function list_size_real - function list_size_char(this) result(size) class(ListChar) :: this diff --git a/src/matrix_header.F90 b/src/matrix_header.F90 deleted file mode 100644 index b2cd3b1291..0000000000 --- a/src/matrix_header.F90 +++ /dev/null @@ -1,390 +0,0 @@ -module matrix_header - - implicit none - private - - type, public :: Matrix - integer :: n ! number of rows/cols in matrix - integer :: nnz ! number of nonzeros in matrix - integer :: n_count ! counter for length of matrix - integer :: nz_count ! counter for number of non zeros - integer, allocatable :: row(:) ! csr row vector - integer, allocatable :: col(:) ! column vector - real(8), allocatable :: val(:) ! matrix value vector - contains - procedure :: create => matrix_create - procedure :: destroy => matrix_destroy - procedure :: add_value => matrix_add_value - procedure :: new_row => matrix_new_row - procedure :: assemble => matrix_assemble - procedure :: get_row => matrix_get_row - procedure :: get_col => matrix_get_col - procedure :: vector_multiply => matrix_vector_multiply - procedure :: search_indices => matrix_search_indices - procedure :: write => matrix_write - procedure :: copy => matrix_copy - procedure :: transpose => matrix_transpose - end type matrix - -contains - -!=============================================================================== -! MATRIX_CREATE allocates CSR vectors -!=============================================================================== - - subroutine matrix_create(self, n, nnz) - - integer, intent(in) :: n ! dimension of matrix - integer, intent(in) :: nnz ! number of nonzeros - class(Matrix), intent(inout) :: self ! matrix instance - - ! Preallocate vectors - if (.not.allocated(self % row)) allocate(self % row(n+1)) - if (.not.allocated(self % col)) allocate(self % col(nnz)) - if (.not.allocated(self % val)) allocate(self % val(nnz)) - - ! Set counters to 1 - self % n_count = 1 - self % nz_count = 1 - - ! Set n and nnz - self % n = n - self % nnz = nnz - - end subroutine matrix_create - -!=============================================================================== -! MATRIX_DESTROY deallocates all space associated with a matrix -!=============================================================================== - - subroutine matrix_destroy(self) - - class(Matrix), intent(inout) :: self ! matrix instance - - if (allocated(self % row)) deallocate(self % row) - if (allocated(self % col)) deallocate(self % col) - if (allocated(self % val)) deallocate(self % val) - - end subroutine matrix_destroy - -!=============================================================================== -! MATRIX_ADD_VALUE adds a value to the matrix -!=============================================================================== - - subroutine matrix_add_value(self, col, val) - - integer, intent(in) :: col ! col location in matrix - real(8), intent(in) :: val ! value to store in matrix - class(Matrix), intent(inout) :: self ! matrix instance - - ! Record the data - self % col(self % nz_count) = col - self % val(self % nz_count) = val - - ! Increment the number of nonzeros currently stored - self % nz_count = self % nz_count + 1 - - end subroutine matrix_add_value - -!=============================================================================== -! MATRIX_NEW_ROW adds a new row by saving the column position of the new row -!=============================================================================== - - subroutine matrix_new_row(self) - - class(Matrix), intent(inout) :: self ! matrix instance - - ! Record the current number of nonzeros - self % row(self % n_count) = self % nz_count - - ! Increment the current row that we are on - self % n_count = self % n_count + 1 - - end subroutine matrix_new_row - -!=============================================================================== -! MATRIX_ASSEMBLE main rountine to sort all the columns in CSR matrix -!=============================================================================== - - subroutine matrix_assemble(self) - - class(Matrix), intent(inout) :: self ! matrix instance - - integer :: i - integer :: first - integer :: last - - ! Loop around row vector - do i = 1, self % n - - ! Get bounds - first = self % row(i) - last = self % row(i+1) - 1 - - ! Sort a row - call sort_csr(self % col, self % val, first, last) - - end do - - end subroutine matrix_assemble - -!=============================================================================== -! SORT_CSR main routine that performs a sort on a CSR col/val subvector -!=============================================================================== - - recursive subroutine sort_csr(col, val, first, last) - - integer :: col(:) ! column vector to sort - integer :: first ! first value in sort - integer :: last ! last value in sort - real(8) :: val(:) ! value vector to be sorted like columns - - integer :: mid ! midpoint value - - if (first < last) then - call split(col, val, first, last, mid) - call sort_csr(col, val, first, mid-1) - call sort_csr(col, val, mid+1, last) - end if - - end subroutine sort_csr - -!=============================================================================== -! SPLIT bisects the search space for the sorting routine -!=============================================================================== - - subroutine split(col, val, low, high, mid) - - integer :: col(:) ! column vector to sort - integer :: low ! low index to sort - integer :: high ! high index to sort - integer :: mid ! middle of sort - real(8) :: val(:) ! value vector to be sorted like columns - - integer :: left ! contains left value in sort - integer :: right ! contains right value in sort - integer :: iswap ! temporary interger swap - integer :: pivot ! pivotting variable for columns - real(8) :: rswap ! temporary real swap - real(8) :: val0 ! pivot for value vector - - left = low - right = high - pivot = col(low) - val0 = val(low) - - ! Repeat the following while left and right havent met - do while (left < right) - - ! Scan right to left to find element < pivot - do while (left < right .and. col(right) >= pivot) - right = right - 1 - end do - - ! Scan left to right to find element > pivot - do while (left < right .and. col(left) <= pivot) - left = left + 1 - end do - - ! If left and right havent met, exchange the items - if (left < right) then - iswap = col(left) - col(left) = col(right) - col(right) = iswap - rswap = val(left) - val(left) = val(right) - val(right) = rswap - end if - - end do - - ! Switch the element in split position with pivot - col(low) = col(right) - col(right) = pivot - mid = right - val(low) = val(right) - val(right) = val0 - - end subroutine split - -!=============================================================================== -! MATRIX_GET_ROW is a method to get row -!=============================================================================== - - function matrix_get_row(self, i) result(row) - - class(Matrix), intent(in) :: self ! matrix instance - integer, intent(in) :: i ! row to get - integer :: row ! index in col where row begins - - row = self % row(i) - - end function matrix_get_row - -!=============================================================================== -! MATRIX_GET_COL is a method to get column -!=============================================================================== - - function matrix_get_col(self, i) result(col) - - class(Matrix), intent(in) :: self ! matrix instance - integer, intent(in) :: i ! index from row vector - integer :: col ! the actual column - - col = self % col(i) - - end function matrix_get_col - -!=============================================================================== -! MATRIX_TRANSPOSE transposes a sparse matrix -!=============================================================================== - - function matrix_transpose(self) result(mat) - - class(Matrix), intent(in) :: self ! matrix instance - type(Matrix) :: mat ! transposed matrix - - integer :: i, j ! loop indices for row/column - integer :: i_transpose ! loop index for row in transposed matrix - integer :: first, last ! indices in col() array - - call mat % create(self % n, self % nnz) - do i_transpose = 1, mat % n - ! Set up row in transposed matrix - call mat % new_row() - - ROW: do i = 1, self % n - ! Get range of columns for row i - first = self % row(i) - last = self % row(i+1) - 1 - - COL: do j = first, last - if (self % col(j) == i_transpose) then - ! If column in original matrix matches row in transposed, add value - call mat % add_value(i, self % val(j)) - elseif (self % col(j) > i_transpose) then - exit COL - end if - end do COL - end do ROW - end do - - call mat % new_row() - - end function matrix_transpose - -!=============================================================================== -! MATRIX_VECTOR_MULTIPLY allow a vector to multiply the matrix -!=============================================================================== - - subroutine matrix_vector_multiply(self, vec_in, vec_out) - - use constants, only: ZERO - use vector_header, only: Vector - - class(Matrix), intent(in) :: self ! matrix instance - type(Vector), intent(in) :: vec_in ! vector to multiply matrix against - type(Vector), intent(inout) :: vec_out ! resulting vector - - integer :: i ! row iteration counter - integer :: j ! column iteration counter - - ! Begin loop around rows - ROWS: do i = 1, self % n - - ! Initialize target location in vector - vec_out % val(i) = ZERO - - ! Begin loop around columns - COLS: do j = self % get_row(i), self % get_row(i + 1) - 1 - - vec_out % val(i) = vec_out % val(i) + self % val(j) * & - vec_in % val(self % get_col(j)) - - end do COLS - - end do ROWS - - end subroutine matrix_vector_multiply - -!=============================================================================== -! MATRIX_SEARCH_INDICES searches for an index in column corresponding to a row -!=============================================================================== - - subroutine matrix_search_indices(self, row, col, idx, found) - - class(Matrix), intent(inout) :: self - integer, intent(in) :: row - integer, intent(in) :: col - integer, intent(out) :: idx - logical, intent(out) :: found - - integer :: j - - found = .false. - - COLS: do j = self % get_row(row), self % get_row(row + 1) - 1 - - if (self % get_col(j) == col) then - idx = j - found = .true. - exit - end if - - end do COLS - - end subroutine matrix_search_indices - -!=============================================================================== -! MATRIX_WRITE writes a matrix to file -!=============================================================================== - - subroutine matrix_write(self, filename) - - character(*), intent(in) :: filename - class(Matrix), intent(inout) :: self - - integer :: unit_ - integer :: i - integer :: j - - open(newunit=unit_, file=filename) - - do i = 1, self % n - do j = self % get_row(i), self % get_row(i + 1) - 1 - write(unit_,*) i, self % get_col(j), self % val(j) - end do - end do - - close(unit_) - - end subroutine matrix_write - -!=============================================================================== -! MATRIX_COPY copies a matrix -!=============================================================================== - - subroutine matrix_copy(self, mattocopy) - - class(Matrix), intent(inout) :: self - type(Matrix), intent(in) :: mattocopy - - ! Set n and nnz - self % n_count = mattocopy % n_count - self % nz_count = mattocopy % nz_count - self % n = mattocopy % n - self % nnz = mattocopy % nnz - - ! Allocate vectors - if (.not.allocated(self % row)) allocate(self % row(self % n + 1)) - if (.not.allocated(self % col)) allocate(self % col(self % nnz)) - if (.not.allocated(self % val)) allocate(self % val(self % nnz)) - - ! Copy over data - self % row = mattocopy % row - self % col = mattocopy % col - self % val = mattocopy % val - - end subroutine matrix_copy - -end module matrix_header diff --git a/src/output.F90 b/src/output.F90 index d2130a87a9..f0b3fe6749 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -37,6 +37,11 @@ module output integer(C_INT), value :: i real(C_DOUBLE) :: h end function + + subroutine print_particle(p) bind(C) + import Particle + type(Particle), intent(in) :: p + end subroutine end interface contains @@ -207,77 +212,6 @@ contains end subroutine print_usage -!=============================================================================== -! PRINT_PARTICLE displays the attributes of a particle -!=============================================================================== - - subroutine print_particle(p) - type(Particle), intent(in) :: p - - integer :: i ! index for coordinate levels - type(Cell), pointer :: c - class(Lattice), pointer :: l - - ! display type of particle - select case (p % type) - case (NEUTRON) - write(ou,*) 'Neutron ' // to_str(p % id) - case (PHOTON) - write(ou,*) 'Photon ' // to_str(p % id) - case (ELECTRON) - write(ou,*) 'Electron ' // to_str(p % id) - case default - write(ou,*) 'Unknown Particle ' // to_str(p % id) - end select - - ! loop through each level of universes - do i = 1, p % n_coord - ! Print level - write(ou,*) ' Level ' // trim(to_str(i - 1)) - - ! Print cell for this level - if (p % coord(i) % cell /= C_NONE) then - c => cells(p % coord(i) % cell + 1) - write(ou,*) ' Cell = ' // trim(to_str(c % id())) - end if - - ! Print universe for this level - if (p % coord(i) % universe /= NONE) then - write(ou,*) ' Universe = ' & - // trim(to_str(universe_id(p % coord(i) % universe))) - end if - - ! Print information on lattice - if (p % coord(i) % lattice /= NONE) then - l => lattices(p % coord(i) % lattice) - write(ou,*) ' Lattice = ' // trim(to_str(l % id())) - write(ou,*) ' Lattice position = (' // trim(to_str(& - p % coord(i) % lattice_x)) // ',' // trim(to_str(& - p % coord(i) % lattice_y)) // ')' - end if - - ! Print local coordinates - write(ou,'(1X,A,3ES12.4)') ' xyz = ', p % coord(i) % xyz - write(ou,'(1X,A,3ES12.4)') ' uvw = ', p % coord(i) % uvw - end do - - ! Print surface - if (p % surface /= ERROR_INT) then - write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%id(), p % surface)) - end if - - ! Display weight, energy, grid index, and interpolation factor - write(ou,*) ' Weight = ' // to_str(p % wgt) - if (run_CE) then - write(ou,*) ' Energy = ' // to_str(p % E) - else - write(ou,*) ' Energy Group = ' // to_str(p % g) - end if - write(ou,*) ' Delayed Group = ' // to_str(p % delayed_group) - write(ou,*) - - end subroutine print_particle - !=============================================================================== ! PRINT_COLUMNS displays a header listing what physical values will displayed ! below them diff --git a/src/output.cpp b/src/output.cpp index 1b197a29e5..16b96d4be2 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -7,15 +7,20 @@ #include #include -#include "openmc/cell.h" -#include "openmc/geometry.h" -#include "openmc/message_passing.h" #include "openmc/capi.h" -#include "openmc/settings.h" +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/geometry.h" +#include "openmc/lattice.h" +#include "openmc/message_passing.h" #include "openmc/plot.h" +#include "openmc/settings.h" +#include "openmc/surface.h" namespace openmc { +//============================================================================== + void header(const char* msg, int level) { // Determine how many times to repeat the '=' character. @@ -40,6 +45,8 @@ header(const char* msg, int level) { } } +//============================================================================== + std::string time_stamp() { int base_year = 1990; @@ -54,12 +61,75 @@ std::string time_stamp() //============================================================================== -//=============================================================================== -// PRINT_PLOT displays selected options for plotting -//=============================================================================== +extern "C" void print_particle(Particle* p) +{ + // Display particle type and ID. + switch (p->type) { + case static_cast(ParticleType::neutron): + std::cout << "Neutron "; + break; + case static_cast(ParticleType::photon): + std::cout << "Photon "; + break; + case static_cast(ParticleType::electron): + std::cout << "Electron "; + break; + case static_cast(ParticleType::positron): + std::cout << "Positron "; + break; + default: + std::cout << "Unknown Particle "; + } + std::cout << p->id << "\n"; -void print_plot() { + // Display particle geometry hierarchy. + for (auto i = 0; i < p->n_coord; i++) { + std::cout << " Level " << i << "\n"; + if (p->coord[i].cell != C_NONE) { + const Cell& c {*model::cells[p->coord[i].cell]}; + std::cout << " Cell = " << c.id_ << "\n"; + } + + if (p->coord[i].universe != C_NONE) { + const Universe& u {*model::universes[p->coord[i].universe]}; + std::cout << " Universe = " << u.id_ << "\n"; + } + + if (p->coord[i].lattice != F90_NONE) { + const Lattice& lat {*model::lattices[p->coord[i].lattice]}; + std::cout << " Lattice = " << lat.id_ << "\n"; + std::cout << " Lattice position = (" << p->coord[i].lattice_x + << "," << p->coord[i].lattice_y << "," + << p->coord[i].lattice_z << ")\n"; + } + + std::cout << " xyz = " << p->coord[i].xyz[0] << " " + << p->coord[i].xyz[1] << " " << p->coord[i].xyz[2] << "\n"; + std::cout << " uvw = " << p->coord[i].uvw[0] << " " + << p->coord[i].uvw[1] << " " << p->coord[i].uvw[2] << "\n"; + } + + // Display miscellaneous info. + if (p->surface != ERROR_INT) { + const Surface& surf {*model::surfaces[std::abs(p->surface)-1]}; + std::cout << " Surface = " << std::copysign(surf.id_, p->surface) << "\n"; + } + std::cout << " Weight = " << p->wgt << "\n"; + if (settings::run_CE) { + std::cout << " Energy = " << p->E << "\n"; + } else { + std::cout << " Energy Group = " << p->g << "\n"; + } + std::cout << " Delayed Group = " << p->delayed_group << "\n"; + + std::cout << "\n"; +} + +//============================================================================== + +void print_plot() +{ header("PLOTTING SUMMARY", 5); for (auto pl : model::plots) { @@ -126,14 +196,17 @@ void print_plot() { } } +//============================================================================== + void -print_overlap_check() { -#ifdef OPENMC_MPI - std::vector temp(model::overlap_check_count); - int err = MPI_Reduce(temp.data(), model::overlap_check_count.data(), - model::overlap_check_count.size(), MPI_INT64_T, MPI_SUM, 0, - mpi::intracomm); -#endif +print_overlap_check() +{ + #ifdef OPENMC_MPI + std::vector temp(model::overlap_check_count); + int err = MPI_Reduce(temp.data(), model::overlap_check_count.data(), + model::overlap_check_count.size(), MPI_INT64_T, + MPI_SUM, 0, mpi::intracomm); + #endif if (mpi::master) { header("cell overlap check summary", 1); diff --git a/src/set_header.F90 b/src/set_header.F90 index 5180d5b142..a5f57c1d52 100644 --- a/src/set_header.F90 +++ b/src/set_header.F90 @@ -21,18 +21,6 @@ module set_header ! basic tasks. !=============================================================================== - type :: SetInt - private - type(ListInt) :: elements - contains - procedure :: add => set_add_int - procedure :: clear => set_clear_int - procedure :: contains => set_contains_int - procedure :: get_item => set_get_item_int - procedure :: remove => set_remove_int - procedure :: size => set_size_int - end type SetInt - type :: SetChar private type(ListChar) :: elements @@ -51,16 +39,6 @@ contains ! SET_ADD adds an item to a set if it is not already present in the set !=============================================================================== - subroutine set_add_int(this, data) - class(SetInt) :: this - integer :: data - - if (.not. this % elements % contains(data)) then - call this % elements % append(data) - end if - - end subroutine set_add_int - subroutine set_add_char(this, data) class(SetChar) :: this character(*) :: data @@ -75,13 +53,6 @@ contains ! SET_CLEAR removes all items in a set !=============================================================================== - subroutine set_clear_int(this) - class(SetInt) :: this - - call this % elements % clear() - - end subroutine set_clear_int - subroutine set_clear_char(this) class(SetChar) :: this @@ -93,15 +64,6 @@ contains ! SET_CONTAINS determines if a specified item is in a set !=============================================================================== - function set_contains_int(this, data) result(in_set) - class(SetInt) :: this - integer :: data - logical :: in_set - - in_set = this % elements % contains(data) - - end function set_contains_int - function set_contains_char(this, data) result(in_set) class(SetChar) :: this character(*) :: data @@ -115,16 +77,6 @@ contains ! SET_GET_ITEM returns the i-th item in the set !=============================================================================== - function set_get_item_int(this, i_list) result(data) - - class(SetInt) :: this - integer :: i_list - integer :: data - - data = this % elements % get_item(i_list) - - end function set_get_item_int - function set_get_item_char(this, i_list) result(data) class(SetChar) :: this @@ -140,15 +92,6 @@ contains ! no action is taken. !=============================================================================== - subroutine set_remove_int(this, data) - - class(SetInt) :: this - integer :: data - - call this % elements % remove(data) - - end subroutine set_remove_int - subroutine set_remove_char(this, data) class(SetChar) :: this @@ -162,15 +105,6 @@ contains ! SET_SIZE returns the number of elements in the set !=============================================================================== - function set_size_int(this) result(size) - - class(SetInt) :: this - integer :: size - - size = this % elements % size() - - end function set_size_int - function set_size_char(this) result(size) class(SetChar) :: this diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index b8d5d60c0f..a1fc01dab2 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -10,7 +10,6 @@ module tally_filter use tally_filter_delayedgroup use tally_filter_distribcell use tally_filter_energy - use tally_filter_legendre use tally_filter_mesh use tally_filter_meshsurface use tally_filter_particle diff --git a/src/tallies/tally_filter_header.F90 b/src/tallies/tally_filter_header.F90 index e79b2fa368..0d2ada1684 100644 --- a/src/tallies/tally_filter_header.F90 +++ b/src/tallies/tally_filter_header.F90 @@ -93,6 +93,9 @@ module tally_filter_header type, public, extends(TallyFilter) :: EnergyFunctionFilter end type + type, public, extends(TallyFilter) :: LegendreFilter + end type + type, public, extends(TallyFilter) :: MaterialFilter end type diff --git a/src/tallies/tally_filter_legendre.F90 b/src/tallies/tally_filter_legendre.F90 deleted file mode 100644 index e5eb34bc05..0000000000 --- a/src/tallies/tally_filter_legendre.F90 +++ /dev/null @@ -1,21 +0,0 @@ -module tally_filter_legendre - - use, intrinsic :: ISO_C_BINDING - - use tally_filter_header - - implicit none - - interface - function openmc_legendre_filter_set_order(index, order) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value :: index - integer(C_INT), value :: order - integer(C_INT) :: err - end function - end interface - - type, extends(TallyFilter) :: LegendreFilter - end type LegendreFilter - -end module tally_filter_legendre