diff --git a/.gitignore b/.gitignore index e758f64e43..5749870dfb 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,9 @@ # OpenMC executable src/openmc +# Inputs generated from Python API +examples/python/**/*.xml + # emacs backups *~ @@ -41,3 +44,6 @@ results_error.dat # Data downloaded from NNDC data/nndc + +.idea +.idea/* \ No newline at end of file diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index bf6f93a3c7..d98804c867 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -174,11 +174,15 @@ should be performed. It has the following attributes/sub-elements: ------------------------- The ```` element determines the treatment of the energy grid during -a simulation. The valid options are "nuclide" and "logarithm". Setting this -element to "nuclide" will cause OpenMC to use a nuclide's energy grid when -determining what points to interpolate between for determining cross sections -(i.e. non-unionized energy grid). Setting this element to "logarithm" causes -OpenMC to use a logarithmic mapping technique described in LA-UR-14-24530_. +a simulation. The valid options are "nuclide", "logarithm", and +"material-union". Setting this element to "nuclide" will cause OpenMC to use a +nuclide's energy grid when determining what points to interpolate between for +determining cross sections (i.e. non-unionized energy grid). Setting this +element to "logarithm" causes OpenMC to use a logarithmic mapping technique +described in LA-UR-14-24530_. Setting this element to "material-union" will +cause OpenMC to create energy grids that are unionized material-by-material and +use these grids when determining the energy-cross section pairs to interpolate +cross section values between. *Default*: logarithm diff --git a/examples/python/basic/build-xml.py b/examples/python/basic/build-xml.py index 58825ee494..312a05a8db 100644 --- a/examples/python/basic/build-xml.py +++ b/examples/python/basic/build-xml.py @@ -33,7 +33,7 @@ fuel.add_nuclide(u235, 1.) # Instantiate a MaterialsFile, register all Materials, and export to XML materials_file = openmc.MaterialsFile() -materials_file.set_default_xs('71c') +materials_file.default_xs = '71c' materials_file.add_materials([moderator, fuel]) materials_file.export_to_xml() @@ -46,7 +46,7 @@ materials_file.export_to_xml() surf1 = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=7, name='surf 1') surf2 = openmc.ZCylinder(surface_id=2, x0=0, y0=0, R=9, name='surf 2') surf3 = openmc.ZCylinder(surface_id=3, x0=0, y0=0, R=11, name='surf 3') -surf3.set_boundary_type('vacuum') +surf3.boundary_type = 'vacuum' # Instantiate Cells cell1 = openmc.Cell(cell_id=1, name='cell 1') @@ -62,14 +62,14 @@ cell4.add_surface(surface=surf2, halfspace=+1) cell4.add_surface(surface=surf3, halfspace=-1) # Register Materials with Cells -cell2.set_fill(fuel) -cell3.set_fill(moderator) -cell4.set_fill(moderator) +cell2.fill = fuel +cell3.fill = moderator +cell4.fill = moderator # Instantiate Universes universe1 = openmc.Universe(universe_id=37) root = openmc.Universe(universe_id=0, name='root universe') -cell1.set_fill(universe1) +cell1.fill = universe1 # Register Cells with Universes universe1.add_cells([cell2, cell3]) @@ -77,11 +77,11 @@ root.add_cells([cell1, cell4]) # Instantiate a Geometry and register the root Universe geometry = openmc.Geometry() -geometry.set_root_universe(root) +geometry.root_universe = root # Instantiate a GeometryFile, register Geometry, and export to XML geometry_file = openmc.GeometryFile() -geometry_file.set_geometry(geometry) +geometry_file.geometry = geometry geometry_file.export_to_xml() @@ -91,9 +91,9 @@ geometry_file.export_to_xml() # Instantiate a SettingsFile, set all runtime parameters, and export to XML settings_file = openmc.SettingsFile() -settings_file.set_batches(batches) -settings_file.set_inactive(inactive) -settings_file.set_particles(particles) +settings_file.batches = batches +settings_file.inactive = inactive +settings_file.particles = particles settings_file.set_source_space('box', [-4, -4, -4, 4, 4, 4]) settings_file.export_to_xml() diff --git a/examples/python/lattice/hexagonal/build-xml.py b/examples/python/lattice/hexagonal/build-xml.py index 6e0a58bc1b..9456b899b8 100644 --- a/examples/python/lattice/hexagonal/build-xml.py +++ b/examples/python/lattice/hexagonal/build-xml.py @@ -38,7 +38,7 @@ iron.add_nuclide(fe56, 1.) # Instantiate a MaterialsFile, register all Materials, and export to XML materials_file = openmc.MaterialsFile() -materials_file.set_default_xs('71c') +materials_file.default_xs = '71c' materials_file.add_materials([moderator, fuel, iron]) materials_file.export_to_xml() @@ -54,10 +54,10 @@ bottom = openmc.YPlane(surface_id=3, y0=-4, name='bottom') top = openmc.YPlane(surface_id=4, y0=4, name='top') fuel_surf = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4) -left.set_boundary_type('vacuum') -right.set_boundary_type('vacuum') -top.set_boundary_type('vacuum') -bottom.set_boundary_type('vacuum') +left.boundary_type = 'vacuum' +right.boundary_type = 'vacuum' +top.boundary_type = 'vacuum' +bottom.boundary_type = 'vacuum' # Instantiate Cells cell1 = openmc.Cell(cell_id=1, name='Cell 1') @@ -78,11 +78,11 @@ cell5.add_surface(fuel_surf, halfspace=-1) cell6.add_surface(fuel_surf, halfspace=+1) # Register Materials with Cells -cell2.set_fill(fuel) -cell3.set_fill(moderator) -cell4.set_fill(moderator) -cell5.set_fill(iron) -cell6.set_fill(moderator) +cell2.fill = fuel +cell3.fill = moderator +cell4.fill = moderator +cell5.fill = iron +cell6.fill = moderator # Instantiate Universe univ1 = openmc.Universe(universe_id=1) @@ -98,24 +98,24 @@ root.add_cell(cell1) # Instantiate a Lattice lattice = openmc.HexLattice(lattice_id=5) -lattice.set_center([0., 0., 0.]) -lattice.set_pitch([1., 2.]) -lattice.set_universes([ - [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ], - [ [univ2] + [univ1]*11, [univ2] + [univ1]*5, [univ1] ], - [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ]]) -lattice.set_outer(univ2) +lattice.center = [0., 0., 0.] +lattice.pitch = [1., 2.] +lattice.universes = \ + [ [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ], + [ [univ2] + [univ1]*11, [univ2] + [univ1]*5, [univ1] ], + [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ] ] +lattice.outer = univ2 # Fill Cell with the Lattice -cell1.set_fill(lattice) +cell1.fill = lattice # Instantiate a Geometry and register the root Universe geometry = openmc.Geometry() -geometry.set_root_universe(root) +geometry.root_universe = root # Instantiate a GeometryFile, register Geometry, and export to XML geometry_file = openmc.GeometryFile() -geometry_file.set_geometry(geometry) +geometry_file.geometry = geometry geometry_file.export_to_xml() @@ -125,9 +125,9 @@ geometry_file.export_to_xml() # Instantiate a SettingsFile, set all runtime parameters, and export to XML settings_file = openmc.SettingsFile() -settings_file.set_batches(batches) -settings_file.set_inactive(inactive) -settings_file.set_particles(particles) +settings_file.batches = batches +settings_file.inactive = inactive +settings_file.particles = particles settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1]) settings_file.export_to_xml() @@ -137,19 +137,19 @@ settings_file.export_to_xml() ############################################################################### plot_xy = openmc.Plot(plot_id=1) -plot_xy.set_filename('plot_xy') -plot_xy.set_origin([0, 0, 0]) -plot_xy.set_width([6, 6]) -plot_xy.set_pixels([400, 400]) -plot_xy.set_color('mat') +plot_xy.filename = 'plot_xy' +plot_xy.origin = [0, 0, 0] +plot_xy.width = [6, 6] +plot_xy.pixels = [400, 400] +plot_xy.color = 'mat' plot_yz = openmc.Plot(plot_id=2) -plot_yz.set_filename('plot_yz') -plot_yz.set_basis('yz') -plot_yz.set_origin([0, 0, 0]) -plot_yz.set_width([8, 8]) -plot_yz.set_pixels([400, 400]) -plot_yz.set_color('mat') +plot_yz.filename = 'plot_yz' +plot_yz.basis = 'yz' +plot_yz.origin = [0, 0, 0] +plot_yz.width = [8, 8] +plot_yz.pixels = [400, 400] +plot_yz.color = 'mat' # Instantiate a PlotsFile, add Plot, and export to XML plot_file = openmc.PlotsFile() diff --git a/examples/python/lattice/nested/build-xml.py b/examples/python/lattice/nested/build-xml.py index df59bf0870..ce67665429 100644 --- a/examples/python/lattice/nested/build-xml.py +++ b/examples/python/lattice/nested/build-xml.py @@ -33,7 +33,7 @@ moderator.add_s_alpha_beta('HH2O', '71t') # Instantiate a MaterialsFile, register all Materials, and export to XML materials_file = openmc.MaterialsFile() -materials_file.set_default_xs('71c') +materials_file.default_xs = '71c' materials_file.add_materials([moderator, fuel]) materials_file.export_to_xml() @@ -51,10 +51,10 @@ fuel1 = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4) fuel2 = openmc.ZCylinder(surface_id=6, x0=0, y0=0, R=0.3) fuel3 = openmc.ZCylinder(surface_id=7, x0=0, y0=0, R=0.2) -left.set_boundary_type('vacuum') -right.set_boundary_type('vacuum') -top.set_boundary_type('vacuum') -bottom.set_boundary_type('vacuum') +left.boundary_type = 'vacuum' +right.boundary_type = 'vacuum' +top.boundary_type = 'vacuum' +bottom.boundary_type = 'vacuum' # Instantiate Cells cell1 = openmc.Cell(cell_id=1, name='Cell 1') @@ -83,12 +83,12 @@ cell7.add_surface(fuel3, halfspace=-1) cell8.add_surface(fuel3, halfspace=+1) # Register Materials with Cells -cell3.set_fill(fuel) -cell4.set_fill(moderator) -cell5.set_fill(fuel) -cell6.set_fill(moderator) -cell7.set_fill(fuel) -cell8.set_fill(moderator) +cell3.fill = fuel +cell4.fill = moderator +cell5.fill = fuel +cell6.fill = moderator +cell7.fill = fuel +cell8.fill = moderator # Instantiate Universe univ1 = openmc.Universe(universe_id=1) @@ -106,30 +106,30 @@ univ4.add_cell(cell2) # Instantiate nested Lattices lattice1 = openmc.RectLattice(lattice_id=4, name='4x4 assembly') -lattice1.set_dimension([2, 2]) -lattice1.set_lower_left([-1., -1.]) -lattice1.set_pitch([1., 1.]) -lattice1.set_universes([[univ1, univ2], - [univ2, univ3]]) +lattice1.dimension = [2, 2] +lattice1.lower_left = [-1., -1.] +lattice1.pitch = [1., 1.] +lattice1.universes = [[univ1, univ2], + [univ2, univ3]] lattice2 = openmc.RectLattice(lattice_id=6, name='4x4 core') -lattice2.set_dimension([2, 2]) -lattice2.set_lower_left([-2., -2.]) -lattice2.set_pitch([2., 2.]) -lattice2.set_universes([[univ4, univ4], - [univ4, univ4]]) +lattice2.dimension = [2, 2] +lattice2.lower_left = [-2., -2.] +lattice2.pitch = [2., 2.] +lattice2.universes = [[univ4, univ4], + [univ4, univ4]] # Fill Cell with the Lattice -cell1.set_fill(lattice2) -cell2.set_fill(lattice1) +cell1.fill = lattice2 +cell2.fill = lattice1 # Instantiate a Geometry and register the root Universe geometry = openmc.Geometry() -geometry.set_root_universe(root) +geometry.root_universe = root # Instantiate a GeometryFile, register Geometry, and export to XML geometry_file = openmc.GeometryFile() -geometry_file.set_geometry(geometry) +geometry_file.geometry = geometry geometry_file.export_to_xml() @@ -139,9 +139,9 @@ geometry_file.export_to_xml() # Instantiate a SettingsFile, set all runtime parameters, and export to XML settings_file = openmc.SettingsFile() -settings_file.set_batches(batches) -settings_file.set_inactive(inactive) -settings_file.set_particles(particles) +settings_file.batches = batches +settings_file.inactive = inactive +settings_file.particles = particles settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1]) settings_file.export_to_xml() @@ -151,10 +151,10 @@ settings_file.export_to_xml() ############################################################################### plot = openmc.Plot(plot_id=1) -plot.set_origin([0, 0, 0]) -plot.set_width([4, 4]) -plot.set_pixels([400, 400]) -plot.set_color('mat') +plot.origin = [0, 0, 0] +plot.width = [4, 4] +plot.pixels = [400, 400] +plot.color = 'mat' # Instantiate a PlotsFile, add Plot, and export to XML plot_file = openmc.PlotsFile() @@ -168,14 +168,14 @@ plot_file.export_to_xml() # Instantiate a tally mesh mesh = openmc.Mesh(mesh_id=1) -mesh.set_type('rectangular') -mesh.set_dimension([4, 4]) -mesh.set_lower_left([-2, -2]) -mesh.set_width([1, 1]) +mesh.type = 'rectangular' +mesh.dimension = [4, 4] +mesh.lower_left = [-2, -2] +mesh.width = [1, 1] # Instantiate tally Filter mesh_filter = openmc.Filter() -mesh_filter.set_mesh(mesh) +mesh_filter.mesh = mesh # Instantiate the Tally tally = openmc.Tally(tally_id=1) diff --git a/examples/python/lattice/simple/build-xml.py b/examples/python/lattice/simple/build-xml.py index 83c403abd6..98ee1fa669 100644 --- a/examples/python/lattice/simple/build-xml.py +++ b/examples/python/lattice/simple/build-xml.py @@ -33,7 +33,7 @@ moderator.add_s_alpha_beta('HH2O', '71t') # Instantiate a MaterialsFile, register all Materials, and export to XML materials_file = openmc.MaterialsFile() -materials_file.set_default_xs('71c') +materials_file.default_xs = '71c' materials_file.add_materials([moderator, fuel]) materials_file.export_to_xml() @@ -51,10 +51,10 @@ fuel1 = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4) fuel2 = openmc.ZCylinder(surface_id=6, x0=0, y0=0, R=0.3) fuel3 = openmc.ZCylinder(surface_id=7, x0=0, y0=0, R=0.2) -left.set_boundary_type('vacuum') -right.set_boundary_type('vacuum') -top.set_boundary_type('vacuum') -bottom.set_boundary_type('vacuum') +left.boundary_type = 'vacuum' +right.boundary_type = 'vacuum' +top.boundary_type = 'vacuum' +bottom.boundary_type = 'vacuum' # Instantiate Cells cell1 = openmc.Cell(cell_id=1, name='Cell 1') @@ -78,12 +78,12 @@ cell6.add_surface(fuel3, halfspace=-1) cell7.add_surface(fuel3, halfspace=+1) # Register Materials with Cells -cell2.set_fill(fuel) -cell3.set_fill(moderator) -cell4.set_fill(fuel) -cell5.set_fill(moderator) -cell6.set_fill(fuel) -cell7.set_fill(moderator) +cell2.fill = fuel +cell3.fill = moderator +cell4.fill = fuel +cell5.fill = moderator +cell6.fill = fuel +cell7.fill = moderator # Instantiate Universe univ1 = openmc.Universe(universe_id=1) @@ -99,24 +99,24 @@ root.add_cell(cell1) # Instantiate a Lattice lattice = openmc.RectLattice(lattice_id=5) -lattice.set_dimension([4, 4]) -lattice.set_lower_left([-2., -2.]) -lattice.set_pitch([1., 1.]) -lattice.set_universes([[univ1, univ2, univ1, univ2], - [univ2, univ3, univ2, univ3], - [univ1, univ2, univ1, univ2], - [univ2, univ3, univ2, univ3]]) +lattice.dimension = [4, 4] +lattice.lower_left = [-2., -2.] +lattice.pitch = [1., 1.] +lattice.universes = [[univ1, univ2, univ1, univ2], + [univ2, univ3, univ2, univ3], + [univ1, univ2, univ1, univ2], + [univ2, univ3, univ2, univ3]] # Fill Cell with the Lattice -cell1.set_fill(lattice) +cell1.fill = lattice # Instantiate a Geometry and register the root Universe geometry = openmc.Geometry() -geometry.set_root_universe(root) +geometry.root_universe = root # Instantiate a GeometryFile, register Geometry, and export to XML geometry_file = openmc.GeometryFile() -geometry_file.set_geometry(geometry) +geometry_file.geometry = geometry geometry_file.export_to_xml() @@ -126,9 +126,9 @@ geometry_file.export_to_xml() # Instantiate a SettingsFile, set all runtime parameters, and export to XML settings_file = openmc.SettingsFile() -settings_file.set_batches(batches) -settings_file.set_inactive(inactive) -settings_file.set_particles(particles) +settings_file.batches = batches +settings_file.inactive = inactive +settings_file.particles = particles settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1]) settings_file.export_to_xml() @@ -138,10 +138,10 @@ settings_file.export_to_xml() ############################################################################### plot = openmc.Plot(plot_id=1) -plot.set_origin([0, 0, 0]) -plot.set_width([4, 4]) -plot.set_pixels([400, 400]) -plot.set_color('mat') +plot.origin = [0, 0, 0] +plot.width = [4, 4] +plot.pixels = [400, 400] +plot.color = 'mat' # Instantiate a PlotsFile, add Plot, and export to XML plot_file = openmc.PlotsFile() @@ -155,14 +155,14 @@ plot_file.export_to_xml() # Instantiate a tally mesh mesh = openmc.Mesh(mesh_id=1) -mesh.set_type('rectangular') -mesh.set_dimension([4, 4]) -mesh.set_lower_left([-2, -2]) -mesh.set_width([1, 1]) +mesh.type = 'rectangular' +mesh.dimension = [4, 4] +mesh.lower_left = [-2, -2] +mesh.width = [1, 1] # Instantiate tally Filter mesh_filter = openmc.Filter() -mesh_filter.set_mesh(mesh) +mesh_filter.mesh = mesh # Instantiate the Tally tally = openmc.Tally(tally_id=1) diff --git a/examples/python/pincell/build-xml.py b/examples/python/pincell/build-xml.py index 4da66fe459..c79fec7fe9 100644 --- a/examples/python/pincell/build-xml.py +++ b/examples/python/pincell/build-xml.py @@ -103,7 +103,7 @@ borated_water.add_s_alpha_beta('HH2O', '71t') # Instantiate a MaterialsFile, register all Materials, and export to XML materials_file = openmc.MaterialsFile() -materials_file.set_default_xs('71c') +materials_file.default_xs = '71c' materials_file.add_materials([uo2, helium, zircaloy, borated_water]) materials_file.export_to_xml() @@ -121,10 +121,10 @@ right = openmc.XPlane(surface_id=5, x0=0.62992, name='right') bottom = openmc.YPlane(surface_id=6, y0=-0.62992, name='bottom') top = openmc.YPlane(surface_id=7, y0=0.62992, name='top') -left.set_boundary_type('reflective') -right.set_boundary_type('reflective') -top.set_boundary_type('reflective') -bottom.set_boundary_type('reflective') +left.boundary_type = 'reflective' +right.boundary_type = 'reflective' +top.boundary_type = 'reflective' +bottom.boundary_type = 'reflective' # Instantiate Cells fuel = openmc.Cell(cell_id=1, name='cell 1') @@ -145,10 +145,10 @@ water.add_surface(bottom, halfspace=+1) water.add_surface(top, halfspace=-1) # Register Materials with Cells -fuel.set_fill(uo2) -gap.set_fill(helium) -clad.set_fill(zircaloy) -water.set_fill(borated_water) +fuel.fill = uo2 +gap.fill = helium +clad.fill = zircaloy +water.fill = borated_water # Instantiate Universe root = openmc.Universe(universe_id=0, name='root universe') @@ -158,11 +158,11 @@ root.add_cells([fuel, gap, clad, water]) # Instantiate a Geometry and register the root Universe geometry = openmc.Geometry() -geometry.set_root_universe(root) +geometry.root_universe = root # Instantiate a GeometryFile, register Geometry, and export to XML geometry_file = openmc.GeometryFile() -geometry_file.set_geometry(geometry) +geometry_file.geometry = geometry geometry_file.export_to_xml() @@ -172,14 +172,14 @@ geometry_file.export_to_xml() # Instantiate a SettingsFile, set all runtime parameters, and export to XML settings_file = openmc.SettingsFile() -settings_file.set_batches(batches) -settings_file.set_inactive(inactive) -settings_file.set_particles(particles) +settings_file.batches = batches +settings_file.inactive = inactive +settings_file.particles = particles settings_file.set_source_space('box', [-0.62992, -0.62992, -1, \ 0.62992, 0.62992, 1]) -settings_file.set_entropy_lower_left([-0.39218, -0.39218, -1.e50]) -settings_file.set_entropy_upper_right([0.39218, 0.39218, 1.e50]) -settings_file.set_entropy_dimension([10, 10, 1]) +settings_file.entropy_lower_left = [-0.39218, -0.39218, -1.e50] +settings_file.entropy_upper_right = [0.39218, 0.39218, 1.e50] +settings_file.entropy_dimension = [10, 10, 1] settings_file.export_to_xml() @@ -189,15 +189,15 @@ settings_file.export_to_xml() # Instantiate a tally mesh mesh = openmc.Mesh(mesh_id=1) -mesh.set_type('rectangular') -mesh.set_dimension([100, 100, 1]) -mesh.set_lower_left([-0.62992, -0.62992, -1.e50]) -mesh.set_upper_right([0.62992, 0.62992, 1.e50]) +mesh.type = 'rectangular' +mesh.dimension = [100, 100, 1] +mesh.lower_left = [-0.62992, -0.62992, -1.e50] +mesh.upper_right = [0.62992, 0.62992, 1.e50] # Instantiate some tally Filters energy_filter = openmc.Filter(type='energy', bins=[0., 4.e-6, 20.]) mesh_filter = openmc.Filter() -mesh_filter.set_mesh(mesh) +mesh_filter.mesh = mesh # Instantiate the Tally tally = openmc.Tally(tally_id=1) diff --git a/examples/python/reflective/build-xml.py b/examples/python/reflective/build-xml.py index 4e09eeb972..9a182903e1 100644 --- a/examples/python/reflective/build-xml.py +++ b/examples/python/reflective/build-xml.py @@ -25,7 +25,7 @@ fuel.add_nuclide(u235, 1.) # Instantiate a MaterialsFile, register Material, and export to XML materials_file = openmc.MaterialsFile() -materials_file.set_default_xs('71c') +materials_file.default_xs = '71c' materials_file.add_material(fuel) materials_file.export_to_xml() @@ -42,12 +42,12 @@ surf4 = openmc.YPlane(surface_id=4, y0=+1, name='surf 4') surf5 = openmc.ZPlane(surface_id=5, z0=-1, name='surf 5') surf6 = openmc.ZPlane(surface_id=6, z0=+1, name='surf 6') -surf1.set_boundary_type('vacuum') -surf2.set_boundary_type('vacuum') -surf3.set_boundary_type('reflective') -surf4.set_boundary_type('reflective') -surf5.set_boundary_type('reflective') -surf6.set_boundary_type('reflective') +surf1.boundary_type = 'vacuum' +surf2.boundary_type = 'vacuum' +surf3.boundary_type = 'reflective' +surf4.boundary_type = 'reflective' +surf5.boundary_type = 'reflective' +surf6.boundary_type = 'reflective' # Instantiate Cell cell = openmc.Cell(cell_id=1, name='cell 1') @@ -61,7 +61,7 @@ cell.add_surface(surface=surf5, halfspace=+1) cell.add_surface(surface=surf6, halfspace=-1) # Register Material with Cell -cell.set_fill(fuel) +cell.fill = fuel # Instantiate Universes root = openmc.Universe(universe_id=0, name='root universe') @@ -71,11 +71,11 @@ root.add_cell(cell) # Instantiate a Geometry and register the root Universe geometry = openmc.Geometry() -geometry.set_root_universe(root) +geometry.root_universe = root # Instantiate a GeometryFile, register Geometry, and export to XML geometry_file = openmc.GeometryFile() -geometry_file.set_geometry(geometry) +geometry_file.geometry = geometry geometry_file.export_to_xml() @@ -85,8 +85,8 @@ geometry_file.export_to_xml() # Instantiate a SettingsFile, set all runtime parameters, and export to XML settings_file = openmc.SettingsFile() -settings_file.set_batches(batches) -settings_file.set_inactive(inactive) -settings_file.set_particles(particles) +settings_file.batches = batches +settings_file.inactive = inactive +settings_file.particles = particles settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1]) settings_file.export_to_xml() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b3e1a88779..89a4a60e12 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -37,7 +37,7 @@ endif() set(MPI_ENABLED FALSE) set(HDF5_ENABLED FALSE) -if($ENV{FC} MATCHES "mpi.*") +if($ENV{FC} MATCHES "mpi[^/]*$") message("-- Detected MPI wrapper: $ENV{FC}") add_definitions(-DMPI) set(MPI_ENABLED TRUE) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 1f34a2311b..7753f4734a 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -103,7 +103,7 @@ module ace_header ! Energy grid information integer :: n_grid ! # of nuclide grid points - integer, allocatable :: grid_index(:) ! union grid pointers / log grid mapping + integer, allocatable :: grid_index(:) ! log grid mapping indices real(8), allocatable :: energy(:) ! energy values corresponding to xs ! Microscopic cross sections @@ -372,9 +372,6 @@ module ace_header integer :: i ! Loop counter - if (allocated(this % grid_index)) & - deallocate(this % grid_index) - if (allocated(this % energy)) & deallocate(this % energy, this % total, this % elastic, & & this % fission, this % nu_fission, this % absorption) diff --git a/src/constants.F90 b/src/constants.F90 index bac301e36a..85976c00ac 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -368,8 +368,9 @@ module constants ! Energy grid methods integer, parameter :: & - GRID_NUCLIDE = 1, & ! non-unionized energy grid - GRID_LOGARITHM = 2 ! logarithmic mapping + GRID_NUCLIDE = 1, & ! unique energy grid for each nuclide + GRID_MAT_UNION = 2, & ! material union grids with pointers + GRID_LOGARITHM = 3 ! lethargy mapping ! Running modes integer, parameter :: & diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1672405e8d..dce4d75f32 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -15,6 +15,9 @@ module cross_section implicit none save + integer :: union_grid_index +!$omp threadprivate(union_grid_index) + contains !=============================================================================== @@ -36,11 +39,11 @@ contains !$omp threadprivate(mat) ! Set all material macroscopic cross sections to zero - material_xs % total = ZERO - material_xs % elastic = ZERO - material_xs % absorption = ZERO - material_xs % fission = ZERO - material_xs % nu_fission = ZERO + material_xs % total = ZERO + material_xs % elastic = ZERO + material_xs % absorption = ZERO + material_xs % fission = ZERO + material_xs % nu_fission = ZERO material_xs % kappa_fission = ZERO ! Exit subroutine if material is void @@ -48,6 +51,10 @@ contains mat => materials(p % material) + ! Find energy index on global or material unionized grid + if (grid_method == GRID_MAT_UNION) & + call find_energy_index(p % E, p % material) + ! Determine if this material has S(a,b) tables check_sab = (mat % n_sab > 0) @@ -88,9 +95,9 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i) end if ! ======================================================================== @@ -131,24 +138,32 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat) 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 - - integer :: i_grid ! index on nuclide energy grid - integer :: i_low, i_high ! bounding indices from logarithmic mapping - integer :: u ! index into logarithmic mapping array + integer, intent(in) :: i_mat ! index into materials array + integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material + integer :: i_grid ! index on nuclide energy grid + integer :: i_low ! lower logarithmic mapping index + integer :: i_high ! upper logarithmic mapping index + integer :: u ! index into logarithmic mapping array + real(8), intent(in) :: E ! energy real(8) :: f ! interp factor on nuclide energy grid - type(Nuclide), pointer, save :: nuc => null() -!$omp threadprivate(nuc) + type(Nuclide), pointer, save :: nuc => null() + type(Material), pointer, save :: mat => null() +!$omp threadprivate(nuc, mat) - ! Set pointer to nuclide + ! Set pointer to nuclide and material nuc => nuclides(i_nuclide) + mat => materials(i_mat) ! Determine index on nuclide energy grid select case (grid_method) + case (GRID_MAT_UNION) + + i_grid = mat % nuclide_grid_index(i_nuc_mat, union_grid_index) + 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 @@ -173,7 +188,7 @@ contains ! Perform binary search on the nuclide energy grid in order to determine ! which points to interpolate between - if (E < nuc % energy(1)) then + if (E <= nuc % energy(1)) then i_grid = 1 elseif (E > nuc % energy(nuc % n_grid)) then i_grid = nuc % n_grid - 1 @@ -506,6 +521,32 @@ contains end subroutine calculate_urr_xs +!=============================================================================== +! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain +! energy +!=============================================================================== + + subroutine find_energy_index(E, i_mat) + + real(8), intent(in) :: E ! energy of particle + integer, intent(in) :: i_mat ! material index + type(Material), pointer, save :: mat => null() ! pointer to current material +!$omp threadprivate(mat) + + mat => materials(i_mat) + + ! if the energy is outside of energy grid range, set to first or last + ! index. Otherwise, do a binary search through the union energy grid. + if (E <= mat % e_grid(1)) then + union_grid_index = 1 + elseif (E > mat % e_grid(mat % n_grid)) then + union_grid_index = mat % n_grid - 1 + else + union_grid_index = binary_search(mat % e_grid, mat % n_grid, E) + end if + + end subroutine find_energy_index + !=============================================================================== ! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section ! for a given nuclide at the trial relative energy used in resonance scattering diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 234a856c1e..676875a5b4 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -1,6 +1,9 @@ module energy_grid use global + use list_header, only: ListReal + use output, only: write_message + implicit none @@ -10,6 +13,53 @@ module energy_grid contains +!=============================================================================== +! UNIONIZED_GRID creates a unionized energy grid, for the entire problem or for +! each material, composed of the grids from each nuclide in the entire problem, +! or each material, respectively. Right now, the grid for each nuclide is added +! into a linked list one at a time with an effective insertion sort. Could be +! done with a hash for all energy points and then a quicksort at the end (what +! hash function to use?) +!=============================================================================== + + subroutine unionized_grid() + + integer :: i ! index in nuclides array + integer :: j ! index in materials array + type(ListReal), pointer, save :: list => null() + type(Nuclide), pointer, save :: nuc => null() + type(Material), pointer, save :: mat => null() + !$omp threadprivate(list, nuc, mat) + + call write_message("Creating unionized energy grid...", 5) + + ! add grid points for each nuclide in the material + do j = 1, n_materials + mat => materials(j) + do i = 1, mat % n_nuclides + nuc => nuclides(mat % nuclide(i)) + call add_grid_points(list, nuc % energy) + end do + + ! set size of unionized material energy grid + mat % n_grid = list % size() + + ! create allocated array from linked list + allocate(mat % e_grid(mat % n_grid)) + do i = 1, mat % n_grid + mat % e_grid(i) = list % get_item(i) + end do + + ! delete linked list and dictionary + call list % clear() + deallocate(list) + end do + + ! Set pointers to unionized energy grid for each nuclide + call grid_pointers() + + end subroutine unionized_grid + !=============================================================================== ! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding ! indices on a nuclide energy grid @@ -59,4 +109,109 @@ contains end subroutine logarithmic_grid +!=============================================================================== +! ADD_GRID_POINTS adds energy points from the 'energy' array into a linked list +! of points already stored from previous arrays. +!=============================================================================== + + subroutine add_grid_points(list, energy) + + type(ListReal), pointer :: list + real(8), intent(in) :: energy(:) + + integer :: i ! index in energy array + integer :: n ! size of energy array + integer :: current ! current index + real(8) :: E ! actual energy value + + i = 1 + n = size(energy) + + ! If the original list is empty, we need to allocate the first element and + ! store first energy point + if (.not. associated(list)) then + allocate(list) + do i = 1, n + call list % append(energy(i)) + end do + return + end if + + ! Set current index to beginning of the list + current = 1 + + do while (i <= n) + E = energy(i) + + ! If we've reached the end of the grid energy list, add the remaining + ! energy points to the end + if (current > list % size()) then + ! Finish remaining energies + do while (i <= n) + call list % append(energy(i)) + i = i + 1 + end do + exit + end if + + if (E < list % get_item(current)) then + + ! Insert new energy in this position + call list % insert(current, E) + + ! Advance index in linked list and in new energy grid + i = i + 1 + current = current + 1 + + elseif (E == list % get_item(current)) then + ! Found the exact same energy, no need to store duplicates so just + ! skip and move to next index + i = i + 1 + current = current + 1 + else + current = current + 1 + end if + + end do + + end subroutine add_grid_points + +!=============================================================================== +! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link +! each point on the nuclide energy grid to one on a unionized energy grid +!=============================================================================== + + subroutine grid_pointers() + + integer :: i ! loop index for nuclides + integer :: j ! loop index for nuclide energy grid + integer :: k ! loop index for materials + integer :: index_e ! index on union energy grid + real(8) :: union_energy ! energy on union grid + real(8) :: energy ! energy on nuclide grid + type(Nuclide), pointer :: nuc => null() + type(Material), pointer :: mat => null() + + do k = 1, n_materials + mat => materials(k) + allocate(mat % nuclide_grid_index(mat % n_nuclides, mat % n_grid)) + do i = 1, mat % n_nuclides + nuc => nuclides(mat % nuclide(i)) + + index_e = 1 + energy = nuc % energy(index_e) + + do j = 1, mat % n_grid + union_energy = mat % e_grid(j) + if (union_energy >= energy .and. index_e < nuc % n_grid) then + index_e = index_e + 1 + energy = nuc % energy(index_e) + end if + mat % nuclide_grid_index(i,j) = index_e - 1 + end do + end do + end do + + end subroutine grid_pointers + end module energy_grid diff --git a/src/global.F90 b/src/global.F90 index 004c43746e..06be1c992f 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -218,6 +218,7 @@ module global type(Timer) :: time_total ! timer for total run type(Timer) :: time_initialize ! timer for initialization type(Timer) :: time_read_xs ! timer for reading cross sections + type(Timer) :: time_unionize ! timer for material xs-energy grid union type(Timer) :: time_bank ! timer for fission bank synchronization type(Timer) :: time_bank_sample ! timer for fission bank sampling type(Timer) :: time_bank_sendrecv ! timer for fission bank SEND/RECV diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 125fdaa0f0..2252c3e763 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -38,7 +38,6 @@ module hdf5_interface module procedure hdf5_write_integer_4Darray module procedure hdf5_write_long module procedure hdf5_write_string - module procedure hdf5_write_string_1Darray #ifdef MPI module procedure hdf5_write_double_parallel module procedure hdf5_write_double_1Darray_parallel @@ -52,7 +51,6 @@ module hdf5_interface module procedure hdf5_write_integer_4Darray_parallel module procedure hdf5_write_long_parallel module procedure hdf5_write_string_parallel -! module procedure hdf5_write_string_1Darray_parallel #endif end interface hdf5_write_data @@ -70,7 +68,6 @@ module hdf5_interface module procedure hdf5_read_integer_4Darray module procedure hdf5_read_long module procedure hdf5_read_string - module procedure hdf5_read_string_1Darray #ifdef MPI module procedure hdf5_read_double_parallel module procedure hdf5_read_double_1Darray_parallel @@ -84,7 +81,6 @@ module hdf5_interface module procedure hdf5_read_integer_4Darray_parallel module procedure hdf5_read_long_parallel module procedure hdf5_read_string_parallel -! module procedure hdf5_read_string_1Darray_parallel #endif end interface hdf5_read_data @@ -210,7 +206,7 @@ contains logical :: status ! does the group exist ! Check if group exists - call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err) + call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err) ! Either create or open group if (status) then @@ -261,7 +257,7 @@ contains integer(HID_T), intent(in) :: group ! name of group character(*), intent(in) :: name ! name of data - integer, intent(inout) :: buffer ! read data to here + integer, intent(inout) :: buffer ! read data to here integer :: buffer_copy(1) ! need an array for read @@ -465,7 +461,7 @@ contains integer(HID_T), intent(in) :: group ! name of group character(*), intent(in) :: name ! name of data - real(8), intent(inout) :: buffer ! read data to here + real(8), intent(inout) :: buffer ! read data to here real(8) :: buffer_copy(1) ! need an array for read @@ -730,7 +726,7 @@ contains call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err) ! Set up dimesnions of string to write - dims2 = (/length, 1/) ! full array of strings to write + dims2 = (/length, 1/) ! full array of strings to write dims1(1) = length ! length of string ! Copy over string buffer to a rank 1 array @@ -764,7 +760,7 @@ contains ! character(len=length, kind=c_char), pointer :: chr_ptr ! f_ptr = c_loc(buf_ptr(1)) ! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist) -! call c_f_pointer(buf_ptr(1), chr_ptr) +! call c_f_pointer(buf_ptr(1), chr_ptr) ! buffer = chr_ptr ! nullify(chr_ptr) @@ -790,76 +786,6 @@ contains end subroutine hdf5_read_string -!=============================================================================== -! HDF5_WRITE_STRING_1DARRAY writes string data -!=============================================================================== - - subroutine hdf5_write_string_1Darray(group, name, buffer, length, rank) - - integer(HID_T), intent(in) :: group ! name of group - character(*), intent(in) :: name ! name of data - integer, intent(in) :: length ! length of strings - integer, intent(in) :: rank ! number of strings - character(*), intent(in) :: buffer(length,rank) ! data to write - - integer :: dims ! number of strings - - ! Insert null character at end of string when writing - call h5tset_strpad_f(H5T_STRING, H5T_STR_NULLPAD_F, hdf5_err) - - ! Create the dataspace and dataset - call h5screate_simple_f(1, dims1, dspace, hdf5_err) - call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err) - - ! Set up dimesnions of string to write - dims2 = (/length, rank/) ! full array of strings to write - dims = 1 ! length of string - - ! Write the variable dataset - !call h5ltmake_dataset_string_f(dset, H5T_STRING, buffer, dims2, dims1, hdf5_err, & - ! mem_space_id=dspace) - - !call h5ltmake_dataset_f(group, name, dims, dims2, H5T_STRING, buffer, hdf5_err) - - ! Close all - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - - end subroutine hdf5_write_string_1Darray - -!=============================================================================== -! HDF5_READ_STRING_1DARRAY reads string data -!=============================================================================== - - subroutine hdf5_read_string_1Darray(group, name, buffer, length, rank) - - integer(HID_T), intent(in) :: group ! name of group - character(*), intent(in) :: name ! name of data - integer, intent(in) :: length ! length of strings - integer, intent(in) :: rank ! number of strings - character(*), intent(inout) :: buffer(length,rank) ! read data to here - - ! Open dataset - call h5dopen_f(group, name, dset, hdf5_err) - - ! Get dataspace to read - call h5dget_space_f(dset, dspace, hdf5_err) - - ! Set dimensions - dims2 = (/length, rank/) - dims1(1) = length - - ! Read in the data - !call h5ltread_dataset_string_f(dset, H5T_STRING, buffer, dims2, dims1, hdf5_err, & - ! mem_space_id=dspace, xfer_prp = plist) - - !call h5ltread_dataset_f(group, name, dims, dims2, H5T_STRING, buffer, hdf5_err) - - ! Close dataset - call h5dclose_f(dset, hdf5_err) - - end subroutine hdf5_read_string_1Darray - !=============================================================================== ! HDF5_WRITE_ATTRIBUTE_STRING writes a string attribute to a variables !=============================================================================== @@ -912,7 +838,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -990,7 +916,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1070,7 +996,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1151,7 +1077,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1233,7 +1159,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1312,7 +1238,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1390,7 +1316,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1470,7 +1396,7 @@ contains f_ptr = c_loc(buffer(1,1)) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1551,7 +1477,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1633,7 +1559,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1713,7 +1639,7 @@ contains f_ptr = c_loc(buffer) call h5dwrite_f(dset, long_type, f_ptr, hdf5_err, xfer_prp=plist) - ! Close all + ! Close all call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) call h5pclose_f(plist, hdf5_err) @@ -1800,7 +1726,7 @@ contains call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err) ! Set up dimesnions of string to write - dims2 = (/length, 1/) ! full array of strings to write + dims2 = (/length, 1/) ! full array of strings to write dims1(1) = length ! length of string ! Copy over string buffer to a rank 1 array @@ -1836,7 +1762,7 @@ contains ! character(len=length, kind=c_char), pointer :: chr_ptr ! f_ptr = c_loc(buf_ptr(1)) ! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist) -! call c_f_pointer(buf_ptr(1), chr_ptr) +! call c_f_pointer(buf_ptr(1), chr_ptr) ! buffer = chr_ptr ! nullify(chr_ptr) diff --git a/src/initialize.F90 b/src/initialize.F90 index e39331675c..51dbf5ed52 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -4,7 +4,7 @@ module initialize use bank_header, only: Bank use constants use dict_header, only: DictIntInt, ElemKeyValueII - use energy_grid, only: logarithmic_grid, grid_method + use energy_grid, only: logarithmic_grid, grid_method, unionized_grid use error, only: fatal_error, warning use geometry, only: neighbor_lists use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& @@ -109,10 +109,17 @@ contains ! Create linked lists for multiple instances of the same nuclide call same_nuclide_list() - ! Construct logarithmic energy grid for cross-sections - if (grid_method == GRID_LOGARITHM) then + ! Construct unionized or log energy grid for cross-sections + select case (grid_method) + case (GRID_NUCLIDE) + continue + case (GRID_MAT_UNION) + call time_unionize % start() + call unionized_grid() + call time_unionize % stop() + case (GRID_LOGARITHM) call logarithmic_grid() - end if + end select ! Allocate and setup tally stride, matching_bins, and tally maps call configure_tallies() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 7f5587ea1b..9c02f61b45 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -213,8 +213,11 @@ contains select case (trim(temp_str)) case ('nuclide') grid_method = GRID_NUCLIDE - case ('union') - call fatal_error("Union energy grid is no longer supported.") + case ('material-union', 'union') + grid_method = GRID_MAT_UNION + if (trim(temp_str) == 'union') & + call warning('Energy grids will be unionized by material. Global& + & energy grid unionization is no longer an allowed option.') case ('logarithm', 'logarithmic', 'log') grid_method = GRID_LOGARITHM case default diff --git a/src/material_header.F90 b/src/material_header.F90 index 61aa5e5586..4c94ce98f2 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -14,6 +14,13 @@ module material_header real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm + ! Energy grid information + integer :: n_grid ! # of union material grid points + real(8), allocatable :: e_grid(:) ! union material grid energies + + ! Unionized energy grid information + integer, allocatable :: nuclide_grid_index(:,:) ! nuclide e_grid pointers + ! S(a,b) data references integer :: n_sab = 0 ! number of S(a,b) tables integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide diff --git a/src/output.F90 b/src/output.F90 index 185a96d5bc..aff46bbfcc 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1624,22 +1624,26 @@ contains ! write global tallies if (n_realizations > 1) then - write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) & - % sum, global_tallies(K_COLLISION) % sum_sq - write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) & - % sum, global_tallies(K_TRACKLENGTH) % sum_sq - write(ou,102) "k-effective (Absorption)", global_tallies(K_ABSORPTION) & - % sum, global_tallies(K_ABSORPTION) % sum_sq - if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined + if (run_mode == MODE_EIGENVALUE) then + write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) & + % sum, global_tallies(K_COLLISION) % sum_sq + write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) & + % sum, global_tallies(K_TRACKLENGTH) % sum_sq + write(ou,102) "k-effective (Absorption)", global_tallies(K_ABSORPTION) & + % sum, global_tallies(K_ABSORPTION) % sum_sq + if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined + end if write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, & global_tallies(LEAKAGE) % sum_sq else if (master) call warning("Could not compute uncertainties -- only one & &active batch simulated!") - write(ou,103) "k-effective (Collision)", global_tallies(K_COLLISION) % sum - write(ou,103) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum - write(ou,103) "k-effective (Absorption)", global_tallies(K_ABSORPTION) % sum + if (run_mode == MODE_EIGENVALUE) then + write(ou,103) "k-effective (Collision)", global_tallies(K_COLLISION) % sum + write(ou,103) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum + write(ou,103) "k-effective (Absorption)", global_tallies(K_ABSORPTION) % sum + end if write(ou,103) "Leakage Fraction", global_tallies(LEAKAGE) % sum end if write(ou,*) diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index a7b92e1738..092bf9983c 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -27,7 +27,7 @@ element settings { (element weight_avg { xsd:double } | attribute weight_avg { xsd:double })? }? & - element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" ) }? & + element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" | "material-union" | "union" ) }? & element entropy { (element dimension { list { xsd:int+ } } | diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 02dba3688e..10770f80cc 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -106,6 +106,8 @@ log logarithm logarithmic + material-union + union diff --git a/src/state_point.F90 b/src/state_point.F90 index c3edc3bff9..e30be0ea31 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -39,17 +39,13 @@ contains subroutine write_state_point() character(MAX_FILE_LEN) :: filename - integer :: i, j, k, m - integer :: n_x, n_y, n_z - character(11), allocatable :: name_array(:) + integer :: i, j, k integer, allocatable :: id_array(:) integer, allocatable :: key_array(:) - type(StructuredMesh), pointer :: mesh => null() - type(TallyObject), pointer :: tally => null() - type(ElemKeyValueII), pointer :: current => null() - type(ElemKeyValueII), pointer :: next => null() - type(ElemKeyValueCI), pointer :: cur_nuclide => null() - type(ElemKeyValueCI), pointer :: next_nuclide => null() + type(StructuredMesh), pointer :: mesh + type(TallyObject), pointer :: tally + type(ElemKeyValueII), pointer :: current + type(ElemKeyValueII), pointer :: next character(8) :: moment_name ! name of moment (e.g, P3) integer :: n_order ! loop index for moment orders integer :: nm_order ! loop index for Ynm moment orders @@ -306,14 +302,14 @@ contains case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) moment_name = 'P' // to_str(tally % moment_order(k)) call sp % write_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & + group="tallies/tally " // trim(to_str(tally % id)) // & "/moments") k = k + 1 case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) do n_order = 0, tally % moment_order(k) moment_name = 'P' // trim(to_str(n_order)) call sp % write_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & + group="tallies/tally " // trim(to_str(tally % id)) // & "/moments") k = k + 1 end do @@ -325,7 +321,7 @@ contains trim(to_str(nm_order)) call sp % write_data(moment_name, "order" // & trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & + group="tallies/tally " // trim(to_str(tally % id)) // & "/moments") k = k + 1 end do @@ -504,9 +500,9 @@ contains real(8) :: dummy ! temporary receive buffer for non-root reduces #endif integer, allocatable :: id_array(:) - type(ElemKeyValueII), pointer :: current => null() - type(ElemKeyValueII), pointer :: next => null() - type(TallyObject), pointer :: tally => null() + type(ElemKeyValueII), pointer :: current + type(ElemKeyValueII), pointer :: next + type(TallyObject), pointer :: tally type(TallyResult), allocatable :: tallyresult_temp(:,:) ! ========================================================================== @@ -661,14 +657,10 @@ contains integer, allocatable :: key_array(:) integer :: curr_key integer, allocatable :: temp_array(:) - integer, allocatable :: temp_array3D(:,:,:) - integer, allocatable :: temp_array4D(:,:,:,:) logical :: source_present - real(8) :: l real(8) :: real_array(3) - real(8), allocatable :: temp_real_array(:) - type(StructuredMesh), pointer :: mesh => null() - type(TallyObject), pointer :: tally => null() + type(StructuredMesh), pointer :: mesh + type(TallyObject), pointer :: tally integer :: n_order ! loop index for moment orders integer :: nm_order ! loop index for Ynm moment orders character(8) :: moment_name ! name of moment (e.g, P3, Y-1,1) diff --git a/src/tally.F90 b/src/tally.F90 index c96f1ecd6e..ded646d4da 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -26,6 +26,490 @@ module tally contains +!=============================================================================== +! SCORE_GENERAL adds scores to the tally array for the given filter and nuclide. +! This will work for either analog or tracklength tallies. Note that +! atom_density and flux are not used for analog tallies. +!=============================================================================== + + subroutine score_general(p, t, start_index, filter_index, i_nuclide, & + atom_density, flux) + type(Particle), intent(in) :: p + type(TallyObject), pointer, intent(inout) :: t + integer, intent(in) :: start_index + integer, intent(in) :: i_nuclide + integer, intent(in) :: filter_index ! for % results + real(8), intent(in) :: flux ! flux estimate + real(8), intent(in) :: atom_density ! atom/b-cm + + integer :: i ! loop index for scoring bins + integer :: l ! loop index for nuclides in material + integer :: m ! loop index for reactions + integer :: n ! loop index for legendre order + integer :: num_nm ! Number of N,M orders in harmonic + integer :: q ! loop index for scoring bins + integer :: i_nuc ! index in nuclides array (from material) + integer :: i_energy ! index in nuclide energy grid + integer :: score_bin ! scoring bin, e.g. SCORE_FLUX + integer :: score_index ! scoring bin index + real(8) :: atom_density_ ! atom/b-cm + real(8) :: f ! interpolation factor + real(8) :: score ! analog tally score + real(8) :: macro_total ! material macro total xs + real(8) :: macro_scatt ! material macro scatt xs + real(8) :: uvw(3) ! particle direction + type(Material), pointer, save :: mat => null() + type(Reaction), pointer, save :: rxn => null() +!$omp threadprivate(mat, rxn) + + i = 0 + SCORE_LOOP: do q = 1, t % n_user_score_bins + i = i + 1 + + ! determine what type of score bin + score_bin = t % score_bins(i) + + ! determine scoring bin index + score_index = start_index + i + + !######################################################################### + ! Determine appropirate scoring value. + + select case(score_bin) + + + case (SCORE_FLUX, SCORE_FLUX_YN) + if (t % estimator == ESTIMATOR_ANALOG) then + ! All events score to a flux bin. We actually use a collision + ! estimator since there is no way to count 'events' exactly for + ! the flux + if (survival_biasing) then + ! We need to account for the fact that some weight was already + ! absorbed + score = p % last_wgt + p % absorb_wgt + else + score = p % last_wgt + end if + score = score / material_xs % total + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + ! For flux, we need no cross section + score = flux + end if + + + case (SCORE_TOTAL, SCORE_TOTAL_YN) + if (t % estimator == ESTIMATOR_ANALOG) then + ! All events will score to the total reaction rate. We can just + ! use the weight of the particle entering the collision as the + ! score + if (survival_biasing) then + ! We need to account for the fact that some weight was already + ! absorbed + score = p % last_wgt + p % absorb_wgt + else + score = p % last_wgt + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % total * atom_density * flux + else + score = material_xs % total * flux + end if + end if + + + case (SCORE_SCATTER, SCORE_SCATTER_N) + if (t % estimator == ESTIMATOR_ANALOG) then + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + ! Since only scattering events make it here, again we can use + ! the weight entering the collision as the estimator for the + ! reaction rate + score = p % last_wgt + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + ! Note SCORE_SCATTER_N not available for tracklength. + if (i_nuclide > 0) then + score = (micro_xs(i_nuclide) % total & + - micro_xs(i_nuclide) % absorption) * atom_density * flux + else + score = (material_xs % total - material_xs % absorption) * flux + end if + end if + + + case (SCORE_SCATTER_PN, SCORE_SCATTER_YN) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) then + i = i + t % moment_order(i) + cycle SCORE_LOOP + end if + ! Since only scattering events make it here, again we can use + ! the weight entering the collision as the estimator for the + ! reaction rate + score = p % last_wgt + + + case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + ! For scattering production, we need to use the post-collision + ! weight as the estimate for the number of neutrons exiting a + ! reaction with neutrons in the exit channel + score = p % wgt + + + case (SCORE_NU_SCATTER_PN, SCORE_NU_SCATTER_YN) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) then + i = i + t % moment_order(i) + cycle SCORE_LOOP + end if + ! For scattering production, we need to use the post-collision + ! weight as the estimate for the number of neutrons exiting a + ! reaction with neutrons in the exit channel + score = p % wgt + + + case (SCORE_TRANSPORT) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + ! get material macros + macro_total = material_xs % total + macro_scatt = material_xs % total - material_xs % absorption + ! Score total rate - p1 scatter rate Note estimator needs to be + ! adjusted since tallying is only occuring when a scatter has + ! happened. Effectively this means multiplying the estimator by + ! total/scatter macro + score = (macro_total - p % mu * macro_scatt) * (ONE / macro_scatt) + + + case (SCORE_N_1N) + ! Only analog estimators are available. + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + ! Skip any events where weight of particle changed + if (p % wgt /= p % last_wgt) cycle SCORE_LOOP + ! All events that reach this point are (n,1n) reactions + score = p % last_wgt + + + case (SCORE_ABSORPTION) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No absorption events actually occur if survival biasing is on -- + ! just use weight absorbed in survival biasing + score = p % absorb_wgt + else + ! Skip any event where the particle wasn't absorbed + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP + ! All fission and absorption events will contribute here, so we + ! can just use the particle's weight entering the collision + score = p % last_wgt + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % absorption * atom_density * flux + else + score = material_xs % absorption * flux + end if + end if + + + case (SCORE_FISSION) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! fission + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if + else + ! Skip any non-absorption events + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP + ! All fission events will contribute, so again we can use + ! particle's weight entering the collision as the estimate for the + ! fission reaction rate + score = p % last_wgt * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % fission * atom_density * flux + else + score = material_xs % fission * flux + end if + end if + + + case (SCORE_NU_FISSION) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing .or. p % fission) then + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + ! Normally, we only need to make contributions to one scoring + ! bin. However, in the case of fission, since multiple fission + ! neutrons were emitted with different energies, multiple + ! outgoing energy bins may have been scored to. The following + ! logic treats this special case and results to multiple bins + call score_fission_eout(p, t, score_index) + cycle SCORE_LOOP + end if + end if + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! nu-fission + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * micro_xs(p % event_nuclide) % & + nu_fission / micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if + else + ! Skip any non-fission events + if (.not. p % fission) cycle SCORE_LOOP + ! If there is no outgoing energy filter, than we only need to + ! score to one bin. For the score to be 'analog', we need to + ! score the number of particles that were banked in the fission + ! bank. Since this was weighted by 1/keff, we multiply by keff + ! to get the proper score. + score = keff * p % wgt_bank + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % nu_fission * atom_density * flux + else + score = material_xs % nu_fission * flux + end if + end if + + + case (SCORE_KAPPA_FISSION) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! fission scale by kappa-fission + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * & + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if + else + ! Skip any non-absorption events + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP + ! All fission events will contribute, so again we can use + ! particle's weight entering the collision as the estimate for + ! the fission energy production rate + score = p % last_wgt * & + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption + end if + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux + else + score = material_xs % kappa_fission * flux + end if + end if + + + case (SCORE_EVENTS) + ! Simply count number of scoring events + score = ONE + + + case default + if (t % estimator == ESTIMATOR_ANALOG) then + ! Any other score is assumed to be a MT number. Thus, we just need + ! to check if it matches the MT number of the event + if (p % event_MT /= score_bin) cycle SCORE_LOOP + score = p % last_wgt + + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + ! Any other cross section has to be calculated on-the-fly. For + ! cross sections that are used often (e.g. n2n, ngamma, etc. for + ! depletion), it might make sense to optimize this section or + ! pre-calculate cross sections + if (score_bin > 1) then + ! Set default score + score = ZERO + + if (i_nuclide > 0) then + ! TODO: The following search for the matching reaction could + ! be replaced by adding a dictionary on each Nuclide instance + ! of the form {MT: i_reaction, ...} + REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction + ! Get pointer to reaction + rxn => nuclides(i_nuclide) % reactions(m) + ! Check if this is the desired MT + if (score_bin == rxn % MT) then + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor + if (i_energy >= rxn % threshold) then + score = ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density * flux + end if + exit REACTION_LOOP + end if + end do REACTION_LOOP + + else + ! Get pointer to current material + mat => materials(p % material) + do l = 1, mat % n_nuclides + ! Get atom density + atom_density_ = mat % atom_density(l) + ! Get index in nuclides array + i_nuc = mat % nuclide(l) + ! TODO: The following search for the matching reaction could + ! be replaced by adding a dictionary on each Nuclide + ! instance of the form {MT: i_reaction, ...} + do m = 1, nuclides(i_nuc) % n_reaction + ! Get pointer to reaction + rxn => nuclides(i_nuc) % reactions(m) + ! Check if this is the desired MT + if (score_bin == rxn % MT) then + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuc) % index_grid + f = micro_xs(i_nuc) % interp_factor + if (i_energy >= rxn % threshold) then + score = score + ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density_ * flux + end if + exit + end if + end do + end do + end if + + else + call fatal_error("Invalid score type on tally " & + // to_str(t % id) // ".") + end if + end if + + + end select + + !######################################################################### + ! Expand score if necessary and add to tally results. + + select case(score_bin) + + + case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) + ! Find the scattering order for a singly requested moment, and + ! store its moment contribution. + if (t % moment_order(i) == 1) then + score = score * p % mu ! avoid function call overhead + else + score = score * calc_pn(t % moment_order(i), p % mu) + endif +!$omp atomic + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score + + + case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN) + score_index = score_index - 1 + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(i) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) & + % value = t & + % results(score_index: score_index + num_nm - 1, filter_index)& + % value & + + score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw) +!$omp end critical + end do + i = i + (t % moment_order(i) + 1)**2 - 1 + + + case(SCORE_FLUX_YN, SCORE_TOTAL_YN) + score_index = score_index - 1 + num_nm = 1 + if (t % estimator == ESTIMATOR_ANALOG) then + uvw = p % last_uvw + else if (t % estimator == ESTIMATOR_TRACKLENGTH) then + uvw = p % coord0 % uvw + end if + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(i) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) & + % value = t & + % results(score_index: score_index + num_nm - 1, filter_index)& + % value & + + score * calc_rn(n, uvw) +!$omp end critical + end do + i = i + (t % moment_order(i) + 1)**2 - 1 + + + case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) + score_index = score_index - 1 + ! Find the scattering order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(i) + ! determine scoring bin index + score_index = score_index + 1 + + ! get the score and tally it +!$omp atomic + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value & + + score * calc_pn(n, p % mu) + end do + i = i + t % moment_order(i) + + + case default +!$omp atomic + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score + + + end select + end do SCORE_LOOP + end subroutine score_general + !=============================================================================== ! SCORE_ANALOG_TALLY keeps track of how many events occur in a specified cell, ! energy range, etc. Note that since these are "analog" tallies, they are only @@ -38,22 +522,13 @@ contains integer :: i integer :: i_tally - integer :: j ! loop index for scoring bins integer :: k ! loop index for nuclide bins - integer :: n ! loop index for legendre order - integer :: num_nm ! Number of N,M orders in harmonic - integer :: l ! scoring bin loop index, allowing for changing ! position during the loop integer :: filter_index ! single index for single bin - integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: i_nuclide ! index in nuclides array - integer :: score_index ! scoring bin index - real(8) :: score ! analog tally score real(8) :: last_wgt ! pre-collision particle weight real(8) :: wgt ! post-collision particle weight real(8) :: mu ! cosine of angle of collision - real(8) :: macro_total ! material macro total xs - real(8) :: macro_scatt ! material macro scatt xs logical :: found_bin ! scoring bin found? type(TallyObject), pointer, save :: t => null() !$omp threadprivate(t) @@ -124,427 +599,8 @@ contains end if ! Determine score for each bin - j = 0 - SCORE_LOOP: do l = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) - - ! determine scoring bin index - score_index = (k - 1)*t % n_score_bins + j - - select case (score_bin) - case (SCORE_FLUX) - ! All events score to a flux bin. We actually use a collision - ! estimator since there is no way to count 'events' exactly for - ! the flux - - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - - score = score / material_xs % total - case (SCORE_FLUX_YN) - ! All events score to a flux bin. We actually use a collision - ! estimator since there is no way to count 'events' exactly for - ! the flux - - score_index = score_index - 1 - - ! get the score - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - - score = score / material_xs % total - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % last_uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - ! All events will score to the total reaction rate. We can just - ! use the weight of the particle entering the collision as the - ! score - - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - - case (SCORE_TOTAL_YN) - ! All events will score to the total reaction rate. We can just - ! use the weight of the particle entering the collision as the - ! score - - score_index = score_index - 1 - - ! get the score - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % last_uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! Since only scattering events make it here, again we can use - ! the weight entering the collision as the estimator for the - ! reaction rate - - score = last_wgt - - case (SCORE_NU_SCATTER) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! For scattering production, we need to use the post-collision - ! weight as the estimate for the number of neutrons exiting a - ! reaction with neutrons in the exit channel - - score = wgt - - case (SCORE_SCATTER_N) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! Find the scattering order for a singly requested moment, and - ! store its moment contribution. - - if (t % moment_order(j) == 1) then - score = last_wgt * mu ! avoid function call overhead - else - score = last_wgt * calc_pn(t % moment_order(j), mu) - endif - - case (SCORE_SCATTER_PN) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) - cycle SCORE_LOOP - end if - score_index = score_index - 1 - ! Find the scattering order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + 1 - ! get the score and tally it - score = last_wgt * calc_pn(n, mu) - -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - end do - j = j + t % moment_order(j) - cycle SCORE_LOOP - - case (SCORE_SCATTER_YN) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) - cycle SCORE_LOOP - end if - score_index = score_index - 1 - - ! Calculate the number of moments from t % moment_order - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - ! get the score of the scattering moment - score = last_wgt * calc_pn(n, mu) - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % last_uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_NU_SCATTER_N) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! Find the scattering order for a singly requested moment, and - ! store its moment contribution. - - if (t % moment_order(j) == 1) then - score = wgt * mu ! avoid function call overhead - else - score = wgt * calc_pn(t % moment_order(j), mu) - endif - - case (SCORE_NU_SCATTER_PN) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) - cycle SCORE_LOOP - end if - score_index = score_index - 1 - ! Find the scattering order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + 1 - ! get the score and tally it - score = wgt * calc_pn(n, mu) - -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - end do - j = j + t % moment_order(j) - cycle SCORE_LOOP - - case (SCORE_NU_SCATTER_YN) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) then - j = j + t % moment_order(j) - cycle SCORE_LOOP - end if - score_index = score_index - 1 - - ! Calculate the number of moments from t % moment_order - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - ! get the score of the scattering moment - score = wgt * calc_pn(n, mu) - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % last_uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TRANSPORT) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! get material macros - macro_total = material_xs % total - macro_scatt = material_xs % total - material_xs % absorption - - ! Score total rate - p1 scatter rate Note estimator needs to be - ! adjusted since tallying is only occuring when a scatter has - ! happened. Effectively this means multiplying the estimator by - ! total/scatter macro - score = (macro_total - mu*macro_scatt)*(ONE/macro_scatt) - - case (SCORE_N_1N) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! Skip any events where weight of particle changed - if (wgt /= last_wgt) cycle SCORE_LOOP - - ! All events that reach this point are (n,1n) reactions - score = last_wgt - - case (SCORE_ABSORPTION) - if (survival_biasing) then - ! No absorption events actually occur if survival biasing is on -- - ! just use weight absorbed in survival biasing - - score = p % absorb_wgt - - else - ! Skip any event where the particle wasn't absorbed - if (p % event == EVENT_SCATTER) cycle SCORE_LOOP - - ! All fission and absorption events will contribute here, so we - ! can just use the particle's weight entering the collision - - score = last_wgt - end if - - case (SCORE_FISSION) - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! fission - - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption - else - score = ZERO - end if - - else - ! Skip any non-absorption events - if (p % event == EVENT_SCATTER) cycle SCORE_LOOP - - ! All fission events will contribute, so again we can use - ! particle's weight entering the collision as the estimate for the - ! fission reaction rate - - score = last_wgt * micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption - end if - - case (SCORE_NU_FISSION) - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! nu-fission - - if (t % find_filter(FILTER_ENERGYOUT) > 0) then - ! Normally, we only need to make contributions to one scoring - ! bin. However, in the case of fission, since multiple fission - ! neutrons were emitted with different energies, multiple - ! outgoing energy bins may have been scored to. The following - ! logic treats this special case and results to multiple bins - - call score_fission_eout(p, t, score_index) - cycle SCORE_LOOP - - else - - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * micro_xs(p % event_nuclide) % & - nu_fission / micro_xs(p % event_nuclide) % absorption - else - score = ZERO - end if - end if - - else - ! Skip any non-fission events - if (.not. p % fission) cycle SCORE_LOOP - - if (t % find_filter(FILTER_ENERGYOUT) > 0) then - ! Normally, we only need to make contributions to one scoring - ! bin. However, in the case of fission, since multiple fission - ! neutrons were emitted with different energies, multiple - ! outgoing energy bins may have been scored to. The following - ! logic treats this special case and results to multiple bins - - call score_fission_eout(p, t, score_index) - cycle SCORE_LOOP - - else - ! If there is no outgoing energy filter, than we only need to - ! score to one bin. For the score to be 'analog', we need to - ! score the number of particles that were banked in the fission - ! bank. Since this was weighted by 1/keff, we multiply by keff - ! to get the proper score. - - score = keff * p % wgt_bank - - end if - end if - - case (SCORE_KAPPA_FISSION) - if (survival_biasing) then - ! No fission events occur if survival biasing is on -- need to - ! calculate fraction of absorptions that would have resulted in - ! fission scale by kappa-fission - - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * & - micro_xs(p % event_nuclide) % kappa_fission / & - micro_xs(p % event_nuclide) % absorption - else - score = ZERO - end if - - else - ! Skip any non-absorption events - if (p % event == EVENT_SCATTER) cycle SCORE_LOOP - - ! All fission events will contribute, so again we can use - ! particle's weight entering the collision as the estimate for - ! the fission energy production rate - score = last_wgt * & - micro_xs(p % event_nuclide) % kappa_fission / & - micro_xs(p % event_nuclide) % absorption - end if - case (SCORE_EVENTS) - ! Simply count number of scoring events - score = ONE - - case default - ! Any other score is assumed to be a MT number. Thus, we just need - ! to check if it matches the MT number of the event - if (p % event_MT /= score_bin) cycle SCORE_LOOP - - score = last_wgt - - end select - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do SCORE_LOOP + call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & + i_nuclide, ZERO, ZERO) end do NUCLIDE_LOOP @@ -640,26 +696,14 @@ contains integer :: i_tally integer :: j ! loop index for scoring bins integer :: k ! loop index for nuclide bins - integer :: l ! loop index for nuclides in material - integer :: m ! loop index for reactions - integer :: n ! loop index for legendre order - integer :: num_nm ! Number of N,M orders in harmonic - integer :: q ! loop index for scoring bins integer :: filter_index ! single index for single bin integer :: i_nuclide ! index in nuclides array (from bins) - integer :: i_nuc ! index in nuclides array (from material) - integer :: i_energy ! index in nuclide energy grid - integer :: score_bin ! scoring type, e.g. SCORE_FLUX - integer :: score_index ! scoring bin index - real(8) :: f ! interpolation factor real(8) :: flux ! tracklength estimate of flux - real(8) :: score ! actual score (e.g., flux*xs) real(8) :: atom_density ! atom density of single nuclide in atom/b-cm logical :: found_bin ! scoring bin found? type(TallyObject), pointer, save :: t => null() type(Material), pointer, save :: mat => null() - type(Reaction), pointer, save :: rxn => null() -!$omp threadprivate(t, mat, rxn) +!$omp threadprivate(t, mat) ! Determine track-length estimate of flux flux = p % wgt * distance @@ -731,291 +775,8 @@ contains end if ! Determine score for each bin - j = 0 - SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) - - ! determine scoring bin index - score_index = (k - 1)*t % n_score_bins + j - - if (i_nuclide > 0) then - ! ================================================================ - ! DETERMINE NUCLIDE CROSS SECTION - - select case(score_bin) - case (SCORE_FLUX) - ! For flux, we need no cross section - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - ! For flux, we need no cross section - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - ! Total cross section is pre-calculated - score = micro_xs(i_nuclide) % total * & - atom_density * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - ! Scattering cross section is pre-calculated - score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux - - case (SCORE_ABSORPTION) - ! Absorption cross section is pre-calculated - score = micro_xs(i_nuclide) % absorption * & - atom_density * flux - - case (SCORE_FISSION) - ! Fission cross section is pre-calculated - score = micro_xs(i_nuclide) % fission * & - atom_density * flux - - case (SCORE_NU_FISSION) - ! Nu-fission cross section is pre-calculated - score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux - - case (SCORE_KAPPA_FISSION) - score = micro_xs(i_nuclide) % kappa_fission * & - atom_density * flux - - case (SCORE_EVENTS) - ! For number of events, just score unity - score = ONE - - case default - ! Any other cross section has to be calculated on-the-fly. For - ! cross sections that are used often (e.g. n2n, ngamma, etc. for - ! depletion), it might make sense to optimize this section or - ! pre-calculate cross sections - - if (score_bin > 1) then - ! Set default score - score = ZERO - - ! TODO: The following search for the matching reaction could - ! be replaced by adding a dictionary on each Nuclide instance - ! of the form {MT: i_reaction, ...} - - REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction - ! Get pointer to reaction - rxn => nuclides(i_nuclide) % reactions(m) - - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - - if (i_energy >= rxn % threshold) then - score = ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - - exit REACTION_LOOP - end if - end do REACTION_LOOP - - else - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end if - end select - - else - ! ================================================================ - ! DETERMINE MATERIAL CROSS SECTION - - select case(score_bin) - case (SCORE_FLUX) - ! For flux, we need no cross section - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - ! For flux, we need no cross section - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - ! Total cross section is pre-calculated - score = material_xs % total * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - ! Total cross section is pre-calculated - score = material_xs % total * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - ! Scattering cross section is pre-calculated - score = (material_xs % total - material_xs % absorption) * flux - - case (SCORE_ABSORPTION) - ! Absorption cross section is pre-calculated - score = material_xs % absorption * flux - - case (SCORE_FISSION) - ! Fission cross section is pre-calculated - score = material_xs % fission * flux - - case (SCORE_NU_FISSION) - ! Nu-fission cross section is pre-calculated - score = material_xs % nu_fission * flux - - case (SCORE_KAPPA_FISSION) - score = material_xs % kappa_fission * flux - - case (SCORE_EVENTS) - ! For number of events, just score unity - score = ONE - - case default - ! Any other cross section has to be calculated on-the-fly. This - ! is somewhat costly since it requires a loop over each nuclide - ! in a material and each reaction in the nuclide - - if (score_bin > 1) then - ! Set default score - score = ZERO - - ! Get pointer to current material - mat => materials(p % material) - - do l = 1, mat % n_nuclides - ! Get atom density - atom_density = mat % atom_density(l) - - ! Get index in nuclides array - i_nuc = mat % nuclide(l) - - ! TODO: The following search for the matching reaction could - ! be replaced by adding a dictionary on each Nuclide - ! instance of the form {MT: i_reaction, ...} - - do m = 1, nuclides(i_nuc) % n_reaction - ! Get pointer to reaction - rxn => nuclides(i_nuc) % reactions(m) - - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuc) % index_grid - f = micro_xs(i_nuc) % interp_factor - - if (i_energy >= rxn % threshold) then - score = score + ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - - exit - end if - end do - - end do - - else - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end if - end select - end if - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do SCORE_LOOP + call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & + i_nuclide, atom_density, flux) end do NUCLIDE_BIN_LOOP end if @@ -1047,22 +808,11 @@ contains integer, intent(in) :: filter_index integer :: i ! loop index for nuclides in material - integer :: j ! loop index for scoring bin types - integer :: m ! loop index for reactions in nuclide - integer :: n ! loop index for legendre order - integer :: num_nm ! Number of N,M orders in harmonic - integer :: q ! loop index for scoring bins integer :: i_nuclide ! index in nuclides array - integer :: score_bin ! type of score, e.g. SCORE_FLUX - integer :: score_index ! scoring bin index - integer :: i_energy ! index in nuclide energy grid - real(8) :: f ! interpolation factor - real(8) :: score ! actual scoring tally value real(8) :: atom_density ! atom density of single nuclide in atom/b-cm type(TallyObject), pointer, save :: t => null() type(Material), pointer, save :: mat => null() - type(Reaction), pointer, save :: rxn => null() -!$omp threadprivate(t, mat, rxn) +!$omp threadprivate(t, mat) ! Get pointer to tally t => tallies(i_tally) @@ -1081,290 +831,21 @@ contains i_nuclide = mat % nuclide(i) atom_density = mat % atom_density(i) - ! Loop over score types for each bin - j = 0 - SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) - - ! Determine scoring bin index based on what the index of the nuclide - ! is in the nuclides array - score_index = (i_nuclide - 1)*t % n_score_bins + j - - ! Determine macroscopic nuclide cross section - select case(score_bin) - case (SCORE_FLUX) - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - ! For flux, we need no cross section - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - score = micro_xs(i_nuclide) % total * atom_density * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - score = micro_xs(i_nuclide) % total * atom_density * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * atom_density * flux - - case (SCORE_ABSORPTION) - score = micro_xs(i_nuclide) % absorption * atom_density * flux - - case (SCORE_FISSION) - score = micro_xs(i_nuclide) % fission * atom_density * flux - - case (SCORE_NU_FISSION) - score = micro_xs(i_nuclide) % nu_fission * atom_density * flux - - case (SCORE_KAPPA_FISSION) - score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux - - case (SCORE_EVENTS) - score = ONE - - case default - ! Any other cross section has to be calculated on-the-fly. For cross - ! sections that are used often (e.g. n2n, ngamma, etc. for depletion), - ! it might make sense to optimize this section or pre-calculate cross - ! sections - - if (score_bin > 1) then - ! Set default score - score = ZERO - - ! TODO: The following search for the matching reaction could be - ! replaced by adding a dictionary on each Nuclide instance of the - ! form {MT: i_reaction, ...} - - REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction - ! Get pointer to reaction - rxn => nuclides(i_nuclide) % reactions(m) - - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - - if (i_energy >= rxn % threshold) then - score = ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - - exit REACTION_LOOP - end if - end do REACTION_LOOP - - else - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end if - end select - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do SCORE_LOOP + ! Determine score for each bin + call score_general(p, t, (i_nuclide-1)*t % n_score_bins, filter_index, & + i_nuclide, atom_density, flux) end do NUCLIDE_LOOP ! ========================================================================== ! SCORE TOTAL MATERIAL REACTION RATES - ! Loop over score types for each bin - j = 0 - MATERIAL_SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) + i_nuclide = -1 + atom_density = ZERO - ! Determine scoring bin index based on what the index of the nuclide - ! is in the nuclides array - score_index = n_nuclides_total*t % n_score_bins + j - - ! Determine macroscopic material cross section - select case(score_bin) - case (SCORE_FLUX) - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - ! For flux, we need no cross section - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle MATERIAL_SCORE_LOOP - - case (SCORE_TOTAL) - score = material_xs % total * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - ! Total cross section is pre-calculated - score = material_xs % total * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle MATERIAL_SCORE_LOOP - - case (SCORE_SCATTER) - score = (material_xs % total - material_xs % absorption) * flux - - case (SCORE_ABSORPTION) - score = material_xs % absorption * flux - - case (SCORE_FISSION) - score = material_xs % fission * flux - - case (SCORE_NU_FISSION) - score = material_xs % nu_fission * flux - - case (SCORE_KAPPA_FISSION) - score = material_xs % kappa_fission * flux - - case (SCORE_EVENTS) - score = ONE - - case default - ! Any other cross section has to be calculated on-the-fly. This is - ! somewhat costly since it requires a loop over each nuclide in a - ! material and each reaction in the nuclide - - if (score_bin > 1) then - ! Set default score - score = ZERO - - ! Get pointer to current material - mat => materials(p % material) - - do i = 1, mat % n_nuclides - ! Get atom density - atom_density = mat % atom_density(i) - - ! Get index in nuclides array - i_nuclide = mat % nuclide(i) - - ! TODO: The following search for the matching reaction could - ! be replaced by adding a dictionary on each Nuclide - ! instance of the form {MT: i_reaction, ...} - - do m = 1, nuclides(i_nuclide) % n_reaction - ! Get pointer to reaction - rxn => nuclides(i_nuclide) % reactions(m) - - ! Check if this is the desired MT - if (score_bin == rxn % MT) then - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - - if (i_energy >= rxn % threshold) then - score = score + ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - - exit - end if - end do - - end do - - else - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end if - end select - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do MATERIAL_SCORE_LOOP + ! Determine score for each bin + call score_general(p, t, n_nuclides_total*t % n_score_bins, filter_index, & + i_nuclide, atom_density, flux) end subroutine score_all_nuclides @@ -1384,21 +865,15 @@ contains integer :: j ! loop index for direction integer :: k ! loop index for mesh cell crossings integer :: b ! loop index for nuclide bins - integer :: n ! loop index for legendre order - integer :: num_nm ! Number of N,M orders in harmonic - integer :: q ! loop index for scoring bins integer :: ijk0(3) ! indices of starting coordinates integer :: ijk1(3) ! indices of ending coordinates integer :: ijk_cross(3) ! indices of mesh cell crossed integer :: n_cross ! number of surface crossings integer :: filter_index ! single index for single bin - integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: i_nuclide ! index in nuclides array - integer :: score_index ! scoring bin index integer :: i_filter_mesh ! index of mesh filter in filters array real(8) :: atom_density ! density of individual nuclide in atom/b-cm real(8) :: flux ! tracklength estimate of flux - real(8) :: score ! actual score (e.g., flux*xs) real(8) :: uvw(3) ! cosine of angle of particle real(8) :: xyz0(3) ! starting/intermediate coordinates real(8) :: xyz1(3) ! ending coordinates of particle @@ -1612,179 +1087,8 @@ contains end if ! Determine score for each bin - j = 0 - SCORE_LOOP: do q = 1, t % n_user_score_bins - j = j + 1 - ! determine what type of score bin - score_bin = t % score_bins(j) - - ! Determine scoring bin index - score_index = (b - 1)*t % n_score_bins + j - - if (i_nuclide > 0) then - ! Determine macroscopic nuclide cross section - select case(score_bin) - case (SCORE_FLUX) - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - score = micro_xs(i_nuclide) % total * & - atom_density * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - ! Total cross section is pre-calculated - score = micro_xs(i_nuclide) % total * & - atom_density * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux - case (SCORE_ABSORPTION) - score = micro_xs(i_nuclide) % absorption * & - atom_density * flux - case (SCORE_FISSION) - score = micro_xs(i_nuclide) % fission * & - atom_density * flux - case (SCORE_NU_FISSION) - score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux - case (SCORE_KAPPA_FISSION) - score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux - case (SCORE_EVENTS) - score = ONE - case default - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end select - - else - ! Determine macroscopic material cross section - select case(score_bin) - case (SCORE_FLUX) - score = flux - - case (SCORE_FLUX_YN) - score_index = score_index - 1 - - score = flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_TOTAL) - score = material_xs % total * flux - - case (SCORE_TOTAL_YN) - score_index = score_index - 1 - - ! Total cross section is pre-calculated - score = material_xs % total * flux - - num_nm = 1 - ! Find the order for a collection of requested moments - ! and store the moment contribution of each - do n = 0, t % moment_order(j) - ! determine scoring bin index - score_index = score_index + num_nm - ! Update number of total n,m bins for this n (m = [-n: n]) - num_nm = 2 * n + 1 - - ! multiply score by the angular flux moments and store -!$omp critical - t % results(score_index: score_index + num_nm - 1, filter_index) % value = & - t % results(score_index: score_index + num_nm - 1, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) -!$omp end critical - end do - j = j + (t % moment_order(j) + 1)**2 - 1 - cycle SCORE_LOOP - - case (SCORE_SCATTER) - score = (material_xs % total - material_xs % absorption) * flux - case (SCORE_ABSORPTION) - score = material_xs % absorption * flux - case (SCORE_FISSION) - score = material_xs % fission * flux - case (SCORE_NU_FISSION) - score = material_xs % nu_fission * flux - case (SCORE_KAPPA_FISSION) - score = material_xs % kappa_fission * flux - case (SCORE_EVENTS) - score = ONE - case default - call fatal_error("Invalid score type on tally " & - &// to_str(t % id) // ".") - end select - end if - - ! Add score to tally -!$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score - - end do SCORE_LOOP + call score_general(p, t, (b-1)*t % n_score_bins, filter_index, & + i_nuclide, atom_density, flux) end do NUCLIDE_BIN_LOOP end if diff --git a/src/utils/openmc/__init__.py b/src/utils/openmc/__init__.py index 2f3a9b6de4..42ee1eebfa 100644 --- a/src/utils/openmc/__init__.py +++ b/src/utils/openmc/__init__.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.element import * from openmc.geometry import * from openmc.nuclide import * diff --git a/src/utils/openmc/checkvalue.py b/src/utils/openmc/checkvalue.py index 679b5c2e42..503e7e24d5 100644 --- a/src/utils/openmc/checkvalue.py +++ b/src/utils/openmc/checkvalue.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - import numpy as np @@ -13,4 +11,3 @@ def is_float(val): def is_string(val): return isinstance(val, (str, np.str)) - diff --git a/src/utils/openmc/clean_xml.py b/src/utils/openmc/clean_xml.py index cac5833dc3..619475fa07 100644 --- a/src/utils/openmc/clean_xml.py +++ b/src/utils/openmc/clean_xml.py @@ -1,12 +1,10 @@ -#!/usr/bin/env python - def sort_xml_elements(tree): # Retrieve all children of the root XML node in the tree elements = tree.getchildren() # Initialize empty lists for the sorted and comment elements - sorted_elements = list() + sorted_elements = [] # Initialize an empty set of tags (e.g., Surface, Cell, and Lattice) tags = set() @@ -16,7 +14,7 @@ def sort_xml_elements(tree): tags.add(element.tag) # Initialize an empty list for the comment elements - comment_elements = list() + comment_elements = [] # Find the comment elements and record their ordering within the # tree using a precedence with respect to the subsequent nodes @@ -40,7 +38,7 @@ def sort_xml_elements(tree): continue # Initialize an empty list of tuples to sort (id, element) - tagged_data = list() + tagged_data = [] # Retrieve the IDs for each of the elements for element in tagged_elements: @@ -68,11 +66,11 @@ def sort_xml_elements(tree): def clean_xml_indentation(element, level=0): - ''' + """ copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint it basically walks your tree and adds spaces and newlines so the tree is printed in a nice way - ''' + """ i = "\n" + level*" " @@ -92,4 +90,4 @@ def clean_xml_indentation(element, level=0): else: if level and (not element.tail or not element.tail.strip()): - element.tail = i \ No newline at end of file + element.tail = i diff --git a/src/utils/openmc/cmfd.py b/src/utils/openmc/cmfd.py index 4b7802ba3d..6e7eb162a8 100644 --- a/src/utils/openmc/cmfd.py +++ b/src/utils/openmc/cmfd.py @@ -1,9 +1,9 @@ -#!/usr/bin/env python +from xml.etree import ElementTree as ET + +import numpy as np from openmc.checkvalue import * from openmc.clean_xml import * -from xml.etree import ElementTree as ET -import numpy as np class CMFDMesh(object): @@ -19,7 +19,43 @@ class CMFDMesh(object): self._map = None - def set_lower_left(self, lower_left): + @property + def lower_left(self): + return self._lower_left + + + @property + def upper_right(self): + return self._upper_right + + + @property + def dimension(self): + return self._dimension + + + @property + def width(self): + return self._width + + + @property + def energy(self): + return self._energy + + + @property + def albedo(self): + return self._albedo + + + @property + def map(self): + return self._mape + + + @lower_left.setter + def lower_left(self, lower_left): if not isinstance(lower_left, (tuple, list, np.ndarray)): msg = 'Unable to set CMFD Mesh with lower_left {0} which is ' \ @@ -41,7 +77,8 @@ class CMFDMesh(object): self._lower_left = lower_left - def set_upper_right(self, upper_right): + @upper_right.setter + def upper_right(self, upper_right): if not isinstance(upper_right, (tuple, list, np.ndarray)): msg = 'Unable to set CMFD Mesh with upper_right {0} which is ' \ @@ -63,7 +100,8 @@ class CMFDMesh(object): self._upper_right = upper_right - def set_dimension(self, dimension): + @dimension.setter + def dimension(self, dimension): if not isinstance(dimension, (tuple, list, np.ndarray)): msg = 'Unable to set CMFD Mesh with dimension {0} which is ' \ @@ -85,7 +123,8 @@ class CMFDMesh(object): self._dimension = dimension - def set_width(self, width): + @width.setter + def width(self, width): if not width is None: @@ -109,7 +148,8 @@ class CMFDMesh(object): self._width = width - def set_energy(self, energy): + @energy.setter + def energy(self, energy): if not isinstance(energy, (tuple, list, np.ndarray)): msg = 'Unable to set CMFD Mesh energy to {0} which is not ' \ @@ -131,7 +171,8 @@ class CMFDMesh(object): self._energy = energy - def set_albedo(self, albedo): + @albedo.setter + def albedo(self, albedo): if not isinstance(albedo, (tuple, list, np.ndarray)): msg = 'Unable to set CMFD Mesh albedo to {0} which is not ' \ @@ -158,7 +199,8 @@ class CMFDMesh(object): self._albedo = albedo - def set_map(self, map): + @map.setter + def map(self, map): if not isinstance(map, (tuple, list, np.ndarray)): msg = 'Unable to set CMFD Mesh map to {0} which is not ' \ @@ -278,7 +320,83 @@ class CMFDFile(object): self._cmfd_mesh_element = None - def set_active_flush(self, active_flush): + @property + def active_flush(self): + return self._active_flush + + + @property + def begin(self): + return self._begin + + + @property + def display(self): + return self._display + + + @property + def feedback(self): + return self._feedback + + + @property + def inactive(self): + return self._inactive + + + @property + def inactive_flush(self): + return self._inactive_flush + + + @property + def ksp_monitor(self): + return self._ksp_monitor + + + @property + def cmfd_mesh(self): + return self._cmfd_mesh + + + @property + def norm(self): + return self._norm + + + @property + def num_flushes(self): + return self._num_flushes + + + @property + def power_monitor(self): + return self._power_monitor + + + @property + def run_adjoint(self): + return self._run_adjoint + + + @property + def snes_monitor(self): + return self._snes_monitor + + + @property + def solver(self): + return self._solver + + + @property + def write_matrices(self): + return self._write_matrices + + + @active_flush.setter + def active_flush(self, active_flush): if not is_integer(active_flush): msg = 'Unable to set CMFD active flush batch to a non-integer ' \ @@ -293,7 +411,8 @@ class CMFDFile(object): self._active_flush = active_flush - def set_begin(self, begin): + @begin.setter + def begin(self, begin): if not is_integer(begin): msg = 'Unable to set CMFD begin batch to a non-integer ' \ @@ -308,7 +427,8 @@ class CMFDFile(object): self._begin = begin - def set_display(self, display): + @display.setter + def display(self, display): if not is_string(display): msg = 'Unable to set CMFD display to a non-string ' \ @@ -323,7 +443,8 @@ class CMFDFile(object): self._display = display - def set_feedback(self, feedback): + @feedback.setter + def feedback(self, feedback): if not isinstance(feedback, bool): msg = 'Unable to set CMFD feedback to {0} which is ' \ @@ -333,7 +454,8 @@ class CMFDFile(object): self._feedback = feedback - def set_inactive(self, inactive): + @inactive.setter + def inactive(self, inactive): if not isinstance(inactive, bool): msg = 'Unable to set CMFD inactive batch to {0} which is ' \ @@ -343,7 +465,8 @@ class CMFDFile(object): self._inactive = inactive - def set_inactive_flush(self, inactive_flush): + @inactive_flush.setter + def inactive_flush(self, inactive_flush): if not is_integer(inactive_flush): msg = 'Unable to set CMFD inactive flush batch to {0} which is ' \ @@ -358,7 +481,8 @@ class CMFDFile(object): self._inactive_flush = inactive_flush - def set_ksp_monitor(self, ksp_monitor): + @ksp_monitor.setter + def ksp_monitor(self, ksp_monitor): if not isinstance(ksp_monitor, bool): msg = 'Unable to set CMFD ksp monitor to {0} which is a ' \ @@ -368,7 +492,8 @@ class CMFDFile(object): self._ksp_monitor = ksp_monitor - def set_mesh(self, mesh): + @cmfd_mesh.setter + def cmfd_mesh(self, mesh): if not isinstance(mesh, CMFDMesh): msg = 'Unable to set CMFD mesh to {0} which is not a ' \ @@ -378,7 +503,8 @@ class CMFDFile(object): self._mesh = mesh - def set_norm(self, norm): + @norm.setter + def norm(self, norm): if not is_integer(norm) and not is_float(norm): msg = 'Unable to set the CMFD norm to {0} which is not ' \ @@ -388,8 +514,8 @@ class CMFDFile(object): self._norm = norm - - def set_num_flushes(self, num_flushes): + @snes_monitor.setter + def snum_flushes(self, num_flushes): if not is_integer(num_flushes): msg = 'Unable to set the CMFD number of flushes to {0} ' \ @@ -404,7 +530,8 @@ class CMFDFile(object): self._num_flushes = num_flushes - def set_power_monitor(self, power_monitor): + @power_monitor.setter + def power_monitor(self, power_monitor): if not isinstance(power_monitor, bool): msg = 'Unable to set CMFD power monitor to {0} which is a ' \ @@ -414,7 +541,8 @@ class CMFDFile(object): self._power_monitor = power_monitor - def set_run_adjoint(self, run_adjoint): + @run_adjoint.setter + def run_adjoint(self, run_adjoint): if not isinstance(run_adjoint, bool): msg = 'Unable to set CMFD run adjoint to {0} which is a ' \ @@ -424,7 +552,8 @@ class CMFDFile(object): self._run_adjoint = run_adjoint - def set_snes_monitor(self, snes_monitor): + @snes_monitor.setter + def snes_monitor(self, snes_monitor): if not isinstance(snes_monitor, bool): msg = 'Unable to set CMFD snes monitor to {0} which is a ' \ @@ -434,7 +563,8 @@ class CMFDFile(object): self._snes_monitor = snes_monitor - def set_solver(self, solver): + @solver.setter + def solver(self, solver): if not solver in ['power', 'jfnk']: msg = 'Unable to set CMFD solver to {0} which is not ' \ @@ -444,7 +574,8 @@ class CMFDFile(object): self._solver = solver - def set_write_matrices(self, write_matrices): + @write_matrices.setter + def write_matrices(self, write_matrices): if not isinstance(write_matrices, bool): msg = 'Unable to set CMFD write matrices to {0} which is a ' \ @@ -454,14 +585,14 @@ class CMFDFile(object): self._write_matrices = write_matrices - def create_active_flush_subelement(self): + def get_active_flush_subelement(self): if not self._active_flush is None: element = ET.SubElement(self._cmfd_file, "active_flush") element.text = '{0}'.format(str(self._active_flush)) - def create_begin_subelement(self): + def get_begin_subelement(self): if not self._begin is None: element = ET.SubElement(self._cmfd_file, "begin") @@ -583,4 +714,4 @@ class CMFDFile(object): # Write the XML Tree to the cmfd.xml file tree = ET.ElementTree(self._cmfd_file) tree.write("cmfd.xml", xml_declaration=True, - encoding='utf-8', method="xml") \ No newline at end of file + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/constants.py b/src/utils/openmc/constants.py index d23adfd24d..7c4a7aed41 100644 --- a/src/utils/openmc/constants.py +++ b/src/utils/openmc/constants.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - """Dictionaries of integer-to-string mappings from openmc/src/constants.F90""" SURFACE_TYPES = {1: 'x-plane', @@ -49,16 +47,19 @@ SCORE_TYPES = {-1: 'flux', -8: 'nu-scatter-pn', -9: 'transport', -10: 'n1n', - -11: 'absorption', - -12: 'fission', - -13: 'nu-fission', - -14: 'kappa-fission', - -15: 'current', - -16: 'flux-yn', - -17: 'total-yn', - -18: 'scatter-yn', - -19: 'nu-scatter-yn', - -20: 'events', + -11: 'n2n', + -12: 'n3n', + -13: 'n4n', + -14: 'absorption', + -15: 'fission', + -16: 'nu-fission', + -17: 'kappa-fission', + -18: 'current', + -19: 'flux-yn', + -20: 'total-yn', + -21: 'scatter-yn', + -22: 'nu-scatter-yn', + -23: 'events', 1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)', @@ -122,4 +123,4 @@ SCORE_TYPES.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)}) SCORE_TYPES.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)}) SCORE_TYPES.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)}) SCORE_TYPES.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)}) -SCORE_TYPES.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)}) \ No newline at end of file +SCORE_TYPES.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)}) diff --git a/src/utils/openmc/element.py b/src/utils/openmc/element.py index 745e2ad6bc..b6d5775ee8 100644 --- a/src/utils/openmc/element.py +++ b/src/utils/openmc/element.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.checkvalue import * class Element(object): @@ -11,10 +9,10 @@ class Element(object): self._xs = None # Set the Material class attributes - self.set_name(name) + self.name = name if not xs is None: - self.set_xs(xs) + self.xs = xs def __eq__(self, element2): @@ -36,23 +34,24 @@ class Element(object): def __hash__(self): - hashable = list() + hashable = [] hashable.append(self._name) hashable.append(self._xs) return hash(tuple(hashable)) - def set_name(self, name): - - if not is_string(name): - msg = 'Unable to set name for Element with a non-string ' \ - 'value {0}'.format(name) - raise ValueError(msg) - - self._name = name + @property + def xs(self): + return self._xs - def set_xs(self, xs): + @property + def name(self): + return self._name + + + @xs.setter + def xs(self, xs): if not is_string(xs): msg = 'Unable to set cross-section identifier xs for Element ' \ @@ -62,6 +61,17 @@ class Element(object): self._xs = xs + @name.setter + def name(self, name): + + if not is_string(name): + msg = 'Unable to set name for Element with a non-string ' \ + 'value {0}'.format(name) + raise ValueError(msg) + + self._name = name + + def __repr__(self): string = 'Element - {0}\n'.format(self._name) diff --git a/src/utils/openmc/executor.py b/src/utils/openmc/executor.py index cb09ea8515..8f828ba87b 100644 --- a/src/utils/openmc/executor.py +++ b/src/utils/openmc/executor.py @@ -1,11 +1,7 @@ -#!/usr/bin/env python - -from openmc.checkvalue import * import subprocess import os - -FNULL = open(os.devnull, 'w') +from openmc.checkvalue import * class Executor(object): @@ -15,7 +11,13 @@ class Executor(object): self._working_directory = '.' - def set_working_directory(self, working_directory): + @property + def working_directory(self): + return self._working_directory + + + @working_directory.setter + def working_directory(self, working_directory): if not is_string(working_directory): msg = 'Unable to set Executor\'s working directory to {0} ' \ @@ -36,7 +38,8 @@ class Executor(object): subprocess.check_call('openmc -p', shell=True, cwd=self._working_directory) else: - subprocess.check_call('openmc -p', shell=True, stdout=FNULL, + subprocess.check_call('openmc -p', shell=True, + stdout=open(os.devnull, 'w'), cwd=self._working_directory) @@ -71,5 +74,6 @@ class Executor(object): subprocess.check_call(command, shell=True, cwd=self._working_directory) else: - subprocess.check_call(command, shell=True, stdout=FNULL, - cwd=self._working_directory) \ No newline at end of file + subprocess.check_call(command, shell=True, + stdout=open(os.devnull, 'w'), + cwd=self._working_directory) diff --git a/src/utils/openmc/geometry.py b/src/utils/openmc/geometry.py index 59fb158ee2..84fb6005b0 100644 --- a/src/utils/openmc/geometry.py +++ b/src/utils/openmc/geometry.py @@ -1,8 +1,8 @@ -#!/usr/bin/env python +from xml.etree import ElementTree as ET import openmc from openmc.clean_xml import * -from xml.etree import ElementTree as ET + def reset_auto_ids(): openmc.reset_auto_material_id() @@ -17,7 +17,29 @@ class Geometry(object): # Initialize Geometry class attributes self._root_universe = None - self._offsets = dict() + self._offsets = {} + + + @property + def root_universe(self): + return self._root_universe + + + @root_universe.setter + def root_universe(self, root_universe): + + if not isinstance(root_universe, openmc.Universe): + msg = 'Unable to add root Universe {0} to Geometry since ' \ + 'it is not a Universe'.format(root_universe) + raise ValueError(msg) + + elif root_universe._id != 0: + msg = 'Unable to add root Universe {0} to Geometry since ' \ + 'it has ID={1} instead of ' \ + 'ID=0'.format(root_universe, root_universe._id) + raise ValueError(msg) + + self._root_universe = root_universe def get_offset(self, path, filter_offset): @@ -63,7 +85,7 @@ class Geometry(object): def get_all_nuclides(self): - nuclides = dict() + nuclides = {} materials = self.get_all_materials() for material in materials: @@ -111,22 +133,6 @@ class Geometry(object): return list(material_universes) - def set_root_universe(self, root_universe): - - if not isinstance(root_universe, openmc.Universe): - msg = 'Unable to add root Universe {0} to Geometry since ' \ - 'it is not a Universe'.format(root_universe) - raise ValueError(msg) - - elif root_universe._id != 0: - msg = 'Unable to add root Universe {0} to Geometry since ' \ - 'it has ID={1} instead of ' \ - 'ID=0'.format(root_universe, root_universe._id) - raise ValueError(msg) - - self._root_universe = root_universe - - class GeometryFile(object): @@ -137,7 +143,13 @@ class GeometryFile(object): self._geometry_file = ET.Element("geometry") - def set_geometry(self, geometry): + @property + def geometry(self): + return self._geometry + + + @geometry.setter + def geometry(self, geometry): if not isinstance(geometry, Geometry): msg = 'Unable to set the Geometry to {0} for the GeometryFile ' \ @@ -159,4 +171,4 @@ class GeometryFile(object): # Write the XML Tree to the materials.xml file tree = ET.ElementTree(self._geometry_file) tree.write("geometry.xml", xml_declaration=True, - encoding='utf-8', method="xml") \ No newline at end of file + encoding='utf-8', method="xml") diff --git a/src/utils/openmc/material.py b/src/utils/openmc/material.py index 92c78bfdb1..7e066bfce8 100644 --- a/src/utils/openmc/material.py +++ b/src/utils/openmc/material.py @@ -1,18 +1,15 @@ -#!/usr/bin/env python - +from collections import MappingView +from copy import deepcopy import warnings +from xml.etree import ElementTree as ET import openmc from openmc.checkvalue import * from openmc.clean_xml import * -from xml.etree import ElementTree as ET -from collections import MappingView -from copy import deepcopy -import numpy as np # A list of all IDs for all Materials created -MATERIAL_IDS = list() +MATERIAL_IDS = [] # A static variable for auto-generated Material IDs AUTO_MATERIAL_ID = 10000 @@ -20,7 +17,7 @@ AUTO_MATERIAL_ID = 10000 def reset_auto_material_id(): global AUTO_MATERIAL_ID, MATERIAL_IDS AUTO_MATERIAL_ID = 10000 - MATERIAL_IDS = list() + MATERIAL_IDS = [] # Units for density supported by OpenMC @@ -41,23 +38,23 @@ class Material(object): def __init__(self, material_id=None, name=''): # Initialize class attributes - self._id = None - self._name = '' + self.id = material_id + self.name = name self._density = None self._density_units = '' # A dictionary of Nuclides # Keys - Nuclide names # Values - tuple (nuclide, percent, percent type) - self._nuclides = dict() + self._nuclides = {} # A dictionary of Elements # Keys - Element names # Values - tuple (element, percent, percent type) - self._elements = dict() + self._elements = {} # If specified, a list of tuples of (table name, xs identifier) - self._sab = list() + self._sab = [] # If true, the material will be initialized as distributed self._convert_to_distrib_comps = False @@ -65,21 +62,47 @@ class Material(object): # If specified, this file will be used instead of composition values self._distrib_otf_file = None - # Set the Material class attributes - self.set_id(material_id) - self.set_name(name) + + @property + def id(self): + return self._id - def set_id(self, material_id=None): + @property + def name(self): + return self._name - global MATERIAL_IDS + + @property + def density(self): + return self._density + + + @property + def density_units(self): + return self._density_units + + + @property + def convert_to_distrib_comps(self): + return self._convert_to_distrib_comps + + + @property + def distrib_otf_file(self): + return self._distrib_otf_file + + + @id.setter + def id(self, material_id): + + global AUTO_MATERIAL_ID, MATERIAL_IDS # If the Material already has an ID, remove it from global list - if not self._id is None: + if hasattr(self, '_id') and not self._id is None: MATERIAL_IDS.remove(self._id) if material_id is None: - global AUTO_MATERIAL_ID self._id = AUTO_MATERIAL_ID MATERIAL_IDS.append(AUTO_MATERIAL_ID) AUTO_MATERIAL_ID += 1 @@ -105,7 +128,8 @@ class Material(object): MATERIAL_IDS.append(material_id) - def set_name(self, name): + @name.setter + def name(self, name): if not is_string(name): msg = 'Unable to set name for Material ID={0} with a non-string ' \ @@ -138,6 +162,31 @@ class Material(object): self._density_units = units + @distrib_otf_file.setter + def distrib_otf_file(self, filename): + + # TODO: remove this when distributed materials are merged + warnings.warn('This feature is not yet implemented in a release ' \ + 'version of openmc') + + if not is_string(filename) and not filename is None: + msg = 'Unable to add OTF material file to Material ID={0} with a ' \ + 'non-string name {1}'.format(self._id, filename) + raise ValueError(msg) + + self._distrib_otf_file = filename + + + @convert_to_distrib_comps.setter + def convert_to_distrib_comps(self): + + # TODO: remove this when distributed materials are merged + warnings.warn('This feature is not yet implemented in a release ' \ + 'version of openmc') + + self._convert_to_distrib_comps = True + + def add_nuclide(self, nuclide, percent, percent_type='ao'): if not isinstance(nuclide, openmc.Nuclide): @@ -218,32 +267,9 @@ class Material(object): self._sab.append((name, xs)) - def set_otf_mat_file(self, name): - - # TODO: remove this when distributed materials are merged - warnings.warn('This feature is not yet implemented in a release ' \ - 'version of openmc') - - if not is_string(name): - msg = 'Unable to add OTF material file to Material ID={0} with a ' \ - 'non-string name {1}'.format(self._id, name) - raise ValueError(msg) - - self._distrib_otf_file = name - - - def set_as_distrib_comp(self): - - # TODO: remove this when distributed materials are merged - warnings.warn('This feature is not yet implemented in a release ' \ - 'version of openmc') - - self._convert_to_distrib_comps = True - - def get_all_nuclides(self): - nuclides = dict() + nuclides = {} for nuclide_name, nuclide_tuple in self._nuclides.items(): nuclide = nuclide_tuple[0] @@ -262,13 +288,13 @@ class Material(object): string += '{0: <16}{1}{2}'.format('\tDensity', '=\t', self._density) string += ' [{0}]\n'.format(self._density_units) - string += '{0: <16}'.format('\tS(a,b) Tables') + '\n' + string += '{0: <16}\n'.format('\tS(a,b) Tables') for sab in self._sab: string += '{0: <16}{1}[{2}{3}]\n'.format('\tS(a,b)', '=\t', sab[0], sab[1]) - string += '{0: <16}'.format('\tNuclides') + '\n' + string += '{0: <16}\n'.format('\tNuclides') for nuclide in self._nuclides: percent = self._nuclides[nuclide][1] @@ -320,7 +346,7 @@ class Material(object): def get_nuclides_xml(self, nuclides, distrib=False): - xml_elements = list() + xml_elements = [] for nuclide in nuclides.values(): xml_elements.append(self.get_nuclide_xml(nuclide, distrib)) @@ -330,7 +356,7 @@ class Material(object): def get_elements_xml(self, elements, distrib=False): - xml_elements = list() + xml_elements = [] for element in elements.values(): xml_elements.append(self.get_element_xml(element, distrib)) @@ -368,7 +394,7 @@ class Material(object): else: subelement = ET.SubElement(element, "compositions") - + comps = [] allnucs = self._nuclides.values() + self._elements.values() dist_per_type = allnucs[0][2] @@ -419,11 +445,26 @@ class MaterialsFile(object): def __init__(self): # Initialize MaterialsFile class attributes - self._materials = list() + self._materials = [] self._default_xs = None self._materials_file = ET.Element("materials") + @property + def default_xs(self): + return self._default_xs + + + @default_xs.setter + def default_xs(self, xs): + + if not is_string(xs): + msg = 'Unable to set default xs to a non-string value'.format(xs) + raise ValueError(msg) + + self._default_xs = xs + + def add_material(self, material): if not isinstance(material, Material): @@ -455,15 +496,6 @@ class MaterialsFile(object): self._materials.remove(material) - def set_default_xs(self, xs): - - if not is_string(xs): - msg = 'Unable to set default xs to a non-string value'.format(xs) - raise ValueError(msg) - - self._default_xs = xs - - def create_material_subelements(self): subelement = ET.SubElement(self._materials_file, "default_xs") diff --git a/src/utils/openmc/nuclide.py b/src/utils/openmc/nuclide.py index 8603c833b8..c2287b23d9 100644 --- a/src/utils/openmc/nuclide.py +++ b/src/utils/openmc/nuclide.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - from openmc.checkvalue import * @@ -13,10 +11,10 @@ class Nuclide(object): self._zaid = None # Set the Material class attributes - self.set_name(name) + self.name = name if not xs is None: - self.set_xs(xs) + self.xs = xs def __eq__(self, nuclide2): @@ -38,13 +36,26 @@ class Nuclide(object): def __hash__(self): - hashable = list() - hashable.append(self._name) - hashable.append(self._xs) - return hash(tuple(hashable)) + return hash((self._name, self._xs)) - def set_name(self, name): + @property + def name(self): + return self._name + + + @property + def xs(self): + return self._xs + + + @property + def zaid(self): + return self._zaid + + + @name.setter + def name(self, name): if not is_string(name): msg = 'Unable to set name for Nuclide with a non-string ' \ @@ -54,7 +65,8 @@ class Nuclide(object): self._name = name - def set_xs(self, xs): + @xs.setter + def xs(self, xs): if not is_string(xs): msg = 'Unable to set cross-section identifier xs for Nuclide ' \ @@ -64,7 +76,8 @@ class Nuclide(object): self._xs = xs - def set_zaid(self, zaid): + @zaid.setter + def zaid(self, zaid): if not is_integer(zaid): msg = 'Unable to set zaid for Nuclide ' \ diff --git a/src/utils/openmc/opencg_compatible.py b/src/utils/openmc/opencg_compatible.py index 1c65d285cf..320529f68a 100644 --- a/src/utils/openmc/opencg_compatible.py +++ b/src/utils/openmc/opencg_compatible.py @@ -1,60 +1,60 @@ -#!/usr/bin/env python +import copy + +import numpy as np +import opencg import openmc -import opencg -import copy -import numpy as np # A dictionary of all OpenMC Materials created # Keys - Material IDs # Values - Materials -OPENMC_MATERIALS = dict() +OPENMC_MATERIALS = {} # A dictionary of all OpenCG Materials created # Keys - Material IDs # Values - Materials -OPENCG_MATERIALS = dict() +OPENCG_MATERIALS = {} # A dictionary of all OpenMC Surfaces created # Keys - Surface IDs # Values - Surfaces -OPENMC_SURFACES = dict() +OPENMC_SURFACES = {} # A dictionary of all OpenCG Surfaces created # Keys - Surface IDs # Values - Surfaces -OPENCG_SURFACES = dict() +OPENCG_SURFACES = {} # A dictionary of all OpenMC Cells created # Keys - Cell IDs # Values - Cells -OPENMC_CELLS = dict() +OPENMC_CELLS = {} # A dictionary of all OpenCG Cells created # Keys - Cell IDs # Values - Cells -OPENCG_CELLS = dict() +OPENCG_CELLS = {} # A dictionary of all OpenMC Universes created # Keys - Universes IDs # Values - Universes -OPENMC_UNIVERSES = dict() +OPENMC_UNIVERSES = {} # A dictionary of all OpenCG Universes created # Keys - Universes IDs # Values - Universes -OPENCG_UNIVERSES = dict() +OPENCG_UNIVERSES = {} # A dictionary of all OpenMC Lattices created # Keys - Lattice IDs # Values - Lattices -OPENMC_LATTICES = dict() +OPENMC_LATTICES = {} # A dictionary of all OpenCG Lattices created # Keys - Lattice IDs # Values - Lattices -OPENCG_LATTICES = dict() +OPENCG_LATTICES = {} @@ -144,7 +144,7 @@ def get_opencg_surface(openmc_surface): name = openmc_surface._name # Correct for OpenMC's syntax for Surfaces dividing Cells - boundary = openmc_surface._bc_type + boundary = openmc_surface._boundary_type if boundary == 'transmission': boundary = 'interface' @@ -410,7 +410,7 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace): raise ValueError(msg) # Initialize an empty list for the new compatible cells - compatible_cells = list() + compatible_cells = [] # SquarePrism Surfaces if opencg_surface._type in ['x-squareprism', @@ -562,15 +562,15 @@ def get_openmc_cell(opencg_cell): fill = opencg_cell._fill if (opencg_cell._type == 'universe'): - openmc_cell.set_fill(get_openmc_universe(fill)) + openmc_cell.fill = get_openmc_universe(fill) elif (opencg_cell._type == 'lattice'): - openmc_cell.set_fill(get_openmc_lattice(fill)) + openmc_cell.fill = get_openmc_lattice(fill) else: - openmc_cell.set_fill(get_openmc_material(fill)) + openmc_cell.fill = get_openmc_material(fill) if opencg_cell._rotation: rotation = np.asarray(opencg_cell._rotation, dtype=np.int) - openmc_cell.set_rotation(rotation) + openmc_cell.rotation = rotation if opencg_cell._translation: translation = np.asarray(opencg_cell._translation, dtype=np.float64) @@ -679,11 +679,21 @@ def get_opencg_lattice(openmc_lattice): return OPENCG_LATTICES[lattice_id] # Create an OpenCG Lattice to represent this OpenMC Lattice - name = openmc_lattice._name - dimension = openmc_lattice._dimension - pitch = openmc_lattice._pitch - lower_left = openmc_lattice._lower_left - universes = openmc_lattice._universes + name = openmc_lattice.name + dimension = openmc_lattice.dimension + pitch = openmc_lattice.pitch + lower_left = openmc_lattice.lower_left + universes = openmc_lattice.universes + + if len(pitch) == 2: + new_pitch = np.ones(3, dtype=np.float64) + new_pitch[:2] = pitch + pitch = new_pitch + + if len(lower_left) == 2: + new_lower_left = np.ones(3, dtype=np.float64) + new_lower_left[:2] = lower_left + lower_left = new_lower_left # Initialize an empty array for the OpenCG nested Universes in this Lattice universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), \ @@ -765,10 +775,10 @@ def get_openmc_lattice(opencg_lattice): np.array(dimension, dtype=np.float64))) / -2.0 openmc_lattice = openmc.RectLattice(lattice_id=lattice_id) - openmc_lattice.set_dimension(dimension) - openmc_lattice.set_pitch(width) - openmc_lattice.set_universes(universe_array) - openmc_lattice.set_lower_left(lower_left) + openmc_lattice.dimension = dimension + openmc_lattice.pitch = width + openmc_lattice.universes = universe_array + openmc_lattice.lower_left = lower_left # Add the OpenMC Lattice to the global collection of all OpenMC Lattices OPENMC_LATTICES[lattice_id] = openmc_lattice @@ -815,6 +825,7 @@ def get_openmc_geometry(opencg_geometry): # Deep copy the goemetry since it may be modified to make all Surfaces # compatible with OpenMC's specifications + opencg_geometry.assignAutoIds() opencg_geometry = copy.deepcopy(opencg_geometry) # Update Cell bounding boxes in Geometry @@ -842,6 +853,6 @@ def get_openmc_geometry(opencg_geometry): openmc_root_universe = get_openmc_universe(opencg_root_universe) openmc_geometry = openmc.Geometry() - openmc_geometry.set_root_universe(openmc_root_universe) + openmc_geometry.root_universe = openmc_root_universe return openmc_geometry diff --git a/src/utils/openmc/plots.py b/src/utils/openmc/plots.py index 551b2e76f4..1468b11b0d 100644 --- a/src/utils/openmc/plots.py +++ b/src/utils/openmc/plots.py @@ -1,9 +1,9 @@ -#!/usr/bin/env python +from xml.etree import ElementTree as ET + +import numpy as np from openmc.checkvalue import * from openmc.clean_xml import * -from xml.etree import ElementTree as ET -import numpy as np # A static variable for auto-generated Plot IDs @@ -22,8 +22,8 @@ class Plot(object): def __init__(self, plot_id=None, name=''): # Initialize Plot class attributes - self._id = None - self._name = '' + self.id = plot_id + self.name = name self._width = [4.0, 4.0] self._pixels = [1000, 1000] self._origin = [0., 0., 0.] @@ -36,11 +36,74 @@ class Plot(object): self._mask_background = None self._col_spec = None - self.set_id(plot_id) - self.set_name(name) + + @property + def id(self): + return self._id - def set_id(self, plot_id=None): + @property + def name(self): + return self._name + + + @property + def width(self): + return self._width + + + @property + def pixels(self): + return self._pixels + + + @property + def origin(self): + return self._origin + + + @property + def filename(self): + return self._filename + + + @property + def color(self): + return self._color + + + @property + def type(self): + return self._type + + + @property + def basis(self): + return self._basis + + + @property + def background(self): + return self._background + + + @property + def mask_componenets(self): + return self._mask_components + + + @property + def mask_background(self): + return self._mask_background + + + @property + def col_spec(self): + return self._col_spec + + + @id.setter + def id(self, plot_id): if plot_id is None: global AUTO_PLOT_ID @@ -54,25 +117,27 @@ class Plot(object): elif plot_id < 0: msg = 'Unable to set Plot ID to {0} since it must be a ' \ - 'non-negative integer'.format(plot_id) + 'non-negative integer'.format(plot_id) raise ValueError(msg) else: self._id = plot_id - def set_name(self, name): + @name.setter + def name(self, name): if not is_string(name): msg = 'Unable to set name for Plot ID={0} with a non-string ' \ - 'value {1}'.format(self._id, name) + 'value {1}'.format(self._id, name) raise ValueError(msg) else: self._name = name - def set_width(self, width): + @width.setter + def width(self, width): if not isinstance(width, (tuple, list, np.ndarray)): msg = 'Unable to create Plot ID={0} with width {1} which is not ' \ @@ -94,7 +159,8 @@ class Plot(object): self._width = width - def set_origin(self, origin): + @origin.setter + def origin(self, origin): if not isinstance(origin, (tuple, list, np.ndarray)): msg = 'Unable to create Plot ID={0} with origin {1} which is not ' \ @@ -117,7 +183,8 @@ class Plot(object): self._origin = origin - def set_pixels(self, pixels): + @pixels.setter + def pixels(self, pixels): if not isinstance(pixels, (tuple, list, np.ndarray)): msg = 'Unable to create Plot ID={0} with pixels {1} which is not ' \ @@ -144,7 +211,8 @@ class Plot(object): self._pixels = pixels - def set_filename(self, filename): + @filename.setter + def filename(self, filename): if not is_string(filename): msg = 'Unable to create Plot ID={0} with filename {1} which is ' \ @@ -154,7 +222,8 @@ class Plot(object): self._filename = filename - def set_color(self, color): + @color.setter + def color(self, color): if not is_string(color): msg = 'Unable to create Plot ID={0} with color {1} which is not ' \ @@ -169,7 +238,8 @@ class Plot(object): self._color = color - def set_type(self, type): + @type.setter + def type(self, type): if not is_string(type): msg = 'Unable to create Plot ID={0} with type {1} which is not ' \ @@ -184,7 +254,8 @@ class Plot(object): self._type = type - def set_basis(self, basis): + @basis.setter + def basis(self, basis): if not is_string(basis): msg = 'Unable to create Plot ID={0} with basis {1} which is not ' \ @@ -199,7 +270,8 @@ class Plot(object): self._basis = basis - def set_background(self, background): + @background.setter + def background(self, background): if not isinstance(background, (tuple, list, np.ndarray)): msg = 'Unable to create Plot ID={0} with background {1} ' \ @@ -228,7 +300,8 @@ class Plot(object): self._background = background - def set_col_spec(self, col_spec): + @col_spec.setter + def col_spec(self, col_spec): if not isinstance(col_spec, dict): msg= 'Unable to create Plot ID={0} with col_spec parameter {1} ' \ @@ -263,7 +336,8 @@ class Plot(object): self._col_spec = col_spec - def set_mask_components(self, mask_components): + @mask_componenets.setter + def mask_components(self, mask_components): if not isinstance(mask_components, (list, tuple, np.ndarray)): msg = 'Unable to create Plot ID={0} with mask components {1} ' \ @@ -285,7 +359,8 @@ class Plot(object): self._mask_components = mask_components - def set_mask_background(self, mask_background): + @mask_background.setter + def mask_background(self, mask_background): if not isinstance(mask_background, (list, tuple, np.ndarray)): msg = 'Unable to create Plot ID={0} with mask background {1} ' \ @@ -400,7 +475,7 @@ class PlotsFile(object): def __init__(self): # Initialize PlotsFile class attributes - self._plots = list() + self._plots = [] self._plots_file = ET.Element("plots") diff --git a/src/utils/openmc/settings.py b/src/utils/openmc/settings.py index 6acb720cbe..ce054487bc 100644 --- a/src/utils/openmc/settings.py +++ b/src/utils/openmc/settings.py @@ -1,12 +1,11 @@ -#!/usr/bin/env python - +import collections import warnings +from xml.etree import ElementTree as ET + +import numpy as np from openmc.checkvalue import * from openmc.clean_xml import * -from xml.etree import ElementTree as ET -import numpy as np -import multiprocessing class SettingsFile(object): @@ -77,13 +76,244 @@ class SettingsFile(object): self._dd_mesh_upper_right = None self._dd_nodemap = None self._dd_allow_leakage = False + self._dd_count_interactions = False self._settings_file = ET.Element("settings") self._eigenvalue_element = None self._source_element = None - def set_batches(self, batches): + @property + def batches(self): + return self._batches + + + @property + def generations_per_batch(self): + return self._generations_per_batch + + + @property + def inactive(self): + return self._inactive + + + @property + def particles(self): + return self._particles + + + @property + def source_file(self): + return self._source_file + + + @property + def source_space_type(self): + return self._source_space_type + + + @property + def source_space_params(self): + return self._source_space_params + + + @property + def source_angle_type(self): + return self._source_angle_type + + + @property + def source_angle_params(self): + return self._source_angle_params + + + @property + def source_energy_type(self): + return self._source_energy_type + + + @property + def source_energy_params(self): + return self._source_energy_params + + + @property + def confidence_intervals(self): + return self._confidence_intervals + + + @property + def cross_sections(self): + return self._cross_sections + + + @property + def energy_grid(self): + return self._energy_grid + + + @property + def ptables(self): + return self._ptables + + + @property + def run_cmfd(self): + return self._run_cmfd + + + @property + def seed(self): + return self._seed + + + @property + def survival_biasing(self): + return self._survival_biasing + + + @property + def entropy_dimension(self): + return self._entropy_dimension + + + @property + def entropy_lower_left(self): + return self._entropy_lower_left + + + @property + def entropy_upper_right(self): + return self._entropy_upper_right + + + @property + def output(self): + return self._output + + + @property + def output_path(self): + return self._output_path + + + @property + def statepoint_batches(self): + return self._statepoint_batches + + + @property + def statepoint_interval(self): + return self._statepoint_interval + + @property + def sourcepoint_batches(self): + return self._sourcepoint_interval + + + @property + def sourcepoint_interval(self): + return self._sourcepoint_interval + + + @property + def sourcepoint_separate(self): + return self._sourcepoint_separate + + + @property + def sourcepoint_write(self): + return self._sourcepoint_write + + + @property + def sourcepoint_overwrite(self): + return self._sourcepoint_overwrite + + + @property + def threads(self): + return self._threads + + + @property + def no_reduce(self): + return self._no_reduce + + + @property + def verbosity(self): + return self._verbosity + + + @property + def trace(self): + return self._trace + + + @property + def track(self): + return self._track + + + @property + def weight(self): + return self._weight + + + @property + def weight_avg(self): + return self._weight_avg + + + @property + def ufs_dimension(self): + return self._ufs_dimension + + + @property + def ufs_lower_left(self): + return self._ufs_lower_left + + + @property + def ufs_upper_right(self): + return self._ufs_upper_right + + + @property + def dd_mesh_dimension(self): + return self._dd_mesh_dimension + + + @property + def dd_mesh_lower_left(self): + return self._dd_mesh_lower_left + + + @property + def dd_mesh_upper_right(self): + return self._dd_mesh_upper_right + + + @property + def dd_nodemap(self): + return self._dd_nodemap + + + @property + def dd_allow_leakage(self): + return self._dd_allow_leakage + + + @property + def dd_count_interactions(self): + return self._dd_count_interactions + + + @batches.setter + def batches(self, batches): if not is_integer(batches): msg = 'Unable to set batches to a non-integer ' \ @@ -98,7 +328,8 @@ class SettingsFile(object): self._batches = batches - def set_generations_per_batch(self, generations_per_batch): + @generations_per_batch.setter + def generations_per_batch(self, generations_per_batch): if not is_integer(generations_per_batch): msg = 'Unable to set generations per batch to a non-integer ' \ @@ -113,7 +344,8 @@ class SettingsFile(object): self._generations_per_batch = generations_per_batch - def set_inactive(self, inactive): + @inactive.setter + def inactive(self, inactive): if not is_integer(inactive): msg = 'Unable to set inactive batches to a non-integer ' \ @@ -128,7 +360,8 @@ class SettingsFile(object): self._inactive = inactive - def set_particles(self, particles): + @particles.setter + def particles(self, particles): if not is_integer(particles): msg = 'Unable to set particles to a non-integer ' \ @@ -143,7 +376,8 @@ class SettingsFile(object): self._particles = particles - def set_source_file(self, source_file): + @source_file.setter + def source_file(self, source_file): if not is_string(source_file): msg = 'Unable to set source file to a non-string ' \ @@ -153,16 +387,16 @@ class SettingsFile(object): self._source_file = source_file - def set_source_space(self, type, params): + def set_source_space(self, stype, params): - if not is_string(type): + if not is_string(stype): msg = 'Unable to set source space type to a non-string ' \ - 'value {0}'.format(type) + 'value {0}'.format(stype) raise ValueError(msg) - elif not type in ['box', 'point']: + elif not stype in ['box', 'point']: msg = 'Unable to set source space type to {0} since it is not ' \ - 'box or point'.format(type) + 'box or point'.format(stype) raise ValueError(msg) elif not isinstance(params, (tuple, list, np.ndarray)): @@ -182,20 +416,20 @@ class SettingsFile(object): 'is not an integer or floating point value'.format(param) raise ValueError(msg) - self._source_space_type = type + self._source_space_type = stype self._source_space_params = params - def set_source_angle(self, type, params=[]): + def set_source_angle(self, stype, params=[]): - if not is_string(type): + if not is_string(stype): msg = 'Unable to set source angle type to a non-string ' \ - 'value {0}'.format(type) + 'value {0}'.format(stype) raise ValueError(msg) - elif not type in ['isotropic', 'monodirectional']: + elif not stype in ['isotropic', 'monodirectional']: msg = 'Unable to set source angle type to {0} since it is not ' \ - 'isotropic or monodirectional'.format(type) + 'isotropic or monodirectional'.format(stype) raise ValueError(msg) elif not isinstance(params, (tuple, list, np.ndarray)): @@ -203,12 +437,12 @@ class SettingsFile(object): 'not a Python list/tuple or NumPy array'.format(params) raise ValueError(msg) - elif type is 'isotropic' and not params is None: + elif stype == 'isotropic' and not params is None: msg = 'Unable to set source angle parameters since they are not ' \ 'it is not supported for isotropic type sources' raise ValueError(msg) - elif type is 'monodirectional' and len(params) != 3: + elif stype == 'monodirectional' and len(params) != 3: msg = 'Unable to set source angle parameters to {0} ' \ 'since 3 parameters are required for monodirectional ' \ 'sources'.format(params) @@ -221,20 +455,20 @@ class SettingsFile(object): 'is not an integer or floating point value'.format(param) raise ValueError(msg) - self._source_angle_type = type + self._source_angle_type = stype self._source_angle_params = params - def set_source_energy(self, type, params=[]): + def set_source_energy(self, stype, params=[]): - if not is_string(type): + if not is_string(stype): msg = 'Unable to set source energy type to a non-string ' \ - 'value {0}'.format(type) + 'value {0}'.format(stype) raise ValueError(msg) - elif not type in ['monoenergetic', 'watt', 'maxwell']: + elif not stype in ['monoenergetic', 'watt', 'maxwell']: msg = 'Unable to set source energy type to {0} since it is not ' \ - 'monoenergetic, watt or maxwell'.format(type) + 'monoenergetic, watt or maxwell'.format(stype) raise ValueError(msg) elif not isinstance(params, (tuple, list, np.ndarray)): @@ -242,19 +476,19 @@ class SettingsFile(object): 'is not a Python list/tuple or NumPy array'.format(params) raise ValueError(msg) - elif type is 'monoenergetic' and not len(params) != 1: + elif stype == 'monoenergetic' and not len(params) != 1: msg = 'Unable to set source energy parameters to {0} ' \ 'since 1 paramater is required for monenergetic ' \ 'sources'.format(params) raise ValueError(msg) - elif type is 'watt' and len(params) != 2: + elif stype == 'watt' and len(params) != 2: msg = 'Unable to set source energy parameters to {0} ' \ 'since 2 parameters are required for monoenergetic ' \ 'sources'.format(params) raise ValueError(msg) - elif type is 'maxwell' and len(params) != 2: + elif stype == 'maxwell' and len(params) != 2: msg = 'Unable to set source energy parameters to {0} since 1 ' \ 'parameter is required for maxwell sources'.format(params) raise ValueError(msg) @@ -267,11 +501,12 @@ class SettingsFile(object): 'value'.format(param) raise ValueError(msg) - self._source_energy_type = type + self._source_energy_type = stype self._source_energy_params = params - def set_output(self, output): + @output.setter + def output(self, output): if not isinstance(output, dict): msg = 'Unable to set output to {0} which is not a Python ' \ @@ -295,7 +530,8 @@ class SettingsFile(object): self._output = output - def set_output_path(self, output_path): + @output_path.setter + def output_path(self, output_path): if not is_string(output_path): msg = 'Unable to set output path to non-string ' \ @@ -305,7 +541,8 @@ class SettingsFile(object): self._output_path = output_path - def set_verbosity(self, verbosity): + @verbosity.setter + def verbosity(self, verbosity): if not is_integer(verbosity): msg = 'Unable to set verbosity to non-integer ' \ @@ -320,7 +557,8 @@ class SettingsFile(object): self._verbosity = verbosity - def set_statepoint_batches(self, batches): + @statepoint_batches.setter + def statepoint_batches(self, batches): if not isinstance(batches, (tuple, list, np.ndarray)): msg = 'Unable to set statepoint batches to {0} which is not a ' \ @@ -342,7 +580,8 @@ class SettingsFile(object): self._statepoint_batches = batches - def set_statepoint_interval(self, interval): + @statepoint_interval.setter + def statepoint_interval(self, interval): if not is_integer(interval): msg = 'Unable to set statepoint interval to non-integer ' \ @@ -352,7 +591,8 @@ class SettingsFile(object): self._statepoint_interval = interval - def set_sourcepoint_batches(self, batches): + @sourcepoint_batches.setter + def sourcepoint_batches(self, batches): if not isinstance(batches, (tuple, list, np.ndarray)): msg = 'Unable to set sourcepoint batches to {0} which is ' \ @@ -374,7 +614,8 @@ class SettingsFile(object): self._sourcepoint_batches = batches - def set_sourcepoint_interval(self, interval): + @sourcepoint_interval.setter + def sourcepoint_interval(self, interval): if not is_integer(interval): msg = 'Unable to set sourcepoint interval to non-integer ' \ @@ -384,7 +625,8 @@ class SettingsFile(object): self._sourcepoint_interval = interval - def set_sourcepoint_separate(self, source_separate): + @sourcepoint_separate.setter + def sourcepoint_separate(self, source_separate): if not isinstance(source_separate, (bool, np.bool)): msg = 'Unable to set sourcepoint separate to non-boolean ' \ @@ -394,7 +636,8 @@ class SettingsFile(object): self._sourcepoint_separate = source_separate - def set_sourcepoint_write(self, source_write): + @sourcepoint_write.setter + def sourcepoint_write(self, source_write): if not isinstance(source_write, (bool, np.bool)): msg = 'Unable to set sourcepoint write to non-boolean ' \ @@ -404,7 +647,8 @@ class SettingsFile(object): self._sourcepoint_write = source_write - def set_sourcepoint_overwrite(self, source_overwrite): + @sourcepoint_overwrite.setter + def sourcepoint_overwrite(self, source_overwrite): if not isinstance(source_overwrite, (bool, np.bool)): msg = 'Unable to set sourcepoint overwrite to non-boolean ' \ @@ -414,7 +658,8 @@ class SettingsFile(object): self._sourcepoint_overwrite = source_overwrite - def set_confidence_intervals(self, confidence_intervals): + @confidence_intervals.setter + def confidence_intervals(self, confidence_intervals): if not isinstance(confidence_intervals, (bool, np.bool)): msg = 'Unable to set confidence interval to non-boolean ' \ @@ -424,7 +669,8 @@ class SettingsFile(object): self._confidence_intervals = confidence_intervals - def set_cross_sections(self, cross_sections): + @cross_sections.setter + def cross_sections(self, cross_sections): if not is_string(cross_sections): msg = 'Unable to set cross sections to non-string ' \ @@ -434,17 +680,19 @@ class SettingsFile(object): self._cross_sections = cross_sections - def set_energy_grid(self, energy_grid): + @energy_grid.setter + def energy_grid(self, energy_grid): - if not energy_grid in ['union', 'nuclide']: + if not energy_grid in ['nuclide', 'logarithm', 'material-union']: msg = 'Unable to set energy grid to {0} which is neither ' \ - 'union nor nuclide'.format(energy_grid) + 'nuclide, logarithm, nor material-union'.format(energy_grid) raise ValueError(msg) self._energy_grid = energy_grid - def set_ptables(self, ptables): + @ptables.setter + def ptables(self, ptables): if not isinstance(ptables, (bool, np.bool)): msg = 'Unable to set ptables to non-boolean ' \ @@ -454,7 +702,8 @@ class SettingsFile(object): self._ptables = ptables - def set_run_cmfd(self, run_cmfd): + @run_cmfd.setter + def run_cmfd(self, run_cmfd): if not isinstance(run_cmfd, (bool, np.bool)): msg = 'Unable to set run_cmfd to non-boolean ' \ @@ -464,7 +713,8 @@ class SettingsFile(object): self._run_cmfd = run_cmfd - def set_seed(self, seed): + @seed.setter + def seed(self, seed): if not is_integer(seed): msg = 'Unable to set seed to non-integer value {0}'.format(seed) @@ -477,7 +727,8 @@ class SettingsFile(object): self._seed = seed - def set_survival_biasing(self, survival_biasing): + @survival_biasing.setter + def survival_biasing(self, survival_biasing): if not isinstance(survival_biasing, (bool, np.bool)): msg = 'Unable to set survival biasing to non-boolean ' \ @@ -487,7 +738,8 @@ class SettingsFile(object): self._survival_biasing = survival_biasing - def set_weight(self, weight, weight_avg): + @weight.setter + def weight(self, weight, weight_avg): if not is_float(weight): msg = 'Unable to set weight cutoff to non-floating point ' \ @@ -513,7 +765,8 @@ class SettingsFile(object): self._weight_avg = weight_avg - def set_entropy_dimension(self, dimension): + @entropy_dimension.setter + def entropy_dimension(self, dimension): if not isinstance(dimension, (tuple, list)): msg = 'Unable to set entropy mesh dimension to {0} which is ' \ @@ -535,7 +788,8 @@ class SettingsFile(object): self._entropy_dimension = dimension - def set_entropy_lower_left(self, lower_left): + @entropy_lower_left.setter + def entropy_lower_left(self, lower_left): if not isinstance(lower_left, (tuple, list)): msg = 'Unable to set entropy mesh lower left corner to {0} which ' \ @@ -558,7 +812,8 @@ class SettingsFile(object): self._entropy_lower_left = lower_left - def set_entropy_upper_right(self, upper_right): + @entropy_upper_right.setter + def entropy_upper_right(self, upper_right): if not isinstance(upper_right, (tuple, list)): msg = 'Unable to set entropy mesh upper right corner to {0} ' \ @@ -581,7 +836,8 @@ class SettingsFile(object): self._entropy_upper_right = upper_right - def set_no_reduce(self, no_reduce): + @no_reduce.setter + def no_reduce(self, no_reduce): if not isinstance(no_reduce, (bool, np.bool)): msg = 'Unable to set the no_reduce to a non-boolean ' \ @@ -591,7 +847,8 @@ class SettingsFile(object): self._no_reduce = no_reduce - def set_threads(self, threads): + @threads.setter + def threads(self, threads): if not is_integer(threads): msg = 'Unable to set the threads to a non-integer ' \ @@ -606,7 +863,8 @@ class SettingsFile(object): self._threads = threads - def set_trace(self, trace): + @trace.setter + def trace(self, trace): if not isinstance(trace, (list, tuple)): msg = 'Unable to set the trace to {0} which is not a Python ' \ @@ -636,7 +894,8 @@ class SettingsFile(object): self._trace = trace - def set_track(self, track): + @track.setter + def track(self, track): if not isinstance(track, (list, tuple)): msg = 'Unable to set the track to {0} which is not a Python ' \ @@ -667,7 +926,8 @@ class SettingsFile(object): self._track = track - def set_ufs_dimension(self, dimension): + @ufs_dimension.setter + def ufs_dimension(self, dimension): if not is_integer(dimension) and not is_float(dimension): msg = 'Unable to set UFS dimension to non-integer or ' \ @@ -682,7 +942,8 @@ class SettingsFile(object): self._ufs_dimension = dimension - def set_ufs_lower_left(self, lower_left): + @ufs_lower_left.setter + def ufs_lower_left(self, lower_left): if not isinstance(lower_left, (tuple, list, np.ndarray)): msg = 'Unable to set UFS mesh lower left corner to {0} which is ' \ @@ -697,7 +958,8 @@ class SettingsFile(object): self._ufs_lower_left = lower_left - def set_ufs_upper_right(self, upper_right): + @ufs_upper_right.setter + def ufs_upper_right(self, upper_right): if not isinstance(upper_right, tuple) and \ not isinstance(upper_right, list): @@ -713,7 +975,8 @@ class SettingsFile(object): self._ufs_upper_right = upper_right - def set_dd_mesh_dimension(self, dimension): + @dd_mesh_dimension.setter + def dd_mesh_dimension(self, dimension): # TODO: remove this when domain decomposition is merged warnings.warn('This feature is not yet implemented in a release ' \ @@ -733,7 +996,8 @@ class SettingsFile(object): self._dd_mesh_dimension = dimension - def set_dd_mesh_lower_left(self, lower_left): + @dd_mesh_lower_left.setter + def dd_mesh_lower_left(self, lower_left): # TODO: remove this when domain decomposition is merged warnings.warn('This feature is not yet implemented in a release ' \ @@ -752,7 +1016,8 @@ class SettingsFile(object): self._dd_mesh_lower_left = lower_left - def set_dd_mesh_upper_right(self, upper_right): + @dd_mesh_upper_right.setter + def dd_mesh_upper_right(self, upper_right): # TODO: remove this when domain decomposition is merged warnings.warn('This feature is not yet implemented in a release ' \ @@ -772,7 +1037,8 @@ class SettingsFile(object): self._dd_mesh_upper_right = upper_right - def set_dd_nodemap(self, nodemap): + @dd_nodemap.setter + def dd_nodemap(self, nodemap): # TODO: remove this when domain decomposition is merged warnings.warn('This feature is not yet implemented in a release ' \ @@ -781,7 +1047,7 @@ class SettingsFile(object): if not isinstance(nodemap, tuple) and \ not isinstance(nodemap, list): msg = 'Unable to set DD nodemap {0} which is ' \ - 'not a Python tuple or list'.format(dimension) + 'not a Python tuple or list'.format(nodemap) raise ValueError(msg) nodemap = np.array(nodemap).flatten() @@ -798,21 +1064,39 @@ class SettingsFile(object): 'mesh'.format(len(nodemap)) raise ValueError(msg) - self._dd_mesh_dimension = dimension + self._dd_nodemap = nodemap - def set_dd_allow_leakage(self, allow): + + @dd_allow_leakage.setter + def dd_allow_leakage(self, allow): # TODO: remove this when domain decomposition is merged warnings.warn('This feature is not yet implemented in a release ' \ 'version of openmc') - if not type(allow) == bool: + if not isinstance(allow, bool): msg = 'Unable to set DD allow_leakage {0} which is ' \ - 'not a Python bool'.format(dimension) + 'not a Python bool'.format(allow) raise ValueError(msg) self._dd_allow_leakage = allow + + @dd_count_interactions.setter + def dd_count_interactions(self, interactions): + + # TODO: remove this when domain decomposition is merged + warnings.warn('This feature is not yet implemented in a release ' \ + 'version of openmc') + + if not isinstance(interactions, bool): + msg = 'Unable to set DD count_interactions {0} which is ' \ + 'not a Python bool'.format(interactions) + raise ValueError(msg) + + self._dd_count_interactions = interactions + + def create_eigenvalue_subelement(self): self.create_particles_subelement() @@ -1182,14 +1466,20 @@ class SettingsFile(object): if not self._dd_nodemap is None: subelement = ET.SubElement(element, "nodemap") - subsubelement.text = ' '.join([str(n) for n in self._dd_nodemap]) + subelement.text = ' '.join([str(n) for n in self._dd_nodemap]) subelement = ET.SubElement(element, "allow_leakage") if self._dd_allow_leakage: subelement.text = 'true' else: subelement.text = 'false' - + + subelement = ET.SubElement(element, "count_interactions") + if self._dd_count_interactions: + subelement.text = 'true' + else: + subelement.text = 'false' + def export_to_xml(self): diff --git a/src/utils/openmc/statepoint.py b/src/utils/openmc/statepoint.py index 7c4f90ae38..73ec85fc2c 100644 --- a/src/utils/openmc/statepoint.py +++ b/src/utils/openmc/statepoint.py @@ -1,13 +1,13 @@ -#!/usr/bin/env python +import copy +import struct -import struct, copy import numpy as np import scipy.stats + import openmc from openmc.constants import * - class SourceSite(object): def __init__(self): @@ -28,6 +28,26 @@ class SourceSite(object): return string + @property + def weight(self): + return self._weight + + + @property + def xyz(self): + return self._xyz + + + @property + def uvw(self): + return self._uvw + + + @property + def E(self): + return self._E + + class StatePoint(object): @@ -188,7 +208,7 @@ class StatePoint(object): # Initialize dictionaries for the Meshes # Keys - Mesh IDs # Values - Mesh objects - self._meshes = dict() + self._meshes = {} # Read the number of Meshes self._n_meshes = self._get_int(path='tallies/meshes/n_meshes')[0] @@ -205,8 +225,8 @@ class StatePoint(object): path='tallies/meshes/keys') else: - self._mesh_keys = list() - self._mesh_ids = list() + self._mesh_keys = [] + self._mesh_ids = [] # Build dictionary of Meshes base = 'tallies/meshes/mesh ' @@ -236,13 +256,13 @@ class StatePoint(object): # Create the Mesh and assign properties to it mesh = openmc.Mesh(mesh_id) - mesh.set_dimension(dimension) - mesh.set_width(width) - mesh.set_lower_left(lower_left) - mesh.set_upper_right(upper_right) + mesh.dimension = dimension + mesh.width = width + mesh.lower_left = lower_left + mesh.upper_right = upper_right #FIXME: Set the mesh type to 'rectangular' by default - mesh.set_type('rectangular') + mesh.type = 'rectangular' # Add mesh to the global dictionary of all Meshes self._meshes[mesh_id] = mesh @@ -253,7 +273,7 @@ class StatePoint(object): # Initialize dictionaries for the Tallies # Keys - Tally IDs # Values - Tally objects - self._tallies = dict() + self._tallies = {} # Read the number of tallies self._n_tallies = self._get_int(path='/tallies/n_tallies')[0] @@ -270,8 +290,8 @@ class StatePoint(object): self._n_tallies, path='tallies/keys') else: - self._tally_keys = list() - self._tally_ids = list() + self._tally_keys = [] + self._tally_ids = [] base = 'tallies/tally ' @@ -300,8 +320,8 @@ class StatePoint(object): # Create Tally object and assign basic properties tally = openmc.Tally(tally_key, label) - tally.set_estimator(ESTIMATOR_TYPES[estimator_type]) - tally.set_num_realizations(n_realizations) + tally.estimator = ESTIMATOR_TYPES[estimator_type] + tally.num_realizations = n_realizations # Read the number of Filters n_filters = self._get_int( @@ -343,11 +363,11 @@ class StatePoint(object): # Create Filter object filter = openmc.Filter(FILTER_TYPES[filter_type], bins) - filter.set_stride(stride) - filter.set_num_bins(n_bins) + filter.stride = stride + filter.num_bins = n_bins if FILTER_TYPES[filter_type] == 'mesh': - filter.set_mesh(self._meshes[bins]) + filter.mesh = self._meshes[bins] # Add Filter to the Tally tally.add_filter(filter) @@ -370,7 +390,7 @@ class StatePoint(object): n_score_bins = self._get_int( path='{0}{1}/n_score_bins'.format(base, tally_key))[0] - tally.set_num_score_bins(n_score_bins) + tally.num_score_bins = n_score_bins scores = [SCORE_TYPES[j] for j in self._get_int( n_score_bins, path='{0}{1}/score_bins'.format(base, tally_key))] @@ -378,12 +398,12 @@ class StatePoint(object): path='{0}{1}/n_user_score_bins'.format(base, tally_key))[0] # Read scattering moment order strings (e.g., P3, Y-1,2, etc.) - moments = list() + moments = [] subbase = '{0}{1}/moments/'.format(base, tally_key) # Extract the moment order string for each score for k in range(len(scores)): - moment = self._get_string(8, + moment = self._get_string(8, path='{0}order{1}'.format(subbase, k+1)) moment = moment.lstrip('[\'') moment = moment.rstrip('\']') @@ -440,7 +460,7 @@ class StatePoint(object): tally = self._tallies[tally_key] # Compute the total number of bins for this Tally - num_tot_bins = tally.get_num_bins() + num_tot_bins = tally.num_bins # Extract Tally data from the file if self._hdf5: @@ -457,9 +477,9 @@ class StatePoint(object): nonzero = lambda val: 1 if not val else val # Reshape the results arrays - new_shape = (nonzero(tally.get_num_filter_bins()), - nonzero(tally.get_num_nuclides()), - nonzero(tally.get_num_score_bins())) + new_shape = (nonzero(tally.num_filter_bins), + nonzero(tally.num_nuclides), + nonzero(tally.num_score_bins)) sum = np.reshape(sum, new_shape) sum_sq = np.reshape(sum_sq, new_shape) @@ -666,28 +686,28 @@ class StatePoint(object): for filter in tally._filters: if filter._type == 'surface': - surface_ids = list() + surface_ids = [] for bin in filter._bins: surface_ids.append(summary.surfaces[bin]._id) - filter.set_bin_edges(surface_ids) + filter.bins = surface_ids if filter._type in ['cell', 'distribcell']: - distribcell_ids = list() + distribcell_ids = [] for bin in filter._bins: distribcell_ids.append(summary.cells[bin]._id) - filter.set_bin_edges(distribcell_ids) + filter.bins = distribcell_ids if filter._type == 'universe': - universe_ids = list() + universe_ids = [] for bin in filter._bins: universe_ids.append(summary.universes[bin]._id) - filter.set_bin_edges(universe_ids) + filter.bins = universe_ids if filter._type == 'material': - material_ids = list() + material_ids = [] for bin in filter._bins: material_ids.append(summary.materials[bin]._id) - filter.set_bin_edges(material_ids) + filter.bins = material_ids self._with_summary = True diff --git a/src/utils/openmc/summary.py b/src/utils/openmc/summary.py index cde2067213..2e558fe726 100644 --- a/src/utils/openmc/summary.py +++ b/src/utils/openmc/summary.py @@ -1,8 +1,7 @@ -#!/usr/bin/env python +import numpy as np import openmc from openmc.opencg_compatible import get_opencg_geometry -import numpy as np try: import h5py @@ -65,7 +64,7 @@ class Summary(object): # Initialize dictionary for each Nuclide # Keys - Nuclide ZAIDs # Values - Nuclide objects - self.nuclides = dict() + self.nuclides = {} for key in self._f['nuclides'].keys(): @@ -87,7 +86,7 @@ class Summary(object): self.nuclides[zaid] = openmc.Element(name=name, xs=xs) else: self.nuclides[zaid] = openmc.Nuclide(name=name, xs=xs) - self.nuclides[zaid].set_zaid(zaid) + self.nuclides[zaid].zaid = zaid def _read_materials(self): @@ -97,7 +96,7 @@ class Summary(object): # Initialize dictionary for each Material # Keys - Material keys # Values - Material objects - self.materials = dict() + self.materials = {} for key in self._f['materials'].keys(): @@ -112,8 +111,8 @@ class Summary(object): nuclides = self._f['materials'][key]['nuclides'][...] n_sab = self._f['materials'][key]['n_sab'][0] - sab_names = list() - sab_xs = list() + sab_names = [] + sab_xs = [] # Read the names of the S(a,b) tables for this Material for i in range(1, n_sab+1): @@ -157,7 +156,7 @@ class Summary(object): # Initialize dictionary for each Surface # Keys - Surface keys # Values - Surfacee objects - self.surfaces = dict() + self.surfaces = {} for key in self._f['geometry/surfaces'].keys(): @@ -241,7 +240,7 @@ class Summary(object): # Initialize dictionary for each Cell # Keys - Cell keys # Values - Cell objects - self.cells = dict() + self.cells = {} # Initialize dictionary for each Cell's fill # (e.g., Material, Universe or Lattice ID) @@ -249,7 +248,7 @@ class Summary(object): # the corresponding objects # Keys - Cell keys # Values - Filling Material, Universe or Lattice ID - self._cell_fills = dict() + self._cell_fills = {} for key in self._f['geometry/cells'].keys(): @@ -271,7 +270,7 @@ class Summary(object): if 'surfaces' in self._f['geometry/cells'][key].keys(): surfaces = self._f['geometry/cells'][key]['surfaces'][...] else: - surfaces = list() + surfaces = [] # Create this Cell cell = openmc.Cell(cell_id=cell_id, name=name) @@ -282,14 +281,14 @@ class Summary(object): translation = \ self._f['geometry/cells'][key]['translation'][...] translation = np.asarray(translation, dtype=np.float64) - cell.set_translation(translation) + cell.translation = translation rotated = self._f['geometry/cells'][key]['rotated'][0] if rotated: rotation = \ self._f['geometry/cells'][key]['rotation'][...] rotation = np.asarray(rotation, dtype=np.int) - cell.set_rotation(rotation) + cell.rotation = rotation # Store Cell fill information for after Universe/Lattice creation self._cell_fills[index] = (fill_type, fill) @@ -313,7 +312,7 @@ class Summary(object): # Initialize dictionary for each Universe # Keys - Universe keys # Values - Universe objects - self.universes = dict() + self.universes = {} for key in self._f['geometry/universes'].keys(): @@ -343,7 +342,7 @@ class Summary(object): # Initialize lattices for each Lattice # Keys - Lattice keys # Values - Lattice objects - self.lattices = dict() + self.lattices = {} for key in self._f['geometry/lattices'].keys(): @@ -356,9 +355,9 @@ class Summary(object): lattice_type = self._f['geometry/lattices'][key]['type'][...][0] if lattice_type == 'rectangular': - dimension = self._f['geometry/lattices'][key]['n_cells'][...] + dimension = self._f['geometry/lattices'][key]['dimension'][...] lower_left = \ - self._f['geometry/lattices'][key]['lower_left'][...] + self._f['geometry/lattices'][key]['lower_left'][...] pitch = self._f['geometry/lattices'][key]['pitch'][...] outer = self._f['geometry/lattices'][key]['outer'][0] @@ -369,13 +368,14 @@ class Summary(object): # Create the Lattice lattice = openmc.RectLattice(lattice_id=lattice_id, name=name) - lattice.set_dimension(tuple(dimension)) - lattice.set_lower_left(lower_left) - lattice.set_pitch(pitch) + lattice.dimension = tuple(dimension) + lattice.lower_left = lower_left + lattice.pitch = pitch + # If the Universe specified outer the Lattice is not void (-22) if outer != -22: - lattice.set_outer(self.universes[outer]) + lattice.outer = self.universes[outer] # Build array of Universe pointers for the Lattice universes = \ @@ -392,7 +392,7 @@ class Summary(object): universes = np.transpose(universes, (1,0,2)) universes.shape = shape universes = universes[:,::-1,:] - lattice.set_universes(universes) + lattice.universes = universes # Add the Lattice to the global dictionary of all Lattices self.lattices[index] = lattice @@ -408,15 +408,15 @@ class Summary(object): self._f['geometry/lattices'][key]['universes'][...] # Create the Lattice - lattice = openmc.HexLattice(lattice_id=lattice_id, name=name) - lattice.set_num_rings(n_rings) - lattice.set_num_axial(n_axial) - lattice.set_center(center) - lattice.set_pitch(pitch) + lattice = openmc.HexLattice(lattice_id=lattice_id, name) + lattice.num_rings = n_rings + lattice.num_axial = n_axial + lattice.center = center + lattice.pitch = pitch # If the Universe specified outer the Lattice is not void (-22) if outer != -22: - lattice.set_outer(self.universes[outer]) + lattice.outer = self.universes[outer] # Build array of Universe pointers for the Lattice universes = \ @@ -457,11 +457,11 @@ class Summary(object): fill = self.get_lattice_by_id(fill_id) # Set the fill for the Cell - self.cells[cell_key].set_fill(fill) + self.cells[cell_key].fill = fill # Set the root universe for the Geometry root_universe = self.get_universe_by_id(0) - self.openmc_geometry.set_root_universe(root_universe) + self.openmc_geometry.root_universe = root_universe def make_opencg_geometry(self): diff --git a/src/utils/openmc/surface.py b/src/utils/openmc/surface.py index 5894cf2503..f6067a5fdf 100644 --- a/src/utils/openmc/surface.py +++ b/src/utils/openmc/surface.py @@ -1,8 +1,7 @@ -#!/usr/bin/env python +from xml.etree import ElementTree as ET from openmc.checkvalue import * from openmc.constants import BC_TYPES -from xml.etree import ElementTree as ET # A static variable for auto-generated Surface IDs @@ -16,29 +15,51 @@ def reset_auto_surface_id(): class Surface(object): - def __init__(self, surface_id=None, bc_type='transmission', name=''): + def __init__(self, surface_id=None, boundary_type='transmission', name=''): # Initialize class attributes - self._id = None - self._name = '' + self.id = surface_id + self.name = name self._type = '' - self._bc_type = '' + self.boundary_type = boundary_type # A dictionary of the quadratic surface coefficients # Key - coefficeint name # Value - coefficient value - self._coeffs = dict() + self._coeffs = {} # An ordered list of the coefficient names to export to XML in the # proper order - self._coeff_keys = list() - - self.set_id(surface_id) - self.set_boundary_type(bc_type) - self.set_name(name) + self._coeff_keys = [] - def set_id(self, surface_id=None): + @property + def id(self): + return self._id + + + @property + def name(self): + return self._name + + + @property + def type(self): + return self._type + + + @property + def boundary_type(self): + return self._boundary_type + + + @property + def coeffs(self): + return self._coeffs + + + @id.setter + def id(self, surface_id): if surface_id is None: global AUTO_SURFACE_ID @@ -60,7 +81,8 @@ class Surface(object): self._id = surface_id - def set_name(self, name): + @name.setter + def name(self, name): if not is_string(name): msg = 'Unable to set name for Surface ID={0} with a non-string ' \ @@ -71,21 +93,22 @@ class Surface(object): self._name = name - def set_boundary_type(self, bc_type): + @boundary_type.setter + def boundary_type(self, boundary_type): - if not is_string(bc_type): + if not is_string(boundary_type): msg = 'Unable to set boundary type for Surface ID={0} with a ' \ - 'non-string value {1}'.format(self._id, bc_type) + 'non-string value {1}'.format(self._id, boundary_type) raise ValueError(msg) - elif not bc_type in BC_TYPES.values(): + elif not boundary_type in BC_TYPES.values(): msg = 'Unable to set boundary type for Surface ID={0} to ' \ '{1} which is not trasmission, vacuum or ' \ - 'reflective'.format(bc_type) + 'reflective'.format(boundary_type) raise ValueError(msg) else: - self._bc_type = bc_type + self._boundary_type = boundary_type def __repr__(self): @@ -94,7 +117,7 @@ class Surface(object): string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type) - string += '{0: <16}{1}{2}\n'.format('\tBoundary', '=\t', self._bc_type) + string += '{0: <16}{1}{2}\n'.format('\tBoundary', '=\t', self._boundary_type) coeffs = '{0: <16}'.format('\tCoefficients') + '\n' @@ -115,7 +138,7 @@ class Surface(object): element.set("label", str(self._name)) element.set("type", self._type) - element.set("boundary", self._bc_type) + element.set("boundary", self._boundary_type) coeffs = '' @@ -130,33 +153,50 @@ class Surface(object): class Plane(Surface): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', A=None, B=None, C=None, D=None, name='',): # Initialize Plane class attributes - super(Plane, self).__init__(surface_id, bc_type, name=name) + super(Plane, self).__init__(surface_id, boundary_type, name=name) - self._A = None - self._B = None - self._C = None - self._D = None self._type = 'plane' self._coeff_keys = ['A', 'B', 'C', 'D'] if not A is None: - self.set_A(A) + self.a = A if not B is None: - self.set_B(B) + self.b = B if not C is None: - self.set_C(C) + self.c = C if not D is None: - self.set_D(D) + self.d = D - def set_A(self, A): + @property + def a(self): + return self.coeffs['A'] + + + @property + def b(self): + return self.coeffs['B'] + + + @property + def c(self): + return self.coeffs['C'] + + + @property + def d(self): + return self.coeffs['D'] + + + @a.setter + def a(self, A): if not is_integer(A) and not is_float(A): msg = 'Unable to set A coefficient for Plane ID={0} to a ' \ @@ -166,7 +206,8 @@ class Plane(Surface): self._coeffs['A'] = A - def set_B(self, B): + @b.setter + def b(self, B): if not is_integer(B) and not is_float(B): msg = 'Unable to set B coefficient for Plane ID={0} to a ' \ @@ -176,7 +217,8 @@ class Plane(Surface): self._coeffs['B'] = B - def set_C(self, C): + @c.setter + def c(self, C): if not is_integer(C) and not is_float(C): msg = 'Unable to set C coefficient for Plane ID={0} to a ' \ @@ -186,7 +228,8 @@ class Plane(Surface): self._coeffs['C'] = C - def set_D(self, D): + @d.setter + def d(self, D): if not is_integer(D) and not is_float(D): msg = 'Unable to set D coefficient for Plane ID={0} to a ' \ @@ -199,21 +242,26 @@ class Plane(Surface): class XPlane(Plane): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', x0=None, name=''): # Initialize XPlane class attributes - super(XPlane, self).__init__(surface_id, bc_type, name=name) + super(XPlane, self).__init__(surface_id, boundary_type, name=name) - self._x0 = None self._type = 'x-plane' self._coeff_keys = ['x0'] if not x0 is None: - self.set_X0(x0) + self.x0 = x0 - def set_X0(self, x0): + @property + def x0(self): + return self.coeff['x0'] + + + @x0.setter + def x0(self, x0): if not is_integer(x0) and not is_float(x0): msg = 'Unable to set x0 coefficient for XPlane ID={0} to a ' \ @@ -226,21 +274,26 @@ class XPlane(Plane): class YPlane(Plane): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', y0=None, name=''): # Initialize YPlane class attributes - super(YPlane, self).__init__(surface_id, bc_type, name=name) + super(YPlane, self).__init__(surface_id, boundary_type, name=name) - self._y0 = None self._type = 'y-plane' self._coeff_keys = ['y0'] if not y0 is None: - self.set_Y0(y0) + self.y0 = y0 - def set_Y0(self, y0): + @property + def y0(self): + return self.coeffs['y0'] + + + @y0.setter + def y0(self, y0): if not is_integer(y0) and not is_float(y0): msg = 'Unable to set y0 coefficient for XPlane ID={0} to a ' \ @@ -253,21 +306,26 @@ class YPlane(Plane): class ZPlane(Plane): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', z0=None, name=''): # Initialize ZPlane class attributes - super(ZPlane, self).__init__(surface_id, bc_type, name=name) + super(ZPlane, self).__init__(surface_id, boundary_type, name=name) - self._z0 = None self._type = 'z-plane' self._coeff_keys = ['z0'] if not z0 is None: - self.set_Z0(z0) + self.z0 = z0 - def set_Z0(self, z0): + @property + def z0(self): + return self.coeffs['z0'] + + + @z0.setter + def z0(self, z0): if not is_integer(z0) and not is_float(z0): msg = 'Unable to set z0 coefficient for ZPlane ID={0} to a ' \ @@ -280,20 +338,25 @@ class ZPlane(Plane): class Cylinder(Surface): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', R=None, name=''): # Initialize Cylinder class attributes - super(Cylinder, self).__init__(surface_id, bc_type, name=name) + super(Cylinder, self).__init__(surface_id, boundary_type, name=name) - self._R = None self._coeff_keys = ['R'] if not R is None: - self.set_R(R) + self.r = R - def set_R(self, R): + @property + def r(self): + return self.coeffs['R'] + + + @r.setter + def r(self, R): if not is_integer(R) and not is_float(R): msg = 'Unable to set R coefficient for Cylinder ID={0} to a ' \ @@ -306,25 +369,34 @@ class Cylinder(Surface): class XCylinder(Cylinder): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', y0=None, z0=None, R=None, name=''): # Initialize XCylinder class attributes - super(XCylinder, self).__init__(surface_id, bc_type, R, name=name) + super(XCylinder, self).__init__(surface_id, boundary_type, R, name=name) self._type = 'x-cylinder' - self._y0 = None - self._z0 = None self._coeff_keys = ['y0', 'z0', 'R'] if not y0 is None: - self.set_Y0(y0) + self.y0 = y0 if not z0 is None: - self.set_Z0(z0) + self.z0 = z0 - def set_Y0(self, y0): + @property + def y0(self): + return self.coeffs['y0'] + + + @property + def z0(self): + return self.coeffs['z0'] + + + @y0.setter + def y0(self, y0): if not is_integer(y0) and not is_float(y0): msg = 'Unable to set y0 coefficient for XCylinder ID={0} to a ' \ @@ -334,7 +406,8 @@ class XCylinder(Cylinder): self._coeffs['y0'] = y0 - def set_Z0(self, z0): + @z0.setter + def z0(self, z0): if not is_integer(z0) and not is_float(z0): msg = 'Unable to set z0 coefficient for XCylinder ID={0} to a ' \ @@ -347,25 +420,34 @@ class XCylinder(Cylinder): class YCylinder(Cylinder): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', x0=None, z0=None, R=None, name=''): # Initialize YCylinder class attributes - super(YCylinder, self).__init__(surface_id, bc_type, R, name=name) + super(YCylinder, self).__init__(surface_id, boundary_type, R, name=name) self._type = 'y-cylinder' - self._x0 = None - self._z0 = None self._coeff_keys = ['x0', 'z0', 'R'] if not x0 is None: - self.set_X0(x0) + self.x0 = x0 if not z0 is None: - self.set_Z0(z0) + self.z0 = z0 - def set_X0(self, x0): + @property + def x0(self): + return self.coeffs['x0'] + + + @property + def z0(self): + return self.coeffs['z0'] + + + @x0.setter + def x0(self, x0): if not is_integer(x0) and not is_float(x0): msg = 'Unable to set x0 coefficient for YCylinder ID={0} to a ' \ @@ -375,7 +457,8 @@ class YCylinder(Cylinder): self._coeffs['x0'] = x0 - def set_Z0(self, z0): + @z0.setter + def z0(self, z0): if not is_integer(z0) and not is_float(z0): msg = 'Unable to set z0 coefficient for YCylinder ID={0} to a ' \ @@ -387,26 +470,35 @@ class YCylinder(Cylinder): class ZCylinder(Cylinder): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', x0=None, y0=None, R=None, name=''): # Initialize ZCylinder class attributes # Initialize YPlane class attributes - super(ZCylinder, self).__init__(surface_id, bc_type, R, name=name) + super(ZCylinder, self).__init__(surface_id, boundary_type, R, name=name) self._type = 'z-cylinder' - self._x0 = None - self._y0 = None self._coeff_keys = ['x0', 'y0', 'R'] if not x0 is None: - self.set_X0(x0) + self.x0 = x0 if not y0 is None: - self.set_Y0(y0) + self.y0 = y0 - def set_X0(self, x0): + @property + def x0(self): + return self.coeffs['x0'] + + + @property + def y0(self): + return self.coeffs['y0'] + + + @x0.setter + def x0(self, x0): if not is_integer(x0) and not is_float(x0): msg = 'Unable to set x0 coefficient for ZCylinder ID={0} to a ' \ @@ -416,7 +508,8 @@ class ZCylinder(Cylinder): self._coeffs['x0'] = x0 - def set_Y0(self, y0): + @y0.setter + def y0(self, y0): if not is_integer(y0) and not is_float(y0): msg = 'Unable to set y0 coefficient for ZCylinder ID={0} to a ' \ @@ -429,33 +522,50 @@ class ZCylinder(Cylinder): class Sphere(Surface): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', x0=None, y0=None, z0=None, R=None, name=''): # Initialize Sphere class attributes - super(Sphere, self).__init__(surface_id, bc_type, name=name) + super(Sphere, self).__init__(surface_id, boundary_type, name=name) self._type = 'sphere' - self._x0 = None - self._y0 = None - self._z0 = None - self._R = None self._coeff_keys = ['x0', 'y0', 'z0', 'R'] if not x0 is None: - self.set_X0(x0) + self.x0 = x0 if not y0 is None: - self.set_Y0(y0) + self.y0 = y0 if not z0 is None: - self.set_Z0(z0) + self.z0 = z0 if not R is None: - self.set_Z0(R) + self.r = R - def set_X0(self, x0): + @property + def x0(self): + return self.coeffs['x0'] + + + @property + def y0(self): + return self.coeffs['y0'] + + + @property + def z0(self): + return self.coeffs['z0'] + + + @property + def r(self): + return self.coeffs['R'] + + + @x0.setter + def x0(self, x0): if not is_integer(x0) and not is_float(x0): msg = 'Unable to set x0 coefficient for Sphere ID={0} to a ' \ @@ -465,7 +575,8 @@ class Sphere(Surface): self._coeffs['x0'] = x0 - def set_Y0(self, y0): + @y0.setter + def y0(self, y0): if not is_integer(y0) and not is_float(y0): msg = 'Unable to set y0 coefficient for Sphere ID={0} to a ' \ @@ -475,7 +586,8 @@ class Sphere(Surface): self._coeffs['y0'] = y0 - def set_Z0(self, z0): + @z0.setter + def z0(self, z0): if not is_integer(z0) and not is_float(z0): msg = 'Unable to set z0 coefficient for Sphere ID={0} to a ' \ @@ -485,7 +597,8 @@ class Sphere(Surface): self._coeffs['z0'] = z0 - def set_R(self, R): + @r.setter + def r(self, R): if not is_integer(R) and not is_float(R): msg = 'Unable to set R coefficient for Sphere ID={0} to a ' \ @@ -498,32 +611,49 @@ class Sphere(Surface): class Cone(Surface): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', x0=None, y0=None, z0=None, R2=None, name=''): # Initialize Cone class attributes - super(Cone, self).__init__(surface_id, bc_type, name=name) + super(Cone, self).__init__(surface_id, boundary_type, name=name) - self._x0 = None - self._y0 = None - self._z0 = None - self._R2 = None self._coeff_keys = ['x0', 'y0', 'z0', 'R2'] if not x0 is None: - self.set_X0(x0) + self.x0 = x0 if not y0 is None: - self.set_Y0(y0) + self.y0 = y0 if not z0 is None: - self.set_Z0(z0) + self.z0 = z0 if not R2 is None: - self.set_Z0(R2) + self.r2 = R2 - def set_X0(self, x0): + @property + def x0(self): + return self.coeffs['x0'] + + + @property + def y0(self): + return self.coeffs['y0'] + + + @property + def z0(self): + return self.coeffs['z0'] + + + @property + def r2(self): + return self.coeffs['r2'] + + + @x0.setter + def x0(self, x0): if not is_integer(x0) and not is_float(x0): msg = 'Unable to set x0 coefficient for Cone ID={0} to a ' \ @@ -533,7 +663,8 @@ class Cone(Surface): self._coeffs['x0'] = x0 - def set_Y0(self, y0): + @y0.setter + def y0(self, y0): if not is_integer(y0) and not is_float(y0): msg = 'Unable to set y0 coefficient for Cone ID={0} to a ' \ @@ -543,7 +674,8 @@ class Cone(Surface): self._coeffs['y0'] = y0 - def set_Z0(self, z0): + @z0.setter + def z0(self, z0): if not is_integer(z0) and not is_float(z0): msg = 'Unable to set z0 coefficient for Cone ID={0} to a ' \ @@ -553,7 +685,8 @@ class Cone(Surface): self._coeffs['z0'] = z0 - def set_R2(self, R2): + @r2.setter + def r2(self, R2): if not is_integer(R2) and not is_float(R2): msg = 'Unable to set R^2 coefficient for Cone ID={0} to a ' \ @@ -566,11 +699,12 @@ class Cone(Surface): class XCone(Cone): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', x0=None, y0=None, z0=None, R2=None, name=''): # Initialize XCone class attributes - super(XCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name) + super(XCone, self).__init__(surface_id, boundary_type, x0, y0, + z0, R2, name=name) self._type = 'x-cone' @@ -578,11 +712,12 @@ class XCone(Cone): class YCone(Cone): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', x0=None, y0=None, z0=None, R2=None, name=''): # Initialize YCone class attributes - super(YCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name) + super(YCone, self).__init__(surface_id, boundary_type, x0, y0, z0, + R2, name=name) self._type = 'y-cone' @@ -590,10 +725,11 @@ class YCone(Cone): class ZCone(Cone): - def __init__(self, surface_id=None, bc_type='transmission', + def __init__(self, surface_id=None, boundary_type='transmission', x0=None, y0=None, z0=None, R2=None, name=''): # Initialize ZCone class attributes - super(ZCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name) + super(ZCone, self).__init__(surface_id, boundary_type, x0, y0, z0, + R2, name=name) self._type = 'z-cone' diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 4dfb77f4ee..810611b26c 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -1,12 +1,11 @@ -#!/usr/bin/env python +import copy +import os +from xml.etree import ElementTree as ET from openmc import Nuclide from openmc.clean_xml import * from openmc.checkvalue import * from openmc.constants import * -from xml.etree import ElementTree as ET -import numpy as np -import os, copy # "Static" variables for auto-generated Tally and Mesh IDs @@ -29,21 +28,12 @@ class Filter(object): # Initialize Filter class attributes def __init__(self, type=None, bins=None): - self._type = None - self._bins = None + self.type = type + self.bins = bins self._mesh = None self._offset = -1 self._stride = None - # FIXME - self._num_bins = None - - if not type is None: - self.set_type(type) - - if not bins is None: - self.set_bin_edges(bins) - def __eq__(self, filter2): @@ -64,7 +54,7 @@ class Filter(object): def __hash__(self): - hashable = list() + hashable = [] hashable.append(self._type) hashable.append(self._bins) return hash(tuple(hashable)) @@ -94,9 +84,43 @@ class Filter(object): return existing - def set_type(self, type): + @property + def type(self): + return self._type - if not type in FILTER_TYPES.values(): + + @property + def bins(self): + return self._bins + + + @property + def num_bins(self): + return self._num_bins + + + @property + def mesh(self): + return self._mesh + + + @property + def offset(self): + return self._offset + + + @property + def stride(self): + return self._stride + + + @type.setter + def type(self, type): + + if type is None: + self._type = type + + elif not type in FILTER_TYPES.values(): msg = 'Unable to set Filter type to {0} since it is not one ' \ 'of the supported types'.format(type) raise ValueError(msg) @@ -104,10 +128,14 @@ class Filter(object): self._type = type - def set_bin_edges(self, bins): + @bins.setter + def bins(self, bins): - if self._type is None: - msg = 'Unable to set bin edges for Filter to {0} since ' \ + if bins is None: + self.num_bins = 0 + + elif self._type is None: + msg = 'Unable to set bins for Filter to {0} since ' \ 'the Filter type has not yet been set'.format(bins) raise ValueError(msg) @@ -125,12 +153,12 @@ class Filter(object): for edge in bins: if not is_integer(edge): - msg = 'Unable to add bin edge {0} to a {1} Filter since ' \ + msg = 'Unable to add bin {0} to a {1} Filter since ' \ 'it is a non-integer'.format(edge, self._type) raise ValueError(msg) elif edge < 0: - msg = 'Unable to add bin edge {0} to a {1} Filter since ' \ + msg = 'Unable to add bin {0} to a {1} Filter since ' \ 'it is a negative integer'.format(edge, self._type) raise ValueError(msg) @@ -164,17 +192,17 @@ class Filter(object): elif self._type == 'mesh': if not len(bins) == 1: - msg = 'Unable to add bin edges {0} to a mesh Filter since ' \ + msg = 'Unable to add bins {0} to a mesh Filter since ' \ 'only a single mesh can be used per tally'.format(bins) raise ValueError(msg) elif not is_integer(bins[0]): - msg = 'Unable to add bin edge {0} to mesh Filter since it ' \ + msg = 'Unable to add bin {0} to mesh Filter since it ' \ 'is a non-integer'.format(bins[0]) raise ValueError(msg) elif bins[0] < 0: - msg = 'Unable to add bin edge {0} to mesh Filter since it ' \ + msg = 'Unable to add bin {0} to mesh Filter since it ' \ 'is a negative integer'.format(bins[0]) raise ValueError(msg) @@ -183,7 +211,8 @@ class Filter(object): # FIXME - def set_num_bins(self, num_bins): + @num_bins.setter + def num_bins(self, num_bins): if not is_integer(num_bins) or num_bins < 0: msg = 'Unable to set the number of bins {0} for a {1} Filter ' \ @@ -194,7 +223,8 @@ class Filter(object): self._num_bins = num_bins - def set_mesh(self, mesh): + @mesh.setter + def mesh(self, mesh): if not isinstance(mesh, Mesh): msg = 'Unable to set Mesh to {0} for Filter since it is not a ' \ @@ -202,11 +232,12 @@ class Filter(object): raise ValueError(msg) self._mesh = mesh - self.set_type('mesh') - self.set_bin_edges(self._mesh._id) + self.type = 'mesh' + self.bins = self._mesh._id - def set_offset(self, offset): + @offset.setter + def offset(self, offset): if not is_integer(offset): msg = 'Unable to set offset {0} for a {1} Filter since it is a ' \ @@ -216,7 +247,8 @@ class Filter(object): self._offset = offset - def set_stride(self, stride): + @stride.setter + def stride(self, stride): if not is_integer(stride): msg = 'Unable to set stride {0} for a {1} Filter since it is a ' \ @@ -231,21 +263,6 @@ class Filter(object): self._stride = stride - def get_num_bins(self): - - # FIXME - #if self._type == 'mesh': - # num_bins = self._mesh.getNumMeshCells() - #elif self._type == 'energy' or self._type == 'energyout': - # num_bins = len(self._bins) - 1 - #else: - # num_bins = len(self._bins) - - #return num_bins - - return self._num_bins - - def get_bin_index(self, bin): try: @@ -274,17 +291,14 @@ class Mesh(object): def __init__(self, mesh_id=None, name=''): # Initialize Mesh class attributes - self._id = None - self._name = '' + self.id = mesh_id + self.name = name self._type = 'rectangular' self._dimension = None self._lower_left = None self._upper_right = None self._width = None - self.set_id(mesh_id) - self.set_name(name) - def __deepcopy__(self, memo): @@ -311,7 +325,48 @@ class Mesh(object): return existing - def set_id(self, mesh_id=None): + @property + def id(self): + return self._id + + + @property + def name(self): + return self._name + + + @property + def type(self): + return self._type + + + @property + def dimension(self): + return self._dimension + + + @property + def lower_left(self): + return self._lower_left + + + @property + def upper_right(self): + return self._upper_right + + + @property + def width(self): + return self._width + + + @property + def num_mesh_cells(self): + return np.prod(self._dimension) + + + @id.setter + def id(self, mesh_id): if mesh_id is None: global AUTO_MESH_ID @@ -332,7 +387,8 @@ class Mesh(object): self._id = mesh_id - def set_name(self, name): + @name.setter + def name(self, name): if not is_string(name): msg = 'Unable to set name for Mesh ID={0} with a non-string ' \ @@ -343,7 +399,8 @@ class Mesh(object): self._name = name - def set_type(self, type): + @type.setter + def type(self, type): if not is_string(type): msg = 'Unable to set Mesh ID={0} for type {1} which is not ' \ @@ -359,7 +416,8 @@ class Mesh(object): self._type = type - def set_dimension(self, dimension): + @dimension.setter + def dimension(self, dimension): if not isinstance(dimension, (tuple, list, np.ndarray)): msg = 'Unable to set Mesh ID={0} with dimension {1} which is ' \ @@ -382,7 +440,8 @@ class Mesh(object): self._dimension = dimension - def set_lower_left(self, lower_left): + @lower_left.setter + def lower_left(self, lower_left): if not isinstance(lower_left, (tuple, list, np.ndarray)): msg = 'Unable to set Mesh ID={0} with lower_left {1} which is ' \ @@ -406,7 +465,8 @@ class Mesh(object): self._lower_left = lower_left - def set_upper_right(self, upper_right): + @upper_right.setter + def upper_right(self, upper_right): if not isinstance(upper_right, (tuple, list, np.ndarray)): msg = 'Unable to set Mesh ID={0} with upper_right {1} which ' \ @@ -430,7 +490,8 @@ class Mesh(object): self._upper_right = upper_right - def set_width(self, width): + @width.setter + def width(self, width): if not width is None: @@ -469,10 +530,6 @@ class Mesh(object): return string - def get_num_mesh_cells(self): - return np.prod(self._dimension) - - def get_mesh_xml(self): element = ET.Element("mesh") @@ -530,11 +587,11 @@ class Tally(object): def __init__(self, tally_id=None, label=''): # Initialize Tally class attributes - self._id = None - self._label = None - self._filters = list() - self._nuclides = list() - self._scores = list() + self.id = tally_id + self.label = label + self._filters = [] + self._nuclides = [] + self._scores = [] self._estimator = None self._num_score_bins = 0 @@ -545,9 +602,6 @@ class Tally(object): self._mean = None self._std_dev = None - self.set_id(tally_id) - self.set_label(label) - def __deepcopy__(self, memo): @@ -567,15 +621,15 @@ class Tally(object): clone._mean = copy.deepcopy(self._mean, memo) clone._std_dev = copy.deepcopy(self._std_dev, memo) - clone._filters = list() + clone._filters = [] for filter in self._filters: clone.add_filter(copy.deepcopy(filter, memo)) - clone._nuclides = list() + clone._nuclides = [] for nuclide in self._nuclides: clone.add_nuclide(copy.deepcopy(nuclide, memo)) - clone._scores = list() + clone._scores = [] for score in self._scores: clone.add_score(score) @@ -612,7 +666,7 @@ class Tally(object): def __hash__(self): - hashable = list() + hashable = [] for filter in self._filters: hashable.append((filter._type, tuple(filter._bins))) @@ -639,7 +693,102 @@ class Tally(object): new_tally._std_dev = np.sqrt(self._std_dev**2 + other._std_dev**2) - def set_estimator(self, estimator): + @property + def id(self): + return self._id + + + @property + def label(self): + return self._label + + + @property + def filters(self): + return self._filters + + + @property + def num_filters(self): + return len(self._filters) + + + @property + def nuclides(self): + return self._nuclides + + + @property + def num_nuclides(self): + return len(self._nuclides) + + + @property + def scores(self): + return self._scores + + + @property + def num_scores(self): + return len(self._scores) + + + @property + def num_score_bins(self): + return self._num_score_bins + + + @property + def num_filter_bins(self): + + num_bins = 1 + + for filter in self._filters: + num_bins *= filter.num_bins + + return num_bins + + + @property + def num_bins(self): + num_bins = self.num_filter_bins + num_bins *= self.num_nuclides + num_bins *= self.num_score_bins + return num_bins + + + @property + def estimator(self): + return self._estimator + + + @property + def num_realizations(self): + return self._num_realizations + + + @property + def sum(self): + return self._sum + + + @property + def sum_sq(self): + return self._sum_sq + + + @property + def mean(self): + return self._mean + + + @property + def std_dev(self): + return self._std_dev + + + @estimator.setter + def estimator(self, estimator): if not estimator in ['analog', 'tracklength']: msg = 'Unable to set the estimator for Tally ID={0} to {1} since ' \ 'it is not a valid estimator type'.format(self._id, estimator) @@ -648,7 +797,8 @@ class Tally(object): self._estimator = estimator - def set_id(self, tally_id=None): + @id.setter + def id(self, tally_id): if tally_id is None: global AUTO_TALLY_ID @@ -669,7 +819,8 @@ class Tally(object): self._id = tally_id - def set_label(self, label=None): + @label.setter + def label(self, label): if not is_string(label): msg = 'Unable to set name for Tally ID={0} with a non-string ' \ @@ -710,11 +861,13 @@ class Tally(object): self._scores.append(score) - def set_num_score_bins(self, num_score_bins): + @num_score_bins.setter + def num_score_bins(self, num_score_bins): self._num_score_bins = num_score_bins - def set_num_realizations(self, num_realizations): + @num_realizations.setter + def num_realizations(self, num_realizations): if not is_integer(num_realizations): msg = 'Unable to set the number of realizations to {0} for ' \ @@ -803,7 +956,10 @@ class Tally(object): string += '{0: <16}{1}'.format('\tNuclides', '=\t') for nuclide in self._nuclides: - string += '{0} '.format(nuclide._name) + if isinstance(nuclide, Nuclide): + string += '{0} '.format(nuclide._name) + else: + string += '{0} '.format(nuclide) string += '\n' @@ -813,39 +969,6 @@ class Tally(object): return string - def get_num_filters(self): - return len(self._filters) - - - def get_num_filter_bins(self): - - num_bins = 1 - - for filter in self._filters: - num_bins *= filter.get_num_bins() - - return num_bins - - - def get_num_nuclides(self): - return len(self._nuclides) - - - def get_num_scores(self): - return len(self._scores) - - - def get_num_score_bins(self): - return self._num_score_bins - - - def get_num_bins(self): - num_bins = self.get_num_filter_bins() - num_bins *= self.get_num_nuclides() - num_bins *= self.get_num_score_bins() - return num_bins - - def get_tally_xml(self): element = ET.Element("tally") @@ -876,7 +999,10 @@ class Tally(object): nuclides = '' for nuclide in self._nuclides: - nuclides += '{0} '.format(nuclide._name) + if isinstance(nuclide, Nuclide): + nuclides += '{0} '.format(nuclide._name) + else: + nuclides += '{0} '.format(nuclide) subelement = ET.SubElement(element, "nuclides") subelement.text = nuclides.rstrip(' ') @@ -1099,7 +1225,7 @@ class Tally(object): tally_group.create_dataset('scores', data=np.array(self._scores)) # Add a string array of the nuclides to the HDF5 group - nuclides = list() + nuclides = [] for nuclide in self._nuclides: nuclides.append(nuclide._name) @@ -1134,10 +1260,10 @@ class Tally(object): if os.path.exists(filename) and append: tally_results = pickle.load(file(filename, 'rb')) else: - tally_results = dict() + tally_results = {} # Create a nested dictionary within the file for this particular Tally - tally_results['Tally-{0}'.format(self._id)] = dict() + tally_results['Tally-{0}'.format(self._id)] = {} tally_group = tally_results['Tally-{0}'.format(self._id)] # Add basic Tally data to the nested dictionary @@ -1147,7 +1273,7 @@ class Tally(object): tally_group['scores'] = np.array(self._scores) # Add a string array of the nuclides to the HDF5 group - nuclides = list() + nuclides = [] for nuclide in self._nuclides: nuclides.append(nuclide._name) @@ -1155,7 +1281,7 @@ class Tally(object): tally_group['nuclides']= np.array(nuclides) # Create a nested dictionary for the Filters - tally_group['filters'] = dict() + tally_group['filters'] = {} filter_group = tally_group['filters'] for filter in self._filters: @@ -1176,8 +1302,8 @@ class TalliesFile(object): def __init__(self): # Initialize TalliesFile class attributes - self._tallies = list() - self._meshes = list() + self._tallies = [] + self._meshes = [] self._tallies_file = ET.Element("tallies") diff --git a/src/utils/openmc/universe.py b/src/utils/openmc/universe.py index 77d541e865..810a8b5e70 100644 --- a/src/utils/openmc/universe.py +++ b/src/utils/openmc/universe.py @@ -1,11 +1,9 @@ -#!/usr/bin/env python +import abc +from collections import OrderedDict +from xml.etree import ElementTree as ET import openmc from openmc.checkvalue import * -from xml.etree import ElementTree as ET -from collections import OrderedDict -import numpy as np -import abc ################################################################################ @@ -15,6 +13,10 @@ import abc # A static variable for auto-generated Cell IDs AUTO_CELL_ID = 10000 +# A dictionary for storing IDs of cell elements that have already been written, +# used to optimize the writing process +WRITTEN_IDS = {} + def reset_auto_cell_id(): global AUTO_CELL_ID AUTO_CELL_ID = 10000 @@ -25,48 +27,58 @@ class Cell(object): def __init__(self, cell_id=None, name=''): # Initialize Cell class attributes - self._id = None - self._name = None + self.id = cell_id + self.name= name self._fill = None self._type = None - self._surfaces = dict() + self._surfaces = {} self._rotation = None self._translation = None self._offset = None - self.set_id(cell_id) - self.set_name(name) + + @property + def id(self): + return self._id - def get_offset(self, path, filter_offset): - - # Get the current element and remove it from the list - cell_id = path[0] - path = path[1:] - - # If the Cell is filled by a Material - if self._type == 'normal': - offset = 0 - - # If the Cell is filled by a Universe - elif self._type == 'fill': - offset = self._offset[filter_offset-1] - offset += self._fill.get_offset(path, filter_offset) - - # If the Cell is filled by a Lattice - else: - offset = self._fill.get_offset(path, filter_offset) - - return offset - - # Make a recursive call to the Universe filling this Cell - offset = self._cells[cell_id].get_offset(path, filter_offset) - - # Return the offset computed at all nested Universe levels - return offset + @property + def name(self): + return self._name - def set_id(self, cell_id=None): + @property + def fill(self): + return self._fill + + + @property + def type(self): + return self._fill + + + @property + def surfaces(self): + return self._surfaces + + + @property + def rotation(self): + return self._rotation + + + @property + def translation(self): + return self._translation + + + @property + def offset(self): + return self._offset + + + @id.setter + def id(self, cell_id): if cell_id is None: global AUTO_CELL_ID @@ -87,7 +99,8 @@ class Cell(object): self._id = cell_id - def set_name(self, name): + @name.setter + def name(self, name): if not isinstance(name, str): msg = 'Unable to set name for Cell ID={0} with a non-string ' \ @@ -98,7 +111,8 @@ class Cell(object): self._name = name - def set_fill(self, fill): + @fill.setter + def fill(self, fill): if not isinstance(fill, (openmc.Material, Universe, Lattice)) \ and fill != 'void': @@ -118,6 +132,67 @@ class Cell(object): self._type = 'normal' + @rotation.setter + def rotation(self, rotation): + + if not isinstance(rotation, (tuple, list, np.ndarray)): + msg = 'Unable to add rotation {0} to Cell ID={1} since ' \ + 'it is not a Python tuple/list or NumPy ' \ + 'array'.format(rotation, self._id) + raise ValueError(msg) + + elif len(rotation) != 3: + msg = 'Unable to add rotation {0} to Cell ID={1} since ' \ + 'it does not contain 3 values'.format(rotation, self._id) + raise ValueError(msg) + + for axis in rotation: + + if not is_integer(axis) and not is_float(axis): + msg = 'Unable to add rotation {0} to Cell ID={1} since ' \ + 'it is not an integer or floating point ' \ + 'value'.format(axis, self._id) + raise ValueError(msg) + + self._rotation = rotation + + + @translation.setter + def translation(self, translation): + + if not isinstance(translation, (tuple, list, np.ndarray)): + msg = 'Unable to add translation {0} to Cell ID={1} since ' \ + 'it is not a Python tuple/list or NumPy ' \ + 'array'.format(translation, self._id) + raise ValueError(msg) + + elif len(translation) != 3: + msg = 'Unable to add translation {0} to Cell ID={1} since ' \ + 'it does not contain 3 values'.format(translation, self._id) + raise ValueError(msg) + + for axis in translation: + if not is_integer(axis) and not is_float(axis): + msg = 'Unable to add translation {0} to Cell ID={1} since ' \ + 'it is not an integer or floating point ' \ + 'value'.format(axis, self._id) + raise ValueError(msg) + + self._translation = translation + + + @offset.setter + def offset(self, offset): + + if not isinstance(offset, (tuple, list, np.ndarray)): + msg = 'Unable to set offset {0} to Cell ID={1} since ' \ + 'it is not a Python tuple/list or NumPy ' \ + 'array'.format(offset, self._id) + raise ValueError(msg) + + self._offset = offset + + def add_surface(self, surface, halfspace): if not isinstance(surface, openmc.Surface): @@ -165,66 +240,37 @@ class Cell(object): self.add_surface(surface) - def set_rotation(self, rotation): + def get_offset(self, path, filter_offset): - if not isinstance(rotation, (tuple, list, np.ndarray)): - msg = 'Unable to add rotation {0} to Cell ID={1} since ' \ - 'it is not a Python tuple/list or NumPy ' \ - 'array'.format(rotation, self._id) - raise ValueError(msg) + # Get the current element and remove it from the list + cell_id = path[0] + path = path[1:] - elif len(rotation) != 3: - msg = 'Unable to add rotation {0} to Cell ID={1} since ' \ - 'it does not contain 3 values'.format(rotation, self._id) - raise ValueError(msg) + # If the Cell is filled by a Material + if self._type == 'normal': + offset = 0 - for axis in rotation: + # If the Cell is filled by a Universe + elif self._type == 'fill': + offset = self._offset[filter_offset-1] + offset += self._fill.get_offset(path, filter_offset) - if not is_integer(axis) and not is_float(axis): - msg = 'Unable to add rotation {0} to Cell ID={1} since ' \ - 'it is not an integer or floating point ' \ - 'value'.format(axis, self._id) - raise ValueError(msg) + # If the Cell is filled by a Lattice + else: + offset = self._fill.get_offset(path, filter_offset) - self._rotation = rotation + return offset + # Make a recursive call to the Universe filling this Cell + offset = self._cells[cell_id].get_offset(path, filter_offset) - def set_translation(self, translation): - - if not isinstance(translation, (tuple, list, np.ndarray)): - msg = 'Unable to add translation {0} to Cell ID={1} since ' \ - 'it is not a Python tuple/list or NumPy ' \ - 'array'.format(translation, self._id) - raise ValueError(msg) - - elif len(translation) != 3: - msg = 'Unable to add translation {0} to Cell ID={1} since ' \ - 'it does not contain 3 values'.format(translation, self._id) - raise ValueError(msg) - - for axis in translation: - if not is_integer(axis) and not is_float(axis): - msg = 'Unable to add translation {0} to Cell ID={1} since ' \ - 'it is not an integer or floating point ' \ - 'value'.format(axis, self._id) - raise ValueError(msg) - - self._translation = translation - - def set_offset(self, offset): - - if not isinstance(offset, (tuple, list, np.ndarray)): - msg = 'Unable to set offset {0} to Cell ID={1} since ' \ - 'it is not a Python tuple/list or NumPy ' \ - 'array'.format(offset, self._id) - raise ValueError(msg) - - self._offset = offset + # Return the offset computed at all nested Universe levels + return offset def get_all_nuclides(self): - nuclides = dict() + nuclides = {} if self._type != 'void': nuclides.update(self._fill.get_all_nuclides()) @@ -234,7 +280,7 @@ class Cell(object): def get_all_cells(self): - cells = dict() + cells = {} if self._type == 'fill' or self._type == 'lattice': cells.update(self._fill.get_all_cells()) @@ -244,7 +290,7 @@ class Cell(object): def get_all_universes(self): - universes = dict() + universes = {} if self._type == 'fill': universes[self._fill._id] = self._fill @@ -372,23 +418,36 @@ class Universe(object): def __init__(self, universe_id=None, name=''): # Initialize Cell class attributes - self._id = None - self._name = None + self.id = universe_id + self.name = name # Keys - Cell IDs # Values - Cells - self._cells = dict() + self._cells = {} # Keys - Cell IDs # Values - Offsets self._cell_offsets = OrderedDict() self._num_regions = 0 - self.set_id(universe_id) - self.set_name(name) + + @property + def id(self): + return self._id - def set_id(self, universe_id=None): + @property + def name(self): + return self._name + + + @property + def cells(self): + return self._cells + + + @id.setter + def id(self, universe_id): if universe_id is None: global AUTO_UNIVERSE_ID @@ -410,7 +469,8 @@ class Universe(object): self._id = universe_id - def set_name(self, name): + @name.setter + def name(self, name): if not is_string(name): msg = 'Unable to set name for Universe ID={0} with a non-string ' \ @@ -480,7 +540,7 @@ class Universe(object): def get_all_nuclides(self): - nuclides = dict() + nuclides = {} # Append all Nuclides in each Cell in the Universe to the dictionary for cell_id, cell in self._cells.items(): @@ -491,7 +551,7 @@ class Universe(object): def get_all_cells(self): - cells = dict() + cells = {} # Add this Universe's cells to the dictionary cells.update(self._cells) @@ -508,7 +568,7 @@ class Universe(object): # Get all Cells in this Universe cells = self.get_all_cells() - universes = dict() + universes = {} # Append all Universes containing each Cell to the dictionary for cell_id, cell in cells.items(): @@ -534,12 +594,9 @@ class Universe(object): # Iterate over all Cells for cell_id, cell in self._cells.items(): - # Determine if XML element already contains subelement for this Cell - path = './cell[@id=\'{0}\']'.format(cell_id) - test = xml_element.find(path) - - # If the element does not contain the Cell subelement, then add it - if test is None: + # If the cell was not already written, write it + if not cell_id in WRITTEN_IDS: + WRITTEN_IDS[cell_id] = None # Create XML subelement for this Cell cell_subelement = cell.create_xml_subelement(xml_element) @@ -563,17 +620,40 @@ class Lattice(object): def __init__(self, lattice_id=None, name=''): # Initialize Lattice class attributes - self._id = None - self._name = None + self.id = lattice_id + self.name = name self._pitch = None self._outer = None self._universes = None - self.set_id(lattice_id) - self.set_name(name) + + @property + def id(self): + return self._id - def set_id(self, lattice_id=None): + @property + def name(self): + return self._name + + + @property + def pitch(self): + return self._pitch + + + @property + def outer(self): + return self._outer + + + @property + def universes(self): + return self._universes + + + @id.setter + def id(self, lattice_id): if lattice_id is None: global AUTO_UNIVERSE_ID @@ -594,7 +674,8 @@ class Lattice(object): self._id = lattice_id - def set_name(self, name): + @name.setter + def name(self, name): if not is_string(name): msg = 'Unable to set name for Lattice ID={0} with a non-string ' \ @@ -605,7 +686,8 @@ class Lattice(object): self._name = name - def set_outer(self, outer): + @outer.setter + def outer(self, outer): if not isinstance(outer, Universe): msg = 'Unable to set Lattice ID={0} outer universe to {1} ' \ @@ -615,7 +697,8 @@ class Lattice(object): self._outer = outer - def set_universes(self, universes): + @universes.setter + def universes(self, universes): if not isinstance(universes, (tuple, list, np.ndarray)): msg = 'Unable to set Lattice ID={0} universes to {1} since ' \ @@ -629,7 +712,7 @@ class Lattice(object): def get_unique_universes(self): unique_universes = np.unique(self._universes.ravel()) - universes = dict() + universes = {} for universe in unique_universes: universes[universe._id] = universe @@ -639,7 +722,7 @@ class Lattice(object): def get_all_nuclides(self): - nuclides = dict() + nuclides = {} # Get all unique Universes contained in each of the lattice cells unique_universes = self.get_unique_universes() @@ -653,7 +736,7 @@ class Lattice(object): def get_all_cells(self): - cells = dict() + cells = {} unique_universes = self.get_unique_universes() for universe_id, universe in unique_universes.items(): @@ -666,7 +749,7 @@ class Lattice(object): # Initialize a dictionary of all Universes contained by the Lattice # in each nested Universe level - all_universes = dict() + all_universes = {} # Get all unique Universes contained in each of the lattice cells unique_universes = self.get_unique_universes() @@ -700,7 +783,23 @@ class RectLattice(Lattice): self._offsets = None - def set_dimension(self, dimension): + @property + def dimension(self): + return self._dimension + + + @property + def lower_left(self): + return self._lower_left + + + @property + def offsets(self): + return self._offsets + + + @dimension.setter + def dimension(self, dimension): if not isinstance(dimension, (tuple, list, np.ndarray)): msg = 'Unable to set RectLattice ID={0} dimension to {1} since ' \ @@ -731,7 +830,8 @@ class RectLattice(Lattice): self._dimension = dimension - def set_lower_left(self, lower_left): + @lower_left.setter + def lower_left(self, lower_left): if not isinstance(lower_left, (tuple, list, np.ndarray)): msg = 'Unable to set RectLattice ID={0} lower_left to {1} since ' \ @@ -757,7 +857,8 @@ class RectLattice(Lattice): self._lower_left = lower_left - def set_offsets(self, offsets): + @offsets.setter + def offsets(self, offsets): if not isinstance(offsets, (tuple, list, np.ndarray)): msg = 'Unable to set Lattice ID={0} offsets to {1} since ' \ @@ -768,7 +869,8 @@ class RectLattice(Lattice): self._offsets = offsets - def set_pitch(self, pitch): + @Lattice.pitch.setter + def pitch(self, pitch): if not isinstance(pitch, (tuple, list, np.ndarray)): msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \ @@ -985,7 +1087,23 @@ class HexLattice(Lattice): self._center = None - def set_num_rings(self, num_rings): + @property + def num_rings(self): + return self._num_rings + + + @property + def num_axial(self): + return self._num_axial + + + @property + def center(self): + return self._center + + + @num_rings.setter + def num_rings(self, num_rings): if not is_integer(num_rings) and num_rings < 1: msg = 'Unable to set HexLattice ID={0} number of rings to {1} ' \ @@ -995,7 +1113,8 @@ class HexLattice(Lattice): self._num_rings = num_rings - def set_num_axial(self, num_axial): + @num_axial.setter + def num_axial(self, num_axial): if not is_integer(num_axial) and num_axial < 1: msg = 'Unable to set HexLattice ID={0} number of axial to {1} ' \ @@ -1005,7 +1124,8 @@ class HexLattice(Lattice): self._num_axial = num_axial - def set_center(self, center): + @center.setter + def center(self, center): if not isinstance(center, (tuple, list, np.ndarray)): msg = 'Unable to set HexLattice ID={0} dimension to {1} since ' \ @@ -1030,7 +1150,8 @@ class HexLattice(Lattice): self._center = center - def set_pitch(self, pitch): + @Lattice.pitch.setter + def pitch(self, pitch): if not isinstance(pitch, (tuple, list, np.ndarray)): msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \ @@ -1060,8 +1181,11 @@ class HexLattice(Lattice): self._pitch = pitch - def set_universes(self, universes): - super(HexLattice, self).set_universes(universes) + @Lattice.universes.setter + def universes(self, universes): + + # Call Lattice.universes parent class setter property + Lattice.universes.fset(self, universes) # NOTE: This routine assumes that the user creates a "ragged" list of # lists, where each sub-list corresponds to one ring of Universes. @@ -1084,7 +1208,7 @@ class HexLattice(Lattice): # Set the number of axial positions. if n_dims == 3: - self.set_num_axial(self._universes.shape[0]) + self.num_axial = self._universes.shape[0] else: self._num_axial = None @@ -1092,7 +1216,7 @@ class HexLattice(Lattice): # Set the number of rings and make sure this number is consistent for # all axial positions. if n_dims == 3: - self.set_num_rings(len(self._universes[0])) + self.num_rings = len(self._universes[0]) for rings in self._universes: if len(rings) != self._num_rings: msg = 'HexLattice ID={0:d} has an inconsistent number of ' \ @@ -1100,7 +1224,7 @@ class HexLattice(Lattice): raise ValueError(msg) else: - self.set_num_rings(self._universes.shape[0]) + self.num_rings = self._universes.shape[0] # Make sure there are the correct number of elements in each ring. if n_dims == 3: diff --git a/src/utils/setup.py b/src/utils/setup.py index ad07cc6b4e..72682aea9d 100644 --- a/src/utils/setup.py +++ b/src/utils/setup.py @@ -3,7 +3,7 @@ from distutils.core import setup setup(name='openmc', - version='0.6.0', + version='0.6.2', description='OpenMC Python API', author='Will Boyd', author_email='wbinventor@gmail.com', diff --git a/tests/run_tests.py b/tests/run_tests.py index 531225d9ad..42fb846d88 100755 --- a/tests/run_tests.py +++ b/tests/run_tests.py @@ -179,6 +179,10 @@ class Test(object): # Runs cmake when in non-script mode def run_cmake(self): os.environ['FC'] = self.fc + if self.petsc: + os.environ['PETSC_DIR'] = PETSC_DIR + if self.mpi: + os.environ['MPI_DIR'] = MPI_DIR build_opts = self.build_opts.split() self.cmake += build_opts rc = call(self.cmake) diff --git a/tests/test_score_total_yn/results_true.dat b/tests/test_score_total_yn/results_true.dat index b78df10085..36d7040bd9 100644 --- a/tests/test_score_total_yn/results_true.dat +++ b/tests/test_score_total_yn/results_true.dat @@ -60,6 +60,106 @@ tally 2: 0.000000E+00 0.000000E+00 0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +8.795501E-01 +1.600341E-01 +-1.530661E-02 +5.215322E-04 +1.571095E-02 +9.166057E-04 +7.226989E-03 +3.683499E-04 +-2.456165E-02 +1.571441E-04 +1.387449E-02 +5.600648E-04 +-2.186870E-02 +2.081164E-04 +-1.106814E-02 +1.065144E-04 +-4.579593E-03 +2.008225E-04 +7.573253E-03 +2.493075E-04 +-1.505016E-02 +1.451400E-04 +1.503792E-02 +9.539029E-05 +-1.183668E-02 +1.741393E-04 +4.060912E-03 +5.890591E-05 +1.789747E-02 +8.239749E-05 +-2.168438E-02 +1.555757E-04 +-1.045826E-02 +8.108402E-05 +1.227413E-02 +2.502871E-04 +-2.806122E-02 +1.886120E-04 +-4.918718E-03 +2.129038E-04 +-1.014729E-02 +6.914145E-05 +-5.873760E-03 +2.229738E-04 +6.666510E-03 +1.394147E-04 +-3.245528E-03 +1.550876E-04 +-1.037659E-02 +2.163837E-04 1.517577E+01 4.747271E+01 -2.463504E-01 @@ -110,6 +210,56 @@ tally 2: 2.383142E-02 -1.463402E-01 1.526988E-02 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 3.151504E+00 2.051857E+00 -4.923938E-02 @@ -160,6 +310,56 @@ tally 2: 1.674419E-04 -1.682320E-02 8.389254E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 4.536316E+01 4.258781E+02 -6.747275E-01 diff --git a/tests/test_score_total_yn/tallies.xml b/tests/test_score_total_yn/tallies.xml index df6de1c134..ca5dcd2dcb 100644 --- a/tests/test_score_total_yn/tallies.xml +++ b/tests/test_score_total_yn/tallies.xml @@ -9,6 +9,7 @@ total-y4 + U-235 total diff --git a/tests/test_union_energy_grids/geometry.xml b/tests/test_union_energy_grids/geometry.xml new file mode 100644 index 0000000000..612e46132e --- /dev/null +++ b/tests/test_union_energy_grids/geometry.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/tests/test_union_energy_grids/materials.xml b/tests/test_union_energy_grids/materials.xml new file mode 100644 index 0000000000..45ccc65540 --- /dev/null +++ b/tests/test_union_energy_grids/materials.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/tests/test_union_energy_grids/results.py b/tests/test_union_energy_grids/results.py new file mode 100644 index 0000000000..b616e2874b --- /dev/null +++ b/tests/test_union_energy_grids/results.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python + +import sys + +sys.path.insert(0, '../../src/utils') +from openmc.statepoint import StatePoint + +# read in statepoint file +if len(sys.argv) > 1: + print(sys.argv) + sp = StatePoint(sys.argv[1]) +else: + sp = StatePoint('statepoint.10.binary') +sp.read_results() + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_union_energy_grids/results_true.dat b/tests/test_union_energy_grids/results_true.dat new file mode 100644 index 0000000000..05cda2e980 --- /dev/null +++ b/tests/test_union_energy_grids/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +3.215828E-01 2.966835E-03 diff --git a/tests/test_union_energy_grids/settings.xml b/tests/test_union_energy_grids/settings.xml new file mode 100644 index 0000000000..1eb22241cc --- /dev/null +++ b/tests/test_union_energy_grids/settings.xml @@ -0,0 +1,18 @@ + + + + union + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + diff --git a/tests/test_union_energy_grids/test_union_energy_grids.py b/tests/test_union_energy_grids/test_union_energy_grids.py new file mode 100644 index 0000000000..6fdbf87459 --- /dev/null +++ b/tests/test_union_energy_grids/test_union_energy_grids.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + try: + test_run() + test_created_statepoint() + test_results() + finally: + teardown()