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

Fission and level scattering added. Ability to sample secondary energy distributions added (incomplete). Limited MPI capability.

Browse source
This commit is contained in:
Paul Romano 2011-03-18 15:45:07 +00:00
parent 1c42639ef5
commit 6222890a8e
12 changed files with 1411 additions and 252 deletions
Download patch file Download diff file Expand all files Collapse all files

16
src/Makefile
View file

@ -11,6 +11,7 @@ modules = ace.f90 \
geometry.f90 \
global.f90 \
mcnp_random.f90 \
mpi_routines.f90 \
output.f90 \
physics.f90 \
search.f90 \
@ -23,20 +24,24 @@ test_objects = $(modules:.f90=.o) $(test:.f90=.o)
#--------------------------------------------------------------------
# Compiler Options
F90 = gfortran
F90FLAGS = -g
F90 = ifort
F90FLAGS = -g -fpp
MPI = /opt/mpich2/1.3.3-gcc
MPIF90 = $(MPI)/bin/mpif90
MPI_COMPILE_FLAGS = $(shell mpif90 -f90='' -compile_info)
MPI_LINK_FLAGS = $(shell mpif90 -f90='' -link_info)
#--------------------------------------------------------------------
# Targets
all: $(program)
$(program): $(main_objects)
$(F90) $(F90FLAGS) $(main_objects) -o $@
$(MPIF90) -f90=$(F90) $(MPI_LINK_FLAGS) $(F90FLAGS) $(main_objects) -o $@
clean:
@rm -f *.o *.mod $(program)
neat:
@rm -f *.o *.mod
unittest: $(test_objects)
$(F90) $(F90FLAGS) $(test_objects) -o $@
$(F90) $(F90FLAGS) $(MPI_COMPILE_FLAGS) $(test_objects) -o $@
#--------------------------------------------------------------------
# Rules & misc
@ -59,7 +64,8 @@ global.o: types.o
main.o: global.o fileio.o output.o geometry.o mcnp_random.o \
source.o physics.o cross_section.o data_structures.o \
ace.o energy_grid.o
output.o: global.o types.o data_structures.o
mpi_routines.o: global.o output.o types.o
output.o: global.o types.o data_structures.o endf.o
physics.o: types.o global.o mcnp_random.o geometry.o output.o \
search.o endf.o
search.o: output.o

302
src/ace.f90
View file

@ -38,7 +38,7 @@ contains
integer :: n
integer :: index_continuous
integer :: index_thermal
type(DictionaryCI), pointer :: temp_dict
type(DictionaryCI), pointer :: temp_dict => null()
call dict_create(ace_dict)
@ -206,8 +206,9 @@ contains
end do
call parse_ESZ(table)
call parse_MTR(table)
call read_esz(table)
call read_nu_data(table)
call read_reactions(table)
call read_angular_dist(table)
call read_energy_dist(table)
@ -220,14 +221,14 @@ contains
end subroutine read_ACE_continuous
!=====================================================================
! PARSE_ESZ - reads through the ESZ block. This block contains the
! READ_ESZ - reads through the ESZ block. This block contains the
! energy grid, total xs, absorption xs, elastic scattering xs, and
! heating numbers.
!=====================================================================
subroutine parse_ESZ(table)
subroutine read_esz(table)
type(AceContinuous), pointer, intent(inout) :: table
type(AceContinuous), pointer :: table
integer :: NE
@ -254,19 +255,170 @@ contains
table%sigma_el = get_real(NE)
table%heating = get_real(NE)
end subroutine parse_ESZ
end subroutine read_esz
!=====================================================================
! PARSE_MTR - Get the list of reaction MTs for this cross-section
! table. The MT values are somewhat arbitrary. Also read in Q-values,
! neutron multiplicities, and cross-sections.
! =====================================================================
! READ_NU_DATA reads data given on the number of neutrons emitted from
! fission as a function of the incoming energy of a neutron. This data
! may be broken down into prompt and delayed neutrons emitted as well.
!=====================================================================
subroutine parse_MTR(table)
subroutine read_nu_data(table)
type(AceContinuous), pointer :: table
type(AceReaction), pointer :: rxn
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 :: JXS24
integer :: DEC1
integer :: DEC2
JXS2 = JXS(2)
JXS24 = JXS(24)
if (JXS2 == 0) then
! =============================================================
! No prompt/total nu data is present
table % nu_t_type = NU_NONE
table % nu_p_type = NU_NONE
elseif (XSS(JXS2) > 0) then
! =============================================================
! Prompt or total nu data is present
KNU = JXS2
LNU = int(XSS(KNU))
if (LNU == 1) then
! Polynomial data
table % nu_t_type = NU_POLYNOMIAL
table % nu_p_type = NU_NONE
! allocate determine how many coefficients for polynomial
NC = int(XSS(KNU+1))
length = NC + 1
elseif (LNU == 2) then
! Tabular data
table % nu_t_type = NU_TABULAR
table % nu_p_type = NU_NONE
! determine number of interpolation regions and number of
! energies
NR = int(XSS(KNU+1))
NE = int(XSS(KNU+2+2*NR))
length = 2 + 2*NR + 2*NE
end if
! allocate space for nu data storage
allocate(table % nu_t_data(length))
! read data -- for polynomial, this is the number of
! coefficients and the coefficients themselves, and for
! tabular, this is interpolation data and tabular E/nu
XSS_index = KNU + 1
table % nu_t_data = get_real(length)
elseif (XSS(JXS2) < 0) then
! =============================================================
! Prompt and total nu data is present -- read prompt data first
KNU = JXS2 + 1
LNU = XSS(KNU)
if (LNU == 1) then
! Polynomial data
table % nu_p_type = NU_POLYNOMIAL
! allocate determine how many coefficients for polynomial
NC = XSS(KNU+1)
length = NC + 1
elseif (LNU == 2) then
! Tabular data
table % nu_p_type = NU_TABULAR
! determine number of interpolation regions and number of
! energies
NR = XSS(KNU+1)
NE = XSS(KNU+2+2*NR)
length = 2 + 2*NR + 2*NE
end if
! allocate space for nu data storage
allocate(table % nu_p_data(length))
! read data
XSS_index = KNU + 1
table % nu_p_data = get_real(length)
! Now read total nu data
KNU = JXS2 + abs(XSS(JXS2)) + 1
LNU = XSS(KNU)
if (LNU == 1) then
! Polynomial data
table % nu_t_type = NU_POLYNOMIAL
! allocate determine how many coefficients for polynomial
NC = int(XSS(KNU+1))
length = NC + 1
elseif (LNU == 2) then
! Tabular data
table % nu_t_type = NU_TABULAR
! determine number of interpolation regions and number of
! energies
NR = int(XSS(KNU+1))
NE = int(XSS(KNU+2+2*NR))
length = 2 + 2*NR + 2*NE
end if
! allocate space for nu data storage
allocate(table % nu_t_data(length))
! read data
XSS_index = KNU + 1
table % nu_t_data = get_real(length)
end if
if (JXS24 > 0) then
! =============================================================
! Delayed nu data is present
table % nu_d_type = NU_TABULAR
KNU = JXS24
! determine size of tabular delayed nu data
NR = int(XSS(KNU+1))
NE = int(XSS(KNU+2+2*NR))
length = 2 + 2*NR + 2*NE
! allocate space for delayed nu data
allocate(table % nu_d_data(length))
! read delayed nu data
XSS_index = KNU + 1
table % nu_d_data = get_real(length)
! TODO: Read secondary energy distribution
! TODO: Read precursor data
else
table % nu_d_type = NU_NONE
end if
end subroutine read_nu_data
!=====================================================================
! READ_REACTIONS - Get the list of reaction MTs for this cross-section
! table. The MT values are somewhat arbitrary. Also read in Q-values,
! neutron multiplicities, and cross-sections.
!=====================================================================
subroutine read_reactions(table)
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
@ -291,7 +443,7 @@ contains
! Store elastic scattering cross-section on reaction one
rxn => table%reactions(1)
rxn%MT = 2
rxn%Q_value = 0.0_8
rxn%Q_value = ZERO
rxn%TY = 1
rxn%IE = 1
allocate(rxn%sigma(table%n_grid))
@ -318,7 +470,7 @@ contains
rxn % has_energy_dist = .false.
end do
end subroutine parse_MTR
end subroutine read_reactions
!=====================================================================
! READ_ANGULAR_DIST parses the angular distribution for each reaction
@ -328,7 +480,7 @@ contains
type(AceContinuous), pointer :: table
type(AceReaction), pointer :: rxn
type(AceReaction), pointer :: rxn => null()
integer :: JXS8
integer :: JXS9
integer :: NMT
@ -337,7 +489,7 @@ contains
integer :: NP
integer :: LC
integer :: i, j
integer :: size
integer :: length
JXS8 = JXS(8)
JXS9 = JXS(9)
@ -363,49 +515,49 @@ contains
! allocate space for incoming energies and locations
NE = XSS(JXS9 + LOCB - 1)
rxn % adist_n_energy = NE
allocate(rxn % adist_energy(NE))
allocate(rxn % adist_type(NE))
allocate(rxn % adist_location(NE))
rxn % adist % n_energy = NE
allocate(rxn % adist % energy(NE))
allocate(rxn % adist % type(NE))
allocate(rxn % adist % location(NE))
! read incoming energy grid and location of tables
XSS_index = JXS9 + LOCB
rxn % adist_energy = get_real(NE)
rxn % adist_location = get_int(NE)
rxn % adist % energy = get_real(NE)
rxn % adist % location = get_int(NE)
! determine dize of data block
size = 0
length = 0
do j = 1, NE
LC = rxn % adist_location(j)
LC = rxn % adist % location(j)
if (LC == 0) then
! isotropic
rxn % adist_type(j) = ANGLE_ISOTROPIC
rxn % adist % type(j) = ANGLE_ISOTROPIC
elseif (LC > 0) then
! 32 equiprobable bins
rxn % adist_type(j) = ANGLE_32_EQUI
size = size + 33
rxn % adist % type(j) = ANGLE_32_EQUI
length = length + 33
elseif (LC < 0) then
! tabular distribution
rxn % adist_type(j) = ANGLE_TABULAR
rxn % adist % type(j) = ANGLE_TABULAR
NP = XSS(JXS9 + abs(LC))
size = size + 2 + 3*NP
length = length + 2 + 3*NP
end if
end do
! allocate angular distribution data and read
allocate(rxn % adist_data(size))
allocate(rxn % adist % data(length))
! read angular distribution -- currently this does not actually
! parse the angular distribution tables for each incoming
! energy, that must be done on-the-fly
LC = rxn % adist_location(1)
LC = rxn % adist % location(1)
XSS_index = JXS9 + abs(LC) - 1
rxn % adist_data = get_real(size)
rxn % adist % data = get_real(length)
! change location pointers since they are currently relative to
! JXS(9)
LC = abs(rxn % adist_location(1))
rxn % adist_location = abs(rxn % adist_location) - LC + 1
LC = abs(rxn % adist % location(1))
rxn % adist % location = abs(rxn % adist % location) - LC
end do
@ -419,7 +571,7 @@ contains
type(AceContinuous), pointer :: table
type(AceReaction), pointer :: rxn
type(AceReaction), pointer :: rxn => null()
integer :: LOCC
integer :: LED ! location of LDLW block
integer :: LDIS ! location of DLW block
@ -431,7 +583,7 @@ contains
integer :: NP2
integer :: NRa, NEa, NRb, NEb
integer :: IDAT
integer :: start, size, size_interp_data
integer :: start, length, length_interp_data
integer :: i, j, k
LED = JXS(10)
@ -440,6 +592,7 @@ contains
! Loop over all reactions
do i = 1, NXS(5)
rxn => table % reactions(i+1) ! skip over elastic scattering
rxn % has_energy_dist = .true.
! find location of energy distribution data
LOCC = XSS(LED + i - 1)
@ -448,80 +601,87 @@ contains
LAW = XSS(LDIS + LOCC)
IDAT = XSS(LDIS + LOCC + 1)
NR = XSS(LDIS + LOCC + 2)
rxn % edist % law = LAW
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(rxn % edist_nbt(NR))
allocate(rxn % edist_int(NR))
allocate(rxn % edist % nbt(NR))
allocate(rxn % edist % int(NR))
end if
! read ENDF interpolation parameters
XSS_index = LDIS + LOCC + 3
rxn % edist_nbt = get_real(NR)
rxn % edist_int = get_real(NR)
rxn % edist % nbt = get_real(NR)
rxn % edist % int = get_real(NR)
! allocate space for law validity data
NE = XSS(LDIS + LOCC + 3 + 2*NR)
allocate(rxn % edist_energy(NE))
allocate(rxn % edist_pvalid(NE))
allocate(rxn % edist % energy(NE))
allocate(rxn % edist % pvalid(NE))
size_interp_data = 5 + 2*(NR + NE)
length_interp_data = 5 + 2*(NR + NE)
! read law validity data
XSS_index = LDIS + LOCC + 4 + 2*NR
rxn % edist_energy = get_real(NE)
rxn % edist_pvalid = get_real(NE)
rxn % edist % energy = get_real(NE)
rxn % edist % pvalid = get_real(NE)
! Set index to beginning of IDAT array
start = LDIS + IDAT - 2
! Determine size of LDAT array based on which secondary energy
! law it is
size = 0
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)
size = 3 + 2*NR + NE + 3*NET*NE
length = 3 + 2*NR + NE + 3*NET*NE
case (2)
! Discrete photon energy
size = 2
length = 2
case (3)
! Level scattering
size = 2
length = 2
case (4)
! Continuous tabular distribution
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
size = size + 2 + 2*NR + 2*NE
length = length + 2 + 2*NR + 2*NE
do j = 1,NE
NP = XSS(start + size + 2)
size = size + 2 + 3*NP
! 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)
size = 3 + 2*NR + 2*NE + NET
length = 3 + 2*NR + 2*NE + NET
case (7)
! Maxwell fission spectrum
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
size = 3 + 2*NR + 2*NE
length = 3 + 2*NR + 2*NE
case (9)
! Evaporation spectrum
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
size = 3 + 2*NR + 2*NE
length = 3 + 2*NR + 2*NE
case (11)
! Watt spectrum
@ -529,53 +689,57 @@ contains
NEa = XSS(start + 2 + 2*NRa)
NRb = XSS(start + 3 + 2*(NRa+NEa))
NEb = XSS(start + 4 + 2*(NRa+NEa+NRb))
size = 5 + 2*(NRa + NEa + NRb + NEb)
length = 5 + 2*(NRa + NEa + NRb + NEb)
case (44)
! Kalbach-Mann correlated scattering
NR = XSS(start + 1)
NE = XSS(start + 2 + 2*NR)
size = size + 2 + 2*NR + 2*NE
length = length + 2 + 2*NR + 2*NE
do j = 1,NE
NP = XSS(start + size + 2)
size = size + 2 + 5*NP
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)
size = size + 2 + 2*NR + 2*NE
length = length + 2 + 2*NR + 2*NE
do j = 1,NE
! outgoing energy distribution
NP = XSS(start + size + 2)
size = size + 2 + 4*NP
NP = XSS(start + length + 2)
length = length + 2 + 4*NP
do k = 1, NP
! outgoing angle distribution
NP2 = XSS(start + size + 2)
size = size + 2 + 3*NP2
NP2 = XSS(start + length + 2)
length = length + 2 + 3*NP2
end do
end do
case (66)
! N-body phase space distribution
size = 2
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)
size = 4 + 2*(NR + NE + NMU)
length = 4 + 2*(NR + NE + NMU)
end select
! allocate secondary energy distribution array
allocate(rxn % edist_data(size))
allocate(rxn % edist % data(length))
! read secondary energy distribution
XSS_index = start + 1
rxn % edist_data = get_real(size)
rxn % edist % data = get_real(length)
end do

