mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 14:35:27 -04:00
Sampling of reactions, several important bugs fixed.
This commit is contained in:
Paul Romano
parent
2610260f56
commit
e6af8d3544
10 changed files with 422 additions and 103 deletions
19
ChangeLog
19
ChangeLog
|
|
@ -1,3 +1,22 @@
|
|||
2011-02-13 Paul Romano <romano7@mit.edu>
|
||||
|
||||
* ace.f90: Fixed xs_continuous allocation again -- cross sections
|
||||
were not being assigned correctly to materials. Added elastic
|
||||
scattering in list of reactions. Changed name of rxn %
|
||||
energy_index to IE.
|
||||
* endf.f90: New module. Includes subroutine reaction_name which
|
||||
returns the name of a reaction for a given ENDF MT value.
|
||||
* geometry.f90: Fixed bug in dist_to_boundary (the 'w' in uvw was
|
||||
actually being set to 'z', not 'w'). Changed 0's and 1's to ZEROs
|
||||
and ONEs.
|
||||
* global.f90: Added ZERO and ONE.
|
||||
* physics.f90: Can now sample reactions once nuclide is
|
||||
selected. Added subroutines elastic_scatter, level_inelastic, and
|
||||
n_gamma (level_inelastic is just a stub), and significantly
|
||||
changed elastic scattering (currently based on stationary target).
|
||||
* types.f90: Added n_reaction attribute on AceContinuous, changed
|
||||
AceReaction % energy_index to IE.
|
||||
|
||||
2011-02-09 Paul Romano <romano7@mit.edu>
|
||||
|
||||
* ace.f90: Fixed bug in the way xs_continuous was allocated. If,
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@ test = unittest.f90
|
|||
modules = ace.f90 \
|
||||
cross_section.f90 \
|
||||
data_structures.f90 \
|
||||
endf.f90 \
|
||||
energy_grid.f90 \
|
||||
fileio.f90 \
|
||||
geometry.f90 \
|
||||
|
|
@ -50,6 +51,7 @@ unittest: $(test_objects)
|
|||
ace.o: global.o output.o string.o fileio.o string.o
|
||||
cross_section.o: global.o string.o data_structures.o output.o
|
||||
data_structures.o: global.o
|
||||
endf.o: global.o
|
||||
energy_grid.o: global.o output.o data_structures.o
|
||||
fileio.o: types.o global.o string.o output.o data_structures.o
|
||||
geometry.o: types.o global.o output.o string.o data_structures.o
|
||||
|
|
@ -59,7 +61,7 @@ main.o: global.o fileio.o output.o geometry.o mcnp_random.o \
|
|||
ace.o energy_grid.o
|
||||
output.o: global.o
|
||||
physics.o: types.o global.o mcnp_random.o geometry.o output.o \
|
||||
search.o
|
||||
search.o endf.o
|
||||
search.o: output.o
|
||||
source.o: global.o mcnp_random.o
|
||||
string.o: global.o output.o
|
||||
|
|
|
|||
58
src/ace.f90
58
src/ace.f90
|
|
@ -6,7 +6,7 @@ module ace
|
|||
use fileio, only: read_line, read_data, skip_lines
|
||||
use string, only: split_string, str_to_real
|
||||
use data_structures, only: dict_create, dict_add_key, dict_has_key, &
|
||||
& dict_get_key
|
||||
& dict_get_key, dict_delete
|
||||
|
||||
integer :: NXS(16)
|
||||
integer :: JXS(32)
|
||||
|
|
@ -37,7 +37,7 @@ contains
|
|||
integer :: n
|
||||
integer :: index_continuous
|
||||
integer :: index_thermal
|
||||
type(ListKeyValueCI), pointer :: elem
|
||||
type(DictionaryCI), pointer :: temp_dict
|
||||
|
||||
call dict_create(ace_dict)
|
||||
|
||||
|
|
@ -78,14 +78,31 @@ contains
|
|||
allocate(xs_thermal(n_thermal))
|
||||
|
||||
! loop over all nuclides in xsdata
|
||||
call dict_create(temp_dict)
|
||||
|
||||
index_continuous = 0
|
||||
do i = 1, size(xsdatas)
|
||||
key = xsdatas(i)%alias
|
||||
if (dict_has_key(ace_dict, key)) then
|
||||
index_continuous = index_continuous + 1
|
||||
call read_ACE_continuous(index_continuous, i)
|
||||
end if
|
||||
index_thermal = 0
|
||||
do i = 1, n_materials
|
||||
mat => materials(i)
|
||||
do j = 1, mat%n_isotopes
|
||||
index = mat%isotopes(j)
|
||||
key = xsdatas(index)%id
|
||||
n = len_trim(key)
|
||||
call lower_case(key)
|
||||
select case (key(n:n))
|
||||
case ('c')
|
||||
if (.not. dict_has_key(temp_dict, key)) then
|
||||
index_continuous = index_continuous + 1
|
||||
call read_ACE_continuous(index_continuous, index)
|
||||
end if
|
||||
case ('t')
|
||||
n_thermal = n_thermal + 1
|
||||
end select
|
||||
end do
|
||||
end do
|
||||
|
||||
! delete dictionary
|
||||
call dict_delete(temp_dict)
|
||||
|
||||
end subroutine read_xs
|
||||
|
||||
|
|
@ -222,7 +239,11 @@ contains
|
|||
allocate(table%sigma_el(NE))
|
||||
allocate(table%heating(NE))
|
||||
|
||||
! read data from XSS
|
||||
! 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
|
||||
! generation.
|
||||
XSS_index = 1
|
||||
table%energy = get_real(NE)
|
||||
table%sigma_t = get_real(NE)
|
||||
|
|
@ -259,11 +280,22 @@ contains
|
|||
JXS7 = JXS(7)
|
||||
NMT = NXS(4)
|
||||
|
||||
! allocate array of reactions
|
||||
allocate(table%reactions(NMT))
|
||||
! allocate array of reactions. Add one since we need to include an
|
||||
! elastic scattering channel
|
||||
table%n_reaction = NMT + 1
|
||||
allocate(table%reactions(NMT+1))
|
||||
|
||||
! Store elastic scattering cross-section on reaction one
|
||||
rxn => table%reactions(1)
|
||||
rxn%MT = 2
|
||||
rxn%Q_value = 0.0_8
|
||||
rxn%TY = 1
|
||||
rxn%IE = 1
|
||||
allocate(rxn%sigma(table%n_grid))
|
||||
rxn%sigma = table%sigma_el
|
||||
|
||||
do i = 1, NMT
|
||||
rxn => table%reactions(i)
|
||||
rxn => table%reactions(i+1)
|
||||
|
||||
! read MT number, Q-value, and neutrons produced
|
||||
rxn%MT = XSS(LMT+i-1)
|
||||
|
|
@ -272,7 +304,7 @@ contains
|
|||
|
||||
! read cross section values
|
||||
LOCA = XSS(LXS+i-1)
|
||||
rxn%energy_index = XSS(JXS7 + LOCA - 1)
|
||||
rxn%IE = XSS(JXS7 + LOCA - 1)
|
||||
NE = XSS(JXS7 + LOCA)
|
||||
allocate(rxn%sigma(NE))
|
||||
XSS_index = JXS7 + LOCA + 1
|
||||
|
|
|
|||
147
src/endf.f90
Normal file
147
src/endf.f90
Normal file
|
|
@ -0,0 +1,147 @@
|
|||
module endf
|
||||
|
||||
use global, only: int_to_str
|
||||
|
||||
contains
|
||||
|
||||
!=====================================================================
|
||||
! REACTION_NAME gives the name of the reaction for a given MT value
|
||||
!=====================================================================
|
||||
|
||||
function reaction_name(MT) result(string)
|
||||
|
||||
integer, intent(in) :: MT
|
||||
character(20) :: string
|
||||
|
||||
select case (MT)
|
||||
case (1)
|
||||
string = '(n,total)'
|
||||
case (2)
|
||||
string = '(elastic)'
|
||||
case (11)
|
||||
string = '(n,2nd)'
|
||||
case (16)
|
||||
string = '(n,2n)'
|
||||
case (17)
|
||||
string = '(n,3n)'
|
||||
case (18)
|
||||
string = '(n,fission)'
|
||||
case (19)
|
||||
string = '(n,f)'
|
||||
case (20)
|
||||
string = '(n,2nf)'
|
||||
case (21)
|
||||
string = '(n,2nf)'
|
||||
case (22)
|
||||
string = '(n,na)'
|
||||
case (23)
|
||||
string = '(n,n3a)'
|
||||
case (24)
|
||||
string = '(n,2na)'
|
||||
case (25)
|
||||
string = '(n,3na)'
|
||||
case (28)
|
||||
string = '(n,np)'
|
||||
case (29)
|
||||
string = '(n,n2a)'
|
||||
case (30)
|
||||
string = '(n,2n2a)'
|
||||
case (32)
|
||||
string = '(n,nd)'
|
||||
case (33)
|
||||
string = '(n,nt)'
|
||||
case (34)
|
||||
string = '(n,n He-3)'
|
||||
case (35)
|
||||
string = '(n,nd3a)'
|
||||
case (36)
|
||||
string = '(n,nt2a)'
|
||||
case (37)
|
||||
string = '(n,4n)'
|
||||
case (38)
|
||||
string = '(n,3nf)'
|
||||
case (41)
|
||||
string = '(n,2np)'
|
||||
case (42)
|
||||
string = '(n,3np)'
|
||||
case (44)
|
||||
string = '(n,2np)'
|
||||
case (45)
|
||||
string = '(n,npa)'
|
||||
case (51 : 90)
|
||||
string = '(n,n' // trim(int_to_str(MT-50)) // ')'
|
||||
case (91)
|
||||
string = '(n,nc)'
|
||||
case (102)
|
||||
string = '(n,gamma)'
|
||||
case (103)
|
||||
string = '(n,p)'
|
||||
case (104)
|
||||
string = '(n,d)'
|
||||
case (105)
|
||||
string = '(n,t)'
|
||||
case (106)
|
||||
string = '(n,3He)'
|
||||
case (107)
|
||||
string = '(n,a)'
|
||||
case (108)
|
||||
string = '(n,2a)'
|
||||
case (109)
|
||||
string = '(n,3a)'
|
||||
case (111)
|
||||
string = '(n,2p)'
|
||||
case (112)
|
||||
string = '(n,pa)'
|
||||
case (113)
|
||||
string = '(n,t2a)'
|
||||
case (114)
|
||||
string = '(n,d2a)'
|
||||
case (115)
|
||||
string = '(n,pd)'
|
||||
case (116)
|
||||
string = '(n,pt)'
|
||||
case (117)
|
||||
string = '(n,da)'
|
||||
case (201)
|
||||
string = '(n,Xn)'
|
||||
case (202)
|
||||
string = '(n,Xgamma)'
|
||||
case (203)
|
||||
string = '(n,Xp)'
|
||||
case (204)
|
||||
string = '(n,Xd)'
|
||||
case (205)
|
||||
string = '(n,Xt)'
|
||||
case (206)
|
||||
string = '(n,X3He)'
|
||||
case (207)
|
||||
string = '(n,Xa)'
|
||||
case (444)
|
||||
string = '(damage)'
|
||||
case (600 : 648)
|
||||
string = '(n,p' // trim(int_to_str(MT-600)) // ')'
|
||||
case (649)
|
||||
string = '(n,pc)'
|
||||
case (650 : 698)
|
||||
string = '(n,d' // trim(int_to_str(MT-650)) // ')'
|
||||
case (699)
|
||||
string = '(n,dc)'
|
||||
case (700 : 748)
|
||||
string = '(n,t' // trim(int_to_str(MT-600)) // ')'
|
||||
case (749)
|
||||
string = '(n,tc)'
|
||||
case (750 : 798)
|
||||
string = '(n,3He' // trim(int_to_str(MT-650)) // ')'
|
||||
case (799)
|
||||
string = '(n,3Hec)'
|
||||
case (800 : 848)
|
||||
string = '(n,a' // trim(int_to_str(MT-800)) // ')'
|
||||
case (849)
|
||||
string = '(n,tc)'
|
||||
case default
|
||||
string = 'MT=' // trim(int_to_str(MT))
|
||||
end select
|
||||
|
||||
end function reaction_name
|
||||
|
||||
end module endf
|
||||
|
|
@ -87,7 +87,7 @@ contains
|
|||
|
||||
subroutine find_cell(neut)
|
||||
|
||||
type(Neutron), pointer, intent(inout) :: neut
|
||||
type(Neutron), pointer :: neut
|
||||
|
||||
type(Cell), pointer :: this_cell
|
||||
logical :: found_cell
|
||||
|
|
@ -123,7 +123,7 @@ contains
|
|||
|
||||
subroutine cross_boundary(neut)
|
||||
|
||||
type(Neutron), pointer, intent(in) :: neut
|
||||
type(Neutron), pointer :: neut
|
||||
|
||||
type(Surface), pointer :: surf
|
||||
type(Cell), pointer :: c
|
||||
|
|
@ -155,6 +155,8 @@ contains
|
|||
c => cells(index_cell)
|
||||
if (cell_contains(c, neut)) then
|
||||
neut%cell = index_cell
|
||||
cCell => cells(index_cell)
|
||||
cMaterial => materials(cCell%material)
|
||||
return
|
||||
end if
|
||||
end do
|
||||
|
|
@ -166,6 +168,8 @@ contains
|
|||
c => cells(index_cell)
|
||||
if (cell_contains(c, neut)) then
|
||||
neut%cell = index_cell
|
||||
cCell => cells(index_cell)
|
||||
cMaterial => materials(cCell%material)
|
||||
return
|
||||
end if
|
||||
end do
|
||||
|
|
@ -177,6 +181,8 @@ contains
|
|||
c => cells(i)
|
||||
if (cell_contains(c, neut)) then
|
||||
neut%cell = i
|
||||
cCell => cells(i)
|
||||
cMaterial => materials(cCell%material)
|
||||
return
|
||||
end if
|
||||
end do
|
||||
|
|
@ -196,7 +202,7 @@ contains
|
|||
|
||||
subroutine dist_to_boundary(neut, dist, surf, other_cell)
|
||||
|
||||
type(Neutron), intent(in) :: neut
|
||||
type(Neutron), pointer :: neut
|
||||
real(8), intent(out) :: dist
|
||||
integer, intent(out) :: surf
|
||||
integer, optional, intent(in) :: other_cell
|
||||
|
|
@ -224,18 +230,18 @@ contains
|
|||
|
||||
current_surf = neut%surface
|
||||
|
||||
x = neut%xyz(1)
|
||||
y = neut%xyz(2)
|
||||
z = neut%xyz(3)
|
||||
u = neut%uvw(1)
|
||||
v = neut%uvw(2)
|
||||
z = neut%uvw(3)
|
||||
|
||||
dist = INFINITY
|
||||
n_boundaries = size(cell_p%boundary_list)
|
||||
allocate(expression(n_boundaries))
|
||||
expression = cell_p%boundary_list
|
||||
do i = 1, n_boundaries
|
||||
x = neut%xyz(1)
|
||||
y = neut%xyz(2)
|
||||
z = neut%xyz(3)
|
||||
u = neut%uvw(1)
|
||||
v = neut%uvw(2)
|
||||
w = neut%uvw(3)
|
||||
|
||||
! check for coincident surfaec
|
||||
index_surf = expression(i)
|
||||
if (index_surf == current_surf) cycle
|
||||
|
|
@ -247,30 +253,30 @@ contains
|
|||
surf_p => surfaces(index_surf)
|
||||
select case (surf_p%type)
|
||||
case (SURF_PX)
|
||||
if (u == 0.0) then
|
||||
if (u == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
x0 = surf_p%coeffs(1)
|
||||
d = (x0 - x)/u
|
||||
if (d < 0) d = INFINITY
|
||||
if (d < ZERO) d = INFINITY
|
||||
end if
|
||||
|
||||
case (SURF_PY)
|
||||
if (v == 0.0) then
|
||||
if (v == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
y0 = surf_p%coeffs(1)
|
||||
d = (y0 - y)/v
|
||||
if (d < 0) d = INFINITY
|
||||
if (d < ZERO) d = INFINITY
|
||||
end if
|
||||
|
||||
case (SURF_PZ)
|
||||
if (w == 0.0) then
|
||||
if (w == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
z0 = surf_p%coeffs(1)
|
||||
d = (z0 - z)/w
|
||||
if (d < 0.0) d = INFINITY
|
||||
if (d < ZERO) d = INFINITY
|
||||
end if
|
||||
|
||||
case (SURF_PLANE)
|
||||
|
|
@ -280,16 +286,16 @@ contains
|
|||
D = surf_p%coeffs(4)
|
||||
|
||||
tmp = A*u + B*v + C*w
|
||||
if (tmp == 0.0) then
|
||||
if (tmp == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = -(A*x + B*y + C*w - D)/tmp
|
||||
if (d < 0.0) d = INFINITY
|
||||
if (d < ZERO) d = INFINITY
|
||||
end if
|
||||
|
||||
case (SURF_CYL_X)
|
||||
a = 1.0 - u**2 ! v^2 + w^2
|
||||
if (a == 0.0) then
|
||||
a = ONE - u**2 ! v^2 + w^2
|
||||
if (a == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
y0 = surf_p%coeffs(1)
|
||||
|
|
@ -302,7 +308,7 @@ contains
|
|||
c = y**2 + z**2 - r**2
|
||||
quad = k**2 - a*c
|
||||
|
||||
if (c < 0) then
|
||||
if (c < ZERO) then
|
||||
! particle is inside the cylinder, thus one distance
|
||||
! must be negative and one must be positive. The
|
||||
! positive distance will be the one with negative sign
|
||||
|
|
@ -317,14 +323,14 @@ contains
|
|||
! positive sign on sqrt(quad)
|
||||
|
||||
d = -(k + sqrt(quad))/a
|
||||
if (d < 0) d = INFINITY
|
||||
if (d < ZERO) d = INFINITY
|
||||
|
||||
end if
|
||||
end if
|
||||
|
||||
case (SURF_CYL_Y)
|
||||
a = 1.0 - v**2 ! u^2 + w^2
|
||||
if (a == 0.0) then
|
||||
a = ONE - v**2 ! u^2 + w^2
|
||||
if (a == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
x0 = surf_p%coeffs(1)
|
||||
|
|
@ -337,7 +343,7 @@ contains
|
|||
c = x**2 + z**2 - r**2
|
||||
quad = k**2 - a*c
|
||||
|
||||
if (c < 0) then
|
||||
if (c < ZERO) then
|
||||
! particle is inside the cylinder, thus one distance
|
||||
! must be negative and one must be positive. The
|
||||
! positive distance will be the one with negative sign
|
||||
|
|
@ -352,14 +358,14 @@ contains
|
|||
! positive sign on sqrt(quad)
|
||||
|
||||
d = -(k + sqrt(quad))/a
|
||||
if (d < 0) d = INFINITY
|
||||
if (d < ZERO) d = INFINITY
|
||||
|
||||
end if
|
||||
end if
|
||||
|
||||
case (SURF_CYL_Z)
|
||||
a = 1.0 - w**2 ! u^2 + v^2
|
||||
if (a == 0.0) then
|
||||
a = ONE - w**2 ! u^2 + v^2
|
||||
if (a == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
x0 = surf_p%coeffs(1)
|
||||
|
|
@ -372,12 +378,12 @@ contains
|
|||
c = x**2 + y**2 - r**2
|
||||
quad = k**2 - a*c
|
||||
|
||||
if (quad < 0) then
|
||||
if (quad < ZERO) then
|
||||
! no intersection with cylinder
|
||||
|
||||
d = INFINITY
|
||||
|
||||
elseif (c < 0) then
|
||||
elseif (c < ZERO) then
|
||||
! particle is inside the cylinder, thus one distance
|
||||
! must be negative and one must be positive. The
|
||||
! positive distance will be the one with negative sign
|
||||
|
|
@ -392,7 +398,7 @@ contains
|
|||
! positive sign on sqrt(quad)
|
||||
|
||||
d = -(k + sqrt(quad))/a
|
||||
if (d < 0) d = INFINITY
|
||||
if (d < ZERO) d = INFINITY
|
||||
|
||||
end if
|
||||
end if
|
||||
|
|
@ -410,12 +416,12 @@ contains
|
|||
c = x**2 + y**2 + z**2 - r**2
|
||||
quad = k**2 - c
|
||||
|
||||
if (quad < 0) then
|
||||
if (quad < ZERO) then
|
||||
! no intersection with sphere
|
||||
|
||||
d = INFINITY
|
||||
|
||||
elseif (c < 0) then
|
||||
elseif (c < ZERO) then
|
||||
! particle is inside the sphere, thus one distance
|
||||
! must be negative and one must be positive. The
|
||||
! positive distance will be the one with negative sign
|
||||
|
|
@ -430,7 +436,7 @@ contains
|
|||
! positive sign on sqrt(quad)
|
||||
|
||||
d = -(k + sqrt(quad))
|
||||
if (d < 0) d = INFINITY
|
||||
if (d < ZERO) d = INFINITY
|
||||
|
||||
end if
|
||||
|
||||
|
|
|
|||
|
|
@ -56,13 +56,15 @@ module global
|
|||
|
||||
! Physical constants
|
||||
real(8), parameter :: &
|
||||
& PI = 2.*acos(0.0), & ! pi
|
||||
& MASS_NEUTRON = 1.0086649156, & ! mass of a neutron
|
||||
& MASS_PROTON = 1.00727646677, & ! mass of a proton
|
||||
& AMU = 1.66053873e-27, & ! 1 amu in kg
|
||||
& N_AVOGADRO = 0.602214179, & ! Avogadro's number in 10^24/mol
|
||||
& K_BOLTZMANN = 8.617342e-5, & ! Boltzmann constant in eV/K
|
||||
& INFINITY = huge(0.0_8) ! positive infinity
|
||||
& PI = 2.0_8*acos(0.0_8), & ! pi
|
||||
& MASS_NEUTRON = 1.0086649156, & ! mass of a neutron
|
||||
& MASS_PROTON = 1.00727646677, & ! mass of a proton
|
||||
& AMU = 1.66053873e-27, & ! 1 amu in kg
|
||||
& N_AVOGADRO = 0.602214179, & ! Avogadro's number in 10^24/mol
|
||||
& 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
|
||||
|
||||
! Boundary conditions
|
||||
integer, parameter :: &
|
||||
|
|
|
|||
|
|
@ -1,17 +1,20 @@
|
|||
# test input file
|
||||
|
||||
cell 100 40 -1
|
||||
cell 200 40 1 -2
|
||||
cell 200 41 1 -2
|
||||
|
||||
surface 1 sph 0 0 0 3
|
||||
surface 2 sph 0 0 0 5
|
||||
|
||||
material 40 -20.0 &
|
||||
3007.03c 1.0 &
|
||||
material 40 -4.5 &
|
||||
92238.03c 1.0
|
||||
|
||||
source box -3 -3 -3 3 3 3
|
||||
material 41 -1.0 &
|
||||
1001.03c 2.0 &
|
||||
8016.03c 1.0
|
||||
|
||||
source box -1 -1 -1 1 1 1
|
||||
|
||||
xs_library endfb7
|
||||
xs_data /opt/serpent/xsdata/endfb7
|
||||
criticality 1 1 50
|
||||
criticality 1 1 15
|
||||
|
|
|
|||
|
|
@ -20,7 +20,7 @@ program main
|
|||
|
||||
! Print the OpenMC title and version/date/time information
|
||||
call title()
|
||||
verbosity = 9
|
||||
verbosity = 10
|
||||
|
||||
! Initialize random number generator
|
||||
call RN_init_problem( 3, 0_8, 0_8, 0_8, 0 )
|
||||
|
|
@ -68,7 +68,9 @@ program main
|
|||
surfaces(2)%bc = BC_VACUUM
|
||||
call run_problem()
|
||||
|
||||
|
||||
! deallocate arrays
|
||||
call free_memory()
|
||||
|
||||
contains
|
||||
|
||||
!=====================================================================
|
||||
|
|
@ -115,8 +117,6 @@ contains
|
|||
|
||||
! print run time
|
||||
|
||||
call free_memory()
|
||||
|
||||
end subroutine run_problem
|
||||
|
||||
end program main
|
||||
|
|
|
|||
165
src/physics.f90
165
src/physics.f90
|
|
@ -6,6 +6,7 @@ module physics
|
|||
use mcnp_random, only: rang
|
||||
use output, only: error, message
|
||||
use search, only: binary_search
|
||||
use endf, only: reaction_name
|
||||
|
||||
implicit none
|
||||
|
||||
|
|
@ -29,6 +30,7 @@ contains
|
|||
real(8) :: distance
|
||||
real(8) :: Sigma ! total cross-section
|
||||
real(8) :: f ! interpolation factor
|
||||
real(8) :: tmp(3)
|
||||
integer :: IE ! index on energy grid
|
||||
|
||||
! determine what cell the particle is in
|
||||
|
|
@ -91,18 +93,14 @@ contains
|
|||
|
||||
type(Neutron), pointer :: neut
|
||||
|
||||
|
||||
type(AceContinuous), pointer :: table
|
||||
type(AceReaction), pointer :: rxn
|
||||
real(8) :: r1
|
||||
real(8) :: phi ! azimuthal angle
|
||||
real(8) :: mu ! cosine of polar angle
|
||||
character(250) :: msg
|
||||
integer :: cell_num
|
||||
integer :: i,j
|
||||
integer :: n_isotopes
|
||||
integer :: IE
|
||||
|
||||
real(8) :: f, Sigma
|
||||
real(8) :: f, Sigma, total
|
||||
real(8) :: density, density_i
|
||||
real(8) :: p
|
||||
real(8), allocatable :: Sigma_t(:)
|
||||
|
|
@ -127,40 +125,149 @@ contains
|
|||
Sigma_t(i) = Sigma
|
||||
end do
|
||||
|
||||
! normalize to create a discrete pdf
|
||||
Sigma_t = Sigma_t / sum(Sigma_t)
|
||||
|
||||
! sample nuclide
|
||||
r1 = rang()
|
||||
p = 0.0_8
|
||||
total = sum(Sigma_t)
|
||||
do i = 1, n_isotopes
|
||||
p = p + Sigma_t(i)
|
||||
p = p + Sigma_t(i) / total
|
||||
if (r1 < p) exit
|
||||
end do
|
||||
table => xs_continuous(cMaterial%table(i))
|
||||
! print *, 'sampled nuclide ', i
|
||||
|
||||
! select collision type
|
||||
r1 = rang()
|
||||
if (r1 <= 0.5) then
|
||||
! scatter
|
||||
phi = 2.*pi*rang()
|
||||
mu = 2.*rang() - 1
|
||||
neut%uvw(1) = mu
|
||||
neut%uvw(2) = sqrt(1. - mu**2) * cos(phi)
|
||||
neut%uvw(3) = sqrt(1. - mu**2) * sin(phi)
|
||||
else
|
||||
neut%alive = .false.
|
||||
if (verbosity >= 10) then
|
||||
cell_num = cells(neut%cell)%uid
|
||||
msg = " Absorbed in cell " // trim(int_to_str(cell_num))
|
||||
call message(msg, 10)
|
||||
end if
|
||||
return
|
||||
! Get table, total xs, interpolation factor
|
||||
table => xs_continuous(cMaterial%table(i))
|
||||
Sigma = Sigma_t(i)
|
||||
IE = table%grid_index(neut%IE)
|
||||
f = (neut%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
|
||||
density = cMaterial%atom_percent(i)*density
|
||||
|
||||
! free memory
|
||||
deallocate(Sigma_t)
|
||||
|
||||
! sample reaction channel
|
||||
r1 = rang()*Sigma
|
||||
p = 0.0_8
|
||||
do i = 1, table%n_reaction
|
||||
rxn => table%reactions(i)
|
||||
if (rxn%MT >= 200) cycle
|
||||
if (IE < rxn%IE) cycle
|
||||
p = p + density * (f*rxn%sigma(IE-rxn%IE+1) + (1-f)*(rxn%sigma(IE-rxn%IE+2)))
|
||||
if (r1 < p) exit
|
||||
end do
|
||||
if (verbosity >= 10) then
|
||||
msg = " " // trim(reaction_name(rxn%MT)) // " with nuclide " // &
|
||||
& trim(table%name)
|
||||
call message(msg, 10)
|
||||
end if
|
||||
|
||||
! call appropriate subroutine
|
||||
select case (rxn%MT)
|
||||
case (2)
|
||||
call elastic_scatter(neut, table%awr)
|
||||
case (102)
|
||||
call n_gamma(neut)
|
||||
case default
|
||||
call elastic_scatter(neut, table%awr)
|
||||
end select
|
||||
|
||||
end subroutine collision
|
||||
|
||||
!=====================================================================
|
||||
! ELASTIC_SCATTER
|
||||
!=====================================================================
|
||||
|
||||
subroutine elastic_scatter(neut, awr)
|
||||
|
||||
type(Neutron), pointer :: neut
|
||||
real(8), intent(in) :: awr
|
||||
|
||||
real(8) :: phi ! azimuthal angle
|
||||
real(8) :: mu ! cosine of polar angle
|
||||
real(8) :: vx, vy, vz
|
||||
real(8) :: vcx, vcy ,vcz
|
||||
real(8) :: vel
|
||||
real(8) :: u, v, w
|
||||
real(8) :: E
|
||||
integer :: IE
|
||||
|
||||
vel = sqrt(neut%E)
|
||||
|
||||
vx = vel*neut%uvw(1)
|
||||
vy = vel*neut%uvw(2)
|
||||
vz = vel*neut%uvw(3)
|
||||
|
||||
|
||||
vcx = vx/(awr + 1.0_8)
|
||||
vcy = vy/(awr + 1.0_8)
|
||||
vcz = vz/(awr + 1.0_8)
|
||||
|
||||
! Transform to CM frame
|
||||
vx = vx - vcx
|
||||
vy = vy - vcy
|
||||
vz = vz - vcz
|
||||
|
||||
vel = sqrt(vx*vx + vy*vy + vz*vz)
|
||||
|
||||
! Select isotropic direcion -- this is only valid for s-wave
|
||||
! scattering
|
||||
phi = 2.0_8*PI*rang()
|
||||
mu = 2.0_8*rang() - 1.0_8
|
||||
u = mu
|
||||
v = sqrt(1.0_8 - mu**2) * cos(phi)
|
||||
w = sqrt(1.0_8 - mu**2) * sin(phi)
|
||||
|
||||
vx = u*vel
|
||||
vy = v*vel
|
||||
vz = w*vel
|
||||
|
||||
! Transform back to LAB frame
|
||||
vx = vx + vcx
|
||||
vy = vy + vcy
|
||||
vz = vz + vcz
|
||||
|
||||
E = vx*vx + vy*vy + vz*vz
|
||||
vel = sqrt(E)
|
||||
|
||||
neut%E = E
|
||||
neut%uvw(1) = vx/vel
|
||||
neut%uvw(2) = vy/vel
|
||||
neut%uvw(3) = vz/vel
|
||||
|
||||
! find energy index, interpolation factor
|
||||
IE = binary_search(e_grid, n_grid, E)
|
||||
neut%IE = IE
|
||||
neut%interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE))
|
||||
|
||||
end subroutine elastic_scatter
|
||||
|
||||
!=====================================================================
|
||||
! LEVEL_INELASTIC
|
||||
!=====================================================================
|
||||
|
||||
subroutine level_inelastic
|
||||
|
||||
end subroutine level_inelastic
|
||||
|
||||
!=====================================================================
|
||||
! N_GAMMA
|
||||
!=====================================================================
|
||||
|
||||
subroutine n_gamma(neut)
|
||||
|
||||
type(Neutron), pointer :: neut
|
||||
|
||||
integer :: cell_num
|
||||
character(250) :: msg
|
||||
|
||||
neut%alive = .false.
|
||||
if (verbosity >= 10) then
|
||||
cell_num = cells(neut%cell)%uid
|
||||
msg = " Absorbed in cell " // trim(int_to_str(cell_num))
|
||||
call message(msg, 10)
|
||||
end if
|
||||
|
||||
end subroutine n_gamma
|
||||
|
||||
!=====================================================================
|
||||
! MAXWELL_SPECTRUM samples an energy from the Maxwell fission
|
||||
! distribution based on a rejection sampling scheme. This is described
|
||||
|
|
|
|||
|
|
@ -107,7 +107,7 @@ module types
|
|||
integer :: MT
|
||||
real(8) :: Q_value
|
||||
real(8) :: TY
|
||||
integer :: energy_index
|
||||
integer :: IE
|
||||
real(8), allocatable :: sigma(:)
|
||||
real(8), allocatable :: ang_cos(:,:)
|
||||
real(8), allocatable :: ang_pdf(:,:)
|
||||
|
|
@ -145,7 +145,8 @@ module types
|
|||
real(8), allocatable :: nu_d_precursor_const(:,:)
|
||||
real(8), allocatable :: nu_d_precursor_energy(:,:)
|
||||
real(8), allocatable :: nu_d_precursor_prob(:,:)
|
||||
type(AceReaction), pointer :: reactions(:)
|
||||
integer :: n_reaction
|
||||
type(AceReaction), pointer :: reactions(:) => null()
|
||||
|
||||
end type AceContinuous
|
||||
|
||||
|
|
@ -214,7 +215,7 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type ListKeyValueCI
|
||||
type(ListKeyValueCI), pointer :: next
|
||||
type(ListKeyValueCI), pointer :: next => null()
|
||||
type(KeyValueCI) :: data
|
||||
end type ListKeyValueCI
|
||||
|
||||
|
|
@ -224,7 +225,7 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type ListKeyValueII
|
||||
type(ListKeyValueII), pointer :: next
|
||||
type(ListKeyValueII), pointer :: next => null()
|
||||
type(KeyValueII) :: data
|
||||
end type ListKeyValueII
|
||||
|
||||
|
|
@ -234,7 +235,7 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type ListReal
|
||||
type(ListReal), pointer :: next
|
||||
type(ListReal), pointer :: next => null()
|
||||
real(8) :: data
|
||||
end type ListReal
|
||||
|
||||
|
|
@ -243,7 +244,7 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type ListInt
|
||||
type(ListInt), pointer :: next
|
||||
type(ListInt), pointer :: next => null()
|
||||
integer :: data
|
||||
end type ListInt
|
||||
|
||||
|
|
@ -253,7 +254,7 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type HashListCI
|
||||
type(ListKeyValueCI), pointer :: list
|
||||
type(ListKeyValueCI), pointer :: list => null()
|
||||
end type HashListCI
|
||||
|
||||
!=====================================================================
|
||||
|
|
@ -262,7 +263,7 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type HashListII
|
||||
type(ListKeyValueII), pointer :: list
|
||||
type(ListKeyValueII), pointer :: list => null()
|
||||
end type HashListII
|
||||
|
||||
!=====================================================================
|
||||
|
|
@ -271,7 +272,7 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type DictionaryCI
|
||||
type(HashListCI), pointer :: table(:)
|
||||
type(HashListCI), pointer :: table(:) => null()
|
||||
end type DictionaryCI
|
||||
|
||||
!=====================================================================
|
||||
|
|
@ -280,7 +281,7 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type DictionaryII
|
||||
type(HashListII), pointer :: table(:)
|
||||
type(HashListII), pointer :: table(:) => null()
|
||||
end type DictionaryII
|
||||
|
||||
end module types
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue