From 8e4ce676ff56359cf962e016d0c8c3e5013319d5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 12 Jul 2022 13:50:48 -0500 Subject: [PATCH 01/83] Make sure weight window lower/upper bounds are flattened before writing to XML --- openmc/weight_windows.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 2e5b7d8e0..257a17e9a 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -291,10 +291,10 @@ class WeightWindows(IDManagerMixin): subelement.text = ' '.join(str(e) for e in self.energy_bounds) subelement = ET.SubElement(element, 'lower_ww_bounds') - subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds) + subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds.ravel()) subelement = ET.SubElement(element, 'upper_ww_bounds') - subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds) + subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds.ravel()) subelement = ET.SubElement(element, 'survival_ratio') subelement.text = str(self.survival_ratio) From cd60109a77dd7a3e6a492ef734684af272fb5cce Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 13 Jul 2022 19:14:14 +0100 Subject: [PATCH 02/83] added bounding_region --- openmc/universe.py | 43 +++++++++++++++++++ .../regression_tests/dagmc/universes/test.py | 14 ++++++ 2 files changed, 57 insertions(+) diff --git a/openmc/universe.py b/openmc/universe.py index b65263ef9..2fd931506 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -713,6 +713,49 @@ class DAGMCUniverse(UniverseBase): dagmc_element.set('filename', self.filename) xml_element.append(dagmc_element) + def bounding_region(self, bounded_type='box', boundary_type='vacuum'): + """Creates a either a spherical or box shaped bounding region around + the DAGMC geometry. + Parameters + ---------- + bounded_type : str + The type of bounding surface(s) to use when constructing the region. + Options include a single spherical surface (sphere) or a rectangle + made from siz planes (box). + boundary_type : str + Boundary condition that defines the behavior for particles hitting + the surface. Defaults to vacuum boundary condition. Passed into + Returns + ------- + openmc.Region + Region instance + """ + + bounding_box = self.bounding_box + + if bounded_type == 'sphere': + import math + bounding_box_center = (bounding_box[0] + bounding_box[1])/2 + radius = math.dist(bounding_box[0], bounding_box[1]) + bounding_surface = openmc.Sphere(x0=bounding_box_center[0], y0=bounding_box_center[1], z0=bounding_box_center[2], r=radius) + + return -bounding_surface + + if bounded_type == 'box': + # defines plane surfaces for all six faces of the bounding box + lower_x = openmc.XPlane(bounding_box[0][0], boundary_type=boundary_type) + upper_x = openmc.XPlane(bounding_box[1][0], boundary_type=boundary_type) + lower_y = openmc.YPlane(bounding_box[0][1], boundary_type=boundary_type) + upper_y = openmc.YPlane(bounding_box[1][1], boundary_type=boundary_type) + lower_z = openmc.ZPlane(bounding_box[0][2], boundary_type=boundary_type) + upper_z = openmc.ZPlane(bounding_box[1][2], boundary_type=boundary_type) + + return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z + + def bounded_universe(self, **kwargs): + bounding_cell = openmc.Cell(fill=self, region=self.bounding_region(**kwargs)) + return openmc.Universe(cells=[bounding_cell]) + @classmethod def from_hdf5(cls, group): """Create DAGMC universe from HDF5 group diff --git a/tests/regression_tests/dagmc/universes/test.py b/tests/regression_tests/dagmc/universes/test.py index 6726b4b44..f6964857d 100644 --- a/tests/regression_tests/dagmc/universes/test.py +++ b/tests/regression_tests/dagmc/universes/test.py @@ -51,6 +51,20 @@ class DAGMCUniverseTest(PyAPITestHarness): assert bounding_box[0].tolist() == [-25., -25., -25.] assert bounding_box[1].tolist() == [25., 25., 25.] + # checks that the bounding region is six surfaces each with a vacuum boundary type + b_region = pincell_univ.bounding_region(bounded_type='box', boundary_type='vacuum') + assert isinstance(b_region, openmc.Region) + assert len(b_region.get_surfaces()) == 6 + for surface in list(b_region.get_surfaces().values()): + assert surface.boundary_type == 'vacuum' + + # checks that the bounding region is a single surface with a reflective boundary type + b_region = pincell_univ.bounding_region(bounded_type='sphere', boundary_type='reflective') + assert isinstance(b_region, openmc.Region) + assert len(b_region.get_surfaces()) == 1 + for surface in list(b_region.get_surfaces().values()): + assert surface.boundary_type == 'reflective' + # create a 2 x 2 lattice using the DAGMC pincell pitch = np.asarray((24.0, 24.0)) lattice = openmc.RectLattice() From 227c4e330fa2a339afa9cb3e5fb5f536683f3689 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 14 Jul 2022 09:36:17 +0100 Subject: [PATCH 03/83] improved doc strings added checks --- openmc/universe.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/openmc/universe.py b/openmc/universe.py index 2fd931506..06807bab3 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -12,9 +12,12 @@ import numpy as np import openmc import openmc.checkvalue as cv + from ._xml import get_text +from .checkvalue import check_type, check_value from .mixin import IDManagerMixin from .plots import _SVG_COLORS +from .surface import _BOUNDARY_TYPES class UniverseBase(ABC, IDManagerMixin): @@ -724,13 +727,17 @@ class DAGMCUniverse(UniverseBase): made from siz planes (box). boundary_type : str Boundary condition that defines the behavior for particles hitting - the surface. Defaults to vacuum boundary condition. Passed into + the surface. Defaults to vacuum boundary condition. Passed into the + surface construction. Returns ------- openmc.Region Region instance """ + check_type('boundary type', boundary_type, str) + check_value('boundary type', boundary_type, _BOUNDARY_TYPES) + bounding_box = self.bounding_box if bounded_type == 'sphere': @@ -753,6 +760,10 @@ class DAGMCUniverse(UniverseBase): return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z def bounded_universe(self, **kwargs): + """Returns an openmc.Universe filled with the DAGMCUniverse and bounded + with a cell. Defaults to a box cell with a vacuum surface. kwargs are + passed directly to DAGMCUniverse.bounding_region()""" + bounding_cell = openmc.Cell(fill=self, region=self.bounding_region(**kwargs)) return openmc.Universe(cells=[bounding_cell]) From 28e203c8286acd83ccd97ce23d4c8e19b4ceb559 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 14 Jul 2022 15:10:27 +0100 Subject: [PATCH 04/83] typos fixed in review by @pshriwise Co-authored-by: Patrick Shriwise --- openmc/universe.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 06807bab3..7ecf489ce 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -724,7 +724,7 @@ class DAGMCUniverse(UniverseBase): bounded_type : str The type of bounding surface(s) to use when constructing the region. Options include a single spherical surface (sphere) or a rectangle - made from siz planes (box). + made from six planes (box). boundary_type : str Boundary condition that defines the behavior for particles hitting the surface. Defaults to vacuum boundary condition. Passed into the @@ -760,7 +760,7 @@ class DAGMCUniverse(UniverseBase): return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z def bounded_universe(self, **kwargs): - """Returns an openmc.Universe filled with the DAGMCUniverse and bounded + """Returns an openmc.Universe filled with this DAGMCUniverse and bounded with a cell. Defaults to a box cell with a vacuum surface. kwargs are passed directly to DAGMCUniverse.bounding_region()""" From 839b862d5645be785fd1f79dc6c3ba3f19e39b9e Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 14 Jul 2022 15:15:04 +0100 Subject: [PATCH 05/83] added type value check for bounded type --- openmc/universe.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/universe.py b/openmc/universe.py index 7ecf489ce..e46d4d027 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -737,6 +737,8 @@ class DAGMCUniverse(UniverseBase): check_type('boundary type', boundary_type, str) check_value('boundary type', boundary_type, _BOUNDARY_TYPES) + check_type('bounded type', bounded_type, str) + check_value('bounded type', bounded_type, ('box', 'sphere')) bounding_box = self.bounding_box From 4c62517b1aa4cbc4b1c742c97e39473551be4682 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 14 Jul 2022 16:50:33 +0100 Subject: [PATCH 06/83] added test to check bounded_universe is usable in model (should fail) --- tests/regression_tests/dagmc/universes/test.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/regression_tests/dagmc/universes/test.py b/tests/regression_tests/dagmc/universes/test.py index f6964857d..7abd80ef7 100644 --- a/tests/regression_tests/dagmc/universes/test.py +++ b/tests/regression_tests/dagmc/universes/test.py @@ -46,6 +46,13 @@ class DAGMCUniverseTest(PyAPITestHarness): # create the DAGMC universe pincell_univ = openmc.DAGMCUniverse(filename='dagmc.h5m', auto_geom_ids=True) + # creates another DAGMC universe, this time with within a bounded cell + bound_pincell_cell = openmc.DAGMCUniverse(filename='dagmc.h5m').bounded_universe() + # uses the bound_dag_cell as the root argument to test the type checks in openmc.Geometry + bound_pincell_geometry = openmc.Geometry(root=[bound_pincell_cell]) + # assigns the bound_dag_geometry to the model to test the type checks in model.Geometry setter + model.Geometry = bound_pincell_geometry + # checks that the bounding box is calculated correctly bounding_box = pincell_univ.bounding_box assert bounding_box[0].tolist() == [-25., -25., -25.] From 3d091a6cfbdb24242886e7a8d076852b8b085d4d Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 14 Jul 2022 16:51:22 +0100 Subject: [PATCH 07/83] returning Cell instead of Universe --- openmc/universe.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index e46d4d027..9931724fb 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -762,12 +762,18 @@ class DAGMCUniverse(UniverseBase): return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z def bounded_universe(self, **kwargs): - """Returns an openmc.Universe filled with this DAGMCUniverse and bounded - with a cell. Defaults to a box cell with a vacuum surface. kwargs are - passed directly to DAGMCUniverse.bounding_region()""" + """Returns an openmc.Cell that bounds the DAGMCUniverse and is filled + with the DAGMCUniverse. Defaults to a box cell with a vacuum surface. + kwargs are passed directly to DAGMCUniverse.bounding_region() + + Returns + ------- + openmc.Cell + cell filled with the DAGMCUniverse + """ bounding_cell = openmc.Cell(fill=self, region=self.bounding_region(**kwargs)) - return openmc.Universe(cells=[bounding_cell]) + return bounding_cell @classmethod def from_hdf5(cls, group): From c3e60da1b9f2c04545200e8d5f16360304c166ca Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 14 Jul 2022 16:55:26 +0100 Subject: [PATCH 08/83] testing and returing openmc.Universe --- openmc/universe.py | 12 ++++++------ tests/regression_tests/dagmc/universes/test.py | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 9931724fb..811b27b31 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -762,18 +762,18 @@ class DAGMCUniverse(UniverseBase): return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z def bounded_universe(self, **kwargs): - """Returns an openmc.Cell that bounds the DAGMCUniverse and is filled - with the DAGMCUniverse. Defaults to a box cell with a vacuum surface. - kwargs are passed directly to DAGMCUniverse.bounding_region() + """Returns an openmc.Universe filled with this DAGMCUniverse and bounded + with a cell. Defaults to a box cell with a vacuum surface. kwargs are + passed directly to DAGMCUniverse.bounding_region() Returns ------- - openmc.Cell - cell filled with the DAGMCUniverse + openmc.Universe + Universe instance """ bounding_cell = openmc.Cell(fill=self, region=self.bounding_region(**kwargs)) - return bounding_cell + return openmc.Universe(cells=[bounding_cell]) @classmethod def from_hdf5(cls, group): diff --git a/tests/regression_tests/dagmc/universes/test.py b/tests/regression_tests/dagmc/universes/test.py index 7abd80ef7..f2eb3882f 100644 --- a/tests/regression_tests/dagmc/universes/test.py +++ b/tests/regression_tests/dagmc/universes/test.py @@ -47,9 +47,9 @@ class DAGMCUniverseTest(PyAPITestHarness): pincell_univ = openmc.DAGMCUniverse(filename='dagmc.h5m', auto_geom_ids=True) # creates another DAGMC universe, this time with within a bounded cell - bound_pincell_cell = openmc.DAGMCUniverse(filename='dagmc.h5m').bounded_universe() + bound_pincell_universe = openmc.DAGMCUniverse(filename='dagmc.h5m').bounded_universe() # uses the bound_dag_cell as the root argument to test the type checks in openmc.Geometry - bound_pincell_geometry = openmc.Geometry(root=[bound_pincell_cell]) + bound_pincell_geometry = openmc.Geometry(root=bound_pincell_universe) # assigns the bound_dag_geometry to the model to test the type checks in model.Geometry setter model.Geometry = bound_pincell_geometry From 34e4a5efd12b8adc21ae9c3bb83fc073fc95df71 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 14 Jul 2022 18:06:07 +0100 Subject: [PATCH 09/83] setting auto_geom_ids to true for bounded universe --- openmc/universe.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 811b27b31..989b50c73 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -761,10 +761,18 @@ class DAGMCUniverse(UniverseBase): return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z - def bounded_universe(self, **kwargs): + def bounded_universe(self, auto_geom_ids=True, **kwargs): """Returns an openmc.Universe filled with this DAGMCUniverse and bounded - with a cell. Defaults to a box cell with a vacuum surface. kwargs are - passed directly to DAGMCUniverse.bounding_region() + with a cell. Defaults to a box cell with a vacuum surface however this + can be changed using the kwargs which are passed directly to + DAGMCUniverse.bounding_region(). + + Parameters + ---------- + auto_geom_ids : bool + Set IDs automatically on initialization (True) or report overlaps + in ID space between CSG and DAGMC (False). Defaults to True to avoid + overlapping ID numbers between the CSG and DAGMC geometry ID numbers Returns ------- @@ -772,6 +780,7 @@ class DAGMCUniverse(UniverseBase): Universe instance """ + self.auto_geom_ids = auto_geom_ids bounding_cell = openmc.Cell(fill=self, region=self.bounding_region(**kwargs)) return openmc.Universe(cells=[bounding_cell]) From 953b35e742d4526b4bf4c5429d876f21506ec5c3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Jul 2022 07:21:08 -0500 Subject: [PATCH 10/83] Add openmc_sample_external_source function and Python binding --- openmc/lib/core.py | 39 ++++++++++++++++++++++++++++++++++++++- openmc/lib/math.py | 11 +++++------ src/source.cpp | 10 ++++++++++ 3 files changed, 53 insertions(+), 7 deletions(-) diff --git a/openmc/lib/core.py b/openmc/lib/core.py index de5f4adf2..3f4a301cd 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -1,8 +1,10 @@ from contextlib import contextmanager from ctypes import (c_bool, c_int, c_int32, c_int64, c_double, c_char_p, - c_char, POINTER, Structure, c_void_p, create_string_buffer) + c_char, POINTER, Structure, c_void_p, create_string_buffer, + c_uint64) import sys import os +from random import getrandbits import numpy as np from numpy.ctypeslib import as_array @@ -10,6 +12,7 @@ from numpy.ctypeslib import as_array from . import _dll from .error import _error_handler import openmc.lib +import openmc class _SourceSite(Structure): @@ -95,6 +98,9 @@ _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double), POINTER(c_double)] _dll.openmc_global_bounding_box.restype = c_int _dll.openmc_global_bounding_box.errcheck = _error_handler +_dll.openmc_sample_external_source.argtypes = [POINTER(c_uint64), POINTER(_SourceSite)] +_dll.openmc_sample_external_source.restype = c_int +_dll.openmc_sample_external_source.errcheck = _error_handler def global_bounding_box(): @@ -415,6 +421,37 @@ def run(output=True): _dll.openmc_run() +def sample_external_source(prn_seed=None): + """Sample external source + + .. versionadded:: 0.13.1 + + Parameters + ---------- + prn_seed : int + PRNG seed; if None, one will be generated randomly + + Returns + ------- + openmc.SourceParticle + Sampled source particle + + """ + if prn_seed is None: + prn_seed = getrandbits(63) + + # Call into C API to sample source + site = _SourceSite() + _dll.openmc_sample_external_source(c_uint64(prn_seed), site) + + # Convert to SourceParticle and return + return openmc.SourceParticle( + r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt, + delayed_group=site.delayed_group, surf_id=site.surf_id, + particle=openmc.ParticleType(site.particle) + ) + + def simulation_init(): """Initialize simulation""" _dll.openmc_simulation_init() diff --git a/openmc/lib/math.py b/openmc/lib/math.py index bdb09b2d1..d70c20a9f 100644 --- a/openmc/lib/math.py +++ b/openmc/lib/math.py @@ -1,12 +1,11 @@ from ctypes import c_int, c_double, POINTER, c_uint64 +from random import getrandbits import numpy as np from numpy.ctypeslib import ndpointer from . import _dll -from random import getrandbits - _dll.t_percentile.restype = c_double _dll.t_percentile.argtypes = [c_double, c_int] @@ -240,10 +239,10 @@ def maxwell_spectrum(T, prn_seed=None): Sampled outgoing energy """ - + if prn_seed is None: prn_seed = getrandbits(63) - + return _dll.maxwell_spectrum(T, c_uint64(prn_seed)) @@ -265,7 +264,7 @@ def watt_spectrum(a, b, prn_seed=None): Sampled outgoing energy """ - + if prn_seed is None: prn_seed = getrandbits(63) @@ -290,7 +289,7 @@ def normal_variate(mean_value, std_dev, prn_seed=None): Sampled outgoing normally distributed value """ - + if prn_seed is None: prn_seed = getrandbits(63) diff --git a/src/source.cpp b/src/source.cpp index 11ede5fc4..8b634ec45 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -399,4 +399,14 @@ void free_memory_source() model::external_sources.clear(); } +//============================================================================== +// C API +//============================================================================== + +extern "C" int openmc_sample_external_source(uint64_t* seed, SourceSite* site) +{ + *site = sample_external_source(seed); + return 0; +} + } // namespace openmc From b1267f05448d80bce3069779f03de68749c44297 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Jul 2022 13:23:07 -0500 Subject: [PATCH 11/83] Allow openmc_sample_external_source to produce multiple samples --- docs/source/pythonapi/capi.rst | 1 + include/openmc/capi.h | 1 + openmc/lib/core.py | 33 ++++++++++++++++++++------------- src/source.cpp | 13 +++++++++++-- 4 files changed, 33 insertions(+), 15 deletions(-) diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 03be25abf..882e3c71b 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -29,6 +29,7 @@ Functions reset run run_in_memory + sample_external_source simulation_init simulation_finalize source_bank diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 0929a11f7..69d3ff1f6 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -121,6 +121,7 @@ int openmc_regular_mesh_set_params(int32_t index, int n, const double* ll, int openmc_reset(); int openmc_reset_timers(); int openmc_run(); +int openmc_sample_external_source(size_t n, uint64_t* seed, void* sites); void openmc_set_seed(int64_t new_seed); int openmc_set_n_batches( int32_t n_batches, bool set_max_batches, bool add_statepoint_batch); diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 3f4a301cd..3a437a295 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -1,7 +1,7 @@ from contextlib import contextmanager from ctypes import (c_bool, c_int, c_int32, c_int64, c_double, c_char_p, c_char, POINTER, Structure, c_void_p, create_string_buffer, - c_uint64) + c_uint64, c_size_t) import sys import os from random import getrandbits @@ -98,7 +98,7 @@ _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double), POINTER(c_double)] _dll.openmc_global_bounding_box.restype = c_int _dll.openmc_global_bounding_box.errcheck = _error_handler -_dll.openmc_sample_external_source.argtypes = [POINTER(c_uint64), POINTER(_SourceSite)] +_dll.openmc_sample_external_source.argtypes = [c_size_t, POINTER(c_uint64), POINTER(_SourceSite)] _dll.openmc_sample_external_source.restype = c_int _dll.openmc_sample_external_source.errcheck = _error_handler @@ -421,35 +421,42 @@ def run(output=True): _dll.openmc_run() -def sample_external_source(prn_seed=None): +def sample_external_source(n_samples=1, prn_seed=None): """Sample external source .. versionadded:: 0.13.1 Parameters ---------- + n_samples : int + Number of samples prn_seed : int PRNG seed; if None, one will be generated randomly Returns ------- - openmc.SourceParticle - Sampled source particle + list of openmc.SourceParticle + List of samples source particles """ + if n_samples <= 0: + raise ValueError("Number of samples must be positive") if prn_seed is None: prn_seed = getrandbits(63) # Call into C API to sample source - site = _SourceSite() - _dll.openmc_sample_external_source(c_uint64(prn_seed), site) + sites_array = (_SourceSite * n_samples)() + _dll.openmc_sample_external_source(c_size_t(n_samples), c_uint64(prn_seed), sites_array) - # Convert to SourceParticle and return - return openmc.SourceParticle( - r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt, - delayed_group=site.delayed_group, surf_id=site.surf_id, - particle=openmc.ParticleType(site.particle) - ) + # Convert to list of SourceParticle and return + return [ + openmc.SourceParticle( + r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt, + delayed_group=site.delayed_group, surf_id=site.surf_id, + particle=openmc.ParticleType(site.particle) + ) + for site in sites_array + ] def simulation_init(): diff --git a/src/source.cpp b/src/source.cpp index 8b634ec45..3716720df 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -403,9 +403,18 @@ void free_memory_source() // C API //============================================================================== -extern "C" int openmc_sample_external_source(uint64_t* seed, SourceSite* site) +extern "C" int openmc_sample_external_source( + size_t n, uint64_t* seed, void* sites) { - *site = sample_external_source(seed); + if (!sites || !seed) { + set_errmsg("Received null pointer."); + return OPENMC_E_INVALID_ARGUMENT; + } + + auto sites_array = static_cast(sites); + for (size_t i = 0; i < n; ++i) { + sites_array[i] = sample_external_source(seed); + } return 0; } From ba5f6f58fad0b28a2c69d4cc63d7dff0ac069067 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Jul 2022 19:57:58 -0500 Subject: [PATCH 12/83] Add test for openmc.lib.sample_external_source --- tests/unit_tests/test_lib.py | 41 ++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index a0f18f851..f1f84bff2 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -807,3 +807,44 @@ def test_cell_rotation(pincell_model_w_univ, mpi_intracomm): cell.rotation = (180., 0., 0.) assert cell.rotation == pytest.approx([180., 0., 0.]) openmc.lib.finalize() + + +def test_sample_external_source(run_in_tmpdir, mpi_intracomm): + # Define a simple model and export + mat = openmc.Material() + mat.add_nuclide('U235', 1.0e-2) + sph = openmc.Sphere(r=100.0, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-sph) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.settings.source = openmc.Source( + space=openmc.stats.Box([-5., -5., -5.], [5., 5., 5.]), + angle=openmc.stats.Monodirectional((0., 0., 1.)), + energy=openmc.stats.Discrete([1.0e5], [1.0]) + ) + model.settings.particles = 1000 + model.settings.batches = 10 + model.export_to_xml() + + # Sample some particles and make sure they match specified source + openmc.lib.init() + particles = openmc.lib.sample_external_source(10, prn_seed=3) + for p in particles: + assert -5. < p.r[0] < 5. + assert -5. < p.r[1] < 5. + assert -5. < p.r[2] < 5. + assert p.u[0] == 0.0 + assert p.u[1] == 0.0 + assert p.u[2] == 1.0 + assert p.E == 1.0e5 + + # Using the same seed should produce the same particles + other_particles = openmc.lib.sample_external_source(10, prn_seed=3) + for p1, p2 in zip(particles, other_particles): + assert p1.r == p2.r + assert p1.u == p2.u + assert p1.E == p2.E + assert p1.time == p2.time + assert p1.wgt == p2.wgt + + openmc.lib.finalize() From 38df4f6928cdef6ef92196b92e79565d1e42c078 Mon Sep 17 00:00:00 2001 From: Richard Morrison Date: Sun, 17 Jul 2022 20:54:02 +0100 Subject: [PATCH 13/83] Fixed std::cout sync bug in output.cpp --- src/output.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/output.cpp b/src/output.cpp index ece475de0..0bd9693a8 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -93,7 +93,7 @@ void title() // Write number of OpenMP threads fmt::print(" OpenMP Threads | {}\n", omp_get_max_threads()); #endif - std::cout << std::endl; + fmt::print("\n"); } //============================================================================== @@ -133,7 +133,7 @@ void header(const char* msg, int level) // Print header based on verbosity level. if (settings::verbosity >= level) - std::cout << '\n' << out << "\n" << std::endl; + fmt::print("\n{}\n\n", out); } //============================================================================== @@ -379,7 +379,7 @@ void print_generation() if (n > 1) { fmt::print(" {:8.5f} +/-{:8.5f}", simulation::keff, simulation::keff_std); } - std::cout << std::endl; + fmt::print("\n"); } //============================================================================== @@ -546,7 +546,7 @@ void print_results() fmt::print(" Leakage Fraction = {:.5f}\n", gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n); } - std::cout << std::endl; + fmt::print("\n"); } //============================================================================== From a42ce82396ed3f3e40718894a3208b3b3bbfc6a3 Mon Sep 17 00:00:00 2001 From: Richard Morrison Date: Sun, 17 Jul 2022 21:11:17 +0100 Subject: [PATCH 14/83] ARG & ENV scoping fix in Dockerfile --- Dockerfile | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/Dockerfile b/Dockerfile index 341390c9a..c83460f35 100644 --- a/Dockerfile +++ b/Dockerfile @@ -27,10 +27,6 @@ ARG compile_cores=1 # Set default value of HOME to /root ENV HOME=/root -# OpenMC variables -ARG openmc_branch=master -ENV OPENMC_REPO='https://github.com/openmc-dev/openmc' - # Embree variables ENV EMBREE_TAG='v3.12.2' ENV EMBREE_REPO='https://github.com/embree/embree' @@ -174,6 +170,17 @@ RUN if [ "$build_libmesh" = "on" ]; then \ FROM dependencies AS build +ENV HOME=/root + +ARG openmc_branch=master +ENV OPENMC_REPO='https://github.com/openmc-dev/openmc' + +ARG build_dagmc +ARG build_libmesh + +ENV DAGMC_INSTALL_DIR=$HOME/DAGMC/ +ENV LIBMESH_INSTALL_DIR=$HOME/LIBMESH + # clone and install openmc RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ && git clone --shallow-submodules --recurse-submodules --single-branch -b ${openmc_branch} --depth=1 ${OPENMC_REPO} \ @@ -211,5 +218,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ FROM build AS release +ENV HOME=/root + # Download cross sections (NNDC and WMP) and ENDF data needed by test suite RUN ${HOME}/OpenMC/openmc/tools/ci/download-xs.sh From d8425d2669b71efbda64a3fbdaf645cf8468e291 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Jul 2022 06:48:51 -0500 Subject: [PATCH 15/83] Apply @shimwell suggestions from code review Co-authored-by: Jonathan Shimwell --- tests/unit_tests/test_lib.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index f1f84bff2..b57242983 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -829,6 +829,7 @@ def test_sample_external_source(run_in_tmpdir, mpi_intracomm): # Sample some particles and make sure they match specified source openmc.lib.init() particles = openmc.lib.sample_external_source(10, prn_seed=3) + assert len(particles) == 10 for p in particles: assert -5. < p.r[0] < 5. assert -5. < p.r[1] < 5. @@ -840,6 +841,7 @@ def test_sample_external_source(run_in_tmpdir, mpi_intracomm): # Using the same seed should produce the same particles other_particles = openmc.lib.sample_external_source(10, prn_seed=3) + assert len(other_particles) == 10 for p1, p2 in zip(particles, other_particles): assert p1.r == p2.r assert p1.u == p2.u From 476fdd1bf7169ac08d17cc529ae9fe1f8bbabcc8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Jul 2022 06:49:14 -0500 Subject: [PATCH 16/83] Improve description of prn_seed argument in openmc.lib functions --- openmc/lib/core.py | 3 ++- openmc/lib/math.py | 12 ++++++++---- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 3a437a295..479e414ae 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -431,7 +431,8 @@ def sample_external_source(n_samples=1, prn_seed=None): n_samples : int Number of samples prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- diff --git a/openmc/lib/math.py b/openmc/lib/math.py index d70c20a9f..0f61450bf 100644 --- a/openmc/lib/math.py +++ b/openmc/lib/math.py @@ -199,7 +199,8 @@ def rotate_angle(uvw0, mu, phi, prn_seed=None): phi : float Azimuthal angle; if None, one will be sampled uniformly prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- @@ -231,7 +232,8 @@ def maxwell_spectrum(T, prn_seed=None): T : float Spectrum parameter prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- @@ -256,7 +258,8 @@ def watt_spectrum(a, b, prn_seed=None): b : float Spectrum parameter b prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- @@ -281,7 +284,8 @@ def normal_variate(mean_value, std_dev, prn_seed=None): std_dev : float Standard deviation of the normal distribution prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- From dc7851dda262980943eb5d0c4461862248a6809e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Jul 2022 15:30:52 -0500 Subject: [PATCH 17/83] Apply @pshriwise suggestions from code review Co-authored-by: Patrick Shriwise --- openmc/weight_windows.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 257a17e9a..c5b7a00b6 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -291,10 +291,10 @@ class WeightWindows(IDManagerMixin): subelement.text = ' '.join(str(e) for e in self.energy_bounds) subelement = ET.SubElement(element, 'lower_ww_bounds') - subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds.ravel()) + subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds.ravel('F')) subelement = ET.SubElement(element, 'upper_ww_bounds') - subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds.ravel()) + subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds.ravel('F')) subelement = ET.SubElement(element, 'survival_ratio') subelement.text = str(self.survival_ratio) From fbb15496db978693cb5df9e963d69853b6f1ee9f Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 18 Jul 2022 21:33:37 +0100 Subject: [PATCH 18/83] boundary_type review suggestion from @pshriwise --- openmc/universe.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/openmc/universe.py b/openmc/universe.py index 989b50c73..959717b6a 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -746,7 +746,13 @@ class DAGMCUniverse(UniverseBase): import math bounding_box_center = (bounding_box[0] + bounding_box[1])/2 radius = math.dist(bounding_box[0], bounding_box[1]) - bounding_surface = openmc.Sphere(x0=bounding_box_center[0], y0=bounding_box_center[1], z0=bounding_box_center[2], r=radius) + bounding_surface = openmc.Sphere( + x0=bounding_box_center[0], + y0=bounding_box_center[1], + z0=bounding_box_center[2], + boundary_type=boundary_type, + r=radius, + ) return -bounding_surface From ae4556adda4112565f99e1c5b9153a143a322a41 Mon Sep 17 00:00:00 2001 From: Richard Morrison Date: Mon, 18 Jul 2022 23:51:14 +0100 Subject: [PATCH 19/83] flushing stdout --- src/output.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/output.cpp b/src/output.cpp index 0bd9693a8..56fea4db9 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -2,6 +2,7 @@ #include // for transform, max #include // for strlen +#include // for stdout #include // for time, localtime #include #include // for setw, setprecision, put_time @@ -94,6 +95,7 @@ void title() fmt::print(" OpenMP Threads | {}\n", omp_get_max_threads()); #endif fmt::print("\n"); + std::fflush(stdout); } //============================================================================== @@ -134,6 +136,7 @@ void header(const char* msg, int level) // Print header based on verbosity level. if (settings::verbosity >= level) fmt::print("\n{}\n\n", out); + std::fflush(stdout); } //============================================================================== @@ -380,6 +383,7 @@ void print_generation() fmt::print(" {:8.5f} +/-{:8.5f}", simulation::keff, simulation::keff_std); } fmt::print("\n"); + std::fflush(stdout); } //============================================================================== @@ -547,6 +551,7 @@ void print_results() gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n); } fmt::print("\n"); + std::fflush(stdout); } //============================================================================== From e72f97c4ed25e9b7c271bdd4f106afc89a4dd2ae Mon Sep 17 00:00:00 2001 From: Richard Morrison Date: Tue, 19 Jul 2022 00:04:28 +0100 Subject: [PATCH 20/83] Dockerfile ARGS to global scope --- Dockerfile | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/Dockerfile b/Dockerfile index c83460f35..158fed43d 100644 --- a/Dockerfile +++ b/Dockerfile @@ -15,14 +15,16 @@ # sudo docker run image_name:tag_name or ID with no tag sudo docker run ID number -FROM debian:bullseye-slim AS dependencies + +# global ARG as these ARGS are used in multiple stages +# By default one core is used to compile +ARG compile_cores=1 # By default this Dockerfile builds OpenMC without DAGMC and LIBMESH support ARG build_dagmc=off ARG build_libmesh=off -# By default one core is used to compile -ARG compile_cores=1 +FROM debian:bullseye-slim AS dependencies # Set default value of HOME to /root ENV HOME=/root @@ -56,7 +58,6 @@ ENV NJOY_REPO='https://github.com/njoy/NJOY2016' # Setup environment variables for Docker image ENV LD_LIBRARY_PATH=${DAGMC_INSTALL_DIR}/lib:$LD_LIBRARY_PATH \ - OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml \ OPENMC_ENDF_DATA=/root/endf-b-vii.1 \ DEBIAN_FRONTEND=noninteractive @@ -175,9 +176,6 @@ ENV HOME=/root ARG openmc_branch=master ENV OPENMC_REPO='https://github.com/openmc-dev/openmc' -ARG build_dagmc -ARG build_libmesh - ENV DAGMC_INSTALL_DIR=$HOME/DAGMC/ ENV LIBMESH_INSTALL_DIR=$HOME/LIBMESH @@ -219,6 +217,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ FROM build AS release ENV HOME=/root +ENV OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml # Download cross sections (NNDC and WMP) and ENDF data needed by test suite RUN ${HOME}/OpenMC/openmc/tools/ci/download-xs.sh From c8fb5eb2bdbacfc53d4280030ecae99cf25f0157 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 19 Jul 2022 11:43:28 +0100 Subject: [PATCH 21/83] updated comparison input --- .../dagmc/universes/inputs_true.dat | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/regression_tests/dagmc/universes/inputs_true.dat b/tests/regression_tests/dagmc/universes/inputs_true.dat index 4443f9a35..1a5ca7486 100644 --- a/tests/regression_tests/dagmc/universes/inputs_true.dat +++ b/tests/regression_tests/dagmc/universes/inputs_true.dat @@ -1,8 +1,8 @@ - + - + 24.0 24.0 2 2 -24.0 -24.0 @@ -10,12 +10,12 @@ 9 9 9 9 - - - - - - + + + + + + From b144f1a0fae08d0069d3bdbb1f24d71e6cff78a5 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 20 Jul 2022 13:22:01 +0100 Subject: [PATCH 22/83] added return types to openmc.Settings --- openmc/settings.py | 94 +++++++++++++++++++++++----------------------- 1 file changed, 47 insertions(+), 47 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 2d8121659..83a08fd77 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -298,191 +298,191 @@ class Settings: self._max_tracks = None @property - def run_mode(self): + def run_mode(self) -> str: return self._run_mode.value @property - def batches(self): + def batches(self) -> int: return self._batches @property - def generations_per_batch(self): + def generations_per_batch(self) -> int: return self._generations_per_batch @property - def inactive(self): + def inactive(self) -> int: return self._inactive @property - def max_lost_particles(self): + def max_lost_particles(self) -> int: return self._max_lost_particles @property - def rel_max_lost_particles(self): + def rel_max_lost_particles(self) -> int: return self._rel_max_lost_particles @property - def particles(self): + def particles(self) -> int: return self._particles @property - def keff_trigger(self): + def keff_trigger(self) -> dict: return self._keff_trigger @property - def energy_mode(self): + def energy_mode(self) -> str: return self._energy_mode @property - def max_order(self): + def max_order(self) -> int: return self._max_order @property - def source(self): + def source(self) -> Union[Source, typing.Iterable[Source]]: return self._source @property - def confidence_intervals(self): + def confidence_intervals(self) -> bool: return self._confidence_intervals @property - def electron_treatment(self): + def electron_treatment(self) -> str: return self._electron_treatment @property - def ptables(self): + def ptables(self) -> bool: return self._ptables @property - def photon_transport(self): + def photon_transport(self) -> bool: return self._photon_transport @property - def seed(self): + def seed(self) -> int: return self._seed @property - def survival_biasing(self): + def survival_biasing(self) -> bool: return self._survival_biasing @property - def entropy_mesh(self): + def entropy_mesh(self) -> RegularMesh: return self._entropy_mesh @property - def trigger_active(self): + def trigger_active(self) -> bool: return self._trigger_active @property - def trigger_max_batches(self): + def trigger_max_batches(self) -> int: return self._trigger_max_batches @property - def trigger_batch_interval(self): + def trigger_batch_interval(self) -> int: return self._trigger_batch_interval @property - def output(self): + def output(self) -> dict: return self._output @property - def sourcepoint(self): + def sourcepoint(self) -> dict: return self._sourcepoint @property - def statepoint(self): + def statepoint(self) -> dict: return self._statepoint @property - def surf_source_read(self): + def surf_source_read(self) -> dict: return self._surf_source_read @property - def surf_source_write(self): + def surf_source_write(self) -> dict: return self._surf_source_write @property - def no_reduce(self): + def no_reduce(self) -> bool: return self._no_reduce @property - def verbosity(self): + def verbosity(self) -> int: return self._verbosity @property - def tabular_legendre(self): + def tabular_legendre(self) -> dict: return self._tabular_legendre @property - def temperature(self): + def temperature(self) -> dict: return self._temperature @property - def trace(self): + def trace(self) -> typing.Iterable: return self._trace @property - def track(self): + def track(self) -> typing.Iterable[int]: return self._track @property - def cutoff(self): + def cutoff(self) -> dict: return self._cutoff @property - def ufs_mesh(self): + def ufs_mesh(self) -> RegularMesh: return self._ufs_mesh @property - def resonance_scattering(self): + def resonance_scattering(self) -> dict: return self._resonance_scattering @property - def volume_calculations(self): + def volume_calculations(self) -> Union[VolumeCalculation, typing.Iterable[VolumeCalculation]]: return self._volume_calculations @property - def create_fission_neutrons(self): + def create_fission_neutrons(self) -> bool: return self._create_fission_neutrons @property - def delayed_photon_scaling(self): + def delayed_photon_scaling(self) -> bool: return self._delayed_photon_scaling @property - def material_cell_offsets(self): + def material_cell_offsets(self) -> bool: return self._material_cell_offsets @property - def log_grid_bins(self): + def log_grid_bins(self) -> int: return self._log_grid_bins @property - def event_based(self): + def event_based(self) -> bool: return self._event_based @property - def max_particles_in_flight(self): + def max_particles_in_flight(self) -> int: return self._max_particles_in_flight @property - def write_initial_source(self): + def write_initial_source(self) -> bool: return self._write_initial_source @property - def weight_windows(self): + def weight_windows(self) -> Union[WeightWindows, typing.Iterable[WeightWindows]]: return self._weight_windows @property - def weight_windows_on(self): + def weight_windows_on(self) -> bool: return self._weight_windows_on @property - def max_splits(self): + def max_splits(self) -> int: return self._max_splits @property - def max_tracks(self): + def max_tracks(self) -> int: return self._max_tracks @run_mode.setter From 66206f8b147075e00af37ab0fe867e1a0ada055e Mon Sep 17 00:00:00 2001 From: Richard Morrison Date: Wed, 20 Jul 2022 20:31:12 +0100 Subject: [PATCH 23/83] ARG reference, use mpicxx --- Dockerfile | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Dockerfile b/Dockerfile index 158fed43d..70bccc041 100644 --- a/Dockerfile +++ b/Dockerfile @@ -176,6 +176,9 @@ ENV HOME=/root ARG openmc_branch=master ENV OPENMC_REPO='https://github.com/openmc-dev/openmc' +ARG build_dagmc +ARG build_libmesh + ENV DAGMC_INSTALL_DIR=$HOME/DAGMC/ ENV LIBMESH_INSTALL_DIR=$HOME/LIBMESH @@ -185,6 +188,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ && mkdir build && cd build ; \ if [ ${build_dagmc} = "on" ] && [ ${build_libmesh} = "on" ]; then \ cmake ../openmc \ + -DCMAKE_CXX_COMPILER=mpicxx \ -DOPENMC_USE_MPI=on \ -DHDF5_PREFER_PARALLEL=on \ -DOPENMC_USE_DAGMC=on \ @@ -193,6 +197,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ fi ; \ if [ ${build_dagmc} = "on" ] && [ ${build_libmesh} = "off" ]; then \ cmake ../openmc \ + -DCMAKE_CXX_COMPILER=mpicxx \ -DOPENMC_USE_MPI=on \ -DHDF5_PREFER_PARALLEL=on \ -DOPENMC_USE_DAGMC=ON \ @@ -200,6 +205,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ fi ; \ if [ ${build_dagmc} = "off" ] && [ ${build_libmesh} = "on" ]; then \ cmake ../openmc \ + -DCMAKE_CXX_COMPILER=mpicxx \ -DOPENMC_USE_MPI=on \ -DHDF5_PREFER_PARALLEL=on \ -DOPENMC_USE_LIBMESH=on \ @@ -207,6 +213,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ fi ; \ if [ ${build_dagmc} = "off" ] && [ ${build_libmesh} = "off" ]; then \ cmake ../openmc \ + -DCMAKE_CXX_COMPILER=mpicxx \ -DOPENMC_USE_MPI=on \ -DHDF5_PREFER_PARALLEL=on ; \ fi ; \ From 56b593d0b5632592260807f84288564968046f66 Mon Sep 17 00:00:00 2001 From: Richard Morrison Date: Wed, 20 Jul 2022 20:40:17 +0100 Subject: [PATCH 24/83] compile_cores ARG included in build stages --- Dockerfile | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Dockerfile b/Dockerfile index 70bccc041..d408b9de5 100644 --- a/Dockerfile +++ b/Dockerfile @@ -26,6 +26,10 @@ ARG build_libmesh=off FROM debian:bullseye-slim AS dependencies +ARG compile_cores +ARG build_dagmc +ARG build_libmesh + # Set default value of HOME to /root ENV HOME=/root @@ -176,6 +180,7 @@ ENV HOME=/root ARG openmc_branch=master ENV OPENMC_REPO='https://github.com/openmc-dev/openmc' +ARG compile_cores ARG build_dagmc ARG build_libmesh From 4df20f2d6c796cb76af8d8b9eaef73ff272ea60b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Jul 2022 22:32:03 -0500 Subject: [PATCH 25/83] Account for neutron wave number in SLBW competitve rho. Closes #907 --- openmc/data/reconstruct.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/data/reconstruct.pyx b/openmc/data/reconstruct.pyx index f63a155b1..cd0bbc38b 100644 --- a/openmc/data/reconstruct.pyx +++ b/openmc/data/reconstruct.pyx @@ -299,8 +299,8 @@ def reconstruct_slbw(slbw, double E): # Determine shift and penetration at modified energy if slbw._competitive[i]: Ex = E + slbw.q_value[l]*(A + 1)/A - rhoc = slbw.channel_radius[l](Ex) - rhochat = slbw.scattering_radius[l](Ex) + rhoc = k*slbw.channel_radius[l](Ex) + rhochat = k*slbw.scattering_radius[l](Ex) P_c, S_c = penetration_shift(l, rhoc) if Ex < 0: P_c = 0 From a27d0e079fe1d69f227e21d55b2d72e8e3e176f5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Jul 2022 22:43:20 -0500 Subject: [PATCH 26/83] Allow inplace argument on Region.translate --- openmc/region.py | 13 ++++++++----- openmc/surface.py | 19 +++++++++---------- tests/unit_tests/test_region.py | 14 ++++++++++++++ 3 files changed, 31 insertions(+), 15 deletions(-) diff --git a/openmc/region.py b/openmc/region.py index 56926187d..4e74a08ad 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -259,13 +259,16 @@ class Region(ABC): clone[:] = [n.clone(memo) for n in self] return clone - def translate(self, vector, memo=None): + def translate(self, vector, inplace=False, memo=None): """Translate region in given direction Parameters ---------- vector : iterable of float Direction in which region should be translated + inplace : bool + Whether or not to return a region based on new surfaces or one based + on the original surfaces that have been modified. memo : dict or None Dictionary used for memoization. This parameter is used internally and should not be specified by the user. @@ -279,7 +282,7 @@ class Region(ABC): if memo is None: memo = {} - return type(self)(n.translate(vector, memo) for n in self) + return type(self)(n.translate(vector, inplace, memo) for n in self) def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, memo=None): @@ -308,7 +311,7 @@ class Region(ABC): :math:`\psi` about z. This corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. - inplace : boolean + inplace : bool Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. memo : dict or None @@ -622,10 +625,10 @@ class Complement(Region): clone.node = self.node.clone(memo) return clone - def translate(self, vector, memo=None): + def translate(self, vector, inplace=False, memo=None): if memo is None: memo = {} - return type(self)(self.node.translate(vector, memo)) + return type(self)(self.node.translate(vector, inplace, memo)) def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, memo=None): diff --git a/openmc/surface.py b/openmc/surface.py index b3bcebb11..2996897a4 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -336,9 +336,9 @@ class Surface(IDManagerMixin, ABC): ---------- vector : iterable of float Direction in which surface should be translated - inplace : boolean + inplace : bool Whether or not to return a new instance of this Surface or to - modify the coefficients of this Surface. Defaults to False + modify the coefficients of this Surface. Returns ------- @@ -374,7 +374,7 @@ class Surface(IDManagerMixin, ABC): :math:`\psi` about z. This corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. - inplace : boolean + inplace : bool Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. @@ -568,9 +568,9 @@ class PlaneMixin: ---------- vector : iterable of float Direction in which surface should be translated - inplace : boolean + inplace : bool Whether or not to return a new instance of a Plane or to modify the - coefficients of this plane. Defaults to False + coefficients of this plane. Returns ------- @@ -1012,9 +1012,8 @@ class QuadricMixin: ---------- vector : iterable of float Direction in which surface should be translated - inplace : boolean + inplace : bool Whether to return a clone of the Surface or the Surface itself. - Defaults to False Returns ------- @@ -2563,7 +2562,7 @@ class Halfspace(Region): clone.surface = self.surface.clone(memo) return clone - def translate(self, vector, memo=None): + def translate(self, vector, inplace=False, memo=None): """Translate half-space in given direction Parameters @@ -2585,7 +2584,7 @@ class Halfspace(Region): # If translated surface not in memo, add it key = (self.surface, tuple(vector)) if key not in memo: - memo[key] = self.surface.translate(vector) + memo[key] = self.surface.translate(vector, inplace) # Return translated half-space return type(self)(memo[key], self.side) @@ -2617,7 +2616,7 @@ class Halfspace(Region): :math:`\psi` about z. This corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. - inplace : boolean + inplace : bool Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. memo : dict or None diff --git a/tests/unit_tests/test_region.py b/tests/unit_tests/test_region.py index 8537e3b80..086ce6401 100644 --- a/tests/unit_tests/test_region.py +++ b/tests/unit_tests/test_region.py @@ -208,3 +208,17 @@ def test_from_expression(reset): # Opening parenthesis immediately after halfspace r = openmc.Region.from_expression('1(2|-3)', surfs) assert str(r) == '(1 (2 | -3))' + + +def test_translate_inplace(): + sph = openmc.Sphere() + x = openmc.XPlane() + region = -sph & +x + + # Translating a region should produce new surfaces + region2 = region.translate((0.5, -6.7, 3.9), inplace=False) + assert str(region) != str(region2) + + # Translating a region in-place should *not* produce new surfaces + region3 = region.translate((0.5, -6.7, 3.9), inplace=True) + assert str(region) == str(region3) From 8bf727022f18d87878d3837b663d33bae0aa4b7b Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 21 Jul 2022 09:12:08 +0100 Subject: [PATCH 27/83] return types corrections by @paulromano Co-authored-by: Paul Romano --- openmc/settings.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 83a08fd77..5b5c265f6 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -318,7 +318,7 @@ class Settings: return self._max_lost_particles @property - def rel_max_lost_particles(self) -> int: + def rel_max_lost_particles(self) -> float: return self._rel_max_lost_particles @property @@ -338,7 +338,7 @@ class Settings: return self._max_order @property - def source(self) -> Union[Source, typing.Iterable[Source]]: + def source(self) -> typing.List[Source]: return self._source @property @@ -422,7 +422,7 @@ class Settings: return self._trace @property - def track(self) -> typing.Iterable[int]: + def track(self) -> typing.Iterable[typing.Iterable[int]]: return self._track @property @@ -438,7 +438,7 @@ class Settings: return self._resonance_scattering @property - def volume_calculations(self) -> Union[VolumeCalculation, typing.Iterable[VolumeCalculation]]: + def volume_calculations(self) -> typing.List[VolumeCalculation]: return self._volume_calculations @property @@ -470,7 +470,7 @@ class Settings: return self._write_initial_source @property - def weight_windows(self) -> Union[WeightWindows, typing.Iterable[WeightWindows]]: + def weight_windows(self) -> typing.List[WeightWindows]: return self._weight_windows @property From 035b0df3d85e22693d8206cc662efdee1b7557e7 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 21 Jul 2022 09:14:15 +0100 Subject: [PATCH 28/83] corrected rel_max_lost_particles type --- openmc/settings.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 5b5c265f6..47c61e6ae 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -75,7 +75,7 @@ class Settings: Maximum number of lost particles .. versionadded:: 0.12 - rel_max_lost_particles : int + rel_max_lost_particles : float Maximum number of lost particles, relative to the total number of particles .. versionadded:: 0.12 @@ -517,7 +517,7 @@ class Settings: self._max_lost_particles = max_lost_particles @rel_max_lost_particles.setter - def rel_max_lost_particles(self, rel_max_lost_particles: int): + def rel_max_lost_particles(self, rel_max_lost_particles: float): cv.check_type('rel_max_lost_particles', rel_max_lost_particles, Real) cv.check_greater_than('rel_max_lost_particles', rel_max_lost_particles, 0) cv.check_less_than('rel_max_lost_particles', rel_max_lost_particles, 1) From ada2493a6b901aba7e404da2ca90f14ea9c166f6 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 21 Jul 2022 12:05:31 +0100 Subject: [PATCH 29/83] changed auto_geom_ids approach --- openmc/universe.py | 10 +--------- src/dagmc.cpp | 4 +++- 2 files changed, 4 insertions(+), 10 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 959717b6a..5bcbfde46 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -767,26 +767,18 @@ class DAGMCUniverse(UniverseBase): return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z - def bounded_universe(self, auto_geom_ids=True, **kwargs): + def bounded_universe(self, **kwargs): """Returns an openmc.Universe filled with this DAGMCUniverse and bounded with a cell. Defaults to a box cell with a vacuum surface however this can be changed using the kwargs which are passed directly to DAGMCUniverse.bounding_region(). - Parameters - ---------- - auto_geom_ids : bool - Set IDs automatically on initialization (True) or report overlaps - in ID space between CSG and DAGMC (False). Defaults to True to avoid - overlapping ID numbers between the CSG and DAGMC geometry ID numbers - Returns ------- openmc.Universe Universe instance """ - self.auto_geom_ids = auto_geom_ids bounding_cell = openmc.Cell(fill=self, region=self.bounding_region(**kwargs)) return openmc.Universe(cells=[bounding_cell]) diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 8ce5aabda..c9a6e6ad2 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -183,7 +183,9 @@ void DAGUniverse::init_geometry() } else { warning(fmt::format("DAGMC Cell IDs: {}", dagmc_ids_for_dim(3))); fatal_error(fmt::format("Cell ID {} exists in both DAGMC Universe {} " - "and the CSG geometry.", + "and the CSG geometry. Setting auto_geom_ids + to True when initiating the DAGMC Universe may + resolve this issue", c->id_, this->id_)); } From 0256fe90d47a642ace88b5c6168fde50f832920b Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 21 Jul 2022 16:11:57 +0100 Subject: [PATCH 30/83] corrected end of line --- src/dagmc.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dagmc.cpp b/src/dagmc.cpp index c9a6e6ad2..dc3bc2b2c 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -183,9 +183,9 @@ void DAGUniverse::init_geometry() } else { warning(fmt::format("DAGMC Cell IDs: {}", dagmc_ids_for_dim(3))); fatal_error(fmt::format("Cell ID {} exists in both DAGMC Universe {} " - "and the CSG geometry. Setting auto_geom_ids - to True when initiating the DAGMC Universe may - resolve this issue", + "and the CSG geometry. Setting auto_geom_ids " + "to True when initiating the DAGMC Universe may " + "resolve this issue", c->id_, this->id_)); } From d024a3e00842a6535d8d0711a4131e4bee1ee841 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Oct 2021 07:26:43 -0500 Subject: [PATCH 31/83] Implementation of C++ Sum1D function --- include/openmc/endf.h | 20 ++++++++++++++++++++ src/endf.cpp | 26 ++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 7beb8e452..e580874b4 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -126,6 +126,26 @@ private: debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1] }; +//============================================================================== +//! Sum of multiple 1D functions +//============================================================================== + +class Sum1D : public Function1D { +public: + // Constructors + explicit Sum1D(hid_t group); + + //! Evaluate each function and sum results + //! \param[in] x independent variable + //! \return Function evaluated at x + double operator()(double E) const override; + + const unique_ptr& functions(int i) const { return functions_[i]; } + +private: + vector> functions_; //!< individual functions +}; + //! Read 1D function from HDF5 dataset //! \param[in] group HDF5 group containing dataset //! \param[in] name Name of dataset diff --git a/src/endf.cpp b/src/endf.cpp index b42c8641d..60db3efa4 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -271,4 +271,30 @@ double IncoherentElasticXS::operator()(double E) const return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W)); } +//============================================================================== +// Sum1D implementation +//============================================================================== + +Sum1D::Sum1D(hid_t group) +{ + // Get number of functions + int n; + read_attribute(group, "n", n); + + // Get each function + for (int i = 0; i < n; ++i) { + auto dset_name = fmt::format("func_{}", i + 1); + functions_.push_back(read_function(group, dset_name.c_str())); + } +} + +double Sum1D::operator()(double x) const +{ + double result = 0.0; + for (auto& func : functions_) { + result += (*func)(x); + } + return result; +} + } // namespace openmc From 8bb2002d8f8bebd58f11d84d88fc88df20552e8d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Oct 2021 07:27:56 -0500 Subject: [PATCH 32/83] Add open_object, close_object functions in HDF5 interface --- include/openmc/hdf5_interface.h | 2 ++ src/endf.cpp | 16 ++++++++-------- src/hdf5_interface.cpp | 12 ++++++++++++ 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index e9d9b4f96..0092c08f8 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -61,6 +61,8 @@ void ensure_exists(hid_t obj_id, const char* name, bool attribute = false); vector group_names(hid_t group_id); vector object_shape(hid_t obj_id); std::string object_name(hid_t obj_id); +hid_t open_object(hid_t group_id, const std::string& name); +void close_object(hid_t obj_id); //============================================================================== // Fortran compatibility functions diff --git a/src/endf.cpp b/src/endf.cpp index 60db3efa4..3ccfdbd21 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -90,23 +90,23 @@ bool is_inelastic_scatter(int mt) unique_ptr read_function(hid_t group, const char* name) { - hid_t dset = open_dataset(group, name); + hid_t obj_id = open_object(group, name); std::string func_type; - read_attribute(dset, "type", func_type); + read_attribute(obj_id, "type", func_type); unique_ptr func; if (func_type == "Tabulated1D") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "Polynomial") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "CoherentElastic") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "IncoherentElastic") { - func = make_unique(dset); + func = make_unique(obj_id); } else { throw std::runtime_error {"Unknown function type " + func_type + - " for dataset " + object_name(dset)}; + " for dataset " + object_name(obj_id)}; } - close_dataset(dset); + close_object(obj_id); return func; } diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index d53ad3879..27b750b93 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -118,6 +118,12 @@ void close_group(hid_t group_id) fatal_error("Failed to close group"); } +void close_object(hid_t obj_id) +{ + if (H5Oclose(obj_id) < 0) + fatal_error("Failed to close object"); +} + int dataset_ndims(hid_t dset) { hid_t dspace = H5Dget_space(dset); @@ -394,6 +400,12 @@ hid_t open_group(hid_t group_id, const char* name) return H5Gopen(group_id, name, H5P_DEFAULT); } +hid_t open_object(hid_t group_id, const std::string& name) +{ + ensure_exists(group_id, name.c_str()); + return H5Oopen(group_id, name.c_str(), H5P_DEFAULT); +} + void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer) { hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT); From 8cc815d4f2ebf56460f0f0f6f9195bf46afc968a Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 21 Jul 2022 22:24:17 -0500 Subject: [PATCH 33/83] Add numpy arrays to a list so they persist to the time of the VTK write --- openmc/mesh.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 1fe5429d0..d08213f77 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -248,8 +248,10 @@ class StructuredMesh(MeshBase): # create VTK arrays for each of # the data sets + datasets_out = [] for label, dataset in datasets.items(): dataset = np.asarray(dataset).flatten() + datasets_out.append(dataset) if volume_normalization: dataset /= self.volumes.flatten() @@ -692,7 +694,7 @@ class RegularMesh(StructuredMesh): root_cell.fill = lattice return root_cell, cells - + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh @@ -1605,7 +1607,7 @@ class UnstructuredMesh(MeshBase): Raises ------ RuntimeError - when the size of a dataset doesn't match the number of cells + when the size of a dataset doesn't match the number of cells """ import vtk From ff5490e3ba432209a95cb75ce07d21b167f90efd Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 22 Jul 2022 08:01:28 +0100 Subject: [PATCH 34/83] User able to set bounding cell id --- openmc/universe.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 5bcbfde46..fe41674d5 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -767,19 +767,25 @@ class DAGMCUniverse(UniverseBase): return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z - def bounded_universe(self, **kwargs): + def bounded_universe(self, bounding_cell_id=10000, **kwargs): """Returns an openmc.Universe filled with this DAGMCUniverse and bounded with a cell. Defaults to a box cell with a vacuum surface however this can be changed using the kwargs which are passed directly to DAGMCUniverse.bounding_region(). + Parameters + ---------- + bounding_cell_id : int + The cell ID number to use for the bounding cell, defaults to 1000 to reduce + the chance of overlapping ID numbers with the DAGMC geometry. + Returns ------- openmc.Universe Universe instance """ - bounding_cell = openmc.Cell(fill=self, region=self.bounding_region(**kwargs)) + bounding_cell = openmc.Cell(fill=self, id=bounding_cell_id, region=self.bounding_region(**kwargs)) return openmc.Universe(cells=[bounding_cell]) @classmethod From f0cda31e73311a2673bc7fb0dac698f396cf2933 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 22 Jul 2022 10:45:15 +0100 Subject: [PATCH 35/83] corrected arg name --- openmc/universe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/universe.py b/openmc/universe.py index fe41674d5..78f75bbfc 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -785,7 +785,7 @@ class DAGMCUniverse(UniverseBase): Universe instance """ - bounding_cell = openmc.Cell(fill=self, id=bounding_cell_id, region=self.bounding_region(**kwargs)) + bounding_cell = openmc.Cell(fill=self, cell_id =bounding_cell_id, region=self.bounding_region(**kwargs)) return openmc.Universe(cells=[bounding_cell]) @classmethod From e453901c84e57eca983a52124e7dc903b06dfb00 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Oct 2021 07:28:35 -0500 Subject: [PATCH 36/83] Implement to/from_hdf5 methods for Python Sum class --- openmc/data/function.py | 39 +++++++++++++++++++++++++++++++++++++++ src/endf.cpp | 2 ++ 2 files changed, 41 insertions(+) diff --git a/openmc/data/function.py b/openmc/data/function.py index b5aa2117d..7a73987df 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -578,6 +578,45 @@ class Sum(EqualityMixin): cv.check_type('functions', functions, Iterable, Callable) self._functions = functions + def to_hdf5(self, group, name='xy'): + """Write sum of functions to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + name : str + Name of the dataset to create + + """ + sum_group = group.create_group(name) + sum_group.attrs['type'] = np.string_(type(self).__name__) + sum_group.attrs['n'] = len(self.functions) + for i, f in enumerate(self.functions): + f.to_hdf5(sum_group, f'func_{i+1}') + + @classmethod + def from_hdf5(cls, group): + """Generate sum of functions from an HDF5 group + + Parameters + ---------- + group : h5py.Group + Group to read from + + Returns + ------- + openmc.data.Sum + Functions read from the group + + """ + n = group.attrs['n'] + functions = [ + Function1D.from_hdf5(group[f'func_{i+1}']) + for i in range(n) + ] + return cls(functions) + class Regions1D(EqualityMixin): r"""Piecewise composition of multiple functions. diff --git a/src/endf.cpp b/src/endf.cpp index 3ccfdbd21..c0c1d2e7e 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -102,6 +102,8 @@ unique_ptr read_function(hid_t group, const char* name) func = make_unique(obj_id); } else if (func_type == "IncoherentElastic") { func = make_unique(obj_id); + } else if (func_type == "Sum") { + func = make_unique(obj_id); } else { throw std::runtime_error {"Unknown function type " + func_type + " for dataset " + object_name(obj_id)}; From 70c16d4b04a60606742f6b1c6825bfd04adcad6d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Oct 2021 07:30:56 -0500 Subject: [PATCH 37/83] Implement mixed coherent/incoheret thermal elastic --- include/openmc/secondary_thermal.h | 28 +++++++++ openmc/data/angle_energy.py | 2 + openmc/data/thermal.py | 91 ++++++++++++++++++++--------- openmc/data/thermal_angle_energy.py | 57 ++++++++++++++++++ src/secondary_thermal.cpp | 35 +++++++++++ src/thermal.cpp | 25 +++++--- 6 files changed, 203 insertions(+), 35 deletions(-) diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 84ed68f51..81de5c451 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -150,6 +150,34 @@ private: //!< each incident energy }; +//============================================================================== +//! Mixed coherent/incoherent elastic angle-energy distribution +//============================================================================== + +class MixedElasticAE : public AngleEnergy { +public: + //! Construct from HDF5 file + // + //! \param[in] group HDF5 group + explicit MixedElasticAE( + hid_t group, const CoherentElasticXS& coh_xs, const Tabulated1D& incoh_xs); + + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction + //! \param[inout] seed Pseudorandom number seed pointer + void sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const override; + +private: + CoherentElasticAE coherent_dist_; //!< Coherent distribution + unique_ptr incoherent_dist_; //!< Incoherent distribution + + const CoherentElasticXS& coherent_xs_; //!< Ref. to coherent XS + const Tabulated1D& incoherent_xs_; //!< Ref. to incoherent XS +}; + } // namespace openmc #endif // OPENMC_SECONDARY_THERMAL_H diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index 6009dc748..b8fde5478 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -44,6 +44,8 @@ class AngleEnergy(EqualityMixin, ABC): return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group) elif dist_type == 'incoherent_inelastic': return openmc.data.IncoherentInelasticAE.from_hdf5(group) + elif dist_type == 'mixed_elastic': + return openmc.data.MixedElasticAE.from_hdf5(group) @staticmethod def from_ace(ace, location_dist, location_start, rx=None): diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 4e0a30b64..6cf5a4e9e 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -19,12 +19,12 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, isotopes from .ace import Table, get_table, Library from .angle_energy import AngleEnergy -from .function import Tabulated1D, Function1D +from .function import Tabulated1D, Function1D, Sum from .njoy import make_ace_thermal from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE, IncoherentElasticAEDiscrete, IncoherentInelasticAEDiscrete, - IncoherentInelasticAE) + IncoherentInelasticAE, MixedElasticAE) _THERMAL_NAMES = { @@ -694,29 +694,53 @@ class ThermalScattering(EqualityMixin): # Incoherent/coherent elastic scattering cross section idx = ace.jxs[4] - n_mu = ace.nxs[6] + 1 if idx != 0: - n_energy = int(ace.xss[idx]) - energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV - P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] - - if ace.nxs[5] == 4: + if ace.nxs[5] in (4, 5): # Coherent elastic - xs = CoherentElastic(energy, P*EV_PER_MEV) - distribution = CoherentElasticAE(xs) + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV + P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] + coherent_xs = CoherentElastic(energy, P*EV_PER_MEV) + coherent_dist = CoherentElasticAE(xs) # Coherent elastic shouldn't have angular distributions listed + n_mu = ace.nxs[6] + 1 assert n_mu == 0 - else: - # Incoherent elastic - xs = Tabulated1D(energy, P) + + if ace.nxs[5] in (3, 5): + # Incoherent elastic scattering -- first determine if both + # incoherent and coherent are present (mixed) + mixed = (ace.nxs[5] == 5) + + # Get cross section values + idx = ace.jxs[7] if mixed else ace.jxs[4] + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV + values = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] + + incoherent_xs = Tabulated1D(energy, values) # Angular distribution + n_mu = (ace.nxs[8] if mixed else ace.nxs[6]) + 1 assert n_mu > 0 - idx = ace.jxs[6] + idx = ace.jxs[9] if mixed else ace.jxs[6] mu_out = ace.xss[idx:idx + n_energy * n_mu] mu_out.shape = (n_energy, n_mu) - distribution = IncoherentElasticAEDiscrete(mu_out) + incoherent_dist = IncoherentElasticAEDiscrete(mu_out) + + if ace.nxs[5] == 3: + xs = incoherent_xs + dist = incoherent_dist + elif ace.nxs[5] == 4: + xs = coherent_xs + dist = coherent_dist + else: + # Create mixed cross section -- note that coherent must come + # first due to assumption on C++ side + xs = Sum([coherent_xs, incoherent_xs]) + + # Create mixed distribution + distribution = MixedElasticAE(coherent_dist, incoherent_dist) table.elastic = ThermalScatteringReaction({T: xs}, {T: distribution}) @@ -802,7 +826,7 @@ class ThermalScattering(EqualityMixin): # Replace ACE data with ENDF data rx, rx_endf = data.elastic, data_endf.elastic for t in temperatures: - if isinstance(rx_endf.xs[t], IncoherentElastic): + if isinstance(rx_endf.xs[t], (IncoherentElastic, Sum)): rx.xs[t] = rx_endf.xs[t] rx.distribution[t] = rx_endf.distribution[t] @@ -832,20 +856,14 @@ class ThermalScattering(EqualityMixin): # Read coherent/incoherent elastic data elastic = None if (7, 2) in ev.section: - xs = {} - distribution = {} - - file_obj = StringIO(ev.section[7, 2]) - lhtr = endf.get_head_record(file_obj)[2] - if lhtr == 1: - # coherent elastic - + # Define helper functions to avoid duplication + def get_coherent_elastic(file_obj): # Get structure factor at first temperature params, S = endf.get_tab1_record(file_obj) strT = _temperature_str(params[0]) n_temps = params[2] bragg_edges = S.x - xs[strT] = CoherentElastic(bragg_edges, S.y) + xs = {strT: CoherentElastic(bragg_edges, S.y)} distribution = {strT: CoherentElasticAE(xs[strT])} # Get structure factor for subsequent temperatures @@ -854,15 +872,34 @@ class ThermalScattering(EqualityMixin): strT = _temperature_str(params[0]) xs[strT] = CoherentElastic(bragg_edges, S) distribution[strT] = CoherentElasticAE(xs[strT]) + return xs, distribution - elif lhtr == 2: - # incoherent elastic + def get_incoherent_elastic(file_obj): params, W = endf.get_tab1_record(file_obj) bound_xs = params[0] + xs = {} + distribution = {} for T, debye_waller in zip(W.x, W.y): strT = _temperature_str(T) xs[strT] = IncoherentElastic(bound_xs, debye_waller) distribution[strT] = IncoherentElasticAE(debye_waller) + return xs, distribution + + file_obj = StringIO(ev.section[7, 2]) + lhtr = endf.get_head_record(file_obj)[2] + if lhtr == 1: + # coherent elastic + xs, distribution = get_coherent_elastic(file_obj) + elif lhtr == 2: + # incoherent elastic + xs, distribution = get_incoherent_elastic(file_obj) + elif lhtr == 3: + # mixed coherent / incoherent elastic + xs_c, dist_c = get_coherent_elastic(file_obj) + xs_i, dist_i = get_incoherent_elastic(file_obj) + assert sorted(xs_c) == sorted(xs_i) + xs = {T: Sum([xs_c[T], xs_i[T]]) for T in xs_c} + distribution = {T: MixedElasticAE(dist_c[T], dist_i[T]) for T in dist_c} elastic = ThermalScatteringReaction(xs, distribution) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index 05393e7bb..6ef874600 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -210,3 +210,60 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): class IncoherentInelasticAE(CorrelatedAngleEnergy): _name = 'incoherent_inelastic' + + +class MixedElasticAE(AngleEnergy): + """Secondary distribution for mixed coherent/incoherent thermal elastic + + Parameters + ---------- + coherent : AngleEnergy + Secondary distribution for coherent elastic scattering + incoherent : AngleEnergy + Secondary distribution for incoherent elastic scattering + + Attributes + ---------- + coherent : AngleEnergy + Secondary distribution for coherent elastic scattering + incoherent : AngleEnergy + Secondary distribution for incoherent elastic scattering + + """ + def __init__(self, coherent, incoherent): + self.coherent = coherent + self.incoherent = incoherent + + def to_hdf5(self, group): + """Write mixed elastic distribution to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + + """ + group.attrs['type'] = np.string_('mixed_elastic') + coherent_group = group.create_group('coherent') + self.coherent.to_hdf5(coherent_group) + incoherent_group = group.create_group('incoherent') + self.incoherent.to_hdf5(incoherent_group) + + @classmethod + def from_hdf5(cls, group): + """Generate mixed thermal elastic distribution from HDF5 data + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.MixedElasticAE + Mixed thermal elastic distribution + + """ + coherent = AngleEnergy.from_hdf5(group['coherent']) + incoherent = AngleEnergy.from_hdf5(group['incoherent']) + return cls(coherent, incoherent) diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 3677f6ed3..e5e8274da 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -332,4 +332,39 @@ void IncoherentInelasticAE::sample( mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5); } +//============================================================================== +// MixedElasticAE implementation +//============================================================================== + +MixedElasticAE::MixedElasticAE( + hid_t group, const CoherentElasticXS& coh_xs, const Tabulated1D& incoh_xs) + : coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs) +{ + // Read incoherent elastic distribution + hid_t incoherent_group = open_group(group, "incoherent"); + std::string temp; + read_attribute(incoherent_group, "type", temp); + if (temp == "incoherent_elastic") { + incoherent_dist_ = make_unique(incoherent_group); + } else if (temp == "incoherent_elastic_discrete") { + incoherent_dist_ = + make_unique(incoherent_group, incoh_xs.x()); + } + close_group(incoherent_group); +} + +void MixedElasticAE::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + // Evaluate coherent and incoherent elastic cross sections + double xs_coh = coherent_xs_(E_in); + double xs_incoh = incoherent_xs_(E_in); + + if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) { + coherent_dist_.sample(E_in, E_out, mu, seed); + } else { + incoherent_dist_->sample(E_in, E_out, mu, seed); + } +} + } // namespace openmc diff --git a/src/thermal.cpp b/src/thermal.cpp index 1101d5485..303786915 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -210,14 +210,23 @@ ThermalData::ThermalData(hid_t group) if (temp == "coherent_elastic") { auto xs = dynamic_cast(elastic_.xs.get()); elastic_.distribution = make_unique(*xs); - } else { - if (temp == "incoherent_elastic") { - elastic_.distribution = make_unique(dgroup); - } else if (temp == "incoherent_elastic_discrete") { - auto xs = dynamic_cast(elastic_.xs.get()); - elastic_.distribution = - make_unique(dgroup, xs->x()); - } + } else if (temp == "incoherent_elastic") { + elastic_.distribution = make_unique(dgroup); + } else if (temp == "incoherent_elastic_discrete") { + auto xs = dynamic_cast(elastic_.xs.get()); + elastic_.distribution = + make_unique(dgroup, xs->x()); + } else if (temp == "mixed_elastic") { + // Get coherent/incoherent cross sections + auto mixed_xs = dynamic_cast(elastic_.xs.get()); + const auto& coh_xs = + dynamic_cast(mixed_xs->functions(0).get()); + const auto& incoh_xs = + dynamic_cast(mixed_xs->functions(1).get()); + + // Create mixed elastic distribution + elastic_.distribution = + make_unique(dgroup, *coh_xs, *incoh_xs); } close_group(elastic_group); From edf023adb9d48b5b8078ad3ca4886b0176457902 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Oct 2021 17:05:49 -0500 Subject: [PATCH 38/83] Fix HDF5 export of mixed elastic distribution --- openmc/data/thermal_angle_energy.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index 6ef874600..d7dfe92f9 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -33,17 +33,22 @@ class CoherentElasticAE(AngleEnergy): def __init__(self, coherent_xs): self.coherent_xs = coherent_xs - def to_hdf5(self, group): + def to_hdf5(self, group, xs_group=None): """Write coherent elastic distribution to an HDF5 group Parameters ---------- group : h5py.Group HDF5 group to write to + xs_group : h5py.Group, optional + Group containing 'xs' dataset """ group.attrs['type'] = np.string_('coherent_elastic') - group['coherent_xs'] = group.parent['xs'] + if xs_group is not None: + group['coherent_xs'] = xs_group['xs'] + else: + group['coherent_xs'] = group.parent['xs'] class IncoherentElasticAE(AngleEnergy): @@ -245,7 +250,7 @@ class MixedElasticAE(AngleEnergy): """ group.attrs['type'] = np.string_('mixed_elastic') coherent_group = group.create_group('coherent') - self.coherent.to_hdf5(coherent_group) + self.coherent.to_hdf5(coherent_group, group.parent) incoherent_group = group.create_group('incoherent') self.incoherent.to_hdf5(incoherent_group) From 4e611f59b9ac60839fa0224a91be929572f0dbc3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Oct 2021 06:53:20 -0500 Subject: [PATCH 39/83] Make incoherent XS polymorphic in MixedElasticAE --- include/openmc/secondary_thermal.h | 4 ++-- src/secondary_thermal.cpp | 5 +++-- src/thermal.cpp | 3 +-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 81de5c451..5b18902af 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -160,7 +160,7 @@ public: // //! \param[in] group HDF5 group explicit MixedElasticAE( - hid_t group, const CoherentElasticXS& coh_xs, const Tabulated1D& incoh_xs); + hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs); //! Sample distribution for an angle and energy //! \param[in] E_in Incoming energy in [eV] @@ -175,7 +175,7 @@ private: unique_ptr incoherent_dist_; //!< Incoherent distribution const CoherentElasticXS& coherent_xs_; //!< Ref. to coherent XS - const Tabulated1D& incoherent_xs_; //!< Ref. to incoherent XS + const Function1D& incoherent_xs_; //!< Polymorphic ref. to incoherent XS }; } // namespace openmc diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index e5e8274da..0b8e1ab42 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -337,7 +337,7 @@ void IncoherentInelasticAE::sample( //============================================================================== MixedElasticAE::MixedElasticAE( - hid_t group, const CoherentElasticXS& coh_xs, const Tabulated1D& incoh_xs) + hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs) : coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs) { // Read incoherent elastic distribution @@ -347,8 +347,9 @@ MixedElasticAE::MixedElasticAE( if (temp == "incoherent_elastic") { incoherent_dist_ = make_unique(incoherent_group); } else if (temp == "incoherent_elastic_discrete") { + auto xs = dynamic_cast(&incoh_xs); incoherent_dist_ = - make_unique(incoherent_group, incoh_xs.x()); + make_unique(incoherent_group, xs->x()); } close_group(incoherent_group); } diff --git a/src/thermal.cpp b/src/thermal.cpp index 303786915..1a9592d0a 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -221,8 +221,7 @@ ThermalData::ThermalData(hid_t group) auto mixed_xs = dynamic_cast(elastic_.xs.get()); const auto& coh_xs = dynamic_cast(mixed_xs->functions(0).get()); - const auto& incoh_xs = - dynamic_cast(mixed_xs->functions(1).get()); + const auto& incoh_xs = mixed_xs->functions(1).get(); // Create mixed elastic distribution elastic_.distribution = From e8294d52a0162b484a0961195dc330ef58b5f051 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 21 Jul 2022 14:37:34 -0500 Subject: [PATCH 40/83] Update documentation for mixed thermal elastic --- docs/source/io_formats/nuclear_data.rst | 23 +++++++++++++++++++++++ docs/source/pythonapi/data.rst | 1 + openmc/data/function.py | 4 ++++ openmc/data/thermal_angle_energy.py | 5 +++++ 4 files changed, 33 insertions(+) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index ef61604cb..8174108e1 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -339,6 +339,16 @@ Incoherent elastic scattering [eV\ :math:`^{-1}`]. :Attributes: - **type** (*char[]*) -- 'IncoherentElastic' +Sum of functions +---------------- + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "Sum" + - **n** (*int*) -- Number of functions +:Datasets: + - ***func_** (:ref:`function <1d_functions>`) -- Dataset for the + i-th function (indexing starts at 1) + .. _angle_energy: -------------------------- @@ -501,6 +511,19 @@ equiprobable bins. - **skewed** (*int8_t*) -- Whether discrete angles are equi-probable (0) or have a skewed distribution (1). +Mixed Elastic +------------- + +This angle-energy distribution is used when an evaluation specifies both +coherent and incoherent elastic thermal neutron scattering. + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "mixed_elastic" +:Groups: - **coherent** -- Distribution for coherent elastic scattering. The + format is given in :ref:`angle_energy`. + - **incoherent** -- Distribution for incoherent elastic scattering. + The format is given in :ref:`angle_energy`. + .. _energy_distribution: -------------------- diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 287738774..c74769be6 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -116,6 +116,7 @@ Angle-Energy Distributions IncoherentElasticAE IncoherentElasticAEDiscrete IncoherentInelasticAEDiscrete + MixedElasticAE Resonance Data -------------- diff --git a/openmc/data/function.py b/openmc/data/function.py index 7a73987df..751f6cdb5 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -581,6 +581,8 @@ class Sum(EqualityMixin): def to_hdf5(self, group, name='xy'): """Write sum of functions to an HDF5 group + .. versionadded:: 0.13.1 + Parameters ---------- group : h5py.Group @@ -599,6 +601,8 @@ class Sum(EqualityMixin): def from_hdf5(cls, group): """Generate sum of functions from an HDF5 group + .. versionadded:: 0.13.1 + Parameters ---------- group : h5py.Group diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index d7dfe92f9..b7ca349f1 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -36,6 +36,9 @@ class CoherentElasticAE(AngleEnergy): def to_hdf5(self, group, xs_group=None): """Write coherent elastic distribution to an HDF5 group + .. versionchanged:: 0.13.1 + The *xs_group* argument was added. + Parameters ---------- group : h5py.Group @@ -220,6 +223,8 @@ class IncoherentInelasticAE(CorrelatedAngleEnergy): class MixedElasticAE(AngleEnergy): """Secondary distribution for mixed coherent/incoherent thermal elastic + .. versionadded:: 0.13.1 + Parameters ---------- coherent : AngleEnergy From 0f03710e367094db1a822f00db7a1e857ea89dcc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 21 Jul 2022 16:09:49 -0500 Subject: [PATCH 41/83] Fix hdf5 write/read for coherent and mixed elastic --- openmc/data/function.py | 2 +- openmc/data/thermal.py | 2 +- openmc/data/thermal_angle_energy.py | 35 ++++++++++++++++++++--------- 3 files changed, 26 insertions(+), 13 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index 751f6cdb5..b0390d19c 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -544,7 +544,7 @@ class Combination(EqualityMixin): self._operations = operations -class Sum(EqualityMixin): +class Sum(Function1D): """Sum of multiple functions. This class allows you to create a callable object which represents the sum diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 6cf5a4e9e..3d5f2df64 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -701,7 +701,7 @@ class ThermalScattering(EqualityMixin): energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] coherent_xs = CoherentElastic(energy, P*EV_PER_MEV) - coherent_dist = CoherentElasticAE(xs) + coherent_dist = CoherentElasticAE(coherent_xs) # Coherent elastic shouldn't have angular distributions listed n_mu = ace.nxs[6] + 1 diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index b7ca349f1..1f37863ea 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -2,6 +2,7 @@ import numpy as np from .angle_energy import AngleEnergy from .correlated import CorrelatedAngleEnergy +import openmc.data class CoherentElasticAE(AngleEnergy): @@ -33,25 +34,37 @@ class CoherentElasticAE(AngleEnergy): def __init__(self, coherent_xs): self.coherent_xs = coherent_xs - def to_hdf5(self, group, xs_group=None): + def to_hdf5(self, group): """Write coherent elastic distribution to an HDF5 group - .. versionchanged:: 0.13.1 - The *xs_group* argument was added. - Parameters ---------- group : h5py.Group HDF5 group to write to - xs_group : h5py.Group, optional - Group containing 'xs' dataset """ group.attrs['type'] = np.string_('coherent_elastic') - if xs_group is not None: - group['coherent_xs'] = xs_group['xs'] - else: - group['coherent_xs'] = group.parent['xs'] + self.coherent_xs.to_hdf5(group, 'coherent_xs') + + @classmethod + def from_hdf5(cls, group): + """Generate coherent elastic distribution from HDF5 data + + .. versionadded:: 0.13.1 + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.CoherentElasticAE + Coherent elastic distribution + + """ + coherent_xs = openmc.data.CoherentElastic.from_hdf5(group['coherent_xs']) + return cls(coherent_xs) class IncoherentElasticAE(AngleEnergy): @@ -255,7 +268,7 @@ class MixedElasticAE(AngleEnergy): """ group.attrs['type'] = np.string_('mixed_elastic') coherent_group = group.create_group('coherent') - self.coherent.to_hdf5(coherent_group, group.parent) + self.coherent.to_hdf5(coherent_group) incoherent_group = group.create_group('incoherent') self.incoherent.to_hdf5(incoherent_group) From 425dba5d25f6a3d737838f9cb408dd57bf45c24e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Jul 2022 07:46:34 -0500 Subject: [PATCH 42/83] Add test for mixed elastic thermal scattering --- openmc/mixin.py | 7 +- tests/unit_tests/test_data_thermal.py | 105 ++++++++++++++++++++++++++ 2 files changed, 110 insertions(+), 2 deletions(-) diff --git a/openmc/mixin.py b/openmc/mixin.py index 516162464..31c26ec76 100644 --- a/openmc/mixin.py +++ b/openmc/mixin.py @@ -14,8 +14,11 @@ class EqualityMixin: def __eq__(self, other): if isinstance(other, type(self)): for key, value in self.__dict__.items(): - if not np.array_equal(value, other.__dict__.get(key)): - return False + if isinstance(value, np.ndarray): + if not np.array_equal(value, other.__dict__.get(key)): + return False + else: + return value == other.__dict__.get(key) else: return False diff --git a/tests/unit_tests/test_data_thermal.py b/tests/unit_tests/test_data_thermal.py index dc4d24628..61ae61b24 100644 --- a/tests/unit_tests/test_data_thermal.py +++ b/tests/unit_tests/test_data_thermal.py @@ -262,3 +262,108 @@ def test_get_thermal_name(): # Names that don't remotely match anything assert f('boogie_monster') == 'c_boogie_monster' + + +@pytest.fixture +def fake_mixed_elastic(): + fake_tsl = openmc.data.ThermalScattering("c_D_in_7LiD", 1.9968, 4.9, [0.0253]) + fake_tsl.nuclides = ['H2'] + + # Create elastic reaction + bragg_edges = [0.00370672, 0.00494229, 0.00988458, 0.01359131, 0.01482688, + 0.01976918, 0.02347589, 0.02471147, 0.02965376, 0.03336048, + 0.03953834, 0.04324506, 0.04448063, 0.04942292, 0.05312964, + 0.05436522, 0.05930751, 0.06301423, 0.0642498 , 0.06919209, + 0.07289881, 0.07907667, 0.08278339, 0.08401896, 0.08896126, + 0.09266798, 0.09390355, 0.09884584, 0.1025526 , 0.1037882 , + 0.1087305 , 0.1124372 , 0.1186151 , 0.1223218 , 0.1235574 , + 0.1284997 , 0.1322064 , 0.133442 , 0.142091 , 0.1433266 , + 0.1482688 , 0.1519756 , 0.1581534 , 0.1618601 , 0.1630957 , + 0.168038 , 0.1717447 , 0.1729803 , 0.1779226 , 0.1816293 , + 0.1828649 , 0.1878072 , 0.1915139 , 0.1976918 , 0.2026341 , + 0.2075763 , 0.2125186 , 0.2174609 , 0.2224032 , 0.2273455 , + 0.2421724 , 0.2471147 , 0.252057 , 0.2569993 , 0.2619415 , + 0.2668838 , 0.2767684 , 0.2817107 , 0.2915953 , 0.3064222 , + 0.3261913 , 0.366965] + factors = [0.00375735, 0.01386287, 0.02595574, 0.02992438, 0.03549502, + 0.03855745, 0.04058831, 0.04986305, 0.05703106, 0.05855471, + 0.06078031, 0.06212291, 0.06656602, 0.06930339, 0.0697072 , + 0.07201456, 0.07263853, 0.07313129, 0.07465531, 0.07714482, + 0.07759976, 0.077809 , 0.07790282, 0.07927957, 0.08013058, + 0.08026637, 0.08073475, 0.08112202, 0.08123039, 0.08187171, + 0.08213756, 0.08218236, 0.08236572, 0.08240729, 0.08259795, + 0.08297893, 0.08300455, 0.08314566, 0.08315611, 0.08337715, + 0.08350026, 0.08350663, 0.08352815, 0.08353776, 0.0836098 , + 0.08367017, 0.08367361, 0.0837242 , 0.08375069, 0.08375227, + 0.08377006, 0.08381488, 0.08381644, 0.08382698, 0.08386266, + 0.08387756, 0.08388445, 0.08388974, 0.08390341, 0.08391088, + 0.08391695, 0.08392361, 0.08392684, 0.08392818, 0.08393161, + 0.08393546, 0.08393685, 0.08393801, 0.08393976, 0.08394167, + 0.08394288, 0.08394398] + coherent_xs = openmc.data.CoherentElastic(bragg_edges, factors) + incoherent_xs = openmc.data.Tabulated1D([0.00370672, 0.00370672], [0.00370672, 0.00370672]) + elastic_xs = {'294K': openmc.data.Sum((coherent_xs, incoherent_xs))} + coherent_dist = openmc.data.CoherentElasticAE(coherent_xs) + incoherent_dist = openmc.data.IncoherentElasticAEDiscrete([ + [-0.6, -0.18, 0.18, 0.6], [-0.6, -0.18, 0.18, 0.6] + ]) + elastic_dist = {'294K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist)} + fake_tsl.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist) + + # Create inelastic reaction + inelastic_xs = {'294K': openmc.data.Tabulated1D([1.0e-5, 4.9], [13.4, 3.35])} + breakpoints = [3] + interpolation = [2] + energy = [1.0e-5, 4.3e-2, 4.9] + energy_out = [ + openmc.data.Tabular([0.0002, 0.067, 0.146, 0.366], [0.25, 0.25, 0.25, 0.25]), + openmc.data.Tabular([0.0001, 0.009, 0.137, 0.277], [0.25, 0.25, 0.25, 0.25]), + openmc.data.Tabular([0.0579, 4.555, 4.803, 4.874], [0.25, 0.25, 0.25, 0.25]), + ] + for eout in energy_out: + eout.normalize() + eout.c = eout.cdf() + discrete = openmc.stats.Discrete([-0.9, -0.6, -0.3, -0.1, 0.1, 0.3, 0.6, 0.9], [1/8]*8) + discrete.c = discrete.cdf()[1:] + mu = [[discrete]*4]*3 + inelastic_dist = {'294K': openmc.data.IncoherentInelasticAE( + breakpoints, interpolation, energy, energy_out, mu)} + inelastic = openmc.data.ThermalScatteringReaction(inelastic_xs, inelastic_dist) + fake_tsl.inelastic = inelastic + + return fake_tsl + + +def test_mixed_elastic(fake_mixed_elastic, run_in_tmpdir): + # Write data to HDF5 and then read back + original = fake_mixed_elastic + original.export_to_hdf5('c_D_in_7LiD.h5') + copy = openmc.data.ThermalScattering.from_hdf5('c_D_in_7LiD.h5') + + # Make sure data did not change as a result of HDF5 writing/reading + assert original == copy + + # Create modified cross_sections.xml file that includes the above data + xs = openmc.data.DataLibrary.from_xml() + xs.register_file('c_D_in_7LiD.h5') + xs.export_to_xml('cross_sections_mixed.xml') + + # Create a minimal model that includes the new data and run it + mat = openmc.Material() + mat.add_nuclide('H2', 1.0) + mat.add_nuclide('Li7', 1.0) + mat.set_density('g/cm3', 1.0) + mat.add_s_alpha_beta('c_D_in_7LiD') + sph = openmc.Sphere(r=10.0, boundary_type="vacuum") + cell = openmc.Cell(fill=mat, region=-sph) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.materials = openmc.Materials([mat]) + model.materials.cross_sections = "cross_sections_mixed.xml" + model.settings.particles = 1000 + model.settings.batches = 10 + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.Source( + energy=openmc.stats.Discrete([3.0], [1.0]) # 3 eV source + ) + model.run() From 3b33fb830624d8c0e7959e6326214f52ae988eb0 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 22 Jul 2022 07:55:22 -0500 Subject: [PATCH 43/83] Adding comment as suggested by @paulromano --- openmc/mesh.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index d08213f77..4af701892 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -248,6 +248,10 @@ class StructuredMesh(MeshBase): # create VTK arrays for each of # the data sets + + # maintain a list of the datasets as added + # to the VTK arrays to ensure they persist + # in memory until the file is written datasets_out = [] for label, dataset in datasets.items(): dataset = np.asarray(dataset).flatten() From 717115938df14c8561da1052dd2e849bdc26b4bb Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Mon, 25 Jul 2022 16:25:29 +0100 Subject: [PATCH 44/83] updating fill id to match --- tests/regression_tests/dagmc/universes/inputs_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/dagmc/universes/inputs_true.dat b/tests/regression_tests/dagmc/universes/inputs_true.dat index 1a5ca7486..2165a78b6 100644 --- a/tests/regression_tests/dagmc/universes/inputs_true.dat +++ b/tests/regression_tests/dagmc/universes/inputs_true.dat @@ -1,6 +1,6 @@ - + 24.0 24.0 From b93ab95ba14e7af55442cc9d9ec154e1d34b48ed Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 31 May 2022 08:33:43 -0500 Subject: [PATCH 45/83] Adding more mesh interrogation options for the libmesh Finishing methods for connectivity and coordinates. Writing vertices and connectivity to statepoint file Loading vertices and connectivity from statepoint. Correcting string repr Correcting connectivity length Adding method to write the mesh elements to VTK with data applied. Updating hdf5 output to include element types Adding support for hex elements when writing unstructured meshes to VTK Adding simple check for VTK writing if the module is present Removing centroids from the statepoint file and Python UM class Updating test check for vtk Adding warning for skipped elements. Correcting element type Adding warning for skipped elements. Using an enum to indicate element types for readability Updating to element types on the Python side as well Handling integer data applied to VTK files. Doc updates for Python API UM class Incrementing statepoint version number Refactor of unstructured mesh tests to extract model Updating inputs for floating point surface coefficients Adding test for hexes and refactoring comparison funcs Updating reference mesh files Adding reference file for the hexes test case Passing test for hex mesh Adding inputs for the hexes test case. Adding hex test meshes. Skipping hex mesh test if not built with libmesh Adding small VTK write tests for unstructured mesh. Allowing file path to be a pathlib path. Adding skips if libmesh or dagmc not enabled Adding a few comments to test file Changing where conversion to str happens for mesh filename. Setting output to false. Removing VTK check from unstructured mesh regression test Removnig VTK test files for regression test -- too large Adding __init__.py file for pytest --- include/openmc/constants.h | 2 +- include/openmc/mesh.h | 39 +++ openmc/mesh.py | 239 +++++++++----- openmc/statepoint.py | 2 +- src/mesh.cpp | 134 +++++++- .../unstructured_mesh/inputs_true.dat | 91 ++++++ .../unstructured_mesh/inputs_true0.dat | 12 +- .../unstructured_mesh/inputs_true1.dat | 12 +- .../unstructured_mesh/inputs_true10.dat | 12 +- .../unstructured_mesh/inputs_true11.dat | 12 +- .../unstructured_mesh/inputs_true12.dat | 12 +- .../unstructured_mesh/inputs_true13.dat | 12 +- .../unstructured_mesh/inputs_true14.dat | 12 +- .../unstructured_mesh/inputs_true15.dat | 12 +- .../unstructured_mesh/inputs_true2.dat | 12 +- .../unstructured_mesh/inputs_true3.dat | 12 +- .../unstructured_mesh/inputs_true8.dat | 12 +- .../unstructured_mesh/inputs_true9.dat | 12 +- .../unstructured_mesh/test.py | 213 ++++++++----- .../unstructured_mesh/test_mesh_hexes.e | 1 + .../unstructured_mesh/test_mesh_hexes.exo | Bin 0 -> 79668 bytes tests/unit_tests/mesh_to_vtk/__init__.py | 0 tests/unit_tests/mesh_to_vtk/hexes.exo | Bin 0 -> 4700 bytes .../mesh_to_vtk/libmesh_hexes_ref.vtk | 76 +++++ .../mesh_to_vtk/libmesh_tets_ref.vtk | 292 ++++++++++++++++++ .../unit_tests/mesh_to_vtk/moab_tets_ref.vtk | 1 + tests/unit_tests/mesh_to_vtk/test.py | 70 +++++ tests/unit_tests/mesh_to_vtk/tets.exo | Bin 0 -> 7708 bytes tests/unit_tests/mesh_to_vtk/umesh.vtk | 292 ++++++++++++++++++ 29 files changed, 1344 insertions(+), 252 deletions(-) create mode 100644 tests/regression_tests/unstructured_mesh/inputs_true.dat create mode 120000 tests/regression_tests/unstructured_mesh/test_mesh_hexes.e create mode 100644 tests/regression_tests/unstructured_mesh/test_mesh_hexes.exo create mode 100644 tests/unit_tests/mesh_to_vtk/__init__.py create mode 100644 tests/unit_tests/mesh_to_vtk/hexes.exo create mode 100644 tests/unit_tests/mesh_to_vtk/libmesh_hexes_ref.vtk create mode 100644 tests/unit_tests/mesh_to_vtk/libmesh_tets_ref.vtk create mode 120000 tests/unit_tests/mesh_to_vtk/moab_tets_ref.vtk create mode 100644 tests/unit_tests/mesh_to_vtk/test.py create mode 100644 tests/unit_tests/mesh_to_vtk/tets.exo create mode 100644 tests/unit_tests/mesh_to_vtk/umesh.vtk diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 73f96ff95..a7362ea9e 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -24,7 +24,7 @@ using double_4dvec = vector>>>; constexpr int HDF5_VERSION[] {3, 0}; // Version numbers for binary files -constexpr array VERSION_STATEPOINT {17, 0}; +constexpr array VERSION_STATEPOINT {18, 0}; constexpr array VERSION_PARTICLE_RESTART {2, 0}; constexpr array VERSION_TRACK {3, 0}; constexpr array VERSION_SUMMARY {6, 0}; diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 21c58862d..1ff20dc11 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -36,6 +36,12 @@ namespace openmc { +//============================================================================== +// Constants +//============================================================================== + +enum class ElementType { UNSUPPORTED=-1, LINEAR_TET, LINEAR_HEX }; + //============================================================================== // Global variables //============================================================================== @@ -485,6 +491,23 @@ public: //! \return The centroid of the bin virtual Position centroid(int bin) const = 0; + //! Get the number of vertices in the mesh + // + //! \return Number of vertices + virtual int n_vertices() const = 0; + + //! Retrieve a vertex of the mesh + // + //! \param[in] vertex ID + //! \return vertex coordinates + virtual Position vertex(int id) const = 0; + + //! Retrieve connectivity of a mesh element + // + //! \param[in] element ID + //! \return element connectivity as IDs of the vertices + virtual std::vector connectivity(int id) const = 0; + //! Get the volume of a mesh bin // //! \param[in] bin Bin to return the volume for @@ -557,6 +580,12 @@ public: Position centroid(int bin) const override; + int n_vertices() const override; + + Position vertex(int id) const override; + + std::vector connectivity(int id) const override; + double volume(int bin) const override; private: @@ -621,6 +650,9 @@ private: //! \return MOAB EntityHandle of tet moab::EntityHandle get_ent_handle_from_bin(int bin) const; + //! Get a vertex index into the global range from a handle + int get_vert_idx_from_handle(moab::EntityHandle vert) const; + //! Get the bin for a given mesh cell index // //! \param[in] idx Index of the mesh cell. @@ -655,6 +687,7 @@ private: // Data members moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh + moab::Range verts_; //!< Range of vertex EntityHandle's in the mesh moab::EntityHandle tetset_; //!< EntitySet containing all tetrahedra moab::EntityHandle kdtree_root_; //!< Root of the MOAB KDTree std::shared_ptr mbi_; //!< MOAB instance @@ -701,6 +734,12 @@ public: Position centroid(int bin) const override; + int n_vertices() const override; + + Position vertex(int id) const override; + + std::vector connectivity(int id) const override; + double volume(int bin) const override; private: diff --git a/openmc/mesh.py b/openmc/mesh.py index 4af701892..601d0ac94 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2,6 +2,7 @@ from abc import ABC, abstractmethod from collections.abc import Iterable from math import pi from numbers import Real, Integral +from pathlib import Path import warnings from xml.etree import ElementTree as ET @@ -1440,7 +1441,7 @@ class UnstructuredMesh(MeshBase): Parameters ---------- - filename : str + filename : str or pathlib.Path Location of the unstructured mesh file library : {'moab', 'libmesh'} Mesh library used for the unstructured mesh tally @@ -1468,20 +1469,41 @@ class UnstructuredMesh(MeshBase): be generated for this mesh volumes : Iterable of float Volumes of the unstructured mesh elements + centroids : np.narray (3, n_elements) + Centroids of the mesh elements + + vertices : np.ndarray (3,) + Coordinates of the mesh vertices + + .. versionadded:: 0.13 + + connectivity : np.ndarray (8, n_elements) + Connectivity of the elements + + .. versionadded:: 0.13 + + element_types : Iterable of integers + Mesh element types + + .. versionadded:: 0.13 total_volume : float Volume of the unstructured mesh in total - centroids : Iterable of tuple - An iterable of element centroid coordinates, e.g. [(0.0, 0.0, 0.0), - (1.0, 1.0, 1.0), ...] """ + + _UNSUPPORTED_ELEM = -1 + _LINEAR_TET = 0 + _LINEAR_HEX = 1 + def __init__(self, filename, library, mesh_id=None, name='', length_multiplier=1.0): super().__init__(mesh_id, name) self.filename = filename self._volumes = None - self._centroids = None + self._n_elements = None + self._conectivity = None + self._vertices = None self.library = library - self._output = True + self._output = False self.length_multiplier = length_multiplier @property @@ -1490,7 +1512,7 @@ class UnstructuredMesh(MeshBase): @filename.setter def filename(self, filename): - cv.check_type('Unstructured Mesh filename', filename, str) + cv.check_type('Unstructured Mesh filename', filename, (str, Path)) self._filename = filename @property @@ -1546,12 +1568,33 @@ class UnstructuredMesh(MeshBase): def centroids(self): return self._centroids + @property + def vertices(self): + return self._vertices + + @property + def connectivity(self): + return self._connectivity + + @property + def element_types(self): + return self._element_types + + @property + def centroids(self): + return np.array([self.centroid(i) for i in range(self.n_elements)]) + @property def n_elements(self): - if self._centroids is None: + if self._n_elements is None: raise RuntimeError("No information about this mesh has " "been loaded from a statepoint file.") - return len(self._centroids) + return self._n_elements + + @n_elements.setter + def n_elements(self, val): + cv.check_type('Number of elements', val, Integral) + self._n_elements = val @centroids.setter @@ -1579,23 +1622,27 @@ class UnstructuredMesh(MeshBase): def n_dimension(self): return 3 - @property - def vertices(self): - raise NotImplementedError("Vertices for UnstructuredMesh objects are " - "not yet available") - def __repr__(self): string = super().__repr__() string += '{: <16}=\t{}\n'.format('\tFilename', self.filename) - string += '{: <16}=\t{}\n'.format('\tMesh Library', self.mesh_lib) + string += '{: <16}=\t{}\n'.format('\tMesh Library', self.library) if self.length_multiplier != 1.0: string += '{: <16}=\t{}\n'.format('\tLength multiplier', self.length_multiplier) return string - def write_data_to_vtk(self, filename, datasets, volume_normalization=True): - """Map data to the unstructured mesh element centroids - to create a VTK point-cloud dataset. + def centroid(self, bin): + """Return the vertex averaged centroid of an element + """ + conn = self.connectivity[bin] + coords = self.vertices[conn] + return coords.mean(axis=0) + + def write_vtk_mesh(self, **kwargs): + """Map data to unstructured VTK mesh elements. + + .. deprecated:: 0.13 + Use :func:`UnstructuredMesh.write_data_to_vtk` instead. Parameters ---------- @@ -1607,83 +1654,105 @@ class UnstructuredMesh(MeshBase): volume_normalization : bool Whether or not to normalize the data by the volume of the mesh elements - - Raises - ------ - RuntimeError - when the size of a dataset doesn't match the number of cells """ + warnings.warn( + "The 'UnstructuredMesh.write_vtk_mesh' method has been renamed " + "to 'write_data_to_vtk' and will be removed in a future version " + " of OpenMC.", FutureWarning + ) + self.write_data_to_vtk(**kwargs) + def write_data_to_vtk(self, filename=None, datasets=None, volume_normalization=True): + """Map data to unstructured VTK mesh elements. + + Parameters + ---------- + filename : str or pathlib.Path + Name of the VTK file to write + datasets : dict + Dictionary whose keys are the data labels + and values are numpy appropriately sized arrays + of the data + volume_normalization : bool + Whether or not to normalize the data by the + volume of the mesh elements + """ import vtk - from vtk.util import numpy_support as vtk_npsup + from vtk.util import numpy_support as nps - if self.centroids is None: - raise RuntimeError("No centroid information is present on this " - "unstructured mesh. Please load this " - "information from a relevant statepoint file.") + if self.connectivity is None or self.vertices is None: + raise RuntimeError('This mesh has not been ' + 'loaded from a statepoint file.') - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") + if filename is None: + filename = f'mesh_{self.id}.vtk' - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) + writer = vtk.vtkUnstructuredGridWriter() + + writer.SetFileName(str(filename)) + + grid = vtk.vtkUnstructuredGrid() + + vtk_pnts = vtk.vtkPoints() + vtk_pnts.SetData(nps.numpy_to_vtk(self.vertices)) + grid.SetPoints(vtk_pnts) + + n_skipped = 0 + elems = [] + for elem_type, conn in zip(self.element_types, self.connectivity): + if elem_type == self._LINEAR_TET: + elem = vtk.vtkTetra() + elif elem_type == self._LINEAR_HEX: + elem = vtk.vtkHexahedron() + elif elem_type == self._UNSUPPORTED_ELEM: + n_skipped += 1 else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - cv.check_type('label', label, str) + raise RuntimeError(f'Invalid element type {elem_type} found') + for i, c in enumerate(conn): + if c == -1: + break + elem.GetPointIds().SetId(i, c) + elems.append(elem) - # create data arrays for the cells/points - cell_dim = 1 - vertices = vtk.vtkCellArray() - points = vtk.vtkPoints() + if n_skipped > 0: + warnings.warn(f'{n_skipped} elements were not written because ' + 'they are not of type linear tet/hex') - for centroid in self.centroids: - # create a point for each centroid - point_id = points.InsertNextPoint(centroid * self.length_multiplier) - # create a cell of type "Vertex" for each point - cell_id = vertices.InsertNextCell(cell_dim, (point_id,)) + for elem in elems: + grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds()) - # create a VTK data object - poly_data = vtk.vtkPolyData() - poly_data.SetPoints(points) - poly_data.SetVerts(vertices) - - # strange VTK nuance: - # data must be held in some container - # until the vtk file is written - data_holder = [] - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() + # check that datasets are the correct size + if datasets is not None: + for name, data in datasets.items(): + if data.shape != (self.dimension,): + raise ValueError(f'Cannot apply dataset {name} with ' + f'shape {data.shape} to mesh {self.id} ' + f'with dimensions {self.dimension}') if volume_normalization: - dataset /= self.volumes.flatten() + for name, data in datasets.items(): + if np.issubdtype(data.dtype, np.integer): + warnings.warn(f'Integer data set {name} will ' + 'not be volume-normalized.') + continue + data /= self.volumes - array = vtk.vtkDoubleArray() - array.SetName(label) - array.SetNumberOfComponents(1) - array.SetArray(vtk_npsup.numpy_to_vtk(dataset), - dataset.size, - True) + # add data to the mesh + for name, data in datasets.items(): + arr = vtk.vtkDoubleArray() + arr.SetName(name) + arr.SetNumberOfTuples(data.size) - data_holder.append(dataset) - poly_data.GetPointData().AddArray(array) + for i in range(data.size): + arr.SetTuple1(i, data.flat[i]) + grid.GetCellData().AddArray(arr) - # set filename - if not filename.endswith(".vtk"): - filename += ".vtk" + if vtk.VTK_MAJOR_VERSION == 5: + grid.update() + writer.SetInput(grid) + else: + writer.SetInputData(grid) - writer = vtk.vtkGenericDataObjectWriter() - writer.SetFileName(filename) - writer.SetInputData(poly_data) writer.Write() @classmethod @@ -1694,10 +1763,14 @@ class UnstructuredMesh(MeshBase): mesh = cls(filename, library, mesh_id=mesh_id) vol_data = group['volumes'][()] - centroids = group['centroids'][()] mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],)) - mesh.centroids = np.reshape(centroids, (vol_data.shape[0], 3)) - mesh.size = mesh.volumes.size + mesh.n_elements = mesh.volumes.size + + vertices = group['vertices'][()] + mesh._vertices = vertices.reshape((-1, 3)) + connectvity = group['connectivity'][()] + mesh._connectivity = connectvity.reshape((-1, 8)) + mesh._element_types = group['element_types'][()] if 'length_multiplier' in group: mesh.length_multiplier = group['length_multiplier'][()] @@ -1719,7 +1792,7 @@ class UnstructuredMesh(MeshBase): element.set("type", "unstructured") element.set("library", self._library) subelement = ET.SubElement(element, "filename") - subelement.text = self.filename + subelement.text = str(self.filename) if self._length_multiplier != 1.0: element.set("length_multiplier", str(self.length_multiplier)) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 630675816..138c0c718 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -11,7 +11,7 @@ from uncertainties import ufloat import openmc import openmc.checkvalue as cv -_VERSION_STATEPOINT = 17 +_VERSION_STATEPOINT = 18 class StatePoint: diff --git a/src/mesh.cpp b/src/mesh.cpp index e8f79aa37..e701f954e 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -220,21 +220,56 @@ void UnstructuredMesh::to_hdf5(hid_t group) const write_dataset(mesh_group, "type", mesh_type); write_dataset(mesh_group, "filename", filename_); write_dataset(mesh_group, "library", this->library()); - // write volume of each element - vector tet_vols; - xt::xtensor centroids({static_cast(this->n_bins()), 3}); - for (int i = 0; i < this->n_bins(); i++) { - tet_vols.emplace_back(this->volume(i)); - auto c = this->centroid(i); - xt::view(centroids, i, xt::all()) = xt::xarray({c.x, c.y, c.z}); - } - - write_dataset(mesh_group, "volumes", tet_vols); - write_dataset(mesh_group, "centroids", centroids); if (specified_length_multiplier_) write_dataset(mesh_group, "length_multiplier", length_multiplier_); + // write vertex coordinates + xt::xtensor vertices({static_cast(this->n_vertices()), 3}); + for (int i = 0; i < this->n_vertices(); i++) { + auto v = this->vertex(i); + xt::view(vertices, i, xt::all()) = xt::xarray({v.x, v.y, v.z}); + } + write_dataset(mesh_group, "vertices", vertices); + + int num_elem_skipped = 0; + + // write element types and connectivity + vector volumes; + xt::xtensor connectivity ({static_cast(this->n_bins()), 8}); + xt::xtensor elem_types ({static_cast(this->n_bins()), 1}); + for (int i = 0; i < this->n_bins(); i++) { + auto conn = this->connectivity(i); + + volumes.emplace_back(this->volume(i)); + + // write linear tet element + if (conn.size() == 4) { + xt::view(elem_types, i, xt::all()) = static_cast(ElementType::LINEAR_TET); + xt::view(connectivity, i, xt::all()) = xt::xarray({conn[0], conn[1], conn[2], conn[3], + -1, -1, -1, -1}); + // write linear hex element + } else if (conn.size() == 8) { + xt::view(elem_types, i, xt::all()) = static_cast(ElementType::LINEAR_HEX); + xt::view(connectivity, i, xt::all()) = xt::xarray({conn[0], conn[1], conn[2], conn[3], + conn[4], conn[5], conn[6], conn[7]}); + } else { + num_elem_skipped++; + xt::view(elem_types, i, xt::all()) = static_cast(ElementType::UNSUPPORTED); + xt::view(connectivity, i, xt::all()) = xt::xarray({-1, -1, -1, -1, -1, -1, -1, -1}); + } + } + + // warn users that some elements were skipped + if (num_elem_skipped > 0) { + warning(fmt::format("The connectivity of {} elements on mesh {} were not written " + "because they are not of type linear tet/hex.", num_elem_skipped, this->id_)); + } + + write_dataset(mesh_group, "volumes", volumes); + write_dataset(mesh_group, "connectivity", connectivity); + write_dataset(mesh_group, "element_types", elem_types); + close_group(mesh_group); } @@ -1775,6 +1810,14 @@ void MOABMesh::initialize() filename_); } + // set member range of vertices + int vertex_dim = 0; + rval = mbi_->get_entities_by_dimension(0, vertex_dim, verts_); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get all vertex handles"); + } + + // make an entity set for all tetrahedra // this is used for convenience later in output rval = mbi_->create_meshset(moab::MESHSET_SET, tetset_); @@ -2124,6 +2167,14 @@ std::pair, vector> MOABMesh::plot( return {}; } +int MOABMesh::get_vert_idx_from_handle(moab::EntityHandle vert) const { + int idx = vert - verts_[0]; + if (idx >= n_vertices()) { + fatal_error(fmt::format("Invalid vertex idx {} (# vertices {})", idx, n_vertices())); + } + return idx; +} + int MOABMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { int bin = eh - ehs_[0]; @@ -2191,6 +2242,46 @@ Position MOABMesh::centroid(int bin) const return {centroid[0], centroid[1], centroid[2]}; } +int MOABMesh::n_vertices() const { + return verts_.size(); +} + +Position MOABMesh::vertex(int id) const { + + moab::ErrorCode rval; + + moab::EntityHandle vert = verts_[id]; + + moab::CartVect coords; + rval = mbi_->get_coords(&vert, 1, coords.array()); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get the coordinates of a vertex."); + } + + return {coords[0], coords[1], coords[2]}; +} + +std::vector MOABMesh::connectivity(int bin) const { + moab::ErrorCode rval; + + auto tet = get_ent_handle_from_bin(bin); + + // look up the tet connectivity + vector conn; + rval = mbi_->get_connectivity(&tet, 1, conn); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get connectivity of a mesh element."); + return {}; + } + + std::vector verts(4); + for (int i = 0; i < verts.size(); i++) { + verts[i] = get_vert_idx_from_handle(conn[i]); + } + + return verts; +} + std::pair MOABMesh::get_score_tags( std::string score) const { @@ -2394,6 +2485,27 @@ Position LibMesh::centroid(int bin) const return {centroid(0), centroid(1), centroid(2)}; } +int LibMesh::n_vertices() const +{ + return m_->n_nodes(); +} + +Position LibMesh::vertex(int vertex_id) const +{ + const auto node_ref = m_->node_ref(vertex_id); + return {node_ref(0), node_ref(1), node_ref(2)}; +} + +std::vector LibMesh::connectivity(int elem_id) const +{ + std::vector conn; + const auto* elem_ptr = m_->elem_ptr(elem_id); + for (int i = 0; i < elem_ptr->n_nodes(); i++) { + conn.push_back(elem_ptr->node_id(i)); + } + return conn; +} + std::string LibMesh::library() const { return mesh_lib_type; diff --git a/tests/regression_tests/unstructured_mesh/inputs_true.dat b/tests/regression_tests/unstructured_mesh/inputs_true.dat new file mode 100644 index 000000000..e556782ad --- /dev/null +++ b/tests/regression_tests/unstructured_mesh/inputs_true.dat @@ -0,0 +1,91 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_hexes.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true0.dat b/tests/regression_tests/unstructured_mesh/inputs_true0.dat index 1ed60d507..2e2594cde 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true0.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true0.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true1.dat b/tests/regression_tests/unstructured_mesh/inputs_true1.dat index 67dde56a9..84a3b186b 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true1.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true1.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true10.dat b/tests/regression_tests/unstructured_mesh/inputs_true10.dat index 9ca47a356..e9272a333 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true10.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true10.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true11.dat b/tests/regression_tests/unstructured_mesh/inputs_true11.dat index 683e1fade..0b89d280c 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true11.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true11.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true12.dat b/tests/regression_tests/unstructured_mesh/inputs_true12.dat index 8ce9f3cf2..b3673a254 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true12.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true12.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true13.dat b/tests/regression_tests/unstructured_mesh/inputs_true13.dat index bca0bee4c..c466396b6 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true13.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true13.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true14.dat b/tests/regression_tests/unstructured_mesh/inputs_true14.dat index 226331ba8..8c6a5120a 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true14.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true14.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true15.dat b/tests/regression_tests/unstructured_mesh/inputs_true15.dat index a6a084165..1b9abbd83 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true15.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true15.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true2.dat b/tests/regression_tests/unstructured_mesh/inputs_true2.dat index 4b442c757..7136d485a 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true2.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true2.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true3.dat b/tests/regression_tests/unstructured_mesh/inputs_true3.dat index 52e53498e..23900a060 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true3.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true3.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true8.dat b/tests/regression_tests/unstructured_mesh/inputs_true8.dat index e484c95a2..a0cbbcae7 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true8.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true8.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true9.dat b/tests/regression_tests/unstructured_mesh/inputs_true9.dat index 5e83d71bd..930bd8579 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true9.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true9.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/test.py b/tests/regression_tests/unstructured_mesh/test.py index daebe7948..090a47599 100644 --- a/tests/regression_tests/unstructured_mesh/test.py +++ b/tests/regression_tests/unstructured_mesh/test.py @@ -1,6 +1,8 @@ +import filecmp import glob from itertools import product import os +import warnings import openmc import openmc.lib @@ -9,15 +11,29 @@ import numpy as np import pytest from tests.testing_harness import PyAPITestHarness -TETS_PER_VOXEL = 12 - class UnstructuredMeshTest(PyAPITestHarness): - def __init__(self, statepoint_name, model, inputs_true, holes): + ELEM_PER_VOXEL = 12 + + def __init__(self, + statepoint_name, + model, + inputs_true='inputs_true.dat', + holes=False, + scale_factor=10.0): super().__init__(statepoint_name, model, inputs_true) self.holes = holes # holes in the test mesh + self.scale_bounding_cell(scale_factor) + + def scale_bounding_cell(self, scale_factor): + geometry = self._model.geometry + for surface in geometry.get_all_surfaces().values(): + if surface.boundary_type != 'vacuum': + continue + for coeff in surface._coefficients: + surface._coefficients[coeff] *= scale_factor def _compare_results(self): with openmc.StatePoint(self._sp_name) as sp: @@ -28,32 +44,27 @@ class UnstructuredMeshTest(PyAPITestHarness): flt = tally.find_filter(openmc.MeshFilter) if isinstance(flt.mesh, openmc.RegularMesh): - reg_mesh_data, reg_mesh_std_dev = self.get_mesh_tally_data(tally) + reg_mesh_data = self.get_mesh_tally_data(tally) if self.holes: reg_mesh_data = np.delete(reg_mesh_data, self.holes) - reg_mesh_std_dev = np.delete(reg_mesh_std_dev, self.holes) else: umesh_tally = tally - unstructured_data, unstructured_std_dev = self.get_mesh_tally_data(tally, True) + unstructured_data = self.get_mesh_tally_data(tally, True) - # we expect these results to be the same to within at least ten - # decimal places - decimals = 10 if umesh_tally.estimator == 'collision' else 8 - np.testing.assert_array_almost_equal(unstructured_data, - reg_mesh_data, - decimals) + # we expect these results to be the same to within at least ten + # decimal places + decimals = 10 if umesh_tally.estimator == 'collision' else 8 + np.testing.assert_array_almost_equal(np.sort(unstructured_data), + np.sort(reg_mesh_data), + decimals) - @staticmethod - def get_mesh_tally_data(tally, structured=False): + def get_mesh_tally_data(self, tally, structured=False): data = tally.get_reshaped_data(value='mean') - std_dev = tally.get_reshaped_data(value='std_dev') if structured: - data.shape = (data.size // TETS_PER_VOXEL, TETS_PER_VOXEL) - std_dev.shape = (std_dev.size // TETS_PER_VOXEL, TETS_PER_VOXEL) + data = data.reshape((-1, self.ELEM_PER_VOXEL)) else: data.shape = (data.size, 1) - std_dev.shape = (std_dev.size, 1) - return np.sum(data, axis=1), np.sum(std_dev, axis=1) + return np.sum(data, axis=1) def _cleanup(self): super()._cleanup() @@ -64,35 +75,11 @@ class UnstructuredMeshTest(PyAPITestHarness): os.remove(f) -param_values = (['libmesh', 'moab'], # mesh libraries - ['collision', 'tracklength'], # estimators - [True, False], # geometry outside of the mesh - [(333, 90, 77), None]) # location of holes in the mesh -test_cases = [] -for i, (lib, estimator, ext_geom, holes) in enumerate(product(*param_values)): - test_cases.append({'library' : lib, - 'estimator' : estimator, - 'external_geom' : ext_geom, - 'holes' : holes, - 'inputs_true' : 'inputs_true{}.dat'.format(i)}) - - -@pytest.mark.parametrize("test_opts", test_cases) -def test_unstructured_mesh(test_opts): - +@pytest.fixture +def model(): openmc.reset_auto_ids() - # skip the test if the library is not enabled - if test_opts['library'] == 'moab' and not openmc.lib._dagmc_enabled(): - pytest.skip("DAGMC (and MOAB) mesh not enbaled in this build.") - - if test_opts['library'] == 'libmesh' and not openmc.lib._libmesh_enabled(): - pytest.skip("LibMesh is not enabled in this build.") - - # skip the tracklength test for libmesh - if test_opts['library'] == 'libmesh' and \ - test_opts['estimator'] == 'tracklength': - pytest.skip("Tracklength tallies are not supported using libmesh.") + model = openmc.Model() ### Materials ### materials = openmc.Materials() @@ -113,7 +100,7 @@ def test_unstructured_mesh(test_opts): water_mat.set_density("atom/b-cm", 0.07416) materials.append(water_mat) - materials.export_to_xml() + model.materials = materials ### Geometry ### fuel_min_x = openmc.XPlane(-5.0, name="minimum x") @@ -149,29 +136,26 @@ def test_unstructured_mesh(test_opts): +clad_min_z & -clad_max_z) clad_cell.fill = zirc_mat - if test_opts['external_geom']: - bounds = (15, 15, 15) - else: - bounds = (10, 10, 10) - - water_min_x = openmc.XPlane(x0=-bounds[0], + # set bounding cell dimension to one + # this will be updated later according to the test case parameters + water_min_x = openmc.XPlane(x0=-1.0, name="minimum x", boundary_type='vacuum') - water_max_x = openmc.XPlane(x0=bounds[0], + water_max_x = openmc.XPlane(x0=1.0, name="maximum x", boundary_type='vacuum') - water_min_y = openmc.YPlane(y0=-bounds[1], + water_min_y = openmc.YPlane(y0=-1.0, name="minimum y", boundary_type='vacuum') - water_max_y = openmc.YPlane(y0=bounds[1], + water_max_y = openmc.YPlane(y0=1.0, name="maximum y", boundary_type='vacuum') - water_min_z = openmc.ZPlane(z0=-bounds[2], + water_min_z = openmc.ZPlane(z0=-1.0, name="minimum z", boundary_type='vacuum') - water_max_z = openmc.ZPlane(z0=bounds[2], + water_max_z = openmc.ZPlane(z0=1.0, name="maximum z", boundary_type='vacuum') @@ -185,9 +169,9 @@ def test_unstructured_mesh(test_opts): water_cell.fill = water_mat # create a containing universe - geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell]) + model.geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell]) - ### Tallies ### + ### Reference Tally ### # create meshes and mesh filters regular_mesh = openmc.RegularMesh() @@ -196,29 +180,11 @@ def test_unstructured_mesh(test_opts): regular_mesh.upper_right = (10.0, 10.0, 10.0) regular_mesh_filter = openmc.MeshFilter(mesh=regular_mesh) - - if test_opts['holes']: - mesh_filename = "test_mesh_tets_w_holes.e" - else: - mesh_filename = "test_mesh_tets.e" - - uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_opts['library']) - uscd_filter = openmc.MeshFilter(mesh=uscd_mesh) - - # create tallies - tallies = openmc.Tallies() - regular_mesh_tally = openmc.Tally(name="regular mesh tally") regular_mesh_tally.filters = [regular_mesh_filter] regular_mesh_tally.scores = ['flux'] - regular_mesh_tally.estimator = test_opts['estimator'] - tallies.append(regular_mesh_tally) - uscd_tally = openmc.Tally(name="unstructured mesh tally") - uscd_tally.filters = [uscd_filter] - uscd_tally.scores = ['flux'] - uscd_tally.estimator = test_opts['estimator'] - tallies.append(uscd_tally) + model.tallies = openmc.Tallies([regular_mesh_tally]) ### Settings ### settings = openmc.Settings() @@ -236,13 +202,92 @@ def test_unstructured_mesh(test_opts): source = openmc.Source(space=space, energy=energy) settings.source = source - model = openmc.model.Model(geometry=geometry, - materials=materials, - tallies=tallies, - settings=settings) + model.settings = settings + + return model + + +param_values = (['libmesh', 'moab'], # mesh libraries + ['collision', 'tracklength'], # estimators + [True, False], # geometry outside of the mesh + [(333, 90, 77), None]) # location of holes in the mesh +test_cases = [] +for i, (lib, estimator, ext_geom, holes) in enumerate(product(*param_values)): + test_cases.append({'library' : lib, + 'estimator' : estimator, + 'external_geom' : ext_geom, + 'holes' : holes, + 'inputs_true' : 'inputs_true{}.dat'.format(i)}) + + +@pytest.mark.parametrize("test_opts", test_cases) +def test_unstructured_mesh_tets(model, test_opts): + # skip the test if the library is not enabled + if test_opts['library'] == 'moab' and not openmc.lib._dagmc_enabled(): + pytest.skip("DAGMC (and MOAB) mesh not enbaled in this build.") + + if test_opts['library'] == 'libmesh' and not openmc.lib._libmesh_enabled(): + pytest.skip("LibMesh is not enabled in this build.") + + # skip the tracklength test for libmesh + if test_opts['library'] == 'libmesh' and \ + test_opts['estimator'] == 'tracklength': + pytest.skip("Tracklength tallies are not supported using libmesh.") + + if test_opts['holes']: + mesh_filename = "test_mesh_tets_w_holes.e" + else: + mesh_filename = "test_mesh_tets.e" + + # add reference mesh tally + regular_mesh_tally = model.tallies[0] + regular_mesh_tally.estimator = test_opts['estimator'] + + # add analagous unstructured mesh tally + uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_opts['library']) + uscd_filter = openmc.MeshFilter(mesh=uscd_mesh) + + # create tallies + uscd_tally = openmc.Tally(name="unstructured mesh tally") + uscd_tally.filters = [uscd_filter] + uscd_tally.scores = ['flux'] + uscd_tally.estimator = test_opts['estimator'] + model.tallies.append(uscd_tally) + + # modify model geometry according to test opts + if test_opts['external_geom']: + scale_factor = 15.0 + else: + scale_factor = 10.0 harness = UnstructuredMeshTest('statepoint.10.h5', model, test_opts['inputs_true'], - test_opts['holes']) + test_opts['holes'], + scale_factor) harness.main() + + +@pytest.mark.skipif(not openmc.lib._libmesh_enabled(), + reason='LibMesh is not enabled in this build.') +def test_unstructured_mesh_hexes(model): + regular_mesh_tally = model.tallies[0] + regular_mesh_tally.estimator = 'collision' + + # add analagous unstructured mesh tally + uscd_mesh = openmc.UnstructuredMesh('test_mesh_hexes.e', 'libmesh') + uscd_filter = openmc.MeshFilter(mesh=uscd_mesh) + + # create tallies + uscd_tally = openmc.Tally(name="unstructured mesh tally") + uscd_tally.filters = [uscd_filter] + uscd_tally.scores = ['flux'] + uscd_tally.estimator = 'collision' + model.tallies.append(uscd_tally) + + harness = UnstructuredMeshTest('statepoint.10.h5', + model) + harness.ELEM_PER_VOXEL = 1 + + harness.main() + diff --git a/tests/regression_tests/unstructured_mesh/test_mesh_hexes.e b/tests/regression_tests/unstructured_mesh/test_mesh_hexes.e new file mode 120000 index 000000000..421bf6a89 --- /dev/null +++ b/tests/regression_tests/unstructured_mesh/test_mesh_hexes.e @@ -0,0 +1 @@ +test_mesh_hexes.exo \ No newline at end of file diff --git a/tests/regression_tests/unstructured_mesh/test_mesh_hexes.exo b/tests/regression_tests/unstructured_mesh/test_mesh_hexes.exo new file mode 100644 index 0000000000000000000000000000000000000000..03682c27564b6857b4397ecded768766e948377b GIT binary patch literal 79668 zcmeFa1(;;j-EUnjGq}4>IYAgDJIr!SdvJdC*=59hge4eCq1HAw#`;8ShWsX{z-NI%M!+tRn@J5~g(=y62$7 z2JdY*ZD9Mp!;UzN+wCy9+kW`)!-fpq&$uxUObL^~6nh>rWYCd=4?BFwu%RqA`a#QD zjblG%oTQKYXWjRpVS5f|ZG#ROHe~4VK}QWcY@b1g4>=Y$>G^yOZr5b{4mk)%IT$Z9 zpW)y=j~;TcesDU74cu?|{@mSsd=0Uvlu#!xmY! zW9gNaT4wpBmR)+;Wx2ibs!K1k>Iz9F>+!vw3MNnQ?@{{?tDeN?jpycOqd)i>xp$@I z`R`#Tl{BI`iC4g>=2(fy1pmgZqNH!dR=WBt;;?6vGMqKSnj{~ zFsvV)G9Tmq7C&s;B3I{y_750y_n36gzvwPGyxBc=V%;6xR$qgg>U-mR5U)30$K-?7 zA9Bz^+wXbM5&4a>uIhnug|_e|K|Ea!Ktoq(!GZb9XfdL;reww zzF+R|v98HH1_up4=8*CBp?({#x8urm{K54J`BvkMw0-&b!g)TJeX_KVGx`}3KR(y# zX`!VxZAW^ZgZK1f!WYbXn?Ky1W<*+QXley*ZEdjP4`v`|&-uXJYFg z(`)_l_&oN0{jvSgoAB!!*VR`Z_s2Nl$B%9Lj66SY?IS$(%*Qof&xG@_%{)JU{rYmh z_K}zCe!0G$&)fU`DEsB>o$&Mj``bUU_2;~=_kV1!_orN6Uant%x&386vB&4U^jp54 zb>7qO=!tG`xqW%LwRpel_E+zJ6Mj7t%j@SW zx3~0PzMm6bu3vx7%jHqc{d>aOTW`;V^ZNDl+}<%0J>JUgDPLc?ocs0jm5(pCuY7zg z&)P^oT`|tm?a$b*JdUfz?rI#xq7Be=Nu2V|(Ix-{nO=drK6-0oPekH=}N)Zysc}R^f}XAjboeJ)ARl@&FzZ&^XrN2i!okLEXQ^lj_U-ntjtIpOt{?`J$N)?1IUo}Swmo??vctzTEq=Z(Dg80#&^o{x|HEB33{ z|H4nKHcxe`~eeTIJob zT#RkS<6^xn)&2de$B(O)qu*SAb$?6Ii}ggWz3Rm{Nj$C^qZjM%U)&$tGilu4$MZF+ z`@gGt&F!l0A5)Fr#N(^6JMws3HOBVDcE@tRDv#GWX{@KZKl(Ao>x+Jj<^9d&SYI`U z|A{>=9v@@$V>`lQEJt3xA8pn1j;hAkzuK$iZ^KV)Uol1=&lh9)e6gM3Gq$ULwf?w2 zmLre*V>!lHf2=p&ueiS)<8_ALc%FEC*^lM&`O4>s<1W0#`eJPS=jT5L{$t?(eFL4< zOPI*mUga@126m^Pan=28edOKM{jD)pj~i2sqvQT+>|ZRmR?CfQY$=u_kL`*)mg8~# z<9VycPg<-m))TKQw!5Wz+;>$!)))QuD*tDd$NjC<{o$v1f7OfE8QUBFVvIZ<7uy?S zY*(y5wkx(Lwm;Svd2H9%YCGcjBabnjH}Xk}elNdO9*>Va#^&>uUaI9W)#GD(apkZ#`CpT>usyX{?!=k zi}gpprCM)y(Tm5oHlKIG&)2VbT=_h)9M2c)jqNDM*!~#fanTDO_4ZctwWjCr9rx$$ zj`hUz$9C0Y&we~lyl=5w`j7S2>*={4#Qo8WG48L|n;#e37q2&#qaTltF+7!HU+XP< z6Mz3F{JzKY#~9lg`$KGhyI(+FUMG4++X*r?$FKzIfcY=63Y#*Ll9b+}^h0abv6Hdi&z_ z^c-7?$Cc}C6nX8T=X%OsY;SBwqxnAM$CcY#_w(xuZ}s-~^by<9*ZSl6V!iRc#d^!U zUSEEmSdRB89v|zAG4gsjy!71O=yf-JMlZ(Puk&&_dgU0;6ZdyC(og4rZhyuUb0Fc^K!kho_JiW zKh_ss#unQXKFYkc=*9Ztd15&pH@fJ@{V~RN)Q=Cp`k_EXNqVc)nPU?JVE7`u?8#MXaw}Puw5d8+kd# z{@YsZf3Y3$xL8kpoaE<^^+aC#%6Z@Y@Hu}Y{dE3k{Er6={AWM-fAe@~{@qJ;9Czj8 zzxUs*lzH_z(R}=bm&?Di>3g{&`zzL0_a~Oec9+}PbGcmK#QyH9{=1(EKY!nO^sDQi zuIl=yD=&|##`wK|{2jsQ;{NzMgXq=2_mB1FG5#(f))W00BaijRdcQBWD}FCt|Nb)1 zFWbASd2?4af3KDs#Tbu^F&Q#q*WFpJ{Y8@=rC! z@2A?U`p)9s>ij)Slu7LZyXi(_c6xvH=nOj{50U>+j_n=_i(B zjK`0w9v|Q1#^d9A#^}fM$5?+~koVX49wEMGi|;L>7x%|pkoWg`xucPO{<}UOn#bw<`FX1A>e1C$ zAGhUkoR`b%%KADrZ+E?1y{_i-<-GnJnDF}Qf6o!mpU)Gm#r0qur!kJM^48*bkK;Ca z^>u1K57)7{(Hrq{YG*97uSL1-+ztk%9{ICV;s-%cR+ES7Jui|{5!Jf zkL*rAV^%*ouPFW-$M5PF$2~pAW?midk*;1W56tWTz8d5Ay3OC?cIF4h_0OnkIr8|u zbZd2gOSPVHo!M6$hcU+Q1$+MVzE~$^}yg$^}Clfx7)YULMe$eya0x zeO%}Jm+SLs&dYyiQhpBfomU>3$9d(eSdQZ` z#yB63s*cyVKaS5h@7Cw(>V2>JaXiQI8pmUthsPB2!TRrSW1bktd(U}a^u~4Nx3@gc zmdA0Nk7FJ+$;j^f)5qiUTF3iuc|4ZO<5knmIIhe42UgG9Sq0VDT8wetjWN!wlp*1x~5d{xJH%qODPUHv{i#`qkld-*soe~%RB<(OZ@80Y2sIGynEKB_bBACN@E z@mgQ6m-lxx(oc1~M=!>Bd>qgJZP)9~Lxc9%`thTq=pYU>h{w?$T`s?F5=XF2l z<>y4%kAJsQ{yU%gc+BrZTt7#@(LBz39>?{+2a4md+CQrEYkhs4?~n6!d7Q`fb$lMQ zSC5b5vK-_58^?Kh9LL{vM;@O8ab4Z`@A?_ooqpo`g*YzTiu>dEjJ!VY<_(MM@7B(0 zLz~BC9DjL?^K$&%NBMUvaU93_I?ls!Jjd^^TRZdn8pn0a&+GGVxWAfAQEnE_)uwy~m&<-^Pd-k|<2IJ-^KtH#|2`~waUPF;d|z50r{!^*m*aXlj^lV-EcbOB$8i{A zoR9l=HPTNU$C3B79QU{6hy8#5?-N}Uer#vHyR{nQc#QLHjQy)Tj^jzH<#=4=F~)gz zO!NF(eQs9AcX^&3+dS_3I&a74LG?aY=h6DSTizeXaoiuvG3NWr`EmT-w?3cd{jfaW z#(6jHk1-w>`Tyb{QS^0+RS>*Ktmk$#%Tak*Ui=rf;iKaRipJX*ak)%iBApUdm) zcwBwl=EueFVd8k4*l`@kYuS(EJ02fnoLA%b6!GuGn)i1lON}^wo9FG3-RbB1Dv$oS zSk4do_x~Hm?drt};{I4uXLX|(+Xv?DiTQ6FzcI%7dQ|ne))=eD#rd>8kM=wt$9Xya z&ays^<9N+ue7{wnm*cqXImY#LjB$L%aUR#-u^heco9Fcb=>?a^WBgr59G9^iy|};r zJClw^`f1Ly=dYIki^t}1RIRr<565vlvGZ%4SFgLdzH0lM%Q^4;{d4)e6I;&jR~(ly z_CD|C&)M?(u-?ai^?BU94zBi_uAbM~@jX{uKi6Z=e*C-VIIiQkjO*dnuKYH}Se}=o z7w6r`qaXLT=hf8zjqjf)HZQIpG2vUC`R+LH#@N5OznuS0{Cu1})yH4<8|UZx`>XQ& z+jBXN_vn|`(XkxIZH(o-xIBL2IF2#Sqj5Zz$7w9b7{_mUT$juLx97+I?r~h6k0-WV z|L@;)q#dT5AJ@lu`MD6waoqJho_in1WuBMgd>z-D zuV2+G|J`UD*Kz;YuB4%PK929NV=Vh|ULThq_FwgDzwtQ`>nV@hI1axnuCwd!ukwBn$9IhJxLA(k_`B*n z8pmsSogKZ%qaXLj{5pPr{eS;B?#TN{aXu}N)AD>==M#Q>-Jft?KR%Az`K#@SejLX! z#_`(Ma@=2z`Em8}ALrGcV|>pQV;tY*7{_zmU-$Fo$MG0r9LI5;9p~G~o4>E>%8!q{ zd3;x^{VyLnt2Gv5c^sGDbH(qW;`dao)#GD~<1~)PI8N*DwX*jzujkD@&*$a&xIEv+ zaT({`GB2;IW4Zo&n7se>{$4BJU;o~?qmh29<2}y5aU9p*gGH~}pF97%#(19U`04dL zaes`l9o6wseE!v5s^g;iexrV$^8MxhP=9`BJ$?Rn_3G~&x~fKuBfHbjxN00#TrZdB z%}I*iL&y9%&fBBAI`R`V#{qf$&E?8N^YgA)?yMf)^qH4in#W}7IU3qA&#`qkG&z1J-dxMyN)W3)5$PX;HyV@V??_=`o9oyXAc)#-)=f~#z zIxy+!j`MPU68^jTh+6cx&OFBTQ+#f>b>{79jWI8>KRlanb>$(> zPqE*|_Qw7j$K|Nz{*bq$zJ3qC<@R=@C#a9No_@>ijq_3S?@sgPvAtdSRyEY$S5@y% z^*Pg4?GLfNu|L%NakZUs-1V{lj%n@>J^hyFz0z-cb9+bh*xs0rmwx;8+8?U-x9R!c z+upwRhrYJAxgW>&#yGM&{gj{IvEBVEzxC&H&+Vlc6LP(6Qr9_UaU;G*s@5Cl zk^1k(vcKl%T6MgQD%Kmz43P>qpSBZj_an6KDIYLUt&LMt2R5X*JHnq zzmJK}iROM>ZEt=4t$a0qKa^M49^b1~+Z%tU8t463Z|$KYKd|W`j)Ryd=eRWwjq>v@ z{En_ZCu$F|-^PAC;r*e1{5ylL`~>xW*w_9r;o~5_m+RMAttSus@4x38;WxIo*S|OF z<9&%Y}@#OsXnWsLE9 z<2)4aQ?KppitX%3Pf?#I<9*BfL-Tyt>$u3<8_yTV#sBVp);x|!=9e9R?-l!ReEu}Q z&+BzuWS?;z8t3sb#dT->_n-Oo#s6-7y}bj{2Gr-Xc%Lc{1M~KE#~U)TTCP4v<2bIb ztHSSqRIJyVpFj5B*3P`b^7EnoeOzpBM=mzsx4ho^dcEiN#=mbF6Q92$(`Iyk-D`XE zdM0)p$LCp`zgzOwG>^A9j%yEjy%YO9i|d5e&b-3b_*}2H_nT^a<9aaOx90ag)&9`i zL*74I26m^P=J%M@dgJ#X_1`}Y$i=4Lp6_2A$My5XagxW{Pjy_x>xnVHyk6(yxc-gv zR~*Ox=hu_5U3veF^ILnb>w%ufaeQ7(?0t*#W_&)h$Ls6LH|p0D`)~GJu6I=PJd}H5 zs`Fd?KDJzMzv}m1<>yTx#D3PxZ}q%!T^*mleeJ*HdYgU+=J%<2e(Uvm z;&pU1(oeCy<>z8~z25tM3s05zUe{68^E8i(c%O>*C2vjhJQSaw_1~*j>*;fQW4-xv zGCp^j+dDFQiR+U1_bKuJBUGM;YQK5C@pm6_J}%Eg@q7CC+^N5(EBD)2j{P)mZ+zd@ z+}^G%B+ke6yl!G1@_OUCq}Tpld5G_o;yCVoy;D8UfBE|suYW*#L*n;f<>y;`{~X6l z{d^Ul#rqVWU%jqRi~X#x&&lvxUiZc4Z#++2@06dv z^?p|D4+E2)?s0MbTfLr`7slsgT(8%kli6$Ov)mv073a6u-g3R)7ar=*hdADP9>?YW z8|TgX@6Yr0mB(A@x7^<5^`p)Ndxr0t$8TOw`FnNRZuRPx3JXC*9=JnSeI?|Tc*Qq^UZ+y@HT`_AcJD0rP@_4JaH@^>Y{o7jk&6fFZSH2bBN5%0LpX<%* zzg~W;_cy-p8<<`mpX=2MN9N_2$HsmZzqf4uz0*YI$FZOFecc!P*?|0nO}}xxcH|dW zJx}SSJRiq#+`NwM%nxtwkH!47wJYc0FUIf~pX2p@mAA9G-!`|mBd@44w!5d_xbB;@ z_&#-#UdLNrfBC-E$8oH;=jUYqUgyo6$NLuN<9faMdFt~;UT<8lw|3^2(;lA>UDa~6 zz2$n#&z<_Xi1q$^*OT@4B^~(*n|}Y@?X5Pn8e-lO$6Fj1_3>6cPxF1tpAWGg#d)*k zi3zVau76wlm>-w>Yg|9Hb~e(_q{Vt;e~9mA`kH5rYHn}O{UO%duUKzfug7}heB9iR z`#ipS?q}t7Nvya0oUGRy$3u18^zl6LcQNrf7yje*v<&P{KV!1C`oCS(&5>1J?Vs_! z$M5yxyfdn+BR`-yR@+zXXXX1A+tu2cAKsUT_}qx|+jqr!%j=z3Z|rA{;&>}Rk7EC= zkGBEoh1cif3GWYW#d^!@hu-V$$S<_m&*JkTKG)0hae2O|KR0{s2j%BZT<^sBxVhf0 zycw~*vH$k<`AbM@xHZo=Iv`O)*IIk@j2Oh z-aX;^L^{g_2z%8_4eN0ynVfo-Cl%^9sxDo%npH^3*(zi|y@<4Xwtm7^~OQ{5y`0 z9Q+^pjd^!`9{2t}^Sg$33kk;f~mod*)-r8!;aZ>Fg@5Yn*nA7J7zP% z%wWfC7Vs0!j@fK5JJ>Os1Lg!fW^=*ZV8?6#bbuYRd0<|!V>Tan1=;~-3&4V4$7~^3 z80?rW0*itjv&CR>u%k1*ECEY`&h)YrEDdIQSqAKw>1A24W2Tqoz>b+-mIpg#dRYPN znCWFjuw$l|mB5afURDM>W_no#?3n3gRj^~Gm({?ISqH2RcFgp$2H5EUy{rj#%=9u4 z?3n4L6YQAjZY{85rn|MF6LhARbzoi4nO@d|^}$Rp8-N`%y=(|}%=EGm*fG<~#$d-x zFPnfJGrepIcFgp$8Q3w?%jRImOfOr29W%Xb33kl%vK81d)63Ri$4oEVfE_cvYzubG z^s*h;G3$iw!H${kb^tq_pt~Kxj?VP56YLB+)5|WfE12nJH?U)-m)*gRnO+8g9W%Y` z0d~yvvM1Ov)5~742bk$)Z?I#gmwmvFnO+8i9W%Y`3wF%(vLDzn)64!~$4oCnz>b+- z4gfo5dN~m6nCaypuw$mXgTaoO?uLRLo#|y590EGi%b{=>nCazkuw$l|;b6y1FGqkK zvpwKQuw$l|qrlD{pqHb;j+tJL0Xt@TITq}g>E$@EW2Tqm!H$_;P5?V*dN~p7nCayt zuw$l|lfjOeUQPi!W_mdlP6jjGod$Nybay(~(V1S(fHOg7dN~Wu1~a{!19r^xaxU00 z)602a$4oEhgB>%yTmW{=^l~BCG1JRMV8=`^7lR!$y<7rz%=B_8*fG<~WnjllFPDQI zGre2^cFazOE5VMLUakT=Cxc$D20LcDy9Vr->F!#vqcgqy0Ima_>E(L(A(-jq2C!qM zml0scOfNTr9W%Y$1S7yqFE@i7GrimbcFgp0E7&p9%WYuCOfR>C9W%Y$0d~yv@*}Wg zrk6Xxj+tKW0y}1Uxf|@5>E#}w>qzj+tJ54R*RfFYkaIGrhbEcFgqh9@sI{-TPq2Om`oE9i8drH}G4~ znO;7GkHAbXAA=n;y?g?8%=Gdp*fG<~@4$|kUVaaD%=GdZ*fG<~=U~T7FJFKiGrjx) z?3n50OR!_6mp_6XGrjx??3n50E3jjxmp_9YGrjx;?3n50uVBYacYgyrX1e=3*y#uT zVG@`W{sI4lui+aon*yc;J7(X)cVNeC8kiRBn0*hUz>e7rFeBJ88x3Q?j@c|QE7&m` z3**3!*&Hw@*fHzQ)27rh8vq?($E*cftB%=xFh6uuSsP<}?v(S57Bif*@Vf1k^Np0@ zp#@iVdd)XlJD?2(%?$%^no%_smHy!wX>UCl01Uf0x-SoLr&Nq@* z-SNKLDd!tHOYYbgT4o~GneH}(^{PyFGlTn!+48V4&`Dl($9`j{oNpwrx?^9mQ_eS% zSKZB)JLP;MdDY$Qxl_(Jl2_faU)m|>8_BEg*jMe8^Nr+HcXQ=V^13#?9btR!Gt=GN zK&PB<=u9sy?Ds7ltM1q@?3DA3 zBGisU#Of4-42 z-K_+4%K1j}s=JkQr<`vjuew_$cgp!j@~XR4bElkdB(J(#EqBWKM)ImVKkMj}^Nr+H zcWdNM<77Ig!qspk_nGO=&q6xod?RJL8<;!gd_!k?Y3YQuc>Gy#E{uTdt4w!m=T14_ zNM3cfPVSWRjpS8#>*h{5-$-6{$LEHfa=wwg>TdnqDd!u>tL`?)opQdByz0)+YC7e7 zBYD-`M!8eYHTc`YDd!u>tM0bRopQdByy|Y-+$rZ9 z$*bEol=F?`Rd>7PPNR!;JPq%_+uUcSJ3iCxl=F?0>26T&l=F?`Rd;*jPC4Ju zS#q~0<6f-)MR)~1fe))pcY6b!a=wwg>TaLhDd!u>tL_HpPC4I5UUj!`?v(S5fvwqMY?3kShCxIQaNnld2V|Fr}0(8=|CrW9mVY_VG1zQ-5Ee9&C$kidn&MFrn@tNPMV{QE#@-W2Tqsz>b+-&ILPWdYK;VnCb33uw$mX8NiO2?#>50 zI@8^ZFcaWRFD)0qgE&{;W2Tom!H$_;t^hk`dYKFCnCaz8uw$mXxxtQ^ z?ydqmI@8Ml=m4GRb+-76v=YgPC4713PAVc@FHD>1A`UW2Tqq!H$_; zwg5Y3dU*ltnCWFpuw$l|7r~C1UbX@|W_o!E?3n3gYp`RcmzTkgnO?R5J7#)$1?-sV zWm~Xgrk7X2j+ySZ13PBVz-wSfXL{Klb^x8}<#l)i%=EG&*fG<~n_$OGFFS!9v*+O# zV8=`^JA)lFz5EjFnCWE~uw$l|x4@2>UUmgLW_tM**fG<~ZeYhuFK>e#GrjB%cFgqh zYp`RcmqB31OfT<%9W&kS0d~xE_b%AcnO^pUy+CJrc@N$PGrjB$cFgqh0oXCq%RXSo zOfSCyJ7#(r40g=)@>{TDrk8!ej+tIQ1UqJW*$?cP>E$D^W2Tq=!H$_;J_b8xdKm(C z%=GdJ*fG<~0bs{WFQ0-PGu<5scFc75JFs&Q91KHY82lbSgU{g$FgpYe1v_ScfG@$0 z*-Bpd~H%>Dv@1v_R(!!cmT>~HXQuw!;C90zvH z{sI34J7&kj31G+UYxoB2n4JhGfgQ7N;XANnb~2m-cFewqQD8@BdN~zN1D)w*G>ic= zy_^nq%=9uA?3n5046tLSmvLanOfP4G9W%XjgB>%yoCS8w@Y32~v}2~1v%!we@Y32+ zbF>*Tr9YrbKIyVfakr<`x3Om|b}PC4I5 zUUfHB?v(S5Tb5& zDd!u>tL|pcopQdByy|X_+$rZ9I@3$*oG@2*HvW7gWxAU?cgp!j@~XQ5xl_(Jl2_ez zk z^Np10Zo%9s=Nmdp?iONPn8)7?x59mJZ}>q`M#^-zeD0LZvl=F?` zRd=i8PC4I5UUj!>?v(S5*P*3-$-6{w{GsVbg_+3!7K1G z_nGN#J)l$0H&Uj%^>e43ZzQj}+aPz!`G(H)(z+pR#QJeI{(K{4y4yH+%K1j}s=G~c zr<`vjue#ebcgp!j@~XSda;KbcB(J*LJa@|ZM)In=Epn%vZzQj}+cI~``9|`pyRCAk zoNpwry4yN;%K3)QlDlmfx5eG-@Fsi!?^T)ZwgWold?R_)-S)Xt&Nq@*-R+P&<$NP~ z)!mM{Q_eS%SKaNDJLP;MdDY#{xl_(Jl2_gBk~`&mBYD-`uDMgrHuQ_eSZmfQ_zJR& z(Z1n!8?=I%?v4jKX^!?Cx7)#vneI*iI%$qJirf9bj+yRG1UhMsHk#Z0VG_WZURqCr zlgY`(pKqj0cc5_GBF!LhW2Tp>!H${k&H_7TdYK07nCb3puw$mXX~B+}?#=-_ zI@8^BFg@TbxjUEfJTTMC3}8pEx;r22nCWFkuw#ao)(gOnnO19^1W2To&z>b+-W&=BBdbt$rnCWG9uw$l|%fODA?&bhH zX1co^?C4A{bHZGpGre2^SAv;d<_0@vdbtYhnCWEz*fG<~)nLa=FCAdVOfT1f9W%Ym z19r^xaxK^~)62YI$4oCj06S)SnGfuk>E$}GW2Tq+!H$_;t_M42dRYMMnCaz*V8={% z3xXXp-Q56obf%YuU}4aiUPi!;V5XNvz>b+-ZUQ@IdRY|gnCazauw$l|#lVi4UTy(9 zW_np1?3n50RS7FSmmoGrcSYcFgp02iP&w%hF)SOfNqI zJ7&6D2JD#W?oP0yGrcSe%Yn}Hau?hUW_np3?3n509&7Jp^`irk6EfP0*QM9)?H2OfLh$j+tH_1v_SX=>$7wdU*`&nCWFL zuw$l|pMo7Ty{rv(%=Gd&*fG<~I$+04FHe9SGrg<}cFgqhGq7W(m-WDonO>d*J7#)W zAMBXvWhB@!)7=JO$4qx!U`J!H$_;HUm3mdU+1)nCWG6uw$l|=fRGdUbX-`W_o!6?3n3gOR!_6mlwf~ znO?R6J7#)$3GA5ZZfme()&(zv9i8cA8`u_frk7XXRWQ@bc3{U$FRy_eGrepNcFgqh zI@mGO%MM`2OfPSM9W%Y`2zJc$@+R0Z)5}g^$4oE306S)S*%|DZ>E)MT$4oD~fE_cv zyajg5^s+11G1JShz>b;jb^|+Rx_cY!><)uq57-la4e!9a@E(}$1$%=Xv-jZxuw%9l z3b+-{tkA`^l~iNG1JRGz>b+- zjsrVpdif{VG1JTOV8=`^UxOVpy_^7c%=GdN*fG=HiD1V}ci(~?o$2KyI2m-Nm+#sd?R_)-7L9N&Nq@*-OZXi<$ObDdTE;tX3x&XpKqj0cXQ-UIp0WL zbvI}3l=F?`Rd;jcPC4I5UUfHj?v(S5t8yZzQj}TR3;h z`9|`pyG3%RoNpwrx?41N%K1j}s=LKTbo{Dd!u>tL|3HopQdByy|Y{+$rZ9$*b;G$(?e(p|j*}RmRou_b@yP&%jev zrn}XFPC4I5UUj!d?v(S5I8_BEg*2TccKDd!tH(@WcWuzq$n{(K{4y4xUk%K1j}s=Ezyr<`vj zue#ePcgp!j@~XRybElkdB(J*LBzMaBM)In=O>?K5ZzQj}+bnm=`9|`pyUlZ_oNpwr zy4xam%K1jymfXJ;Yz^Ds?;Us_K7ebxWIN?fIp1hA z+ZlGropQdBvR%n`%bjw*(Pp+g49cByzLBy$$o9;ga=y`KwioOTbjtZg%Jv}}oIB-w zLubj|zKr{2XXDQ|Ql`88flfK!NM3a}BzMaBM)In=19GSJBi8W={2l(reP+5l5a^Wi zjg;x`pxi0v8_BEg4$hr&zLC7@ZfNe5^Nr+Hcf)e0oNpwrx;rFy%K3)Q^wM@H9LDy4 z1%H82@LiSZ?(p0x=NrkZ?uO@1Ip0WLb$3MWl=F?`Rd+|`PC4I5UUhd=?v(S5u7HG zhkjtDyOV%Unxl>3_9S4(Om`;(ois-q%k4?Qj+yRG0Xk`pHjdkqfgLm5oeFf)9Icz% zlY<>I-JJ$>(p;p$?J2;HneI*pJ7z5~CD<|3-5FrVtQDpLJ7&5&6YR7DD@oVI(=g&p zFKuVR*$4oENfgLm5oeOr%^fEozG1J|7V8=`^Gk_g4-JK71%=9uN z*fG=H1z^WaFEfE1Gre30cFgoLGuScH%SB+vOfR#59W%XL40gRtM`wDO z4Q2eaW_pnO?30J7#*B8|;|r%ybbuW*y<7u!%=9u3*fG<~wP43gFY|&OGrjx(?3n3pKCok^yX(M? z&h#=rEC4#w%k}U>Fw@I|V8=`^H-H^8y(|QF%=9t>?3n3gVX$MSmm9&3nO+tFJ7#*h z3GA5ZWl^wWrk9(+j+tH-13PAVxdrT)>1A=SW2To|!H(IyumspK)5~pO$4qxif*mv6 z-41qirkAB)Y0#No?tmYG*@CbP*fG<~onXgIFUx`*GrimecFgp$9N00_%iUndOfSoW z9kZL@9$7wdU*`&nCWFLuw$l|pMo7Ty{rv(%=Gd&*fG<~I$+04FHe9S zGu^EVcFc75Gq9sGy{rf8gU%=FR)cFgp$A=oj~%Tr*-OfMUO9W%W= z4R*}*vN6~()5|kp$4oDqfE_cvJPUTr)`CsJj+tJ54tC7+vKiPh)5~*U$4oDqgB>%y zJP&rvbhicAG1J`(U}sC%3buxA;6-=|UWQk|Y+Kk4?3ld@uYn!2?O_M7WA-|{0d~xG zgq^^S*_-eSuw%9}>;iVoehF`Z9kX3wH?U*&D|j31nC%XOz>e9k;T^DJwg>D9cFf*| z_rQ+XUa&XVF?%0A06S*;z+kXr_8a&u*wL9@_J#dGXL|V%J_0kn><@O#^zt#-G1JQs zuw$l|Pr#0uUJd{|W_tM)?3n50K(J$`m*0UMGrb%HcFgqhd$41smxIBMnO;5vJ7#(r z3Ub;j4gou6y88pz(V1Qjg~LE+difIm2xfXY9PF6sE+L0$4oCrf*mux`~~co>E$S}W2TqCf*mux91V8N^zt{b zW2Tp5z>b+-{tkA`^l~iNG1JRGz>b;jjsrVpy89>C(V1S3hZ8_&diff@0W-av2zJc$ z@-5gg)5}R<$4oEZfgLlwoD6o%^zuE}G1JQ_V8=`^qri@tUQPu&W_lS7cFgp08rU(@ z%NVd@rkB&fj+tJ@f*muxoB?*s^fC_YnCb3Juw$mXZm^>>y_^MSgU;~M-hdWvo9X2o zu%k1)w72F?ula_~@Y3FvJH6%`I>SplPhclG)64mAUe5H=-Y<80%{O$0m-ha-(`&w= zGrY7F`=q(kYrdf~ytGf2JH6%`I>TN2tL~=GopQdByy|Y6+$rZ9$*bTagoDd!tH(@XozFiUne{(K{4x|=n3%K1j}s=L{8r<`vj zuezH(cgp!j@~XQza;KbcB(J)gGk41QM)In=xpJqRZzQj}n>%;P`9|`py8*dV&Nq@* z-F4(nIp5G(ayJj-y!g8wZh$-Cjw;jLd_bq1ZzQj}n?HBjM$owt?uNU#&rEj<0G)Ea zkuu#am^f_q)c~9=1w`^NM3cfRPL1XjpS8#OXp5G-$-6{w@mJo^Nr+HcgyBZ zIp0WLb+=sZl=F?`Rd>thPC4I5UUj!Z?v(S5tL|3KopQdByy|YX+$rZ9$*b;G&z*9&^PC4I5+2&+hRzJLP;MW!sbOkUQmkLubj|j*L5HXUW~i-2Sl2bhk6mN!c5Wp8%cYRd>7OPC4I5 zUUj!??v(S5wZzQj}+aq_%`9|`pyFGKKoNpwr zy4x#v+CIu&(sl74xesT0Y2O?6sWRRD3BIf{-3`v2a=wwg>TciMDd!u>tM2y8opQdB zyy|ZM+$rZ9$*b;$^Nr+HcL(QAIp5G( zayOK5Saz1&ear2yt4w!?0G*V5#rPf2NnUk#XzrBrjpS8#hviN=-$-6{cX;lU^Nr+H zcf)h1oNpwrx;r9w%K1j}s=Fg|r<`vjuev)bcgp!j@~XR|bElkd=u9u|$H1}K+4%E~ zldN&h#=f z%mO;o%SCW8nCWF!uw$l|OTdnq?q&l!W_r04?3n3gcCcfnm&?G8nO^1qJ7#*h9PF6s zWlpeTrk5+gj+tKO0y}1Uxf1M{>1A%PW2Tp@z>b+-27nziy<81;%yicQcFc5l4cO6{ zUgm*$L1%io7JdL`dYKRGnCay@uw$l|`N58vUaki_W_now?3n50hhWD{FAIVlGrimZ zcFgp$5ZE!(%LuSzrk91mj+tI=1UqJWSp@8u>E$M{W2Toy!H$_;ZU#GMx?2qFnCb2o zu%okiU~yOibf%YE;WjYS%aUNnOfR>C9W%Ww1$NBzatGKk)63Fe$4oCj0y}1USqAKw z>E%waW2Tp7!H$_;?gBezdRY$anB53>gB>%yEDv_f^l}f_G1JQmV8=`^_ktZW-K_|A z%yjo-u%k1*tOP5A&h&C0`~=MOvI^KS)64x}$4oD)f*muxJOFmg^s*Y*G1JR~V8=`^ ztAiagy*vbV%=EGb*fG<~!(hiuFKdDwGrc?lcFgoL5bT)grk75zW2TqKz>b;j z)&e_by89{ESsT`Ybzwbt9G-xm!INONK5PJX%tk^N*fHA>HUc|lPr=h*$82NR1nihS z1J8mTvrS<$uw(Xfcn<8CZ4O(29kb`*1+ZhbC2R$D%wB|-z>eA0unpKTdl_B1A)&2Xv;F58yXorkBBB$4oE31v_SX*%$1X z>E%PPW2Tq=z>b+-J_0*tdf6ZBnCazXuw$l|Az;T$FQ0%NGrb%DcFgqhDcCX7%Yk6W zOfSC!J7#)02<({Y<@aF6Om_!^9W&j126l9|Hw=YgpfkOE4qt$oUJd~}W_tMp*fG<~ zpd!Jaypy=I@8No7zbv0ITP%d z>7^U&nCay#uw#aoeho%DW_mdr?C1E%+eV}`qaQ{+x^rkBfdCuO>uGIz@PM)In=sdA^B zZzQj}n>u&O`9|`pyJ>Q#oNpwrx|=q4%K1j}s=MiOr<`vjuezH)cgp#O&XT(s7-!@) z>r5}#!L?PUyP1GaIp0WLbvJYFl=F?`Rd=)GPC4I5UUfHX?v(S50T}SSe^Nr+H zck|>S6XVB<(Hzy^G)z)>6KSqdYM&M_!p~a zG+KBhpLq>F*Bb4>=Tn2vpa!2a4L*k(e5N$`%xLhr(BRl_aC|p7h8rBG4UWYI$6JGA zs=;y7;P`29d^9))8tmr{_T>iqWP^RO!G70Z-)gY`G&sH*>?@59;Qem!o;G;j8oW0R z{5Npiz+(ef4fc};`%;5_ron#E;Jt6~{x*0I8@x{qU!1QWr;X)d1y~VQf|X$vSQS=- z)nN@-69z&jtOaYsIcsJ z>I4ufD1*c0}Gyb72o8p!FbocXL*Xzu9EQUYa3mZB zN5e62EF1^N!wGOAoCGJsDR3&B2B*Ura3-7uXTv#gE}RGF!v%05Tm%=xC2%QR2A9JX za3x#?SHm@ME&KqkgX`gka085h8{sCn8E%1F;WoG(?tmY`op2Z24fnvk@ME|SeggNy z1Mna`1P{X_@F+Y6KZVEP3HTX22_vBko`R?08F&_c4$r~!@B+LDFTu<33cL!h!Rzn_ zya~U6U&34PD|j1z4e!9a@E*JmAHZ+mx9}l+1RujE@G1Nbeh;6)=kNvm0ltJk!k^$P z_%r+k{tADCzr#P^pYSz&1K+}T@I8!z(J%(a!Z_&8hki>7V=J^lJM@G8FbPZwlfmRL z1xyK3!PGDfObgS&^e_X=2s6RVFbm8Iv%%~z2h0g`!Q3zaI$$1{7v_WcVF6eW7J`Lg z5m*!!gT-M9SQ3_krC}LZ7M6qMVFg$bR)Upb6<8HkgVkXTSQ7?9C#(f)!#c1otOx7E z2CyM)1RKL9uqkW?o5L2cC2R#-!#1!jYzN!J4zMHa1Uthnuq*5ayTc&Z1NMZyU~kw5 z2E)FvAM6i9-~c!f4uXSWC=7!`;7~XW4u|1z1RM!R!O?IG91F+6@o)m12q(eGa0;9X zr@`rP2Am0J!P#&QoD1i{`EUVT2p7S{a0y%rm%-(51zZVN!PRgLTnj&d>)?9$A>05X z;6}I!ZiZXnR=5prhdba$a3|aacf&n!FZ>wpgP*|t@Blmr55dFm2s{dp!B63FcmjR~ zPr^v(f~VkVcm|$@pTl$TJiGue!b|WnyaKPnYw$X}0dK-D;Fs_g{0iQNU&A}_F1!ct z!w2vi_$_=0AHm1)3498_gWtnv@Hu<|e}FIHkMJk>3jPd#fxp7v;P3Dc_$Pb~-@v!< z9efX?U^I+@u`mw0^J%}eg|QXdpdI=_f0zU&g~?!Ym;$DRsbFfD2BwASV0xGVW`vnw zW|#$Lh1pg7xnOP>039$7%nS3u{ICEl2n)f&um~&)i^1Zs1S|;ZeiUa&Xp1A}2-*bnxHA#ea32nWHzFcgNt zA#f-h28Y9NI0BA@qu^*b29AZ};CMIzPK1--WH<#*h11}4I0Mdvv*2tv2hN4_;C#3M zE`*EVVz>k@h0EY_xB{+(tKe$52Cju4z;$pv{19${5pW~i1UJJia4Xyfx5FLqBe)ap zg1g}!xEFp5_rXu#es};Lgoof^cmy7W$Ka>%I6MJAgC}7mbiq^bG&}>(!q4G3cphGW z7vUv%8D4=`;Wc<2-hemZ7w}7X3w{M}!>{2Tco*J-_u&Kh4g3~9gpc53_yj(M-@)(U zGx!|7fIq;O@JIL)d0o-80cM1mU}l&FW`)^cc9;X^gt=gD7yunG z56lbm!ThiQEC>t1!mtP|3X8$wummg#OTp5x3@i)F!Sb*ItOzT?%CHKo3ai2Dum-FN z1ECYvg0*2CSQplV^;r>gU)T@!haqqP90&)&!7vnt!69%c90rHOa5w^vgrneSI0lY| zBy_=3@H9LF&%)2)Id~pkfEVEmk>0t(#5oUs!VHTJbW`o&b4ww_>g1KP;bih0?FU$w?!ve4% zECdU~BCsed28+WIup}%6OT#j-EG!4h!wRq>tOP5=DzGZ72CKswuqF(IPFM@phIL?F zSP#~R4PZmq2sVaIU{lx(His=>OV|pwhHYS5*bcUb9biY;33i5EU{}};c85W*2kZ%Z z!QQYB42FGSKiD6JzyWX|90Ui$P#6Y>z@cy$91g?b2sjdsf}`OWI2MkBRPd+zhwCt#BLM4tKzh;7+&;?uL8dUidNG2S0)P;Q@FM9)gGA5qK0HgP+3V z@C5t}o`jLm1y8}#@C-Z)KZobwd3XU{gqPrDcm-aC*Wh({1KxyRz%St~_!YbjzlL|< zU3d@PhY#R4@LTv0K7xbz^1SnYz|w%mar9U z4coxBupMjx6YLDTz^<)uq57-m-g1uoM7!3Quey~3bfdk+`I0z1gp)d>% zfkWXiI2?w<5pX0N1xLd%a4Z}L$HNJ5BAf&#!zpkooCc@E8E__?1!u!Ka4wt&=feeX zAzTC(!zFMjTn3lJ6>ue71y{p0a4q}*u7m60hj0UofE(c^xEXGNTj4gi9qxc1!JTjy z+zt1@z3^kW4}Jpo!vpXjJOmHJBk(9Z20w+z;R*N|JP9MA3!Z|f;Td=qeh$yU^Y8+^ z2rt3Q@Cv*Nufgl^2D}NsfM3E}@GE#5ehu%yyYL>o4;OB$POvlV0=vR)usaNb zJz!7R3-*S6U@+_p`@#M&1P*`$;UG8|hQcs71P+D6;BXiYN5GMA6dVo5z_D-~91kbJ ziEt8}45z@Ua2lKrXTX_o7Mu;|z`1Z9oDUbkg>VsE441&Aa2Z?;ZeiUa&Xp1A}2-*bnxHA#ea32nWHz zFcgNtA#f-h28Y9NI0BA@qu^*b29AZ};CMIzPK1--WH<#*h11}4I0Mdvv*2tv2hN4_ z;C#3ME`*EVVz>k@h0EY_xB{+(tKe$52Cju4z;$pv{19${5pW~i1UJJia4Xyfx5FLq zBe)apg1g}!xEFp5_rXu#es};Lgoof^cmy7W$Ka>%I6MJAgC}7mbiq^bG&}>(!q4G3 zcphGW7vUv%8D4=`;Wc<2-hemZ7w}7X3w{M}!>{2Tco*J-_u&Kh4g3~9gpc53_yj(M z-@)(UGx!|7Xym!{KR93GBj(I4{WyMiVVU>dYyNbd)SvryW%+BCan>>k%e%4s4Y%>v zGAYZuv-~Z$aoI8%%Y#_{j@x)`nVjW4SpJ^dIBuDORT{^}L$q z!GB@-G;YtqeJ3&Ud|q3(p3I0Z&$(V>d*1ar+jFni+Ma*C-u4{qHMi$sue&`Lr@8f? zS!NsjKE-ph^L(%CJx6mt_0Bas-#g#%oS(h}_!G;#$9@mw zHGtnoIG^!az&VZA1I}x_CU9=!b%FC6uMM2zw2olpz4aQwxsKNf>3fAgvOI{}vv42# zkKb{$?7?mJA?H(GLpZ1MI>LFC*AmXHyq<7==EUy5Z|julvmRalKRY7X;Sja*Xjl?8gB-ipTIc4&n(lPxF%; z^WZ7OzvKP!dC|YPpSXWxmg04DBKE|})=0@UW1SRxXJu=pq#yA&RP1Y({T+3DkE@-! zn3G4n$NbEiY9Hff{2@2@o$jTq1!vrnb8f|6$Xav3r*c){-8iUiex>nT12s_#wNVFk zQ4i);l3$#g#GYc=*-4UjMgL;&v7&$f?+hgwpGsS`SF?Mpv_}VYL??7c7j#u)JmWki z_A<-PRFd(JbCuZhEIV6C=1ZKf#NKGx8A~$9VjUNIsAX%pWM0PK(y@PC_V@JhJxuqW zwQ=3y-s8I>lS`}xv&dxxM>Zol@)^aE(P)mG#&BdcmLo6kPj|ud5$na+GcH>*CeKr4 zyxMUb;~e7`;}UsCj*(B?_cdxff3bdyJ?OGEWb(YmIx_aE%hr;~^BwES*wZdsQzqXB zv9650@3OUJ^8FF(%h)3?TVp2QIkC=+z4WrRX4)$6LuwoHVH|<)j{mil$}ri}rLjZ_3uRX|tM{YR#*% zwQbs@W}#YhuWXH*HmX^v*8D76>!uBATszJ-%GSJTy&7{Z&Opl6zGpFU*+EX z8u#JXxi7!L{rFAp&u{SnewzpKJFamDR%^e@7x8;|A0OZ-j=}GM!plYGKdU(BAzXnM ziYiV(4b(&})K+86iaI=$>!Kb`M13@XdoIlFqLG>v+?X%slhFjHpedT6Ia;8ln&sSz zFX7hkd#7lFwrGd;=zxxDmT@P(lslsfy5cl=E{bmOJQUs4EM?C^oS7A#g`x*!UC6kQ zZz0=4u7ym~5|(M4trapYF- zb>Vp}MqngHsdJh->gAd>LQCSMfD9^ZDyMh`)hv;#>GOzJu@Ld-%SZdHe&8`~M++glq9*T!-s% z18!6^mv7?8=w{r4TX7qHg4^*^L}qjN4vrk}#9g=>KgTcdOZ*DIRx?}P8@WDj!e(s2 zRy>SHunpVQ`~#h}jbZT{?Vj4kvG^?;$KrQ9lz-30vG@ZU$KsD{42wT;jNzX-#_=z% zq0i|-?Z0v_{u};|d*Jhm37CjUn5R6aT`$)y(Aoa38)O4`2r3 z=W(6*T=!qdE5*-coceN%ON?WTd(4xVSBtP1OVq@5mU2H{hUHj+l~{$LLe_<>3t1PkPWoENI}PA%*p408iCx%@J=m+p{T2IoAn(TkJc`HgI1b_o zJgG(=nf|2{Ssycf%=9tS$4nnH{mb+%)3-DMU9^XJ?^D`cwTEjz&6n{rJc1ALNIuM? z_y~{YXL$@i$7A_<*SH-qE-&zD{32e$%XkH^;x)XEH`Lt6Z}RE<7T(4?co*;CeSCnU zYNqos?#3UcRMo6z$WcK7KdheM1@oj@gU{fasD;|7gSw~(8CC16nZ^ycJD-GxXoSW% z8BK5snyQ(~&G<}ijuvQ%R%ne=(FSeROyPEX7Pm(SbVMg~hU}{{ugWP+W*MX&?Eb3J z@90l_PP{*U-V?plOp^Bkb`RM?n6sHVo0+qjIh&cU*)o``X(D8lWSG1EB&Q_HTno!7Wf{eQ<#Nj9l*=h^gSL=S-X8Kvx3G+oG0f$Z%PE&rE~i{hxr}o83vy75) z&gGQLDVI|&r(8~XZ}dT5H8*iTHW%{#7=VE|8y~?rI2Y%sxseC4Ig_7{3veL@<01^f zP+Y9$2EK%4n_r5L;$s+w;kXPVFjCF+Jc?zXkH#2`#pSpH<8UReQgfZWXRtBIW?~i| P#B9vLT+G9KHP`+Z(-MX5QbNYKP1I*?1Ya12`Ge*gb&M-m4x!e2?;53;{X(gR4uZ+j<%3?g|xzo zPdV^|_y^|5fdfa5rivfHkt6T3v#mJm+9a+5Dpj7UcY9{~-R_>A9@VUF-?@9jbFC_I zN+|b&$dB6n!1GFIDoNN6{817N7`Zkdjr;yq*k^Q7s4$M=tzhKz*Coo#1-$^$OG4F+ z{mtHYDGB3?*<2Vg-AL!*p2_U_gJ38dsiLg;)W-L1e;9O9E;eSiH9VH7i9c;OOdyeF)3s%yi1!+!lZ{o zHTLR8XS^9ESL)rkAJhk$P=M{kvhW_jS8BcAEnS|1|ZjTiMsR%6g{OoIUuo=d@=jsP5~Xawm?W zpp%ea?66E8;6UBXe*W1vpNKc}c>xYF;rBJ+b%B`h`-X5qK!@Lp0vx;NzNtPObohNsfJ039tqE|@;rEgN z2OWOj7HR@IV!JHBK}T%w2ypWGdRKGr3GWM6gnYixT@^kM;E>CM`9k-h@R0x~pD%PD z3o8Pge7?}tg;fDgK3~LyU%m~RLOx&U@Vh3!$>$3lep>>Ze7?}(_nH7FpD%R8#(x0U zZh!UZy)V@%Z2CxTtTc$DpgsPUEg(Gcfm^^SY5r+v!0O-X-s;`XGx%5;#xPA@$BaP`}h2^-kW2pM6^th?~0SoKaiUb4_3$sdws$ nI%nTlFKcCeoKfnG^T&B$PdR_=73Y-NvV2if=Fns3$_Zh#p063hAeJ&Lv~{Vt=b1w z+lW-8Xqk>>pgv}_)WURHJFT1Q)X|wRGG?ZLolamnZN>TrjCAU#P5S-(_IG~9*fjm0 z{$@V8_uO;tx#ym9?{9vam8({d@ueTXnsFjG8jgix^{ru_F9S_pGSV6jC6eJbr1X;$ z>u3!%MOu-@iag|4yeXVW>uc2HG8T@8F*+>#X*tvoz1uOu`2=??5@WdB8OOH3^%rUj zx66m4$ke#wV^@8sJ>2Nn;v=&;-qB8rp4bz~_DF1pqo1w6<3)}bSKk&1-4kw4MB*`& zSt4VFZ_fp_IU|f0HAmz1NqGyk#Urs~s59Q)6iP()C^yb}te-3Lnj=x=D6JTxLqDap zzAMs7ew5RY0^E^o(XcnaSpD&GzNAW#*C#tjA$wsiPgXLLjIy8}YPF=Xqal(kDQk(h zhRfO#E$xxcNFrR83?~z%;jZ|Cn&QRPW#v_66^kn>YKki=YbvX1Dl44nPJQ#$^*c4~ zY>B&$SZQ@dC|j`jadBw)S`%g)>{~KtD5mww%3; zUAIaj>A5z)S8{e^T$%22bj;z5*`wGpPN#)T+O+L%XXL-wW*_lv_l#zn%TZ?!gi1Zg zyY1KZ;Sa5jM59~kqaEquV8i-6>RZ}5aeBT9LACz$ol?%z7nsGUe;ztv4{PmB z3)9Z4&2`pkttVe>*sY;@PJc{`)AJQIsU7dwa~ z^oR8`Qiu1yTXJTuFLkcxokg2ZtWI4*@AIX4U;CipXHOIb9Q|xxz{a2Ypy$!S=kh%} z8}EPc$Uk5H^Ul=aUrc#*Q1*TG4T;*jUhw9fdTHuUt9CU7PUiIXoq3@rb?$KYGnMz% zruw(^6@K)Wr_5i#nLo>Ca_U%W^49kj{?+7kbUt6;Y)<8YzMu4)fAj0fZ*@&8sO!3D zvdtf@_5AUU3t=5&pjuW zzrfk!z#l%G_gvvci_`SO<)G_t zc&oVHdA@XTuZ>%LCR;n&c#rGfS}$L=zxH##8vJe0=k3>#mkyJqcgHPHo(J^G9d%?o|}W6$1A#vY%nfxp>;xQU1o>btIbpu#iu|IVC(vVu9(lgF-j(r7pE2~v#6}FPhdpymk6a=a z2>jmQ_8Bvowdu1ye(*^k~ z$c$s3AumcN`n^$eep*{FMoO;%~7i7x+ZB^~bt0 ze`u%jxl7`B^F|2edRfu6jxfAU6* z+$Yu~@5ok1;>ABY)+0`2?6EBt*b^T*ix+!jt2g%CRo)Bi89Uan__0UFIC_i6@@M{7 z)ADDrB9nX00X1QN)R%nQIJM_2FB9a+Wcr+G?gBX?kK~a2)ClrJo>-rCEFaWq=ounz z_Gft`Hf)&7JaqUn{e1O_fpbP(?7XoiHL^HuKiCp0>626JEuXgkIT|A;)OC`WE2aqYMm{Xh_+WmOAQ!~U zn%FYW@`oMkTOP`DP56b%}*M!IwAH;;7dQ!(MVfALrWc<^& z^I_+Vm{<=zI{YEqe(B?%dQeL{59rASW7wD<;>VtRkOS&UKFAHZw7ilV#>g>x-Y48? z{P9jz@|nOcx*-T|e8x>0ZP!#Sgt)P)*S|7yY6;XYCKo5aoH7O`9eL{O{{wPK}MCF;a# zu|}*F>%@9-tN6UwAU2BI#23UD#qHt_K{+>zE#gb!%i=3ytGH8qRcsS?iICVX>P3TS z6ip&5n#B&$A|hg^xLZU;tB8rXXcN0cyGV$n=n(gaPSGWHi#_6N;_KpGai6$fd_#Ow zd`o;=>=oY;4~PfFL*l#Qd!k!BEcS^<#D4L8@u)Z;9uo(}55ytyxOhVRP#hK+VvHCo zGDVii7C9nUF+miFiDHtNET)J;F;!e6rimgkUCa>wKWA6>IpdyjhqzF} A;s5{u literal 0 HcmV?d00001 diff --git a/tests/unit_tests/mesh_to_vtk/umesh.vtk b/tests/unit_tests/mesh_to_vtk/umesh.vtk new file mode 100644 index 000000000..90f226121 --- /dev/null +++ b/tests/unit_tests/mesh_to_vtk/umesh.vtk @@ -0,0 +1,292 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 58 double +-0.02593964576 -1 -1.1195739536 -0.40239958217 -0.40962166746 -1.9256035383 -1 0.02593964576 -1.1195739536 +-0.2248833639 0.23189144433 -1.2841109811 0.02593964576 1 -1.1195739536 1 -0.02593964576 -1.1195739536 +-0.042928712675 0.066425810853 -0.54270728017 -1 1 -1.5 -1 1 -0.5 +-0.49265312381 0.49332495053 -1.9656359509 -0.02608137991 -1 1.121521147 0.30281888252 -0.30293622513 1.8067963894 +1 -0.02608137991 1.121521147 -0.23236342389 0.26314805583 1.2867556783 0.074337770194 0.074337770194 2.5 +-0.40005770996 -0.36909660956 1.8768871762 -0.4896743493 0.49451064621 1.952829049 -1 1 1.5 +0.02608137991 1 1.121521147 -1 1 0.5 -0.038567136385 0.14628887118 0.55111724874 +-1 0.02608137991 1.121521147 0 -1 -2.5 0.37543954218 -0.36967061971 -1.8691105401 +1 -1 -2.5 0.074337770194 -0.074337770194 -2.5 1 0 -2.5 +0.4135797166 0.38351746262 -1.9274856063 1 -1 -1.5 -1 0 -2.5 +-1 -1 -1.5 -1 -1 -2.5 -1 0 2.5 +-1 1 2.5 0 1 2.5 0.46519182209 0.49356920409 1.95286052 +0 -1 2.5 -1 -1 1.5 -1 -1 2.5 +1 -1 2.5 1 -1 1.5 1 0 2.5 +1 1 2.5 1 1 1.5 0 1 -2.5 +-1 1 -2.5 1 -1 -0.5 1 1 -0.5 +0.038669824604 1 0.0046021677516 -1 -1 0.5 -0.038669824604 -1 0.0046021677516 +1 1 0.5 -1 0.038669824604 0.0046021677516 1 -0.038669824604 0.0046021677516 +1 1 -2.5 1 1 -1.5 1 -1 0.5 +-1 -1 -0.5 +CELLS 155 616 +OFFSETS vtktypeint64 +0 4 8 12 16 20 24 28 32 +36 40 44 48 52 56 60 64 68 +72 76 80 84 88 92 96 100 104 +108 112 116 120 124 128 132 136 140 +144 148 152 156 160 164 168 172 176 +180 184 188 192 196 200 204 208 212 +216 220 224 228 232 236 240 244 248 +252 256 260 264 268 272 276 280 284 +288 292 296 300 304 308 312 316 320 +324 328 332 336 340 344 348 352 356 +360 364 368 372 376 380 384 388 392 +396 400 404 408 412 416 420 424 428 +432 436 440 444 448 452 456 460 464 +468 472 476 480 484 488 492 496 500 +504 508 512 516 520 524 528 532 536 +540 544 548 552 556 560 564 568 572 +576 580 584 588 592 596 600 604 608 +612 616 +CONNECTIVITY vtktypeint64 +0 1 2 3 4 3 5 6 7 +4 8 3 7 3 2 9 8 3 +4 6 10 11 12 13 14 13 15 +16 17 13 18 16 19 18 13 20 +10 13 21 15 22 1 0 23 22 +24 25 23 25 1 3 9 26 23 +5 27 24 28 26 23 0 1 3 +23 29 30 31 1 22 31 30 1 +2 3 1 9 29 2 1 9 22 +0 28 23 32 15 21 16 14 13 +11 15 33 17 34 16 33 34 32 +16 34 16 18 35 10 11 13 15 +18 13 12 35 36 11 10 15 36 +37 38 15 32 38 37 15 32 21 +17 16 36 10 37 15 36 39 40 +11 41 40 39 11 41 12 11 35 +17 21 13 16 10 13 12 20 12 +13 11 35 14 13 16 35 34 17 +18 16 42 34 43 35 41 12 40 +11 14 41 11 35 42 41 14 35 +0 3 5 23 0 5 3 6 44 +4 9 27 29 2 30 1 44 45 +29 9 29 45 7 9 0 46 5 +6 47 48 4 6 0 3 2 6 +4 5 3 27 49 10 50 20 8 +2 3 6 51 18 48 20 49 50 +52 20 19 13 21 20 51 53 12 +20 19 48 18 20 25 3 1 23 +7 4 3 9 25 3 23 27 44 +4 7 9 25 23 26 27 54 55 +44 27 54 26 55 27 18 12 13 +20 14 16 34 35 19 21 52 20 +30 2 0 1 22 25 29 1 5 +23 3 27 25 26 44 27 22 29 +31 1 54 44 26 27 25 44 9 +27 22 28 24 23 26 5 55 27 +25 24 26 23 22 25 1 23 25 +44 29 9 25 29 1 9 48 6 +53 20 8 4 48 6 10 40 56 +12 10 40 12 11 17 19 18 13 +30 57 0 2 25 9 3 27 7 +8 2 3 26 28 5 23 29 7 +2 9 28 0 5 23 37 10 49 +21 49 50 57 52 56 12 53 20 +10 12 56 20 43 18 12 35 32 +37 21 15 14 11 13 35 19 52 +8 48 19 52 48 20 41 39 36 +11 50 6 52 20 46 53 5 6 +4 3 9 27 44 7 45 9 55 +47 4 5 28 0 46 5 8 52 +2 6 57 2 52 6 50 56 53 +20 10 56 50 20 57 50 0 6 +37 10 21 15 33 32 17 16 34 +18 43 35 57 52 50 6 51 48 +47 53 18 16 13 35 50 56 46 +53 14 36 32 15 14 11 36 15 +17 21 19 13 10 21 13 20 51 +12 18 20 50 53 46 6 14 15 +32 16 42 43 41 35 14 41 36 +11 49 52 21 20 36 40 10 11 +47 53 48 6 49 21 10 20 21 +15 13 16 51 48 53 20 55 5 +4 27 44 55 4 27 47 4 5 +6 8 48 52 6 47 5 53 6 +50 53 6 20 57 0 2 6 50 +46 0 6 52 6 48 20 42 14 +34 35 43 18 51 12 41 43 12 +35 22 30 0 1 14 32 34 16 +36 38 32 15 +CELL_TYPES 154 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 + +CELL_DATA 154 +FIELD FieldData 1 +ids 1 154 double +0 1 2 3 4 5 6 7 8 +9 10 11 12 13 14 15 16 17 +18 19 20 21 22 23 24 25 26 +27 28 29 30 31 32 33 34 35 +36 37 38 39 40 41 42 43 44 +45 46 47 48 49 50 51 52 53 +54 55 56 57 58 59 60 61 62 +63 64 65 66 67 68 69 70 71 +72 73 74 75 76 77 78 79 80 +81 82 83 84 85 86 87 88 89 +90 91 92 93 94 95 96 97 98 +99 100 101 102 103 104 105 106 107 +108 109 110 111 112 113 114 115 116 +117 118 119 120 121 122 123 124 125 +126 127 128 129 130 131 132 133 134 +135 136 137 138 139 140 141 142 143 +144 145 146 147 148 149 150 151 152 +153 From 46e98156ac105926d3c2f0b158fc12da51f3e807 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 18 Jul 2022 16:49:44 -0500 Subject: [PATCH 46/83] Adding notice about MOAB meshes --- src/state_point.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/state_point.cpp b/src/state_point.cpp index 9ea045915..a2d8b8881 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -812,12 +812,18 @@ void write_unstructured_mesh_results() auto umesh = dynamic_cast(model::meshes[mesh_idx].get()); + if (!umesh) continue; if (!umesh->output_) continue; + if (umesh->library() == "moab" && !mpi::master) { + warning(fmt::format("Output for a MOAB mesh (mesh {}) was requested but will not be written. Please use the Python API to generated the desired VTK tetrahedral mesh.", umesh->id_)); + continue; + } + // if this tally has more than one filter, print // warning and skip writing the mesh if (tally->filters().size() > 1) { @@ -889,12 +895,10 @@ void write_unstructured_mesh_results() std::string filename = fmt::format("tally_{0}.{1:0{2}}", tally->id_, simulation::current_batch, batch_width); - if (umesh->library() == "moab" && !mpi::master) - continue; - // Write the unstructured mesh and data to file umesh->write(filename); + // remove score data added for this mesh write umesh->remove_scores(); } } From 2f8d2fe003a7713672e8b61be16da453ddab882a Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 20 Jul 2022 09:33:14 -0500 Subject: [PATCH 47/83] Empty commit attempting to fix CI. From 782c2411906041c8db1c01bb65b3ad1c6692d1af Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 21 Jul 2022 05:38:58 -0500 Subject: [PATCH 48/83] Always continue, but only write warning on mpi master rank --- src/state_point.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/state_point.cpp b/src/state_point.cpp index a2d8b8881..b4200ffe2 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -819,7 +819,8 @@ void write_unstructured_mesh_results() if (!umesh->output_) continue; - if (umesh->library() == "moab" && !mpi::master) { + if (umesh->library() == "moab") { + if (mpi::master) warning(fmt::format("Output for a MOAB mesh (mesh {}) was requested but will not be written. Please use the Python API to generated the desired VTK tetrahedral mesh.", umesh->id_)); continue; } From f4bb731687acfe003080e9614c117795d69c71aa Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 21 Jul 2022 22:29:51 -0500 Subject: [PATCH 49/83] Keep outgoing datasets in a list until VTK file is written --- openmc/mesh.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 601d0ac94..d88f6e363 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1722,6 +1722,7 @@ class UnstructuredMesh(MeshBase): grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds()) # check that datasets are the correct size + datasets_out = [] if datasets is not None: for name, data in datasets.items(): if data.shape != (self.dimension,): @@ -1739,6 +1740,7 @@ class UnstructuredMesh(MeshBase): # add data to the mesh for name, data in datasets.items(): + datasets_out.append(data) arr = vtk.vtkDoubleArray() arr.SetName(name) arr.SetNumberOfTuples(data.size) From 0118067e2f7cc1fd37f50591f998b2f97ce897ae Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Jul 2022 16:00:50 -0500 Subject: [PATCH 50/83] Fix for IncoherentElasticAE.from_hdf5 --- openmc/data/thermal_angle_energy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index 1f37863ea..17a560092 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -122,7 +122,7 @@ class IncoherentElasticAE(AngleEnergy): Incoherent elastic distribution """ - return cls(group['debye_waller']) + return cls(group['debye_waller'][()]) class IncoherentElasticAEDiscrete(AngleEnergy): From 91000f3ea0dac86e4129c22b75937bec85f64e75 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Fri, 1 Jul 2022 16:51:39 -0500 Subject: [PATCH 51/83] refactored mesh.cpp to allow mesh pointer or filename and reduced copied code --- include/openmc/mesh.h | 9 ++++++--- src/mesh.cpp | 27 ++++++++++++++++++++------- 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 21c58862d..4d5ebc46e 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -53,7 +53,7 @@ extern vector> meshes; #ifdef LIBMESH namespace settings { -// used when creating new libMesh::Mesh instances +// used when creating new libMesh::MeshBase instances extern unique_ptr libmesh_init; extern const libMesh::Parallel::Communicator* libmesh_comm; } // namespace settings @@ -671,7 +671,8 @@ class LibMesh : public UnstructuredMesh { public: // Constructors LibMesh(pugi::xml_node node); - LibMesh(const std::string& filename, double length_multiplier = 1.0); + LibMesh(const std::string & filename, double length_multiplier = 1.0); + LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier = 1.0); static const std::string mesh_lib_type; @@ -705,6 +706,7 @@ public: private: void initialize() override; + void set_mesh_pointer_from_filename(const std::string& filename); // Methods @@ -715,7 +717,8 @@ private: int get_bin_from_element(const libMesh::Elem* elem) const; // Data members - unique_ptr m_; //!< pointer to the libMesh mesh instance + unique_ptr unique_m_; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC + libMesh::MeshBase * m_; //!< pointer to libMesh MeshBase instance vector> pl_; //!< per-thread point locators unique_ptr diff --git a/src/mesh.cpp b/src/mesh.cpp index e8f79aa37..68675b969 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2327,13 +2327,31 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) initialize(); } -LibMesh::LibMesh(const std::string& filename, double length_multiplier) +// create the mesh from a pointer to a libMesh Mesh +LibMesh::LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier) { - filename_ = filename; + m_ = &input_mesh; set_length_multiplier(length_multiplier); initialize(); } +// create the mesh from an input file +LibMesh::LibMesh(const std::string& filename, double length_multiplier) +{ + set_mesh_pointer_from_filename(filename); + set_length_multiplier(length_multiplier); + initialize(); +} + +void LibMesh::set_mesh_pointer_from_filename(const std::string& filename) +{ + filename_ = filename; + unique_m_ = make_unique(*settings::libmesh_comm, n_dimension_); + m_ = unique_m_.get(); + m_->read(filename_); +} + +// intialize from mesh file void LibMesh::initialize() { if (!settings::libmesh_comm) { @@ -2344,15 +2362,10 @@ void LibMesh::initialize() // assuming that unstructured meshes used in OpenMC are 3D n_dimension_ = 3; - m_ = make_unique(*settings::libmesh_comm, n_dimension_); - m_->read(filename_); - if (specified_length_multiplier_) { libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); } - m_->prepare_for_use(); - // ensure that the loaded mesh is 3 dimensional if (m_->mesh_dimension() != n_dimension_) { fatal_error(fmt::format("Mesh file {} specified for use in an unstructured " From 23646c3c7cf8beec9de07a02ea8a9caa38bf7086 Mon Sep 17 00:00:00 2001 From: David Andrs Date: Sat, 23 Jul 2022 08:56:12 -0600 Subject: [PATCH 52/83] Allow OpenMC to compile against fmt v9 --- include/openmc/position.h | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/include/openmc/position.h b/include/openmc/position.h index 5ab0774f4..4fde420cc 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -5,6 +5,7 @@ #include #include // for out_of_range +#include "fmt/format.h" #include "openmc/array.h" #include "openmc/vector.h" @@ -210,4 +211,14 @@ using Direction = Position; } // namespace openmc +template<> +struct fmt::formatter : formatter { + template + auto format(const openmc::Position& pos, FormatContext& ctx) + { + return fmt::formatter::format( + fmt::format("({}, {}, {})", pos.x, pos.y, pos.z), ctx); + } +}; + #endif // OPENMC_POSITION_H From 7fc9d20a3011d6aec405f6c22669da365502de1e Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 12:20:19 -0500 Subject: [PATCH 53/83] added back accidentally deleted line --- src/mesh.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh.cpp b/src/mesh.cpp index 68675b969..e4ede6a04 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2366,6 +2366,8 @@ void LibMesh::initialize() libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); } + m_->prepare_for_use() + // ensure that the loaded mesh is 3 dimensional if (m_->mesh_dimension() != n_dimension_) { fatal_error(fmt::format("Mesh file {} specified for use in an unstructured " From deda82d6319e628d4ca549e741730a6c1821386e Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 12:20:47 -0500 Subject: [PATCH 54/83] added clarification to comment --- include/openmc/mesh.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 4d5ebc46e..ce512c767 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -718,7 +718,7 @@ private: // Data members unique_ptr unique_m_; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC - libMesh::MeshBase * m_; //!< pointer to libMesh MeshBase instance + libMesh::MeshBase * m_; //!< pointer to libMesh MeshBase instance, always set during intialization vector> pl_; //!< per-thread point locators unique_ptr From 6a9be64d7112ebb65cb671076385c564c9462747 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 12:22:59 -0500 Subject: [PATCH 55/83] left off semicolon --- src/mesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index e4ede6a04..44580d81b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2366,7 +2366,7 @@ void LibMesh::initialize() libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); } - m_->prepare_for_use() + m_->prepare_for_use(); // ensure that the loaded mesh is 3 dimensional if (m_->mesh_dimension() != n_dimension_) { From fd8a8207082522a7a454931c49453ee5f0c833e9 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 12:40:13 -0500 Subject: [PATCH 56/83] added public mesh_ptr() accessor. moved m->prepare_for_use() to filename constructor since it only happens there. for consistency, added theh length_multiplier scale modification to constructors --- include/openmc/mesh.h | 2 ++ src/mesh.cpp | 18 ++++++++++++------ 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index ce512c767..cd2e56bb0 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -704,6 +704,8 @@ public: double volume(int bin) const override; + libMesh::MeshBase * mesh_ptr() const; + private: void initialize() override; void set_mesh_pointer_from_filename(const std::string& filename); diff --git a/src/mesh.cpp b/src/mesh.cpp index 44580d81b..5c2c8ea25 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2332,6 +2332,9 @@ LibMesh::LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier) { m_ = &input_mesh; set_length_multiplier(length_multiplier); + if (specified_length_multiplier_) { + libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); + } initialize(); } @@ -2340,6 +2343,10 @@ LibMesh::LibMesh(const std::string& filename, double length_multiplier) { set_mesh_pointer_from_filename(filename); set_length_multiplier(length_multiplier); + if (specified_length_multiplier_) { + libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); + } + m_->prepare_for_use(); initialize(); } @@ -2362,12 +2369,6 @@ void LibMesh::initialize() // assuming that unstructured meshes used in OpenMC are 3D n_dimension_ = 3; - if (specified_length_multiplier_) { - libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); - } - - m_->prepare_for_use(); - // ensure that the loaded mesh is 3 dimensional if (m_->mesh_dimension() != n_dimension_) { fatal_error(fmt::format("Mesh file {} specified for use in an unstructured " @@ -2565,6 +2566,11 @@ double LibMesh::volume(int bin) const return m_->elem_ref(bin).volume(); } +libMesh::MeshBase * mesh_ptr() const +{ + return m_; +} + #endif // LIBMESH //============================================================================== From 045bc34ec4353ed0ff94448198e61027a26efe66 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 12:59:43 -0500 Subject: [PATCH 57/83] moved accessor to one liner in header --- include/openmc/mesh.h | 4 ++-- src/mesh.cpp | 5 ----- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index cd2e56bb0..6d2c2f029 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -704,7 +704,7 @@ public: double volume(int bin) const override; - libMesh::MeshBase * mesh_ptr() const; + libMesh::MeshBase* mesh_ptr() const { return m_; }; private: void initialize() override; @@ -720,7 +720,7 @@ private: // Data members unique_ptr unique_m_; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC - libMesh::MeshBase * m_; //!< pointer to libMesh MeshBase instance, always set during intialization + libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set during intialization vector> pl_; //!< per-thread point locators unique_ptr diff --git a/src/mesh.cpp b/src/mesh.cpp index 5c2c8ea25..d849e5d55 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2566,11 +2566,6 @@ double LibMesh::volume(int bin) const return m_->elem_ref(bin).volume(); } -libMesh::MeshBase * mesh_ptr() const -{ - return m_; -} - #endif // LIBMESH //============================================================================== From 5161da8031ce97d2c5b19a98614d9d5fc2a52f27 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 16:39:54 -0500 Subject: [PATCH 58/83] retrireve filename and length_multiplier from XML file and call methods similar to filename case for XML constructor --- include/openmc/mesh.h | 2 +- src/mesh.cpp | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 6d2c2f029..0ace6dfea 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -7,13 +7,13 @@ #include #include "hdf5.h" -#include "pugixml.hpp" #include "xtensor/xtensor.hpp" #include "openmc/memory.h" // for unique_ptr #include "openmc/particle.h" #include "openmc/position.h" #include "openmc/vector.h" +#include "openmc/xml_interface.h" #ifdef DAGMC #include "moab/AdaptiveKDTree.hpp" diff --git a/src/mesh.cpp b/src/mesh.cpp index d849e5d55..771d6f81f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2324,6 +2324,14 @@ const std::string LibMesh::mesh_lib_type = "libmesh"; LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) { + std::string& filename = node.get_node_value(node, "filename"); + double length_multiplier = std::stod(node.get_node_value(node, "length_multiplier")); + set_mesh_pointer_from_filename(filename); + set_length_multiplier(length_multiplier); + if (specified_length_multiplier_) { + libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); + } + m_->prepare_for_use(); initialize(); } From 3c2d79a24198b87a1ed4a6473c920c02da773374 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 17:01:50 -0500 Subject: [PATCH 59/83] added back in header (despite redundancy) due to use in mesh files --- include/openmc/mesh.h | 1 + 1 file changed, 1 insertion(+) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 0ace6dfea..13138a739 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -8,6 +8,7 @@ #include "hdf5.h" #include "xtensor/xtensor.hpp" +#include "pugixml.hpp" #include "openmc/memory.h" // for unique_ptr #include "openmc/particle.h" From 5f1f0827aed488f1fde8876fcf28961c2ebd7d31 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 17:55:43 -0500 Subject: [PATCH 60/83] not supposed to call method on node object, other places use it from no object. following those examples --- src/mesh.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 771d6f81f..e2c1e50ec 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2324,8 +2324,8 @@ const std::string LibMesh::mesh_lib_type = "libmesh"; LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) { - std::string& filename = node.get_node_value(node, "filename"); - double length_multiplier = std::stod(node.get_node_value(node, "length_multiplier")); + std::string& filename = get_node_value(node, "filename"); + double length_multiplier = std::stod(get_node_value(node, "length_multiplier")); set_mesh_pointer_from_filename(filename); set_length_multiplier(length_multiplier); if (specified_length_multiplier_) { From 35c3ef9e2dbfccfe212267a02a6a72a2a1849821 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Tue, 26 Jul 2022 18:18:26 -0500 Subject: [PATCH 61/83] filenanme shouldn't be a reference since the RHS initializes with a non-const rvalue, causing an error --- src/mesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index e2c1e50ec..ca19a02b6 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2324,7 +2324,7 @@ const std::string LibMesh::mesh_lib_type = "libmesh"; LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) { - std::string& filename = get_node_value(node, "filename"); + std::string filename = get_node_value(node, "filename"); double length_multiplier = std::stod(get_node_value(node, "length_multiplier")); set_mesh_pointer_from_filename(filename); set_length_multiplier(length_multiplier); From 1e4371549c6a3478a79d064933c2f8edd451e528 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 27 Jul 2022 09:41:56 -0500 Subject: [PATCH 62/83] Returning unstructured mesh dimension as a tuple. Improving error message format --- openmc/mesh.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index d88f6e363..9c837fd8e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1616,7 +1616,7 @@ class UnstructuredMesh(MeshBase): @property def dimension(self): - return self.n_elements + return (self.n_elements,) @property def n_dimension(self): @@ -1725,15 +1725,15 @@ class UnstructuredMesh(MeshBase): datasets_out = [] if datasets is not None: for name, data in datasets.items(): - if data.shape != (self.dimension,): - raise ValueError(f'Cannot apply dataset {name} with ' + if data.shape != self.dimension: + raise ValueError(f'Cannot apply dataset "{name}" with ' f'shape {data.shape} to mesh {self.id} ' f'with dimensions {self.dimension}') if volume_normalization: for name, data in datasets.items(): if np.issubdtype(data.dtype, np.integer): - warnings.warn(f'Integer data set {name} will ' + warnings.warn(f'Integer data set "{name}" will ' 'not be volume-normalized.') continue data /= self.volumes From 3010c5e9ad27cff6c4f3a54e8d202f0ef796ef5e Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Wed, 27 Jul 2022 12:48:30 -0500 Subject: [PATCH 63/83] moved m->prepare_for_use() back to initialize, but only occurs if unique_m exists. XML constructor can access member variables from UnstucturedMesh constructor --- src/mesh.cpp | 26 +++++++++++--------------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index ca19a02b6..8e9d4bcec 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2324,14 +2324,8 @@ const std::string LibMesh::mesh_lib_type = "libmesh"; LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) { - std::string filename = get_node_value(node, "filename"); - double length_multiplier = std::stod(get_node_value(node, "length_multiplier")); - set_mesh_pointer_from_filename(filename); - set_length_multiplier(length_multiplier); - if (specified_length_multiplier_) { - libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); - } - m_->prepare_for_use(); + set_mesh_pointer_from_filename(filename_); + set_length_multiplier(length_multiplier_); initialize(); } @@ -2340,9 +2334,6 @@ LibMesh::LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier) { m_ = &input_mesh; set_length_multiplier(length_multiplier); - if (specified_length_multiplier_) { - libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); - } initialize(); } @@ -2351,10 +2342,6 @@ LibMesh::LibMesh(const std::string& filename, double length_multiplier) { set_mesh_pointer_from_filename(filename); set_length_multiplier(length_multiplier); - if (specified_length_multiplier_) { - libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); - } - m_->prepare_for_use(); initialize(); } @@ -2377,6 +2364,15 @@ void LibMesh::initialize() // assuming that unstructured meshes used in OpenMC are 3D n_dimension_ = 3; + if (specified_length_multiplier_) { + libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); + } + // if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh. + // Otherwise assume that it is prepared by it's owning application + if (unique_m_) { + m_->prepare_for_use(); + } + // ensure that the loaded mesh is 3 dimensional if (m_->mesh_dimension() != n_dimension_) { fatal_error(fmt::format("Mesh file {} specified for use in an unstructured " From 7c215543ebe7f3090683096b93047082b325179a Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Wed, 27 Jul 2022 12:53:19 -0500 Subject: [PATCH 64/83] added comment to XML constructor. initialized unique_m_ to nullptr to help some compilers --- include/openmc/mesh.h | 2 +- src/mesh.cpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 13138a739..02bd3bc5a 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -720,7 +720,7 @@ private: int get_bin_from_element(const libMesh::Elem* elem) const; // Data members - unique_ptr unique_m_; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC + unique_ptr unique_m_ = nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set during intialization vector> pl_; //!< per-thread point locators diff --git a/src/mesh.cpp b/src/mesh.cpp index 8e9d4bcec..dfa443b7f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2324,6 +2324,7 @@ const std::string LibMesh::mesh_lib_type = "libmesh"; LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) { + // filename_ and length_multiplier_ will already be set by the UnstructuredMesh constructor set_mesh_pointer_from_filename(filename_); set_length_multiplier(length_multiplier_); initialize(); From b3d0f20a0010dd930e4f403848d72164f32b294e Mon Sep 17 00:00:00 2001 From: yardasol Date: Wed, 27 Jul 2022 13:15:10 -0500 Subject: [PATCH 65/83] added add_nuclides method --- openmc/material.py | 14 ++++++++++++++ tests/unit_tests/test_material.py | 28 +++++++++++++++++++++++++++- 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index a837b1f17..cd181539c 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -403,6 +403,20 @@ class Material(IDManagerMixin): self._nuclides.append(NuclideTuple(nuclide, percent, percent_type)) + def add_nuclides(self, nuclides: dict): + """ Add multiple nuclides to a material + + Parameters + ---------- + nuclides : dict of str to tuple + Dictionary mapping nuclide names to a tuple containing their + atom or weight percent. + + """ + + for nuclide, (percent, percent_type) in nuclides.items(): + self.add_nuclide(nuclide, percent, percent_type) + def remove_nuclide(self, nuclide: str): """Remove a nuclide from the material diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index fd0e4f624..8dde30485 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -25,6 +25,32 @@ def test_add_nuclide(): with pytest.raises(ValueError): m.add_nuclide('H1', 1.0, 'oa') +def test_add_nuclides(): + """Test adding multipe nuclides at once""" + m = openmc.Material() + nuclides = {'H1': (2.0, 'ao'), + 'O16': (1.0, 'wo')} + m.add_nuclides(nuclides) + + # wt-% + m = openmc.Material() + nuclides = {'H1': (2.0, 'wo'), + 'O16': (1.0, 'wo')} + m.add_nuclides(nuclides) + + # mixed + m = openmc.Material() + nuclides = {'H1': (2.0, 'ao'), + 'O16': (1.0, 'wo')} + m.add_nuclides(nuclides) + + with pytest.raises(TypeError): + m.add_nuclides({'H1': ('1.0', 'ao')}) + with pytest.raises(TypeError): + m.add_nuclides({1.0: ('H1', 'wo')}) + with pytest.raises(ValueError): + m.add_nuclides({'H1': (1.0, 'oa')}) + def test_remove_nuclide(): """Test removing nuclides.""" @@ -441,7 +467,7 @@ def test_activity_of_tritium(): m1.add_nuclide("H3", 1) m1.set_density('g/cm3', 1) m1.volume = 1 - assert pytest.approx(m1.activity) == 3.559778e14 + assert pytest.approx(m1.activity) == 3.559778e14 def test_activity_of_metastable(): From 9ae0c6bd55e045681b1cef213fd6081f9b694534 Mon Sep 17 00:00:00 2001 From: Olek <45364492+yardasol@users.noreply.github.com> Date: Thu, 28 Jul 2022 08:50:48 -0500 Subject: [PATCH 66/83] add `add_elements` Co-authored-by: Jonathan Shimwell --- openmc/material.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/openmc/material.py b/openmc/material.py index cd181539c..40b495883 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -417,6 +417,18 @@ class Material(IDManagerMixin): for nuclide, (percent, percent_type) in nuclides.items(): self.add_nuclide(nuclide, percent, percent_type) + def add_elements(self, elements: dict): + """ Add multiple elements to a material + + Parameters + ---------- + elements : dict of str to tuple + Dictionary mapping element names to a tuple containing their + atom or weight percent. + """ + + for element, (percent, percent_type) in elements.items(): + self.add_element(element, percent, percent_type) def remove_nuclide(self, nuclide: str): """Remove a nuclide from the material From 3374ec12f5a305e83b4d853208bd1912be03c6b5 Mon Sep 17 00:00:00 2001 From: Lewis Gross <43077972+lewisgross1296@users.noreply.github.com> Date: Thu, 28 Jul 2022 09:39:56 -0500 Subject: [PATCH 67/83] Update src/mesh.cpp grammar Co-authored-by: Patrick Shriwise --- src/mesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index dfa443b7f..a8d57eef1 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2369,7 +2369,7 @@ void LibMesh::initialize() libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); } // if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh. - // Otherwise assume that it is prepared by it's owning application + // Otherwise assume that it is prepared by its owning application if (unique_m_) { m_->prepare_for_use(); } From 623658870d3a51320465aabed254acd64c7671f6 Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 28 Jul 2022 10:58:08 -0500 Subject: [PATCH 68/83] merged add_nuclides and add_elements into one function --- openmc/material.py | 45 ++++++++++++++---------- tests/unit_tests/test_material.py | 57 +++++++++++++++++-------------- 2 files changed, 59 insertions(+), 43 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 40b495883..1463c7fa1 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -403,32 +403,41 @@ class Material(IDManagerMixin): self._nuclides.append(NuclideTuple(nuclide, percent, percent_type)) - def add_nuclides(self, nuclides: dict): - """ Add multiple nuclides to a material + def add_elements_or_nuclides(self, components: dict): + """ Add multiple elements or nuclides to a material Parameters ---------- - nuclides : dict of str to tuple - Dictionary mapping nuclide names to a tuple containing their - atom or weight percent. + components : dict of str to tuple + Dictionary mapping element or nuclide names to a tuple containing + their atom or weight percent, as well as any other arguments. + + Examples + -------- + >> mat = openmc.Material() + >> components = {'Li': (1.0, {'enrichment': 60.0, 'enrichment_target': 'Li7'}), + 'Fl': 1.0, + 'Be6': (0.5, 'wo')} + >> mat.add_elements_or_nuclides(components) """ - for nuclide, (percent, percent_type) in nuclides.items(): - self.add_nuclide(nuclide, percent, percent_type) + for component, args in components.items(): + cv.check_type('component', component, str) + percent_type = 'ao' + if isinstance(args, float): + args = (args,) - def add_elements(self, elements: dict): - """ Add multiple elements to a material + ## check if nuclide + if str.isdigit(component[-1]): + self.add_nuclide(component, *args) + else: # is element + kwargs = {} + if isinstance(args[-1], dict): + kwargs = args[-1] + args = args[:-1] + self.add_element(component, *args, **kwargs) - Parameters - ---------- - elements : dict of str to tuple - Dictionary mapping element names to a tuple containing their - atom or weight percent. - """ - - for element, (percent, percent_type) in elements.items(): - self.add_element(element, percent, percent_type) def remove_nuclide(self, nuclide: str): """Remove a nuclide from the material diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 8dde30485..69ca87f2e 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -25,31 +25,39 @@ def test_add_nuclide(): with pytest.raises(ValueError): m.add_nuclide('H1', 1.0, 'oa') -def test_add_nuclides(): - """Test adding multipe nuclides at once""" +def test_add_elements_or_nuclides(): + """Test adding multipe elements or nuclides at once""" m = openmc.Material() - nuclides = {'H1': (2.0, 'ao'), - 'O16': (1.0, 'wo')} - m.add_nuclides(nuclides) - - # wt-% - m = openmc.Material() - nuclides = {'H1': (2.0, 'wo'), - 'O16': (1.0, 'wo')} - m.add_nuclides(nuclides) - - # mixed - m = openmc.Material() - nuclides = {'H1': (2.0, 'ao'), - 'O16': (1.0, 'wo')} - m.add_nuclides(nuclides) - - with pytest.raises(TypeError): - m.add_nuclides({'H1': ('1.0', 'ao')}) - with pytest.raises(TypeError): - m.add_nuclides({1.0: ('H1', 'wo')}) + components = {'H1': 2.0, + 'O16': (1.0, 'wo'), + 'Zr': 1.0, + 'O': (1.0, 'wo'), + 'U': (1.0, {'enrichment': 4.5}), + 'Li': (1.0, 'ao', {'enrichment': 60.0, 'enrichment_target': 'Li7'}), + 'H': (1.0, 'ao', {'enrichment': 50.0, 'enrichment_target': 'H2', 'enrichment_type': 'wo'})} + m.add_elements_or_nuclides(components) with pytest.raises(ValueError): - m.add_nuclides({'H1': (1.0, 'oa')}) + m.add_elements_or_nuclides({'U': (1.0, 'ao', {'enrichment': 100.0})}) + with pytest.raises(ValueError): + m.add_elements_or_nuclides({'Pu': (1.0, 'ao', {'enrichment': 3.0})}) + with pytest.raises(ValueError): + m.add_elements_or_nuclides({'U': (1.0, 'ao', {'enrichment': 70.0, 'enrichment_target': 'U235'})}) + with pytest.raises(ValueError): + m.add_elements_or_nuclides({'He': (1.0, 'ao', {'enrichment': 17.0, 'enrichment_target': 'He6'})}) + with pytest.raises(ValueError): + m.add_elements_or_nuclides({'li': (1.0, 'ao')}) # should fail as 1st char is lowercase + with pytest.raises(ValueError): + m.add_elements_or_nuclides({'LI': (1.0, 'ao')}) # should fail as 2nd char is uppercase + with pytest.raises(ValueError): + m.add_elements_or_nuclides({'Xx': (1.0, 'ao')}) # should fail as Xx is not an element + with pytest.raises(ValueError): + m.add_elements_or_nuclides({'n': (1.0, 'ao')}) # check to avoid n for neutron being accepted + with pytest.raises(TypeError): + m.add_elements_or_nuclides({'H1': ('1.0', 'ao')}) + with pytest.raises(TypeError): + m.add_elements_or_nuclides({1.0: ('H1', 'wo')}) + with pytest.raises(ValueError): + m.add_elements_or_nuclides({'H1': (1.0, 'oa')}) def test_remove_nuclide(): @@ -75,7 +83,7 @@ def test_remove_elements(): assert m.nuclides[0].percent == 1.0 -def test_elements(): +def test_add_element(): """Test adding elements.""" m = openmc.Material() m.add_element('Zr', 1.0) @@ -100,7 +108,6 @@ def test_elements(): with pytest.raises(ValueError): m.add_element('n', 1.0) # check to avoid n for neutron being accepted - def test_elements_by_name(): """Test adding elements by name""" m = openmc.Material() From 9317724a5ba3ff2b6a8ea15eb34a1e3628828c41 Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 28 Jul 2022 16:03:42 -0500 Subject: [PATCH 69/83] fix examples formatting --- openmc/material.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 1463c7fa1..8fee0585c 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -414,11 +414,11 @@ class Material(IDManagerMixin): Examples -------- - >> mat = openmc.Material() - >> components = {'Li': (1.0, {'enrichment': 60.0, 'enrichment_target': 'Li7'}), - 'Fl': 1.0, - 'Be6': (0.5, 'wo')} - >> mat.add_elements_or_nuclides(components) + >>> mat = openmc.Material() + >>> components = {'Li': (1.0, {'enrichment': 60.0, 'enrichment_target': 'Li7'}), + >>> 'Fl': 1.0, + >>> 'Be6': (0.5, 'wo')} + >>> mat.add_elements_or_nuclides(components) """ From 6351b49ea2acf13f90914ef756c6fbf89a296d0a Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 28 Jul 2022 16:06:19 -0500 Subject: [PATCH 70/83] add reference to add_elements_or_nuclides in the docstring for Material --- openmc/material.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index 8fee0585c..14bec12f3 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -34,7 +34,9 @@ class Material(IDManagerMixin): To create a material, one should create an instance of this class, add nuclides or elements with :meth:`Material.add_nuclide` or :meth:`Material.add_element`, respectively, and set the total material - density with :meth:`Material.set_density()`. The material can then be + density with :meth:`Material.set_density()`. Alternatively, you can + use :meth:`Material.add_elements_or_nuclides()` to pass a dictionary + containing all the component information. The material can then be assigned to a cell using the :attr:`Cell.fill` attribute. Parameters From 25ab9d28ce0c45f886639681aa9233eb1a2dfac1 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 29 Jul 2022 00:47:27 +0100 Subject: [PATCH 71/83] Apply suggestions from @paulromano Co-authored-by: Paul Romano --- openmc/mesh.py | 6 +++--- src/state_point.cpp | 1 - tests/unit_tests/mesh_to_vtk/test.py | 11 +++++------ 3 files changed, 8 insertions(+), 10 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 9c837fd8e..211c4075d 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1485,7 +1485,7 @@ class UnstructuredMesh(MeshBase): element_types : Iterable of integers Mesh element types - .. versionadded:: 0.13 + .. versionadded:: 0.13.1 total_volume : float Volume of the unstructured mesh in total """ @@ -1770,8 +1770,8 @@ class UnstructuredMesh(MeshBase): vertices = group['vertices'][()] mesh._vertices = vertices.reshape((-1, 3)) - connectvity = group['connectivity'][()] - mesh._connectivity = connectvity.reshape((-1, 8)) + connectivity = group['connectivity'][()] + mesh._connectivity = connectivity.reshape((-1, 8)) mesh._element_types = group['element_types'][()] if 'length_multiplier' in group: diff --git a/src/state_point.cpp b/src/state_point.cpp index b4200ffe2..11c93447d 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -812,7 +812,6 @@ void write_unstructured_mesh_results() auto umesh = dynamic_cast(model::meshes[mesh_idx].get()); - if (!umesh) continue; diff --git a/tests/unit_tests/mesh_to_vtk/test.py b/tests/unit_tests/mesh_to_vtk/test.py index 4ac875002..d9f050226 100644 --- a/tests/unit_tests/mesh_to_vtk/test.py +++ b/tests/unit_tests/mesh_to_vtk/test.py @@ -17,11 +17,10 @@ def ids_func(param): test_params = (['libmesh', 'moab'], ['tets', 'hexes']) -test_cases = [] -for library, elem_type in product(*test_params): - test_case = {'library' : library, - 'elem_type' : elem_type} - test_cases.append(test_case) +test_cases = [ + {'library' : library, 'elem_type' : elem_type} + for library, elem_type in product(*test_params) +] @pytest.mark.parametrize("test_opts", test_cases, ids=ids_func) def test_unstructured_mesh_to_vtk(run_in_tmpdir, request, test_opts): @@ -60,7 +59,7 @@ def test_unstructured_mesh_to_vtk(run_in_tmpdir, request, test_opts): with openmc.StatePoint(sp_file) as sp: umesh = sp.meshes[umesh.id] - test_data = {'ids' : np.arange(umesh.n_elements)} + test_data = {'ids': np.arange(umesh.n_elements)} umesh.write_data_to_vtk('umesh.vtk', datasets=test_data, volume_normalization=False) From 4564607eff50e69c77071e52ff71242d1d0a6083 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 18:48:11 -0500 Subject: [PATCH 72/83] Addressing comments from @paulromano --- docs/source/io_formats/statepoint.rst | 6 ++++-- openmc/mesh.py | 19 +++++++++++++------ src/mesh.cpp | 2 +- 3 files changed, 18 insertions(+), 9 deletions(-) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index f96d7181d..26dbd9c2e 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -80,8 +80,10 @@ The current version of the statepoint file format is 17.0. dimension. - **Unstructured Mesh Only:** - **volumes** (*double[]*) -- Volume of each mesh cell. - - **centroids** (*double[]*) -- Location of the mesh cell - centroids. + - **vertices** (*double[]*) -- Location of the mesh vertices. + - **connectivity** (*int[]*) -- Connectivity array for the mesh + cells. + - **element_types** (*int[]*) -- Mesh element types. **/tallies/filters/** diff --git a/openmc/mesh.py b/openmc/mesh.py index 211c4075d..f2aced9f8 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1477,7 +1477,7 @@ class UnstructuredMesh(MeshBase): .. versionadded:: 0.13 - connectivity : np.ndarray (8, n_elements) + connectivity : np.ndarray (n_elements, 8) Connectivity of the elements .. versionadded:: 0.13 @@ -1633,6 +1633,17 @@ class UnstructuredMesh(MeshBase): def centroid(self, bin): """Return the vertex averaged centroid of an element + + Parameters + ---------- + bin : int + Bin ID for the returned centroid + + Returns + ------- + numpy.ndarray + x, y, z values of the element centroid + """ conn = self.connectivity[bin] coords = self.vertices[conn] @@ -1749,11 +1760,7 @@ class UnstructuredMesh(MeshBase): arr.SetTuple1(i, data.flat[i]) grid.GetCellData().AddArray(arr) - if vtk.VTK_MAJOR_VERSION == 5: - grid.update() - writer.SetInput(grid) - else: - writer.SetInputData(grid) + writer.SetInputData(grid) writer.Write() diff --git a/src/mesh.cpp b/src/mesh.cpp index e701f954e..df5c7f2b3 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -256,7 +256,7 @@ void UnstructuredMesh::to_hdf5(hid_t group) const } else { num_elem_skipped++; xt::view(elem_types, i, xt::all()) = static_cast(ElementType::UNSUPPORTED); - xt::view(connectivity, i, xt::all()) = xt::xarray({-1, -1, -1, -1, -1, -1, -1, -1}); + xt::view(connectivity, i, xt::all()) = -1; } } From 8e747cc8747410326da04fbd4fcebc587a4a0081 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 29 Jul 2022 00:51:25 +0100 Subject: [PATCH 73/83] Apply suggestion from @paulromano Co-authored-by: Paul Romano --- openmc/mesh.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 211c4075d..d7ea4bc11 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1469,19 +1469,17 @@ class UnstructuredMesh(MeshBase): be generated for this mesh volumes : Iterable of float Volumes of the unstructured mesh elements - centroids : np.narray (3, n_elements) - Centroids of the mesh elements + centroids : numpy.ndarray + Centroids of the mesh elements with array shape (3, n_elements) - vertices : np.ndarray (3,) - Coordinates of the mesh vertices + vertices : numpy.ndarray + Coordinates of the mesh vertices with array shape (3, n_elements) - .. versionadded:: 0.13 - - connectivity : np.ndarray (8, n_elements) - Connectivity of the elements - - .. versionadded:: 0.13 + .. versionadded:: 0.13.1 + connectivity : numpy.ndarray + Connectivity of the elements with array shape (8, n_elements) + .. versionadded:: 0.13.1 element_types : Iterable of integers Mesh element types From 044c28d6519b76ef72d8c4b79f48f0faeff78fd7 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 18:50:21 -0500 Subject: [PATCH 74/83] Applying clang-format --- src/mesh.cpp | 8 ++++++-- src/state_point.cpp | 3 ++- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index df5c7f2b3..aa35f3064 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -262,8 +262,12 @@ void UnstructuredMesh::to_hdf5(hid_t group) const // warn users that some elements were skipped if (num_elem_skipped > 0) { - warning(fmt::format("The connectivity of {} elements on mesh {} were not written " - "because they are not of type linear tet/hex.", num_elem_skipped, this->id_)); + warning(fmt::format( + "The connectivity of {} elements + on mesh {} were not written + " + "because they are not of type linear tet/hex.", + num_elem_skipped, this->id_)); } write_dataset(mesh_group, "volumes", volumes); diff --git a/src/state_point.cpp b/src/state_point.cpp index 11c93447d..3b217a658 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -820,7 +820,8 @@ void write_unstructured_mesh_results() if (umesh->library() == "moab") { if (mpi::master) - warning(fmt::format("Output for a MOAB mesh (mesh {}) was requested but will not be written. Please use the Python API to generated the desired VTK tetrahedral mesh.", umesh->id_)); + warning(fmt::format("Output for a MOAB mesh (mesh {}) was + requested but will not be written. Please use the Python API to generated the desired VTK tetrahedral mesh.", umesh->id_)); continue; } From 2edd390d30d98c64282268034a469e0e607b0df3 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 18:56:19 -0500 Subject: [PATCH 75/83] Newline at end of test file --- tests/unit_tests/mesh_to_vtk/test.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/mesh_to_vtk/test.py b/tests/unit_tests/mesh_to_vtk/test.py index d9f050226..c433913ee 100644 --- a/tests/unit_tests/mesh_to_vtk/test.py +++ b/tests/unit_tests/mesh_to_vtk/test.py @@ -66,4 +66,5 @@ def test_unstructured_mesh_to_vtk(run_in_tmpdir, request, test_opts): # compare file content with reference file ref_file = Path(f"{test_opts['library']}_{test_opts['elem_type']}_ref.vtk") - assert filecmp.cmp('umesh.vtk', request.node.path.parent / ref_file) \ No newline at end of file + assert filecmp.cmp('umesh.vtk', request.node.path.parent / ref_file) + \ No newline at end of file From d5e963220601df0739449172cd249b24cf66c538 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 19:19:39 -0500 Subject: [PATCH 76/83] Updating statepoint doc a bit --- docs/source/io_formats/statepoint.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 26dbd9c2e..28a6849fd 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -80,7 +80,7 @@ The current version of the statepoint file format is 17.0. dimension. - **Unstructured Mesh Only:** - **volumes** (*double[]*) -- Volume of each mesh cell. - - **vertices** (*double[]*) -- Location of the mesh vertices. + - **vertices** (*double[]*) -- x, y, z values of the mesh vertices. - **connectivity** (*int[]*) -- Connectivity array for the mesh cells. - **element_types** (*int[]*) -- Mesh element types. From bdf6de5e516bc18a38f786cb09981905c74379fd Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 19:21:17 -0500 Subject: [PATCH 77/83] EOL for test --- tests/unit_tests/mesh_to_vtk/test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/unit_tests/mesh_to_vtk/test.py b/tests/unit_tests/mesh_to_vtk/test.py index c433913ee..5f6cd6a41 100644 --- a/tests/unit_tests/mesh_to_vtk/test.py +++ b/tests/unit_tests/mesh_to_vtk/test.py @@ -67,4 +67,3 @@ def test_unstructured_mesh_to_vtk(run_in_tmpdir, request, test_opts): # compare file content with reference file ref_file = Path(f"{test_opts['library']}_{test_opts['elem_type']}_ref.vtk") assert filecmp.cmp('umesh.vtk', request.node.path.parent / ref_file) - \ No newline at end of file From 7b083c7b0ce68b1a615ef3df521f6d253a1c20bc Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 29 Jul 2022 02:03:03 -0500 Subject: [PATCH 78/83] Fixing output message syntax --- src/mesh.cpp | 8 +++----- src/state_point.cpp | 7 +++++-- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index aa35f3064..a8ee8e73f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -262,11 +262,9 @@ void UnstructuredMesh::to_hdf5(hid_t group) const // warn users that some elements were skipped if (num_elem_skipped > 0) { - warning(fmt::format( - "The connectivity of {} elements - on mesh {} were not written - " - "because they are not of type linear tet/hex.", + warning(fmt::format("The connectivity of {} elements " + "on mesh {} were not written " + "because they are not of type linear tet/hex.", num_elem_skipped, this->id_)); } diff --git a/src/state_point.cpp b/src/state_point.cpp index 3b217a658..648d1a249 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -820,8 +820,11 @@ void write_unstructured_mesh_results() if (umesh->library() == "moab") { if (mpi::master) - warning(fmt::format("Output for a MOAB mesh (mesh {}) was - requested but will not be written. Please use the Python API to generated the desired VTK tetrahedral mesh.", umesh->id_)); + warning(fmt::format( + "Output for a MOAB mesh (mesh {}) was " + "requested but will not be written. Please use the Python " + "API to generated the desired VTK tetrahedral mesh.", + umesh->id_)); continue; } From 5bd2a14a2cd021594b52fe90a9475e145f1838a5 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 29 Jul 2022 15:30:35 -0500 Subject: [PATCH 79/83] Adding missing datasets to the statepoint docs --- docs/source/io_formats/statepoint.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 28a6849fd..c1598f2c3 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -79,6 +79,10 @@ The current version of the statepoint file format is 17.0. - **width** (*double[]*) -- Width of each mesh cell in each dimension. - **Unstructured Mesh Only:** + - **filename** (*char[]*) -- Name of the mesh file. + - **library** (*char[]*) -- Mesh library used to represent the + mesh ("moab" or "libmesh"). + - **length_multiplier** (*double*) Scaling factor applied to the mesh. - **volumes** (*double[]*) -- Volume of each mesh cell. - **vertices** (*double[]*) -- x, y, z values of the mesh vertices. - **connectivity** (*int[]*) -- Connectivity array for the mesh From 647863c7ace4c6c5736ae75ba04c82b258897b41 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 29 Jul 2022 15:33:29 -0500 Subject: [PATCH 80/83] Removing centroids setter --- openmc/mesh.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 2b405bf6a..bb559a383 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1594,13 +1594,6 @@ class UnstructuredMesh(MeshBase): cv.check_type('Number of elements', val, Integral) self._n_elements = val - - @centroids.setter - def centroids(self, centroids): - cv.check_type("Unstructured mesh centroids", centroids, - Iterable, Real) - self._centroids = centroids - @property def length_multiplier(self): return self._length_multiplier From a3e2fcdcca22184f87ecf8ebbdddecf770360057 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 29 Jul 2022 15:34:23 -0500 Subject: [PATCH 81/83] Correcting dimensions of attributes in um docstrings --- openmc/mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index bb559a383..daff4dedc 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1470,14 +1470,14 @@ class UnstructuredMesh(MeshBase): volumes : Iterable of float Volumes of the unstructured mesh elements centroids : numpy.ndarray - Centroids of the mesh elements with array shape (3, n_elements) + Centroids of the mesh elements with array shape (n_elements, 3) vertices : numpy.ndarray - Coordinates of the mesh vertices with array shape (3, n_elements) + Coordinates of the mesh vertices with array shape (n_elements, 3) .. versionadded:: 0.13.1 connectivity : numpy.ndarray - Connectivity of the elements with array shape (8, n_elements) + Connectivity of the elements with array shape (n_elements, 8) .. versionadded:: 0.13.1 element_types : Iterable of integers From effc35b8eed3c33c63a56fdead9120143f5551ed Mon Sep 17 00:00:00 2001 From: yardasol Date: Fri, 29 Jul 2022 18:21:06 -0500 Subject: [PATCH 82/83] add_elements_or_nuclides -> add_components - changed the assumptions the function makes about the 'components' parameter - added 'percent_type' parameter - passed changes to unit test --- openmc/material.py | 48 ++++++++++++++++++----------- tests/unit_tests/test_material.py | 50 +++++++++++++++++++------------ 2 files changed, 62 insertions(+), 36 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 14bec12f3..9ba7360ff 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -35,7 +35,7 @@ class Material(IDManagerMixin): nuclides or elements with :meth:`Material.add_nuclide` or :meth:`Material.add_element`, respectively, and set the total material density with :meth:`Material.set_density()`. Alternatively, you can - use :meth:`Material.add_elements_or_nuclides()` to pass a dictionary + use :meth:`Material.add_components()` to pass a dictionary containing all the component information. The material can then be assigned to a cell using the :attr:`Cell.fill` attribute. @@ -405,41 +405,55 @@ class Material(IDManagerMixin): self._nuclides.append(NuclideTuple(nuclide, percent, percent_type)) - def add_elements_or_nuclides(self, components: dict): + def add_components(self, components: dict, percent_type: str = 'ao'): """ Add multiple elements or nuclides to a material Parameters ---------- - components : dict of str to tuple - Dictionary mapping element or nuclide names to a tuple containing - their atom or weight percent, as well as any other arguments. + components : dict of str to float or dict + Dictionary mapping element or nuclide names to their atom or weight + percent. To specify enrichment of an element, the entry of + ``components`` for that element must instead be a dictionary + containing the keyword arguments as well as a value for + ``'percent'`` + percent_type : {'ao', 'wo'} + 'ao' for atom percent and 'wo' for weight percent Examples -------- >>> mat = openmc.Material() - >>> components = {'Li': (1.0, {'enrichment': 60.0, 'enrichment_target': 'Li7'}), - >>> 'Fl': 1.0, - >>> 'Be6': (0.5, 'wo')} - >>> mat.add_elements_or_nuclides(components) + >>> components = {'Li': {'percent': 1.0, + >>> 'enrichment': 60.0, + >>> 'enrichment_target': 'Li7'}, + >>> 'Fl': 1.0, + >>> 'Be6': 0.5} + >>> mat.add_components(components) """ - for component, args in components.items(): + for component, params in components.items(): cv.check_type('component', component, str) - percent_type = 'ao' - if isinstance(args, float): - args = (args,) + if isinstance(params, float): + params = {'percent': params} + + else: + cv.check_type('params', params, dict) + if 'percent' not in params: + raise ValueError("An entry in the dictionary does not have " + "a required key: 'percent'") + + percent = params.pop('percent') + args = (percent, percent_type) ## check if nuclide if str.isdigit(component[-1]): self.add_nuclide(component, *args) else: # is element - kwargs = {} - if isinstance(args[-1], dict): - kwargs = args[-1] - args = args[:-1] + kwargs = params self.add_element(component, *args, **kwargs) + params['percent'] = percent + def remove_nuclide(self, nuclide: str): """Remove a nuclide from the material diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 69ca87f2e..2c5cc2874 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -25,39 +25,51 @@ def test_add_nuclide(): with pytest.raises(ValueError): m.add_nuclide('H1', 1.0, 'oa') -def test_add_elements_or_nuclides(): +def test_add_components(): """Test adding multipe elements or nuclides at once""" m = openmc.Material() components = {'H1': 2.0, - 'O16': (1.0, 'wo'), - 'Zr': 1.0, - 'O': (1.0, 'wo'), - 'U': (1.0, {'enrichment': 4.5}), - 'Li': (1.0, 'ao', {'enrichment': 60.0, 'enrichment_target': 'Li7'}), - 'H': (1.0, 'ao', {'enrichment': 50.0, 'enrichment_target': 'H2', 'enrichment_type': 'wo'})} - m.add_elements_or_nuclides(components) + 'O16': 1.0, + 'Zr': 1.0, + 'O': 1.0, + 'U': {'percent': 1.0, + 'enrichment': 4.5}, + 'Li': {'percent': 1.0, + 'enrichment': 60.0, + 'enrichment_target': 'Li7'}, + 'H': {'percent': 1.0, + 'enrichment': 50.0, + 'enrichment_target': 'H2', + 'enrichment_type': 'wo'}} + m.add_components(components) with pytest.raises(ValueError): - m.add_elements_or_nuclides({'U': (1.0, 'ao', {'enrichment': 100.0})}) + m.add_components({'U': {'percent': 1.0, + 'enrichment': 100.0}}) with pytest.raises(ValueError): - m.add_elements_or_nuclides({'Pu': (1.0, 'ao', {'enrichment': 3.0})}) + m.add_components({'Pu': {'percent': 1.0, + 'enrichment': 3.0}}) with pytest.raises(ValueError): - m.add_elements_or_nuclides({'U': (1.0, 'ao', {'enrichment': 70.0, 'enrichment_target': 'U235'})}) + m.add_components({'U': {'percent': 1.0, + 'enrichment': 70.0, + 'enrichment_target':'U235'}}) with pytest.raises(ValueError): - m.add_elements_or_nuclides({'He': (1.0, 'ao', {'enrichment': 17.0, 'enrichment_target': 'He6'})}) + m.add_components({'He': {'percent': 1.0, + 'enrichment': 17.0, + 'enrichment_target': 'He6'}}) with pytest.raises(ValueError): - m.add_elements_or_nuclides({'li': (1.0, 'ao')}) # should fail as 1st char is lowercase + m.add_components({'li': 1.0}) # should fail as 1st char is lowercase with pytest.raises(ValueError): - m.add_elements_or_nuclides({'LI': (1.0, 'ao')}) # should fail as 2nd char is uppercase + m.add_components({'LI': 1.0}) # should fail as 2nd char is uppercase with pytest.raises(ValueError): - m.add_elements_or_nuclides({'Xx': (1.0, 'ao')}) # should fail as Xx is not an element + m.add_components({'Xx': 1.0}) # should fail as Xx is not an element with pytest.raises(ValueError): - m.add_elements_or_nuclides({'n': (1.0, 'ao')}) # check to avoid n for neutron being accepted + m.add_components({'n': 1.0}) # check to avoid n for neutron being accepted with pytest.raises(TypeError): - m.add_elements_or_nuclides({'H1': ('1.0', 'ao')}) + m.add_components({'H1': '1.0'}) with pytest.raises(TypeError): - m.add_elements_or_nuclides({1.0: ('H1', 'wo')}) + m.add_components({1.0: 'H1'}, percent_type = 'wo') with pytest.raises(ValueError): - m.add_elements_or_nuclides({'H1': (1.0, 'oa')}) + m.add_components({'H1': 1.0}, percent_type = 'oa') def test_remove_nuclide(): From 0dd8f70939d8103e11b205b559c2e71a254e4f22 Mon Sep 17 00:00:00 2001 From: yardasol Date: Mon, 1 Aug 2022 09:22:10 -0500 Subject: [PATCH 83/83] add paulromano's 2nd round of suggestions --- openmc/material.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 9ba7360ff..6c2641680 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -408,6 +408,8 @@ class Material(IDManagerMixin): def add_components(self, components: dict, percent_type: str = 'ao'): """ Add multiple elements or nuclides to a material + .. versionadded:: 0.13.1 + Parameters ---------- components : dict of str to float or dict @@ -442,17 +444,14 @@ class Material(IDManagerMixin): raise ValueError("An entry in the dictionary does not have " "a required key: 'percent'") - percent = params.pop('percent') - args = (percent, percent_type) + params['percent_type'] = percent_type ## check if nuclide if str.isdigit(component[-1]): - self.add_nuclide(component, *args) + self.add_nuclide(component, **params) else: # is element kwargs = params - self.add_element(component, *args, **kwargs) - - params['percent'] = percent + self.add_element(component, **params) def remove_nuclide(self, nuclide: str): """Remove a nuclide from the material