7
src/fileio.f90
View file

@ -271,6 +271,13 @@ contains
case ('xs_data')
path_xsdata = words(2)
case ('verbosity')
verbosity = str_to_int(words(2))
if (verbosity == ERROR_CODE) then
msg = "Invalid verbosity: " // words(2)
call error(msg)
end if
case default
! do nothing
end select

7
src/geometry.f90
View file

@ -107,8 +107,6 @@ contains
c => cells(univ % cells(i))
if (cell_contains(c, p)) then
p%cell = dict_get_key(cell_dict, c % uid)
! If this cell contains a universe of lattice, search for
! the particle in that universe/lattice
if (c % fill > 0) then
@ -125,6 +123,11 @@ contains
found = .true.
cCell => c
cMaterial => materials(cCell%material)
cUniverse => univ
! set particle attributes
p%cell = dict_get_key(cell_dict, c % uid)
p%universe = dict_get_key(universe_dict, univ % uid)
exit
end if
end if

36
src/global.f90
View file

@ -56,6 +56,18 @@ module global
! Source and fission bank
type(Particle), allocatable, target :: source_bank(:)
type(Bank), allocatable, target :: fission_bank(:)
integer :: n_bank ! # of sites in fission bank
integer :: bank_first ! index of first particle in bank
integer :: bank_last ! index of last particle in bank
! cycle keff
real(8) :: keff
! Parallel processing variables
integer :: n_procs ! number of processes
integer :: rank ! rank of process
logical :: master ! master process?
logical :: mpi_enabled ! is MPI in use and initialized?
! Paths to input file, cross section data, etc
character(100) :: &
@ -72,7 +84,8 @@ module global
& K_BOLTZMANN = 8.617342e-5, & ! Boltzmann constant in eV/K
& INFINITY = huge(0.0_8), & ! positive infinity
& ZERO = 0.0_8, &
& ONE = 1.0_8
& ONE = 1.0_8, &
& TWO = 2.0_8
! Boundary conditions
integer, parameter :: &
@ -146,19 +159,28 @@ module global
& PHOTON = 2, &
& ELECTRON = 3
! Fission neutron emission (nu) type
integer, parameter :: &
& NU_NONE = 0, & ! No nu values (non-fissionable)
& NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
& NU_TABULAR = 2 ! Nu values given by tabular distribution
! Integer code for read error -- better hope this number is never
! used in an input file!
integer, parameter :: &
& ERROR_CODE = -huge(0)
integer :: verbosity = 5
! The verbosity controls how much information will be printed to the
! screen and in logs
integer :: verbosity
integer, parameter :: max_words = 500
integer, parameter :: max_line = 250
! Versioning numbers
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 2
integer, parameter :: VERSION_RELEASE = 0
integer, parameter :: VERSION_RELEASE = 1
contains
@ -174,6 +196,9 @@ contains
! Default number of particles
n_particles = 10000
! Default verbosity
verbosity = 5
end subroutine set_defaults
!=====================================================================
@ -183,6 +208,8 @@ contains
subroutine free_memory()
integer :: ierr
! Deallocate cells, surfaces, materials
if (allocated(cells)) deallocate(cells)
if (allocated(surfaces)) deallocate(surfaces)
@ -195,6 +222,9 @@ contains
if (allocated(fission_bank)) deallocate(fission_bank)
if (allocated(source_bank)) deallocate(source_bank)
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(ierr)
! End program
stop

