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

Merge branch 'opencg-underscores' into mgxs

Browse source
This commit is contained in:
Will Boyd 2015-10-03 14:35:29 -04:00
commit 5f3f58f3ab
4 changed files with 65 additions and 92 deletions
Download patch file Download diff file Expand all files Collapse all files

82
openmc/opencg_compatible.py
View file

@ -213,42 +213,42 @@ def get_opencg_surface(openmc_surface):
opencg_surface = None opencg_surface = None
if openmc_surface.type == 'plane': if openmc_surface.type == 'plane':
A = openmc_surface.coeffs['A'] A = openmc_surface.a
B = openmc_surface.coeffs['B'] B = openmc_surface.b
C = openmc_surface.coeffs['C'] C = openmc_surface.c
D = openmc_surface.coeffs['D'] D = openmc_surface.d
opencg_surface = opencg.Plane(surface_id, name, boundary, A, B, C, D) opencg_surface = opencg.Plane(surface_id, name, boundary, A, B, C, D)
elif openmc_surface.type == 'x-plane': elif openmc_surface.type == 'x-plane':
x0 = openmc_surface.coeffs['x0'] x0 = openmc_surface.y0
opencg_surface = opencg.XPlane(surface_id, name, boundary, x0) opencg_surface = opencg.XPlane(surface_id, name, boundary, x0)
elif openmc_surface.type == 'y-plane': elif openmc_surface.type == 'y-plane':
y0 = openmc_surface.coeffs['y0'] y0 = openmc_surface.y0
opencg_surface = opencg.YPlane(surface_id, name, boundary, y0) opencg_surface = opencg.YPlane(surface_id, name, boundary, y0)
elif openmc_surface.type == 'z-plane': elif openmc_surface.type == 'z-plane':
z0 = openmc_surface.coeffs['z0'] z0 = openmc_surface.z0
opencg_surface = opencg.ZPlane(surface_id, name, boundary, z0) opencg_surface = opencg.ZPlane(surface_id, name, boundary, z0)
elif openmc_surface.type == 'x-cylinder': elif openmc_surface.type == 'x-cylinder':
y0 = openmc_surface.coeffs['y0'] y0 = openmc_surface.y0
z0 = openmc_surface.coeffs['z0'] z0 = openmc_surface.z0
R = openmc_surface.coeffs['R'] R = openmc_surface.r
opencg_surface = opencg.XCylinder(surface_id, name, opencg_surface = opencg.XCylinder(surface_id, name,
boundary, y0, z0, R) boundary, y0, z0, R)
elif openmc_surface.type == 'y-cylinder': elif openmc_surface.type == 'y-cylinder':
x0 = openmc_surface.coeffs['x0'] x0 = openmc_surface.x0
z0 = openmc_surface.coeffs['z0'] z0 = openmc_surface.z0
R = openmc_surface.coeffs['R'] R = openmc_surface.r
opencg_surface = opencg.YCylinder(surface_id, name, opencg_surface = opencg.YCylinder(surface_id, name,
boundary, x0, z0, R) boundary, x0, z0, R)
elif openmc_surface.type == 'z-cylinder': elif openmc_surface.type == 'z-cylinder':
x0 = openmc_surface.coeffs['x0'] x0 = openmc_surface.x0
y0 = openmc_surface.coeffs['y0'] y0 = openmc_surface.y0
R = openmc_surface.coeffs['R'] R = openmc_surface.r
opencg_surface = opencg.ZCylinder(surface_id, name, opencg_surface = opencg.ZCylinder(surface_id, name,
boundary, x0, y0, R) boundary, x0, y0, R)
@ -297,40 +297,40 @@ def get_openmc_surface(opencg_surface):
boundary = 'transmission' boundary = 'transmission'
if opencg_surface.type == 'plane': if opencg_surface.type == 'plane':
A = opencg_surface.coeffs['A'] A = opencg_surface.a
B = opencg_surface.coeffs['B'] B = opencg_surface.b
C = opencg_surface.coeffs['C'] C = opencg_surface.c
D = opencg_surface.coeffs['D'] D = opencg_surface.d
openmc_surface = openmc.Plane(surface_id, boundary, A, B, C, D, name) openmc_surface = openmc.Plane(surface_id, boundary, A, B, C, D, name)
elif opencg_surface.type == 'x-plane': elif opencg_surface.type == 'x-plane':
x0 = opencg_surface.coeffs['x0'] x0 = opencg_surface.x0
openmc_surface = openmc.XPlane(surface_id, boundary, x0, name) openmc_surface = openmc.XPlane(surface_id, boundary, x0, name)
elif opencg_surface.type == 'y-plane': elif opencg_surface.type == 'y-plane':
y0 = opencg_surface.coeffs['y0'] y0 = opencg_surface.y0
openmc_surface = openmc.YPlane(surface_id, boundary, y0, name) openmc_surface = openmc.YPlane(surface_id, boundary, y0, name)
elif opencg_surface.type == 'z-plane': elif opencg_surface.type == 'z-plane':
z0 = opencg_surface.coeffs['z0'] z0 = opencg_surface.z0
openmc_surface = openmc.ZPlane(surface_id, boundary, z0, name) openmc_surface = openmc.ZPlane(surface_id, boundary, z0, name)
elif opencg_surface.type == 'x-cylinder': elif opencg_surface.type == 'x-cylinder':
y0 = opencg_surface.coeffs['y0'] y0 = opencg_surface.y0
z0 = opencg_surface.coeffs['z0'] z0 = opencg_surface.z0
R = opencg_surface.coeffs['R'] R = opencg_surface.r
openmc_surface = openmc.XCylinder(surface_id, boundary, y0, z0, R, name) openmc_surface = openmc.XCylinder(surface_id, boundary, y0, z0, R, name)
elif opencg_surface.type == 'y-cylinder': elif opencg_surface.type == 'y-cylinder':
x0 = opencg_surface.coeffs['x0'] x0 = opencg_surface.x0
z0 = opencg_surface.coeffs['z0'] z0 = opencg_surface.z0
R = opencg_surface.coeffs['R'] R = opencg_surface.r
openmc_surface = openmc.YCylinder(surface_id, boundary, x0, z0, R, name) openmc_surface = openmc.YCylinder(surface_id, boundary, x0, z0, R, name)
elif opencg_surface.type == 'z-cylinder': elif opencg_surface.type == 'z-cylinder':
x0 = opencg_surface.coeffs['x0'] x0 = opencg_surface.x0
y0 = opencg_surface.coeffs['y0'] y0 = opencg_surface.y0
R = opencg_surface.coeffs['R'] R = opencg_surface.r
openmc_surface = openmc.ZCylinder(surface_id, boundary, x0, y0, R, name) openmc_surface = openmc.ZCylinder(surface_id, boundary, x0, y0, R, name)
else: else:
@ -384,9 +384,9 @@ def get_compatible_opencg_surfaces(opencg_surface):
boundary = opencg_surface.boundary_type boundary = opencg_surface.boundary_type
if opencg_surface.type == 'x-squareprism': if opencg_surface.type == 'x-squareprism':
y0 = opencg_surface.coeffs['y0'] y0 = opencg_surface.y0
z0 = opencg_surface.coeffs['z0'] z0 = opencg_surface.z0
R = opencg_surface.coeffs['R'] R = opencg_surface.r
# Create a list of the four planes we need # Create a list of the four planes we need
left = opencg.YPlane(name=name, boundary=boundary, y0=y0-R) left = opencg.YPlane(name=name, boundary=boundary, y0=y0-R)
@ -396,9 +396,9 @@ def get_compatible_opencg_surfaces(opencg_surface):
surfaces = [left, right, bottom, top] surfaces = [left, right, bottom, top]
elif opencg_surface.type == 'y-squareprism': elif opencg_surface.type == 'y-squareprism':
x0 = opencg_surface.coeffs['x0'] x0 = opencg_surface.x0
z0 = opencg_surface.coeffs['z0'] z0 = opencg_surface.z0
R = opencg_surface.coeffs['R'] R = opencg_surface.r
# Create a list of the four planes we need # Create a list of the four planes we need
left = opencg.XPlane(name=name, boundary=boundary, x0=x0-R) left = opencg.XPlane(name=name, boundary=boundary, x0=x0-R)
@ -408,9 +408,9 @@ def get_compatible_opencg_surfaces(opencg_surface):
surfaces = [left, right, bottom, top] surfaces = [left, right, bottom, top]
elif opencg_surface.type == 'z-squareprism': elif opencg_surface.type == 'z-squareprism':
x0 = opencg_surface.coeffs['x0'] x0 = opencg_surface.x0['x0']
y0 = opencg_surface.coeffs['y0'] y0 = opencg_surface.y0['y0']
R = opencg_surface.coeffs['R'] R = opencg_surface.r['R']
# Create a list of the four planes we need # Create a list of the four planes we need
left = opencg.XPlane(name=name, boundary=boundary, x0=x0-R) left = opencg.XPlane(name=name, boundary=boundary, x0=x0-R)

