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

Merge pull request #1211 from pshriwise/hex_indices

Browse source 
Hex lattices indices update and lattice interior cell coincidence
This commit is contained in:
Sterling Harper 2019-04-03 23:49:08 -04:00 • committed by GitHub
commit 4f724b8905
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
7 changed files with 293 additions and 36 deletions
Download patch file Download diff file Expand all files Collapse all files

8
include/openmc/geometry.h
View file

@ -35,6 +35,14 @@ struct BoundaryInfo {
std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
};
//==============================================================================
//! Check two distances by coincidence tolerance
//==============================================================================
inline bool coincident(double d1, double d2) {
return std::abs(d1 - d2) < FP_COINCIDENT;
}
//==============================================================================
//! Check for overlapping cells at a particle's position.
//==============================================================================

2
src/geometry.cpp
View file

@ -396,7 +396,7 @@ BoundaryInfo distance_to_boundary(Particle* p)
// a higher level then we need to make sure that the higher level boundary
// is selected. This logic must consider floating point precision.
double& d = info.distance;
if (d_surf < d_lat) {
if (d_surf < d_lat - FP_COINCIDENT) {
if (d == INFINITY || (d - d_surf)/d >= FP_REL_PRECISION) {
d = d_surf;

86
src/lattice.cpp
View file

@ -6,6 +6,7 @@
#include "openmc/cell.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/string_utils.h"
@ -291,7 +292,7 @@ RectLattice::get_indices(Position r, Direction u) const
double ix_ {(r.x - lower_left_.x) / pitch_.x};
long ix_close {std::lround(ix_)};
int ix;
if (std::abs(ix_ - ix_close) < FP_COINCIDENT) {
if (coincident(ix_, ix_close)) {
ix = (u.x > 0) ? ix_close : ix_close - 1;
} else {
ix = std::floor(ix_);
@ -301,7 +302,7 @@ RectLattice::get_indices(Position r, Direction u) const
double iy_ {(r.y - lower_left_.y) / pitch_.y};
long iy_close {std::lround(iy_)};
int iy;
if (std::abs(iy_ - iy_close) < FP_COINCIDENT) {
if (coincident(iy_, iy_close)) {
iy = (u.y > 0) ? iy_close : iy_close - 1;
} else {
iy = std::floor(iy_);
@ -312,7 +313,7 @@ RectLattice::get_indices(Position r, Direction u) const
if (is_3d_) {
double iz_ {(r.z - lower_left_.z) / pitch_.z};
long iz_close {std::lround(iz_)};
if (std::abs(iz_ - iz_close) < FP_COINCIDENT) {
if (coincident(iz_, iz_close)) {
iz = (u.z > 0) ? iz_close : iz_close - 1;
} else {
iz = std::floor(iz_);
@ -710,69 +711,84 @@ const
std::array<int, 3>
HexLattice::get_indices(Position r, Direction u) const
{
// The implementation for HexLattice currently doesn't use direction
// information. As a result, we move the position slightly forward to
// determine what lattice index the particle is most likely to be in.
r += TINY_BIT * u;
// Offset the xyz by the lattice center.
Position r_o {r.x - center_.x, r.y - center_.y, r.z};
if (is_3d_) {r_o.z -= center_.z;}
// Index the z direction.
std::array<int, 3> out;
// Index the z direction, accounting for coincidence
int iz = 0;
if (is_3d_) {
out[2] = std::floor(r_o.z / pitch_[1] + 0.5 * n_axial_);
} else {
out[2] = 0;
double iz_ {r_o.z / pitch_[1] + 0.5 * n_axial_};
long iz_close {std::lround(iz_)};
if (coincident(iz_, iz_close)) {
iz = (u.z > 0) ? iz_close : iz_close - 1;
} else {
iz = std::floor(iz_);
}
}
// Convert coordinates into skewed bases. The (x, alpha) basis is used to
// find the index of the global coordinates to within 4 cells.
double alpha = r_o.y - r_o.x / std::sqrt(3.0);
out[0] = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0]));
out[1] = std::floor(alpha / pitch_[0]);
int ix = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0]));
int ia = std::floor(alpha / pitch_[0]);
// Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
// the array is offset so that the indices never go below 0).
out[0] += n_rings_-1;
out[1] += n_rings_-1;
ix += n_rings_-1;
ia += n_rings_-1;
// Calculate the (squared) distance between the particle and the centers of
// the four possible cells. Regular hexagonal tiles form a Voronoi
// tessellation so the xyz should be in the hexagonal cell that it is closest
// to the center of. This method is used over a method that uses the
// remainders of the floor divisions above because it provides better finite
// precision performance. Squared distances are used becasue they are more
// precision performance. Squared distances are used because they are more
// computationally efficient than normal distances.
int k {1};
int k_min {1};
// COINCIDENCE CHECK
// if a distance to center, d, is within the coincidence tolerance of the
// current minimum distance, d_min, the particle is on an edge or vertex.
// In this case, the dot product of the position vector and direction vector
// for the current indices, dp, and the dot product for the currently selected
// indices, dp_min, are compared. The cell which the particle is moving into
// is kept (i.e. the cell with the lowest dot product as the vectors will be
// completely opposed if the particle is moving directly toward the center of
// the cell).
int ix_chg {};
int ia_chg {};
double d_min {INFTY};
double dp_min {INFTY};
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 2; j++) {
const std::array<int, 3> i_xyz {out[0] + j, out[1] + i, 0};
// get local coordinates
const std::array<int, 3> i_xyz {ix + j, ia + i, 0};
Position r_t = get_local_position(r, i_xyz);
// calculate distance
double d = r_t.x*r_t.x + r_t.y*r_t.y;
if (d < d_min) {
// check for coincidence
bool on_boundary = coincident(d, d_min);
if (d < d_min || on_boundary) {
// normalize r_t and find dot product
r_t /= std::sqrt(d);
double dp = u.x * r_t.x + u.y * r_t.y;
// do not update values if particle is on a
// boundary and not moving into this cell
if (on_boundary && dp > dp_min) continue;
// update values
d_min = d;
k_min = k;
ix_chg = j;
ia_chg = i;
dp_min = dp;
}
k++;
}
}
// Select the minimum squared distance which corresponds to the cell the
// coordinates are in.
if (k_min == 2) {
++out[0];
} else if (k_min == 3) {
++out[1];
} else if (k_min == 4) {
++out[0];
++out[1];
}
// update outgoing indices
ix += ix_chg;
ia += ia_chg;
return out;
return {ix, ia, iz};
}
//==============================================================================

0
tests/regression_tests/lattice_hex_coincident/__init__.py Normal file
View file

81
tests/regression_tests/lattice_hex_coincident/inputs_true.dat Normal file
View file

@ -0,0 +1,81 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="13" material="13" region="-9 10 -11" universe="9" />
<cell id="14" material="14" region="9" universe="9" />
<cell id="15" material="16" name="coolant" region="-12 10 -11" universe="10" />
<cell id="16" material="17" name="zirconium_shell" region="12 -13 10 -11" universe="10" />
<cell id="17" material="15" name="lead_shell" region="13 -14 10 -11" universe="10" />
<cell id="18" material="14" name="matrix coolant surround" region="14 10 -11" universe="10" />
<cell id="19" material="14" universe="11" />
<cell fill="12" id="20" name="container cell" region="-15 16 -17 18 19 -20 10 -11" universe="13" />
<hex_lattice id="12" n_rings="2" name="regular fuel assembly">
<pitch>1.4</pitch>
<outer>11</outer>
<center>0.0 0.0</center>
<universes>
10
10 10
9
10 10
10</universes>
</hex_lattice>
<surface coeffs="0.0 0.0 0.7" id="9" type="z-cylinder" />
<surface boundary="reflective" coeffs="0.0" id="10" type="z-plane" />
<surface boundary="reflective" coeffs="10.0" id="11" type="z-plane" />
<surface coeffs="0.0 0.0 0.293" id="12" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.35" id="13" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.352" id="14" type="z-cylinder" />
<surface boundary="reflective" coeffs="1.4" id="15" type="y-plane" />
<surface boundary="reflective" coeffs="-1.4" id="16" type="y-plane" />
<surface boundary="reflective" coeffs="1.7320508075688772 1.0 0.0 2.8" id="17" type="plane" />
<surface boundary="reflective" coeffs="-1.7320508075688772 1.0 0.0 -2.8" id="18" type="plane" />
<surface boundary="reflective" coeffs="1.7320508075688772 1.0 0.0 -2.8" id="19" type="plane" />
<surface boundary="reflective" coeffs="-1.7320508075688772 1.0 0.0 2.8" id="20" type="plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="13">
<density units="sum" />
<nuclide ao="0.0049817" name="U235" />
</material>
<material id="14">
<density units="atom/b-cm" value="0.087742" />
<nuclide ao="1.0" name="C0" />
<sab name="c_Graphite" />
</material>
<material id="15" name="Lead">
<density units="g/cm3" value="10.32" />
<nuclide ao="0.014" name="Pb204" />
<nuclide ao="0.241" name="Pb206" />
<nuclide ao="0.221" name="Pb207" />
<nuclide ao="0.524" name="Pb208" />
</material>
<material id="16">
<density units="atom/b-cm" value="0.00054464" />
<nuclide ao="1.0" name="He4" />
</material>
<material id="17" name="Zirc4">
<density units="sum" />
<nuclide ao="0.02217" name="Zr90" />
<nuclide ao="0.004781" name="Zr91" />
<nuclide ao="0.007228" name="Zr92" />
<nuclide ao="0.007169" name="Zr94" />
<nuclide ao="0.001131" name="Zr96" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>-0.9899494936611666 -0.9899494936611666 0.0 0.9899494936611666 0.9899494936611666 10.0</parameters>
</space>
</source>
<output>
<summary>false</summary>
</output>
<seed>22</seed>
</settings>

2
tests/regression_tests/lattice_hex_coincident/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
1.741370E+00 1.384609E-03

150
tests/regression_tests/lattice_hex_coincident/test.py Normal file
View file

@ -0,0 +1,150 @@
import numpy as np
import openmc
from tests.testing_harness import PyAPITestHarness
class HexLatticeCoincidentTestHarness(PyAPITestHarness):
def _build_inputs(self):
materials = openmc.Materials()
fuel_mat = openmc.Material()
fuel_mat.add_nuclide('U235', 4.9817E-03, 'ao')
materials.append(fuel_mat)
matrix = openmc.Material()
matrix.set_density('atom/b-cm', 8.7742E-02)
matrix.add_element('C', 1.0, 'ao')
matrix.add_s_alpha_beta('c_Graphite')
materials.append(matrix)
lead = openmc.Material(name="Lead")
lead.set_density('g/cm3', 10.32)
lead.add_nuclide('Pb204', 0.014, 'ao')
lead.add_nuclide('Pb206', 0.241, 'ao')
lead.add_nuclide('Pb207', 0.221, 'ao')
lead.add_nuclide('Pb208', 0.524, 'ao')
materials.append(lead)
coolant = openmc.Material()
coolant.set_density('atom/b-cm', 5.4464E-04)
coolant.add_nuclide('He4', 1.0, 'ao')
materials.append(coolant)
zirc = openmc.Material(name="Zirc4")
zirc.add_nuclide('Zr90', 2.217E-02, 'ao')
zirc.add_nuclide('Zr91', 4.781E-03, 'ao')
zirc.add_nuclide('Zr92', 7.228E-03, 'ao')
zirc.add_nuclide('Zr94', 7.169E-03, 'ao')
zirc.add_nuclide('Zr96', 1.131E-03, 'ao')
materials.append(zirc)
materials.export_to_xml()
### Geometry ###
pin_rad = 0.7 # cm
assembly_pitch = 1.4 # cm
cool_rad = 0.293 # cm
zirc_clad_thickness = 0.057 # cm
zirc_ir = cool_rad # cm
zirc_or = cool_rad + zirc_clad_thickness # cm
lead_thickness = 0.002 # cm
lead_ir = zirc_or # cm
lead_or = zirc_or + lead_thickness # cm
cyl = openmc.ZCylinder(x0=0., y0=0., r=pin_rad)
fuel_btm = openmc.ZPlane(z0=0.0, boundary_type = 'reflective')
fuel_top = openmc.ZPlane(z0=10.0, boundary_type = 'reflective')
region = -cyl & +fuel_btm & -fuel_top
container = openmc.Cell(region=region)
container.fill = fuel_mat
fuel_outside = openmc.Cell()
fuel_outside.region = +cyl
fuel_outside.fill = matrix
fuel_ch_univ = openmc.Universe(cells=[container, fuel_outside])
# Coolant Channel
cool_outer = openmc.ZCylinder(x0=0.0, y0=0.0, r=cool_rad)
zirc_outer = openmc.ZCylinder(x0=0.0, y0=0.0, r=zirc_or)
lead_outer = openmc.ZCylinder(x0=0.0, y0=0.0, r=lead_or)
coolant_ch = openmc.Cell(name="coolant")
coolant_ch.region = -cool_outer & +fuel_btm & -fuel_top
coolant_ch.fill = coolant
zirc_shell = openmc.Cell(name="zirconium_shell")
zirc_shell.region = +cool_outer & -zirc_outer & +fuel_btm & -fuel_top
zirc_shell.fill = zirc
lead_shell = openmc.Cell(name="lead_shell")
lead_shell.region = +zirc_outer & -lead_outer & +fuel_btm & -fuel_top
lead_shell.fill = lead
coolant_matrix = openmc.Cell(name="matrix coolant surround")
coolant_matrix.region = +lead_outer & +fuel_btm & -fuel_top
coolant_matrix.fill = matrix
coolant_channel = [coolant_ch, zirc_shell, lead_shell, coolant_matrix]
coolant_univ = openmc.Universe(name="coolant universe")
coolant_univ.add_cells(coolant_channel)
half_width = assembly_pitch # cm
edge_length = (2./np.sqrt(3.0)) * half_width
inf_mat = openmc.Cell()
inf_mat.fill = matrix
inf_mat_univ = openmc.Universe(cells=[inf_mat,])
# a hex surface for the core to go inside of
hexprism = openmc.model.get_hexagonal_prism(edge_length=edge_length,
origin=(0.0, 0.0),
boundary_type = 'reflective',
orientation='x')
pincell_only_lattice = openmc.HexLattice(name="regular fuel assembly")
pincell_only_lattice.center = (0., 0.)
pincell_only_lattice.pitch = (assembly_pitch,)
pincell_only_lattice.outer = inf_mat_univ
# setup hex rings
ring0 = [fuel_ch_univ]
ring1 = [coolant_univ] * 6
pincell_only_lattice.universes = [ring1, ring0]
pincell_only_cell = openmc.Cell(name="container cell")
pincell_only_cell.region = hexprism & +fuel_btm & -fuel_top
pincell_only_cell.fill = pincell_only_lattice
root_univ = openmc.Universe(name="root universe", cells=[pincell_only_cell,])
geom = openmc.Geometry(root_univ)
geom.export_to_xml()
### Settings ###
settings = openmc.Settings()
settings.run_mode = 'eigenvalue'
source = openmc.Source()
corner_dist = np.sqrt(2) * pin_rad
ll = [-corner_dist, -corner_dist, 0.0]
ur = [corner_dist, corner_dist, 10.0]
source.space = openmc.stats.Box(ll, ur)
source.strength = 1.0
settings.source = source
settings.output = {'summary' : False}
settings.batches = 10
settings.inactive = 5
settings.particles = 1000
settings.seed = 22
settings.export_to_xml()
def test_lattice_hex_coincident_surf():
harness = HexLatticeCoincidentTestHarness('statepoint.10.h5')
harness.main()