3
src/input_sample
View file

@ -20,4 +20,5 @@ source box -4 -4 -4 4 4 4
xs_library endfb7
xs_data /opt/serpent/xsdata/endfb7
criticality 1 1 15
criticality 1 1 100
verbosity 8

23
src/main.f90
View file

@ -12,6 +12,7 @@ program main
use cross_section, only: read_xsdata, material_total_xs
use ace, only: read_xs
use energy_grid, only: unionized_grid, original_indices
use mpi_routines
implicit none
@ -19,12 +20,16 @@ program main
character(250) :: msg
type(Universe), pointer :: univ
! Print the OpenMC title and version/date/time information
call title()
verbosity = 10
! Setup MPI
call setup_mpi()
! Initialize random number generator
call RN_init_problem( 3, 0_8, 0_8, 0_8, 0 )
! Print the OpenMC title and version/date/time information
if (master) call title()
! Initialize random number generator. The first argument corresponds
! to which random number generator to use- in this case one of the
! L'Ecuyer 63-bit RNGs.
call RN_init_problem(3, 0_8, 0_8, 0_8, 0)
! Set default values for settings
call set_defaults()
@ -65,8 +70,10 @@ program main
! calculate total material cross-sections for sampling path lenghts
call material_total_xs()
call echo_input()
call print_summary()
if (master) then
! call echo_input()
! call print_summary()
end if
! create source particles
call init_source()
@ -95,6 +102,7 @@ contains
CYCLE_LOOP: do i_cycle = 1, n_cycles
! Set all tallies to zero
n_bank = 0
HISTORY_LOOP: do j = 1, n_particles
@ -115,6 +123,7 @@ contains
end do HISTORY_LOOP
! Collect results and statistics
print *, n_bank
! print cycle information

