diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index 000801ca2a..a6f90818f9 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -213,42 +213,42 @@ def get_opencg_surface(openmc_surface): opencg_surface = None if openmc_surface.type == 'plane': - A = openmc_surface.coeffs['A'] - B = openmc_surface.coeffs['B'] - C = openmc_surface.coeffs['C'] - D = openmc_surface.coeffs['D'] + A = openmc_surface.a + B = openmc_surface.b + C = openmc_surface.c + D = openmc_surface.d opencg_surface = opencg.Plane(surface_id, name, boundary, A, B, C, D) elif openmc_surface.type == 'x-plane': - x0 = openmc_surface.coeffs['x0'] + x0 = openmc_surface.y0 opencg_surface = opencg.XPlane(surface_id, name, boundary, x0) elif openmc_surface.type == 'y-plane': - y0 = openmc_surface.coeffs['y0'] + y0 = openmc_surface.y0 opencg_surface = opencg.YPlane(surface_id, name, boundary, y0) elif openmc_surface.type == 'z-plane': - z0 = openmc_surface.coeffs['z0'] + z0 = openmc_surface.z0 opencg_surface = opencg.ZPlane(surface_id, name, boundary, z0) elif openmc_surface.type == 'x-cylinder': - y0 = openmc_surface.coeffs['y0'] - z0 = openmc_surface.coeffs['z0'] - R = openmc_surface.coeffs['R'] + y0 = openmc_surface.y0 + z0 = openmc_surface.z0 + R = openmc_surface.r opencg_surface = opencg.XCylinder(surface_id, name, boundary, y0, z0, R) elif openmc_surface.type == 'y-cylinder': - x0 = openmc_surface.coeffs['x0'] - z0 = openmc_surface.coeffs['z0'] - R = openmc_surface.coeffs['R'] + x0 = openmc_surface.x0 + z0 = openmc_surface.z0 + R = openmc_surface.r opencg_surface = opencg.YCylinder(surface_id, name, boundary, x0, z0, R) elif openmc_surface.type == 'z-cylinder': - x0 = openmc_surface.coeffs['x0'] - y0 = openmc_surface.coeffs['y0'] - R = openmc_surface.coeffs['R'] + x0 = openmc_surface.x0 + y0 = openmc_surface.y0 + R = openmc_surface.r opencg_surface = opencg.ZCylinder(surface_id, name, boundary, x0, y0, R) @@ -297,40 +297,40 @@ def get_openmc_surface(opencg_surface): boundary = 'transmission' if opencg_surface.type == 'plane': - A = opencg_surface.coeffs['A'] - B = opencg_surface.coeffs['B'] - C = opencg_surface.coeffs['C'] - D = opencg_surface.coeffs['D'] + A = opencg_surface.a + B = opencg_surface.b + C = opencg_surface.c + D = opencg_surface.d openmc_surface = openmc.Plane(surface_id, boundary, A, B, C, D, name) elif opencg_surface.type == 'x-plane': - x0 = opencg_surface.coeffs['x0'] + x0 = opencg_surface.x0 openmc_surface = openmc.XPlane(surface_id, boundary, x0, name) elif opencg_surface.type == 'y-plane': - y0 = opencg_surface.coeffs['y0'] + y0 = opencg_surface.y0 openmc_surface = openmc.YPlane(surface_id, boundary, y0, name) elif opencg_surface.type == 'z-plane': - z0 = opencg_surface.coeffs['z0'] + z0 = opencg_surface.z0 openmc_surface = openmc.ZPlane(surface_id, boundary, z0, name) elif opencg_surface.type == 'x-cylinder': - y0 = opencg_surface.coeffs['y0'] - z0 = opencg_surface.coeffs['z0'] - R = opencg_surface.coeffs['R'] + y0 = opencg_surface.y0 + z0 = opencg_surface.z0 + R = opencg_surface.r openmc_surface = openmc.XCylinder(surface_id, boundary, y0, z0, R, name) elif opencg_surface.type == 'y-cylinder': - x0 = opencg_surface.coeffs['x0'] - z0 = opencg_surface.coeffs['z0'] - R = opencg_surface.coeffs['R'] + x0 = opencg_surface.x0 + z0 = opencg_surface.z0 + R = opencg_surface.r openmc_surface = openmc.YCylinder(surface_id, boundary, x0, z0, R, name) elif opencg_surface.type == 'z-cylinder': - x0 = opencg_surface.coeffs['x0'] - y0 = opencg_surface.coeffs['y0'] - R = opencg_surface.coeffs['R'] + x0 = opencg_surface.x0 + y0 = opencg_surface.y0 + R = opencg_surface.r openmc_surface = openmc.ZCylinder(surface_id, boundary, x0, y0, R, name) else: @@ -384,9 +384,9 @@ def get_compatible_opencg_surfaces(opencg_surface): boundary = opencg_surface.boundary_type if opencg_surface.type == 'x-squareprism': - y0 = opencg_surface.coeffs['y0'] - z0 = opencg_surface.coeffs['z0'] - R = opencg_surface.coeffs['R'] + y0 = opencg_surface.y0 + z0 = opencg_surface.z0 + R = opencg_surface.r # Create a list of the four planes we need 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] elif opencg_surface.type == 'y-squareprism': - x0 = opencg_surface.coeffs['x0'] - z0 = opencg_surface.coeffs['z0'] - R = opencg_surface.coeffs['R'] + x0 = opencg_surface.x0 + z0 = opencg_surface.z0 + R = opencg_surface.r # Create a list of the four planes we need 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] elif opencg_surface.type == 'z-squareprism': - x0 = opencg_surface.coeffs['x0'] - y0 = opencg_surface.coeffs['y0'] - R = opencg_surface.coeffs['R'] + x0 = opencg_surface.x0['x0'] + y0 = opencg_surface.y0['y0'] + R = opencg_surface.r['R'] # Create a list of the four planes we need left = opencg.XPlane(name=name, boundary=boundary, x0=x0-R) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b937b03a15..4d8fb2f0fb 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -33,6 +33,7 @@ contains integer :: i_nuclide ! index into nuclides array integer :: i_sab ! index into sab_tables array integer :: j ! index in mat % i_sab_nuclides + integer :: u ! index into logarithmic mapping array real(8) :: atom_density ! atom density of a nuclide logical :: check_sab ! should we check for S(a,b) table? type(Material), pointer :: mat ! current material @@ -50,9 +51,13 @@ contains mat => materials(p % material) - ! Find energy index on global or material unionized grid - if (grid_method == GRID_MAT_UNION) & - call find_energy_index(p % E, p % material) + ! Find energy index on energy grid + u = 0 + 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 check_sab = (mat % n_sab > 0) @@ -94,9 +99,9 @@ contains ! Calculate microscopic cross section for this nuclide 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 - 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 ! ======================================================================== @@ -137,16 +142,16 @@ contains ! 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_sab ! index into sab_tables 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) :: u ! index into logarithmic mapping array integer :: i_grid ! index on nuclide energy grid integer :: i_low ! lower 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) :: f ! interp factor on nuclide energy grid type(Nuclide), pointer :: nuc @@ -173,7 +178,6 @@ contains else ! Determine bounding indices based on which equal log-spaced interval ! the energy is in - u = int(log(E/1.0e-11_8)/log_spacing) i_low = nuc % grid_index(u) i_high = nuc % grid_index(u + 1) + 1 diff --git a/src/search.F90 b/src/search.F90 index dab7fa67ca..d38dfb986e 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -28,7 +28,6 @@ contains integer :: L integer :: R integer :: n_iteration - real(8) :: testval L = 1 R = n @@ -39,22 +38,11 @@ contains n_iteration = 0 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 array_index = L + (R - L)/2 - testval = array(array_index) - if (val >= testval) then + if (val >= array(array_index)) then L = array_index - elseif (val < testval) then + else R = array_index end if @@ -80,7 +68,6 @@ contains integer :: L integer :: R integer :: n_iteration - real(8) :: testval L = 1 R = n @@ -91,22 +78,11 @@ contains n_iteration = 0 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 array_index = L + (R - L)/2 - testval = array(array_index) - if (val >= testval) then + if (val >= array(array_index)) then L = array_index - elseif (val < testval) then + else R = array_index end if @@ -132,7 +108,6 @@ contains integer :: L integer :: R integer :: n_iteration - real(8) :: testval L = 1 R = n @@ -143,22 +118,11 @@ contains n_iteration = 0 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 array_index = L + (R - L)/2 - testval = array(array_index) - if (val >= testval) then + if (val >= array(array_index)) then L = array_index - elseif (val < testval) then + else R = array_index end if diff --git a/tests/run_tests.py b/tests/run_tests.py index d3b79aa3b1..70ea4c3dc9 100755 --- a/tests/run_tests.py +++ b/tests/run_tests.py @@ -126,7 +126,12 @@ class Test(object): # Check for 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: self.fc = FC