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

Merge remote-tracking branch 'upstream/develop' into capi_sourcepoint

Browse source
This commit is contained in:
Sterling Harper 2018-08-14 11:06:50 -04:00
commit ec67e25111
29 changed files with 981 additions and 551 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -395,6 +395,7 @@ add_library(libopenmc SHARED
src/mgxs_interface.cpp
src/particle.cpp
src/plot.cpp
src/position.cpp
src/pugixml/pugixml_c.cpp
src/random_lcg.cpp
src/scattdata.cpp

1
docs/source/pythonapi/base.rst
View file

@ -122,6 +122,7 @@ Constructing Tallies
openmc.SpatialLegendreFilter
openmc.SphericalHarmonicsFilter
openmc.ZernikeFilter
openmc.ZernikeRadialFilter
openmc.ParticleFilter
openmc.Mesh
openmc.Trigger

9
openmc/capi/filter.py
View file

@ -20,7 +20,7 @@ __all__ = ['Filter', 'AzimuthalFilter', 'CellFilter',
'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MeshFilter',
'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SphericalHarmonicsFilter',
'SpatialLegendreFilter', 'SurfaceFilter',
'UniverseFilter', 'ZernikeFilter', 'filters']
'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters']
# Tally functions
_dll.openmc_cell_filter_get_bins.argtypes = [
@ -360,6 +360,10 @@ class ZernikeFilter(Filter):
_dll.openmc_zernike_filter_set_order(self._index, order)
class ZernikeRadialFilter(ZernikeFilter):
filter_type = 'zernikeradial'
_FILTER_TYPE_MAP = {
'azimuthal': AzimuthalFilter,
'cell': CellFilter,
@ -380,7 +384,8 @@ _FILTER_TYPE_MAP = {
'spatiallegendre': SpatialLegendreFilter,
'surface': SurfaceFilter,
'universe': UniverseFilter,
'zernike': ZernikeFilter
'zernike': ZernikeFilter,
'zernikeradial': ZernikeRadialFilter
}

29
openmc/capi/math.py
View file

@ -21,6 +21,9 @@ _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
_dll.calc_zn_c.restype = None
_dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
_dll.calc_zn_rad_c.restype = None
_dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
_dll.rotate_angle_c.restype = None
_dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
POINTER(c_double)]
@ -155,6 +158,32 @@ def calc_zn(n, rho, phi):
return zn
def calc_zn_rad(n, rho):
""" Calculate the even orders in n-th order modified Zernike polynomial
moment with no azimuthal dependency (m=0) for a given radial location in
the unit disk. The normalization of the polynomials is such that the
integral of Z_pq*Z_pq over the unit disk is exactly pi.
Parameters
----------
n : int
Maximum order
rho : float
Radial location in the unit disk
Returns
-------
numpy.ndarray
Corresponding resulting list of coefficients
"""
num_bins = n // 2 + 1
zn_rad = np.zeros(num_bins, dtype=np.float64)
_dll.calc_zn_rad_c(n, rho, zn_rad)
return zn_rad
def rotate_angle(uvw0, mu, phi=None):
""" Rotates direction cosines through a polar angle whose cosine is
mu and through an azimuthal angle sampled uniformly.

2
openmc/filter.py
View file

@ -22,7 +22,7 @@ _FILTER_TYPES = (
'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre',
'sphericalharmonics', 'zernike', 'particle'
'sphericalharmonics', 'zernike', 'zernikeradial', 'particle'
)
_CURRENT_NAMES = (

72
openmc/filter_expansion.py
View file

@ -234,8 +234,8 @@ class SphericalHarmonicsFilter(ExpansionFilter):
r"""Score spherical harmonic expansion moments up to specified order.
This filter allows you to obtain real spherical harmonic moments of either
the particle's direction or the cosine of the scattering angle. Specifying a
filter with order :math:`\ell` tallies moments for all orders from 0 to
the particle's direction or the cosine of the scattering angle. Specifying
a filter with order :math:`\ell` tallies moments for all orders from 0 to
:math:`\ell`.
Parameters
@ -342,11 +342,11 @@ class ZernikeFilter(ExpansionFilter):
\frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}.
With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk
is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is
2 if :math:`m` equals 0 and 1 otherwise.
is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where
:math:`\epsilon_m` is 2 if :math:`m` equals 0 and 1 otherwise.
Specifying a filter with order N tallies moments for all :math:`n` from 0 to
N and each value of :math:`m`. The ordering of the Zernike polynomial
Specifying a filter with order N tallies moments for all :math:`n` from 0
to N and each value of :math:`m`. The ordering of the Zernike polynomial
moments follows the ANSI Z80.28 standard, where the one-dimensional index
:math:`j` corresponds to the :math:`n` and :math:`m` by
@ -463,3 +463,63 @@ class ZernikeFilter(ExpansionFilter):
subelement.text = str(self.r)
return element
class ZernikeRadialFilter(ZernikeFilter):
r"""Score the :math:`m = 0` (radial variation only) Zernike moments up to
specified order.
The Zernike polynomials are defined the same as in :class:`ZernikeFilter`.
.. math::
Z_n^{0}(\rho, \theta) = R_n^{0}(\rho)
where the radial polynomials are
.. math::
R_n^{0}(\rho) = \sum\limits_{k=0}^{n/2} \frac{(-1)^k (n-k)!}{k! ((
\frac{n}{2} - k)!)^{2}} \rho^{n-2k}.
With this definition, the integral of :math:`(Z_n^0)^2` over the unit disk
is :math:`\frac{\pi}{n+1}`.
If there is only radial dependency, the polynomials are integrated over
the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta)
= R_n^{0}(\rho)`. Note that :math:`n` could only be even orders.
Therefore, for a radial Zernike polynomials up to order of :math:`n`,
there are :math:`\frac{n}{2} + 1` terms in total. The indexing is from the
lowest even order (0) to highest even order.
Parameters
----------
order : int
Maximum radial Zernike polynomial order
x : float
x-coordinate of center of circle for normalization
y : float
y-coordinate of center of circle for normalization
r : int or None
Radius of circle for normalization
Attributes
----------
order : int
Maximum radial Zernike polynomial order
x : float
x-coordinate of center of circle for normalization
y : float
y-coordinate of center of circle for normalization
r : int or None
Radius of circle for normalization
id : int
Unique identifier for the filter
num_bins : int
The number of filter bins
"""
@ExpansionFilter.order.setter
def order(self, order):
ExpansionFilter.order.__set__(self, order)
self.bins = ['Z{},0'.format(n) for n in range(0, order+1, 2)]

34
src/cell.cpp
View file

@ -279,33 +279,31 @@ Cell::Cell(pugi::xml_node cell_node)
//==============================================================================
bool
Cell::contains(const double xyz[3], const double uvw[3],
int32_t on_surface) const
Cell::contains(Position r, Direction u, int32_t on_surface) const
{
if (simple) {
return contains_simple(xyz, uvw, on_surface);
return contains_simple(r, u, on_surface);
} else {
return contains_complex(xyz, uvw, on_surface);
return contains_complex(r, u, on_surface);
}
}
//==============================================================================
std::pair<double, int32_t>
Cell::distance(const double xyz[3], const double uvw[3],
int32_t on_surface) const
Cell::distance(Position r, Direction u, int32_t on_surface) const
{
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
for (int32_t token : rpn) {
// Ignore this token if it corresponds to an operator rather than a region.
if (token >= OP_UNION) {continue;}
if (token >= OP_UNION) continue;
// Calculate the distance to this surface.
// Note the off-by-one indexing
bool coincident {token == on_surface};
double d {surfaces_c[abs(token)-1]->distance(xyz, uvw, coincident)};
double d {surfaces_c[abs(token)-1]->distance(r, u, coincident)};
// Check if this distance is the new minimum.
if (d < min_dist) {
@ -356,8 +354,7 @@ Cell::to_hdf5(hid_t cell_group) const
//==============================================================================
bool
Cell::contains_simple(const double xyz[3], const double uvw[3],
int32_t on_surface) const
Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
{
for (int32_t token : rpn) {
if (token < OP_UNION) {
@ -370,7 +367,7 @@ Cell::contains_simple(const double xyz[3], const double uvw[3],
return false;
} else {
// Note the off-by-one indexing
bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);
bool sense = surfaces_c[abs(token)-1]->sense(r, u);
if (sense != (token > 0)) {return false;}
}
}
@ -381,8 +378,7 @@ Cell::contains_simple(const double xyz[3], const double uvw[3],
//==============================================================================
bool
Cell::contains_complex(const double xyz[3], const double uvw[3],
int32_t on_surface) const
Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
{
// Make a stack of booleans. We don't know how big it needs to be, but we do
// know that rpn.size() is an upper-bound.
@ -413,7 +409,7 @@ Cell::contains_complex(const double xyz[3], const double uvw[3],
stack[i_stack] = false;
} else {
// Note the off-by-one indexing
bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);;
bool sense = surfaces_c[abs(token)-1]->sense(r, u);;
stack[i_stack] = (sense == (token > 0));
}
}
@ -494,12 +490,18 @@ extern "C" {
bool cell_simple(Cell *c) {return c->simple;}
bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface)
{return c->contains(xyz, uvw, on_surface);}
{
Position r {xyz};
Direction u {uvw};
return c->contains(r, u, on_surface);
}
void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface,
double *min_dist, int32_t *i_surf)
{
std::pair<double, int32_t> out = c->distance(xyz, uvw, on_surface);
Position r {xyz};
Direction u {uvw};
std::pair<double, int32_t> out = c->distance(r, u, on_surface);
*min_dist = out.first;
*i_surf = out.second;
}

26
src/cell.h
View file

@ -1,5 +1,5 @@
#ifndef CELL_H
#define CELL_H
#ifndef OPENMC_CELL_H
#define OPENMC_CELL_H
#include <cstdint>
#include <string>
@ -9,6 +9,8 @@
#include "hdf5.h"
#include "pugixml.hpp"
#include "position.h"
namespace openmc {
@ -90,29 +92,27 @@ public:
//! provides a performance benefit for the common case. In
//! contains_complex, we evaluate the RPN expression using a stack, similar to
//! how a RPN calculator would work.
//! @param xyz[3] The 3D Cartesian coordinate to check.
//! @param uvw[3] A direction used to "break ties" the coordinates are very
//! \param r The 3D Cartesian coordinate to check.
//! \param u A direction used to "break ties" the coordinates are very
//! close to a surface.
//! @param on_surface The signed index of a surface that the coordinate is
//! \param on_surface The signed index of a surface that the coordinate is
//! known to be on. This index takes precedence over surface sense
//! calculations.
bool
contains(const double xyz[3], const double uvw[3], int32_t on_surface) const;
contains(Position r, Direction u, int32_t on_surface) const;
//! Find the oncoming boundary of this cell.
std::pair<double, int32_t>
distance(const double xyz[3], const double uvw[3], int32_t on_surface) const;
distance(Position r, Direction u, int32_t on_surface) const;
//! \brief Write cell information to an HDF5 group.
//! @param group_id An HDF5 group id.
//! \param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
protected:
bool contains_simple(const double xyz[3], const double uvw[3],
int32_t on_surface) const;
bool contains_complex(const double xyz[3], const double uvw[3],
int32_t on_surface) const;
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
};
} // namespace openmc
#endif // CELL_H
#endif // OPENMC_CELL_H

6
src/constants.F90
View file

@ -357,7 +357,7 @@ module constants
integer, parameter :: NO_BIN_FOUND = -1
! Tally filter and map types
integer, parameter :: N_FILTER_TYPES = 21
integer, parameter :: N_FILTER_TYPES = 22
integer, parameter :: &
FILTER_UNIVERSE = 1, &
FILTER_MATERIAL = 2, &
@ -379,7 +379,9 @@ module constants
FILTER_SPH_HARMONICS = 18, &
FILTER_SPTL_LEGENDRE = 19, &
FILTER_ZERNIKE = 20, &
FILTER_PARTICLE = 21
FILTER_ZERNIKE_RADIAL = 21, &
FILTER_PARTICLE = 22
! Mesh types
integer, parameter :: &

6
src/eigenvalue.F90
View file

@ -380,7 +380,11 @@ contains
! Determine overall generation and number of active generations
i = overall_generation()
n = i - n_inactive*gen_per_batch
if (current_batch > n_inactive) then
n = gen_per_batch*n_realizations + current_gen
else
n = 0
end if
if (n <= 0) then
! For inactive generations, use current generation k as estimate for next

10
src/geometry.h
View file

@ -1,6 +1,10 @@
#ifndef GEOMETRY_H
#define GEOMETRY_H
#ifndef OPENMC_GEOMETRY_H
#define OPENMC_GEOMETRY_H
namespace openmc {
extern "C" int openmc_root_universe;
#endif // GEOMETRY_H
} // namespace openmc
#endif // OPENMC_GEOMETRY_H

15
src/hdf5_interface.h
View file

@ -1,5 +1,5 @@
#ifndef HDF5_INTERFACE_H
#define HDF5_INTERFACE_H
#ifndef OPENMC_HDF5_INTERFACE_H
#define OPENMC_HDF5_INTERFACE_H
#include "hdf5.h"
#include "hdf5_hl.h"
@ -11,6 +11,8 @@
#include <vector>
#include <complex.h>
#include "position.h"
namespace openmc {
@ -181,5 +183,12 @@ write_dataset(hid_t obj_id, const char* name, const std::array<T, N>& buffer)
write_dataset(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data(), false);
}
inline void
write_dataset(hid_t obj_id, const char* name, Position r)
{
std::array<double, 3> buffer {r.x, r.y, r.z};
write_dataset(obj_id, name, buffer);
}
} // namespace openmc
#endif //HDF5_INTERFACE_H
#endif // OPENMC_HDF5_INTERFACE_H

158
src/lattice.cpp
View file

@ -223,26 +223,23 @@ RectLattice::are_valid_indices(const int i_xyz[3]) const
//==============================================================================
std::pair<double, std::array<int, 3>>
RectLattice::distance(const double xyz[3], const double uvw[3],
const int i_xyz[3]) const
RectLattice::distance(Position r, Direction u, const int i_xyz[3]) const
{
// Get short aliases to the coordinates.
double x {xyz[0]};
double y {xyz[1]};
double z {xyz[2]};
double u {uvw[0]};
double v {uvw[1]};
double x = r.x;
double y = r.y;
double z = r.z;
// Determine the oncoming edge.
double x0 {copysign(0.5 * pitch[0], u)};
double y0 {copysign(0.5 * pitch[1], v)};
double x0 {copysign(0.5 * pitch[0], u.x)};
double y0 {copysign(0.5 * pitch[1], u.y)};
// Left and right sides
double d {INFTY};
std::array<int, 3> lattice_trans;
if ((std::abs(x - x0) > FP_PRECISION) && u != 0) {
d = (x0 - x) / u;
if (u > 0) {
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};
@ -250,11 +247,11 @@ RectLattice::distance(const double xyz[3], const double uvw[3],
}
// Front and back sides
if ((std::abs(y - y0) > FP_PRECISION) && v != 0) {
double this_d = (y0 - y) / v;
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 (v > 0) {
if (u.y > 0) {
lattice_trans = {0, 1, 0};
} else {
lattice_trans = {0, -1, 0};
@ -264,13 +261,12 @@ RectLattice::distance(const double xyz[3], const double uvw[3],
// Top and bottom sides
if (is_3d) {
double w {uvw[2]};
double z0 {copysign(0.5 * pitch[2], w)};
if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
double this_d = (z0 - z) / w;
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 (w > 0) {
if (u.z > 0) {
lattice_trans = {0, 0, 1};
} else {
lattice_trans = {0, 0, -1};
@ -285,13 +281,13 @@ RectLattice::distance(const double xyz[3], const double uvw[3],
//==============================================================================
std::array<int, 3>
RectLattice::get_indices(const double xyz[3]) const
RectLattice::get_indices(Position r) const
{
int ix {static_cast<int>(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))-1};
int iy {static_cast<int>(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))-1};
int ix {static_cast<int>(std::ceil((r.x - lower_left.x) / pitch.x))-1};
int iy {static_cast<int>(std::ceil((r.y - lower_left.y) / pitch.y))-1};
int iz;
if (is_3d) {
iz = static_cast<int>(std::ceil((xyz[2] - lower_left[2]) / pitch[2]))-1;
iz = static_cast<int>(std::ceil((r.z - lower_left.z) / pitch.z))-1;
} else {
iz = 0;
}
@ -300,18 +296,15 @@ RectLattice::get_indices(const double xyz[3]) const
//==============================================================================
std::array<double, 3>
RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
Position
RectLattice::get_local_position(Position r, const int i_xyz[3]) const
{
std::array<double, 3> local_xyz;
local_xyz[0] = global_xyz[0] - (lower_left[0] + (i_xyz[0] + 0.5)*pitch[0]);
local_xyz[1] = global_xyz[1] - (lower_left[1] + (i_xyz[1] + 0.5)*pitch[1]);
r.x -= (lower_left.x + (i_xyz[0] + 0.5)*pitch.x);
r.y -= (lower_left.y + (i_xyz[1] + 0.5)*pitch.y);
if (is_3d) {
local_xyz[2] = global_xyz[2] - (lower_left[2] + (i_xyz[2] + 0.5)*pitch[2]);
} else {
local_xyz[2] = global_xyz[2];
r.z -= (lower_left.z + (i_xyz[2] + 0.5)*pitch.z);
}
return local_xyz;
return r;
}
//==============================================================================
@ -583,12 +576,11 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const
//==============================================================================
std::pair<double, std::array<int, 3>>
HexLattice::distance(const double xyz[3], const double uvw[3],
const int i_xyz[3]) const
HexLattice::distance(Position r, Direction u, const int i_xyz[3]) const
{
// Compute the direction on the hexagonal basis.
double beta_dir = uvw[0] * std::sqrt(3.0) / 2.0 + uvw[1] / 2.0;
double gamma_dir = uvw[0] * std::sqrt(3.0) / 2.0 - uvw[1] / 2.0;
double beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0;
double gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0;
// Note that hexagonal lattice distance calculations are performed
// using the particle's coordinates relative to the neighbor lattice
@ -600,15 +592,15 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
double d {INFTY};
std::array<int, 3> lattice_trans;
double edge = -copysign(0.5*pitch[0], beta_dir); // Oncoming edge
std::array<double, 3> xyz_t;
Position r_t;
if (beta_dir > 0) {
const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]};
xyz_t = get_local_xyz(xyz, i_xyz_t);
r_t = get_local_position(r, i_xyz_t);
} else {
const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]};
xyz_t = get_local_xyz(xyz, i_xyz_t);
r_t = get_local_position(r, i_xyz_t);
}
double beta = xyz_t[0] * std::sqrt(3.0) / 2.0 + xyz_t[1] / 2.0;
double beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;
if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) {
d = (edge - beta) / beta_dir;
if (beta_dir > 0) {
@ -622,12 +614,12 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
edge = -copysign(0.5*pitch[0], gamma_dir);
if (gamma_dir > 0) {
const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]};
xyz_t = get_local_xyz(xyz, i_xyz_t);
r_t = get_local_position(r, i_xyz_t);
} else {
const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]};
xyz_t = get_local_xyz(xyz, i_xyz_t);
r_t = get_local_position(r, i_xyz_t);
}
double gamma = xyz_t[0] * std::sqrt(3.0) / 2.0 - xyz_t[1] / 2.0;
double gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;
if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) {
double this_d = (edge - gamma) / gamma_dir;
if (this_d < d) {
@ -641,18 +633,18 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
}
// Upper and lower sides.
edge = -copysign(0.5*pitch[0], uvw[1]);
if (uvw[1] > 0) {
edge = -copysign(0.5*pitch[0], u.y);
if (u.y > 0) {
const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]};
xyz_t = get_local_xyz(xyz, i_xyz_t);
r_t = get_local_position(r, i_xyz_t);
} else {
const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]};
xyz_t = get_local_xyz(xyz, i_xyz_t);
r_t = get_local_position(r, i_xyz_t);
}
if ((std::abs(xyz_t[1] - edge) > FP_PRECISION) && uvw[1] != 0) {
double this_d = (edge - xyz_t[1]) / uvw[1];
if ((std::abs(r_t.y - edge) > FP_PRECISION) && u.y != 0) {
double this_d = (edge - r_t.y) / u.y;
if (this_d < d) {
if (uvw[1] > 0) {
if (u.y > 0) {
lattice_trans = {0, 1, 0};
} else {
lattice_trans = {0, -1, 0};
@ -663,14 +655,13 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
// Top and bottom sides
if (is_3d) {
double z {xyz[2]};
double w {uvw[2]};
double z0 {copysign(0.5 * pitch[1], w)};
if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
double this_d = (z0 - z) / w;
double z = r.z;
double z0 {copysign(0.5 * pitch[1], 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 (w > 0) {
if (u.z > 0) {
lattice_trans = {0, 0, 1};
} else {
lattice_trans = {0, 0, -1};
@ -686,24 +677,24 @@ HexLattice::distance(const double xyz[3], const double uvw[3],
//==============================================================================
std::array<int, 3>
HexLattice::get_indices(const double xyz[3]) const
HexLattice::get_indices(Position r) const
{
// Offset the xyz by the lattice center.
double xyz_o[3] {xyz[0] - center[0], xyz[1] - center[1], xyz[2]};
if (is_3d) {xyz_o[2] -= center[2];}
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;
if (is_3d) {
out[2] = static_cast<int>(std::ceil(xyz_o[2] / pitch[1] + 0.5 * n_axial))-1;
out[2] = static_cast<int>(std::ceil(r_o.z / pitch[1] + 0.5 * n_axial))-1;
} else {
out[2] = 0;
}
// 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 = xyz_o[1] - xyz_o[0] / std::sqrt(3.0);
out[0] = static_cast<int>(std::floor(xyz_o[0]
double alpha = r_o.y - r_o.x / std::sqrt(3.0);
out[0] = static_cast<int>(std::floor(r_o.x
/ (0.5*std::sqrt(3.0) * pitch[0])));
out[1] = static_cast<int>(std::floor(alpha / pitch[0]));
@ -725,8 +716,8 @@ HexLattice::get_indices(const double xyz[3]) const
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 2; j++) {
int i_xyz[3] {out[0] + j, out[1] + i, 0};
std::array<double, 3> xyz_t = get_local_xyz(xyz, i_xyz);
double d = xyz_t[0]*xyz_t[0] + xyz_t[1]*xyz_t[1];
Position r_t = get_local_position(r, i_xyz);
double d = r_t.x*r_t.x + r_t.y*r_t.y;
if (d < d_min) {
d_min = d;
k_min = k;
@ -751,26 +742,19 @@ HexLattice::get_indices(const double xyz[3]) const
//==============================================================================
std::array<double, 3>
HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
Position
HexLattice::get_local_position(Position r, const int i_xyz[3]) const
{
std::array<double, 3> local_xyz;
// x_l = x_g - (center + pitch_x*cos(30)*index_x)
local_xyz[0] = global_xyz[0] - (center[0]
+ std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]);
r.x -= (center.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]);
// y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
local_xyz[1] = global_xyz[1] - (center[1]
+ (i_xyz[1] - n_rings + 1) * pitch[0]
+ (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0);
r.y -= (center.y + (i_xyz[1] - n_rings + 1) * pitch[0]
+ (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0);
if (is_3d) {
local_xyz[2] = global_xyz[2] - center[2]
+ (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1];
} else {
local_xyz[2] = global_xyz[2];
r.z -= center.z - (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1];
}
return local_xyz;
return r;
}
//==============================================================================
@ -908,7 +892,9 @@ extern "C" {
void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3],
const int i_xyz[3], double *d, int lattice_trans[3])
{
std::pair<double, std::array<int, 3>> ld {lat->distance(xyz, uvw, i_xyz)};
Position r {xyz};
Direction u {uvw};
std::pair<double, std::array<int, 3>> ld {lat->distance(r, u, i_xyz)};
*d = ld.first;
lattice_trans[0] = ld.second[0];
lattice_trans[1] = ld.second[1];
@ -917,7 +903,8 @@ extern "C" {
void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3])
{
std::array<int, 3> inds = lat->get_indices(xyz);
Position r {xyz};
std::array<int, 3> inds = lat->get_indices(r);
i_xyz[0] = inds[0];
i_xyz[1] = inds[1];
i_xyz[2] = inds[2];
@ -926,10 +913,11 @@ extern "C" {
void lattice_get_local_xyz(Lattice *lat, const double global_xyz[3],
const int i_xyz[3], double local_xyz[3])
{
std::array<double, 3> xyz = lat->get_local_xyz(global_xyz, i_xyz);
local_xyz[0] = xyz[0];
local_xyz[1] = xyz[1];
local_xyz[2] = xyz[2];
Position global {global_xyz};
Position local = lat->get_local_position(global, i_xyz);
local_xyz[0] = local.x;
local_xyz[1] = local.y;
local_xyz[2] = local.z;
}
int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3])

80
src/lattice.h
View file

@ -1,5 +1,5 @@
#ifndef LATTICE_H
#define LATTICE_H
#ifndef OPENMC_LATTICE_H
#define OPENMC_LATTICE_H
#include <array>
#include <cstdint>
@ -7,10 +7,12 @@
#include <unordered_map>
#include <vector>
#include "constants.h"
#include "hdf5.h"
#include "pugixml.hpp"
#include "constants.h"
#include "position.h"
namespace openmc {
@ -69,54 +71,54 @@ public:
int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map);
//! \brief Check lattice indices.
//! @param i_xyz[3] The indices for a lattice tile.
//! @return true if the given indices fit within the lattice bounds. False
//! \param i_xyz[3] The indices for a lattice tile.
//! \return true if the given indices fit within the lattice bounds. False
//! otherwise.
virtual bool are_valid_indices(const int i_xyz[3]) const = 0;
//! \brief Find the next lattice surface crossing
//! @param xyz[3] A 3D Cartesian coordinate.
//! @param uvw[3] A 3D Cartesian direction.
//! @param i_xyz[3] The indices for a lattice tile.
//! @return The distance to the next crossing and an array indicating how the
//! \param r A 3D Cartesian coordinate.
//! \param u A 3D Cartesian direction.
//! \param i_xyz[3] The indices for a lattice tile.
//! \return The distance to the next crossing and an array indicating how the
//! lattice indices would change after crossing that boundary.
virtual std::pair<double, std::array<int, 3>>
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const
distance(Position r, Direction u, const int i_xyz[3]) const
= 0;
//! \brief Find the lattice tile indices for a given point.
//! @param xyz[3] A 3D Cartesian coordinate.
//! @return An array containing the indices of a lattice tile.
virtual std::array<int, 3> get_indices(const double xyz[3]) const = 0;
//! \param r A 3D Cartesian coordinate.
//! \return An array containing the indices of a lattice tile.
virtual std::array<int, 3> get_indices(Position r) const = 0;
//! \brief Get coordinates local to a lattice tile.
//! @param global_xyz[3] A 3D Cartesian coordinate.
//! @param i_xyz[3] The indices for a lattice tile.
//! @return Local 3D Cartesian coordinates.
virtual std::array<double, 3>
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const = 0;
//! \param r A 3D Cartesian coordinate.
//! \param i_xyz[3] The indices for a lattice tile.
//! \return Local 3D Cartesian coordinates.
virtual Position
get_local_position(Position r, const int i_xyz[3]) const = 0;
//! \brief Check flattened lattice index.
//! @param indx The index for a lattice tile.
//! @return true if the given index fit within the lattice bounds. False
//! \param indx The index for a lattice tile.
//! \return true if the given index fit within the lattice bounds. False
//! otherwise.
virtual bool is_valid_index(int indx) const
{return (indx >= 0) && (indx < universes.size());}
//! \brief Get the distribcell offset for a lattice tile.
//! @param The map index for the target cell.
//! @param i_xyz[3] The indices for a lattice tile.
//! @return Distribcell offset i.e. the largest instance number for the target
//! \param The map index for the target cell.
//! \param i_xyz[3] The indices for a lattice tile.
//! \return Distribcell offset i.e. the largest instance number for the target
//! cell found in the geometry tree under this lattice tile.
virtual int32_t& offset(int map, const int i_xyz[3]) = 0;
//! \brief Convert an array index to a useful human-readable string.
//! @param indx The index for a lattice tile.
//! @return A string representing the lattice tile.
//! \param indx The index for a lattice tile.
//! \return A string representing the lattice tile.
virtual std::string index_to_string(int indx) const = 0;
//! \brief Write lattice information to an HDF5 group.
//! @param group_id An HDF5 group id.
//! \param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
protected:
@ -193,12 +195,12 @@ public:
bool are_valid_indices(const int i_xyz[3]) const;
std::pair<double, std::array<int, 3>>
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const;
distance(Position r, Direction u, const int i_xyz[3]) const;
std::array<int, 3> get_indices(const double xyz[3]) const;
std::array<int, 3> get_indices(Position r) const;
std::array<double, 3>
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;
Position
get_local_position(Position r, const int i_xyz[3]) const;
int32_t& offset(int map, const int i_xyz[3]);
@ -207,9 +209,9 @@ public:
void to_hdf5_inner(hid_t group_id) const;
private:
std::array<int, 3> n_cells; //!< Number of cells along each axis
std::array<double, 3> lower_left; //!< Global lower-left corner of the lattice
std::array<double, 3> pitch; //!< Lattice tile width along each axis
std::array<int, 3> n_cells; //!< Number of cells along each axis
Position lower_left; //!< Global lower-left corner of the lattice
Position pitch; //!< Lattice tile width along each axis
// Convenience aliases
int &nx {n_cells[0]};
@ -233,12 +235,12 @@ public:
bool are_valid_indices(const int i_xyz[3]) const;
std::pair<double, std::array<int, 3>>
distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const;
distance(Position r, Direction u, const int i_xyz[3]) const;
std::array<int, 3> get_indices(const double xyz[3]) const;
std::array<int, 3> get_indices(Position r) const;
std::array<double, 3>
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;
Position
get_local_position(Position r, const int i_xyz[3]) const;
bool is_valid_index(int indx) const;
@ -251,9 +253,9 @@ public:
private:
int n_rings; //!< Number of radial tile positions
int n_axial; //!< Number of axial tile positions
std::array<double, 3> center; //!< Global center of lattice
Position center; //!< Global center of lattice
std::array<double, 2> pitch; //!< Lattice tile width and height
};
} // namespace openmc
#endif // LATTICE_H
#endif // OPENMC_LATTICE_H

9
src/math.F90
View file

@ -11,6 +11,7 @@ module math
public :: calc_pn
public :: calc_rn
public :: calc_zn
public :: calc_zn_rad
public :: evaluate_legendre
public :: rotate_angle
public :: maxwell_spectrum
@ -68,6 +69,14 @@ module math
real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2)
end subroutine calc_zn
pure subroutine calc_zn_rad(n, rho, zn_rad) bind(C, name='calc_zn_rad_c')
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), value, intent(in) :: rho
real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1)
end subroutine calc_zn_rad
subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c')
use ISO_C_BINDING
implicit none

27
src/math_functions.cpp
View file

@ -587,6 +587,33 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) {
}
void calc_zn_rad_c(int n, double rho, double zn_rad[]) {
// Calculate R_p0(rho) as Zn_p0(rho)
// Set up the array of the coefficients
double q = 0;
// R_00 is always 1
zn_rad[0] = 1;
// Fill in the rest of the array (Eq 3.8 and Eq 3.10 in Chong)
for (int p = 2; p <= n; p += 2) {
int index = int(p/2);
if (p == 2) {
// Setting up R_22 to calculate R_20 (Eq 3.10 in Chong)
double R_22 = rho * rho;
zn_rad[index] = 2 * R_22 - zn_rad[0];
} else {
double k1 = ((p + q) * (p - q) * (p - 2)) / 2.;
double k2 = 2 * p * (p - 1) * (p - 2);
double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2);
double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.;
zn_rad[index] =
((k2 * rho * rho + k3) * zn_rad[index-1] + k4 * zn_rad[index-2]) / k1;
}
}
}
void rotate_angle_c(double uvw[3], double mu, double* phi) {
// Copy original directional cosines

18
src/math_functions.h
View file

@ -91,6 +91,24 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
//==============================================================================
//! Calculate only the even order components of n-th order modified Zernike
//! polynomial moment with azimuthal dependency m = 0 for a given radial (rho)
//! location on the unit disk.
//!
//! Since m = 0, n could only be even orders. Z_q0 = R_q0
//!
//! See calc_zn_c for methodology.
//!
//! @param n The maximum order requested
//! @param rho The radial parameter to specify location on the unit disk
//! @param phi The angle parameter to specify location on the unit disk
//! @param zn_rad The requested moments of order 0 to n (inclusive)
//! evaluated at rho and phi when m = 0.
//==============================================================================
extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]);
//==============================================================================
//! Rotate the direction cosines through a polar angle whose cosine is mu and
//! through an azimuthal angle sampled uniformly.

8
src/output.F90
View file

@ -324,7 +324,11 @@ contains
! Determine overall generation and number of active generations
i = overall_generation()
n = i - n_inactive*gen_per_batch
if (current_batch > n_inactive) then
n = gen_per_batch*n_realizations + current_gen
else
n = 0
end if
! write out information about batch and generation
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
@ -357,7 +361,7 @@ contains
! Determine overall generation and number of active generations
i = current_batch*gen_per_batch
n = i - n_inactive*gen_per_batch
n = n_realizations*gen_per_batch
! write out information batch and option independent output
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &

63
src/position.cpp Normal file
View file

@ -0,0 +1,63 @@
#include "position.h"
namespace openmc {
//==============================================================================
// Position implementation
//==============================================================================
Position&
Position::operator+=(Position other)
{
x += other.x;
y += other.y;
z += other.z;
return *this;
}
Position&
Position::operator+=(double v)
{
x += v;
y += v;
z += v;
return *this;
}
Position&
Position::operator-=(Position other)
{
x -= other.x;
y -= other.y;
z -= other.z;
return *this;
}
Position&
Position::operator-=(double v)
{
x -= v;
y -= v;
z -= v;
return *this;
}
Position&
Position::operator*=(Position other)
{
x *= other.x;
y *= other.y;
z *= other.z;
return *this;
}
Position&
Position::operator*=(double v)
{
x *= v;
y *= v;
z *= v;
return *this;
}
} // namespace openmc

74
src/position.h Normal file
View file

@ -0,0 +1,74 @@
#ifndef OPENMC_POSITION_H
#define OPENMC_POSITION_H
namespace openmc {
//==============================================================================
//! Type representing a position in Cartesian coordinates
//==============================================================================
struct Position {
// Constructors
Position() = default;
Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { };
Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
// Unary operators
Position& operator+=(Position);
Position& operator+=(double);
Position& operator-=(Position);
Position& operator-=(double);
Position& operator*=(Position);
Position& operator*=(double);
const double& operator[](int i) const {
switch (i) {
case 0: return x;
case 1: return y;
case 2: return z;
}
}
double& operator[](int i) {
switch (i) {
case 0: return x;
case 1: return y;
case 2: return z;
}
}
// Other member functions
//! Dot product of two vectors
//! \param[in] other Vector to take dot product with
//! \result Resulting dot product
inline double dot(Position other) {
return x*other.x + y*other.y + z*other.z;
}
// Data members
double x = 0.;
double y = 0.;
double z = 0.;
};
// Binary operators
inline Position operator+(Position a, Position b) { return a += b; }
inline Position operator+(Position a, double b) { return a += b; }
inline Position operator+(double a, Position b) { return b += a; }
inline Position operator-(Position a, Position b) { return a -= b; }
inline Position operator-(Position a, double b) { return a -= b; }
inline Position operator-(double a, Position b) { return b -= a; }
inline Position operator*(Position a, Position b) { return a *= b; }
inline Position operator*(Position a, double b) { return a *= b; }
inline Position operator*(double a, Position b) { return b *= a; }
//==============================================================================
//! Type representing a vector direction in Cartesian coordinates
//==============================================================================
using Direction = Position;
} // namespace openmc
#endif // OPENMC_POSITION_H

20
src/simulation.F90
View file

@ -397,6 +397,13 @@ contains
end do
end if
! Increment n_realizations as would ordinarily be done in finalize_batch
if (reduce_tallies) then
n_realizations = n_realizations + 1
else
n_realizations = n_realizations + n_procs
end if
! Write message at end
if (current_batch == restart_batch) then
call write_message("Resuming simulation...", 6)
@ -460,6 +467,19 @@ contains
call initialize_source()
end if
! In restart, set the batch to begin from in order to reproduce the correct
! average keff (used in sampling the fission bank). Use n_realizations from
! the statepoint rather than n_inactive in case openmc_reset was called in
! the previous run.
if (restart_run) then
if (reduce_tallies) then
current_batch = restart_batch - n_realizations
else
current_batch = restart_batch - n_realizations*n_procs
end if
n_realizations = 0
end if
! Display header
if (master) then
if (run_mode == MODE_FIXEDSOURCE) then

7
src/state_point.F90
View file

@ -754,6 +754,9 @@ contains
end if
end if
! Read number of realizations for global tallies
call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.)
! Check to make sure source bank is present
if (path_source_point == path_state_point .and. .not. source_present) then
call fatal_error("Source bank must be contained in statepoint restart &
@ -767,10 +770,6 @@ contains
#else
if (master) then
#endif
! Read number of realizations for global tallies
call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.)
! Read global tally data
call read_dataset(global_tallies, file_id, "global_tallies")

483
src/surface.cpp
View file

@ -3,6 +3,7 @@
#include <array>
#include <cmath>
#include <sstream>
#include <utility>
#include "error.h"
#include "hdf5_interface.h"
@ -177,38 +178,33 @@ Surface::Surface(pugi::xml_node surf_node)
}
bool
Surface::sense(const double xyz[3], const double uvw[3]) const
Surface::sense(Position r, Direction u) const
{
// Evaluate the surface equation at the particle's coordinates to determine
// which side the particle is on.
const double f = evaluate(xyz);
const double f = evaluate(r);
// Check which side of surface the point is on.
if (std::abs(f) < FP_COINCIDENT) {
// Particle may be coincident with this surface. To determine the sense, we
// look at the direction of the particle relative to the surface normal (by
// default in the positive direction) via their dot product.
double norm[3];
normal(xyz, norm);
return uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0;
return u.dot(normal(r)) > 0.0;
}
return f > 0.0;
}
void
Surface::reflect(const double xyz[3], double uvw[3]) const
Direction
Surface::reflect(Position r, Direction u) const
{
// Determine projection of direction onto normal and squared magnitude of
// normal.
double norm[3];
normal(xyz, norm);
const double projection = norm[0]*uvw[0] + norm[1]*uvw[1] + norm[2]*uvw[2];
const double magnitude = norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2];
Direction n = normal(r);
const double projection = n.dot(u);
const double magnitude = n.dot(n);
// Reflect direction according to normal.
uvw[0] -= 2.0 * projection / magnitude * norm[0];
uvw[1] -= 2.0 * projection / magnitude * norm[1];
uvw[2] -= 2.0 * projection / magnitude * norm[2];
return u -= (2.0 * projection / magnitude) * n;
}
void
@ -261,33 +257,15 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node)
// The template parameter indicates the axis normal to the plane.
template<int i> double
axis_aligned_plane_evaluate(const double xyz[3], double offset)
axis_aligned_plane_distance(Position r, Direction u, bool coincident, double offset)
{
return xyz[i] - offset;
}
// The template parameter indicates the axis normal to the plane.
template<int i> double
axis_aligned_plane_distance(const double xyz[3], const double uvw[3],
bool coincident, double offset)
{
const double f = offset - xyz[i];
if (coincident or std::abs(f) < FP_COINCIDENT or uvw[i] == 0.0) return INFTY;
const double d = f / uvw[i];
const double f = offset - r[i];
if (coincident or std::abs(f) < FP_COINCIDENT or u[i] == 0.0) return INFTY;
const double d = f / u[i];
if (d < 0.0) return INFTY;
return d;
}
// The first template parameter indicates the axis normal to the plane. The
// other two parameters indicate the other two axes.
template<int i1, int i2, int i3> void
axis_aligned_plane_normal(const double xyz[3], double uvw[3])
{
uvw[i1] = 1.0;
uvw[i2] = 0.0;
uvw[i3] = 0.0;
}
//==============================================================================
// SurfaceXPlane implementation
//==============================================================================
@ -298,20 +276,19 @@ SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node)
read_coeffs(surf_node, id, x0);
}
inline double SurfaceXPlane::evaluate(const double xyz[3]) const
double SurfaceXPlane::evaluate(Position r) const
{
return axis_aligned_plane_evaluate<0>(xyz, x0);
return r.x - x0;
}
inline double SurfaceXPlane::distance(const double xyz[3], const double uvw[3],
bool coincident) const
double SurfaceXPlane::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_plane_distance<0>(xyz, uvw, coincident, x0);
return axis_aligned_plane_distance<0>(r, u, coincident, x0);
}
inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceXPlane::normal(Position r) const
{
axis_aligned_plane_normal<0, 1, 2>(xyz, uvw);
return {1., 0., 0.};
}
void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
@ -321,26 +298,24 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const
bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Direction& u) const
{
double other_norm[3];
other->normal(xyz, other_norm);
if (other_norm[0] == 1 and other_norm[1] == 0 and other_norm[2] == 0) {
xyz[0] = x0;
Direction other_n = other->normal(r);
if (other_n.x == 1 and other_n.y == 0 and other_n.z == 0) {
r.x = x0;
return false;
} else {
// Assume the partner is an YPlane (the only supported partner). Use the
// evaluate function to find y0, then adjust xyz and uvw for rotational
// evaluate function to find y0, then adjust position/Direction for rotational
// symmetry.
double xyz_test[3] {0, 0, 0};
double y0 = -other->evaluate(xyz_test);
xyz[1] = xyz[0] - x0 + y0;
xyz[0] = x0;
double y0 = -other->evaluate({0., 0., 0.});
r.y = r.x - x0 + y0;
r.x = x0;
double u = uvw[0];
uvw[0] = -uvw[1];
uvw[1] = u;
double ux = u.x;
u.x = -u.y;
u.y = ux;
return true;
}
@ -349,8 +324,7 @@ bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
BoundingBox
SurfaceXPlane::bounding_box() const
{
BoundingBox out {x0, x0, -INFTY, INFTY, -INFTY, INFTY};
return out;
return {x0, x0, -INFTY, INFTY, -INFTY, INFTY};
}
//==============================================================================
@ -363,20 +337,19 @@ SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node)
read_coeffs(surf_node, id, y0);
}
inline double SurfaceYPlane::evaluate(const double xyz[3]) const
double SurfaceYPlane::evaluate(Position r) const
{
return axis_aligned_plane_evaluate<1>(xyz, y0);
return r.y - y0;
}
inline double SurfaceYPlane::distance(const double xyz[3], const double uvw[3],
bool coincident) const
double SurfaceYPlane::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_plane_distance<1>(xyz, uvw, coincident, y0);
return axis_aligned_plane_distance<1>(r, u, coincident, y0);
}
inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceYPlane::normal(Position r) const
{
axis_aligned_plane_normal<1, 0, 2>(xyz, uvw);
return {0., 1., 0.};
}
void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
@ -386,27 +359,25 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const
bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Direction& u) const
{
double other_norm[3];
other->normal(xyz, other_norm);
if (other_norm[0] == 0 and other_norm[1] == 1 and other_norm[2] == 0) {
Direction other_n = other->normal(r);
if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) {
// The periodic partner is also aligned along y. Just change the y coord.
xyz[1] = y0;
r.y = y0;
return false;
} else {
// Assume the partner is an XPlane (the only supported partner). Use the
// evaluate function to find x0, then adjust xyz and uvw for rotational
// evaluate function to find x0, then adjust position/Direction for rotational
// symmetry.
double xyz_test[3] {0, 0, 0};
double x0 = -other->evaluate(xyz_test);
xyz[0] = xyz[1] - y0 + x0;
xyz[1] = y0;
double x0 = -other->evaluate({0., 0., 0.});
r.x = r.y - y0 + x0;
r.y = y0;
double u = uvw[0];
uvw[0] = uvw[1];
uvw[1] = -u;
double ux = u.x;
u.x = u.y;
u.y = -ux;
return true;
}
@ -415,8 +386,7 @@ bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
BoundingBox
SurfaceYPlane::bounding_box() const
{
BoundingBox out {-INFTY, INFTY, y0, y0, -INFTY, INFTY};
return out;
return {-INFTY, INFTY, y0, y0, -INFTY, INFTY};
}
//==============================================================================
@ -429,20 +399,19 @@ SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node)
read_coeffs(surf_node, id, z0);
}
inline double SurfaceZPlane::evaluate(const double xyz[3]) const
double SurfaceZPlane::evaluate(Position r) const
{
return axis_aligned_plane_evaluate<2>(xyz, z0);
return r.z - z0;
}
inline double SurfaceZPlane::distance(const double xyz[3], const double uvw[3],
bool coincident) const
double SurfaceZPlane::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_plane_distance<2>(xyz, uvw, coincident, z0);
return axis_aligned_plane_distance<2>(r, u, coincident, z0);
}
inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceZPlane::normal(Position r) const
{
axis_aligned_plane_normal<2, 0, 1>(xyz, uvw);
return {0., 0., 1.};
}
void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
@ -452,19 +421,18 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const
bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r,
Direction& u) const
{
// Assume the other plane is aligned along z. Just change the z coord.
xyz[2] = z0;
r.z = z0;
return false;
}
BoundingBox
SurfaceZPlane::bounding_box() const
{
BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, z0, z0};
return out;
return {-INFTY, INFTY, -INFTY, INFTY, z0, z0};
}
//==============================================================================
@ -478,17 +446,16 @@ SurfacePlane::SurfacePlane(pugi::xml_node surf_node)
}
double
SurfacePlane::evaluate(const double xyz[3]) const
SurfacePlane::evaluate(Position r) const
{
return A*xyz[0] + B*xyz[1] + C*xyz[2] - D;
return A*r.x + B*r.y + C*r.z - D;
}
double
SurfacePlane::distance(const double xyz[3], const double uvw[3],
bool coincident) const
SurfacePlane::distance(Position r, Direction u, bool coincident) const
{
const double f = A*xyz[0] + B*xyz[1] + C*xyz[2] - D;
const double projection = A*uvw[0] + B*uvw[1] + C*uvw[2];
const double f = A*r.x + B*r.y + C*r.z - D;
const double projection = A*u.x + B*u.y + C*u.z;
if (coincident or std::abs(f) < FP_COINCIDENT or projection == 0.0) {
return INFTY;
} else {
@ -498,12 +465,10 @@ SurfacePlane::distance(const double xyz[3], const double uvw[3],
}
}
void
SurfacePlane::normal(const double xyz[3], double uvw[3]) const
Direction
SurfacePlane::normal(Position r) const
{
uvw[0] = A;
uvw[1] = B;
uvw[2] = C;
return {A, B, C};
}
void SurfacePlane::to_hdf5_inner(hid_t group_id) const
@ -513,18 +478,18 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const
bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r,
Direction& u) const
{
// This function assumes the other plane shares this plane's normal direction.
// Determine the distance to intersection.
double d = evaluate(xyz) / (A*A + B*B + C*C);
double d = evaluate(r) / (A*A + B*B + C*C);
// Move the particle that distance along the normal vector.
xyz[0] -= d * A;
xyz[1] -= d * B;
xyz[2] -= d * C;
r.x -= d * A;
r.y -= d * B;
r.z -= d * C;
return false;
}
@ -532,8 +497,7 @@ bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3],
BoundingBox
SurfacePlane::bounding_box() const
{
BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
return out;
return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
}
//==============================================================================
@ -544,28 +508,28 @@ SurfacePlane::bounding_box() const
// cylinder. offset1 and offset2 should correspond with i1 and i2,
// respectively.
template<int i1, int i2> double
axis_aligned_cylinder_evaluate(const double xyz[3], double offset1,
axis_aligned_cylinder_evaluate(Position r, double offset1,
double offset2, double radius)
{
const double xyz1 = xyz[i1] - offset1;
const double xyz2 = xyz[i2] - offset2;
return xyz1*xyz1 + xyz2*xyz2 - radius*radius;
const double r1 = r[i1] - offset1;
const double r2 = r[i2] - offset2;
return r1*r1 + r2*r2 - radius*radius;
}
// The first template parameter indicates which axis the cylinder is aligned to.
// The other two parameters indicate the other two axes. offset1 and offset2
// should correspond with i2 and i3, respectively.
template<int i1, int i2, int i3> double
axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3],
axis_aligned_cylinder_distance(Position r, Direction u,
bool coincident, double offset1, double offset2, double radius)
{
const double a = 1.0 - uvw[i1]*uvw[i1]; // u^2 + v^2
const double a = 1.0 - u[i1]*u[i1]; // u^2 + v^2
if (a == 0.0) return INFTY;
const double xyz2 = xyz[i2] - offset1;
const double xyz3 = xyz[i3] - offset2;
const double k = xyz2 * uvw[i2] + xyz3 * uvw[i3];
const double c = xyz2*xyz2 + xyz3*xyz3 - radius*radius;
const double r2 = r[i2] - offset1;
const double r3 = r[i3] - offset2;
const double k = r2 * u[i2] + r3 * u[i3];
const double c = r2*r2 + r3*r3 - radius*radius;
const double quad = k*k - a*c;
if (quad < 0.0) {
@ -601,13 +565,14 @@ axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3],
// The first template parameter indicates which axis the cylinder is aligned to.
// The other two parameters indicate the other two axes. offset1 and offset2
// should correspond with i2 and i3, respectively.
template<int i1, int i2, int i3> void
axis_aligned_cylinder_normal(const double xyz[3], double uvw[3], double offset1,
double offset2)
template<int i1, int i2, int i3> Direction
axis_aligned_cylinder_normal(Position r, double offset1, double offset2)
{
uvw[i2] = 2.0 * (xyz[i2] - offset1);
uvw[i3] = 2.0 * (xyz[i3] - offset2);
uvw[i1] = 0.0;
Direction u;
u[i2] = 2.0 * (r[i2] - offset1);
u[i3] = 2.0 * (r[i3] - offset2);
u[i1] = 0.0;
return u;
}
//==============================================================================
@ -617,31 +582,30 @@ axis_aligned_cylinder_normal(const double xyz[3], double uvw[3], double offset1,
SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node)
: Surface(surf_node)
{
read_coeffs(surf_node, id, y0, z0, r);
read_coeffs(surf_node, id, y0, z0, radius);
}
inline double SurfaceXCylinder::evaluate(const double xyz[3]) const
double SurfaceXCylinder::evaluate(Position r) const
{
return axis_aligned_cylinder_evaluate<1, 2>(xyz, y0, z0, r);
return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, radius);
}
inline double SurfaceXCylinder::distance(const double xyz[3],
const double uvw[3], bool coincident) const
double SurfaceXCylinder::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cylinder_distance<0, 1, 2>(xyz, uvw, coincident, y0, z0,
r);
return axis_aligned_cylinder_distance<0, 1, 2>(r, u, coincident, y0, z0,
radius);
}
inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceXCylinder::normal(Position r) const
{
axis_aligned_cylinder_normal<0, 1, 2>(xyz, uvw, y0, z0);
return axis_aligned_cylinder_normal<0, 1, 2>(r, y0, z0);
}
void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-cylinder", false);
std::array<double, 3> coeffs {{y0, z0, r}};
std::array<double, 3> coeffs {{y0, z0, radius}};
write_dataset(group_id, "coefficients", coeffs);
}
@ -652,30 +616,29 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node)
: Surface(surf_node)
{
read_coeffs(surf_node, id, x0, z0, r);
read_coeffs(surf_node, id, x0, z0, radius);
}
inline double SurfaceYCylinder::evaluate(const double xyz[3]) const
double SurfaceYCylinder::evaluate(Position r) const
{
return axis_aligned_cylinder_evaluate<0, 2>(xyz, x0, z0, r);
return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, radius);
}
inline double SurfaceYCylinder::distance(const double xyz[3],
const double uvw[3], bool coincident) const
double SurfaceYCylinder::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cylinder_distance<1, 0, 2>(xyz, uvw, coincident, x0, z0,
r);
return axis_aligned_cylinder_distance<1, 0, 2>(r, u, coincident, x0, z0,
radius);
}
inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceYCylinder::normal(Position r) const
{
axis_aligned_cylinder_normal<1, 0, 2>(xyz, uvw, x0, z0);
return axis_aligned_cylinder_normal<1, 0, 2>(r, x0, z0);
}
void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "y-cylinder", false);
std::array<double, 3> coeffs {{x0, z0, r}};
std::array<double, 3> coeffs {{x0, z0, radius}};
write_dataset(group_id, "coefficients", coeffs);
}
@ -686,30 +649,29 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node)
: Surface(surf_node)
{
read_coeffs(surf_node, id, x0, y0, r);
read_coeffs(surf_node, id, x0, y0, radius);
}
inline double SurfaceZCylinder::evaluate(const double xyz[3]) const
double SurfaceZCylinder::evaluate(Position r) const
{
return axis_aligned_cylinder_evaluate<0, 1>(xyz, x0, y0, r);
return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, radius);
}
inline double SurfaceZCylinder::distance(const double xyz[3],
const double uvw[3], bool coincident) const
double SurfaceZCylinder::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cylinder_distance<2, 0, 1>(xyz, uvw, coincident, x0, y0,
r);
return axis_aligned_cylinder_distance<2, 0, 1>(r, u, coincident, x0, y0,
radius);
}
inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceZCylinder::normal(Position r) const
{
axis_aligned_cylinder_normal<2, 0, 1>(xyz, uvw, x0, y0);
return axis_aligned_cylinder_normal<2, 0, 1>(r, x0, y0);
}
void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "z-cylinder", false);
std::array<double, 3> coeffs {{x0, y0, r}};
std::array<double, 3> coeffs {{x0, y0, radius}};
write_dataset(group_id, "coefficients", coeffs);
}
@ -720,25 +682,24 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node)
: Surface(surf_node)
{
read_coeffs(surf_node, id, x0, y0, z0, r);
read_coeffs(surf_node, id, x0, y0, z0, radius);
}
double SurfaceSphere::evaluate(const double xyz[3]) const
double SurfaceSphere::evaluate(Position r) const
{
const double x = xyz[0] - x0;
const double y = xyz[1] - y0;
const double z = xyz[2] - z0;
return x*x + y*y + z*z - r*r;
const double x = r.x - x0;
const double y = r.y - y0;
const double z = r.z - z0;
return x*x + y*y + z*z - radius*radius;
}
double SurfaceSphere::distance(const double xyz[3], const double uvw[3],
bool coincident) const
double SurfaceSphere::distance(Position r, Direction u, bool coincident) const
{
const double x = xyz[0] - x0;
const double y = xyz[1] - y0;
const double z = xyz[2] - z0;
const double k = x*uvw[0] + y*uvw[1] + z*uvw[2];
const double c = x*x + y*y + z*z - r*r;
const double x = r.x - x0;
const double y = r.y - y0;
const double z = r.z - z0;
const double k = x*u.x + y*u.y + z*u.z;
const double c = x*x + y*y + z*z - radius*radius;
const double quad = k*k - c;
if (quad < 0.0) {
@ -770,17 +731,15 @@ double SurfaceSphere::distance(const double xyz[3], const double uvw[3],
}
}
inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceSphere::normal(Position r) const
{
uvw[0] = 2.0 * (xyz[0] - x0);
uvw[1] = 2.0 * (xyz[1] - y0);
uvw[2] = 2.0 * (xyz[2] - z0);
return {2.0*(r.x - x0), 2.0*(r.y - y0), 2.0*(r.z - z0)};
}
void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "sphere", false);
std::array<double, 4> coeffs {{x0, y0, z0, r}};
std::array<double, 4> coeffs {{x0, y0, z0, radius}};
write_dataset(group_id, "coefficients", coeffs);
}
@ -792,30 +751,30 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
// The other two parameters indicate the other two axes. offset1, offset2,
// and offset3 should correspond with i1, i2, and i3, respectively.
template<int i1, int i2, int i3> double
axis_aligned_cone_evaluate(const double xyz[3], double offset1,
axis_aligned_cone_evaluate(Position r, double offset1,
double offset2, double offset3, double radius_sq)
{
const double xyz1 = xyz[i1] - offset1;
const double xyz2 = xyz[i2] - offset2;
const double xyz3 = xyz[i3] - offset3;
return xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1;
const double r1 = r[i1] - offset1;
const double r2 = r[i2] - offset2;
const double r3 = r[i3] - offset3;
return r2*r2 + r3*r3 - radius_sq*r1*r1;
}
// The first template parameter indicates which axis the cone is aligned to.
// The other two parameters indicate the other two axes. offset1, offset2,
// and offset3 should correspond with i1, i2, and i3, respectively.
template<int i1, int i2, int i3> double
axis_aligned_cone_distance(const double xyz[3], const double uvw[3],
axis_aligned_cone_distance(Position r, Direction u,
bool coincident, double offset1, double offset2, double offset3,
double radius_sq)
{
const double xyz1 = xyz[i1] - offset1;
const double xyz2 = xyz[i2] - offset2;
const double xyz3 = xyz[i3] - offset3;
const double a = uvw[i2]*uvw[i2] + uvw[i3]*uvw[i3]
- radius_sq*uvw[i1]*uvw[i1];
const double k = xyz2*uvw[i2] + xyz3*uvw[i3] - radius_sq*xyz1*uvw[i1];
const double c = xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1;
const double r1 = r[i1] - offset1;
const double r2 = r[i2] - offset2;
const double r3 = r[i3] - offset3;
const double a = u[i2]*u[i2] + u[i3]*u[i3]
- radius_sq*u[i1]*u[i1];
const double k = r2*u[i2] + r3*u[i3] - radius_sq*r1*u[i1];
const double c = r2*r2 + r3*r3 - radius_sq*r1*r1;
double quad = k*k - a*c;
double d;
@ -858,13 +817,15 @@ axis_aligned_cone_distance(const double xyz[3], const double uvw[3],
// The first template parameter indicates which axis the cone is aligned to.
// The other two parameters indicate the other two axes. offset1, offset2,
// and offset3 should correspond with i1, i2, and i3, respectively.
template<int i1, int i2, int i3> void
axis_aligned_cone_normal(const double xyz[3], double uvw[3], double offset1,
double offset2, double offset3, double radius_sq)
template<int i1, int i2, int i3> Direction
axis_aligned_cone_normal(Position r, double offset1, double offset2,
double offset3, double radius_sq)
{
uvw[i1] = -2.0 * radius_sq * (xyz[i1] - offset1);
uvw[i2] = 2.0 * (xyz[i2] - offset2);
uvw[i3] = 2.0 * (xyz[i3] - offset3);
Direction u;
u[i1] = -2.0 * radius_sq * (r[i1] - offset1);
u[i2] = 2.0 * (r[i2] - offset2);
u[i3] = 2.0 * (r[i3] - offset3);
return u;
}
//==============================================================================
@ -874,30 +835,29 @@ axis_aligned_cone_normal(const double xyz[3], double uvw[3], double offset1,
SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node)
: Surface(surf_node)
{
read_coeffs(surf_node, id, x0, y0, z0, r_sq);
read_coeffs(surf_node, id, x0, y0, z0, radius_sq);
}
inline double SurfaceXCone::evaluate(const double xyz[3]) const
double SurfaceXCone::evaluate(Position r) const
{
return axis_aligned_cone_evaluate<0, 1, 2>(xyz, x0, y0, z0, r_sq);
return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, radius_sq);
}
inline double SurfaceXCone::distance(const double xyz[3],
const double uvw[3], bool coincident) const
double SurfaceXCone::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cone_distance<0, 1, 2>(xyz, uvw, coincident, x0, y0, z0,
r_sq);
return axis_aligned_cone_distance<0, 1, 2>(r, u, coincident, x0, y0, z0,
radius_sq);
}
inline void SurfaceXCone::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceXCone::normal(Position r) const
{
axis_aligned_cone_normal<0, 1, 2>(xyz, uvw, x0, y0, z0, r_sq);
return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, radius_sq);
}
void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
std::array<double, 4> coeffs {{x0, y0, z0, radius_sq}};
write_dataset(group_id, "coefficients", coeffs);
}
@ -908,30 +868,29 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node)
: Surface(surf_node)
{
read_coeffs(surf_node, id, x0, y0, z0, r_sq);
read_coeffs(surf_node, id, x0, y0, z0, radius_sq);
}
inline double SurfaceYCone::evaluate(const double xyz[3]) const
double SurfaceYCone::evaluate(Position r) const
{
return axis_aligned_cone_evaluate<1, 0, 2>(xyz, y0, x0, z0, r_sq);
return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, radius_sq);
}
inline double SurfaceYCone::distance(const double xyz[3],
const double uvw[3], bool coincident) const
double SurfaceYCone::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cone_distance<1, 0, 2>(xyz, uvw, coincident, y0, x0, z0,
r_sq);
return axis_aligned_cone_distance<1, 0, 2>(r, u, coincident, y0, x0, z0,
radius_sq);
}
inline void SurfaceYCone::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceYCone::normal(Position r) const
{
axis_aligned_cone_normal<1, 0, 2>(xyz, uvw, y0, x0, z0, r_sq);
return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, radius_sq);
}
void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "y-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
std::array<double, 4> coeffs {{x0, y0, z0, radius_sq}};
write_dataset(group_id, "coefficients", coeffs);
}
@ -942,30 +901,29 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node)
: Surface(surf_node)
{
read_coeffs(surf_node, id, x0, y0, z0, r_sq);
read_coeffs(surf_node, id, x0, y0, z0, radius_sq);
}
inline double SurfaceZCone::evaluate(const double xyz[3]) const
double SurfaceZCone::evaluate(Position r) const
{
return axis_aligned_cone_evaluate<2, 0, 1>(xyz, z0, x0, y0, r_sq);
return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, radius_sq);
}
inline double SurfaceZCone::distance(const double xyz[3],
const double uvw[3], bool coincident) const
double SurfaceZCone::distance(Position r, Direction u, bool coincident) const
{
return axis_aligned_cone_distance<2, 0, 1>(xyz, uvw, coincident, z0, x0, y0,
r_sq);
return axis_aligned_cone_distance<2, 0, 1>(r, u, coincident, z0, x0, y0,
radius_sq);
}
inline void SurfaceZCone::normal(const double xyz[3], double uvw[3]) const
Direction SurfaceZCone::normal(Position r) const
{
axis_aligned_cone_normal<2, 0, 1>(xyz, uvw, z0, x0, y0, r_sq);
return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, radius_sq);
}
void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "z-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
std::array<double, 4> coeffs {{x0, y0, z0, radius_sq}};
write_dataset(group_id, "coefficients", coeffs);
}
@ -980,26 +938,25 @@ SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node)
}
double
SurfaceQuadric::evaluate(const double xyz[3]) const
SurfaceQuadric::evaluate(Position r) const
{
const double &x = xyz[0];
const double &y = xyz[1];
const double &z = xyz[2];
const double x = r.x;
const double y = r.y;
const double z = r.z;
return x*(A*x + D*y + G) +
y*(B*y + E*z + H) +
z*(C*z + F*x + J) + K;
}
double
SurfaceQuadric::distance(const double xyz[3],
const double uvw[3], bool coincident) const
SurfaceQuadric::distance(Position r, Direction ang, bool coincident) const
{
const double &x = xyz[0];
const double &y = xyz[1];
const double &z = xyz[2];
const double &u = uvw[0];
const double &v = uvw[1];
const double &w = uvw[2];
const double &x = r.x;
const double &y = r.y;
const double &z = r.z;
const double &u = ang.x;
const double &v = ang.y;
const double &w = ang.z;
const double a = A*u*u + B*v*v + C*w*w + D*u*v + E*v*w + F*u*w;
const double k = (A*u*x + B*v*y + C*w*z + 0.5*(D*(u*y + v*x) +
@ -1045,15 +1002,15 @@ SurfaceQuadric::distance(const double xyz[3],
return d;
}
void
SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const
Direction
SurfaceQuadric::normal(Position r) const
{
const double &x = xyz[0];
const double &y = xyz[1];
const double &z = xyz[2];
uvw[0] = 2.0*A*x + D*y + F*z + G;
uvw[1] = 2.0*B*y + D*x + E*z + H;
uvw[2] = 2.0*C*z + E*y + F*x + J;
const double &x = r.x;
const double &y = r.y;
const double &z = r.z;
return {2.0*A*x + D*y + F*z + G,
2.0*B*y + D*x + E*z + H,
2.0*C*z + E*y + F*x + J};
}
void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
@ -1263,10 +1220,24 @@ extern "C" {
int surface_bc(Surface *surf) {return surf->bc;}
void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
{surf->reflect(xyz, uvw);}
{
Position r {xyz};
Direction u {uvw};
u = surf->reflect(r, u);
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
}
void surface_normal(Surface *surf, double xyz[3], double uvw[3])
{return surf->normal(xyz, uvw);}
{
Position r {xyz};
Direction u = surf->normal(r);
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
}
void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);}
@ -1275,7 +1246,21 @@ extern "C" {
bool
surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
double uvw[3])
{return surf->periodic_translate(other, xyz, uvw);}
{
Position r {xyz};
Direction u {uvw};
bool rotational = surf->periodic_translate(other, r, u);
// Copy back to arrays
xyz[0] = r.x;
xyz[1] = r.y;
xyz[2] = r.z;
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
return rotational;
}
void free_memory_surfaces_c()
{

189
src/surface.h
View file

@ -1,5 +1,5 @@
#ifndef SURFACE_H
#define SURFACE_H
#ifndef OPENMC_SURFACE_H
#define OPENMC_SURFACE_H
#include <map>
#include <limits> // For numeric_limits
@ -9,6 +9,7 @@
#include "pugixml.hpp"
#include "constants.h"
#include "position.h"
namespace openmc {
@ -65,42 +66,40 @@ public:
virtual ~Surface() {}
//! Determine which side of a surface a point lies on.
//! @param xyz[3] The 3D Cartesian coordinate of a point.
//! @param uvw[3] A direction used to "break ties" and pick a sense when the
//! \param r The 3D Cartesian coordinate of a point.
//! \param u A direction used to "break ties" and pick a sense when the
//! point is very close to the surface.
//! @return true if the point is on the "positive" side of the surface and
//! \return true if the point is on the "positive" side of the surface and
//! false otherwise.
bool sense(const double xyz[3], const double uvw[3]) const;
bool sense(Position r, Direction u) const;
//! Determine the direction of a ray reflected from the surface.
//! @param xyz[3] The point at which the ray is incident.
//! @param uvw[3] A direction. This is both an input and an output parameter.
//! It specifies the icident direction on input and the reflected direction
//! on output.
void reflect(const double xyz[3], double uvw[3]) const;
//! \param[in] r The point at which the ray is incident.
//! \param[in] u Incident direction of the ray
//! \return Outgoing direction of the ray
Direction 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
//! function evaluates that mathematical function.
//! @param xyz[3] A 3D Cartesian coordinate.
virtual double evaluate(const double xyz[3]) const = 0;
//! \param r A 3D Cartesian coordinate.
virtual double evaluate(Position r) const = 0;
//! Compute the distance between a point and the surface along a ray.
//! @param xyz[3] A 3D Cartesian coordinate.
//! @param uvw[3] The direction of the ray.
//! @param coincident A hint to the code that the given point should lie
//! \param r A 3D Cartesian coordinate.
//! \param u The direction of the ray.
//! \param coincident A hint to the code that the given point should lie
//! exactly on the surface.
virtual double distance(const double xyz[3], const double uvw[3],
bool coincident) const = 0;
virtual double distance(Position r, Direction u, bool coincident) const = 0;
//! Compute the local outward normal direction of the surface.
//! @param xyz[3] A 3D Cartesian coordinate.
//! @param uvw[3] This output argument provides the normal.
virtual void normal(const double xyz[3], double uvw[3]) const = 0;
//! \param r A 3D Cartesian coordinate.
//! \return Normal direction
virtual Direction normal(Position r) const = 0;
//! Write all information needed to reconstruct the surface to an HDF5 group.
//! @param group_id An HDF5 group id.
//! \param group_id An HDF5 group id.
//TODO: this probably needs to include i_periodic for PeriodicSurface
void to_hdf5(hid_t group_id) const;
@ -124,15 +123,15 @@ public:
explicit PeriodicSurface(pugi::xml_node surf_node);
//! Translate a particle onto this surface from a periodic partner surface.
//! @param other A pointer to the partner surface in this periodic BC.
//! @param xyz[3] A point on the partner surface that will be translated onto
//! \param other A pointer to the partner surface in this periodic BC.
//! \param r A point on the partner surface that will be translated onto
//! this surface.
//! @param uvw[3] A direction that will be rotated for systems with rotational
//! \param u A direction that will be rotated for systems with rotational
//! periodicity.
//! @return true if this surface and its partner make a rotationally-periodic
//! \return true if this surface and its partner make a rotationally-periodic
//! boundary condition.
virtual bool periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const = 0;
virtual bool periodic_translate(const PeriodicSurface *other, Position& r,
Direction& u) const = 0;
//! Get the bounding box for this surface.
virtual BoundingBox bounding_box() const = 0;
@ -149,12 +148,11 @@ class SurfaceXPlane : public PeriodicSurface
double x0;
public:
explicit SurfaceXPlane(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3], bool coincident)
const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
BoundingBox bounding_box() const;
};
@ -170,12 +168,11 @@ class SurfaceYPlane : public PeriodicSurface
double y0;
public:
explicit SurfaceYPlane(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
BoundingBox bounding_box() const;
};
@ -191,12 +188,11 @@ class SurfaceZPlane : public PeriodicSurface
double z0;
public:
explicit SurfaceZPlane(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
BoundingBox bounding_box() const;
};
@ -212,12 +208,11 @@ class SurfacePlane : public PeriodicSurface
double A, B, C, D;
public:
explicit SurfacePlane(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3], bool coincident)
const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
const;
BoundingBox bounding_box() const;
};
@ -231,13 +226,12 @@ public:
class SurfaceXCylinder : public Surface
{
double y0, z0, r;
double y0, z0, radius;
public:
explicit SurfaceXCylinder(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -250,13 +244,12 @@ public:
class SurfaceYCylinder : public Surface
{
double x0, z0, r;
double x0, z0, radius;
public:
explicit SurfaceYCylinder(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -269,13 +262,12 @@ public:
class SurfaceZCylinder : public Surface
{
double x0, y0, r;
double x0, y0, radius;
public:
explicit SurfaceZCylinder(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -288,13 +280,12 @@ public:
class SurfaceSphere : public Surface
{
double x0, y0, z0, r;
double x0, y0, z0, radius;
public:
explicit SurfaceSphere(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -307,13 +298,12 @@ public:
class SurfaceXCone : public Surface
{
double x0, y0, z0, r_sq;
double x0, y0, z0, radius_sq;
public:
explicit SurfaceXCone(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -326,13 +316,12 @@ public:
class SurfaceYCone : public Surface
{
double x0, y0, z0, r_sq;
double x0, y0, z0, radius_sq;
public:
explicit SurfaceYCone(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -345,13 +334,12 @@ public:
class SurfaceZCone : public Surface
{
double x0, y0, z0, r_sq;
double x0, y0, z0, radius_sq;
public:
explicit SurfaceZCone(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -367,10 +355,9 @@ class SurfaceQuadric : public Surface
double A, B, C, D, E, F, G, H, J, K;
public:
explicit SurfaceQuadric(pugi::xml_node surf_node);
double evaluate(const double xyz[3]) const;
double distance(const double xyz[3], const double uvw[3],
bool coincident) const;
void normal(const double xyz[3], double uvw[3]) const;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
};
@ -379,20 +366,20 @@ public:
//==============================================================================
extern "C" {
Surface* surface_pointer(int surf_ind);
int surface_id(Surface *surf);
int surface_bc(Surface *surf);
bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
double surface_distance(Surface *surf, double xyz[3], double uvw[3],
bool coincident);
void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
void surface_to_hdf5(Surface *surf, hid_t group);
int surface_i_periodic(PeriodicSurface *surf);
bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
double xyz[3], double uvw[3]);
void free_memory_surfaces_c();
Surface* surface_pointer(int surf_ind);
int surface_id(Surface *surf);
int surface_bc(Surface *surf);
bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
double surface_distance(Surface *surf, double xyz[3], double uvw[3],
bool coincident);
void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
void surface_to_hdf5(Surface *surf, hid_t group);
int surface_i_periodic(PeriodicSurface *surf);
bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
double xyz[3], double uvw[3]);
void free_memory_surfaces_c();
}
} // namespace openmc
#endif // SURFACE_H
#endif // OPENMC_SURFACE_H

4
src/tallies/tally_filter.F90
View file

@ -91,6 +91,8 @@ contains
type_ = 'zernike'
type is (ParticleFilter)
type_ = 'particle'
type is (ZernikeRadialFilter)
type_ = 'zernikeradial'
end select
! Convert Fortran string to null-terminated C string. We assume the
@ -164,6 +166,8 @@ contains
allocate(UniverseFilter :: filters(index) % obj)
case ('zernike')
allocate(ZernikeFilter :: filters(index) % obj)
case ('zernikeradial')
allocate(ZernikeRadialFilter :: filters(index) % obj)
case default
err = E_UNASSIGNED
call set_errmsg("Unknown filter type: " // trim(type_))

142
src/tallies/tally_filter_zernike.F90
View file

@ -5,7 +5,7 @@ module tally_filter_zernike
use constants
use error
use hdf5_interface
use math, only: calc_zn
use math, only: calc_zn, calc_zn_rad
use particle_header, only: Particle
use string, only: to_str
use tally_filter_header
@ -13,6 +13,10 @@ module tally_filter_zernike
implicit none
private
public :: openmc_zernike_filter_get_order
public :: openmc_zernike_filter_get_params
public :: openmc_zernike_filter_set_order
public :: openmc_zernike_filter_set_params
!===============================================================================
! ZERNIKEFILTER gives Zernike polynomial moments of a particle's position
@ -24,19 +28,43 @@ module tally_filter_zernike
real(8) :: y
real(8) :: r
contains
procedure :: from_xml
procedure :: get_all_bins
procedure :: to_statepoint
procedure :: text_label
procedure :: calc_n_bins => calc_n_bins_zn
procedure :: from_xml => from_xml_zn
procedure :: get_all_bins => get_all_bins_zn
procedure :: to_statepoint => to_statepoint_zn
procedure :: text_label => text_label_zn
end type ZernikeFilter
!===============================================================================
! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a
! particle's position
!===============================================================================
type, public, extends(ZernikeFilter) :: ZernikeRadialFilter
contains
procedure :: calc_n_bins => calc_n_bins_zn_rad
! Inherit from_xml from ZernikeFilter
procedure :: get_all_bins => get_all_bins_zn_rad
procedure :: to_statepoint => to_statepoint_zn_rad
procedure :: text_label => text_label_zn_rad
end type ZernikeRadialFilter
contains
!===============================================================================
! ZernikeFilter methods
!===============================================================================
subroutine from_xml(this, node)
function calc_n_bins_zn(this) result(n_bins)
class(ZernikeFilter), intent(in) :: this
integer :: n
integer :: n_bins
n = this % order
n_bins = ((n+1) * (n+2))/2
end function calc_n_bins_zn
subroutine from_xml_zn(this, node)
class(ZernikeFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
@ -50,10 +78,10 @@ contains
! Get specified order
call get_node_value(node, "order", n)
this % order = n
this % n_bins = ((n + 1)*(n + 2))/2
end subroutine from_xml
this % n_bins = this % calc_n_bins()
end subroutine from_xml_zn
subroutine get_all_bins(this, p, estimator, match)
subroutine get_all_bins_zn(this, p, estimator, match)
class(ZernikeFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
@ -67,18 +95,20 @@ contains
x = p % coord(1) % xyz(1) - this % x
y = p % coord(1) % xyz(2) - this % y
r = sqrt(x*x + y*y)/this % r
theta = atan2(y, x)
if (r <= ONE) then
theta = atan2(y, x)
! Get moments for Zernike polynomial orders 0..n
call calc_zn(this % order, r, theta, zn)
! Get moments for Zernike polynomial orders 0..n
call calc_zn(this % order, r, theta, zn)
do i = 1, this % n_bins
call match % bins % push_back(i)
call match % weights % push_back(zn(i))
end do
endif
end subroutine get_all_bins_zn
do i = 1, this % n_bins
call match % bins % push_back(i)
call match % weights % push_back(zn(i))
end do
end subroutine get_all_bins
subroutine to_statepoint(this, filter_group)
subroutine to_statepoint_zn(this, filter_group)
class(ZernikeFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
@ -88,9 +118,9 @@ contains
call write_dataset(filter_group, "x", this % x)
call write_dataset(filter_group, "y", this % y)
call write_dataset(filter_group, "r", this % r)
end subroutine to_statepoint
end subroutine to_statepoint_zn
function text_label(this, bin) result(label)
function text_label_zn(this, bin) result(label)
class(ZernikeFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
@ -108,7 +138,66 @@ contains
exit
end if
end do
end function text_label
end function text_label_zn
!===============================================================================
! ZernikeRadialFilter methods
!===============================================================================
function calc_n_bins_zn_rad(this) result(num_n_bins)
class(ZernikeRadialFilter), intent(in) :: this
integer :: n
integer :: num_n_bins
n = this % order
num_n_bins = n/2 + 1
end function calc_n_bins_zn_rad
subroutine get_all_bins_zn_rad(this, p, estimator, match)
class(ZernikeRadialFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer :: i
real(8) :: x, y, r
real(C_DOUBLE) :: zn_rad(this % n_bins)
! Determine normalized (r,theta) positions
x = p % coord(1) % xyz(1) - this % x
y = p % coord(1) % xyz(2) - this % y
r = sqrt(x*x + y*y)/this % r
if (r <= ONE) then
! Get moments for even order Zernike polynomial orders 0..n
call calc_zn_rad(this % order, r, zn_rad)
do i = 1, this % n_bins
call match % bins % push_back(i)
call match % weights % push_back(zn_rad(i))
end do
endif
end subroutine get_all_bins_zn_rad
subroutine to_statepoint_zn_rad(this, filter_group)
class(ZernikeRadialFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
call write_dataset(filter_group, "type", "zernikeradial")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "order", this % order)
call write_dataset(filter_group, "x", this % x)
call write_dataset(filter_group, "y", this % y)
call write_dataset(filter_group, "r", this % r)
end subroutine to_statepoint_zn_rad
function text_label_zn_rad(this, bin) result(label)
class(ZernikeRadialFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
label = "Zernike expansion, Z" // trim(to_str(2*(bin-1))) // ",0"
end function text_label_zn_rad
!===============================================================================
! C API FUNCTIONS
@ -125,6 +214,8 @@ contains
select type (f => filters(index) % obj)
type is (ZernikeFilter)
order = f % order
type is (ZernikeRadialFilter)
order = f % order
class default
err = E_INVALID_TYPE
call set_errmsg("Not a Zernike filter.")
@ -144,7 +235,7 @@ contains
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (ZernikeFilter)
class is (ZernikeFilter)
x = f % x
y = f % y
r = f % r
@ -168,6 +259,9 @@ contains
type is (ZernikeFilter)
f % order = order
f % n_bins = ((order + 1)*(order + 2))/2
type is (ZernikeRadialFilter)
f % order = order
f % n_bins = order/2 + 1
class default
err = E_INVALID_TYPE
call set_errmsg("Not a Zernike filter.")
@ -187,7 +281,7 @@ contains
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (ZernikeFilter)
class is (ZernikeFilter)
if (present(x)) f % x = x
if (present(y)) f % y = y
if (present(r)) f % r = r

3
src/tallies/tally_header.F90
View file

@ -340,6 +340,9 @@ contains
type is (ZernikeFilter)
j = FILTER_ZERNIKE
this % estimator = ESTIMATOR_COLLISION
type is (ZernikeRadialFilter)
j = FILTER_ZERNIKE_RADIAL
this % estimator = ESTIMATOR_COLLISION
type is (ParticleFilter)
j = FILTER_PARTICLE
end select

19
src/xml_interface.F90
View file

@ -261,6 +261,7 @@ contains
! Get value of text node/attribute
str = node_value_string(node, name)
call whitespace_to_blanks(str)
! Read numbers into array
read(UNIT=str, FMT=*, IOSTAT=stat) array
@ -283,6 +284,7 @@ contains
! Get value of text node/attribute
str = node_value_string(node, name)
call whitespace_to_blanks(str)
! Read numbers into array
read(UNIT=str, FMT=*, IOSTAT=stat) array
@ -333,4 +335,21 @@ contains
end if
end subroutine get_node_array_string
!===============================================================================
! WHITESPACE_TO_BLANKS converts all whitespace to blanks
!===============================================================================
subroutine whitespace_to_blanks(str)
character(len=*, kind=C_CHAR), intent(inout) :: str
integer :: i
do i = 1, len(str)
select case (str(i:i))
case (C_NEW_LINE, C_HORIZONTAL_TAB, C_CARRIAGE_RETURN)
str(i:i) = ' '
end select
end do
end subroutine
end module xml_interface

17
tests/unit_tests/test_filters.py
View file

@ -100,6 +100,23 @@ def test_zernike():
assert elem.attrib['type'] == 'zernike'
assert elem.find('order').text == str(n)
def test_zernike_radial():
n = 4
f = openmc.ZernikeRadialFilter(n, 0., 0., 1.)
assert f.order == n
assert f.bins[0] == 'Z0,0'
assert f.bins[-1] == 'Z{},0'.format(n)
assert len(f.bins) == n//2 + 1
# Make sure __repr__ works
repr(f)
# to_xml_element()
elem = f.to_xml_element()
assert elem.tag == 'filter'
assert elem.attrib['type'] == 'zernikeradial'
assert elem.find('order').text == str(n)
def test_first_moment(run_in_tmpdir, box_model):
plain_tally = openmc.Tally()