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

Merge 6125408721 into 61c8a59cff

Browse source
This commit is contained in:
pferney05 2026-07-20 18:24:36 +03:00 • committed by GitHub
commit f54d67bb0f
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
63 changed files with 5933 additions and 7 deletions
Download patch file Download diff file Expand all files Collapse all files

2
CMakeLists.txt
View file

@ -380,6 +380,8 @@ list(APPEND libopenmc_SOURCES
src/geometry_aux.cpp
src/hdf5_interface.cpp
src/ifp.cpp
src/implicit.cpp
src/implicit_solvers.cpp
src/initialize.cpp
src/lattice.cpp
src/material.cpp

35
README.md
View file

@ -1,3 +1,38 @@
# For Implicit Function development
This branch aims to develop a framework to create any surface defined by an implicit function. A dedicated solver and architecture is proposed.
Developer cross sections : https://anl.box.com/shared/static/teaup95cqv8s9nn56hfn7ku8mmelr95p.xz
## Quick install
```sh
conda create -n openmc-IF compilers=1.9.0 cmake hdf5 python libpng
git clone --recurse-submodules https://github.com/pferney05/openmc-dev.git
cd openmc
git checkout ImplicitFunction
mkdir build
cd build
cmake -DOPENMC_ENABLE_STRICT_FP=on -DOPENMC_BUILD_TESTS=ON -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX ..
make -j 12
make install
cd ..
python -m pip install .[test]
```
To get git-clang-format:
```sh
conda install clang-format=18.1.8 wget
curl -L \
https://raw.githubusercontent.com/llvm/llvm-project/llvmorg-18.1.8/clang/tools/clang-format/git-clang-format \
-o "$CONDA_PREFIX/bin/git-clang-format"
chmod +x "$CONDA_PREFIX/bin/git-clang-format"
./tools/dev/install-commit-hooks.sh
```
## Local Notes:
- TPMS: max event particles, to be checked
- TPMS: Particle could not be located after crossing a boundary of lattice
- SolidRayTracing, aliasing effect, to be checked
# OpenMC Monte Carlo Particle Transport Code
[![License](https://img.shields.io/badge/license-MIT-green)](https://docs.openmc.org/en/latest/license.html)

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

@ -84,6 +84,8 @@ Building geometry
openmc.XTorus
openmc.YTorus
openmc.ZTorus
openmc.ImplicitSurface
openmc.TPMS
openmc.Halfspace
openmc.Intersection
openmc.Union
@ -107,6 +109,15 @@ Many of the above classes are derived from several abstract classes:
openmc.Region
openmc.Lattice
Implicit surfaces are defined with a
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.ImplicitFunction
.. _pythonapi_tallies:
Constructing Tallies

32
docs/source/pythonapi/implicit.rst Normal file
View file

@ -0,0 +1,32 @@
.. _pythonapi_implicit:
--------------------------------------------
:mod:`openmc.implicit` -- Implicit Functions
--------------------------------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.implicit.ImplicitFunction
openmc.implicit.X
openmc.implicit.Y
openmc.implicit.Z
openmc.implicit.Constant
openmc.implicit.Add
openmc.implicit.Sub
openmc.implicit.Neg
openmc.implicit.Scale
openmc.implicit.Mul
openmc.implicit.Div
openmc.implicit.Pow
openmc.implicit.Sin
openmc.implicit.Cos
openmc.implicit.Sqrt
openmc.implicit.Exp
openmc.implicit.Log
openmc.implicit.Abs
openmc.implicit.Min
openmc.implicit.Max
openmc.implicit.Cached

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

@ -45,6 +45,7 @@ or class.
model
examples
deplete
implicit
mgxs
stats
data

314
docs/source/usersguide/geometry.rst
View file

@ -157,6 +157,320 @@ work::
(array([-0.8660254, -inf, -inf]),
array([ 0.8660254, inf, inf]))
Implicit Surface with custom equations
--------------------------------------
In addition to the algebraic surfaces listed above, OpenMC supports
**implicit surfaces** defined by an arbitrary smooth function
:math:`f(x, y, z) = c`. This makes it possible to model geometries that
cannot be expressed as polynomials - most notably Triply Periodic Minimal
Surfaces (`TPMS <https://www.sciencedirect.com/science/article/pii/S014919702300330X>`_).
The surface is the level set:
.. math::
f\!\left(\mathbf{R}(\mathbf{r} - \mathbf{r}_0)\right) = c
where :math:`\mathbf{r}_0` is a translation vector, :math:`\mathbf{R}` is a
rotation matrix, and :math:`c` is the isovalue. The negative half-space
(:math:`f < c`) is the "inside" and the positive half-space (:math:`f > c`)
is the "outside". This formulation separates the surface geometry from the
coordinate transform, simplifying translation and rotation without modifying
the function itself.
Defining a custom function
~~~~~~~~~~~~~~~~~~~~~~~~~~
Functions are built as expression trees using nodes from the
:mod:`openmc.implicit` module. Each node represents either a terminal value
or an operation, and they compose with natural Python syntax::
from openmc.implicit import X, Y, Z, Sin, Cos, Cached
import numpy as np
# Sphere: f(x,y,z) = x² + y² + z² = R²
sphere_func = X()**2 + Y()**2 + Z()**2
sphere = openmc.ImplicitSurface(function=sphere_func, isovalue=100.)
# Gyroid TPMS with pitch L
L = 1.0
cx = Cached(2 * np.pi * X() / L)
cy = Cached(2 * np.pi * Y() / L)
cz = Cached(2 * np.pi * Z() / L)
gyroid_func = Sin(cx)*Cos(cz) + Sin(cy)*Cos(cx) + Sin(cz)*Cos(cy)
gyroid = openmc.ImplicitSurface(function=gyroid_func, isovalue=0.)
Alternatively, common TPMS families can be constructed directly from a pitch
and isovalue using the :class:`openmc.TPMS` convenience class::
gyroid = openmc.TPMS.from_pitch_isovalue('gyroid', pitch=1.0, isovalue=0.)
prim = openmc.TPMS.from_pitch_isovalue('primitive', pitch=1.0, isovalue=0.)
diamond = openmc.TPMS.from_pitch_isovalue('diamond', pitch=1.0, isovalue=0.)
Available node types
~~~~~~~~~~~~~~~~~~~~
The following node types are available in :mod:`openmc.implicit`:
.. table:: Terminal nodes
+--------------------+---------------------------+
| Node | Value |
+====================+===========================+
| ``X()`` | :math:`x` |
+--------------------+---------------------------+
| ``Y()`` | :math:`y` |
+--------------------+---------------------------+
| ``Z()`` | :math:`z` |
+--------------------+---------------------------+
| ``Constant(v)`` | :math:`v` |
+--------------------+---------------------------+
.. table:: Arithmetic nodes
+--------------------+----------------------------------+
| Node | Value |
+====================+==================================+
| ``f + g`` | :math:`f + g` |
+--------------------+----------------------------------+
| ``f - g`` | :math:`f - g` |
+--------------------+----------------------------------+
| ``k * f`` | :math:`k \cdot f` (scalar) |
+--------------------+----------------------------------+
| ``f * g`` | :math:`f \cdot g` |
+--------------------+----------------------------------+
| ``f / g`` | :math:`f / g` :math:`g \ne 0` |
+--------------------+----------------------------------+
| ``f ** n`` | :math:`f^n` (positive int only) |
+--------------------+----------------------------------+
| ``Min(f, g)`` | :math:`\min(f, g)` |
+--------------------+----------------------------------+
| ``Max(f, g)`` | :math:`\max(f, g)` |
+--------------------+----------------------------------+
.. table:: Transcendental nodes
+--------------------+-------------------------------+
| Node | Value |
+====================+===============================+
| ``Sin(f)`` | :math:`\sin(f)` |
+--------------------+-------------------------------+
| ``Cos(f)`` | :math:`\cos(f)` |
+--------------------+-------------------------------+
| ``Sqrt(f)`` | :math:`\sqrt{f}`, :math:`f>0` |
+--------------------+-------------------------------+
| ``Exp(f)`` | :math:`e^f` |
+--------------------+-------------------------------+
| ``Log(f)`` | :math:`\ln f`, :math:`f>0` |
+--------------------+-------------------------------+
| ``Abs(f)`` | :math:`|f|` |
+--------------------+-------------------------------+
.. note::
``Pow`` only accepts strictly positive integer exponents. Use
``Div(Constant(1.), f)`` for :math:`f^{-1}` and ``Sqrt(f)`` for
:math:`f^{0.5}`.
Sub-expression caching
~~~~~~~~~~~~~~~~~~~~~~
When the same sub-expression appears multiple times in a function - as
``cx``, ``cy``, ``cz`` do in the gyroid example above - wrapping it in
:class:`~openmc.implicit.Cached` tells the C++ solver to compute it once and
reuse the result within each geometry step. This can significantly reduce
the number of function evaluations for complex TPMS expressions.
The **same Python object** must be reused at every occurrence for sharing to
take effect. Constructing a new ``Cached(...)`` at each use site produces
independent caches with no benefit::
# Correct - one object, two references
cx = Cached(2 * np.pi * X())
f = Sin(cx) + Cos(cx)
# Wrong - two independent caches, no sharing
f = Sin(Cached(2 * np.pi * X())) + Cos(Cached(2 * np.pi * X()))
Finite region requirement
~~~~~~~~~~~~~~~~~~~~~~~~~
The implicit surface solver needs an upper bound on how far to search along
a ray. This bound is provided automatically by the nearest analytical surface
in the same cell region. **Every cell containing an implicit surface must
therefore also reference at least one finite analytical surface** (plane,
sphere, cylinder, etc.) to enclose it::
R = 10.0
outer = openmc.Sphere(r=R * 1.1, boundary_type='vacuum')
impl = openmc.ImplicitSurface(function=gyroid_func, isovalue=0.)
fuel_cell = openmc.Cell(region=-impl & -outer, fill=fuel)
void_cell = openmc.Cell(region=+impl & -outer)
A fatal error is raised at runtime if an implicit surface is placed in a
region with no finite analytical boundary.
Bounding region and Lipschitz computation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The size of the analytical bounding region matters beyond simply excluding
geometrically invalid points. Even when the implicit surface is well-defined
everywhere inside the cell, the solver may fail because it computes Lipschitz
bounds **over the entire ray interval** - from the particle's current position
to the nearest analytical boundary. If that interval passes through a region
where a node's domain condition is violated, the Lipschitz computation throws
a domain error before the solver begins marching.
Consider the following example. The function involves ``Sqrt(x² + y² + z²)``,
which requires a strictly positive argument. The domain condition is satisfied
at every point in the shell :math:`1 \leq r \leq 5`, so the geometry appears
valid::
# Bounding shell
inner = openmc.Sphere(r=1.0)
outer = openmc.Sphere(r=5.0)
shell = +inner & -outer
# Vacuum boundary
world = openmc.Sphere(r=20.0, boundary_type='vacuum')
# Implicit surface
r = Sqrt(X()**2 + Y()**2 + Z()**2)
surf = ImplicitSurface(function=r, isovalue=5.)
# Cells
fuel_cell = openmc.Cell(region=-surf & shell, fill=material)
void_cell = openmc.Cell(region=+surf & shell)
outer_cell = openmc.Cell(region= ~shell & -world)
geometry = openmc.Geometry([fuel_cell, void_cell, outer_cell])
This crashes with::
std::domain_error: Sqrt::compute_lipschitz: argument reaches zero or
negative on ray between r=(1.057, 0.571, -1.830) and r=(-4.235, -2.657,
0.047) in expression X**2 + Y**2 + Z**2
The ray interval spans from the particle's current position to the nearest
bounding surface. Even though both endpoints lie within the shell, the ray is
a straight line whose closest approach to the origin can be much smaller than
either endpoint's radius - potentially passing through the region :math:`r = 0`
where ``Sqrt`` derivative is undefined.
The fix is to subdivide the geometry into cells whose bounding regions
guarantee that the function is well-defined over every possible ray interval
within each cell.
Ray-surface intersection solvers
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Finding where a neutron ray crosses an implicit surface is not trivial - unlike
algebraic surfaces, there is no closed-form intersection formula. OpenMC uses
a **Lipschitz-marching** algorithm that guarantees no root is missed.
The key insight is the Lipschitz condition: for a smooth function :math:`f`,
the rate of change along a ray is bounded by a constant :math:`L` (the
Lipschitz constant):
.. math::
|f(t + \delta) - f(t)| \leq L \cdot \delta
This means that if the function currently has value :math:`f(t)`, it cannot
reach zero before advancing at least :math:`|f(t)| / L` along the ray. The
solver exploits this to take large, guaranteed-safe steps, only slowing down
as it approaches a root.
Two solvers are available:
``naive``
Computes :math:`L` once over the full ray interval and marches with steps
of size :math:`|f| / L`. When a sign change is detected, bisection
refines the root location. Simple and robust, but can be slow when
:math:`L` is large relative to the actual function variation (e.g. for
highly oscillatory TPMS over long rays).
``fast``
Subdivides the ray interval into a stack of sub-intervals and recomputes
tight Lipschitz bounds on each one. Sub-intervals where :math:`f` has
constant sign are discarded immediately; only intervals containing a root
are refined further. Significantly faster than ``naive`` for TPMS
geometries where the global :math:`L` is large but local bounds are tight.
The solver is selected via :attr:`openmc.Settings.implicit`::
settings.implicit = {'name': 'fast'} # default
Tuning solver parameters
~~~~~~~~~~~~~~~~~~~~~~~~
Several parameters control the accuracy and robustness of the solver and can
be set through :attr:`openmc.Settings.implicit`.
``atol`` - geometric tolerance
The solver stops as soon as the estimated distance to the root is smaller
than ``atol`` [cm]. Decreasing ``atol`` gives a more accurate root
location at the cost of more function evaluations. The default value of
``1e-9`` cm is appropriate for most
applications::
settings.implicit = {'atol': 1e-9}
``ftol`` - function value tolerance
A root is also accepted when the function value :math:`|f - c|` falls
below ``ftol``. This can be useful when the surface passes through a
nearly-flat region where the gradient is small and the geometric distance
to the root is hard to estimate from :math:`L` alone. Defaults to
``1e-9``::
settings.implicit = {'ftol': 1e-9}
``maxiter`` - maximum solver iterations
Safety cap on the number of marching steps. If the solver reaches
``maxiter`` without finding a root, it skips the surface and emits a
warning. This prevents an infinite loop when the function has very
small gradients over a long ray segment. The default of ``1 000 000``
is rarely reached in practice but can be lowered for faster failure
detection during debugging::
settings.implicit = {'maxiter': 10000}
If you see the warning
.. code-block:: none
WARNING: NaiveLipschitz reached max iterations.
it might mean that the solver steps are very small relative to the length
of the bounding interval because of a large Lipschitz constant. Switching
to ``'fast'`` (which recomputes tighter per-interval bounds) may eliminate
the issue.
``margin`` - coincident surface margin
After the solver finds a root, the intersection point is nudged forward
by ``margin`` [cm] to ensure the particle lands unambiguously on the
destination side of the surface. Without this nudge, floating-point
rounding can place the particle within the tolerance band of the surface,
causing OpenMC's ``sense()`` function to invoke the surface normal for a
tiebreak - which may fail for functions with restricted gradient domains.
The default of ``1e-6`` cm is sufficient for all tested geometries::
settings.implicit = {'margin': 1e-6}
Increasing ``margin`` slightly can resolve rare geometry errors where a
particle appears to bounce back across a surface it just crossed.
All parameters can be combined in a single assignment::
settings.implicit = {
'name': 'fast',
'atol': 1e-9,
'ftol': 1e-9,
'maxiter': 1000000,
'margin': 1e-6,
}
Boundary Conditions
-------------------

2
include/openmc/constants.h
View file

@ -373,7 +373,7 @@ enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT };
//==============================================================================
// Geometry Constants
enum class GeometryType { CSG, DAG };
enum class GeometryType { CSG, DAG, IMP };
// a surface token cannot be zero due to the unsigned nature of zero for integer
// representations. This value represents no surface.

550
include/openmc/implicit.h Normal file
View file

@ -0,0 +1,550 @@
#ifndef OPENMC_IMPLICIT_H
#define OPENMC_IMPLICIT_H
#include <atomic>
#include <cmath>
#include <cstdint>
#include <limits> // For numeric_limits
#include <memory>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
#include "pugixml.hpp"
#include "openmc/position.h"
namespace openmc {
//==============================================================================
// Type aliases
//==============================================================================
// Gradient shares the memory layout of Position (x, y, z doubles).
// Using a named alias makes intent explicit at every call site.
using Gradient = Position;
//==============================================================================
// StepCache
//
// Thread-local epoch counter. Increment step_cache.epoch once at the
// start of each geometry step. All Cached nodes across the DAG are
// lazily invalidated at zero cost — no traversal needed.
//==============================================================================
struct CacheEntry {
uint64_t epoch {std::numeric_limits<uint64_t>::max()};
Position r {0.0, 0.0, 0.0};
double eval {0.0};
Gradient grad {0.0, 0.0, 0.0};
bool grad_valid {false};
// Interval cache — keyed on (u, t0, t1); r is covered by r + epoch.
Direction u {0.0, 0.0, 1.0};
double t0 {0.0}, t1 {0.0};
double L {1.0};
bool L_valid {false};
std::pair<double, double> min_max {0.0, 0.0};
bool min_max_valid {false};
};
struct StepCache {
uint64_t epoch {0};
// Per-thread cache indexed by Cached::slot_.
// Grows to exactly one entry per Cached node in the model — bounded and tiny.
std::vector<CacheEntry> node_cache;
};
extern thread_local StepCache step_cache;
//==============================================================================
// Implicit
//
// Abstract base class for expression DAG nodes representing a smooth
// 3D function f(x, y, z). Concrete subclasses are defined in the
// openmc::implicit namespace below.
//
// The DAG is constructed in Python (implicit_function.py) and serialised
// to XML. On the C++ side it is reconstructed via from_xml_element() and
// evaluated during implicit surface geometry queries.
//==============================================================================
class Implicit {
public:
virtual ~Implicit() = default;
virtual std::string expression() const = 0;
//! Evaluate f(r).
virtual double evaluate(Position r) const = 0;
//! Gradient df/dx, df/dy, df/dz.
//! Computed analytically — no finite differences.
virtual Gradient gradient(Position r) const = 0;
//! Evaluate f along the ray: f(r + t*u).
double along_ray(double t, Position r, Direction u) const;
//! Lipschitz constant of f on the ray segment [t0, t1]:
//! |f(r + a*u) - f(r + b*u)| <= L * |a - b| for all a,b in [t0,t1].
//! Computed analytically via interval arithmetic propagation.
virtual double compute_lipschitz(
Position r, Direction u, double t0, double t1) const = 0;
//! Tight lower and upper bounds of f on the ray segment [t0, t1].
//! Returned as {f_min, f_max}.
//! Computed jointly because many nodes (Mul, Div, Pow) need both bounds
//! together to propagate the interval correctly.
virtual std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const = 0;
// create a xml element out of the node
virtual pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const = 0;
// Convenience wrapper — serialises the whole DAG to an XML string
std::string to_xml_string() const;
//! Reconstruct an Implicit DAG from an XML element produced by
//! Python's to_xml_element().
//!
//! \param node The XML element to parse.
//! \param cache_map Maps numeric ids to already-constructed shared nodes.
//! Populated when a <to_cache> element is encountered;
//! looked up when a <from_cache> element is encountered.
//! Passed by reference through the recursion.
static std::shared_ptr<Implicit> from_xml_element(pugi::xml_node node,
std::unordered_map<int, std::shared_ptr<Implicit>>& cache_map);
//! Convenience overload for the root call — creates an empty cache_map.
static std::shared_ptr<Implicit> from_xml_element(pugi::xml_node node);
};
//==============================================================================
// Concrete node types (openmc::implicit namespace)
//
// All implementations are in implicit.cpp.
// Nodes follow the expression tree structure of the Python counterparts in
// implicit_function.py. The openmc::implicit sub-namespace avoids name
// collisions with common identifiers (X, Y, Z, Log, Exp, etc.).
//==============================================================================
namespace implicit {
// ============================================================================
// Terminal nodes
// ============================================================================
class X final : public Implicit {
public:
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
};
class Y final : public Implicit {
public:
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
};
class Z final : public Implicit {
public:
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
};
class Constant final : public Implicit {
public:
explicit Constant(double value) : value_(value) {}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
double value_;
};
// ============================================================================
// Arithmetic nodes
// ============================================================================
class Add final : public Implicit {
public:
explicit Add(std::shared_ptr<Implicit> f, std::shared_ptr<Implicit> g)
: f_(std::move(f)), g_(std::move(g))
{}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_, g_;
};
class Sub final : public Implicit {
public:
explicit Sub(std::shared_ptr<Implicit> f, std::shared_ptr<Implicit> g)
: f_(std::move(f)), g_(std::move(g))
{}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_, g_;
};
class Mul final : public Implicit {
public:
explicit Mul(std::shared_ptr<Implicit> f, std::shared_ptr<Implicit> g)
: f_(std::move(f)), g_(std::move(g))
{}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_, g_;
};
class Div final : public Implicit {
public:
explicit Div(std::shared_ptr<Implicit> f, std::shared_ptr<Implicit> g)
: f_(std::move(f)), g_(std::move(g))
{}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_, g_;
};
class Scale final : public Implicit {
public:
explicit Scale(std::shared_ptr<Implicit> f, double k)
: f_(std::move(f)), k_(k)
{}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_;
double k_;
};
class Neg final : public Implicit {
public:
explicit Neg(std::shared_ptr<Implicit> f) : f_(std::move(f)) {}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_;
};
class Pow final : public Implicit {
public:
explicit Pow(std::shared_ptr<Implicit> f, int exp)
: f_(std::move(f)), exp_(exp)
{
if (exp_ <= 0)
throw std::domain_error(
"Pow: exponent must be a strictly positive integer, got " +
std::to_string(exp_));
}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_;
int exp_;
};
// ============================================================================
// Transcendental nodes
// ============================================================================
class Sin final : public Implicit {
public:
explicit Sin(std::shared_ptr<Implicit> arg) : arg_(std::move(arg)) {}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> arg_;
};
class Cos final : public Implicit {
public:
explicit Cos(std::shared_ptr<Implicit> arg) : arg_(std::move(arg)) {}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> arg_;
};
class Sqrt final : public Implicit {
public:
explicit Sqrt(std::shared_ptr<Implicit> arg) : arg_(std::move(arg)) {}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> arg_;
};
class Exp final : public Implicit {
public:
explicit Exp(std::shared_ptr<Implicit> arg) : arg_(std::move(arg)) {}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> arg_;
};
class Log final : public Implicit {
public:
explicit Log(std::shared_ptr<Implicit> arg) : arg_(std::move(arg)) {}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> arg_;
};
class Abs final : public Implicit {
public:
explicit Abs(std::shared_ptr<Implicit> arg) : arg_(std::move(arg)) {}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> arg_;
};
class Min final : public Implicit {
public:
explicit Min(std::shared_ptr<Implicit> f, std::shared_ptr<Implicit> g)
: f_(std::move(f)), g_(std::move(g))
{}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_, g_;
};
class Max final : public Implicit {
public:
explicit Max(std::shared_ptr<Implicit> f, std::shared_ptr<Implicit> g)
: f_(std::move(f)), g_(std::move(g))
{}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
std::shared_ptr<Implicit> f_, g_;
};
// ============================================================================
// Cached node
//
// Memoises evaluate() and gradient() for the duration of one geometry step
// AND one position. A cache hit requires BOTH a matching epoch AND a
// matching position.
//
// Epoch-only invalidation is INCORRECT: the ray-tracing solver evaluates f
// at many different positions within a single step, so the position check
// is essential for correctness.
//
// compute_lipschitz and compute_f_min_max are cached per ray query,
// keyed on (epoch, r, u, t0, t1). Within one solver call the three
// branches of a shared subtree hit the cache on the 2nd and 3rd traversal.
//
// Thread safety: Cached has no mutable members. All cache data lives in
// thread_local step_cache.node_cache, indexed by slot. Each thread maintains
// a fully independent cache — no shared mutable state, no data race.
//
// The map grows to at most one entry per Cached node in the model and is
// never cleared — entries are simply overwritten when the epoch or position
// changes.
//
// WARNING: the SAME shared_ptr instance must be used wherever this
// subexpression appears in the DAG. Constructing two separate Cached nodes
// wrapping the same expression gives two independent caches and no XML
// sharing. See the Python-side docstring in implicit_function.py for the
// full pitfall and correct usage pattern.
// ============================================================================
class Cached final : public Implicit {
public:
explicit Cached(std::shared_ptr<Implicit> child)
: child_(std::move(child)), slot_(next_slot_++)
{}
std::string expression() const override;
double evaluate(Position r) const override;
Gradient gradient(Position r) const override;
double compute_lipschitz(
Position r, Direction u, double t0, double t1) const override;
std::pair<double, double> compute_f_min_max(
Position r, Direction u, double t0, double t1) const override;
pugi::xml_node to_xml_element(pugi::xml_node parent,
std::unordered_map<const Implicit*, int>& cache_map) const override;
private:
//! Return this node's cache entry.
//! Grows the thread-local cache vector on first access.
CacheEntry& get_cache_entry() const;
//! Recompute eval if the epoch or position changed.
CacheEntry& refresh(Position r) const;
//! Like refresh, but additionally invalidates the interval cache
//! (L, min_max) if the ray parameters (u, t0, t1) changed.
CacheEntry& refresh_interval(
Position r, Direction u, double t0, double t1) const;
//! Global counter assigning a unique slot to each Cached node.
//! Only incremented during model construction (single-threaded).
static std::atomic<int> next_slot_;
std::shared_ptr<Implicit> child_;
int slot_;
};
} // namespace implicit
} // namespace openmc
#endif // OPENMC_IMPLICIT_H

84
include/openmc/implicit_solvers.h Normal file
View file

@ -0,0 +1,84 @@
#ifndef OPENMC_IMPLICIT_SOLVER_H
#define OPENMC_IMPLICIT_SOLVER_H
#include "openmc/implicit.h"
#include "openmc/position.h"
namespace openmc {
//==============================================================================
// ImplicitSolver
//
// Abstract base class for ray-implicit surface intersection solvers.
//
// Contract:
// - func has already been shifted by isovalue (SurfaceImplicit::evaluate
// returns f(r) - c, so the solver always looks for f = 0).
// - [t0, t1] is already clipped to the bounding box by the caller.
// - Returns the smallest positive root in [t0, t1], or INFTY if none found.
//==============================================================================
class ImplicitSolver {
public:
ImplicitSolver(double atol = 1e-9, double ftol = 1e-9, int max_iter = 1000000)
: atol_(atol), ftol_(ftol), max_iter_(max_iter)
{}
virtual ~ImplicitSolver() = default;
//! Find the smallest root of func(r + t*u) = 0 in [t0, t1].
virtual double solve(const Implicit& func, Position r, Direction u, double t0,
double t1, double isovalue = 0.0) const = 0;
// access atol outside the solver
double get_atol() const { return atol_; }
static std::unique_ptr<ImplicitSolver> create(
const std::string& name, int max_iter, double atol, double ftol);
protected:
double atol_;
double ftol_;
int max_iter_;
};
//==============================================================================
// NaiveLipschitz
//
// Lipschitz-marching solver. At each step, the Lipschitz constant L gives a
// guaranteed safe skip distance: if |f(t)| / L > remaining interval, no root
// is possible and we advance. When a sign change is detected, bisection
// finds the root to within atol.
//
// The Lipschitz constant is computed once over [t0, t1] at the start and
// reused throughout the march. This is the "naive" part — a more
// sophisticated solver would recompute L on each sub-interval.
//
// Parameters:
// atol Geometric tolerance for root detection (metres). A value of
// 1e-8 (roughly one atomic radius) is a sensible default.
// max_iter Safety cap on the number of marching steps. If reached, the
// solver returns INFTY and the surface is treated as not hit.
//==============================================================================
class NaiveLipschitz : public ImplicitSolver {
public:
using ImplicitSolver::ImplicitSolver;
double solve(const Implicit& function, Position r, Direction u, double t0,
double t1, double isovalue) const override;
private:
//! Bisect on [ta, tb] given that f(ta) and f(tb) have opposite signs.
//! Returns the root to within atol_.
double bisect(const Implicit& func, Position r, Direction u, double ta,
double tb, double fa, double fb, double isovalue) const;
};
class FastLipschitz : public ImplicitSolver {
public:
using ImplicitSolver::ImplicitSolver;
double solve(const Implicit& function, Position r, Direction u, double t0,
double t1, double isovalue) const override;
};
} // namespace openmc
#endif // OPENMC_IMPLICIT_SOLVER_H

7
include/openmc/settings.h
View file

@ -204,6 +204,13 @@ extern "C" int verbosity; //!< How verbose to make output
extern double weight_cutoff; //!< Weight cutoff for Russian roulette
extern double weight_survive; //!< Survival weight after Russian roulette
// implicit solvers
extern int implicit_maxiter;
extern std::string implicit_solver;
extern double implicit_atol;
extern double implicit_ftol;
extern double implicit_margin;
} // namespace settings
//==============================================================================

29
include/openmc/surface.h
View file

@ -12,6 +12,8 @@
#include "openmc/boundary_condition.h"
#include "openmc/bounding_box.h"
#include "openmc/constants.h"
#include "openmc/implicit.h"
#include "openmc/implicit_solvers.h"
#include "openmc/memory.h" // for unique_ptr
#include "openmc/particle.h"
#include "openmc/position.h"
@ -375,6 +377,33 @@ public:
double x0_, y0_, z0_, A_, B_, C_;
};
//==============================================================================
//! An implicit surface described by a user specified function
//==============================================================================
class SurfaceImplicit : public Surface {
public:
explicit SurfaceImplicit(pugi::xml_node surf_node);
double evaluate(Position r) const override;
double distance_finite(
Position r, Direction u, bool coincident, double distance_max) const;
double distance(Position r, Direction u, bool coincident) const override;
Direction normal(Position r) const override;
void to_hdf5_inner(hid_t group_id) const override;
GeometryType geom_type() const override { return GeometryType::IMP; }
private:
//! Apply the surface transform: r_local = R * (r - r0)
Position transform(Position r) const;
//! Apply rotation only (no translation) for directions: u_local = R * u
Direction transform_dir(Direction u) const;
double x0_, y0_, z0_, A_, B_, C_, D_, E_, F_, G_, H_, I_;
double isovalue_;
std::shared_ptr<Implicit> function_;
std::unique_ptr<ImplicitSolver> solver_;
};
//==============================================================================
// Non-member functions
//==============================================================================

