mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Merge remote-tracking branch 'upstream/develop' into transportoperator-subclass-notransport
This commit is contained in:
yardasol
commit
b742258954
67 changed files with 2353 additions and 428 deletions
36
Dockerfile
36
Dockerfile
|
|
@ -15,22 +15,24 @@
|
|||
|
||||
# sudo docker run image_name:tag_name or ID with no tag sudo docker run ID number
|
||||
|
||||
FROM debian:bullseye-slim AS dependencies
|
||||
|
||||
# global ARG as these ARGS are used in multiple stages
|
||||
# By default one core is used to compile
|
||||
ARG compile_cores=1
|
||||
|
||||
# By default this Dockerfile builds OpenMC without DAGMC and LIBMESH support
|
||||
ARG build_dagmc=off
|
||||
ARG build_libmesh=off
|
||||
|
||||
# By default one core is used to compile
|
||||
ARG compile_cores=1
|
||||
FROM debian:bullseye-slim AS dependencies
|
||||
|
||||
ARG compile_cores
|
||||
ARG build_dagmc
|
||||
ARG build_libmesh
|
||||
|
||||
# Set default value of HOME to /root
|
||||
ENV HOME=/root
|
||||
|
||||
# OpenMC variables
|
||||
ARG openmc_branch=master
|
||||
ENV OPENMC_REPO='https://github.com/openmc-dev/openmc'
|
||||
|
||||
# Embree variables
|
||||
ENV EMBREE_TAG='v3.12.2'
|
||||
ENV EMBREE_REPO='https://github.com/embree/embree'
|
||||
|
|
@ -60,7 +62,6 @@ ENV NJOY_REPO='https://github.com/njoy/NJOY2016'
|
|||
|
||||
# Setup environment variables for Docker image
|
||||
ENV LD_LIBRARY_PATH=${DAGMC_INSTALL_DIR}/lib:$LD_LIBRARY_PATH \
|
||||
OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml \
|
||||
OPENMC_ENDF_DATA=/root/endf-b-vii.1 \
|
||||
DEBIAN_FRONTEND=noninteractive
|
||||
|
||||
|
|
@ -174,12 +175,25 @@ RUN if [ "$build_libmesh" = "on" ]; then \
|
|||
|
||||
FROM dependencies AS build
|
||||
|
||||
ENV HOME=/root
|
||||
|
||||
ARG openmc_branch=master
|
||||
ENV OPENMC_REPO='https://github.com/openmc-dev/openmc'
|
||||
|
||||
ARG compile_cores
|
||||
ARG build_dagmc
|
||||
ARG build_libmesh
|
||||
|
||||
ENV DAGMC_INSTALL_DIR=$HOME/DAGMC/
|
||||
ENV LIBMESH_INSTALL_DIR=$HOME/LIBMESH
|
||||
|
||||
# clone and install openmc
|
||||
RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
|
||||
&& git clone --shallow-submodules --recurse-submodules --single-branch -b ${openmc_branch} --depth=1 ${OPENMC_REPO} \
|
||||
&& mkdir build && cd build ; \
|
||||
if [ ${build_dagmc} = "on" ] && [ ${build_libmesh} = "on" ]; then \
|
||||
cmake ../openmc \
|
||||
-DCMAKE_CXX_COMPILER=mpicxx \
|
||||
-DOPENMC_USE_MPI=on \
|
||||
-DHDF5_PREFER_PARALLEL=on \
|
||||
-DOPENMC_USE_DAGMC=on \
|
||||
|
|
@ -188,6 +202,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
|
|||
fi ; \
|
||||
if [ ${build_dagmc} = "on" ] && [ ${build_libmesh} = "off" ]; then \
|
||||
cmake ../openmc \
|
||||
-DCMAKE_CXX_COMPILER=mpicxx \
|
||||
-DOPENMC_USE_MPI=on \
|
||||
-DHDF5_PREFER_PARALLEL=on \
|
||||
-DOPENMC_USE_DAGMC=ON \
|
||||
|
|
@ -195,6 +210,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
|
|||
fi ; \
|
||||
if [ ${build_dagmc} = "off" ] && [ ${build_libmesh} = "on" ]; then \
|
||||
cmake ../openmc \
|
||||
-DCMAKE_CXX_COMPILER=mpicxx \
|
||||
-DOPENMC_USE_MPI=on \
|
||||
-DHDF5_PREFER_PARALLEL=on \
|
||||
-DOPENMC_USE_LIBMESH=on \
|
||||
|
|
@ -202,6 +218,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
|
|||
fi ; \
|
||||
if [ ${build_dagmc} = "off" ] && [ ${build_libmesh} = "off" ]; then \
|
||||
cmake ../openmc \
|
||||
-DCMAKE_CXX_COMPILER=mpicxx \
|
||||
-DOPENMC_USE_MPI=on \
|
||||
-DHDF5_PREFER_PARALLEL=on ; \
|
||||
fi ; \
|
||||
|
|
@ -211,5 +228,8 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
|
|||
|
||||
FROM build AS release
|
||||
|
||||
ENV HOME=/root
|
||||
ENV OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml
|
||||
|
||||
# Download cross sections (NNDC and WMP) and ENDF data needed by test suite
|
||||
RUN ${HOME}/OpenMC/openmc/tools/ci/download-xs.sh
|
||||
|
|
|
|||
|
|
@ -339,6 +339,16 @@ Incoherent elastic scattering
|
|||
[eV\ :math:`^{-1}`].
|
||||
:Attributes: - **type** (*char[]*) -- 'IncoherentElastic'
|
||||
|
||||
Sum of functions
|
||||
----------------
|
||||
|
||||
:Object type: Group
|
||||
:Attributes: - **type** (*char[]*) -- "Sum"
|
||||
- **n** (*int*) -- Number of functions
|
||||
:Datasets:
|
||||
- ***func_<i>** (:ref:`function <1d_functions>`) -- Dataset for the
|
||||
i-th function (indexing starts at 1)
|
||||
|
||||
.. _angle_energy:
|
||||
|
||||
--------------------------
|
||||
|
|
@ -501,6 +511,19 @@ equiprobable bins.
|
|||
- **skewed** (*int8_t*) -- Whether discrete angles are equi-probable
|
||||
(0) or have a skewed distribution (1).
|
||||
|
||||
Mixed Elastic
|
||||
-------------
|
||||
|
||||
This angle-energy distribution is used when an evaluation specifies both
|
||||
coherent and incoherent elastic thermal neutron scattering.
|
||||
|
||||
:Object type: Group
|
||||
:Attributes: - **type** (*char[]*) -- "mixed_elastic"
|
||||
:Groups: - **coherent** -- Distribution for coherent elastic scattering. The
|
||||
format is given in :ref:`angle_energy`.
|
||||
- **incoherent** -- Distribution for incoherent elastic scattering.
|
||||
The format is given in :ref:`angle_energy`.
|
||||
|
||||
.. _energy_distribution:
|
||||
|
||||
--------------------
|
||||
|
|
|
|||
|
|
@ -79,9 +79,15 @@ The current version of the statepoint file format is 17.0.
|
|||
- **width** (*double[]*) -- Width of each mesh cell in each
|
||||
dimension.
|
||||
- **Unstructured Mesh Only:**
|
||||
- **filename** (*char[]*) -- Name of the mesh file.
|
||||
- **library** (*char[]*) -- Mesh library used to represent the
|
||||
mesh ("moab" or "libmesh").
|
||||
- **length_multiplier** (*double*) Scaling factor applied to the mesh.
|
||||
- **volumes** (*double[]*) -- Volume of each mesh cell.
|
||||
- **centroids** (*double[]*) -- Location of the mesh cell
|
||||
centroids.
|
||||
- **vertices** (*double[]*) -- x, y, z values of the mesh vertices.
|
||||
- **connectivity** (*int[]*) -- Connectivity array for the mesh
|
||||
cells.
|
||||
- **element_types** (*int[]*) -- Mesh element types.
|
||||
|
||||
**/tallies/filters/**
|
||||
|
||||
|
|
|
|||
|
|
@ -29,6 +29,7 @@ Functions
|
|||
reset
|
||||
run
|
||||
run_in_memory
|
||||
sample_external_source
|
||||
simulation_init
|
||||
simulation_finalize
|
||||
source_bank
|
||||
|
|
|
|||
|
|
@ -116,6 +116,7 @@ Angle-Energy Distributions
|
|||
IncoherentElasticAE
|
||||
IncoherentElasticAEDiscrete
|
||||
IncoherentInelasticAEDiscrete
|
||||
MixedElasticAE
|
||||
|
||||
Resonance Data
|
||||
--------------
|
||||
|
|
|
|||
|
|
@ -121,6 +121,7 @@ int openmc_regular_mesh_set_params(int32_t index, int n, const double* ll,
|
|||
int openmc_reset();
|
||||
int openmc_reset_timers();
|
||||
int openmc_run();
|
||||
int openmc_sample_external_source(size_t n, uint64_t* seed, void* sites);
|
||||
void openmc_set_seed(int64_t new_seed);
|
||||
int openmc_set_n_batches(
|
||||
int32_t n_batches, bool set_max_batches, bool add_statepoint_batch);
|
||||
|
|
|
|||
|
|
@ -24,7 +24,7 @@ using double_4dvec = vector<vector<vector<vector<double>>>>;
|
|||
constexpr int HDF5_VERSION[] {3, 0};
|
||||
|
||||
// Version numbers for binary files
|
||||
constexpr array<int, 2> VERSION_STATEPOINT {17, 0};
|
||||
constexpr array<int, 2> VERSION_STATEPOINT {18, 0};
|
||||
constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
|
||||
constexpr array<int, 2> VERSION_TRACK {3, 0};
|
||||
constexpr array<int, 2> VERSION_SUMMARY {6, 0};
|
||||
|
|
|
|||
|
|
@ -126,6 +126,26 @@ private:
|
|||
debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1]
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! Sum of multiple 1D functions
|
||||
//==============================================================================
|
||||
|
||||
class Sum1D : public Function1D {
|
||||
public:
|
||||
// Constructors
|
||||
explicit Sum1D(hid_t group);
|
||||
|
||||
//! Evaluate each function and sum results
|
||||
//! \param[in] x independent variable
|
||||
//! \return Function evaluated at x
|
||||
double operator()(double E) const override;
|
||||
|
||||
const unique_ptr<Function1D>& functions(int i) const { return functions_[i]; }
|
||||
|
||||
private:
|
||||
vector<unique_ptr<Function1D>> functions_; //!< individual functions
|
||||
};
|
||||
|
||||
//! Read 1D function from HDF5 dataset
|
||||
//! \param[in] group HDF5 group containing dataset
|
||||
//! \param[in] name Name of dataset
|
||||
|
|
|
|||
|
|
@ -61,6 +61,8 @@ void ensure_exists(hid_t obj_id, const char* name, bool attribute = false);
|
|||
vector<std::string> group_names(hid_t group_id);
|
||||
vector<hsize_t> object_shape(hid_t obj_id);
|
||||
std::string object_name(hid_t obj_id);
|
||||
hid_t open_object(hid_t group_id, const std::string& name);
|
||||
void close_object(hid_t obj_id);
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility functions
|
||||
|
|
|
|||
|
|
@ -7,13 +7,14 @@
|
|||
#include <unordered_map>
|
||||
|
||||
#include "hdf5.h"
|
||||
#include "pugixml.hpp"
|
||||
#include "xtensor/xtensor.hpp"
|
||||
#include "pugixml.hpp"
|
||||
|
||||
#include "openmc/memory.h" // for unique_ptr
|
||||
#include "openmc/particle.h"
|
||||
#include "openmc/position.h"
|
||||
#include "openmc/vector.h"
|
||||
#include "openmc/xml_interface.h"
|
||||
|
||||
#ifdef DAGMC
|
||||
#include "moab/AdaptiveKDTree.hpp"
|
||||
|
|
@ -36,6 +37,12 @@
|
|||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Constants
|
||||
//==============================================================================
|
||||
|
||||
enum class ElementType { UNSUPPORTED=-1, LINEAR_TET, LINEAR_HEX };
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
|
@ -53,7 +60,7 @@ extern vector<unique_ptr<Mesh>> meshes;
|
|||
|
||||
#ifdef LIBMESH
|
||||
namespace settings {
|
||||
// used when creating new libMesh::Mesh instances
|
||||
// used when creating new libMesh::MeshBase instances
|
||||
extern unique_ptr<libMesh::LibMeshInit> libmesh_init;
|
||||
extern const libMesh::Parallel::Communicator* libmesh_comm;
|
||||
} // namespace settings
|
||||
|
|
@ -485,6 +492,23 @@ public:
|
|||
//! \return The centroid of the bin
|
||||
virtual Position centroid(int bin) const = 0;
|
||||
|
||||
//! Get the number of vertices in the mesh
|
||||
//
|
||||
//! \return Number of vertices
|
||||
virtual int n_vertices() const = 0;
|
||||
|
||||
//! Retrieve a vertex of the mesh
|
||||
//
|
||||
//! \param[in] vertex ID
|
||||
//! \return vertex coordinates
|
||||
virtual Position vertex(int id) const = 0;
|
||||
|
||||
//! Retrieve connectivity of a mesh element
|
||||
//
|
||||
//! \param[in] element ID
|
||||
//! \return element connectivity as IDs of the vertices
|
||||
virtual std::vector<int> connectivity(int id) const = 0;
|
||||
|
||||
//! Get the volume of a mesh bin
|
||||
//
|
||||
//! \param[in] bin Bin to return the volume for
|
||||
|
|
@ -557,6 +581,12 @@ public:
|
|||
|
||||
Position centroid(int bin) const override;
|
||||
|
||||
int n_vertices() const override;
|
||||
|
||||
Position vertex(int id) const override;
|
||||
|
||||
std::vector<int> connectivity(int id) const override;
|
||||
|
||||
double volume(int bin) const override;
|
||||
|
||||
private:
|
||||
|
|
@ -621,6 +651,9 @@ private:
|
|||
//! \return MOAB EntityHandle of tet
|
||||
moab::EntityHandle get_ent_handle_from_bin(int bin) const;
|
||||
|
||||
//! Get a vertex index into the global range from a handle
|
||||
int get_vert_idx_from_handle(moab::EntityHandle vert) const;
|
||||
|
||||
//! Get the bin for a given mesh cell index
|
||||
//
|
||||
//! \param[in] idx Index of the mesh cell.
|
||||
|
|
@ -655,6 +688,7 @@ private:
|
|||
|
||||
// Data members
|
||||
moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh
|
||||
moab::Range verts_; //!< Range of vertex EntityHandle's in the mesh
|
||||
moab::EntityHandle tetset_; //!< EntitySet containing all tetrahedra
|
||||
moab::EntityHandle kdtree_root_; //!< Root of the MOAB KDTree
|
||||
std::shared_ptr<moab::Interface> mbi_; //!< MOAB instance
|
||||
|
|
@ -671,7 +705,8 @@ class LibMesh : public UnstructuredMesh {
|
|||
public:
|
||||
// Constructors
|
||||
LibMesh(pugi::xml_node node);
|
||||
LibMesh(const std::string& filename, double length_multiplier = 1.0);
|
||||
LibMesh(const std::string & filename, double length_multiplier = 1.0);
|
||||
LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier = 1.0);
|
||||
|
||||
static const std::string mesh_lib_type;
|
||||
|
||||
|
|
@ -701,10 +736,19 @@ public:
|
|||
|
||||
Position centroid(int bin) const override;
|
||||
|
||||
int n_vertices() const override;
|
||||
|
||||
Position vertex(int id) const override;
|
||||
|
||||
std::vector<int> connectivity(int id) const override;
|
||||
|
||||
double volume(int bin) const override;
|
||||
|
||||
libMesh::MeshBase* mesh_ptr() const { return m_; };
|
||||
|
||||
private:
|
||||
void initialize() override;
|
||||
void set_mesh_pointer_from_filename(const std::string& filename);
|
||||
|
||||
// Methods
|
||||
|
||||
|
|
@ -715,7 +759,8 @@ private:
|
|||
int get_bin_from_element(const libMesh::Elem* elem) const;
|
||||
|
||||
// Data members
|
||||
unique_ptr<libMesh::Mesh> m_; //!< pointer to the libMesh mesh instance
|
||||
unique_ptr<libMesh::MeshBase> unique_m_ = nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC
|
||||
libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set during intialization
|
||||
vector<unique_ptr<libMesh::PointLocatorBase>>
|
||||
pl_; //!< per-thread point locators
|
||||
unique_ptr<libMesh::EquationSystems>
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@
|
|||
#include <iostream>
|
||||
#include <stdexcept> // for out_of_range
|
||||
|
||||
#include "fmt/format.h"
|
||||
#include "openmc/array.h"
|
||||
#include "openmc/vector.h"
|
||||
|
||||
|
|
@ -210,4 +211,14 @@ using Direction = Position;
|
|||
|
||||
} // namespace openmc
|
||||
|
||||
template<>
|
||||
struct fmt::formatter<openmc::Position> : formatter<std::string> {
|
||||
template<typename FormatContext>
|
||||
auto format(const openmc::Position& pos, FormatContext& ctx)
|
||||
{
|
||||
return fmt::formatter<std::string>::format(
|
||||
fmt::format("({}, {}, {})", pos.x, pos.y, pos.z), ctx);
|
||||
}
|
||||
};
|
||||
|
||||
#endif // OPENMC_POSITION_H
|
||||
|
|
|
|||
|
|
@ -150,6 +150,34 @@ private:
|
|||
//!< each incident energy
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! Mixed coherent/incoherent elastic angle-energy distribution
|
||||
//==============================================================================
|
||||
|
||||
class MixedElasticAE : public AngleEnergy {
|
||||
public:
|
||||
//! Construct from HDF5 file
|
||||
//
|
||||
//! \param[in] group HDF5 group
|
||||
explicit MixedElasticAE(
|
||||
hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs);
|
||||
|
||||
//! Sample distribution for an angle and energy
|
||||
//! \param[in] E_in Incoming energy in [eV]
|
||||
//! \param[out] E_out Outgoing energy in [eV]
|
||||
//! \param[out] mu Outgoing cosine with respect to current direction
|
||||
//! \param[inout] seed Pseudorandom number seed pointer
|
||||
void sample(
|
||||
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
|
||||
|
||||
private:
|
||||
CoherentElasticAE coherent_dist_; //!< Coherent distribution
|
||||
unique_ptr<AngleEnergy> incoherent_dist_; //!< Incoherent distribution
|
||||
|
||||
const CoherentElasticXS& coherent_xs_; //!< Ref. to coherent XS
|
||||
const Function1D& incoherent_xs_; //!< Polymorphic ref. to incoherent XS
|
||||
};
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_SECONDARY_THERMAL_H
|
||||
|
|
|
|||
|
|
@ -44,6 +44,8 @@ class AngleEnergy(EqualityMixin, ABC):
|
|||
return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group)
|
||||
elif dist_type == 'incoherent_inelastic':
|
||||
return openmc.data.IncoherentInelasticAE.from_hdf5(group)
|
||||
elif dist_type == 'mixed_elastic':
|
||||
return openmc.data.MixedElasticAE.from_hdf5(group)
|
||||
|
||||
@staticmethod
|
||||
def from_ace(ace, location_dist, location_start, rx=None):
|
||||
|
|
|
|||
|
|
@ -544,7 +544,7 @@ class Combination(EqualityMixin):
|
|||
self._operations = operations
|
||||
|
||||
|
||||
class Sum(EqualityMixin):
|
||||
class Sum(Function1D):
|
||||
"""Sum of multiple functions.
|
||||
|
||||
This class allows you to create a callable object which represents the sum
|
||||
|
|
@ -578,6 +578,49 @@ class Sum(EqualityMixin):
|
|||
cv.check_type('functions', functions, Iterable, Callable)
|
||||
self._functions = functions
|
||||
|
||||
def to_hdf5(self, group, name='xy'):
|
||||
"""Write sum of functions to an HDF5 group
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to write to
|
||||
name : str
|
||||
Name of the dataset to create
|
||||
|
||||
"""
|
||||
sum_group = group.create_group(name)
|
||||
sum_group.attrs['type'] = np.string_(type(self).__name__)
|
||||
sum_group.attrs['n'] = len(self.functions)
|
||||
for i, f in enumerate(self.functions):
|
||||
f.to_hdf5(sum_group, f'func_{i+1}')
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, group):
|
||||
"""Generate sum of functions from an HDF5 group
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
Group to read from
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.data.Sum
|
||||
Functions read from the group
|
||||
|
||||
"""
|
||||
n = group.attrs['n']
|
||||
functions = [
|
||||
Function1D.from_hdf5(group[f'func_{i+1}'])
|
||||
for i in range(n)
|
||||
]
|
||||
return cls(functions)
|
||||
|
||||
|
||||
class Regions1D(EqualityMixin):
|
||||
r"""Piecewise composition of multiple functions.
|
||||
|
|
|
|||
|
|
@ -299,8 +299,8 @@ def reconstruct_slbw(slbw, double E):
|
|||
# Determine shift and penetration at modified energy
|
||||
if slbw._competitive[i]:
|
||||
Ex = E + slbw.q_value[l]*(A + 1)/A
|
||||
rhoc = slbw.channel_radius[l](Ex)
|
||||
rhochat = slbw.scattering_radius[l](Ex)
|
||||
rhoc = k*slbw.channel_radius[l](Ex)
|
||||
rhochat = k*slbw.scattering_radius[l](Ex)
|
||||
P_c, S_c = penetration_shift(l, rhoc)
|
||||
if Ex < 0:
|
||||
P_c = 0
|
||||
|
|
|
|||
|
|
@ -19,12 +19,12 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf
|
|||
from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, isotopes
|
||||
from .ace import Table, get_table, Library
|
||||
from .angle_energy import AngleEnergy
|
||||
from .function import Tabulated1D, Function1D
|
||||
from .function import Tabulated1D, Function1D, Sum
|
||||
from .njoy import make_ace_thermal
|
||||
from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE,
|
||||
IncoherentElasticAEDiscrete,
|
||||
IncoherentInelasticAEDiscrete,
|
||||
IncoherentInelasticAE)
|
||||
IncoherentInelasticAE, MixedElasticAE)
|
||||
|
||||
|
||||
_THERMAL_NAMES = {
|
||||
|
|
@ -694,29 +694,53 @@ class ThermalScattering(EqualityMixin):
|
|||
|
||||
# Incoherent/coherent elastic scattering cross section
|
||||
idx = ace.jxs[4]
|
||||
n_mu = ace.nxs[6] + 1
|
||||
if idx != 0:
|
||||
n_energy = int(ace.xss[idx])
|
||||
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
|
||||
P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
|
||||
|
||||
if ace.nxs[5] == 4:
|
||||
if ace.nxs[5] in (4, 5):
|
||||
# Coherent elastic
|
||||
xs = CoherentElastic(energy, P*EV_PER_MEV)
|
||||
distribution = CoherentElasticAE(xs)
|
||||
n_energy = int(ace.xss[idx])
|
||||
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
|
||||
P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
|
||||
coherent_xs = CoherentElastic(energy, P*EV_PER_MEV)
|
||||
coherent_dist = CoherentElasticAE(coherent_xs)
|
||||
|
||||
# Coherent elastic shouldn't have angular distributions listed
|
||||
n_mu = ace.nxs[6] + 1
|
||||
assert n_mu == 0
|
||||
else:
|
||||
# Incoherent elastic
|
||||
xs = Tabulated1D(energy, P)
|
||||
|
||||
if ace.nxs[5] in (3, 5):
|
||||
# Incoherent elastic scattering -- first determine if both
|
||||
# incoherent and coherent are present (mixed)
|
||||
mixed = (ace.nxs[5] == 5)
|
||||
|
||||
# Get cross section values
|
||||
idx = ace.jxs[7] if mixed else ace.jxs[4]
|
||||
n_energy = int(ace.xss[idx])
|
||||
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
|
||||
values = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
|
||||
|
||||
incoherent_xs = Tabulated1D(energy, values)
|
||||
|
||||
# Angular distribution
|
||||
n_mu = (ace.nxs[8] if mixed else ace.nxs[6]) + 1
|
||||
assert n_mu > 0
|
||||
idx = ace.jxs[6]
|
||||
idx = ace.jxs[9] if mixed else ace.jxs[6]
|
||||
mu_out = ace.xss[idx:idx + n_energy * n_mu]
|
||||
mu_out.shape = (n_energy, n_mu)
|
||||
distribution = IncoherentElasticAEDiscrete(mu_out)
|
||||
incoherent_dist = IncoherentElasticAEDiscrete(mu_out)
|
||||
|
||||
if ace.nxs[5] == 3:
|
||||
xs = incoherent_xs
|
||||
dist = incoherent_dist
|
||||
elif ace.nxs[5] == 4:
|
||||
xs = coherent_xs
|
||||
dist = coherent_dist
|
||||
else:
|
||||
# Create mixed cross section -- note that coherent must come
|
||||
# first due to assumption on C++ side
|
||||
xs = Sum([coherent_xs, incoherent_xs])
|
||||
|
||||
# Create mixed distribution
|
||||
distribution = MixedElasticAE(coherent_dist, incoherent_dist)
|
||||
|
||||
table.elastic = ThermalScatteringReaction({T: xs}, {T: distribution})
|
||||
|
||||
|
|
@ -802,7 +826,7 @@ class ThermalScattering(EqualityMixin):
|
|||
# Replace ACE data with ENDF data
|
||||
rx, rx_endf = data.elastic, data_endf.elastic
|
||||
for t in temperatures:
|
||||
if isinstance(rx_endf.xs[t], IncoherentElastic):
|
||||
if isinstance(rx_endf.xs[t], (IncoherentElastic, Sum)):
|
||||
rx.xs[t] = rx_endf.xs[t]
|
||||
rx.distribution[t] = rx_endf.distribution[t]
|
||||
|
||||
|
|
@ -832,20 +856,14 @@ class ThermalScattering(EqualityMixin):
|
|||
# Read coherent/incoherent elastic data
|
||||
elastic = None
|
||||
if (7, 2) in ev.section:
|
||||
xs = {}
|
||||
distribution = {}
|
||||
|
||||
file_obj = StringIO(ev.section[7, 2])
|
||||
lhtr = endf.get_head_record(file_obj)[2]
|
||||
if lhtr == 1:
|
||||
# coherent elastic
|
||||
|
||||
# Define helper functions to avoid duplication
|
||||
def get_coherent_elastic(file_obj):
|
||||
# Get structure factor at first temperature
|
||||
params, S = endf.get_tab1_record(file_obj)
|
||||
strT = _temperature_str(params[0])
|
||||
n_temps = params[2]
|
||||
bragg_edges = S.x
|
||||
xs[strT] = CoherentElastic(bragg_edges, S.y)
|
||||
xs = {strT: CoherentElastic(bragg_edges, S.y)}
|
||||
distribution = {strT: CoherentElasticAE(xs[strT])}
|
||||
|
||||
# Get structure factor for subsequent temperatures
|
||||
|
|
@ -854,15 +872,34 @@ class ThermalScattering(EqualityMixin):
|
|||
strT = _temperature_str(params[0])
|
||||
xs[strT] = CoherentElastic(bragg_edges, S)
|
||||
distribution[strT] = CoherentElasticAE(xs[strT])
|
||||
return xs, distribution
|
||||
|
||||
elif lhtr == 2:
|
||||
# incoherent elastic
|
||||
def get_incoherent_elastic(file_obj):
|
||||
params, W = endf.get_tab1_record(file_obj)
|
||||
bound_xs = params[0]
|
||||
xs = {}
|
||||
distribution = {}
|
||||
for T, debye_waller in zip(W.x, W.y):
|
||||
strT = _temperature_str(T)
|
||||
xs[strT] = IncoherentElastic(bound_xs, debye_waller)
|
||||
distribution[strT] = IncoherentElasticAE(debye_waller)
|
||||
return xs, distribution
|
||||
|
||||
file_obj = StringIO(ev.section[7, 2])
|
||||
lhtr = endf.get_head_record(file_obj)[2]
|
||||
if lhtr == 1:
|
||||
# coherent elastic
|
||||
xs, distribution = get_coherent_elastic(file_obj)
|
||||
elif lhtr == 2:
|
||||
# incoherent elastic
|
||||
xs, distribution = get_incoherent_elastic(file_obj)
|
||||
elif lhtr == 3:
|
||||
# mixed coherent / incoherent elastic
|
||||
xs_c, dist_c = get_coherent_elastic(file_obj)
|
||||
xs_i, dist_i = get_incoherent_elastic(file_obj)
|
||||
assert sorted(xs_c) == sorted(xs_i)
|
||||
xs = {T: Sum([xs_c[T], xs_i[T]]) for T in xs_c}
|
||||
distribution = {T: MixedElasticAE(dist_c[T], dist_i[T]) for T in dist_c}
|
||||
|
||||
elastic = ThermalScatteringReaction(xs, distribution)
|
||||
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@ import numpy as np
|
|||
|
||||
from .angle_energy import AngleEnergy
|
||||
from .correlated import CorrelatedAngleEnergy
|
||||
import openmc.data
|
||||
|
||||
|
||||
class CoherentElasticAE(AngleEnergy):
|
||||
|
|
@ -43,7 +44,27 @@ class CoherentElasticAE(AngleEnergy):
|
|||
|
||||
"""
|
||||
group.attrs['type'] = np.string_('coherent_elastic')
|
||||
group['coherent_xs'] = group.parent['xs']
|
||||
self.coherent_xs.to_hdf5(group, 'coherent_xs')
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, group):
|
||||
"""Generate coherent elastic distribution from HDF5 data
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to read from
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.data.CoherentElasticAE
|
||||
Coherent elastic distribution
|
||||
|
||||
"""
|
||||
coherent_xs = openmc.data.CoherentElastic.from_hdf5(group['coherent_xs'])
|
||||
return cls(coherent_xs)
|
||||
|
||||
|
||||
class IncoherentElasticAE(AngleEnergy):
|
||||
|
|
@ -101,7 +122,7 @@ class IncoherentElasticAE(AngleEnergy):
|
|||
Incoherent elastic distribution
|
||||
|
||||
"""
|
||||
return cls(group['debye_waller'])
|
||||
return cls(group['debye_waller'][()])
|
||||
|
||||
|
||||
class IncoherentElasticAEDiscrete(AngleEnergy):
|
||||
|
|
@ -210,3 +231,62 @@ class IncoherentInelasticAEDiscrete(AngleEnergy):
|
|||
|
||||
class IncoherentInelasticAE(CorrelatedAngleEnergy):
|
||||
_name = 'incoherent_inelastic'
|
||||
|
||||
|
||||
class MixedElasticAE(AngleEnergy):
|
||||
"""Secondary distribution for mixed coherent/incoherent thermal elastic
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
coherent : AngleEnergy
|
||||
Secondary distribution for coherent elastic scattering
|
||||
incoherent : AngleEnergy
|
||||
Secondary distribution for incoherent elastic scattering
|
||||
|
||||
Attributes
|
||||
----------
|
||||
coherent : AngleEnergy
|
||||
Secondary distribution for coherent elastic scattering
|
||||
incoherent : AngleEnergy
|
||||
Secondary distribution for incoherent elastic scattering
|
||||
|
||||
"""
|
||||
def __init__(self, coherent, incoherent):
|
||||
self.coherent = coherent
|
||||
self.incoherent = incoherent
|
||||
|
||||
def to_hdf5(self, group):
|
||||
"""Write mixed elastic distribution to an HDF5 group
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to write to
|
||||
|
||||
"""
|
||||
group.attrs['type'] = np.string_('mixed_elastic')
|
||||
coherent_group = group.create_group('coherent')
|
||||
self.coherent.to_hdf5(coherent_group)
|
||||
incoherent_group = group.create_group('incoherent')
|
||||
self.incoherent.to_hdf5(incoherent_group)
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, group):
|
||||
"""Generate mixed thermal elastic distribution from HDF5 data
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to read from
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.data.MixedElasticAE
|
||||
Mixed thermal elastic distribution
|
||||
|
||||
"""
|
||||
coherent = AngleEnergy.from_hdf5(group['coherent'])
|
||||
incoherent = AngleEnergy.from_hdf5(group['incoherent'])
|
||||
return cls(coherent, incoherent)
|
||||
|
|
|
|||
|
|
@ -1,8 +1,10 @@
|
|||
from contextlib import contextmanager
|
||||
from ctypes import (c_bool, c_int, c_int32, c_int64, c_double, c_char_p,
|
||||
c_char, POINTER, Structure, c_void_p, create_string_buffer)
|
||||
c_char, POINTER, Structure, c_void_p, create_string_buffer,
|
||||
c_uint64, c_size_t)
|
||||
import sys
|
||||
import os
|
||||
from random import getrandbits
|
||||
|
||||
import numpy as np
|
||||
from numpy.ctypeslib import as_array
|
||||
|
|
@ -10,6 +12,7 @@ from numpy.ctypeslib import as_array
|
|||
from . import _dll
|
||||
from .error import _error_handler
|
||||
import openmc.lib
|
||||
import openmc
|
||||
|
||||
|
||||
class _SourceSite(Structure):
|
||||
|
|
@ -95,6 +98,9 @@ _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double),
|
|||
POINTER(c_double)]
|
||||
_dll.openmc_global_bounding_box.restype = c_int
|
||||
_dll.openmc_global_bounding_box.errcheck = _error_handler
|
||||
_dll.openmc_sample_external_source.argtypes = [c_size_t, POINTER(c_uint64), POINTER(_SourceSite)]
|
||||
_dll.openmc_sample_external_source.restype = c_int
|
||||
_dll.openmc_sample_external_source.errcheck = _error_handler
|
||||
|
||||
|
||||
def global_bounding_box():
|
||||
|
|
@ -415,6 +421,45 @@ def run(output=True):
|
|||
_dll.openmc_run()
|
||||
|
||||
|
||||
def sample_external_source(n_samples=1, prn_seed=None):
|
||||
"""Sample external source
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n_samples : int
|
||||
Number of samples
|
||||
prn_seed : int
|
||||
Pseudorandom number generator (PRNG) seed; if None, one will be
|
||||
generated randomly.
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of openmc.SourceParticle
|
||||
List of samples source particles
|
||||
|
||||
"""
|
||||
if n_samples <= 0:
|
||||
raise ValueError("Number of samples must be positive")
|
||||
if prn_seed is None:
|
||||
prn_seed = getrandbits(63)
|
||||
|
||||
# Call into C API to sample source
|
||||
sites_array = (_SourceSite * n_samples)()
|
||||
_dll.openmc_sample_external_source(c_size_t(n_samples), c_uint64(prn_seed), sites_array)
|
||||
|
||||
# Convert to list of SourceParticle and return
|
||||
return [
|
||||
openmc.SourceParticle(
|
||||
r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt,
|
||||
delayed_group=site.delayed_group, surf_id=site.surf_id,
|
||||
particle=openmc.ParticleType(site.particle)
|
||||
)
|
||||
for site in sites_array
|
||||
]
|
||||
|
||||
|
||||
def simulation_init():
|
||||
"""Initialize simulation"""
|
||||
_dll.openmc_simulation_init()
|
||||
|
|
|
|||
|
|
@ -1,12 +1,11 @@
|
|||
from ctypes import c_int, c_double, POINTER, c_uint64
|
||||
from random import getrandbits
|
||||
|
||||
import numpy as np
|
||||
from numpy.ctypeslib import ndpointer
|
||||
|
||||
from . import _dll
|
||||
|
||||
from random import getrandbits
|
||||
|
||||
|
||||
_dll.t_percentile.restype = c_double
|
||||
_dll.t_percentile.argtypes = [c_double, c_int]
|
||||
|
|
@ -200,7 +199,8 @@ def rotate_angle(uvw0, mu, phi, prn_seed=None):
|
|||
phi : float
|
||||
Azimuthal angle; if None, one will be sampled uniformly
|
||||
prn_seed : int
|
||||
PRNG seed; if None, one will be generated randomly
|
||||
Pseudorandom number generator (PRNG) seed; if None, one will be
|
||||
generated randomly.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -232,7 +232,8 @@ def maxwell_spectrum(T, prn_seed=None):
|
|||
T : float
|
||||
Spectrum parameter
|
||||
prn_seed : int
|
||||
PRNG seed; if None, one will be generated randomly
|
||||
Pseudorandom number generator (PRNG) seed; if None, one will be
|
||||
generated randomly.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -240,10 +241,10 @@ def maxwell_spectrum(T, prn_seed=None):
|
|||
Sampled outgoing energy
|
||||
|
||||
"""
|
||||
|
||||
|
||||
if prn_seed is None:
|
||||
prn_seed = getrandbits(63)
|
||||
|
||||
|
||||
return _dll.maxwell_spectrum(T, c_uint64(prn_seed))
|
||||
|
||||
|
||||
|
|
@ -257,7 +258,8 @@ def watt_spectrum(a, b, prn_seed=None):
|
|||
b : float
|
||||
Spectrum parameter b
|
||||
prn_seed : int
|
||||
PRNG seed; if None, one will be generated randomly
|
||||
Pseudorandom number generator (PRNG) seed; if None, one will be
|
||||
generated randomly.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -265,7 +267,7 @@ def watt_spectrum(a, b, prn_seed=None):
|
|||
Sampled outgoing energy
|
||||
|
||||
"""
|
||||
|
||||
|
||||
if prn_seed is None:
|
||||
prn_seed = getrandbits(63)
|
||||
|
||||
|
|
@ -282,7 +284,8 @@ def normal_variate(mean_value, std_dev, prn_seed=None):
|
|||
std_dev : float
|
||||
Standard deviation of the normal distribution
|
||||
prn_seed : int
|
||||
PRNG seed; if None, one will be generated randomly
|
||||
Pseudorandom number generator (PRNG) seed; if None, one will be
|
||||
generated randomly.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -290,7 +293,7 @@ def normal_variate(mean_value, std_dev, prn_seed=None):
|
|||
Sampled outgoing normally distributed value
|
||||
|
||||
"""
|
||||
|
||||
|
||||
if prn_seed is None:
|
||||
prn_seed = getrandbits(63)
|
||||
|
||||
|
|
|
|||
|
|
@ -34,7 +34,9 @@ class Material(IDManagerMixin):
|
|||
To create a material, one should create an instance of this class, add
|
||||
nuclides or elements with :meth:`Material.add_nuclide` or
|
||||
:meth:`Material.add_element`, respectively, and set the total material
|
||||
density with :meth:`Material.set_density()`. The material can then be
|
||||
density with :meth:`Material.set_density()`. Alternatively, you can
|
||||
use :meth:`Material.add_components()` to pass a dictionary
|
||||
containing all the component information. The material can then be
|
||||
assigned to a cell using the :attr:`Cell.fill` attribute.
|
||||
|
||||
Parameters
|
||||
|
|
@ -403,6 +405,54 @@ class Material(IDManagerMixin):
|
|||
|
||||
self._nuclides.append(NuclideTuple(nuclide, percent, percent_type))
|
||||
|
||||
def add_components(self, components: dict, percent_type: str = 'ao'):
|
||||
""" Add multiple elements or nuclides to a material
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
components : dict of str to float or dict
|
||||
Dictionary mapping element or nuclide names to their atom or weight
|
||||
percent. To specify enrichment of an element, the entry of
|
||||
``components`` for that element must instead be a dictionary
|
||||
containing the keyword arguments as well as a value for
|
||||
``'percent'``
|
||||
percent_type : {'ao', 'wo'}
|
||||
'ao' for atom percent and 'wo' for weight percent
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> mat = openmc.Material()
|
||||
>>> components = {'Li': {'percent': 1.0,
|
||||
>>> 'enrichment': 60.0,
|
||||
>>> 'enrichment_target': 'Li7'},
|
||||
>>> 'Fl': 1.0,
|
||||
>>> 'Be6': 0.5}
|
||||
>>> mat.add_components(components)
|
||||
|
||||
"""
|
||||
|
||||
for component, params in components.items():
|
||||
cv.check_type('component', component, str)
|
||||
if isinstance(params, float):
|
||||
params = {'percent': params}
|
||||
|
||||
else:
|
||||
cv.check_type('params', params, dict)
|
||||
if 'percent' not in params:
|
||||
raise ValueError("An entry in the dictionary does not have "
|
||||
"a required key: 'percent'")
|
||||
|
||||
params['percent_type'] = percent_type
|
||||
|
||||
## check if nuclide
|
||||
if str.isdigit(component[-1]):
|
||||
self.add_nuclide(component, **params)
|
||||
else: # is element
|
||||
kwargs = params
|
||||
self.add_element(component, **params)
|
||||
|
||||
def remove_nuclide(self, nuclide: str):
|
||||
"""Remove a nuclide from the material
|
||||
|
||||
|
|
|
|||
261
openmc/mesh.py
261
openmc/mesh.py
|
|
@ -2,6 +2,7 @@ from abc import ABC, abstractmethod
|
|||
from collections.abc import Iterable
|
||||
from math import pi
|
||||
from numbers import Real, Integral
|
||||
from pathlib import Path
|
||||
import warnings
|
||||
from xml.etree import ElementTree as ET
|
||||
|
||||
|
|
@ -248,8 +249,14 @@ class StructuredMesh(MeshBase):
|
|||
|
||||
# create VTK arrays for each of
|
||||
# the data sets
|
||||
|
||||
# maintain a list of the datasets as added
|
||||
# to the VTK arrays to ensure they persist
|
||||
# in memory until the file is written
|
||||
datasets_out = []
|
||||
for label, dataset in datasets.items():
|
||||
dataset = np.asarray(dataset).flatten()
|
||||
datasets_out.append(dataset)
|
||||
|
||||
if volume_normalization:
|
||||
dataset /= self.volumes.flatten()
|
||||
|
|
@ -692,7 +699,7 @@ class RegularMesh(StructuredMesh):
|
|||
root_cell.fill = lattice
|
||||
|
||||
return root_cell, cells
|
||||
|
||||
|
||||
def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
|
||||
"""Creates a VTK object of the mesh
|
||||
|
||||
|
|
@ -1434,7 +1441,7 @@ class UnstructuredMesh(MeshBase):
|
|||
|
||||
Parameters
|
||||
----------
|
||||
filename : str
|
||||
filename : str or pathlib.Path
|
||||
Location of the unstructured mesh file
|
||||
library : {'moab', 'libmesh'}
|
||||
Mesh library used for the unstructured mesh tally
|
||||
|
|
@ -1462,20 +1469,39 @@ class UnstructuredMesh(MeshBase):
|
|||
be generated for this mesh
|
||||
volumes : Iterable of float
|
||||
Volumes of the unstructured mesh elements
|
||||
centroids : numpy.ndarray
|
||||
Centroids of the mesh elements with array shape (n_elements, 3)
|
||||
|
||||
vertices : numpy.ndarray
|
||||
Coordinates of the mesh vertices with array shape (n_elements, 3)
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
connectivity : numpy.ndarray
|
||||
Connectivity of the elements with array shape (n_elements, 8)
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
element_types : Iterable of integers
|
||||
Mesh element types
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
total_volume : float
|
||||
Volume of the unstructured mesh in total
|
||||
centroids : Iterable of tuple
|
||||
An iterable of element centroid coordinates, e.g. [(0.0, 0.0, 0.0),
|
||||
(1.0, 1.0, 1.0), ...]
|
||||
"""
|
||||
|
||||
_UNSUPPORTED_ELEM = -1
|
||||
_LINEAR_TET = 0
|
||||
_LINEAR_HEX = 1
|
||||
|
||||
def __init__(self, filename, library, mesh_id=None, name='',
|
||||
length_multiplier=1.0):
|
||||
super().__init__(mesh_id, name)
|
||||
self.filename = filename
|
||||
self._volumes = None
|
||||
self._centroids = None
|
||||
self._n_elements = None
|
||||
self._conectivity = None
|
||||
self._vertices = None
|
||||
self.library = library
|
||||
self._output = True
|
||||
self._output = False
|
||||
self.length_multiplier = length_multiplier
|
||||
|
||||
@property
|
||||
|
|
@ -1484,7 +1510,7 @@ class UnstructuredMesh(MeshBase):
|
|||
|
||||
@filename.setter
|
||||
def filename(self, filename):
|
||||
cv.check_type('Unstructured Mesh filename', filename, str)
|
||||
cv.check_type('Unstructured Mesh filename', filename, (str, Path))
|
||||
self._filename = filename
|
||||
|
||||
@property
|
||||
|
|
@ -1540,19 +1566,33 @@ class UnstructuredMesh(MeshBase):
|
|||
def centroids(self):
|
||||
return self._centroids
|
||||
|
||||
@property
|
||||
def vertices(self):
|
||||
return self._vertices
|
||||
|
||||
@property
|
||||
def connectivity(self):
|
||||
return self._connectivity
|
||||
|
||||
@property
|
||||
def element_types(self):
|
||||
return self._element_types
|
||||
|
||||
@property
|
||||
def centroids(self):
|
||||
return np.array([self.centroid(i) for i in range(self.n_elements)])
|
||||
|
||||
@property
|
||||
def n_elements(self):
|
||||
if self._centroids is None:
|
||||
if self._n_elements is None:
|
||||
raise RuntimeError("No information about this mesh has "
|
||||
"been loaded from a statepoint file.")
|
||||
return len(self._centroids)
|
||||
return self._n_elements
|
||||
|
||||
|
||||
@centroids.setter
|
||||
def centroids(self, centroids):
|
||||
cv.check_type("Unstructured mesh centroids", centroids,
|
||||
Iterable, Real)
|
||||
self._centroids = centroids
|
||||
@n_elements.setter
|
||||
def n_elements(self, val):
|
||||
cv.check_type('Number of elements', val, Integral)
|
||||
self._n_elements = val
|
||||
|
||||
@property
|
||||
def length_multiplier(self):
|
||||
|
|
@ -1567,29 +1607,44 @@ class UnstructuredMesh(MeshBase):
|
|||
|
||||
@property
|
||||
def dimension(self):
|
||||
return self.n_elements
|
||||
return (self.n_elements,)
|
||||
|
||||
@property
|
||||
def n_dimension(self):
|
||||
return 3
|
||||
|
||||
@property
|
||||
def vertices(self):
|
||||
raise NotImplementedError("Vertices for UnstructuredMesh objects are "
|
||||
"not yet available")
|
||||
|
||||
def __repr__(self):
|
||||
string = super().__repr__()
|
||||
string += '{: <16}=\t{}\n'.format('\tFilename', self.filename)
|
||||
string += '{: <16}=\t{}\n'.format('\tMesh Library', self.mesh_lib)
|
||||
string += '{: <16}=\t{}\n'.format('\tMesh Library', self.library)
|
||||
if self.length_multiplier != 1.0:
|
||||
string += '{: <16}=\t{}\n'.format('\tLength multiplier',
|
||||
self.length_multiplier)
|
||||
return string
|
||||
|
||||
def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
|
||||
"""Map data to the unstructured mesh element centroids
|
||||
to create a VTK point-cloud dataset.
|
||||
def centroid(self, bin):
|
||||
"""Return the vertex averaged centroid of an element
|
||||
|
||||
Parameters
|
||||
----------
|
||||
bin : int
|
||||
Bin ID for the returned centroid
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.ndarray
|
||||
x, y, z values of the element centroid
|
||||
|
||||
"""
|
||||
conn = self.connectivity[bin]
|
||||
coords = self.vertices[conn]
|
||||
return coords.mean(axis=0)
|
||||
|
||||
def write_vtk_mesh(self, **kwargs):
|
||||
"""Map data to unstructured VTK mesh elements.
|
||||
|
||||
.. deprecated:: 0.13
|
||||
Use :func:`UnstructuredMesh.write_data_to_vtk` instead.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
|
@ -1601,83 +1656,103 @@ class UnstructuredMesh(MeshBase):
|
|||
volume_normalization : bool
|
||||
Whether or not to normalize the data by the
|
||||
volume of the mesh elements
|
||||
|
||||
Raises
|
||||
------
|
||||
RuntimeError
|
||||
when the size of a dataset doesn't match the number of cells
|
||||
"""
|
||||
warnings.warn(
|
||||
"The 'UnstructuredMesh.write_vtk_mesh' method has been renamed "
|
||||
"to 'write_data_to_vtk' and will be removed in a future version "
|
||||
" of OpenMC.", FutureWarning
|
||||
)
|
||||
self.write_data_to_vtk(**kwargs)
|
||||
|
||||
def write_data_to_vtk(self, filename=None, datasets=None, volume_normalization=True):
|
||||
"""Map data to unstructured VTK mesh elements.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str or pathlib.Path
|
||||
Name of the VTK file to write
|
||||
datasets : dict
|
||||
Dictionary whose keys are the data labels
|
||||
and values are numpy appropriately sized arrays
|
||||
of the data
|
||||
volume_normalization : bool
|
||||
Whether or not to normalize the data by the
|
||||
volume of the mesh elements
|
||||
"""
|
||||
import vtk
|
||||
from vtk.util import numpy_support as vtk_npsup
|
||||
from vtk.util import numpy_support as nps
|
||||
|
||||
if self.centroids is None:
|
||||
raise RuntimeError("No centroid information is present on this "
|
||||
"unstructured mesh. Please load this "
|
||||
"information from a relevant statepoint file.")
|
||||
if self.connectivity is None or self.vertices is None:
|
||||
raise RuntimeError('This mesh has not been '
|
||||
'loaded from a statepoint file.')
|
||||
|
||||
if self.volumes is None and volume_normalization:
|
||||
raise RuntimeError("No volume data is present on this "
|
||||
"unstructured mesh. Please load the "
|
||||
" mesh information from a statepoint file.")
|
||||
if filename is None:
|
||||
filename = f'mesh_{self.id}.vtk'
|
||||
|
||||
# check that the data sets are appropriately sized
|
||||
errmsg = "The size of the dataset {} should be equal to the number of cells"
|
||||
for label, dataset in datasets.items():
|
||||
if isinstance(dataset, np.ndarray):
|
||||
if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]:
|
||||
raise RuntimeError(errmsg.format(label))
|
||||
writer = vtk.vtkUnstructuredGridWriter()
|
||||
|
||||
writer.SetFileName(str(filename))
|
||||
|
||||
grid = vtk.vtkUnstructuredGrid()
|
||||
|
||||
vtk_pnts = vtk.vtkPoints()
|
||||
vtk_pnts.SetData(nps.numpy_to_vtk(self.vertices))
|
||||
grid.SetPoints(vtk_pnts)
|
||||
|
||||
n_skipped = 0
|
||||
elems = []
|
||||
for elem_type, conn in zip(self.element_types, self.connectivity):
|
||||
if elem_type == self._LINEAR_TET:
|
||||
elem = vtk.vtkTetra()
|
||||
elif elem_type == self._LINEAR_HEX:
|
||||
elem = vtk.vtkHexahedron()
|
||||
elif elem_type == self._UNSUPPORTED_ELEM:
|
||||
n_skipped += 1
|
||||
else:
|
||||
if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]:
|
||||
raise RuntimeError(errmsg.format(label))
|
||||
cv.check_type('label', label, str)
|
||||
raise RuntimeError(f'Invalid element type {elem_type} found')
|
||||
for i, c in enumerate(conn):
|
||||
if c == -1:
|
||||
break
|
||||
elem.GetPointIds().SetId(i, c)
|
||||
elems.append(elem)
|
||||
|
||||
# create data arrays for the cells/points
|
||||
cell_dim = 1
|
||||
vertices = vtk.vtkCellArray()
|
||||
points = vtk.vtkPoints()
|
||||
if n_skipped > 0:
|
||||
warnings.warn(f'{n_skipped} elements were not written because '
|
||||
'they are not of type linear tet/hex')
|
||||
|
||||
for centroid in self.centroids:
|
||||
# create a point for each centroid
|
||||
point_id = points.InsertNextPoint(centroid * self.length_multiplier)
|
||||
# create a cell of type "Vertex" for each point
|
||||
cell_id = vertices.InsertNextCell(cell_dim, (point_id,))
|
||||
for elem in elems:
|
||||
grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())
|
||||
|
||||
# create a VTK data object
|
||||
poly_data = vtk.vtkPolyData()
|
||||
poly_data.SetPoints(points)
|
||||
poly_data.SetVerts(vertices)
|
||||
|
||||
# strange VTK nuance:
|
||||
# data must be held in some container
|
||||
# until the vtk file is written
|
||||
data_holder = []
|
||||
|
||||
# create VTK arrays for each of
|
||||
# the data sets
|
||||
for label, dataset in datasets.items():
|
||||
dataset = np.asarray(dataset).flatten()
|
||||
# check that datasets are the correct size
|
||||
datasets_out = []
|
||||
if datasets is not None:
|
||||
for name, data in datasets.items():
|
||||
if data.shape != self.dimension:
|
||||
raise ValueError(f'Cannot apply dataset "{name}" with '
|
||||
f'shape {data.shape} to mesh {self.id} '
|
||||
f'with dimensions {self.dimension}')
|
||||
|
||||
if volume_normalization:
|
||||
dataset /= self.volumes.flatten()
|
||||
for name, data in datasets.items():
|
||||
if np.issubdtype(data.dtype, np.integer):
|
||||
warnings.warn(f'Integer data set "{name}" will '
|
||||
'not be volume-normalized.')
|
||||
continue
|
||||
data /= self.volumes
|
||||
|
||||
array = vtk.vtkDoubleArray()
|
||||
array.SetName(label)
|
||||
array.SetNumberOfComponents(1)
|
||||
array.SetArray(vtk_npsup.numpy_to_vtk(dataset),
|
||||
dataset.size,
|
||||
True)
|
||||
# add data to the mesh
|
||||
for name, data in datasets.items():
|
||||
datasets_out.append(data)
|
||||
arr = vtk.vtkDoubleArray()
|
||||
arr.SetName(name)
|
||||
arr.SetNumberOfTuples(data.size)
|
||||
|
||||
data_holder.append(dataset)
|
||||
poly_data.GetPointData().AddArray(array)
|
||||
for i in range(data.size):
|
||||
arr.SetTuple1(i, data.flat[i])
|
||||
grid.GetCellData().AddArray(arr)
|
||||
|
||||
# set filename
|
||||
if not filename.endswith(".vtk"):
|
||||
filename += ".vtk"
|
||||
writer.SetInputData(grid)
|
||||
|
||||
writer = vtk.vtkGenericDataObjectWriter()
|
||||
writer.SetFileName(filename)
|
||||
writer.SetInputData(poly_data)
|
||||
writer.Write()
|
||||
|
||||
@classmethod
|
||||
|
|
@ -1688,10 +1763,14 @@ class UnstructuredMesh(MeshBase):
|
|||
|
||||
mesh = cls(filename, library, mesh_id=mesh_id)
|
||||
vol_data = group['volumes'][()]
|
||||
centroids = group['centroids'][()]
|
||||
mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
|
||||
mesh.centroids = np.reshape(centroids, (vol_data.shape[0], 3))
|
||||
mesh.size = mesh.volumes.size
|
||||
mesh.n_elements = mesh.volumes.size
|
||||
|
||||
vertices = group['vertices'][()]
|
||||
mesh._vertices = vertices.reshape((-1, 3))
|
||||
connectivity = group['connectivity'][()]
|
||||
mesh._connectivity = connectivity.reshape((-1, 8))
|
||||
mesh._element_types = group['element_types'][()]
|
||||
|
||||
if 'length_multiplier' in group:
|
||||
mesh.length_multiplier = group['length_multiplier'][()]
|
||||
|
|
@ -1713,7 +1792,7 @@ class UnstructuredMesh(MeshBase):
|
|||
element.set("type", "unstructured")
|
||||
element.set("library", self._library)
|
||||
subelement = ET.SubElement(element, "filename")
|
||||
subelement.text = self.filename
|
||||
subelement.text = str(self.filename)
|
||||
|
||||
if self._length_multiplier != 1.0:
|
||||
element.set("length_multiplier", str(self.length_multiplier))
|
||||
|
|
|
|||
|
|
@ -14,8 +14,11 @@ class EqualityMixin:
|
|||
def __eq__(self, other):
|
||||
if isinstance(other, type(self)):
|
||||
for key, value in self.__dict__.items():
|
||||
if not np.array_equal(value, other.__dict__.get(key)):
|
||||
return False
|
||||
if isinstance(value, np.ndarray):
|
||||
if not np.array_equal(value, other.__dict__.get(key)):
|
||||
return False
|
||||
else:
|
||||
return value == other.__dict__.get(key)
|
||||
else:
|
||||
return False
|
||||
|
||||
|
|
|
|||
|
|
@ -259,13 +259,16 @@ class Region(ABC):
|
|||
clone[:] = [n.clone(memo) for n in self]
|
||||
return clone
|
||||
|
||||
def translate(self, vector, memo=None):
|
||||
def translate(self, vector, inplace=False, memo=None):
|
||||
"""Translate region in given direction
|
||||
|
||||
Parameters
|
||||
----------
|
||||
vector : iterable of float
|
||||
Direction in which region should be translated
|
||||
inplace : bool
|
||||
Whether or not to return a region based on new surfaces or one based
|
||||
on the original surfaces that have been modified.
|
||||
memo : dict or None
|
||||
Dictionary used for memoization. This parameter is used internally
|
||||
and should not be specified by the user.
|
||||
|
|
@ -279,7 +282,7 @@ class Region(ABC):
|
|||
|
||||
if memo is None:
|
||||
memo = {}
|
||||
return type(self)(n.translate(vector, memo) for n in self)
|
||||
return type(self)(n.translate(vector, inplace, memo) for n in self)
|
||||
|
||||
def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False,
|
||||
memo=None):
|
||||
|
|
@ -308,7 +311,7 @@ class Region(ABC):
|
|||
:math:`\psi` about z. This corresponds to an x-y-z extrinsic
|
||||
rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan
|
||||
angles :math:`(\phi, \theta, \psi)`.
|
||||
inplace : boolean
|
||||
inplace : bool
|
||||
Whether or not to return a new instance of Surface or to modify the
|
||||
coefficients of this Surface in place. Defaults to False.
|
||||
memo : dict or None
|
||||
|
|
@ -622,10 +625,10 @@ class Complement(Region):
|
|||
clone.node = self.node.clone(memo)
|
||||
return clone
|
||||
|
||||
def translate(self, vector, memo=None):
|
||||
def translate(self, vector, inplace=False, memo=None):
|
||||
if memo is None:
|
||||
memo = {}
|
||||
return type(self)(self.node.translate(vector, memo))
|
||||
return type(self)(self.node.translate(vector, inplace, memo))
|
||||
|
||||
def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False,
|
||||
memo=None):
|
||||
|
|
|
|||
|
|
@ -75,7 +75,7 @@ class Settings:
|
|||
Maximum number of lost particles
|
||||
|
||||
.. versionadded:: 0.12
|
||||
rel_max_lost_particles : int
|
||||
rel_max_lost_particles : float
|
||||
Maximum number of lost particles, relative to the total number of particles
|
||||
|
||||
.. versionadded:: 0.12
|
||||
|
|
@ -298,191 +298,191 @@ class Settings:
|
|||
self._max_tracks = None
|
||||
|
||||
@property
|
||||
def run_mode(self):
|
||||
def run_mode(self) -> str:
|
||||
return self._run_mode.value
|
||||
|
||||
@property
|
||||
def batches(self):
|
||||
def batches(self) -> int:
|
||||
return self._batches
|
||||
|
||||
@property
|
||||
def generations_per_batch(self):
|
||||
def generations_per_batch(self) -> int:
|
||||
return self._generations_per_batch
|
||||
|
||||
@property
|
||||
def inactive(self):
|
||||
def inactive(self) -> int:
|
||||
return self._inactive
|
||||
|
||||
@property
|
||||
def max_lost_particles(self):
|
||||
def max_lost_particles(self) -> int:
|
||||
return self._max_lost_particles
|
||||
|
||||
@property
|
||||
def rel_max_lost_particles(self):
|
||||
def rel_max_lost_particles(self) -> float:
|
||||
return self._rel_max_lost_particles
|
||||
|
||||
@property
|
||||
def particles(self):
|
||||
def particles(self) -> int:
|
||||
return self._particles
|
||||
|
||||
@property
|
||||
def keff_trigger(self):
|
||||
def keff_trigger(self) -> dict:
|
||||
return self._keff_trigger
|
||||
|
||||
@property
|
||||
def energy_mode(self):
|
||||
def energy_mode(self) -> str:
|
||||
return self._energy_mode
|
||||
|
||||
@property
|
||||
def max_order(self):
|
||||
def max_order(self) -> int:
|
||||
return self._max_order
|
||||
|
||||
@property
|
||||
def source(self):
|
||||
def source(self) -> typing.List[Source]:
|
||||
return self._source
|
||||
|
||||
@property
|
||||
def confidence_intervals(self):
|
||||
def confidence_intervals(self) -> bool:
|
||||
return self._confidence_intervals
|
||||
|
||||
@property
|
||||
def electron_treatment(self):
|
||||
def electron_treatment(self) -> str:
|
||||
return self._electron_treatment
|
||||
|
||||
@property
|
||||
def ptables(self):
|
||||
def ptables(self) -> bool:
|
||||
return self._ptables
|
||||
|
||||
@property
|
||||
def photon_transport(self):
|
||||
def photon_transport(self) -> bool:
|
||||
return self._photon_transport
|
||||
|
||||
@property
|
||||
def seed(self):
|
||||
def seed(self) -> int:
|
||||
return self._seed
|
||||
|
||||
@property
|
||||
def survival_biasing(self):
|
||||
def survival_biasing(self) -> bool:
|
||||
return self._survival_biasing
|
||||
|
||||
@property
|
||||
def entropy_mesh(self):
|
||||
def entropy_mesh(self) -> RegularMesh:
|
||||
return self._entropy_mesh
|
||||
|
||||
@property
|
||||
def trigger_active(self):
|
||||
def trigger_active(self) -> bool:
|
||||
return self._trigger_active
|
||||
|
||||
@property
|
||||
def trigger_max_batches(self):
|
||||
def trigger_max_batches(self) -> int:
|
||||
return self._trigger_max_batches
|
||||
|
||||
@property
|
||||
def trigger_batch_interval(self):
|
||||
def trigger_batch_interval(self) -> int:
|
||||
return self._trigger_batch_interval
|
||||
|
||||
@property
|
||||
def output(self):
|
||||
def output(self) -> dict:
|
||||
return self._output
|
||||
|
||||
@property
|
||||
def sourcepoint(self):
|
||||
def sourcepoint(self) -> dict:
|
||||
return self._sourcepoint
|
||||
|
||||
@property
|
||||
def statepoint(self):
|
||||
def statepoint(self) -> dict:
|
||||
return self._statepoint
|
||||
|
||||
@property
|
||||
def surf_source_read(self):
|
||||
def surf_source_read(self) -> dict:
|
||||
return self._surf_source_read
|
||||
|
||||
@property
|
||||
def surf_source_write(self):
|
||||
def surf_source_write(self) -> dict:
|
||||
return self._surf_source_write
|
||||
|
||||
@property
|
||||
def no_reduce(self):
|
||||
def no_reduce(self) -> bool:
|
||||
return self._no_reduce
|
||||
|
||||
@property
|
||||
def verbosity(self):
|
||||
def verbosity(self) -> int:
|
||||
return self._verbosity
|
||||
|
||||
@property
|
||||
def tabular_legendre(self):
|
||||
def tabular_legendre(self) -> dict:
|
||||
return self._tabular_legendre
|
||||
|
||||
@property
|
||||
def temperature(self):
|
||||
def temperature(self) -> dict:
|
||||
return self._temperature
|
||||
|
||||
@property
|
||||
def trace(self):
|
||||
def trace(self) -> typing.Iterable:
|
||||
return self._trace
|
||||
|
||||
@property
|
||||
def track(self):
|
||||
def track(self) -> typing.Iterable[typing.Iterable[int]]:
|
||||
return self._track
|
||||
|
||||
@property
|
||||
def cutoff(self):
|
||||
def cutoff(self) -> dict:
|
||||
return self._cutoff
|
||||
|
||||
@property
|
||||
def ufs_mesh(self):
|
||||
def ufs_mesh(self) -> RegularMesh:
|
||||
return self._ufs_mesh
|
||||
|
||||
@property
|
||||
def resonance_scattering(self):
|
||||
def resonance_scattering(self) -> dict:
|
||||
return self._resonance_scattering
|
||||
|
||||
@property
|
||||
def volume_calculations(self):
|
||||
def volume_calculations(self) -> typing.List[VolumeCalculation]:
|
||||
return self._volume_calculations
|
||||
|
||||
@property
|
||||
def create_fission_neutrons(self):
|
||||
def create_fission_neutrons(self) -> bool:
|
||||
return self._create_fission_neutrons
|
||||
|
||||
@property
|
||||
def delayed_photon_scaling(self):
|
||||
def delayed_photon_scaling(self) -> bool:
|
||||
return self._delayed_photon_scaling
|
||||
|
||||
@property
|
||||
def material_cell_offsets(self):
|
||||
def material_cell_offsets(self) -> bool:
|
||||
return self._material_cell_offsets
|
||||
|
||||
@property
|
||||
def log_grid_bins(self):
|
||||
def log_grid_bins(self) -> int:
|
||||
return self._log_grid_bins
|
||||
|
||||
@property
|
||||
def event_based(self):
|
||||
def event_based(self) -> bool:
|
||||
return self._event_based
|
||||
|
||||
@property
|
||||
def max_particles_in_flight(self):
|
||||
def max_particles_in_flight(self) -> int:
|
||||
return self._max_particles_in_flight
|
||||
|
||||
@property
|
||||
def write_initial_source(self):
|
||||
def write_initial_source(self) -> bool:
|
||||
return self._write_initial_source
|
||||
|
||||
@property
|
||||
def weight_windows(self):
|
||||
def weight_windows(self) -> typing.List[WeightWindows]:
|
||||
return self._weight_windows
|
||||
|
||||
@property
|
||||
def weight_windows_on(self):
|
||||
def weight_windows_on(self) -> bool:
|
||||
return self._weight_windows_on
|
||||
|
||||
@property
|
||||
def max_splits(self):
|
||||
def max_splits(self) -> int:
|
||||
return self._max_splits
|
||||
|
||||
@property
|
||||
def max_tracks(self):
|
||||
def max_tracks(self) -> int:
|
||||
return self._max_tracks
|
||||
|
||||
@run_mode.setter
|
||||
|
|
@ -517,7 +517,7 @@ class Settings:
|
|||
self._max_lost_particles = max_lost_particles
|
||||
|
||||
@rel_max_lost_particles.setter
|
||||
def rel_max_lost_particles(self, rel_max_lost_particles: int):
|
||||
def rel_max_lost_particles(self, rel_max_lost_particles: float):
|
||||
cv.check_type('rel_max_lost_particles', rel_max_lost_particles, Real)
|
||||
cv.check_greater_than('rel_max_lost_particles', rel_max_lost_particles, 0)
|
||||
cv.check_less_than('rel_max_lost_particles', rel_max_lost_particles, 1)
|
||||
|
|
|
|||
|
|
@ -11,7 +11,7 @@ from uncertainties import ufloat
|
|||
import openmc
|
||||
import openmc.checkvalue as cv
|
||||
|
||||
_VERSION_STATEPOINT = 17
|
||||
_VERSION_STATEPOINT = 18
|
||||
|
||||
|
||||
class StatePoint:
|
||||
|
|
|
|||
|
|
@ -336,9 +336,9 @@ class Surface(IDManagerMixin, ABC):
|
|||
----------
|
||||
vector : iterable of float
|
||||
Direction in which surface should be translated
|
||||
inplace : boolean
|
||||
inplace : bool
|
||||
Whether or not to return a new instance of this Surface or to
|
||||
modify the coefficients of this Surface. Defaults to False
|
||||
modify the coefficients of this Surface.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -374,7 +374,7 @@ class Surface(IDManagerMixin, ABC):
|
|||
:math:`\psi` about z. This corresponds to an x-y-z extrinsic
|
||||
rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan
|
||||
angles :math:`(\phi, \theta, \psi)`.
|
||||
inplace : boolean
|
||||
inplace : bool
|
||||
Whether or not to return a new instance of Surface or to modify the
|
||||
coefficients of this Surface in place. Defaults to False.
|
||||
|
||||
|
|
@ -568,9 +568,9 @@ class PlaneMixin:
|
|||
----------
|
||||
vector : iterable of float
|
||||
Direction in which surface should be translated
|
||||
inplace : boolean
|
||||
inplace : bool
|
||||
Whether or not to return a new instance of a Plane or to modify the
|
||||
coefficients of this plane. Defaults to False
|
||||
coefficients of this plane.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -1012,9 +1012,8 @@ class QuadricMixin:
|
|||
----------
|
||||
vector : iterable of float
|
||||
Direction in which surface should be translated
|
||||
inplace : boolean
|
||||
inplace : bool
|
||||
Whether to return a clone of the Surface or the Surface itself.
|
||||
Defaults to False
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -2563,7 +2562,7 @@ class Halfspace(Region):
|
|||
clone.surface = self.surface.clone(memo)
|
||||
return clone
|
||||
|
||||
def translate(self, vector, memo=None):
|
||||
def translate(self, vector, inplace=False, memo=None):
|
||||
"""Translate half-space in given direction
|
||||
|
||||
Parameters
|
||||
|
|
@ -2585,7 +2584,7 @@ class Halfspace(Region):
|
|||
# If translated surface not in memo, add it
|
||||
key = (self.surface, tuple(vector))
|
||||
if key not in memo:
|
||||
memo[key] = self.surface.translate(vector)
|
||||
memo[key] = self.surface.translate(vector, inplace)
|
||||
|
||||
# Return translated half-space
|
||||
return type(self)(memo[key], self.side)
|
||||
|
|
@ -2617,7 +2616,7 @@ class Halfspace(Region):
|
|||
:math:`\psi` about z. This corresponds to an x-y-z extrinsic
|
||||
rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan
|
||||
angles :math:`(\phi, \theta, \psi)`.
|
||||
inplace : boolean
|
||||
inplace : bool
|
||||
Whether or not to return a new instance of Surface or to modify the
|
||||
coefficients of this Surface in place. Defaults to False.
|
||||
memo : dict or None
|
||||
|
|
|
|||
|
|
@ -12,9 +12,12 @@ import numpy as np
|
|||
|
||||
import openmc
|
||||
import openmc.checkvalue as cv
|
||||
|
||||
from ._xml import get_text
|
||||
from .checkvalue import check_type, check_value
|
||||
from .mixin import IDManagerMixin
|
||||
from .plots import _SVG_COLORS
|
||||
from .surface import _BOUNDARY_TYPES
|
||||
|
||||
|
||||
class UniverseBase(ABC, IDManagerMixin):
|
||||
|
|
@ -713,6 +716,78 @@ class DAGMCUniverse(UniverseBase):
|
|||
dagmc_element.set('filename', self.filename)
|
||||
xml_element.append(dagmc_element)
|
||||
|
||||
def bounding_region(self, bounded_type='box', boundary_type='vacuum'):
|
||||
"""Creates a either a spherical or box shaped bounding region around
|
||||
the DAGMC geometry.
|
||||
Parameters
|
||||
----------
|
||||
bounded_type : str
|
||||
The type of bounding surface(s) to use when constructing the region.
|
||||
Options include a single spherical surface (sphere) or a rectangle
|
||||
made from six planes (box).
|
||||
boundary_type : str
|
||||
Boundary condition that defines the behavior for particles hitting
|
||||
the surface. Defaults to vacuum boundary condition. Passed into the
|
||||
surface construction.
|
||||
Returns
|
||||
-------
|
||||
openmc.Region
|
||||
Region instance
|
||||
"""
|
||||
|
||||
check_type('boundary type', boundary_type, str)
|
||||
check_value('boundary type', boundary_type, _BOUNDARY_TYPES)
|
||||
check_type('bounded type', bounded_type, str)
|
||||
check_value('bounded type', bounded_type, ('box', 'sphere'))
|
||||
|
||||
bounding_box = self.bounding_box
|
||||
|
||||
if bounded_type == 'sphere':
|
||||
import math
|
||||
bounding_box_center = (bounding_box[0] + bounding_box[1])/2
|
||||
radius = math.dist(bounding_box[0], bounding_box[1])
|
||||
bounding_surface = openmc.Sphere(
|
||||
x0=bounding_box_center[0],
|
||||
y0=bounding_box_center[1],
|
||||
z0=bounding_box_center[2],
|
||||
boundary_type=boundary_type,
|
||||
r=radius,
|
||||
)
|
||||
|
||||
return -bounding_surface
|
||||
|
||||
if bounded_type == 'box':
|
||||
# defines plane surfaces for all six faces of the bounding box
|
||||
lower_x = openmc.XPlane(bounding_box[0][0], boundary_type=boundary_type)
|
||||
upper_x = openmc.XPlane(bounding_box[1][0], boundary_type=boundary_type)
|
||||
lower_y = openmc.YPlane(bounding_box[0][1], boundary_type=boundary_type)
|
||||
upper_y = openmc.YPlane(bounding_box[1][1], boundary_type=boundary_type)
|
||||
lower_z = openmc.ZPlane(bounding_box[0][2], boundary_type=boundary_type)
|
||||
upper_z = openmc.ZPlane(bounding_box[1][2], boundary_type=boundary_type)
|
||||
|
||||
return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z
|
||||
|
||||
def bounded_universe(self, bounding_cell_id=10000, **kwargs):
|
||||
"""Returns an openmc.Universe filled with this DAGMCUniverse and bounded
|
||||
with a cell. Defaults to a box cell with a vacuum surface however this
|
||||
can be changed using the kwargs which are passed directly to
|
||||
DAGMCUniverse.bounding_region().
|
||||
|
||||
Parameters
|
||||
----------
|
||||
bounding_cell_id : int
|
||||
The cell ID number to use for the bounding cell, defaults to 1000 to reduce
|
||||
the chance of overlapping ID numbers with the DAGMC geometry.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.Universe
|
||||
Universe instance
|
||||
"""
|
||||
|
||||
bounding_cell = openmc.Cell(fill=self, cell_id =bounding_cell_id, region=self.bounding_region(**kwargs))
|
||||
return openmc.Universe(cells=[bounding_cell])
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, group):
|
||||
"""Create DAGMC universe from HDF5 group
|
||||
|
|
|
|||
|
|
@ -291,10 +291,10 @@ class WeightWindows(IDManagerMixin):
|
|||
subelement.text = ' '.join(str(e) for e in self.energy_bounds)
|
||||
|
||||
subelement = ET.SubElement(element, 'lower_ww_bounds')
|
||||
subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds)
|
||||
subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds.ravel('F'))
|
||||
|
||||
subelement = ET.SubElement(element, 'upper_ww_bounds')
|
||||
subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds)
|
||||
subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds.ravel('F'))
|
||||
|
||||
subelement = ET.SubElement(element, 'survival_ratio')
|
||||
subelement.text = str(self.survival_ratio)
|
||||
|
|
|
|||
|
|
@ -183,7 +183,9 @@ void DAGUniverse::init_geometry()
|
|||
} else {
|
||||
warning(fmt::format("DAGMC Cell IDs: {}", dagmc_ids_for_dim(3)));
|
||||
fatal_error(fmt::format("Cell ID {} exists in both DAGMC Universe {} "
|
||||
"and the CSG geometry.",
|
||||
"and the CSG geometry. Setting auto_geom_ids "
|
||||
"to True when initiating the DAGMC Universe may "
|
||||
"resolve this issue",
|
||||
c->id_, this->id_));
|
||||
}
|
||||
|
||||
|
|
|
|||
44
src/endf.cpp
44
src/endf.cpp
|
|
@ -90,23 +90,25 @@ bool is_inelastic_scatter(int mt)
|
|||
|
||||
unique_ptr<Function1D> read_function(hid_t group, const char* name)
|
||||
{
|
||||
hid_t dset = open_dataset(group, name);
|
||||
hid_t obj_id = open_object(group, name);
|
||||
std::string func_type;
|
||||
read_attribute(dset, "type", func_type);
|
||||
read_attribute(obj_id, "type", func_type);
|
||||
unique_ptr<Function1D> func;
|
||||
if (func_type == "Tabulated1D") {
|
||||
func = make_unique<Tabulated1D>(dset);
|
||||
func = make_unique<Tabulated1D>(obj_id);
|
||||
} else if (func_type == "Polynomial") {
|
||||
func = make_unique<Polynomial>(dset);
|
||||
func = make_unique<Polynomial>(obj_id);
|
||||
} else if (func_type == "CoherentElastic") {
|
||||
func = make_unique<CoherentElasticXS>(dset);
|
||||
func = make_unique<CoherentElasticXS>(obj_id);
|
||||
} else if (func_type == "IncoherentElastic") {
|
||||
func = make_unique<IncoherentElasticXS>(dset);
|
||||
func = make_unique<IncoherentElasticXS>(obj_id);
|
||||
} else if (func_type == "Sum") {
|
||||
func = make_unique<Sum1D>(obj_id);
|
||||
} else {
|
||||
throw std::runtime_error {"Unknown function type " + func_type +
|
||||
" for dataset " + object_name(dset)};
|
||||
" for dataset " + object_name(obj_id)};
|
||||
}
|
||||
close_dataset(dset);
|
||||
close_object(obj_id);
|
||||
return func;
|
||||
}
|
||||
|
||||
|
|
@ -271,4 +273,30 @@ double IncoherentElasticXS::operator()(double E) const
|
|||
return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W));
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Sum1D implementation
|
||||
//==============================================================================
|
||||
|
||||
Sum1D::Sum1D(hid_t group)
|
||||
{
|
||||
// Get number of functions
|
||||
int n;
|
||||
read_attribute(group, "n", n);
|
||||
|
||||
// Get each function
|
||||
for (int i = 0; i < n; ++i) {
|
||||
auto dset_name = fmt::format("func_{}", i + 1);
|
||||
functions_.push_back(read_function(group, dset_name.c_str()));
|
||||
}
|
||||
}
|
||||
|
||||
double Sum1D::operator()(double x) const
|
||||
{
|
||||
double result = 0.0;
|
||||
for (auto& func : functions_) {
|
||||
result += (*func)(x);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -118,6 +118,12 @@ void close_group(hid_t group_id)
|
|||
fatal_error("Failed to close group");
|
||||
}
|
||||
|
||||
void close_object(hid_t obj_id)
|
||||
{
|
||||
if (H5Oclose(obj_id) < 0)
|
||||
fatal_error("Failed to close object");
|
||||
}
|
||||
|
||||
int dataset_ndims(hid_t dset)
|
||||
{
|
||||
hid_t dspace = H5Dget_space(dset);
|
||||
|
|
@ -394,6 +400,12 @@ hid_t open_group(hid_t group_id, const char* name)
|
|||
return H5Gopen(group_id, name, H5P_DEFAULT);
|
||||
}
|
||||
|
||||
hid_t open_object(hid_t group_id, const std::string& name)
|
||||
{
|
||||
ensure_exists(group_id, name.c_str());
|
||||
return H5Oopen(group_id, name.c_str(), H5P_DEFAULT);
|
||||
}
|
||||
|
||||
void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer)
|
||||
{
|
||||
hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT);
|
||||
|
|
|
|||
171
src/mesh.cpp
171
src/mesh.cpp
|
|
@ -220,21 +220,58 @@ void UnstructuredMesh::to_hdf5(hid_t group) const
|
|||
write_dataset(mesh_group, "type", mesh_type);
|
||||
write_dataset(mesh_group, "filename", filename_);
|
||||
write_dataset(mesh_group, "library", this->library());
|
||||
// write volume of each element
|
||||
vector<double> tet_vols;
|
||||
xt::xtensor<double, 2> centroids({static_cast<size_t>(this->n_bins()), 3});
|
||||
for (int i = 0; i < this->n_bins(); i++) {
|
||||
tet_vols.emplace_back(this->volume(i));
|
||||
auto c = this->centroid(i);
|
||||
xt::view(centroids, i, xt::all()) = xt::xarray<double>({c.x, c.y, c.z});
|
||||
}
|
||||
|
||||
write_dataset(mesh_group, "volumes", tet_vols);
|
||||
write_dataset(mesh_group, "centroids", centroids);
|
||||
|
||||
if (specified_length_multiplier_)
|
||||
write_dataset(mesh_group, "length_multiplier", length_multiplier_);
|
||||
|
||||
// write vertex coordinates
|
||||
xt::xtensor<double, 2> vertices({static_cast<size_t>(this->n_vertices()), 3});
|
||||
for (int i = 0; i < this->n_vertices(); i++) {
|
||||
auto v = this->vertex(i);
|
||||
xt::view(vertices, i, xt::all()) = xt::xarray<double>({v.x, v.y, v.z});
|
||||
}
|
||||
write_dataset(mesh_group, "vertices", vertices);
|
||||
|
||||
int num_elem_skipped = 0;
|
||||
|
||||
// write element types and connectivity
|
||||
vector<double> volumes;
|
||||
xt::xtensor<int, 2> connectivity ({static_cast<size_t>(this->n_bins()), 8});
|
||||
xt::xtensor<int, 2> elem_types ({static_cast<size_t>(this->n_bins()), 1});
|
||||
for (int i = 0; i < this->n_bins(); i++) {
|
||||
auto conn = this->connectivity(i);
|
||||
|
||||
volumes.emplace_back(this->volume(i));
|
||||
|
||||
// write linear tet element
|
||||
if (conn.size() == 4) {
|
||||
xt::view(elem_types, i, xt::all()) = static_cast<int>(ElementType::LINEAR_TET);
|
||||
xt::view(connectivity, i, xt::all()) = xt::xarray<int>({conn[0], conn[1], conn[2], conn[3],
|
||||
-1, -1, -1, -1});
|
||||
// write linear hex element
|
||||
} else if (conn.size() == 8) {
|
||||
xt::view(elem_types, i, xt::all()) = static_cast<int>(ElementType::LINEAR_HEX);
|
||||
xt::view(connectivity, i, xt::all()) = xt::xarray<int>({conn[0], conn[1], conn[2], conn[3],
|
||||
conn[4], conn[5], conn[6], conn[7]});
|
||||
} else {
|
||||
num_elem_skipped++;
|
||||
xt::view(elem_types, i, xt::all()) = static_cast<int>(ElementType::UNSUPPORTED);
|
||||
xt::view(connectivity, i, xt::all()) = -1;
|
||||
}
|
||||
}
|
||||
|
||||
// warn users that some elements were skipped
|
||||
if (num_elem_skipped > 0) {
|
||||
warning(fmt::format("The connectivity of {} elements "
|
||||
"on mesh {} were not written "
|
||||
"because they are not of type linear tet/hex.",
|
||||
num_elem_skipped, this->id_));
|
||||
}
|
||||
|
||||
write_dataset(mesh_group, "volumes", volumes);
|
||||
write_dataset(mesh_group, "connectivity", connectivity);
|
||||
write_dataset(mesh_group, "element_types", elem_types);
|
||||
|
||||
close_group(mesh_group);
|
||||
}
|
||||
|
||||
|
|
@ -1775,6 +1812,14 @@ void MOABMesh::initialize()
|
|||
filename_);
|
||||
}
|
||||
|
||||
// set member range of vertices
|
||||
int vertex_dim = 0;
|
||||
rval = mbi_->get_entities_by_dimension(0, vertex_dim, verts_);
|
||||
if (rval != moab::MB_SUCCESS) {
|
||||
fatal_error("Failed to get all vertex handles");
|
||||
}
|
||||
|
||||
|
||||
// make an entity set for all tetrahedra
|
||||
// this is used for convenience later in output
|
||||
rval = mbi_->create_meshset(moab::MESHSET_SET, tetset_);
|
||||
|
|
@ -2124,6 +2169,14 @@ std::pair<vector<double>, vector<double>> MOABMesh::plot(
|
|||
return {};
|
||||
}
|
||||
|
||||
int MOABMesh::get_vert_idx_from_handle(moab::EntityHandle vert) const {
|
||||
int idx = vert - verts_[0];
|
||||
if (idx >= n_vertices()) {
|
||||
fatal_error(fmt::format("Invalid vertex idx {} (# vertices {})", idx, n_vertices()));
|
||||
}
|
||||
return idx;
|
||||
}
|
||||
|
||||
int MOABMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const
|
||||
{
|
||||
int bin = eh - ehs_[0];
|
||||
|
|
@ -2191,6 +2244,46 @@ Position MOABMesh::centroid(int bin) const
|
|||
return {centroid[0], centroid[1], centroid[2]};
|
||||
}
|
||||
|
||||
int MOABMesh::n_vertices() const {
|
||||
return verts_.size();
|
||||
}
|
||||
|
||||
Position MOABMesh::vertex(int id) const {
|
||||
|
||||
moab::ErrorCode rval;
|
||||
|
||||
moab::EntityHandle vert = verts_[id];
|
||||
|
||||
moab::CartVect coords;
|
||||
rval = mbi_->get_coords(&vert, 1, coords.array());
|
||||
if (rval != moab::MB_SUCCESS) {
|
||||
fatal_error("Failed to get the coordinates of a vertex.");
|
||||
}
|
||||
|
||||
return {coords[0], coords[1], coords[2]};
|
||||
}
|
||||
|
||||
std::vector<int> MOABMesh::connectivity(int bin) const {
|
||||
moab::ErrorCode rval;
|
||||
|
||||
auto tet = get_ent_handle_from_bin(bin);
|
||||
|
||||
// look up the tet connectivity
|
||||
vector<moab::EntityHandle> conn;
|
||||
rval = mbi_->get_connectivity(&tet, 1, conn);
|
||||
if (rval != moab::MB_SUCCESS) {
|
||||
fatal_error("Failed to get connectivity of a mesh element.");
|
||||
return {};
|
||||
}
|
||||
|
||||
std::vector<int> verts(4);
|
||||
for (int i = 0; i < verts.size(); i++) {
|
||||
verts[i] = get_vert_idx_from_handle(conn[i]);
|
||||
}
|
||||
|
||||
return verts;
|
||||
}
|
||||
|
||||
std::pair<moab::Tag, moab::Tag> MOABMesh::get_score_tags(
|
||||
std::string score) const
|
||||
{
|
||||
|
|
@ -2324,16 +2417,37 @@ const std::string LibMesh::mesh_lib_type = "libmesh";
|
|||
|
||||
LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
|
||||
{
|
||||
// filename_ and length_multiplier_ will already be set by the UnstructuredMesh constructor
|
||||
set_mesh_pointer_from_filename(filename_);
|
||||
set_length_multiplier(length_multiplier_);
|
||||
initialize();
|
||||
}
|
||||
|
||||
LibMesh::LibMesh(const std::string& filename, double length_multiplier)
|
||||
// create the mesh from a pointer to a libMesh Mesh
|
||||
LibMesh::LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier)
|
||||
{
|
||||
filename_ = filename;
|
||||
m_ = &input_mesh;
|
||||
set_length_multiplier(length_multiplier);
|
||||
initialize();
|
||||
}
|
||||
|
||||
// create the mesh from an input file
|
||||
LibMesh::LibMesh(const std::string& filename, double length_multiplier)
|
||||
{
|
||||
set_mesh_pointer_from_filename(filename);
|
||||
set_length_multiplier(length_multiplier);
|
||||
initialize();
|
||||
}
|
||||
|
||||
void LibMesh::set_mesh_pointer_from_filename(const std::string& filename)
|
||||
{
|
||||
filename_ = filename;
|
||||
unique_m_ = make_unique<libMesh::Mesh>(*settings::libmesh_comm, n_dimension_);
|
||||
m_ = unique_m_.get();
|
||||
m_->read(filename_);
|
||||
}
|
||||
|
||||
// intialize from mesh file
|
||||
void LibMesh::initialize()
|
||||
{
|
||||
if (!settings::libmesh_comm) {
|
||||
|
|
@ -2344,14 +2458,14 @@ void LibMesh::initialize()
|
|||
// assuming that unstructured meshes used in OpenMC are 3D
|
||||
n_dimension_ = 3;
|
||||
|
||||
m_ = make_unique<libMesh::Mesh>(*settings::libmesh_comm, n_dimension_);
|
||||
m_->read(filename_);
|
||||
|
||||
if (specified_length_multiplier_) {
|
||||
libMesh::MeshTools::Modification::scale(*m_, length_multiplier_);
|
||||
}
|
||||
|
||||
m_->prepare_for_use();
|
||||
// if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh.
|
||||
// Otherwise assume that it is prepared by its owning application
|
||||
if (unique_m_) {
|
||||
m_->prepare_for_use();
|
||||
}
|
||||
|
||||
// ensure that the loaded mesh is 3 dimensional
|
||||
if (m_->mesh_dimension() != n_dimension_) {
|
||||
|
|
@ -2394,6 +2508,27 @@ Position LibMesh::centroid(int bin) const
|
|||
return {centroid(0), centroid(1), centroid(2)};
|
||||
}
|
||||
|
||||
int LibMesh::n_vertices() const
|
||||
{
|
||||
return m_->n_nodes();
|
||||
}
|
||||
|
||||
Position LibMesh::vertex(int vertex_id) const
|
||||
{
|
||||
const auto node_ref = m_->node_ref(vertex_id);
|
||||
return {node_ref(0), node_ref(1), node_ref(2)};
|
||||
}
|
||||
|
||||
std::vector<int> LibMesh::connectivity(int elem_id) const
|
||||
{
|
||||
std::vector<int> conn;
|
||||
const auto* elem_ptr = m_->elem_ptr(elem_id);
|
||||
for (int i = 0; i < elem_ptr->n_nodes(); i++) {
|
||||
conn.push_back(elem_ptr->node_id(i));
|
||||
}
|
||||
return conn;
|
||||
}
|
||||
|
||||
std::string LibMesh::library() const
|
||||
{
|
||||
return mesh_lib_type;
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@
|
|||
|
||||
#include <algorithm> // for transform, max
|
||||
#include <cstring> // for strlen
|
||||
#include <cstdio> // for stdout
|
||||
#include <ctime> // for time, localtime
|
||||
#include <fstream>
|
||||
#include <iomanip> // for setw, setprecision, put_time
|
||||
|
|
@ -93,7 +94,8 @@ void title()
|
|||
// Write number of OpenMP threads
|
||||
fmt::print(" OpenMP Threads | {}\n", omp_get_max_threads());
|
||||
#endif
|
||||
std::cout << std::endl;
|
||||
fmt::print("\n");
|
||||
std::fflush(stdout);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -133,7 +135,8 @@ void header(const char* msg, int level)
|
|||
|
||||
// Print header based on verbosity level.
|
||||
if (settings::verbosity >= level)
|
||||
std::cout << '\n' << out << "\n" << std::endl;
|
||||
fmt::print("\n{}\n\n", out);
|
||||
std::fflush(stdout);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -379,7 +382,8 @@ void print_generation()
|
|||
if (n > 1) {
|
||||
fmt::print(" {:8.5f} +/-{:8.5f}", simulation::keff, simulation::keff_std);
|
||||
}
|
||||
std::cout << std::endl;
|
||||
fmt::print("\n");
|
||||
std::fflush(stdout);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -546,7 +550,8 @@ void print_results()
|
|||
fmt::print(" Leakage Fraction = {:.5f}\n",
|
||||
gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n);
|
||||
}
|
||||
std::cout << std::endl;
|
||||
fmt::print("\n");
|
||||
std::fflush(stdout);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -332,4 +332,40 @@ void IncoherentInelasticAE::sample(
|
|||
mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// MixedElasticAE implementation
|
||||
//==============================================================================
|
||||
|
||||
MixedElasticAE::MixedElasticAE(
|
||||
hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs)
|
||||
: coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs)
|
||||
{
|
||||
// Read incoherent elastic distribution
|
||||
hid_t incoherent_group = open_group(group, "incoherent");
|
||||
std::string temp;
|
||||
read_attribute(incoherent_group, "type", temp);
|
||||
if (temp == "incoherent_elastic") {
|
||||
incoherent_dist_ = make_unique<IncoherentElasticAE>(incoherent_group);
|
||||
} else if (temp == "incoherent_elastic_discrete") {
|
||||
auto xs = dynamic_cast<const Tabulated1D*>(&incoh_xs);
|
||||
incoherent_dist_ =
|
||||
make_unique<IncoherentElasticAEDiscrete>(incoherent_group, xs->x());
|
||||
}
|
||||
close_group(incoherent_group);
|
||||
}
|
||||
|
||||
void MixedElasticAE::sample(
|
||||
double E_in, double& E_out, double& mu, uint64_t* seed) const
|
||||
{
|
||||
// Evaluate coherent and incoherent elastic cross sections
|
||||
double xs_coh = coherent_xs_(E_in);
|
||||
double xs_incoh = incoherent_xs_(E_in);
|
||||
|
||||
if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {
|
||||
coherent_dist_.sample(E_in, E_out, mu, seed);
|
||||
} else {
|
||||
incoherent_dist_->sample(E_in, E_out, mu, seed);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -399,4 +399,23 @@ void free_memory_source()
|
|||
model::external_sources.clear();
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// C API
|
||||
//==============================================================================
|
||||
|
||||
extern "C" int openmc_sample_external_source(
|
||||
size_t n, uint64_t* seed, void* sites)
|
||||
{
|
||||
if (!sites || !seed) {
|
||||
set_errmsg("Received null pointer.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
||||
auto sites_array = static_cast<SourceSite*>(sites);
|
||||
for (size_t i = 0; i < n; ++i) {
|
||||
sites_array[i] = sample_external_source(seed);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -818,6 +818,16 @@ void write_unstructured_mesh_results()
|
|||
if (!umesh->output_)
|
||||
continue;
|
||||
|
||||
if (umesh->library() == "moab") {
|
||||
if (mpi::master)
|
||||
warning(fmt::format(
|
||||
"Output for a MOAB mesh (mesh {}) was "
|
||||
"requested but will not be written. Please use the Python "
|
||||
"API to generated the desired VTK tetrahedral mesh.",
|
||||
umesh->id_));
|
||||
continue;
|
||||
}
|
||||
|
||||
// if this tally has more than one filter, print
|
||||
// warning and skip writing the mesh
|
||||
if (tally->filters().size() > 1) {
|
||||
|
|
@ -889,12 +899,10 @@ void write_unstructured_mesh_results()
|
|||
std::string filename = fmt::format("tally_{0}.{1:0{2}}", tally->id_,
|
||||
simulation::current_batch, batch_width);
|
||||
|
||||
if (umesh->library() == "moab" && !mpi::master)
|
||||
continue;
|
||||
|
||||
// Write the unstructured mesh and data to file
|
||||
umesh->write(filename);
|
||||
|
||||
// remove score data added for this mesh write
|
||||
umesh->remove_scores();
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -210,14 +210,22 @@ ThermalData::ThermalData(hid_t group)
|
|||
if (temp == "coherent_elastic") {
|
||||
auto xs = dynamic_cast<CoherentElasticXS*>(elastic_.xs.get());
|
||||
elastic_.distribution = make_unique<CoherentElasticAE>(*xs);
|
||||
} else {
|
||||
if (temp == "incoherent_elastic") {
|
||||
elastic_.distribution = make_unique<IncoherentElasticAE>(dgroup);
|
||||
} else if (temp == "incoherent_elastic_discrete") {
|
||||
auto xs = dynamic_cast<Tabulated1D*>(elastic_.xs.get());
|
||||
elastic_.distribution =
|
||||
make_unique<IncoherentElasticAEDiscrete>(dgroup, xs->x());
|
||||
}
|
||||
} else if (temp == "incoherent_elastic") {
|
||||
elastic_.distribution = make_unique<IncoherentElasticAE>(dgroup);
|
||||
} else if (temp == "incoherent_elastic_discrete") {
|
||||
auto xs = dynamic_cast<Tabulated1D*>(elastic_.xs.get());
|
||||
elastic_.distribution =
|
||||
make_unique<IncoherentElasticAEDiscrete>(dgroup, xs->x());
|
||||
} else if (temp == "mixed_elastic") {
|
||||
// Get coherent/incoherent cross sections
|
||||
auto mixed_xs = dynamic_cast<Sum1D*>(elastic_.xs.get());
|
||||
const auto& coh_xs =
|
||||
dynamic_cast<const CoherentElasticXS*>(mixed_xs->functions(0).get());
|
||||
const auto& incoh_xs = mixed_xs->functions(1).get();
|
||||
|
||||
// Create mixed elastic distribution
|
||||
elastic_.distribution =
|
||||
make_unique<MixedElasticAE>(dgroup, *coh_xs, *incoh_xs);
|
||||
}
|
||||
|
||||
close_group(elastic_group);
|
||||
|
|
|
|||
|
|
@ -1,8 +1,8 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<geometry>
|
||||
<cell fill="10" id="13" region="9 -10 11 -12 13 -14" universe="11" />
|
||||
<cell fill="12" id="13" region="22 -23 24 -25 26 -27" universe="13" />
|
||||
<dagmc_universe auto_geom_ids="true" filename="dagmc.h5m" id="9" />
|
||||
<lattice id="10">
|
||||
<lattice id="12">
|
||||
<pitch>24.0 24.0</pitch>
|
||||
<dimension>2 2</dimension>
|
||||
<lower_left>-24.0 -24.0</lower_left>
|
||||
|
|
@ -10,12 +10,12 @@
|
|||
9 9
|
||||
9 9 </universes>
|
||||
</lattice>
|
||||
<surface boundary="reflective" coeffs="-24.0" id="9" name="left" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="24.0" id="10" name="right" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="-24.0" id="11" name="front" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="24.0" id="12" name="back" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="-10.0" id="13" name="bottom" type="z-plane" />
|
||||
<surface boundary="reflective" coeffs="10.0" id="14" name="top" type="z-plane" />
|
||||
<surface boundary="reflective" coeffs="-24.0" id="22" name="left" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="24.0" id="23" name="right" type="x-plane" />
|
||||
<surface boundary="reflective" coeffs="-24.0" id="24" name="front" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="24.0" id="25" name="back" type="y-plane" />
|
||||
<surface boundary="reflective" coeffs="-10.0" id="26" name="bottom" type="z-plane" />
|
||||
<surface boundary="reflective" coeffs="10.0" id="27" name="top" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -46,11 +46,32 @@ class DAGMCUniverseTest(PyAPITestHarness):
|
|||
# create the DAGMC universe
|
||||
pincell_univ = openmc.DAGMCUniverse(filename='dagmc.h5m', auto_geom_ids=True)
|
||||
|
||||
# creates another DAGMC universe, this time with within a bounded cell
|
||||
bound_pincell_universe = openmc.DAGMCUniverse(filename='dagmc.h5m').bounded_universe()
|
||||
# uses the bound_dag_cell as the root argument to test the type checks in openmc.Geometry
|
||||
bound_pincell_geometry = openmc.Geometry(root=bound_pincell_universe)
|
||||
# assigns the bound_dag_geometry to the model to test the type checks in model.Geometry setter
|
||||
model.Geometry = bound_pincell_geometry
|
||||
|
||||
# checks that the bounding box is calculated correctly
|
||||
bounding_box = pincell_univ.bounding_box
|
||||
assert bounding_box[0].tolist() == [-25., -25., -25.]
|
||||
assert bounding_box[1].tolist() == [25., 25., 25.]
|
||||
|
||||
# checks that the bounding region is six surfaces each with a vacuum boundary type
|
||||
b_region = pincell_univ.bounding_region(bounded_type='box', boundary_type='vacuum')
|
||||
assert isinstance(b_region, openmc.Region)
|
||||
assert len(b_region.get_surfaces()) == 6
|
||||
for surface in list(b_region.get_surfaces().values()):
|
||||
assert surface.boundary_type == 'vacuum'
|
||||
|
||||
# checks that the bounding region is a single surface with a reflective boundary type
|
||||
b_region = pincell_univ.bounding_region(bounded_type='sphere', boundary_type='reflective')
|
||||
assert isinstance(b_region, openmc.Region)
|
||||
assert len(b_region.get_surfaces()) == 1
|
||||
for surface in list(b_region.get_surfaces().values()):
|
||||
assert surface.boundary_type == 'reflective'
|
||||
|
||||
# create a 2 x 2 lattice using the DAGMC pincell
|
||||
pitch = np.asarray((24.0, 24.0))
|
||||
lattice = openmc.RectLattice()
|
||||
|
|
|
|||
91
tests/regression_tests/unstructured_mesh/inputs_true.dat
Normal file
91
tests/regression_tests/unstructured_mesh/inputs_true.dat
Normal file
|
|
@ -0,0 +1,91 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<geometry>
|
||||
<cell id="1" material="1" name="fuel" region="1 -2 3 -4 5 -6" universe="1" />
|
||||
<cell id="2" material="2" name="clad" region="(-1 | 2 | -3 | 4 | -5 | 6) (7 -8 9 -10 11 -12)" universe="1" />
|
||||
<cell id="3" material="3" name="water" region="(-7 | 8 | -9 | 10 | -11 | 12) (13 -14 15 -16 17 -18)" universe="1" />
|
||||
<surface coeffs="-5.0" id="1" name="minimum x" type="x-plane" />
|
||||
<surface coeffs="5.0" id="2" name="maximum x" type="x-plane" />
|
||||
<surface coeffs="-5.0" id="3" name="minimum y" type="y-plane" />
|
||||
<surface coeffs="5.0" id="4" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-5.0" id="5" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="5.0" id="6" name="maximum z" type="z-plane" />
|
||||
<surface coeffs="-6.0" id="7" name="minimum x" type="x-plane" />
|
||||
<surface coeffs="6.0" id="8" name="maximum x" type="x-plane" />
|
||||
<surface coeffs="-6.0" id="9" name="minimum y" type="y-plane" />
|
||||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
<material depletable="true" id="1" name="fuel">
|
||||
<density units="g/cc" value="4.5" />
|
||||
<nuclide ao="1.0" name="U235" />
|
||||
</material>
|
||||
<material id="2" name="zircaloy">
|
||||
<density units="g/cc" value="5.77" />
|
||||
<nuclide ao="0.5145" name="Zr90" />
|
||||
<nuclide ao="0.1122" name="Zr91" />
|
||||
<nuclide ao="0.1715" name="Zr92" />
|
||||
<nuclide ao="0.1738" name="Zr94" />
|
||||
<nuclide ao="0.028" name="Zr96" />
|
||||
</material>
|
||||
<material id="3" name="water">
|
||||
<density units="atom/b-cm" value="0.07416" />
|
||||
<nuclide ao="2.0" name="H1" />
|
||||
<nuclide ao="1.0" name="O16" />
|
||||
</material>
|
||||
</materials>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<settings>
|
||||
<run_mode>fixed source</run_mode>
|
||||
<particles>1000</particles>
|
||||
<batches>10</batches>
|
||||
<source strength="1.0">
|
||||
<space origin="0.0 0.0 0.0" type="spherical">
|
||||
<r parameters="0.0 0.0" type="uniform" />
|
||||
<cos_theta type="discrete">
|
||||
<parameters>1.0 1.0</parameters>
|
||||
</cos_theta>
|
||||
<phi type="discrete">
|
||||
<parameters>0.0 1.0</parameters>
|
||||
</phi>
|
||||
</space>
|
||||
<energy type="discrete">
|
||||
<parameters>15000000.0 1.0</parameters>
|
||||
</energy>
|
||||
</source>
|
||||
</settings>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<tallies>
|
||||
<mesh id="1">
|
||||
<dimension>10 10 10</dimension>
|
||||
<lower_left>-10.0 -10.0 -10.0</lower_left>
|
||||
<upper_right>10.0 10.0 10.0</upper_right>
|
||||
</mesh>
|
||||
<mesh id="2" library="libmesh" type="unstructured">
|
||||
<filename>test_mesh_hexes.e</filename>
|
||||
</mesh>
|
||||
<filter id="1" type="mesh">
|
||||
<bins>1</bins>
|
||||
</filter>
|
||||
<filter id="2" type="mesh">
|
||||
<bins>2</bins>
|
||||
</filter>
|
||||
<tally id="1" name="regular mesh tally">
|
||||
<filters>1</filters>
|
||||
<scores>flux</scores>
|
||||
<estimator>collision</estimator>
|
||||
</tally>
|
||||
<tally id="2" name="unstructured mesh tally">
|
||||
<filters>2</filters>
|
||||
<scores>flux</scores>
|
||||
<estimator>collision</estimator>
|
||||
</tally>
|
||||
</tallies>
|
||||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-10.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="10.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -15,12 +15,12 @@
|
|||
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
|
||||
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
|
||||
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15" id="18" name="maximum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="13" name="minimum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="14" name="maximum x" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="15" name="minimum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="16" name="maximum y" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-15.0" id="17" name="minimum z" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="15.0" id="18" name="maximum z" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
|
|
|
|||
|
|
@ -1,6 +1,8 @@
|
|||
import filecmp
|
||||
import glob
|
||||
from itertools import product
|
||||
import os
|
||||
import warnings
|
||||
|
||||
import openmc
|
||||
import openmc.lib
|
||||
|
|
@ -9,15 +11,29 @@ import numpy as np
|
|||
import pytest
|
||||
from tests.testing_harness import PyAPITestHarness
|
||||
|
||||
TETS_PER_VOXEL = 12
|
||||
|
||||
|
||||
class UnstructuredMeshTest(PyAPITestHarness):
|
||||
|
||||
def __init__(self, statepoint_name, model, inputs_true, holes):
|
||||
ELEM_PER_VOXEL = 12
|
||||
|
||||
def __init__(self,
|
||||
statepoint_name,
|
||||
model,
|
||||
inputs_true='inputs_true.dat',
|
||||
holes=False,
|
||||
scale_factor=10.0):
|
||||
|
||||
super().__init__(statepoint_name, model, inputs_true)
|
||||
self.holes = holes # holes in the test mesh
|
||||
self.scale_bounding_cell(scale_factor)
|
||||
|
||||
def scale_bounding_cell(self, scale_factor):
|
||||
geometry = self._model.geometry
|
||||
for surface in geometry.get_all_surfaces().values():
|
||||
if surface.boundary_type != 'vacuum':
|
||||
continue
|
||||
for coeff in surface._coefficients:
|
||||
surface._coefficients[coeff] *= scale_factor
|
||||
|
||||
def _compare_results(self):
|
||||
with openmc.StatePoint(self._sp_name) as sp:
|
||||
|
|
@ -28,32 +44,27 @@ class UnstructuredMeshTest(PyAPITestHarness):
|
|||
flt = tally.find_filter(openmc.MeshFilter)
|
||||
|
||||
if isinstance(flt.mesh, openmc.RegularMesh):
|
||||
reg_mesh_data, reg_mesh_std_dev = self.get_mesh_tally_data(tally)
|
||||
reg_mesh_data = self.get_mesh_tally_data(tally)
|
||||
if self.holes:
|
||||
reg_mesh_data = np.delete(reg_mesh_data, self.holes)
|
||||
reg_mesh_std_dev = np.delete(reg_mesh_std_dev, self.holes)
|
||||
else:
|
||||
umesh_tally = tally
|
||||
unstructured_data, unstructured_std_dev = self.get_mesh_tally_data(tally, True)
|
||||
unstructured_data = self.get_mesh_tally_data(tally, True)
|
||||
|
||||
# we expect these results to be the same to within at least ten
|
||||
# decimal places
|
||||
decimals = 10 if umesh_tally.estimator == 'collision' else 8
|
||||
np.testing.assert_array_almost_equal(unstructured_data,
|
||||
reg_mesh_data,
|
||||
decimals)
|
||||
# we expect these results to be the same to within at least ten
|
||||
# decimal places
|
||||
decimals = 10 if umesh_tally.estimator == 'collision' else 8
|
||||
np.testing.assert_array_almost_equal(np.sort(unstructured_data),
|
||||
np.sort(reg_mesh_data),
|
||||
decimals)
|
||||
|
||||
@staticmethod
|
||||
def get_mesh_tally_data(tally, structured=False):
|
||||
def get_mesh_tally_data(self, tally, structured=False):
|
||||
data = tally.get_reshaped_data(value='mean')
|
||||
std_dev = tally.get_reshaped_data(value='std_dev')
|
||||
if structured:
|
||||
data.shape = (data.size // TETS_PER_VOXEL, TETS_PER_VOXEL)
|
||||
std_dev.shape = (std_dev.size // TETS_PER_VOXEL, TETS_PER_VOXEL)
|
||||
data = data.reshape((-1, self.ELEM_PER_VOXEL))
|
||||
else:
|
||||
data.shape = (data.size, 1)
|
||||
std_dev.shape = (std_dev.size, 1)
|
||||
return np.sum(data, axis=1), np.sum(std_dev, axis=1)
|
||||
return np.sum(data, axis=1)
|
||||
|
||||
def _cleanup(self):
|
||||
super()._cleanup()
|
||||
|
|
@ -64,35 +75,11 @@ class UnstructuredMeshTest(PyAPITestHarness):
|
|||
os.remove(f)
|
||||
|
||||
|
||||
param_values = (['libmesh', 'moab'], # mesh libraries
|
||||
['collision', 'tracklength'], # estimators
|
||||
[True, False], # geometry outside of the mesh
|
||||
[(333, 90, 77), None]) # location of holes in the mesh
|
||||
test_cases = []
|
||||
for i, (lib, estimator, ext_geom, holes) in enumerate(product(*param_values)):
|
||||
test_cases.append({'library' : lib,
|
||||
'estimator' : estimator,
|
||||
'external_geom' : ext_geom,
|
||||
'holes' : holes,
|
||||
'inputs_true' : 'inputs_true{}.dat'.format(i)})
|
||||
|
||||
|
||||
@pytest.mark.parametrize("test_opts", test_cases)
|
||||
def test_unstructured_mesh(test_opts):
|
||||
|
||||
@pytest.fixture
|
||||
def model():
|
||||
openmc.reset_auto_ids()
|
||||
|
||||
# skip the test if the library is not enabled
|
||||
if test_opts['library'] == 'moab' and not openmc.lib._dagmc_enabled():
|
||||
pytest.skip("DAGMC (and MOAB) mesh not enbaled in this build.")
|
||||
|
||||
if test_opts['library'] == 'libmesh' and not openmc.lib._libmesh_enabled():
|
||||
pytest.skip("LibMesh is not enabled in this build.")
|
||||
|
||||
# skip the tracklength test for libmesh
|
||||
if test_opts['library'] == 'libmesh' and \
|
||||
test_opts['estimator'] == 'tracklength':
|
||||
pytest.skip("Tracklength tallies are not supported using libmesh.")
|
||||
model = openmc.Model()
|
||||
|
||||
### Materials ###
|
||||
materials = openmc.Materials()
|
||||
|
|
@ -113,7 +100,7 @@ def test_unstructured_mesh(test_opts):
|
|||
water_mat.set_density("atom/b-cm", 0.07416)
|
||||
materials.append(water_mat)
|
||||
|
||||
materials.export_to_xml()
|
||||
model.materials = materials
|
||||
|
||||
### Geometry ###
|
||||
fuel_min_x = openmc.XPlane(-5.0, name="minimum x")
|
||||
|
|
@ -149,29 +136,26 @@ def test_unstructured_mesh(test_opts):
|
|||
+clad_min_z & -clad_max_z)
|
||||
clad_cell.fill = zirc_mat
|
||||
|
||||
if test_opts['external_geom']:
|
||||
bounds = (15, 15, 15)
|
||||
else:
|
||||
bounds = (10, 10, 10)
|
||||
|
||||
water_min_x = openmc.XPlane(x0=-bounds[0],
|
||||
# set bounding cell dimension to one
|
||||
# this will be updated later according to the test case parameters
|
||||
water_min_x = openmc.XPlane(x0=-1.0,
|
||||
name="minimum x",
|
||||
boundary_type='vacuum')
|
||||
water_max_x = openmc.XPlane(x0=bounds[0],
|
||||
water_max_x = openmc.XPlane(x0=1.0,
|
||||
name="maximum x",
|
||||
boundary_type='vacuum')
|
||||
|
||||
water_min_y = openmc.YPlane(y0=-bounds[1],
|
||||
water_min_y = openmc.YPlane(y0=-1.0,
|
||||
name="minimum y",
|
||||
boundary_type='vacuum')
|
||||
water_max_y = openmc.YPlane(y0=bounds[1],
|
||||
water_max_y = openmc.YPlane(y0=1.0,
|
||||
name="maximum y",
|
||||
boundary_type='vacuum')
|
||||
|
||||
water_min_z = openmc.ZPlane(z0=-bounds[2],
|
||||
water_min_z = openmc.ZPlane(z0=-1.0,
|
||||
name="minimum z",
|
||||
boundary_type='vacuum')
|
||||
water_max_z = openmc.ZPlane(z0=bounds[2],
|
||||
water_max_z = openmc.ZPlane(z0=1.0,
|
||||
name="maximum z",
|
||||
boundary_type='vacuum')
|
||||
|
||||
|
|
@ -185,9 +169,9 @@ def test_unstructured_mesh(test_opts):
|
|||
water_cell.fill = water_mat
|
||||
|
||||
# create a containing universe
|
||||
geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell])
|
||||
model.geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell])
|
||||
|
||||
### Tallies ###
|
||||
### Reference Tally ###
|
||||
|
||||
# create meshes and mesh filters
|
||||
regular_mesh = openmc.RegularMesh()
|
||||
|
|
@ -196,29 +180,11 @@ def test_unstructured_mesh(test_opts):
|
|||
regular_mesh.upper_right = (10.0, 10.0, 10.0)
|
||||
|
||||
regular_mesh_filter = openmc.MeshFilter(mesh=regular_mesh)
|
||||
|
||||
if test_opts['holes']:
|
||||
mesh_filename = "test_mesh_tets_w_holes.e"
|
||||
else:
|
||||
mesh_filename = "test_mesh_tets.e"
|
||||
|
||||
uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_opts['library'])
|
||||
uscd_filter = openmc.MeshFilter(mesh=uscd_mesh)
|
||||
|
||||
# create tallies
|
||||
tallies = openmc.Tallies()
|
||||
|
||||
regular_mesh_tally = openmc.Tally(name="regular mesh tally")
|
||||
regular_mesh_tally.filters = [regular_mesh_filter]
|
||||
regular_mesh_tally.scores = ['flux']
|
||||
regular_mesh_tally.estimator = test_opts['estimator']
|
||||
tallies.append(regular_mesh_tally)
|
||||
|
||||
uscd_tally = openmc.Tally(name="unstructured mesh tally")
|
||||
uscd_tally.filters = [uscd_filter]
|
||||
uscd_tally.scores = ['flux']
|
||||
uscd_tally.estimator = test_opts['estimator']
|
||||
tallies.append(uscd_tally)
|
||||
model.tallies = openmc.Tallies([regular_mesh_tally])
|
||||
|
||||
### Settings ###
|
||||
settings = openmc.Settings()
|
||||
|
|
@ -236,13 +202,92 @@ def test_unstructured_mesh(test_opts):
|
|||
source = openmc.Source(space=space, energy=energy)
|
||||
settings.source = source
|
||||
|
||||
model = openmc.model.Model(geometry=geometry,
|
||||
materials=materials,
|
||||
tallies=tallies,
|
||||
settings=settings)
|
||||
model.settings = settings
|
||||
|
||||
return model
|
||||
|
||||
|
||||
param_values = (['libmesh', 'moab'], # mesh libraries
|
||||
['collision', 'tracklength'], # estimators
|
||||
[True, False], # geometry outside of the mesh
|
||||
[(333, 90, 77), None]) # location of holes in the mesh
|
||||
test_cases = []
|
||||
for i, (lib, estimator, ext_geom, holes) in enumerate(product(*param_values)):
|
||||
test_cases.append({'library' : lib,
|
||||
'estimator' : estimator,
|
||||
'external_geom' : ext_geom,
|
||||
'holes' : holes,
|
||||
'inputs_true' : 'inputs_true{}.dat'.format(i)})
|
||||
|
||||
|
||||
@pytest.mark.parametrize("test_opts", test_cases)
|
||||
def test_unstructured_mesh_tets(model, test_opts):
|
||||
# skip the test if the library is not enabled
|
||||
if test_opts['library'] == 'moab' and not openmc.lib._dagmc_enabled():
|
||||
pytest.skip("DAGMC (and MOAB) mesh not enbaled in this build.")
|
||||
|
||||
if test_opts['library'] == 'libmesh' and not openmc.lib._libmesh_enabled():
|
||||
pytest.skip("LibMesh is not enabled in this build.")
|
||||
|
||||
# skip the tracklength test for libmesh
|
||||
if test_opts['library'] == 'libmesh' and \
|
||||
test_opts['estimator'] == 'tracklength':
|
||||
pytest.skip("Tracklength tallies are not supported using libmesh.")
|
||||
|
||||
if test_opts['holes']:
|
||||
mesh_filename = "test_mesh_tets_w_holes.e"
|
||||
else:
|
||||
mesh_filename = "test_mesh_tets.e"
|
||||
|
||||
# add reference mesh tally
|
||||
regular_mesh_tally = model.tallies[0]
|
||||
regular_mesh_tally.estimator = test_opts['estimator']
|
||||
|
||||
# add analagous unstructured mesh tally
|
||||
uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_opts['library'])
|
||||
uscd_filter = openmc.MeshFilter(mesh=uscd_mesh)
|
||||
|
||||
# create tallies
|
||||
uscd_tally = openmc.Tally(name="unstructured mesh tally")
|
||||
uscd_tally.filters = [uscd_filter]
|
||||
uscd_tally.scores = ['flux']
|
||||
uscd_tally.estimator = test_opts['estimator']
|
||||
model.tallies.append(uscd_tally)
|
||||
|
||||
# modify model geometry according to test opts
|
||||
if test_opts['external_geom']:
|
||||
scale_factor = 15.0
|
||||
else:
|
||||
scale_factor = 10.0
|
||||
|
||||
harness = UnstructuredMeshTest('statepoint.10.h5',
|
||||
model,
|
||||
test_opts['inputs_true'],
|
||||
test_opts['holes'])
|
||||
test_opts['holes'],
|
||||
scale_factor)
|
||||
harness.main()
|
||||
|
||||
|
||||
@pytest.mark.skipif(not openmc.lib._libmesh_enabled(),
|
||||
reason='LibMesh is not enabled in this build.')
|
||||
def test_unstructured_mesh_hexes(model):
|
||||
regular_mesh_tally = model.tallies[0]
|
||||
regular_mesh_tally.estimator = 'collision'
|
||||
|
||||
# add analagous unstructured mesh tally
|
||||
uscd_mesh = openmc.UnstructuredMesh('test_mesh_hexes.e', 'libmesh')
|
||||
uscd_filter = openmc.MeshFilter(mesh=uscd_mesh)
|
||||
|
||||
# create tallies
|
||||
uscd_tally = openmc.Tally(name="unstructured mesh tally")
|
||||
uscd_tally.filters = [uscd_filter]
|
||||
uscd_tally.scores = ['flux']
|
||||
uscd_tally.estimator = 'collision'
|
||||
model.tallies.append(uscd_tally)
|
||||
|
||||
harness = UnstructuredMeshTest('statepoint.10.h5',
|
||||
model)
|
||||
harness.ELEM_PER_VOXEL = 1
|
||||
|
||||
harness.main()
|
||||
|
||||
|
|
|
|||
1
tests/regression_tests/unstructured_mesh/test_mesh_hexes.e
Symbolic link
1
tests/regression_tests/unstructured_mesh/test_mesh_hexes.e
Symbolic link
|
|
@ -0,0 +1 @@
|
|||
test_mesh_hexes.exo
|
||||
BIN
tests/regression_tests/unstructured_mesh/test_mesh_hexes.exo
Normal file
BIN
tests/regression_tests/unstructured_mesh/test_mesh_hexes.exo
Normal file
Binary file not shown.
0
tests/unit_tests/mesh_to_vtk/__init__.py
Normal file
0
tests/unit_tests/mesh_to_vtk/__init__.py
Normal file
BIN
tests/unit_tests/mesh_to_vtk/hexes.exo
Normal file
BIN
tests/unit_tests/mesh_to_vtk/hexes.exo
Normal file
Binary file not shown.
76
tests/unit_tests/mesh_to_vtk/libmesh_hexes_ref.vtk
Normal file
76
tests/unit_tests/mesh_to_vtk/libmesh_hexes_ref.vtk
Normal file
|
|
@ -0,0 +1,76 @@
|
|||
# vtk DataFile Version 5.1
|
||||
vtk output
|
||||
ASCII
|
||||
DATASET UNSTRUCTURED_GRID
|
||||
POINTS 54 double
|
||||
-1 -1 2.5 -1 -1 1.5 -1 0 1.5
|
||||
-1 0 2.5 0 -1 2.5 0 -1 1.5
|
||||
0 0 1.5 0 0 2.5 -1 -1 0.5
|
||||
-1 0 0.5 0 -1 0.5 0 0 0.5
|
||||
-1 -1 -0.5 -1 0 -0.5 0 -1 -0.5
|
||||
0 0 -0.5 -1 -1 -1.5 -1 0 -1.5
|
||||
0 -1 -1.5 0 0 -1.5 -1 -1 -2.5
|
||||
-1 0 -2.5 0 -1 -2.5 0 0 -2.5
|
||||
-1 1 1.5 -1 1 2.5 0 1 1.5
|
||||
0 1 2.5 -1 1 0.5 0 1 0.5
|
||||
-1 1 -0.5 0 1 -0.5 -1 1 -1.5
|
||||
0 1 -1.5 -1 1 -2.5 0 1 -2.5
|
||||
1 -1 2.5 1 -1 1.5 1 0 1.5
|
||||
1 0 2.5 1 -1 0.5 1 0 0.5
|
||||
1 -1 -0.5 1 0 -0.5 1 -1 -1.5
|
||||
1 0 -1.5 1 -1 -2.5 1 0 -2.5
|
||||
1 1 1.5 1 1 2.5 1 1 0.5
|
||||
1 1 -0.5 1 1 -1.5 1 1 -2.5
|
||||
|
||||
CELLS 21 160
|
||||
OFFSETS vtktypeint64
|
||||
0 8 16 24 32 40 48 56 64
|
||||
72 80 88 96 104 112 120 128 136
|
||||
144 152 160
|
||||
CONNECTIVITY vtktypeint64
|
||||
0 1 2 3 4 5 6 7 1
|
||||
8 9 2 5 10 11 6 8 12
|
||||
13 9 10 14 15 11 12 16 17
|
||||
13 14 18 19 15 16 20 21 17
|
||||
18 22 23 19 3 2 24 25 7
|
||||
6 26 27 2 9 28 24 6 11
|
||||
29 26 9 13 30 28 11 15 31
|
||||
29 13 17 32 30 15 19 33 31
|
||||
17 21 34 32 19 23 35 33 4
|
||||
5 6 7 36 37 38 39 5 10
|
||||
11 6 37 40 41 38 10 14 15
|
||||
11 40 42 43 41 14 18 19 15
|
||||
42 44 45 43 18 22 23 19 44
|
||||
46 47 45 7 6 26 27 39 38
|
||||
48 49 6 11 29 26 38 41 50
|
||||
48 11 15 31 29 41 43 51 50
|
||||
15 19 33 31 43 45 52 51 19
|
||||
23 35 33 45 47 53 52
|
||||
CELL_TYPES 20
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
12
|
||||
|
||||
CELL_DATA 20
|
||||
FIELD FieldData 1
|
||||
ids 1 20 double
|
||||
0 1 2 3 4 5 6 7 8
|
||||
9 10 11 12 13 14 15 16 17
|
||||
18 19
|
||||
292
tests/unit_tests/mesh_to_vtk/libmesh_tets_ref.vtk
Normal file
292
tests/unit_tests/mesh_to_vtk/libmesh_tets_ref.vtk
Normal file
|
|
@ -0,0 +1,292 @@
|
|||
# vtk DataFile Version 5.1
|
||||
vtk output
|
||||
ASCII
|
||||
DATASET UNSTRUCTURED_GRID
|
||||
POINTS 58 double
|
||||
-0.02593964576 -1 -1.1195739536 -0.40239958217 -0.40962166746 -1.9256035383 -1 0.02593964576 -1.1195739536
|
||||
-0.2248833639 0.23189144433 -1.2841109811 0.02593964576 1 -1.1195739536 1 -0.02593964576 -1.1195739536
|
||||
-0.042928712675 0.066425810853 -0.54270728017 -1 1 -1.5 -1 1 -0.5
|
||||
-0.49265312381 0.49332495053 -1.9656359509 -0.02608137991 -1 1.121521147 0.30281888252 -0.30293622513 1.8067963894
|
||||
1 -0.02608137991 1.121521147 -0.23236342389 0.26314805583 1.2867556783 0.074337770194 0.074337770194 2.5
|
||||
-0.40005770996 -0.36909660956 1.8768871762 -0.4896743493 0.49451064621 1.952829049 -1 1 1.5
|
||||
0.02608137991 1 1.121521147 -1 1 0.5 -0.038567136385 0.14628887118 0.55111724874
|
||||
-1 0.02608137991 1.121521147 0 -1 -2.5 0.37543954218 -0.36967061971 -1.8691105401
|
||||
1 -1 -2.5 0.074337770194 -0.074337770194 -2.5 1 0 -2.5
|
||||
0.4135797166 0.38351746262 -1.9274856063 1 -1 -1.5 -1 0 -2.5
|
||||
-1 -1 -1.5 -1 -1 -2.5 -1 0 2.5
|
||||
-1 1 2.5 0 1 2.5 0.46519182209 0.49356920409 1.95286052
|
||||
0 -1 2.5 -1 -1 1.5 -1 -1 2.5
|
||||
1 -1 2.5 1 -1 1.5 1 0 2.5
|
||||
1 1 2.5 1 1 1.5 0 1 -2.5
|
||||
-1 1 -2.5 1 -1 -0.5 1 1 -0.5
|
||||
0.038669824604 1 0.0046021677516 -1 -1 0.5 -0.038669824604 -1 0.0046021677516
|
||||
1 1 0.5 -1 0.038669824604 0.0046021677516 1 -0.038669824604 0.0046021677516
|
||||
1 1 -2.5 1 1 -1.5 1 -1 0.5
|
||||
-1 -1 -0.5
|
||||
CELLS 155 616
|
||||
OFFSETS vtktypeint64
|
||||
0 4 8 12 16 20 24 28 32
|
||||
36 40 44 48 52 56 60 64 68
|
||||
72 76 80 84 88 92 96 100 104
|
||||
108 112 116 120 124 128 132 136 140
|
||||
144 148 152 156 160 164 168 172 176
|
||||
180 184 188 192 196 200 204 208 212
|
||||
216 220 224 228 232 236 240 244 248
|
||||
252 256 260 264 268 272 276 280 284
|
||||
288 292 296 300 304 308 312 316 320
|
||||
324 328 332 336 340 344 348 352 356
|
||||
360 364 368 372 376 380 384 388 392
|
||||
396 400 404 408 412 416 420 424 428
|
||||
432 436 440 444 448 452 456 460 464
|
||||
468 472 476 480 484 488 492 496 500
|
||||
504 508 512 516 520 524 528 532 536
|
||||
540 544 548 552 556 560 564 568 572
|
||||
576 580 584 588 592 596 600 604 608
|
||||
612 616
|
||||
CONNECTIVITY vtktypeint64
|
||||
0 1 2 3 4 3 5 6 7
|
||||
4 8 3 7 3 2 9 8 3
|
||||
4 6 10 11 12 13 14 13 15
|
||||
16 17 13 18 16 19 18 13 20
|
||||
10 13 21 15 22 1 0 23 22
|
||||
24 25 23 25 1 3 9 26 23
|
||||
5 27 24 28 26 23 0 1 3
|
||||
23 29 30 31 1 22 31 30 1
|
||||
2 3 1 9 29 2 1 9 22
|
||||
0 28 23 32 15 21 16 14 13
|
||||
11 15 33 17 34 16 33 34 32
|
||||
16 34 16 18 35 10 11 13 15
|
||||
18 13 12 35 36 11 10 15 36
|
||||
37 38 15 32 38 37 15 32 21
|
||||
17 16 36 10 37 15 36 39 40
|
||||
11 41 40 39 11 41 12 11 35
|
||||
17 21 13 16 10 13 12 20 12
|
||||
13 11 35 14 13 16 35 34 17
|
||||
18 16 42 34 43 35 41 12 40
|
||||
11 14 41 11 35 42 41 14 35
|
||||
0 3 5 23 0 5 3 6 44
|
||||
4 9 27 29 2 30 1 44 45
|
||||
29 9 29 45 7 9 0 46 5
|
||||
6 47 48 4 6 0 3 2 6
|
||||
4 5 3 27 49 10 50 20 8
|
||||
2 3 6 51 18 48 20 49 50
|
||||
52 20 19 13 21 20 51 53 12
|
||||
20 19 48 18 20 25 3 1 23
|
||||
7 4 3 9 25 3 23 27 44
|
||||
4 7 9 25 23 26 27 54 55
|
||||
44 27 54 26 55 27 18 12 13
|
||||
20 14 16 34 35 19 21 52 20
|
||||
30 2 0 1 22 25 29 1 5
|
||||
23 3 27 25 26 44 27 22 29
|
||||
31 1 54 44 26 27 25 44 9
|
||||
27 22 28 24 23 26 5 55 27
|
||||
25 24 26 23 22 25 1 23 25
|
||||
44 29 9 25 29 1 9 48 6
|
||||
53 20 8 4 48 6 10 40 56
|
||||
12 10 40 12 11 17 19 18 13
|
||||
30 57 0 2 25 9 3 27 7
|
||||
8 2 3 26 28 5 23 29 7
|
||||
2 9 28 0 5 23 37 10 49
|
||||
21 49 50 57 52 56 12 53 20
|
||||
10 12 56 20 43 18 12 35 32
|
||||
37 21 15 14 11 13 35 19 52
|
||||
8 48 19 52 48 20 41 39 36
|
||||
11 50 6 52 20 46 53 5 6
|
||||
4 3 9 27 44 7 45 9 55
|
||||
47 4 5 28 0 46 5 8 52
|
||||
2 6 57 2 52 6 50 56 53
|
||||
20 10 56 50 20 57 50 0 6
|
||||
37 10 21 15 33 32 17 16 34
|
||||
18 43 35 57 52 50 6 51 48
|
||||
47 53 18 16 13 35 50 56 46
|
||||
53 14 36 32 15 14 11 36 15
|
||||
17 21 19 13 10 21 13 20 51
|
||||
12 18 20 50 53 46 6 14 15
|
||||
32 16 42 43 41 35 14 41 36
|
||||
11 49 52 21 20 36 40 10 11
|
||||
47 53 48 6 49 21 10 20 21
|
||||
15 13 16 51 48 53 20 55 5
|
||||
4 27 44 55 4 27 47 4 5
|
||||
6 8 48 52 6 47 5 53 6
|
||||
50 53 6 20 57 0 2 6 50
|
||||
46 0 6 52 6 48 20 42 14
|
||||
34 35 43 18 51 12 41 43 12
|
||||
35 22 30 0 1 14 32 34 16
|
||||
36 38 32 15
|
||||
CELL_TYPES 154
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
|
||||
CELL_DATA 154
|
||||
FIELD FieldData 1
|
||||
ids 1 154 double
|
||||
0 1 2 3 4 5 6 7 8
|
||||
9 10 11 12 13 14 15 16 17
|
||||
18 19 20 21 22 23 24 25 26
|
||||
27 28 29 30 31 32 33 34 35
|
||||
36 37 38 39 40 41 42 43 44
|
||||
45 46 47 48 49 50 51 52 53
|
||||
54 55 56 57 58 59 60 61 62
|
||||
63 64 65 66 67 68 69 70 71
|
||||
72 73 74 75 76 77 78 79 80
|
||||
81 82 83 84 85 86 87 88 89
|
||||
90 91 92 93 94 95 96 97 98
|
||||
99 100 101 102 103 104 105 106 107
|
||||
108 109 110 111 112 113 114 115 116
|
||||
117 118 119 120 121 122 123 124 125
|
||||
126 127 128 129 130 131 132 133 134
|
||||
135 136 137 138 139 140 141 142 143
|
||||
144 145 146 147 148 149 150 151 152
|
||||
153
|
||||
1
tests/unit_tests/mesh_to_vtk/moab_tets_ref.vtk
Symbolic link
1
tests/unit_tests/mesh_to_vtk/moab_tets_ref.vtk
Symbolic link
|
|
@ -0,0 +1 @@
|
|||
libmesh_tets_ref.vtk
|
||||
69
tests/unit_tests/mesh_to_vtk/test.py
Normal file
69
tests/unit_tests/mesh_to_vtk/test.py
Normal file
|
|
@ -0,0 +1,69 @@
|
|||
import filecmp
|
||||
from itertools import product
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
|
||||
import openmc
|
||||
import openmc.lib
|
||||
import pytest
|
||||
|
||||
pytest.importorskip('vtk')
|
||||
|
||||
|
||||
def ids_func(param):
|
||||
return f"{param['library']}_{param['elem_type']}"
|
||||
|
||||
test_params = (['libmesh', 'moab'],
|
||||
['tets', 'hexes'])
|
||||
|
||||
test_cases = [
|
||||
{'library' : library, 'elem_type' : elem_type}
|
||||
for library, elem_type in product(*test_params)
|
||||
]
|
||||
|
||||
@pytest.mark.parametrize("test_opts", test_cases, ids=ids_func)
|
||||
def test_unstructured_mesh_to_vtk(run_in_tmpdir, request, test_opts):
|
||||
|
||||
if test_opts['library'] == 'moab' and test_opts['elem_type'] == 'hexes':
|
||||
pytest.skip('Hexes are not supported with MOAB')
|
||||
|
||||
if test_opts['library'] == 'libmesh' and not openmc.lib._libmesh_enabled():
|
||||
pytest.skip('LibMesh is not enabled in this build.')
|
||||
|
||||
if test_opts['library'] == 'moab' and not openmc.lib._dagmc_enabled():
|
||||
pytest.skip('DAGMC (and MOAB) mesh not enabled in this build.')
|
||||
|
||||
# pull in a simple model -- just need to create the statepoint file
|
||||
openmc.reset_auto_ids()
|
||||
model = openmc.examples.pwr_pin_cell()
|
||||
|
||||
if test_opts['elem_type'] == 'tets':
|
||||
filename = Path('tets.exo')
|
||||
else:
|
||||
filename = Path('hexes.exo')
|
||||
|
||||
# create a basic tally using the unstructured mesh
|
||||
umesh = openmc.UnstructuredMesh(request.node.path.parent / filename,
|
||||
test_opts['library'])
|
||||
umesh.output = False
|
||||
mesh_filter = openmc.MeshFilter(umesh)
|
||||
tally = openmc.Tally()
|
||||
tally.filters = [mesh_filter]
|
||||
tally.scores = ['flux']
|
||||
tally.estimator = 'collision'
|
||||
model.tallies = openmc.Tallies([tally])
|
||||
sp_file = model.run()
|
||||
|
||||
# check VTK output after reading mesh from statepoint file
|
||||
with openmc.StatePoint(sp_file) as sp:
|
||||
umesh = sp.meshes[umesh.id]
|
||||
|
||||
test_data = {'ids': np.arange(umesh.n_elements)}
|
||||
umesh.write_data_to_vtk('umesh.vtk',
|
||||
datasets=test_data,
|
||||
volume_normalization=False)
|
||||
|
||||
# compare file content with reference file
|
||||
ref_file = Path(f"{test_opts['library']}_{test_opts['elem_type']}_ref.vtk")
|
||||
assert filecmp.cmp('umesh.vtk', request.node.path.parent / ref_file)
|
||||
BIN
tests/unit_tests/mesh_to_vtk/tets.exo
Normal file
BIN
tests/unit_tests/mesh_to_vtk/tets.exo
Normal file
Binary file not shown.
292
tests/unit_tests/mesh_to_vtk/umesh.vtk
Normal file
292
tests/unit_tests/mesh_to_vtk/umesh.vtk
Normal file
|
|
@ -0,0 +1,292 @@
|
|||
# vtk DataFile Version 5.1
|
||||
vtk output
|
||||
ASCII
|
||||
DATASET UNSTRUCTURED_GRID
|
||||
POINTS 58 double
|
||||
-0.02593964576 -1 -1.1195739536 -0.40239958217 -0.40962166746 -1.9256035383 -1 0.02593964576 -1.1195739536
|
||||
-0.2248833639 0.23189144433 -1.2841109811 0.02593964576 1 -1.1195739536 1 -0.02593964576 -1.1195739536
|
||||
-0.042928712675 0.066425810853 -0.54270728017 -1 1 -1.5 -1 1 -0.5
|
||||
-0.49265312381 0.49332495053 -1.9656359509 -0.02608137991 -1 1.121521147 0.30281888252 -0.30293622513 1.8067963894
|
||||
1 -0.02608137991 1.121521147 -0.23236342389 0.26314805583 1.2867556783 0.074337770194 0.074337770194 2.5
|
||||
-0.40005770996 -0.36909660956 1.8768871762 -0.4896743493 0.49451064621 1.952829049 -1 1 1.5
|
||||
0.02608137991 1 1.121521147 -1 1 0.5 -0.038567136385 0.14628887118 0.55111724874
|
||||
-1 0.02608137991 1.121521147 0 -1 -2.5 0.37543954218 -0.36967061971 -1.8691105401
|
||||
1 -1 -2.5 0.074337770194 -0.074337770194 -2.5 1 0 -2.5
|
||||
0.4135797166 0.38351746262 -1.9274856063 1 -1 -1.5 -1 0 -2.5
|
||||
-1 -1 -1.5 -1 -1 -2.5 -1 0 2.5
|
||||
-1 1 2.5 0 1 2.5 0.46519182209 0.49356920409 1.95286052
|
||||
0 -1 2.5 -1 -1 1.5 -1 -1 2.5
|
||||
1 -1 2.5 1 -1 1.5 1 0 2.5
|
||||
1 1 2.5 1 1 1.5 0 1 -2.5
|
||||
-1 1 -2.5 1 -1 -0.5 1 1 -0.5
|
||||
0.038669824604 1 0.0046021677516 -1 -1 0.5 -0.038669824604 -1 0.0046021677516
|
||||
1 1 0.5 -1 0.038669824604 0.0046021677516 1 -0.038669824604 0.0046021677516
|
||||
1 1 -2.5 1 1 -1.5 1 -1 0.5
|
||||
-1 -1 -0.5
|
||||
CELLS 155 616
|
||||
OFFSETS vtktypeint64
|
||||
0 4 8 12 16 20 24 28 32
|
||||
36 40 44 48 52 56 60 64 68
|
||||
72 76 80 84 88 92 96 100 104
|
||||
108 112 116 120 124 128 132 136 140
|
||||
144 148 152 156 160 164 168 172 176
|
||||
180 184 188 192 196 200 204 208 212
|
||||
216 220 224 228 232 236 240 244 248
|
||||
252 256 260 264 268 272 276 280 284
|
||||
288 292 296 300 304 308 312 316 320
|
||||
324 328 332 336 340 344 348 352 356
|
||||
360 364 368 372 376 380 384 388 392
|
||||
396 400 404 408 412 416 420 424 428
|
||||
432 436 440 444 448 452 456 460 464
|
||||
468 472 476 480 484 488 492 496 500
|
||||
504 508 512 516 520 524 528 532 536
|
||||
540 544 548 552 556 560 564 568 572
|
||||
576 580 584 588 592 596 600 604 608
|
||||
612 616
|
||||
CONNECTIVITY vtktypeint64
|
||||
0 1 2 3 4 3 5 6 7
|
||||
4 8 3 7 3 2 9 8 3
|
||||
4 6 10 11 12 13 14 13 15
|
||||
16 17 13 18 16 19 18 13 20
|
||||
10 13 21 15 22 1 0 23 22
|
||||
24 25 23 25 1 3 9 26 23
|
||||
5 27 24 28 26 23 0 1 3
|
||||
23 29 30 31 1 22 31 30 1
|
||||
2 3 1 9 29 2 1 9 22
|
||||
0 28 23 32 15 21 16 14 13
|
||||
11 15 33 17 34 16 33 34 32
|
||||
16 34 16 18 35 10 11 13 15
|
||||
18 13 12 35 36 11 10 15 36
|
||||
37 38 15 32 38 37 15 32 21
|
||||
17 16 36 10 37 15 36 39 40
|
||||
11 41 40 39 11 41 12 11 35
|
||||
17 21 13 16 10 13 12 20 12
|
||||
13 11 35 14 13 16 35 34 17
|
||||
18 16 42 34 43 35 41 12 40
|
||||
11 14 41 11 35 42 41 14 35
|
||||
0 3 5 23 0 5 3 6 44
|
||||
4 9 27 29 2 30 1 44 45
|
||||
29 9 29 45 7 9 0 46 5
|
||||
6 47 48 4 6 0 3 2 6
|
||||
4 5 3 27 49 10 50 20 8
|
||||
2 3 6 51 18 48 20 49 50
|
||||
52 20 19 13 21 20 51 53 12
|
||||
20 19 48 18 20 25 3 1 23
|
||||
7 4 3 9 25 3 23 27 44
|
||||
4 7 9 25 23 26 27 54 55
|
||||
44 27 54 26 55 27 18 12 13
|
||||
20 14 16 34 35 19 21 52 20
|
||||
30 2 0 1 22 25 29 1 5
|
||||
23 3 27 25 26 44 27 22 29
|
||||
31 1 54 44 26 27 25 44 9
|
||||
27 22 28 24 23 26 5 55 27
|
||||
25 24 26 23 22 25 1 23 25
|
||||
44 29 9 25 29 1 9 48 6
|
||||
53 20 8 4 48 6 10 40 56
|
||||
12 10 40 12 11 17 19 18 13
|
||||
30 57 0 2 25 9 3 27 7
|
||||
8 2 3 26 28 5 23 29 7
|
||||
2 9 28 0 5 23 37 10 49
|
||||
21 49 50 57 52 56 12 53 20
|
||||
10 12 56 20 43 18 12 35 32
|
||||
37 21 15 14 11 13 35 19 52
|
||||
8 48 19 52 48 20 41 39 36
|
||||
11 50 6 52 20 46 53 5 6
|
||||
4 3 9 27 44 7 45 9 55
|
||||
47 4 5 28 0 46 5 8 52
|
||||
2 6 57 2 52 6 50 56 53
|
||||
20 10 56 50 20 57 50 0 6
|
||||
37 10 21 15 33 32 17 16 34
|
||||
18 43 35 57 52 50 6 51 48
|
||||
47 53 18 16 13 35 50 56 46
|
||||
53 14 36 32 15 14 11 36 15
|
||||
17 21 19 13 10 21 13 20 51
|
||||
12 18 20 50 53 46 6 14 15
|
||||
32 16 42 43 41 35 14 41 36
|
||||
11 49 52 21 20 36 40 10 11
|
||||
47 53 48 6 49 21 10 20 21
|
||||
15 13 16 51 48 53 20 55 5
|
||||
4 27 44 55 4 27 47 4 5
|
||||
6 8 48 52 6 47 5 53 6
|
||||
50 53 6 20 57 0 2 6 50
|
||||
46 0 6 52 6 48 20 42 14
|
||||
34 35 43 18 51 12 41 43 12
|
||||
35 22 30 0 1 14 32 34 16
|
||||
36 38 32 15
|
||||
CELL_TYPES 154
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
|
||||
CELL_DATA 154
|
||||
FIELD FieldData 1
|
||||
ids 1 154 double
|
||||
0 1 2 3 4 5 6 7 8
|
||||
9 10 11 12 13 14 15 16 17
|
||||
18 19 20 21 22 23 24 25 26
|
||||
27 28 29 30 31 32 33 34 35
|
||||
36 37 38 39 40 41 42 43 44
|
||||
45 46 47 48 49 50 51 52 53
|
||||
54 55 56 57 58 59 60 61 62
|
||||
63 64 65 66 67 68 69 70 71
|
||||
72 73 74 75 76 77 78 79 80
|
||||
81 82 83 84 85 86 87 88 89
|
||||
90 91 92 93 94 95 96 97 98
|
||||
99 100 101 102 103 104 105 106 107
|
||||
108 109 110 111 112 113 114 115 116
|
||||
117 118 119 120 121 122 123 124 125
|
||||
126 127 128 129 130 131 132 133 134
|
||||
135 136 137 138 139 140 141 142 143
|
||||
144 145 146 147 148 149 150 151 152
|
||||
153
|
||||
|
|
@ -262,3 +262,108 @@ def test_get_thermal_name():
|
|||
|
||||
# Names that don't remotely match anything
|
||||
assert f('boogie_monster') == 'c_boogie_monster'
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def fake_mixed_elastic():
|
||||
fake_tsl = openmc.data.ThermalScattering("c_D_in_7LiD", 1.9968, 4.9, [0.0253])
|
||||
fake_tsl.nuclides = ['H2']
|
||||
|
||||
# Create elastic reaction
|
||||
bragg_edges = [0.00370672, 0.00494229, 0.00988458, 0.01359131, 0.01482688,
|
||||
0.01976918, 0.02347589, 0.02471147, 0.02965376, 0.03336048,
|
||||
0.03953834, 0.04324506, 0.04448063, 0.04942292, 0.05312964,
|
||||
0.05436522, 0.05930751, 0.06301423, 0.0642498 , 0.06919209,
|
||||
0.07289881, 0.07907667, 0.08278339, 0.08401896, 0.08896126,
|
||||
0.09266798, 0.09390355, 0.09884584, 0.1025526 , 0.1037882 ,
|
||||
0.1087305 , 0.1124372 , 0.1186151 , 0.1223218 , 0.1235574 ,
|
||||
0.1284997 , 0.1322064 , 0.133442 , 0.142091 , 0.1433266 ,
|
||||
0.1482688 , 0.1519756 , 0.1581534 , 0.1618601 , 0.1630957 ,
|
||||
0.168038 , 0.1717447 , 0.1729803 , 0.1779226 , 0.1816293 ,
|
||||
0.1828649 , 0.1878072 , 0.1915139 , 0.1976918 , 0.2026341 ,
|
||||
0.2075763 , 0.2125186 , 0.2174609 , 0.2224032 , 0.2273455 ,
|
||||
0.2421724 , 0.2471147 , 0.252057 , 0.2569993 , 0.2619415 ,
|
||||
0.2668838 , 0.2767684 , 0.2817107 , 0.2915953 , 0.3064222 ,
|
||||
0.3261913 , 0.366965]
|
||||
factors = [0.00375735, 0.01386287, 0.02595574, 0.02992438, 0.03549502,
|
||||
0.03855745, 0.04058831, 0.04986305, 0.05703106, 0.05855471,
|
||||
0.06078031, 0.06212291, 0.06656602, 0.06930339, 0.0697072 ,
|
||||
0.07201456, 0.07263853, 0.07313129, 0.07465531, 0.07714482,
|
||||
0.07759976, 0.077809 , 0.07790282, 0.07927957, 0.08013058,
|
||||
0.08026637, 0.08073475, 0.08112202, 0.08123039, 0.08187171,
|
||||
0.08213756, 0.08218236, 0.08236572, 0.08240729, 0.08259795,
|
||||
0.08297893, 0.08300455, 0.08314566, 0.08315611, 0.08337715,
|
||||
0.08350026, 0.08350663, 0.08352815, 0.08353776, 0.0836098 ,
|
||||
0.08367017, 0.08367361, 0.0837242 , 0.08375069, 0.08375227,
|
||||
0.08377006, 0.08381488, 0.08381644, 0.08382698, 0.08386266,
|
||||
0.08387756, 0.08388445, 0.08388974, 0.08390341, 0.08391088,
|
||||
0.08391695, 0.08392361, 0.08392684, 0.08392818, 0.08393161,
|
||||
0.08393546, 0.08393685, 0.08393801, 0.08393976, 0.08394167,
|
||||
0.08394288, 0.08394398]
|
||||
coherent_xs = openmc.data.CoherentElastic(bragg_edges, factors)
|
||||
incoherent_xs = openmc.data.Tabulated1D([0.00370672, 0.00370672], [0.00370672, 0.00370672])
|
||||
elastic_xs = {'294K': openmc.data.Sum((coherent_xs, incoherent_xs))}
|
||||
coherent_dist = openmc.data.CoherentElasticAE(coherent_xs)
|
||||
incoherent_dist = openmc.data.IncoherentElasticAEDiscrete([
|
||||
[-0.6, -0.18, 0.18, 0.6], [-0.6, -0.18, 0.18, 0.6]
|
||||
])
|
||||
elastic_dist = {'294K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist)}
|
||||
fake_tsl.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist)
|
||||
|
||||
# Create inelastic reaction
|
||||
inelastic_xs = {'294K': openmc.data.Tabulated1D([1.0e-5, 4.9], [13.4, 3.35])}
|
||||
breakpoints = [3]
|
||||
interpolation = [2]
|
||||
energy = [1.0e-5, 4.3e-2, 4.9]
|
||||
energy_out = [
|
||||
openmc.data.Tabular([0.0002, 0.067, 0.146, 0.366], [0.25, 0.25, 0.25, 0.25]),
|
||||
openmc.data.Tabular([0.0001, 0.009, 0.137, 0.277], [0.25, 0.25, 0.25, 0.25]),
|
||||
openmc.data.Tabular([0.0579, 4.555, 4.803, 4.874], [0.25, 0.25, 0.25, 0.25]),
|
||||
]
|
||||
for eout in energy_out:
|
||||
eout.normalize()
|
||||
eout.c = eout.cdf()
|
||||
discrete = openmc.stats.Discrete([-0.9, -0.6, -0.3, -0.1, 0.1, 0.3, 0.6, 0.9], [1/8]*8)
|
||||
discrete.c = discrete.cdf()[1:]
|
||||
mu = [[discrete]*4]*3
|
||||
inelastic_dist = {'294K': openmc.data.IncoherentInelasticAE(
|
||||
breakpoints, interpolation, energy, energy_out, mu)}
|
||||
inelastic = openmc.data.ThermalScatteringReaction(inelastic_xs, inelastic_dist)
|
||||
fake_tsl.inelastic = inelastic
|
||||
|
||||
return fake_tsl
|
||||
|
||||
|
||||
def test_mixed_elastic(fake_mixed_elastic, run_in_tmpdir):
|
||||
# Write data to HDF5 and then read back
|
||||
original = fake_mixed_elastic
|
||||
original.export_to_hdf5('c_D_in_7LiD.h5')
|
||||
copy = openmc.data.ThermalScattering.from_hdf5('c_D_in_7LiD.h5')
|
||||
|
||||
# Make sure data did not change as a result of HDF5 writing/reading
|
||||
assert original == copy
|
||||
|
||||
# Create modified cross_sections.xml file that includes the above data
|
||||
xs = openmc.data.DataLibrary.from_xml()
|
||||
xs.register_file('c_D_in_7LiD.h5')
|
||||
xs.export_to_xml('cross_sections_mixed.xml')
|
||||
|
||||
# Create a minimal model that includes the new data and run it
|
||||
mat = openmc.Material()
|
||||
mat.add_nuclide('H2', 1.0)
|
||||
mat.add_nuclide('Li7', 1.0)
|
||||
mat.set_density('g/cm3', 1.0)
|
||||
mat.add_s_alpha_beta('c_D_in_7LiD')
|
||||
sph = openmc.Sphere(r=10.0, boundary_type="vacuum")
|
||||
cell = openmc.Cell(fill=mat, region=-sph)
|
||||
model = openmc.Model()
|
||||
model.geometry = openmc.Geometry([cell])
|
||||
model.materials = openmc.Materials([mat])
|
||||
model.materials.cross_sections = "cross_sections_mixed.xml"
|
||||
model.settings.particles = 1000
|
||||
model.settings.batches = 10
|
||||
model.settings.run_mode = 'fixed source'
|
||||
model.settings.source = openmc.Source(
|
||||
energy=openmc.stats.Discrete([3.0], [1.0]) # 3 eV source
|
||||
)
|
||||
model.run()
|
||||
|
|
|
|||
|
|
@ -807,3 +807,46 @@ def test_cell_rotation(pincell_model_w_univ, mpi_intracomm):
|
|||
cell.rotation = (180., 0., 0.)
|
||||
assert cell.rotation == pytest.approx([180., 0., 0.])
|
||||
openmc.lib.finalize()
|
||||
|
||||
|
||||
def test_sample_external_source(run_in_tmpdir, mpi_intracomm):
|
||||
# Define a simple model and export
|
||||
mat = openmc.Material()
|
||||
mat.add_nuclide('U235', 1.0e-2)
|
||||
sph = openmc.Sphere(r=100.0, boundary_type='vacuum')
|
||||
cell = openmc.Cell(fill=mat, region=-sph)
|
||||
model = openmc.Model()
|
||||
model.geometry = openmc.Geometry([cell])
|
||||
model.settings.source = openmc.Source(
|
||||
space=openmc.stats.Box([-5., -5., -5.], [5., 5., 5.]),
|
||||
angle=openmc.stats.Monodirectional((0., 0., 1.)),
|
||||
energy=openmc.stats.Discrete([1.0e5], [1.0])
|
||||
)
|
||||
model.settings.particles = 1000
|
||||
model.settings.batches = 10
|
||||
model.export_to_xml()
|
||||
|
||||
# Sample some particles and make sure they match specified source
|
||||
openmc.lib.init()
|
||||
particles = openmc.lib.sample_external_source(10, prn_seed=3)
|
||||
assert len(particles) == 10
|
||||
for p in particles:
|
||||
assert -5. < p.r[0] < 5.
|
||||
assert -5. < p.r[1] < 5.
|
||||
assert -5. < p.r[2] < 5.
|
||||
assert p.u[0] == 0.0
|
||||
assert p.u[1] == 0.0
|
||||
assert p.u[2] == 1.0
|
||||
assert p.E == 1.0e5
|
||||
|
||||
# Using the same seed should produce the same particles
|
||||
other_particles = openmc.lib.sample_external_source(10, prn_seed=3)
|
||||
assert len(other_particles) == 10
|
||||
for p1, p2 in zip(particles, other_particles):
|
||||
assert p1.r == p2.r
|
||||
assert p1.u == p2.u
|
||||
assert p1.E == p2.E
|
||||
assert p1.time == p2.time
|
||||
assert p1.wgt == p2.wgt
|
||||
|
||||
openmc.lib.finalize()
|
||||
|
|
|
|||
|
|
@ -25,6 +25,52 @@ def test_add_nuclide():
|
|||
with pytest.raises(ValueError):
|
||||
m.add_nuclide('H1', 1.0, 'oa')
|
||||
|
||||
def test_add_components():
|
||||
"""Test adding multipe elements or nuclides at once"""
|
||||
m = openmc.Material()
|
||||
components = {'H1': 2.0,
|
||||
'O16': 1.0,
|
||||
'Zr': 1.0,
|
||||
'O': 1.0,
|
||||
'U': {'percent': 1.0,
|
||||
'enrichment': 4.5},
|
||||
'Li': {'percent': 1.0,
|
||||
'enrichment': 60.0,
|
||||
'enrichment_target': 'Li7'},
|
||||
'H': {'percent': 1.0,
|
||||
'enrichment': 50.0,
|
||||
'enrichment_target': 'H2',
|
||||
'enrichment_type': 'wo'}}
|
||||
m.add_components(components)
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'U': {'percent': 1.0,
|
||||
'enrichment': 100.0}})
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'Pu': {'percent': 1.0,
|
||||
'enrichment': 3.0}})
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'U': {'percent': 1.0,
|
||||
'enrichment': 70.0,
|
||||
'enrichment_target':'U235'}})
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'He': {'percent': 1.0,
|
||||
'enrichment': 17.0,
|
||||
'enrichment_target': 'He6'}})
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'li': 1.0}) # should fail as 1st char is lowercase
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'LI': 1.0}) # should fail as 2nd char is uppercase
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'Xx': 1.0}) # should fail as Xx is not an element
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'n': 1.0}) # check to avoid n for neutron being accepted
|
||||
with pytest.raises(TypeError):
|
||||
m.add_components({'H1': '1.0'})
|
||||
with pytest.raises(TypeError):
|
||||
m.add_components({1.0: 'H1'}, percent_type = 'wo')
|
||||
with pytest.raises(ValueError):
|
||||
m.add_components({'H1': 1.0}, percent_type = 'oa')
|
||||
|
||||
|
||||
def test_remove_nuclide():
|
||||
"""Test removing nuclides."""
|
||||
|
|
@ -49,7 +95,7 @@ def test_remove_elements():
|
|||
assert m.nuclides[0].percent == 1.0
|
||||
|
||||
|
||||
def test_elements():
|
||||
def test_add_element():
|
||||
"""Test adding elements."""
|
||||
m = openmc.Material()
|
||||
m.add_element('Zr', 1.0)
|
||||
|
|
@ -74,7 +120,6 @@ def test_elements():
|
|||
with pytest.raises(ValueError):
|
||||
m.add_element('n', 1.0) # check to avoid n for neutron being accepted
|
||||
|
||||
|
||||
def test_elements_by_name():
|
||||
"""Test adding elements by name"""
|
||||
m = openmc.Material()
|
||||
|
|
@ -441,7 +486,7 @@ def test_activity_of_tritium():
|
|||
m1.add_nuclide("H3", 1)
|
||||
m1.set_density('g/cm3', 1)
|
||||
m1.volume = 1
|
||||
assert pytest.approx(m1.activity) == 3.559778e14
|
||||
assert pytest.approx(m1.activity) == 3.559778e14
|
||||
|
||||
|
||||
def test_activity_of_metastable():
|
||||
|
|
|
|||
|
|
@ -208,3 +208,17 @@ def test_from_expression(reset):
|
|||
# Opening parenthesis immediately after halfspace
|
||||
r = openmc.Region.from_expression('1(2|-3)', surfs)
|
||||
assert str(r) == '(1 (2 | -3))'
|
||||
|
||||
|
||||
def test_translate_inplace():
|
||||
sph = openmc.Sphere()
|
||||
x = openmc.XPlane()
|
||||
region = -sph & +x
|
||||
|
||||
# Translating a region should produce new surfaces
|
||||
region2 = region.translate((0.5, -6.7, 3.9), inplace=False)
|
||||
assert str(region) != str(region2)
|
||||
|
||||
# Translating a region in-place should *not* produce new surfaces
|
||||
region3 = region.translate((0.5, -6.7, 3.9), inplace=True)
|
||||
assert str(region) == str(region3)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue