Merge pull request #1022 from paulromano/cmake-updates

Improvements to CMakeLists.txt
This commit is contained in:
Sterling Harper 2018-07-02 14:31:13 -04:00 committed by GitHub
commit 499f2f2494
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
17 changed files with 110 additions and 415 deletions

View file

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
project(openmc Fortran C CXX)
# Setup output directories
@ -9,14 +9,6 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
# Set module path
set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules)
#===============================================================================
# Architecture specific definitions
#===============================================================================
if (${UNIX})
add_definitions(-DUNIX)
endif()
#===============================================================================
# Command line options
#===============================================================================
@ -27,10 +19,7 @@ option(debug "Compile with debug flags" OFF)
option(optimize "Turn on all compiler optimization flags" OFF)
option(coverage "Compile with coverage analysis flags" OFF)
option(mpif08 "Use Fortran 2008 MPI interface" OFF)
# Maximum number of nested coordinates levels
set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels")
add_definitions(-DMAX_COORD=${maxcoord})
#===============================================================================
# MPI for distributed-memory parallelism
@ -39,17 +28,12 @@ add_definitions(-DMAX_COORD=${maxcoord})
set(MPI_ENABLED FALSE)
if($ENV{FC} MATCHES "(mpi[^/]*|ftn)$")
message("-- Detected MPI wrapper: $ENV{FC}")
add_definitions(-DOPENMC_MPI)
set(MPI_ENABLED TRUE)
# Get directory containing MPI wrapper
get_filename_component(MPI_DIR $ENV{FC} DIRECTORY)
endif()
# Check for Fortran 2008 MPI interface
if(MPI_ENABLED AND mpif08)
message("-- Using Fortran 2008 MPI bindings")
add_definitions(-DOPENMC_MPIF08)
endif()
#===============================================================================
@ -77,7 +61,6 @@ if(HDF5_IS_PARALLEL)
if(NOT MPI_ENABLED)
message(FATAL_ERROR "Parallel HDF5 must be used with MPI.")
endif()
add_definitions(-DPHDF5)
message("-- Using parallel HDF5")
endif()
@ -85,19 +68,15 @@ endif()
# Set compile/link flags based on which compiler is being used
#===============================================================================
# Support for Fortran in FindOpenMP was added in CMake 3.1. To support lower
# versions, we manually add the flags. However, at some point in time, the
# manual logic can be removed in favor of the block below
#if(NOT (CMAKE_VERSION VERSION_LESS 3.1))
# if(openmp)
# find_package(OpenMP)
# if(OPENMP_FOUND)
# list(APPEND f90flags ${OpenMP_Fortran_FLAGS})
# list(APPEND ldflags ${OpenMP_Fortran_FLAGS})
# endif()
# endif()
#endif()
if(openmp)
# Requires CMake 3.1+
find_package(OpenMP)
if(OPENMP_FOUND)
list(APPEND f90flags ${OpenMP_Fortran_FLAGS})
list(APPEND cxxflags ${OpenMP_CXX_FLAGS})
list(APPEND ldflags ${OpenMP_Fortran_FLAGS})
endif()
endif()
set(CMAKE_POSITION_INDEPENDENT_CODE ON)
@ -125,10 +104,6 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
list(REMOVE_ITEM f90flags -O2)
list(APPEND f90flags -O3)
endif()
if(openmp)
list(APPEND f90flags -fopenmp)
list(APPEND ldflags -fopenmp)
endif()
if(coverage)
list(APPEND f90flags -coverage)
list(APPEND ldflags -coverage)
@ -149,10 +124,6 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel)
if(optimize)
list(APPEND f90flags -O3)
endif()
if(openmp)
list(APPEND f90flags -qopenmp)
list(APPEND ldflags -qopenmp)
endif()
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL PGI)
# PGI Fortran compiler options
@ -187,10 +158,6 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL XL)
list(REMOVE_ITEM f90flags -O2)
list(APPEND f90flags -O3)
endif()
if(openmp)
list(APPEND f90flags -qsmp=omp)
list(APPEND ldflags -qsmp=omp)
endif()
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Cray)
# Cray Fortran compiler options
@ -204,7 +171,7 @@ endif()
if(CMAKE_C_COMPILER_ID STREQUAL GNU)
# GCC compiler options
list(APPEND cflags -std=c99 -O2)
list(APPEND cflags -O2)
if(debug)
list(REMOVE_ITEM cflags -O2)
list(APPEND cflags -g -Wall -pedantic -fbounds-check)
@ -222,7 +189,6 @@ if(CMAKE_C_COMPILER_ID STREQUAL GNU)
elseif(CMAKE_C_COMPILER_ID STREQUAL Intel)
# Intel compiler options
list(APPEND cflags -std=c99)
if(debug)
list(APPEND cflags -g -w3 -ftrapuv -fp-stack-check -O0)
endif()
@ -235,7 +201,6 @@ elseif(CMAKE_C_COMPILER_ID STREQUAL Intel)
elseif(CMAKE_C_COMPILER_ID MATCHES Clang)
# Clang options
list(APPEND cflags -std=c99)
if(debug)
list(APPEND cflags -g -O0 -ftrapv)
endif()
@ -257,9 +222,6 @@ if(optimize)
list(REMOVE_ITEM cxxflags -O2)
list(APPEND cxxflags -O3)
endif()
if(openmp)
list(APPEND cxxflags -fopenmp)
endif()
# Show flags being used
message(STATUS "Fortran flags: ${f90flags}")
@ -267,26 +229,12 @@ message(STATUS "C flags: ${cflags}")
message(STATUS "C++ flags: ${cxxflags}")
message(STATUS "Linker flags: ${ldflags}")
#===============================================================================
# git SHA1 hash
#===============================================================================
execute_process(COMMAND git rev-parse HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
RESULT_VARIABLE GIT_SHA1_SUCCESS
OUTPUT_VARIABLE GIT_SHA1
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
if(GIT_SHA1_SUCCESS EQUAL 0)
add_definitions(-DGIT_SHA1="${GIT_SHA1}")
endif()
#===============================================================================
# pugixml library
#===============================================================================
add_library(pugixml src/pugixml/pugixml_c.cpp src/pugixml/pugixml.cpp)
add_library(pugixml_fortran src/pugixml/pugixml_f.F90)
target_link_libraries(pugixml_fortran pugixml)
add_library(pugixml vendor/pugixml/pugixml.cpp)
target_include_directories(pugixml PUBLIC vendor/pugixml/)
#===============================================================================
# RPATH information
@ -295,7 +243,7 @@ target_link_libraries(pugixml_fortran pugixml)
# This block of code ensures that dynamic libraries can be found via the RPATH
# whether the executable is the original one from the build directory or the
# installed one in CMAKE_INSTALL_PREFIX. Ref:
# https://cmake.org/Wiki/CMake_RPATH_handling#Always_full_RPATH
# https://gitlab.kitware.com/cmake/community/wikis/doc/cmake/RPATH-handling
# use, i.e. don't skip the full RPATH for the build tree
set(CMAKE_SKIP_BUILD_RPATH FALSE)
@ -317,17 +265,20 @@ if("${isSystemDir}" STREQUAL "-1")
endif()
#===============================================================================
# Build faddeeva library
# faddeeva library
#===============================================================================
add_library(faddeeva STATIC src/faddeeva/Faddeeva.c)
add_library(faddeeva STATIC vendor/faddeeva/Faddeeva.c)
target_compile_options(faddeeva PRIVATE ${cflags})
set_target_properties(faddeeva PROPERTIES
C_STANDARD 99
C_STANDARD_REQUIRED ON)
#===============================================================================
# List source files. Define the libopenmc and the OpenMC executable
# libopenmc
#===============================================================================
set(program "openmc")
set(LIBOPENMC_FORTRAN_SRC
add_library(libopenmc SHARED
src/algorithm.F90
src/angle_distribution.F90
src/angleenergy_header.F90
@ -344,7 +295,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/dict_header.F90
src/distribution_multivariate.F90
src/distribution_univariate.F90
src/doppler.F90
src/eigenvalue.F90
src/endf.F90
src/endf_header.F90
@ -376,6 +326,7 @@ set(LIBOPENMC_FORTRAN_SRC
src/plot_header.F90
src/product_header.F90
src/progress_header.F90
src/pugixml/pugixml_f.F90
src/random_lcg.F90
src/reaction_header.F90
src/relaxng
@ -429,8 +380,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/tallies/tally_header.F90
src/tallies/trigger.F90
src/tallies/trigger_header.F90
)
set(LIBOPENMC_CXX_SRC
src/cell.cpp
src/initialize.cpp
src/finalize.cpp
@ -442,6 +391,7 @@ set(LIBOPENMC_CXX_SRC
src/mgxs.cpp
src/mgxs_interface.cpp
src/plot.cpp
src/pugixml/pugixml_c.cpp
src/random_lcg.cpp
src/scattdata.cpp
src/simulation.cpp
@ -449,48 +399,60 @@ set(LIBOPENMC_CXX_SRC
src/string_functions.cpp
src/surface.cpp
src/xml_interface.cpp
src/xsdata.cpp
src/pugixml/pugixml.cpp)
add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC})
src/xsdata.cpp)
set_target_properties(libopenmc PROPERTIES
OUTPUT_NAME openmc
PUBLIC_HEADER include/openmc.h)
add_executable(${program} src/main.cpp)
PUBLIC_HEADER include/openmc.h
LINKER_LANGUAGE Fortran)
#===============================================================================
# Add compiler/linker flags
#===============================================================================
set_property(TARGET ${program} libopenmc pugixml_fortran
PROPERTY LINKER_LANGUAGE Fortran)
target_include_directories(libopenmc PUBLIC include ${HDF5_INCLUDE_DIRS})
# The executable and the faddeeva package use only one language. They can be
# set via target_compile_options which accepts a list.
target_compile_options(${program} PUBLIC ${cxxflags})
target_compile_options(faddeeva PRIVATE ${cflags})
target_include_directories(libopenmc
PUBLIC include
PRIVATE ${HDF5_INCLUDE_DIRS})
# The libopenmc library has both F90 and C++ so the compile flags must be set
# file-by-file via set_source_file_properties. The compile flags must first be
# converted from lists to strings.
string(REPLACE ";" " " f90flags "${f90flags}")
string(REPLACE ";" " " cxxflags "${cxxflags}")
set_source_files_properties(${LIBOPENMC_FORTRAN_SRC} PROPERTIES COMPILE_FLAGS
${f90flags})
set_source_files_properties(${LIBOPENMC_CXX_SRC} PROPERTIES COMPILE_FLAGS
${cxxflags})
# differently depending on the language. The $<COMPILE_LANGUAGE> generator
# expression was added in CMake 3.3
target_compile_options(libopenmc PRIVATE
$<$<COMPILE_LANGUAGE:Fortran>:${f90flags}>
$<$<COMPILE_LANGUAGE:CXX>:${cxxflags}>)
# Add HDF5 library directories to link line with -L
foreach(LIBDIR ${HDF5_LIBRARY_DIRS})
list(APPEND ldflags "-L${LIBDIR}")
endforeach()
target_compile_definitions(libopenmc PRIVATE -DMAX_COORD=${maxcoord})
if (UNIX)
# Used in progress_header.F90 for calling check_isatty
target_compile_definitions(libopenmc PRIVATE -DUNIX)
endif()
if (HDF5_IS_PARALLEL)
target_compile_definitions(libopenmc PRIVATE -DPHDF5)
endif()
if (MPI_ENABLED)
target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI)
if (mpif08)
target_compile_definitions(libopenmc PRIVATE -DOPENMC_MPIF08)
endif()
endif()
# Set git SHA1 hash as a compile definition
execute_process(COMMAND git rev-parse HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
RESULT_VARIABLE GIT_SHA1_SUCCESS
OUTPUT_VARIABLE GIT_SHA1
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
if(GIT_SHA1_SUCCESS EQUAL 0)
target_compile_definitions(libopenmc PRIVATE -DGIT_SHA1="${GIT_SHA1}")
endif()
# target_link_libraries treats any arguments starting with - but not -l as
# linker flags. Thus, we can pass both linker flags and libraries together.
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml_fortran
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml
faddeeva)
target_link_libraries(${program} ${ldflags} libopenmc)
#===============================================================================
# openmc executable
#===============================================================================
add_executable(openmc src/main.cpp)
target_compile_options(openmc PRIVATE ${cxxflags})
target_link_libraries(openmc libopenmc)
#===============================================================================
# Python package
@ -506,7 +468,7 @@ add_custom_command(TARGET libopenmc POST_BUILD
# Install executable, scripts, manpage, license
#===============================================================================
install(TARGETS ${program} libopenmc
install(TARGETS openmc libopenmc
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
@ -514,4 +476,4 @@ install(TARGETS ${program} libopenmc
)
install(DIRECTORY src/relaxng DESTINATION share/openmc)
install(FILES man/man1/openmc.1 DESTINATION share/man/man1)
install(FILES LICENSE DESTINATION "share/doc/${program}" RENAME copyright)
install(FILES LICENSE DESTINATION "share/doc/openmc" RENAME copyright)

View file

@ -213,7 +213,7 @@ The following tables show all valid scores:
+----------------------+---------------------------------------------------+
|nu-fission |Total production of neutrons due to fission. |
+----------------------+---------------------------------------------------+
|nu-scatter, |This score is similar in functionality to the |
|nu-scatter |This score is similar in functionality to the |
| |``scatter`` score except the total production of |
| |neutrons due to scattering is scored vice simply |
| |the scattering rate. This accounts for |

View file

@ -134,6 +134,7 @@ extern "C" {
extern int32_t n_lattices;
extern int32_t n_materials;
extern int32_t n_meshes;
extern int n_nuclides;
extern int64_t n_particles;
extern int32_t n_plots;
extern int32_t n_realizations;

View file

@ -24,40 +24,44 @@ from openmc.stats import Discrete, Tabular
_THERMAL_NAMES = {
'c_Al27': ('al', 'al27'),
'c_Be': ('be', 'be-metal'),
'c_BeO': ('beo'),
'c_Be_in_BeO': ('bebeo', 'be-o', 'be/o'),
'c_Al27': ('al', 'al27', 'al-27'),
'c_Be': ('be', 'be-metal', 'be-met'),
'c_BeO': ('beo',),
'c_Be_in_BeO': ('bebeo', 'be-beo', 'be-o', 'be/o'),
'c_C6H6': ('benz', 'c6h6'),
'c_C_in_SiC': ('csic',),
'c_Ca_in_CaH2': ('cah'),
'c_D_in_D2O': ('dd2o', 'hwtr', 'hw'),
'c_Fe56': ('fe', 'fe56'),
'c_C_in_SiC': ('csic', 'c-sic'),
'c_Ca_in_CaH2': ('cah',),
'c_D_in_D2O': ('dd2o', 'd-d2o', 'hwtr', 'hw'),
'c_Fe56': ('fe', 'fe56', 'fe-56'),
'c_Graphite': ('graph', 'grph', 'gr'),
'c_H_in_CaH2': ('hcah2'),
'c_H_in_CH2': ('hch2', 'poly', 'pol'),
'c_Graphite_10p': ('grph10',),
'c_Graphite_30p': ('grph30',),
'c_H_in_CaH2': ('hcah2',),
'c_H_in_CH2': ('hch2', 'poly', 'pol', 'h-poly'),
'c_H_in_CH4_liquid': ('lch4', 'lmeth'),
'c_H_in_CH4_solid': ('sch4', 'smeth'),
'c_H_in_H2O': ('hh2o', 'lwtr', 'lw'),
'c_H_in_H2O_solid': ('hice',),
'c_H_in_C5O2H8': ('lucite', 'c5o2h8'),
'c_H_in_YH2': ('hyh2'),
'c_H_in_ZrH': ('hzrh', 'h-zr', 'h/zr', 'hzr'),
'c_H_in_H2O': ('hh2o', 'h-h2o', 'lwtr', 'lw'),
'c_H_in_H2O_solid': ('hice', 'h-ice'),
'c_H_in_C5O2H8': ('lucite', 'c5o2h8', 'h-luci'),
'c_H_in_YH2': ('hyh2', 'h-yh2'),
'c_H_in_ZrH': ('hzrh', 'h-zrh', 'h-zr', 'h/zr', 'hzr'),
'c_Mg24': ('mg', 'mg24'),
'c_O_in_BeO': ('obeo', 'o-be', 'o/be'),
'c_O_in_D2O': ('od2o'),
'c_O_in_H2O_ice': ('oice'),
'c_O_in_UO2': ('ouo2', 'o2-u', 'o2/u'),
'c_ortho_D': ('orthod', 'dortho'),
'c_ortho_H': ('orthoh', 'hortho'),
'c_Si_in_SiC': ('sisic'),
'c_O_in_BeO': ('obeo', 'o-beo', 'o-be', 'o/be'),
'c_O_in_D2O': ('od2o', 'o-d2o'),
'c_O_in_H2O_ice': ('oice', 'o-ice'),
'c_O_in_UO2': ('ouo2', 'o-uo2', 'o2-u', 'o2/u'),
'c_N_in_UN': ('n-un',),
'c_ortho_D': ('orthod', 'orthoD', 'dortho'),
'c_ortho_H': ('orthoh', 'orthoH', 'hortho'),
'c_Si_in_SiC': ('sisic', 'si-sic'),
'c_SiO2_alpha': ('sio2', 'sio2a'),
'c_SiO2_beta': ('sio2b'),
'c_para_D': ('parad', 'dpara'),
'c_para_H': ('parah', 'hpara'),
'c_U_in_UO2': ('uuo2', 'u-o2', 'u/o2'),
'c_Y_in_YH2': ('yyh2'),
'c_Zr_in_ZrH': ('zrzrh', 'zr-h', 'zr/h')
'c_SiO2_beta': ('sio2b',),
'c_para_D': ('parad', 'paraD', 'dpara'),
'c_para_H': ('parah', 'paraH', 'hpara'),
'c_U_in_UN': ('u-un',),
'c_U_in_UO2': ('uuo2', 'u-uo2', 'u-o2', 'u/o2'),
'c_Y_in_YH2': ('yyh2', 'y-yh2'),
'c_Zr_in_ZrH': ('zrzrh', 'zr-zrh', 'zr-h', 'zr/h')
}

View file

@ -54,7 +54,7 @@ class LegendreFilter(ExpansionFilter):
r"""Score Legendre expansion moments up to specified order.
This filter allows scores to be multiplied by Legendre polynomials of the
change in particle angle ($\mu$) up to a user-specified order.
change in particle angle (:math:`\mu`) up to a user-specified order.
Parameters
----------

View file

@ -7,7 +7,7 @@
#include <vector>
#include "hdf5.h"
#include "pugixml/pugixml.hpp"
#include "pugixml.hpp"
namespace openmc {
@ -50,7 +50,7 @@ public:
//! A geometry primitive that links surfaces, universes, and materials
//==============================================================================
class Cell
class Cell
{
public:
int32_t id; //!< Unique ID

View file

@ -1,216 +0,0 @@
module doppler
use constants, only: ZERO, ONE, PI, K_BOLTZMANN
implicit none
real(8), parameter :: sqrt_pi_inv = ONE / sqrt(PI)
contains
!===============================================================================
! BROADEN takes a microscopic cross section at a temperature T_1 and Doppler
! broadens it to a higher temperature T_2 based on a method originally developed
! by Cullen and Weisbin (see "Exact Doppler Broadening of Tabulated Cross
! Sections," Nucl. Sci. Eng. 60, 199-229 (1976)). The only difference here is
! the F functions are evaluated based on complementary error functions rather
! than error functions as is done in the BROADR module of NJOY.
!===============================================================================
subroutine broaden(energy, xs, A_target, T, sigmaNew)
real(8), intent(in) :: energy(:) ! energy grid
real(8), intent(in) :: xs(:) ! unbroadened cross section
integer, intent(in) :: A_target ! mass number of target
real(8), intent(in) :: T ! temperature (difference)
real(8), intent(out) :: sigmaNew(:) ! broadened cross section
integer :: i, k ! loop indices
integer :: n ! number of energy points
real(8) :: F_a(0:4) ! F(a) functions as per C&W
real(8) :: F_b(0:4) ! F(b) functions as per C&W
real(8) :: H(0:4) ! H functions as per C&W
real(8), allocatable :: x(:) ! proportional to relative velocity
real(8) :: y ! proportional to neutron velocity
real(8) :: y_sq ! y**2
real(8) :: y_inv ! 1/y
real(8) :: y_inv_sq ! 1/y**2
real(8) :: alpha ! constant equal to A/kT
real(8) :: slope ! slope of xs between adjacent points
real(8) :: Ak, Bk ! coefficients at each point
real(8) :: a, b ! values of x(k)-y and x(k+1)-y
real(8) :: sigma ! broadened cross section at one point
! Determine alpha parameter -- have to convert k to eV/K
alpha = A_target/(K_BOLTZMANN * T)
! Allocate memory for x and assign values
n = size(energy)
allocate(x(n))
x = sqrt(alpha * energy)
! Loop over incoming neutron energies
ENERGY_NEUTRON: do i = 1, n
sigma = ZERO
y = x(i)
y_sq = y*y
y_inv = ONE / y
y_inv_sq = y_inv / y
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4
k = i
a = ZERO
call calculate_F(F_a, a)
do while (a >= -4.0 .and. k > 1)
! Move to next point
F_b = F_a
k = k - 1
a = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k+1) - xs(k)) / (x(k+1)**2 - x(k)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! =======================================================================
! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE
if (k == 1 .and. a >= -4.0) then
! Since x = 0, this implies that a = -y
F_b = F_a
a = -y
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma + xs(k)*x(k)*(y_inv_sq*H(1) + y_inv*H(0))
end if
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4
k = i
b = ZERO
call calculate_F(F_b, b)
do while (b <= 4.0 .and. k < n)
! Move to next point
F_a = F_b
k = k + 1
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_b, b)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! =======================================================================
! EXTEND CROSS SECTION TO INFINITY ASSUMING CONSTANT SHAPE
if (k == n .and. b <= 4.0) then
! Calculate F function at last energy point
a = x(k) - y
call calculate_F(F_a, a)
! Add contribution to broadened cross section
sigma = sigma + xs(k) * (y_inv_sq*F_a(2) + 2.0*y_inv*F_a(1) + F_a(0))
end if
! =======================================================================
! EVALUATE SECOND TERM FROM x(k) + y = 0 to +4
if (y <= 4.0) then
! Swap signs on y
y = -y
y_inv = -y_inv
k = 1
! Calculate a and b based on 0 and x(1)
a = -y
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
call calculate_F(F_b, b)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma - xs(k) * x(k) * (y_inv_sq*H(1) + y_inv*H(0))
! Now progress forward doing the remainder of the second term
do while (b <= 4.0)
! Move to next point
F_a = F_b
k = k + 1
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_b, b)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) &
+ y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma - Ak*(xs(k) - slope*x(k)**2) - slope*Bk
end do
end if
! Set broadened cross section
sigmaNew(i) = sigma
end do ENERGY_NEUTRON
end subroutine broaden
!===============================================================================
! CALCULATE_F evaluates the function:
!
! F(n,a) = 1/sqrt(pi)*int(z^n*exp(-z^2), z = a to infinity)
!
! The five values returned in a vector correspond to the integral for n = 0
! through 4. These functions are called over and over during the Doppler
! broadening routine.
!===============================================================================
subroutine calculate_F(F, a)
real(8), intent(inout) :: F(0:4)
real(8), intent(in) :: a
#ifndef NO_F2008
F(0) = 0.5*erfc(a)
#endif
F(1) = 0.5*sqrt_pi_inv*exp(-a*a)
F(2) = 0.5*F(0) + a*F(1)
F(3) = F(1)*(1.0 + a*a)
F(4) = 0.75*F(0) + F(1)*a*(1.5 + a*a)
end subroutine calculate_F
end module doppler

View file

@ -37,7 +37,6 @@ module endf_header
contains
procedure :: from_hdf5 => polynomial_from_hdf5
procedure :: evaluate => polynomial_evaluate
procedure :: from_ace => polynomial_from_ace
end type Polynomial
!===============================================================================
@ -52,7 +51,6 @@ module endf_header
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
contains
procedure :: from_ace => tabulated1d_from_ace
procedure :: from_hdf5 => tabulated1d_from_hdf5
procedure :: evaluate => tabulated1d_evaluate
end type Tabulated1D
@ -63,24 +61,6 @@ contains
! Polynomial implementation
!===============================================================================
subroutine polynomial_from_ace(this, xss, idx)
class(Polynomial), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nc ! number of coefficients (order - 1)
! Clear space
if (allocated(this % coef)) deallocate(this % coef)
! Determine number of coefficients
nc = nint(xss(idx))
! Allocate space for and read coefficients
allocate(this % coef(nc))
this % coef(:) = xss(idx + 1 : idx + nc)
end subroutine polynomial_from_ace
subroutine polynomial_from_hdf5(this, dset_id)
class(Polynomial), intent(inout) :: this
integer(HID_T), intent(in) :: dset_id
@ -111,42 +91,6 @@ contains
! Tabulated1D implementation
!===============================================================================
subroutine tabulated1d_from_ace(this, xss, idx)
class(Tabulated1D), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nr, ne
! Clear space
if (allocated(this % nbt)) deallocate(this % nbt)
if (allocated(this % int)) deallocate(this % int)
if (allocated(this % x)) deallocate(this % x)
if (allocated(this % y)) deallocate(this % y)
! Determine number of regions
nr = nint(xss(idx))
this % n_regions = nr
! Read interpolation region data
if (nr > 0) then
allocate(this % nbt(nr))
allocate(this % int(nr))
this % nbt(:) = nint(xss(idx + 1 : idx + nr))
this % int(:) = nint(xss(idx + nr + 1 : idx + 2*nr))
end if
! Determine number of pairs
ne = int(XSS(idx + 2*nr + 1))
this % n_pairs = ne
! Read (x,y) pairs
allocate(this % x(ne))
allocate(this % y(ne))
this % x(:) = xss(idx + 2*nr + 2 : idx + 2*nr + 1 + ne)
this % y(:) = xss(idx + 2*nr + 2 + ne : idx + 2*nr + 1 + 2*ne)
end subroutine tabulated1d_from_ace
subroutine tabulated1d_from_hdf5(this, dset_id)
class(Tabulated1D), intent(inout) :: this
integer(HID_T), intent(in) :: dset_id

View file

@ -9,7 +9,7 @@
#include "constants.h"
#include "hdf5.h"
#include "pugixml/pugixml.hpp"
#include "pugixml.hpp"
namespace openmc {

View file

@ -6,7 +6,7 @@
#include <string>
#include "hdf5.h"
#include "pugixml/pugixml.hpp"
#include "pugixml.hpp"
#include "constants.h"

View file

@ -4,7 +4,7 @@
#include <string>
#include <vector>
#include "pugixml/pugixml.hpp"
#include "pugixml.hpp"
namespace openmc {