56
src/mpi_routines.f90 Normal file
View file

@ -0,0 +1,56 @@
module mpi_routines
use mpi
use global
use output, only: message, error
implicit none
integer :: MPI_bank ! MPI datatype for fission bank
integer(8) :: bank_index ! Fission bank site unique identifier
contains
!=====================================================================
! SETUP_MPI initilizes the Message Passing Interface (MPI) and
! determines the number of processors the problem is being run with as
! well as the rank of each processor.
!=====================================================================
subroutine setup_mpi()
integer :: ierr
character(250) :: msg
! Initialize MPI
call MPI_INIT(ierr)
if (ierr /= MPI_SUCCESS) then
msg = "Failed to initialize MPI."
call error(msg)
end if
mpi_enabled = .true.
! Determine number of processors
call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, ierr)
if (ierr /= MPI_SUCCESS) then
msg = "Could not determine number of processors."
call error(msg)
end if
! Determine rank of each processor
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
if (ierr /= MPI_SUCCESS) then
msg = "Could not determine MPI rank."
call error(msg)
end if
! Determine master
if (rank == 0) then
master = .true.
else
master = .false.
end if
end subroutine setup_mpi
end module mpi_routines

144
src/output.f90
View file

@ -4,9 +4,13 @@ module output
use global
use types, only: Cell, Universe, Surface
use data_structures, only: dict_get_key
use endf, only: reaction_name
implicit none
integer :: ou = OUTPUT_UNIT
integer :: eu = ERROR_UNIT
contains
!=====================================================================
@ -19,9 +23,6 @@ module output
character(10) :: date
character(8) :: time
integer :: ou
ou = OUTPUT_UNIT
write(ou,*)
write(ou,*) ' .d88888b. 888b d888 .d8888b. '
@ -61,9 +62,6 @@ module output
subroutine echo_input()
character(32) :: string
integer :: ou
ou = OUTPUT_UNIT
! Display problem summary
write(ou,*) '============================================='
@ -104,18 +102,17 @@ module output
character(*), intent(in) :: msg
integer, intent(in) :: level
integer :: ou
integer :: n_lines
integer :: i
! Only allow master to print to screen
if (.not. master) return
if ( level <= verbosity ) then
ou = OUTPUT_UNIT
n_lines = (len_trim(msg)-1)/79 + 1
do i = 1, n_lines
write(ou, fmt='(1X,A79)') msg(79*(i-1)+1:79*i)
end do
end if
end subroutine message
@ -130,9 +127,10 @@ module output
character(*), intent(in) :: msg
integer :: n_lines
integer :: i, ou
integer :: i
ou = OUTPUT_UNIT
! Only allow master to print to screen
if (.not. master) return
write(ou, fmt='(1X,A9)', advance='no') 'WARNING: '
@ -158,9 +156,10 @@ module output
character(*), intent(in) :: msg
integer :: n_lines
integer :: i, eu
integer :: i
eu = ERROR_UNIT
! Only allow master to print to screen
if (.not. master) return
write(eu, fmt='(1X,A7)', advance='no') 'ERROR: '
@ -183,7 +182,7 @@ module output
! execution and returns it in a readable format
!=====================================================================
subroutine get_today( today_date, today_time )
subroutine get_today(today_date, today_time)
character(10), intent(out) :: today_date
character(8), intent(out) :: today_time
@ -220,6 +219,91 @@ module output
end subroutine get_today
!=====================================================================
! PRINT_PARTICLE displays the attributes of a particle
!=====================================================================
subroutine print_particle(p)
type(Particle), pointer :: p
integer :: i
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
character(250) :: string
select case (p % type)
case (NEUTRON)
write(ou,*) 'Neutron ' // int_to_str(p % uid)
case (PHOTON)
write(ou,*) 'Photon ' // int_to_str(p % uid)
case (ELECTRON)
write(ou,*) 'Electron ' // int_to_str(p % uid)
case default
write(ou,*) 'Unknown Particle ' // int_to_str(p % uid)
end select
write(ou,100) 'x = ', p % xyz(1)
write(ou,100) 'y = ', p % xyz(2)
write(ou,100) 'z = ', p % xyz(3)
write(ou,100) 'u = ', p % uvw(1)
write(ou,100) 'v = ', p % uvw(2)
write(ou,100) 'w = ', p % uvw(3)
write(ou,100) 'Weight = ', p % wgt
write(ou,100) 'Energy = ', p % E
write(ou,*) ' IE = ' // int_to_str(p % IE)
write(ou,100) 'Interpolation factor = ', p % interp
if (p % cell > 0) then
c => cells(p % cell)
write(ou,*) ' Cell = ' // int_to_str(c % uid)
else
write(ou,*) ' Cell not determined'
end if
if (p % surface > 0) then
s => surfaces(p % surface)
write(ou,*) ' Surface = ' // int_to_str(s % uid)
else
write(ou,*) ' Surface = None'
end if
u => universes(p % universe)
write(ou,*) ' Universe = ' // int_to_str(u % uid)
write(ou,*)
! Format for a single real
100 format (5X,A,G10.3)
nullify(c)
nullify(s)
nullify(u)
end subroutine print_particle
!=====================================================================
! PRINT_REACTION displays the attributes of a reaction
!=====================================================================
subroutine print_reaction(rxn)
type(AceReaction), pointer :: rxn
write(ou,*) 'Reaction ' // reaction_name(rxn % MT)
write(ou,*) ' MT = ' // int_to_str(rxn % MT)
write(ou,100) 'Q-value = ', rxn % Q_value
write(ou,*) ' TY = ' // int_to_str(rxn % TY)
write(ou,*) ' Starting index = ' // int_to_str(rxn % IE)
if (rxn % has_energy_dist) then
write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law)
end if
write(ou,*)
! Format for a single real
100 format (5X,A,G10.3)
end subroutine print_reaction
!=====================================================================
! PRINT_CELL displays the attributes of a cell
!=====================================================================
@ -228,15 +312,12 @@ module output
type(Cell), pointer :: c
integer :: ou
integer :: temp
integer :: i
type(Universe), pointer :: u => null()
type(Material), pointer :: m => null()
character(250) :: string
ou = OUTPUT_UNIT
write(ou,*) 'Cell ' // int_to_str(c % uid)
temp = dict_get_key(cell_dict, c % uid)
write(ou,*) ' Array Index = ' // int_to_str(temp)
@ -287,12 +368,9 @@ module output
type(Universe), pointer :: univ
integer :: ou
integer :: i
character(250) :: string
type(Cell), pointer :: c
ou = OUTPUT_UNIT
type(Cell), pointer :: c => null()
write(ou,*) 'Universe ' // int_to_str(univ % uid)
write(ou,*) ' Level = ' // int_to_str(univ % level)
@ -304,6 +382,8 @@ module output
write(ou,*) ' Cells =' // trim(string)
write(ou,*)
nullify(c)
end subroutine print_universe
!=====================================================================
@ -314,12 +394,9 @@ module output
type(Surface), pointer :: surf
integer :: ou
integer :: i
character(80) :: string
ou = OUTPUT_UNIT
write(ou,*) 'Surface ' // int_to_str(surf % uid)
select case (surf % type)
case (SURF_PX)
@ -354,7 +431,7 @@ module output
string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i))
end do
end if
write(ou,*) ' Positive Neighbors = ', string
write(ou,*) ' Positive Neighbors = ', trim(string)
string = ""
if (allocated(surf % neighbor_neg)) then
@ -362,7 +439,7 @@ module output
string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i))
end do
end if
write(ou,*) ' Negative Neighbors =', string
write(ou,*) ' Negative Neighbors =', trim(string)
end subroutine print_surface
@ -374,14 +451,11 @@ module output
type(Material), pointer :: mat
integer :: ou
integer :: i
integer :: n_lines
type(AceContinuous), pointer :: table
character(250) :: string
ou = OUTPUT_UNIT
write(ou,*) 'Material ' // int_to_str(mat % uid)
! Make string of all isotopes
string = ""
@ -403,6 +477,8 @@ module output
& ' atom/b-cm'
write(ou,*)
nullify(table)
end subroutine print_material
!=====================================================================
@ -418,9 +494,6 @@ module output
type(Universe), pointer :: u
type(Material), pointer :: m
integer :: i
integer :: ou
ou = OUTPUT_UNIT
! print summary of cells
write(ou,*) '============================================='
@ -462,6 +535,11 @@ module output
call print_material(m)
end do
nullify(s)
nullify(c)
nullify(u)
nullify(m)
end subroutine print_summary
end module output

