mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #861 from paulromano/multipole-tally-fix
Fix bug tallying reaction rates with multipole
This commit is contained in:
Colin Josey
GitHub
commit
9d22f149d7
8 changed files with 147 additions and 125 deletions
|
|
@ -235,12 +235,6 @@ endif()
|
|||
#===============================================================================
|
||||
|
||||
add_library(pugixml src/pugixml/pugixml_c.cpp src/pugixml/pugixml.cpp)
|
||||
if(CMAKE_VERSION VERSION_LESS 3.1)
|
||||
target_compile_options(pugixml PRIVATE -std=c++11)
|
||||
else()
|
||||
set_property(TARGET pugixml PROPERTY CXX_STANDARD 11)
|
||||
endif()
|
||||
|
||||
add_library(pugixml_fortran src/pugixml/pugixml_f.F90)
|
||||
target_link_libraries(pugixml_fortran pugixml)
|
||||
|
||||
|
|
|
|||
|
|
@ -997,7 +997,10 @@ class Settings(object):
|
|||
for key, value in sorted(self.temperature.items()):
|
||||
element = ET.SubElement(root,
|
||||
"temperature_{}".format(key))
|
||||
element.text = str(value)
|
||||
if isinstance(value, bool):
|
||||
element.text = str(value).lower()
|
||||
else:
|
||||
element.text = str(value)
|
||||
|
||||
def _create_threads_subelement(self, root):
|
||||
if self._threads is not None:
|
||||
|
|
|
|||
|
|
@ -44,7 +44,7 @@
|
|||
// #define PUGIXML_HEADER_ONLY
|
||||
|
||||
// Uncomment this to enable long long support
|
||||
// #define PUGIXML_HAS_LONG_LONG
|
||||
#define PUGIXML_HAS_LONG_LONG
|
||||
|
||||
#endif
|
||||
|
||||
|
|
|
|||
|
|
@ -1169,16 +1169,28 @@ contains
|
|||
! Retrieve temperature and energy grid index and interpolation
|
||||
! factor
|
||||
i_temp = micro_xs(i_nuclide) % index_temp
|
||||
i_energy = micro_xs(i_nuclide) % index_grid
|
||||
f = micro_xs(i_nuclide) % interp_factor
|
||||
if (i_temp > 0) then
|
||||
i_energy = micro_xs(i_nuclide) % index_grid
|
||||
f = micro_xs(i_nuclide) % interp_factor
|
||||
|
||||
associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp))
|
||||
if (i_energy >= xs % threshold) then
|
||||
score = ((ONE - f) * xs % value(i_energy - &
|
||||
xs % threshold + 1) + f * xs % value(i_energy - &
|
||||
xs % threshold + 2)) * atom_density * flux
|
||||
associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp))
|
||||
if (i_energy >= xs % threshold) then
|
||||
score = ((ONE - f) * xs % value(i_energy - &
|
||||
xs % threshold + 1) + f * xs % value(i_energy - &
|
||||
xs % threshold + 2)) * atom_density * flux
|
||||
end if
|
||||
end associate
|
||||
else
|
||||
! This block is reached if multipole is turned on and we're in
|
||||
! the resolved range. For (n,gamma), use absorption -
|
||||
! fission. For everything else, assume it's zero.
|
||||
if (score_bin == N_GAMMA) then
|
||||
score = (micro_xs(i_nuclide) % absorption - &
|
||||
micro_xs(i_nuclide) % fission) * atom_density * flux
|
||||
else
|
||||
score = ZERO
|
||||
end if
|
||||
end associate
|
||||
end if
|
||||
end if
|
||||
|
||||
else
|
||||
|
|
@ -1195,16 +1207,28 @@ contains
|
|||
! Retrieve temperature and energy grid index and interpolation
|
||||
! factor
|
||||
i_temp = micro_xs(i_nuc) % index_temp
|
||||
i_energy = micro_xs(i_nuc) % index_grid
|
||||
f = micro_xs(i_nuc) % interp_factor
|
||||
if (i_temp > 0) then
|
||||
i_energy = micro_xs(i_nuc) % index_grid
|
||||
f = micro_xs(i_nuc) % interp_factor
|
||||
|
||||
associate (xs => nuclides(i_nuc) % reactions(m) % xs(i_temp))
|
||||
if (i_energy >= xs % threshold) then
|
||||
score = score + ((ONE - f) * xs % value(i_energy - &
|
||||
xs % threshold + 1) + f * xs % value(i_energy - &
|
||||
xs % threshold + 2)) * atom_density_ * flux
|
||||
associate (xs => nuclides(i_nuc) % reactions(m) % xs(i_temp))
|
||||
if (i_energy >= xs % threshold) then
|
||||
score = score + ((ONE - f) * xs % value(i_energy - &
|
||||
xs % threshold + 1) + f * xs % value(i_energy - &
|
||||
xs % threshold + 2)) * atom_density_ * flux
|
||||
end if
|
||||
end associate
|
||||
else
|
||||
! This block is reached if multipole is turned on and we're
|
||||
! in the resolved range. For (n,gamma), use absorption -
|
||||
! fission. For everything else, assume it's zero.
|
||||
if (score_bin == N_GAMMA) then
|
||||
score = (micro_xs(i_nuc) % absorption - micro_xs(i_nuc) &
|
||||
% fission) * atom_density_ * flux
|
||||
else
|
||||
score = ZERO
|
||||
end if
|
||||
end associate
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
|
|
|
|||
|
|
@ -306,7 +306,7 @@
|
|||
<parameters>-160 -160 -183 160 160 183</parameters>
|
||||
</space>
|
||||
</source>
|
||||
<temperature_multipole>True</temperature_multipole>
|
||||
<temperature_multipole>true</temperature_multipole>
|
||||
</settings>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<tallies>
|
||||
|
|
|
|||
|
|
@ -43,19 +43,13 @@
|
|||
<parameters>-1 -1 -1 1 1 1</parameters>
|
||||
</space>
|
||||
</source>
|
||||
<temperature_multipole>True</temperature_multipole>
|
||||
<temperature_multipole>true</temperature_multipole>
|
||||
<temperature_tolerance>1000</temperature_tolerance>
|
||||
</settings>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<plots>
|
||||
<plot basis="xy" color_by="cell" filename="cellplot" id="1" type="slice">
|
||||
<origin>0 0 0</origin>
|
||||
<width>7 7</width>
|
||||
<pixels>400 400</pixels>
|
||||
</plot>
|
||||
<plot basis="xy" color_by="material" filename="matplot" id="2" type="slice">
|
||||
<origin>0 0 0</origin>
|
||||
<width>7 7</width>
|
||||
<pixels>400 400</pixels>
|
||||
</plot>
|
||||
</plots>
|
||||
<tallies>
|
||||
<tally id="10000">
|
||||
<nuclides>U235 O16 total</nuclides>
|
||||
<scores>total fission (n,gamma) elastic (n,p)</scores>
|
||||
</tally>
|
||||
</tallies>
|
||||
|
|
|
|||
|
|
@ -1,5 +1,36 @@
|
|||
k-combined:
|
||||
1.363786E+00 1.103929E-02
|
||||
tally 1:
|
||||
3.960375E+00
|
||||
3.138356E+00
|
||||
2.875799E+00
|
||||
1.655329E+00
|
||||
5.521904E-01
|
||||
6.105083E-02
|
||||
4.617974E-01
|
||||
4.274130E-02
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
2.254165E+01
|
||||
1.016444E+02
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
6.839351E-04
|
||||
9.359599E-08
|
||||
2.251829E+01
|
||||
1.014341E+02
|
||||
4.903575E-05
|
||||
1.199780E-09
|
||||
3.595002E+02
|
||||
2.586079E+04
|
||||
2.875799E+00
|
||||
1.655329E+00
|
||||
2.174747E+00
|
||||
9.465589E-01
|
||||
3.543564E+02
|
||||
2.512619E+04
|
||||
4.903575E-05
|
||||
1.199780E-09
|
||||
Cell
|
||||
ID = 11
|
||||
Name =
|
||||
|
|
|
|||
|
|
@ -4,94 +4,69 @@ import sys
|
|||
sys.path.insert(0, os.pardir)
|
||||
from testing_harness import TestHarness, PyAPITestHarness
|
||||
import openmc
|
||||
import openmc.model
|
||||
|
||||
|
||||
def make_model():
|
||||
model = openmc.model.Model()
|
||||
|
||||
# Materials
|
||||
moderator = openmc.Material(material_id=1)
|
||||
moderator.set_density('g/cc', 1.0)
|
||||
moderator.add_nuclide('H1', 2.0)
|
||||
moderator.add_nuclide('O16', 1.0)
|
||||
moderator.add_s_alpha_beta('c_H_in_H2O')
|
||||
|
||||
dense_fuel = openmc.Material(material_id=2)
|
||||
dense_fuel.set_density('g/cc', 4.5)
|
||||
dense_fuel.add_nuclide('U235', 1.0)
|
||||
|
||||
model.materials += [moderator, dense_fuel]
|
||||
|
||||
# Geometry
|
||||
c1 = openmc.Cell(cell_id=1, fill=moderator)
|
||||
mod_univ = openmc.Universe(universe_id=1, cells=(c1,))
|
||||
|
||||
r0 = openmc.ZCylinder(R=0.3)
|
||||
c11 = openmc.Cell(cell_id=11, fill=dense_fuel, region=-r0)
|
||||
c11.temperature = [500, 0, 700, 800]
|
||||
c12 = openmc.Cell(cell_id=12, fill=moderator, region=+r0)
|
||||
fuel_univ = openmc.Universe(universe_id=11, 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, fill=lat, region=+x0 & -x1 & +y0 & -y1)
|
||||
model.geometry.root_universe = openmc.Universe(universe_id=0, cells=(c101,))
|
||||
|
||||
# Settings
|
||||
model.settings.batches = 5
|
||||
model.settings.inactive = 0
|
||||
model.settings.particles = 1000
|
||||
model.settings.source = openmc.Source(space=openmc.stats.Box(
|
||||
[-1, -1, -1], [1, 1, 1]))
|
||||
model.settings.temperature = {'tolerance': 1000, 'multipole': True}
|
||||
|
||||
# Tallies
|
||||
tally = openmc.Tally()
|
||||
tally.nuclides = ['U235', 'O16', 'total']
|
||||
tally.scores = ['total', 'fission', '(n,gamma)', 'elastic', '(n,p)']
|
||||
model.tallies.append(tally)
|
||||
|
||||
return model
|
||||
|
||||
|
||||
class MultipoleTestHarness(PyAPITestHarness):
|
||||
def _build_inputs(self):
|
||||
####################
|
||||
# Materials
|
||||
####################
|
||||
|
||||
moderator = openmc.Material(material_id=1)
|
||||
moderator.set_density('g/cc', 1.0)
|
||||
moderator.add_nuclide('H1', 2.0)
|
||||
moderator.add_nuclide('O16', 1.0)
|
||||
moderator.add_s_alpha_beta('c_H_in_H2O')
|
||||
|
||||
dense_fuel = openmc.Material(material_id=2)
|
||||
dense_fuel.set_density('g/cc', 4.5)
|
||||
dense_fuel.add_nuclide('U235', 1.0)
|
||||
|
||||
mats_file = openmc.Materials([moderator, dense_fuel])
|
||||
mats_file.export_to_xml()
|
||||
|
||||
####################
|
||||
# Geometry
|
||||
####################
|
||||
|
||||
c1 = openmc.Cell(cell_id=1, fill=moderator)
|
||||
mod_univ = openmc.Universe(universe_id=1, cells=(c1,))
|
||||
|
||||
r0 = openmc.ZCylinder(R=0.3)
|
||||
c11 = openmc.Cell(cell_id=11, fill=dense_fuel, region=-r0)
|
||||
c11.temperature = [500, 0, 700, 800]
|
||||
c12 = openmc.Cell(cell_id=12, fill=moderator, region=+r0)
|
||||
fuel_univ = openmc.Universe(universe_id=11, 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, fill=lat, region=+x0 & -x1 & +y0 & -y1)
|
||||
root_univ = openmc.Universe(universe_id=0, cells=(c101,))
|
||||
|
||||
geometry = openmc.Geometry(root_univ)
|
||||
geometry.export_to_xml()
|
||||
|
||||
####################
|
||||
# Settings
|
||||
####################
|
||||
|
||||
sets_file = openmc.Settings()
|
||||
sets_file.batches = 5
|
||||
sets_file.inactive = 0
|
||||
sets_file.particles = 1000
|
||||
sets_file.source = openmc.Source(space=openmc.stats.Box(
|
||||
[-1, -1, -1], [1, 1, 1]))
|
||||
sets_file.temperature = {'tolerance': 1000, 'multipole': True}
|
||||
sets_file.export_to_xml()
|
||||
|
||||
####################
|
||||
# Plots
|
||||
####################
|
||||
|
||||
plot1 = openmc.Plot(plot_id=1)
|
||||
plot1.basis = 'xy'
|
||||
plot1.color_by = 'cell'
|
||||
plot1.filename = 'cellplot'
|
||||
plot1.origin = (0, 0, 0)
|
||||
plot1.width = (7, 7)
|
||||
plot1.pixels = (400, 400)
|
||||
|
||||
plot2 = openmc.Plot(plot_id=2)
|
||||
plot2.basis = 'xy'
|
||||
plot2.color_by = 'material'
|
||||
plot2.filename = 'matplot'
|
||||
plot2.origin = (0, 0, 0)
|
||||
plot2.width = (7, 7)
|
||||
plot2.pixels = (400, 400)
|
||||
|
||||
plots = openmc.Plots([plot1, plot2])
|
||||
plots.export_to_xml()
|
||||
|
||||
def execute_test(self):
|
||||
if not 'OPENMC_MULTIPOLE_LIBRARY' in os.environ:
|
||||
raise RuntimeError("The 'OPENMC_MULTIPOLE_LIBRARY' environment "
|
||||
|
|
@ -107,5 +82,6 @@ class MultipoleTestHarness(PyAPITestHarness):
|
|||
|
||||
|
||||
if __name__ == '__main__':
|
||||
harness = MultipoleTestHarness('statepoint.5.h5')
|
||||
model = make_model()
|
||||
harness = MultipoleTestHarness('statepoint.5.h5', model)
|
||||
harness.main()
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue