crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Remove some dead F90 code

Browse source
This commit is contained in:
Sterling Harper 2018-12-24 15:15:36 -05:00
parent e23be47c57
commit a7fdb386ad
15 changed files with 120 additions and 1169 deletions
Download patch file Download diff file Expand all files Collapse all files

2
CMakeLists.txt
View file

@ -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

3
include/openmc/constants.h
View file

@ -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 <numeric_limits> when F90
// interop is gone
// Interpolation rules
enum class Interpolation {

13
include/openmc/output.h
View file

@ -6,6 +6,8 @@
#include <string>
#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();

14
src/cell.cpp
View file

@ -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_;}

53
src/geometry_header.F90
View file

@ -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

43
src/input_xml.F90
View file

@ -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 <lattice>
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

10
src/lattice.cpp
View file

@ -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

491
src/list_header.F90
View file

@ -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

390
src/matrix_header.F90
View file

@ -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

76
src/output.F90
View file

@ -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

103
src/output.cpp
View file

@ -7,15 +7,20 @@
#include <sstream>
#include <ctime>
#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<int>(ParticleType::neutron):
std::cout << "Neutron ";
break;
case static_cast<int>(ParticleType::photon):
std::cout << "Photon ";
break;
case static_cast<int>(ParticleType::electron):
std::cout << "Electron ";
break;
case static_cast<int>(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<int64_t> 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<int64_t> 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);

66
src/set_header.F90
View file

@ -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

1
src/tallies/tally_filter.F90
View file

@ -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

3
src/tallies/tally_filter_header.F90
View file

@ -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

21
src/tallies/tally_filter_legendre.F90
View file

@ -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