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

Sampling of delayed neutrons. Added many comments.

Browse source
This commit is contained in:
Paul Romano 2011-03-25 18:31:34 +00:00
parent 948f2bab02
commit 583897eb53
7 changed files with 850 additions and 616 deletions
Download patch file Download diff file Expand all files Collapse all files

15
ChangeLog
View file

@ -1,3 +1,18 @@
2011-03-25 Paul Romano <romano7@mit.edu>
* ace.f90: Added ability to read delayed neutron energy
distributions and precursor yields/constants. Created function
length_energy_dist from inner section of read_energy_dist so that
it could used for delayed neutrons as well. Added many comments
where they were previously missing.
* fileio.f90: Added many comments.
* global.f90: Added lattice types
* physics.f90: Added sampling of delayed neutrons (delayed/prompt,
precursor group, outgoing energy). Added many comments.
* types.f90: Moved all delayed neutron precursor data into one
array. Made delayed neutron energy distribution allocatable (since
there are 6 groups, each with a dist). Added many comments.
2011-03-24 Paul Romano <romano7@mit.edu>
* geometry.f90: Fixed significant bug in how quadratic surfaces

511
src/ace.f90
View file

@ -9,10 +9,10 @@ module ace
& dict_get_key, dict_delete
use endf, only: reaction_name
integer :: NXS(16)
integer :: JXS(32)
real(8), allocatable :: XSS(:)
integer :: XSS_index
integer :: NXS(16) ! Descriptors for ACE XSS tables
integer :: JXS(32) ! Pointers into ACE XSS tables
real(8), allocatable :: XSS(:) ! Cross section data
integer :: XSS_index ! current index in XSS data
private :: NXS
private :: JXS
@ -27,18 +27,19 @@ contains
subroutine read_xs()
integer :: i ! index in materials array
integer :: j ! index over isotopes in material
integer :: index ! index in xsdatas array
integer :: n ! length of key
integer :: index_continuous ! index in xs_continuous
integer :: index_thermal ! index in xs_thermal
character(10) :: key ! name of isotope, e.g. 92235.03c
character(250) :: msg ! output/error message
type(Material), pointer :: mat => null()
type(xsData), pointer :: iso => null()
type(AceContinuous), pointer :: ace_cont => null()
type(AceThermal), pointer :: ace_thermal => null()
integer :: i, j
integer :: index
character(10) :: key
character(250) :: msg
integer :: n
integer :: index_continuous
integer :: index_thermal
type(DictionaryCI), pointer :: temp_dict => null()
type(DictionaryCI), pointer :: temp_dict => null()
call dict_create(ace_dict)
@ -114,24 +115,24 @@ contains
subroutine read_ACE_continuous(index_table, index)
integer, intent(in) :: index_table
integer, intent(in) :: index
integer, intent(in) :: index_table ! index in xs_continuous array
integer, intent(in) :: index ! index in xsdatas array
integer :: in = 7 ! unit to read from
integer :: ioError ! error status for file access
integer :: words_per_line ! number of words per line (data)
integer :: lines ! number of lines (data
integer :: n ! number of data values
real(8) :: kT ! ACE table temperature
logical :: file_exists ! does ACE library exist?
logical :: found_xs ! did we find table in library?
character(7) :: readable ! is ACE library readable?
character(250) :: msg ! output/error message
character(250) :: line ! single line to read
character(32) :: words(max_words) ! words on a line
character(100) :: filename ! name of ACE library file
character(10) :: tablename ! name of cross section table
type(AceContinuous), pointer :: table => null()
integer :: i
integer :: in = 7
integer :: ioError
integer :: words_per_line
integer :: lines
integer :: n
logical :: file_exists
logical :: found_xs
character(7) :: readable
character(250) :: msg, line
character(32) :: words(max_words)
character(100) :: filename
character(10) :: tablename
real(8) :: kT
filename = xsdatas(index)%path
tablename = xsdatas(index)%id
@ -230,7 +231,8 @@ contains
type(AceContinuous), pointer :: table
integer :: NE
integer :: NE ! number of energy points for total and elastic
! cross sections
! determine number of energy points
NE = NXS(3)
@ -246,7 +248,7 @@ contains
! read data from XSS -- right now the total, absorption and
! elastic scattering are read in to these special arrays, but in
! reality, it should be necessary to only store elastic scattering
! and possible total cross-section for total material xs
! and possibly total cross-section for total material xs
! generation.
XSS_index = 1
table%energy = get_real(NE)
@ -267,30 +269,36 @@ contains
type(AceContinuous), pointer :: table
integer :: JXS2
integer :: KNU ! Location for nu data
integer :: LNU ! Type of nu data (polynomial or tabular)
integer :: NC ! Number of polynomial coefficients
integer :: NR ! Number of interpolation regions
integer :: NE ! Number of energies
integer :: length ! Length of data block to allocate
integer :: JXS2 ! location for fission nu data
integer :: JXS24 ! location for delayed neutron data
integer :: KNU ! location for nu data
integer :: LNU ! type of nu data (polynomial or tabular)
integer :: NC ! number of polynomial coefficients
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies
integer :: NPCR ! number of delayed neutron precursor groups
integer :: LED ! location of energy distribution locators
integer :: LDIS ! location of all energy distributions
integer :: LOCC ! location of energy distributions for given MT
integer :: LNW ! location of next energy distribution if multiple
integer :: LAW ! secondary energy distribution law
integer :: IDAT ! location of first energy distribution for given MT
integer :: loc ! locator
integer :: length ! length of data to allocate
integer :: length_interp_data ! length of interpolation data
integer :: JXS24
integer :: DEC1
integer :: DEC2
JXS2 = JXS(2)
JXS2 = JXS(2)
JXS24 = JXS(24)
if (JXS2 == 0) then
! =============================================================
! No prompt/total nu data is present
! NO PROMPT/TOTAL NU DATA
table % nu_t_type = NU_NONE
table % nu_p_type = NU_NONE
elseif (XSS(JXS2) > 0) then
! =============================================================
! Prompt or total nu data is present
! PROMPT OR TOTAL NU DATA
KNU = JXS2
LNU = int(XSS(KNU))
if (LNU == 1) then
@ -324,7 +332,7 @@ contains
elseif (XSS(JXS2) < 0) then
! =============================================================
! Prompt and total nu data is present -- read prompt data first
! PROMPT AND TOTAL NU DATA -- read prompt data first
KNU = JXS2 + 1
LNU = XSS(KNU)
if (LNU == 1) then
@ -383,7 +391,7 @@ contains
if (JXS24 > 0) then
! =============================================================
! Delayed nu data is present
! DELAYED NU DATA
table % nu_d_type = NU_TABULAR
KNU = JXS24
@ -399,9 +407,85 @@ contains
! read delayed nu data
XSS_index = KNU + 1
table % nu_d_data = get_real(length)
! =============================================================
! DELAYED NEUTRON ENERGY DISTRIBUTION
! Allocate space for secondary energy distribution
NPCR = NXS(8)
table % n_precursor = NPCR
allocate(table % nu_d_edist(NPCR))
LED = JXS(26)
LDIS = JXS(27)
! Loop over all delayed neutron precursor groups
do i = 1, NPCR
! find location of energy distribution data
LOCC = XSS(LED + i - 1)
LNW = XSS(LDIS + LOCC - 1)
LAW = XSS(LDIS + LOCC)
IDAT = XSS(LDIS + LOCC + 1)
NR = XSS(LDIS + LOCC + 2)
table % nu_d_edist(i) % law = LAW
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(table % nu_d_edist(i) % nbt(NR))
allocate(table % nu_d_edist(i) % int(NR))
end if
! read ENDF interpolation parameters
XSS_index = LDIS + LOCC + 3
table % nu_d_edist(i) % nbt = get_real(NR)
table % nu_d_edist(i) % int = get_real(NR)
! allocate space for law validity data
NE = XSS(LDIS + LOCC + 3 + 2*NR)
allocate(table % nu_d_edist(i) % energy(NE))
allocate(table % nu_d_edist(i) % pvalid(NE))
length_interp_data = 5 + 2*(NR + NE)
! read law validity data
XSS_index = LDIS + LOCC + 4 + 2*NR
table % nu_d_edist(i) % energy = get_real(NE)
table % nu_d_edist(i) % pvalid = get_real(NE)
! Set index to beginning of IDAT array
loc = LDIS + IDAT - 2
! determine length of energy distribution
length = length_energy_dist(loc, LAW, LOCC, length_interp_data)
! allocate secondary energy distribution array
allocate(table % nu_d_edist(i) % data(length))
! read secondary energy distribution
XSS_index = loc + 1
table % nu_d_edist(i) % data = get_real(length)
end do
! =============================================================
! DELAYED NEUTRON PRECUSOR YIELDS AND CONSTANTS
! determine length of all precursor constants/yields/interp data
length = 0
loc = JXS(25)
do i = 1, NPCR
NR = XSS(loc + length + 1)
NE = XSS(loc + length + 2 + 2*NR)
length = length + 3 + 2*NR + 2*NE
end do
! allocate space for precusor data
allocate(table % nu_d_precursor_data(length))
! read delayed neutron precursor data
XSS_index = loc
table % nu_d_precursor_data = get_real(length)
! TODO: Read secondary energy distribution
! TODO: Read precursor data
else
table % nu_d_type = NU_NONE
end if
@ -418,22 +502,22 @@ contains
type(AceContinuous), pointer :: table
type(AceReaction), pointer :: rxn => null()
integer :: LMT ! index of MT list in XSS
integer :: NMT ! Number of reactions
integer :: JXS4 ! index of Q values in XSS
integer :: JXS5 ! index of neutron multiplicities in XSS
integer :: JXS7 ! index of reactions cross-sections in XSS
integer :: LXS !
integer :: LOCA !
integer :: LXS ! location of cross-section locators
integer :: LOCA ! location of cross-section for given MT
integer :: NE ! number of energies for reaction
type(AceReaction), pointer :: rxn => null()
LMT = JXS(3)
LMT = JXS(3)
JXS4 = JXS(4)
JXS5 = JXS(5)
LXS = JXS(6)
JXS7 = JXS(7)
NMT = NXS(4)
NMT = NXS(4)
! allocate array of reactions. Add one since we need to include an
! elastic scattering channel
@ -474,22 +558,23 @@ contains
!=====================================================================
! READ_ANGULAR_DIST parses the angular distribution for each reaction
! with secondary neutrons
!=====================================================================
subroutine read_angular_dist(table)
type(AceContinuous), pointer :: table
integer :: JXS8 ! location of angular distribution locators
integer :: JXS9 ! location of angular distributions
integer :: LOCB ! location of angular distribution for given MT
integer :: NE ! number of incoming energies
integer :: NP ! number of points for cosine distribution
integer :: LC ! locator
integer :: i ! index in reactions array
integer :: j ! index over incoming energies
integer :: length ! length of data array to allocate
type(AceReaction), pointer :: rxn => null()
integer :: JXS8
integer :: JXS9
integer :: NMT
integer :: LOCB
integer :: NE
integer :: NP
integer :: LC
integer :: i, j
integer :: length
JXS8 = JXS(8)
JXS9 = JXS(9)
@ -564,29 +649,29 @@ contains
end subroutine read_angular_dist
!=====================================================================
! READ_ENERGY_DIST
! READ_ENERGY_DIST parses the secondary energy distribution for each
! reaction with seconary neutrons (except elastic scattering)
!=====================================================================
subroutine read_energy_dist(table)
type(AceContinuous), pointer :: table
integer :: LED ! location of energy distribution locators
integer :: LDIS ! location of all energy distributions
integer :: LOCC ! location of energy distributions for given MT
integer :: LNW ! location of next energy distribution if multiple
integer :: LAW ! secondary energy distribution law
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
integer :: IDAT ! location of first energy distribution for given MT
integer :: loc ! locator
integer :: length ! length of data to allocate
integer :: length_interp_data ! length of interpolation data
integer :: i, j, k, l ! indices
type(AceReaction), pointer :: rxn => null()
integer :: LOCC
integer :: LED ! location of LDLW block
integer :: LDIS ! location of DLW block
integer :: LNW ! location of next law
integer :: LAW
integer :: NR, NE
integer :: NMU
integer :: NP
integer :: NP2
integer :: NRa, NEa, NRb, NEb
integer :: IDAT
integer :: start, length, length_interp_data
integer :: i, j, k, l
LED = JXS(10)
LED = JXS(10)
LDIS = JXS(11)
! Loop over all reactions
@ -627,135 +712,163 @@ contains
rxn % edist % pvalid = get_real(NE)
! Set index to beginning of IDAT array
start = LDIS + IDAT - 2
loc = LDIS + IDAT - 2
! Determine size of LDAT array based on which secondary energy
! law it is
length = 0
select case (LAW)
case (1)
! Tabular equiprobable energy bins
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
NET = XSS(start + 3 + 2*NR + NE)
length = 3 + 2*NR + NE + 3*NET*NE
case (2)
! Discrete photon energy
length = 2
case (3)
! Level scattering
length = 2
case (4)
! Continuous tabular distribution
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
length = length + 2 + 2*NR + 2*NE
do j = 1,NE
! determine length
NP = XSS(start + length + 2)
length = length + 2 + 3*NP
! adjust location for this block
k = start + 2 + 2*NR + NE + j
XSS(k) = XSS(k) - LOCC - length_interp_data
end do
k = start + 2 + 2*NR + NE
case (5)
! General evaporation spectrum
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
NET = XSS(start + 3 + 2*NR + 2*NE)
length = 3 + 2*NR + 2*NE + NET
case (7)
! Maxwell fission spectrum
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
length = 3 + 2*NR + 2*NE
case (9)
! Evaporation spectrum
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
length = 3 + 2*NR + 2*NE
case (11)
! Watt spectrum
NRa = XSS(start + 1)
NEa = XSS(start + 2 + 2*NRa)
NRb = XSS(start + 3 + 2*(NRa+NEa))
NEb = XSS(start + 4 + 2*(NRa+NEa+NRb))
length = 5 + 2*(NRa + NEa + NRb + NEb)
case (44)
! Kalbach-Mann correlated scattering
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
length = length + 2 + 2*NR + 2*NE
do j = 1,NE
NP = XSS(start + length + 2)
length = length + 2 + 5*NP
! adjust location for this block
k = start + 2 + 2*NR + NE + j
XSS(k) = XSS(k) - LOCC - length_interp_data
end do
case (61)
! Correlated energy and angle distribution
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
length = length + 2 + 2*NR + 2*NE
do j = 1,NE
! outgoing energy distribution
NP = XSS(start + length + 2)
! adjust locators for angular distribution
do k = 1, NP
l = start + length + 2 + 3*NP + k
if (XSS(l) /= 0) XSS(l) = XSS(l) - LOCC - length_interp_data
end do
length = length + 2 + 4*NP
do k = 1, NP
! outgoing angle distribution
NP2 = XSS(start + length + 2)
length = length + 2 + 3*NP2
end do
! adjust locators for energy distribution
k = start + 2 + 2*NR + NE + j
XSS(k) = XSS(k) - LOCC - length_interp_data
end do
case (66)
! N-body phase space distribution
length = 2
case (67)
! Laboratory energy-angle law
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
NMU = XSS(start + 4 + 2*NR + 2*NE)
length = 4 + 2*(NR + NE + NMU)
end select
! determine length of energy distribution
length = length_energy_dist(loc, LAW, LOCC, length_interp_data)
! allocate secondary energy distribution array
allocate(rxn % edist % data(length))
! read secondary energy distribution
XSS_index = start + 1
XSS_index = loc + 1
rxn % edist % data = get_real(length)
end do
end subroutine read_energy_dist
!=====================================================================
! LENGTH_ENERGY_DIST determines how many values are contained in an
! LDAT energy distribution array based on the secondary energy law and
! location in XSS
!=====================================================================
function length_energy_dist(loc, law, LOCC, lid) result(length)
integer, intent(in) :: loc ! location in XSS array
integer, intent(in) :: law ! energy distribution law
integer, intent(in) :: LOCC ! location of energy distribution
integer, intent(in) :: lid ! length of interpolation data
integer :: length ! length of energy distribution (LDAT)
integer :: i,j,k
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
integer :: NP ! number of points in outgoing energy distribution
integer :: NMU ! number of points in outgoing cosine distribution
integer :: NRa ! number of interpolation regions for Watt 'a'
integer :: NEa ! number of energies for Watt 'a'
integer :: NRb ! number of interpolation regions for Watt 'b'
integer :: NEb ! number of energies for Watt 'b'
! initialize length
length = 0
select case (law)
case (1)
! Tabular equiprobable energy bins
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NP = XSS(loc + 3 + 2*NR + NE)
length = 3 + 2*NR + NE + 3*NP*NE
case (2)
! Discrete photon energy
length = 2
case (3)
! Level scattering
length = 2
case (4)
! Continuous tabular distribution
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! determine length
NP = XSS(loc + length + 2)
length = length + 2 + 3*NP
! adjust location for this block
j = loc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
case (5)
! General evaporation spectrum
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NP = XSS(loc + 3 + 2*NR + 2*NE)
length = 3 + 2*NR + 2*NE + NP
case (7)
! Maxwell fission spectrum
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
length = 3 + 2*NR + 2*NE
case (9)
! Evaporation spectrum
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
length = 3 + 2*NR + 2*NE
case (11)
! Watt spectrum
NRa = XSS(loc + 1)
NEa = XSS(loc + 2 + 2*NRa)
NRb = XSS(loc + 3 + 2*(NRa+NEa))
NEb = XSS(loc + 4 + 2*(NRa+NEa+NRb))
length = 5 + 2*(NRa + NEa + NRb + NEb)
case (44)
! Kalbach-Mann correlated scattering
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
NP = XSS(loc + length + 2)
length = length + 2 + 5*NP
! adjust location for this block
j = loc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
case (61)
! Correlated energy and angle distribution
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! outgoing energy distribution
NP = XSS(loc + length + 2)
! adjust locators for angular distribution
do j = 1, NP
k = loc + length + 2 + 3*NP + j
if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid
end do
length = length + 2 + 4*NP
do j = 1, NP
! outgoing angle distribution -- NMU here is actually
! referred to as NP in the MCNP documentation
NMU = XSS(loc + length + 2)
length = length + 2 + 3*NMU
end do
! adjust locators for energy distribution
j = loc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
case (66)
! N-body phase space distribution
length = 2
case (67)
! Laboratory energy-angle law
NR = XSS(loc + 1)
NE = XSS(loc + 2 + 2*NR)
NMU = XSS(loc + 4 + 2*NR + 2*NE)
length = 4 + 2*(NR + NE + NMU)
end select
end function length_energy_dist
!=====================================================================
! GET_INT returns an array of integers read from the current position
! in the XSS array
@ -763,8 +876,8 @@ contains
function get_int(n_values) result(array)
integer, intent(in) :: n_values
integer :: array(n_values)
integer, intent(in) :: n_values ! number of values to read
integer :: array(n_values) ! array of values
array = int(XSS(XSS_index:XSS_index + n_values - 1))
XSS_index = XSS_index + n_values
@ -778,8 +891,8 @@ contains
function get_real(n_values) result(array)
integer, intent(in) :: n_values
real(8) :: array(n_values)
integer, intent(in) :: n_values ! number of values to read
real(8) :: array(n_values) ! array of values
array = XSS(XSS_index:XSS_index + n_values - 1)
XSS_index = XSS_index + n_values

2
src/endf.f90
View file

@ -18,6 +18,8 @@ contains
string = '(n,total)'
case (2)
string = '(elastic)'
case (4)
string = '(n,level)'
case (11)
string = '(n,2nd)'
case (16)

404
src/fileio.f90
View file

@ -36,13 +36,13 @@ contains
subroutine read_command_line()
integer :: argc
character(250) :: msg
logical :: file_exists
character(7) :: readable
integer :: argc ! number of command line arguments
logical :: file_exists ! does specified input file exist?
character(250) :: msg ! error message
character(7) :: readable ! is input file readable?
argc = COMMAND_ARGUMENT_COUNT()
if ( argc > 0 ) then
if (argc > 0) then
call GET_COMMAND_ARGUMENT(1, path_input)
else
msg = "No input file specified!"
@ -75,16 +75,17 @@ contains
character(*), intent(in) :: filename
character(250) :: line, msg
character(32) :: words(max_words)
integer :: in = 7
integer :: ioError
integer :: n
integer :: count
integer :: universe_num
type(ListKeyValueII), pointer :: key_list
type(Universe), pointer :: univ
integer :: index
integer :: in = 7 ! unit # for input file
integer :: index ! index in universes array
integer :: ioError ! error status for file access
integer :: n ! number of words on a line
integer :: count ! number of cells in a universe
integer :: universe_num ! user-specified universe #
character(250) :: line ! a line of words in input file
character(250) :: msg ! output/error message
character(32) :: words(max_words) ! words on a line
type(ListKeyValueII), pointer :: key_list => null()
type(Universe), pointer :: univ => null()
msg = "First pass through input file..."
call message(msg, 5)
@ -206,17 +207,18 @@ contains
character(*), intent(in) :: filename
character(250) :: line, msg
character(32) :: words(max_words)
integer :: in = 7
integer :: ioError
integer :: n, i
integer :: index_cell
integer :: index_universe
integer :: index_surface
integer :: index_lattice
integer :: index_material
integer :: index_source
integer :: in = 7 ! unit # for input file
integer :: ioError ! error status for file access
integer :: n ! number of words on a line
integer :: index_cell ! index in cells array
integer :: index_universe ! index in universes array
integer :: index_surface ! index in surfaces array
integer :: index_lattice ! index in lattices array
integer :: index_material ! index in materials array
integer :: index_source ! index in source array (?)
character(250) :: line ! a line of words
character(250) :: msg ! output/error message
character(32) :: words(max_words) ! words on a single line
msg = "Second pass through input file..."
call message(msg, 5)
@ -237,7 +239,7 @@ contains
do
call get_next_line(in, words, n, ioError)
if ( ioError /= 0 ) exit
if (ioError /= 0) exit
! skip comments
if (n==0 .or. words(1)(1:1) == '#') cycle
@ -294,7 +296,7 @@ contains
end subroutine read_input
!=====================================================================
! ADJUST_INDICES changes the boundary_list values for each cell and
! ADJUST_INDICES changes the values for 'surfaces' for each cell and
! the material index assigned to each to the indices in the surfaces
! and material array rather than the unique IDs assigned to each
! surface and material
@ -302,12 +304,12 @@ contains
subroutine adjust_indices()
type(Cell), pointer :: c => null()
integer :: i, j
integer :: index
integer :: surf_num
character(250) :: msg
integer :: i ! index in cells array
integer :: j ! index over surface list
integer :: index ! index in surfaces/materials array
integer :: surf_num ! user-specified surface number
character(250) :: msg ! output/error message
type(Cell), pointer :: c => null() ! cell pointer
do i = 1, n_cells
! adjust boundary list
@ -327,9 +329,7 @@ contains
end do
! adjust material indices
if (c % fill /= 0) then
! cell is filled with universe -- do nothing
else
if (c % fill == 0) then
index = dict_get_key(material_dict, c%material)
if (index == DICT_NULL) then
msg = "Could not find material " // trim(int_to_str(c%material)) // &
@ -349,12 +349,12 @@ contains
recursive subroutine build_universe(univ, parent, level)
type(Universe), pointer :: univ
integer, intent(in) :: parent
integer, intent(in) :: level
type(Universe), pointer :: univ ! univese pointer
integer, intent(in) :: parent ! cell containing universe
integer, intent(in) :: level ! level of universe
integer :: i
integer :: index
integer :: i ! index for cells in universe
integer :: index ! index in cells array
type(Cell), pointer :: c => null()
type(Universe), pointer :: subuniverse => null()
@ -367,8 +367,8 @@ contains
c => cells(index)
c%parent = parent
! if this cell is filled with another universe, recursive call
! this subroutine
! if this cell is filled with another universe, recursively
! call this subroutine
if (c % fill > 0) then
subuniverse => universes(c % fill)
call build_universe(subuniverse, index, level + 1)
@ -383,16 +383,16 @@ contains
subroutine read_cell(index, words, n_words)
integer, intent(in) :: index
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
integer, intent(in) :: index ! index in cells array
character(*), intent(in) :: words(n_words) ! words on cell card entry
integer, intent(in) :: n_words ! number of words
integer :: ioError
integer :: i
integer :: universe_num
integer :: n_surfaces
character(250) :: msg
character(32) :: word
integer :: ioError ! error status for file access
integer :: i ! index for surface list in a cell
integer :: universe_num ! user-specified universe number
integer :: n_surfaces ! number of surfaces in a cell
character(250) :: msg ! output/error message
character(32) :: word ! single word
type(Cell), pointer :: c => null()
c => cells(index)
@ -417,7 +417,7 @@ contains
if (trim(words(4)) == 'fill') then
c % type = CELL_FILL
c % material = 0
n_surfaces = n_words - 5
n_surfaces = n_words - 5
! find universe
universe_num = str_to_int(words(5))
@ -471,69 +471,69 @@ contains
subroutine read_surface(index, words, n_words)
integer, intent(in) :: index
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
integer, intent(in) :: index ! index in surfaces array
character(*), intent(in) :: words(n_words) ! words in surface card entry
integer, intent(in) :: n_words ! number of words
integer :: ioError
integer :: i
integer :: coeffs_reqd
character(250) :: msg
character(32) :: word
type(Surface), pointer :: this_surface => null()
integer :: ioError ! error status for file access
integer :: i ! index for surface coefficients
integer :: coeffs_reqd ! number of coefficients are required
character(250) :: msg ! output/error message
character(32) :: word ! single word
type(Surface), pointer :: surf => null()
this_surface => surfaces(index)
surf => surfaces(index)
! Read surface identifier
read(words(2), FMT='(I8)', IOSTAT=ioError) this_surface%uid
read(words(2), FMT='(I8)', IOSTAT=ioError) surf % uid
if (ioError > 0) then
msg = "Invalid surface name: " // words(2)
call error(msg)
end if
call dict_add_key(surface_dict, this_surface%uid, index)
call dict_add_key(surface_dict, surf % uid, index)
! Read surface type
word = words(3)
call lower_case(word)
select case (trim(word))
case ('px')
this_surface%type = SURF_PX
surf % type = SURF_PX
coeffs_reqd = 1
case ('py')
this_surface%type = SURF_PY
surf % type = SURF_PY
coeffs_reqd = 1
case ('pz')
this_surface%type = SURF_PZ
surf % type = SURF_PZ
coeffs_reqd = 1
case ('plane')
this_surface%type = SURF_PLANE
surf % type = SURF_PLANE
coeffs_reqd = 4
case ('cylx')
this_surface%type = SURF_CYL_X
surf % type = SURF_CYL_X
coeffs_reqd = 3
case ('cyly')
this_surface%type = SURF_CYL_Y
surf % type = SURF_CYL_Y
coeffs_reqd = 3
case ('cylz')
this_surface%type = SURF_CYL_Z
surf % type = SURF_CYL_Z
coeffs_reqd = 3
case ('sph')
this_surface%type = SURF_SPHERE
surf % type = SURF_SPHERE
coeffs_reqd = 4
case ('boxx')
this_surface%type = SURF_BOX_X
surf % type = SURF_BOX_X
coeffs_reqd = 4
case ('boxy')
this_surface%type = SURF_BOX_Y
surf % type = SURF_BOX_Y
coeffs_reqd = 4
case ('boxz')
this_surface%type = SURF_BOX_Z
surf % type = SURF_BOX_Z
coeffs_reqd = 4
case ('box')
this_surface%type = SURF_BOX
surf % type = SURF_BOX
coeffs_reqd = 6
case ('gq')
this_surface%type = SURF_GQ
surf % type = SURF_GQ
coeffs_reqd = 10
case default
msg = "Invalid surface type: " // words(3)
@ -547,9 +547,9 @@ contains
end if
! Read list of surfaces
allocate(this_surface%coeffs(n_words-3))
allocate(surf%coeffs(n_words-3))
do i = 1, n_words-3
this_surface%coeffs(i) = str_to_real(words(i+3))
surf % coeffs(i) = str_to_real(words(i+3))
end do
end subroutine read_surface
@ -560,9 +560,9 @@ contains
subroutine read_lattice(index, words, n_words)
integer, intent(in) :: index
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
integer, intent(in) :: index
character(*), intent(in) :: words(n_words)
integer, intent(in) :: n_words
end subroutine read_lattice
@ -572,21 +572,21 @@ contains
subroutine read_source(words, n_words)
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
character(*), intent(in) :: words(n_words) ! words on source entry
integer, intent(in) :: n_words ! number of words
character(250) :: msg
character(32) :: word
integer :: ioError
integer :: values_reqd
integer :: i
integer :: i ! index in values list
integer :: ioError ! error status for file access
integer :: values_reqd ! number of values required to specify source
character(250) :: msg ! output/error message
character(32) :: word ! single word
! Read source type
word = words(2)
call lower_case(word)
select case (trim(word))
case ('box')
external_source%type = SRC_BOX
external_source % type = SRC_BOX
values_reqd = 6
case default
msg = "Invalid source type: " // words(2)
@ -599,10 +599,10 @@ contains
call error(msg)
end if
! Read list of surfaces
! Read values
allocate(external_source%values(n_words-2))
do i = 1, n_words-2
external_source%values(i) = str_to_real(words(i+2))
external_source % values(i) = str_to_real(words(i+2))
end do
end subroutine read_source
@ -614,27 +614,29 @@ contains
subroutine read_material(index, words, n_words)
integer, intent(in) :: index
character(*), intent(in) :: words(n_words)
integer, intent(in) :: n_words
integer, intent(in) :: index ! index in materials array
character(*), intent(in) :: words(n_words) ! words on material entry
integer, intent(in) :: n_words ! number of words
character(100) :: msg
integer :: i
integer :: ioError
integer :: n_isotopes
integer :: i ! index over isotopes
integer :: ioError ! error status for file access
integer :: n_isotopes ! number of isotopes in material
character(250) :: msg ! output/error message
type(Material), pointer :: mat => null()
! Check for correct number of arguments
if (mod(n_words,2) == 0 .or. n_words < 5) then
msg = "Invalid number of arguments for material: " // words(2)
call error(msg)
end if
! Determine and set number of isotopes
n_isotopes = (n_words-3)/2
mat => materials(index)
mat%n_isotopes = n_isotopes
mat % n_isotopes = n_isotopes
! Read surface identifier
read(words(2), FMT='(I8)', IOSTAT=ioError) mat%uid
read(words(2), FMT='(I8)', IOSTAT=ioError) mat % uid
if (ioError > 0) then
msg = "Invalid surface name: " // words(2)
call error(msg)
@ -642,7 +644,7 @@ contains
call dict_add_key(material_dict, mat%uid, index)
! Read atom density
mat%atom_density = str_to_real(words(3))
mat % atom_density = str_to_real(words(3))
! allocate isotope and density list
allocate(mat%names(n_isotopes))
@ -652,8 +654,8 @@ contains
! read names and percentage
do i = 1, n_isotopes
mat%names(i) = words(2*i+2)
mat%atom_percent(i) = str_to_real(words(2*i+3))
mat % names(i) = words(2*i+2)
mat % atom_percent(i) = str_to_real(words(2*i+3))
end do
end subroutine read_material
@ -665,18 +667,19 @@ contains
subroutine normalize_ao()
integer :: index ! index used for several purposes
integer :: i ! index in materials array
integer :: j ! index over isotopes in material
real(8) :: sum_percent !
real(8) :: awr ! atomic weight ratio
real(8) :: w ! weight percent
real(8) :: x ! atom percent
logical :: percent_in_atom ! isotopes specified in atom percent?
logical :: density_in_atom ! density specified in atom/b-cm?
character(10) :: key ! name of isotopes, e.g. 92235.03c
character(100) :: msg ! output/error message
type(xsData), pointer :: iso => null()
type(Material), pointer :: mat => null()
integer :: index
integer :: i, j
real(8) :: sum_percent
real(8) :: awr ! atomic weight ratio
real(8) :: w ! weight percent
real(8) :: x ! atom percent
logical :: percent_in_atom
logical :: density_in_atom
character(10) :: key
character(100) :: msg
! first find the index in the xsdata array for each isotope in
! each material
@ -685,31 +688,31 @@ contains
! Check to make sure either all atom percents or all weight
! percents are given
if (.not. (all(mat%atom_percent > 0) .or. &
& all(mat%atom_percent < 0))) then
if (.not. (all(mat%atom_percent > ZERO) .or. &
& all(mat%atom_percent < ZERO))) then
msg = "Cannot mix atom and weight percents in material " // &
& int_to_str(mat%uid)
call error(msg)
end if
percent_in_atom = (mat%atom_percent(1) > 0)
density_in_atom = (mat%atom_density > 0)
percent_in_atom = (mat%atom_percent(1) > ZERO)
density_in_atom = (mat%atom_density > ZERO)
sum_percent = 0
do j = 1, mat%n_isotopes
sum_percent = ZERO
do j = 1, mat % n_isotopes
! Set indices for isotopes
key = mat%names(j)
key = mat % names(j)
index = dict_get_key(xsdata_dict, key)
mat%isotopes(j) = index
mat % isotopes(j) = index
! determine atomic weight ratio
awr = xsdatas(index)%awr
awr = xsdatas(index) % awr
! if given weight percent, convert all values so that they
! are divided by awr. thus, when a sum is done over the
! values, it's actually sum(w/awr)
if (.not. percent_in_atom) then
mat%atom_percent(j) = -mat%atom_percent(j) / awr
mat % atom_percent(j) = -mat % atom_percent(j) / awr
end if
end do
@ -717,20 +720,20 @@ contains
! this is straightforward. if given weight percents, the value
! is w/awr and is divided by sum(w/awr)
sum_percent = sum(mat%atom_percent)
mat%atom_percent = mat%atom_percent / sum_percent
mat % atom_percent = mat % atom_percent / sum_percent
! Change density in g/cm^3 to atom/b-cm. Since all values are
! now in atom percent, the sum needs to be re-evaluated as
! 1/sum(x*awr)
if (.not. density_in_atom) then
sum_percent = 0
do j = 1, mat%n_isotopes
index = mat%isotopes(j)
awr = xsdatas(index)%awr
x = mat%atom_percent(j)
sum_percent = ZERO
do j = 1, mat % n_isotopes
index = mat % isotopes(j)
awr = xsdatas(index) % awr
x = mat % atom_percent(j)
sum_percent = sum_percent + x*awr
end do
sum_percent = 1.0_8 / sum_percent
sum_percent = ONE / sum_percent
mat%atom_density = -mat%atom_density * N_AVOGADRO &
& / MASS_NEUTRON * sum_percent
end if
@ -743,10 +746,10 @@ contains
! specifies what cross section library should be used by default
!=====================================================================
subroutine read_xs_library( words, n_words )
subroutine read_xs_library(words, n_words)
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
character(*), intent(in) :: words(n_words)
integer, intent(in) :: n_words
end subroutine read_xs_library
@ -757,34 +760,33 @@ contains
! cycle.
!=====================================================================
subroutine read_criticality( words, n_words )
subroutine read_criticality(words, n_words)
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
character(*), intent(in) :: words(n_words) ! words on criticality card
integer, intent(in) :: n_words ! number of words
character(250) :: msg
integer :: ioError
character(250) :: msg ! output/error message
! Set problem type to criticality
problem_type = PROB_CRITICALITY
! Read number of cycles
read(words(2), FMT='(I8)', IOSTAT=ioError) n_cycles
if (ioError > 0) then
n_cycles = str_to_int(words(2))
if (n_cycles == ERROR_CODE) then
msg = "Invalid number of cycles: " // words(2)
call error(msg)
end if
! Read number of inactive cycles
read(words(3), FMT='(I8)', IOSTAT=ioError) n_inactive
if (ioError > 0) then
n_inactive = str_to_int(words(3))
if (n_inactive == ERROR_CODE) then
msg = "Invalid number of inactive cycles: " // words(2)
call error(msg)
end if
! Read number of particles
read(words(4), FMT='(I8)', IOSTAT=ioError) n_particles
if (ioError > 0) then
n_particles = str_to_int(words(4))
if (n_particles == ERROR_CODE) then
msg = "Invalid number of particles: " // words(2)
call error(msg)
end if
@ -797,9 +799,9 @@ contains
subroutine read_line(unit, line, ioError)
integer, intent(in) :: unit
character(max_line), intent(out) :: line
integer, intent(out) :: ioError
integer, intent(in) :: unit ! unit to read from
character(max_line), intent(out) :: line ! line to return
integer, intent(out) :: ioError ! error status
read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) line
@ -814,30 +816,41 @@ contains
subroutine get_next_line(unit, words, n, ioError)
integer, intent(in) :: unit
character(*), intent(out) :: words(max_words)
integer, intent(out) :: n
integer, intent(out) :: ioError
integer, intent(in) :: unit ! unit to read from
character(*), intent(out) :: words(max_words) ! words read
integer, intent(out) :: n ! number of words
integer, intent(out) :: ioError ! error status
character(250) :: line
character(32) :: local_words(max_words)
integer :: index
character(250) :: line ! single line
character(32) :: local_words(max_words) ! words on one line
integer :: index ! index of words
index = 0
do
do
! read line from file
read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) line
! if we're at the end of the file, return
if (ioError /= 0) return
! split a single line into words
call split_string_wl(line, local_words, n)
! if there are no words, we're done
if (n == 0) exit
! Check whether there is a continuation line
if (local_words(n) == '&') then
words(index+1:index+n-1) = local_words(1:n-1)
index = index + n-1
index = index + n - 1
else
words(index+1:index+n) = local_words(1:n)
index = index + n
exit
end if
end do
! set total number of words
n = index
end subroutine get_next_line
@ -848,12 +861,12 @@ contains
subroutine skip_lines(unit, n_lines, ioError)
integer, intent(in) :: unit
integer, intent(in) :: n_lines
integer, intent(out) :: ioError
integer, intent(in) :: unit ! unit to read from
integer, intent(in) :: n_lines ! number of lines to skip
integer, intent(out) :: ioError ! error status
integer :: i
character(max_line) :: tmp
integer :: i ! index for number of lines
character(max_line) :: tmp ! single line
do i = 1, n_lines
read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) tmp
@ -862,63 +875,56 @@ contains
end subroutine skip_lines
!=====================================================================
! READ_DATA_INT reads integer data into an arrya from a file open
! READ_DATA_INT reads integer data into an array from a file open
!=====================================================================
subroutine read_data_int(in, array, n, lines, words_per_line)
subroutine read_data_int(unit, array, n, lines, words_per_line)
integer, intent(in) :: in
integer, intent(out) :: array(n)
integer, intent(in) :: n
integer, intent(in) :: lines
integer, intent(in) :: words_per_line
integer, intent(in) :: unit ! unit to read from
integer, intent(out) :: array(n) ! ints read from file
integer, intent(in) :: n ! total number of ints
integer, intent(in) :: lines ! total number of lines
integer, intent(in) :: words_per_line ! number of words per line
integer :: i, j
character(250) :: line
character(32) :: words(max_words)
integer :: ioError
integer :: n_words
integer :: index
integer :: n_last
integer :: i ! line index
integer :: loc ! locator for array
integer :: ioError ! error status
index = 1
loc = 0
do i = 1, lines
if (i == lines) then
read(in,*) array(index:n)
read(UNIT=unit,FMT=*) array(loc+1:n)
else
read(in,*) array(index:index+words_per_line-1)
index = index + words_per_line
read(UNIT=unit,FMT=*) array(loc+1:loc+words_per_line)
loc = loc + words_per_line
end if
end do
end subroutine read_data_int
!=====================================================================
! READ_DATA_REAL reads real(8) data into an arrya from a file open
! READ_DATA_REAL reads real(8) data into an array from a file open
!=====================================================================
subroutine read_data_real(in, array, n, lines, words_per_line)
subroutine read_data_real(unit, array, n, lines, words_per_line)
integer, intent(in) :: in
real(8), intent(out) :: array(n)
integer, intent(in) :: n
integer, intent(in) :: lines
integer, intent(in) :: words_per_line
integer, intent(in) :: unit ! unit to read from
real(8), intent(out) :: array(n) ! real(8)s read from file
integer, intent(in) :: n ! total number of ints
integer, intent(in) :: lines ! total number of lines
integer, intent(in) :: words_per_line ! number of words per line
integer :: i, j
character(250) :: line
character(32) :: words(max_words)
integer :: ioError
integer :: n_words
integer :: index
integer :: i ! line index
integer :: loc ! locator for array
integer :: ioError ! error status
index = 1
loc = 0
do i = 1, lines
if (i == lines) then
read(in,*) array(index:n)
read(UNIT=unit,FMT=*) array(loc+1:n)
else
read(in,*) array(index:index+words_per_line-1)
index = index + words_per_line
read(UNIT=unit,FMT=*) array(loc+1:loc+words_per_line)
loc = loc + words_per_line
end if
end do

