From 693314d097c79c1f596f336eb69e4a212072e5df Mon Sep 17 00:00:00 2001 From: Yuvraj Jain <20016119+nutcasev15@users.noreply.github.com> Date: Mon, 30 Oct 2023 20:51:31 +0000 Subject: [PATCH] Leaky and Albedo Boundary Conditions Implementation (#2724) Co-authored-by: Paul Romano --- docs/source/io_formats/summary.rst | 8 +- docs/source/usersguide/geometry.rst | 32 +++- include/openmc/boundary_condition.h | 39 ++++ include/openmc/surface.h | 8 +- openmc/model/funcs.py | 93 ++++++---- openmc/surface.py | 168 ++++++++++++++++-- src/boundary_condition.cpp | 12 ++ src/dagmc.cpp | 6 +- src/surface.cpp | 53 ++++-- tests/regression_tests/albedo_box/__init__.py | 0 .../regression_tests/albedo_box/geometry.xml | 13 ++ .../regression_tests/albedo_box/materials.xml | 9 + .../albedo_box/results_true.dat | 2 + .../regression_tests/albedo_box/settings.xml | 9 + tests/regression_tests/albedo_box/test.py | 6 + 15 files changed, 376 insertions(+), 82 deletions(-) create mode 100755 tests/regression_tests/albedo_box/__init__.py create mode 100644 tests/regression_tests/albedo_box/geometry.xml create mode 100644 tests/regression_tests/albedo_box/materials.xml create mode 100644 tests/regression_tests/albedo_box/results_true.dat create mode 100644 tests/regression_tests/albedo_box/settings.xml create mode 100755 tests/regression_tests/albedo_box/test.py diff --git a/docs/source/io_formats/summary.rst b/docs/source/io_formats/summary.rst index c140d98ae3..7d3ab94d9f 100644 --- a/docs/source/io_formats/summary.rst +++ b/docs/source/io_formats/summary.rst @@ -60,9 +60,11 @@ The current version of the summary file format is 6.0. - **coefficients** (*double[]*) -- Array of coefficients that define the surface. See :ref:`surface_element` for what coefficients are defined for each surface type. - - **boundary_condition** (*char[]*) -- Boundary condition applied to - the surface. Can be 'transmission', 'vacuum', 'reflective', or - 'periodic'. + - **boundary_type** (*char[]*) -- Boundary condition applied to + the surface. Can be 'transmission', 'vacuum', 'reflective', + 'periodic', or 'white'. + - **albedo** (*double*) -- Boundary albedo as a positive multiplier + of particle weight. If absent, it is assumed to be 1.0. - **geom_type** (*char[]*) -- Type of geometry used to create the cell. Either 'csg' or 'dagmc'. diff --git a/docs/source/usersguide/geometry.rst b/docs/source/usersguide/geometry.rst index 2382609a48..379642996c 100644 --- a/docs/source/usersguide/geometry.rst +++ b/docs/source/usersguide/geometry.rst @@ -172,13 +172,17 @@ surface. To specify a vacuum boundary condition, simply change the outer_surface = openmc.Sphere(r=100.0) outer_surface.boundary_type = 'vacuum' -Reflective and periodic boundary conditions can be set with the strings -'reflective' and 'periodic'. Vacuum and reflective boundary conditions can be -applied to any type of surface. Periodic boundary conditions can be applied to -pairs of planar surfaces. If there are only two periodic surfaces they will be -matched automatically. Otherwise it is necessary to specify pairs explicitly -using the :attr:`Surface.periodic_surface` attribute as in the following -example:: +Reflective, periodic, and white boundary conditions can be set with the +strings 'reflective', 'periodic', and 'white' respectively. +Vacuum, reflective and white boundary conditions can be applied to any +type of surface. The 'white' boundary condition supports diffuse particle +reflection in contrast to specular reflection provided by the 'reflective' +boundary condition. + +Periodic boundary conditions can be applied to pairs of planar surfaces. +If there are only two periodic surfaces they will be matched automatically. +Otherwise it is necessary to specify pairs explicitly using the +:attr:`Surface.periodic_surface` attribute as in the following example:: p1 = openmc.Plane(a=0.3, b=5.0, d=1.0, boundary_type='periodic') p2 = openmc.Plane(a=0.3, b=5.0, d=-1.0, boundary_type='periodic') @@ -196,6 +200,20 @@ lies in the first quadrant of the Cartesian grid. If the geometry instead lies in the fourth quadrant, the :class:`YPlane` must be replaced by a :class:`Plane` with the normal vector pointing in the :math:`-y` direction. +Additionally, 'reflective', 'periodic', and 'white' boundary conditions have +an albedo parameter that can be used to modify the importance of particles +that encounter the boundary. The albedo value specifies the ratio between +the particle's importance after interaction with the boundary to its initial +importance. The following example creates a reflective planar surface which +reduces the reflected particles' importance by 33.3%:: + + x1 = openmc.XPlane(1.0, boundary_type='reflective', albedo=0.667) + + # This is equivalent + x1 = openmc.XPlane(1.0) + x1.boundary_type = 'reflective' + x1.albedo = 0.667 + .. _usersguide_cells: ----- diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index 0f16fa3b31..fde860dbca 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -1,7 +1,10 @@ #ifndef OPENMC_BOUNDARY_CONDITION_H #define OPENMC_BOUNDARY_CONDITION_H +#include "openmc/hdf5_interface.h" +#include "openmc/particle.h" #include "openmc/position.h" +#include namespace openmc { @@ -22,8 +25,44 @@ public: //! \param surf The specific surface on the boundary the particle struck. virtual void handle_particle(Particle& p, const Surface& surf) const = 0; + //! Modify the incident particle's weight according to the boundary's albedo. + //! \param p The particle that struck the boundary. This function calculates + //! the reduction in the incident particle's weight as it interacts + //! with a boundary. The lost weight is tallied before the remaining weight + //! is reassigned to the incident particle. Implementations of the + //! handle_particle function typically call this method in its body. + //! \param surf The specific surface on the boundary the particle struck. + void handle_albedo(Particle& p, const Surface& surf) const + { + if (!has_albedo()) + return; + double initial_wgt = p.wgt(); + // Treat the lost weight fraction as leakage, similar to VacuumBC. + // This ensures the lost weight is tallied properly. + p.wgt() *= (1.0 - albedo_); + p.cross_vacuum_bc(surf); + p.wgt() = initial_wgt * albedo_; + }; + //! Return a string classification of this BC. virtual std::string type() const = 0; + + //! Write albedo data of this BC to hdf5. + void to_hdf5(hid_t surf_group) const + { + if (has_albedo()) { + write_string(surf_group, "albedo", fmt::format("{}", albedo_), false); + } + }; + + //! Set albedo of this BC. + void set_albedo(double albedo) { albedo_ = albedo; } + + //! Return if this BC has an albedo. + bool has_albedo() const { return (albedo_ > 0.0); } + +private: + double albedo_ = -1.0; }; //============================================================================== diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 4c97a05145..3bb84c6e00 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -83,10 +83,10 @@ struct BoundingBox { class Surface { public: - int id_; //!< Unique ID - std::string name_; //!< User-defined name - std::shared_ptr bc_ {nullptr}; //!< Boundary condition - GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC) + int id_; //!< Unique ID + std::string name_; //!< User-defined name + unique_ptr bc_; //!< Boundary condition + GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC) bool surf_source_ {false}; //!< Activate source banking for the surface? explicit Surface(pugi::xml_node surf_node); diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index cbf3566515..d185c13df6 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -109,14 +109,16 @@ def borated_water(boron_ppm, temperature=293., pressure=0.1013, temp_unit='K', # Define function to create a plane on given axis -def _plane(axis, name, value, boundary_type='transmission'): +def _plane(axis, name, value, boundary_type='transmission', albedo=1.): cls = getattr(openmc, f'{axis.upper()}Plane') return cls(value, name=f'{name} {axis}', - boundary_type=boundary_type) + boundary_type=boundary_type, + albedo=albedo) def rectangular_prism(width, height, axis='z', origin=(0., 0.), - boundary_type='transmission', corner_radius=0.): + boundary_type='transmission', albedo=1., + corner_radius=0.): """Get an infinite rectangular prism from four planar surfaces. .. versionchanged:: 0.11 @@ -138,9 +140,14 @@ def rectangular_prism(width, height, axis='z', origin=(0., 0.), Origin of the prism. The two floats correspond to (y,z), (x,z) or (x,y) for prisms parallel to the x, y or z axis, respectively. Defaults to (0., 0.). - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surfaces comprising the rectangular prism (default is 'transmission'). + albedo : float, optional + Albedo of the prism's surfaces as a ratio of particle weight after + interaction with the surface to the initial weight. Values must be + positive. Only applicable if the boundary type is 'reflective', + 'periodic', or 'white'. corner_radius: float Prism corner radius in units of cm. Defaults to 0. @@ -153,6 +160,7 @@ def rectangular_prism(width, height, axis='z', origin=(0., 0.), check_type('width', width, Real) check_type('height', height, Real) + check_type('albedo', albedo, Real) check_type('corner_radius', corner_radius, Real) check_value('axis', axis, ['x', 'y', 'z']) check_type('origin', origin, Iterable, Real) @@ -167,15 +175,14 @@ def rectangular_prism(width, height, axis='z', origin=(0., 0.), # Get cylinder class corresponding to given axis cyl = getattr(openmc, f'{axis.upper()}Cylinder') + # Create container for boundary arguments + bc_args = {'boundary_type': boundary_type, 'albedo': albedo} + # Create rectangular region - min_x1 = _plane(x1, 'minimum', -width/2 + origin[0], - boundary_type=boundary_type) - max_x1 = _plane(x1, 'maximum', width/2 + origin[0], - boundary_type=boundary_type) - min_x2 = _plane(x2, 'minimum', -height/2 + origin[1], - boundary_type=boundary_type) - max_x2 = _plane(x2, 'maximum', height/2 + origin[1], - boundary_type=boundary_type) + min_x1 = _plane(x1, 'minimum', -width/2 + origin[0], **bc_args) + max_x1 = _plane(x1, 'maximum', width/2 + origin[0], **bc_args) + min_x2 = _plane(x2, 'minimum', -height/2 + origin[1], **bc_args) + max_x2 = _plane(x2, 'maximum', height/2 + origin[1], **bc_args) if boundary_type == 'periodic': min_x1.periodic_surface = max_x1 min_x2.periodic_surface = max_x2 @@ -187,7 +194,7 @@ def rectangular_prism(width, height, axis='z', origin=(0., 0.), raise ValueError('Periodic boundary conditions not permitted when ' 'rounded corners are used.') - args = {'r': corner_radius, 'boundary_type': boundary_type} + args = {'r': corner_radius, 'boundary_type': boundary_type, 'albedo' : albedo} args[x1 + '0'] = origin[0] - width/2 + corner_radius args[x2 + '0'] = origin[1] - height/2 + corner_radius @@ -210,13 +217,13 @@ def rectangular_prism(width, height, axis='z', origin=(0., 0.), x1_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args) x1_min = _plane(x1, 'min', -width/2 + origin[0] + corner_radius, - boundary_type=boundary_type) + **bc_args) x1_max = _plane(x1, 'max', width/2 + origin[0] - corner_radius, - boundary_type=boundary_type) + **bc_args) x2_min = _plane(x2, 'min', -height/2 + origin[1] + corner_radius, - boundary_type=boundary_type) + **bc_args) x2_max = _plane(x2, 'max', height/2 + origin[1] - corner_radius, - boundary_type=boundary_type) + **bc_args) corners = (+x1_min_x2_min & -x1_min & -x2_min) | \ (+x1_min_x2_max & -x1_min & +x2_max) | \ @@ -236,7 +243,7 @@ def get_rectangular_prism(*args, **kwargs): def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), - boundary_type='transmission', corner_radius=0.): + boundary_type='transmission', albedo=1., corner_radius=0.): """Create a hexagon region from six surface planes. .. versionchanged:: 0.11 @@ -253,9 +260,14 @@ def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), parallel to the y-axis. origin: Iterable of two floats Origin of the prism. Defaults to (0., 0.). - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surfaces comprising the hexagonal prism (default is 'transmission'). + albedo : float, optional + Albedo of the prism's surfaces as a ratio of particle weight after + interaction with the surface to the initial weight. Values must be + positive. Only applicable if the boundary type is 'reflective', + 'periodic', or 'white'. corner_radius: float Prism corner radius in units of cm. Defaults to 0. @@ -266,25 +278,35 @@ def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), """ + check_type('edge_length', edge_length, Real) + check_type('albedo', albedo, Real) + check_type('corner_radius', corner_radius, Real) + check_value('orientation', orientation, ['x', 'y']) + check_type('origin', origin, Iterable, Real) + l = edge_length x, y = origin + + # Create container for boundary arguments + bc_args = {'boundary_type': boundary_type, 'albedo' : albedo} + if orientation == 'y': - right = openmc.XPlane(x + sqrt(3.)/2*l, boundary_type=boundary_type) - left = openmc.XPlane(x - sqrt(3.)/2*l, boundary_type=boundary_type) + right = openmc.XPlane(x + sqrt(3.)/2*l, **bc_args) + left = openmc.XPlane(x - sqrt(3.)/2*l, **bc_args) c = sqrt(3.)/3. # y = -x/sqrt(3) + a - upper_right = Plane(a=c, b=1., d=l+x*c+y, boundary_type=boundary_type) + upper_right = Plane(a=c, b=1., d=l+x*c+y, **bc_args) # y = x/sqrt(3) + a - upper_left = Plane(a=-c, b=1., d=l-x*c+y, boundary_type=boundary_type) + upper_left = Plane(a=-c, b=1., d=l-x*c+y, **bc_args) # y = x/sqrt(3) - a - lower_right = Plane(a=-c, b=1., d=-l-x*c+y, boundary_type=boundary_type) + lower_right = Plane(a=-c, b=1., d=-l-x*c+y, **bc_args) # y = -x/sqrt(3) - a - lower_left = Plane(a=c, b=1., d=-l+x*c+y, boundary_type=boundary_type) + lower_left = Plane(a=c, b=1., d=-l+x*c+y, **bc_args) prism = -right & +left & -upper_right & -upper_left & \ +lower_right & +lower_left @@ -295,22 +317,21 @@ def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), lower_right.periodic_surface = upper_left elif orientation == 'x': - top = openmc.YPlane(y0=y + sqrt(3.)/2*l, boundary_type=boundary_type) - bottom = openmc.YPlane(y0=y - sqrt(3.)/2*l, boundary_type=boundary_type) + top = openmc.YPlane(y0=y + sqrt(3.)/2*l, **bc_args) + bottom = openmc.YPlane(y0=y - sqrt(3.)/2*l, **bc_args) c = sqrt(3.) # y = -sqrt(3)*(x - a) - upper_right = Plane(a=c, b=1., d=c*l+x*c+y, boundary_type=boundary_type) + upper_right = Plane(a=c, b=1., d=c*l+x*c+y, **bc_args) # y = sqrt(3)*(x + a) - lower_right = Plane(a=-c, b=1., d=-c*l-x*c+y, - boundary_type=boundary_type) + lower_right = Plane(a=-c, b=1., d=-c*l-x*c+y, **bc_args) # y = -sqrt(3)*(x + a) - lower_left = Plane(a=c, b=1., d=-c*l+x*c+y, boundary_type=boundary_type) + lower_left = Plane(a=c, b=1., d=-c*l+x*c+y, **bc_args) # y = sqrt(3)*(x + a) - upper_left = Plane(a=-c, b=1., d=c*l-x*c+y, boundary_type=boundary_type) + upper_left = Plane(a=-c, b=1., d=c*l-x*c+y, **bc_args) prism = -top & +bottom & -upper_right & +lower_right & \ +lower_left & -upper_left @@ -329,11 +350,9 @@ def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), c = sqrt(3.)/2 t = l - corner_radius/c - # Cylinder with corner radius and boundary type pre-applied - cyl1 = partial(openmc.ZCylinder, r=corner_radius, - boundary_type=boundary_type) - cyl2 = partial(openmc.ZCylinder, r=corner_radius/(2*c), - boundary_type=boundary_type) + # Cylinder with corner radius and boundary conditions pre-applied + cyl1 = partial(openmc.ZCylinder, r=corner_radius, **bc_args) + cyl2 = partial(openmc.ZCylinder, r=corner_radius/(2*c), **bc_args) if orientation == 'x': x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t) diff --git a/openmc/surface.py b/openmc/surface.py index dc1068b229..a684f83fea 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -8,13 +8,14 @@ from warnings import warn, catch_warnings, simplefilter import numpy as np -from .checkvalue import check_type, check_value, check_length +from .checkvalue import check_type, check_value, check_length, check_greater_than from .mixin import IDManagerMixin, IDWarning from .region import Region, Intersection, Union from .bounding_box import BoundingBox _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic', 'white'] +_ALBEDO_BOUNDARIES = ['reflective', 'periodic', 'white'] _WARNING_UPPER = """\ "{}(...) accepts an argument named '{}', not '{}'. Future versions of OpenMC \ @@ -123,6 +124,10 @@ class Surface(IDManagerMixin, ABC): freely pass through the surface. Note that periodic boundary conditions can only be applied to x-, y-, and z-planes, and only axis-aligned periodicity is supported. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the surface. If not specified, the name will be the empty string. @@ -132,6 +137,8 @@ class Surface(IDManagerMixin, ABC): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -147,10 +154,12 @@ class Surface(IDManagerMixin, ABC): used_ids = set() _atol = 1.e-12 - def __init__(self, surface_id=None, boundary_type='transmission', name=''): + def __init__(self, surface_id=None, boundary_type='transmission', + albedo=1., name=''): self.id = surface_id self.name = name self.boundary_type = boundary_type + self.albedo = albedo # A dictionary of the quadratic surface coefficients # Key - coefficient name @@ -165,15 +174,20 @@ class Surface(IDManagerMixin, ABC): def __repr__(self): string = 'Surface\n' - string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) - string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) - string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type) - string += '{0: <16}{1}{2}\n'.format('\tBoundary', '=\t', self._boundary_type) + string += '{0: <20}{1}{2}\n'.format('\tID', '=\t', self._id) + string += '{0: <20}{1}{2}\n'.format('\tName', '=\t', self._name) + string += '{0: <20}{1}{2}\n'.format('\tType', '=\t', self._type) + string += '{0: <20}{1}{2}\n'.format('\tBoundary', '=\t', + self._boundary_type) + if (self._boundary_type in _ALBEDO_BOUNDARIES and + not math.isclose(self._albedo, 1.0)): + string += '{0: <20}{1}{2}\n'.format('\tBoundary Albedo', '=\t', + self._albedo) - coefficients = '{0: <16}'.format('\tCoefficients') + '\n' + coefficients = '{0: <20}'.format('\tCoefficients') + '\n' for coeff in self._coefficients: - coefficients += '{0: <16}{1}{2}\n'.format( + coefficients += '{0: <20}{1}{2}\n'.format( coeff, '=\t', self._coefficients[coeff]) string += coefficients @@ -206,6 +220,16 @@ class Surface(IDManagerMixin, ABC): check_value('boundary type', boundary_type, _BOUNDARY_TYPES) self._boundary_type = boundary_type + @property + def albedo(self): + return self._albedo + + @albedo.setter + def albedo(self, albedo): + check_type('albedo', albedo, Real) + check_greater_than('albedo', albedo, 0.0) + self._albedo = float(albedo) + @property def coefficients(self): return self._coefficients @@ -402,6 +426,9 @@ class Surface(IDManagerMixin, ABC): element.set("type", self._type) if self.boundary_type != 'transmission': element.set("boundary", self.boundary_type) + if (self.boundary_type in _ALBEDO_BOUNDARIES and + not math.isclose(self.albedo, 1.0)): + element.set("albedo", str(self.albedo)) element.set("coeffs", ' '.join([str(self._coefficients.setdefault(key, 0.0)) for key in self._coeff_keys])) @@ -427,10 +454,12 @@ class Surface(IDManagerMixin, ABC): surf_type = elem.get('type') cls = _SURFACE_CLASSES[surf_type] - # Determine ID, boundary type, coefficients + # Determine ID, boundary type, boundary albedo, coefficients kwargs = {} kwargs['surface_id'] = int(elem.get('id')) kwargs['boundary_type'] = elem.get('boundary', 'transmission') + if kwargs['boundary_type'] in _ALBEDO_BOUNDARIES: + kwargs['albedo'] = float(elem.get('albedo', 1.0)) kwargs['name'] = elem.get('name') coeffs = [float(x) for x in elem.get('coeffs').split()] kwargs.update(dict(zip(cls._coeff_keys, coeffs))) @@ -462,8 +491,13 @@ class Surface(IDManagerMixin, ABC): name = group['name'][()].decode() if 'name' in group else '' bc = group['boundary_type'][()].decode() + if 'albedo' in group: + bc_alb = float(group['albedo'][()].decode()) + else: + bc_alb = 1.0 coeffs = group['coefficients'][...] - kwargs = {'boundary_type': bc, 'name': name, 'surface_id': surface_id} + kwargs = {'boundary_type': bc, 'albedo': bc_alb, 'name': name, + 'surface_id': surface_id} surf_type = group['type'][()].decode() cls = _SURFACE_CLASSES[surf_type] @@ -607,7 +641,9 @@ class PlaneMixin: # Compute new rotated coefficients a, b, c a, b, c = Rmat @ [a, b, c] - kwargs = {'boundary_type': surf.boundary_type, 'name': surf.name} + kwargs = {'boundary_type': surf.boundary_type, + 'albedo': surf.albedo, + 'name': surf.name} if inplace: kwargs['surface_id'] = surf.id @@ -651,6 +687,10 @@ class Plane(PlaneMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the plane. If not specified, the name will be the empty string. surface_id : int, optional @@ -670,6 +710,8 @@ class Plane(PlaneMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight periodic_surface : openmc.Surface If a periodic boundary condition is used, the surface with which this one is periodic with @@ -774,6 +816,10 @@ class XPlane(PlaneMixin, Surface): surface. Defaults to transmissive boundary condition where particles freely pass through the surface. Only axis-aligned periodicity is supported, i.e., x-planes can only be paired with x-planes. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the plane. If not specified, the name will be the empty string. surface_id : int, optional @@ -787,6 +833,8 @@ class XPlane(PlaneMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight periodic_surface : openmc.Surface If a periodic boundary condition is used, the surface with which this one is periodic with @@ -832,7 +880,11 @@ class YPlane(PlaneMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. Only axis-aligned periodicity is - supported, i.e., x-planes can only be paired with x-planes. + supported, i.e., y-planes can only be paired with y-planes. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the plane. If not specified, the name will be the empty string. surface_id : int, optional @@ -846,6 +898,8 @@ class YPlane(PlaneMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight periodic_surface : openmc.Surface If a periodic boundary condition is used, the surface with which this one is periodic with @@ -891,7 +945,11 @@ class ZPlane(PlaneMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. Only axis-aligned periodicity is - supported, i.e., x-planes can only be paired with x-planes. + supported, i.e., z-planes can only be paired with z-planes. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the plane. If not specified, the name will be the empty string. surface_id : int, optional @@ -905,6 +963,8 @@ class ZPlane(PlaneMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight periodic_surface : openmc.Surface If a periodic boundary condition is used, the surface with which this one is periodic with @@ -1076,7 +1136,8 @@ class QuadricMixin: else: base_cls = type(tsurf)._virtual_base # Copy necessary surface attributes to new kwargs dictionary - kwargs = {'boundary_type': tsurf.boundary_type, 'name': tsurf.name} + kwargs = {'boundary_type': tsurf.boundary_type, + 'albedo': tsurf.albedo, 'name': tsurf.name} if inplace: kwargs['surface_id'] = tsurf.id kwargs.update({k: getattr(tsurf, k) for k in base_cls._coeff_keys}) @@ -1133,6 +1194,10 @@ class Cylinder(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the cylinder. If not specified, the name will be the empty string. @@ -1159,6 +1224,8 @@ class Cylinder(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -1280,8 +1347,8 @@ class Cylinder(QuadricMixin, Surface): # since the C++ layer doesn't support Cylinders right now with catch_warnings(): simplefilter('ignore', IDWarning) - kwargs = {'boundary_type': self.boundary_type, 'name': self.name, - 'surface_id': self.id} + kwargs = {'boundary_type': self.boundary_type, 'albedo': self.albedo, + 'name': self.name, 'surface_id': self.id} quad_rep = Quadric(*self._get_base_coeffs(), **kwargs) return quad_rep.to_xml_element() @@ -1302,6 +1369,10 @@ class XCylinder(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the cylinder. If not specified, the name will be the empty string. @@ -1320,6 +1391,8 @@ class XCylinder(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -1394,6 +1467,10 @@ class YCylinder(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the cylinder. If not specified, the name will be the empty string. @@ -1412,6 +1489,8 @@ class YCylinder(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -1486,6 +1565,10 @@ class ZCylinder(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the cylinder. If not specified, the name will be the empty string. @@ -1504,6 +1587,8 @@ class ZCylinder(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -1579,6 +1664,10 @@ class Sphere(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the sphere. If not specified, the name will be the empty string. surface_id : int, optional @@ -1598,6 +1687,8 @@ class Sphere(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -1683,6 +1774,11 @@ class Cone(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. + name : str Name of the cone. If not specified, the name will be the empty string. @@ -1705,6 +1801,8 @@ class Cone(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -1790,7 +1888,9 @@ class Cone(QuadricMixin, Surface): # since the C++ layer doesn't support Cones right now with catch_warnings(): simplefilter('ignore', IDWarning) - kwargs = {'boundary_type': self.boundary_type, 'name': self.name, + kwargs = {'boundary_type': self.boundary_type, + 'albedo': self.albedo, + 'name': self.name, 'surface_id': self.id} quad_rep = Quadric(*self._get_base_coeffs(), **kwargs) return quad_rep.to_xml_element() @@ -1814,6 +1914,10 @@ class XCone(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the cone. If not specified, the name will be the empty string. surface_id : int, optional @@ -1833,6 +1937,8 @@ class XCone(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -1903,6 +2009,10 @@ class YCone(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the cone. If not specified, the name will be the empty string. surface_id : int, optional @@ -1922,6 +2032,8 @@ class YCone(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -1992,6 +2104,10 @@ class ZCone(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the cone. If not specified, the name will be the empty string. surface_id : int, optional @@ -2011,6 +2127,8 @@ class ZCone(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -2075,6 +2193,10 @@ class Quadric(QuadricMixin, Surface): Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. + albedo : float, optional + Albedo of the surfaces as a ratio of particle weight after interaction + with the surface to the initial weight. Values must be positive. Only + applicable if the boundary type is 'reflective', 'periodic', or 'white'. name : str, optional Name of the surface. If not specified, the name will be the empty string. surface_id : int, optional @@ -2088,6 +2210,8 @@ class Quadric(QuadricMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -2183,7 +2307,9 @@ class TorusMixin: # Create rotated torus kwargs = { - 'boundary_type': surf.boundary_type, 'name': surf.name, + 'boundary_type': surf.boundary_type, + 'albedo': surf.albedo, + 'name': surf.name, 'a': surf.a, 'b': surf.b, 'c': surf.c } if inplace: @@ -2236,6 +2362,8 @@ class XTorus(TorusMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -2309,6 +2437,8 @@ class YTorus(TorusMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int @@ -2382,6 +2512,8 @@ class ZTorus(TorusMixin, Surface): boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. + albedo : float + Boundary albedo as a positive multiplier of particle weight coefficients : dict Dictionary of surface coefficients id : int diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 5b842399e8..42a8e27a16 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -28,6 +28,9 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const Direction u = surf.reflect(p.r(), p.u(), &p); u /= u.norm(); + // Handle the effects of the surface albedo on the particle's weight. + BoundaryCondition::handle_albedo(p, surf); + p.cross_reflective_bc(surf, u); } @@ -40,6 +43,9 @@ void WhiteBC::handle_particle(Particle& p, const Surface& surf) const Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed()); u /= u.norm(); + // Handle the effects of the surface albedo on the particle's weight. + BoundaryCondition::handle_albedo(p, surf); + p.cross_reflective_bc(surf, u); } @@ -130,6 +136,9 @@ void TranslationalPeriodicBC::handle_particle( "hitting a surface, but that surface is not recognized by the BC."); } + // Handle the effects of the surface albedo on the particle's weight. + BoundaryCondition::handle_albedo(p, surf); + // Pass the new location and surface to the particle. p.cross_periodic_bc(surf, new_r, p.u(), new_surface); } @@ -264,6 +273,9 @@ void RotationalPeriodicBC::handle_particle( Direction new_u = { cos_theta * u.x - sin_theta * u.y, sin_theta * u.x + cos_theta * u.y, u.z}; + // Handle the effects of the surface albedo on the particle's weight. + BoundaryCondition::handle_albedo(p, surf); + // Pass the new location, direction, and surface to the particle. p.cross_periodic_bc(surf, new_r, new_u, new_surface); } diff --git a/src/dagmc.cpp b/src/dagmc.cpp index e245d4475d..4ffe4e6526 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -289,10 +289,10 @@ void DAGUniverse::init_geometry() bc_value == "transmission") { // set to transmission by default (nullptr) } else if (bc_value == "vacuum") { - s->bc_ = std::make_shared(); + s->bc_ = make_unique(); } else if (bc_value == "reflective" || bc_value == "reflect" || bc_value == "reflecting") { - s->bc_ = std::make_shared(); + s->bc_ = make_unique(); } else if (bc_value == "periodic") { fatal_error("Periodic boundary condition not supported in DAGMC."); } else { @@ -310,7 +310,7 @@ void DAGUniverse::init_geometry() // if this surface belongs to the graveyard if (graveyard && parent_vols.find(graveyard) != parent_vols.end()) { // set graveyard surface BC's to vacuum - s->bc_ = std::make_shared(); + s->bc_ = make_unique(); } // add to global array and map diff --git a/src/surface.cpp b/src/surface.cpp index 01e2f275eb..394a4b3663 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -80,19 +80,36 @@ Surface::Surface(pugi::xml_node surf_node) if (surf_bc == "transmission" || surf_bc == "transmit" || surf_bc.empty()) { // Leave the bc_ a nullptr } else if (surf_bc == "vacuum") { - bc_ = std::make_shared(); + bc_ = make_unique(); } else if (surf_bc == "reflective" || surf_bc == "reflect" || surf_bc == "reflecting") { - bc_ = std::make_shared(); + bc_ = make_unique(); } else if (surf_bc == "white") { - bc_ = std::make_shared(); + bc_ = make_unique(); } else if (surf_bc == "periodic") { - // periodic BC's are handled separately + // Periodic BCs are handled separately } else { fatal_error(fmt::format("Unknown boundary condition \"{}\" specified " "on surface {}", surf_bc, id_)); } + + if (check_for_node(surf_node, "albedo") && bc_) { + double surf_alb = std::stod(get_node_value(surf_node, "albedo")); + + if (surf_alb < 0.0) + fatal_error(fmt::format("Surface {} has an albedo of {}. " + "Albedo values must be positive.", + id_, surf_alb)); + + if (surf_alb > 1.0) + warning(fmt::format("Surface {} has an albedo of {}. " + "Albedos greater than 1 may cause " + "unphysical behaviour.", + id_, surf_alb)); + + bc_->set_albedo(surf_alb); + } } } @@ -154,6 +171,7 @@ void Surface::to_hdf5(hid_t group_id) const if (bc_) { write_string(surf_group, "boundary_type", bc_->type(), false); + bc_->to_hdf5(surf_group); } else { write_string(surf_group, "boundary_type", "transmission", false); } @@ -1156,9 +1174,10 @@ void read_surfaces(pugi::xml_node node) } // Loop over XML surface elements and populate the array. Keep track of - // periodic surfaces. + // periodic surfaces and their albedos. model::surfaces.reserve(n_surfaces); std::set> periodic_pairs; + std::unordered_map albedo_map; { pugi::xml_node surf_node; int i_surf; @@ -1221,6 +1240,12 @@ void read_surfaces(pugi::xml_node node) if (check_for_node(surf_node, "boundary")) { std::string surf_bc = get_node_value(surf_node, "boundary", true, true); if (surf_bc == "periodic") { + // Check for surface albedo. Skip sanity check as it is already done + // in the Surface class's constructor. + if (check_for_node(surf_node, "albedo")) { + albedo_map[model::surfaces.back()->id_] = + std::stod(get_node_value(surf_node, "albedo")); + } if (check_for_node(surf_node, "periodic_surface_id")) { int i_periodic = std::stoi(get_node_value(surf_node, "periodic_surface_id")); @@ -1284,7 +1309,7 @@ void read_surfaces(pugi::xml_node node) periodic_pairs.erase(second_unresolved); } - // Assign the periodic boundary conditions + // Assign the periodic boundary conditions with albedos for (auto periodic_pair : periodic_pairs) { int i_surf = model::surface_map[periodic_pair.first]; int j_surf = model::surface_map[periodic_pair.second]; @@ -1302,11 +1327,19 @@ void read_surfaces(pugi::xml_node node) // planes are parallel which indicates a translational periodic boundary // condition. Otherwise, it is a rotational periodic BC. if (std::abs(1.0 - dot_prod) < FP_PRECISION) { - surf1.bc_ = std::make_shared(i_surf, j_surf); - surf2.bc_ = surf1.bc_; + surf1.bc_ = make_unique(i_surf, j_surf); + surf2.bc_ = make_unique(i_surf, j_surf); } else { - surf1.bc_ = std::make_shared(i_surf, j_surf); - surf2.bc_ = surf1.bc_; + surf1.bc_ = make_unique(i_surf, j_surf); + surf2.bc_ = make_unique(i_surf, j_surf); + } + + // If albedo data is present in albedo map, set the boundary albedo. + if (albedo_map.count(surf1.id_)) { + surf1.bc_->set_albedo(albedo_map[surf1.id_]); + } + if (albedo_map.count(surf2.id_)) { + surf2.bc_->set_albedo(albedo_map[surf2.id_]); } } } diff --git a/tests/regression_tests/albedo_box/__init__.py b/tests/regression_tests/albedo_box/__init__.py new file mode 100755 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/albedo_box/geometry.xml b/tests/regression_tests/albedo_box/geometry.xml new file mode 100644 index 0000000000..7d0e9b5f5d --- /dev/null +++ b/tests/regression_tests/albedo_box/geometry.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/tests/regression_tests/albedo_box/materials.xml b/tests/regression_tests/albedo_box/materials.xml new file mode 100644 index 0000000000..2472a74717 --- /dev/null +++ b/tests/regression_tests/albedo_box/materials.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/regression_tests/albedo_box/results_true.dat b/tests/regression_tests/albedo_box/results_true.dat new file mode 100644 index 0000000000..0f571f2459 --- /dev/null +++ b/tests/regression_tests/albedo_box/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.590800E+00 4.251788E-03 diff --git a/tests/regression_tests/albedo_box/settings.xml b/tests/regression_tests/albedo_box/settings.xml new file mode 100644 index 0000000000..cc66683057 --- /dev/null +++ b/tests/regression_tests/albedo_box/settings.xml @@ -0,0 +1,9 @@ + + + + eigenvalue + 10 + 5 + 1000 + + diff --git a/tests/regression_tests/albedo_box/test.py b/tests/regression_tests/albedo_box/test.py new file mode 100755 index 0000000000..179f58e5b3 --- /dev/null +++ b/tests/regression_tests/albedo_box/test.py @@ -0,0 +1,6 @@ +from tests.testing_harness import TestHarness + + +def test_albedo_box(): + harness = TestHarness('statepoint.10.h5') + harness.main()