diff --git a/CMakeLists.txt b/CMakeLists.txt index 96547ebbf6..2c625517a7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -308,7 +308,6 @@ add_library(libopenmc SHARED src/dagmc_header.F90 src/constants.F90 src/dict_header.F90 - src/eigenvalue.F90 src/error.F90 src/geometry.F90 src/geometry_header.F90 @@ -316,8 +315,6 @@ add_library(libopenmc SHARED src/initialize.F90 src/input_xml.F90 src/material_header.F90 - src/math.F90 - src/mesh_header.F90 src/message_passing.F90 src/mgxs_interface.F90 src/nuclide_header.F90 diff --git a/src/constants.F90 b/src/constants.F90 index e68f5a6c34..9061bb8004 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -7,16 +7,8 @@ module constants ! ============================================================================ ! VERSIONING NUMBERS - ! OpenMC major, minor, and release numbers - integer, parameter :: VERSION_MAJOR = 0 - integer, parameter :: VERSION_MINOR = 10 - integer, parameter :: VERSION_RELEASE = 0 - integer, parameter :: & - VERSION(3) = [VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE] - ! Version numbers for binary files - integer, parameter :: VERSION_TRACK(2) = [2, 0] - integer, parameter :: VERSION_VOLUME(2) = [1, 0] + integer, parameter :: VERSION_TRACK(2) = [2, 0] ! ============================================================================ ! ADJUSTABLE PARAMETERS @@ -24,15 +16,6 @@ module constants ! NOTE: This is the only section of the constants module that should ever be ! adjusted. Modifying constants in other sections may cause the code to fail. - ! Significance level for confidence intervals - real(8), parameter :: CONFIDENCE_LEVEL = 0.95_8 - - ! Used for surface current tallies - real(8), parameter :: TINY_BIT = 1e-8_8 - - ! Maximum number of collisions/crossings - integer, parameter :: MAX_EVENTS = 1000000 - ! Maximum number of secondary particles created integer, parameter :: MAX_SECONDARY = 1000 @@ -50,24 +33,12 @@ module constants real(8), parameter :: & PI = 3.1415926535898_8, & ! pi - MASS_NEUTRON = 1.00866491588_8, & ! mass of a neutron in amu - MASS_NEUTRON_EV = 939.5654133e6_8, & ! mass of a neutron in eV/c^2 - C_LIGHT = 2.99792458e8_8, & ! speed of light in m/s - K_BOLTZMANN = 8.6173303e-5_8, & ! Boltzmann constant in eV/K INFINITY = huge(0.0_8), & ! positive infinity ZERO = 0.0_8, & HALF = 0.5_8, & ONE = 1.0_8, & TWO = 2.0_8 - ! Electron subshell labels - character(3), parameter :: SUBSHELLS(39) = [ & - 'K ', 'L1 ', 'L2 ', 'L3 ', 'M1 ', 'M2 ', 'M3 ', 'M4 ', 'M5 ', & - 'N1 ', 'N2 ', 'N3 ', 'N4 ', 'N5 ', 'N6 ', 'N7 ', 'O1 ', 'O2 ', & - 'O3 ', 'O4 ', 'O5 ', 'O6 ', 'O7 ', 'O8 ', 'O9 ', 'P1 ', 'P2 ', & - 'P3 ', 'P4 ', 'P5 ', 'P6 ', 'P7 ', 'P8 ', 'P9 ', 'P10', 'P11', & - 'Q1 ', 'Q2 ', 'Q3 '] - ! ============================================================================ ! GEOMETRY-RELATED CONSTANTS @@ -99,55 +70,11 @@ module constants ELECTRON = 2, & POSITRON = 3 - ! Reaction types - integer, parameter :: & - TOTAL_XS = 1, ELASTIC = 2, N_NONELASTIC = 3, N_LEVEL = 4, MISC = 5, & - N_2ND = 11, N_2N = 16, N_3N = 17, N_FISSION = 18, N_F = 19, & - N_NF = 20, N_2NF = 21, N_NA = 22, N_N3A = 23, N_2NA = 24, & - N_3NA = 25, N_NP = 28, N_N2A = 29, N_2N2A = 30, N_ND = 32, & - N_NT = 33, N_N3HE = 34, N_ND2A = 35, N_NT2A = 36, N_4N = 37, & - N_3NF = 38, N_2NP = 41, N_3NP = 42, N_N2P = 44, N_NPA = 45, & - N_N1 = 51, N_N40 = 90, N_NC = 91, N_DISAPPEAR = 101, N_GAMMA = 102, & - N_P = 103, N_D = 104, N_T = 105, N_3HE = 106, N_A = 107, & - N_2A = 108, N_3A = 109, N_2P = 111, N_PA = 112, N_T2A = 113, & - N_D2A = 114, N_PD = 115, N_PT = 116, N_DA = 117, N_5N = 152, & - N_6N = 153, N_2NT = 154, N_TA = 155, N_4NP = 156, N_3ND = 157, & - N_NDA = 158, N_2NPA = 159, N_7N = 160, N_8N = 161, N_5NP = 162, & - N_6NP = 163, N_7NP = 164, N_4NA = 165, N_5NA = 166, N_6NA = 167, & - N_7NA = 168, N_4ND = 169, N_5ND = 170, N_6ND = 171, N_3NT = 172, & - N_4NT = 173, N_5NT = 174, N_6NT = 175, N_2N3HE = 176, N_3N3HE = 177, & - N_4N3HE = 178, N_3N2P = 179, N_3N3A = 180, N_3NPA = 181, N_DT = 182, & - N_NPD = 183, N_NPT = 184, N_NDT = 185, N_NP3HE = 186, N_ND3HE = 187, & - N_NT3HE = 188, N_NTA = 189, N_2N2P = 190, N_P3HE = 191, N_D3HE = 192, & - N_3HEA = 193, N_4N2P = 194, N_4N2A = 195, N_4NPA = 196, N_3P = 197, & - N_N3P = 198, N_3N2PA = 199, N_5N2P = 200, N_P0 = 600, N_PC = 649, & - N_D0 = 650, N_DC = 699, N_T0 = 700, N_TC = 749, N_3HE0 = 750, & - N_3HEC = 799, N_A0 = 800, N_AC = 849, N_2N0 = 875, N_2NC = 891, & - COHERENT = 502, INCOHERENT = 504, PHOTOELECTRIC = 522, & - PAIR_PROD_ELEC = 515, PAIR_PROD = 516, PAIR_PROD_NUC = 517 - - ! Depletion reactions - integer, parameter :: DEPLETION_RX(6) = [N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N] - ! MGXS Table Types integer, parameter :: & MGXS_ISOTROPIC = 1, & ! Isotropically Weighted Data MGXS_ANGLE = 2 ! Data by Angular Bins - ! Secondary particle emission type - integer, parameter :: & - EMISSION_PROMPT = 1, & ! Prompt emission of secondary particle - EMISSION_DELAYED = 2, & ! Delayed emission of secondary particle - EMISSION_TOTAL = 3 ! Yield represents total emission (prompt + delayed) - - ! Maximum number of partial fission reactions - integer, parameter :: PARTIAL_FISSION_MAX = 4 - - ! Temperature treatment method - integer, parameter :: & - TEMPERATURE_NEAREST = 1, & - TEMPERATURE_INTERPOLATION = 2 - ! ============================================================================ ! TALLY-RELATED CONSTANTS @@ -206,21 +133,6 @@ module constants FILTER_MATERIAL = 2, & FILTER_CELL = 3 - ! Tally surface current directions - integer, parameter :: & - OUT_LEFT = 1, & ! x min - IN_LEFT = 2, & ! x min - OUT_RIGHT = 3, & ! x max - IN_RIGHT = 4, & ! x max - OUT_BACK = 5, & ! y min - IN_BACK = 6, & ! y min - OUT_FRONT = 7, & ! y max - IN_FRONT = 8, & ! y max - OUT_BOTTOM = 9, & ! z min - IN_BOTTOM = 10, & ! z min - OUT_TOP = 11, & ! z max - IN_TOP = 12 ! z max - ! Global tally parameters integer, parameter :: N_GLOBAL_TALLIES = 4 integer, parameter :: & @@ -235,25 +147,6 @@ module constants DIFF_NUCLIDE_DENSITY = 2, & DIFF_TEMPERATURE = 3 - - ! Mgxs::get_xs enumerated types - integer(C_INT), parameter :: & - MG_GET_XS_TOTAL = 0, & - MG_GET_XS_ABSORPTION = 1, & - MG_GET_XS_INVERSE_VELOCITY = 2, & - MG_GET_XS_DECAY_RATE = 3, & - MG_GET_XS_SCATTER = 4, & - MG_GET_XS_SCATTER_MULT = 5, & - MG_GET_XS_SCATTER_FMU_MULT = 6, & - MG_GET_XS_SCATTER_FMU = 7, & - MG_GET_XS_FISSION = 8, & - MG_GET_XS_KAPPA_FISSION = 9, & - MG_GET_XS_PROMPT_NU_FISSION = 10, & - MG_GET_XS_DELAYED_NU_FISSION = 11, & - MG_GET_XS_NU_FISSION = 12, & - MG_GET_XS_CHI_PROMPT = 13, & - MG_GET_XS_CHI_DELAYED = 14 - ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 deleted file mode 100644 index 4de3ae8a67..0000000000 --- a/src/eigenvalue.F90 +++ /dev/null @@ -1,15 +0,0 @@ -module eigenvalue - - use, intrinsic :: ISO_C_BINDING - - implicit none - - interface - function openmc_get_keff(k_combined) result(err) bind(C) - import C_INT, C_DOUBLE - real(C_DOUBLE), intent(out) :: k_combined(2) - integer(C_INT) :: err - end function - end interface - -end module eigenvalue diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 94a1d61905..574d2949b3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -11,12 +11,10 @@ module input_xml #endif use hdf5_interface use material_header - use mesh_header use message_passing use mgxs_interface use nuclide_header use photon_header - use random_lcg, only: prn use settings use stl_vector, only: VectorInt, VectorReal, VectorChar use string, only: to_lower, to_str, str_to_int, & @@ -424,6 +422,11 @@ contains import C_PTR type(C_PTR) :: node_ptr end subroutine + + subroutine read_meshes(node_ptr) bind(C) + import C_PTR + type(C_PTR) :: node_ptr + end subroutine end interface ! Check if tallies.xml exists diff --git a/src/math.F90 b/src/math.F90 deleted file mode 100644 index 81856c93cf..0000000000 --- a/src/math.F90 +++ /dev/null @@ -1,20 +0,0 @@ -module math - - use, intrinsic :: ISO_C_BINDING - - implicit none - private - public :: t_percentile - - interface - - pure function t_percentile(p, df) bind(C) result(t) - use ISO_C_BINDING - implicit none - real(C_DOUBLE), value, intent(in) :: p - integer(C_INT), value, intent(in) :: df - real(C_DOUBLE) :: t - end function t_percentile - end interface - -end module math diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 deleted file mode 100644 index a0c14077a8..0000000000 --- a/src/mesh_header.F90 +++ /dev/null @@ -1,296 +0,0 @@ -module mesh_header - - use, intrinsic :: ISO_C_BINDING - - implicit none - -!=============================================================================== -! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by -! congruent squares or cubes -!=============================================================================== - - type :: RegularMesh - type(C_PTR) :: ptr - contains - procedure :: id => regular_id - procedure :: volume_frac => regular_volume_frac - procedure :: n_dimension => regular_n_dimension - procedure :: dimension => regular_dimension - procedure :: lower_left => regular_lower_left - procedure :: upper_right => regular_upper_right - procedure :: width => regular_width - - procedure :: get_bin => regular_get_bin - procedure :: get_indices => regular_get_indices - procedure :: get_bin_from_indices => regular_get_bin_from_indices - procedure :: get_indices_from_bin => regular_get_indices_from_bin - end type RegularMesh - - interface - function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value, intent(in) :: n - integer(C_INT32_T), optional, intent(out) :: index_start - integer(C_INT32_T), optional, intent(out) :: index_end - integer(C_INT) :: err - end function openmc_extend_meshes - - function openmc_get_mesh_index(id, index) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value :: id - integer(C_INT32_T), intent(out) :: index - integer(C_INT) :: err - end function openmc_get_mesh_index - - function openmc_mesh_get_id(index, id) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value :: index - integer(C_INT32_T), intent(out) :: id - integer(C_INT) :: err - end function openmc_mesh_get_id - - function openmc_mesh_set_id(index, id) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), value, intent(in) :: id - integer(C_INT) :: err - end function openmc_mesh_set_id - - function openmc_mesh_get_dimension(index, dims, n) result(err) bind(C) - import C_INT32_T, C_PTR, C_INT - integer(C_INT32_T), value, intent(in) :: index - type(C_PTR), intent(out) :: dims - integer(C_INT), intent(out) :: n - integer(C_INT) :: err - end function openmc_mesh_get_dimension - - function openmc_mesh_set_dimension(index, n, dims) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT), value, intent(in) :: n - integer(C_INT), intent(in) :: dims(n) - integer(C_INT) :: err - end function openmc_mesh_set_dimension - - function openmc_mesh_get_params(index, ll, ur, width, n) result(err) bind(C) - import C_INT32_T, C_PTR, C_INT - integer(C_INT32_T), value, intent(in) :: index - type(C_PTR), intent(out) :: ll - type(C_PTR), intent(out) :: ur - type(C_PTR), intent(out) :: width - integer(C_INT), intent(out) :: n - integer(C_INT) :: err - end function openmc_mesh_get_params - - function openmc_mesh_set_params(index, n, ll, ur, width) result(err) bind(C) - import C_INT32_T, C_INT, C_DOUBLE - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), intent(in), optional :: ll(n) - real(C_DOUBLE), intent(in), optional :: ur(n) - real(C_DOUBLE), intent(in), optional :: width(n) - integer(C_INT) :: err - end function openmc_mesh_set_params - - function mesh_id(ptr) result(id) bind(C) - import C_PTR, C_INT32_T - type(C_PTR), value :: ptr - integer(C_INT32_T) :: id - end function - - function mesh_volume_frac(ptr) result(volume_frac) bind(C) - import C_PTR, C_DOUBLE - type(C_PTR), value :: ptr - real(C_DOUBLE) :: volume_frac - end function - - function mesh_n_dimension(ptr) result(n) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT) :: n - end function - - function mesh_dimension(ptr, i) result(d) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT), value :: i - integer(C_INT) :: d - end function - - function mesh_lower_left(ptr, i) result(ll) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: i - real(C_DOUBLE) :: ll - end function - - function mesh_upper_right(ptr, i) result(ur) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: i - real(C_DOUBLE) :: ur - end function - - function mesh_width(ptr, i) result(w) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: i - real(C_DOUBLE) :: w - end function - - pure function mesh_get_bin(ptr, xyz) result(bin) bind(C) - import C_PTR, C_DOUBLE, C_INT - type(C_PTR), intent(in), value :: ptr - real(C_DOUBLE), intent(in) :: xyz(*) - integer(C_INT) :: bin - end function - - pure function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: ptr - integer(C_INT), intent(in) :: ijk(*) - integer(C_INT) :: bin - end function - - pure subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C) - import C_PTR, C_DOUBLE, C_INT, C_BOOL - type(C_PTR), intent(in), value :: ptr - real(C_DOUBLE), intent(in) :: xyz(*) - integer(C_INT), intent(out) :: ijk(*) - logical(C_BOOL), intent(out) :: in_mesh - end subroutine - - pure subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: ptr - integer(C_INT), intent(in), value :: bin - integer(C_INT), intent(out) :: ijk(*) - end subroutine - - function mesh_ptr(i) result(ptr) bind(C) - import C_INT, C_PTR - integer(C_INT), value :: i - type(C_PTR) :: ptr - end function - - subroutine read_meshes(node_ptr) bind(C) - import C_PTR - type(C_PTR) :: node_ptr - end subroutine - - function n_meshes() result(n) bind(C) - import C_INT - integer(C_INT) :: n - end function - end interface - -contains - - function meshes(i) result(m) - integer, intent(in) :: i - type(RegularMesh) :: m - - m % ptr = mesh_ptr(i) - end function - - function regular_id(this) result(id) - class(RegularMesh), intent(in) :: this - integer(C_INT32_T) :: id - id = mesh_id(this % ptr) - end function - - function regular_volume_frac(this) result(volume_frac) - class(RegularMesh), intent(in) :: this - real(C_DOUBLE) :: volume_frac - volume_frac = mesh_volume_frac(this % ptr) - end function - - function regular_n_dimension(this) result(n) - class(RegularMesh), intent(in) :: this - integer(C_INT) :: n - n = mesh_n_dimension(this % ptr) - end function - - function regular_dimension(this, i) result(d) - class(RegularMesh), intent(in) :: this - integer(C_INT), intent(in) :: i - integer(C_INT) :: d - d = mesh_dimension(this % ptr, i) - end function - - function regular_lower_left(this, i) result(ll) - class(RegularMesh), intent(in) :: this - integer(C_INT), intent(in) :: i - real(C_DOUBLE) :: ll - ll = mesh_lower_left(this % ptr, i) - end function - - function regular_upper_right(this, i) result(ur) - class(RegularMesh), intent(in) :: this - integer(C_INT), intent(in) :: i - real(C_DOUBLE) :: ur - ur = mesh_upper_right(this % ptr, i) - end function - - function regular_width(this, i) result(w) - class(RegularMesh), intent(in) :: this - integer(C_INT), intent(in) :: i - real(C_DOUBLE) :: w - w = mesh_width(this % ptr, i) - end function - -!=============================================================================== -! GET_MESH_BIN determines the tally bin for a particle in a structured mesh -!=============================================================================== - - pure subroutine regular_get_bin(this, xyz, bin) - class(RegularMesh), intent(in) :: this - real(8), intent(in) :: xyz(:) ! coordinates - integer, intent(out) :: bin ! tally bin - - bin = mesh_get_bin(this % ptr, xyz) - end subroutine - -!=============================================================================== -! GET_MESH_INDICES determines the indices of a particle in a structured mesh -!=============================================================================== - - pure subroutine regular_get_indices(this, xyz, ijk, in_mesh) - class(RegularMesh), intent(in) :: this - real(8), intent(in) :: xyz(:) ! coordinates to check - integer, intent(out) :: ijk(:) ! indices in mesh - logical, intent(out) :: in_mesh ! were given coords in mesh? - - logical(C_BOOL) :: in_mesh_ - - call mesh_get_indices(this % ptr, xyz, ijk, in_mesh_) - in_mesh = in_mesh_ - end subroutine regular_get_indices - -!=============================================================================== -! MESH_INDICES_TO_BIN maps (i), (i,j), or (i,j,k) indices to a single bin number -! for use in a TallyObject results array -!=============================================================================== - - pure function regular_get_bin_from_indices(this, ijk) result(bin) - class(RegularMesh), intent(in) :: this - integer, intent(in) :: ijk(:) - integer :: bin - - bin = mesh_get_bin_from_indices(this % ptr, ijk) - end function regular_get_bin_from_indices - -!=============================================================================== -! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to -! (i), (i,j), or (i,j,k) indices -!=============================================================================== - - pure subroutine regular_get_indices_from_bin(this, bin, ijk) - class(RegularMesh), intent(in) :: this - integer, intent(in) :: bin - integer, intent(out) :: ijk(:) - - call mesh_get_indices_from_bin(this % ptr, bin, ijk) - end subroutine regular_get_indices_from_bin - -end module mesh_header diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index b0a3f64442..57a283c650 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -31,7 +31,7 @@ module nuclide_header ! Cross sections for depletion reactions (note that these are not stored in ! macroscopic cache) - real(C_DOUBLE) :: reaction(size(DEPLETION_RX)) + real(C_DOUBLE) :: reaction(6) ! Indicies and factors needed to compute cross sections from the data tables integer(C_INT) :: index_grid ! Index on nuclide energy grid