15
openmc/geometry.py
View file

@ -701,9 +701,18 @@ class Geometry:
key = (surf._type, surf._boundary_type) + coeffs
redundancies[key].append(surf)
redundant_surfaces = {replace.id: keep
for keep, *redundant in redundancies.values()
for replace in redundant}
redundant_surfaces = {}
for group in redundancies.values():
kept = [group[0]]
for surf in group[1:]:
equivalent = next(
(keep for keep in kept if surf.is_equal(keep)),
None
)
if equivalent is not None:
redundant_surfaces[surf.id] = equivalent
else:
kept.append(surf)
if redundant_surfaces:
# Iterate through all cells contained in the geometry

544
openmc/implicit.py Normal file
View file

@ -0,0 +1,544 @@
from __future__ import annotations
from abc import ABC, abstractmethod
import lxml.etree as ET
import numpy as np
# ------------------------------------------------------------------
# Helper functions
# ------------------------------------------------------------------
def _to_function(f: int | float | ImplicitFunction) -> ImplicitFunction:
return Constant(float(f)) if isinstance(f, (int, float)) else f
# ------------------------------------------------------------------
# Main ImplicitFunction Abstract class
# ------------------------------------------------------------------
class ImplicitFunction(ABC):
"""Abstract base class for nodes in an implicit function expression tree.
Subclasses represent either terminal values (coordinates, constants) or
operations that combine child nodes. The tree is evaluated by calling
:meth:`evaluate` and serialised to XML via :meth:`to_xml_element`.
Operator overloading allows trees to be built with natural Python syntax::
f = Sin(2 * X()) * Cos(Z()) + Constant(1.)
"""
@abstractmethod
def __repr__(self) -> str: ...
@abstractmethod
def evaluate(self, point) -> float: ...
@abstractmethod
def to_xml_element(self, _cached: list[int] | None = None) -> ET.Element: ...
@classmethod
def from_xml_element(cls, element: ET.Element, _cached: dict | None = None) -> ImplicitFunction:
if _cached is None: _cached = {}
tag = element.tag
attrib = element.attrib
# Handle cache reference before parsing children (no children to parse)
if tag == "from_cache":
return _cached[int(attrib["id"])]
# Recursively parse children for all other tags
children = [cls.from_xml_element(child, _cached) for child in element]
# Register a new cached node and return it
if tag == "to_cache":
node = Cached(children[0])
_cached[int(attrib["id"])] = node
return node
dispatch = {
"x": lambda: X(),
"y": lambda: Y(),
"z": lambda: Z(),
"constant": lambda: Constant(float(attrib["value"])),
"add": lambda: Add(*children),
"neg": lambda: Neg(*children),
"sub": lambda: Sub(*children),
"scale": lambda: Scale(children[0], float(attrib["value"])),
"mul": lambda: Mul(*children),
"div": lambda: Div(*children),
"pow": lambda: Pow(children[0], int(attrib["value"])),
"sin": lambda: Sin(*children),
"cos": lambda: Cos(*children),
"sqrt": lambda: Sqrt(*children),
"exp": lambda: Exp(*children),
"log": lambda: Log(*children),
"abs": lambda: Abs(*children),
"max": lambda: Max(*children),
"min": lambda: Min(*children),
}
if tag not in dispatch:
raise ValueError(f"Unknown tag '{tag}'")
return dispatch[tag]()
# ------------------------------------------------------------------
# Operator overloading - enables natural Python expression syntax
# e.g. Sin(X()) * Cos(Z()) + Constant(1)
# ------------------------------------------------------------------
def __add__(self, other: float | ImplicitFunction) -> ImplicitFunction:
return Add(_to_function(self), _to_function(other))
def __radd__(self, other: float | ImplicitFunction) -> ImplicitFunction:
return Add(_to_function(other), _to_function(self))
def __neg__(self) -> ImplicitFunction:
return Neg(_to_function(self))
def __sub__(self, other: float | ImplicitFunction) -> ImplicitFunction:
return Sub(_to_function(self), _to_function(other))
def __rsub__(self, other: float | ImplicitFunction) -> ImplicitFunction:
return Sub(_to_function(other), _to_function(self))
def __truediv__(self, other: float | ImplicitFunction) -> ImplicitFunction:
return Div(_to_function(self), _to_function(other))
def __rtruediv__(self, other: float | ImplicitFunction) -> ImplicitFunction:
return Div(_to_function(other), _to_function(self))
def __pow__(self, exp: int) -> ImplicitFunction:
if not isinstance(exp, (int)):
raise TypeError(f"Pow exponent must be an int, got {type(exp)}")
if exp <= 0.:
raise TypeError(f"Pow exponent must be strictly positive, got {exp}")
return Pow(_to_function(self), exp)
def __mul__(self, other: float | ImplicitFunction) -> ImplicitFunction:
if isinstance(other, (int, float)):
return Scale(self, float(other))
return Mul(self, other)
def __rmul__(self, other: float | ImplicitFunction) -> ImplicitFunction:
if isinstance(other, (int, float)):
return Scale(self, float(other))
# ---------------------------------------------------------------------------
# Terminal nodes
# ---------------------------------------------------------------------------
class X(ImplicitFunction):
"""The x-coordinate: f(x, y, z) = x."""
def __repr__(self): return "X"
def evaluate(self, point): return point[0]
def to_xml_element(self, _cached=None): return ET.Element("x")
class Y(ImplicitFunction):
"""The y-coordinate: f(x, y, z) = y."""
def __repr__(self): return "Y"
def evaluate(self, point): return point[1]
def to_xml_element(self, _cached=None): return ET.Element("y")
class Z(ImplicitFunction):
"""The z-coordinate: f(x, y, z) = z."""
def __repr__(self): return "Z"
def evaluate(self, point): return point[2]
def to_xml_element(self, _cached=None): return ET.Element("z")
class Constant(ImplicitFunction):
"""A scalar constant: f(x, y, z) = value.
Parameters
----------
value : float
The constant value.
"""
def __repr__(self): return f"{self.value}"
def __init__(self, value: float) -> None:
self.value = float(value)
def evaluate(self, point): return self.value
def to_xml_element(self, _cached=None):
element = ET.Element("constant")
element.set("value", str(self.value))
return element
# ---------------------------------------------------------------------------
# Operations
# ---------------------------------------------------------------------------
class Add(ImplicitFunction):
"""Element-wise sum: f(x, y, z) = f(x,y,z) + g(x,y,z).
Parameters
----------
f, g : ImplicitFunction
Operands.
"""
def __repr__(self): return f"{self.f} + {self.g}"
def __init__(self, f:ImplicitFunction, g:ImplicitFunction):
self.f = f
self.g = g
def evaluate(self, point): return self.f.evaluate(point) + self.g.evaluate(point)
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("add")
element.append(self.f.to_xml_element(_cached))
element.append(self.g.to_xml_element(_cached))
return element
class Neg(ImplicitFunction):
"""Negation: f(x, y, z) = -f(x, y, z).
Parameters
----------
f : ImplicitFunction
Operand.
"""
def __repr__(self): return f"-{self.f}"
def __init__(self, f:ImplicitFunction):
self.f = f
def evaluate(self, point): return -self.f.evaluate(point)
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("neg")
element.append(self.f.to_xml_element(_cached))
return element
class Sub(ImplicitFunction):
"""Element-wise difference: h(x, y, z) = f(x,y,z) - g(x,y,z).
Parameters
----------
f, g : ImplicitFunction
Operands.
"""
def __repr__(self): return f"{self.f} - {self.g}"
def __init__(self, f:ImplicitFunction, g:ImplicitFunction):
self.f = f
self.g = g
def evaluate(self, point): return self.f.evaluate(point) - self.g.evaluate(point)
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("sub")
element.append(self.f.to_xml_element(_cached))
element.append(self.g.to_xml_element(_cached))
return element
class Scale(ImplicitFunction):
"""Multiplication by a scalar constant: h = scalar * f.
Prefer this over ``Mul(Constant(k), f)`` when one factor is a plain
number - the Lipschitz constant and interval bounds are tighter.
Parameters
----------
f : ImplicitFunction
Function to scale.
scalar : float
Scalar multiplier.
"""
def __repr__(self): return f"{self.scalar} * {self.f}"
def __init__(self, f:ImplicitFunction, scalar: float):
self.f = f
self.scalar = float(scalar)
def evaluate(self, point): return self.scalar * self.f.evaluate(point)
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("scale")
element.set("value", str(self.scalar))
element.append(self.f.to_xml_element(_cached))
return element
class Mul(ImplicitFunction):
"""Point-wise product of two functions: h = f * g.
Use :class:`Scale` instead when one factor is a plain number.
Parameters
----------
f, g : ImplicitFunction
Operands.
"""
def __repr__(self): return f"{self.f} * {self.g}"
def __init__(self, f:ImplicitFunction, g:ImplicitFunction):
self.f = f
self.g = g
def evaluate(self, point): return self.f.evaluate(point) * self.g.evaluate(point)
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("mul")
element.append(self.f.to_xml_element(_cached))
element.append(self.g.to_xml_element(_cached))
return element
class Div(ImplicitFunction):
"""Point-wise quotient: h = f / g.
Parameters
----------
f, g : ImplicitFunction
Numerator and denominator.
Raises
------
ValueError
If ``g`` evaluates to zero at the query point.
"""
def __repr__(self): return f"{self.f} / {self.g}"
def __init__(self, f:ImplicitFunction, g:ImplicitFunction):
self.f = f
self.g = g
def evaluate(self, point):
v = self.g.evaluate(point)
if v == 0.: raise ValueError(f"Denominator value cannot be zero.")
return self.f.evaluate(point) / v
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("div")
element.append(self.f.to_xml_element(_cached))
element.append(self.g.to_xml_element(_cached))
return element
class Pow(ImplicitFunction):
"""Integer power: h = f ** exp.
Parameters
----------
f : ImplicitFunction
Base.
exp : int
Exponent - must be a strictly positive integer.
Use ``Div(Constant(1.), f)`` for exp = -1 and
``Sqrt(f)`` for exp = 0.5.
Raises
------
TypeError
If ``exp`` is not a strictly positive integer.
"""
def __repr__(self): return f"{self.f} ** {self.exp}"
def __init__(self, f:ImplicitFunction, exp: int):
if not isinstance(exp, int) or exp <= 0:
raise TypeError(
f"Pow exponent must be a strictly positive integer, got {exp!r}. "
f"For exp=-1 use Div, for exp=0.5 use Sqrt.")
self.f = f
self.exp = exp
def evaluate(self, point): return self.f.evaluate(point) ** self.exp
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("pow")
element.set("value", str(self.exp))
element.append(self.f.to_xml_element(_cached))
return element
# ---------------------------------------------------------------------------
# Functions
# ---------------------------------------------------------------------------
class Sin(ImplicitFunction):
"""Sine: h(x, y, z) = sin(arg(x, y, z)).
Parameters
----------
arg : ImplicitFunction
Argument in radians.
"""
def __repr__(self): return f"Sin({self.arg})"
def __init__(self, arg:ImplicitFunction):
self.arg = arg
def evaluate(self, point): return np.sin(self.arg.evaluate(point))
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("sin")
element.append(self.arg.to_xml_element(_cached))
return element
class Cos(ImplicitFunction):
"""Cosine: h(x, y, z) = cos(arg(x, y, z)).
Parameters
----------
arg : ImplicitFunction
Argument in radians.
"""
def __repr__(self): return f"Cos({self.arg})"
def __init__(self, arg:ImplicitFunction):
self.arg = arg
def evaluate(self, point): return np.cos(self.arg.evaluate(point))
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("cos")
element.append(self.arg.to_xml_element(_cached))
return element
class Sqrt(ImplicitFunction):
"""Square root: h(x, y, z) = sqrt(arg(x, y, z)).
Parameters
----------
arg : ImplicitFunction
Argument - must be non-negative at every evaluation point.
Raises
------
ValueError
If ``arg`` evaluates to a negative value.
"""
def __repr__(self): return f"Sqrt({self.arg})"
def __init__(self, arg:ImplicitFunction):
self.arg = arg
def evaluate(self, point):
v = self.arg.evaluate(point)
if v < 0.: raise ValueError(f"Sqrt input value cannot be negative but is: '{v}'.")
return np.sqrt(v)
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("sqrt")
element.append(self.arg.to_xml_element(_cached))
return element
class Exp(ImplicitFunction):
"""Natural exponential: h(x, y, z) = exp(arg(x, y, z)).
Parameters
----------
arg : ImplicitFunction
Exponent.
"""
def __repr__(self): return f"Exp({self.arg})"
def __init__(self, arg:ImplicitFunction):
self.arg = arg
def evaluate(self, point): return np.exp(self.arg.evaluate(point))
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("exp")
element.append(self.arg.to_xml_element(_cached))
return element
class Log(ImplicitFunction):
"""Natural logarithm: h(x, y, z) = log(arg(x, y, z)).
Parameters
----------
arg : ImplicitFunction
Argument - must be strictly positive at every evaluation point.
Raises
------
ValueError
If ``arg`` evaluates to a non-positive value.
"""
def __repr__(self): return f"Log({self.arg})"
def __init__(self, arg:ImplicitFunction):
self.arg = arg
def evaluate(self, point):
v = self.arg.evaluate(point)
if v < 0.: raise ValueError(f"Log input value cannot be negative but is: '{v}'.")
return np.log(v)
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("log")
element.append(self.arg.to_xml_element(_cached))
return element
class Abs(ImplicitFunction):
"""Absolute value: h(x, y, z) = Abs(arg(x, y, z)).
Parameters
----------
arg : ImplicitFunction
Argument.
"""
def __repr__(self): return f"|{self.arg}|"
def __init__(self, arg:ImplicitFunction):
self.arg = arg
def evaluate(self, point): return np.abs(self.arg.evaluate(point))
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("abs")
element.append(self.arg.to_xml_element(_cached))
return element
class Min(ImplicitFunction):
"""Point-wise minimum: h(x,y,z) = min(f(x,y,z), g(x,y,z)).
Parameters
----------
f, g : ImplicitFunction
Operands.
"""
def __repr__(self): return f"Min({self.f}, {self.g})"
def __init__(self, f: ImplicitFunction, g: ImplicitFunction):
self.f = f
self.g = g
def evaluate(self, point):
return np.min((self.f.evaluate(point), self.g.evaluate(point)))
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("min")
element.append(self.f.to_xml_element(_cached))
element.append(self.g.to_xml_element(_cached))
return element
class Max(ImplicitFunction):
"""Point-wise maximum: h(x,y,z) = max(f(x,y,z), g(x,y,z)).
Parameters
----------
f, g : ImplicitFunction
Operands.
"""
def __repr__(self): return f"Max({self.f}, {self.g})"
def __init__(self, f:ImplicitFunction, g:ImplicitFunction):
self.f = f
self.g = g
def evaluate(self, point):
return np.max((self.f.evaluate(point), self.g.evaluate(point)))
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
element = ET.Element("max")
element.append(self.f.to_xml_element(_cached))
element.append(self.g.to_xml_element(_cached))
return element
# ---------------------------------------------------------------------------
# Cache
# ---------------------------------------------------------------------------
class Cached(ImplicitFunction):
"""
Marks a subexpression for memoisation in C++.
The Python object identity (id()) is used to detect shared nodes during
XML serialisation. This means the SAME Python object must be reused
wherever you want sharing to occur - do NOT construct a new Cached()
at each use site.
Correct - one object, two references:
cx = Cached(2 * np.pi * X())
f = Sin(cx) * Cos(cx) # cx serialised once as <to_cache>
Wrong - two objects, same expression:
f = Sin(Cached(2 * np.pi * X())) * Cos(Cached(2 * np.pi * X()))
"""
def __repr__(self): return f"@[{self.f}]"
def __init__(self, f:ImplicitFunction):
self.f = f
def evaluate(self, point): return self.f.evaluate(point)
def to_xml_element(self, _cached=None):
if _cached is None: _cached = []
key = id(self)
if key in _cached:
node = ET.Element("from_cache")
node.set("id", str(_cached.index(key)))
return node
else:
_cached.append(key)
xml_id = len(_cached) - 1
node = ET.Element("to_cache")
node.set("id", str(xml_id))
node.append(self.f.to_xml_element(_cached))
return node

