mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Simplify rectangular lattice crossing and correctly handle corner checks (#3703)
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
GuySten
GitHub
parent
60ddafa9b3
commit
818fd11b18
5 changed files with 358 additions and 38 deletions
|
|
@ -260,46 +260,24 @@ std::pair<double, array<int, 3>> RectLattice::distance(
|
|||
// Determine the oncoming edge.
|
||||
double x0 {copysign(0.5 * pitch_[0], u.x)};
|
||||
double y0 {copysign(0.5 * pitch_[1], u.y)};
|
||||
double z0;
|
||||
|
||||
// Left and right sides
|
||||
double d {INFTY};
|
||||
array<int, 3> lattice_trans;
|
||||
if ((std::abs(x - x0) > FP_PRECISION) && u.x != 0) {
|
||||
d = (x0 - x) / u.x;
|
||||
if (u.x > 0) {
|
||||
lattice_trans = {1, 0, 0};
|
||||
} else {
|
||||
lattice_trans = {-1, 0, 0};
|
||||
}
|
||||
}
|
||||
|
||||
// Front and back sides
|
||||
if ((std::abs(y - y0) > FP_PRECISION) && u.y != 0) {
|
||||
double this_d = (y0 - y) / u.y;
|
||||
if (this_d < d) {
|
||||
d = this_d;
|
||||
if (u.y > 0) {
|
||||
lattice_trans = {0, 1, 0};
|
||||
} else {
|
||||
lattice_trans = {0, -1, 0};
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Top and bottom sides
|
||||
double d = std::min(
|
||||
u.x != 0.0 ? (x0 - x) / u.x : INFTY, u.y != 0.0 ? (y0 - y) / u.y : INFTY);
|
||||
if (is_3d_) {
|
||||
double z0 {copysign(0.5 * pitch_[2], u.z)};
|
||||
if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) {
|
||||
double this_d = (z0 - z) / u.z;
|
||||
if (this_d < d) {
|
||||
d = this_d;
|
||||
if (u.z > 0) {
|
||||
lattice_trans = {0, 0, 1};
|
||||
} else {
|
||||
lattice_trans = {0, 0, -1};
|
||||
}
|
||||
}
|
||||
}
|
||||
z0 = copysign(0.5 * pitch_[2], u.z);
|
||||
d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY);
|
||||
}
|
||||
|
||||
// Determine which lattice boundaries are being crossed
|
||||
array<int, 3> lattice_trans = {0, 0, 0};
|
||||
if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION)
|
||||
lattice_trans[0] = copysign(1, u.x);
|
||||
if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION)
|
||||
lattice_trans[1] = copysign(1, u.y);
|
||||
if (is_3d_) {
|
||||
if (u.z != 0.0 && std::abs(z + u.z * d - z0) < FP_PRECISION)
|
||||
lattice_trans[2] = copysign(1, u.z);
|
||||
}
|
||||
|
||||
return {d, lattice_trans};
|
||||
|
|
|
|||
|
|
@ -0,0 +1,58 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<model>
|
||||
<materials>
|
||||
<material id="1">
|
||||
<density value="0.001" units="g/cm3"/>
|
||||
<nuclide name="N14" ao="1.0"/>
|
||||
</material>
|
||||
<material id="2">
|
||||
<density value="7.0" units="g/cm3"/>
|
||||
<nuclide name="Fe56" ao="1.0"/>
|
||||
</material>
|
||||
</materials>
|
||||
<geometry>
|
||||
<cell id="1" material="2" universe="1"/>
|
||||
<cell id="2" material="1" universe="2"/>
|
||||
<cell id="3" material="1" region="-5 (-1 | 2 | -3 | 4)" universe="4"/>
|
||||
<cell id="4" fill="3" region="1 -2 3 -4" universe="4"/>
|
||||
<lattice id="3">
|
||||
<pitch>10.0 10.0</pitch>
|
||||
<dimension>2 2</dimension>
|
||||
<lower_left>-10.0 -10.0</lower_left>
|
||||
<universes>
|
||||
1 2
|
||||
2 1 </universes>
|
||||
</lattice>
|
||||
<surface id="1" name="minimum x" type="x-plane" coeffs="-10.0"/>
|
||||
<surface id="2" name="maximum x" type="x-plane" coeffs="10.0"/>
|
||||
<surface id="3" name="minimum y" type="y-plane" coeffs="-10.0"/>
|
||||
<surface id="4" name="maximum y" type="y-plane" coeffs="10.0"/>
|
||||
<surface id="5" type="z-cylinder" boundary="vacuum" coeffs="0.0 0.0 20.0"/>
|
||||
</geometry>
|
||||
<settings>
|
||||
<run_mode>fixed source</run_mode>
|
||||
<particles>1000</particles>
|
||||
<batches>10</batches>
|
||||
<source type="independent" strength="1.0" particle="neutron">
|
||||
<space type="point">
|
||||
<parameters>-0.7071067811865476 -0.7071067811865475 0.0</parameters>
|
||||
</space>
|
||||
<angle type="monodirectional" reference_uvw="0.7071067811865476 0.7071067811865475 0.0"/>
|
||||
</source>
|
||||
</settings>
|
||||
<tallies>
|
||||
<mesh id="1">
|
||||
<dimension>10 10</dimension>
|
||||
<lower_left>-20.0 -20.0</lower_left>
|
||||
<upper_right>20.0 20.0</upper_right>
|
||||
</mesh>
|
||||
<filter id="1" type="mesh">
|
||||
<bins>1</bins>
|
||||
</filter>
|
||||
<tally id="1">
|
||||
<filters>1</filters>
|
||||
<scores>flux</scores>
|
||||
<estimator>tracklength</estimator>
|
||||
</tally>
|
||||
</tallies>
|
||||
</model>
|
||||
201
tests/regression_tests/lattice_corner_crossing/results_true.dat
Normal file
201
tests/regression_tests/lattice_corner_crossing/results_true.dat
Normal file
|
|
@ -0,0 +1,201 @@
|
|||
tally 1:
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
2.648775E-03
|
||||
7.016009E-06
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
2.981847E-03
|
||||
8.891414E-06
|
||||
3.331677E-03
|
||||
1.110007E-05
|
||||
6.742237E-03
|
||||
4.545776E-05
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.388206E-03
|
||||
1.925636E-05
|
||||
2.032104E-03
|
||||
4.129446E-06
|
||||
3.014946E-03
|
||||
9.089898E-06
|
||||
7.058968E-03
|
||||
4.982903E-05
|
||||
9.013188E-04
|
||||
8.123755E-07
|
||||
7.501461E-04
|
||||
5.627191E-07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
5.134898E-03
|
||||
2.636718E-05
|
||||
1.207302E-03
|
||||
1.457579E-06
|
||||
3.045925E-03
|
||||
9.277657E-06
|
||||
4.132399E-03
|
||||
1.707672E-05
|
||||
1.055271E-02
|
||||
5.829711E-05
|
||||
4.014909E-03
|
||||
1.611949E-05
|
||||
2.698215E-03
|
||||
7.280363E-06
|
||||
1.751215E-02
|
||||
1.073504E-04
|
||||
1.356908E-02
|
||||
9.589157E-05
|
||||
5.451466E-03
|
||||
2.971848E-05
|
||||
3.165676E-04
|
||||
1.002151E-07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
3.724938E-03
|
||||
1.387516E-05
|
||||
2.812863E-03
|
||||
7.278837E-06
|
||||
1.000700E+01
|
||||
1.001402E+01
|
||||
2.345904E-02
|
||||
2.127641E-04
|
||||
2.451653E-02
|
||||
1.750725E-04
|
||||
8.695904E-03
|
||||
5.553981E-05
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.947565E-03
|
||||
2.447840E-05
|
||||
1.725721E-02
|
||||
8.878803E-05
|
||||
5.656444E+01
|
||||
3.199538E+02
|
||||
2.159961E-02
|
||||
1.034364E-04
|
||||
3.091063E-02
|
||||
5.543317E-04
|
||||
4.232209E-03
|
||||
1.791159E-05
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.716170E-03
|
||||
2.224226E-05
|
||||
1.321855E-02
|
||||
5.288265E-05
|
||||
5.309576E-02
|
||||
7.984865E-04
|
||||
5.656179E+01
|
||||
3.199238E+02
|
||||
2.399311E-02
|
||||
2.936670E-04
|
||||
4.680753E-02
|
||||
9.526509E-04
|
||||
1.074699E-02
|
||||
1.154978E-04
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.526425E-03
|
||||
2.048853E-05
|
||||
3.983003E-03
|
||||
9.553596E-06
|
||||
1.726789E-02
|
||||
9.579848E-05
|
||||
3.355605E-02
|
||||
3.331438E-04
|
||||
3.578422E-02
|
||||
3.560242E-04
|
||||
5.649209E+01
|
||||
3.191358E+02
|
||||
2.591495E-02
|
||||
4.169753E-04
|
||||
9.484744E-03
|
||||
8.065715E-05
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.232587E-03
|
||||
9.113263E-06
|
||||
5.787378E-03
|
||||
2.194940E-05
|
||||
6.193626E-03
|
||||
2.357423E-05
|
||||
4.552787E-03
|
||||
6.987360E-06
|
||||
9.641308E-03
|
||||
3.506339E-05
|
||||
1.303892E-02
|
||||
4.378408E-05
|
||||
1.831744E-02
|
||||
8.836683E-05
|
||||
3.024725E+01
|
||||
9.148969E+01
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
1.488204E-03
|
||||
2.214752E-06
|
||||
5.067487E-03
|
||||
1.520987E-05
|
||||
7.550927E-03
|
||||
2.993991E-05
|
||||
3.215559E-03
|
||||
1.033982E-05
|
||||
4.258525E-03
|
||||
1.205064E-05
|
||||
2.115487E-03
|
||||
4.475286E-06
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
83
tests/regression_tests/lattice_corner_crossing/test.py
Normal file
83
tests/regression_tests/lattice_corner_crossing/test.py
Normal file
|
|
@ -0,0 +1,83 @@
|
|||
"""
|
||||
This test is designed to ensure that we account for potential corner crossings
|
||||
in floating point precision.
|
||||
|
||||
"""
|
||||
|
||||
from math import pi, cos, sin
|
||||
|
||||
import openmc
|
||||
import pytest
|
||||
|
||||
from tests.testing_harness import PyAPITestHarness
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def model():
|
||||
model = openmc.Model()
|
||||
|
||||
# Length of lattice on each side
|
||||
lat_size = 20.0
|
||||
|
||||
# Angle that we're crossing the corner at
|
||||
phi = pi/4.0
|
||||
|
||||
air = openmc.Material()
|
||||
air.set_density('g/cm3', 0.001)
|
||||
air.add_nuclide('N14', 1.0)
|
||||
metal = openmc.Material()
|
||||
metal.set_density('g/cm3', 7.0)
|
||||
metal.add_nuclide('Fe56', 1.0)
|
||||
|
||||
metal_cell = openmc.Cell(fill=metal)
|
||||
metal_uni = openmc.Universe(cells=[metal_cell])
|
||||
|
||||
air_cell = openmc.Cell(fill=air)
|
||||
air_uni = openmc.Universe(cells=[air_cell])
|
||||
|
||||
# Define a checkerboard lattice
|
||||
lattice = openmc.RectLattice()
|
||||
lattice.lower_left = (-lat_size/2.0, -lat_size/2.0)
|
||||
lattice.pitch = (lat_size/2, lat_size/2)
|
||||
lattice.universes = [
|
||||
[metal_uni, air_uni],
|
||||
[air_uni, metal_uni]
|
||||
]
|
||||
|
||||
box = openmc.model.RectangularPrism(lat_size, lat_size)
|
||||
cyl = openmc.ZCylinder(r=lat_size, boundary_type='vacuum')
|
||||
outside_lattice = openmc.Cell(region=-cyl & +box, fill=air)
|
||||
inside_lattice = openmc.Cell(region=-box, fill=lattice)
|
||||
|
||||
model.geometry = openmc.Geometry([outside_lattice, inside_lattice])
|
||||
|
||||
# Set all runtime parameters
|
||||
model.settings.run_mode = 'fixed source'
|
||||
model.settings.batches = 10
|
||||
model.settings.particles = 1000
|
||||
|
||||
# Define a source located outside the lattice and pointing straight into its
|
||||
# corner at 45 degrees
|
||||
model.settings.source = openmc.IndependentSource(
|
||||
space=openmc.stats.Point((-cos(phi), -sin(phi), 0.0)),
|
||||
angle=openmc.stats.Monodirectional((cos(phi), sin(phi), 0.0))
|
||||
)
|
||||
|
||||
# Create a mesh tally
|
||||
mesh = openmc.RegularMesh()
|
||||
mesh.dimension = (10, 10)
|
||||
mesh.lower_left = (-lat_size, -lat_size)
|
||||
mesh.upper_right = (lat_size, lat_size)
|
||||
mesh_filter = openmc.MeshFilter(mesh)
|
||||
tally = openmc.Tally(tally_id=1)
|
||||
tally.filters = [mesh_filter]
|
||||
tally.scores = ['flux']
|
||||
tally.estimator = 'tracklength'
|
||||
model.tallies = [tally]
|
||||
|
||||
return model
|
||||
|
||||
|
||||
def test_lattice_corner_crossing(model):
|
||||
harness = PyAPITestHarness('statepoint.10.h5', model)
|
||||
harness.main()
|
||||
Loading…
Add table
Add a link
Reference in a new issue