916
src/physics.f90
View file

File diff suppressed because it is too large Load diff

51
src/source.f90
View file

@ -20,6 +20,7 @@ contains
real(8) :: phi ! azimuthal angle
real(8) :: mu ! cosine of polar angle
real(8) :: p_min(3), p_max(3)
type(Particle), pointer :: p => null()
! Allocate fission and source banks
allocate( source_bank(n_particles) )
@ -34,24 +35,28 @@ contains
! Initialize first cycle source bank
do i = 1, n_particles
call RN_init_particle(int(i,8))
p => source_bank(i)
! position
r = (/ (rang(), j = 1,3) /)
source_bank(i)%uid = i
source_bank(i)%xyz = p_min + r*(p_max - p_min)
p % uid = i
p % xyz = p_min + r*(p_max - p_min)
! angle
phi = 2.0_8*PI*rang()
mu = 2.0_8*rang() - 1.0_8
source_bank(i)%uvw(1) = mu
source_bank(i)%uvw(2) = sqrt(1. - mu**2) * cos(phi)
source_bank(i)%uvw(3) = sqrt(1. - mu**2) * sin(phi)
phi = TWO*PI*rang()
mu = TWO*rang() - ONE
p % uvw(1) = mu
p % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
p % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! cell
source_bank(i)%cell = 0
! set defaults
p % type = NEUTRON
p % cell = 0
p % surface = 0
p % universe = 0
p % wgt = ONE
p % alive = .true.
! set particle to be alive
source_bank(i)%alive = .true.
end do
! Reset source index
@ -81,4 +86,28 @@ contains
end function get_source_particle
!=====================================================================
! ADD_BANK_SITES
!=====================================================================
subroutine add_bank_sites(p, table, n)
type(Particle), pointer :: p
type(AceContinuous), pointer :: table
integer, intent(in) :: n
integer :: i
if (n == 0 .or. n_bank == 3*n_particles) return
do i = n_bank + 1, min(n_bank + n, 3*n_particles)
! Copy particle data
fission_bank(i)%uid = p%uid
fission_bank(i)%xyz = p%xyz
! TODO: Sample angle and energy from secondary distribution
end do
n_bank = min(n_bank + n, 3*n_particles)
end subroutine add_bank_sites
end module source

