From 5f394b281614fda5a9aa5a3695db4513c53e5cb7 Mon Sep 17 00:00:00 2001 From: rockfool Date: Mon, 9 Sep 2019 16:22:25 -0400 Subject: [PATCH] add white boundary condition --- include/openmc/surface.h | 3 ++ openmc/surface.py | 56 ++++++++++++++++++------------------ src/particle.cpp | 50 ++++++++++++++++++++++++++++++++ src/surface.cpp | 61 +++++++++++++++++++++++++++++++++++++++- 4 files changed, 141 insertions(+), 29 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index d61ae07cef..5272dab3fb 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -25,6 +25,7 @@ extern "C" const int BC_TRANSMIT; extern "C" const int BC_VACUUM; extern "C" const int BC_REFLECT; extern "C" const int BC_PERIODIC; +extern "C" const int BC_WHITE; //============================================================================== // Global variables @@ -115,6 +116,8 @@ public: //! \return Outgoing direction of the ray virtual Direction reflect(Position r, Direction u) const; + virtual Direction diffuse_reflect(Position r, Direction u) const; + //! Evaluate the equation describing the surface. //! //! Surfaces can be described by some function f(x, y, z) = 0. This member diff --git a/openmc/surface.py b/openmc/surface.py index dc251485d2..ec5b833048 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -12,7 +12,7 @@ from openmc.region import Region, Intersection, Union from openmc.mixin import IDManagerMixin -_BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic'] +_BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic', 'white'] _WARNING_UPPER = """\ "{}(...) accepts an argument named '{}', not '{}'. Future versions of OpenMC \ @@ -32,7 +32,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. Note that periodic boundary conditions @@ -535,7 +535,7 @@ class XPlane(Plane): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional 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 @@ -549,7 +549,7 @@ class XPlane(Plane): ---------- x0 : float Location of the plane - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. periodic_surface : openmc.Surface @@ -660,7 +660,7 @@ class YPlane(Plane): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional 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 @@ -674,7 +674,7 @@ class YPlane(Plane): ---------- y0 : float Location of the plane - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. periodic_surface : openmc.Surface @@ -786,7 +786,7 @@ class ZPlane(Plane): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional 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 @@ -800,7 +800,7 @@ class ZPlane(Plane): ---------- z0 : float Location of the plane - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. periodic_surface : openmc.Surface @@ -912,7 +912,7 @@ class Cylinder(Surface): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -926,7 +926,7 @@ class Cylinder(Surface): ---------- r : float Radius of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -963,7 +963,7 @@ class XCylinder(Cylinder): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -983,7 +983,7 @@ class XCylinder(Cylinder): y-coordinate of the center of the cylinder z0 : float z-coordinate of the center of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1110,7 +1110,7 @@ class YCylinder(Cylinder): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1130,7 +1130,7 @@ class YCylinder(Cylinder): x-coordinate of the center of the cylinder z0 : float z-coordinate of the center of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1257,7 +1257,7 @@ class ZCylinder(Cylinder): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1277,7 +1277,7 @@ class ZCylinder(Cylinder): x-coordinate of the center of the cylinder y0 : float y-coordinate of the center of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1403,7 +1403,7 @@ class Sphere(Surface): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1428,7 +1428,7 @@ class Sphere(Surface): z-coordinate of the center of the sphere r : float Radius of the sphere - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1577,7 +1577,7 @@ class Cone(Surface): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1602,7 +1602,7 @@ class Cone(Surface): z-coordinate of the apex r2 : float Parameter related to the aperature - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1698,7 +1698,7 @@ class XCone(Cone): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1723,7 +1723,7 @@ class XCone(Cone): z-coordinate of the apex r2 : float Parameter related to the aperature - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1769,7 +1769,7 @@ class YCone(Cone): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1794,7 +1794,7 @@ class YCone(Cone): z-coordinate of the apex r2 : float Parameter related to the aperature - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1840,7 +1840,7 @@ class ZCone(Cone): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1865,7 +1865,7 @@ class ZCone(Cone): z-coordinate of the apex r2 : float Parameter related to the aperature - boundary_type : {'transmission, 'vacuum', 'reflective'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1911,7 +1911,7 @@ class Quadric(Surface): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1924,7 +1924,7 @@ class Quadric(Surface): ---------- a, b, c, d, e, f, g, h, j, k : float coefficients for the surface - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict diff --git a/src/particle.cpp b/src/particle.cpp index 1fb1eadb09..b6a3f5dee1 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -490,6 +490,56 @@ Particle::cross_surface() } return; + } else if (surf->bc_ == BC_WHITE && settings::run_mode != RUN_MODE_PLOTTING) { + // ======================================================================= + // WHITE BOUNDARY + + // Do not handle white boundary conditions on lower universes + if (n_coord_ != 1) { + this->mark_as_lost("Cannot diffuse reflect particle " + std::to_string(id_) + + " off surface in a lower universe."); + return; + } + + if (!model::active_surface_tallies.empty()) { + score_surface_tally(this, model::active_surface_tallies); + } + + if (!model::active_meshsurf_tallies.empty()) { + Position r {this->r()}; + this->r() -= TINY_BIT * this->u(); + score_surface_tally(this, model::active_meshsurf_tallies); + this->r() = r; + } + + // Diffuse Reflect particle off surface + Direction u = surf->diffuse_reflect(this->r(), this->u()); + + // Make sure new particle direction is normalized + this->u() = u / u.norm(); + + // Reassign particle's cell and surface + coord_[0].cell = cell_last_[n_coord_last_ - 1]; + surface_ = -surface_; + + if (!settings::dagmc) { + n_coord_ = 1; + if (!find_cell(this, true)) { + this->mark_as_lost("Couldn't find particle after diffuse reflecting from surface " + + std::to_string(surf->id_) + "."); + return; + } + } + + // Set previous coordinate going slightly past surface crossing + // r_last_current_ = this->r() + TINY_BIT*this->u(); + + // Diagnostic message + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Diffuse reflected from surface " + std::to_string(surf->id_)); + } + return; + } else if (surf->bc_ == BC_PERIODIC && settings::run_mode != RUN_MODE_PLOTTING) { // ======================================================================= // PERIODIC BOUNDARY diff --git a/src/surface.cpp b/src/surface.cpp index 14d077cfff..f759e7752b 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -11,6 +11,7 @@ #include "openmc/settings.h" #include "openmc/string_utils.h" #include "openmc/xml_interface.h" +#include "openmc/random_lcg.h" namespace openmc { @@ -22,7 +23,7 @@ extern "C" const int BC_TRANSMIT {0}; extern "C" const int BC_VACUUM {1}; extern "C" const int BC_REFLECT {2}; extern "C" const int BC_PERIODIC {3}; - +extern "C" const int BC_WHITE {4}; //============================================================================== // Global variables //============================================================================== @@ -147,6 +148,8 @@ Surface::Surface(pugi::xml_node surf_node) } else if (surf_bc == "reflective" || surf_bc == "reflect" || surf_bc == "reflecting") { bc_ = BC_REFLECT; + } else if (surf_bc == "white" || surf_bc == "diffuse") { + bc_ = BC_WHITE; } else if (surf_bc == "periodic") { bc_ = BC_PERIODIC; } else { @@ -192,6 +195,59 @@ Surface::reflect(Position r, Direction u) const return u -= (2.0 * projection / magnitude) * n; } +Direction +Surface::diffuse_reflect(Position r, Direction u) const +{ + // Diffuse reflect direction according to the normal. + // cosine distribution + + Direction n = normal(r); + n = n/n.norm(); + const double projection = std::max(-1.0, std::min(1.0, n.dot(u))); + + // sample from inverse function, u=sqrt(rand) since p(u)=2u, so F(u)=u^2 + const double cosine = (projection>=0.0) ? + -std::sqrt(prn()) : std::sqrt(prn()); + + // handle special case for very large cosine + if(std::fabs(cosine)>1.0001){ + std::cout<<"cosine="<<"\n"; + exit(-1); + }else if(std::fabs(cosine)>=1.000){ + return u = cosine * n; + } + + // sample azimuthal distribution uniformly + double t1, t2, rr, ss, tt; + for(;;){ + t1 = 2.0 * prn() - 1.0; + t2 = 2.0 * prn() - 1.0; + rr = t1 * t1 + t2 * t2; + if(rr<=1.0) break; + } + rr = std::sqrt((1.0 - cosine * cosine)/rr); + t1 = t1 * rr; + t2 = t2 * rr; + + // determine the direction off surface relative to norm + if (std::fabs(n[2]<=0.9)) { // for the general normal + tt = (std::sqrt(n[0]*n[0]+n[1]*n[1])); + ss = 1.0/tt; + u[0] = n[0] * cosine + (t1*n[0]*n[2]-t2*n[1])*ss; + u[1] = n[1] * cosine + (t1*n[1]*n[2]+t2*n[0])*ss; + u[2] = n[2] * cosine - t1*tt; + } else { // for the normal is almost along with z-axis + tt = (std::sqrt(n[0]*n[0]+n[2]*n[2])); + ss = 1.0/tt; + u[0] = n[0] * cosine + (t1*n[0]*n[1]+t2*n[2])*ss; + u[1] = n[1] * cosine - t1*tt; + u[2] = n[2] * cosine + (t1*n[1]*n[2]-t2*n[0])*ss; + } + + // normalize the direction + return u = u/u.norm(); +} + CSGSurface::CSGSurface() : Surface{} {}; CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface{surf_node} {}; @@ -213,6 +269,9 @@ CSGSurface::to_hdf5(hid_t group_id) const case BC_REFLECT : write_string(surf_group, "boundary_type", "reflective", false); break; + case BC_WHITE : + write_string(surf_group, "boundary_type", "white", false); + break; case BC_PERIODIC : write_string(surf_group, "boundary_type", "periodic", false); break;