crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

add white boundary condition

Browse source
This commit is contained in:
rockfool 2019-09-09 16:22:25 -04:00
parent 3040d4d3bb
commit 5f394b2816
4 changed files with 141 additions and 29 deletions
Download patch file Download diff file Expand all files Collapse all files

3
include/openmc/surface.h
View file

@ -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

56
openmc/surface.py
View file

@ -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

50
src/particle.cpp
View file

@ -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

61
src/surface.cpp
View file

@ -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;