20
src/cross_section.F90
View file

@ -33,6 +33,7 @@ contains
integer :: i_nuclide ! index into nuclides array integer :: i_nuclide ! index into nuclides array
integer :: i_sab ! index into sab_tables array integer :: i_sab ! index into sab_tables array
integer :: j ! index in mat % i_sab_nuclides integer :: j ! index in mat % i_sab_nuclides
integer :: u ! index into logarithmic mapping array
real(8) :: atom_density ! atom density of a nuclide real(8) :: atom_density ! atom density of a nuclide
logical :: check_sab ! should we check for S(a,b) table? logical :: check_sab ! should we check for S(a,b) table?
type(Material), pointer :: mat ! current material type(Material), pointer :: mat ! current material
@ -50,9 +51,13 @@ contains
mat => materials(p % material) mat => materials(p % material)
! Find energy index on global or material unionized grid ! Find energy index on energy grid
if (grid_method == GRID_MAT_UNION) & u = 0
call find_energy_index(p % E, p % material) if (grid_method == GRID_MAT_UNION) then
call find_energy_index(p % E, p % material)
else if (grid_method == GRID_LOGARITHM) then
u = int(log(p % E/1.0e-11_8)/log_spacing)
end if
! Determine if this material has S(a,b) tables ! Determine if this material has S(a,b) tables
check_sab = (mat % n_sab > 0) check_sab = (mat % n_sab > 0)
@ -94,9 +99,9 @@ contains
! Calculate microscopic cross section for this nuclide ! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E) then if (p % E /= micro_xs(i_nuclide) % last_E) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i) call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i) call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u)
end if end if
! ======================================================================== ! ========================================================================
@ -137,16 +142,16 @@ contains
! given index in the nuclides array at the energy of the given particle ! given index in the nuclides array at the energy of the given particle
!=============================================================================== !===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat) subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, u)
integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array integer, intent(in) :: i_sab ! index into sab_tables array
integer, intent(in) :: i_mat ! index into materials array integer, intent(in) :: i_mat ! index into materials array
integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material
integer, intent(in) :: u ! index into logarithmic mapping array
integer :: i_grid ! index on nuclide energy grid integer :: i_grid ! index on nuclide energy grid
integer :: i_low ! lower logarithmic mapping index integer :: i_low ! lower logarithmic mapping index
integer :: i_high ! upper logarithmic mapping index integer :: i_high ! upper logarithmic mapping index
integer :: u ! index into logarithmic mapping array
real(8), intent(in) :: E ! energy real(8), intent(in) :: E ! energy
real(8) :: f ! interp factor on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide), pointer :: nuc type(Nuclide), pointer :: nuc
@ -173,7 +178,6 @@ contains
else else
! Determine bounding indices based on which equal log-spaced interval ! Determine bounding indices based on which equal log-spaced interval
! the energy is in ! the energy is in
u = int(log(E/1.0e-11_8)/log_spacing)
i_low = nuc % grid_index(u) i_low = nuc % grid_index(u)
i_high = nuc % grid_index(u + 1) + 1 i_high = nuc % grid_index(u + 1) + 1

48
src/search.F90
View file

@ -28,7 +28,6 @@ contains
integer :: L integer :: L
integer :: R integer :: R
integer :: n_iteration integer :: n_iteration
real(8) :: testval
L = 1 L = 1
R = n R = n
@ -39,22 +38,11 @@ contains
n_iteration = 0 n_iteration = 0
do while (R - L > 1) do while (R - L > 1)
! Check boundaries
if (val > array(L) .and. val < array(L+1)) then
array_index = L
return
elseif (val > array(R-1) .and. val < array(R)) then
array_index = R - 1
return
end if
! Find values at midpoint ! Find values at midpoint
array_index = L + (R - L)/2 array_index = L + (R - L)/2
testval = array(array_index) if (val >= array(array_index)) then
if (val >= testval) then
L = array_index L = array_index
elseif (val < testval) then else
R = array_index R = array_index
end if end if
@ -80,7 +68,6 @@ contains
integer :: L integer :: L
integer :: R integer :: R
integer :: n_iteration integer :: n_iteration
real(8) :: testval
L = 1 L = 1
R = n R = n
@ -91,22 +78,11 @@ contains
n_iteration = 0 n_iteration = 0
do while (R - L > 1) do while (R - L > 1)
! Check boundaries
if (val > array(L) .and. val < array(L+1)) then
array_index = L
return
elseif (val > array(R-1) .and. val < array(R)) then
array_index = R - 1
return
end if
! Find values at midpoint ! Find values at midpoint
array_index = L + (R - L)/2 array_index = L + (R - L)/2
testval = array(array_index) if (val >= array(array_index)) then
if (val >= testval) then
L = array_index L = array_index
elseif (val < testval) then else
R = array_index R = array_index
end if end if
@ -132,7 +108,6 @@ contains
integer :: L integer :: L
integer :: R integer :: R
integer :: n_iteration integer :: n_iteration
real(8) :: testval
L = 1 L = 1
R = n R = n
@ -143,22 +118,11 @@ contains
n_iteration = 0 n_iteration = 0
do while (R - L > 1) do while (R - L > 1)
! Check boundaries
if (val > array(L) .and. val < array(L+1)) then
array_index = L
return
elseif (val > array(R-1) .and. val < array(R)) then
array_index = R - 1
return
end if
! Find values at midpoint ! Find values at midpoint
array_index = L + (R - L)/2 array_index = L + (R - L)/2
testval = array(array_index) if (val >= array(array_index)) then
if (val >= testval) then
L = array_index L = array_index
elseif (val < testval) then else
R = array_index R = array_index
end if end if

7
tests/run_tests.py
View file

@ -126,7 +126,12 @@ class Test(object):
# Check for MPI # Check for MPI
if self.mpi: if self.mpi:
self.fc = os.path.join(MPI_DIR, 'bin', 'mpifort') if os.path.exists(os.path.join(MPI_DIR, 'bin', 'mpifort')):
self.fc = os.path.join(MPI_DIR, 'bin', 'mpifort')
elif os.path.exists(os.path.join(MPI_DIR, 'bin', 'mpif90')):
self.fc = os.path.join(MPI_DIR, 'bin', 'mpif90')
else:
raise RuntimeError('Cannot find an MPI Fortran compiler')
else: else:
self.fc = FC self.fc = FC