mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 06:25:30 -04:00
Fix a bug in rotational periodic boundary conditions (#3692)
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
GuySten
GitHub
parent
818fd11b18
commit
c7d7fa4613
16 changed files with 297 additions and 73 deletions
|
|
@ -152,6 +152,8 @@ public:
|
|||
protected:
|
||||
//! Angle about the axis by which particle coordinates will be rotated
|
||||
double angle_;
|
||||
//! Do we need to flip surfaces senses when applying the transformation?
|
||||
bool flip_sense_;
|
||||
//! Ensure that choice of axes is right handed. axis_1_idx_ corresponds to the
|
||||
//! independent axis and axis_2_idx_ corresponds to the dependent axis in the
|
||||
//! 2D plane perpendicular to the planes' axis of rotation
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@
|
|||
#define OPENMC_SURFACE_H
|
||||
|
||||
#include <limits> // For numeric_limits
|
||||
#include <set>
|
||||
#include <string>
|
||||
#include <unordered_map>
|
||||
|
||||
|
|
@ -378,7 +379,39 @@ public:
|
|||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
void read_surfaces(pugi::xml_node node);
|
||||
//! Read surface definitions from XML and populate the global surfaces vector.
|
||||
//!
|
||||
//! This function parses surface elements from the XML input, creates the
|
||||
//! appropriate surface objects, and identifies periodic surfaces along with
|
||||
//! their albedo values and sense information.
|
||||
//!
|
||||
//! \param node XML node containing surface definitions
|
||||
//! \param[out] periodic_pairs Set of surface ID pairs representing periodic
|
||||
//! boundary conditions
|
||||
//! \param[out] albedo_map Map of surface IDs to albedo values for periodic
|
||||
//! surfaces
|
||||
//! \param[out] periodic_sense_map Map of surface IDs to their sense values
|
||||
//! (used to determine orientation for periodic BCs)
|
||||
void read_surfaces(pugi::xml_node node,
|
||||
std::set<std::pair<int, int>>& periodic_pairs,
|
||||
std::unordered_map<int, double>& albedo_map,
|
||||
std::unordered_map<int, int>& periodic_sense_map);
|
||||
|
||||
//! Resolve periodic surface pairs and assign boundary conditions.
|
||||
//!
|
||||
//! This function completes the setup of periodic boundary conditions by
|
||||
//! resolving unpaired periodic surfaces, determining whether each pair
|
||||
//! represents translational or rotational periodicity based on surface
|
||||
//! normals, and assigning the appropriate boundary condition objects.
|
||||
//!
|
||||
//! \param[inout] periodic_pairs Set of surface ID pairs representing periodic
|
||||
//! boundary conditions; unpaired entries are resolved
|
||||
//! \param albedo_map Map of surface IDs to albedo values for periodic surfaces
|
||||
//! \param periodic_sense_map Map of surface IDs to their sense values (used to
|
||||
//! determine orientation for periodic BCs)
|
||||
void prepare_boundary_conditions(std::set<std::pair<int, int>>& periodic_pairs,
|
||||
std::unordered_map<int, double>& albedo_map,
|
||||
std::unordered_map<int, int>& periodic_sense_map);
|
||||
|
||||
void free_memory_surfaces();
|
||||
|
||||
|
|
|
|||
|
|
@ -160,7 +160,7 @@ void TranslationalPeriodicBC::handle_particle(
|
|||
|
||||
RotationalPeriodicBC::RotationalPeriodicBC(
|
||||
int i_surf, int j_surf, PeriodicAxis axis)
|
||||
: PeriodicBC(i_surf, j_surf)
|
||||
: PeriodicBC(std::abs(i_surf) - 1, std::abs(j_surf) - 1)
|
||||
{
|
||||
Surface& surf1 {*model::surfaces[i_surf_]};
|
||||
Surface& surf2 {*model::surfaces[j_surf_]};
|
||||
|
|
@ -173,14 +173,9 @@ RotationalPeriodicBC::RotationalPeriodicBC(
|
|||
axis_2_idx_ = 2; // z component dependent
|
||||
break;
|
||||
case y:
|
||||
// for a right handed coordinate system, z should be the independent axis
|
||||
// but this would cause the y-rotation case to be different than the other
|
||||
// two. using a left handed coordinate system and a negative rotation the
|
||||
// compute angle and rotation matrix behavior mimics that of the x and z
|
||||
// cases
|
||||
zero_axis_idx_ = 1; // y component of plane must be zero
|
||||
axis_1_idx_ = 0; // x component independent
|
||||
axis_2_idx_ = 2; // z component dependent
|
||||
axis_1_idx_ = 2; // z component independent
|
||||
axis_2_idx_ = 0; // x component dependent
|
||||
break;
|
||||
case z:
|
||||
zero_axis_idx_ = 2; // z component of plane must be zero
|
||||
|
|
@ -192,10 +187,16 @@ RotationalPeriodicBC::RotationalPeriodicBC(
|
|||
fmt::format("You've specified an axis that is not x, y, or z."));
|
||||
}
|
||||
|
||||
Direction ax = {0.0, 0.0, 0.0};
|
||||
ax[zero_axis_idx_] = 1.0;
|
||||
|
||||
auto i_sign = std::copysign(1, i_surf);
|
||||
auto j_sign = -std::copysign(1, j_surf);
|
||||
|
||||
// Compute the surface normal vectors and make sure they are perpendicular
|
||||
// to the correct axis
|
||||
Direction norm1 = surf1.normal({0, 0, 0});
|
||||
Direction norm2 = surf2.normal({0, 0, 0});
|
||||
Direction norm1 = i_sign * surf1.normal({0, 0, 0});
|
||||
Direction norm2 = j_sign * surf2.normal({0, 0, 0});
|
||||
// Make sure both surfaces intersect the origin
|
||||
if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {
|
||||
throw std::invalid_argument(fmt::format(
|
||||
|
|
@ -212,8 +213,15 @@ RotationalPeriodicBC::RotationalPeriodicBC(
|
|||
surf2.id_));
|
||||
}
|
||||
|
||||
angle_ = compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_],
|
||||
norm2[axis_2_idx_], norm2[axis_1_idx_]);
|
||||
// Compute the signed rotation angle about the periodic axis. Note that
|
||||
// (n1×n2)·a = |n1||n2|sin(θ) and n1·n2 = |n1||n2|cos(θ), where a is the axis
|
||||
// of rotation.
|
||||
auto c = norm1.cross(norm2);
|
||||
angle_ = std::atan2(c.dot(ax), norm1.dot(norm2));
|
||||
|
||||
// If the normals point in the same general direction, the surface sense
|
||||
// should change when crossing the boundary
|
||||
flip_sense_ = (i_sign * j_sign > 0.0);
|
||||
|
||||
// Warn the user if the angle does not evenly divide a circle
|
||||
double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
|
||||
|
|
@ -225,47 +233,18 @@ RotationalPeriodicBC::RotationalPeriodicBC(
|
|||
}
|
||||
}
|
||||
|
||||
double RotationalPeriodicBC::compute_periodic_rotation(
|
||||
double rise_1, double run_1, double rise_2, double run_2) const
|
||||
{
|
||||
// Compute the BC rotation angle. Here it is assumed that both surface
|
||||
// normal vectors point inwards---towards the valid geometry region.
|
||||
// Consequently, the rotation angle is not the difference between the two
|
||||
// normals, but is instead the difference between one normal and one
|
||||
// anti-normal. (An incident ray on one surface must be an outgoing ray on
|
||||
// the other surface after rotation hence the anti-normal.)
|
||||
double theta1 = std::atan2(rise_1, run_1);
|
||||
double theta2 = std::atan2(rise_2, run_2) + PI;
|
||||
return theta2 - theta1;
|
||||
}
|
||||
|
||||
void RotationalPeriodicBC::handle_particle(
|
||||
Particle& p, const Surface& surf) const
|
||||
{
|
||||
int i_particle_surf = p.surface_index();
|
||||
int new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1;
|
||||
if (flip_sense_)
|
||||
new_surface = -new_surface;
|
||||
|
||||
// Figure out which of the two BC surfaces were struck to figure out if a
|
||||
// forward or backward rotation is required. Specify the other surface as
|
||||
// the particle's new surface.
|
||||
double theta;
|
||||
int new_surface;
|
||||
if (i_particle_surf == i_surf_) {
|
||||
theta = angle_;
|
||||
new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1;
|
||||
} else if (i_particle_surf == j_surf_) {
|
||||
theta = -angle_;
|
||||
new_surface = p.surface() > 0 ? -(i_surf_ + 1) : i_surf_ + 1;
|
||||
} else {
|
||||
throw std::runtime_error(
|
||||
"Called BoundaryCondition::handle_particle after "
|
||||
"hitting a surface, but that surface is not recognized by the BC.");
|
||||
}
|
||||
|
||||
// Rotate the particle's position and direction about the z-axis.
|
||||
// Rotate the particle's position and direction.
|
||||
Position r = p.r();
|
||||
Direction u = p.u();
|
||||
double cos_theta = std::cos(theta);
|
||||
double sin_theta = std::sin(theta);
|
||||
double cos_theta = std::cos(angle_);
|
||||
double sin_theta = std::sin(angle_);
|
||||
|
||||
Position new_r;
|
||||
new_r[zero_axis_idx_] = r[zero_axis_idx_];
|
||||
|
|
|
|||
|
|
@ -55,8 +55,13 @@ void read_geometry_xml()
|
|||
void read_geometry_xml(pugi::xml_node root)
|
||||
{
|
||||
// Read surfaces, cells, lattice
|
||||
read_surfaces(root);
|
||||
std::set<std::pair<int, int>> periodic_pairs;
|
||||
std::unordered_map<int, double> albedo_map;
|
||||
std::unordered_map<int, int> periodic_sense_map;
|
||||
|
||||
read_surfaces(root, periodic_pairs, albedo_map, periodic_sense_map);
|
||||
read_cells(root);
|
||||
prepare_boundary_conditions(periodic_pairs, albedo_map, periodic_sense_map);
|
||||
read_lattices(root);
|
||||
|
||||
// Check to make sure a boundary condition was applied to at least one
|
||||
|
|
|
|||
|
|
@ -744,9 +744,7 @@ void Particle::cross_periodic_bc(
|
|||
if (!neighbor_list_find_cell(*this)) {
|
||||
mark_as_lost("Couldn't find particle after hitting periodic "
|
||||
"boundary on surface " +
|
||||
std::to_string(surf.id_) +
|
||||
". The normal vector "
|
||||
"of one periodic surface may need to be reversed.");
|
||||
std::to_string(surf.id_) + ".");
|
||||
return;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@
|
|||
#include <fmt/core.h>
|
||||
|
||||
#include "openmc/array.h"
|
||||
#include "openmc/cell.h"
|
||||
#include "openmc/container_util.h"
|
||||
#include "openmc/error.h"
|
||||
#include "openmc/external/quartic_solver.h"
|
||||
|
|
@ -1169,7 +1170,10 @@ Direction SurfaceZTorus::normal(Position r) const
|
|||
|
||||
//==============================================================================
|
||||
|
||||
void read_surfaces(pugi::xml_node node)
|
||||
void read_surfaces(pugi::xml_node node,
|
||||
std::set<std::pair<int, int>>& periodic_pairs,
|
||||
std::unordered_map<int, double>& albedo_map,
|
||||
std::unordered_map<int, int>& periodic_sense_map)
|
||||
{
|
||||
// Count the number of surfaces
|
||||
int n_surfaces = 0;
|
||||
|
|
@ -1180,8 +1184,6 @@ void read_surfaces(pugi::xml_node node)
|
|||
// Loop over XML surface elements and populate the array. Keep track of
|
||||
// periodic surfaces and their albedos.
|
||||
model::surfaces.reserve(n_surfaces);
|
||||
std::set<std::pair<int, int>> periodic_pairs;
|
||||
std::unordered_map<int, double> albedo_map;
|
||||
{
|
||||
pugi::xml_node surf_node;
|
||||
int i_surf;
|
||||
|
|
@ -1244,6 +1246,7 @@ 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") {
|
||||
periodic_sense_map[model::surfaces.back()->id_] = 0;
|
||||
// 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")) {
|
||||
|
|
@ -1275,6 +1278,28 @@ void read_surfaces(pugi::xml_node node)
|
|||
fmt::format("Two or more surfaces use the same unique ID: {}", id));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void prepare_boundary_conditions(std::set<std::pair<int, int>>& periodic_pairs,
|
||||
std::unordered_map<int, double>& albedo_map,
|
||||
std::unordered_map<int, int>& periodic_sense_map)
|
||||
{
|
||||
// Fill the senses map for periodic surfaces
|
||||
auto n_periodic = periodic_sense_map.size();
|
||||
for (const auto& cell : model::cells) {
|
||||
if (n_periodic == 0)
|
||||
break; // Early exit once all periodic surfaces found
|
||||
|
||||
for (auto s : cell->surfaces()) {
|
||||
auto surf_idx = std::abs(s) - 1;
|
||||
auto id = model::surfaces[surf_idx]->id_;
|
||||
|
||||
if (periodic_sense_map.count(id)) {
|
||||
periodic_sense_map[id] = std::copysign(1, s);
|
||||
--n_periodic;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Resolve unpaired periodic surfaces. A lambda function is used with
|
||||
// std::find_if to identify the unpaired surfaces.
|
||||
|
|
@ -1370,8 +1395,12 @@ void read_surfaces(pugi::xml_node node)
|
|||
"indicates that the two planes are not periodic about the X, Y, or Z "
|
||||
"axis, which is not supported."));
|
||||
}
|
||||
surf1.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf, axis);
|
||||
surf2.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf, axis);
|
||||
auto i_sign = periodic_sense_map[periodic_pair.first];
|
||||
auto j_sign = periodic_sense_map[periodic_pair.second];
|
||||
surf1.bc_ = make_unique<RotationalPeriodicBC>(
|
||||
i_sign * (i_surf + 1), j_sign * (j_surf + 1), axis);
|
||||
surf2.bc_ = make_unique<RotationalPeriodicBC>(
|
||||
j_sign * (j_surf + 1), i_sign * (i_surf + 1), axis);
|
||||
}
|
||||
|
||||
// If albedo data is present in albedo map, set the boundary albedo.
|
||||
|
|
|
|||
|
|
@ -0,0 +1,34 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<model>
|
||||
<materials>
|
||||
<material id="5">
|
||||
<density value="1.0" units="g/cc"/>
|
||||
<nuclide name="H1" ao="2.0"/>
|
||||
<nuclide name="O16" ao="1.0"/>
|
||||
<sab name="c_H_in_H2O"/>
|
||||
</material>
|
||||
<material id="6" depletable="true">
|
||||
<density value="4.5" units="g/cc"/>
|
||||
<nuclide name="U235" ao="1.0"/>
|
||||
</material>
|
||||
</materials>
|
||||
<geometry>
|
||||
<cell id="1" material="5" region="9 -10 -11 12" universe="0"/>
|
||||
<cell id="2" material="6" region="9 -10 -12" universe="0"/>
|
||||
<surface id="9" type="plane" boundary="periodic" coeffs="0.4999999999999999 0.8660254037844387 0.0 0.0"/>
|
||||
<surface id="10" type="plane" boundary="periodic" coeffs="-0.4999999999999999 0.8660254037844387 0.0 0.0"/>
|
||||
<surface id="11" type="x-plane" boundary="reflective" coeffs="5.0"/>
|
||||
<surface id="12" type="z-cylinder" coeffs="2.598076211353316 1.4999999999999998 2.0"/>
|
||||
</geometry>
|
||||
<settings>
|
||||
<run_mode>eigenvalue</run_mode>
|
||||
<particles>1000</particles>
|
||||
<batches>4</batches>
|
||||
<inactive>0</inactive>
|
||||
<source type="independent" strength="1.0" particle="neutron">
|
||||
<space type="box">
|
||||
<parameters>0 0 0 5 5 0</parameters>
|
||||
</space>
|
||||
</source>
|
||||
</settings>
|
||||
</model>
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
k-combined:
|
||||
1.848492E+00 2.933785E-03
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<model>
|
||||
<materials>
|
||||
<material id="3">
|
||||
<density value="1.0" units="g/cc"/>
|
||||
<nuclide name="H1" ao="2.0"/>
|
||||
<nuclide name="O16" ao="1.0"/>
|
||||
<sab name="c_H_in_H2O"/>
|
||||
</material>
|
||||
<material id="4" depletable="true">
|
||||
<density value="4.5" units="g/cc"/>
|
||||
<nuclide name="U235" ao="1.0"/>
|
||||
</material>
|
||||
</materials>
|
||||
<geometry>
|
||||
<cell id="1" material="3" region="-5 6 -7 8" universe="0"/>
|
||||
<cell id="2" material="4" region="-5 6 -8" universe="0"/>
|
||||
<surface id="5" type="plane" boundary="periodic" coeffs="-0.4999999999999999 -0.8660254037844387 0.0 0.0"/>
|
||||
<surface id="6" type="plane" boundary="periodic" coeffs="0.4999999999999999 -0.8660254037844387 0.0 0.0"/>
|
||||
<surface id="7" type="x-plane" boundary="reflective" coeffs="5.0"/>
|
||||
<surface id="8" type="z-cylinder" coeffs="2.598076211353316 1.4999999999999998 2.0"/>
|
||||
</geometry>
|
||||
<settings>
|
||||
<run_mode>eigenvalue</run_mode>
|
||||
<particles>1000</particles>
|
||||
<batches>4</batches>
|
||||
<inactive>0</inactive>
|
||||
<source type="independent" strength="1.0" particle="neutron">
|
||||
<space type="box">
|
||||
<parameters>0 0 0 5 5 0</parameters>
|
||||
</space>
|
||||
</source>
|
||||
</settings>
|
||||
</model>
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
k-combined:
|
||||
1.848492E+00 2.933785E-03
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<model>
|
||||
<materials>
|
||||
<material id="7">
|
||||
<density value="1.0" units="g/cc"/>
|
||||
<nuclide name="H1" ao="2.0"/>
|
||||
<nuclide name="O16" ao="1.0"/>
|
||||
<sab name="c_H_in_H2O"/>
|
||||
</material>
|
||||
<material id="8" depletable="true">
|
||||
<density value="4.5" units="g/cc"/>
|
||||
<nuclide name="U235" ao="1.0"/>
|
||||
</material>
|
||||
</materials>
|
||||
<geometry>
|
||||
<cell id="1" material="7" region="-13 -14 -15 16" universe="0"/>
|
||||
<cell id="2" material="8" region="-13 -14 -16" universe="0"/>
|
||||
<surface id="13" type="plane" boundary="periodic" coeffs="-0.4999999999999999 -0.8660254037844387 0.0 0.0"/>
|
||||
<surface id="14" type="plane" boundary="periodic" coeffs="-0.4999999999999999 0.8660254037844387 0.0 0.0"/>
|
||||
<surface id="15" type="x-plane" boundary="reflective" coeffs="5.0"/>
|
||||
<surface id="16" type="z-cylinder" coeffs="2.598076211353316 1.4999999999999998 2.0"/>
|
||||
</geometry>
|
||||
<settings>
|
||||
<run_mode>eigenvalue</run_mode>
|
||||
<particles>1000</particles>
|
||||
<batches>4</batches>
|
||||
<inactive>0</inactive>
|
||||
<source type="independent" strength="1.0" particle="neutron">
|
||||
<space type="box">
|
||||
<parameters>0 0 0 5 5 0</parameters>
|
||||
</space>
|
||||
</source>
|
||||
</settings>
|
||||
</model>
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
k-combined:
|
||||
1.848492E+00 2.933785E-03
|
||||
|
|
@ -1,14 +1,16 @@
|
|||
from math import sin, cos, pi
|
||||
|
||||
import openmc
|
||||
from openmc.utility_funcs import change_directory
|
||||
import pytest
|
||||
|
||||
from tests.testing_harness import PyAPITestHarness
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def model():
|
||||
model = openmc.model.Model()
|
||||
@pytest.mark.parametrize("flip1", [False, True])
|
||||
@pytest.mark.parametrize("flip2", [False, True])
|
||||
def test_periodic(flip1, flip2):
|
||||
model = openmc.Model()
|
||||
|
||||
# Define materials
|
||||
water = openmc.Material()
|
||||
|
|
@ -27,31 +29,52 @@ def model():
|
|||
# answers.
|
||||
theta1 = (-1/6 + 1/2) * pi
|
||||
theta2 = (1/6 - 1/2) * pi
|
||||
plane1 = openmc.Plane(a=cos(theta1), b=sin(theta1), boundary_type='periodic')
|
||||
plane2 = openmc.Plane(a=cos(theta2), b=sin(theta2), boundary_type='periodic')
|
||||
if flip1:
|
||||
plane1 = openmc.Plane(a=-cos(theta1), b=-sin(theta1), boundary_type='periodic')
|
||||
else:
|
||||
plane1 = openmc.Plane(a=cos(theta1), b=sin(theta1), boundary_type='periodic')
|
||||
if flip2:
|
||||
plane2 = openmc.Plane(a=-cos(theta2), b=-sin(theta2), boundary_type='periodic')
|
||||
else:
|
||||
plane2 = openmc.Plane(a=cos(theta2), b=sin(theta2), boundary_type='periodic')
|
||||
|
||||
x_max = openmc.XPlane(5., boundary_type='reflective')
|
||||
|
||||
z_cyl = openmc.ZCylinder(x0=3*cos(pi/6), y0=3*sin(pi/6), r=2.0)
|
||||
|
||||
outside_cyl = openmc.Cell(1, fill=water, region=(
|
||||
+plane1 & +plane2 & -x_max & +z_cyl))
|
||||
inside_cyl = openmc.Cell(2, fill=fuel, region=(
|
||||
+plane1 & +plane2 & -z_cyl))
|
||||
match (flip1, flip2):
|
||||
case (False, False):
|
||||
outside_cyl = openmc.Cell(1, fill=water, region=(
|
||||
+plane1 & +plane2 & -x_max & +z_cyl))
|
||||
inside_cyl = openmc.Cell(2, fill=fuel, region=(
|
||||
+plane1 & +plane2 & -z_cyl))
|
||||
case (False, True):
|
||||
outside_cyl = openmc.Cell(1, fill=water, region=(
|
||||
+plane1 & -plane2 & -x_max & +z_cyl))
|
||||
inside_cyl = openmc.Cell(2, fill=fuel, region=(
|
||||
+plane1 & -plane2 & -z_cyl))
|
||||
case (True, False):
|
||||
outside_cyl = openmc.Cell(1, fill=water, region=(
|
||||
-plane1 & +plane2 & -x_max & +z_cyl))
|
||||
inside_cyl = openmc.Cell(2, fill=fuel, region=(
|
||||
-plane1 & +plane2 & -z_cyl))
|
||||
case (True, True):
|
||||
outside_cyl = openmc.Cell(1, fill=water, region=(
|
||||
-plane1 & -plane2 & -x_max & +z_cyl))
|
||||
inside_cyl = openmc.Cell(2, fill=fuel, region=(
|
||||
-plane1 & -plane2 & -z_cyl))
|
||||
root_universe = openmc.Universe(0, cells=(outside_cyl, inside_cyl))
|
||||
model.geometry = openmc.Geometry(root_universe)
|
||||
|
||||
# Define settings
|
||||
model.settings = openmc.Settings()
|
||||
model.settings.particles = 1000
|
||||
model.settings.batches = 4
|
||||
model.settings.inactive = 0
|
||||
model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
|
||||
(0, 0, 0), (5, 5, 0))
|
||||
model.settings.source = openmc.IndependentSource(
|
||||
space=openmc.stats.Box((0, 0, 0), (5, 5, 0))
|
||||
)
|
||||
return model
|
||||
|
||||
with change_directory(f'{flip1}-{flip2}'):
|
||||
harness = PyAPITestHarness('statepoint.4.h5', model)
|
||||
harness.main()
|
||||
|
||||
def test_periodic(model):
|
||||
harness = PyAPITestHarness('statepoint.4.h5', model)
|
||||
harness.main()
|
||||
|
|
|
|||
47
tests/unit_tests/test_periodic_bc.py
Normal file
47
tests/unit_tests/test_periodic_bc.py
Normal file
|
|
@ -0,0 +1,47 @@
|
|||
from math import cos, sin, radians
|
||||
import random
|
||||
|
||||
import openmc
|
||||
import pytest
|
||||
|
||||
|
||||
@pytest.mark.parametrize("angle", [30., 45., 60., 90., 120.])
|
||||
def test_rotational_periodic_bc(angle):
|
||||
# Pick random starting angle
|
||||
start = random.uniform(0., 360.)
|
||||
degrees = angle
|
||||
ang1 = radians(start)
|
||||
ang2 = radians(start + degrees)
|
||||
|
||||
# Define three points on each plane and then randomly shuffle them
|
||||
p1_points = [(0., 0., 0.), (cos(ang1), sin(ang1), 0.), (0., 0., 1.)]
|
||||
p2_points = [(0., 0., 0.), (cos(ang2), sin(ang2), 0.), (0., 0., 1.)]
|
||||
random.shuffle(p1_points)
|
||||
random.shuffle(p2_points)
|
||||
|
||||
# Create periodic planes and a cylinder
|
||||
p1 = openmc.Plane.from_points(*p1_points, boundary_type='periodic')
|
||||
p2 = openmc.Plane.from_points(*p2_points, boundary_type='periodic')
|
||||
p1.periodic_surface = p2
|
||||
zcyl = openmc.ZCylinder(r=5., boundary_type='vacuum')
|
||||
|
||||
# Figure out which side of planes to use based on a point in the middle
|
||||
ang_mid = radians(start + degrees/2.)
|
||||
mid_point = (cos(ang_mid), sin(ang_mid), 0.)
|
||||
r1 = -p1 if mid_point in -p1 else +p1
|
||||
r2 = -p2 if mid_point in -p2 else +p2
|
||||
|
||||
# Create one cell bounded by the two planes and the cylinder
|
||||
mat = openmc.Material(density=1.0, density_units='g/cm3', components={'U235': 1.0})
|
||||
cell = openmc.Cell(fill=mat, region=r1 & r2 & -zcyl)
|
||||
|
||||
# Make the model complete
|
||||
model = openmc.Model()
|
||||
model.geometry = openmc.Geometry([cell])
|
||||
model.settings.source = openmc.IndependentSource(space=openmc.stats.Point(mid_point))
|
||||
model.settings.particles = 1000
|
||||
model.settings.batches = 10
|
||||
model.settings.inactive = 5
|
||||
|
||||
# Run the model
|
||||
model.run()
|
||||
Loading…
Add table
Add a link
Reference in a new issue