Skip to content

Commit

Permalink
Leaky and Albedo Boundary Conditions Implementation (#2724)
Browse files Browse the repository at this point in the history 
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
  • Loading branch information
@nutcasev15 @paulromano
nutcasev15 and paulromano committed Oct 30, 2023
1 parent 2c1e304 commit 693314d
Show file tree
Hide file tree
Showing 15 changed files with 377 additions and 83 deletions.
8 changes: 5 additions & 3 deletions docs/source/io_formats/summary.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'.

Expand Down
32 changes: 25 additions & 7 deletions docs/source/usersguide/geometry.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Expand All @@ -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:

-----
Expand Down
39 changes: 39 additions & 0 deletions include/openmc/boundary_condition.h
View file
Original file line number Diff line number Diff line change
@@ -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 <fmt/core.h>

namespace openmc {

Expand All @@ -22,8 +25,44 @@ class BoundaryCondition {
//! \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;
};

//==============================================================================
Expand Down
8 changes: 4 additions & 4 deletions include/openmc/surface.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,10 +83,10 @@ struct BoundingBox {

class Surface {
public:
int id_; //!< Unique ID
std::string name_; //!< User-defined name
std::shared_ptr<BoundaryCondition> 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<BoundaryCondition> 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);
Expand Down
93 changes: 56 additions & 37 deletions openmc/model/funcs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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) | \
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down

0 comments on commit 693314d

Please sign in to comment.