9
src/global.f90
View file

@ -107,6 +107,11 @@ module global
& CELL_FILL = 2, & ! Cell filled by a separate universe
& CELL_LATTICE = 3 ! Cell filled with a lattice
! Lattice types
integer, parameter :: &
& LATTICE_RECT = 1, &
& LATTICE_HEX = 2
! array index of universe 0
integer :: BASE_UNIVERSE
@ -247,13 +252,13 @@ contains
integer, intent(in) :: num
character(10) :: str
write ( str, '(I10)' ) num
write (str, '(I10)') num
str = adjustl(str)
end function int_to_str
!=====================================================================
! INT_TO_STR converts a string to an integer.
! STR_TO_INT converts a string to an integer.
!=====================================================================
function str_to_int(str) result(num)

516
src/physics.f90
View file

File diff suppressed because it is too large Load diff

9
src/types.f90
View file

@ -159,7 +159,7 @@ module types
integer :: n_interp ! # of interpolation regions
integer, allocatable :: nbt(:) ! ENDF interpolation parameters
integer, allocatable :: int(:) ! ''
integer :: n_energy ! # of incoming energies
integer :: n_energy ! # of energies for law validity
real(8), allocatable :: energy(:) ! energy grid for law validity
real(8), allocatable :: pvalid(:) ! probability of law validity
real(8), allocatable :: data(:) ! energy distribution data
@ -211,11 +211,10 @@ module types
! Delayed fission neutron emission
integer :: nu_d_type
integer :: n_precursor
real(8), allocatable :: nu_d_data(:)
type(AceDistEnergy) :: nu_d_edist
real(8), allocatable :: nu_d_precursor_const(:,:)
real(8), allocatable :: nu_d_precursor_energy(:,:)
real(8), allocatable :: nu_d_precursor_prob(:,:)
real(8), allocatable :: nu_d_precursor_data(:)
type(AceDistEnergy), allocatable :: nu_d_edist(:)
! Reactions
integer :: n_reaction