53
openmc/settings.py
View file

@ -500,6 +500,7 @@ class Settings:
self._use_decay_photons = None
self._random_ray = {}
self._implicit = {}
for key, value in kwargs.items():
setattr(self, key, value)
@ -1491,6 +1492,37 @@ class Settings:
cv.check_greater_than('free gas threshold', free_gas_threshold, 0.0)
self._free_gas_threshold = free_gas_threshold
@property
def implicit(self) -> dict:
return self._implicit
@implicit.setter
def implicit(self, implicit: dict):
if not isinstance(implicit, Mapping):
raise ValueError(f'Unable to set implicit from "{implicit}" '
'which is not a dict.')
for key, value in implicit.items():
if key == 'maxiter':
cv.check_type('maxiter', value, Integral)
cv.check_greater_than('maxiter', value, 0)
elif key == 'atol':
cv.check_type('atol', value, Real)
cv.check_greater_than('atol', value, 0.0)
elif key == 'ftol':
cv.check_type('ftol', value, Real)
cv.check_greater_than('ftol', value, 0.0)
elif key == 'name':
cv.check_type('name', value, str)
cv.check_value('name', value, ['fast', 'naive'])
elif key == 'margin':
cv.check_type('margin', value, Real)
cv.check_greater_than('margin', value, 0.0, True)
else:
raise ValueError(f'Unable to set implicit to "{key}" which is '
'unsupported by OpenMC')
self._implicit = implicit
def _create_run_mode_subelement(self, root):
elem = ET.SubElement(root, "run_mode")
elem.text = self._run_mode.value
@ -2062,6 +2094,13 @@ class Settings:
element = ET.SubElement(root, "free_gas_threshold")
element.text = str(self._free_gas_threshold)
def _create_implicit_solvers_subelement(self, root):
if self._implicit:
element = ET.SubElement(root, "implicit")
for key, value in self._implicit.items():
subelement = ET.SubElement(element, key)
subelement.text = str(value)
def _eigenvalue_from_xml_element(self, root):
elem = root.find('eigenvalue')
if elem is not None:
@ -2560,6 +2599,18 @@ class Settings:
if text is not None:
self.free_gas_threshold = float(text)
def _implicit_solvers_from_xml_element(self, root):
elem = root.find('implicit')
if elem is not None:
self.implicit = {}
for child in elem:
if child.tag in ('atol', 'ftol', 'margin'):
self.implicit[child.tag] = float(child.text)
elif child.tag == 'maxiter':
self.implicit[child.tag] = int(child.text)
elif child.tag == 'name':
self.implicit[child.tag] = str(child.text)
def to_xml_element(self, mesh_memo=None):
"""Create a 'settings' element to be written to an XML file.
@ -2637,6 +2688,7 @@ class Settings:
self._create_use_decay_photons_subelement(element)
self._create_source_rejection_fraction_subelement(element)
self._create_free_gas_threshold_subelement(element)
self._create_implicit_solvers_subelement(element)
# Clean the indentation in the file to be user-readable
clean_indentation(element)
@ -2755,6 +2807,7 @@ class Settings:
settings._use_decay_photons_from_xml_element(elem)
settings._source_rejection_fraction_from_xml_element(elem)
settings._free_gas_threshold_from_xml_element(elem)
settings._implicit_solvers_from_xml_element(elem)
return settings

371
openmc/surface.py
View file

@ -13,6 +13,8 @@ from .checkvalue import check_type, check_value, check_length, check_greater_tha
from .mixin import IDManagerMixin, IDWarning
from .region import Region, Intersection, Union
from .bounding_box import BoundingBox
from . import implicit
from .implicit import ImplicitFunction
from ._xml import get_elem_list, get_text
@ -470,6 +472,10 @@ class Surface(IDManagerMixin, ABC):
coeffs = get_elem_list(elem, "coeffs", float)
kwargs.update(dict(zip(cls._coeff_keys, coeffs)))
if surf_type == "implicit":
kwargs['function'] = ImplicitFunction.from_xml_element(elem.find("function")[0])
kwargs['isovalue'] = float(elem.get("isovalue"))
return cls(**kwargs)
@staticmethod
@ -508,6 +514,16 @@ class Surface(IDManagerMixin, ABC):
surf_type = group['type'][()].decode()
cls = _SURFACE_CLASSES[surf_type]
if surf_type == 'implicit':
xml_str = group['function_xml'][()].decode()
func_elem = ET.fromstring(xml_str)
func = ImplicitFunction.from_xml_element(func_elem[0])
isovalue = float(group['isovalue'][()])
kwargs.update(dict(zip(cls._coeff_keys, coeffs)))
kwargs['function'] = func
kwargs['isovalue'] = isovalue
return cls(**kwargs)
return cls(*coeffs, **kwargs)
@ -2582,6 +2598,359 @@ class ZTorus(TorusMixin, Surface):
elif side == '+':
return BoundingBox.infinite()
class ImplicitSurface(Surface):
_type = 'implicit'
_coeff_keys = ('x0', 'y0', 'z0', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i')
def __init__(self, function:ImplicitFunction, isovalue:float=0., x0=0., y0=0., z0=0., a=1., b=0., c=0., d=0., e=1., f=0., g=0., h=0., i=1., **kwargs):
# Create the surface
super().__init__(**kwargs)
for key, val in zip(self._coeff_keys, (x0, y0, z0, a, b, c, d, e, f, g, h, i)):
setattr(self, key, val)
self.function = function
self.isovalue = float(isovalue)
# Check the implicit surface
if not self.boundary_type == "transmission":
raise ValueError(f"ImplicitSurface boundary type must be 'transmission' but is '{self.boundary_type}'.")
if not self._is_valid_rotation():
raise ValueError(f"Coefficients a,...,i must form a valid rotation matrix.")
if not isinstance(self.function, ImplicitFunction):
raise TypeError(f"func expected type is an ImplicitFunction, but a '{type(function)}' type was given.")
x0 = SurfaceCoefficient('x0')
y0 = SurfaceCoefficient('y0')
z0 = SurfaceCoefficient('z0')
a = SurfaceCoefficient('a')
b = SurfaceCoefficient('b')
c = SurfaceCoefficient('c')
d = SurfaceCoefficient('d')
e = SurfaceCoefficient('e')
f = SurfaceCoefficient('f')
g = SurfaceCoefficient('g')
h = SurfaceCoefficient('h')
i = SurfaceCoefficient('i')
def __repr__(self):
stringlines = super().__repr__().split('\n')
fline = '\n{0: <20}{1}{2}\n'.format('\tFunction', '=\t', self.function)
isoline = '{0: <20}{1}{2}\n'.format('\tIsovalue', '=\t', self.isovalue)
string = "\n".join(stringlines[:4]) + fline + isoline + "\n".join(stringlines[4:])
return string
def _is_valid_rotation(self):
Rmat = self.get_rotation_matrix()
if not np.allclose(Rmat @ Rmat.T, np.identity(3), rtol=0., atol=self._atol): return False
if not np.isclose(np.linalg.det(Rmat), 1.0, rtol=0., atol=self._atol): return False
return True
def get_rotation_matrix(self):
return np.array([[self.a,self.b,self.c],[self.d,self.e,self.f],[self.g,self.h,self.i]])
def bounding_box(self, side):
return BoundingBox.infinite()
def clone(self, memo=None):
if memo is None:
memo = {}
# If no memoize'd clone exists, instantiate one
if self not in memo:
clone = deepcopy(self)
clone.function = self.function
clone.id = None
# Memoize the clone
memo[self] = clone
return memo[self]
def normalize(self, coeffs=None):
if coeffs is None:
coeffs = self._get_base_coeffs()
coeffs = np.asarray(coeffs)
return tuple([c for c in coeffs])
def is_equal(self, other: ImplicitSurface):
coeffs1 = self._get_base_coeffs()
coeffs2 = other._get_base_coeffs()
if not np.allclose(coeffs1, coeffs2, rtol=0., atol=self._atol): return False
if not np.isclose(self.isovalue, other.isovalue, rtol=0., atol=self._atol): return False
if not self.function is other.function: return False
return True
def _get_base_coeffs(self):
return self.x0, self.y0, self.z0, self.a, self.b, self.c, self.d, self.e, self.f, self.g, self.h, self.i
def evaluate(self, point):
Rmat = self.get_rotation_matrix()
point = np.asarray(point) # handles tuple or array
translated = np.array([point[0] - self.x0,
point[1] - self.y0,
point[2] - self.z0]) # (3, ...) for any batch shape
newpoint = np.einsum('ij,j...->i...', Rmat, translated) # rotate, preserving batch dims
return self.function.evaluate(newpoint) - self.isovalue
def translate(self, vector, inplace=False):
if np.allclose(vector, 0., rtol=0., atol=self._atol):
return self if inplace else self.clone()
x0, y0, z0 = self._get_base_coeffs()[:3]
x0 += vector[0]
y0 += vector[1]
z0 += vector[2]
surf = self if inplace else self.clone()
setattr(surf, surf._coeff_keys[0], x0)
setattr(surf, surf._coeff_keys[1], y0)
setattr(surf, surf._coeff_keys[2], z0)
return surf
def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False):
pivot = np.asarray(pivot)
rotation = np.asarray(rotation, dtype=float)
# Allow rotation matrix to be passed in directly, otherwise build it
if rotation.ndim == 2:
check_length('surface rotation', rotation.ravel(), 9)
Rmat = rotation
else:
Rmat = get_rotation_matrix(rotation, order=order)
# Translate surface to pivot
surf = self.translate(-pivot, inplace=inplace)
x0, y0, z0, a, b, c, d, e, f, g, h, i = surf._get_base_coeffs()
# Compute new rotated coefficients a, b, c
newR = surf.get_rotation_matrix() @ Rmat.T
x0, y0, z0 = Rmat @ np.array([x0, y0, z0])
a, b, c = newR[0,:]
d, e, f = newR[1,:]
g, h, i = newR[2,:]
kwargs = {'boundary_type': surf.boundary_type,
'albedo': surf.albedo,
'name': surf.name}
if inplace:
kwargs['surface_id'] = surf.id
surf = ImplicitSurface(surf.function, surf.isovalue, x0, y0, z0, a, b, c, d, e, f, g, h, i, **kwargs)
return surf.translate(pivot, inplace=inplace)
def to_xml_element(self):
root = super().to_xml_element()
root.set("isovalue", str(self.isovalue))
fnode = ET.Element("function")
cached_list = []
fnode.append(self.function.to_xml_element(cached_list))
# Check cached nodes
from_cache_ids = {int(e.get("id")) for e in fnode.iter("from_cache")}
for i in range(len(cached_list)):
if i not in from_cache_ids:
warn(f"Cached node id={i} has no <from_cache> reference and is only used once. Did you forget to reuse it?", UserWarning)
root.append(fnode)
return root
@staticmethod
def from_xml_element(elem):
return Surface.from_xml_element(elem)
@staticmethod
def from_hdf5(group):
return Surface.from_hdf5(group)
"""An implicit surface defined by a user-specified function f(x, y, z) = c.
The surface is the set of points where ``function(R @ (r - r0)) = isovalue``,
where ``r0 = (x0, y0, z0)`` is the translation vector and ``R`` is the
3x3 rotation matrix encoded by coefficients ``a`` through ``i``.
Unlike algebraic surfaces (planes, spheres, quadrics), the function is an
arbitrary smooth expression built from :mod:`openmc.implicit` nodes and
evaluated at runtime by the C++ NaiveLipschitz or FastLipschitz solver.
This makes ImplicitSurface suitable for TPMS geometries (gyroid, Schwartz-P,
diamond) and any other smooth level-set surface.
Parameters
----------
function : ImplicitFunction
Expression tree representing f(x, y, z). Built from nodes in
:mod:`openmc.implicit` (``X()``, ``Sin``, ``Cos``, etc.).
isovalue : float, optional
Level-set value c such that the surface is f = c. Defaults to 0.
x0, y0, z0 : float, optional
Translation of the surface origin in world coordinates. Defaults to 0.
a, b, c, d, e, f, g, h, i : float, optional
Coefficients of the 3x3 rotation matrix R, stored row-major::
R = [[a, b, c],
[d, e, f],
[g, h, i]]
Must form a valid rotation matrix (orthogonal, det = +1).
Defaults to the identity matrix.
boundary_type : str, optional
Must be ``'transmission'`` (the only supported boundary condition).
surface_id : int, optional
Unique identifier. Assigned automatically if not specified.
name : str, optional
Human-readable label.
Raises
------
ValueError
If ``boundary_type`` is not ``'transmission'``.
ValueError
If the a-i coefficients do not form a valid rotation matrix.
TypeError
If ``function`` is not an :class:`~openmc.implicit.ImplicitFunction`.
Notes
-----
**Finite region requirement.** The C++ solver computes the distance to the
implicit surface using the distance to surrounding analytical surfaces as an
upper bound. Every cell containing an ImplicitSurface must therefore also
reference at least one finite analytical surface (plane, sphere, cylinder,
etc.) so that ``Region::distance`` can establish a bound for the solver.
A fatal error is raised at runtime if this condition is not met.
**Caching.** Sub-expressions that appear more than once should be wrapped
in :class:`~openmc.implicit.Cached` to avoid redundant evaluations in the
C++ solver. See the :class:`~openmc.implicit.Cached` docstring for the
correct usage pattern.
**Transform convention.** The surface evaluates the function in local
coordinates ``r_local = R @ (r - r0)``. Translating by ``v`` adds ``v``
to ``r0``. Rotating by ``Rmat`` updates ``R ← R @ Rmat^T`` and
``r0 ← Rmat @ r0``.
Examples
--------
A sphere of radius 5 defined implicitly:
>>> from openmc.implicit import X, Y, Z
>>> func = X()**2 + Y()**2 + Z()**2
>>> sphere = ImplicitSurface(function=func, isovalue=25.)
A Schwartz-P TPMS with pitch 1 cm:
>>> from openmc.surface import TPMS
>>> tpms = TPMS.from_pitch_isovalue('primitive', pitch=1.0, isovalue=0.)
"""
class TPMS(ImplicitSurface):
"""A Triply Periodic Minimal Surface (TPMS) implicit surface.
Convenience subclass of :class:`ImplicitSurface` that constructs the
expression tree for common TPMS families from a pitch length and isovalue.
The resulting surface is periodic in all three Cartesian directions with
the given pitch.
Do not instantiate directly — use the factory method
:meth:`from_pitch_isovalue`.
Notes
-----
TPMS geometries are widely used in nuclear fuel design for their high
surface-area-to-volume ratio, mechanical isotropy, and tuneable porosity.
The isovalue controls the volume fraction: at isovalue = 0 the surface
divides space into two equal-volume phases; positive values shift the
balance toward the positive half-space.
The gyroid and diamond expressions use :class:`~openmc.implicit.Cached`
nodes for the scaled coordinates ``2π x / pitch``, ``2π y / pitch``,
``2π z / pitch``, since each appears in two trigonometric sub-expressions.
Examples
--------
>>> gyroid = TPMS.from_pitch_isovalue('gyroid', pitch=1.0, isovalue=0.)
>>> prim = TPMS.from_pitch_isovalue('primitive', pitch=0.5, isovalue=0.3)
>>> diamond = TPMS.from_pitch_isovalue('diamond', pitch=1.0, isovalue=0.,
... surface_id=5)
"""
@classmethod
def from_pitch_isovalue(cls, tpms: str, pitch: float, isovalue: float,
**kwargs) -> 'TPMS':
"""Construct a TPMS surface from a pitch length and isovalue.
Parameters
----------
tpms : str
Name of the TPMS family. Case-insensitive. Supported values:
+-----------------------+------------------------------------------+
| Name | Equation |
+=======================+==========================================+
| ``'primitive'``, | cos(x') + cos(y') + cos(z') = c |
| ``'schwarz_p'`` | |
+-----------------------+------------------------------------------+
| ``'gyroid'``, | sin(x')cos(z') + sin(y')cos(x') |
| ``'schoen-g'`` | + sin(z')cos(y') = c |
+-----------------------+------------------------------------------+
| ``'diamond'``, | sin(x')cos(y'-z') |
| ``'schwarz_d'`` | + sin(y'+z')cos(x') = c |
+-----------------------+------------------------------------------+
where ``x' = 2π x / pitch``, and similarly for y' and z'.
pitch : float
Spatial period of the surface in [cm]. All three Cartesian
directions share the same pitch.
isovalue : float
Level-set value c. At ``isovalue = 0`` the surface is a true
minimal surface dividing space into two equal-volume phases.
Increasing the isovalue increases the volume fraction of the
negative half-space (the ``-surf`` region).
**kwargs
Additional keyword arguments passed to :class:`ImplicitSurface`
(e.g. ``surface_id``, ``name``, transform coefficients).
Returns
-------
TPMS
Constructed TPMS surface ready for use in cell definitions.
Raises
------
NotImplementedError
If ``tpms`` is not one of the supported names.
Examples
--------
>>> surf = TPMS.from_pitch_isovalue('gyroid', pitch=1.0, isovalue=0.)
>>> surf.evaluate((0., 0., 0.)) # gyroid passes through the origin
0.0
"""
@classmethod
def from_pitch_isovalue(cls, tpms:str, pitch:float, isovalue:float, **kwargs):
# Shortcuts
X, Y, Z = implicit.X, implicit.Y, implicit.Z
Cos, Sin = implicit.Cos, implicit.Sin
Cached = implicit.Cached
# Get cached or uncached variables, for performance improvement
def _get_xyz(cached=False):
x = 2 * np.pi * X() / pitch
y = 2 * np.pi * Y() / pitch
z = 2 * np.pi * Z() / pitch
if cached:
return Cached(x), Cached(y), Cached(z)
else:
return x, y, z
# Choice of TPMS
if tpms.lower() in ["primitive", "schwarz_p"]:
x, y, z = _get_xyz(False)
func = Cos(x) + Cos(y) + Cos(z)
elif tpms.lower() in ["gyroid", "schoen-g"]:
x, y, z = _get_xyz(True)
func = Sin(x)*Cos(z) + Sin(y)*Cos(x) + Sin(z)*Cos(y)
elif tpms.lower() in ["diamond", "schwarz_d"]:
x, y, z = _get_xyz(True)
func = Sin(x)*Cos(y - z) + Sin(y + z)*Cos(x)
else:
raise NotImplementedError(f"The TPMS named '{tpms.lower()}' is not implemented.")
return cls(func, isovalue, **kwargs)
class Halfspace(Region):
"""A positive or negative half-space region.
@ -2840,3 +3209,5 @@ YCone._virtual_base = Cone
ZCone._virtual_base = Cone
Sphere._virtual_base = Sphere
Quadric._virtual_base = Quadric
ImplicitSurface._virtual_base = ImplicitSurface
TPMS._virtual_base = TPMS

67
src/cell.cpp
View file

@ -22,6 +22,7 @@
#include "openmc/material.h"
#include "openmc/nuclide.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -949,6 +950,9 @@ std::pair<double, int32_t> Region::distance(
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
// Implicit surfaces are deferred until we know min_dist from the
// analytical surfaces, which bounds the solver interval.
std::vector<int32_t> implicit_tokens;
for (int32_t token : expression_) {
// Ignore this token if it corresponds to an operator rather than a region.
if (token >= OP_UNION)
@ -956,8 +960,16 @@ std::pair<double, int32_t> Region::distance(
// Calculate the distance to this surface.
// Note the off-by-one indexing
Surface* surf = model::surfaces[std::abs(token) - 1].get();
// Defer implicit surfaces to pass 2.
if (surf->geom_type() == GeometryType::IMP) {
implicit_tokens.push_back(token);
continue;
}
bool coincident {std::abs(token) == std::abs(on_surface)};
double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
double d {surf->distance(r, u, coincident)};
// Check if this distance is the new minimum.
if (d < min_dist) {
@ -968,6 +980,59 @@ std::pair<double, int32_t> Region::distance(
}
}
// Finite region check
if (!implicit_tokens.empty() && min_dist == INFTY) {
// Ensure we are actually in the region: False errors sometimes comes with
// SolidRayTracing
bool isInCell = true;
for (int32_t token : expression_) {
if (token >= OP_UNION)
continue;
Surface* surf = model::surfaces[std::abs(token) - 1].get();
if (surf->geom_type() != GeometryType::CSG)
continue;
double f = surf->evaluate(r);
bool in_halfspace = (token < 0) ? (f < 0.) : (f > 0.);
if (!in_halfspace) {
isInCell = false;
break;
}
}
// If we are actually in the region: throw error, if not, return min dist.
if (isInCell) {
fatal_error(
"An implicit surface belongs to a region with no finite analytical "
"boundary. Implicit surfaces must be enclosed in a finite region "
"defined by standard surfaces (planes, spheres, cylinders, etc.)."
"r=(" +
std::to_string(r.x) + ", " + std::to_string(r.y) + ", " +
std::to_string(r.z) +
")"
"u=(" +
std::to_string(u.x) + ", " + std::to_string(u.y) + ", " +
std::to_string(u.z) + ")");
} else {
return {min_dist, i_surf};
}
}
// Implicit surfaces treatment
for (int32_t token : implicit_tokens) {
auto* surf =
static_cast<SurfaceImplicit*>(model::surfaces[std::abs(token) - 1].get());
bool coincident {std::abs(token) == std::abs(on_surface)};
double d {surf->distance_finite(r, u, coincident, min_dist)};
if (d < min_dist) {
if (min_dist - d >= FP_PRECISION * min_dist) {
min_dist = d;
i_surf = -token;
}
}
}
return {min_dist, i_surf};
}

2
src/geometry.cpp
View file

@ -390,6 +390,8 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
double d_surf = INFINITY;
int32_t level_surf_cross;
array<int, 3> level_lat_trans {};
++step_cache.epoch; // New geometry step: invalidate implicit surface
// subexpression cache.
// Loop over each coordinate level.
for (int i = 0; i < p.n_coord(); i++) {

1129
src/implicit.cpp Normal file
View file

File diff suppressed because it is too large Load diff

173
src/implicit_solvers.cpp Normal file
View file

@ -0,0 +1,173 @@
#include "openmc/implicit_solvers.h"
#include <fmt/core.h>
#include "openmc/constants.h"
#include "openmc/error.h"
namespace openmc {
std::unique_ptr<ImplicitSolver> ImplicitSolver::create(
const std::string& name, int max_iter, double atol, double ftol)
{
if (name == "naive")
return std::make_unique<NaiveLipschitz>(atol, ftol, max_iter);
if (name == "fast")
return std::make_unique<FastLipschitz>(atol, ftol, max_iter);
fatal_error("make_implicit_solver: unknown solver name '" + name +
"'. "
"Valid options are: 'naive', 'simple', 'fast'.");
return nullptr; // unreachable — suppresses compiler warning
}
double NaiveLipschitz::solve(const Implicit& function, Position r, Direction u,
double t0, double t1, double isovalue) const
{
const double L = function.compute_lipschitz(r, u, t0, t1);
// Constant function — root only if already within tolerance at t0.
if (L == 0.0) {
double f0 = function.along_ray(t0, r, u) - isovalue;
return (std::abs(f0) < ftol_) ? t0 : INFTY;
}
double t_curr = t0;
double f_curr = function.along_ray(t_curr, r, u) - isovalue;
for (int i = 0; i < max_iter_; ++i) {
// Safe Lipschitz step: f can only change by L*step over this distance,
// so no root is possible before t_curr + |f_curr| / L.
const double step = std::abs(f_curr) / L;
const double t_next = t_curr + step;
// Stepped past the end of the interval — check endpoint.
if (t_next >= t1) {
double f1 = function.along_ray(t1, r, u) - isovalue;
if (f_curr * f1 < 0.0)
return bisect(function, r, u, t_curr, t1, f_curr, f1, isovalue);
return INFTY;
}
const double f_next = function.along_ray(t_next, r, u) - isovalue;
if (std::abs(f_next) < ftol_)
return t_next;
// Sign change detected — root is bracketed in [t_curr, t_next].
if (f_curr * f_next < 0.0)
return bisect(function, r, u, t_curr, t_next, f_curr, f_next, isovalue);
t_curr = t_next;
f_curr = f_next;
}
// Max iterations reached — skip this surface.
warning(fmt::format("NaiveLipschitz reached max iterations ({})."
" The surface will be skipped. Results could be biased.",
max_iter_));
return INFTY;
}
double NaiveLipschitz::bisect(const Implicit& function, Position r, Direction u,
double ta, double tb, double fa, double fb, double isovalue) const
{
// Invariant: fa and fb have opposite signs, root is in (ta, tb).
// tb is always on the destination side of the crossing.
// We return tb rather than the midpoint so that the returned position
// is unambiguously on the destination side — this prevents sense()
// from placing the particle back in the cell it just left.
while ((tb - ta) > atol_) {
const double tm = 0.5 * (ta + tb);
const double fm = function.along_ray(tm, r, u) - isovalue;
if (fa * fm < 0.0) {
tb = tm;
fb = fm;
} else {
ta = tm;
fa = fm;
}
}
return tb;
}
double FastLipschitz::solve(const Implicit& function, Position r, Direction u,
double t0, double t1, double isovalue) const
{
// Stack of intervals to explore. Stored as (t0, t1) pairs.
// Because the stack is LIFO, pushing the right half first ensures
// the left half (smaller t) is always explored first, so the
// solver finds the smallest root.
std::vector<std::pair<double, double>> stack;
stack.reserve(64);
stack.push_back({t0, t1});
int n_iter = 0;
while (!stack.empty()) {
auto [t0, t1] = stack.back();
stack.pop_back();
// Interval has collapsed to within tolerance — root is at t0.
if (t1 - t0 < atol_)
return t1;
// Compute tight Lipschitz bound and function range over [t0, t1].
const double L = function.compute_lipschitz(r, u, t0, t1);
if (L == 0.0) {
// Constant on this interval — root only if already on surface.
double f0 = function.along_ray(t0, r, u) - isovalue;
if (std::abs(f0) < atol_)
return t0;
continue;
}
auto [fmin_raw, fmax_raw] = function.compute_f_min_max(r, u, t0, t1);
const double fmin = fmin_raw - isovalue;
const double fmax = fmax_raw - isovalue;
// Both bounds have the same sign — no root in this interval.
if (fmin * fmax > 0.0)
continue;
// Endpoint values (needed to compute safe sub-interval).
const double f0 = function.along_ray(t0, r, u) - isovalue;
const double f1 = function.along_ray(t1, r, u) - isovalue;
// Safe start: root cannot exist before tSafeStart.
// Derivation: |f(t)| >= |f0| - L*(t-t0), so f cannot cross zero
// until |f0| - L*(t-t0) <= 0 => t >= t0 + |f0|/L.
// Subtract L*atol_ to account for the tolerance band.
const double tSafeStart = t0 + std::max(0.0, std::abs(f0) - L * atol_) / L;
if (tSafeStart - t0 < atol_)
return tSafeStart; // root is within atol of t0
// Safe end: root cannot exist after tSafeEnd (symmetric argument).
const double tSafeEnd = t1 - std::max(0.0, std::abs(f1) - L * atol_) / L;
if (t1 - tSafeEnd < atol_)
return t1; // root is within atol of t1
// Safe interval has collapsed — root is at tSafeStart.
if (tSafeEnd - tSafeStart < atol_)
return tSafeEnd;
if (++n_iter >= max_iter_) {
warning(
fmt::format("FastLipschitz reached max iterations ({})."
" The surface will be skipped. Results could be biased.",
max_iter_));
return INFTY;
}
// Split the safe interval and push both halves.
// Right is pushed first so that left is popped and explored first.
const double tSplit = 0.5 * (tSafeStart + tSafeEnd);
stack.push_back({tSplit, tSafeEnd}); // right — explored second
stack.push_back({tSafeStart, tSplit}); // left — explored first
}
return INFTY;
}
} // namespace openmc

28
src/settings.cpp
View file

@ -154,6 +154,13 @@ int verbosity {-1};
double weight_cutoff {0.25};
double weight_survive {1.0};
// implicit surface
int implicit_maxiter {1000000};
std::string implicit_solver {"fast"};
double implicit_atol {1.e-9};
double implicit_ftol {1.e-9};
double implicit_margin {1.e-6};
} // namespace settings
//==============================================================================
@ -1343,6 +1350,27 @@ void read_settings_xml(pugi::xml_node root)
settings::use_shared_secondary_bank = true;
}
}
// Implicit surfaces settings
if (check_for_node(root, "implicit")) {
xml_node implicit_node = root.child("implicit");
if (check_for_node(implicit_node, "maxiter")) {
implicit_maxiter = std::stoi(get_node_value(implicit_node, "maxiter"));
}
if (check_for_node(implicit_node, "name")) {
implicit_solver = get_node_value(implicit_node, "name");
if (implicit_solver != "naive" && implicit_solver != "fast")
fatal_error("Unrecognized implicit solver name: " + implicit_solver);
}
if (check_for_node(implicit_node, "atol")) {
implicit_atol = std::stod(get_node_value(implicit_node, "atol"));
}
if (check_for_node(implicit_node, "ftol")) {
implicit_ftol = std::stod(get_node_value(implicit_node, "ftol"));
}
if (check_for_node(implicit_node, "margin")) {
implicit_margin = std::stod(get_node_value(implicit_node, "margin"));
}
}
}
void free_memory_settings()

78
src/surface.cpp
View file

@ -167,8 +167,10 @@ void Surface::to_hdf5(hid_t group_id) const
if (geom_type() == GeometryType::DAG) {
write_string(surf_group, "geom_type", "dagmc", false);
} else if (geom_type() == GeometryType::CSG) {
write_string(surf_group, "geom_type", "csg", false);
} else if (geom_type() == GeometryType::CSG ||
geom_type() == GeometryType::IMP) {
write_string(surf_group, "geom_type",
geom_type() == GeometryType::IMP ? "imp" : "csg", false);
if (bc_) {
write_string(surf_group, "boundary_type", bc_->type(), false);
@ -1168,6 +1170,75 @@ Direction SurfaceZTorus::normal(Position r) const
return n / n.norm();
}
//==============================================================================
// SurfaceImplicit implementation
//==============================================================================
SurfaceImplicit::SurfaceImplicit(pugi::xml_node surf_node) : Surface(surf_node)
{
read_coeffs(surf_node, id_,
{&x0_, &y0_, &z0_, &A_, &B_, &C_, &D_, &E_, &F_, &G_, &H_, &I_});
isovalue_ = surf_node.attribute("isovalue").as_double();
auto func_node = surf_node.child("function");
if (!func_node)
fatal_error(fmt::format("Surface {} missing <function> element.", id_));
function_ = Implicit::from_xml_element(func_node.first_child());
solver_ = ImplicitSolver::create(settings::implicit_solver,
settings::implicit_maxiter, settings::implicit_atol,
settings::implicit_ftol);
}
void SurfaceImplicit::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "implicit", false);
std::array<double, 12> coeffs {
{x0_, y0_, z0_, A_, B_, C_, D_, E_, F_, G_, H_, I_}};
write_dataset(group_id, "coefficients", coeffs);
write_dataset(group_id, "isovalue", isovalue_);
write_string(group_id, "function_xml", function_->to_xml_string(), false);
}
Position SurfaceImplicit::transform(Position r) const
{
double dx = r.x - x0_, dy = r.y - y0_, dz = r.z - z0_;
return Position(A_ * dx + B_ * dy + C_ * dz, D_ * dx + E_ * dy + F_ * dz,
G_ * dx + H_ * dy + I_ * dz);
}
Direction SurfaceImplicit::transform_dir(Direction u) const
{
// Rotation only — no translation for directions
return Direction(A_ * u.x + B_ * u.y + C_ * u.z,
D_ * u.x + E_ * u.y + F_ * u.z, G_ * u.x + H_ * u.y + I_ * u.z);
}
double SurfaceImplicit::evaluate(Position r) const
{
return function_->evaluate(transform(r)) - isovalue_;
}
double SurfaceImplicit::distance_finite(
Position r, Direction u, bool coincident, double distance_max) const
{
double f0 = function_->evaluate(r);
double t0 = (std::abs(f0) <= settings::implicit_ftol || coincident)
? settings::implicit_margin
: 0.0;
Position r_tr = transform(r);
Direction u_tr = transform_dir(u);
return solver_->solve(*function_, r_tr, u_tr, t0, distance_max, isovalue_) +
settings::implicit_margin;
}
double SurfaceImplicit::distance(Position r, Direction u, bool coincident) const
{
fatal_error("SurfaceImplicit::distance: shouldn't be called.");
}
Direction SurfaceImplicit::normal(Position r) const
{
Position r_tr = transform(r);
Gradient g = function_->gradient(r_tr);
double nx = A_ * g.x + D_ * g.y + G_ * g.z;
double ny = B_ * g.x + E_ * g.y + H_ * g.z;
double nz = C_ * g.x + F_ * g.y + I_ * g.z;
double len = std::sqrt(nx * nx + ny * ny + nz * nz);
return {nx / len, ny / len, nz / len};
}
//==============================================================================
void read_surfaces(pugi::xml_node node,
@ -1238,6 +1309,9 @@ void read_surfaces(pugi::xml_node node,
} else if (surf_type == "z-torus") {
model::surfaces.push_back(std::make_unique<SurfaceZTorus>(surf_node));
} else if (surf_type == "implicit") {
model::surfaces.push_back(std::make_unique<SurfaceImplicit>(surf_node));
} else {
fatal_error(fmt::format("Invalid surface type, \"{}\"", surf_type));
}

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

0
tests/regression_tests/implicit/diamond/__init__.py Normal file
View file

80
tests/regression_tests/implicit/diamond/inputs_true.dat Normal file
View file

@ -0,0 +1,80 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="5.0"/>
<nuclide name="U238" ao="95.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-7 1 -2 3 -4 5 -6" universe="1"/>
<cell id="2" material="void" region="7 1 -2 3 -4 5 -6" universe="1"/>
<surface id="1" type="x-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="3" type="y-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="5" type="z-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="7" type="implicit" coeffs="0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" isovalue="0.0">
<function>
<add>
<mul>
<sin>
<to_cache id="0">
<div>
<scale value="6.283185307179586">
<x/>
</scale>
<constant value="8.0"/>
</div>
</to_cache>
</sin>
<cos>
<sub>
<to_cache id="1">
<div>
<scale value="6.283185307179586">
<y/>
</scale>
<constant value="8.0"/>
</div>
</to_cache>
<to_cache id="2">
<div>
<scale value="6.283185307179586">
<z/>
</scale>
<constant value="8.0"/>
</div>
</to_cache>
</sub>
</cos>
</mul>
<mul>
<sin>
<add>
<from_cache id="1"/>
<from_cache id="2"/>
</add>
</sin>
<cos>
<from_cache id="0"/>
</cos>
</mul>
</add>
</function>
</surface>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/diamond/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
2.372756E-01 1.317173E-03

47
tests/regression_tests/implicit/diamond/test.py Normal file
View file

@ -0,0 +1,47 @@
import openmc
import pytest
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def implicit_diamond_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 5.0)
material.add_nuclide('U238', 95.0)
material.set_density('g/cm3', 16.0)
# Diamond
x0 = openmc.XPlane(-16.0, boundary_type="vacuum")
x1 = openmc.XPlane(+16.0, boundary_type="vacuum")
y0 = openmc.YPlane(-16.0, boundary_type="vacuum")
y1 = openmc.YPlane(+16.0, boundary_type="vacuum")
z0 = openmc.ZPlane(-16.0, boundary_type="vacuum")
z1 = openmc.ZPlane(+16.0, boundary_type="vacuum")
box = +x0 & -x1 & +y0 & -y1 & +z0 & -z1
impl = openmc.TPMS.from_pitch_isovalue("diamond", 8., 0.)
fuel_cell = openmc.Cell(region=-impl & box, fill=material)
void_cell = openmc.Cell(region=+impl & box)
geometry = openmc.Geometry([fuel_cell, void_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_implicit_diamond(implicit_diamond_model):
harness = PyAPITestHarness('statepoint.20.h5', model=implicit_diamond_model)
harness.main()

0
tests/regression_tests/implicit/fuel_graded/__init__.py Normal file
View file

72
tests/regression_tests/implicit/fuel_graded/inputs_true.dat Normal file
View file

@ -0,0 +1,72 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-7 1 -2 3 -4 5 -6" universe="1"/>
<cell id="2" material="void" region="7 1 -2 3 -4 5 -6" universe="1"/>
<surface id="1" type="x-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="3" type="y-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="5" type="z-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="7" type="implicit" coeffs="0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" isovalue="0.0">
<function>
<add>
<add>
<cos>
<mul>
<scale value="6.283185307179586">
<x/>
</scale>
<to_cache id="0">
<div>
<add>
<pow value="2">
<z/>
</pow>
<constant value="8.0"/>
</add>
<constant value="64.0"/>
</div>
</to_cache>
</mul>
</cos>
<cos>
<mul>
<scale value="6.283185307179586">
<y/>
</scale>
<from_cache id="0"/>
</mul>
</cos>
</add>
<cos>
<mul>
<scale value="6.283185307179586">
<z/>
</scale>
<from_cache id="0"/>
</mul>
</cos>
</add>
</function>
</surface>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/fuel_graded/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
9.475677E-01 7.264105E-03

52
tests/regression_tests/implicit/fuel_graded/test.py Normal file
View file

@ -0,0 +1,52 @@
import openmc
import pytest
import numpy as np
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z, Cos, Sin, Cached
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def implicit_fuel_graded_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 1.0)
material.set_density('g/cm3', 16.0)
# Box
x0 = openmc.XPlane(-16., boundary_type="vacuum")
x1 = openmc.XPlane(+16., boundary_type="vacuum")
y0 = openmc.YPlane(-16., boundary_type="vacuum")
y1 = openmc.YPlane(+16., boundary_type="vacuum")
z0 = openmc.ZPlane(-16., boundary_type="vacuum")
z1 = openmc.ZPlane(+16., boundary_type="vacuum")
box = +x0 & -x1 & +y0 & -y1 & +z0 & -z1
# Fuel grading
invPitch = Cached((Z() ** 2 + 8.) / 64.)
x = 2 * np.pi * X() * invPitch
y = 2 * np.pi * Y() * invPitch
z = 2 * np.pi * Z() * invPitch
func = Cos(x) + Cos(y) + Cos(z)
impl = ImplicitSurface(function=func)
fuel_cell = openmc.Cell(region=-impl & box, fill=material)
void_cell = openmc.Cell(region=+impl & box)
geometry = openmc.Geometry([fuel_cell, void_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_implicit_fuel_graded(implicit_fuel_graded_model):
harness = PyAPITestHarness('statepoint.20.h5', model=implicit_fuel_graded_model)
harness.main()

0
tests/regression_tests/implicit/fuel_graded_no_cache/__init__.py Normal file
View file

86
tests/regression_tests/implicit/fuel_graded_no_cache/inputs_true.dat Normal file
View file

@ -0,0 +1,86 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-7 1 -2 3 -4 5 -6" universe="1"/>
<cell id="2" material="void" region="7 1 -2 3 -4 5 -6" universe="1"/>
<surface id="1" type="x-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="3" type="y-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="5" type="z-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="7" type="implicit" coeffs="0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" isovalue="0.0">
<function>
<add>
<add>
<cos>
<mul>
<scale value="6.283185307179586">
<x/>
</scale>
<div>
<add>
<pow value="2">
<z/>
</pow>
<constant value="8.0"/>
</add>
<constant value="64.0"/>
</div>
</mul>
</cos>
<cos>
<mul>
<scale value="6.283185307179586">
<y/>
</scale>
<div>
<add>
<pow value="2">
<z/>
</pow>
<constant value="8.0"/>
</add>
<constant value="64.0"/>
</div>
</mul>
</cos>
</add>
<cos>
<mul>
<scale value="6.283185307179586">
<z/>
</scale>
<div>
<add>
<pow value="2">
<z/>
</pow>
<constant value="8.0"/>
</add>
<constant value="64.0"/>
</div>
</mul>
</cos>
</add>
</function>
</surface>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/fuel_graded_no_cache/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
9.475677E-01 7.264105E-03

52
tests/regression_tests/implicit/fuel_graded_no_cache/test.py Normal file
View file

@ -0,0 +1,52 @@
import openmc
import pytest
import numpy as np
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z, Cos, Sin, Cached
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def implicit_fuel_graded_nocache_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 1.0)
material.set_density('g/cm3', 16.0)
# Box
x0 = openmc.XPlane(-16., boundary_type="vacuum")
x1 = openmc.XPlane(+16., boundary_type="vacuum")
y0 = openmc.YPlane(-16., boundary_type="vacuum")
y1 = openmc.YPlane(+16., boundary_type="vacuum")
z0 = openmc.ZPlane(-16., boundary_type="vacuum")
z1 = openmc.ZPlane(+16., boundary_type="vacuum")
box = +x0 & -x1 & +y0 & -y1 & +z0 & -z1
# Fuel grading
invPitch = (Z() ** 2 + 8.) / 64.
x = 2 * np.pi * X() * invPitch
y = 2 * np.pi * Y() * invPitch
z = 2 * np.pi * Z() * invPitch
func = Cos(x) + Cos(y) + Cos(z)
impl = ImplicitSurface(function=func)
fuel_cell = openmc.Cell(region=-impl & box, fill=material)
void_cell = openmc.Cell(region=+impl & box)
geometry = openmc.Geometry([fuel_cell, void_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_implicit_fuel_graded_nocache(implicit_fuel_graded_nocache_model):
harness = PyAPITestHarness('statepoint.20.h5', model=implicit_fuel_graded_nocache_model)
harness.main()

0
tests/regression_tests/implicit/gyroid/__init__.py Normal file
View file

84
tests/regression_tests/implicit/gyroid/inputs_true.dat Normal file
View file

@ -0,0 +1,84 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="5.0"/>
<nuclide name="U238" ao="95.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-7 1 -2 3 -4 5 -6" universe="1"/>
<cell id="2" material="void" region="7 1 -2 3 -4 5 -6" universe="1"/>
<surface id="1" type="x-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="3" type="y-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="5" type="z-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="7" type="implicit" coeffs="0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" isovalue="0.0">
<function>
<add>
<add>
<mul>
<sin>
<to_cache id="0">
<div>
<scale value="6.283185307179586">
<x/>
</scale>
<constant value="4.0"/>
</div>
</to_cache>
</sin>
<cos>
<to_cache id="1">
<div>
<scale value="6.283185307179586">
<z/>
</scale>
<constant value="4.0"/>
</div>
</to_cache>
</cos>
</mul>
<mul>
<sin>
<to_cache id="2">
<div>
<scale value="6.283185307179586">
<y/>
</scale>
<constant value="4.0"/>
</div>
</to_cache>
</sin>
<cos>
<from_cache id="0"/>
</cos>
</mul>
</add>
<mul>
<sin>
<from_cache id="1"/>
</sin>
<cos>
<from_cache id="2"/>
</cos>
</mul>
</add>
</function>
</surface>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/gyroid/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
2.356259E-01 3.813307E-03

47
tests/regression_tests/implicit/gyroid/test.py Normal file
View file

@ -0,0 +1,47 @@
import openmc
import pytest
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def implicit_gyroid_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 5.0)
material.add_nuclide('U238', 95.0)
material.set_density('g/cm3', 16.0)
# Gyroid
x0 = openmc.XPlane(-16.0, boundary_type="vacuum")
x1 = openmc.XPlane(+16.0, boundary_type="vacuum")
y0 = openmc.YPlane(-16.0, boundary_type="vacuum")
y1 = openmc.YPlane(+16.0, boundary_type="vacuum")
z0 = openmc.ZPlane(-16.0, boundary_type="vacuum")
z1 = openmc.ZPlane(+16.0, boundary_type="vacuum")
box = +x0 & -x1 & +y0 & -y1 & +z0 & -z1
impl = openmc.TPMS.from_pitch_isovalue("gyroid", 4., 0.)
fuel_cell = openmc.Cell(region=-impl & box, fill=material)
void_cell = openmc.Cell(region=+impl & box)
geometry = openmc.Geometry([fuel_cell, void_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_implicit_gyroid(implicit_gyroid_model):
harness = PyAPITestHarness('statepoint.20.h5', model=implicit_gyroid_model)
harness.main()

0
tests/regression_tests/implicit/primitive/__init__.py Normal file
View file

63
tests/regression_tests/implicit/primitive/inputs_true.dat Normal file
View file

@ -0,0 +1,63 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="5.0"/>
<nuclide name="U238" ao="95.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-7 1 -2 3 -4 5 -6" universe="1"/>
<cell id="2" material="void" region="7 1 -2 3 -4 5 -6" universe="1"/>
<surface id="1" type="x-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="3" type="y-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="5" type="z-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="7" type="implicit" coeffs="0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" isovalue="0.0">
<function>
<add>
<add>
<cos>
<div>
<scale value="6.283185307179586">
<x/>
</scale>
<constant value="4.0"/>
</div>
</cos>
<cos>
<div>
<scale value="6.283185307179586">
<y/>
</scale>
<constant value="4.0"/>
</div>
</cos>
</add>
<cos>
<div>
<scale value="6.283185307179586">
<z/>
</scale>
<constant value="4.0"/>
</div>
</cos>
</add>
</function>
</surface>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/primitive/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
2.370020E-01 3.809965E-03

47
tests/regression_tests/implicit/primitive/test.py Normal file
View file

@ -0,0 +1,47 @@
import openmc
import pytest
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def implicit_primitive_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 5.0)
material.add_nuclide('U238', 95.0)
material.set_density('g/cm3', 16.0)
# Primitive
x0 = openmc.XPlane(-16.0, boundary_type="vacuum")
x1 = openmc.XPlane(+16.0, boundary_type="vacuum")
y0 = openmc.YPlane(-16.0, boundary_type="vacuum")
y1 = openmc.YPlane(+16.0, boundary_type="vacuum")
z0 = openmc.ZPlane(-16.0, boundary_type="vacuum")
z1 = openmc.ZPlane(+16.0, boundary_type="vacuum")
box = +x0 & -x1 & +y0 & -y1 & +z0 & -z1
impl = openmc.TPMS.from_pitch_isovalue("primitive", 4., 0.)
fuel_cell = openmc.Cell(region=-impl & box, fill=material)
void_cell = openmc.Cell(region=+impl & box)
geometry = openmc.Geometry([fuel_cell, void_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_implicit_primitive(implicit_primitive_model):
harness = PyAPITestHarness('statepoint.20.h5', model=implicit_primitive_model)
harness.main()

0
tests/regression_tests/implicit/sphere/__init__.py Normal file
View file

42
tests/regression_tests/implicit/sphere/inputs_true.dat Normal file
View file

@ -0,0 +1,42 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-2 -1" universe="1"/>
<cell id="2" material="void" region="2 -1" universe="1"/>
<surface id="1" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 11.0"/>
<surface id="2" type="implicit" coeffs="0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" isovalue="100.0">
<function>
<add>
<add>
<pow value="2">
<x/>
</pow>
<pow value="2">
<y/>
</pow>
</add>
<pow value="2">
<z/>
</pow>
</add>
</function>
</surface>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/sphere/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
1.009025E+00 5.384490E-03

42
tests/regression_tests/implicit/sphere/test.py Normal file
View file

@ -0,0 +1,42 @@
import openmc
import pytest
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def implicit_sphere_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 1.0)
material.set_density('g/cm3', 16.0)
# Geometry: implicit sphere enclosed in an analytic vacuum sphere.
# The outer analytic sphere provides the finite region required by
# the two-pass Region::distance algorithm.
R = 10.0
outer = openmc.Sphere(r=R * 1.1, boundary_type='vacuum')
impl = ImplicitSurface(function=X()**2 + Y()**2 + Z()**2, isovalue=R**2)
fuel_cell = openmc.Cell(region=-impl & -outer, fill=material)
void_cell = openmc.Cell(region=+impl & -outer)
geometry = openmc.Geometry([fuel_cell, void_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_implicit_sphere(implicit_sphere_model):
harness = PyAPITestHarness('statepoint.20.h5', model=implicit_sphere_model)
harness.main()

0
tests/regression_tests/implicit/sphere_classic/__init__.py Normal file
View file

26
tests/regression_tests/implicit/sphere_classic/inputs_true.dat Normal file
View file

@ -0,0 +1,26 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-2 -1" universe="1"/>
<cell id="2" material="void" region="2 -1" universe="1"/>
<surface id="1" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 11.0"/>
<surface id="2" type="sphere" coeffs="0.0 0.0 0.0 10.0"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/sphere_classic/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
1.009025E+00 5.384496E-03

39
tests/regression_tests/implicit/sphere_classic/test.py Normal file
View file

@ -0,0 +1,39 @@
import openmc
import pytest
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def classic_sphere_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 1.0)
material.set_density('g/cm3', 16.0)
R = 10.0
outer = openmc.Sphere(r=R * 1.1, boundary_type='vacuum')
impl = openmc.Sphere(r=R)
fuel_cell = openmc.Cell(region=-impl & -outer, fill=material)
void_cell = openmc.Cell(region=+impl & -outer)
geometry = openmc.Geometry([fuel_cell, void_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_classic_sphere(classic_sphere_model):
harness = PyAPITestHarness('statepoint.20.h5', model=classic_sphere_model)
harness.main()

0
tests/regression_tests/implicit/square_circle_fit/__init__.py Normal file
View file

121
tests/regression_tests/implicit/square_circle_fit/inputs_true.dat Normal file
View file

@ -0,0 +1,121 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-7 -8" universe="1"/>
<cell id="2" material="void" region="7 -8" universe="1"/>
<cell id="3" material="void" region="1 -2 3 -4 5 -6 8" universe="1"/>
<surface id="1" type="x-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="3" type="y-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="5" type="z-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="7" type="implicit" coeffs="0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" isovalue="0.0">
<function>
<add>
<add>
<cos>
<div>
<mul>
<scale value="6.283185307179586">
<to_cache id="0">
<scale value="0.125">
<x/>
</scale>
</to_cache>
</scale>
<to_cache id="1">
<sqrt>
<add>
<add>
<add>
<constant value="1.0"/>
<pow value="2">
<from_cache id="0"/>
</pow>
</add>
<pow value="2">
<to_cache id="2">
<scale value="0.125">
<y/>
</scale>
</to_cache>
</pow>
</add>
<pow value="2">
<to_cache id="3">
<scale value="0.125">
<z/>
</scale>
</to_cache>
</pow>
</add>
</sqrt>
</to_cache>
</mul>
<to_cache id="4">
<add>
<constant value="1.0"/>
<max>
<max>
<abs>
<from_cache id="0"/>
</abs>
<abs>
<from_cache id="2"/>
</abs>
</max>
<abs>
<from_cache id="3"/>
</abs>
</max>
</add>
</to_cache>
</div>
</cos>
<cos>
<div>
<mul>
<scale value="6.283185307179586">
<from_cache id="2"/>
</scale>
<from_cache id="1"/>
</mul>
<from_cache id="4"/>
</div>
</cos>
</add>
<cos>
<div>
<mul>
<scale value="6.283185307179586">
<from_cache id="3"/>
</scale>
<from_cache id="1"/>
</mul>
<from_cache id="4"/>
</div>
</cos>
</add>
</function>
</surface>
<surface id="8" type="sphere" coeffs="0.0 0.0 0.0 16.0"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/square_circle_fit/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
8.192616E-01 4.877463E-03

58
tests/regression_tests/implicit/square_circle_fit/test.py Normal file
View file

@ -0,0 +1,58 @@
import openmc
import pytest
import numpy as np
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z, Cos, Abs, Max, Sqrt, Cached
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def implicit_square_circle_fit_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 1.0)
material.set_density('g/cm3', 16.0)
# Square Circle Fit
x0 = openmc.XPlane(-16., boundary_type="vacuum")
x1 = openmc.XPlane(+16., boundary_type="vacuum")
y0 = openmc.YPlane(-16., boundary_type="vacuum")
y1 = openmc.YPlane(+16., boundary_type="vacuum")
z0 = openmc.ZPlane(-16., boundary_type="vacuum")
z1 = openmc.ZPlane(+16., boundary_type="vacuum")
box = +x0 & -x1 & +y0 & -y1 & +z0 & -z1
# Square Circle Fit
L = 8.0
scale = 1.0/L
xp, yp, zp = Cached(scale * X()), Cached(scale * Y()), Cached(scale * Z())
r = Cached(Sqrt(1. + xp**2 + yp**2 + zp**2))
c = Cached(1. + Max(Max(Abs(xp), Abs(yp)), Abs(zp)))
x = 2 * np.pi * xp * r / c
y = 2 * np.pi * yp * r / c
z = 2 * np.pi * zp * r / c
func = Cos(x) + Cos(y) + Cos(z)
impl = ImplicitSurface(function=func)
sphere = openmc.Sphere(r=2*L)
fuel_cell = openmc.Cell(region=-impl & -sphere, fill=material)
void_cell = openmc.Cell(region=+impl & -sphere)
outer_cell = openmc.Cell(region=box & +sphere)
geometry = openmc.Geometry([fuel_cell, void_cell, outer_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_implicit_square_circle_fit(implicit_square_circle_fit_model):
harness = PyAPITestHarness('statepoint.20.h5', model=implicit_square_circle_fit_model)
harness.main()

0
tests/regression_tests/implicit/squircle/__init__.py Normal file
View file

163
tests/regression_tests/implicit/squircle/inputs_true.dat Normal file
View file

@ -0,0 +1,163 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-7 -8" universe="1"/>
<cell id="2" material="void" region="7 -8" universe="1"/>
<cell id="3" material="void" region="1 -2 3 -4 5 -6 8" universe="1"/>
<surface id="1" type="x-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="3" type="y-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="5" type="z-plane" boundary="vacuum" coeffs="-16.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="16.0"/>
<surface id="7" type="implicit" coeffs="0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0" isovalue="0.0">
<function>
<add>
<add>
<cos>
<mul>
<scale value="6.283185307179586">
<to_cache id="0">
<scale value="0.0625">
<x/>
</scale>
</to_cache>
</scale>
<sqrt>
<add>
<sub>
<sub>
<constant value="1.0"/>
<scale value="0.5">
<pow value="2">
<to_cache id="1">
<scale value="0.0625">
<y/>
</scale>
</to_cache>
</pow>
</scale>
</sub>
<scale value="0.5">
<pow value="2">
<to_cache id="2">
<scale value="0.0625">
<z/>
</scale>
</to_cache>
</pow>
</scale>
</sub>
<div>
<mul>
<pow value="2">
<from_cache id="1"/>
</pow>
<pow value="2">
<from_cache id="2"/>
</pow>
</mul>
<constant value="3.0"/>
</div>
</add>
</sqrt>
</mul>
</cos>
<cos>
<mul>
<scale value="6.283185307179586">
<from_cache id="1"/>
</scale>
<sqrt>
<add>
<sub>
<sub>
<constant value="1.0"/>
<scale value="0.5">
<pow value="2">
<from_cache id="0"/>
</pow>
</scale>
</sub>
<scale value="0.5">
<pow value="2">
<from_cache id="2"/>
</pow>
</scale>
</sub>
<div>
<mul>
<pow value="2">
<from_cache id="0"/>
</pow>
<pow value="2">
<from_cache id="2"/>
</pow>
</mul>
<constant value="3.0"/>
</div>
</add>
</sqrt>
</mul>
</cos>
</add>
<cos>
<mul>
<scale value="6.283185307179586">
<from_cache id="2"/>
</scale>
<sqrt>
<add>
<sub>
<sub>
<constant value="1.0"/>
<scale value="0.5">
<pow value="2">
<from_cache id="0"/>
</pow>
</scale>
</sub>
<scale value="0.5">
<pow value="2">
<from_cache id="1"/>
</pow>
</scale>
</sub>
<div>
<mul>
<pow value="2">
<from_cache id="0"/>
</pow>
<pow value="2">
<from_cache id="1"/>
</pow>
</mul>
<constant value="3.0"/>
</div>
</add>
</sqrt>
</mul>
</cos>
</add>
</function>
</surface>
<surface id="8" type="sphere" coeffs="0.0 0.0 0.0 16.0"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>
</model>

2
tests/regression_tests/implicit/squircle/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
9.246343E-01 5.728620E-03

56
tests/regression_tests/implicit/squircle/test.py Normal file
View file

@ -0,0 +1,56 @@
import openmc
import pytest
import numpy as np
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z, Cos, Sqrt, Cached
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def implicit_squircle_model():
# Material
material = openmc.Material()
material.add_nuclide('U235', 1.0)
material.set_density('g/cm3', 16.0)
# Box
x0 = openmc.XPlane(-16., boundary_type="vacuum")
x1 = openmc.XPlane(+16., boundary_type="vacuum")
y0 = openmc.YPlane(-16., boundary_type="vacuum")
y1 = openmc.YPlane(+16., boundary_type="vacuum")
z0 = openmc.ZPlane(-16., boundary_type="vacuum")
z1 = openmc.ZPlane(+16., boundary_type="vacuum")
box = +x0 & -x1 & +y0 & -y1 & +z0 & -z1
# Squircle
L = 8.0
scale = 0.5/L
xp, yp, zp = Cached(scale * X()), Cached(scale * Y()), Cached(scale * Z())
x = 2 * np.pi * xp * Sqrt(1 - 0.5 * yp**2 - 0.5 * zp**2 + yp**2*zp**2/3)
y = 2 * np.pi * yp * Sqrt(1 - 0.5 * xp**2 - 0.5 * zp**2 + xp**2*zp**2/3)
z = 2 * np.pi * zp * Sqrt(1 - 0.5 * xp**2 - 0.5 * yp**2 + xp**2*yp**2/3)
func = Cos(x) + Cos(y) + Cos(z)
impl = ImplicitSurface(function=func)
sphere = openmc.Sphere(r=2*L)
fuel_cell = openmc.Cell(region=-impl & -sphere, fill=material)
void_cell = openmc.Cell(region=+impl & -sphere)
outer_cell = openmc.Cell(region=box & +sphere)
geometry = openmc.Geometry([fuel_cell, void_cell, outer_cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
model = openmc.Model(settings=settings, geometry=geometry)
model.settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0., 0., 0.)),
)
return model
def test_implicit_squircle(implicit_squircle_model):
harness = PyAPITestHarness('statepoint.20.h5', model=implicit_squircle_model)
harness.main()

7
tests/unit_tests/settings.xml Normal file
View file

@ -0,0 +1,7 @@
<?xml version='1.0' encoding='UTF-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<implicit>
<maxiter>0</maxiter>
</implicit>
</settings>

262
tests/unit_tests/test_implicit.py Normal file
View file

@ -0,0 +1,262 @@
"""
Temporary pytest suite for openmc/implicit_function.py
"""
import pytest
import numpy as np
import lxml.etree as ET
from openmc.implicit import (
ImplicitFunction,
X, Y, Z, Constant,
Add, Sub, Mul, Div, Scale, Neg, Pow,
Sin, Cos, Sqrt, Exp, Log, Abs, Min, Max,
Cached,
)
# Canonical test point used throughout
PT = (1.0, 2.0, 3.0)
# ── helpers ────────────────────────────────────────────────────────────────
def xml_roundtrip(func: ImplicitFunction, point=PT) -> float:
"""Serialise → deserialise → evaluate. Must match the original."""
elem = func.to_xml_element([])
restored = ImplicitFunction.from_xml_element(elem)
return restored.evaluate(point)
# ── terminals ──────────────────────────────────────────────────────────────
def test_x(): assert X().evaluate(PT) == 1.0
def test_y(): assert Y().evaluate(PT) == 2.0
def test_z(): assert Z().evaluate(PT) == 3.0
def test_constant(): assert Constant(5.0).evaluate(PT) == 5.0
def test_constant_int(): assert Constant(3).evaluate(PT) == 3.0
# ── arithmetic ─────────────────────────────────────────────────────────────
def test_add(): assert Add(X(), Y()).evaluate(PT) == 3.0
def test_sub(): assert Sub(X(), Y()).evaluate(PT) == -1.0
def test_mul(): assert Mul(X(), Y()).evaluate(PT) == 2.0
def test_div(): assert Div(Y(), X()).evaluate(PT) == 2.0
def test_scale(): assert Scale(X(), 3.0).evaluate(PT) == 3.0
def test_neg(): assert Neg(X()).evaluate(PT) == -1.0
def test_pow_int(): assert Pow(Y(), 2).evaluate(PT) == 4.0
# ── transcendentals ────────────────────────────────────────────────────────
def test_sin(): assert Sin(Constant(np.pi / 2)).evaluate(PT) == pytest.approx(1.0)
def test_cos(): assert Cos(Constant(0.0)).evaluate(PT) == pytest.approx(1.0)
def test_sqrt(): assert Sqrt(Constant(4.0)).evaluate(PT) == pytest.approx(2.0)
def test_exp(): assert Exp(Constant(0.0)).evaluate(PT) == pytest.approx(1.0)
def test_log(): assert Log(Constant(np.e)).evaluate(PT) == pytest.approx(1.0)
def test_abs_pos(): assert Abs(X()).evaluate(PT) == 1.0
def test_abs_neg(): assert Abs(Neg(X())).evaluate(PT) == 1.0
# ── domain guards ──────────────────────────────────────────────────────────
def test_div_zero():
with pytest.raises(ValueError, match="zero"):
Div(X(), Constant(0.0)).evaluate(PT)
def test_sqrt_negative():
with pytest.raises(ValueError, match="negative"):
Sqrt(Constant(-1.0)).evaluate(PT)
def test_log_negative():
with pytest.raises(ValueError, match="negative"):
Log(Constant(-1.0)).evaluate(PT)
def test_pow_float():
with pytest.raises(TypeError):
X() ** 2.0
with pytest.raises(TypeError):
Pow(X(), 0.5)
def test_pow_neg():
with pytest.raises(TypeError):
X() ** -2
def test_pow_non_scalar_via_operator():
with pytest.raises(TypeError):
X() ** Y()
# ── operator overloading ───────────────────────────────────────────────────
def test_add(): assert (X() + Y()).evaluate(PT) == 3.0
def test_radd(): assert (1.0 + X()).evaluate(PT) == 2.0
def test_sub(): assert (X() - Y()).evaluate(PT) == -1.0
def test_rsub(): assert (3.0 - X()).evaluate(PT) == 2.0
def test_mul_func(): assert (X() * Y()).evaluate(PT) == 2.0
def test_mul_scalar(): assert (X() * 3.0).evaluate(PT) == 3.0
def test_rmul_scalar(): assert (3.0 * X()).evaluate(PT) == 3.0
def test_div(): assert (Y() / X()).evaluate(PT) == 2.0
def test_rdiv(): assert (2.0 / Y()).evaluate(PT) == 1.0
def test_neg(): assert (-X()).evaluate(PT) == -1.0
def test_pow(): assert (Y() ** 2).evaluate(PT) == 4.0
def test_scalar_mul_produces_Scale(): assert isinstance(2 * X(), Scale)
def test_func_mul_produces_Mul(): assert isinstance(X() * Y(), Mul)
def test_nested():
f = X()**2 + Y()**2 + Z()**2
assert f.evaluate(PT) == pytest.approx(14.0)
# ── XML tag correctness ────────────────────────────────────────────────────
@pytest.mark.parametrize("node,expected_tag", [
(X(), "x"),
(Y(), "y"),
(Z(), "z"),
(Constant(1.0), "constant"),
(Add(X(), Y()), "add"),
(Sub(X(), Y()), "sub"),
(Mul(X(), Y()), "mul"),
(Div(X(), Y()), "div"),
(Scale(X(), 2.0), "scale"),
(Neg(X()), "neg"),
(Pow(X(), 2), "pow"),
(Sin(X()), "sin"),
(Cos(X()), "cos"),
(Sqrt(X()), "sqrt"),
(Exp(X()), "exp"),
(Log(X()), "log"),
(Abs(X()), "abs"),
(Min(X(), Y()), "min"),
(Max(X(), Y()), "max"),
])
def test_tag(node, expected_tag):
assert node.to_xml_element([]).tag == expected_tag
def test_constant_value_attr():
assert float(Constant(3.14).to_xml_element([]).get("value")) == pytest.approx(3.14)
def test_scale_value_attr():
assert float(Scale(X(), 2.5).to_xml_element([]).get("value")) == pytest.approx(2.5)
def test_pow_value_attr():
assert float(Pow(X(), 3).to_xml_element([]).get("value")) == pytest.approx(3.0)
def test_binary_node_has_two_children():
for binary in [Add, Sub, Mul, Div]:
elem = binary(X(), Y()).to_xml_element([])
assert len(elem) == 2
def test_unary_node_has_one_child():
for unary in [Neg, Sin, Cos, Abs, Exp, Log, Sqrt]:
elem = unary(X()).to_xml_element([])
assert len(elem) == 1
for unary in [Scale, Pow]:
elem = unary(X(), 1).to_xml_element([])
assert len(elem) == 1
def test_wrong_tag():
bad_elem = ET.Element("tanh") # not in the dispatch table
with pytest.raises(ValueError, match="Unknown tag 'tanh'"):
ImplicitFunction.from_xml_element(bad_elem)
# ── Cached serialisation ───────────────────────────────────────────────────
def test_single_use_emits_to_cache():
elem = Cached(X()).to_xml_element([])
assert elem.tag == "to_cache"
assert elem.get("id") == "0"
def test_to_cache_wraps_child():
elem = Cached(X()).to_xml_element([])
assert len(elem) == 1
assert elem[0].tag == "x"
def test_shared_node_first_to_cache_then_from_cache():
node = Cached(X())
_cached = []
first = node.to_xml_element(_cached)
second = node.to_xml_element(_cached)
assert first.tag == "to_cache"
assert second.tag == "from_cache"
assert first.get("id") == second.get("id") == "0"
def test_distinct_cached_nodes_get_distinct_ids():
node1 = Cached(X())
node2 = Cached(Y())
_cached = []
ex = node1.to_xml_element(_cached)
ey = node2.to_xml_element(_cached)
assert ex.get("id") == "0"
assert ey.get("id") == "1"
def test_cached_in_expression():
cx = Cached(X())
elem = (cx + cx).to_xml_element([])
assert elem.tag == "add"
assert elem[0].tag == "to_cache"
assert elem[1].tag == "from_cache"
assert elem[0].get("id") == elem[1].get("id")
# ── XML round-trip ─────────────────────────────────────────────────────────
@pytest.mark.parametrize("func,expected", [
(X(), 1.0),
(Y(), 2.0),
(Z(), 3.0),
(Constant(7.0), 7.0),
(Add(X(), Y()), 3.0),
(Sub(Z(), X()), 2.0),
(Mul(X(), Y()), 2.0),
(Div(Y(), X()), 2.0),
(Scale(Z(), 2.0), 6.0),
(Neg(X()), -1.0),
(Pow(Y(), 2), 4.0),
(Sin(Constant(0.0)), 0.0),
(Cos(Constant(0.0)), 1.0),
(Sqrt(Constant(9.0)), 3.0),
(Exp(Constant(0.0)), 1.0),
(Log(Constant(1.0)), 0.0),
(Abs(Neg(Y())), 2.0),
(X()**2 + Y()**2 + Z()**2, 14.0),
(Min(Y(), X()), 1.0),
(Max(Y(), X()), 2.0),
])
def test_roundtrip(func, expected):
assert xml_roundtrip(func) == pytest.approx(expected)
def test_roundtrip_cached_single():
cx = Cached(X())
f = Sin(cx) + Cos(cx)
assert xml_roundtrip(f) == pytest.approx(np.sin(1.0) + np.cos(1.0))
def test_roundtrip_cached_shared():
cx = Cached(2.0 * X())
f = cx * cx # (2*1)^2 = 4
assert xml_roundtrip(f) == pytest.approx(4.0)
def test_roundtrip_deeply_nested():
f = Sin(Cos(X() + Constant(np.pi)))
assert xml_roundtrip(f) == pytest.approx(f.evaluate(PT))
# ── TPMS-like integration smoke test ──────────────────────────────────────
@pytest.fixture
def gyroid():
L = 1.0
k = 2 * np.pi / L
cx = Cached(k * X())
cy = Cached(k * Y())
cz = Cached(k * Z())
return Sin(cx)*Cos(cz) + Sin(cy)*Cos(cx) + Sin(cz)*Cos(cy)
def test_gyroid_evaluates(gyroid):
# value at origin should be 0 for standard gyroid
assert gyroid.evaluate((0.0, 0.0, 0.0)) == pytest.approx(0.0)
def test_gyroid_roundtrip(gyroid):
pt = (0.1, 0.2, 0.3)
assert xml_roundtrip(gyroid, pt) == pytest.approx(gyroid.evaluate(pt))

173
tests/unit_tests/test_implicit_redundant.py Normal file
View file

@ -0,0 +1,173 @@
"""
Tests for Geometry.remove_redundant_surfaces with ImplicitSurface.
The fix ensures that is_equal() is used within each coefficient-key bucket,
so that ImplicitSurface objects with identical transforms but different
isovalue or function are NOT incorrectly declared redundant.
"""
import numpy as np
import pytest
import openmc
from openmc.surface import ImplicitSurface
from openmc.implicit import X, Y, Z, Sin, Cos, Cached
# ==============================================================================
# Helpers
# ==============================================================================
def make_geometry(*cells):
"""Wrap cells in a minimal geometry."""
universe = openmc.Universe(cells=list(cells))
return openmc.Geometry(universe)
def make_cell(region, fill=None):
mat = openmc.Material()
mat.add_nuclide('U235', 1.0)
mat.set_density('g/cm3', 16.0)
return openmc.Cell(region=region, fill=fill or mat)
# ==============================================================================
# Standard surfaces — existing behaviour preserved
# ==============================================================================
def test_identical_spheres_are_redundant():
"""Two spheres with the same radius are redundant."""
s1 = openmc.Sphere(r=5.0, boundary_type='vacuum')
s2 = openmc.Sphere(r=5.0, boundary_type='vacuum')
c1 = make_cell(-s1)
c2 = make_cell(-s2)
geom = make_geometry(c1, c2)
redundant = geom.remove_redundant_surfaces()
assert len(redundant) == 1
assert s2.id in redundant
assert redundant[s2.id] is s1
def test_different_spheres_are_not_redundant():
"""Two spheres with different radii are not redundant."""
s1 = openmc.Sphere(r=5.0, boundary_type='vacuum')
s2 = openmc.Sphere(r=6.0, boundary_type='vacuum')
c1 = make_cell(-s1)
c2 = make_cell(-s2)
geom = make_geometry(c1, c2)
redundant = geom.remove_redundant_surfaces()
assert len(redundant) == 0
def test_no_redundant_surfaces():
"""Geometry with all unique surfaces returns empty dict."""
s = openmc.Sphere(r=5.0, boundary_type='vacuum')
geom = make_geometry(make_cell(-s))
assert geom.remove_redundant_surfaces() == {}
# ==============================================================================
# ImplicitSurface
# ==============================================================================
def test_same_function_same_isovalue_is_redundant():
"""Two implicit surfaces sharing the same function object and isovalue
are genuinely redundant."""
func = X()**2 + Y()**2 + Z()**2
outer = openmc.Sphere(r=11.0, boundary_type='vacuum')
s1 = ImplicitSurface(function=func, isovalue=25.)
s2 = ImplicitSurface(function=func, isovalue=25.)
c1 = make_cell(-s1 & -outer)
c2 = make_cell(-s2 & -outer)
geom = make_geometry(c1, c2)
redundant = geom.remove_redundant_surfaces()
assert len(redundant) == 1
assert s2.id in redundant
assert redundant[s2.id] is s1
def test_same_function_different_isovalue_is_not_redundant():
"""Same function, different isovalue — different surfaces.
This was the bug: both surfaces had identical _coefficients keys
but different isovalues, yet were incorrectly declared redundant."""
func = X()**2 + Y()**2 + Z()**2
outer = openmc.Sphere(r=11.0, boundary_type='vacuum')
s1 = ImplicitSurface(function=func, isovalue=25.) # r = 5
s2 = ImplicitSurface(function=func, isovalue=100.) # r = 10
c1 = make_cell(-s1 & -outer)
c2 = make_cell(-s2 & -outer)
geom = make_geometry(c1, c2)
redundant = geom.remove_redundant_surfaces()
assert len(redundant) == 0
def test_different_function_objects_not_redundant():
"""Equivalent but distinct function objects are not redundant —
is_equal uses Python object identity for the function."""
outer = openmc.Sphere(r=11.0, boundary_type='vacuum')
s1 = ImplicitSurface(function=X()**2 + Y()**2 + Z()**2, isovalue=25.)
s2 = ImplicitSurface(function=X()**2 + Y()**2 + Z()**2, isovalue=25.)
c1 = make_cell(-s1 & -outer)
c2 = make_cell(-s2 & -outer)
geom = make_geometry(c1, c2)
redundant = geom.remove_redundant_surfaces()
assert len(redundant) == 0
def test_different_transforms_not_redundant():
"""Same function and isovalue but different x0 — not redundant."""
func = X()**2 + Y()**2 + Z()**2
outer = openmc.Sphere(r=11.0, boundary_type='vacuum')
s1 = ImplicitSurface(function=func, isovalue=25., x0=0.)
s2 = ImplicitSurface(function=func, isovalue=25., x0=1.)
c1 = make_cell(-s1 & -outer)
c2 = make_cell(-s2 & -outer)
geom = make_geometry(c1, c2)
redundant = geom.remove_redundant_surfaces()
assert len(redundant) == 0
def test_tpms_same_function_object_is_redundant():
"""Two TPMS surfaces sharing the same underlying function object
and isovalue are correctly identified as redundant."""
func = Sin(X()) * Cos(Z()) + Sin(Y()) * Cos(X()) + Sin(Z()) * Cos(Y())
outer = openmc.Sphere(r=11.0, boundary_type='vacuum')
s1 = ImplicitSurface(function=func, isovalue=0.)
s2 = ImplicitSurface(function=func, isovalue=0.)
c1 = make_cell(-s1 & -outer)
c2 = make_cell(-s2 & -outer)
geom = make_geometry(c1, c2)
redundant = geom.remove_redundant_surfaces()
assert len(redundant) == 1
assert redundant[s2.id] is s1
# ==============================================================================
# Mixed geometry
# ==============================================================================
def test_standard_and_implicit_both_handled():
"""In a geometry with both standard and implicit surfaces,
redundancy detection works correctly for both types independently."""
# Standard redundant pair
sphere_a1 = openmc.Sphere(r=5.0)
sphere_a2 = openmc.Sphere(r=5.0)
# Standard non-redundant pair
sphere_b = openmc.Sphere(r=8.0)
# Implicit non-redundant pair (different isovalue)
func = X()**2 + Y()**2 + Z()**2
outer = openmc.Sphere(r=11.0, boundary_type='vacuum')
impl1 = ImplicitSurface(function=func, isovalue=25.)
impl2 = ImplicitSurface(function=func, isovalue=64.)
c1 = make_cell(-sphere_a1 & -outer)
c2 = make_cell(-sphere_a2 & -outer)
c3 = make_cell(-sphere_b & -outer)
c4 = make_cell(-impl1 & -outer)
c5 = make_cell(-impl2 & -outer)
geom = make_geometry(c1, c2, c3, c4, c5)
redundant = geom.remove_redundant_surfaces()
# Only sphere_a2 is redundant
assert len(redundant) == 1
assert sphere_a2.id in redundant
assert impl1.id not in redundant
assert impl2.id not in redundant

213
tests/unit_tests/test_implicit_settings.py Normal file
View file

@ -0,0 +1,213 @@
"""
Unit tests for implicit surface solver settings.
Tests:
- Default values
- Valid assignments
- Invalid assignments raise errors
- XML round-trip preserves all values
- ImplicitSurface uses the solver specified in Settings
"""
import pytest
from collections.abc import Mapping
import openmc
# ==============================================================================
# Default value
# ==============================================================================
def test_implicit_default_is_none_or_dict():
"""Before assignment, implicit should be None or an empty dict."""
s = openmc.Settings()
assert s.implicit is None or isinstance(s.implicit, Mapping)
# ==============================================================================
# Valid assignments
# ==============================================================================
def test_empty_dict():
"""Empty dict is valid — all values remain at C++ defaults."""
s = openmc.Settings()
s.implicit = {}
assert s.implicit == {}
@pytest.mark.parametrize("name", ["naive", "fast"])
def test_valid_solver_names(name):
s = openmc.Settings()
s.implicit = {'name': name}
assert s.implicit['name'] == name
def test_set_atol():
s = openmc.Settings()
s.implicit = {'atol': 1e-10}
assert s.implicit['atol'] == pytest.approx(1e-10)
def test_set_ftol():
s = openmc.Settings()
s.implicit = {'ftol': 1e-6}
assert s.implicit['ftol'] == pytest.approx(1e-6)
def test_set_maxiter():
s = openmc.Settings()
s.implicit = {'maxiter': 5000}
assert s.implicit['maxiter'] == 5000
def test_set_margin():
s = openmc.Settings()
s.implicit = {'margin': 2e-8}
assert s.implicit['margin'] == pytest.approx(2e-8)
def test_set_all_fields():
s = openmc.Settings()
s.implicit = {
'name': 'fast',
'atol': 1e-9,
'ftol': 1e-6,
'maxiter': 500,
'margin': 3e-8,
}
assert s.implicit['name'] == 'fast'
assert s.implicit['atol'] == pytest.approx(1e-9)
assert s.implicit['ftol'] == pytest.approx(1e-6)
assert s.implicit['maxiter'] == 500
assert s.implicit['margin'] == pytest.approx(3e-8)
# ==============================================================================
# Invalid assignments
# ==============================================================================
def test_not_a_dict():
"""Setting implicit to a non-Mapping raises ValueError."""
s = openmc.Settings()
with pytest.raises(ValueError):
s.implicit = 'fast'
def test_not_a_dict_list():
s = openmc.Settings()
with pytest.raises(ValueError):
s.implicit = ['fast', 1e-8]
def test_unknown_key():
"""Unknown keys raise ValueError."""
s = openmc.Settings()
with pytest.raises(ValueError):
s.implicit = {'unknown_key': 1}
def test_atol_negative():
s = openmc.Settings()
with pytest.raises((ValueError, Exception)):
s.implicit = {'atol': -1e-8}
def test_atol_zero():
s = openmc.Settings()
with pytest.raises((ValueError, Exception)):
s.implicit = {'atol': 0.0}
def test_atol_wrong_type():
s = openmc.Settings()
with pytest.raises((ValueError, TypeError, Exception)):
s.implicit = {'atol': 'small'}
def test_maxiter_negative():
s = openmc.Settings()
with pytest.raises((ValueError, Exception)):
s.implicit = {'maxiter': -1}
def test_maxiter_zero():
s = openmc.Settings()
with pytest.raises((ValueError, Exception)):
s.implicit = {'maxiter': 0}
def test_maxiter_wrong_type():
s = openmc.Settings()
with pytest.raises((ValueError, TypeError, Exception)):
s.implicit = {'maxiter': 1.5}
def test_ftol_negative():
s = openmc.Settings()
with pytest.raises((ValueError, Exception)):
s.implicit = {'ftol': -1.0}
def test_margin_negative():
s = openmc.Settings()
with pytest.raises((ValueError, Exception)):
s.implicit = {'margin': -1e-7}
def test_name_wrong_type():
s = openmc.Settings()
with pytest.raises((ValueError, TypeError, Exception)):
s.implicit = {'name': 42}
def test_name_wrong_name():
s = openmc.Settings()
with pytest.raises((ValueError, TypeError, Exception)):
s.implicit = {'name': 'screugneugneu'}
# ==============================================================================
# XML round-trip
# ==============================================================================
def _roundtrip(settings):
elem = settings.to_xml_element()
restored = openmc.Settings()
return restored.from_xml_element(elem)
def test_name_preserved():
s = openmc.Settings()
s.implicit = {'name': 'naive'}
r = _roundtrip(s)
assert r.implicit['name'] == 'naive'
def test_atol_preserved():
s = openmc.Settings()
s.implicit = {'atol': 1e-10}
r = _roundtrip(s)
assert r.implicit['atol'] == pytest.approx(1e-10)
def test_maxiter_preserved():
s = openmc.Settings()
s.implicit = {'maxiter': 2000}
r = _roundtrip(s)
assert r.implicit['maxiter'] == 2000
def test_ftol_preserved():
s = openmc.Settings()
s.implicit = {'ftol': 1e-6}
r = _roundtrip(s)
assert r.implicit['ftol'] == pytest.approx(1e-6)
def test_margin_preserved():
s = openmc.Settings()
s.implicit = {'margin': 2e-7}
r = _roundtrip(s)
assert r.implicit['margin'] == pytest.approx(2e-7)
def test_implicit_element_present_in_xml():
"""<implicit> element must appear in settings XML."""
s = openmc.Settings()
s.implicit = {'name': 'fast', 'atol': 1e-9}
elem = s.to_xml_element()
assert elem.find('implicit') is not None
def test_full_roundtrip():
s = openmc.Settings()
s.implicit = {
'name': 'naive',
'atol': 1e-9,
'ftol': 1e-6,
'maxiter': 300,
'margin': 5e-8,
}
r = _roundtrip(s)
assert r.implicit['name'] == 'naive'
assert r.implicit['atol'] == pytest.approx(1e-9)
assert r.implicit['ftol'] == pytest.approx(1e-6)
assert r.implicit['maxiter'] == 300
assert r.implicit['margin'] == pytest.approx(5e-8)

563
tests/unit_tests/test_implicit_surface.py Normal file
View file

@ -0,0 +1,563 @@
"""
Tests for openmc.ImplicitSurface and openmc.TPMS.
Run with: pytest test_implicit_surface.py -v
"""
import warnings
import numpy as np
import h5py
import lxml.etree as ET
import pytest
import openmc.implicit as impl
from openmc.implicit import (
X, Y, Z,
Add, Mul, Scale, Sin, Cos, Cached,
ImplicitFunction,
)
from openmc.surface import ImplicitSurface, TPMS
# ==============================================================================
# Construction
# ==============================================================================
def test_valid_minimal():
"""Only the function is required; all other params default correctly."""
func = X()**2 + Y()**2 + Z()**2
surf = ImplicitSurface(function=func)
assert surf.isovalue == 0.0
assert surf.x0 == 0.0
assert surf.y0 == 0.0
assert surf.z0 == 0.0
assert np.allclose(surf.get_rotation_matrix(), np.eye(3))
def test_valid_full_params():
"""All twelve transform parameters plus isovalue are stored."""
surf = ImplicitSurface(
function=X(), isovalue=3.14,
x0=1., y0=2., z0=3.,
a=1., b=0., c=0.,
d=0., e=1., f=0.,
g=0., h=0., i=1.,
)
assert surf.isovalue == pytest.approx(3.14)
assert surf.x0 == pytest.approx(1.)
assert surf.y0 == pytest.approx(2.)
assert surf.z0 == pytest.approx(3.)
def test_isovalue_stored_as_float():
"""isovalue is always stored as float even when passed as int."""
surf = ImplicitSurface(function=X(), isovalue=5)
assert isinstance(surf.isovalue, float)
def test_default_rotation_is_identity():
"""Default a-i coefficients encode the 3x3 identity."""
surf = ImplicitSurface(function=X())
assert np.allclose(surf.get_rotation_matrix(), np.eye(3))
def test_invalid_boundary_vacuum():
"""Vacuum boundary type raises ValueError."""
with pytest.raises(ValueError, match="transmission"):
ImplicitSurface(function=X(), boundary_type='vacuum')
def test_invalid_boundary_reflective():
"""Reflective boundary type raises ValueError."""
with pytest.raises(ValueError, match="transmission"):
ImplicitSurface(function=X(), boundary_type='reflective')
def test_invalid_rotation_non_orthogonal():
"""Non-orthogonal a-i matrix raises ValueError."""
with pytest.raises(ValueError, match="rotation"):
ImplicitSurface(function=X(),
a=2., b=0., c=0.,
d=0., e=1., f=0.,
g=0., h=0., i=1.)
def test_invalid_rotation_det_minus_one():
"""Orthogonal matrix with det = -1 (reflection) raises ValueError."""
with pytest.raises(ValueError, match="rotation"):
ImplicitSurface(function=X(),
a=-1., b=0., c=0.,
d=0., e=1., f=0.,
g=0., h=0., i=1.)
def test_invalid_function_type_string():
"""Passing a string as function raises TypeError."""
with pytest.raises(TypeError):
ImplicitSurface(function="not a function")
def test_invalid_function_type_none():
"""Passing None as function raises TypeError."""
with pytest.raises(TypeError):
ImplicitSurface(function=None)
def test_repr_does_not_crash():
"""repr() returns a non-empty string that mentions Function and Isovalue."""
surf = ImplicitSurface(function=X()**2 + Y()**2 + Z()**2, isovalue=25.)
s = repr(surf)
assert len(s) > 0
assert "Function" in s
assert "Isovalue" in s
def test_normalize():
surf = ImplicitSurface(
function=X(), isovalue=3.14,
x0=1., y0=2., z0=3.)
coeffs = surf.normalize()
assert coeffs[0]==1.
assert coeffs[1]==2.
assert coeffs[2]==3.
# ==============================================================================
# evaluate
# ==============================================================================
# Sphere: f = x^2 + y^2 + z^2, isovalue = R^2
R = 5.0
@pytest.fixture
def sphere():
func = X()**2 + Y()**2 + Z()**2
return ImplicitSurface(function=func, isovalue=R**2)
def test_interior_negative(sphere):
"""Origin is inside the sphere → evaluate < 0."""
assert sphere.evaluate((0., 0., 0.)) < 0.
def test_exterior_positive(sphere):
"""Point clearly outside sphere → evaluate > 0."""
assert sphere.evaluate((10., 0., 0.)) > 0.
def test_on_surface_zero(sphere):
"""Point exactly on sphere surface → evaluate ≈ 0."""
assert sphere.evaluate((R, 0., 0.)) == pytest.approx(0., abs=1e-10)
assert sphere.evaluate((0., R, 0.)) == pytest.approx(0., abs=1e-10)
assert sphere.evaluate((0., 0., R)) == pytest.approx(0., abs=1e-10)
def test_isovalue_shifts_zero_level():
"""Larger isovalue makes previously-exterior points interior."""
func = X()**2 + Y()**2 + Z()**2
s_small = ImplicitSurface(function=func, isovalue=25.) # r = 5
s_large = ImplicitSurface(function=func, isovalue=100.) # r = 10
pt = (7., 0., 0.) # between the two radii
assert s_small.evaluate(pt) > 0.
assert s_large.evaluate(pt) < 0.
def test_translation_shifts_surface():
"""x0 = 5 moves the sphere centre to (5, 0, 0)."""
func = X()**2 + Y()**2 + Z()**2
surf = ImplicitSurface(function=func, isovalue=16., x0=5.)
assert surf.evaluate((9., 0., 0.)) == pytest.approx(0., abs=1e-10) # on surface
assert surf.evaluate((5., 0., 0.)) < 0. # centre → inside
assert surf.evaluate((5., 0., 0.)) == pytest.approx(-16., abs=1e-10) # centre value
assert surf.evaluate((0., 0., 0.)) > 0. # outside
def test_rotation_applied_correctly():
"""90° z-rotation of the x-plane (f = x) makes f = y in world coords.
Rotation matrix [[0,-1,0],[1,0,0],[0,0,1]] maps r_local.x = -y,
so evaluate(r) = X().evaluate(r_local) = -y.
"""
surf = ImplicitSurface(function=X(), isovalue=0.)
rotated = surf.rotate([0., 0., 90.])
assert rotated.evaluate((0., 1., 0.)) == pytest.approx( 1., abs=1e-10)
assert rotated.evaluate((0., -1., 0.)) == pytest.approx(-1., abs=1e-10)
assert rotated.evaluate((5., 0., 0.)) == pytest.approx( 0., abs=1e-10) # on surface
def test_matrix_rotation_applied_correctly():
"""Same as precedent but with a matrix"""
surf = ImplicitSurface(function=X(), isovalue=0.)
rmat = [[0,-1,0],[1,0,0],[0,0,1]]
rotated = surf.rotate(rmat)
assert rotated.evaluate((0., 1., 0.)) == pytest.approx( 1., abs=1e-10)
assert rotated.evaluate((0., -1., 0.)) == pytest.approx(-1., abs=1e-10)
assert rotated.evaluate((5., 0., 0.)) == pytest.approx( 0., abs=1e-10) # on surface
# ==============================================================================
# translate
# ==============================================================================
def test_updates_origin_components():
"""translate([1, 2, 3]) increments x0, y0, z0 correctly."""
surf = ImplicitSurface(function=X())
t = surf.translate([1., 2., 3.])
assert t.x0 == pytest.approx(1.)
assert t.y0 == pytest.approx(2.)
assert t.z0 == pytest.approx(3.)
def test_returns_new_object_by_default():
"""Default inplace=False returns a different Python object."""
surf = ImplicitSurface(function=X())
t = surf.translate([1., 0., 0.])
assert t is not surf
def test_original_unchanged():
"""Original surface is not modified when inplace=False."""
surf = ImplicitSurface(function=X())
surf.translate([5., 0., 0.])
assert surf.x0 == pytest.approx(0.)
def test_inplace_modifies_same_object():
"""inplace=True modifies and returns the same object."""
surf = ImplicitSurface(function=X())
result = surf.translate([3., 0., 0.], inplace=True)
assert result is surf
assert surf.x0 == pytest.approx(3.)
def test_zero_vector_returns_clone():
"""translate([0,0,0]) returns a clone (new object), not self.
This differs from the base-class behaviour where the same object is
returned for zero translations.
"""
surf = ImplicitSurface(function=X())
result = surf.translate([0., 0., 0.])
assert result is not surf
def test_preserves_function_reference():
"""translate does not copy or modify the ImplicitFunction DAG."""
func = X()**2 + Y()**2 + Z()**2
surf = ImplicitSurface(function=func)
t = surf.translate([1., 0., 0.])
assert t.function is func
def test_evaluate_consistency():
"""After translate([5,0,0]), sphere centre moves accordingly."""
func = X()**2 + Y()**2 + Z()**2
surf = ImplicitSurface(function=func, isovalue=16.)
t = surf.translate([5., 0., 0.])
assert t.evaluate((9., 0., 0.)) == pytest.approx(0., abs=1e-10) # on surface
assert t.evaluate((5., 0., 0.)) == pytest.approx(-16., abs=1e-10) # centre
assert t.evaluate((0., 0., 0.)) > 0. # outside
# ==============================================================================
# rotate
# ==============================================================================
def test_euler_angles_update_rotation_matrix():
"""rotate([0, 0, 90]) applies Rz(90°) to the rotation matrix."""
surf = ImplicitSurface(function=X())
rotated = surf.rotate([0., 0., 90.])
R90z = np.array([[0., -1., 0.], [1., 0., 0.], [0., 0., 1.]])
expected = R90z.T
assert np.allclose(rotated.get_rotation_matrix(), expected, atol=1e-10)
def test_rotation_matrix_passed_directly():
"""A 3x3 ndarray can be passed directly as the rotation."""
R90z = np.array([[0., -1., 0.], [1., 0., 0.], [0., 0., 1.]])
surf = ImplicitSurface(function=X())
rotated = surf.rotate(R90z)
assert np.allclose(rotated.get_rotation_matrix(), R90z.T, atol=1e-10)
@pytest.mark.parametrize("angles", [
[30., 0., 0.], [0., 45., 0.], [0., 0., 90.],
[30., 45., 60.], [15., -30., 75.],
])
def test_result_is_always_valid_rotation(angles):
"""Rotated surface always has an orthogonal matrix with det = +1."""
surf = ImplicitSurface(function=X())
rotated = surf.rotate(angles)
R = rotated.get_rotation_matrix()
assert np.allclose(R @ R.T, np.eye(3), atol=1e-10), \
f"Not orthogonal for angles {angles}"
assert np.isclose(np.linalg.det(R), 1.0, atol=1e-10), \
f"det ≠ 1 for angles {angles}"
def test_preserves_function_reference():
"""rotate does not copy or modify the ImplicitFunction DAG."""
func = X()
surf = ImplicitSurface(function=func)
assert surf.rotate([0., 0., 45.]).function is func
def test_pivot_translates_origin():
"""Rotating about a non-origin pivot moves the surface origin correctly.
Sphere centred at (5,0,0); 90° z-rotation about origin → centre at (0,5,0).
"""
func = X()**2 + Y()**2 + Z()**2
surf = ImplicitSurface(function=func, isovalue=4., x0=5.)
rotated = surf.rotate([0., 0., 90.], pivot=(0., 0., 0.))
assert rotated.x0 == pytest.approx(0., abs=1e-10)
assert rotated.y0 == pytest.approx(5., abs=1e-10)
assert rotated.z0 == pytest.approx(0., abs=1e-10)
def test_rotate_inplace_modifies_same_object():
"""inplace=True modifies and returns the same object."""
surf = ImplicitSurface(function=X(), x0=1.)
result = surf.rotate([0., 0., 90.], inplace=True)
assert result.id == surf.id
assert result.y0 == pytest.approx(1.)
# ==============================================================================
# is_equal and clone
# ==============================================================================
def test_equal_to_itself():
func = X()**2 + Y()**2 + Z()**2
surf = ImplicitSurface(function=func, isovalue=25.)
assert surf.is_equal(surf)
def test_equal_same_function_object_same_params():
"""Two surfaces sharing the same Python function object and params."""
func = X()
s1 = ImplicitSurface(function=func, isovalue=1.)
s2 = ImplicitSurface(function=func, isovalue=1.)
assert s1.is_equal(s2)
def test_not_equal_different_isovalue():
func = X()
s1 = ImplicitSurface(function=func, isovalue=1.)
s2 = ImplicitSurface(function=func, isovalue=2.)
assert not s1.is_equal(s2)
def test_not_equal_different_translation():
func = X()
s1 = ImplicitSurface(function=func, x0=0.)
s2 = ImplicitSurface(function=func, x0=1.)
assert not s1.is_equal(s2)
def test_not_equal_different_function_objects():
"""is_equal uses Python object identity for the function — two
equivalent but distinct objects are NOT equal."""
s1 = ImplicitSurface(function=X())
s2 = ImplicitSurface(function=X())
assert not s1.is_equal(s2)
def test_clone_has_new_id():
surf = ImplicitSurface(function=X())
clone = surf.clone()
assert clone.id != surf.id
def test_clone_shares_function_reference():
"""clone preserves the Python function object, not a deep copy."""
func = X()**2 + Y()**2 + Z()**2
surf = ImplicitSurface(function=func, isovalue=25.)
assert surf.clone().function is func
def test_clone_independent_coefficients():
"""Translating a clone does not affect the original."""
surf = ImplicitSurface(function=X())
clone = surf.clone()
clone.translate([5., 0., 0.], inplace=True)
assert surf.x0 == pytest.approx(0.)
# ==============================================================================
# bounding_box
# ==============================================================================
def test_always_infinite():
"""bounding_box returns an infinite box on both sides — no analytical
bound is available without a user-supplied bbox."""
surf = ImplicitSurface(function=X()**2 + Y()**2 + Z()**2, isovalue=25.)
bb_pos = surf.bounding_box('+')
bb_neg = surf.bounding_box('-')
assert np.all(bb_pos.lower_left == -np.inf)
assert np.all(bb_neg.lower_left == -np.inf)
assert np.all(bb_pos.upper_right == +np.inf)
assert np.all(bb_neg.upper_right == +np.inf)
# ==============================================================================
# XML round-trip
# ==============================================================================
def test_roundtrip_basic():
"""to_xml_element + from_xml_element produces a surface that evaluates
identically."""
func = X()**2 + Y()**2 + Z()**2
surf = ImplicitSurface(function=func, isovalue=25.)
elem = surf.to_xml_element()
restored = ImplicitSurface.from_xml_element(elem)
for pt in [(3., 4., 0.), (0., 0., 0.), (5., 0., 0.)]:
assert restored.evaluate(pt) == pytest.approx(surf.evaluate(pt), abs=1e-10)
def test_roundtrip_isovalue_preserved():
"""isovalue survives the XML round-trip exactly."""
surf = ImplicitSurface(function=X(), isovalue=3.14)
restored = ImplicitSurface.from_xml_element(surf.to_xml_element())
assert restored.isovalue == pytest.approx(3.14)
def test_roundtrip_transform_params_preserved():
"""All 12 transform coefficients survive the round-trip."""
surf = ImplicitSurface(function=X(), isovalue=0.,
x0=1., y0=2., z0=3.,
a=1., b=0., c=0.,
d=0., e=1., f=0.,
g=0., h=0., i=1.)
restored = ImplicitSurface.from_xml_element(surf.to_xml_element())
assert np.allclose(restored._get_base_coeffs(),
surf._get_base_coeffs(), atol=1e-10)
def test_roundtrip_with_cached_nodes():
"""Surface with Cached nodes evaluates identically after round-trip."""
cx = Cached(2. * X())
surf = ImplicitSurface(function=Sin(cx) * Cos(cx), isovalue=0.)
restored = ImplicitSurface.from_xml_element(surf.to_xml_element())
for pt in [(0., 0., 0.), (0.5, 0., 0.), (1., 0., 0.)]:
assert restored.evaluate(pt) == pytest.approx(surf.evaluate(pt), abs=1e-10)
def test_roundtrip_produces_implicit_surface_not_tpms():
"""from_xml_element always produces ImplicitSurface (not TPMS) because
TPMS is not a direct subclass of Surface and is absent from
_SURFACE_CLASSES."""
surf = TPMS.from_pitch_isovalue("gyroid", 1.0, 0.0)
restored = ImplicitSurface.from_xml_element(surf.to_xml_element())
assert type(restored) is ImplicitSurface
def test_single_use_cached_warns():
"""Cached node used only once triggers a UserWarning during serialisation."""
surf = ImplicitSurface(function=Sin(Cached(X())))
with pytest.warns(UserWarning, match="only used once"):
surf.to_xml_element()
def test_reused_cached_no_warning():
"""Cached node used twice (shared) triggers no UserWarning."""
cx = Cached(X())
surf = ImplicitSurface(function=Sin(cx) + Cos(cx))
with warnings.catch_warnings():
warnings.simplefilter("error", UserWarning)
surf.to_xml_element() # must not raise
def test_implicit_surface_from_hdf5():
func = X()**2 + Y()**2 + Z()**2
orig = ImplicitSurface(function=func, isovalue=25.)
# Build the XML string that C++ to_xml_string() produces
fnode = ET.Element("function")
fnode.append(func.to_xml_element([]))
xml_str = ET.tostring(fnode, encoding='unicode').encode()
with h5py.File('test.h5', 'w', driver='core', backing_store=False) as f:
g = f.create_group('surface 1')
g.create_dataset('type', data=b'implicit')
g.create_dataset('boundary_type', data=b'transmission')
g.create_dataset('coefficients',
data=np.array([0.,0.,0.,1.,0.,0.,0.,1.,0.,0.,0.,1.]))
g.create_dataset('isovalue', data=25.)
g.create_dataset('function_xml', data=xml_str)
surf = ImplicitSurface.from_hdf5(g)
assert isinstance(surf, ImplicitSurface)
assert surf.isovalue == pytest.approx(25.)
for pt in [(0.,0.,0.), (3.,4.,0.), (5.,0.,0.), (6.,0.,0.)]:
assert surf.evaluate(pt) == pytest.approx(orig.evaluate(pt), abs=1e-10)
# ==============================================================================
# TPMS
# ==============================================================================
# ── construction ──────────────────────────────────────────────────────────
@pytest.mark.parametrize("name", [
"primitive", "schwarz_p",
"gyroid", "schoen-g",
"diamond", "schwarz_d",
])
def test_valid_names_produce_implicit_surface(name):
assert isinstance(TPMS.from_pitch_isovalue(name, 1.0, 0.0), ImplicitSurface)
def test_unknown_name_raises():
"""Unknown TPMS name raises NotImplementedError."""
with pytest.raises(NotImplementedError):
TPMS.from_pitch_isovalue("screugneugneu", 1.0, 0.0)
def test_isovalue_stored_correctly():
surf = TPMS.from_pitch_isovalue("primitive", 1.0, 0.5)
assert surf.isovalue == pytest.approx(0.5)
# ── primitive (Schwartz-P) ────────────────────────────────────────────────
def test_primitive_at_origin():
"""cos(0)+cos(0)+cos(0) = 3.0, isovalue=0 → evaluate = 3."""
surf = TPMS.from_pitch_isovalue("primitive", 1.0, 0.0)
assert surf.evaluate((0., 0., 0.)) == pytest.approx(3.0, abs=1e-10)
def test_primitive_periodicity():
"""Primitive is periodic with pitch L along each axis."""
L = 1.0
surf = TPMS.from_pitch_isovalue("primitive", L, 0.0)
for pt in [(0.1, 0.2, 0.3), (0.3, 0.4, 0.1)]:
v_orig = surf.evaluate(pt)
v_shift_x = surf.evaluate((pt[0] + L, pt[1], pt[2]))
v_shift_y = surf.evaluate((pt[0], pt[1] + L, pt[2]))
v_shift_z = surf.evaluate((pt[0], pt[1], pt[2] + L))
assert v_orig == pytest.approx(v_shift_x, abs=1e-10)
assert v_orig == pytest.approx(v_shift_y, abs=1e-10)
assert v_orig == pytest.approx(v_shift_z, abs=1e-10)
def test_pitch_scales_spatial_frequency():
"""Doubling the pitch halves the spatial frequency:
evaluate_2L(x, 0, 0) == evaluate_L(x/2, 0, 0)."""
s1 = TPMS.from_pitch_isovalue("primitive", 1.0, 0.0)
s2 = TPMS.from_pitch_isovalue("primitive", 2.0, 0.0)
for x in [0.1, 0.2, 0.35]:
assert s2.evaluate((x, 0., 0.)) == pytest.approx(
s1.evaluate((x / 2., 0., 0.)), abs=1e-10)
def test_isovalue_shifts_evaluate():
"""Increasing isovalue by Δ decreases evaluate by Δ at every point."""
s0 = TPMS.from_pitch_isovalue("primitive", 1.0, 0.0)
s1 = TPMS.from_pitch_isovalue("primitive", 1.0, 1.0)
for pt in [(0., 0., 0.), (0.1, 0.2, 0.3)]:
assert s0.evaluate(pt) - s1.evaluate(pt) == pytest.approx(1.0, abs=1e-10)
# ── gyroid ────────────────────────────────────────────────────────────────
def test_gyroid_at_origin():
"""Gyroid f = sin·cos + sin·cos + sin·cos = 0 at the origin."""
surf = TPMS.from_pitch_isovalue("gyroid", 1.0, 0.0)
assert surf.evaluate((0., 0., 0.)) == pytest.approx(0., abs=1e-10)
def test_gyroid_at_quarter_pitch():
"""Gyroid at (L/4, 0, 0) should be 1.0.
Analytical:
x_scaled = 2π*(L/4)/L = π/2, y_scaled = 0, z_scaled = 0
sin(π/2)*cos(0) + sin(0)*cos(π/2) + sin(0)*cos(0) = 1 + 0 + 0 = 1.
"""
L = 1.0
surf = TPMS.from_pitch_isovalue("gyroid", L, 0.0)
assert surf.evaluate((L / 4., 0., 0.)) == pytest.approx(1.0, abs=1e-10)
def test_gyroid_periodicity():
"""Gyroid is periodic with pitch L along x."""
L = 1.0
surf = TPMS.from_pitch_isovalue("gyroid", L, 0.0)
for pt in [(0.1, 0.2, 0.3), (0.4, 0.1, 0.5)]:
assert surf.evaluate(pt) == pytest.approx(
surf.evaluate((pt[0] + L, pt[1], pt[2])), abs=1e-10)
# ── diamond —──────────────────────────────────────────────────────────────
def test_diamond_at_origin():
"""Diamond f = sin(x)*cos(y-z) + sin(y+z)*cos(x) = 0 at the origin."""
surf = TPMS.from_pitch_isovalue("diamond", 1.0, 0.0)
assert surf.evaluate((0., 0., 0.)) == pytest.approx(0., abs=1e-10)
def test_diamond_at_eighth_pitch():
"""Diamond at (L/8, 0, 0) should equal sin(π/4) = √2/2.
Analytical:
x_scaled = 2π*(L/8)/L = π/4, y_scaled = 0, z_scaled = 0
sin(π/4)*cos(0-0) + sin(0+0)*cos(π/4) = sin(π/4) + 0 = √2/2.
"""
L = 1.0
surf = TPMS.from_pitch_isovalue("diamond", L, 0.0)
expected = np.sin(np.pi / 4.) # ≈ 0.7071
assert surf.evaluate((L / 8., 0., 0.)) == pytest.approx(expected, abs=1e-10)
def test_diamond_periodicity():
"""Diamond is periodic with pitch L along x."""
L = 1.0
surf = TPMS.from_pitch_isovalue("diamond", L, 0.0)
for pt in [(0.1, 0.2, 0.3), (0.4, 0.1, 0.5)]:
assert surf.evaluate(pt) == pytest.approx(
surf.evaluate((pt[0] + L, pt[1], pt[2])), abs=1e-10)