diff --git a/CMakeLists.txt b/CMakeLists.txt index 9fe133a22e3..f7a75a44f44 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -370,10 +370,12 @@ list(APPEND libopenmc_SOURCES src/distribution_energy.cpp src/distribution_multi.cpp src/distribution_spatial.cpp + src/dnp_drift.cpp src/eigenvalue.cpp src/endf.cpp src/error.cpp src/event.cpp + src/field.cpp src/file_utils.cpp src/finalize.cpp src/geometry.cpp @@ -425,6 +427,7 @@ list(APPEND libopenmc_SOURCES src/simulation.cpp src/source.cpp src/state_point.cpp + src/streamline_integrator.cpp src/string_utils.cpp src/summary.cpp src/surface.cpp diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index 1777e795ad7..1ec8de48ad0 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -1030,3 +1030,12 @@ Functions :type scores: const int* :return: Return status (negative if an error occurred) :rtype: int + +.. c:function:: int openmc_temperature_field_set_temperature(int32_t index, double temperature) + + Set the temperature value of a given cell in the temperature field + + :param int32_t index: Index in the tempererature mesh + :param double temperature: Temperature in Kelvin + :return: Return status (negative if an error occurred) + :rtype: int diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index b3911c8b3ba..6b462e36f4c 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -170,6 +170,17 @@ Meshes openmc.SphericalMesh openmc.UnstructuredMesh +Fields +------ + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclassinherit.rst + + openmc.ScalarField + openmc.TemperatureField + Geometry Plotting ----------------- diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index 83af5580574..b33b1421f40 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -191,6 +191,37 @@ attribute, e.g., :attr:`Material.temperature` or :attr:`Cell.temperature` attributes, respectively. +Alternatively, temperatures can be specified using a temperature field composed +of a geometric mesh and a map that associates a temperature value to each mesh +cell. During the simulation, temperatures associated with particles are also +updated every time a temperature field cell surface is crossed. While a particle +is contained inside the temperature mesh, temperatures from the temperature field +take precedence over temperature declared for a cell, a material or globally. + +The following example shows how to specify temperatures using a temperarature +field based on a regular mesh: + +.. code-block:: python + + # Define a mesh (regular mesh) + dim = 5 + mesh = openmc.RegularMesh() + mesh.lower_left = (0., 0., 0.) + mesh.upper_right = (10.0, 10.0, 10.0) + mesh.dimension = (dim, dim, dim) + + # Define temperature values for each cell + temperature_values = [273.0 + i * 10 for i in range(dim**3)] + + # Create a temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Register the temperature field in the settings + settings = openmc.Settings() + settings.temperature_field = temperature_field + +.. note:: Temperature fields are currently limited to structured meshes only. + ----------------- Material Mixtures ----------------- diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 0b425a673dc..f3e1f027e12 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -378,6 +378,14 @@ enum class GeometryType { CSG, DAG }; // representations. This value represents no surface. constexpr int32_t SURFACE_NONE {0}; +//============================================================================== +// EVENT IDENTIFIER IN HISTORY-BASED TRANSPORT + +const int EVENT_UNDEFINED = 0; +const int EVENT_CROSS_SURFACE = 1; +const int EVENT_COLLIDE = 2; +const int EVENT_TIME_CUTOFF = 3; + } // namespace openmc #endif // OPENMC_CONSTANTS_H diff --git a/include/openmc/dnp_drift.h b/include/openmc/dnp_drift.h new file mode 100644 index 00000000000..ca5fed6447b --- /dev/null +++ b/include/openmc/dnp_drift.h @@ -0,0 +1,57 @@ +#ifndef OPENMC_DNP_DRIFT_H +#define OPENMC_DNP_DRIFT_H + +#include "openmc/particle_data.h" +#include "openmc/position.h" + +namespace openmc { + +const int DNP_DRIFT_TRANSPORT_MAX_ITER = + 1000; //! Maximum number of iterations allowed in the DNP transport loop. +const double DNP_DRIFT_DISTANCE_MIN = + 1.0E-13; //! Minimum distance between two integration steps before stopping + //! the DNP transport loop. + +enum class Actions { PLACE_AT_INLET, BLOCK_AT_OUTLET, BLOCK_AT_LOCATION }; + +//! Adjust position linearly. +//! +//! This function is used to adjust the position of a DNP that reached its decay +//! time between two integration steps (i.e., t is greater than decay_time). The +//! final position of the DNP is adjusted linearly so that the time +//! corresponding to this new position is equal to the sampled decay time. +//! +//! \param[inout] y_n Current position +//! \param[in] y_n_minus_1 Previous position +//! \param[in] t Time at current position +//! \param[in] dt Time step +//! \param[in] decay_time Decay time +void _adjust_position(Position& y_n, const Position& y_n_minus_1, double t, + double dt, double decay_time); + +//! Adjust time linearly. +//! +//! This function is used to adjust the time associated with the position of a +//! DNP that needs to be stopped between two integration steps. The position C, +//! where the DNP is stopped, must be located between position A and B. Time is +//! adjusted linearly. +//! +//! \param[inout] t In: time at position B, out: time at position C +//! \param[in] ta Time at position A +//! \param[in] pa Position A +//! \param[in] pb Position B +//! \param[in] pc Position C (between A and B) +void _adjust_time(double& t, double ta, const Position& pa, const Position& pb, + const Position& pc); + +//! Explicitly transport a delayed neutron precursor using streamline +//! integration. +//! +//! \param[inout] site Fission site corresponding to the DNP +//! \param[in] decay_time Sampled decay time (in seconds) +//! \param[in] seed Random number generator seed +bool transport_dnp(SourceSite& site, double decay_time, uint64_t* seed); + +} // namespace openmc + +#endif // OPENMC_DNP_DRIFT_H diff --git a/include/openmc/event.h b/include/openmc/event.h index 2d215a10e46..c97d161ce56 100644 --- a/include/openmc/event.h +++ b/include/openmc/event.h @@ -107,6 +107,10 @@ void process_surface_crossing_events(); //! Execute the collision event for all particles in this event's buffer void process_collision_events(); +//! Execute the temperature mesh crossing event for all particles in this +//! event's buffer +void process_temperature_mesh_crossing_events(); + //! Execute the death event for all particles // //! \param n_particles The number of particles in the particle buffer diff --git a/include/openmc/field.h b/include/openmc/field.h new file mode 100644 index 00000000000..6adeb050556 --- /dev/null +++ b/include/openmc/field.h @@ -0,0 +1,351 @@ +#ifndef OPENMC_FIELD_H +#define OPENMC_FIELD_H + +#include +#include +#include + +#include "openmc/mesh.h" +#include "openmc/position.h" +#include "openmc/vector.h" + +namespace openmc { + +// ----------------------------------------------------------- +// FieldData +// ----------------------------------------------------------- + +template +class Field; + +template +class FieldData { +public: + // Constructor + FieldData(vector values) { values_ = values; } + + //! Get the value for a given index. + // + //! \param[in] idx Index + //! \return Associated value + T evaluate(int idx) const { return values_[idx]; } + + //! Assign a value at a given index. + // + //! \param[in] idx Index + //! \param[in] value Value to store + void assign(int idx, T value) { values_[idx] = value; } + + //! Return the size of the data field. + // + //! \return Number of values + int size() { return values_.size(); } + + // Values accessors + vector& values() { return values_; } + const vector& values() const { return values_; } + +private: + vector values_; //!< Stored data +}; + +// ----------------------------------------------------------- +// Field +// ----------------------------------------------------------- + +enum class FieldMapping { + NODAL, // Nodal representation (values defined for each vertex) + CELL // Cell-based representation (values defined for each cell) +}; + +template +class Field { +public: + // Constructor + Field() = default; + + //! Set the mesh pointer. + //! Returns an error if the mesh pointer is not valid. + // + //! \param[in] value Mesh pointer + void set_mesh(Mesh* value) + { + if (value != nullptr) { + mesh_ = value; + } else { + fatal_error("No mesh found for this field!"); + } + } + + //! Set the mapping type. + //! Returns an error if the mapping type is not valid. + // + //! \param[in] value Mapping type + void set_mapping(std::string value) + { + if (value == "nodal") { + mapping_ = FieldMapping::NODAL; + } else if (value == "cell") { + mapping_ = FieldMapping::CELL; + } else { + fatal_error(fmt::format("Unrecognized mapping type: {}", value)); + } + } + + //! Set data field. + //! Returns an error if the size of the data field is not consistent with its + //! mapping type. + // + //! \param[in] data Data field + void set_data(std::unique_ptr> data) + { + // Values/mesh size consistency check + if (mapping() == FieldMapping::NODAL) { + if (mesh_ptr()->n_vertices() != data->size()) { + fatal_error("The number of bins in the mesh is not consistent with the " + "number of values declared for this field!"); + } + } else if (mapping() == FieldMapping::CELL) { + if (mesh_ptr()->n_bins() != data->size()) { + fatal_error("The number of bins in the mesh is not consistent with the " + "number of values declared for this field!"); + } + } + + data_ = std::move(data); + } + + //! Assign a value to a bin. + // + //! \param[in] bin Bin number + //! \param[in] value Value to store + void assign(int bin, T value) { data().assign(bin, value); } + + //! Return the mesh bin associated with a given position. + // + //! \param[in] r Position + //! \return Bin number + int get_mesh_bin(const Position& r) { return mesh_ptr()->get_bin(r); } + + //! Evaluate the field at a given position with prior knowledge of the current + //! bin. + // + //! \param[in] r Position + //! \param[in] bin Bin number corresponding to the current position + //! \return Value corresponding to the position + T evaluate(const Position& r, int bin) + { + if (bin != C_NONE) { + switch (mapping()) { + case FieldMapping::NODAL: + // TODO: implement other interpolation techniques + return trilinear_interpolation(r, bin); + break; + case FieldMapping::CELL: + return value(bin); + break; + default: + fatal_error("Not implemented for this mapping type!"); + break; + } + } else { + fatal_error("Bin outside the mesh."); + } + } + + //! Evaluate the field at a given position with prior knowledge of the bin and + //! of the previous position. We assume that the transport between the two + //! points is in a straight line, justifying the use of ray-tracing. If the + //! next point is outside the mesh, we use the last bin that was located + //! inside the mesh. + // + //! \param[in] r0 Previous position + //! \param[in] r1 Current position + //! \param[in] bin Bin corresponding to the previous position + //! \return Value corresponding to the current position + T evaluate(const Position& r0, const Position& r1, int bin) + { + int next_bin = mesh_ptr()->get_last_bin_inside_mesh(r0, r1, bin); + return evaluate(r1, next_bin); + } + + //! Interpolate data field at a given position using trilinear interpolation. + //! To avoid extrapolation, any normalized coordinate not in [0.0, 1.0] will + //! trigger an error. + // + //! \param[in] r Position + //! \param[in] bin Bin number corresponding to the position + //! \return Interpolated value + T trilinear_interpolation(const Position& r, int bin) + { + // Normalize coordinates + Position n_r = mesh_ptr()->normalize_position(r); + + // Check that normalized coordinates are in [0.0, 1.0] to avoid + // extrapolation + for (int i = 0; i < 3; i++) { + if ((n_r[i] < 0.0) || (n_r[i] > 1.0)) { + fatal_error( + "Normalized coordinates must be in [0.0, 1.0] for interpolation!"); + } + } + + // Retrieve vertices + vector v = mesh_ptr()->connectivity(bin); + + // Interpolate along x + T c00 = value(v[0]) * (1 - n_r[0]) + value(v[1]) * n_r[0]; + T c01 = value(v[4]) * (1 - n_r[0]) + value(v[5]) * n_r[0]; + T c10 = value(v[2]) * (1 - n_r[0]) + value(v[3]) * n_r[0]; + T c11 = value(v[6]) * (1 - n_r[0]) + value(v[7]) * n_r[0]; + + // Interpolate along y + T c0 = c00 * (1 - n_r[1]) + c10 * n_r[1]; + T c1 = c01 * (1 - n_r[1]) + c11 * n_r[1]; + + // Interpolate along z + T result = c0 * (1 - n_r[2]) + c1 * n_r[2]; + return result; + } + + //! Returns the distance to the next mesh boundary given a particle position + //! and direction. If the particle is initially outside, the distance will + //! correspond to the nearest distance to the outer boundaries of the mesh. + // + //! \param[in] current_bin Current bin number + //! \param[in] r Position of the particle + //! \param[in] u Direction of the particle + //! \param[out] bin_next Next bin number + //! \return The distance in cm to the next mesh boundary + double distance_to_next_boundary( + int current_bin, const Position& r, const Direction& u, int& bin_next) + { + return this->mesh_ptr()->distance_to_next_boundary( + current_bin, r, u, bin_next); + } + + // Mesh pointer accessor + Mesh* mesh_ptr() const + { + if (mesh_ == nullptr) { + fatal_error("No mesh found for this field!"); + } else { + return mesh_; + } + } + + // Mapping accessors + FieldMapping mapping() { return mapping_; } + const FieldMapping mapping() const { return mapping_; } + + // Data field accessor + FieldData& data() const + { + if (data_ == nullptr) { + fatal_error("No data found for this field!"); + } else { + return *data_; + } + } + + // Data field value accessors + const T value(int i) const { return data().evaluate(i); } + const vector values() const { return data().values(); } + +private: + Mesh* mesh_; //!< Pointer to the geometric mesh + FieldMapping mapping_; //!< Relationship between values and mesh + std::unique_ptr> data_; //!< Data associated with the mesh +}; + +// ----------------------------------------------------------- +// TemperatureField +// ----------------------------------------------------------- + +class TemperatureField : public Field { +public: + // Constructors + TemperatureField() : Field() {}; + TemperatureField( + Mesh* mesh_ptr, vector values, std::string mapping = "cell"); + + //! Returns the temperature in Kelvin corresponding to a given bin number + //! relative to the mesh. + // + //! \param[in] r Position of the particle + //! \return Temperature in Kelvin + double get_temperature(int bin); + + //! Returns the square root of the temperature multiplied by the Boltzmann + //! constant in eV for a given bin number relative to the mesh. + // + //! \param[in] bin Bin number + //! \return Sqrt(k_Boltzmann * temperature) in eV + double get_sqrtkT(int bin); +}; + +// ----------------------------------------------------------- +// VelocityField +// ----------------------------------------------------------- + +// Boundary conditions type +enum class BCType { NONE, INLET, OUTLET, WALL }; + +// Boundary conditions map type +using BCMap = std::unordered_map>; + +class VelocityField : public Field { +public: + // Constructors + VelocityField() : Field() {}; + VelocityField( + Mesh* mesh_ptr, vector values, std::string mapping); + + //! Find next bin associated with a given position (r1) knowing the previous + //! position (r0) and the previous bin (bin0). The next bin is evaluated using + //! raytracing. If r1 is outside the mesh, crossed_boundary will indicate + //! the boundary condition associated with the last surface crossed before + //! leaving the mesh. Also, intersection will be the last intersection with + //! the mesh. If r1 is still inside the mesh, crossed_boundary will be NONE + //! and the intersection will coincide with r1. + //! + //! We currently do not check whether the mesh was left during the travel, + //! meaning that a point can leave and reenter the mesh. It is not problematic + //! for convex geometries like regular meshes, but it can be for unstuctured + //! meshes. + // + //! \param[in] r0 First position + //! \param[in] r1 Second position + //! \param[in] bin0 Bin number corresponding to r0 + //! \param[out] crossed_boundary Boundary type of the last crossed surface + //! \param[out] intersection Last intersection with the mesh (or r1) + //! \return Bin number corresponding to r1 + int get_next_bin(const Position& r0, const Position& r1, int bin0, + BCType& crossed_boundary, Position& intersection); + + //! Update the position and the bin by moving the point to a random location + //! on the inlet physical group. + // + //! \param[inout] p Position + //! \param[inout] bin Bin number corresponding to p + //! \param[in] seed Random number generator seed + void randomly_place_on_inlet(Position& p, int& bin, uint64_t* seed); + + //! Retrieve the boundary condition associated with a physical group. + // + //! \param[in] physical_group Physical group + //! \return Boundary condition associated with the physical group + BCType get_boundary_condition(int physical_group); + + // Boundary conditions map accessors + BCMap& bc_map() { return bc_map_; } + const BCMap& bc_map() const { return bc_map_; } + +private: + BCMap bc_map_; //!< Boundary conditions map linking a boundary condition type + //!< to physical group numbers +}; + +} // namespace openmc + +#endif // OPENMC_FIELD_H diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 0d8189caa1d..a337d248de8 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -147,6 +147,9 @@ class MaterialVolumes { //! Base mesh class //============================================================================== +// Physical groups map type +using PGMap = std::unordered_map>; + class Mesh { public: // Constructors and destructor @@ -193,6 +196,43 @@ class Mesh { virtual void surface_bins_crossed( Position r0, Position r1, const Direction& u, vector& bins) const = 0; + virtual void full_raytracing(Position r0, Position r1, + vector& outward_surface_ids, vector& inward_surface_ids, + vector& bins, vector& segment_lengths) const = 0; + + // TODO - document + virtual int get_last_bin_inside_mesh( + Position r0, Position r1, int bin) const = 0; + + // TODO - document + virtual void randomly_place_on_physical_group( + Position& pa, int& cell, uint64_t* seed, vector physical_groups) = 0; + + // TODO - document + virtual Position normalize_position(const Position& r) + { + fatal_error("Not implemented!"); + } + + //! Retrieve connectivity of a mesh element + // + //! \param[in] element ID + //! \return element connectivity as IDs of the vertices + virtual std::vector connectivity(int id) const = 0; + + //! Distance to the next boundary. + //! If the initial position is outside the mesh, the distance + //! will be from the initial position to the external boundary + //! of the mesh if hit. + // + //! \param[in] current_bin Current bin number + //! \param[in] r Position of the particle + //! \param[in] u Direction of the particle + //! \param[out] bin_next Next bin number + //! \return Distance to the next boundary + virtual double distance_to_next_boundary( + int current_bin, Position r, Direction u, int& bin_next) const = 0; + //! Get bin at a given position in space // //! \param[in] r Position to get bin for @@ -205,6 +245,9 @@ class Mesh { //! Get the number of mesh cell surfaces. virtual int n_surface_bins() const = 0; + //! Get the number of unique vertices. + virtual int n_vertices() const = 0; + int32_t id() const { return id_; } const std::string& name() const { return name_; } @@ -289,6 +332,25 @@ class Mesh { int id_ {-1}; //!< Mesh ID std::string name_; //!< User-specified name int n_dimension_ {-1}; //!< Number of dimensions + + // Physical groups map accessors + PGMap& pg_map() { return pg_map_; } + const PGMap& pg_map() const { return pg_map_; } + + //! Returns the physical group associated with a given face ID. + // + //! \param[in] face_id Face ID + //! \return Physical group + int get_physical_group(int face_id) const; + + //! Returns face IDs corresponding to a given physical group. + // + //! \param[in] group Physical group + //! \return Face IDs + const std::vector& get_face_ids(int group) const; + +private: + PGMap pg_map_; //!< Physical groups map linking physical group to face IDs }; class StructuredMesh : public Mesh { @@ -327,12 +389,36 @@ class StructuredMesh : public Mesh { int n_surface_bins() const override; + int n_vertices() const override + { + fatal_error("Not implemented!"); + }; + void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; + void full_raytracing(Position r0, Position r1, + vector& outward_surface_ids, vector& inward_surface_ids, + vector& bins, vector& segment_lengths) const override; + + double distance_to_next_boundary( + int current_bin, Position r, Direction u, int& bin_next) const override; + + int get_last_bin_inside_mesh( + Position r0, Position r1, int bin) const override; + + void randomly_place_on_physical_group(Position& pa, int& cell, uint64_t* seed, + vector physical_groups) override + {} + + std::vector connectivity(int id) const override + { + fatal_error("Not implemented!"); + }; + //! Determine which cell or surface bins were crossed by a particle // //! \param[in] r0 Previous position of the particle @@ -416,6 +502,17 @@ class StructuredMesh : public Mesh { virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const = 0; + //! Find the closest distance from a point to the mesh boundaries that are + //! aligned with the k direction. The point has to be located outside the + //! mesh. + //! + //! \param[in] k direction index of grid surface + //! \param[in] r position, from where to calculate the distance + //! \param[in] u direction of flight + //! \return distance to the mesh boundary + virtual double distance_to_mesh_boundary_from_outside( + int k, const Position& r, const Direction& u) const = 0; + //! Get a label for the mesh bin std::string bin_label(int bin) const override; @@ -468,6 +565,13 @@ class PeriodicStructuredMesh : public StructuredMesh { return r - origin_; }; + double distance_to_mesh_boundary_from_outside( + int k, const Position& r, const Direction& u) const override + { + fatal_error("Not implemented"); + return -1.0; + } + // Data members Position origin_ {0.0, 0.0, 0.0}; //!< Origin of the mesh }; @@ -484,6 +588,8 @@ class RegularMesh : public StructuredMesh { RegularMesh(hid_t group); // Overridden methods + int n_vertices() const override; + int get_index_in_direction(double r, int i) const override; virtual std::string get_mesh_type() const override; @@ -493,6 +599,9 @@ class RegularMesh : public StructuredMesh { MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; + double distance_to_mesh_boundary_from_outside( + int k, const Position& r, const Direction& u) const override; + std::pair, vector> plot( Position plot_ll, Position plot_ur) const override; @@ -523,6 +632,20 @@ class RegularMesh : public StructuredMesh { int set_grid(); + double face_area(int face_id); + + void randomly_place_on_physical_group(Position& pa, int& cell, uint64_t* seed, + vector physical_groups) override; + + // face_id is on the ID system where the inward flag is not taken into account + Position sample_on_face(int face_id, uint64_t* seed); + + int return_vertex_unique_id(MeshIndex ijk, int local_vertex_idx) const; + + std::vector connectivity(int id) const override; + + Position normalize_position(const Position& r) override; + // Data members double volume_frac_; //!< Volume fraction of each mesh element double element_volume_; //!< Volume of each mesh element @@ -546,6 +669,9 @@ class RectilinearMesh : public StructuredMesh { MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; + double distance_to_mesh_boundary_from_outside( + int k, const Position& r, const Direction& u) const override; + std::pair, vector> plot( Position plot_ll, Position plot_ur) const override; @@ -715,6 +841,9 @@ class UnstructuredMesh : public Mesh { UnstructuredMesh(pugi::xml_node node); UnstructuredMesh(hid_t group); + double distance_to_next_boundary( + int current_bin, Position r, Direction u, int& bin_next) const override; + static const std::string mesh_type; virtual std::string get_mesh_type() const override; @@ -723,6 +852,17 @@ class UnstructuredMesh : public Mesh { void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; + void full_raytracing(Position r0, Position r1, + vector& outward_surface_ids, vector& inward_surface_ids, + vector& bins, vector& segment_lengths) const override; + + int get_last_bin_inside_mesh( + Position r0, Position r1, int bin) const override; + + void randomly_place_on_physical_group(Position& pa, int& cell, uint64_t* seed, + vector physical_groups) override + {} + void to_hdf5_inner(hid_t group) const override; std::string bin_label(int bin) const override; @@ -751,23 +891,12 @@ class UnstructuredMesh : public Mesh { //! \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 connectivity(int id) const = 0; - //! Get the library used for this unstructured mesh virtual std::string library() const = 0; diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 5b632bbce53..8496f1e3074 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -265,6 +265,24 @@ class BoundaryInfo { 0, 0, 0}; //!< which way lattice indices will change }; +struct TransportEvent { + int event_type = EVENT_UNDEFINED; //!< Type of transport event (crossing + //!< surface, collision, time cutoff) + bool cross_surface_geometry = + false; //!< if a surface of the geometry has been crossed during the event + bool cross_surface_temperature_field = + false; //!< if a surface of the temperature field has been crossed during + //!< the event + + // resets all information + void clear() + { + event_type = EVENT_UNDEFINED; + cross_surface_geometry = false; + cross_surface_temperature_field = false; + } +}; + /* * Contains all geometry state information for a particle. */ @@ -412,6 +430,12 @@ class GeometryState { const double& density_mult() const { return density_mult_; } double& density_mult_last() { return density_mult_last_; } + // Temperature field related information + int& tf_bin() { return tf_bin_; } + const int& tf_bin() const { return tf_bin_; } + int& tf_bin_next() { return tf_bin_next_; } + const int& tf_bin_next() const { return tf_bin_next_; } + private: int64_t id_ {-1}; //!< Unique ID @@ -437,12 +461,15 @@ class GeometryState { int material_ {-1}; //!< index for current material int material_last_ {-1}; //!< index for last material - double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV - double sqrtkT_last_ {0.0}; //!< last temperature + double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV + double sqrtkT_last_ {-1.0}; //!< last temperature double density_mult_ {1.0}; //!< density multiplier double density_mult_last_ {1.0}; //!< last density multiplier + int tf_bin_ = C_NONE; //!< Current temperature field bin + int tf_bin_next_ = C_NONE; //!< Next temperature field bin + #ifdef OPENMC_DAGMC_ENABLED moab::DagMC::RayHistory history_; Direction last_dir_; @@ -568,6 +595,8 @@ class ParticleData : public GeometryState { int64_t n_progeny_ {0}; + TransportEvent next_event_ = TransportEvent(); + public: //---------------------------------------------------------------------------- // Constructors @@ -784,6 +813,10 @@ class ParticleData : public GeometryState { d = 0; } } + + // Next event in transport + TransportEvent& next_event() { return next_event_; } + const TransportEvent& next_event() const { return next_event_; } }; } // namespace openmc diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 0914a0958b8..66028b82293 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -67,10 +67,12 @@ extern bool create_delayed_neutrons; //!< create delayed fission neutrons? extern "C" bool cmfd_run; //!< is a CMFD run? extern bool delayed_photon_scaling; //!< Scale fission photon yield to include delayed +extern bool dnp_drift_on; //!< Transport delayed neutron precursors? extern "C" bool entropy_on; //!< calculate Shannon entropy? extern "C" bool - event_based; //!< use event-based mode (instead of history-based) -extern bool ifp_on; //!< Use IFP for kinetics parameters? + event_based; //!< use event-based mode (instead of history-based) +extern bool temperature_field_on; //!< Is there a temperature field defined? +extern bool ifp_on; //!< Use IFP for kinetics parameters? extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular? extern bool material_cell_offsets; //!< create material cells offsets? extern "C" bool output_summary; //!< write summary.h5? @@ -203,6 +205,13 @@ extern "C" int verbosity; //!< How verbose to make output extern double weight_cutoff; //!< Weight cutoff for Russian roulette extern double weight_survive; //!< Survival weight after Russian roulette +extern double + dnp_drift_external_travel_time; //!< External travel time outside the modeled + //!< system for delayed neutron precursor + //!< drift +extern bool dnp_drift_recycling_on; //!< Is particle recycling on for delayed + //!< neutron precursor drift? + } // namespace settings //============================================================================== diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b2131..c98536a3c1c 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -4,8 +4,10 @@ #ifndef OPENMC_SIMULATION_H #define OPENMC_SIMULATION_H +#include "openmc/field.h" #include "openmc/mesh.h" #include "openmc/particle.h" +#include "openmc/streamline_integrator.h" #include "openmc/vector.h" #include @@ -46,6 +48,10 @@ extern int64_t work_per_rank; //!< number of particles per MPI rank extern const RegularMesh* entropy_mesh; extern const RegularMesh* ufs_mesh; +extern TemperatureField temperature_field; +extern VelocityField velocity_field; +extern StreamlineIntegrator* streamline_integrator; + extern vector k_generation; extern vector work_index; diff --git a/include/openmc/streamline_integrator.h b/include/openmc/streamline_integrator.h new file mode 100644 index 00000000000..5d465d2bd6f --- /dev/null +++ b/include/openmc/streamline_integrator.h @@ -0,0 +1,46 @@ +#ifndef OPENMC_STREAMLINE_INTEGRATOR_H +#define OPENMC_STREAMLINE_INTEGRATOR_H + +#include "openmc/field.h" +#include "openmc/particle_data.h" + +#include + +namespace openmc { + +//! Streamline integrator +class StreamlineIntegrator { +public: + //! Advance to the next integration step. + //! + //! \param [inout] dt Time step + //! \param [inout] y_n Position + //! \param [in] cell_n Cell containing y_n + //! \param [in] field Pointer to the velocity field + virtual void next_step( + double& t_n, Position& y_n, int cell_n, VelocityField& field) = 0; + + //! Time step accessors + double& dt() { return dt_; } + const double& dt() const { return dt_; } + +private: + double dt_; //!< Time step +}; + +//! Runge Kutta 4 integrator +class RK4StreamlineIntegrator : public StreamlineIntegrator { +public: + //! Constructor + RK4StreamlineIntegrator(double time_step) { dt() = time_step; } + + void next_step( + double& t_n, Position& y_n, int cell_n, VelocityField& field) override; + +private: + std::string method_ = "Runge Kutta 4"; //!< Method name +}; + +} // namespace openmc + +#endif // OPENMC_STREAMLINE_INTEGRATOR_H diff --git a/include/openmc/timer.h b/include/openmc/timer.h index d928aad4560..d4d63ed07fc 100644 --- a/include/openmc/timer.h +++ b/include/openmc/timer.h @@ -30,6 +30,7 @@ extern Timer time_event_calculate_xs; extern Timer time_event_advance_particle; extern Timer time_event_surface_crossing; extern Timer time_event_collision; +extern Timer time_event_temperature_mesh_crossing; extern Timer time_event_death; extern Timer time_update_src; diff --git a/openmc/__init__.py b/openmc/__init__.py index c204929c840..6aba68878be 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -19,6 +19,7 @@ from openmc.source import * from openmc.settings import * from openmc.lattice import * +from openmc.field import * from openmc.filter import * from openmc.filter_expansion import * from openmc.trigger import * diff --git a/openmc/field.py b/openmc/field.py new file mode 100644 index 00000000000..a33d32c695d --- /dev/null +++ b/openmc/field.py @@ -0,0 +1,59 @@ +from abc import ABC, abstractmethod + +import openmc + + +class ScalarField(ABC): + """Scalar field defined on a geometric mesh. + + Attributes + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell + + """ + def __init__(self, mesh, values): + """Initialization. + + Parameters + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell + + """ + self.mesh = mesh + self.values = values + + def __eq__(self, other): + if not isinstance(other, ScalarField): + return NotImplementedError() + return self.mesh.id == other.mesh.id and self.values == other.values + + +class TemperatureField(ScalarField): + """Temperature field. + + Attributes + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell + + """ + def __init__(self, mesh, values): + """Initialization. + + Parameters + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell + + """ + super().__init__(mesh, values) diff --git a/openmc/settings.py b/openmc/settings.py index 6919afca4cc..247e99080c4 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -4,6 +4,7 @@ from math import ceil from numbers import Integral, Real from pathlib import Path +import textwrap import traceback import lxml.etree as ET @@ -15,6 +16,7 @@ from ._xml import clean_indentation, get_elem_list, get_text from .mesh import _read_meshes, RegularMesh, MeshBase from .source import SourceBase, MeshSource, IndependentSource +from .field import TemperatureField from .utility_funcs import input_path from .volume import VolumeCalculation from .weight_windows import WeightWindows, WeightWindowGenerator, WeightWindowsList @@ -328,6 +330,10 @@ class Settings: sections be loaded at all temperatures within the range. 'multipole' is a boolean indicating whether or not the windowed multipole method should be used to evaluate resolved resonance cross sections. + temperature_field : openmc.TemperatureField + Temperature field based on a geometric mesh used to specify temperatures + in a model. Temperatures declared from a mesh take precedence over + all other temperature definition (cell, material, and global). trace : tuple or list Show detailed information about a single particle, indicated by three integers: the batch number, generation number, and particle number @@ -424,6 +430,9 @@ def __init__(self, **kwargs): # Shannon entropy mesh self._entropy_mesh = None + # Temperature field + self._temperature_field = None + # Trigger subelement self._trigger_active = None self._trigger_max_batches = None @@ -776,6 +785,15 @@ def entropy_mesh(self) -> RegularMesh: def entropy_mesh(self, entropy: RegularMesh): cv.check_type('entropy mesh', entropy, RegularMesh) self._entropy_mesh = entropy + + @property + def temperature_field(self) -> TemperatureField: + return self._temperature_field + + @temperature_field.setter + def temperature_field(self, temperature_field: TemperatureField): + cv.check_type('temperature field', temperature_field, TemperatureField) + self._temperature_field = temperature_field @property def trigger_active(self) -> bool: @@ -1744,6 +1762,33 @@ def _create_entropy_mesh_subelement(self, root, mesh_memo=None): if mesh_memo is not None: mesh_memo.add(self.entropy_mesh.id) + def _create_temperature_field_subelement(self, root, mesh_memo=None): + if self.temperature_field is None: + return + + # add mesh ID to this element + element = ET.SubElement(root, "temperature_field") + subelement = ET.SubElement(element, "mesh") + subelement.text = str(self.temperature_field.mesh.id) + + # If this mesh has already been written outside the + # settings element, skip writing it again + if mesh_memo and self.temperature_field.mesh.id in mesh_memo: + return + + # See if a element already exists -- if not, add it + path = f"./mesh[@id='{self.temperature_field.mesh.id}']" + if root.find(path) is None: + root.append(self.temperature_field.mesh.to_xml_element()) + if mesh_memo is not None: + mesh_memo.add(self.temperature_field.mesh.id) + + # Add temperature values + subelement = ET.SubElement(element, "values") + subelement.text = "\n ".join( + textwrap.wrap(" ".join( + [str(i) for i in self.temperature_field.values]), 80)) + def _create_trigger_subelement(self, root): if self._trigger_active is not None: trigger_element = ET.SubElement(root, "trigger") @@ -2254,6 +2299,16 @@ def _entropy_mesh_from_xml_element(self, root, meshes): raise ValueError(f'Could not locate mesh with ID "{mesh_id}"') self.entropy_mesh = meshes[mesh_id] + def _temperature_field_from_xml_element(self, root, meshes): + elem = root.find('temperature_field') + if elem is not None: + mesh_id = int(get_text(elem, 'mesh')) + if mesh_id not in meshes: + raise ValueError(f'Could not locate mesh with ID "{mesh_id}"') + mesh = meshes[mesh_id] + values = [float(x) for x in get_text(elem, 'values').split()] + self.temperature_field = TemperatureField(mesh, values) + def _trigger_from_xml_element(self, root): elem = root.find('trigger') if elem is not None: @@ -2543,6 +2598,7 @@ def to_xml_element(self, mesh_memo=None): self._create_survival_biasing_subelement(element) self._create_cutoff_subelement(element) self._create_entropy_mesh_subelement(element, mesh_memo) + self._create_temperature_field_subelement(element, mesh_memo) self._create_trigger_subelement(element) self._create_no_reduce_subelement(element) self._create_verbosity_subelement(element) @@ -2661,6 +2717,7 @@ def from_xml_element(cls, elem, meshes=None): settings._survival_biasing_from_xml_element(elem) settings._cutoff_from_xml_element(elem) settings._entropy_mesh_from_xml_element(elem, meshes) + settings._temperature_field_from_xml_element(elem, meshes) settings._trigger_from_xml_element(elem) settings._no_reduce_from_xml_element(elem) settings._verbosity_from_xml_element(elem) diff --git a/src/dnp_drift.cpp b/src/dnp_drift.cpp new file mode 100644 index 00000000000..94fc13b03fa --- /dev/null +++ b/src/dnp_drift.cpp @@ -0,0 +1,194 @@ +#include "openmc/dnp_drift.h" +#include "openmc/error.h" +#include "openmc/field.h" +#include "openmc/particle_data.h" +#include "openmc/position.h" +#include "openmc/simulation.h" + +namespace openmc { + +void _adjust_position(Position& y_n, const Position& y_n_minus_1, double t, + double dt, double decay_time) +{ + if (dt <= 0.0) { + fatal_error("The time step dt should be greater than 0."); + } + if (decay_time > t) { + fatal_error("Decay time should be lower or equal to the total time t."); + } + + double excess_time = t - decay_time; + + if (excess_time > dt) { + fatal_error("The excess time cannot be greater than the time step dt."); + } + + double ab = (y_n - y_n_minus_1).norm(); + + if (ab <= 0.0) { + fatal_error( + "The distance between the two points cannot be lower or equal to 0."); + } + + double ac = (dt - excess_time) * ab / dt; + + y_n = y_n_minus_1 + (y_n - y_n_minus_1) * ac / ab; +} + +void _adjust_time(double& t, double ta, const Position& pa, const Position& pb, + const Position& pc) +{ + if (t <= ta) { + fatal_error( + "Time t at point B must be greater than tine ta at point A."); + } + + double ab = (pb - pa).norm(); + double ac = (pc - pa).norm(); + double cb = (pb - pc).norm(); + + if (ab - (ac + cb) > DNP_DRIFT_DISTANCE_MIN) { + fatal_error( + "Point C must be located between point A and point B."); + } + + t -= (t - ta) * cb / ab; +} + +bool transport_dnp(SourceSite& site, double decay_time, uint64_t* seed) +{ + // Initialization + double t_n = 0.0; // Current (n) integration time + double t_n_minus_1; // Previous (n-1) integration time + Position& y_n = site.r; // Current (n) location of the site + Position y_n_minus_1; // Previous (n-1) location of the site + int cell_n; // Current (n) cell + int cell_n_minus_1; // Previous (n-1) cell + int n_iteration = 0; // Total number of performed iterations + double t_before_decay = decay_time; // Time remaining before decay occurs + BCType crossed_boundary = + BCType::NONE; // Boundary type of the last crossed surface + Position intersection; // Intersection with the mesh boundary + + // Localize initial position + cell_n = simulation::velocity_field.get_mesh_bin(y_n); + if (cell_n < 0) { + // Particle can be close to the wall of the mesh so it can be stopped at + // the initial location because we generally put no slip boundary + // conditions on walls which means that the particle would not have been + // able to move even if detected inside the mesh + t_before_decay = 0.; + return true; + } + + // Transport + while (t_n < decay_time) { + + // Increment number of iterations + n_iteration++; + + // Verify that we have not reached the maximum number of iterations + if (n_iteration > DNP_DRIFT_TRANSPORT_MAX_ITER) { + fatal_error("Maximum number of iterations reached during DNP transport!"); + } + + // Save previous states before integration + t_n_minus_1 = t_n; + y_n_minus_1 = y_n; + cell_n_minus_1 = cell_n; + + // Integration + simulation::streamline_integrator->next_step( + t_n, y_n, cell_n, simulation::velocity_field); + + // Find the next cell + cell_n = simulation::velocity_field.get_next_bin( + y_n_minus_1, y_n, cell_n, crossed_boundary, intersection); + + // If the distance between two consecutive points is low, we block the point + // in position + if ((y_n - y_n_minus_1).norm() < DNP_DRIFT_DISTANCE_MIN) { + t_before_decay = 0.; + return true; + } + + // If the particle left the mesh + if (cell_n < 0) { + + // Update time and position + _adjust_time(t_n, t_n_minus_1, y_n_minus_1, y_n, intersection); + y_n = intersection; + + // If decay time occurs before intersection, stop when decay time is + // reached + if (t_n > decay_time) { + _adjust_position(y_n, y_n_minus_1, t_n, t_n - t_n_minus_1, decay_time); + t_n = decay_time; + t_before_decay = 0.; + return true; + } + + // Determine the action to perform based on the crossed boundary + Actions action; + + switch (crossed_boundary) { + case BCType::OUTLET: + if (settings::dnp_drift_recycling_on) { + action = Actions::PLACE_AT_INLET; + } else { + action = Actions::BLOCK_AT_OUTLET; + } + break; + case BCType::INLET: + action = Actions::PLACE_AT_INLET; + break; + case BCType::WALL: + action = Actions::BLOCK_AT_LOCATION; + break; + default: + fatal_error("Not implemented for this boundary type!"); + break; + } + + // Perform action + switch (action) { + case Actions::PLACE_AT_INLET: + if (decay_time > t_n + settings::dnp_drift_external_travel_time) { + // The DNP has time to reenter the modeled part of the system + simulation::velocity_field.randomly_place_on_inlet(y_n, cell_n, seed); + t_n += settings::dnp_drift_external_travel_time; + cell_n_minus_1 = -1; + y_n_minus_1 = Position(); + } else { + // The DNP decays outside the modeled part of the system + t_before_decay = 0.; + return false; + } + break; + case Actions::BLOCK_AT_OUTLET: + t_before_decay = 0.; + return false; // TODO - return true normally + break; + case Actions::BLOCK_AT_LOCATION: + t_before_decay = decay_time - t_n; + return false; // TODO - return true normally + break; + default: + fatal_error("Unrecognized action in DNP transport!"); + break; + } + } + } + + // Adjust time and position to fit exactly + if (t_n > decay_time) { + _adjust_position(y_n, y_n_minus_1, t_n, + simulation::streamline_integrator->dt(), decay_time); + t_n = decay_time; + } + + t_before_decay = 0; + return true; +} + +} // namespace openmc diff --git a/src/event.cpp b/src/event.cpp index f33e132d0af..9aea0deefde 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -115,10 +115,21 @@ void process_advance_particle_events() p.event_advance(); if (!p.alive()) continue; - if (p.collision_distance() > p.boundary().distance()) { + + switch (p.next_event().event_type) { + case EVENT_CROSS_SURFACE: simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx}); - } else { + break; + case EVENT_COLLIDE: simulation::collision_queue.thread_safe_append({p, buffer_idx}); + break; + case EVENT_TIME_CUTOFF: + p.wgt() = 0.0; + break; + default: + fatal_error(fmt::format("Unknown event '{}' in event-based transport!", + p.next_event().event_type)); + break; } } diff --git a/src/field.cpp b/src/field.cpp new file mode 100644 index 00000000000..662aa24165a --- /dev/null +++ b/src/field.cpp @@ -0,0 +1,146 @@ +#include "openmc/field.h" +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/mesh.h" +#include "openmc/simulation.h" +#include "openmc/vector.h" + +namespace openmc { + +// ----------------------------------------------------------- +// TemperatureField implementation +// ----------------------------------------------------------- + +TemperatureField::TemperatureField( + Mesh* mesh_ptr, std::vector values, std::string mapping) +{ + set_mesh(mesh_ptr); + set_mapping(mapping); + + std::unique_ptr> data = + std::make_unique>(values); + set_data(std::move(data)); +} + +double TemperatureField::get_temperature(int bin) +{ + if (bin >= 0 && bin < values().size()) { + return this->value(bin); + } + return -1.0; +} + +double TemperatureField::get_sqrtkT(int bin) +{ + double temperature = get_temperature(bin); + if (temperature >= 0) { + return sqrt(temperature * K_BOLTZMANN); + } + return -1.0; +} + +// ----------------------------------------------------------- +// VelocityField implementation +// ----------------------------------------------------------- + +VelocityField::VelocityField( + Mesh* mesh_ptr, std::vector values, std::string mapping) +{ + set_mesh(mesh_ptr); + set_mapping(mapping); + + std::unique_ptr> data = + std::make_unique>(values); + set_data(std::move(data)); +} + +int VelocityField::get_next_bin(const Position& r0, const Position& r1, + int bin0, BCType& crossed_boundary, Position& intersection) +{ + // TODO - implement early mesh departure detection + + // Initialization + double total_length = 0.0; + vector outward_surface_ids; + vector inward_surface_ids; + vector bins; + vector segment_lengths; + + // Raytracing + mesh_ptr()->full_raytracing( + r0, r1, outward_surface_ids, inward_surface_ids, bins, segment_lengths); + + // Consistency check + if (segment_lengths.size() == 0) { + fatal_error("Inconsistency in raytrace results."); + } + + // If next point outside the mesh + if (outward_surface_ids.size() != inward_surface_ids.size()) { + + // Retrieve ID from the last surface + int surface_id = outward_surface_ids.back(); + + // Calculate total length travelled from r0 to intersection + for (auto l : segment_lengths) { + total_length += l; + } + + // Calculate intersection + intersection = r0 + (r1 - r0) * total_length; + + // Translate surface ID in physical group + int physical_group = mesh_ptr()->get_physical_group(surface_id); + + // Determine crossed_boundary type + crossed_boundary = get_boundary_condition(physical_group); + + // Return bin default value when outside the mesh + return C_NONE; + } + + // Next point inside the mesh + intersection = r1; + crossed_boundary = BCType::NONE; + return bins.back(); +} + +void VelocityField::randomly_place_on_inlet( + Position& p, int& bin, uint64_t* seed) +{ + mesh_ptr()->randomly_place_on_physical_group( + p, bin, seed, bc_map_[BCType::INLET]); +} + +BCType VelocityField::get_boundary_condition(int physical_group) +{ + if (!bc_map_.empty()) { + for (auto& pair : bc_map_) { + std::vector& vec = pair.second; + if (std::find(vec.begin(), vec.end(), physical_group) != vec.end()) { + return pair.first; + } + } + fatal_error("Not found!"); + } else { + fatal_error("Empty map for boundary conditions in the velocity field!"); + } +} + +//============================================================================== +// C API +//============================================================================== + +extern "C" int openmc_temperature_field_set_temperature( + int32_t index, double temperature) +{ + if (index < 0 || index >= simulation::temperature_field.values().size()) { + set_errmsg("Index in temperature field is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + simulation::temperature_field.assign(index, temperature); + return 0; +} + +} // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 98f1125347d..b43370b58a3 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -135,6 +135,7 @@ int openmc_finalize() settings::ssw_max_particles = 0; settings::ssw_max_files = 1; settings::survival_biasing = false; + settings::temperature_field_on = false; settings::temperature_default = 293.6; settings::temperature_method = TemperatureMethod::NEAREST; settings::temperature_multipole = false; @@ -163,6 +164,8 @@ int openmc_finalize() simulation::entropy_mesh = nullptr; simulation::ufs_mesh = nullptr; + simulation::temperature_field = TemperatureField(); + data::energy_max = {INFTY, INFTY, INFTY, INFTY}; data::energy_min = {0.0, 0.0, 0.0, 0.0}; data::temperature_min = 0.0; diff --git a/src/geometry.cpp b/src/geometry.cpp index ddb61385f18..04c0ddc67d1 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -172,11 +172,19 @@ bool find_cell_inner( p.cell_instance() = cell_instance_at_level(p, p.n_coord() - 1); } - // Set the material, temperature and density multiplier. + // Set the material. p.material_last() = p.material(); p.material() = c.material(p.cell_instance()); + + // Set the temperature. p.sqrtkT_last() = p.sqrtkT(); - p.sqrtkT() = c.sqrtkT(p.cell_instance()); + if (settings::temperature_field_on && p.tf_bin() != C_NONE) { + p.sqrtkT() = simulation::temperature_field.get_sqrtkT(p.tf_bin()); + } else { + p.sqrtkT() = c.sqrtkT(p.cell_instance()); + } + + // Set the density multiplier. p.density_mult_last() = p.density_mult(); p.density_mult() = c.density_mult(p.cell_instance()); diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index a740740c1e6..e88d21d4808 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -17,6 +17,7 @@ #include "openmc/lattice.h" #include "openmc/material.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/surface.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_cell_instance.h" @@ -266,6 +267,28 @@ void get_temperatures( } } } + + // Add temperatures from a temperature field. + // We assume that we do not know how geometric cells are impacted by the + // temperature field in advance. If we had access to this information, we + // could limit the declarations of temperature from the temperature field to + // impacted nuclides only. + if (settings::temperature_field_on) { + for (auto t : simulation::temperature_field.values()) { + // Nuclide temperatures + for (size_t i = 0; i < nuc_temps.size(); i++) { + if (!contains(nuc_temps[i], t)) { + nuc_temps[i].push_back(t); + } + } + // Thermal scattering temperatures + for (size_t i = 0; i < thermal_temps.size(); i++) { + if (!contains(thermal_temps[i], t)) { + thermal_temps[i].push_back(t); + } + } + } + } } //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index 5ab7ac3988b..cb0f380617a 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -750,6 +750,35 @@ void Mesh::to_hdf5(hid_t group) const close_group(mesh_group); } +int Mesh::get_physical_group(int face_id) const +{ + if (!pg_map().empty()) { + for (auto& pair : pg_map()) { + const std::vector& vec = pair.second; + if (std::find(vec.begin(), vec.end(), face_id) != vec.end()) { + return pair.first; + } + } + fatal_error("Not found!"); + } else { + fatal_error("Empty map for physical groups!"); + } +} + +const std::vector& Mesh::get_face_ids(int group) const +{ + if (!pg_map().empty()) { + auto it = pg_map().find(group); + if (it != pg_map().end()) { + return it->second; + } else { + fatal_error("Key not found!"); + } + } else { + fatal_error("Empty map for physical groups!"); + } +} + //============================================================================== // Structured Mesh implementation //============================================================================== @@ -872,6 +901,13 @@ UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group) } } +double UnstructuredMesh::distance_to_next_boundary( + int current_bin, Position r, Direction u, int& bin_next) const +{ + fatal_error("Not implemented"); + return -1.0; +} + void UnstructuredMesh::determine_bounds() { double xmin = INFTY; @@ -937,6 +973,19 @@ void UnstructuredMesh::surface_bins_crossed( fatal_error("Unstructured mesh surface tallies are not implemented."); } +void UnstructuredMesh::full_raytracing(Position r0, Position r1, + vector& outward_surface_ids, vector& inward_surface_ids, + vector& bins, vector& segment_lengths) const +{ + fatal_error("Not implemented."); +} + +int UnstructuredMesh::get_last_bin_inside_mesh( + Position r0, Position r1, int bin) const +{ + fatal_error("Not implemented."); +} + std::string UnstructuredMesh::bin_label(int bin) const { return fmt::format("Mesh Index ({})", bin); @@ -1321,6 +1370,150 @@ void StructuredMesh::surface_bins_crossed( raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins)); } +// Get next bin using raytracing +void StructuredMesh::full_raytracing(Position r0, Position r1, + vector& outward_surface_ids, vector& inward_surface_ids, + vector& bins, vector& segment_lengths) const +{ + // TODO: can reconstruct bins and surface_ids instead of explicitly tracking + // them + + // Helper tally class. + struct CompleteAggregator { + CompleteAggregator(const StructuredMesh* _mesh, + vector& _outward_surface_ids, vector& _inward_surface_ids, + vector& _bins, vector& _lengths) + : mesh(_mesh), outward_surface_ids(_outward_surface_ids), + inward_surface_ids(_inward_surface_ids), bins(_bins), lengths(_lengths) + {} + // Returns surface ID without the inward information + // This is different from the other representation with the inward + // information + void surface(const MeshIndex& ijk, int k, bool max, bool inward) const + { + int surface_id = + 2 * mesh->n_dimension_ * mesh->get_bin_from_indices(ijk) + 2 * k; + if (max) + surface_id += 1; + + if (inward) { + inward_surface_ids.push_back(surface_id); + } else { + outward_surface_ids.push_back(surface_id); + } + } + void track(const MeshIndex& ijk, double l) const + { + bins.push_back(mesh->get_bin_from_indices(ijk)); + lengths.push_back(l); + } + + const StructuredMesh* mesh; + vector& outward_surface_ids; + vector& inward_surface_ids; + vector& bins; + vector& lengths; + }; + + // Determine direction + Direction u = (r1 - r0) / (r1 - r0).norm(); + + // Perform the mesh raytrace with the helper class. + raytrace_mesh(r0, r1, u, + CompleteAggregator( + this, outward_surface_ids, inward_surface_ids, bins, segment_lengths)); +} + +double StructuredMesh::distance_to_next_boundary( + int current_bin, Position r, Direction u, int& bin_next) const +{ + Position global_r = r; + Position local_r = local_coords(r); + + double distance; + bool in_mesh; + + // Find the cell indices + StructuredMesh::MeshIndex ijk; + if (current_bin >= 0) { + + ijk = get_indices_from_bin(current_bin); + + // Calculate initial distances to next surfaces in all three dimensions + std::array distances; + for (int k = 0; k < n_dimension_; ++k) { + distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0); + } + + // Find next ijk + auto k_min = + std::min_element(distances.begin(), distances.end()) - distances.begin(); + distance = distances[k_min].distance; + ijk[k_min] = distances[k_min].next_index; + + // If the particle is on a surface, test using the next index + if (distance <= FP_COINCIDENT) { + + // Update distances + for (int k = 0; k < n_dimension_; ++k) { + distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0); + } + + k_min = std::min_element(distances.begin(), distances.end()) - + distances.begin(); + distance = distances[k_min].distance; + ijk[k_min] = distances[k_min].next_index; + } + + // Determine next bin + in_mesh = true; + for (int k = 0; k < n_dimension_; ++k) { + if ((ijk[k] < 1) || (ijk[k] > shape_[k])) { + in_mesh = false; + } + } + if (in_mesh) { + bin_next = get_bin_from_indices(ijk); + } else { + bin_next = C_NONE; + } + + } else { // Outside mesh + + // Calculate distance to mesh from outside + distance = INFTY; + for (int k = 0; k < n_dimension_; k++) { + double d = distance_to_mesh_boundary_from_outside(k, r, u); + if (d < INFTY) { + distance = d; + } + } + + // Determine next bin + if (distance < INFTY) { + ijk = get_indices(global_r + (distance + TINY_BIT) * u, in_mesh); + bin_next = get_bin_from_indices(ijk); + } else { + bin_next = C_NONE; + } + } + + return distance; +} + +int StructuredMesh::get_last_bin_inside_mesh( + Position r0, Position r1, int bin) const +{ + vector bins; + vector lengths; + bins_crossed(r0, r1, (r1 - r0) / (r1 - r0).norm(), bins, lengths); + if (!bins.empty()) { + return bins.back(); + } else { + return bin; + } +} + //============================================================================== // RegularMesh implementation //============================================================================== @@ -1470,9 +1663,22 @@ RegularMesh::RegularMesh(hid_t group) : StructuredMesh {group} } } +int RegularMesh::n_vertices() const +{ + int n_vertices = 1; + for (int i = 0; i < n_dimension_; i++) { + n_vertices *= shape_[i] + 1; + } + return n_vertices; +} + int RegularMesh::get_index_in_direction(double r, int i) const { - return std::ceil((r - lower_left_[i]) / width_[i]); + if (r == lower_left_[i]) { + return 1; + } else { + return std::ceil((r - lower_left_[i]) / width_[i]); + } } const std::string RegularMesh::mesh_type = "regular"; @@ -1513,6 +1719,34 @@ StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary( return d; } +double RegularMesh::distance_to_mesh_boundary_from_outside( + int k, const Position& r, const Direction& u) const +{ + double t; + double distance = INFTY; + + if (u[k] > 0.0) { + t = (lower_left_[k] - r[k]) / u[k]; + } else { + t = (upper_right_[k] - r[k]) / u[k]; + } + + if (t > FP_COINCIDENT) { + bool reenter = true; + for (int i = 0; i < n_dimension_; i++) { + if (i != k) { + double a = r[i] + u[i] * t; + reenter &= (a >= lower_left_[i]); + reenter &= (a <= upper_right_[i]); + } + } + if (reenter) { + distance = t; + } + } + return distance; +} + std::pair, vector> RegularMesh::plot( Position plot_ll, Position plot_ur) const { @@ -1616,6 +1850,192 @@ double RegularMesh::volume(const MeshIndex& ijk) const return element_volume_; } +double RegularMesh::face_area(int face_id) +{ + // Only implemented in 3D + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + int face_local_id = face_id % 6; + int normal_id = face_local_id / 2; + + double area = 0.; + if (normal_id == 0) { + area = width_[1] * width_[2]; + } else if (normal_id == 1) { + area = width_[0] * width_[2]; + } else if (normal_id == 2) { + area = width_[0] * width_[1]; + } + + return area; +} + +void RegularMesh::randomly_place_on_physical_group( + Position& pa, int& cell, uint64_t* seed, vector physical_groups) +{ + // Only implemented in 3D + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + // Retrieve the face IDs from the physical groups + vector face_ids; + for (auto g : physical_groups) { + vector temp = get_face_ids(g); + face_ids.insert(face_ids.end(), temp.begin(), temp.end()); + } + + // Retrieve the face area of each face and calculate the total + double total_area = 0.; + vector areas; + for (auto s : face_ids) { + double area = face_area(s); + areas.push_back(area); + total_area += area; + } + + // Randomly select a face (weighted by surface area) + int selected_face_id; + double random_limit = total_area * prn(seed); + double incremental_area = 0.; + for (int i = 0; i < face_ids.size(); i++) { + incremental_area += areas[i]; + if ((random_limit - incremental_area) <= 1.0E-13) { + selected_face_id = face_ids[i]; + break; + } + } + + // Sample a new position on the selected face + pa = sample_on_face(selected_face_id, seed); + + // Returns cell as well + cell = selected_face_id / (2 * n_dimension_); +} + +Position RegularMesh::sample_on_face(int face_id, uint64_t* seed) +{ + // Only implemented in 3D + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + int bin_id = face_id / (2 * n_dimension_); + MeshIndex ijk = get_indices_from_bin(bin_id); + + int face_local_id = face_id % (2 * n_dimension_); + int normal_id = face_local_id / 2; + int max_value = face_local_id % 2; + + Position p = Position(); + + int dimensions[] = {0, 1, 2}; + + for (int i : dimensions) { + if (i == normal_id) { + if (max_value) { + p[i] = positive_grid_boundary(ijk, i); + } else { + p[i] = negative_grid_boundary(ijk, i); + } + } else { + p[i] = negative_grid_boundary(ijk, i) + width_(i) * prn(seed); + } + } + + return p; +} + +int RegularMesh::return_vertex_unique_id( + MeshIndex ijk, int local_vertex_idx) const +{ + // Only implemented in 3D + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + switch (local_vertex_idx) { + case 0: + ijk[0]--; + ijk[1]--; + ijk[2]--; + break; + case 1: + ijk[1]--; + ijk[2]--; + break; + case 2: + ijk[0]--; + ijk[2]--; + break; + case 3: + ijk[2]--; + break; + case 4: + ijk[0]--; + ijk[1]--; + break; + case 5: + ijk[1]--; + break; + case 6: + ijk[0]--; + break; + case 7: + break; + default: + fatal_error("Local vertex index out of bound!"); + } + + return ijk[0] + ijk[1] * (shape_[0] + 1) + + ijk[2] * (shape_[0] + 1) * (shape_[1] + 1); +} + +std::vector RegularMesh::connectivity(int id) const +{ + // Only implemented in 3D + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + MeshIndex ijk = get_indices_from_bin(id); + + vector v; + for (int i = 0; i < 8; i++) { + v.push_back(return_vertex_unique_id(ijk, i)); + } + + return v; +} + +Position RegularMesh::normalize_position(const Position& r) +{ + // Only implemented in 3D + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + int bin = get_bin(r); + MeshIndex ijk = get_indices_from_bin(bin); + + // Retrieve dimensions + double xmin = negative_grid_boundary(ijk, 0); + double xmax = positive_grid_boundary(ijk, 0); + double ymin = negative_grid_boundary(ijk, 1); + double ymax = positive_grid_boundary(ijk, 1); + double zmin = negative_grid_boundary(ijk, 2); + double zmax = positive_grid_boundary(ijk, 2); + + // Normalize and restrict normalized values between 0.0 and 1.0 + double p_x = std::clamp((r[0] - xmin) / (xmax - xmin), 0.0, 1.0); + double p_y = std::clamp((r[1] - ymin) / (ymax - ymin), 0.0, 1.0); + double p_z = std::clamp((r[2] - zmin) / (zmax - zmin), 0.0, 1.0); + + return Position(p_x, p_y, p_z); +} + //============================================================================== // RectilinearMesh implementation //============================================================================== @@ -1685,6 +2105,34 @@ StructuredMesh::MeshDistance RectilinearMesh::distance_to_grid_boundary( return d; } +double RectilinearMesh::distance_to_mesh_boundary_from_outside( + int k, const Position& r, const Direction& u) const +{ + double t; + double distance = INFTY; + + if (u[k] > 0.0) { + t = (lower_left_[k] - r[k]) / u[k]; + } else { + t = (upper_right_[k] - r[k]) / u[k]; + } + + if (t > FP_COINCIDENT) { + bool reenter = true; + for (int i = 0; i < n_dimension_; i++) { + if (i != k) { + double a = r[i] + u[i] * t; + reenter &= (a >= lower_left_[i]); + reenter &= (a <= upper_right_[i]); + } + } + if (reenter) { + distance = t; + } + } + return distance; +} + int RectilinearMesh::set_grid() { shape_ = {static_cast(grid_[0].size()) - 1, diff --git a/src/output.cpp b/src/output.cpp index c66c313dce5..c74e33e9679 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -453,6 +453,10 @@ void print_runtime() show_time("Advancing", time_event_advance_particle.elapsed(), 2); show_time("Surface crossings", time_event_surface_crossing.elapsed(), 2); show_time("Collisions", time_event_collision.elapsed(), 2); + if (settings::temperature_field_on) { + show_time("Temperature mesh crossings", + time_event_temperature_mesh_crossing.elapsed(), 2); + } show_time("Particle death", time_event_death.elapsed(), 2); } if (settings::run_mode == RunMode::EIGENVALUE) { diff --git a/src/particle.cpp b/src/particle.cpp index 11747e0cca3..71a4b0cc9fb 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -128,6 +128,10 @@ void Particle::from_source(const SourceSite* src) fission() = false; zero_flux_derivs(); lifetime() = 0.0; + if (settings::temperature_field_on) { + tf_bin() = C_NONE; + tf_bin_next() = C_NONE; + } #ifdef OPENMC_DAGMC_ENABLED history().reset(); #endif @@ -184,6 +188,12 @@ void Particle::event_calculate_xs() // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles if (lowest_coord().cell() == C_NONE) { + + // Define temperature field cell + if (settings::temperature_field_on) { + tf_bin() = simulation::temperature_field.get_mesh_bin(r()); + } + if (!exhaustive_find_cell(*this)) { mark_as_lost( "Could not find the cell containing particle " + std::to_string(id())); @@ -237,7 +247,7 @@ void Particle::event_calculate_xs() void Particle::event_advance() { - // Find the distance to the nearest boundary + // Find the distance to the nearest geometry boundary boundary() = distance_to_boundary(*this); // Sample a distance to collision @@ -250,14 +260,42 @@ void Particle::event_advance() collision_distance() = -std::log(prn(current_seed())) / macro_xs().total; } + // Find the distance to the nearest temperature mesh cell surface + double distance_tmesh = INFTY; + if (settings::temperature_field_on) { + distance_tmesh = simulation::temperature_field.distance_to_next_boundary( + tf_bin(), r(), u(), tf_bin_next()); + } + + // Calculate the distance corresponding to the time cutoff double speed = this->speed(); double time_cutoff = settings::time_cutoff[type().transport_index()]; double distance_cutoff = (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; - // Select smaller of the three distances + // Determine minimal distance for cross surface events + double distance_cross_surface = + std::min({boundary().distance(), distance_tmesh}); + + // Determine minimal distance of all events double distance = - std::min({boundary().distance(), collision_distance(), distance_cutoff}); + std::min({distance_cross_surface, collision_distance(), distance_cutoff}); + + // Determine next event + next_event().clear(); + if (distance == distance_cutoff) { + next_event().event_type = EVENT_TIME_CUTOFF; + } else { + if (collision_distance() > distance_cross_surface) { + next_event().event_type = EVENT_CROSS_SURFACE; + next_event().cross_surface_geometry = + (std::abs(distance - boundary().distance()) <= FP_COINCIDENT); + next_event().cross_surface_temperature_field = + (std::abs(distance - distance_tmesh) <= FP_COINCIDENT); + } else { + next_event().event_type = EVENT_COLLIDE; + } + } // Advance particle in space and time this->move_distance(distance); @@ -284,71 +322,95 @@ void Particle::event_advance() if (!model::active_tallies.empty()) { score_track_derivative(*this, distance); } - - // Set particle weight to zero if it hit the time boundary - if (distance == distance_cutoff) { - wgt() = 0.0; - } } void Particle::event_cross_surface() { - // Saving previous cell data - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); - - // Set surface that particle is on and adjust coordinate levels - surface() = boundary().surface(); - n_coord() = boundary().coord_level(); - - if (boundary().lattice_translation()[0] != 0 || - boundary().lattice_translation()[1] != 0 || - boundary().lattice_translation()[2] != 0) { - // Particle crosses lattice boundary - - bool verbose = settings::verbosity >= 10 || trace(); - cross_lattice(*this, boundary(), verbose); - event() = TallyEvent::LATTICE; - - // Score cell to cell partial currents - if (!model::active_surface_tallies.empty()) { - auto& lat {*model::lattices[lowest_coord().lattice()]}; - bool is_valid; - Direction normal = - lat.get_normal(boundary().lattice_translation(), is_valid); - if (is_valid) { + if (next_event().cross_surface_geometry) { + + // Saving previous cell data + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); + + // Set surface that particle is on and adjust coordinate levels + surface() = boundary().surface(); + n_coord() = boundary().coord_level(); + + if (boundary().lattice_translation()[0] != 0 || + boundary().lattice_translation()[1] != 0 || + boundary().lattice_translation()[2] != 0) { + // Particle crosses lattice boundary + + // Update temperature field bin + if (settings::temperature_field_on) { + if (next_event().cross_surface_temperature_field) { + tf_bin() = tf_bin_next(); + } + } + + bool verbose = settings::verbosity >= 10 || trace(); + cross_lattice(*this, boundary(), verbose); + event() = TallyEvent::LATTICE; + + // Score cell to cell partial currents + if (!model::active_surface_tallies.empty()) { + auto& lat {*model::lattices[lowest_coord().lattice()]}; + bool is_valid; + Direction normal = + lat.get_normal(boundary().lattice_translation(), is_valid); + if (is_valid) { + normal /= normal.norm(); + score_surface_tally(*this, model::active_surface_tallies, normal); + } + } + + } else { + + const auto& surf {*model::surfaces[surface_index()].get()}; + + // Particle crosses surface + // If BC, add particle to surface source before crossing surface + if (surf.surf_source_ && surf.bc_) { + add_surf_source_to_bank(*this, surf); + } + this->cross_surface(surf); + // If no BC, add particle to surface source after crossing surface + if (surf.surf_source_ && !surf.bc_) { + add_surf_source_to_bank(*this, surf); + } + if (settings::weight_window_checkpoint_surface) { + apply_weight_windows(*this); + } + event() = TallyEvent::SURFACE; + + // Score cell to cell partial currents + if (!model::active_surface_tallies.empty()) { + Direction normal = surf.normal(r()); normal /= normal.norm(); score_surface_tally(*this, model::active_surface_tallies, normal); } } - } else { + // Update particle temperature from the temperature field + } else if (next_event().cross_surface_temperature_field) { - const auto& surf {*model::surfaces[surface_index()].get()}; + sqrtkT_last() = sqrtkT(); - // Particle crosses surface - // If BC, add particle to surface source before crossing surface - if (surf.surf_source_ && surf.bc_) { - add_surf_source_to_bank(*this, surf); - } - this->cross_surface(surf); - // If no BC, add particle to surface source after crossing surface - if (surf.surf_source_ && !surf.bc_) { - add_surf_source_to_bank(*this, surf); - } - if (settings::weight_window_checkpoint_surface) { - apply_weight_windows(*this); - } - event() = TallyEvent::SURFACE; + tf_bin() = tf_bin_next(); - // Score cell to cell partial currents - if (!model::active_surface_tallies.empty()) { - Direction normal = surf.normal(r()); - normal /= normal.norm(); - score_surface_tally(*this, model::active_surface_tallies, normal); + if (tf_bin() != C_NONE) { + sqrtkT() = simulation::temperature_field.get_sqrtkT(tf_bin()); + } else { + int i_cell = lowest_coord().cell(); + Cell& c {*model::cells[i_cell]}; + sqrtkT() = c.sqrtkT(cell_instance()); } + +#ifdef OPENMC_DAGMC_ENABLED + history().reset(); +#endif } } @@ -476,6 +538,12 @@ void Particle::event_revive_from_secondary() // removed from the pulse-height of this cell. if (lowest_coord().cell() == C_NONE) { bool verbose = settings::verbosity >= 10 || trace(); + + // Define temperature field cell + if (settings::temperature_field_on) { + tf_bin() = simulation::temperature_field.get_mesh_bin(r()); + } + if (!exhaustive_find_cell(*this, verbose)) { mark_as_lost("Could not find the cell containing particle " + std::to_string(id())); @@ -594,6 +662,13 @@ void Particle::cross_surface(const Surface& surf) return; } + // Update temperature field bin after handling boundary conditions + if (settings::temperature_field_on) { + if (next_event().cross_surface_temperature_field) { + tf_bin() = tf_bin_next(); + } + } + // ========================================================================== // SEARCH NEIGHBOR LISTS FOR NEXT CELL @@ -616,7 +691,11 @@ void Particle::cross_surface(const Surface& surf) cell_instance() = cell_instance_at_level(*this, n_coord() - 1); material() = cell->material(cell_instance()); - sqrtkT() = cell->sqrtkT(cell_instance()); + if (settings::temperature_field_on && tf_bin() != C_NONE) { + sqrtkT() = simulation::temperature_field.get_sqrtkT(tf_bin()); + } else { + sqrtkT() = cell->sqrtkT(cell_instance()); + } density_mult() = cell->density_mult(cell_instance()); return; } @@ -718,6 +797,8 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u) coord(0).cell() = cell_last(0); surface() = -surface(); + // Particle's temperature field bin is unchanged + // If a reflective surface is coincident with a lattice or universe // boundary, it is necessary to redetermine the particle's coordinates in // the lower universes. @@ -768,6 +849,11 @@ void Particle::cross_periodic_bc( // Reassign particle's surface surface() = new_surface; + // Reassign particle's temperature field bin + if (settings::temperature_field_on) { + tf_bin() = simulation::temperature_field.get_mesh_bin(r()); + } + // Figure out what cell particle is in now n_coord() = 1; diff --git a/src/physics.cpp b/src/physics.cpp index 106bd1aa2b2..f8f1266c864 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -5,6 +5,7 @@ #include "openmc/chain.h" #include "openmc/constants.h" #include "openmc/distribution_multi.h" +#include "openmc/dnp_drift.h" #include "openmc/eigenvalue.h" #include "openmc/endf.h" #include "openmc/error.h" @@ -205,9 +206,9 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Counter for the number of fission sites successfully stored to the shared // fission bank or the secondary particle bank - int n_sites_stored; + int n_sites_stored = 0; - for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { + for (int i = 0; i < nu; i++) { // Initialize fission site object with particle data SourceSite site; site.r = p.r(); @@ -219,8 +220,19 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Sample delayed group and angle/energy for fission reaction sample_fission_neutron(i_nuclide, rx, &site, p); - // Reject site if it exceeds time cutoff + // If a delayed neutron is sampled if (site.delayed_group > 0) { + + // Explicit transport of Delayed Neutron Precursor (DNP) + if (settings::dnp_drift_on) { + double dnp_decay_time = site.time - p.time(); + if (!transport_dnp(site, dnp_decay_time, p.current_seed())) { + continue; + } + // TODO: update site.time accordingly + } + + // Reject site if it exceeds time cutoff double t_cutoff = settings::time_cutoff[site.particle.transport_index()]; if (site.time > t_cutoff) { continue; @@ -257,6 +269,8 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) p.n_secondaries()++; } + n_sites_stored++; + // Increment the number of neutrons born delayed if (site.delayed_group > 0) { nu_d[site.delayed_group - 1]++; diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index dde5023e44f..a87f63d1a8a 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -288,6 +288,11 @@ void RandomRay::event_advance_ray() boundary() = distance_to_boundary(*this); double distance = boundary().distance(); + // Define next event + next_event().clear(); + next_event().event_type = EVENT_CROSS_SURFACE; + next_event().cross_surface_geometry = true; + if (distance < 0.0) { mark_as_lost("Negative transport distance detected for particle " + std::to_string(id())); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 80cbfc3fe5b..da4c5682e7f 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -235,6 +235,13 @@ void validate_random_ray_inputs() "quality FW-CADIS weight windows. We recommend you use flat source mode " "when generating weight windows with an overlaid mesh tally."); } + + // Warn about the presence of a temperature field (not supported) + /////////////////////////////////////////////////////////////////// + if (settings::temperature_field_on) { + warning("Temperature fields are not supported with the random ray solver. " + "It will be ignored during this simulation."); + } } void print_adjoint_header() diff --git a/src/settings.cpp b/src/settings.cpp index ab9f9a5aafa..6ea465a6111 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -19,6 +19,7 @@ #include "openmc/distribution_spatial.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/field.h" #include "openmc/file_utils.h" #include "openmc/mcpl_interface.h" #include "openmc/mesh.h" @@ -53,6 +54,7 @@ bool confidence_intervals {false}; bool create_delayed_neutrons {true}; bool create_fission_neutrons {true}; bool delayed_photon_scaling {true}; +bool dnp_drift_on {false}; bool entropy_on {false}; bool event_based {false}; bool ifp_on {false}; @@ -76,6 +78,7 @@ bool surf_mcpl_write {false}; bool surf_source_read {false}; bool survival_biasing {false}; bool survival_normalization {false}; +bool temperature_field_on {false}; bool temperature_multipole {false}; bool trigger_on {false}; bool trigger_predict {false}; @@ -153,6 +156,9 @@ int verbosity {-1}; double weight_cutoff {0.25}; double weight_survive {1.0}; +double dnp_drift_external_travel_time {0.0}; +bool dnp_drift_recycling_on {false}; + } // namespace settings //============================================================================== @@ -822,6 +828,58 @@ void read_settings_xml(pugi::xml_node root) "it by specifying its ID in an element."); } } + + // Temperature field + if (check_for_node(root, "temperature_field")) { + temperature_field_on = true; + + // Get pointer to temperature_field node + auto node_tf = root.child("temperature_field"); + + // Mesh parameter + Mesh* tf_mesh_ptr; + if (check_for_node(node_tf, "mesh")) { + int temp = std::stoi(get_node_value(node_tf, "mesh")); + if (model::mesh_map.find(temp) == model::mesh_map.end()) { + throw std::runtime_error(fmt::format( + "Mesh {} specified for the temperature field does not exist.", temp)); + } + tf_mesh_ptr = model::meshes[model::mesh_map.at(temp)].get(); + } else { + throw std::runtime_error( + "A mesh must be given for the temperature field."); + } + + // Values parameter + vector tf_values; + if (check_for_node(node_tf, "values")) { + auto temp = get_node_array(node_tf, "values"); + if (temp.size() != tf_mesh_ptr->n_bins()) { + throw std::runtime_error( + "Inconsistency in the temperature field: the number of " + "values must be equal to the number of bins in the mesh."); + } + for (const auto& b : temp) { + tf_values.push_back(b); + } + } else { + throw std::runtime_error( + "Temperature values must be given for the temperature field."); + } + + // Mapping representation + std::string mapping; + if (check_for_node(node_tf, "mapping")) { + mapping = get_node_value(node_tf, "mapping"); + } else { + mapping = "cell"; + } + + // Instantiate the temperature field + simulation::temperature_field = + TemperatureField(tf_mesh_ptr, tf_values, mapping); + } + // Uniform fission source weighting mesh if (check_for_node(root, "ufs_mesh")) { auto temp = std::stoi(get_node_value(root, "ufs_mesh")); @@ -1173,6 +1231,176 @@ void read_settings_xml(pugi::xml_node root) temperature_range[1] = range.at(1); } + // Explicit transport of Delayed Neutron Precursor (DNP) + if (check_for_node(root, "dnp_drift")) { + dnp_drift_on = true; + auto node_dnp_drift = root.child("dnp_drift"); + + // Mesh + Mesh* mesh_ptr; + if (check_for_node(node_dnp_drift, "field_mesh")) { + int temp = std::stoi(get_node_value(node_dnp_drift, "field_mesh")); + if (model::mesh_map.find(temp) == model::mesh_map.end()) { + fatal_error(fmt::format( + "Mesh {} specified for the velocity field does not exist.", temp)); + } + mesh_ptr = model::meshes[model::mesh_map.at(temp)].get(); + } else { + fatal_error("A mesh must be given for the velocity field."); + } + + // Values + vector vf_values; + if (check_for_node(node_dnp_drift, "field_values")) { + auto temp = get_node_array(node_dnp_drift, "field_values"); + if (temp.size() % 3 != 0) { + fatal_error("The number of values must be a multiple of 3."); + } + for (size_t i = 0; i + 2 < temp.size(); i += 3) { + Direction d = Direction(temp[i], temp[i + 1], temp[i + 2]); + vf_values.push_back(d); + } + } else { + fatal_error("Values must be given for the velocity field."); + } + + // Mapping representation + std::string field_mapping; + if (check_for_node(node_dnp_drift, "field_mapping")) { + field_mapping = get_node_value(node_dnp_drift, "field_mapping"); + } else { + fatal_error( + "A mapping representation must be given for the velocity field."); + } + + // Velocity field + simulation::velocity_field = + VelocityField(mesh_ptr, vf_values, field_mapping); + + // Boundary conditions map + if (check_for_node(node_dnp_drift, "boundary_map")) { + BCMap bc_map; + auto node_boundary = node_dnp_drift.child("boundary_map"); + + if (check_for_node(node_boundary, "inlet")) { + bc_map[BCType::INLET] = get_node_array(node_boundary, "inlet"); + } else { + fatal_error("Inlet boundary conditions must be declared."); + } + + if (check_for_node(node_boundary, "outlet")) { + bc_map[BCType::OUTLET] = get_node_array(node_boundary, "outlet"); + } else { + fatal_error("Outlet boundary conditions must be declared."); + } + + if (check_for_node(node_boundary, "wall")) { + bc_map[BCType::WALL] = get_node_array(node_boundary, "wall"); + } else { + fatal_error("Wall boundary conditions must be declared."); + } + + simulation::velocity_field.bc_map() = bc_map; + + } else { + fatal_error("Boundary conditions must be declared."); + } + + // Integrator + if (check_for_node(node_dnp_drift, "integrator")) { + std::string integration_method = + get_node_value(node_dnp_drift, "integrator"); + + // Runge Kutta 4 + if (integration_method == "RK4") { + + // Time step + double dt; + if (!check_for_node(node_dnp_drift, "integrator_dt")) { + fatal_error("The attribute 'integrator_dt' is not declared in the " + "DNP drift settings."); + } else { + dt = std::stod(get_node_value(node_dnp_drift, "integrator_dt")); + } + + // Instantiate integrator + simulation::streamline_integrator = new RK4StreamlineIntegrator(dt); + + // Undefined integration method + } else { + fatal_error( + fmt::format("Integrator '{}' not implemented", integration_method)); + } + } else { + fatal_error("An integrator should be defined in the DNP drift settings."); + } + + // Recycle precursor when reaching an outlet? + if (check_for_node(node_dnp_drift, "recycling")) { + dnp_drift_recycling_on = get_node_value_bool(node_dnp_drift, "recycling"); + if (dnp_drift_recycling_on) { + if (!check_for_node(node_dnp_drift, "external_travel_time")) { + fatal_error("The external travel time is not declared in " + "the DNP drift settings."); + } else { + dnp_drift_external_travel_time = + std::stod(get_node_value(node_dnp_drift, "external_travel_time")); + } + } + } + } + + // Add physical group information to mesh + if (check_for_node(root, "mesh_physical_group")) { + + auto node_physical_group = root.child("mesh_physical_group"); + + // Mesh pointer + Mesh* mesh_ptr; + if (check_for_node(node_physical_group, "mesh")) { + int temp = std::stoi(get_node_value(node_physical_group, "mesh")); + if (model::mesh_map.find(temp) == model::mesh_map.end()) { + fatal_error(fmt::format( + "Mesh {} specified for the physical groups does not exist.", temp)); + } + mesh_ptr = model::meshes[model::mesh_map.at(temp)].get(); + } else { + fatal_error("A mesh must be given for the velocity field."); + } + + // Face IDs + vector face_ids; + if (check_for_node(node_physical_group, "face_ids")) { + face_ids = get_node_array(node_physical_group, "face_ids"); + } else { + fatal_error("Surface IDs must be declared."); + } + + // Physical groups + vector physical_groups; + if (check_for_node(node_physical_group, "physical_groups")) { + physical_groups = + get_node_array(node_physical_group, "physical_groups"); + } else { + fatal_error("Physical_groups must be declared."); + } + + // Check for consistency + if (face_ids.size() != physical_groups.size()) { + fatal_error( + "The lists of face IDs and physical groups must have the same size!"); + } + + // Create the physical group map + PGMap pg_map; + for (size_t i = 0; i < face_ids.size(); i++) { + pg_map[physical_groups[i]].push_back(face_ids[i]); + } + + // Save the map in the mesh + mesh_ptr->pg_map() = pg_map; + } + // Check for tabular_legendre options if (check_for_node(root, "tabular_legendre")) { // Get pointer to tabular_legendre node diff --git a/src/simulation.cpp b/src/simulation.cpp index 4fad196a604..2f42451b0fa 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -3,10 +3,12 @@ #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/collision_track.h" +#include "openmc/constants.h" #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" #include "openmc/event.h" +#include "openmc/field.h" #include "openmc/geometry_aux.h" #include "openmc/ifp.h" #include "openmc/material.h" @@ -19,6 +21,7 @@ #include "openmc/settings.h" #include "openmc/source.h" #include "openmc/state_point.h" +#include "openmc/streamline_integrator.h" #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" @@ -321,6 +324,10 @@ int64_t work_per_rank; const RegularMesh* entropy_mesh {nullptr}; const RegularMesh* ufs_mesh {nullptr}; +TemperatureField temperature_field; +VelocityField velocity_field; +StreamlineIntegrator* streamline_integrator; + vector k_generation; vector work_index; @@ -810,10 +817,21 @@ void transport_history_based_single_particle(Particle& p) p.event_advance(); } if (p.alive()) { - if (p.collision_distance() > p.boundary().distance()) { + switch (p.next_event().event_type) { + case EVENT_CROSS_SURFACE: p.event_cross_surface(); - } else if (p.alive()) { + break; + case EVENT_COLLIDE: p.event_collide(); + break; + case EVENT_TIME_CUTOFF: + p.wgt() = 0.0; + break; + default: + fatal_error( + fmt::format("Unknown event '{}' in history-based transport!", + p.next_event().event_type)); + break; } } p.event_revive_from_secondary(); diff --git a/src/streamline_integrator.cpp b/src/streamline_integrator.cpp new file mode 100644 index 00000000000..c4ba774f454 --- /dev/null +++ b/src/streamline_integrator.cpp @@ -0,0 +1,32 @@ +#include "openmc/streamline_integrator.h" +#include "openmc/field.h" +#include "openmc/particle_data.h" + +namespace openmc { + +void RK4StreamlineIntegrator::next_step( + double& tn, Position& yn, int cell_n, VelocityField& field) +{ + int datacell_id; + + // Calculate k1 + Direction k1 = field.evaluate(yn, cell_n) * dt(); + + // Calculate k2 + Position p2 = yn + k1 / 2.; + Direction k2 = field.evaluate(yn, p2, cell_n) * dt(); + + // Calculate k3 + Position p3 = yn + k2 / 2.; + Direction k3 = field.evaluate(yn, p3, cell_n) * dt(); + + // Calculate k4 + Position p4 = yn + k3; + Direction k4 = field.evaluate(yn, p4, cell_n) * dt(); + + // Step forward + yn += (k1 + k2 * 2 + k3 * 2 + k4) / 6.; + tn += dt(); +} + +} // namespace openmc diff --git a/src/timer.cpp b/src/timer.cpp index 6d692d4fbf6..12b8f89fce3 100644 --- a/src/timer.cpp +++ b/src/timer.cpp @@ -25,6 +25,7 @@ Timer time_event_calculate_xs; Timer time_event_advance_particle; Timer time_event_surface_crossing; Timer time_event_collision; +Timer time_event_temperature_mesh_crossing; Timer time_event_death; Timer time_update_src; @@ -85,6 +86,7 @@ void reset_timers() simulation::time_event_advance_particle.reset(); simulation::time_event_surface_crossing.reset(); simulation::time_event_collision.reset(); + simulation::time_event_temperature_mesh_crossing.reset(); simulation::time_event_death.reset(); simulation::time_update_src.reset(); } diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index ce7e539ea57..5cba118e2ca 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -1,5 +1,6 @@ set(TEST_NAMES test_distribution + test_field test_file_utils test_tally test_interpolate @@ -14,5 +15,8 @@ set(TEST_NAMES foreach(test ${TEST_NAMES}) add_executable(${test} ${test}.cpp) target_link_libraries(${test} Catch2::Catch2WithMain libopenmc) + if (OPENMC_USE_DAGMC) + target_compile_definitions(${test} PRIVATE OPENMC_DAGMC_ENABLED) + endif() add_test(NAME ${test} COMMAND ${test} WORKING_DIRECTORY ${UNIT_TEST_BIN_OUTPUT_DIR}) endforeach() diff --git a/tests/cpp_unit_tests/test_field.cpp b/tests/cpp_unit_tests/test_field.cpp new file mode 100644 index 00000000000..0e0acbfbbc1 --- /dev/null +++ b/tests/cpp_unit_tests/test_field.cpp @@ -0,0 +1,146 @@ +#include + +#include +#include +#include +#include +#include + +#include "openmc/field.h" +#include "openmc/geometry.h" +#include "openmc/mesh.h" + +using namespace openmc; + +TEST_CASE("Test TemperatureField functions with a regular mesh") +{ + // The XML data as a string + std::string xml_string = R"( + + 2 2 2 + -1 -1 -1 + 1 1 1 + + )"; + + // Create the mesh from a file + pugi::xml_document doc; + pugi::xml_parse_result result = doc.load_string(xml_string.c_str()); + pugi::xml_node root = doc.child("mesh"); + auto mesh = RegularMesh(root); + + // Define some temperature values + vector values = {10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0}; + + // Create a temperature field + TemperatureField temp_field = TemperatureField(&mesh, values); + + // Get temperature + REQUIRE(temp_field.get_temperature(0) == 10.0); + REQUIRE(temp_field.get_temperature(1) == 20.0); + REQUIRE(temp_field.get_temperature(2) == 30.0); + REQUIRE(temp_field.get_temperature(6) == 70.0); + REQUIRE(temp_field.get_temperature(7) == 80.0); + REQUIRE(temp_field.get_temperature(-1) == -1.0); + + // Get sqrtkT + REQUIRE(temp_field.get_sqrtkT(7) == Catch::Approx(0.083029).margin(1.0E-6)); + + // Get mesh bin + REQUIRE(temp_field.get_mesh_bin(Position(0.5, 0.5, 0.5)) == 7); + REQUIRE(temp_field.get_mesh_bin(Position(-0.5, -0.5, -0.5)) == 0); + REQUIRE(temp_field.get_mesh_bin(Position(0.5, -0.5, -0.5)) == 1); + REQUIRE(temp_field.get_mesh_bin(Position(-0.5, -0.5, 0.5)) == 4); + REQUIRE(temp_field.get_mesh_bin(Position(0.0, 0.0, 0.0)) == 0); + REQUIRE(temp_field.get_mesh_bin(Position(2.0, 2.0, 2.0)) == -1); +} + +TEST_CASE("Test settings declaration exceptions for a temperature field", + "[generators]") +{ + auto [input, error] = GENERATE(table({ + {// If the number of values is not equal to the number of bins -> error + R"( + + eigenvalue + 200 + 20 + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + + )", + "Inconsistency in the temperature field: the number of values must be " + "equal to the number of bins in the mesh."}, + {// If the mesh declared is not defined -> error + R"( + + eigenvalue + 200 + 20 + + 2 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + + )", + "Mesh 2 specified for the temperature field does not exist."}, + {// No mesh declared -> error + R"( + + eigenvalue + 200 + 20 + + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + + )", + "A mesh must be given for the temperature field."}, + {// No values declared -> error + R"( + + eigenvalue + 200 + 20 + + 1 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + + )", + "Temperature values must be given for the temperature field."}, + })); + + free_memory_mesh(); + free_memory_settings(); + settings::run_mode = RunMode::UNSET; + + pugi::xml_document doc; + pugi::xml_parse_result result = doc.load_string(input.c_str()); + pugi::xml_node root = doc.child("settings"); + + CAPTURE(input); + REQUIRE_THROWS_WITH(read_settings_xml(root), error); + doc.reset(); +} diff --git a/tests/cpp_unit_tests/test_mesh.cpp b/tests/cpp_unit_tests/test_mesh.cpp index 24c4f77373f..0dc993bfbf3 100644 --- a/tests/cpp_unit_tests/test_mesh.cpp +++ b/tests/cpp_unit_tests/test_mesh.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -255,3 +256,131 @@ TEST_CASE("Test multiple meshes HDF5 roundtrip - spherical") REQUIRE(regular_mesh_hdf5->lower_left() == regular_mesh_xml->lower_left()); REQUIRE(regular_mesh_hdf5->upper_right() == regular_mesh_xml->upper_right()); } + +TEST_CASE("Test distance_to_next_boundary() - regular") +{ + // The XML data as a string + std::string xml_string = R"( + + 2 2 2 + -1 -1 -1 + 1 1 1 + + )"; + + // Create the mesh from a file + pugi::xml_document doc; + pugi::xml_parse_result result = doc.load_string(xml_string.c_str()); + pugi::xml_node root = doc.child("mesh"); + auto mesh = RegularMesh(root); + + int current_bin; + Position r; + Position u; + int next_bin; + double distance; + + // Test inside the mesh + current_bin = 0; + r = Position(0.0, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == 1.0); + REQUIRE(next_bin == -1); + + // Test outside the mesh, going toward the mesh + current_bin = -1; + r = Position(-2.5, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == 1.5); + REQUIRE(next_bin == 0); + + // Test outside the mesh, not going toward the mesh + current_bin = -1; + r = Position(-2.0, 0.0, 0.0); + u = Position(-1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == INFTY); + REQUIRE(next_bin == -1); + + // Test on the mesh boundary, leaving the mesh + current_bin = 1; + r = Position(1.0, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == INFTY); + REQUIRE(next_bin == -1); + + // Test close to the mesh boundary + current_bin = 1; + r = Position(0.99999999999, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == Catch::Approx(0.00000000001).margin(1.0E-12)); + REQUIRE(next_bin == -1); +} + +TEST_CASE("Test distance_to_next_boundary() - rectilinear") +{ + // The XML data as a string + std::string xml_string = R"( + + -1.0 0.5 1.0 + -1.0 0.1 1.0 + -1.0 0.2 1.0 + + )"; + + // Create the mesh from a file + pugi::xml_document doc; + pugi::xml_parse_result result = doc.load_string(xml_string.c_str()); + pugi::xml_node root = doc.child("mesh"); + auto mesh = RectilinearMesh(root); + + int current_bin; + Position r; + Position u; + int next_bin; + double distance; + + // Test inside the mesh + current_bin = 0; + r = Position(0.5, 0.1, 0.2); + u = Position(1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == 0.5); + REQUIRE(next_bin == -1); + + // Test outside the mesh, going toward the mesh + current_bin = -1; + r = Position(-2.5, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == 1.5); + REQUIRE(next_bin == 0); + + // Test outside the mesh, not going toward the mesh + current_bin = -1; + r = Position(-2.0, 0.0, 0.0); + u = Position(-1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == INFTY); + REQUIRE(next_bin == -1); + + // Test on the mesh boundary, leaving the mesh + current_bin = 1; + r = Position(1.0, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == INFTY); + REQUIRE(next_bin == -1); + + // Test close to the mesh boundary + current_bin = 1; + r = Position(0.99999999999, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + distance = mesh.distance_to_next_boundary(current_bin, r, u, next_bin); + REQUIRE(distance == Catch::Approx(0.00000000001).margin(1.0E-12)); + REQUIRE(next_bin == -1); +} diff --git a/tests/regression_tests/dnp_drift/__init__.py b/tests/regression_tests/dnp_drift/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/dnp_drift/regular_mesh/__init__.py b/tests/regression_tests/dnp_drift/regular_mesh/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/dnp_drift/regular_mesh/cell-based/__init__.py b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/dnp_drift/regular_mesh/cell-based/model.xml b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/model.xml new file mode 100644 index 00000000000..6ced13b8fce --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/model.xml @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 4000 + 20 + 2 + + + 0.0 0.0 0.0 10.0 10.0 10.0 + + + true + + + + cell + 1 + 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 + + 1 + 2 + 3 4 + + RK4 + 0.1 + true + 2.0 + + + 1 + 0 24 12 36 7 31 19 43 15 39 21 45 2 26 8 32 4 16 10 22 29 41 35 47 + 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + + + 2 2 2 + 0.0 0.0 0.0 + 10.0 10.0 10.0 + + interpolation + + diff --git a/tests/regression_tests/dnp_drift/regular_mesh/cell-based/results_true.dat b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/results_true.dat new file mode 100644 index 00000000000..16fba2f6b9e --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.459382E+00 1.712447E-03 diff --git a/tests/regression_tests/dnp_drift/regular_mesh/cell-based/test.py b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/test.py new file mode 100755 index 00000000000..5087d3ab711 --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/test.py @@ -0,0 +1,6 @@ +from tests.testing_harness import TestHarness + + +def test_dnp_drift_regular_mesh_cell_based(): + harness = TestHarness('statepoint.20.h5') + harness.main() diff --git a/tests/regression_tests/dnp_drift/regular_mesh/nodal/__init__.py b/tests/regression_tests/dnp_drift/regular_mesh/nodal/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/dnp_drift/regular_mesh/nodal/model.xml b/tests/regression_tests/dnp_drift/regular_mesh/nodal/model.xml new file mode 100644 index 00000000000..1f0a0d02c27 --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/nodal/model.xml @@ -0,0 +1,64 @@ + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 4000 + 20 + 2 + + + 0.0 0.0 0.0 10.0 10.0 10.0 + + + true + + + + nodal + 1 + + 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 + 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 + 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 + + + 1 + 2 + 3 4 + + RK4 + 0.1 + true + 2.0 + + + 1 + 0 24 12 36 7 31 19 43 15 39 21 45 2 26 8 32 4 16 10 22 29 41 35 47 + 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + + + 2 2 2 + 0.0 0.0 0.0 + 10.0 10.0 10.0 + + interpolation + + diff --git a/tests/regression_tests/dnp_drift/regular_mesh/nodal/results_true.dat b/tests/regression_tests/dnp_drift/regular_mesh/nodal/results_true.dat new file mode 100644 index 00000000000..16fba2f6b9e --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/nodal/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.459382E+00 1.712447E-03 diff --git a/tests/regression_tests/dnp_drift/regular_mesh/nodal/test.py b/tests/regression_tests/dnp_drift/regular_mesh/nodal/test.py new file mode 100755 index 00000000000..0aee79dfbfc --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/nodal/test.py @@ -0,0 +1,6 @@ +from tests.testing_harness import TestHarness + + +def test_dnp_drift_regular_mesh_nodal(): + harness = TestHarness('statepoint.20.h5') + harness.main() diff --git a/tests/regression_tests/temperature_field/README.rst b/tests/regression_tests/temperature_field/README.rst new file mode 100644 index 00000000000..b035c9aa601 --- /dev/null +++ b/tests/regression_tests/temperature_field/README.rst @@ -0,0 +1,29 @@ +Temperature field regression tests +================================== + +These regressions tests use two models: +- single_cube: a 5x5x5cm cube, +- nested_cubes: a 5x5x5cm cube surrounded by a 15x15x15cm cube. + +Three types of boundary conditions are tested: +- vacuum, +- reflective, +- periodic (except with DAGMC geometries). + +The temperature field uses a 2x2x2 regular mesh where each cell has its own +temperature value (built by starting with 294K and adding increments of 100K +to assign the next cell's temperature). When modeled, the temperature of the +surrounding cube is set to 250K. + +For DAGMC geometries, the boundary conditions and temperature values are +directly stored in the DAGMC files created with a modified version of +stellarmesh. + +Verifications scripts are made available in **_verification**. They are designed +to compare the results of the regression tests to a strictly equivalent model +where the tracks of particles are expected to be identical. These scripts +should be run using similar settings than what is used with the regression tests. + +.. warning:: The nested cubes verification scripts for dagmc are currently missing + because there are errors with the geometry. These errors are currently under + investigation and the scripts will be made available once the problem is fixed. diff --git a/tests/regression_tests/temperature_field/__init__.py b/tests/regression_tests/temperature_field/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/reflective_bc/single_cube_multicell_reflecting.h5m b/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/reflective_bc/single_cube_multicell_reflecting.h5m new file mode 100644 index 00000000000..1ea62214f59 Binary files /dev/null and b/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/reflective_bc/single_cube_multicell_reflecting.h5m differ diff --git a/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/reflective_bc/verification.py b/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/reflective_bc/verification.py new file mode 100644 index 00000000000..1661dcdf131 --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/reflective_bc/verification.py @@ -0,0 +1,77 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell. +The cube has reflective boundary conditions. + +The entire geometry is defined as a DAGMC geometry which includes the +boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material(name="fuel") +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Read DAGMC universe +dag_univ = openmc.DAGMCUniverse(filename="single_cube_multicell_reflecting.h5m", auto_geom_ids=True) +model.geometry = openmc.Geometry(dag_univ) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box(lower_left, upper_right) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/vacuum_bc/single_cube_multicell_vacuum.h5m b/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/vacuum_bc/single_cube_multicell_vacuum.h5m new file mode 100644 index 00000000000..5abac1fd09f Binary files /dev/null and b/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/vacuum_bc/single_cube_multicell_vacuum.h5m differ diff --git a/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/vacuum_bc/verification.py b/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/vacuum_bc/verification.py new file mode 100644 index 00000000000..cb2f9865d1b --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/dagmc/single_cube/vacuum_bc/verification.py @@ -0,0 +1,77 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell. +The cube has vacuum boundary conditions. + +The entire geometry is defined as a DAGMC geometry which includes the +boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material(name="fuel") +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Read DAGMC universe +dag_univ = openmc.DAGMCUniverse(filename="single_cube_multicell_vacuum.h5m", auto_geom_ids=True) +model.geometry = openmc.Geometry(dag_univ) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box(lower_left, upper_right) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/periodic_bc/verification.py b/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/periodic_bc/verification.py new file mode 100644 index 00000000000..ed1127350de --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/periodic_bc/verification.py @@ -0,0 +1,115 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +which is surrounded by another cube using a lattice. +The surrounding cube has periodic boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Create geometry +box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 +) + +universes = [] +for i, t in enumerate(temperature_values): + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + c.temperature = temperature_values[i] + u.add_cell(c) + universes.append(u) + +lattice = openmc.RectLattice() +lattice.lower_left = (0., 0., 0.) +lattice.pitch = (2.5, 2.5, 2.5) +lattice.universes = [ + [[universes[2], universes[3]], [universes[0], universes[1]]], + [[universes[6], universes[7]], [universes[4], universes[5]]] +] +outer_cell = openmc.Cell(fill=mat) +outer_cell.temperature = 250.0 +lattice.outer = openmc.Universe(cells=[outer_cell]) + +xmin = openmc.XPlane(x0=-5.0, boundary_type="periodic") +xmax = openmc.XPlane(x0=10.0, boundary_type="periodic") +ymin = openmc.YPlane(y0=-5.0, boundary_type="periodic") +ymax = openmc.YPlane(y0=10.0, boundary_type="periodic") +zmin = openmc.ZPlane(z0=-5.0, boundary_type="periodic") +zmax = openmc.ZPlane(z0=10.0, boundary_type="periodic") + +xmin.periodic_surface = xmax +ymin.periodic_surface = ymax +zmin.periodic_surface = zmax + +cell_2 = openmc.Cell(fill=lattice, region=(+xmin & -xmax & +ymin & -ymax & +zmin & -zmax)) + +model.geometry = openmc.Geometry([cell_2]) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/reflective_bc/verification.py b/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/reflective_bc/verification.py new file mode 100644 index 00000000000..2eded430c14 --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/reflective_bc/verification.py @@ -0,0 +1,110 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +which is surrounded by another cube using a lattice. +The surrounding cube has reflective boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Create geometry +box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 +) + +universes = [] +for i, t in enumerate(temperature_values): + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + c.temperature = temperature_values[i] + u.add_cell(c) + universes.append(u) + +lattice = openmc.RectLattice() +lattice.lower_left = (0., 0., 0.) +lattice.pitch = (2.5, 2.5, 2.5) +lattice.universes = [ + [[universes[2], universes[3]], [universes[0], universes[1]]], + [[universes[6], universes[7]], [universes[4], universes[5]]] +] +outer_cell = openmc.Cell(fill=mat) +outer_cell.temperature = 250.0 +lattice.outer = openmc.Universe(cells=[outer_cell]) + +surrounding_box = openmc.model.RectangularParallelepiped( + xmin=-5.0, xmax=10.0, + ymin=-5.0, ymax=10.0, + zmin=-5.0, zmax=10.0, + boundary_type="reflective" +) +cell_2 = openmc.Cell(fill=lattice, region=(-surrounding_box)) + +model.geometry = openmc.Geometry([cell_2]) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/vacuum_bc/verification.py b/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/vacuum_bc/verification.py new file mode 100644 index 00000000000..f441ea6572e --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/lattices/nested_cubes/vacuum_bc/verification.py @@ -0,0 +1,110 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +which is surrounded by another cube using a lattice. +The surrounding cube has vacuum boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Create geometry +box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 +) + +universes = [] +for i, t in enumerate(temperature_values): + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + c.temperature = temperature_values[i] + u.add_cell(c) + universes.append(u) + +lattice = openmc.RectLattice() +lattice.lower_left = (0., 0., 0.) +lattice.pitch = (2.5, 2.5, 2.5) +lattice.universes = [ + [[universes[2], universes[3]], [universes[0], universes[1]]], + [[universes[6], universes[7]], [universes[4], universes[5]]] +] +outer_cell = openmc.Cell(fill=mat) +outer_cell.temperature = 250.0 +lattice.outer = openmc.Universe(cells=[outer_cell]) + +surrounding_box = openmc.model.RectangularParallelepiped( + xmin=-5.0, xmax=10.0, + ymin=-5.0, ymax=10.0, + zmin=-5.0, zmax=10.0, + boundary_type="vacuum" +) +cell_2 = openmc.Cell(fill=lattice, region=(-surrounding_box)) + +model.geometry = openmc.Geometry([cell_2]) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/lattices/single_cube/periodic_bc/verification.py b/tests/regression_tests/temperature_field/_verification/lattices/single_cube/periodic_bc/verification.py new file mode 100644 index 00000000000..abbb3086fdc --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/lattices/single_cube/periodic_bc/verification.py @@ -0,0 +1,114 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +using a lattice with periodic boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Create geometry +box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 +) + +universes = [] +for i, t in enumerate(temperature_values): + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + c.temperature = temperature_values[i] + u.add_cell(c) + universes.append(u) + +lattice = openmc.RectLattice() +lattice.lower_left = (0., 0., 0.) +lattice.pitch = (2.5, 2.5, 2.5) +lattice.universes = [ + [[universes[2], universes[3]], [universes[0], universes[1]]], + [[universes[6], universes[7]], [universes[4], universes[5]]] +] +outer_cell = openmc.Cell(fill=mat) +outer_cell.temperature = 250.0 +lattice.outer = openmc.Universe(cells=[outer_cell]) + +xmin = openmc.XPlane(x0=0.0, boundary_type="periodic") +xmax = openmc.XPlane(x0=5.0, boundary_type="periodic") +ymin = openmc.YPlane(y0=0.0, boundary_type="periodic") +ymax = openmc.YPlane(y0=5.0, boundary_type="periodic") +zmin = openmc.ZPlane(z0=0.0, boundary_type="periodic") +zmax = openmc.ZPlane(z0=5.0, boundary_type="periodic") + +xmin.periodic_surface = xmax +ymin.periodic_surface = ymax +zmin.periodic_surface = zmax + +cell_2 = openmc.Cell(fill=lattice, region=(+xmin & -xmax & +ymin & -ymax & +zmin & -zmax)) + +model.geometry = openmc.Geometry([cell_2]) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box(lower_left, upper_right) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/lattices/single_cube/reflective_bc/verification.py b/tests/regression_tests/temperature_field/_verification/lattices/single_cube/reflective_bc/verification.py new file mode 100644 index 00000000000..1bec0b16fa3 --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/lattices/single_cube/reflective_bc/verification.py @@ -0,0 +1,109 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +using a lattice with reflective boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Create geometry +box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 +) + +universes = [] +for i, t in enumerate(temperature_values): + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + c.temperature = temperature_values[i] + u.add_cell(c) + universes.append(u) + +lattice = openmc.RectLattice() +lattice.lower_left = (0., 0., 0.) +lattice.pitch = (2.5, 2.5, 2.5) +lattice.universes = [ + [[universes[2], universes[3]], [universes[0], universes[1]]], + [[universes[6], universes[7]], [universes[4], universes[5]]] +] +outer_cell = openmc.Cell(fill=mat) +outer_cell.temperature = 250.0 +lattice.outer = openmc.Universe(cells=[outer_cell]) + +surrounding_box = openmc.model.RectangularParallelepiped( + xmin=0.0, xmax=5.0, + ymin=0.0, ymax=5.0, + zmin=0.0, zmax=5.0, + boundary_type="reflective" +) +cell_2 = openmc.Cell(fill=lattice, region=(-surrounding_box)) + +model.geometry = openmc.Geometry([cell_2]) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box(lower_left, upper_right) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/lattices/single_cube/vacuum_bc/verification.py b/tests/regression_tests/temperature_field/_verification/lattices/single_cube/vacuum_bc/verification.py new file mode 100644 index 00000000000..9f23071be86 --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/lattices/single_cube/vacuum_bc/verification.py @@ -0,0 +1,109 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +using a lattice with vacuum boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Create geometry +box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 +) + +universes = [] +for i, t in enumerate(temperature_values): + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + c.temperature = temperature_values[i] + u.add_cell(c) + universes.append(u) + +lattice = openmc.RectLattice() +lattice.lower_left = (0., 0., 0.) +lattice.pitch = (2.5, 2.5, 2.5) +lattice.universes = [ + [[universes[2], universes[3]], [universes[0], universes[1]]], + [[universes[6], universes[7]], [universes[4], universes[5]]] +] +outer_cell = openmc.Cell(fill=mat) +outer_cell.temperature = 250.0 +lattice.outer = openmc.Universe(cells=[outer_cell]) + +surrounding_box = openmc.model.RectangularParallelepiped( + xmin=0.0, xmax=5.0, + ymin=0.0, ymax=5.0, + zmin=0.0, zmax=5.0, + boundary_type="vacuum" +) +cell_2 = openmc.Cell(fill=lattice, region=(-surrounding_box)) + +model.geometry = openmc.Geometry([cell_2]) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box(lower_left, upper_right) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/nested_cubes/periodic_bc/verification.py b/tests/regression_tests/temperature_field/_verification/nested_cubes/periodic_bc/verification.py new file mode 100644 index 00000000000..44c4d1cb120 --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/nested_cubes/periodic_bc/verification.py @@ -0,0 +1,146 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +which is surrounded by another cube. +The surrounding cube has periodic boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Create cells +x_planes = [ + openmc.XPlane(x0=-5., boundary_type="periodic"), + openmc.XPlane(x0=0.), + openmc.XPlane(x0=2.5), + openmc.XPlane(x0=5.), + openmc.XPlane(x0=10., boundary_type="periodic") +] + +y_planes = [ + openmc.YPlane(y0=-5., boundary_type="periodic"), + openmc.YPlane(y0=0.), + openmc.YPlane(y0=2.5), + openmc.YPlane(y0=5.), + openmc.YPlane(y0=10., boundary_type="periodic") +] + +z_planes = [ + openmc.ZPlane(z0=-5., boundary_type="periodic"), + openmc.ZPlane(z0=0.), + openmc.ZPlane(z0=2.5), + openmc.ZPlane(z0=5.), + openmc.ZPlane(z0=10., boundary_type="periodic") +] + +x_planes[0].periodic_surface = x_planes[-1] +y_planes[0].periodic_surface = y_planes[-1] +z_planes[0].periodic_surface = z_planes[-1] + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +cells = [] +i = 0 +for iz in range(dim): + for iy in range(dim): + for ix in range(dim): + print(ix, iy, iz) + region = ( + +x_planes[ix+1] & -x_planes[ix+2] & + +y_planes[iy+1] & -y_planes[iy+2] & + +z_planes[iz+1] & -z_planes[iz+2] + ) + cell = openmc.Cell( + fill=mat, + region=region + ) + + cell.temperature = temperature_values[i] + cells.append(cell) + i += 1 + +# Surrounding cell +box1_region = (+x_planes[0] & -x_planes[1] & +y_planes[0] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box2_region = (+x_planes[-2] & -x_planes[-1] & +y_planes[0] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box3_region = (+x_planes[1] & -x_planes[-2] & +y_planes[0] & -y_planes[1] & +z_planes[0] & -z_planes[-1]) +box4_region = (+x_planes[1] & -x_planes[-2] & +y_planes[-2] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box5_region = (+x_planes[1] & -x_planes[-2] & +y_planes[1] & -y_planes[-2] & +z_planes[0] & -z_planes[1]) +box6_region = (+x_planes[1] & -x_planes[-2] & +y_planes[1] & -y_planes[-2] & +z_planes[-2] & -z_planes[-1]) + +surrounding_region = box1_region | box2_region | box3_region | box4_region | box5_region | box6_region + +cell = openmc.Cell( + fill=mat, + region=surrounding_region +) +cell.temperature = 250.0 + +cells.append(cell) + +root_universe = openmc.Universe(cells=cells) + +# Register geometry +model.geometry = openmc.Geometry(root_universe) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/nested_cubes/reflective_bc/verification.py b/tests/regression_tests/temperature_field/_verification/nested_cubes/reflective_bc/verification.py new file mode 100644 index 00000000000..1c997a122fd --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/nested_cubes/reflective_bc/verification.py @@ -0,0 +1,142 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +which is surrounded by another cube. +The surrounding cube has reflective boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Create cells +x_planes = [ + openmc.XPlane(x0=-5., boundary_type="reflective"), + openmc.XPlane(x0=0.), + openmc.XPlane(x0=2.5), + openmc.XPlane(x0=5.), + openmc.XPlane(x0=10., boundary_type="reflective") +] + +y_planes = [ + openmc.YPlane(y0=-5., boundary_type="reflective"), + openmc.YPlane(y0=0.), + openmc.YPlane(y0=2.5), + openmc.YPlane(y0=5.), + openmc.YPlane(y0=10., boundary_type="reflective") +] + +z_planes = [ + openmc.ZPlane(z0=-5., boundary_type="reflective"), + openmc.ZPlane(z0=0.), + openmc.ZPlane(z0=2.5), + openmc.ZPlane(z0=5.), + openmc.ZPlane(z0=10., boundary_type="reflective") +] + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +cells = [] +i = 0 +for iz in range(dim): + for iy in range(dim): + for ix in range(dim): + print(ix, iy, iz) + region = ( + +x_planes[ix+1] & -x_planes[ix+2] & + +y_planes[iy+1] & -y_planes[iy+2] & + +z_planes[iz+1] & -z_planes[iz+2] + ) + cell = openmc.Cell( + fill=mat, + region=region + ) + + cell.temperature = temperature_values[i] + cells.append(cell) + i += 1 + +# Surrounding cell +box1_region = (+x_planes[0] & -x_planes[1] & +y_planes[0] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box2_region = (+x_planes[-2] & -x_planes[-1] & +y_planes[0] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box3_region = (+x_planes[1] & -x_planes[-2] & +y_planes[0] & -y_planes[1] & +z_planes[0] & -z_planes[-1]) +box4_region = (+x_planes[1] & -x_planes[-2] & +y_planes[-2] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box5_region = (+x_planes[1] & -x_planes[-2] & +y_planes[1] & -y_planes[-2] & +z_planes[0] & -z_planes[1]) +box6_region = (+x_planes[1] & -x_planes[-2] & +y_planes[1] & -y_planes[-2] & +z_planes[-2] & -z_planes[-1]) + +surrounding_region = box1_region | box2_region | box3_region | box4_region | box5_region | box6_region + +cell = openmc.Cell( + fill=mat, + region=surrounding_region +) +cell.temperature = 250.0 + +cells.append(cell) + +root_universe = openmc.Universe(cells=cells) + +# Register geometry +model.geometry = openmc.Geometry(root_universe) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/nested_cubes/vacuum_bc/verification.py b/tests/regression_tests/temperature_field/_verification/nested_cubes/vacuum_bc/verification.py new file mode 100644 index 00000000000..3c6ebcc7a96 --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/nested_cubes/vacuum_bc/verification.py @@ -0,0 +1,142 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell +which is surrounded by another cube. +The surrounding cube has vacuum boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Create cells +x_planes = [ + openmc.XPlane(x0=-5., boundary_type="vacuum"), + openmc.XPlane(x0=0.), + openmc.XPlane(x0=2.5), + openmc.XPlane(x0=5.), + openmc.XPlane(x0=10., boundary_type="vacuum") +] + +y_planes = [ + openmc.YPlane(y0=-5., boundary_type="vacuum"), + openmc.YPlane(y0=0.), + openmc.YPlane(y0=2.5), + openmc.YPlane(y0=5.), + openmc.YPlane(y0=10., boundary_type="vacuum") +] + +z_planes = [ + openmc.ZPlane(z0=-5., boundary_type="vacuum"), + openmc.ZPlane(z0=0.), + openmc.ZPlane(z0=2.5), + openmc.ZPlane(z0=5.), + openmc.ZPlane(z0=10., boundary_type="vacuum") +] + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +cells = [] +i = 0 +for iz in range(dim): + for iy in range(dim): + for ix in range(dim): + print(ix, iy, iz) + region = ( + +x_planes[ix+1] & -x_planes[ix+2] & + +y_planes[iy+1] & -y_planes[iy+2] & + +z_planes[iz+1] & -z_planes[iz+2] + ) + cell = openmc.Cell( + fill=mat, + region=region + ) + + cell.temperature = temperature_values[i] + cells.append(cell) + i += 1 + +# Surrounding cell +box1_region = (+x_planes[0] & -x_planes[1] & +y_planes[0] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box2_region = (+x_planes[-2] & -x_planes[-1] & +y_planes[0] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box3_region = (+x_planes[1] & -x_planes[-2] & +y_planes[0] & -y_planes[1] & +z_planes[0] & -z_planes[-1]) +box4_region = (+x_planes[1] & -x_planes[-2] & +y_planes[-2] & -y_planes[-1] & +z_planes[0] & -z_planes[-1]) +box5_region = (+x_planes[1] & -x_planes[-2] & +y_planes[1] & -y_planes[-2] & +z_planes[0] & -z_planes[1]) +box6_region = (+x_planes[1] & -x_planes[-2] & +y_planes[1] & -y_planes[-2] & +z_planes[-2] & -z_planes[-1]) + +surrounding_region = box1_region | box2_region | box3_region | box4_region | box5_region | box6_region + +cell = openmc.Cell( + fill=mat, + region=surrounding_region +) +cell.temperature = 250.0 + +cells.append(cell) + +root_universe = openmc.Universe(cells=cells) + +# Register geometry +model.geometry = openmc.Geometry(root_universe) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/single_cube/periodic_bc/verification.py b/tests/regression_tests/temperature_field/_verification/single_cube/periodic_bc/verification.py new file mode 100644 index 00000000000..b76fbc36cfa --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/single_cube/periodic_bc/verification.py @@ -0,0 +1,120 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell. +The cube has periodic boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Create cells +x_planes = [ + openmc.XPlane(x0=0., boundary_type="periodic"), + openmc.XPlane(x0=2.5), + openmc.XPlane(x0=5., boundary_type="periodic") +] + +y_planes = [ + openmc.YPlane(y0=0., boundary_type="periodic"), + openmc.YPlane(y0=2.5), + openmc.YPlane(y0=5., boundary_type="periodic") +] + +z_planes = [ + openmc.ZPlane(z0=0., boundary_type="periodic"), + openmc.ZPlane(z0=2.5), + openmc.ZPlane(z0=5., boundary_type="periodic") +] + +x_planes[0].periodic_surface = x_planes[2] +y_planes[0].periodic_surface = y_planes[2] +z_planes[0].periodic_surface = z_planes[2] + +cells = [] +for iz in range(dim): + for iy in range(dim): + for ix in range(dim): + region = ( + +x_planes[ix] & -x_planes[ix+1] & + +y_planes[iy] & -y_planes[iy+1] & + +z_planes[iz] & -z_planes[iz+1] + ) + cell = openmc.Cell( + fill=mat, + region=region + ) + cells.append(cell) + +root_universe = openmc.Universe(cells=cells) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Set the temperature to each CSG cell +for i, c in enumerate(cells): + c.temperature = temperature_values[i] + +# Register geometry +model.geometry = openmc.Geometry(root_universe) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box(lower_left, upper_right) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/single_cube/reflective_bc/verification.py b/tests/regression_tests/temperature_field/_verification/single_cube/reflective_bc/verification.py new file mode 100644 index 00000000000..4c729b13378 --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/single_cube/reflective_bc/verification.py @@ -0,0 +1,116 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell. +The cube has reflective boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Create cells +x_planes = [ + openmc.XPlane(x0=0., boundary_type="reflective"), + openmc.XPlane(x0=2.5), + openmc.XPlane(x0=5., boundary_type="reflective") +] + +y_planes = [ + openmc.YPlane(y0=0., boundary_type="reflective"), + openmc.YPlane(y0=2.5), + openmc.YPlane(y0=5., boundary_type="reflective") +] + +z_planes = [ + openmc.ZPlane(z0=0., boundary_type="reflective"), + openmc.ZPlane(z0=2.5), + openmc.ZPlane(z0=5., boundary_type="reflective") +] + +cells = [] +for iz in range(dim): + for iy in range(dim): + for ix in range(dim): + region = ( + +x_planes[ix] & -x_planes[ix+1] & + +y_planes[iy] & -y_planes[iy+1] & + +z_planes[iz] & -z_planes[iz+1] + ) + cell = openmc.Cell( + fill=mat, + region=region + ) + cells.append(cell) + +root_universe = openmc.Universe(cells=cells) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Set the temperature to each CSG cell +for i, c in enumerate(cells): + c.temperature = temperature_values[i] + +# Register geometry +model.geometry = openmc.Geometry(root_universe) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box(lower_left, upper_right) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/_verification/single_cube/vacuum_bc/verification.py b/tests/regression_tests/temperature_field/_verification/single_cube/vacuum_bc/verification.py new file mode 100644 index 00000000000..7aad432f607 --- /dev/null +++ b/tests/regression_tests/temperature_field/_verification/single_cube/vacuum_bc/verification.py @@ -0,0 +1,116 @@ +"""Runs OpenMC on a 2x2x2 cube with different temperature values for each cell. +The cube has vacuum boundary conditions. + +A "results.dat" file is created and can be compared to the "results_true.dat" +file from the corresponding case. + +""" + +import openmc +import numpy as np + +model = openmc.Model() + +# Material +mat = openmc.Material() +mat.add_nuclide("U235", 0.2) +mat.add_nuclide("U238", 0.8) +mat.add_element("O", 2.0) +mat.add_element("H", 4.0) +mat.set_density("g/cm3", 5.0) +mat.add_s_alpha_beta('c_H_in_H2O') +model.materials = openmc.Materials([mat]) + +# Create mesh for tallying +dim = 2 +lower_left = (0.0, 0.0, 0.0) +upper_right = (5.0, 5.0, 5.0) +mesh = openmc.RegularMesh() +mesh.lower_left = lower_left +mesh.upper_right = upper_right +mesh.dimension = (dim, dim, dim) + +# Create cells +x_planes = [ + openmc.XPlane(x0=0., boundary_type="vacuum"), + openmc.XPlane(x0=2.5), + openmc.XPlane(x0=5., boundary_type="vacuum") +] + +y_planes = [ + openmc.YPlane(y0=0., boundary_type="vacuum"), + openmc.YPlane(y0=2.5), + openmc.YPlane(y0=5., boundary_type="vacuum") +] + +z_planes = [ + openmc.ZPlane(z0=0., boundary_type="vacuum"), + openmc.ZPlane(z0=2.5), + openmc.ZPlane(z0=5., boundary_type="vacuum") +] + +cells = [] +for iz in range(dim): + for iy in range(dim): + for ix in range(dim): + region = ( + +x_planes[ix] & -x_planes[ix+1] & + +y_planes[iy] & -y_planes[iy+1] & + +z_planes[iz] & -z_planes[iz+1] + ) + cell = openmc.Cell( + fill=mat, + region=region + ) + cells.append(cell) + +root_universe = openmc.Universe(cells=cells) + +# Temperature values +temperature_values = [294.0 + i * 100 for i in range(dim**3)] + +# Set the temperature to each CSG cell +for i, c in enumerate(cells): + c.temperature = temperature_values[i] + +# Register geometry +model.geometry = openmc.Geometry(root_universe) + +# Settings +settings = openmc.Settings() +settings.batches = 20 +settings.particles = 200 +spatial_dist = openmc.stats.Box(lower_left, upper_right) +settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) +settings.temperature = {'tolerance': 1000, 'multipole': True} +model.settings = settings + +# Add tallies +mesh_filter = openmc.MeshFilter(mesh) +mesh_tally = openmc.Tally(name="total reaction rate") +mesh_tally.filters = [mesh_filter] +mesh_tally.scores = ["total"] +tallies = openmc.Tallies([mesh_tally]) +model.tallies = tallies + +model.run() + +with openmc.StatePoint(f"statepoint.20.h5") as sp: + + outstr = 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + + for i, tally_ind in enumerate(sp.tallies): + tally = sp.tallies[tally_ind] + results = np.zeros((tally.sum.size * 2, )) + results[0::2] = tally.sum.ravel() + results[1::2] = tally.sum_sq.ravel() + results = ['{0:12.6E}'.format(x) for x in results] + + outstr += 'tally {}:\n'.format(i + 1) + outstr += '\n'.join(results) + '\n' + + with open("results.dat", "w") as res: + res.write(outstr) diff --git a/tests/regression_tests/temperature_field/dagmc/__init__.py b/tests/regression_tests/temperature_field/dagmc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/__init__.py b/tests/regression_tests/temperature_field/dagmc/nested_cubes/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/__init__.py b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/inputs_true.dat b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/inputs_true.dat new file mode 100644 index 00000000000..25c143085d0 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/inputs_true.dat @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/nested_cubes_cell_reflecting.h5m b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/nested_cubes_cell_reflecting.h5m new file mode 100644 index 00000000000..1932d6fe4a2 Binary files /dev/null and b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/nested_cubes_cell_reflecting.h5m differ diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/results_true.dat b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/results_true.dat new file mode 100644 index 00000000000..3afa4fd8117 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.515234E+00 9.683233E-03 +tally 1: +1.480117E+00 +1.283100E-01 +1.694708E+00 +1.654352E-01 +1.463603E+00 +1.204571E-01 +1.548104E+00 +1.354957E-01 +1.746713E+00 +1.729295E-01 +1.397141E+00 +1.116072E-01 +1.522349E+00 +1.334629E-01 +1.661876E+00 +1.587775E-01 diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/test.py b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/test.py new file mode 100644 index 00000000000..b210f8833cb --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/test.py @@ -0,0 +1,78 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +pytestmark = pytest.mark.skipif( + not openmc.lib._dagmc_enabled(), + reason="DAGMC is not enabled.") + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell. This box is surrounded by a larger + box where the temperature field is not defined and has reflective + boundary conditions. + + The entire geometry is defined as a DAGMC geometry which includes the + boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material(name="fuel") + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Read DAGMC universe + dag_univ = openmc.DAGMCUniverse(filename='nested_cubes_cell_reflecting.h5m', auto_geom_ids=True) + model.geometry = openmc.Geometry(dag_univ) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/__init__.py b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/inputs_true.dat b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/inputs_true.dat new file mode 100644 index 00000000000..4844da91eea --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/inputs_true.dat @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/nested_cubes_cell_vacuum.h5m b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/nested_cubes_cell_vacuum.h5m new file mode 100644 index 00000000000..ffc0bc0e583 Binary files /dev/null and b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/nested_cubes_cell_vacuum.h5m differ diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/results_true.dat b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/results_true.dat new file mode 100644 index 00000000000..5bd690224e8 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +3.512624E-01 1.436222E-02 +tally 1: +9.836991E-01 +6.331373E-02 +9.860739E-01 +7.076155E-02 +1.051687E+00 +7.416489E-02 +8.660211E-01 +5.339067E-02 +9.091784E-01 +4.821934E-02 +1.146603E+00 +8.947049E-02 +8.048224E-01 +4.748861E-02 +9.541933E-01 +6.135353E-02 diff --git a/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/test.py b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/test.py new file mode 100644 index 00000000000..63923a935c6 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/test.py @@ -0,0 +1,78 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +pytestmark = pytest.mark.skipif( + not openmc.lib._dagmc_enabled(), + reason="DAGMC is not enabled.") + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell. This box is surrounded by a larger + box where the temperature field is not defined and has vacuum + boundary conditions. + + The entire geometry is defined as a DAGMC geometry which includes the + boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material(name="fuel") + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Read DAGMC universe + dag_univ = openmc.DAGMCUniverse(filename='nested_cubes_cell_vacuum.h5m', auto_geom_ids=True) + model.geometry = openmc.Geometry(dag_univ) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/__init__.py b/tests/regression_tests/temperature_field/dagmc/single_cube/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/__init__.py b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/inputs_true.dat b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/inputs_true.dat new file mode 100644 index 00000000000..90e56de9806 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/inputs_true.dat @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/results_true.dat b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/results_true.dat new file mode 100644 index 00000000000..7e16b60bd4b --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.507597E+00 1.811642E-02 +tally 1: +4.316452E+01 +9.417474E+01 +4.216471E+01 +8.960204E+01 +4.507663E+01 +1.029737E+02 +4.127791E+01 +8.578660E+01 +4.147002E+01 +8.693664E+01 +4.161294E+01 +8.721239E+01 +4.186858E+01 +8.850129E+01 +4.247133E+01 +9.097907E+01 diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/single_cube_cell_reflecting.h5m b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/single_cube_cell_reflecting.h5m new file mode 100644 index 00000000000..b75846f8313 Binary files /dev/null and b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/single_cube_cell_reflecting.h5m differ diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/test.py b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/test.py new file mode 100644 index 00000000000..d9cffcd3e17 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/test.py @@ -0,0 +1,76 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +pytestmark = pytest.mark.skipif( + not openmc.lib._dagmc_enabled(), + reason="DAGMC is not enabled.") + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell - reflective boundary conditions. + + The entire geometry is defined as a DAGMC geometry which includes the + boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material(name="fuel") + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Read DAGMC universe + dag_univ = openmc.DAGMCUniverse(filename='single_cube_cell_reflecting.h5m', auto_geom_ids=True) + model.geometry = openmc.Geometry(dag_univ) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/__init__.py b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/inputs_true.dat b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/inputs_true.dat new file mode 100644 index 00000000000..d94b2dcf5a9 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/inputs_true.dat @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/results_true.dat b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/results_true.dat new file mode 100644 index 00000000000..05ff69dd728 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +5.071104E-02 2.927355E-03 +tally 1: +2.731604E+00 +4.072231E-01 +3.388567E+00 +6.168992E-01 +2.435725E+00 +3.283434E-01 +2.724446E+00 +4.158287E-01 +3.426649E+00 +6.625678E-01 +2.929812E+00 +4.792204E-01 +2.484602E+00 +3.679255E-01 +2.506827E+00 +3.400903E-01 diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/single_cube_cell_vacuum.h5m b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/single_cube_cell_vacuum.h5m new file mode 100644 index 00000000000..2590f485ac9 Binary files /dev/null and b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/single_cube_cell_vacuum.h5m differ diff --git a/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/test.py b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/test.py new file mode 100644 index 00000000000..f7f72b76117 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/test.py @@ -0,0 +1,76 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +pytestmark = pytest.mark.skipif( + not openmc.lib._dagmc_enabled(), + reason="DAGMC is not enabled.") + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell - vacuum boundary conditions. + + The entire geometry is defined as a DAGMC geometry which includes the + boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material(name="fuel") + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Read DAGMC universe + dag_univ = openmc.DAGMCUniverse(filename='single_cube_cell_vacuum.h5m', auto_geom_ids=True) + model.geometry = openmc.Geometry(dag_univ) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/lattices/__init__.py b/tests/regression_tests/temperature_field/lattices/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/__init__.py b/tests/regression_tests/temperature_field/lattices/nested_cubes/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/__init__.py b/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/inputs_true.dat b/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/inputs_true.dat new file mode 100644 index 00000000000..6787b1debad --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/inputs_true.dat @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 3 + 2 2 2 + 0.0 0.0 0.0 + +1 1 +1 1 + +1 1 +1 1 + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/results_true.dat b/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/results_true.dat new file mode 100644 index 00000000000..22829229671 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.521509E+00 1.529817E-02 +tally 1: +1.687083E+00 +1.670239E-01 +1.615093E+00 +1.537921E-01 +1.655584E+00 +1.587550E-01 +1.606709E+00 +1.634004E-01 +1.465197E+00 +1.223359E-01 +1.525386E+00 +1.332914E-01 +1.538143E+00 +1.345132E-01 +1.620989E+00 +1.522083E-01 diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/test.py b/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/test.py new file mode 100644 index 00000000000..23fdb4eb740 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/test.py @@ -0,0 +1,97 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """ + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 + ) + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + u.add_cell(c) + + lattice = openmc.RectLattice() + lattice.lower_left = (0., 0., 0.) + lattice.pitch = (2.5, 2.5, 2.5) + lattice.universes = [ + [[u, u], [u, u]], + [[u, u], [u, u]] + ] + lattice.outer = openmc.Universe(cells=[openmc.Cell(fill=mat)]) + + xmin = openmc.XPlane(x0=-5.0, boundary_type="periodic") + xmax = openmc.XPlane(x0=10.0, boundary_type="periodic") + ymin = openmc.YPlane(y0=-5.0, boundary_type="periodic") + ymax = openmc.YPlane(y0=10.0, boundary_type="periodic") + zmin = openmc.ZPlane(z0=-5.0, boundary_type="periodic") + zmax = openmc.ZPlane(z0=10.0, boundary_type="periodic") + + xmin.periodic_surface = xmax + ymin.periodic_surface = ymax + zmin.periodic_surface = zmax + + cell_2 = openmc.Cell(fill=lattice, region=(+xmin & -xmax & +ymin & -ymax & +zmin & -zmax)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/__init__.py b/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/inputs_true.dat b/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/inputs_true.dat new file mode 100644 index 00000000000..a69be4ad8ca --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/inputs_true.dat @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 3 + 2 2 2 + 0.0 0.0 0.0 + +1 1 +1 1 + +1 1 +1 1 + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/results_true.dat b/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/results_true.dat new file mode 100644 index 00000000000..230771043b7 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.520805E+00 1.313130E-02 +tally 1: +1.620556E+00 +1.493663E-01 +1.646724E+00 +1.566329E-01 +1.504426E+00 +1.323039E-01 +1.519558E+00 +1.371785E-01 +1.319558E+00 +1.104879E-01 +1.454862E+00 +1.313438E-01 +1.380626E+00 +1.192778E-01 +1.409276E+00 +1.269625E-01 diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/test.py b/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/test.py new file mode 100644 index 00000000000..cf767fc4677 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/test.py @@ -0,0 +1,92 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """ + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 + ) + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + u.add_cell(c) + + lattice = openmc.RectLattice() + lattice.lower_left = (0., 0., 0.) + lattice.pitch = (2.5, 2.5, 2.5) + lattice.universes = [ + [[u, u], [u, u]], + [[u, u], [u, u]] + ] + lattice.outer = openmc.Universe(cells=[openmc.Cell(fill=mat)]) + + surrounding_box = openmc.model.RectangularParallelepiped( + xmin=-5.0, xmax=10.0, + ymin=-5.0, ymax=10.0, + zmin=-5.0, zmax=10.0, + boundary_type="reflective" + ) + cell_2 = openmc.Cell(fill=lattice, region=(-surrounding_box)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/__init__.py b/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/inputs_true.dat b/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/inputs_true.dat new file mode 100644 index 00000000000..9c6feb9d053 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/inputs_true.dat @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 3 + 2 2 2 + 0.0 0.0 0.0 + +1 1 +1 1 + +1 1 +1 1 + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/results_true.dat b/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/results_true.dat new file mode 100644 index 00000000000..b61c89f63d9 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +3.809948E-01 1.507099E-02 +tally 1: +1.047605E+00 +6.859683E-02 +9.076171E-01 +5.159802E-02 +9.946558E-01 +6.870688E-02 +1.088986E+00 +7.784666E-02 +8.504478E-01 +4.647715E-02 +1.056518E+00 +6.550817E-02 +1.037332E+00 +6.681476E-02 +8.884792E-01 +5.126623E-02 diff --git a/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/test.py b/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/test.py new file mode 100644 index 00000000000..41c3e5fdbb3 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/test.py @@ -0,0 +1,92 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """ + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 + ) + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + u.add_cell(c) + + lattice = openmc.RectLattice() + lattice.lower_left = (0., 0., 0.) + lattice.pitch = (2.5, 2.5, 2.5) + lattice.universes = [ + [[u, u], [u, u]], + [[u, u], [u, u]] + ] + lattice.outer = openmc.Universe(cells=[openmc.Cell(fill=mat)]) + + surrounding_box = openmc.model.RectangularParallelepiped( + xmin=-5.0, xmax=10.0, + ymin=-5.0, ymax=10.0, + zmin=-5.0, zmax=10.0, + boundary_type="vacuum" + ) + cell_2 = openmc.Cell(fill=lattice, region=(-surrounding_box)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/__init__.py b/tests/regression_tests/temperature_field/lattices/single_cube/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/__init__.py b/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/inputs_true.dat b/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/inputs_true.dat new file mode 100644 index 00000000000..f12851efd38 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/inputs_true.dat @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 3 + 2 2 2 + 0.0 0.0 0.0 + +1 1 +1 1 + +1 1 +1 1 + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/results_true.dat b/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/results_true.dat new file mode 100644 index 00000000000..bfed5a15158 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.511081E+00 1.098748E-02 +tally 1: +4.208340E+01 +8.914019E+01 +4.317567E+01 +9.354453E+01 +4.113208E+01 +8.499925E+01 +4.247578E+01 +9.078685E+01 +4.215567E+01 +8.953567E+01 +4.294495E+01 +9.271687E+01 +4.152245E+01 +8.657952E+01 +4.253211E+01 +9.093756E+01 diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/test.py b/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/test.py new file mode 100644 index 00000000000..de742e855bc --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/test.py @@ -0,0 +1,97 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """ + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 + ) + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + u.add_cell(c) + + lattice = openmc.RectLattice() + lattice.lower_left = (0., 0., 0.) + lattice.pitch = (2.5, 2.5, 2.5) + lattice.universes = [ + [[u, u], [u, u]], + [[u, u], [u, u]] + ] + lattice.outer = openmc.Universe(cells=[openmc.Cell(fill=mat)]) + + xmin = openmc.XPlane(x0=0.0, boundary_type="periodic") + xmax = openmc.XPlane(x0=5.0, boundary_type="periodic") + ymin = openmc.YPlane(y0=0.0, boundary_type="periodic") + ymax = openmc.YPlane(y0=5.0, boundary_type="periodic") + zmin = openmc.ZPlane(z0=0.0, boundary_type="periodic") + zmax = openmc.ZPlane(z0=5.0, boundary_type="periodic") + + xmin.periodic_surface = xmax + ymin.periodic_surface = ymax + zmin.periodic_surface = zmax + + cell_2 = openmc.Cell(fill=lattice, region=(+xmin & -xmax & +ymin & -ymax & +zmin & -zmax)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/__init__.py b/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/inputs_true.dat b/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/inputs_true.dat new file mode 100644 index 00000000000..c459536bf00 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/inputs_true.dat @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 3 + 2 2 2 + 0.0 0.0 0.0 + +1 1 +1 1 + +1 1 +1 1 + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/results_true.dat b/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/results_true.dat new file mode 100644 index 00000000000..7e16b60bd4b --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.507597E+00 1.811642E-02 +tally 1: +4.316452E+01 +9.417474E+01 +4.216471E+01 +8.960204E+01 +4.507663E+01 +1.029737E+02 +4.127791E+01 +8.578660E+01 +4.147002E+01 +8.693664E+01 +4.161294E+01 +8.721239E+01 +4.186858E+01 +8.850129E+01 +4.247133E+01 +9.097907E+01 diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/test.py b/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/test.py new file mode 100644 index 00000000000..aa70e9cc3d3 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/test.py @@ -0,0 +1,92 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """ + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 + ) + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + u.add_cell(c) + + lattice = openmc.RectLattice() + lattice.lower_left = (0., 0., 0.) + lattice.pitch = (2.5, 2.5, 2.5) + lattice.universes = [ + [[u, u], [u, u]], + [[u, u], [u, u]] + ] + lattice.outer = openmc.Universe(cells=[openmc.Cell(fill=mat)]) + + surrounding_box = openmc.model.RectangularParallelepiped( + xmin=0.0, xmax=5.0, + ymin=0.0, ymax=5.0, + zmin=0.0, zmax=5.0, + boundary_type="reflective" + ) + cell_2 = openmc.Cell(fill=lattice, region=(-surrounding_box)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/__init__.py b/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/inputs_true.dat b/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/inputs_true.dat new file mode 100644 index 00000000000..b4b9e40ac34 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/inputs_true.dat @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 3 + 2 2 2 + 0.0 0.0 0.0 + +1 1 +1 1 + +1 1 +1 1 + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/results_true.dat b/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/results_true.dat new file mode 100644 index 00000000000..05ff69dd728 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +5.071104E-02 2.927355E-03 +tally 1: +2.731604E+00 +4.072231E-01 +3.388567E+00 +6.168992E-01 +2.435725E+00 +3.283434E-01 +2.724446E+00 +4.158287E-01 +3.426649E+00 +6.625678E-01 +2.929812E+00 +4.792204E-01 +2.484602E+00 +3.679255E-01 +2.506827E+00 +3.400903E-01 diff --git a/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/test.py b/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/test.py new file mode 100644 index 00000000000..e4c446974de --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/test.py @@ -0,0 +1,92 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """ + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=-1.75, xmax=1.75, + ymin=-1.75, ymax=1.75, + zmin=-1.75, zmax=1.75 + ) + u = openmc.Universe() + c = openmc.Cell(fill=mat, region=-box) + u.add_cell(c) + + lattice = openmc.RectLattice() + lattice.lower_left = (0., 0., 0.) + lattice.pitch = (2.5, 2.5, 2.5) + lattice.universes = [ + [[u, u], [u, u]], + [[u, u], [u, u]] + ] + lattice.outer = openmc.Universe(cells=[openmc.Cell(fill=mat)]) + + surrounding_box = openmc.model.RectangularParallelepiped( + xmin=0.0, xmax=5.0, + ymin=0.0, ymax=5.0, + zmin=0.0, zmax=5.0, + boundary_type="vacuum" + ) + cell_2 = openmc.Cell(fill=lattice, region=(-surrounding_box)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/nested_cubes/__init__.py b/tests/regression_tests/temperature_field/nested_cubes/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/__init__.py b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/inputs_true.dat b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/inputs_true.dat new file mode 100644 index 00000000000..9c3995177be --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/inputs_true.dat @@ -0,0 +1,64 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/results_true.dat b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/results_true.dat new file mode 100644 index 00000000000..e29ba3bde07 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.543179E+00 1.653725E-02 +tally 1: +1.455073E+00 +1.171920E-01 +1.417464E+00 +1.144478E-01 +1.645363E+00 +1.548413E-01 +1.817800E+00 +1.966605E-01 +1.631346E+00 +1.603000E-01 +1.609351E+00 +1.517809E-01 +1.477630E+00 +1.303733E-01 +1.487139E+00 +1.239214E-01 diff --git a/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/test.py b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/test.py new file mode 100644 index 00000000000..4b1a78ef9ea --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/test.py @@ -0,0 +1,90 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell. This box is surrounded by a larger + box where the temperature field is not defined and has vacuum + boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=lower_left[0], xmax=upper_right[0], + ymin=lower_left[1], ymax=upper_right[1], + zmin=lower_left[2], zmax=upper_right[2] + ) + cell_1 = openmc.Cell(fill=mat, region=-box) + + xmin = openmc.XPlane(x0=-5.0, boundary_type="periodic") + xmax = openmc.XPlane(x0=10.0, boundary_type="periodic") + ymin = openmc.YPlane(y0=-5.0, boundary_type="periodic") + ymax = openmc.YPlane(y0=10.0, boundary_type="periodic") + zmin = openmc.ZPlane(z0=-5.0, boundary_type="periodic") + zmax = openmc.ZPlane(z0=10.0, boundary_type="periodic") + + xmin.periodic_surface = xmax + ymin.periodic_surface = ymax + zmin.periodic_surface = zmax + + cell_2 = openmc.Cell(fill=mat, region=((+xmin & -xmax & +ymin & -ymax & +zmin & -zmax) & +box)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_1, cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/__init__.py b/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/inputs_true.dat b/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/inputs_true.dat new file mode 100644 index 00000000000..e13abceaa10 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/inputs_true.dat @@ -0,0 +1,64 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/results_true.dat b/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/results_true.dat new file mode 100644 index 00000000000..22a999d7015 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.538904E+00 1.274036E-02 +tally 1: +1.394072E+00 +1.023965E-01 +1.762378E+00 +1.743081E-01 +1.408695E+00 +1.167464E-01 +1.553015E+00 +1.441475E-01 +1.399127E+00 +1.265146E-01 +1.632604E+00 +1.759006E-01 +1.464540E+00 +1.292863E-01 +1.555127E+00 +1.389518E-01 diff --git a/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/test.py b/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/test.py new file mode 100644 index 00000000000..8438ebe60cf --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/test.py @@ -0,0 +1,85 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell. This box is surrounded by a larger + box where the temperature field is not defined and has reflective + boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=lower_left[0], xmax=upper_right[0], + ymin=lower_left[1], ymax=upper_right[1], + zmin=lower_left[2], zmax=upper_right[2] + ) + cell_1 = openmc.Cell(fill=mat, region=-box) + + surrounding_box = openmc.model.RectangularParallelepiped( + xmin=-5.0, xmax=10.0, + ymin=-5.0, ymax=10.0, + zmin=-5.0, zmax=10.0, + boundary_type="reflective" + ) + cell_2 = openmc.Cell(fill=mat, region=(-surrounding_box & +box)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_1, cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/__init__.py b/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/inputs_true.dat b/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/inputs_true.dat new file mode 100644 index 00000000000..437cb53c46e --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/inputs_true.dat @@ -0,0 +1,64 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/results_true.dat b/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/results_true.dat new file mode 100644 index 00000000000..0a9618e2d94 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +3.814464E-01 1.239404E-02 +tally 1: +1.191064E+00 +8.295499E-02 +1.361519E+00 +1.163894E-01 +1.096751E+00 +7.178199E-02 +9.924760E-01 +6.328229E-02 +9.884224E-01 +6.060763E-02 +1.197255E+00 +7.868653E-02 +1.132349E+00 +7.940930E-02 +9.670501E-01 +6.000441E-02 diff --git a/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/test.py b/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/test.py new file mode 100644 index 00000000000..466e70ce0b9 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/test.py @@ -0,0 +1,85 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell. This box is surrounded by a larger + box where the temperature field is not defined and has vacuum + boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=lower_left[0], xmax=upper_right[0], + ymin=lower_left[1], ymax=upper_right[1], + zmin=lower_left[2], zmax=upper_right[2] + ) + cell_1 = openmc.Cell(fill=mat, region=-box) + + surrounding_box = openmc.model.RectangularParallelepiped( + xmin=-5.0, xmax=10.0, + ymin=-5.0, ymax=10.0, + zmin=-5.0, zmax=10.0, + boundary_type="vacuum" + ) + cell_2 = openmc.Cell(fill=mat, region=(-surrounding_box & +box)) + cell_2.temperature = 250.0 + + model.geometry = openmc.Geometry([cell_1, cell_2]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box((-5.0, -5.0, -5.0), (10.0, 10.0, 10.0)) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/single_cube/__init__.py b/tests/regression_tests/temperature_field/single_cube/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/single_cube/periodic_bc/__init__.py b/tests/regression_tests/temperature_field/single_cube/periodic_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/single_cube/periodic_bc/inputs_true.dat b/tests/regression_tests/temperature_field/single_cube/periodic_bc/inputs_true.dat new file mode 100644 index 00000000000..1d879a3adbd --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/periodic_bc/inputs_true.dat @@ -0,0 +1,57 @@ + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/single_cube/periodic_bc/results_true.dat b/tests/regression_tests/temperature_field/single_cube/periodic_bc/results_true.dat new file mode 100644 index 00000000000..bfed5a15158 --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/periodic_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.511081E+00 1.098748E-02 +tally 1: +4.208340E+01 +8.914019E+01 +4.317567E+01 +9.354453E+01 +4.113208E+01 +8.499925E+01 +4.247578E+01 +9.078685E+01 +4.215567E+01 +8.953567E+01 +4.294495E+01 +9.271687E+01 +4.152245E+01 +8.657952E+01 +4.253211E+01 +9.093756E+01 diff --git a/tests/regression_tests/temperature_field/single_cube/periodic_bc/test.py b/tests/regression_tests/temperature_field/single_cube/periodic_bc/test.py new file mode 100644 index 00000000000..5a8e25e7f52 --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/periodic_bc/test.py @@ -0,0 +1,79 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell - periodic boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + xmin = openmc.XPlane(x0=0.0, boundary_type="periodic") + xmax = openmc.XPlane(x0=5.0, boundary_type="periodic") + ymin = openmc.YPlane(y0=0.0, boundary_type="periodic") + ymax = openmc.YPlane(y0=5.0, boundary_type="periodic") + zmin = openmc.ZPlane(z0=0.0, boundary_type="periodic") + zmax = openmc.ZPlane(z0=5.0, boundary_type="periodic") + + xmin.periodic_surface = xmax + ymin.periodic_surface = ymax + zmin.periodic_surface = zmax + + cell = openmc.Cell(fill=mat, region=(+xmin & -xmax & +ymin & -ymax & +zmin & -zmax)) + model.geometry = openmc.Geometry([cell]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/single_cube/reflective_bc/__init__.py b/tests/regression_tests/temperature_field/single_cube/reflective_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/single_cube/reflective_bc/inputs_true.dat b/tests/regression_tests/temperature_field/single_cube/reflective_bc/inputs_true.dat new file mode 100644 index 00000000000..fba95023626 --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/reflective_bc/inputs_true.dat @@ -0,0 +1,57 @@ + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/single_cube/reflective_bc/results_true.dat b/tests/regression_tests/temperature_field/single_cube/reflective_bc/results_true.dat new file mode 100644 index 00000000000..7e16b60bd4b --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/reflective_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.507597E+00 1.811642E-02 +tally 1: +4.316452E+01 +9.417474E+01 +4.216471E+01 +8.960204E+01 +4.507663E+01 +1.029737E+02 +4.127791E+01 +8.578660E+01 +4.147002E+01 +8.693664E+01 +4.161294E+01 +8.721239E+01 +4.186858E+01 +8.850129E+01 +4.247133E+01 +9.097907E+01 diff --git a/tests/regression_tests/temperature_field/single_cube/reflective_bc/test.py b/tests/regression_tests/temperature_field/single_cube/reflective_bc/test.py new file mode 100644 index 00000000000..0e34b646f55 --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/reflective_bc/test.py @@ -0,0 +1,74 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell - reflective boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=lower_left[0], xmax=upper_right[0], + ymin=lower_left[1], ymax=upper_right[1], + zmin=lower_left[2], zmax=upper_right[2], + boundary_type="reflective" + ) + cell = openmc.Cell(fill=mat, region=-box) + model.geometry = openmc.Geometry([cell]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/temperature_field/single_cube/vacuum_bc/__init__.py b/tests/regression_tests/temperature_field/single_cube/vacuum_bc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/temperature_field/single_cube/vacuum_bc/inputs_true.dat b/tests/regression_tests/temperature_field/single_cube/vacuum_bc/inputs_true.dat new file mode 100644 index 00000000000..63e05c58f8b --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/vacuum_bc/inputs_true.dat @@ -0,0 +1,57 @@ + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + + + 1 + + + 1 + total + + + diff --git a/tests/regression_tests/temperature_field/single_cube/vacuum_bc/results_true.dat b/tests/regression_tests/temperature_field/single_cube/vacuum_bc/results_true.dat new file mode 100644 index 00000000000..05ff69dd728 --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/vacuum_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +5.071104E-02 2.927355E-03 +tally 1: +2.731604E+00 +4.072231E-01 +3.388567E+00 +6.168992E-01 +2.435725E+00 +3.283434E-01 +2.724446E+00 +4.158287E-01 +3.426649E+00 +6.625678E-01 +2.929812E+00 +4.792204E-01 +2.484602E+00 +3.679255E-01 +2.506827E+00 +3.400903E-01 diff --git a/tests/regression_tests/temperature_field/single_cube/vacuum_bc/test.py b/tests/regression_tests/temperature_field/single_cube/vacuum_bc/test.py new file mode 100644 index 00000000000..b3b1f25210a --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/vacuum_bc/test.py @@ -0,0 +1,74 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell - vacuum boundary conditions. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=lower_left[0], xmax=upper_right[0], + ymin=lower_left[1], ymax=upper_right[1], + zmin=lower_left[2], zmax=upper_right[2], + boundary_type="vacuum" + ) + cell = openmc.Cell(fill=mat, region=-box) + model.geometry = openmc.Geometry([cell]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + # Add tallies + mesh_filter = openmc.MeshFilter(mesh) + mesh_tally = openmc.Tally(name="total reaction rate") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["total"] + tallies = openmc.Tallies([mesh_tally]) + model.tallies = tallies + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() diff --git a/tests/regression_tests/weightwindows/survival_biasing/results_true.dat b/tests/regression_tests/weightwindows/survival_biasing/results_true.dat index b457a245c7c..8a4d3368cc4 100644 --- a/tests/regression_tests/weightwindows/survival_biasing/results_true.dat +++ b/tests/regression_tests/weightwindows/survival_biasing/results_true.dat @@ -1 +1 @@ -386e507008ed3c72c6e1a101aafc01cfaaff2c0b10555558d1eab6332441f7b17264b55002d721151bac2f3e7c1a8aac4c5ed243f5270533d171850f985206ed \ No newline at end of file +0d7b17d4e364bda2be21feb2e5b2aae4be69fc8ed634c4f45062e8efc61b7bd24dcf448837692fd603afe3756501fe8821d134f2d6289179618f85178ddb8793 \ No newline at end of file diff --git a/tests/unit_tests/test_temperature_field.py b/tests/unit_tests/test_temperature_field.py new file mode 100644 index 00000000000..a99337eb1c3 --- /dev/null +++ b/tests/unit_tests/test_temperature_field.py @@ -0,0 +1,40 @@ +"""Test the temperature field Python object.""" + +import openmc +import pytest +import numpy as np + + +@pytest.fixture(scope="module") +def mesh(): + """2x2x2 regular mesh""" + dim = 2 + lower_left = (0., 0., 0.) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + return mesh + + +def test_xml_serialization(mesh, run_in_tmpdir): + """Test XMl serialization for a temperature field declaration. + + """ + # Define values + values = [0., 1., 2., 3., 4., 5., 6., 7.] + + # Create field + temperature_field = openmc.TemperatureField(mesh, values) + + # Export settings + settings = openmc.Settings() + settings.temperature_field = temperature_field + settings.export_to_xml() + + # Read settings from xml file + read_settings = openmc.Settings.from_xml() + + # Check consistency + assert read_settings.temperature_field == temperature_field