102
src/types.f90
View file

@ -71,14 +71,15 @@ module types
integer :: type ! Particle type (n, p, e, etc)
real(8) :: xyz(3) ! location
real(8) :: uvw(3) ! directional cosines
real(8) :: wgt ! particle weight
real(8) :: E ! energy
integer :: IE ! index on energy grid
real(8) :: interp ! interpolation factor for energy grid
integer :: cell ! current cell
integer :: universe ! current universe
integer :: surface ! current surface
real(8) :: wgt ! particle weight
logical :: alive ! is particle alive?
integer :: n_coll ! # of collisions
end type Particle
!=====================================================================
@ -89,7 +90,9 @@ module types
type Bank
integer :: uid ! Unique ID
real(8) :: xyz(3) ! Location of bank particle
real(8) :: xyz(3) ! location of bank particle
real(8) :: uvw(3) ! diretional cosines
real(8) :: E ! energy
end type Bank
!=====================================================================
@ -133,6 +136,35 @@ module types
real(8), allocatable :: values(:) ! values for particular source type
end type ExtSource
!=====================================================================
! ACEDISTANGLE contains data for a tabular secondary angle
! distribution whether it be tabular or 32 equiprobable cosine bins
!=====================================================================
type AceDistAngle
integer :: n_energy ! # of incoming energies
real(8), allocatable :: energy(:) ! incoming energy grid
integer, allocatable :: type(:) ! type of distribution
integer, allocatable :: location(:) ! location of each table
real(8), allocatable :: data(:) ! angular distribution data
end type AceDistAngle
!=====================================================================
! ACEDISTENERGY contains data for a secondary energy distribution for
! all scattering laws
!=====================================================================
type AceDistEnergy
integer :: law ! secondary distribution law
integer :: n_interp ! # of interpolation regions
integer, allocatable :: nbt(:) ! ENDF interpolation parameters
integer, allocatable :: int(:) ! ''
integer :: n_energy ! # of incoming energies
real(8), allocatable :: energy(:) ! energy grid for law validity
real(8), allocatable :: pvalid(:) ! probability of law validity
real(8), allocatable :: data(:) ! energy distribution data
end type AceDistEnergy
!=====================================================================
! ACEREACTION contains the cross-section and secondary energy and
! angle distributions for a single reaction in a continuous-energy
@ -142,28 +174,13 @@ module types
type AceReaction
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
real(8) :: TY ! Number of neutrons released
integer :: TY ! Number of neutrons released
integer :: IE ! Starting energy grid index
real(8), allocatable :: sigma(:) ! Cross section values
logical :: has_angle_dist
logical :: has_energy_dist
! Secondary angle distribution
integer :: adist_n_energy ! # of incoming energies
real(8), allocatable :: adist_energy(:) ! incoming energy grid
integer, allocatable :: adist_type(:) ! type of distribution
integer, allocatable :: adist_location(:) ! location of each table
real(8), allocatable :: adist_data(:) ! angular distribution data
! Secondary energy distribution
integer :: edist_law ! secondary distribution law
integer :: edist_n_interp ! # of interpolation regions
integer, allocatable :: edist_nbt(:) ! ENDF interpolation parameters
integer, allocatable :: edist_int(:) ! ''
integer :: edist_n_energy ! # of incoming energies
real(8), allocatable :: edist_energy(:) ! energy grid for law validity
real(8), allocatable :: edist_pvalid(:) ! probability of law validity
real(8), allocatable :: edist_data(:) ! energy distribution data
logical :: has_angle_dist ! Angle distribution present?
logical :: has_energy_dist ! Energy distribution present?
type(AceDistAngle) :: adist ! Secondary angular distribution
type(AceDistEnergy) :: edist ! Secondary energy distribution
end type AceReaction
!=====================================================================
@ -184,19 +201,23 @@ module types
real(8), allocatable :: sigma_el(:)
real(8), allocatable :: heating(:)
integer :: nu_p_type
real(8), allocatable :: nu_p_energy(:)
real(8), allocatable :: nu_p_value(:)
! Total fission neutron emission
integer :: nu_t_type
real(8), allocatable :: nu_t_energy(:)
real(8), allocatable :: nu_t_value(:)
real(8), allocatable :: nu_t_data(:)
real(8), allocatable :: nu_d_energy(:)
real(8), allocatable :: nu_d_value(:)
! Prompt fission neutron emission
integer :: nu_p_type
real(8), allocatable :: nu_p_data(:)
! Delayed fission neutron emission
integer :: nu_d_type
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(:,:)
! Reactions
integer :: n_reaction
type(AceReaction), pointer :: reactions(:) => null()
@ -220,27 +241,6 @@ module types
real(8), allocatable :: elastic_mu_out(:)
end type AceThermal
!=====================================================================
! ACEDISTANGLE contains data for a tabular secondary angle
! distribution whether it be tabular or 32 equiprobable cosine bins
!=====================================================================
!!$ type AceDistAngle
!!$ integer :: n_energy
!!$ real(8), allocatable :: energy
!!$ integer, allocatable :: location
!!$ real(8), allocatable :: data
!!$ end type AceDistAngle
!=====================================================================
! ACEDISTTAB contains data for a tabular secondary energy distribution
! according to MCNP Law 4 (ENDF Law 1)
!=====================================================================
type AceDistEnergy
integer :: law
end type AceDistEnergy
!=====================================================================
! XSDATA contains data read in from a SERPENT xsdata file
!=====================================================================