mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge remote-tracking branch 'cjosey/multipole-final' into diff_tally5
This commit is contained in:
Sterling Harper
commit
9f52c707ff
42 changed files with 4081 additions and 129 deletions
|
|
@ -1,5 +1,5 @@
|
|||
cmake_minimum_required(VERSION 2.8 FATAL_ERROR)
|
||||
project(openmc Fortran)
|
||||
project(openmc Fortran C)
|
||||
|
||||
# Setup output directories
|
||||
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
|
||||
|
|
@ -107,49 +107,62 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
|
|||
message(FATAL_ERROR "gfortran version must be 4.6 or higher")
|
||||
endif()
|
||||
|
||||
# GNU Fortran compiler options
|
||||
# GCC compiler options
|
||||
list(APPEND f90flags -cpp -std=f2008 -fbacktrace)
|
||||
list(APPEND cflags -cpp -std=c99)
|
||||
if(debug)
|
||||
if(NOT (GCC_VERSION VERSION_LESS 4.7))
|
||||
list(APPEND f90flags -Wall)
|
||||
list(APPEND cflags -Wall)
|
||||
endif()
|
||||
list(APPEND f90flags -g -pedantic -fbounds-check
|
||||
-ffpe-trap=invalid,overflow,underflow)
|
||||
list(APPEND cflags -g -pedantic -fbounds-check)
|
||||
list(APPEND ldflags -g)
|
||||
endif()
|
||||
if(profile)
|
||||
list(APPEND f90flags -pg)
|
||||
list(APPEND cflags -pg)
|
||||
list(APPEND ldflags -pg)
|
||||
endif()
|
||||
if(optimize)
|
||||
list(APPEND f90flags -O3 -flto -fuse-linker-plugin)
|
||||
list(APPEND cflags -O3 -flto -fuse-linker-plugin)
|
||||
endif()
|
||||
if(openmp)
|
||||
list(APPEND f90flags -fopenmp)
|
||||
list(APPEND cflags -fopenmp)
|
||||
list(APPEND ldflags -fopenmp)
|
||||
endif()
|
||||
if(coverage)
|
||||
list(APPEND f90flags -coverage)
|
||||
list(APPEND cflags -coverage)
|
||||
list(APPEND ldflags -coverage)
|
||||
endif()
|
||||
|
||||
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel)
|
||||
# Intel Fortran compiler options
|
||||
# Intel compiler options
|
||||
list(APPEND f90flags -fpp -std08 -assume byterecl -traceback)
|
||||
list(APPEND cflags -std=c99)
|
||||
if(debug)
|
||||
list(APPEND f90flags -g -warn -ftrapuv -fp-stack-check
|
||||
"-check all" -fpe0)
|
||||
list(APPEND cflags -g -warn -ftrapuv -fp-stack-check
|
||||
"-check all" -fpe0)
|
||||
list(APPEND ldflags -g)
|
||||
endif()
|
||||
if(profile)
|
||||
list(APPEND f90flags -pg)
|
||||
list(APPEND cflags -pg)
|
||||
list(APPEND ldflags -pg)
|
||||
endif()
|
||||
if(optimize)
|
||||
list(APPEND f90flags -O3)
|
||||
list(APPEND cflags -O3)
|
||||
endif()
|
||||
if(openmp)
|
||||
list(APPEND f90flags -openmp)
|
||||
list(APPEND cflags -openmp)
|
||||
list(APPEND ldflags -openmp)
|
||||
endif()
|
||||
|
||||
|
|
@ -243,6 +256,12 @@ add_subdirectory(src/xml/fox)
|
|||
# which point to directories outside the build tree to the install RPATH
|
||||
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
|
||||
|
||||
#===============================================================================
|
||||
# Build faddeeva library
|
||||
#===============================================================================
|
||||
|
||||
add_library(faddeeva STATIC src/Faddeeva.c)
|
||||
|
||||
#===============================================================================
|
||||
# Build OpenMC executable
|
||||
#===============================================================================
|
||||
|
|
@ -263,11 +282,14 @@ endif()
|
|||
# set compile flags. Note that this sets the COMPILE_OPTIONS property (also
|
||||
# available only in 2.8.12+) rather than the COMPILE_FLAGS property, which is
|
||||
# deprecated. The former can handle lists whereas the latter cannot.
|
||||
if(CMAKE_VERSION VERSION_LESS 2.8.12)
|
||||
if (CMAKE_VERSION VERSION_LESS 2.8.12)
|
||||
string(REPLACE ";" " " f90flags "${f90flags}")
|
||||
string(REPLACE ";" " " cflags "${cflags}")
|
||||
set_property(TARGET ${program} PROPERTY COMPILE_FLAGS "${f90flags}")
|
||||
set_property(TARGET faddeeva PROPERTY COMPILE_FLAGS "${cflags}")
|
||||
else()
|
||||
target_compile_options(${program} PUBLIC ${f90flags})
|
||||
target_compile_options(faddeeva PRIVATE ${cflags})
|
||||
endif()
|
||||
|
||||
# Add HDF5 library directories to link line with -L
|
||||
|
|
@ -277,7 +299,7 @@ endforeach()
|
|||
|
||||
# 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(${program} ${ldflags} ${HDF5_LIBRARIES} fox_dom)
|
||||
target_link_libraries(${program} ${ldflags} ${HDF5_LIBRARIES} fox_dom faddeeva)
|
||||
|
||||
#===============================================================================
|
||||
# Install executable, scripts, manpage, license
|
||||
|
|
|
|||
|
|
@ -154,7 +154,8 @@ html_title = "OpenMC Documentation"
|
|||
# so a file named "default.css" will overwrite the builtin "default.css".
|
||||
html_static_path = ['_static']
|
||||
|
||||
html_context = {'css_files': ['_static/theme_overrides.css']}
|
||||
def setup(app):
|
||||
app.add_stylesheet('theme_overrides.css')
|
||||
|
||||
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
|
||||
# using the given strftime format.
|
||||
|
|
|
|||
|
|
@ -987,6 +987,15 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
|
|||
|
||||
*Default*: A region filling all space.
|
||||
|
||||
:temperature:
|
||||
The temperature of the cell in Kelvin. If windowed-multipole data is
|
||||
avalable, this temperature will be used to Doppler broaden some cross
|
||||
sections in the resolved resonance region. A list of temperatures can be
|
||||
specified for the "distributed temperature" feature. This will give each
|
||||
unique instance of the cell its own temperature.
|
||||
|
||||
*Default*: The temperature of the coldest nuclide in the cell's material(s)
|
||||
|
||||
:rotation:
|
||||
If the cell is filled with a universe, this element specifies the angles in
|
||||
degrees about the x, y, and z axes that the filled universe should be
|
||||
|
|
|
|||
|
|
@ -98,6 +98,10 @@ The current revision of the summary file format is 1.
|
|||
material. The data is an array if the cell uses distributed materials,
|
||||
otherwise it is a scalar.
|
||||
|
||||
**/geometry/cells/cell <uid>/temperature** (*double[]*)
|
||||
|
||||
Temperature of the cell in Kelvin.
|
||||
|
||||
**/geometry/cells/cell <uid>/offset** (*int[]*)
|
||||
|
||||
Offsets used for distribcell tally filter. This dataset is present only if
|
||||
|
|
|
|||
|
|
@ -823,7 +823,8 @@ class AggregateFilter(object):
|
|||
|
||||
"""
|
||||
|
||||
if filter_bin not in self.bins:
|
||||
if filter_bin not in self.bins and \
|
||||
filter_bin != self._aggregate_filter.bins:
|
||||
msg = 'Unable to get the bin index for AggregateFilter since ' \
|
||||
'"{0}" is not one of the bins'.format(filter_bin)
|
||||
raise ValueError(msg)
|
||||
|
|
|
|||
|
|
@ -673,9 +673,11 @@ class Filter(object):
|
|||
|
||||
# Assign entry to Lattice Multi-index column
|
||||
else:
|
||||
# Reverse y index per lattice ordering in OpenCG
|
||||
level_dict[lat_id_key][offset] = coords._lattice._id
|
||||
level_dict[lat_x_key][offset] = coords._lat_x
|
||||
level_dict[lat_y_key][offset] = coords._lat_y
|
||||
level_dict[lat_y_key][offset] = \
|
||||
coords._lattice.dimension[1] - coords._lat_y - 1
|
||||
level_dict[lat_z_key][offset] = coords._lat_z
|
||||
|
||||
# Move to next node in LocalCoords linked list
|
||||
|
|
|
|||
|
|
@ -568,8 +568,9 @@ class Library(object):
|
|||
for domain in self.domains:
|
||||
for mgxs_type in self.mgxs_types:
|
||||
mgxs = subdomain_avg_library.get_mgxs(domain, mgxs_type)
|
||||
avg_mgxs = mgxs.get_subdomain_avg_xs()
|
||||
subdomain_avg_library.all_mgxs[domain.id][mgxs_type] = avg_mgxs
|
||||
if mgxs.domain_type == 'distribcell':
|
||||
avg_mgxs = mgxs.get_subdomain_avg_xs()
|
||||
subdomain_avg_library.all_mgxs[domain.id][mgxs_type] = avg_mgxs
|
||||
|
||||
return subdomain_avg_library
|
||||
|
||||
|
|
|
|||
|
|
@ -264,6 +264,9 @@ class Summary(object):
|
|||
self._f['geometry/cells'][key]['rotation'][...]
|
||||
rotation = np.asarray(rotation, dtype=np.int)
|
||||
cell.rotation = rotation
|
||||
elif fill_type == 'normal':
|
||||
cell.temperature = \
|
||||
self._f['geometry/cells'][key]['temperature'][...]
|
||||
|
||||
# Store Cell fill information for after Universe/Lattice creation
|
||||
self._cell_fills[index] = (fill_type, fill)
|
||||
|
|
@ -367,9 +370,9 @@ class Summary(object):
|
|||
# Set the universes for the lattice
|
||||
lattice.universes = universes
|
||||
|
||||
# Set the distribcell offsets for the lattice
|
||||
if offsets is not None:
|
||||
offsets = np.swapaxes(offsets, 0, 2)
|
||||
lattice.offsets = offsets
|
||||
lattice.offsets = offsets[:, ::-1, :]
|
||||
|
||||
# Add the Lattice to the global dictionary of all Lattices
|
||||
self.lattices[index] = lattice
|
||||
|
|
|
|||
|
|
@ -2785,7 +2785,7 @@ class Tally(object):
|
|||
bin_indices.append(bin_index)
|
||||
num_bins += 1
|
||||
|
||||
find_filter.bins = set(find_filter.bins[bin_indices])
|
||||
find_filter.bins = np.unique(find_filter.bins[bin_indices])
|
||||
find_filter.num_bins = num_bins
|
||||
|
||||
# Update the new tally's filter strides
|
||||
|
|
|
|||
|
|
@ -51,9 +51,13 @@ class Cell(object):
|
|||
name : str
|
||||
Name of the cell
|
||||
fill : Material or Universe or Lattice or 'void' or iterable of Material
|
||||
Indicates what the region of space is filled with
|
||||
Indicates what the region of space is filled with. Multiple materials
|
||||
can be given to give each distributed cell instance a unique material.
|
||||
region : openmc.region.Region
|
||||
Region of space that is assigned to the cell.
|
||||
temperature : float or iterable of float
|
||||
Temperature of the cell in Kelvin. Multiple temperatures can be given
|
||||
to give each distributed cell instance a unique temperature.
|
||||
rotation : ndarray
|
||||
If the cell is filled with a universe, this array specifies the angles
|
||||
in degrees about the x, y, and z axes that the filled universe should be
|
||||
|
|
@ -75,6 +79,7 @@ class Cell(object):
|
|||
self._fill = None
|
||||
self._type = None
|
||||
self._region = None
|
||||
self._temperature = None
|
||||
self._rotation = None
|
||||
self._translation = None
|
||||
self._offsets = None
|
||||
|
|
@ -91,6 +96,8 @@ class Cell(object):
|
|||
return False
|
||||
elif self.region != other.region:
|
||||
return False
|
||||
elif self.temperature != other.temperature:
|
||||
return False
|
||||
elif self.rotation != other.rotation:
|
||||
return False
|
||||
elif self.translation != other.translation:
|
||||
|
|
@ -126,6 +133,10 @@ class Cell(object):
|
|||
|
||||
string += '{0: <16}{1}{2}\n'.format('\tRegion', '=\t', self._region)
|
||||
|
||||
if self.fill_type == 'material':
|
||||
string += '\t{0: <15}=\t{1}\n'.format('Temperature',
|
||||
self.temperature)
|
||||
|
||||
string += '{0: <16}{1}{2}\n'.format('\tRotation', '=\t',
|
||||
self._rotation)
|
||||
string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t',
|
||||
|
|
@ -150,7 +161,7 @@ class Cell(object):
|
|||
|
||||
@property
|
||||
def fill_type(self):
|
||||
if isinstance(self.fill, openmc.Material):
|
||||
if isinstance(self.fill, (openmc.Material, Iterable)):
|
||||
return 'material'
|
||||
elif isinstance(self.fill, openmc.Universe):
|
||||
return 'universe'
|
||||
|
|
@ -163,6 +174,10 @@ class Cell(object):
|
|||
def region(self):
|
||||
return self._region
|
||||
|
||||
@property
|
||||
def temperature(self):
|
||||
return self._temperature
|
||||
|
||||
@property
|
||||
def rotation(self):
|
||||
return self._rotation
|
||||
|
|
@ -251,6 +266,17 @@ class Cell(object):
|
|||
cv.check_type('cell region', region, Region)
|
||||
self._region = region
|
||||
|
||||
@temperature.setter
|
||||
def temperature(self, temperature):
|
||||
cv.check_type('cell temperature', temperature, (Iterable, Real))
|
||||
if isinstance(temperature, Iterable):
|
||||
cv.check_type('cell temperature', temperature, Iterable, Real)
|
||||
for T in temperature:
|
||||
cv.check_greater_than('cell temperature', T, 0.0, True)
|
||||
else:
|
||||
cv.check_greater_than('cell temperature', temperature, 0.0, True)
|
||||
self._temperature = temperature
|
||||
|
||||
@distribcell_index.setter
|
||||
def distribcell_index(self, ind):
|
||||
cv.check_type('distribcell index', ind, Integral)
|
||||
|
|
@ -445,6 +471,13 @@ class Cell(object):
|
|||
# Call the recursive function from the top node
|
||||
create_surface_elements(self.region, xml_element)
|
||||
|
||||
if self.temperature is not None:
|
||||
if isinstance(self.temperature, Iterable):
|
||||
element.set("temperature", ' '.join(
|
||||
[str(t) for t in self.temperature]))
|
||||
else:
|
||||
element.set("temperature", str(self.temperature))
|
||||
|
||||
if self.translation is not None:
|
||||
element.set("translation", ' '.join(map(str, self.translation)))
|
||||
|
||||
|
|
@ -1095,14 +1128,14 @@ class RectLattice(Lattice):
|
|||
|
||||
# For 2D Lattices
|
||||
if len(self._dimension) == 2:
|
||||
offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_index-1]
|
||||
offset = self._offsets[i[3]-1, i[2]-1, i[1]-1, distribcell_index-1]
|
||||
offset += self._universes[i[1]-1][i[2]-1].get_cell_instance(path,
|
||||
distribcell_index)
|
||||
|
||||
# For 3D Lattices
|
||||
else:
|
||||
offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_index-1]
|
||||
offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_cell_instance(
|
||||
offset = self._offsets[i[3]-1, i[2]-1, i[1]-1, distribcell_index-1]
|
||||
offset += self._universes[i[3]-1][i[2]-1][i[1]-1].get_cell_instance(
|
||||
path, distribcell_index)
|
||||
|
||||
return offset
|
||||
|
|
|
|||
2
setup.py
2
setup.py
|
|
@ -32,7 +32,7 @@ kwargs = {'name': 'openmc',
|
|||
if have_setuptools:
|
||||
kwargs.update({
|
||||
# Required dependencies
|
||||
'install_requires': ['numpy', 'h5py', 'matplotlib'],
|
||||
'install_requires': ['numpy>=1.9', 'h5py', 'matplotlib'],
|
||||
|
||||
# Optional dependencies
|
||||
'extras_require': {
|
||||
|
|
|
|||
3
src/Faddeeva.c
Normal file
3
src/Faddeeva.c
Normal file
|
|
@ -0,0 +1,3 @@
|
|||
/* The Faddeeva.cc file contains macros to let it compile as C code
|
||||
(assuming C99 complex-number support), so just #include it. */
|
||||
#include "Faddeeva.cc"
|
||||
2516
src/Faddeeva.cc
Normal file
2516
src/Faddeeva.cc
Normal file
File diff suppressed because it is too large
Load diff
68
src/Faddeeva.h
Normal file
68
src/Faddeeva.h
Normal file
|
|
@ -0,0 +1,68 @@
|
|||
/* Copyright (c) 2012 Massachusetts Institute of Technology
|
||||
*
|
||||
* Permission is hereby granted, free of charge, to any person obtaining
|
||||
* a copy of this software and associated documentation files (the
|
||||
* "Software"), to deal in the Software without restriction, including
|
||||
* without limitation the rights to use, copy, modify, merge, publish,
|
||||
* distribute, sublicense, and/or sell copies of the Software, and to
|
||||
* permit persons to whom the Software is furnished to do so, subject to
|
||||
* the following conditions:
|
||||
*
|
||||
* The above copyright notice and this permission notice shall be
|
||||
* included in all copies or substantial portions of the Software.
|
||||
*
|
||||
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
||||
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
||||
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
||||
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
|
||||
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
|
||||
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
|
||||
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
|
||||
*/
|
||||
|
||||
/* Available at: http://ab-initio.mit.edu/Faddeeva
|
||||
|
||||
Header file for Faddeeva.c; see Faddeeva.cc for more information. */
|
||||
|
||||
#ifndef FADDEEVA_H
|
||||
#define FADDEEVA_H 1
|
||||
|
||||
// Require C99 complex-number support
|
||||
#include <complex.h>
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C"
|
||||
{
|
||||
#endif /* __cplusplus */
|
||||
|
||||
// compute w(z) = exp(-z^2) erfc(-iz) [ Faddeeva / scaled complex error func ]
|
||||
extern double complex Faddeeva_w(double complex z,double relerr);
|
||||
extern double Faddeeva_w_im(double x); // special-case code for Im[w(x)] of real x
|
||||
|
||||
// Various functions that we can compute with the help of w(z)
|
||||
|
||||
// compute erfcx(z) = exp(z^2) erfc(z)
|
||||
extern double complex Faddeeva_erfcx(double complex z, double relerr);
|
||||
extern double Faddeeva_erfcx_re(double x); // special case for real x
|
||||
|
||||
// compute erf(z), the error function of complex arguments
|
||||
extern double complex Faddeeva_erf(double complex z, double relerr);
|
||||
extern double Faddeeva_erf_re(double x); // special case for real x
|
||||
|
||||
// compute erfi(z) = -i erf(iz), the imaginary error function
|
||||
extern double complex Faddeeva_erfi(double complex z, double relerr);
|
||||
extern double Faddeeva_erfi_re(double x); // special case for real x
|
||||
|
||||
// compute erfc(z) = 1 - erf(z), the complementary error function
|
||||
extern double complex Faddeeva_erfc(double complex z, double relerr);
|
||||
extern double Faddeeva_erfc_re(double x); // special case for real x
|
||||
|
||||
// compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z)
|
||||
extern double complex Faddeeva_Dawson(double complex z, double relerr);
|
||||
extern double Faddeeva_Dawson_re(double x); // special case for real x
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif /* __cplusplus */
|
||||
|
||||
#endif // FADDEEVA_H
|
||||
68
src/ace.F90
68
src/ace.F90
|
|
@ -11,6 +11,8 @@ module ace
|
|||
use global
|
||||
use list_header, only: ListInt
|
||||
use material_header, only: Material
|
||||
use multipole, only: multipole_read
|
||||
use multipole_header, only: max_L, max_poles, max_poly
|
||||
use output, only: write_message
|
||||
use set_header, only: SetChar
|
||||
use secondary_header, only: AngleEnergy
|
||||
|
|
@ -111,6 +113,9 @@ contains
|
|||
end do
|
||||
end if
|
||||
|
||||
! Read multipole file into the appropriate entry on the nuclides array
|
||||
call read_multipole_data(i_nuclide)
|
||||
|
||||
! Add name and alias to dictionary
|
||||
call already_read % add(name)
|
||||
call already_read % add(alias)
|
||||
|
|
@ -414,6 +419,69 @@ contains
|
|||
|
||||
end subroutine read_ace_table
|
||||
|
||||
!===============================================================================
|
||||
! READ_MULTIPOLE_DATA checks for the existence of a multipole library in the
|
||||
! directory and loads it using multipole_read
|
||||
!===============================================================================
|
||||
|
||||
subroutine read_multipole_data(i_table)
|
||||
|
||||
integer, intent(in) :: i_table ! index in nuclides/sab_tables
|
||||
|
||||
logical :: file_exists ! does multipole library exist?
|
||||
character(7) :: readable ! is multipole library readable?
|
||||
character(6) :: zaid_string ! String of the ZAID
|
||||
character(9) :: filename ! path to multipole cross section library
|
||||
type(Nuclide), pointer :: nuc => null()
|
||||
|
||||
! For the time being, and I know this is a bit hacky, we just assume
|
||||
! that the file will be zaid.h5.
|
||||
associate (nuc => nuclides(i_table))
|
||||
|
||||
write(zaid_string,'(I6.6)') nuc % zaid
|
||||
filename = zaid_string // ".h5"
|
||||
|
||||
! Check if Multipole library exists and is readable
|
||||
inquire(FILE=filename, EXIST=file_exists, READ=readable)
|
||||
if (.not. file_exists) then
|
||||
nuc % mp_present = .false.
|
||||
return
|
||||
elseif (readable(1:3) == 'NO') then
|
||||
call fatal_error("Multipole library '" // trim(filename) // "' is not readable! &
|
||||
&Change file permissions with chmod command.")
|
||||
end if
|
||||
|
||||
! display message
|
||||
call write_message("Loading Multipole XS table: " // filename, 6)
|
||||
|
||||
allocate(nuc % multipole)
|
||||
|
||||
! Call the read routine
|
||||
call multipole_read(filename, nuc % multipole, i_table)
|
||||
nuc % mp_present = .true.
|
||||
|
||||
! Update the maximum number of poles, l indices, and polynomial order
|
||||
if (nuc % multipole % max_w > max_poles) then
|
||||
max_poles = nuc % multipole % max_w
|
||||
end if
|
||||
|
||||
if (nuc % multipole % num_l > max_L) then
|
||||
max_L = nuc % multipole % num_l
|
||||
end if
|
||||
|
||||
if (nuc % multipole % fit_order + 1 > max_poly) then
|
||||
max_poly = nuc % multipole % fit_order + 1
|
||||
end if
|
||||
|
||||
! Recreate nu-fission tables
|
||||
if (nuc % fissionable) then
|
||||
call generate_nu_fission(nuc)
|
||||
end if
|
||||
|
||||
end associate
|
||||
|
||||
end subroutine read_multipole_data
|
||||
|
||||
!===============================================================================
|
||||
! READ_ESZ - reads through the ESZ block. This block contains the energy grid,
|
||||
! total xs, absorption xs, elastic scattering xs, and heating numbers.
|
||||
|
|
|
|||
|
|
@ -3,6 +3,7 @@ module ace_header
|
|||
use constants, only: MAX_FILE_LEN, ZERO
|
||||
use dict_header, only: DictIntInt
|
||||
use endf_header, only: Tab1
|
||||
use multipole_header, only: MultipoleArray
|
||||
use secondary_header, only: SecondaryDistribution, AngleEnergyContainer
|
||||
use stl_vector, only: VectorInt
|
||||
|
||||
|
|
@ -110,6 +111,10 @@ module ace_header
|
|||
integer :: urr_inelastic
|
||||
type(UrrData), pointer :: urr_data => null()
|
||||
|
||||
! Multipole data
|
||||
logical :: mp_present
|
||||
type(MultipoleArray), allocatable :: multipole
|
||||
|
||||
! Reactions
|
||||
integer :: n_reaction ! # of reactions
|
||||
type(Reaction), allocatable :: reactions(:)
|
||||
|
|
@ -233,6 +238,9 @@ module ace_header
|
|||
! Information for URR probability table use
|
||||
logical :: use_ptable ! in URR range with probability tables?
|
||||
real(8) :: last_prn
|
||||
|
||||
! Information for Doppler broadening
|
||||
real(8) :: last_sqrtkT = ZERO ! last temperature in sqrt(Boltzmann constant * temperature (MeV))
|
||||
end type NuclideMicroXS
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -271,6 +279,10 @@ module ace_header
|
|||
|
||||
if (associated(this % urr_data)) deallocate(this % urr_data)
|
||||
|
||||
if (this % mp_present) then
|
||||
deallocate(this % multipole)
|
||||
end if
|
||||
|
||||
if (allocated(this % reactions)) then
|
||||
do i = 1, size(this % reactions)
|
||||
call this % reactions(i) % clear()
|
||||
|
|
|
|||
|
|
@ -63,6 +63,7 @@ module constants
|
|||
|
||||
real(8), parameter :: &
|
||||
PI = 3.1415926535898_8, & ! pi
|
||||
SQRT_PI = 1.7724538509055_8, & ! square root of pi
|
||||
MASS_NEUTRON = 1.008664916_8, & ! mass of a neutron in amu
|
||||
MASS_NEUTRON_MEV = 939.565379_8, & ! mass of a neutron in MeV/c^2
|
||||
MASS_PROTON = 1.007276466812_8, & ! mass of a proton in amu
|
||||
|
|
@ -77,6 +78,7 @@ module constants
|
|||
TWO = 2.0_8, &
|
||||
THREE = 3.0_8, &
|
||||
FOUR = 4.0_8
|
||||
complex(8), parameter :: ONEI = (ZERO, ONE)
|
||||
|
||||
! ============================================================================
|
||||
! GEOMETRY-RELATED CONSTANTS
|
||||
|
|
|
|||
|
|
@ -1,19 +1,31 @@
|
|||
module cross_section
|
||||
|
||||
use ace_header, only: Nuclide, SAlphaBeta, Reaction, UrrData
|
||||
use ace_header, only: Nuclide, SAlphaBeta, Reaction, UrrData
|
||||
use constants
|
||||
use energy_grid, only: grid_method, log_spacing
|
||||
use error, only: fatal_error
|
||||
use fission, only: nu_total
|
||||
use energy_grid, only: grid_method, log_spacing
|
||||
use error, only: fatal_error
|
||||
use fission, only: nu_total
|
||||
use global
|
||||
use list_header, only: ListElemInt
|
||||
use material_header, only: Material
|
||||
use particle_header, only: Particle
|
||||
use random_lcg, only: prn
|
||||
use search, only: binary_search
|
||||
use list_header, only: ListElemInt
|
||||
use material_header, only: Material
|
||||
use math, only: w, broaden_n_polynomials
|
||||
use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, RM_RF, &
|
||||
MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, FIT_T, FIT_A, FIT_F, &
|
||||
MultipoleArray, max_poly, max_L, max_poles
|
||||
use particle_header, only: Particle
|
||||
use random_lcg, only: prn
|
||||
use search, only: binary_search
|
||||
|
||||
implicit none
|
||||
|
||||
! Allocatable arrays for multipole that are allocated once for speed purposes
|
||||
complex(8), allocatable :: sigT_factor(:)
|
||||
real(8), allocatable :: twophi(:)
|
||||
real(8), allocatable :: broadened_polynomials(:)
|
||||
logical :: mp_already_alloc = .false.
|
||||
|
||||
!$omp threadprivate(sigT_factor, twophi, broadened_polynomials, mp_already_alloc)
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -93,10 +105,13 @@ contains
|
|||
i_nuclide = mat % nuclide(i)
|
||||
|
||||
! 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, i_grid)
|
||||
if (p % E /= micro_xs(i_nuclide) % last_E &
|
||||
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
|
||||
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, &
|
||||
i_grid, p % sqrtkT)
|
||||
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, i_grid)
|
||||
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, &
|
||||
i_grid, p % sqrtkT)
|
||||
end if
|
||||
|
||||
! ========================================================================
|
||||
|
|
@ -133,7 +148,8 @@ 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, i_log_union)
|
||||
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, &
|
||||
i_log_union, sqrtkT)
|
||||
integer, intent(in) :: i_nuclide ! index into nuclides array
|
||||
integer, intent(in) :: i_sab ! index into sab_tables array
|
||||
real(8), intent(in) :: E ! energy
|
||||
|
|
@ -141,11 +157,13 @@ contains
|
|||
integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material
|
||||
integer, intent(in) :: i_log_union ! index into logarithmic mapping array or
|
||||
! material union energy grid
|
||||
real(8), intent(in) :: sqrtkT ! Square root of kT, material dependent
|
||||
|
||||
integer :: i_grid ! index on nuclide energy grid
|
||||
integer :: i_low ! lower logarithmic mapping index
|
||||
integer :: i_high ! upper logarithmic mapping index
|
||||
real(8) :: f ! interp factor on nuclide energy grid
|
||||
real(8) :: sigT, sigA, sigF ! Intermediate multipole variables
|
||||
type(Nuclide), pointer :: nuc
|
||||
type(Material), pointer :: mat
|
||||
|
||||
|
|
@ -153,53 +171,111 @@ contains
|
|||
nuc => nuclides(i_nuclide)
|
||||
mat => materials(i_mat)
|
||||
|
||||
! Determine index on nuclide energy grid
|
||||
select case (grid_method)
|
||||
case (GRID_MAT_UNION)
|
||||
! If MP, don't interpolate, it's all already baked in.
|
||||
if (nuc % mp_present .and. &
|
||||
(E >= nuc % multipole % start_E/1.0e6_8 .and.&
|
||||
E <= nuc % multipole % end_E/1.0e6_8)) then
|
||||
|
||||
i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union)
|
||||
! Call multipole kernel
|
||||
call multipole_eval(nuc % multipole, E, sqrtkT, sigT, sigA, sigF)
|
||||
|
||||
case (GRID_LOGARITHM)
|
||||
! Determine the energy grid index using a logarithmic mapping to reduce
|
||||
! the energy range over which a binary search needs to be performed
|
||||
micro_xs(i_nuclide) % total = sigT
|
||||
micro_xs(i_nuclide) % absorption = sigA
|
||||
micro_xs(i_nuclide) % elastic = sigT - sigA
|
||||
|
||||
if (E < nuc % energy(1)) then
|
||||
i_grid = 1
|
||||
elseif (E > nuc % energy(nuc % n_grid)) then
|
||||
i_grid = nuc % n_grid - 1
|
||||
if (nuc % fissionable) then
|
||||
micro_xs(i_nuclide) % fission = sigF
|
||||
micro_xs(i_nuclide) % nu_fission = sigF * nu_total(nuc, E)
|
||||
else
|
||||
! Determine bounding indices based on which equal log-spaced interval
|
||||
! the energy is in
|
||||
i_low = nuc % grid_index(i_log_union)
|
||||
i_high = nuc % grid_index(i_log_union + 1) + 1
|
||||
|
||||
! Perform binary search over reduced range
|
||||
i_grid = binary_search(nuc % energy(i_low:i_high), &
|
||||
i_high - i_low + 1, E) + i_low - 1
|
||||
micro_xs(i_nuclide) % fission = ZERO
|
||||
micro_xs(i_nuclide) % nu_fission = ZERO
|
||||
end if
|
||||
|
||||
case (GRID_NUCLIDE)
|
||||
! Perform binary search on the nuclide energy grid in order to determine
|
||||
! which points to interpolate between
|
||||
! Ensure these values are set
|
||||
! Note, the only time either is used is in one of 4 places:
|
||||
! 1. physics.F90 - scatter - For inelastic scatter.
|
||||
! 2. physics.F90 - sample_fission - For partial fissions.
|
||||
! 3. tally.F90 - score_general - For tallying on MTxxx reactions.
|
||||
! 4. cross_section.F90 - calculate_urr_xs - For unresolved purposes.
|
||||
! It is worth noting that none of these occur in the resolved
|
||||
! resonance range, so the value here does not matter.
|
||||
micro_xs(i_nuclide) % index_grid = 0
|
||||
micro_xs(i_nuclide) % interp_factor = ZERO
|
||||
else
|
||||
! Determine index on nuclide energy grid
|
||||
select case (grid_method)
|
||||
case (GRID_MAT_UNION)
|
||||
|
||||
if (E <= nuc % energy(1)) then
|
||||
i_grid = 1
|
||||
elseif (E > nuc % energy(nuc % n_grid)) then
|
||||
i_grid = nuc % n_grid - 1
|
||||
else
|
||||
i_grid = binary_search(nuc % energy, nuc % n_grid, E)
|
||||
i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union)
|
||||
|
||||
case (GRID_LOGARITHM)
|
||||
! Determine the energy grid index using a logarithmic mapping to reduce
|
||||
! the energy range over which a binary search needs to be performed
|
||||
|
||||
if (E < nuc % energy(1)) then
|
||||
i_grid = 1
|
||||
elseif (E > nuc % energy(nuc % n_grid)) then
|
||||
i_grid = nuc % n_grid - 1
|
||||
else
|
||||
! Determine bounding indices based on which equal log-spaced interval
|
||||
! the energy is in
|
||||
i_low = nuc % grid_index(i_log_union)
|
||||
i_high = nuc % grid_index(i_log_union + 1) + 1
|
||||
|
||||
! Perform binary search over reduced range
|
||||
i_grid = binary_search(nuc % energy(i_low:i_high), &
|
||||
i_high - i_low + 1, E) + i_low - 1
|
||||
end if
|
||||
|
||||
case (GRID_NUCLIDE)
|
||||
! Perform binary search on the nuclide energy grid in order to determine
|
||||
! which points to interpolate between
|
||||
|
||||
if (E <= nuc % energy(1)) then
|
||||
i_grid = 1
|
||||
elseif (E > nuc % energy(nuc % n_grid)) then
|
||||
i_grid = nuc % n_grid - 1
|
||||
else
|
||||
i_grid = binary_search(nuc % energy, nuc % n_grid, E)
|
||||
end if
|
||||
|
||||
end select
|
||||
|
||||
! check for rare case where two energy points are the same
|
||||
if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1
|
||||
|
||||
! calculate interpolation factor
|
||||
f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid))
|
||||
|
||||
micro_xs(i_nuclide) % index_grid = i_grid
|
||||
micro_xs(i_nuclide) % interp_factor = f
|
||||
|
||||
! Initialize nuclide cross-sections to zero
|
||||
micro_xs(i_nuclide) % fission = ZERO
|
||||
micro_xs(i_nuclide) % nu_fission = ZERO
|
||||
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) &
|
||||
+ f * nuc % total(i_grid+1)
|
||||
|
||||
! Calculate microscopic nuclide elastic cross section
|
||||
micro_xs(i_nuclide) % elastic = (ONE - f) * nuc % elastic(i_grid) &
|
||||
+ f * nuc % elastic(i_grid+1)
|
||||
|
||||
! Calculate microscopic nuclide absorption cross section
|
||||
micro_xs(i_nuclide) % absorption = (ONE - f) * nuc % absorption( &
|
||||
i_grid) + f * nuc % absorption(i_grid+1)
|
||||
|
||||
if (nuc % fissionable) then
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(i_nuclide) % fission = (ONE - f) * nuc % fission(i_grid) &
|
||||
+ f * nuc % fission(i_grid+1)
|
||||
|
||||
! Calculate microscopic nuclide nu-fission cross section
|
||||
micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( &
|
||||
i_grid) + f * nuc % nu_fission(i_grid+1)
|
||||
end if
|
||||
|
||||
end select
|
||||
|
||||
! check for rare case where two energy points are the same
|
||||
if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1
|
||||
|
||||
! calculate interpolation factor
|
||||
f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid))
|
||||
|
||||
micro_xs(i_nuclide) % index_grid = i_grid
|
||||
micro_xs(i_nuclide) % interp_factor = f
|
||||
end if
|
||||
|
||||
! Initialize sab treatment to false
|
||||
micro_xs(i_nuclide) % index_sab = NONE
|
||||
|
|
@ -208,32 +284,6 @@ contains
|
|||
! Initialize URR probability table treatment to false
|
||||
micro_xs(i_nuclide) % use_ptable = .false.
|
||||
|
||||
! Initialize nuclide cross-sections to zero
|
||||
micro_xs(i_nuclide) % fission = ZERO
|
||||
micro_xs(i_nuclide) % nu_fission = ZERO
|
||||
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) &
|
||||
+ f * nuc % total(i_grid+1)
|
||||
|
||||
! Calculate microscopic nuclide elastic cross section
|
||||
micro_xs(i_nuclide) % elastic = (ONE - f) * nuc % elastic(i_grid) &
|
||||
+ f * nuc % elastic(i_grid+1)
|
||||
|
||||
! Calculate microscopic nuclide absorption cross section
|
||||
micro_xs(i_nuclide) % absorption = (ONE - f) * nuc % absorption( &
|
||||
i_grid) + f * nuc % absorption(i_grid+1)
|
||||
|
||||
if (nuc % fissionable) then
|
||||
! Calculate microscopic nuclide total cross section
|
||||
micro_xs(i_nuclide) % fission = (ONE - f) * nuc % fission(i_grid) &
|
||||
+ f * nuc % fission(i_grid+1)
|
||||
|
||||
! Calculate microscopic nuclide nu-fission cross section
|
||||
micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( &
|
||||
i_grid) + f * nuc % nu_fission(i_grid+1)
|
||||
end if
|
||||
|
||||
! If there is S(a,b) data for this nuclide, we need to do a few
|
||||
! things. Since the total cross section was based on non-S(a,b) data, we
|
||||
! need to correct it by subtracting the non-S(a,b) elastic cross section and
|
||||
|
|
@ -253,6 +303,7 @@ contains
|
|||
|
||||
micro_xs(i_nuclide) % last_E = E
|
||||
micro_xs(i_nuclide) % last_index_sab = i_sab
|
||||
micro_xs(i_nuclide) % last_sqrtkT = sqrtkT
|
||||
|
||||
end subroutine calculate_nuclide_xs
|
||||
|
||||
|
|
@ -525,6 +576,192 @@ contains
|
|||
|
||||
end function find_energy_index
|
||||
|
||||
!===============================================================================
|
||||
! MULTIPOLE_EVAL_ALLOCATE allocates fixed-length arrays that vary based on
|
||||
! what nuclides are loaded into the problem
|
||||
!===============================================================================
|
||||
|
||||
subroutine multipole_eval_allocate()
|
||||
allocate(sigT_factor(max_L))
|
||||
allocate(twophi(max_L))
|
||||
allocate(broadened_polynomials(max_poly))
|
||||
|
||||
mp_already_alloc = .true.
|
||||
end subroutine
|
||||
|
||||
!===============================================================================
|
||||
! MULTIPOLE_EVAL evaluates the windowed multipole equations for cross
|
||||
! sections in the resolved resonance regions
|
||||
!===============================================================================
|
||||
|
||||
subroutine multipole_eval(multipole, Emev, sqrtkT_, sigT, sigA, sigF)
|
||||
type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
|
||||
! object to process.
|
||||
real(8), intent(in) :: Emev ! The energy at which to
|
||||
! evaluate the cross section
|
||||
! in MeV
|
||||
real(8), intent(in) :: sqrtkT_ ! The temperature in the form
|
||||
! sqrt(kT (in MeV)), at which
|
||||
! to evaluate the XS.
|
||||
real(8), intent(out) :: sigT ! Total cross section
|
||||
real(8), intent(out) :: sigA ! Absorption cross section
|
||||
real(8), intent(out) :: sigF ! Fission cross section
|
||||
complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic
|
||||
! form
|
||||
complex(8) :: c_temp ! complex temporary variable
|
||||
complex(8) :: w_val ! The faddeeva function evaluated at Z
|
||||
complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
|
||||
real(8) :: sqrtE ! sqrt(E), eV
|
||||
real(8) :: invE ! 1/E, eV
|
||||
real(8) :: dopp ! sqrt(atomic weight ratio / kT)
|
||||
real(8) :: dopp_ecoef ! sqrt(atomic weight ratio * pi / kT) / E
|
||||
real(8) :: temp ! real temporary value
|
||||
real(8) :: E ! energy, eV
|
||||
real(8) :: sqrtkT ! sqrt(kT (in eV))
|
||||
integer :: iP ! index of pole
|
||||
integer :: iC ! index of curvefit
|
||||
integer :: iW ! index of window
|
||||
integer :: startw ! window start pointer (for poles)
|
||||
integer :: startw_1 ! window start pointer - 1
|
||||
integer :: startw_endw ! window start pointer - window end pointer
|
||||
integer :: endw ! window end pointer
|
||||
|
||||
! Convert to eV
|
||||
E = Emev * 1.0e6_8
|
||||
sqrtkT = sqrtkT_ * 1.0e3_8
|
||||
|
||||
sqrtE = sqrt(E)
|
||||
invE = ONE/E
|
||||
|
||||
if(.not. mp_already_alloc) then
|
||||
call multipole_eval_allocate()
|
||||
end if
|
||||
|
||||
! Locate us
|
||||
iW = floor((sqrtE - sqrt(multipole % start_E))/multipole % spacing + ONE)
|
||||
|
||||
startw = multipole % w_start(iW)
|
||||
startw_1 = startw - 1 ! This is an index shift parameter.
|
||||
endw = multipole % w_end(iW)
|
||||
startw_endw = endw - startw + 1
|
||||
|
||||
! Fill in factors
|
||||
if (startw <= endw) then
|
||||
call fill_factors(multipole, sqrtE, sigT_factor, twophi, multipole % num_l)
|
||||
end if
|
||||
|
||||
! Generate some doppler broadening parameters
|
||||
|
||||
! dopp_ecoef is inverse of dopp, divided by E, multiplied by sqrt(pi).
|
||||
dopp = multipole % sqrtAWR / sqrtKT
|
||||
dopp_ecoef = dopp * invE * SQRT_PI
|
||||
|
||||
sigT = ZERO
|
||||
sigA = ZERO
|
||||
sigF = ZERO
|
||||
! Evaluate linefit first
|
||||
|
||||
if(sqrtkT /= 0 .and. multipole % broaden_poly(iW) == 1) then ! Broaden the curvefit.
|
||||
call broaden_n_polynomials(E, dopp, multipole % fit_order + 1, broadened_polynomials)
|
||||
|
||||
do iC = 1, multipole % fit_order+1
|
||||
sigT = sigT + multipole % curvefit(FIT_T, iC, iW)*broadened_polynomials(iC)
|
||||
sigA = sigA + multipole % curvefit(FIT_A, iC, iW)*broadened_polynomials(iC)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + multipole % curvefit(FIT_F, iC, iW)*broadened_polynomials(iC)
|
||||
end if
|
||||
end do
|
||||
else ! Evaluate as if it were a polynomial
|
||||
temp = invE
|
||||
do iC = 1, multipole % fit_order+1
|
||||
|
||||
sigT = sigT + multipole % curvefit(FIT_T, iC, iW)*temp
|
||||
sigA = sigA + multipole % curvefit(FIT_A, iC, iW)*temp
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + multipole % curvefit(FIT_F, iC, iW)*temp
|
||||
end if
|
||||
|
||||
temp = temp * sqrtE
|
||||
end do
|
||||
end if
|
||||
|
||||
! Then get the poles we want and broaden them.
|
||||
|
||||
if (sqrtkT == ZERO) then
|
||||
! If at 0K, use asymptotic form.
|
||||
do iP = startw, endw
|
||||
psi_ki = -ONEI/(multipole % data(MP_EA, iP) - sqrtE)
|
||||
c_temp = psi_ki/E
|
||||
if (multipole % formalism == FORM_MLBW) then
|
||||
sigT = sigT + real(multipole % data(MLBW_RT, iP) * c_temp * &
|
||||
sigT_factor(multipole % l_value(iP))) &
|
||||
+ real(multipole % data(MLBW_RX, iP) * c_temp)
|
||||
sigA = sigA + real(multipole % data(MLBW_RA, iP) * c_temp)
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, iP) * c_temp)
|
||||
else if (multipole % formalism == FORM_RM) then
|
||||
sigT = sigT + real(multipole % data(RM_RT, iP) * c_temp* &
|
||||
sigT_factor(multipole % l_value(iP)))
|
||||
sigA = sigA + real(multipole % data(RM_RA, iP) * c_temp)
|
||||
sigF = sigF + real(multipole % data(RM_RF, iP) * c_temp)
|
||||
end if
|
||||
end do
|
||||
else
|
||||
! At temperature, use Faddeeva function-based form.
|
||||
if(endw >= startw) then
|
||||
do iP = startw, endw
|
||||
Z = (sqrtE - multipole % data(MP_EA, iP)) * dopp
|
||||
w_val = w(Z) * dopp_ecoef
|
||||
|
||||
if (multipole % formalism == FORM_MLBW) then
|
||||
sigT = sigT + real((multipole % data(MLBW_RT, iP) * &
|
||||
sigT_factor(multipole%l_value(iP)) + &
|
||||
multipole % data(MLBW_RX, iP)) * w_val)
|
||||
sigA = sigA + real(multipole % data(MLBW_RA, iP) * w_val)
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, iP) * w_val)
|
||||
else if (multipole % formalism == FORM_RM) then
|
||||
sigT = sigT + real(multipole % data(RM_RT, iP) * w_val * &
|
||||
sigT_factor(multipole % l_value(iP)))
|
||||
sigA = sigA + real(multipole % data(RM_RA, iP) * w_val)
|
||||
sigF = sigF + real(multipole % data(RM_RF, iP) * w_val)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
end subroutine
|
||||
|
||||
!===============================================================================
|
||||
! FILL_FACTORS calculates the value of phi, the hardsphere phase shift factor,
|
||||
! and sigT_factor, a factor inside of the sigT equation not present in the
|
||||
! sigA and sigF equations.
|
||||
!===============================================================================
|
||||
|
||||
subroutine fill_factors(multipole, sqrtE, sigT_factor, twophi, max_L)
|
||||
type(MultipoleArray), intent(in) :: multipole
|
||||
real(8), intent(in) :: sqrtE
|
||||
integer, intent(in) :: max_L
|
||||
complex(8), intent(out) :: sigT_factor(max_L)
|
||||
real(8), intent(out) :: twophi(max_L)
|
||||
|
||||
integer :: iL
|
||||
real(8) :: arg
|
||||
|
||||
do iL = 1, max_L
|
||||
twophi(iL) = multipole % pseudo_k0RS(iL) * sqrtE
|
||||
if (iL == 2) then
|
||||
twophi(iL) = twophi(iL) - atan(twophi(iL))
|
||||
else if (iL == 3) then
|
||||
arg = 3.0_8 * twophi(iL) / (3.0_8 - twophi(iL)**2)
|
||||
twophi(iL) = twophi(iL) - atan(arg)
|
||||
else if (iL == 4) then
|
||||
arg = twophi(iL) * (15.0_8 - twophi(iL)**2) / (15.0_8 - 6.0_8 * twophi(iL)**2)
|
||||
twophi(iL) = twophi(iL) - atan(arg)
|
||||
end if
|
||||
end do
|
||||
|
||||
twophi = 2.0_8 * twophi
|
||||
sigT_factor = cmplx(cos(twophi),-sin(twophi), KIND=8)
|
||||
end subroutine
|
||||
|
||||
!===============================================================================
|
||||
! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
|
||||
! for a given nuclide at the trial relative energy used in resonance scattering
|
||||
|
|
|
|||
|
|
@ -279,6 +279,36 @@ contains
|
|||
p % material = c % material(offset + 1)
|
||||
end if
|
||||
|
||||
! Set the particle temperature
|
||||
if (size(c % sqrtkT) == 1) then
|
||||
! Only one temperature for this cell; assign that one to the particle.
|
||||
p % sqrtkT = c % sqrtkT(1)
|
||||
else
|
||||
! Distributed instances of this cell have different temperatures.
|
||||
! Determine which instance this is and assign the matching temp.
|
||||
distribcell_index = c % distribcell_index
|
||||
offset = 0
|
||||
do k = 1, p % n_coord
|
||||
if (cells(p % coord(k) % cell) % type == CELL_FILL) then
|
||||
offset = offset + cells(p % coord(k) % cell) % &
|
||||
offset(distribcell_index)
|
||||
elseif (cells(p % coord(k) % cell) % type == CELL_LATTICE) then
|
||||
if (lattices(p % coord(k + 1) % lattice) % obj &
|
||||
% are_valid_indices([&
|
||||
p % coord(k + 1) % lattice_x, &
|
||||
p % coord(k + 1) % lattice_y, &
|
||||
p % coord(k + 1) % lattice_z])) then
|
||||
offset = offset + lattices(p % coord(k + 1) % lattice) % obj % &
|
||||
offset(distribcell_index, &
|
||||
p % coord(k + 1) % lattice_x, &
|
||||
p % coord(k + 1) % lattice_y, &
|
||||
p % coord(k + 1) % lattice_z)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
p % sqrtkT = c % sqrtkT(offset + 1)
|
||||
end if
|
||||
|
||||
elseif (c % type == CELL_FILL) then CELL_TYPE
|
||||
! ======================================================================
|
||||
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
|
||||
|
|
|
|||
|
|
@ -139,6 +139,9 @@ module geometry_header
|
|||
! only)
|
||||
integer :: distribcell_index ! Index corresponding to this cell in
|
||||
! distribcell arrays
|
||||
real(8), allocatable :: sqrtkT(:) ! Square root of k_Boltzmann *
|
||||
! temperature in MeV. Multiple for
|
||||
! distribcell
|
||||
|
||||
! Rotation matrix and translation vector
|
||||
real(8), allocatable :: translation(:)
|
||||
|
|
|
|||
|
|
@ -62,6 +62,7 @@ module hdf5_interface
|
|||
module procedure read_string_1D
|
||||
module procedure read_tally_result_1D
|
||||
module procedure read_tally_result_2D
|
||||
module procedure read_complex_2D
|
||||
end interface read_dataset
|
||||
|
||||
public :: write_dataset
|
||||
|
|
@ -1921,4 +1922,93 @@ contains
|
|||
mpio = (driver == H5FD_MPIO_F)
|
||||
end function using_mpio_device
|
||||
|
||||
!===============================================================================
|
||||
! READ_COMPLEX_2D reads double precision complex 2-D array data as output by
|
||||
! the h5py HDF5 python module.
|
||||
!===============================================================================
|
||||
|
||||
subroutine read_complex_2D(group_id, name, buffer, indep)
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
character(*), intent(in) :: name ! name of data
|
||||
complex(8), intent(inout), target :: buffer(:,:) ! data to write
|
||||
logical, intent(in), optional :: indep ! independent I/O
|
||||
|
||||
integer(HSIZE_T) :: dims(2)
|
||||
|
||||
dims(:) = shape(buffer)
|
||||
if (present(indep)) then
|
||||
call read_complex_2D_explicit(group_id, dims, name, buffer, indep)
|
||||
else
|
||||
call read_complex_2D_explicit(group_id, dims, name, buffer)
|
||||
end if
|
||||
end subroutine read_complex_2D
|
||||
|
||||
subroutine read_complex_2D_explicit(group_id, dims, name, buffer, indep)
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
integer(HSIZE_T), intent(in) :: dims(2)
|
||||
character(*), intent(in) :: name ! name of data
|
||||
complex(8), intent(inout), target :: buffer(dims(1),dims(2))
|
||||
logical, intent(in), optional :: indep ! independent I/O
|
||||
|
||||
real(8), target :: buffer_r(dims(1), dims(2))
|
||||
real(8), target :: buffer_i(dims(1), dims(2))
|
||||
|
||||
integer(HSIZE_T) :: i, j
|
||||
|
||||
integer :: hdf5_err
|
||||
integer :: data_xfer_mode
|
||||
#ifdef PHDF5
|
||||
integer(HID_T) :: plist ! property list
|
||||
#endif
|
||||
integer(HID_T) :: dset ! data set handle
|
||||
type(c_ptr) :: f_ptr_r, f_ptr_i
|
||||
|
||||
! Components needed for complex type support
|
||||
integer(HID_T) :: dtype_real
|
||||
integer(HID_T) :: dtype_imag
|
||||
integer(SIZE_T) :: size_double
|
||||
integer :: error
|
||||
|
||||
! Create the complex type
|
||||
call h5tget_size_f(H5T_NATIVE_DOUBLE, size_double, error)
|
||||
|
||||
! Insert the 'r' and 'i' identifiers
|
||||
call h5tcreate_f(H5T_COMPOUND_F, size_double, dtype_real, error)
|
||||
call h5tcreate_f(H5T_COMPOUND_F, size_double, dtype_imag, error)
|
||||
call h5tinsert_f(dtype_real, "r", 0_8, H5T_NATIVE_DOUBLE, error)
|
||||
call h5tinsert_f(dtype_imag, "i", 0_8, H5T_NATIVE_DOUBLE, error)
|
||||
|
||||
! Set up collective vs. independent I/O
|
||||
data_xfer_mode = H5FD_MPIO_COLLECTIVE_F
|
||||
if (present(indep)) then
|
||||
if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F
|
||||
end if
|
||||
|
||||
call h5dopen_f(group_id, trim(name), dset, hdf5_err)
|
||||
f_ptr_r = c_loc(buffer_r)
|
||||
f_ptr_i = c_loc(buffer_i)
|
||||
|
||||
if (using_mpio_device(group_id)) then
|
||||
#ifdef PHDF5
|
||||
call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err)
|
||||
call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err)
|
||||
call h5dread_f(dset, dtype_real, f_ptr_r, hdf5_err, xfer_prp=plist)
|
||||
call h5dread_f(dset, dtype_imag, f_ptr_i, hdf5_err, xfer_prp=plist)
|
||||
call h5pclose_f(plist, hdf5_err)
|
||||
#endif
|
||||
else
|
||||
call h5dread_f(dset, dtype_real, f_ptr_r, hdf5_err)
|
||||
call h5dread_f(dset, dtype_imag, f_ptr_i, hdf5_err)
|
||||
end if
|
||||
|
||||
! Reconstitute the complex numbers
|
||||
do i = 1,dims(1)
|
||||
do j = 1,dims(2)
|
||||
buffer(i,j) = cmplx(buffer_r(i,j), buffer_i(i,j), kind=8)
|
||||
end do
|
||||
end do
|
||||
|
||||
call h5dclose_f(dset, hdf5_err)
|
||||
end subroutine read_complex_2D_explicit
|
||||
|
||||
end module hdf5_interface
|
||||
|
|
|
|||
|
|
@ -120,6 +120,9 @@ contains
|
|||
! Create linked lists for multiple instances of the same nuclide
|
||||
call same_nuclide_list()
|
||||
|
||||
! Set undefined cell temperatures to match the material data.
|
||||
call lookup_material_temperatures()
|
||||
|
||||
! Construct unionized or log energy grid for cross-sections
|
||||
select case (grid_method)
|
||||
case (GRID_NUCLIDE)
|
||||
|
|
@ -976,10 +979,11 @@ contains
|
|||
end do
|
||||
end do
|
||||
|
||||
! We also need distribcell if any distributed materials are present.
|
||||
! We also need distribcell if any distributed materials or distributed
|
||||
! temperatues are present.
|
||||
if (.not. distribcell_active) then
|
||||
do i = 1, n_cells
|
||||
if (size(cells(i) % material) > 1) then
|
||||
if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then
|
||||
distribcell_active = .true.
|
||||
exit
|
||||
end if
|
||||
|
|
@ -1004,8 +1008,8 @@ contains
|
|||
end do
|
||||
end do
|
||||
|
||||
! Make sure the number of materials matches the number of cell instances for
|
||||
! distributed materials.
|
||||
! Make sure the number of materials and temperatures matches the number of
|
||||
! cell instances.
|
||||
do i = 1, n_cells
|
||||
associate (c => cells(i))
|
||||
if (size(c % material) > 1) then
|
||||
|
|
@ -1017,6 +1021,15 @@ contains
|
|||
&equal one or the number of instances.")
|
||||
end if
|
||||
end if
|
||||
if (size(c % sqrtkT) > 1) then
|
||||
if (size(c % sqrtkT) /= c % instances) then
|
||||
call fatal_error("Cell " // trim(to_str(c % id)) // " was &
|
||||
&specified with " // trim(to_str(size(c % sqrtkT))) &
|
||||
// " temperatures but has " // trim(to_str(c % instances)) &
|
||||
// " distributed instances. The number of temperatures must &
|
||||
&equal one or the number of instances.")
|
||||
end if
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
||||
|
|
@ -1058,9 +1071,9 @@ contains
|
|||
end do
|
||||
end do
|
||||
|
||||
! List all cells with multiple (distributed) materials.
|
||||
! List all cells with multiple (distributed) materials or temperatures.
|
||||
do i = 1, n_cells
|
||||
if (size(cells(i) % material) > 1) then
|
||||
if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then
|
||||
call cell_list % add(i)
|
||||
end if
|
||||
end do
|
||||
|
|
@ -1129,4 +1142,57 @@ contains
|
|||
|
||||
end subroutine allocate_offsets
|
||||
|
||||
!===============================================================================
|
||||
! LOOKUP_MATERIAL_TEMPERATURES If any cells have undefined temperatures, try to
|
||||
! find their temperatures from material data.
|
||||
!===============================================================================
|
||||
|
||||
subroutine lookup_material_temperatures()
|
||||
integer :: i, j, k
|
||||
real(8) :: min_temp
|
||||
logical :: warning_given
|
||||
|
||||
warning_given = .false.
|
||||
do i = 1, n_cells
|
||||
! Ignore non-normal cells and cells with defined temperature.
|
||||
if (cells(i) % type /= CELL_NORMAL) cycle
|
||||
if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle
|
||||
|
||||
! Set the number of temperatures equal to the number of materials.
|
||||
deallocate(cells(i) % sqrtkT)
|
||||
allocate(cells(i) % sqrtkT(size(cells(i) % material)))
|
||||
|
||||
! Check each of the cell materials for temperature data.
|
||||
do j = 1, size(cells(i) % material)
|
||||
! Arbitrarily set void regions to 0K.
|
||||
if (cells(i) % material(j) == MATERIAL_VOID) then
|
||||
cells(i) % sqrtkT(j) = ZERO
|
||||
cycle
|
||||
end if
|
||||
|
||||
associate (mat => materials(cells(i) % material(j)))
|
||||
! Find the temperature of the coldest nuclide.
|
||||
min_temp = nuclides(mat % nuclide(1)) % kT
|
||||
do k = 2, mat % n_nuclides
|
||||
! Warn the user if the nuclides don't have identical temperatues.
|
||||
if (nuclides(mat % nuclide(k)) % kT /= min_temp &
|
||||
.and. .not. warning_given) then
|
||||
call warning("OpenMC cannot &
|
||||
&identify the temperature of at least one cell. For the &
|
||||
&purposes of multipole cross section evaluations, all cells &
|
||||
&with unknown temperature will be set to the coldest &
|
||||
&temperature found in the nuclear data for that cell's &
|
||||
&material")
|
||||
warning_given = .true.
|
||||
end if
|
||||
min_temp = min(min_temp, nuclides(mat % nuclide(k)) % kT)
|
||||
end do
|
||||
|
||||
! Set the temperature for this cell instance.
|
||||
cells(i) % sqrtkT(j) = sqrt(min_temp)
|
||||
end associate
|
||||
end do
|
||||
end do
|
||||
end subroutine lookup_material_temperatures
|
||||
|
||||
end module initialize
|
||||
|
|
|
|||
|
|
@ -1298,6 +1298,40 @@ contains
|
|||
call get_node_array(node_cell, "translation", c % translation)
|
||||
end if
|
||||
|
||||
! Read cell temperatures. If the temperature is not specified, set it to
|
||||
! ERROR_REAL for now. During initialization we'll replace ERROR_REAL with
|
||||
! the temperature from the material data.
|
||||
if (check_for_node(node_cell, "temperature")) then
|
||||
n = get_arraysize_double(node_cell, "temperature")
|
||||
if (n > 0) then
|
||||
! Make sure this is a "normal" cell.
|
||||
if (c % material(1) == NONE) call fatal_error("Cell " &
|
||||
// trim(to_str(c % id)) // " was specified with a temperature &
|
||||
&but no material. Temperature specification is only valid for &
|
||||
&cells filled with a material.")
|
||||
|
||||
! Copy in temperatures
|
||||
allocate(c % sqrtkT(n))
|
||||
call get_node_array(node_cell, "temperature", c % sqrtkT)
|
||||
|
||||
! Make sure all temperatues are positive
|
||||
do j = 1, size(c % sqrtkT)
|
||||
if (c % sqrtkT(j) < ZERO) call fatal_error("Cell " &
|
||||
// trim(to_str(c % id)) // " was specified with a negative &
|
||||
&temperature. All cell temperatures must be non-negative.")
|
||||
end do
|
||||
|
||||
! Convert to sqrt(kT)
|
||||
c % sqrtkT(:) = sqrt(K_BOLTZMANN * c % sqrtkT(:))
|
||||
else
|
||||
allocate(c % sqrtkT(1))
|
||||
c % sqrtkT(1) = ERROR_REAL
|
||||
end if
|
||||
else
|
||||
allocate(c % sqrtkT(1))
|
||||
c % sqrtkT = ERROR_REAL
|
||||
end if
|
||||
|
||||
! Add cell to dictionary
|
||||
call cell_dict % add_key(c % id, i)
|
||||
|
||||
|
|
@ -1867,6 +1901,7 @@ contains
|
|||
type(Node), pointer :: doc => null()
|
||||
type(Node), pointer :: node_mat => null()
|
||||
type(Node), pointer :: node_dens => null()
|
||||
type(Node), pointer :: node_temp => null()
|
||||
type(Node), pointer :: node_nuc => null()
|
||||
type(Node), pointer :: node_ele => null()
|
||||
type(Node), pointer :: node_sab => null()
|
||||
|
|
|
|||
103
src/math.F90
103
src/math.F90
|
|
@ -2,9 +2,25 @@ module math
|
|||
|
||||
use constants
|
||||
use random_lcg, only: prn
|
||||
use ISO_C_BINDING
|
||||
|
||||
implicit none
|
||||
|
||||
!===============================================================================
|
||||
! FADDEEVA_W evaluates the scaled complementary error function. This
|
||||
! interfaces with the MIT C library
|
||||
!===============================================================================
|
||||
|
||||
interface
|
||||
function faddeeva_w(z, relerr) bind(C, name='Faddeeva_w') result(w)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
complex(C_DOUBLE_COMPLEX), value :: z
|
||||
real(C_DOUBLE), value :: relerr
|
||||
complex(C_DOUBLE_COMPLEX) :: w
|
||||
end function faddeeva_w
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -659,4 +675,91 @@ contains
|
|||
|
||||
end function watt_spectrum
|
||||
|
||||
!===============================================================================
|
||||
! W acts as a front end to the MIT Faddeeva function, Faddeeva_w.
|
||||
!===============================================================================
|
||||
|
||||
function w(z) result(wv)
|
||||
complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at
|
||||
complex(8) :: wv ! The resulting w(z) value
|
||||
real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT Faddeeva
|
||||
|
||||
! Technically, the value we want is given by the equation:
|
||||
! w(z) = I/Pi * Integrate[Exp[-t^2]/(z-t), {t, -Infinity, Infinity}]
|
||||
! as shown in Equation 63 from Hwang, R. N. "A rigorous pole
|
||||
! representation of multilevel cross sections and its practical
|
||||
! applications." Nuclear Science and Engineering 96.3 (1987): 192-209.
|
||||
!
|
||||
! The MIT Faddeeva function evaluates w(z) = exp(-z^2)erfc(-iz). These
|
||||
! two forms of the Faddeeva function are related by a transformation.
|
||||
!
|
||||
! If we call the integral form w_int, and the function form w_fun:
|
||||
! For imag(z) > 0, w_int(z) = w_fun(z)
|
||||
! For imag(z) < 0, w_int(z) = -conjg(w_fun(conjg(z)))
|
||||
|
||||
relerr = ZERO
|
||||
if (aimag(z) > ZERO) then
|
||||
wv = faddeeva_w(z, relerr)
|
||||
else
|
||||
wv = -conjg(faddeeva_w(conjg(z), relerr))
|
||||
end if
|
||||
|
||||
end function w
|
||||
|
||||
!===============================================================================
|
||||
! BROADEN_N_POLYNOMIALS doppler broadens polynomials of the form
|
||||
! a/En + b/sqrt(En) + c + d sqrt(En) ...
|
||||
! exactly and quickly.
|
||||
!===============================================================================
|
||||
|
||||
subroutine broaden_n_polynomials(En, dopp, n, factors)
|
||||
real(8), intent(in) :: En ! Energy to evaluate at
|
||||
real(8), intent(in) :: dopp ! sqrt(atomic weight ratio / kT), kT given in eV.
|
||||
integer, intent(in) :: n ! number of components to polynomial
|
||||
real(8), intent(out):: factors(n) ! output leading coefficient
|
||||
|
||||
integer :: i
|
||||
|
||||
real(8) :: sqrtE ! Sqrt(energy)
|
||||
real(8) :: beta ! sqrt(atomic weight ratio * E / kT)
|
||||
real(8) :: half_inv_dopp2 ! 0.5 / dopp**2
|
||||
real(8) :: quarter_inv_dopp4 ! 0.25 / dopp**4
|
||||
real(8) :: erfbeta ! error function of beta
|
||||
real(8) :: exp_m_beta2 ! exp(-beta**2)
|
||||
|
||||
sqrtE = sqrt(En)
|
||||
beta = sqrtE * dopp
|
||||
half_inv_dopp2 = HALF / dopp**2
|
||||
quarter_inv_dopp4 = half_inv_dopp2**2
|
||||
|
||||
if (beta > 6.0_8) then
|
||||
! Save time, ERF(6) is 1 to machine precision.
|
||||
! beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon.
|
||||
erfBeta = ONE
|
||||
exp_m_beta2 = ZERO
|
||||
else
|
||||
erfBeta = erf(beta)
|
||||
exp_m_beta2 = exp(-beta**2)
|
||||
end if
|
||||
|
||||
! Assume that, for sure, we'll use a second order (1/E, 1/V, const)
|
||||
! fit, and no less.
|
||||
|
||||
factors(1) = erfbeta / En
|
||||
factors(2) = ONE / sqrtE
|
||||
factors(3) = factors(1) * (half_inv_dopp2 + En) + exp_m_beta2 / (beta * SQRT_PI)
|
||||
|
||||
! Perform recursive broadening of high order components
|
||||
do i = 1, n-3
|
||||
if (i /= 1) then
|
||||
factors(i+3) = -factors(i-1) * (i - ONE) * i * quarter_inv_dopp4 + &
|
||||
factors(i+1) * (En + (ONE + TWO * i) * half_inv_dopp2)
|
||||
else
|
||||
! Although it's mathematically identical, factors(0) will contain
|
||||
! nothing, and we don't want to have to worry about memory.
|
||||
factors(i+3) = factors(i+1)*(En + (ONE + TWO * i) * half_inv_dopp2)
|
||||
end if
|
||||
end do
|
||||
end subroutine broaden_n_polynomials
|
||||
|
||||
end module math
|
||||
|
|
|
|||
183
src/multipole.F90
Normal file
183
src/multipole.F90
Normal file
|
|
@ -0,0 +1,183 @@
|
|||
module multipole
|
||||
|
||||
use constants
|
||||
use global
|
||||
use hdf5
|
||||
use hdf5_interface
|
||||
use multipole_header, only: MultipoleArray, FIT_T, FIT_A, FIT_F, max_L, &
|
||||
max_poles, max_poly, MP_FISS, FORM_MLBW, FORM_RM
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! MULTIPOLE_READ Reads in a multipole HDF5 file with the original API
|
||||
! specification. Subject to change as the library format matures.
|
||||
!===============================================================================
|
||||
|
||||
subroutine multipole_read(filename, multipole, i_table)
|
||||
character(len=*), intent(in) :: filename ! Filename of the multipole library to load
|
||||
type(MultipoleArray), intent(out), target :: multipole ! The object to fill
|
||||
integer, intent(in) :: i_table ! index in nuclides/sab_tables
|
||||
|
||||
type(Nuclide), pointer :: nuc => null()
|
||||
|
||||
integer(HID_T) :: file_id
|
||||
integer(HID_T) :: group_id
|
||||
|
||||
! Intermediate loading components
|
||||
integer :: NMT
|
||||
integer :: i, j, k
|
||||
integer, allocatable :: MT(:)
|
||||
logical :: accumulated_fission
|
||||
character(len=3) :: MT_string
|
||||
character(len=24) :: MT_n ! Takes the form '/nuclide/reactions/MT???'
|
||||
integer :: is_fissionable
|
||||
|
||||
associate (nuc => nuclides(i_table))
|
||||
|
||||
! Open file for reading and move into the /isotope group
|
||||
file_id = file_open(filename, 'r', parallel=.true.)
|
||||
group_id = open_group(file_id, "/nuclide")
|
||||
|
||||
! Load in all the array size scalars
|
||||
call read_dataset(group_id, "length", multipole % length)
|
||||
call read_dataset(group_id, "windows", multipole % windows)
|
||||
call read_dataset(group_id, "num_l", multipole % num_l)
|
||||
call read_dataset(group_id, "fit_order", multipole % fit_order)
|
||||
call read_dataset(group_id, "max_w", multipole % max_w)
|
||||
call read_dataset(group_id, "fissionable", is_fissionable)
|
||||
if (is_fissionable == MP_FISS) then
|
||||
multipole % fissionable = .true.
|
||||
else
|
||||
multipole % fissionable = .false.
|
||||
end if
|
||||
call read_dataset(group_id, "formalism", multipole % formalism)
|
||||
|
||||
call read_dataset(group_id, "spacing", multipole % spacing)
|
||||
call read_dataset(group_id, "sqrtAWR", multipole % sqrtAWR)
|
||||
call read_dataset(group_id, "start_E", multipole % start_E)
|
||||
call read_dataset(group_id, "end_E", multipole % end_E)
|
||||
|
||||
! Allocate the multipole array components
|
||||
call multipole % allocate()
|
||||
|
||||
! Read in arrays
|
||||
call read_dataset(group_id, "data", multipole % data)
|
||||
call read_dataset(group_id, "pseudo_K0RS", multipole % pseudo_k0RS)
|
||||
call read_dataset(group_id, "l_value", multipole % l_value)
|
||||
call read_dataset(group_id, "w_start", multipole % w_start)
|
||||
call read_dataset(group_id, "w_end", multipole % w_end)
|
||||
call read_dataset(group_id, "broaden_poly", multipole % broaden_poly)
|
||||
|
||||
call read_dataset(group_id, "curvefit", multipole % curvefit)
|
||||
|
||||
! Delete ACE pointwise data
|
||||
call read_dataset(group_id, "n_grid", nuc % n_grid)
|
||||
|
||||
deallocate(nuc % energy)
|
||||
deallocate(nuc % total)
|
||||
deallocate(nuc % elastic)
|
||||
deallocate(nuc % fission)
|
||||
deallocate(nuc % nu_fission)
|
||||
deallocate(nuc % absorption)
|
||||
|
||||
allocate(nuc % energy(nuc % n_grid))
|
||||
allocate(nuc % total(nuc % n_grid))
|
||||
allocate(nuc % elastic(nuc % n_grid))
|
||||
allocate(nuc % fission(nuc % n_grid))
|
||||
allocate(nuc % nu_fission(nuc % n_grid))
|
||||
allocate(nuc % absorption(nuc % n_grid))
|
||||
|
||||
nuc % total = ZERO
|
||||
nuc % absorption = ZERO
|
||||
nuc % fission = ZERO
|
||||
|
||||
! Read in new energy axis (converting eV to MeV)
|
||||
call read_dataset(group_id, "energy_points", nuc % energy)
|
||||
nuc % energy = nuc % energy / 1.0D6
|
||||
|
||||
! Get count and list of MT tables
|
||||
call read_dataset(group_id, "MT_count", NMT)
|
||||
allocate(MT(NMT))
|
||||
|
||||
call read_dataset(group_id, "MT_list", MT)
|
||||
|
||||
call close_group(group_id)
|
||||
|
||||
accumulated_fission = .false.
|
||||
|
||||
! Loop over each MT entry and load it into a reaction.
|
||||
do i = 1, NMT
|
||||
write(MT_string, '(I3.3)') MT(i)
|
||||
MT_n = "/nuclide/reactions/MT" // MT_string
|
||||
|
||||
group_id = open_group(file_id, MT_n)
|
||||
|
||||
! Each MT needs to be treated slightly differently.
|
||||
select case (MT(i))
|
||||
case(ELASTIC)
|
||||
call read_dataset(group_id, "MT_sigma", nuc % elastic)
|
||||
nuc % total = nuc % total + nuc % elastic
|
||||
case(N_FISSION)
|
||||
call read_dataset(group_id, "MT_sigma", nuc % fission)
|
||||
nuc % total = nuc % total + nuc % fission
|
||||
nuc % absorption = nuc % absorption + nuc % fission
|
||||
accumulated_fission = .true.
|
||||
case default
|
||||
! Search through all of our secondary reactions
|
||||
do j = 1, nuc % n_reaction
|
||||
if (nuc % reactions(j) % MT == MT(i)) then
|
||||
! Match found
|
||||
|
||||
! Individual Fission components exist, so remove the combined
|
||||
! fission cross section.
|
||||
if ( (MT(i) == N_F .or. MT(i) == N_NF .or. MT(i) == N_2NF &
|
||||
.or. MT(i) == N_3NF) .and. accumulated_fission) then
|
||||
nuc % total = nuc % total - nuc % fission
|
||||
nuc % absorption = nuc % absorption - nuc % fission
|
||||
nuc % fission = 0.0_8
|
||||
accumulated_fission = .false.
|
||||
end if
|
||||
|
||||
deallocate(nuc % reactions(j) % sigma)
|
||||
allocate(nuc % reactions(j) % sigma(nuc % n_grid))
|
||||
|
||||
call read_dataset(group_id, "MT_sigma", nuc % reactions(j) % sigma)
|
||||
call read_dataset(group_id, "Q_value", nuc % reactions(j) % Q_value)
|
||||
call read_dataset(group_id, "threshold", nuc % reactions(j) % threshold)
|
||||
nuc % reactions(j) % threshold = 1 ! TODO: reconsider implications.
|
||||
nuc % reactions(j) % Q_value = nuc % reactions(j) % Q_value / 1.0D6
|
||||
|
||||
! Accumulate total
|
||||
if (MT(i) /= N_LEVEL .and. MT(i) <= N_DA) then
|
||||
nuc % total = nuc % total + nuc % reactions(j) % sigma
|
||||
end if
|
||||
|
||||
! Accumulate absorption
|
||||
if (MT(i) >= N_GAMMA .and. MT(i) <= N_DA) then
|
||||
nuc % absorption = nuc % absorption + nuc % reactions(j) % sigma
|
||||
end if
|
||||
|
||||
! Accumulate fission (if needed)
|
||||
if ( (MT(i) == N_F .or. MT(i) == N_NF .or. MT(i) == N_2NF &
|
||||
.or. MT(i) == N_3NF) ) then
|
||||
nuc % fission = nuc % fission + nuc % reactions(j) % sigma
|
||||
nuc % absorption = nuc % absorption + nuc % reactions(j) % sigma
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end select
|
||||
|
||||
call close_group(group_id)
|
||||
end do
|
||||
|
||||
! Close file
|
||||
call file_close(file_id)
|
||||
|
||||
end associate
|
||||
|
||||
end subroutine multipole_read
|
||||
|
||||
end module multipole
|
||||
133
src/multipole_header.F90
Normal file
133
src/multipole_header.F90
Normal file
|
|
@ -0,0 +1,133 @@
|
|||
module multipole_header
|
||||
|
||||
implicit none
|
||||
|
||||
!========================================================================
|
||||
! Multipole related constants
|
||||
|
||||
! Formalisms
|
||||
integer, parameter :: FORM_MLBW = 2, &
|
||||
FORM_RM = 3, &
|
||||
FORM_RML = 7
|
||||
|
||||
! Constants that determine which value to access
|
||||
integer, parameter :: MP_EA = 1 ! Pole
|
||||
|
||||
! Reich-Moore indices
|
||||
integer, parameter :: RM_RT = 2, & ! Residue total
|
||||
RM_RA = 3, & ! Residue absorption
|
||||
RM_RF = 4 ! Residue fission
|
||||
|
||||
! Multi-level Breit Wigner indices
|
||||
integer, parameter :: MLBW_RT = 2, & ! Residue total
|
||||
MLBW_RX = 3, & ! Residue compettitive
|
||||
MLBW_RA = 4, & ! Residue absorption
|
||||
MLBW_RF = 5 ! Residue fission
|
||||
|
||||
! Polynomial fit indices
|
||||
integer, parameter :: FIT_T = 1, & ! Total
|
||||
FIT_A = 2, & ! Absorption
|
||||
FIT_F = 3 ! Fission
|
||||
|
||||
! Value of 'true' when checking if nuclide is fissionable
|
||||
integer, parameter :: MP_FISS = 1
|
||||
|
||||
! These variables store the maximum value from every nuclide in order
|
||||
! to preallocate some arrays to improve performance.
|
||||
integer :: max_poly ! Maximum number of polynomials we expect
|
||||
integer :: max_poles ! Maximum number of poles in the problem for allocation
|
||||
integer :: max_L ! Maximum L value for allocation
|
||||
|
||||
!===============================================================================
|
||||
! MULTIPOLE contains all the components needed for the windowed multipole
|
||||
! temperature dependent cross section libraries for the resolved resonance
|
||||
! region.
|
||||
!===============================================================================
|
||||
|
||||
type MultipoleArray
|
||||
|
||||
!=========================================================================
|
||||
! Isotope Properties
|
||||
logical :: fissionable = .false. ! Is this isotope fissionable?
|
||||
integer :: length ! Number of poles
|
||||
integer, allocatable :: l_value(:) ! The l index of the pole
|
||||
real(8), allocatable :: pseudo_k0RS(:) ! The value (sqrt(2*mass neutron)/reduced planck constant) * AWR/(AWR + 1) * scattering radius for each l
|
||||
complex(8), allocatable :: data(:,:) ! Contains all of the pole-residue data
|
||||
real(8) :: sqrtAWR ! Square root of the atomic weight ratio
|
||||
|
||||
!=========================================================================
|
||||
! Windows
|
||||
|
||||
integer :: windows ! Number of windows
|
||||
integer :: fit_order ! Order of the fit. 1 linear, 2 quadratic, etc.
|
||||
real(8) :: start_E ! Start energy for the windows
|
||||
real(8) :: end_E ! End energy for the windows
|
||||
real(8) :: spacing ! The actual spacing in sqrt(E) space.
|
||||
! spacing = sqrt(multipole_w%endE - multipole_w%startE)/multipole_w%windows
|
||||
integer, allocatable :: w_start(:) ! Contains the index of the pole at the start of the window
|
||||
integer, allocatable :: w_end(:) ! Contains the index of the pole at the end of the window
|
||||
real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function. (reaction type, coeff index, window index)
|
||||
|
||||
integer, allocatable :: broaden_poly(:) ! if 1, broaden, if 0, don't.
|
||||
|
||||
!=========================================================================
|
||||
! Storage Helpers
|
||||
integer :: num_l
|
||||
integer :: max_w
|
||||
|
||||
integer :: formalism
|
||||
|
||||
contains
|
||||
procedure :: allocate => multipole_allocate ! Allocates Multipole
|
||||
end type MultipoleArray
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! MULTIPOLE_ALLOCATE allocates necessary data for Multipole.
|
||||
!===============================================================================
|
||||
|
||||
subroutine multipole_allocate(multipole)
|
||||
class(MultipoleArray), intent(inout) :: multipole ! Multipole object to allocate.
|
||||
|
||||
! This function assumes length, numL, fissionable, windows, fitorder,
|
||||
! and formalism are known
|
||||
|
||||
! Allocate the pole-residue storage.
|
||||
! MLBW has one more pole than Reich-Moore, and fissionable nuclides
|
||||
! have further one more.
|
||||
if (multipole % formalism == FORM_MLBW) then
|
||||
if (multipole % fissionable) then
|
||||
allocate(multipole % data(5, multipole % length))
|
||||
else
|
||||
allocate(multipole % data(4, multipole % length))
|
||||
end if
|
||||
else if (multipole % formalism == FORM_RM) then
|
||||
if (multipole % fissionable) then
|
||||
allocate(multipole % data(4, multipole % length))
|
||||
else
|
||||
allocate(multipole % data(3, multipole % length))
|
||||
end if
|
||||
end if
|
||||
|
||||
! Allocate the l value table for each pole-residue set.
|
||||
allocate(multipole % l_value(multipole % length))
|
||||
|
||||
! Allocate the table of pseudo_k0RS values at each l.
|
||||
allocate(multipole % pseudo_k0RS(multipole % num_l))
|
||||
|
||||
! Allocate window start, window end
|
||||
allocate(multipole % w_start(multipole % windows))
|
||||
allocate(multipole % w_end(multipole % windows))
|
||||
|
||||
! Allocate broaden_poly
|
||||
allocate(multipole % broaden_poly(multipole % windows))
|
||||
|
||||
! Allocate curvefit
|
||||
if(multipole % fissionable) then
|
||||
allocate(multipole % curvefit(FIT_F, multipole % fit_order+1, multipole % windows))
|
||||
else
|
||||
allocate(multipole % curvefit(FIT_A, multipole % fit_order+1, multipole % windows))
|
||||
end if
|
||||
end subroutine
|
||||
end module multipole_header
|
||||
|
|
@ -1434,7 +1434,7 @@ contains
|
|||
do i = 1, n
|
||||
! If the cell matches the goal and the offset matches final, write to the
|
||||
! geometry stack
|
||||
if (univ % cells(i) == goal .AND. offset == final) then
|
||||
if (univ % cells(i) == goal .and. offset == final) then
|
||||
c => cells(univ % cells(i))
|
||||
path = trim(path) // "->" // to_str(c % id)
|
||||
return
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@ module particle_header
|
|||
|
||||
use bank_header, only: Bank
|
||||
use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, &
|
||||
MAX_DELAYED_GROUPS
|
||||
MAX_DELAYED_GROUPS, ERROR_REAL
|
||||
use error, only: fatal_error
|
||||
use geometry_header, only: BASE_UNIVERSE
|
||||
|
||||
|
|
@ -79,6 +79,9 @@ module particle_header
|
|||
integer :: material ! index for current material
|
||||
integer :: last_material ! index for last material
|
||||
|
||||
! Temperature of the current cell
|
||||
real(8) :: sqrtkT ! sqrt(k_Boltzmann * temperature) in MeV
|
||||
|
||||
! Statistical data
|
||||
integer :: n_collision ! # of collisions
|
||||
|
||||
|
|
@ -86,7 +89,7 @@ module particle_header
|
|||
logical :: write_track = .false.
|
||||
|
||||
! Secondary particles created
|
||||
integer :: n_secondary = 0
|
||||
integer(8) :: n_secondary = 0
|
||||
type(Bank) :: secondary_bank(MAX_SECONDARY)
|
||||
|
||||
contains
|
||||
|
|
@ -124,6 +127,7 @@ contains
|
|||
this % absorb_wgt = ZERO
|
||||
this % n_bank = 0
|
||||
this % wgt_bank = ZERO
|
||||
this % sqrtkT = ERROR_REAL
|
||||
this % n_collision = 0
|
||||
this % fission = .false.
|
||||
this % delayed_group = 0
|
||||
|
|
|
|||
|
|
@ -88,9 +88,14 @@ contains
|
|||
! change when sampling fission sites. The following block handles all
|
||||
! absorption (including fission)
|
||||
|
||||
if (nuc % fissionable .and. run_mode == MODE_EIGENVALUE) then
|
||||
if (nuc % fissionable) then
|
||||
call sample_fission(i_nuclide, i_reaction)
|
||||
call create_fission_sites(p, i_nuclide, i_reaction)
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank)
|
||||
elseif (run_mode == MODE_FIXEDSOURCE) then
|
||||
call create_fission_sites(p, i_nuclide, i_reaction, &
|
||||
p % secondary_bank, p % n_secondary)
|
||||
end if
|
||||
end if
|
||||
|
||||
! If survival biasing is being used, the following subroutine adjusts the
|
||||
|
|
@ -1071,10 +1076,12 @@ contains
|
|||
! neutrons produced from fission and creates appropriate bank sites.
|
||||
!===============================================================================
|
||||
|
||||
subroutine create_fission_sites(p, i_nuclide, i_reaction)
|
||||
subroutine create_fission_sites(p, i_nuclide, i_reaction, bank_array, size_bank)
|
||||
type(Particle), intent(inout) :: p
|
||||
integer, intent(in) :: i_nuclide
|
||||
integer, intent(in) :: i_reaction
|
||||
type(Bank), intent(inout) :: bank_array(:)
|
||||
integer(8), intent(inout) :: size_bank
|
||||
|
||||
integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born
|
||||
integer :: i ! loop index
|
||||
|
|
@ -1123,27 +1130,35 @@ contains
|
|||
nu = int(nu_t) + 1
|
||||
end if
|
||||
|
||||
! Check for fission bank size getting hit
|
||||
if (n_bank + nu > size(fission_bank)) then
|
||||
if (master) call warning("Maximum number of sites in fission bank &
|
||||
&reached. This can result in irreproducible results using different &
|
||||
&numbers of processes/threads.")
|
||||
! Check for bank size getting hit. For fixed source calculations, this is a
|
||||
! fatal error. For eigenvalue calculations, it just means that k-effective
|
||||
! was too high for a single batch.
|
||||
if (size_bank + nu > size(bank_array)) then
|
||||
if (run_mode == MODE_FIXEDSOURCE) then
|
||||
call fatal_error("Secondary particle bank size limit reached. If you &
|
||||
&are running a subcritical multiplication problem, k-effective &
|
||||
&may be too close to one.")
|
||||
else
|
||||
if (master) call warning("Maximum number of sites in fission bank &
|
||||
&reached. This can result in irreproducible results using different &
|
||||
&numbers of processes/threads.")
|
||||
end if
|
||||
end if
|
||||
|
||||
! Bank source neutrons
|
||||
if (nu == 0 .or. n_bank == size(fission_bank)) return
|
||||
if (nu == 0 .or. size_bank == size(bank_array)) return
|
||||
|
||||
! Initialize counter of delayed neutrons encountered for each delayed group
|
||||
! to zero.
|
||||
nu_d(:) = 0
|
||||
|
||||
p % fission = .true. ! Fission neutrons will be banked
|
||||
do i = int(n_bank,4) + 1, int(min(n_bank + nu, int(size(fission_bank),8)),4)
|
||||
do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4)
|
||||
! Bank source neutrons by copying particle data
|
||||
fission_bank(i) % xyz = p % coord(1) % xyz
|
||||
bank_array(i) % xyz = p % coord(1) % xyz
|
||||
|
||||
! Set weight of fission bank site
|
||||
fission_bank(i) % wgt = ONE/weight
|
||||
bank_array(i) % wgt = ONE/weight
|
||||
|
||||
! Sample cosine of angle -- fission neutrons are always emitted
|
||||
! isotropically. Sometimes in ACE data, fission reactions actually have
|
||||
|
|
@ -1153,17 +1168,17 @@ contains
|
|||
|
||||
! Sample azimuthal angle uniformly in [0,2*pi)
|
||||
phi = TWO*PI*prn()
|
||||
fission_bank(i) % uvw(1) = mu
|
||||
fission_bank(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
|
||||
fission_bank(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
|
||||
bank_array(i) % uvw(1) = mu
|
||||
bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
|
||||
bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
|
||||
|
||||
! Sample secondary energy distribution for fission reaction and set energy
|
||||
! in fission bank
|
||||
fission_bank(i) % E = sample_fission_energy(nuc, nuc%reactions(&
|
||||
i_reaction), p)
|
||||
bank_array(i) % E = sample_fission_energy(nuc, &
|
||||
nuc % reactions(i_reaction), p)
|
||||
|
||||
! Set the delayed group of the neutron
|
||||
fission_bank(i) % delayed_group = p % delayed_group
|
||||
bank_array(i) % delayed_group = p % delayed_group
|
||||
|
||||
! Increment the number of neutrons born delayed
|
||||
if (p % delayed_group > 0) then
|
||||
|
|
@ -1172,7 +1187,7 @@ contains
|
|||
end do
|
||||
|
||||
! increment number of bank sites
|
||||
n_bank = min(n_bank + nu, int(size(fission_bank),8))
|
||||
size_bank = min(size_bank + nu, int(size(bank_array),8))
|
||||
|
||||
! Store total and delayed weight banked for analog fission tallies
|
||||
p % n_bank = nu
|
||||
|
|
|
|||
|
|
@ -9,6 +9,8 @@ element geometry {
|
|||
(element material { ( xsd:int | "void" )+ } |
|
||||
attribute material { ( xsd:int | "void" )+ })
|
||||
) &
|
||||
(element temperature { list { xsd:double+ } } |
|
||||
attribute temperature { list { xsd:double+ } } )? &
|
||||
(element region { xsd:string } | attribute region { xsd:string })? &
|
||||
(element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? &
|
||||
(element translation { list { xsd:double+ } } | attribute translation { list { xsd:double+ } })?
|
||||
|
|
|
|||
|
|
@ -64,6 +64,24 @@
|
|||
</attribute>
|
||||
</choice>
|
||||
</choice>
|
||||
<optional>
|
||||
<choice>
|
||||
<element name="temperature">
|
||||
<list>
|
||||
<oneOrMore>
|
||||
<data type="double"/>
|
||||
</oneOrMore>
|
||||
</list>
|
||||
</element>
|
||||
<attribute name="temperature">
|
||||
<list>
|
||||
<oneOrMore>
|
||||
<data type="double"/>
|
||||
</oneOrMore>
|
||||
</list>
|
||||
</attribute>
|
||||
</choice>
|
||||
</optional>
|
||||
<optional>
|
||||
<choice>
|
||||
<element name="region">
|
||||
|
|
|
|||
|
|
@ -108,6 +108,7 @@ contains
|
|||
integer :: i, j, k, m
|
||||
integer, allocatable :: lattice_universes(:,:,:)
|
||||
integer, allocatable :: cell_materials(:)
|
||||
real(8), allocatable :: cell_temperatures(:)
|
||||
integer(HID_T) :: geom_group
|
||||
integer(HID_T) :: cells_group, cell_group
|
||||
integer(HID_T) :: surfaces_group, surface_group
|
||||
|
|
@ -151,6 +152,7 @@ contains
|
|||
select case (c%type)
|
||||
case (CELL_NORMAL)
|
||||
call write_dataset(cell_group, "fill_type", "normal")
|
||||
|
||||
if (size(c % material) == 1) then
|
||||
if (c % material(1) == MATERIAL_VOID) then
|
||||
call write_dataset(cell_group, "material", MATERIAL_VOID)
|
||||
|
|
@ -171,6 +173,12 @@ contains
|
|||
deallocate(cell_materials)
|
||||
end if
|
||||
|
||||
allocate(cell_temperatures(size(c % sqrtkT)))
|
||||
cell_temperatures(:) = c % sqrtkT(:)
|
||||
cell_temperatures(:) = cell_temperatures(:)**2 / K_BOLTZMANN
|
||||
call write_dataset(cell_group, "temperature", cell_temperatures)
|
||||
deallocate(cell_temperatures)
|
||||
|
||||
case (CELL_FILL)
|
||||
call write_dataset(cell_group, "fill_type", "universe")
|
||||
call write_dataset(cell_group, "fill", universes(c%fill)%id)
|
||||
|
|
|
|||
|
|
@ -213,6 +213,7 @@ contains
|
|||
if (p % n_secondary > 0) then
|
||||
call p % initialize_from_source(p % secondary_bank(p % n_secondary))
|
||||
p % n_secondary = p % n_secondary - 1
|
||||
n_event = 0
|
||||
|
||||
! Enter new particle in particle track file
|
||||
if (p % write_track) call add_particle_track()
|
||||
|
|
|
|||
1
tests/test_asymmetric_lattice/inputs_true.dat
Normal file
1
tests/test_asymmetric_lattice/inputs_true.dat
Normal file
|
|
@ -0,0 +1 @@
|
|||
b9b4222c4beea80fe6083590f6b785303d174972d80671fb661bac8e030db6f4a61648240cfad6162799361fc0e08a23c61d31aff844d978528d6dad5b5fbc63
|
||||
1
tests/test_asymmetric_lattice/results_true.dat
Normal file
1
tests/test_asymmetric_lattice/results_true.dat
Normal file
|
|
@ -0,0 +1 @@
|
|||
b5f96919ca474cd1c9c9d0acde3b8aac4a1cf636443c72a38b6c5a4221a8ce3e90182aaef2f664e44b9175ca257a89db2328b63e19388ee0e5006de4b3d92ce6
|
||||
129
tests/test_asymmetric_lattice/test_asymmetric_lattice.py
Normal file
129
tests/test_asymmetric_lattice/test_asymmetric_lattice.py
Normal file
|
|
@ -0,0 +1,129 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import os
|
||||
import sys
|
||||
import glob
|
||||
import hashlib
|
||||
sys.path.insert(0, os.pardir)
|
||||
from testing_harness import PyAPITestHarness
|
||||
import openmc
|
||||
from openmc.source import Source
|
||||
from openmc.stats import Box
|
||||
|
||||
|
||||
class AsymmetricLatticeTestHarness(PyAPITestHarness):
|
||||
|
||||
def _build_inputs(self):
|
||||
"""Build an axis-asymmetric lattice of fuel assemblies"""
|
||||
|
||||
# Build full core geometry from underlying input set
|
||||
self._input_set.build_default_materials_and_geometry()
|
||||
|
||||
# Extract all universes from the full core geometry
|
||||
geometry = self._input_set.geometry.geometry
|
||||
all_univs = geometry.get_all_universes()
|
||||
|
||||
# Extract universes encapsulating fuel and water assemblies
|
||||
water = all_univs[7]
|
||||
fuel = all_univs[8]
|
||||
|
||||
# Construct a 3x3 lattice of fuel assemblies
|
||||
core_lat = openmc.RectLattice(name='3x3 Core Lattice', lattice_id=202)
|
||||
core_lat.dimension = (3, 3)
|
||||
core_lat.lower_left = (-32.13, -32.13)
|
||||
core_lat.pitch = (21.42, 21.42)
|
||||
core_lat.universes = [[fuel, water, water],
|
||||
[fuel, fuel, fuel],
|
||||
[water, water, water]]
|
||||
|
||||
# Create bounding surfaces
|
||||
min_x = openmc.XPlane(x0=-32.13, boundary_type='reflective')
|
||||
max_x = openmc.XPlane(x0=+32.13, boundary_type='reflective')
|
||||
min_y = openmc.YPlane(y0=-32.13, boundary_type='reflective')
|
||||
max_y = openmc.YPlane(y0=+32.13, boundary_type='reflective')
|
||||
min_z = openmc.ZPlane(z0=0, boundary_type='reflective')
|
||||
max_z = openmc.ZPlane(z0=+32.13, boundary_type='reflective')
|
||||
|
||||
# Define root universe
|
||||
root_univ = openmc.Universe(universe_id=0, name='root universe')
|
||||
root_cell = openmc.Cell(cell_id=1)
|
||||
root_cell.region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z
|
||||
root_cell.fill = core_lat
|
||||
root_univ.add_cell(root_cell)
|
||||
|
||||
# Over-ride geometry in the input set with this 3x3 lattice
|
||||
self._input_set.geometry.geometry.root_universe = root_univ
|
||||
|
||||
# Initialize a "distribcell" filter for the fuel pin cell
|
||||
distrib_filter = openmc.Filter(type='distribcell', bins=[27])
|
||||
|
||||
# Initialize the tallies
|
||||
tally = openmc.Tally(name='distribcell tally', tally_id=27)
|
||||
tally.add_filter(distrib_filter)
|
||||
tally.add_score('nu-fission')
|
||||
|
||||
# Initialize the tallies file
|
||||
tallies_file = openmc.TalliesFile()
|
||||
tallies_file.add_tally(tally)
|
||||
|
||||
# Assign the tallies file to the input set
|
||||
self._input_set.tallies = tallies_file
|
||||
|
||||
# Build default settings
|
||||
self._input_set.build_default_settings()
|
||||
|
||||
# Specify summary output and correct source sampling box
|
||||
source = Source(space=Box([-32, -32, 0], [32, 32, 32]))
|
||||
source.space.only_fissionable = True
|
||||
self._input_set.settings.source = source
|
||||
self._input_set.settings.output = {'summary': True}
|
||||
|
||||
# Write input XML files
|
||||
self._input_set.export()
|
||||
|
||||
def _get_results(self, hash_output=True):
|
||||
"""Digest info in statepoint and summary and return as a string."""
|
||||
|
||||
# Read the statepoint file
|
||||
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
|
||||
sp = openmc.StatePoint(statepoint)
|
||||
|
||||
# Read the summary file
|
||||
summary = glob.glob(os.path.join(os.getcwd(), 'summary.h5'))[0]
|
||||
su = openmc.Summary(summary)
|
||||
sp.link_with_summary(su)
|
||||
|
||||
# Extract the tally of interest
|
||||
tally = sp.get_tally(name='distribcell tally')
|
||||
|
||||
# Create a string of all mean, std. dev. values for both tallies
|
||||
outstr = ''
|
||||
outstr += ', '.join(map(str, tally.mean.flatten())) + '\n'
|
||||
outstr += ', '.join(map(str, tally.std_dev.flatten())) + '\n'
|
||||
|
||||
# Extract fuel assembly lattices from the summary
|
||||
all_cells = su.openmc_geometry.get_all_cells()
|
||||
fuel = all_cells[80].fill
|
||||
core = all_cells[1].fill
|
||||
|
||||
# Append a string of lattice distribcell offsets to the string
|
||||
outstr += ', '.join(map(str, fuel.offsets.flatten())) + '\n'
|
||||
outstr += ', '.join(map(str, core.offsets.flatten())) + '\n'
|
||||
|
||||
# Hash the results if necessary
|
||||
if hash_output:
|
||||
sha512 = hashlib.sha512()
|
||||
sha512.update(outstr.encode('utf-8'))
|
||||
outstr = sha512.hexdigest()
|
||||
|
||||
return outstr
|
||||
|
||||
def _cleanup(self):
|
||||
super(AsymmetricLatticeTestHarness, self)._cleanup()
|
||||
f = os.path.join(os.getcwd(), 'tallies.xml')
|
||||
if os.path.exists(f): os.remove(f)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
harness = AsymmetricLatticeTestHarness('statepoint.10.h5', True)
|
||||
harness.main()
|
||||
|
|
@ -5,6 +5,7 @@ Cell
|
|||
Name =
|
||||
Material = [2, 3, void, 2]
|
||||
Region = -10000
|
||||
Temperature = [ 293.60594237 293.60594237 0. 293.60594237]
|
||||
Rotation = None
|
||||
Translation = None
|
||||
Offset = None
|
||||
|
|
|
|||
|
|
@ -4,6 +4,7 @@
|
|||
<material id="1">
|
||||
<density value="7.5" units="g/cc" />
|
||||
<nuclide name="O-16" xs="71c" ao="1.0" />
|
||||
<nuclide name="U-238" xs="71c" ao="0.0001" />
|
||||
</material>
|
||||
|
||||
</materials>
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
tally 1:
|
||||
4.538791E+02
|
||||
2.073271E+04
|
||||
4.563929E+02
|
||||
2.091711E+04
|
||||
leakage:
|
||||
9.830000E+00
|
||||
9.663900E+00
|
||||
9.780000E+00
|
||||
9.566400E+00
|
||||
|
|
|
|||
1
tests/test_multipole/inputs_true.dat
Normal file
1
tests/test_multipole/inputs_true.dat
Normal file
|
|
@ -0,0 +1 @@
|
|||
5c1cec635da5c4c869bdf58f62924a4cb1648e4ddecf0ce71b823cf45178767ba7f2e089b76d36e8616b1b21ffa43b32ab1b17d20bb74120f900b9e3e9ab9bcc
|
||||
12
tests/test_multipole/results_true.dat
Normal file
12
tests/test_multipole/results_true.dat
Normal file
|
|
@ -0,0 +1,12 @@
|
|||
k-combined:
|
||||
1.445285E+00 9.521660E-03
|
||||
Cell
|
||||
ID = 11
|
||||
Name =
|
||||
Material = 2
|
||||
Region = -10000
|
||||
Temperature = [ 500. 0. 700. 800.]
|
||||
Rotation = None
|
||||
Translation = None
|
||||
Offset = None
|
||||
Distribcell index= 1
|
||||
133
tests/test_multipole/test_multipole.py
Normal file
133
tests/test_multipole/test_multipole.py
Normal file
|
|
@ -0,0 +1,133 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import os
|
||||
import sys
|
||||
sys.path.insert(0, os.pardir)
|
||||
from testing_harness import TestHarness, PyAPITestHarness
|
||||
import openmc
|
||||
from openmc.stats import Box
|
||||
from openmc.source import Source
|
||||
|
||||
|
||||
class DistribmatTestHarness(PyAPITestHarness):
|
||||
def _build_inputs(self):
|
||||
####################
|
||||
# Materials
|
||||
####################
|
||||
|
||||
moderator = openmc.Material(material_id=1)
|
||||
moderator.set_density('g/cc', 1.0)
|
||||
moderator.add_nuclide('H-1', 2.0)
|
||||
moderator.add_nuclide('O-16', 1.0)
|
||||
|
||||
dense_fuel = openmc.Material(material_id=2)
|
||||
dense_fuel.set_density('g/cc', 4.5)
|
||||
dense_fuel.add_nuclide('U-235', 1.0)
|
||||
|
||||
mats_file = openmc.MaterialsFile()
|
||||
mats_file.default_xs = '71c'
|
||||
mats_file.add_materials([moderator, dense_fuel])
|
||||
mats_file.export_to_xml()
|
||||
|
||||
|
||||
####################
|
||||
# Geometry
|
||||
####################
|
||||
|
||||
c1 = openmc.Cell(cell_id=1)
|
||||
c1.fill = moderator
|
||||
mod_univ = openmc.Universe(universe_id=1)
|
||||
mod_univ.add_cell(c1)
|
||||
|
||||
r0 = openmc.ZCylinder(R=0.3)
|
||||
c11 = openmc.Cell(cell_id=11)
|
||||
c11.region = -r0
|
||||
c11.fill = dense_fuel
|
||||
c11.temperature = [500, 0, 700, 800]
|
||||
c12 = openmc.Cell(cell_id=12)
|
||||
c12.region = +r0
|
||||
c12.fill = moderator
|
||||
fuel_univ = openmc.Universe(universe_id=11)
|
||||
fuel_univ.add_cells((c11, c12))
|
||||
|
||||
lat = openmc.RectLattice(lattice_id=101)
|
||||
lat.dimension = [2, 2]
|
||||
lat.lower_left = [-2.0, -2.0]
|
||||
lat.pitch = [2.0, 2.0]
|
||||
lat.universes = [[fuel_univ]*2]*2
|
||||
lat.outer = mod_univ
|
||||
|
||||
x0 = openmc.XPlane(x0=-3.0)
|
||||
x1 = openmc.XPlane(x0=3.0)
|
||||
y0 = openmc.YPlane(y0=-3.0)
|
||||
y1 = openmc.YPlane(y0=3.0)
|
||||
for s in [x0, x1, y0, y1]:
|
||||
s.boundary_type = 'reflective'
|
||||
c101 = openmc.Cell(cell_id=101)
|
||||
c101.region = +x0 & -x1 & +y0 & -y1
|
||||
c101.fill = lat
|
||||
root_univ = openmc.Universe(universe_id=0)
|
||||
root_univ.add_cell(c101)
|
||||
|
||||
geometry = openmc.Geometry()
|
||||
geometry.root_universe = root_univ
|
||||
geo_file = openmc.GeometryFile()
|
||||
geo_file.geometry = geometry
|
||||
geo_file.export_to_xml()
|
||||
|
||||
|
||||
####################
|
||||
# Settings
|
||||
####################
|
||||
|
||||
sets_file = openmc.SettingsFile()
|
||||
sets_file.batches = 5
|
||||
sets_file.inactive = 0
|
||||
sets_file.particles = 1000
|
||||
sets_file.source = Source(space=Box([-1, -1, -1], [1, 1, 1]))
|
||||
sets_file.output = {'summary': True}
|
||||
sets_file.export_to_xml()
|
||||
|
||||
|
||||
####################
|
||||
# Plots
|
||||
####################
|
||||
|
||||
plots_file = openmc.PlotsFile()
|
||||
|
||||
plot = openmc.Plot(plot_id=1)
|
||||
plot.basis = 'xy'
|
||||
plot.color = 'cell'
|
||||
plot.filename = 'cellplot'
|
||||
plot.origin = (0, 0, 0)
|
||||
plot.width = (7, 7)
|
||||
plot.pixels = (400, 400)
|
||||
plots_file.add_plot(plot)
|
||||
|
||||
plot = openmc.Plot(plot_id=2)
|
||||
plot.basis = 'xy'
|
||||
plot.color = 'mat'
|
||||
plot.filename = 'matplot'
|
||||
plot.origin = (0, 0, 0)
|
||||
plot.width = (7, 7)
|
||||
plot.pixels = (400, 400)
|
||||
plots_file.add_plot(plot)
|
||||
|
||||
plots_file.export_to_xml()
|
||||
|
||||
def _get_results(self):
|
||||
outstr = super(DistribmatTestHarness, self)._get_results()
|
||||
su = openmc.Summary('summary.h5')
|
||||
outstr += str(su.get_cell_by_id(11))
|
||||
return outstr
|
||||
|
||||
def _cleanup(self):
|
||||
f = os.path.join(os.getcwd(), 'plots.xml')
|
||||
if os.path.exists(f):
|
||||
os.remove(f)
|
||||
super(DistribmatTestHarness, self)._cleanup()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
harness = DistribmatTestHarness('statepoint.5.*')
|
||||
harness.main()
|
||||
Loading…
Add table
Add a link
Reference in a new issue