diff --git a/CMakeLists.txt b/CMakeLists.txt index 2c2cc72e0ae..a59843ef265 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -341,10 +341,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 @@ -399,6 +401,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 87ada334d13..48ee10bd831 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -1027,3 +1027,9 @@ Functions :type scores: const int* :return: Return status (negative if an error occurred) :rtype: int + +.. doxygenfunction:: openmc_temperature_field_set_temperature + +.. doxygenfunction:: openmc_temperature_field_size + +.. doxygenfunction:: openmc_temperature_field_get_value diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 01021e8d22d..e7f3eaba4d4 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -236,6 +236,70 @@ from MF=1, MT=458 on an ENDF evaluation. *Default*: true +.. _dnp_drift_element: + +----------------------- +```` Element +----------------------- + +The ```` element configures the delayed neutron precursor (DNP) drift +feature, which models the transport of precursors along fluid velocity streamlines +in liquid-fueled systems (see :ref:`dnp_drift` for a detailed description of the +underlying methodology). It has the following sub-elements: + + :velocity_field: + The integer ID of a 3D vector field (defined via a :ref:`field_element`) that + provides the fluid velocity used to transport precursor sites. The referenced + field must be of type **velocity** and its values must have three components + per node or element (:math:`v_x, v_y, v_z`). + + :boundary_map: + A group of sub-elements that assigns physical group numbers to boundary + types. Each sub-element contains a whitespace-separated list of integer + physical group numbers: The following sub-elements are recognized: + + :inlet: + Physical group numbers corresponding to inlet boundaries. + + :outlet: + Physical group numbers corresponding to outlet boundaries. + + :wall: + Physical group numbers corresponding to wall boundaries. + + :physical_group_map: + A group of sub-elements that maps individual boundary faces on the velocity + mesh to their physical group number. Because the DNP drift feature is currently + limited to OpenMC regular meshes, this mapping must be provided explicitly + by the user. It contains two sub-elements: + + :face_ids: + A whitespace-separated list of integer face IDs identifying each boundary + face on the mesh. + + :physical_groups: + A whitespace-separated list of integer physical group numbers, one per face, + in the same order as ``face_ids``. Each entry assigns the corresponding face to + a physical group. + + The two lists must have the same length. + + :integrator: + The time integration scheme used to transport precursors along streamlines. + Currently, only **RK4** (fourth-order Runge-Kutta) is supported. + + :integrator_dt: + The time step size (in seconds) used by the integrator. + + :recycling: + Whether to enable precursor recycling. If **true**, precursors that exit through + an outlet can be re-injected at an inlet if they do not decay during the external + transit. If **false**, precursors leaving through an outlet are immediately discarded. + + :external_travel_time: + The average time (in seconds) for a precursor to transit out of the modeled system + before re-entering. Only used when recycling is set to **true**. + -------------------------------- ```` Element -------------------------------- @@ -275,6 +339,51 @@ history-based parallelism. *Default*: false +.. _field_element: + +------------------------------- +```` Element +------------------------------- + +The ```` element describes a field combining a geometric mesh with values +for each node or element depending on the type of mapping. Fields are currently +used for the delayed neutron precursor drift feature (velocity field) and to +declare temperature using a mesh (temperature field). It has the following +attributes/sub-elements: + + :id: + A unique integer that is used to identify the field. + + :type: + The type of field. Currently, **temperature** and **velocity** are the two + field types supported. + + :mesh: + The integer ID used to identify the associated mesh. + + :mapping: + The type of mapping between the mesh and the values. Mapping can be either + **nodal** (values defined at mesh nodes) or **cell** (values defined per mesh + element and assumed constant within each element). + + *Default*: 'cell' + + :values: + A whitespace-separated list of floating point numbers specifying the field + data. Values are listed sequentially, with all components of a given node + or element grouped together before moving to the next (e.g., for a + three-component field: + :math:`[v_{1,1}, v_{1,2}, v_{1,3}, v_{2,1}, v_{2,2}, v_{2,3}, \ldots]`, + where the first index identifies the node/element and the second identifies + the component). The total number of values must equal + :math:`N_{\text{c}} \cdot N_{\text{l}}`, where + :math:`N_{\text{c}}` is the number of components per location + (e.g., 1 for a scalar field, 3 for a vector field) and + :math:`N_{\text{l}}` is the number of mesh nodes (for nodal mapping) + or mesh elements (for cell mapping). For a regular mesh with dimensions + :math:`(n_x, n_y, n_z)`, the number of nodes is :math:`(n_x+1)(n_y+1)(n_z+1)` + and the number of elements is :math:`n_x\cdot n_y\cdot n_z`. + -------------------------------- ```` Element -------------------------------- @@ -1494,6 +1603,18 @@ a material default temperature. *Default*: 293.6 K +.. _temperature_field: + +------------------------------- +```` Element +------------------------------- + +The ```` element specifies the ID of a scalar field (defined via +a :ref:`field_element`) that provides temperatures on a geometric mesh. The referenced +field must be of type **temperature** and its values must be scalar (one component per +node or element) and given in Kelvin. When a temperature field is provided, it is used +to assign temperatures in the model based on the underlying mesh. + .. _temperature_method: -------------------------------- diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 3ffcac69cda..654a262950d 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -172,6 +172,28 @@ Meshes openmc.SphericalMesh openmc.UnstructuredMesh +Fields +------ + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclassinherit.rst + + openmc.FieldBase + openmc.VelocityField + openmc.TemperatureField + +DNP drift +--------- + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + openmc.DNPDrift + Geometry Plotting ----------------- diff --git a/docs/source/usersguide/dnp_drift.rst b/docs/source/usersguide/dnp_drift.rst new file mode 100644 index 00000000000..bc1030e9bec --- /dev/null +++ b/docs/source/usersguide/dnp_drift.rst @@ -0,0 +1,251 @@ +.. _dnp_drift: + +Delayed Neutron Precursor (DNP) Drift +===================================== + +.. note:: + The DNP drift feature in OpenMC is currently limited to regular meshes + (:class:`openmc.RegularMesh`). + +.. warning:: + Because DNP removal is an analog process, it is not taken into account with + collision or track-length processes. One way to correct the averaged + :math:`k_{\text{eff}}` is to normalize it by the ratio of analog and track-length + nu-fission tallies. + +In solid-fueled reactors, delayed neutron precursors (DNPs) are produced and decay at the +same spatial location where the fission event occurred. In liquid-fueled reactors (e.g., +molten salt reactors (MSRs)), the fissile material is dissolved in a carrier fluid that +circulates through the primary loop. As a result, DNPs are carried away from their point +of origin by the flowing fuel and may emit their delayed neutrons outside of the active +core region. The spatial redistribution of precursors alters the delayed neutron source +distribution, which can significantly affect reactor criticality and kinetics. This +phenomenon is commonly referred to as **DNP drift**. + +OpenMC provides capabilities to model DNP drift by explicitly transporting precursors along +fluid velocity streamlines. This feature is designed to support high-fidelity multiphysics +coupling with computational fluid dynamics (CFD) codes such as +`nekRS `_ through the +`Cardinal `_ framework. + +Streamline integration +---------------------- + +In OpenMC, individual delayed neutron precursors are advected along fluid streamlines +until they reach their sampled decay time. The decay time for a precursor belonging to +delayed group :math:`g` is sampled as: + +.. math:: + + t_{\text{decay}} = - \frac{\log(\zeta)}{\lambda_g}, + +where :math:`\zeta` is a randomly generated number and :math:`\lambda_g` is the decay +constant corresponding to the sampled delayed neutron group number :math:`g`. + +Each precursor trajectory is obtained by integrating: + +.. math:: + + \frac{dy}{dt} = v(y,t), + +where :math:`y(t)` is the precursor position and :math:`v(y,t)` is the local fluid velocity +interpolated from the underlying mesh. In the current implementation, the velocity field +is assumed to be static within an OpenMC calculation, so the time dependence is dropped +and :math:`v(y,t)` reduces to :math:`v(y)`. The integration is performed using the classical +fourth-order Runge-Kutta method (RK4). Given a precursor at position :math:`x_n` and +a time step :math:`\Delta t`, the method evaluates four intermediate velocity estimates: + +.. math:: + + k_1 = v(y_n), \\ + k_2 = v\left(y_n + \frac{\Delta t}{2}\,k_1\right), \\ + k_3 = v\left(y_n + \frac{\Delta t}{2}\,k_2\right), \\ + k_4 = v\left(y_n + \Delta t\,k_3\right). + +The position is then advanced as: + +.. math:: + + y_{n+1} = y_n + \frac{\Delta t}{6}\left(k_1 + 2\,k_2 + 2\,k_3 + k_4\right). + +Boundary conditions and recycling +--------------------------------- + +In practice, the velocity mesh used for DNP transport typically covers only a portion +of the full primary loop (e.g., the active core region). Some surfaces of the mesh +therefore represent the inlet and outlet of the modeled open system. +These boundary faces are organized into physical groups, following conventions used +by mesh generators and CFD codes. The user must specify which physical groups correspond +to **inlet**, **outlet**, or **wall** boundaries. Additionally, because the DNP drift feature +is currently limited to OpenMC regular meshes (:class:`openmc.RegularMesh`), the assignment +of individual faces to physical groups must also be provided explicitly by the user. + +During streamline integration, a precursor may encounter one of these boundaries. +Each type is handled as follows: + +- **Outlet**: The precursor exits the modeled domain and is discarded. + +- **Wall**: The precursor is stopped at the position where the wall crossing would have + occurred and decays at that location. + +- **Inlet**: The precursor is repositioned at a location sampled uniformly on the set + of inlet faces and continues its trajectory. + +In a real MSR, however, the fuel salt recirculates: it exits the core, passes through +heat exchangers and piping, and re-enters the core at the inlet. Precursors leaving +the core might re-enter if they survive long enough. The **recycling** option +models this recirculation without requiring the velocity mesh to cover the entire +primary loop. It assumes a fixed external travel time :math:`\tau_{\text{ext}}` +representing the transit duration through the portion of the loop that is not +explicitly represented with the velocity mesh. + +When recycling is enabled, the algorithmic behavior associated with the outlet boundary +condition is modified: instead of being immediately discarded, a precursor that reaches +an outlet has its remaining lifetime :math:`t_{\text{rem}}` compared to +:math:`\tau_{\text{ext}}`: + +- If :math:`t_{\text{rem}} \gt \tau_{\text{ext}}`, the precursor survives the + external loop. It is re-injected at a randomly sampled inlet position with an updated + remaining lifetime of :math:`t_{\text{rem}} - \tau_{\text{ext}}`, and streamline + integration resumes from the new inlet location. + +- If :math:`t_{\text{rem}} \le \tau_{\text{ext}}`, the precursor decays during + the external transit and is discarded. + +This cycle of exit, external transit, and re-injection may repeat multiple times for +long-lived precursors, allowing them to complete several passes through the core +before eventually decaying. + +Geometric reconciliation +------------------------ + +After transport of each DNP, a **reconciliation step** verifies that the +precursor site lies within the OpenMC geometry. Because DNP transport is +performed on a dedicated mesh while particle transport uses the native OpenMC +geometry (CSG or DAGMC), a precursor may be transported to a position that the +geometry routines classify as outside the model. + +Note that this inside/outside classification depends on the **neutron emission direction** +sampled during the original fission event, not the fluid velocity vector. A site physically +located at a boundary surface may be seen as inside or outside the OpenMC geometry depending +on whether this direction points inward or outward relative to the surface normal. + +The reconciliation procedure distinguishes two cases: + +1. **Site classified as inside the geometry** - OpenMC's cell-finding routines successfully + locate the precursor inside a cell (including sites on a boundary surface whose + emission direction points inward). The site is accepted as a valid delayed neutron source + location for the next generation. + +2. **Site not found inside any cell** - The precursor is nudged slightly backward + along its emission direction and the cell search is repeated: + + - **Site on a boundary surface (heading outward)**: If the nudged position is + found inside a cell and the nearest boundary is within the nudge tolerance, + the precursor is determined to be on a surface. The boundary condition associated + with that surface is then applied: + + - *Reflective / White*: the site position and direction are transformed according + to the boundary condition. If the boundary carries an albedo factor, the site + weight is adjusted accordingly. + - *Periodic*: the site is mapped to the corresponding periodic partner surface. + - *Vacuum*: the site is discarded. + + - **Site genuinely outside**: If the nudged position is still outside the geometry, + or the nearest boundary is closer than the nudge tolerance, the site is + discarded. A warning is emitted so that the user can diagnose potential + issues with mesh alignment. + +.. tip:: + A high number of lost sites during reconciliation might indicate a mismatch between + the velocity mesh and the OpenMC geometry, or an integrator time step that is too large. + Reducing ``integrator_dt`` or improving the geometric alignment between the mesh and + the model can mitigate this. + +Usage +----- + +The following example shows how to set up DNP drift on a simple 2x2x2 regular mesh. + +**1. Define the mesh and the velocity field** + +.. code-block:: python + + import numpy as np + import openmc + + # Create a 2x2x2 regular mesh + mesh = openmc.RegularMesh() + mesh.lower_left = (0.0, 0.0, 0.0) + mesh.upper_right = (10.0, 10.0, 10.0) + mesh.dimension = (2, 2, 2) + + # Define a uniform velocity field on the mesh nodes + # For a 2x2x2 mesh there are 3x3x3 = 27 unique nodes, each with 3 components. + velocity_field = openmc.VelocityField( + mesh=mesh, + values=np.array([1.0]*3*27).reshape(-1,3), + mapping="nodal" + ) + +**2. Create the DNPDrift object** + +.. code-block:: python + + dnp_drift = openmc.DNPDrift( + velocity_field=velocity_field, + boundary_map={ + "inlet": [1], + "outlet": [2], + "wall": [3] + }, + physical_group_map={ + "face_ids": [ + 0, 24, 12, 36, 7, 31, 19, 43, 15, 39, 21, 45, 2, + 26, 8,32, 4, 16, 10, 22, 29, 41, 35, 47 + ], + "physical_groups": [ + 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3 + ]}, + integrator="RK4", + integrator_dt=0.1, + recycling=True, + external_travel_time=2.0 + ) + +The parameters are explained below: + +.. list-table:: + :header-rows: 1 + :widths: 30 70 + + * - Parameter + - Description + * - ``velocity_field`` + - The velocity field used to transport precursor sites. + * - ``boundary_map`` + - Dictionary mapping boundary type names (``inlet``, ``outlet``, + ``wall``) to lists of physical group numbers. + * - ``physical_group_map`` + - Dictionary with two keys: ``face_ids`` (integer IDs of each boundary + face on the mesh) and ``physical_groups`` (the corresponding physical + group number for each face). + * - ``integrator`` + - Time integration scheme. Currently, only ``"RK4"`` is supported. + * - ``integrator_dt`` + - Time step size (in seconds) used by the integrator. + * - ``recycling`` + - If ``True``, precursors that exit through an outlet can be re-injected at + an inlet if they do not decay during the external travel. Otherwise, precursors + leaving by an outlet are considered lost. + * - ``external_travel_time`` + - Average time (in seconds) needed for a precursor to re-enter the modeled system + after having left the system by crossing an outlet. Only used when ``recycling=True``. + +**3. Register the DNP drift model in the settings** + +.. code-block:: python + + settings = openmc.Settings() + settings.dnp_drift = dnp_drift diff --git a/docs/source/usersguide/index.rst b/docs/source/usersguide/index.rst index 5f8e0197e70..9d1b1b60e3c 100644 --- a/docs/source/usersguide/index.rst +++ b/docs/source/usersguide/index.rst @@ -23,6 +23,7 @@ essential aspects of using OpenMC to perform simulations. depletion decay_sources kinetics + dnp_drift scripts processing parallel diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index 8d1bad54f35..5a851a900d8 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -265,6 +265,38 @@ attribute, e.g., :attr:`Material.temperature` or :attr:`Cell.temperature` attributes, respectively. +Alternatively, temperatures can be specified using a temperature field composed +of a mesh and a map that associates a temperature value to each mesh element. +During the simulation, temperatures associated with particles are also updated +every time a temperature field mesh surface is crossed. While a particle is +contained inside the temperature mesh, temperatures from the temperature field +take precedence over any temperature declared for a cell, a material, or +globally. + +The following example shows how to specify temperatures using a temperature +field based on a regular mesh: + +.. code-block:: python + + # Define a 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 mesh element + 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/boundary_condition.h b/include/openmc/boundary_condition.h index 89bf7ca8bb3..5e5037ab9c9 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -21,12 +21,38 @@ class BoundaryCondition { public: virtual ~BoundaryCondition() = default; + //! Apply the geometric transformation. + //! \param[in] p Particle + //! \param[in] surf Current surface + //! \param[out] new_r Transformed position + //! \param[out] new_u Transformed direction + //! \param[out] new_surface New surface index + virtual void transform(Particle& p, const Surface& surf, Position& new_r, + Direction& new_u, int& new_surface) const = 0; + + //! Handle the crossing of the boundary condition by the particle + //! \param[inout] p Particle + //! \param[in] surf Current surface + //! \param[in] new_r Transformed position + //! \param[in] new_u Transformed direction + //! \param[in] new_surface New surface index + virtual void cross_bc(Particle& p, const Surface& surf, const Position& new_r, + const Direction& new_u, int new_surface) const = 0; + //! Perform tracking operations for a particle that strikes the boundary. //! \param p The particle that struck the boundary. This class is not meant //! to directly modify anything about the particle, but it will do so //! indirectly by calling the particle's appropriate cross_*_bc function. //! \param surf The specific surface on the boundary the particle struck. - virtual void handle_particle(Particle& p, const Surface& surf) const = 0; + virtual void handle_particle(Particle& p, const Surface& surf) const + { + Position new_r; + Direction new_u; + int new_surface; + transform(p, surf, new_r, new_u, new_surface); + BoundaryCondition::handle_albedo(p, surf); + cross_bc(p, surf, new_r, new_u, new_surface); + } //! Modify the incident particle's weight according to the boundary's albedo. //! \param p The particle that struck the boundary. This function calculates @@ -64,6 +90,9 @@ class BoundaryCondition { //! Return if this BC has an albedo. bool has_albedo() const { return (albedo_ > 0.0); } + //! Return albedo value. + double albedo() const { return albedo_; } + private: double albedo_ = -1.0; }; @@ -74,6 +103,16 @@ class BoundaryCondition { class VacuumBC : public BoundaryCondition { public: + // Not used - handle_particle is fully overriden + void transform(Particle& p, const Surface& surf, Position& new_r, + Direction& new_u, int& new_surface) const override + {} + + // Not used - handle_particle is fully overriden + void cross_bc(Particle& p, const Surface& surf, const Position& new_r, + const Direction& new_u, int new_surface) const override + {} + void handle_particle(Particle& p, const Surface& surf) const override; std::string type() const override { return "vacuum"; } @@ -85,7 +124,11 @@ class VacuumBC : public BoundaryCondition { class ReflectiveBC : public BoundaryCondition { public: - void handle_particle(Particle& p, const Surface& surf) const override; + void transform(Particle& p, const Surface& surf, Position& new_r, + Direction& new_u, int& new_surface) const override; + + void cross_bc(Particle& p, const Surface& surf, const Position& new_r, + const Direction& new_u, int new_surface) const override; std::string type() const override { return "reflective"; } }; @@ -96,7 +139,11 @@ class ReflectiveBC : public BoundaryCondition { class WhiteBC : public BoundaryCondition { public: - void handle_particle(Particle& p, const Surface& surf) const override; + void transform(Particle& p, const Surface& surf, Position& new_r, + Direction& new_u, int& new_surface) const override; + + void cross_bc(Particle& p, const Surface& surf, const Position& new_r, + const Direction& new_u, int new_surface) const override; std::string type() const override { return "white"; } }; @@ -128,7 +175,11 @@ class TranslationalPeriodicBC : public PeriodicBC { public: TranslationalPeriodicBC(int i_surf, int j_surf); - void handle_particle(Particle& p, const Surface& surf) const override; + void transform(Particle& p, const Surface& surf, Position& new_r, + Direction& new_u, int& new_surface) const override; + + void cross_bc(Particle& p, const Surface& surf, const Position& new_r, + const Direction& new_u, int new_surface) const override; protected: //! Vector along which incident particles will be moved @@ -145,9 +196,12 @@ class RotationalPeriodicBC : public PeriodicBC { public: enum PeriodicAxis { x, y, z }; RotationalPeriodicBC(int i_surf, int j_surf, PeriodicAxis axis); - double compute_periodic_rotation( - double rise_1, double run_1, double rise_2, double run_2) const; - void handle_particle(Particle& p, const Surface& surf) const override; + + void transform(Particle& p, const Surface& surf, Position& new_r, + Direction& new_u, int& new_surface) const override; + + void cross_bc(Particle& p, const Surface& surf, const Position& new_r, + const Direction& new_u, int new_surface) const override; protected: //! Angle about the axis by which particle coordinates will be rotated diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 14c1ee26cf2..3036b4237e3 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -272,6 +272,48 @@ int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides); int openmc_tally_set_scores(int32_t index, int n, const char** scores); int openmc_tally_set_type(int32_t index, const char* type); int openmc_tally_set_writable(int32_t index, bool writable); + +//! Set a value at a given index in the temperature field +// +//! \param[in] index Index in the field +//! \param[in] temperature Temperature value to apply +//! \return Status (negative if an error occurred) +int openmc_temperature_field_set_temperature(int32_t index, double temperature); + +//! Return the temperature field size (number of cells) +// +//! \return Number of cell in the temperature field +size_t openmc_temperature_field_size(); + +//! Get the value at a given index in the temperature field +// +//! \param[in] index Index in the field +//! \param[out] temperature Temperature value at the given index +//! \return Status (negative if an error occurred) +int openmc_temperature_field_get_value(int32_t index, double* temperature); + +//! Set a value at a given index in the velocity field +// +//! \param[in] index Index in the field +//! \param[in] vx X-axis velocity component to apply +//! \param[in] vy Y-axis velocity component to apply +//! \param[in] vz Z-axis velocity component to apply +//! \return Status (negative if an error occurred) +int openmc_velocity_field_set_velocity( + int32_t index, double vx, double vy, double vz); + +//! Return the velocity field size (number of cells) +// +//! \return Number of cell in the velocity field +size_t openmc_velocity_field_size(); + +//! Get the value at a given index in the velocity field +// +//! \param[in] index Index in the field +//! \param[out] velocity Velocity value at the given index +//! \return Status (negative if an error occurred) +int openmc_velocity_field_get_value(int32_t index, double velocity[3]); + int openmc_get_weight_windows_index(int32_t id, int32_t* idx); int openmc_weight_windows_get_id(int32_t index, int32_t* id); int openmc_weight_windows_set_id(int32_t index, int32_t id); diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 46ba17c603b..8aeffa7b9d0 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -406,6 +406,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..7b435245f36 --- /dev/null +++ b/include/openmc/dnp_drift.h @@ -0,0 +1,75 @@ +#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 { REENTER, ESCAPE, DECAY_IN_PLACE }; + +//! 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[in,out] 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[in,out] 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[in,out] site Fission site corresponding to the DNP +//! \param[in] decay_time Sampled decay time (in seconds) +//! \param[in] seed Random number generator seed +//! \return True if the site is considered inside the explicitly modeled part of +//! the system. False otherwise. +bool transport_dnp(SourceSite& site, double decay_time, uint64_t* seed); + +//! Reconcile a delayed neutron precursor source site with the model geometry +//! after mesh-based drift. +//! +//! After transporting a delayed neutron precursor on a mesh, the emission site +//! may land outside the model or on a geometry boundary. This function +//! identifies if the site is valid and applies boundary condition corrections +//! as needed. If the function returns true, the site is guaranteed to be usable +//! for the next fission generation. +//! +//! \param[in,out] site Fission site corresponding to the DNP +//! \param[in] particle_id Particle ID used to initialize the RNG seed for the +//! boundary condition sampling (e.g., white BC). +//!\return True if the site is usable for the next generation, false if the site +//! should be discarded. +bool reconcile_precursor_drift(SourceSite& site, int64_t particle_id); + +} // namespace openmc + +#endif // OPENMC_DNP_DRIFT_H diff --git a/include/openmc/event.h b/include/openmc/event.h index 055d5230f29..914bba4ed1d 100644 --- a/include/openmc/event.h +++ b/include/openmc/event.h @@ -109,6 +109,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..072cfd39ddc --- /dev/null +++ b/include/openmc/field.h @@ -0,0 +1,356 @@ +#ifndef OPENMC_FIELD_H +#define OPENMC_FIELD_H + +#include +#include +#include + +#include "fmt/format.h" +#include "pugixml.hpp" + +#include "openmc/mesh.h" +#include "openmc/position.h" +#include "openmc/vector.h" + +namespace openmc { + +class Field; + +namespace model { + +//! Global container of all field instances, indexed by position. +extern vector> fields; + +//! Map from user-facing field ID to index in fields +extern std::unordered_map field_map; + +} // namespace model + +// ----------------------------------------------------------- +// Field +// ----------------------------------------------------------- + +//! Describes how field values are mapped onto a mesh +enum class FieldMapping { + NODAL, // Nodal representation (values defined for each vertex) + CELL // Cell-based representation (values defined for each cell) +}; + +//! Abstract base class for all fields defined on a geometric mesh. +// +//! Provides the mapping type and a non-owning pointer to the associated mesh. +//! Concrete field behaviour is implemented by TypedField and its subclasses. +class Field { +public: + explicit Field(Mesh* mesh, FieldMapping mapping) + : mesh_(mesh), mapping_(mapping) + { + if (mesh_ == nullptr) + fatal_error("Field requires a non-null mesh pointer!"); + } + + virtual ~Field() = default; + + Mesh* mesh_ptr() const { return mesh_; } + FieldMapping mapping() const { return mapping_; } + +protected: + Mesh* mesh_; //!< Non-owning pointer to the geometric mesh + FieldMapping mapping_; //!< Relationship between values and mesh entities + //!< (nodes vs. cells) +}; + +// ----------------------------------------------------------- +// TypedField +// ----------------------------------------------------------- + +//! A field that stores typed data on a mesh with a specified mapping. +//! +//! Provides evaluation (including trilinear interpolation for nodal fields), +//! assignment, and mesh-boundary queries. +//! +//! \tparam T Value type stored per mesh entity (e.g., double for scalar fields, +//! Direction for vector fields). +template +class TypedField : public Field { +public: + //! Construct a fully initialized TypedField. + // + //! \param[in] mesh_ptr Non-owning pointer to the mesh + //! \param[in] values Field values. + //! \param[in] mapping Mapping type: 'nodal' or 'cell' + TypedField(Mesh* mesh_ptr, vector values, FieldMapping mapping) + : Field(mesh_ptr, mapping) + { + int expected = 0; + std::string element_type; + switch (mapping_) { + case FieldMapping::NODAL: + expected = mesh_->n_vertices(); + element_type = "unique vertices"; + break; + case FieldMapping::CELL: + expected = mesh_->n_bins(); + element_type = "elements"; + break; + default: + fatal_error("Logic not implemented!"); + } + + if (static_cast(expected) != values.size()) { + fatal_error(fmt::format( + "The number of {} in the mesh ({}) is not consistent with the " + "number of values declared for this field ({})!", + element_type, expected, values.size())); + } + + data_ = std::move(values); + } + + //! Assign a value to a bin. + // + //! \param[in] bin Bin number + //! \param[in] value Value to store + void assign(int bin, T value) + { + if (bin < 0 || bin >= static_cast(data_.size())) { + fatal_error(fmt::format( + "Bin index {} is out of range [0, {}].", bin, data_.size())); + } + data_[bin] = value; + } + + //! Return the mesh bin associated with a given position. + // + //! \param[in] r Position + //! \return Bin number + int get_bin(const Position& r) const { return mesh_ptr()->get_bin(r); } + + //! Evaluate field at a given position inside the mesh knowing the current + //! bin. + // + //! \param[in] r Position + //! \param[in] bin Bin number corresponding to r + //! \return Value corresponding to r + T evaluate_in_mesh(const Position& r, int bin) const + { + if (bin == C_NONE) + fatal_error("Bin outside the mesh."); + + switch (mapping()) { + case FieldMapping::NODAL: + return trilinear_interpolation(r, bin); + case FieldMapping::CELL: + return value(bin); + default: + fatal_error("Not implemented for this mapping type!"); + break; + } + } + + //! Evaluate field at a given position knowing the previous position and bin. + // + //! The bin used to evaluate the field at r1 can either be: + //! - the bin corresponding to r1, if r1 is inside the mesh, + //! - the bin corresponding to the last crossed bin inside the mesh during ray + //! traversal analysis between r0 and r1, if r1 is outside the mesh, + //! - bin0, if r1 is outside the mesh and no bins were found during ray + //! traversal. + // + //! \param[in] r0 Previous position (inside the mesh) + //! \param[in] r1 Current position + //! \param[in] bin0 Bin corresponding to r0 + //! \return Value corresponding to r1 relative to a clamped bin + T evaluate_clamped(const Position& r0, const Position& r1, int bin0) const + { + int next_bin = mesh_ptr()->get_bin_clamped(r0, r1, bin0); + return evaluate_in_mesh(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. + // + //! The position r and the bin number does not have to be related. + // + //! \param[in] r Position + //! \param[in] bin Bin number + //! \return Interpolated value + T trilinear_interpolation(const Position& r, int bin) const + { + // Normalize coordinates + Position n_r = mesh_ptr()->normalize_coordinates(r, bin); + + // Protect from 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); + if (v.size() != 8) { + fatal_error(fmt::format( + "Trilinear interpolation requires 8 vertices per cell, got {}", + v.size())); + } + + // 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 Distance to the next mesh boundary + double distance_to_next_boundary( + int current_bin, const Position& r, const Direction& u, int& bin_next) const + { + return mesh_ptr()->distance_to_next_boundary(current_bin, r, u, bin_next); + } + + // Data field value accessors + T value(int i) const + { + if (i < 0 || i >= static_cast(data_.size())) { + fatal_error( + fmt::format("Data index {} is out of range [0, {}].", i, data_.size())); + } + return data_[i]; + } + const vector& values() const { return data_; } + +private: + vector data_; //!< Data associated with the mesh +}; + +// ----------------------------------------------------------- +// TemperatureField +// ----------------------------------------------------------- + +class TemperatureField : public TypedField { +public: + // Constructors + TemperatureField(Mesh* mesh_ptr, vector values, + FieldMapping mapping = FieldMapping::CELL) + : TypedField(mesh_ptr, std::move(values), mapping) + {} + + //! Returns the temperature in Kelvin corresponding to a given bin number + //! relative to the mesh. + // + //! \param[in] bin Bin number + //! \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 TypedField { +public: + // Constructors + VelocityField(Mesh* mesh_ptr, vector values, FieldMapping mapping) + : TypedField(mesh_ptr, std::move(values), 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 + 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 +}; + +template +T* get_field(int field_id) +{ + auto it = model::field_map.find(field_id); + if (it == model::field_map.end()) { + fatal_error( + fmt::format("No element with id={} was found.", field_id)); + } + + Field* base = model::fields[it->second].get(); + auto* field = dynamic_cast(base); + if (!field) { + fatal_error(fmt::format("Field id={} is of incompatible type.", field_id)); + } + + return field; +} + +void read_fields(const pugi::xml_node& root); + +} // namespace openmc + +#endif // OPENMC_FIELD_H diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 54c9e7b4195..1f11530b898 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -145,6 +145,9 @@ class MaterialVolumes { //! Base mesh class //============================================================================== +// Physical groups map type +using PGMap = std::unordered_map>; + class Mesh { public: // Constructors and destructor @@ -190,6 +193,80 @@ class Mesh { virtual void surface_bins_crossed( Position r0, Position r1, const Direction& u, vector& bins) const = 0; + //! Determine which bins and surface bins were crossed by a particle. + // + //! Surface bins are identified without the inward flag. + // + //! \param[in] r0 Previous position + //! \param[in] r1 Current position + //! \param[out] leaving_surface_ids Surface bins that were left + //! \param[out] entering_surface_ids Surface bins that were entered + //! \param[out] bins Bins that were crossed + //! \param[out] length_fractions Fraction of tracklength in each bin + virtual void bins_and_surface_bins_crossed(Position r0, Position r1, + vector& leaving_surface_ids, vector& entering_surface_ids, + vector& bins, vector& length_fractions) const + { + fatal_error("Not implemented."); + } + + //! Get bin at a given position r1 if r1 is inside the mesh. If r1 is located + //! outside the mesh, this function returns the nearest valid bin inside the + //! mesh, if bins were crossed during ray traversal between r0 and r1, or + //! bin0. + // + //! \param[in] r0 Previous position + //! \param[in] r1 Current position + //! \param[in] bin0 Bin corresponding to r0 + //! \return If r1 is inside the mesh: bin corresponding to r1, + //! if r1 outside the mesh: last bin inside the mesh when travelling from r0 + //! to r1 or bin0 if no bins were crossed. + virtual int get_bin_clamped(Position r0, Position r1, int bin0) const + { + fatal_error("Not implemented."); + } + + //! Sample a position on mesh faces belonging to one or more physical groups + // + //! \param[out] p Sampled position + //! \param[out] bin Bin number corresponding to r + //! \param[in] seed Random number generator seed + //! \param[in] physical_groups Physical groups considered + virtual void sample_on_physical_groups( + Position& p, int& bin, uint64_t* seed, vector physical_groups) + { + fatal_error("Not implemented!"); + } + + //! Normalize coordinates by scaling them to the dimensions of a given bin + // + //! \param[in] r Position + //! \param[in] bin Bin + //! \return Normalized position + virtual Position normalize_coordinates(const Position& r, int bin) + { + 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 @@ -202,6 +279,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_; } @@ -300,6 +380,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 { @@ -338,12 +437,28 @@ 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 bins_and_surface_bins_crossed(Position r0, Position r1, + vector& leaving_surface_ids, vector& entering_surface_ids, + vector& bins, vector& length_fractions) const override; + + double distance_to_next_boundary( + int current_bin, Position r, Direction u, int& bin_next) const override; + + int get_bin_clamped(Position r0, Position r1, int bin0) const 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 @@ -430,6 +545,21 @@ 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 + { + fatal_error("Not implemented"); + return -1.0; + } + //! Get a label for the mesh bin std::string bin_label(int bin) const override; @@ -508,6 +638,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; @@ -517,6 +649,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; @@ -547,6 +682,36 @@ class RegularMesh : public StructuredMesh { int set_grid(); + //! Calculate the area of a face given its ID number + // + //! \param[in] face_id Face ID number (without inward information) + //! \return Area of the face + double face_area(int face_id); + + void sample_on_physical_groups(Position& pa, int& cell, uint64_t* seed, + vector physical_groups) override; + + //! Sample a position on a face given its ID number + // + //! \param[in] face_id Face ID number (without inward information) + //! \param[in] seed Random number generator seed + //! \return Sampled position + Position sample_on_face(int face_id, uint64_t* seed); + + //! Return the unique ID of a vertex + // + //! The vertex unique ID system ensures that vertices from different bins + //! located at the same position in space share the same unique ID. + // + //! \param[in] ijk Array of mesh indices + //! \param[in] local_verted_idx Local index of a vertex relative to a bin + //! \return Vertex unique ID + int return_vertex_unique_id(MeshIndex ijk, int local_vertex_idx) const; + + std::vector connectivity(int id) const override; + + Position normalize_coordinates(const Position& r, int bin) override; + // Data members double volume_frac_; //!< Volume fraction of each mesh element double element_volume_; //!< Volume of each mesh element @@ -570,6 +735,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; @@ -746,6 +914,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; // Overridden Methods @@ -784,23 +955,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 24b7eb53c0a..25d30cab374 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -269,6 +269,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. */ @@ -416,6 +434,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 @@ -441,12 +465,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_; @@ -574,6 +601,8 @@ class ParticleData : public GeometryState { int64_t n_progeny_ {0}; + TransportEvent next_event_ = TransportEvent(); + public: //---------------------------------------------------------------------------- // Constructors @@ -793,6 +822,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 3ca7e9c3a82..6be7848658b 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -59,6 +59,7 @@ extern bool create_delayed_neutrons; //!< create delayed fission neutrons? extern 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 bool entropy_on; //!< calculate Shannon entropy? extern bool event_based; //!< use event-based mode (instead of history-based) extern bool ifp_delayed_group_on; //!< Store delayed group IFP data? @@ -82,6 +83,7 @@ extern bool surf_source_write; //!< write surface source file? extern bool surf_mcpl_write; //!< write surface mcpl file? extern bool survival_biasing; //!< use survival biasing? extern bool survival_normalization; //!< use survival normalization? +extern bool temperature_field_on; //!< Is there a temperature field defined? extern bool temperature_multipole; //!< use multipole data? extern bool trigger_on; //!< tally triggers enabled? extern bool trigger_predict; //!< predict batches for triggers? @@ -187,6 +189,13 @@ extern 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? + //! Whether Iterated Fission Probability is in use at all. inline bool ifp_on() { diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index b5c26d807d6..a138f5b9b05 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 @@ -44,6 +46,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 std::unique_ptr 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..e094594b0e0 --- /dev/null +++ b/include/openmc/streamline_integrator.h @@ -0,0 +1,49 @@ +#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: + virtual ~StreamlineIntegrator() = default; + + //! 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: + RK4StreamlineIntegrator(double time_step) { dt() = time_step; } + + ~RK4StreamlineIntegrator() override = default; + + 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/openmc/__init__.py b/openmc/__init__.py index c204929c840..4f537b83756 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -19,6 +19,8 @@ from openmc.source import * from openmc.settings import * from openmc.lattice import * +from openmc.field import * +from openmc.dnp_drift import * from openmc.filter import * from openmc.filter_expansion import * from openmc.trigger import * diff --git a/openmc/dnp_drift.py b/openmc/dnp_drift.py new file mode 100644 index 00000000000..a8f9891a39d --- /dev/null +++ b/openmc/dnp_drift.py @@ -0,0 +1,345 @@ +import lxml.etree as ET + +import openmc +import openmc.checkvalue as cv +from .field import VelocityField + + +class DNPDrift: + """Settings for delayed neutron precursor drift. + + This class defines the parameters needed to model the transport of + delayed neutron precursors in a flowing fuel system (e.g., molten + salt reactor). + + Parameters + ---------- + velocity_field : openmc.VelocityField + Velocity field describing the fluid motion. + boundary_map : dict + Mapping of boundary types to physical groups. + Keys: "inlet", "outlet", "wall" (optional). Values: list of int. + "inlet" and "outlet" are required. + physical_group_map : dict + Mapping from mesh face IDs to physical group IDs. Must contain + keys "face_ids" and "physical_groups", each mapping to a list + of int of equal length. + integrator : {'RK4'} + Time integration scheme for precursor transport. + integrator_dt : float + Time step for integration in seconds. + recycling : bool + Whether precursors are recycled from outlet to inlet. + external_travel_time : float + Time for precursors to travel through the external loop in + seconds. Used when recycling is True. Default is 1.0s. + + Attributes + ---------- + velocity_field : openmc.VelocityField + Velocity field describing the fluid motion. + boundary_map : dict + Mapping of boundary types to physical groups. + physical_group_map : dict + Mapping from mesh face IDs to physical group IDs. + integrator : {'RK4'} + Time integration scheme. + integrator_dt : float + Time step for integration in seconds. + recycling : bool + Whether precursors are recycled. + external_travel_time : float + External travel time in seconds. + + """ + def __init__( + self, + velocity_field, + boundary_map, + physical_group_map, + integrator="RK4", + integrator_dt=0.1, + recycling=True, + external_travel_time=1.0, + ): + self.velocity_field = velocity_field + self.boundary_map = boundary_map + self.physical_group_map = physical_group_map + self.integrator = integrator + self.integrator_dt = integrator_dt + self.recycling = recycling + self.external_travel_time = external_travel_time + + @property + def velocity_field(self) -> VelocityField: + return self._velocity_field + + @velocity_field.setter + def velocity_field(self, velocity_field: VelocityField): + cv.check_type("velocity_field", velocity_field, VelocityField) + self._velocity_field = velocity_field + + @property + def boundary_map(self) -> dict[str, list[int]]: + return self._boundary_map + + @boundary_map.setter + def boundary_map(self, boundary_map: dict[str, int | list[int]]): + cv.check_type("boundary_map", boundary_map, dict) + + # Create a copy to avoid mutating the original object + bm = {k: v if isinstance(v, list) else [v] for k, v in boundary_map.items()} + + # Check that inlet and outlet are defined + required_keys = {"inlet", "outlet"} + missing_keys = required_keys - bm.keys() + if missing_keys: + raise ValueError( + f"boundary_map is missing required keys: {missing_keys}" + ) + + # Check entries + valid_keys = {"inlet", "outlet", "wall"} + for key in bm: + if key not in valid_keys: + raise ValueError( + f"Invalid boundary type '{key}'. " + f"Must be one of {valid_keys}." + ) + cv.check_type(f"boundary_map[{key}]", bm[key], list) + for sid in bm[key]: + cv.check_type(f"surface id in boundary_map[{key}]", sid, int) + self._boundary_map = bm + + @property + def physical_group_map(self) -> dict[str, list[int]]: + return self._physical_group_map + + @physical_group_map.setter + def physical_group_map(self, physical_group_map: dict[str, list[int]]): + cv.check_type("physical_group_map", physical_group_map, dict) + + if "face_ids" not in physical_group_map: + raise ValueError("physical_group_map must contain 'face_ids'.") + if "physical_groups" not in physical_group_map: + raise ValueError("physical_group_map must contain 'physical_groups'.") + + face_ids = list(physical_group_map["face_ids"]) + physical_groups = list(physical_group_map["physical_groups"]) + + for fid in face_ids: + cv.check_type("face if in physical_group_map", fid, int) + for pg in physical_groups: + cv.check_type("physical group in physical_group_map", pg, int) + + if len(face_ids) != len(physical_groups): + raise ValueError( + f"physical_group_map: face_ids has {len(face_ids)} entries but " + f"physical_groups has {len(physical_groups)} entries. " + f"They must be identical." + ) + + self._physical_group_map = { + "face_ids": face_ids, + "physical_groups": physical_groups, + } + + @property + def integrator(self) -> str: + return self._integrator + + @integrator.setter + def integrator(self, integrator: str): + cv.check_value("integrator", integrator, ("RK4",)) + self._integrator = integrator + + @property + def integrator_dt(self) -> float: + return self._integrator_dt + + @integrator_dt.setter + def integrator_dt(self, integrator_dt: int | float): + cv.check_type("integrator_dt", integrator_dt, (int, float)) + cv.check_greater_than("integrator_dt", integrator_dt, 0) + self._integrator_dt = float(integrator_dt) + + @property + def recycling(self) -> bool: + return self._recycling + + @recycling.setter + def recycling(self, recycling: bool): + cv.check_type("recycling", recycling, bool) + self._recycling = recycling + + @property + def external_travel_time(self) -> float: + return self._external_travel_time + + @external_travel_time.setter + def external_travel_time(self, external_travel_time: int | float): + cv.check_type("external_travel_time", external_travel_time, (int, float)) + cv.check_greater_than("external_travel_time", external_travel_time, 0, True) + self._external_travel_time = float(external_travel_time) + + def to_xml_element(self): + """Create an XML element for DNP drift settings. + + The velocity field data is written separately as a element. + This element only references it by ID. + + Returns + ------- + element : xml.etree._Element + XML element representing DNP drift settings. + + """ + element = ET.Element("dnp_drift") + + vel_elem = ET.SubElement(element, "velocity_field") + vel_elem.text = str(self.velocity_field.id) + + bm_elem = ET.SubElement(element, "boundary_map") + for btype, surface_ids in self.boundary_map.items(): + bt_elem = ET.SubElement(bm_elem, btype) + bt_elem.text = ' '.join(str(s) for s in surface_ids) + + pgm_elem = ET.SubElement(element, "physical_group_map") + face_ids_elem = ET.SubElement(pgm_elem, "face_ids") + face_ids_elem.text = " ".join( + str(i) for i in self.physical_group_map["face_ids"] + ) + pg_elem = ET.SubElement(pgm_elem, "physical_groups") + pg_elem.text = " ".join( + str(i) for i in self.physical_group_map["physical_groups"] + ) + + integrator_elem = ET.SubElement(element, "integrator") + integrator_elem.text = self.integrator + + dt_elem = ET.SubElement(element, "integrator_dt") + dt_elem.text = str(self.integrator_dt) + + recycling_elem = ET.SubElement(element, "recycling") + recycling_elem.text = str(self.recycling).lower() + + ett_elem = ET.SubElement(element, "external_travel_time") + ett_elem.text = str(self.external_travel_time) + + return element + + @classmethod + def from_xml_element(cls, elem, fields): + """Construct a DNPDrift instance from an XML element. + + Parameters + ---------- + elem : xml.etree._Element + XML element to parse. + fields : dict + Dictionary mapping field IDs (int) to Field instances. + + Returns + ------- + DNPDrift + Reconstructed DNP drift settings. + + Raises + ------ + ValueError + If the referenced velocity field ID is not found or is + not a VelocityField. Also, if the boundary_map is not found. + + """ + # Velocity field (required) + vel_elem = elem.find("velocity_field") + if vel_elem is None: + raise ValueError("Missing required in .") + vel_id = int(vel_elem.text) + if vel_id not in fields: + raise ValueError( + f"DNP drift references velocity field id={vel_id}, " + f"but no element with that id was found." + ) + vel_field = fields[vel_id] + if not isinstance(vel_field, VelocityField): + raise ValueError( + f"Field id={vel_id} referenced by is of " + f"type '{vel_field._field_type}', expected 'velocity'." + ) + kwargs = {"velocity_field": vel_field} + + # Boundary map (required) + bm_elem = elem.find("boundary_map") + if bm_elem is None: + raise ValueError("Missing required in .") + boundary_map = {} + for child in bm_elem: + ids = [int(x) for x in child.text.split()] + boundary_map[child.tag] = ids + kwargs["boundary_map"] = boundary_map + + # Physical group map (required) + pgm_elem = elem.find("physical_group_map") + if pgm_elem is None: + raise ValueError("dnp_drift: missing element.") + face_ids_elem = pgm_elem.find("face_ids") + if face_ids_elem is None or face_ids_elem.text is None: + raise ValueError( + "dnp_drift: missing or empty in the element." + ) + pg_elem = pgm_elem.find("physical_groups") + if pg_elem is None or pg_elem.text is None: + raise ValueError( + "dnp_drift: missing or empty in the element." + ) + kwargs["physical_group_map"] = { + "face_ids": [int(i) for i in face_ids_elem.text.split()], + "physical_groups": [int(i) for i in pg_elem.text.split()] + } + + # Integrator (optional) + integrator = elem.find("integrator") + if integrator is not None: + kwargs["integrator"] = integrator.text + + # Integrator time step (optional) + dt = elem.find("integrator_dt") + if dt is not None: + kwargs["integrator_dt"] = float(dt.text) + + # Recycling (optional) + recycling_elem = elem.find("recycling") + if recycling_elem is not None: + kwargs["recycling"] = recycling_elem.text.lower() == "true" + + # External travel time (optional) + ett_elem = elem.find("external_travel_time") + if ett_elem is not None: + kwargs["external_travel_time"] = float(ett_elem.text) + + return cls(**kwargs) + + def __eq__(self, other): + if type(self) is not type(other): + return NotImplemented + return ( + self.velocity_field == other.velocity_field + and self.boundary_map == other.boundary_map + and self.physical_group_map == other.physical_group_map + and self.integrator == other.integrator + and self.integrator_dt == other.integrator_dt + and self.recycling == other.recycling + and self.external_travel_time == other.external_travel_time + ) + + def __repr__(self): + return ( + f"DNPDrift(velocity_field_id={self.velocity_field.id}, " + f"boundary_map={self.boundary_map}, " + f"physical_group_map=({len(self._physical_group_map['face_ids'])} faces), " + f"integrator='{self.integrator}', " + f"integrator_dt={self.integrator_dt}, " + f"recycling={self.recycling}, " + f"external_travel_time={self.external_travel_time})" + ) diff --git a/openmc/field.py b/openmc/field.py new file mode 100644 index 00000000000..6b4f6c307b9 --- /dev/null +++ b/openmc/field.py @@ -0,0 +1,368 @@ +from abc import ABC +import textwrap + +import lxml.etree as ET +import numpy as np +from numpy.typing import ArrayLike + +import openmc +import openmc.checkvalue as cv +from .mixin import IDManagerMixin +from .mesh import RegularMesh, RectilinearMesh + + +# Registry mapping field type strings (e.g., "temperature", "velocity") +# to their corresponding FieldBase subclasses. Populated automatically +# by __init_subclass__. +_FIELD_TYPES = {} + + +class FieldBase(IDManagerMixin, ABC): + """Abstract base class for fields defined on a geometric mesh. + + While there is no abstract method currently, deriving from ABC enforces + the contract of not instantiating FieldBase objects directly. + + Parameters + ---------- + mesh : openmc.RegularMesh or openmc.RectilinearMesh + Spatial mesh on which the field is defined. + values : array-like + Field values. Shape must be (n_elements,) for scalar fields + or (n_elements, n_components) for vector fields. + mapping : {'nodal', 'cell'} + How field values map to the mesh. Default is 'cell'. + field_id : int + Unique identifier for the field. If None, an ID is automatically assigned. + name : str + User-defined name for the field. Default is ''. + + Attributes + ---------- + mesh : openmc.MeshBase + Spatial mesh associated with the field. + values : numpy.ndarray + Array of field values. + mapping : {'nodal', 'cell'} + How field values map to the mesh. + id : int + Unique identifier for the field. + name : str + User-defined name for the field. + next_id : int + Next auto-assigned field ID. + used_ids : set of int + Set of all field IDs currently in use. + _field_type : str or None + String identifier for the field type. Subclasses must set this + to a unique string to be auto-registered. + _n_components : int or None + Number of components per element (1 for scalar field, 3 for vector field). + Subclasses must define this. + _supported_meshes : tuple of type + Mesh types accepted by this field class. + + """ + # Shared across all FieldBase subclasses so that IDs are globally unique + next_id = 1 + used_ids = set() + + _field_type = None + _n_components = None + _supported_meshes = (RegularMesh, RectilinearMesh) + + def __init__(self, mesh, values, mapping="cell", field_id=None, name=''): + self.mesh = mesh + self.mapping = mapping + self.values = values + self.id = field_id + self.name = name + + def __init_subclass__(cls, **kwargs): + """Auto-register subclasses that define _field_type to _FIELD_TYPES.""" + super().__init_subclass__(**kwargs) + if cls._field_type is not None: + _FIELD_TYPES[cls._field_type] = cls + + @property + def mesh(self) -> RegularMesh | RectilinearMesh: + return self._mesh + + @mesh.setter + def mesh(self, mesh: RegularMesh | RectilinearMesh): + cv.check_type('mesh', mesh, self._supported_meshes) + self._mesh = mesh + + @property + def values(self) -> np.ndarray: + return self._values + + @values.setter + def values(self, values: ArrayLike): + if not hasattr(self, "_mapping"): + raise AttributeError("Mapping must be set before values.") + + # Check that _n_components is defined + if self._n_components is None: + raise NotImplementedError( + f"{type(self).__name__} must define _n_components." + ) + + # Check that mesh is set + if not hasattr(self, "_mesh"): + raise AttributeError("Mesh must be set before values.") + + values = np.asarray(values, dtype=float) + + # Determine expected range shape + if self.mapping == "nodal": + n = self._mesh.n_vertices + elif self.mapping == "cell": + n = self._mesh.n_elements + + if self._n_components == 1: + expected = (n,) + else: + expected = (n, self._n_components) + + if values.shape != expected: + raise ValueError( + f"{type(self).__name__} values has shape {values.shape}, " + f"expected {expected} based on mesh with {n} element(s) " + f"and {self._n_components} component(s)." + ) + + self._validate_data(values) + self._values = values + + @property + def mapping(self) -> str: + return self._mapping + + @mapping.setter + def mapping(self, mapping: str): + cv.check_value("mapping", mapping, ("nodal", "cell")) + self._mapping = mapping + + def _validate_data(self, values): + """Subclass-specific validation. + + Parameters + ---------- + values : numpy.ndarray + Data array to validate. + + """ + pass + + def to_xml_element(self): + """Create an XML element representing this field. + + Returns + ------- + element : lxml.etree._Element + XML element for this field. + + """ + element = ET.Element("field") + + # Field type + element.set("type", self._field_type) + + # ID + element.set("id", str(self.id)) + + # Name + if self.name: + element.set("name", self.name) + + # Mesh + mesh_elem = ET.SubElement(element, "mesh") + mesh_elem.text = str(self.mesh.id) + + # Mapping + mapping_elem = ET.SubElement(element, "mapping") + mapping_elem.text = self.mapping + + # Values + values_elem = ET.SubElement(element, "values") + flat = self.values.flatten(order="C") + + # Assumes four XML block levels for the indent + indent = 8 * ' ' + # f"{v:.15e}" instead of str(v) for consistency? + values_elem.text = f"\n{indent}".join( + textwrap.wrap(" ".join(str(v) for v in flat), 80) + ) + + return element + + @classmethod + def from_xml_element(cls, elem, meshes): + """Construct a field from an XML element. + + Parameters + ---------- + elem : lxml.etree._Element + XML element to parse. + meshes : dict + Dictionary mapping mesh IDs (int) to mesh instances + + Returns + ------- + FieldBase + Reconstructed field instance + + """ + # Resolve subclass + field_type = elem.get("type") + if cls is FieldBase: + if field_type not in _FIELD_TYPES: + raise TypeError( + f"Unknown field type '{field_type}'. " + f"Known types: {list(_FIELD_TYPES.keys())}" + ) + subcls = _FIELD_TYPES[field_type] + else: + subcls = cls + + # Mesh + mesh_elem = elem.find("mesh") + if mesh_elem is None: + raise ValueError(f"Field: missing required sub-element.") + mesh_id = int(mesh_elem.text) + if mesh_id not in meshes: + raise ValueError(f"Mesh with id={mesh_id} not found.") + mesh = meshes[mesh_id] + + # Mapping + mapping_elem = elem.find("mapping") + mapping = mapping_elem.text if mapping_elem is not None else "cell" + + # Values + values_elem = elem.find("values") + if values_elem is None: + raise ValueError(f"Field: missing required sub-element.") + values_text = values_elem.text + flat = np.array([float(v) for v in values_text.split()]) + + if mapping == "nodal": + n = mesh.n_vertices + elif mapping == "cell": + n = mesh.n_elements + else: + raise ValueError(f"Unknown mapping '{mapping}'") + + expected_size = n * subcls._n_components + if flat.size != expected_size: + raise ValueError(f"Expected {expected_size} values, got {flat.size}.") + + if subcls._n_components > 1: + flat = flat.reshape((n, subcls._n_components), order="C") + + # Field ID + field_id = int(elem.get("id")) + + # Name + name = elem.get("name", "") + + kwargs = { + "mesh": mesh, + "mapping": mapping, + "values": flat, + "field_id": field_id, + "name": name + } + return subcls(**kwargs) + + def __eq__(self, other): + if type(self) is not type(other): + return NotImplemented + # Only checks that mesh IDs are consistent + return ( + self.mesh.id == other.mesh.id + and np.array_equal(self.values, other.values) + and self.mapping == other.mapping + ) + + def __repr__(self): + repr_str = f"{type(self).__name__}(id={self.id}" + if self.name: + repr_str += f", name='{self.name}'" + repr_str += f", mapping='{self.mapping}'" + repr_str += f", mesh_id={self.mesh.id}" + repr_str += f", n_elements={self.mesh.n_elements})" + return repr_str + + +class TemperatureField(FieldBase): + """Temperature field defined on a mesh. + + Parameters + ---------- + mesh : openmc.RegularMesh or openmc.RectilinearMesh + Spatial mesh on which the field is defined. + values : array-like + Temperature values in Kelvin. Shape must be (n,) where n is the number of + mesh elements or nodes. All values must be strictly positive. + mapping : {'nodal', 'cell'} + How field values map to the mesh. Default is 'cell'. + field_id : int + Unique identifier for the field. If None, an ID is automatically assigned. + name : str + User-defined name for the field. Default is ''. + + """ + _field_type = "temperature" + _n_components = 1 + + def _validate_data(self, values): + """Validate that all temperatures are strictly positive. + + Parameters + ---------- + values : numpy.ndarray + Temperature data array in Kelvin. + + Raises + ------ + ValueError + If any temperature value is zero or negative. + + """ + if np.any(values <= 0.0): + raise ValueError( + "All temperature values must be strictly positive (Kelvin). " + f"Found minimum value: {values.min():.4e}" + ) + + +class VelocityField(FieldBase): + """Velocity field defined on a mesh. + + Parameters + ---------- + mesh : openmc.RegularMesh or openmc.RectilinearMesh + Spatial mesh on which the field is defined. + values : array-like + Velocity vectors (vx, vy, vz) in cm/s. Shape must be (n, 3) + where n is the number of mesh elements or nodes. + mapping : {'nodal', 'cell'} + How field values map to the mesh. Default is 'cell'. + field_id : int + Unique identifier for the field. If None, an ID is automatically assigned. + name : str + User-defined name for the field. Default is ''. + + """ + _field_type = "velocity" + _n_components = 3 + + @property + def magnitude(self): + """Velocity magnitude at each mesh element (cm/s).""" + return np.linalg.norm(self.values, axis=1) + + def max_speed(self): + """Return the maximum speed in the field (cm/s).""" + return float(self.magnitude.max()) diff --git a/openmc/mesh.py b/openmc/mesh.py index 4299fc5530d..d92b81cca85 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1052,6 +1052,8 @@ class RegularMesh(StructuredMesh): The number of mesh cells in each direction (x, y, z). n_dimension : int Number of mesh dimensions. + n_vertices : int + Number of unique vertices. axis_labels : tuple of str Names of the mesh axes ('x', 'y', 'z'), truncated to the mesh dimensionality. @@ -1097,6 +1099,13 @@ def n_dimension(self): else: return None + @property + def n_vertices(self): + n_vertices = 1 + for d in self._dimension: + n_vertices *= d + 1 + return n_vertices + @property def axis_labels(self): return ('x', 'y', 'z')[:self.n_dimension] @@ -1582,6 +1591,8 @@ class RectilinearMesh(StructuredMesh): The number of mesh cells in each direction (x, y, z). n_dimension : int Number of mesh dimensions (always 3 for a RectilinearMesh). + n_vertices : int + Number of unique vertices. axis_labels : tuple of str Names of the mesh axes ('x', 'y', 'z'). x_grid : numpy.ndarray @@ -1616,6 +1627,13 @@ def dimension(self): def n_dimension(self): return 3 + @property + def n_vertices(self): + n_vertices = 1 + for d in self._dimension: + n_vertices *= d + 1 + return n_vertices + @property def axis_labels(self): return ('x', 'y', 'z') diff --git a/openmc/settings.py b/openmc/settings.py index 5d71df091c5..6065578f783 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -15,6 +15,8 @@ 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 FieldBase, TemperatureField, VelocityField +from .dnp_drift import DNPDrift from .utility_funcs import input_path, set_xml_input_path from .volume import VolumeCalculation from .weight_windows import WeightWindows, WeightWindowGenerator, WeightWindowsList @@ -356,6 +358,14 @@ 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 temperature field take precedence + over all other temperature definition (cell, material, and global). + dnp_drift : openmc.DNPDrift + Settings for delayed neutron precursor (DNP) drift in flowing fuel systems. + When specified, precursors are transported along the velocity field rather + than being treated as stationary. trace : tuple or list Show detailed information about a single particle, indicated by three integers: the batch number, generation number, and particle number @@ -448,6 +458,12 @@ def __init__(self, **kwargs): # Shannon entropy mesh self._entropy_mesh = None + # Temperature field + self._temperature_field = None + + # DNP drift + self._dnp_drift = None + # Trigger subelement self._trigger_active = None self._trigger_max_batches = None @@ -802,6 +818,24 @@ 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 dnp_drift(self) -> DNPDrift: + return self._dnp_drift + + @dnp_drift.setter + def dnp_drift(self, dnp_drift: DNPDrift): + cv.check_type("DNP drift", dnp_drift, DNPDrift) + self._dnp_drift = dnp_drift + @property def trigger_active(self) -> bool: return self._trigger_active @@ -1796,6 +1830,31 @@ def _create_entropy_mesh_subelement(self, root, mesh_memo=None): if mesh_memo is not None: mesh_memo.add(self.entropy_mesh.id) + def _append_mesh_if_needed(self, root, mesh, mesh_memo): + """Append mesh element to root if not already present.""" + if mesh_memo is not None and mesh.id in mesh_memo: + return + path = f"./mesh[@id='{mesh.id}']" + if root.find(path) is None: + root.append(mesh.to_xml_element()) + if mesh_memo is not None: + mesh_memo.add(mesh.id) + + def _create_temperature_field_subelement(self, root, mesh_memo=None): + if self._temperature_field is not None: + temp = self._temperature_field + root.append(temp.to_xml_element()) + ref_elem = ET.SubElement(root, "temperature_field") + ref_elem.text = str(temp.id) + self._append_mesh_if_needed(root, temp.mesh, mesh_memo) + + def _create_dnp_drift_subelement(self, root, mesh_memo=None): + if self._dnp_drift is not None: + vel = self._dnp_drift.velocity_field + root.append(vel.to_xml_element()) + root.append(self._dnp_drift.to_xml_element()) + self._append_mesh_if_needed(root, vel.mesh, mesh_memo) + def _create_trigger_subelement(self, root): if self._trigger_active is not None: trigger_element = ET.SubElement(root, "trigger") @@ -2324,6 +2383,28 @@ 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, fields): + ref_elem = root.find("temperature_field") + if ref_elem is not None: + field_id = int(ref_elem.text) + if field_id not in fields: + raise ValueError( + f"Temperature field references field id={field_id}, " + f"but no element with that id was found." + ) + field = fields[field_id] + if not isinstance(field, TemperatureField): + raise ValueError( + f"Field id={field_id} referenced by " + f"is of type '{field._field_type}', expected 'temperature'." + ) + self.temperature_field = field + + def _dnp_drift_from_xml_element(self, root, fields): + dnp_elem = root.find("dnp_drift") + if dnp_elem is not None: + self.dnp_drift = DNPDrift.from_xml_element(dnp_elem, fields) + def _trigger_from_xml_element(self, root): elem = root.find('trigger') if elem is not None: @@ -2629,6 +2710,8 @@ 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_dnp_drift_subelement(element, mesh_memo) self._create_trigger_subelement(element) self._create_no_reduce_subelement(element) self._create_verbosity_subelement(element) @@ -2689,6 +2772,23 @@ def export_to_xml(self, path: PathLike = 'settings.xml'): tree = ET.ElementTree(root_element) tree.write(str(p), xml_declaration=True, encoding='utf-8') + def _read_fields(root, meshes): + """Parse all elements from XML. + + Parameters + ---------- + root : xml.etree._Element + Root XML Element. + meshes : dict + Dictionnary mapping mesh IDs to mesh instance. + + """ + fields = {} + for field_elem in root.findall("field"): + field = FieldBase.from_xml_element(field_elem, meshes) + fields[field.id] = field + return fields + @classmethod def from_xml_element(cls, elem, meshes=None): """Generate settings from XML element @@ -2714,6 +2814,10 @@ def from_xml_element(cls, elem, meshes=None): meshes.update(settings_meshes) settings = cls() + + # read all fields + fields = cls._read_fields(elem, meshes) + settings._eigenvalue_from_xml_element(elem) settings._run_mode_from_xml_element(elem) settings._particles_from_xml_element(elem) @@ -2748,6 +2852,8 @@ 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, fields) + settings._dnp_drift_from_xml_element(elem, fields) settings._trigger_from_xml_element(elem) settings._no_reduce_from_xml_element(elem) settings._verbosity_from_xml_element(elem) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 9135c1782ab..828ee48d534 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -35,36 +35,38 @@ void VacuumBC::handle_particle(Particle& p, const Surface& surf) const // ReflectiveBC implementation //============================================================================== -void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const +void ReflectiveBC::transform(Particle& p, const Surface& surf, Position& new_r, + Direction& new_u, int& new_surface) const { - Direction u = surf.reflect(p.r(), p.u(), &p); - - // normalize reflected u to ensure no floating point error leads to - // unnormalized directions - u /= u.norm(); - - // Handle the effects of the surface albedo on the particle's weight. - BoundaryCondition::handle_albedo(p, surf); + new_r = p.r(); + new_u = surf.reflect(p.r(), p.u(), &p); + new_u /= new_u.norm(); + new_surface = p.surface(); +} - p.cross_reflective_bc(surf, u); +void ReflectiveBC::cross_bc(Particle& p, const Surface& surf, + const Position& new_r, const Direction& new_u, int new_surface) const +{ + p.cross_reflective_bc(surf, new_u); } //============================================================================== // WhiteBC implementation //============================================================================== -void WhiteBC::handle_particle(Particle& p, const Surface& surf) const +void WhiteBC::transform(Particle& p, const Surface& surf, Position& new_r, + Direction& new_u, int& new_surface) const { - Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed()); - - // normalize outgoing u to ensure no floating point error leads to - // unnormalized directions - u /= u.norm(); - - // Handle the effects of the surface albedo on the particle's weight. - BoundaryCondition::handle_albedo(p, surf); + new_r = p.r(); + new_u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed()); + new_u /= new_u.norm(); + new_surface = p.surface(); +} - p.cross_reflective_bc(surf, u); +void WhiteBC::cross_bc(Particle& p, const Surface& surf, const Position& new_r, + const Direction& new_u, int new_surface) const +{ + p.cross_reflective_bc(surf, new_u); } //============================================================================== @@ -132,17 +134,18 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) translation_ = u * (d2 - d1); } -void TranslationalPeriodicBC::handle_particle( - Particle& p, const Surface& surf) const +void TranslationalPeriodicBC::transform(Particle& p, const Surface& surf, + Position& new_r, Direction& new_u, int& new_surface) const { - auto new_r = p.r() + translation_; - int new_surface = p.surface() > 0 ? j_surf_ + 1 : -(j_surf_ + 1); - - // Handle the effects of the surface albedo on the particle's weight. - BoundaryCondition::handle_albedo(p, surf); + new_r = p.r() + translation_; + new_u = p.u(); + new_surface = p.surface() > 0 ? j_surf_ + 1 : -(j_surf_ + 1); +} - // Pass the new location and surface to the particle. - p.cross_periodic_bc(surf, new_r, p.u(), new_surface); +void TranslationalPeriodicBC::cross_bc(Particle& p, const Surface& surf, + const Position& new_r, const Direction& new_u, int new_surface) const +{ + p.cross_periodic_bc(surf, new_r, new_u, new_surface); } //============================================================================== @@ -224,10 +227,10 @@ RotationalPeriodicBC::RotationalPeriodicBC( } } -void RotationalPeriodicBC::handle_particle( - Particle& p, const Surface& surf) const +void RotationalPeriodicBC::transform(Particle& p, const Surface& surf, + Position& new_r, Direction& new_u, int& new_surface) const { - int new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1; + new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1; if (flip_sense_) new_surface = -new_surface; @@ -237,12 +240,10 @@ void RotationalPeriodicBC::handle_particle( double cos_theta = std::cos(angle_); double sin_theta = std::sin(angle_); - Position new_r; new_r[zero_axis_idx_] = r[zero_axis_idx_]; new_r[axis_1_idx_] = cos_theta * r[axis_1_idx_] - sin_theta * r[axis_2_idx_]; new_r[axis_2_idx_] = sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_]; - Direction new_u; new_u[zero_axis_idx_] = u[zero_axis_idx_]; new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_]; new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_]; @@ -250,11 +251,11 @@ void RotationalPeriodicBC::handle_particle( // normalize new_u to ensure no floating point error leads to unnormalized // directions new_u /= new_u.norm(); +} - // Handle the effects of the surface albedo on the particle's weight. - BoundaryCondition::handle_albedo(p, surf); - - // Pass the new location, direction, and surface to the particle. +void RotationalPeriodicBC::cross_bc(Particle& p, const Surface& surf, + const Position& new_r, const Direction& new_u, int new_surface) const +{ p.cross_periodic_bc(surf, new_r, new_u, new_surface); } diff --git a/src/dnp_drift.cpp b/src/dnp_drift.cpp new file mode 100644 index 00000000000..ccea6629438 --- /dev/null +++ b/src/dnp_drift.cpp @@ -0,0 +1,270 @@ +#include "openmc/dnp_drift.h" +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/field.h" +#include "openmc/geometry.h" +#include "openmc/particle_data.h" +#include "openmc/position.h" +#include "openmc/random_lcg.h" +#include "openmc/simulation.h" +#include "openmc/surface.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_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) { + if (decay_time > t_n + settings::dnp_drift_external_travel_time) { + action = Actions::REENTER; + } else { + action = Actions::ESCAPE; + } + } else { + action = Actions::ESCAPE; + } + break; + case BCType::INLET: + action = Actions::REENTER; + break; + case BCType::WALL: + action = Actions::DECAY_IN_PLACE; + break; + default: + fatal_error("Not implemented for this boundary type!"); + break; + } + + // Perform action + switch (action) { + case Actions::REENTER: + 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(); + break; + case Actions::ESCAPE: + t_before_decay = decay_time - t_n; + return false; + break; + case Actions::DECAY_IN_PLACE: + t_before_decay = 0.; + return true; + 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; +} + +bool reconcile_precursor_drift(SourceSite& site, int64_t particle_id) +{ + // Set up temporary particle at the precursor's position and direction + Particle p; + int64_t particle_seed = compute_transport_seed(particle_id); + init_particle_seeds(particle_seed, p.seeds()); + p.stream() = STREAM_TRACKING; + p.r() = site.r; + p.u() = site.u; + + // Case 1: the particle is seen inside the geometry. + if (exhaustive_find_cell(p)) { + return true; + } + + // Case 2: the particle is seen outside the geometry. This means that the + // particle is genuinely outside the geometry, on a boundary but heading + // outward, or traveling exactly along a boundary plane. + + // Nudge the particle slightly backward to check if it reenters. + p.r() = p.r() - p.u() * TINY_BIT; + + if (!exhaustive_find_cell(p)) { + // The particle is not recoverable: it is either in an undefined region or + // traveling exactly along a surface plane. + warning("DNP lost: unable to locate cell after backward nudge. If this " + "warning is occurring frequently, this might indicate spatial " + "inconsistency between the mesh used for DNP transport and the " + "OpenMC model. Please check both the mesh and the OpenMC model."); + return false; + } + + // Compute the distance from the nudged position to the nearest boundary. + BoundaryInfo boundary = distance_to_boundary(p); + + // If the boundary is closer than the nudge distance, the particle was + // genuinely outside before the nudge and not on a boundary. + if (boundary.distance() < TINY_BIT - FP_COINCIDENT) { + return false; + } + + // The particle is on a boundary heading outward. + // Reset particle to its original state and apply boundary conditions. + p.r() = site.r; + p.u() = site.u; + + const auto& surf = *model::surfaces[std::abs(boundary.surface()) - 1]; + const auto* bc = surf.bc_.get(); + + if (!bc) { + fatal_error( + "Boundary condition not found after crossing an external boundary!"); + } + + const std::string& bc_type = bc->type(); + + // Transmission + if (bc_type == "transmission") { + fatal_error("Potentially crossing an external transmission boundary!"); + } + + // Vacuum + if (bc_type == "vacuum") { + return false; + } + + // Remaining conditions: reflective, white, or periodic + // Apply transformation and albedo + bc->transform(p, surf, site.r, site.u, site.surf_id); + if (bc->has_albedo()) { + site.wgt *= bc->albedo(); + } + + return true; +} + +} // namespace openmc diff --git a/src/event.cpp b/src/event.cpp index 2d436bb9dcf..5b0a604a925 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -119,10 +119,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..57a0bfb4a85 --- /dev/null +++ b/src/field.cpp @@ -0,0 +1,304 @@ +#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 { + +namespace model { + +vector> fields; +std::unordered_map field_map; + +} // namespace model + +// ----------------------------------------------------------- +// TemperatureField implementation +// ----------------------------------------------------------- + +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 +// ----------------------------------------------------------- + +int VelocityField::get_next_bin(const Position& r0, const Position& r1, + int bin0, BCType& crossed_boundary, Position& intersection) +{ + // TODO - implement early mesh departure detection + + // The algorithm used to determine if r1 is outside the mesh assumes that + // r0 is inside the mesh + if (bin0 < 0 || bin0 >= mesh_ptr()->n_bins()) { + fatal_error("r0 must be inside the mesh!"); + } + + double total_length = 0.0; + vector outward_surface_ids; + vector inward_surface_ids; + vector bins; + vector segment_lengths; + + mesh_ptr()->bins_and_surface_bins_crossed( + r0, r1, outward_surface_ids, inward_surface_ids, bins, segment_lengths); + + if (segment_lengths.size() == 0) { + + // Try to move the DNP back a little in case it is right on a mesh boundary + Direction u = r1 - r0; + Position r0_nudged = r0 - (u / u.norm() * TINY_BIT); + + // Reset vectors + outward_surface_ids.clear(); + inward_surface_ids.clear(); + bins.clear(); + segment_lengths.clear(); + + mesh_ptr()->bins_and_surface_bins_crossed(r0_nudged, r1, + outward_surface_ids, inward_surface_ids, bins, segment_lengths); + + if (segment_lengths.size() == 0) { + fatal_error("Inconsistency in raytrace results."); + } + + // Adjust total length + total_length -= TINY_BIT; + } + + // 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()->sample_on_physical_groups(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!"); + } +} + +void read_fields(const pugi::xml_node& root) +{ + for (auto node = root.child("field"); node; + node = node.next_sibling("field")) { + + // Type + if (!node.attribute("type")) { + throw std::runtime_error( + "A element is missing the 'type' attribute"); + } + std::string type = node.attribute("type").value(); + + // ID + if (!node.attribute("id")) { + throw std::runtime_error(fmt::format( + "A element of type '{}' is missing the 'id' attribute.", type)); + } + int id = std::stoi(node.attribute("id").value()); + + // Check for duplicate + if (model::field_map.find(id) != model::field_map.end()) { + throw std::runtime_error(fmt::format("Duplicate field id={}.", id)); + } + + // Mesh + if (!check_for_node(node, "mesh")) { + throw std::runtime_error( + fmt::format("Field id={} must specify a mesh.", id)); + } + int mesh_id = std::stoi(get_node_value(node, "mesh")); + if (model::mesh_map.find(mesh_id) == model::mesh_map.end()) { + throw std::runtime_error(fmt::format( + "Mesh {} specified in field id={} does not exist.", mesh_id, id)); + } + Mesh* mesh_ptr = model::meshes[model::mesh_map.at(mesh_id)].get(); + + // Values + if (!check_for_node(node, "values")) { + throw std::runtime_error( + fmt::format("Field id={} must specify values.", id)); + } + auto raw_values = get_node_array(node, "values"); + + // Mapping (optional, default "cell") + FieldMapping mapping = FieldMapping::CELL; + if (check_for_node(node, "mapping")) { + std::string m = get_node_value(node, "mapping"); + if (m == "cell") { + mapping = FieldMapping::CELL; + } else if (m == "nodal") { + mapping = FieldMapping::NODAL; + } else { + throw std::runtime_error(fmt::format( + "Field id={}: unknown mapping '{}'. Must be 'cell' or 'nodal'.", id, + m)); + } + } + + // Create a field + std::unique_ptr field; + std::string entity_name = + (mapping == FieldMapping::CELL) ? "element" : "node"; + size_t n_entities = (mapping == FieldMapping::CELL) + ? mesh_ptr->n_bins() + : mesh_ptr->n_vertices(); + + // Temperature field + if (type == "temperature") { + if (raw_values.size() != n_entities) { + throw std::runtime_error(fmt::format( + "Field id={}: temperature field has {} values but mesh has {} {}s.", + id, raw_values.size(), n_entities, entity_name)); + } + field = std::make_unique(mesh_ptr, raw_values, mapping); + + // Velocity field + } else if (type == "velocity") { + if (raw_values.size() != 3 * n_entities) { + throw std::runtime_error( + fmt::format("Field id={}: velocity field has {} values but expected " + "{} (3 x {} {}s).", + id, raw_values.size(), 3 * n_entities, n_entities, entity_name)); + } + vector vel_values; + vel_values.reserve(n_entities); + for (size_t i = 0; i < raw_values.size(); i += 3) { + vel_values.emplace_back( + raw_values[i], raw_values[i + 1], raw_values[i + 2]); + } + field = std::make_unique(mesh_ptr, vel_values, mapping); + + } else { + throw std::runtime_error( + fmt::format("Unknown field type '{}' for field id={}.", type, id)); + } + + // Store + model::field_map[id] = model::fields.size(); + model::fields.push_back(std::move(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; +} + +extern "C" size_t openmc_temperature_field_size() +{ + return simulation::temperature_field->values().size(); +} + +extern "C" int openmc_temperature_field_get_value( + 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; + } + + *temperature = simulation::temperature_field->value(index); + return 0; +} + +extern "C" int openmc_velocity_field_set_velocity( + int32_t index, double vx, double vy, double vz) +{ + if (index < 0 || index >= simulation::velocity_field->values().size()) { + set_errmsg("Index in velocity field is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + simulation::velocity_field->assign(index, openmc::Position(vx, vy, vz)); + return 0; +} + +extern "C" size_t openmc_velocity_field_size() +{ + return simulation::velocity_field->values().size(); +} + +extern "C" int openmc_velocity_field_get_value( + int32_t index, double velocity[3]) +{ + if (index < 0 || index >= simulation::velocity_field->values().size()) { + set_errmsg("Index in velocity field is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + Direction v = simulation::velocity_field->value(index); + velocity[0] = v.x; + velocity[1] = v.y; + velocity[2] = v.z; + + return 0; +} + +} // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 515f394c4ee..47600fecf01 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -151,6 +151,7 @@ int openmc_finalize() settings::survival_biasing = false; settings::survival_normalization = false; settings::temperature_default = 293.6; + settings::temperature_field_on = false; settings::temperature_method = TemperatureMethod::NEAREST; settings::temperature_multipole = false; settings::temperature_range = {0.0, 0.0}; @@ -185,6 +186,10 @@ int openmc_finalize() simulation::entropy_mesh = nullptr; simulation::ufs_mesh = nullptr; + simulation::temperature_field = nullptr; + simulation::velocity_field = nullptr; + simulation::streamline_integrator.reset(); + data::energy_max = {INFTY, INFTY, INFTY, INFTY}; data::energy_min = {0.0, 0.0, 0.0, 0.0}; data::temperature_min = INFTY; diff --git a/src/geometry.cpp b/src/geometry.cpp index ecacf0bffbf..6f86c78bd89 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -197,11 +197,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..9d710df882b 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 470044e4129..1726b47c1d9 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -858,6 +858,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 //============================================================================== @@ -1055,6 +1084,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; @@ -1488,6 +1524,149 @@ void StructuredMesh::surface_bins_crossed( raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins)); } +void StructuredMesh::bins_and_surface_bins_crossed(Position r0, Position r1, + vector& leaving_surface_ids, vector& entering_surface_ids, + vector& bins, vector& length_fractions) const +{ + // TODO: can reconstruct bins and surface_ids instead of explicitly tracking + // them + + // Helper tally class. + struct CompleteAggregator { + CompleteAggregator(const StructuredMesh* _mesh, + vector& _leaving_surface_ids, vector& _entering_surface_ids, + vector& _bins, vector& _length_fractions) + : mesh(_mesh), leaving_surface_ids(_leaving_surface_ids), + entering_surface_ids(_entering_surface_ids), bins(_bins), + length_fractions(_length_fractions) + {} + // 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) { + entering_surface_ids.push_back(surface_id); + } else { + leaving_surface_ids.push_back(surface_id); + } + } + void track(const MeshIndex& ijk, double l) const + { + bins.push_back(mesh->get_bin_from_indices(ijk)); + length_fractions.push_back(l); + } + + const StructuredMesh* mesh; + vector& leaving_surface_ids; + vector& entering_surface_ids; + vector& bins; + vector& length_fractions; + }; + + // Determine direction + Direction u = (r1 - r0) / (r1 - r0).norm(); + + // Perform the mesh raytrace with the helper class. + raytrace_mesh(r0, r1, u, + CompleteAggregator( + this, leaving_surface_ids, entering_surface_ids, bins, length_fractions)); +} + +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_bin_clamped(Position r0, Position r1, int bin0) const +{ + vector bins; + vector lengths; + bins_crossed(r0, r1, (r1 - r0) / (r1 - r0).norm(), bins, lengths); + if (!bins.empty()) { + return bins.back(); + } else { + return bin0; + } +} + //============================================================================== // RegularMesh implementation //============================================================================== @@ -1637,6 +1816,15 @@ 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 { if (r <= lower_left_[i]) @@ -1685,6 +1873,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 { @@ -1787,6 +2003,209 @@ double RegularMesh::volume(const MeshIndex& ijk) const return element_volume_; } +double RegularMesh::face_area(int face_id) +{ + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + if ((face_id < 0) || (face_id >= (n_bins() * 2 * n_dimension_))) { + fatal_error("Face ID out of bounds!"); + } + + 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::sample_on_physical_groups( + Position& p, int& bin, 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 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()); + } + + // Remove any duplicates to avoid sampling bias + std::sort(face_ids.begin(), face_ids.end()); + face_ids.erase(std::unique(face_ids.begin(), face_ids.end()), face_ids.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 + p = sample_on_face(selected_face_id, seed); + + // Returns bin as well + bin = selected_face_id / (2 * n_dimension_); +} + +Position RegularMesh::sample_on_face(int face_id, uint64_t* seed) +{ + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + if ((face_id < 0) || (face_id >= (n_bins() * 2 * n_dimension_))) { + fatal_error("Face ID out of bounds!"); + } + + 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 +{ + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + for (int k = 0; k < n_dimension_; ++k) { + if ((ijk[k] < 1) || (ijk[k] > shape_[k])) { + fatal_error("Indices out of bounds!"); + } + } + + 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 bounds!"); + } + + return ijk[0] + ijk[1] * (shape_[0] + 1) + + ijk[2] * (shape_[0] + 1) * (shape_[1] + 1); +} + +std::vector RegularMesh::connectivity(int id) const +{ + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + if ((id < 0) || (id >= n_bins())) { + fatal_error("Element ID out of bounds!"); + } + + 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_coordinates(const Position& r, int bin) +{ + // Only implemented in 3D + if (n_dimension_ != 3) { + fatal_error("Only implemented in 3D!"); + } + + 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 //============================================================================== @@ -1856,6 +2275,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/particle.cpp b/src/particle.cpp index 998b71883a8..4a9ba985cb8 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -160,6 +160,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 @@ -223,6 +227,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_bin(r()); + } + if (!exhaustive_find_cell(*this)) { mark_as_lost( "Could not find the cell containing particle " + std::to_string(id())); @@ -276,7 +286,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 @@ -289,14 +299,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); @@ -324,11 +362,6 @@ void Particle::event_advance() score_track_derivative(*this, distance); } - // Set particle weight to zero if it hit the time boundary - if (distance == distance_cutoff) { - wgt() = 0.0; - } - // Clear surface component if distance is long enough if (distance > TINY_BIT) surface() = SURFACE_NONE; @@ -336,63 +369,92 @@ void Particle::event_advance() 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 - - int i_lattice = coord(boundary().coord_level() - 1).lattice(); - 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[i_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(); + } + } + + int i_lattice = coord(boundary().coord_level() - 1).lattice(); + 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[i_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 } } @@ -520,6 +582,12 @@ void Particle::event_revive_from_secondary(const SourceSite& site) // 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_bin(r()); + } + if (!exhaustive_find_cell(*this, verbose)) { mark_as_lost("Could not find the cell containing particle " + std::to_string(id())); @@ -672,6 +740,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 @@ -694,7 +769,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; } @@ -796,6 +875,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. @@ -846,6 +927,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_bin(r()); + } + // Figure out what cell particle is in now n_coord() = 1; diff --git a/src/physics.cpp b/src/physics.cpp index 71dae94b590..60a72eda164 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" @@ -206,9 +207,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(); @@ -220,8 +221,23 @@ 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(); + + // Reject site if it is longer inside the model + if (!transport_dnp(site, dnp_decay_time, p.current_seed())) + continue; + + // Reject site if it is not usable for the next generation + if (!reconcile_precursor_drift(site, p.id())) + continue; + } + + // Reject site if it exceeds time cutoff double t_cutoff = settings::time_cutoff[site.particle.transport_index()]; if (site.time > t_cutoff) { continue; @@ -261,6 +277,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 2364c090802..ef60cba9a30 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -294,6 +294,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; + #ifdef OPENMC_MPI // Only consumed by calculate_rank_load(), which runs during the load // balancing window; past it these atomics would be pure contention on a diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index db8d2342561..9f937170781 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -276,6 +276,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 openmc_finalize_random_ray() diff --git a/src/settings.cpp b/src/settings.cpp index 5e2da3a52e4..7b73f747d7d 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -20,6 +20,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" @@ -54,6 +55,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_delayed_group_on {false}; @@ -77,6 +79,7 @@ bool surf_source_write {false}; bool surf_mcpl_write {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 //============================================================================== @@ -467,6 +473,9 @@ void read_settings_xml(pugi::xml_node root) // Check for user meshes and allocate read_meshes(root); + // Read all fields + read_fields(root); + // Look for deprecated cross_sections.xml file in settings.xml if (check_for_node(root, "cross_sections")) { warning( @@ -856,6 +865,14 @@ 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; + int field_id = std::stoi(get_node_value(root, "temperature_field")); + simulation::temperature_field = get_field(field_id); + } + // Uniform fission source weighting mesh if (check_for_node(root, "ufs_mesh")) { auto temp = std::stoi(get_node_value(root, "ufs_mesh")); @@ -1207,6 +1224,133 @@ 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"); + + // Velocity field + int field_id = std::stoi(get_node_value(node_dnp_drift, "velocity_field")); + simulation::velocity_field = get_field(field_id); + + // 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."); + } + + // Physical group map + if (check_for_node(node_dnp_drift, "physical_group_map")) { + + auto node_physical_group = node_dnp_drift.child("physical_group_map"); + + // 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."); + } + + // Check for duplicate face IDs + std::set unique_face_ids(face_ids.begin(), face_ids.end()); + if (unique_face_ids.size() != face_ids.size()) { + fatal_error("Duplicate face IDs found in physical groups definition!"); + } + + // 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 + simulation::velocity_field->mesh_ptr()->pg_map() = pg_map; + } + + // 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 = + std::make_unique(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")); + } + } + } + } + // 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 4cef7a8e628..a2aaf858d57 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" @@ -21,6 +23,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" @@ -349,6 +352,10 @@ int64_t work_per_rank; const RegularMesh* entropy_mesh {nullptr}; const RegularMesh* ufs_mesh {nullptr}; +TemperatureField* temperature_field {nullptr}; +VelocityField* velocity_field {nullptr}; +std::unique_ptr streamline_integrator; + vector k_generation; vector work_index; @@ -963,10 +970,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_check_limit_and_revive(); diff --git a/src/streamline_integrator.cpp b/src/streamline_integrator.cpp new file mode 100644 index 00000000000..94d2cea6cc2 --- /dev/null +++ b/src/streamline_integrator.cpp @@ -0,0 +1,21 @@ +#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) +{ + // Intermediate velocity estimates + Direction k1 = field->evaluate_in_mesh(yn, cell_n); + Direction k2 = field->evaluate_clamped(yn, yn + dt() / 2. * k1, cell_n); + Direction k3 = field->evaluate_clamped(yn, yn + dt() / 2. * k2, cell_n); + Direction k4 = field->evaluate_clamped(yn, yn + dt() * k3, cell_n); + + // Step forward + yn += dt() / 6. * (k1 + 2 * k2 + 2 * k3 + k4); + tn += dt(); +} + +} // namespace openmc diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 82bcf319ca0..0f11ff7ab56 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -1,6 +1,8 @@ set(TEST_NAMES test_atomic_mass test_distribution + test_dnp_drift + test_field test_file_utils test_tally test_interpolate @@ -19,6 +21,10 @@ set(TEST_NAMES foreach(test ${TEST_NAMES}) add_executable(${test} ${test}.cpp) target_link_libraries(${test} Catch2::Catch2WithMain libopenmc_test) + target_include_directories(${test} PRIVATE matchers) + 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}) # Catch2 uses exit code 4 when all test cases in a binary are skipped. Treat diff --git a/tests/cpp_unit_tests/matchers/equals_position.hpp b/tests/cpp_unit_tests/matchers/equals_position.hpp new file mode 100644 index 00000000000..61c5c0b0665 --- /dev/null +++ b/tests/cpp_unit_tests/matchers/equals_position.hpp @@ -0,0 +1,38 @@ +#include +#include + +#include + +#include "fmt/format.h" +#include "openmc/position.h" + +using namespace openmc; + +class EqualsPositionMatcher : public Catch::Matchers::MatcherGenericBase { +public: + EqualsPositionMatcher(const Position& target, double margin = 1.0e-10) + : target(target), margin(margin) + {} + + bool match(const Position& other) const + { + return std::abs(other.x - target.x) <= margin && + std::abs(other.y - target.y) <= margin && + std::abs(other.z - target.z) <= margin; + } + + std::string describe() const override + { + return fmt::format("Equals: ({}, {}, {})", target.x, target.y, target.z); + } + +private: + const Position& target; + double margin; +}; + +auto EqualsPosition(const Position& target, double margin) + -> EqualsPositionMatcher +{ + return EqualsPositionMatcher {target, margin}; +} diff --git a/tests/cpp_unit_tests/test_dnp_drift.cpp b/tests/cpp_unit_tests/test_dnp_drift.cpp new file mode 100644 index 00000000000..25c06eb5753 --- /dev/null +++ b/tests/cpp_unit_tests/test_dnp_drift.cpp @@ -0,0 +1,451 @@ +#include +#include +#include + +#include +#include + +#include "openmc/dnp_drift.h" +#include "openmc/mesh.h" +#include "openmc/position.h" +#include "openmc/simulation.h" +#include "openmc/streamline_integrator.h" + +using namespace openmc; + +// ---------------------------------------------------------------------------- +// _adjust_position +// ---------------------------------------------------------------------------- + +TEST_CASE("Test linear position adjustment") +{ + // Typical use case + Position a = Position(-1.0, -1.0, -1.0); + Position b = Position(1.0, 1.0, 1.0); + _adjust_position(a, b, 3.0, 1.0, 2.5); + REQUIRE(a == Position(0.0, 0.0, 0.0)); +} + +// ---------------------------------------------------------------------------- +// _adjust_time +// ---------------------------------------------------------------------------- + +TEST_CASE("Test linear time adjustment") +{ + // Typical use case + double t = 3.0; + Position a = Position(-1.0, -1.0, -1.0); + Position b = Position(1.0, 1.0, 1.0); + Position c = Position(0.0, 0.0, 0.0); + _adjust_time(t, 2.0, a, b, c); + REQUIRE(t == 2.5); +} + +// ---------------------------------------------------------------------------- +// transport_dnp +// ---------------------------------------------------------------------------- + +class DNPTransportModelFixture { +protected: + RegularMesh mesh; + BCMap bc_map; + SourceSite site; + uint64_t seed = 1; + + void set_site(const Position& r, const Direction& u) + { + site.r = r; + site.u = u; + } + + void set_velocity_field( + const Direction& v, FieldMapping type = FieldMapping::CELL) + { + vector velocities(8, v); + simulation::velocity_field = new VelocityField(&mesh, velocities, type); + simulation::velocity_field->bc_map() = bc_map; + } + + void set_integrator(double step_size) + { + simulation::streamline_integrator = + std::make_unique(step_size); + } + +public: + DNPTransportModelFixture() + { + // 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"); + mesh = RegularMesh(root); + + // Add physical groups + PGMap pg_map; + vector face_ids = {0, 24, 12, 36, 7, 31, 19, 43, 15, 39, 21, 45, 2, 26, + 8, 32, 4, 16, 10, 22, 29, 41, 35, 47}; + vector physical_groups = { + 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3}; + for (size_t i = 0; i < face_ids.size(); i++) { + pg_map[physical_groups[i]].push_back(face_ids[i]); + } + mesh.pg_map() = pg_map; + + // Boundary condition map + bc_map[BCType::INLET] = {1}; + bc_map[BCType::OUTLET] = {2}; + bc_map[BCType::WALL] = {3, 4}; + + // Default simulation setup + simulation::streamline_integrator = + std::make_unique(0.8); + set_velocity_field(Direction(1.0, 0.0, 0.0)); + + // Default site + set_site(Position(0.0, 0.0, 0.0), Direction(1.0, 0.0, 0.0)); + } + + ~DNPTransportModelFixture() + { + // Reset all openmc global variables + simulation::streamline_integrator.reset(); + settings::dnp_drift_external_travel_time = 0.0; + settings::dnp_drift_recycling_on = false; + delete simulation::velocity_field; + } +}; + +TEST_CASE_METHOD(DNPTransportModelFixture, + "Test DNP transport - cross outlet with recycling " + "on - decay inside the system") +{ + settings::dnp_drift_external_travel_time = 1.0; + settings::dnp_drift_recycling_on = true; + + bool is_inside = transport_dnp(site, 2.5, &seed); + + REQUIRE(is_inside); + REQUIRE(site.r.x == Catch::Approx(-0.5).margin(1.0E-10)); +} + +TEST_CASE_METHOD(DNPTransportModelFixture, + "Test DNP transport - cross outlet with recycling " + "on - decay outside the system") +{ + settings::dnp_drift_external_travel_time = 100.0; + settings::dnp_drift_recycling_on = true; + + bool is_inside = transport_dnp(site, 2.5, &seed); + + REQUIRE(!is_inside); + REQUIRE(site.r == Position(1.0, 0.0, 0.0)); +} + +TEST_CASE_METHOD(DNPTransportModelFixture, + "Test DNP transport - cross outlet with recycling off") +{ + bool is_inside = transport_dnp(site, 2.5, &seed); + + REQUIRE(!is_inside); + REQUIRE(site.r == Position(1.0, 0.0, 0.0)); +} + +TEST_CASE_METHOD(DNPTransportModelFixture, "Test DNP transport - cross wall") +{ + set_velocity_field(Direction(0.0, 1.0, 0.0)); + + bool is_inside = transport_dnp(site, 10.0, &seed); + + REQUIRE(is_inside); + REQUIRE(site.r == Position(0.0, 1.0, 0.0)); +} + +TEST_CASE_METHOD(DNPTransportModelFixture, "Test DNP transport - cross inlet") +{ + simulation::velocity_field->assign(0, Direction(-1.0, 0.0, 0.0)); + set_site(Position(-0.5, -0.5, -0.5), Direction(1.0, 0.0, 0.0)); + + bool is_inside = transport_dnp(site, 1.0, &seed); + + REQUIRE(is_inside); + REQUIRE(site.r.x == Catch::Approx(-0.5).margin(1.0E-10)); +} + +TEST_CASE_METHOD(DNPTransportModelFixture, + "Test DNP transport - stop exactly on mesh boundary") +{ + settings::dnp_drift_external_travel_time = 1.0; + settings::dnp_drift_recycling_on = true; + set_integrator(1.0); + + bool is_inside = transport_dnp(site, 3.0, &seed); + + REQUIRE(is_inside); + REQUIRE(site.r.x == Catch::Approx(0.0).margin(1.0E-10)); +} + +// ---------------------------------------------------------------------------- +// reconcile_precursor_drift +// ---------------------------------------------------------------------------- + +// TODO: change reflective BC to white BC in surface 3 +const std::string RECONCILE_MODEL_XML = R"xml( + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1 + 100 + + + -5.0 -5.0 -5.0 5.0 5.0 5.0 + + + true + + + + +)xml"; + +TEST_CASE("Test reconcile_precursor_drift") +{ + // RAII guard to ensure cleanup happens + struct Cleanup { + bool active = true; + bool initialized = false; + + Cleanup() = default; + Cleanup(const Cleanup&) = delete; + Cleanup& operator=(const Cleanup&) = delete; + + void release() + { + if (active) { + if (initialized) { + openmc_finalize(); + } + std::filesystem::remove("model.xml"); + active = false; + } + } + + ~Cleanup() { release(); } + } cleanup; + + // Write model file + { + std::ofstream f("model.xml"); + REQUIRE(f.is_open()); + f << RECONCILE_MODEL_XML; + f.flush(); + f.close(); + REQUIRE(f.good()); + } + + // Init OpenMC + const char* argv[] = {"test", nullptr}; + int err = openmc_init(1, const_cast(argv), nullptr); + REQUIRE(err == 0); + cleanup.initialized = true; + + // Initialize particle_ids + int64_t particle_id = 1; + + // Inside cell -> true + { + SourceSite site {}; + site.r = {0.0, 0.0, -2.0}; + site.u = {1.0, 0.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, 0.0, -2.0)); + CHECK(site.u == Direction(1.0, 0.0, 0.0)); + CHECK(site.wgt == 1.0); + } + + // TODO: Cannot test outside model case because fatal_error aborts the process + // rather than throwing an exception + + // Outside model -> should fail (fatal error) + //{ + // SourceSite site {}; + // site.r = {10.0, 0.0, 0.0}; + // site.u = {1.0, 0.0, 0.0}; + // REQUIRE_THROWS(reconcile_precursor_drift(site, particle_id)); + //} + + // ----------------------------------------------------------------------------- + // Particles located at an external boundary. + // ----------------------------------------------------------------------------- + + // Vacuum surface, going outward -> false + { + SourceSite site {}; + site.r = {5.0, 0.0, 0.0}; + site.u = {1.0, 0.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == false); + CHECK(site.r == Position(5.0, 0.0, 0.0)); + CHECK(site.u == Direction(1.0, 0.0, 0.0)); + } + + // Vacuum surface, going inward -> true + { + SourceSite site {}; + site.r = {5.0, 0.0, 0.0}; + site.u = {-1.0, 0.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(5.0, 0.0, 0.0)); + CHECK(site.u == Direction(-1.0, 0.0, 0.0)); + CHECK(site.wgt == 1.0); + } + + // Reflective surface, going outward -> reflect and return true + { + SourceSite site {}; + site.r = {0.0, 5.0, 0.0}; + site.u = {0.0, 1.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, 5.0, 0.0)); + CHECK(site.u == Direction(0.0, -1.0, 0.0)); + CHECK(site.wgt == 0.8); + } + + // Reflective surface, going inward -> no reflection, return true + { + SourceSite site {}; + site.r = {0.0, 5.0, 0.0}; + site.u = {0.0, -1.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, 5.0, 0.0)); + CHECK(site.u == Direction(0.0, -1.0, 0.0)); + CHECK(site.wgt == 1.0); + } + + // Periodic surface, going outward -> true + { + SourceSite site {}; + site.r = {0.0, 0.0, 5.0}; + site.u = {0.0, 0.0, 1.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, 0.0, -5.0)); + CHECK(site.u == Direction(0.0, 0.0, 1.0)); + CHECK(site.wgt == 0.8); + } + + // Periodic surface, going inward -> true + { + SourceSite site {}; + site.r = {0.0, 0.0, 5.0}; + site.u = {0.0, 0.0, -1.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, 0.0, 5.0)); + CHECK(site.u == Direction(0.0, 0.0, -1.0)); + CHECK(site.wgt == 1.0); + } + + // White surface, going outward -> reflect and return true + { + SourceSite site {}; + site.r = {0.0, -5.0, 0.0}; + site.u = {0.0, -1.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, -5.0, 0.0)); + CHECK(site.u.y > 0.0); + CHECK(site.wgt == 0.8); + } + + // White surface, going inward -> no reflection and return true + { + SourceSite site {}; + site.r = {0.0, -5.0, 0.0}; + site.u = {0.0, 1.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, -5.0, 0.0)); + CHECK(site.u == Direction(0.0, 1.0, 0.0)); + CHECK(site.wgt == 1.0); + } + + // ----------------------------------------------------------------------------- + // Particles traveling exactly along external boundary planes. + // ----------------------------------------------------------------------------- + + // Any external surface (reflective in this case). The particle is originally + // seen inside -> true + { + SourceSite site {}; + site.r = {0.0, 5.0, 0.0}; + site.u = {1.0, 0.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, 5.0, 0.0)); + CHECK(site.u == Direction(1.0, 0.0, 0.0)); + CHECK(site.wgt == 1.0); + } + + // Any external surface (white in this case). The particle is originally seen + // outside and the backward nudge cannot recover the particle since it remains + // on the surface -> false + { + SourceSite site {}; + site.r = {0.0, -5.0, 0.0}; + site.u = {1.0, 0.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == false); + CHECK(site.r == Position(0.0, -5.0, 0.0)); + CHECK(site.u == Direction(1.0, 0.0, 0.0)); + CHECK(site.wgt == 1.0); + } + + // Any external surface (periodic in this case). The particle is originally + // seen inside -> true + { + SourceSite site {}; + site.r = {0.0, 0.0, 5.0}; + site.u = {1.0, 0.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == true); + CHECK(site.r == Position(0.0, 0.0, 5.0)); + CHECK(site.u == Direction(1.0, 0.0, 0.0)); + CHECK(site.wgt == 1.0); + } + + // Any external surface (periodic in this case). The particle is originally + // seen outside and the backward nudge cannot recover the particle since it + // remains on the surface -> false + { + SourceSite site {}; + site.r = {0.0, 0.0, -5.0}; + site.u = {1.0, 0.0, 0.0}; + CHECK(reconcile_precursor_drift(site, particle_id) == false); + CHECK(site.r == Position(0.0, 0.0, -5.0)); + CHECK(site.u == Direction(1.0, 0.0, 0.0)); + CHECK(site.wgt == 1.0); + } +} diff --git a/tests/cpp_unit_tests/test_field.cpp b/tests/cpp_unit_tests/test_field.cpp new file mode 100644 index 00000000000..ecceab80304 --- /dev/null +++ b/tests/cpp_unit_tests/test_field.cpp @@ -0,0 +1,338 @@ +#include + +#include +#include +#include +#include +#include + +#include "openmc/field.h" +#include "openmc/geometry.h" +#include "openmc/mesh.h" + +using namespace openmc; + +class RegularMeshFixture { +protected: + RegularMesh mesh; + +public: + RegularMeshFixture() + { + // 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"); + mesh = RegularMesh(root); + } +}; + +TEST_CASE_METHOD(RegularMeshFixture, "Test Field - regular mesh") +{ + auto [mapping, values, outputs] = + GENERATE(table, vector>( + {{FieldMapping::CELL, {10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0}, + {80.0, 10.0, 20.0, 50.0, 10.0, 10.0}}, + {FieldMapping::NODAL, + {10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, + 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, + 32.0, 33.0, 34.0, 35.0, 36.0}, + {29.5, 16.5, 17.5, 25.5, 23.0, 10.0}}})); + + // Create a typed field + TypedField field = TypedField(&mesh, values, mapping); + + // Assign + double saved_value = field.value(2); + field.assign(2, 150.0); + REQUIRE(field.value(2) == 150.0); + field.assign(2, saved_value); // Reassign back + + // Get bin + REQUIRE(field.get_bin(Position(0.5, 0.5, 0.5)) == 7); + REQUIRE(field.get_bin(Position(-0.5, -0.5, -0.5)) == 0); + REQUIRE(field.get_bin(Position(0.5, -0.5, -0.5)) == 1); + REQUIRE(field.get_bin(Position(-0.5, -0.5, 0.5)) == 4); + REQUIRE(field.get_bin(Position(0.0, 0.0, 0.0)) == 0); + REQUIRE(field.get_bin(Position(2.0, 2.0, 2.0)) == -1); + + // Evaluate in mesh (current position, current bin) + REQUIRE(field.evaluate_in_mesh(Position(0.5, 0.5, 0.5), 7) == outputs[0]); + REQUIRE(field.evaluate_in_mesh(Position(-0.5, -0.5, -0.5), 0) == outputs[1]); + REQUIRE(field.evaluate_in_mesh(Position(0.5, -0.5, -0.5), 1) == outputs[2]); + REQUIRE(field.evaluate_in_mesh(Position(-0.5, -0.5, 0.5), 4) == outputs[3]); + REQUIRE(field.evaluate_in_mesh(Position(0.0, 0.0, 0.0), 0) == outputs[4]); + + // Trilinear interpolation + if (mapping == FieldMapping::NODAL) { + REQUIRE( + field.trilinear_interpolation(Position(0.5, 0.5, 0.5), 7) == outputs[0]); + REQUIRE(field.trilinear_interpolation(Position(-0.5, -0.5, -0.5), 0) == + outputs[1]); + REQUIRE(field.trilinear_interpolation(Position(0.5, -0.5, -0.5), 1) == + outputs[2]); + REQUIRE(field.trilinear_interpolation(Position(-0.5, -0.5, 0.5), 4) == + outputs[3]); + REQUIRE( + field.trilinear_interpolation(Position(0.0, 0.0, 0.0), 0) == outputs[4]); + } + + // Evaluate clamped (previous position, current position, previous bin) + REQUIRE(field.evaluate_clamped(Position(0.5, 0.5, 0.5), + Position(-0.5, -0.5, -0.5), 7) == outputs[1]); + REQUIRE(field.evaluate_clamped(Position(-0.5, -0.5, -0.5), + Position(0.5, 0.5, 0.5), 0) == outputs[0]); + REQUIRE(field.evaluate_clamped(Position(0.0, 0.0, 0.0), + Position(-1.0, -1.0, -1.0), 0) == outputs[5]); + + // Distance to next boundary + 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_METHOD( + RegularMeshFixture, "Test TemperatureField - regular mesh - cell-based only") +{ + // 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)); +} + +TEST_CASE_METHOD(RegularMeshFixture, "Test VelocityField - regular mesh") +{ + auto [mapping, values] = GENERATE(table>( + {{FieldMapping::CELL, + {Direction(10.0, 10.0, 10.0), Direction(20.0, 20.0, 20.0), + Direction(30.0, 30.0, 30.0), Direction(40.0, 40.0, 40.0), + Direction(50.0, 50.0, 50.0), Direction(60.0, 60.0, 60.0), + Direction(70.0, 70.0, 70.0), Direction(80.0, 80.0, 80.0)}}, + {FieldMapping::NODAL, + {Direction(10.0, 10.0, 10.0), Direction(11.0, 11.0, 11.0), + Direction(12.0, 12.0, 12.0), Direction(13.0, 13.0, 13.0), + Direction(14.0, 14.0, 14.0), Direction(15.0, 15.0, 15.0), + Direction(16.0, 16.0, 16.0), Direction(17.0, 17.0, 17.0), + Direction(18.0, 18.0, 18.0), Direction(19.0, 19.0, 19.0), + Direction(20.0, 20.0, 20.0), Direction(21.0, 21.0, 21.0), + Direction(22.0, 22.0, 22.0), Direction(23.0, 23.0, 23.0), + Direction(24.0, 24.0, 24.0), Direction(25.0, 25.0, 25.0), + Direction(26.0, 26.0, 26.0), Direction(27.0, 27.0, 27.0), + Direction(28.0, 28.0, 28.0), Direction(29.0, 29.0, 29.0), + Direction(30.0, 30.0, 30.0), Direction(31.0, 31.0, 31.0), + Direction(32.0, 32.0, 32.0), Direction(33.0, 33.0, 33.0), + Direction(34.0, 34.0, 34.0), Direction(35.0, 35.0, 35.0), + Direction(36.0, 36.0, 36.0)}}})); + + // Add physical group map + mesh.pg_map() = {{1, {0, 24, 12, 36}}, {2, {7, 31, 19, 43}}, + {3, {15, 39, 21, 45, 2, 26, 8, 32, 4, 16}}, {4, {10, 22, 29, 41, 35, 47}}}; + + // Create a temperature field + VelocityField velocity_field = VelocityField(&mesh, values, mapping); + + // Add boundary conditions map + velocity_field.bc_map() = { + {BCType::INLET, {1}}, {BCType::OUTLET, {2}}, {BCType::WALL, {3, 4}}}; + + // Get next bin + int bin1; + BCType crossed_boundary; + Position intersection; + + // - Next point is inside the mesh + bin1 = velocity_field.get_next_bin(Position(-0.5, -0.5, -0.5), + Position(0.5, 0.5, 0.5), 0, crossed_boundary, intersection); + CHECK(bin1 == 7); + CHECK(crossed_boundary == BCType::NONE); + CHECK(intersection == Position(0.5, 0.5, 0.5)); + + // - Next point is outside the mesh + bin1 = velocity_field.get_next_bin(Position(-0.5, -0.5, -0.5), + Position(-0.5, -0.5, 1.5), 0, crossed_boundary, intersection); + CHECK(bin1 == -1); + CHECK(crossed_boundary == BCType::WALL); + CHECK(intersection == Position(-0.5, -0.5, 1.0)); + + // - r0 is outside the mesh based on bin0 - should fail + // bin1 = velocity_field.get_next_bin( + // Position(), Position(), -1, crossed_boundary, intersection); + // bin1 = velocity_field.get_next_bin( + // Position(), Position(), 8, crossed_boundary, intersection); + + // Randomly place on inlet + Position p = Position(0., 0., 0.); + int bin = 0; + uint64_t seed = 1; + velocity_field.randomly_place_on_inlet(p, bin, &seed); + CHECK(p.x == -1.0); + CHECK(p.y == Catch::Approx(-0.2470039942).margin(1.0E-10)); + CHECK(p.z == Catch::Approx(0.3301189888).margin(1.0E-10)); + CHECK(bin == 4); + + // Get boundary conditions + CHECK(velocity_field.get_boundary_condition(1) == BCType::INLET); + CHECK(velocity_field.get_boundary_condition(2) == BCType::OUTLET); + CHECK(velocity_field.get_boundary_condition(3) == BCType::WALL); + CHECK(velocity_field.get_boundary_condition(4) == BCType::WALL); +} + +TEST_CASE( + "Test settings declaration exceptions - TemperatureField", "[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 + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 + cell + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + + )", + "Field id=1: temperature field has 7 values but mesh has 8 elements."}, + {// If the mesh declared is not defined -> error + R"( + + eigenvalue + 200 + 20 + 1 + + 2 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + cell + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + + )", + "Mesh 2 specified in field id=1 does not exist."}, + {// No mesh declared -> error + R"( + + eigenvalue + 200 + 20 + 1 + + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + cell + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + + )", + "Field id=1 must specify a mesh."}, + {// No values declared -> error + R"( + + eigenvalue + 200 + 20 + 1 + + 1 + cell + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + + )", + "Field id=1 must specify values."}, + })); + + 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 40fad0ffe76..272e3496471 100644 --- a/tests/cpp_unit_tests/test_mesh.cpp +++ b/tests/cpp_unit_tests/test_mesh.cpp @@ -1,11 +1,15 @@ #include #include #include +#include #include #include #include #include +#include +#include +#include #include #include "openmc/hdf5_interface.h" @@ -16,6 +20,8 @@ #include "xdg/config.h" #endif +#include "matchers/equals_position.hpp" + using namespace openmc; #ifdef OPENMC_XDG_ENABLED @@ -298,81 +304,146 @@ TEST_CASE("Test multiple meshes HDF5 roundtrip - spherical") REQUIRE(regular_mesh_hdf5->upper_right() == regular_mesh_xml->upper_right()); } -TEST_CASE("Test get_index_in_direction - regular") -{ - // The XML data as a string - std::string xml_string = R"( - - 2 1 1 - -2 -2 -2 - 2 2 2 - - )"; - - // Create a pugixml document object - pugi::xml_document doc; - - // Load the XML from the string - pugi::xml_parse_result result = doc.load_string(xml_string.c_str()); - - pugi::xml_node root = doc.child("mesh"); - - auto mesh = RegularMesh(root); - - REQUIRE(mesh.get_index_in_direction(-4.1, 0) == 0); - REQUIRE(mesh.get_index_in_direction(-4.0, 0) == 0); - REQUIRE(mesh.get_index_in_direction(-3.9, 0) == 0); - - REQUIRE(mesh.get_index_in_direction(-2.1, 0) == 0); - REQUIRE(mesh.get_index_in_direction(-2.0, 0) == 1); // lower left - REQUIRE(mesh.get_index_in_direction(-1.9, 0) == 1); - - REQUIRE(mesh.get_index_in_direction(-0.1, 0) == 1); - REQUIRE(mesh.get_index_in_direction(0.0, 0) == 1); - REQUIRE(mesh.get_index_in_direction(0.1, 0) == 2); - - REQUIRE(mesh.get_index_in_direction(1.9, 0) == 2); - REQUIRE(mesh.get_index_in_direction(2.0, 0) == 2); // upper right - REQUIRE(mesh.get_index_in_direction(2.1, 0) == 3); - - REQUIRE(mesh.get_index_in_direction(3.9, 0) == 3); - REQUIRE(mesh.get_index_in_direction(4.0, 0) == 3); - REQUIRE(mesh.get_index_in_direction(4.1, 0) == 3); -} - -TEST_CASE("Test get_index_in_direction - rectilinear") -{ - // The XML data as a string - std::string xml_string = R"( - - -1.0 0.0 2.0 - -1.0 1.0 - -1.0 1.0 +class RegularMeshFixture { +protected: + RegularMesh mesh; + +public: + RegularMeshFixture() + { + std::string xml_string = R"( + + 2 2 2 + -1 -1 -1 + 1 1 1 - )"; - - // Create a pugixml document object - pugi::xml_document doc; - - // Load the XML from the string - pugi::xml_parse_result result = doc.load_string(xml_string.c_str()); - - pugi::xml_node root = doc.child("mesh"); + )"; + + pugi::xml_document doc; + pugi::xml_parse_result result = doc.load_string(xml_string.c_str()); + pugi::xml_node root = doc.child("mesh"); + mesh = RegularMesh(root); + + // Add physical groups + PGMap pg_map; + vector face_ids = {0, 24, 12, 36, 7, 31, 19, 43, 15, 39, 21, 45, 2, 26, + 8, 32, 4, 16, 10, 22, 29, 41, 35, 47}; + vector physical_groups = { + 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3}; + for (size_t i = 0; i < face_ids.size(); i++) { + pg_map[physical_groups[i]].push_back(face_ids[i]); + } + mesh.pg_map() = pg_map; + } +}; + +class RectilinearMeshFixture { +protected: + RectilinearMesh mesh; + +public: + RectilinearMeshFixture() + { + std::string xml_string = R"( + + -1.0 0.0 1.0 + -1.0 0.0 1.0 + -1.0 0.0 1.0 + + )"; + + pugi::xml_document doc; + pugi::xml_parse_result result = doc.load_string(xml_string.c_str()); + pugi::xml_node root = doc.child("mesh"); + mesh = RectilinearMesh(root); + } +}; + +#define GET_DISTANCE_TO_NEXT_BOUNDARY_TESTS \ + int current_bin; \ + Position r; \ + Position u; \ + int next_bin; \ + double distance; \ + \ + SECTION("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); \ + } \ + \ + SECTION("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); \ + } \ + \ + SECTION("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); \ + } \ + \ + SECTION("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); \ + } \ + \ + SECTION("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); \ + } - auto mesh = RectilinearMesh(root); +TEST_CASE_METHOD(RegularMeshFixture, "Test distance_to_next_boundary()") { + GET_DISTANCE_TO_NEXT_BOUNDARY_TESTS} + +TEST_CASE_METHOD(RectilinearMeshFixture, "Test distance_to_next_boundary()") { + GET_DISTANCE_TO_NEXT_BOUNDARY_TESTS} + +#undef GET_DISTANCE_TO_NEXT_BOUNDARY_TESTS + +#define GET_GET_INDEX_IN_DIRECTION_TESTS \ + REQUIRE(mesh.get_index_in_direction(-5.0, 0) == 0); \ + \ + REQUIRE(mesh.get_index_in_direction(-1.1, 0) == 0); \ + REQUIRE(mesh.get_index_in_direction(-1.0, 0) == 1); \ + REQUIRE(mesh.get_index_in_direction(-0.9, 0) == 1); \ + \ + REQUIRE(mesh.get_index_in_direction(0.9, 0) == 2); \ + REQUIRE(mesh.get_index_in_direction(1.0, 0) == 2); \ + REQUIRE(mesh.get_index_in_direction(1.1, 0) == 3); \ + \ + REQUIRE(mesh.get_index_in_direction(5.0, 0) == 3); - REQUIRE(mesh.get_index_in_direction(-5.0, 0) == 0); +TEST_CASE_METHOD(RegularMeshFixture, "Test get_index_in_direction()") { + GET_GET_INDEX_IN_DIRECTION_TESTS} - REQUIRE(mesh.get_index_in_direction(-1.1, 0) == 0); - REQUIRE(mesh.get_index_in_direction(-1.0, 0) == 1); // lower left - REQUIRE(mesh.get_index_in_direction(-0.9, 0) == 1); +TEST_CASE_METHOD(RectilinearMeshFixture, "Test get_index_in_direction()") { + GET_GET_INDEX_IN_DIRECTION_TESTS} - REQUIRE(mesh.get_index_in_direction(1.9, 0) == 2); - REQUIRE(mesh.get_index_in_direction(2.0, 0) == 2); // upper right - REQUIRE(mesh.get_index_in_direction(2.1, 0) == 3); - - REQUIRE(mesh.get_index_in_direction(5.0, 0) == 3); -} +#undef GET_GET_INDEX_IN_DIRECTION_TESTS TEST_CASE("Test regular mesh ray tracing from outside") { @@ -410,3 +481,310 @@ TEST_CASE("Test regular mesh ray tracing from outside") REQUIRE(lengths[0] == Catch::Approx(1.0 / 6.0)); REQUIRE(lengths[1] == Catch::Approx(1.0 / 6.0)); } + +TEST_CASE_METHOD(RegularMeshFixture, "Test bins_and_surface_bins_crossed()") +{ + // Test cases: + // - r0 inside and r1 inside + // - r0 inside and r1 outside + // - r0 outside and r1 inside + // - r0 outside and r1 outside + // - r0 inside and r1 inside: corner check + auto [r0, r1, expected_leaving_surface_ids, expected_entering_surface_ids, + expected_bins, expected_length_fractions] = + GENERATE(table, vector, vector, + vector>({{Position(-0.5, -0.5, -0.5), Position(0.5, -0.5, -0.5), + {1}, {6}, {0, 1}, {0.5, 0.5}}, + {Position(-0.5, -0.5, -0.5), Position(2.0, -0.5, -0.5), {1, 7}, {6}, + {0, 1}, {0.2, 0.4}}, + {Position(-2.0, -0.5, -0.5), Position(0.5, -0.5, -0.5), {1}, {0, 6}, + {0, 1}, {0.4, 0.2}}, + {Position(-2.0, -0.5, -0.5), Position(2.0, -0.5, -0.5), {1, 7}, {0, 6}, + {0, 1}, {0.25, 0.25}}, + {Position(-0.5, -0.5, -0.5), Position(0.5, 0.5, 0.5), {1, 9, 23}, + {6, 20, 46}, {0, 1, 3, 7}, {0.5, 0.0, 0.0, 0.5}}})); + + vector leaving_surface_ids; + vector entering_surface_ids; + vector bins; + vector length_fractions; + + mesh.bins_and_surface_bins_crossed( + r0, r1, leaving_surface_ids, entering_surface_ids, bins, length_fractions); + REQUIRE_THAT( + leaving_surface_ids, Catch::Matchers::Equals(expected_leaving_surface_ids)); + REQUIRE_THAT(entering_surface_ids, + Catch::Matchers::Equals(expected_entering_surface_ids)); + REQUIRE_THAT(bins, Catch::Matchers::Equals(expected_bins)); + REQUIRE(length_fractions.size() == expected_length_fractions.size()); + for (size_t i = 0; i < expected_length_fractions.size(); i++) { + REQUIRE_THAT(length_fractions[i], + Catch::Matchers::WithinAbs(expected_length_fractions[i], 1.0e-10)); + } +} + +TEST_CASE_METHOD(RegularMeshFixture, "Test normalize_coordinates()") +{ + // r inside bin + REQUIRE_THAT(mesh.normalize_coordinates(Position(-0.5, -0.5, -0.5), 0), + EqualsPosition(Position(0.5, 0.5, 0.5), 1.0e-10)); + REQUIRE_THAT(mesh.normalize_coordinates(Position(-0.1, -0.7, -0.2), 0), + EqualsPosition(Position(0.9, 0.3, 0.8), 1.0e-10)); + + // r outside bin + REQUIRE_THAT(mesh.normalize_coordinates(Position(0.5, -0.5, -0.5), 0), + EqualsPosition(Position(1.0, 0.5, 0.5), 1.0e-10)); + REQUIRE_THAT(mesh.normalize_coordinates(Position(-0.5, -0.5, -0.5), 1), + EqualsPosition(Position(0.0, 0.5, 0.5), 1.0e-10)); +} + +TEST_CASE_METHOD(RegularMeshFixture, "Test face_area()") +{ + REQUIRE(mesh.face_area(0) == 1.0); + REQUIRE(mesh.face_area(1) == 1.0); + REQUIRE(mesh.face_area(2) == 1.0); + REQUIRE(mesh.face_area(3) == 1.0); + REQUIRE(mesh.face_area(4) == 1.0); + REQUIRE(mesh.face_area(5) == 1.0); + + // mesh.face_area(-1); // should fail + // mesh.face_area(48); // should fail +} + +TEST_CASE_METHOD(RegularMeshFixture, "Test sample_on_face()") +{ + // Check that sampled points are on the expected surface for different seeds + + uint64_t seed = GENERATE(1, 2); + + REQUIRE(mesh.sample_on_face(0, &seed).x == -1.0); + REQUIRE(mesh.sample_on_face(1, &seed).x == 0.0); + REQUIRE(mesh.sample_on_face(2, &seed).y == -1.0); + REQUIRE(mesh.sample_on_face(3, &seed).y == 0.0); + REQUIRE(mesh.sample_on_face(4, &seed).z == -1.0); + REQUIRE(mesh.sample_on_face(5, &seed).z == 0.0); + + REQUIRE(mesh.sample_on_face(42, &seed).x == 0.0); + REQUIRE(mesh.sample_on_face(43, &seed).x == 1.0); + REQUIRE(mesh.sample_on_face(44, &seed).y == 0.0); + REQUIRE(mesh.sample_on_face(45, &seed).y == 1.0); + REQUIRE(mesh.sample_on_face(46, &seed).z == 0.0); + REQUIRE(mesh.sample_on_face(47, &seed).z == 1.0); + + // mesh.sample_on_face(-1, &seed); // should fail + // mesh.sample_on_face(100, &seed); // should fail +} + +TEST_CASE_METHOD(RegularMeshFixture, "Regression test sample_on_face()") +{ + // Check that the other 2 coordinates are sampled + + uint64_t seed = 1; + REQUIRE_THAT(mesh.sample_on_face(0, &seed), + EqualsPosition(Position(-1.0, -0.2891826401, -0.2470039942), 1.0e-10)); +} + +TEST_CASE_METHOD(RegularMeshFixture, "Test return_vertex_unique_id()") +{ + // Lower left + REQUIRE(mesh.return_vertex_unique_id({1, 1, 1}, 0) == 0); + + // Upper right + REQUIRE(mesh.return_vertex_unique_id({2, 2, 2}, 7) == 26); + + // Middle from lower left bin + REQUIRE(mesh.return_vertex_unique_id({1, 1, 1}, 7) == 13); + + // Middle from upper right bin + REQUIRE(mesh.return_vertex_unique_id({2, 2, 2}, 0) == 13); + + // mesh.return_vertex_unique_id({0, 0, 0}, 0); // should fail + // mesh.return_vertex_unique_id({3, 3, 3}, 0); // should fail + // mesh.return_vertex_unique_id({1, 1, 1}, -1); // should fail + // mesh.return_vertex_unique_id({1, 1, 1}, 8); // should fail +} + +TEST_CASE_METHOD(RegularMeshFixture, "Test connectivity()") +{ + // Lower left element + vector expected = {0, 1, 3, 4, 9, 10, 12, 13}; + REQUIRE(mesh.connectivity(0) == expected); + + // Upper right element + expected = {13, 14, 16, 17, 22, 23, 25, 26}; + REQUIRE(mesh.connectivity(7) == expected); + + // mesh.connectivity(-1); // should fail + // mesh.connectivity(8); // should fail +} + +#define GET_BIN_CLAMPED_TESTS \ + SECTION("r0 and r1 inside mesh - returns bin containing r1") \ + { \ + Position r0 {-0.5, -0.5, -0.5}; \ + Position r1 {0.5, 0.5, -0.5}; \ + int bin0 = 0; \ + REQUIRE(mesh.get_bin_clamped(r0, r1, bin0) == 3); \ + } \ + \ + SECTION( \ + "r0 inside mesh and r1 outside mesh - returns last bin before leaving") \ + { \ + Position r0 {-0.5, -0.5, -0.5}; \ + Position r1 {1.5, 1.5, 1.5}; \ + int bin0 = 0; \ + REQUIRE(mesh.get_bin_clamped(r0, r1, bin0) == 7); \ + } \ + \ + SECTION( \ + "r0 inside mesh and r1 on internal boundary - returns bin containing r1") \ + { \ + Position r0 {-0.5, -0.5, -0.5}; \ + Position r1 {1.0, 1.0, 1.0}; \ + int bin0 = 0; \ + REQUIRE(mesh.get_bin_clamped(r0, r1, bin0) == 7); \ + } \ + \ + SECTION("r1 == r0 - returns bin0 (no movement)") \ + { \ + Position r0 {-0.5, -0.5, -0.5}; \ + Position r1 {-0.5, -0.5, -0.5}; \ + int bin0 = 0; \ + REQUIRE(mesh.get_bin_clamped(r0, r1, bin0) == bin0); \ + } \ + \ + SECTION("r0 and r1 outside mesh with no traversal - returns bin0") \ + { \ + Position r0 {-1.5, -1.5, -1.5}; \ + Position r1 {-1.2, -1.2, -1.2}; \ + int bin0 = -1; \ + REQUIRE(mesh.get_bin_clamped(r0, r1, bin0) == bin0); \ + } \ + \ + SECTION("r0 and r1 outside mesh with traversal - last bin before leaving") \ + { \ + Position r0 {-1.5, -1.5, -1.5}; \ + Position r1 {1.5, 1.5, 1.5}; \ + int bin0 = -1; \ + REQUIRE(mesh.get_bin_clamped(r0, r1, bin0) == 7); \ + } + +TEST_CASE_METHOD(RegularMeshFixture, "Test get_bin_clamped()") { + GET_BIN_CLAMPED_TESTS} + +TEST_CASE_METHOD(RectilinearMeshFixture, "Test get_bin_clamped()") { + GET_BIN_CLAMPED_TESTS} + +#undef GET_BIN_CLAMPED_TESTS + +// TODO: test get_bin_clamped with more mesh types + +TEST_CASE_METHOD(RegularMeshFixture, "Test sample_on_physical_groups()") +{ + const double tol = 1.0e-12; + uint64_t seed = 1; + + SECTION("Accept a single physical group") + { + // Using multiple calls of sample_on_physical_groups(), we verify that each + // sample is on x=-1 (group 1) + std::set expected_bins = {0, 2, 4, 6}; + vector physical_groups = {1}; + + for (int i = 0; i < 1000; ++i) { + Position p; + int bin; + mesh.sample_on_physical_groups(p, bin, &seed, physical_groups); + + REQUIRE_THAT(p.x, Catch::Matchers::WithinAbs(-1.0, tol)); + REQUIRE(expected_bins.count(bin) == 1); + REQUIRE(p.y >= -1.0); + REQUIRE(p.y <= 1.0); + REQUIRE(p.z >= -1.0); + REQUIRE(p.z <= 1.0); + } + } + + SECTION("Accept multiple physical groups") + { + // Using multiple calls of sample_on_physical_groups(), we verify that each + // sample is either on x=-1 (group 1) or x=1 (group 2) + vector physical_groups = {1, 2}; + + bool sampled_x_min = false; + bool sampled_x_max = false; + + for (int i = 0; i < 1000; ++i) { + Position p; + int bin; + mesh.sample_on_physical_groups(p, bin, &seed, physical_groups); + + bool on_x_min = std::abs(p.x - (-1.0)) < tol; + bool on_x_max = std::abs(p.x - 1.0) < tol; + REQUIRE((on_x_min || on_x_max)); + + if (on_x_min) + sampled_x_min = true; + if (on_x_max) + sampled_x_max = true; + + REQUIRE(p.y >= -1.0); + REQUIRE(p.y <= 1.0); + REQUIRE(p.z >= -1.0); + REQUIRE(p.z <= 1.0); + } + + // Both faces should have been sampled at least once + REQUIRE(sampled_x_min); + REQUIRE(sampled_x_max); + } + + SECTION("Statistical test") + { + // Using multiple calls of sample_on_physical_groups(), we verify that each + // of the 4 bins of physical group 1 gets ~25% of the samples, using a loose + // tolerance (between 24% and 26%) + vector physical_groups = {1}; + + std::map bin_counts; + const int n_samples = 10000; + + for (int i = 0; i < n_samples; ++i) { + Position p; + int bin; + mesh.sample_on_physical_groups(p, bin, &seed, physical_groups); + bin_counts[bin]++; + } + + REQUIRE(bin_counts.size() == 4); + for (auto& [bin, count] : bin_counts) { + double fraction = static_cast(count) / n_samples; + CHECK(fraction > 0.24); + CHECK(fraction < 0.26); + } + } + + SECTION("Statistical test - with duplicated faces") + { + // Same test as above but with a duplicate in the physical groups map + vector physical_groups = {1}; + mesh.pg_map()[1].push_back(0); + + std::map bin_counts; + const int n_samples = 10000; + + for (int i = 0; i < n_samples; ++i) { + Position p; + int bin; + mesh.sample_on_physical_groups(p, bin, &seed, physical_groups); + bin_counts[bin]++; + } + + REQUIRE(bin_counts.size() == 4); + for (auto& [bin, count] : bin_counts) { + double fraction = static_cast(count) / n_samples; + CHECK(fraction > 0.24); + CHECK(fraction < 0.26); + } + } +} diff --git a/tests/regression_tests/microxs/__init.py__ b/tests/regression_tests/dnp_drift/__init__.py similarity index 100% rename from tests/regression_tests/microxs/__init.py__ rename to tests/regression_tests/dnp_drift/__init__.py 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/inputs_true.dat b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/inputs_true.dat new file mode 100644 index 00000000000..83863f891e2 --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/inputs_true.dat @@ -0,0 +1,71 @@ + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 4000 + 20 + 2 + + + 0.0 0.0 0.0 10.0 10.0 10.0 + + + true + + + + 1 + cell + 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 + + 1 + 2 + 3 4 + + + 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 + + RK4 + 0.1 + true + 2.0 + + + 2 2 2 + 0.0 0.0 0.0 + 10.0 10.0 10.0 + + + + + 1 + + + 1 + total + + + 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..a802a271245 --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.461752E+00 3.767986E-03 +tally 1: +3.985501E+01 +8.827004E+01 +3.992058E+01 +8.858998E+01 +4.008324E+01 +8.933599E+01 +4.010358E+01 +8.940097E+01 +4.022517E+01 +8.993904E+01 +4.013399E+01 +8.951738E+01 +4.019704E+01 +8.982139E+01 +4.001187E+01 +8.898035E+01 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..54de82b612d --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/cell-based/test.py @@ -0,0 +1,90 @@ +import openmc +import pytest +import numpy as np + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def dnp_drift_model(): + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_nuclide("O16", 3.0) + mat.add_nuclide("H1", 2.0) + mat.set_density("g/cm3", 5.0) + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (10.0, 10.0, 10.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Create the velocity field + velocity_field = openmc.VelocityField( + mesh=mesh, + values=np.array([1.0]*3*8).reshape(-1,3), + mapping="cell" + ) + + # 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]) + + # DNP drift + dnp_drift = openmc.DNPDrift( + velocity_field=velocity_field, + boundary_map={"inlet": [1], "outlet": [2], "wall": [3, 4]}, + physical_group_map={ + "face_ids": [ + 0, 24, 12, 36, 7, 31, 19, 43, 15, 39, 21, 45, 2, + 26, 8,32, 4, 16, 10, 22, 29, 41, 35, 47 + ], + "physical_groups": [ + 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3 + ]}, + integrator="RK4", + integrator_dt=0.1, + recycling=True, + external_travel_time=2.0 + ) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.inactive = 2 + settings.particles = 4000 + settings.dnp_drift = dnp_drift + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": 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_dnp_drift_regular_mesh_cell_based(dnp_drift_model): + harness = PyAPITestHarness('statepoint.20.h5', dnp_drift_model) + 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/inputs_true.dat b/tests/regression_tests/dnp_drift/regular_mesh/nodal/inputs_true.dat new file mode 100644 index 00000000000..8e9ae3029a8 --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/nodal/inputs_true.dat @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 4000 + 20 + 2 + + + 0.0 0.0 0.0 10.0 10.0 10.0 + + + true + + + + 1 + nodal + 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 + + 1 + 2 + 3 4 + + + 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 + + RK4 + 0.1 + true + 2.0 + + + 2 2 2 + 0.0 0.0 0.0 + 10.0 10.0 10.0 + + + + + 1 + + + 1 + total + + + 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..a802a271245 --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/nodal/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.461752E+00 3.767986E-03 +tally 1: +3.985501E+01 +8.827004E+01 +3.992058E+01 +8.858998E+01 +4.008324E+01 +8.933599E+01 +4.010358E+01 +8.940097E+01 +4.022517E+01 +8.993904E+01 +4.013399E+01 +8.951738E+01 +4.019704E+01 +8.982139E+01 +4.001187E+01 +8.898035E+01 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..3294e099ad0 --- /dev/null +++ b/tests/regression_tests/dnp_drift/regular_mesh/nodal/test.py @@ -0,0 +1,90 @@ +import openmc +import pytest +import numpy as np + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def dnp_drift_model(): + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_nuclide("O16", 3.0) + mat.add_nuclide("H1", 2.0) + mat.set_density("g/cm3", 5.0) + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0.0, 0.0, 0.0) + upper_right = (10.0, 10.0, 10.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Create the velocity field + velocity_field = openmc.VelocityField( + mesh=mesh, + values=np.array([1.0]*3*27).reshape(-1,3), + mapping="nodal" + ) + + # 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]) + + # DNP drift + dnp_drift = openmc.DNPDrift( + velocity_field=velocity_field, + boundary_map={"inlet": [1], "outlet": [2], "wall": [3, 4]}, + physical_group_map={ + "face_ids": [ + 0, 24, 12, 36, 7, 31, 19, 43, 15, 39, 21, 45, 2, + 26, 8,32, 4, 16, 10, 22, 29, 41, 35, 47 + ], + "physical_groups": [ + 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3 + ]}, + integrator="RK4", + integrator_dt=0.1, + recycling=True, + external_travel_time=2.0 + ) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.inactive = 2 + settings.particles = 4000 + settings.dnp_drift = dnp_drift + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": 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_dnp_drift_regular_mesh_nodal(dnp_drift_model): + harness = PyAPITestHarness('statepoint.20.h5', dnp_drift_model) + harness.main() diff --git a/tests/regression_tests/microxs/__init__.py b/tests/regression_tests/microxs/__init__.py new file mode 100644 index 00000000000..e69de29bb2d 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/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..e635c8d1353 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/nested_cubes/reflective_bc/inputs_true.dat @@ -0,0 +1,53 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..11ffbedc239 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..1557b23889d --- /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.feature_enabled('dagmc'), + reason="DAGMC CAD geometry 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..eb61a2fb0f2 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/nested_cubes/vacuum_bc/inputs_true.dat @@ -0,0 +1,53 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..2fd6efc48f7 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..6d2a3f26003 --- /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.feature_enabled('dagmc'), + reason="DAGMC CAD geometry 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..38afc34908e --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/single_cube/reflective_bc/inputs_true.dat @@ -0,0 +1,53 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..c10593fd9fd 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..0903b8d6806 --- /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.feature_enabled('dagmc'), + reason="DAGMC CAD geometry 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..0d4bb3004a5 --- /dev/null +++ b/tests/regression_tests/temperature_field/dagmc/single_cube/vacuum_bc/inputs_true.dat @@ -0,0 +1,53 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..5a90098dc7f 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..178540904aa --- /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.feature_enabled('dagmc'), + reason="DAGMC CAD geometry 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..00187d3b438 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/periodic_bc/inputs_true.dat @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + 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 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..74827e855ad --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/reflective_bc/inputs_true.dat @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + 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 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..1c61b6850a4 --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/nested_cubes/vacuum_bc/inputs_true.dat @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + 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 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..9a01200150d --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/periodic_bc/inputs_true.dat @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + 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 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..6bcdf66d05e --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/reflective_bc/inputs_true.dat @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + 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 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..8d0f27a67ae --- /dev/null +++ b/tests/regression_tests/temperature_field/lattices/single_cube/vacuum_bc/inputs_true.dat @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + 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 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..5db88316b54 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/inputs_true.dat @@ -0,0 +1,66 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..1048d71dc76 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/periodic_bc/results_true.dat @@ -0,0 +1,19 @@ +k-combined: +1.524339E+00 1.497787E-02 +tally 1: +1.656787E+00 +1.529189E-01 +1.926381E+00 +2.076675E-01 +1.244409E+00 +8.621588E-02 +1.605885E+00 +1.481062E-01 +1.337474E+00 +1.160469E-01 +1.560340E+00 +1.478962E-01 +1.589180E+00 +1.380172E-01 +1.460992E+00 +1.267146E-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..63b39d2dc39 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/reflective_bc/inputs_true.dat @@ -0,0 +1,66 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..3afa4fd8117 --- /dev/null +++ b/tests/regression_tests/temperature_field/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/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..5865489b8a2 --- /dev/null +++ b/tests/regression_tests/temperature_field/nested_cubes/vacuum_bc/inputs_true.dat @@ -0,0 +1,66 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + -5.0 -5.0 -5.0 10.0 10.0 10.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..5bd690224e8 --- /dev/null +++ b/tests/regression_tests/temperature_field/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/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..4efc0df8662 --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/periodic_bc/inputs_true.dat @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..449aadafeec --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/reflective_bc/inputs_true.dat @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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..b6d9d4a1fe6 --- /dev/null +++ b/tests/regression_tests/temperature_field/single_cube/vacuum_bc/inputs_true.dat @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + cell + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + 1 + + 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/unit_tests/test_dnp_drift.py b/tests/unit_tests/test_dnp_drift.py new file mode 100644 index 00000000000..4ff7ec8e324 --- /dev/null +++ b/tests/unit_tests/test_dnp_drift.py @@ -0,0 +1,136 @@ +"""Test the dnp_drift Python object.""" + +import openmc +import pytest +import numpy as np + + +@pytest.fixture() +def mesh(): + """2x2x2 regular 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) + return mesh + + +def test_xml_serialization(mesh, run_in_tmpdir): + """Test XML serialization.""" + velocity_field = openmc.VelocityField( + mesh=mesh, values=np.tile([0, 0, 10], (8, 1)), mapping="cell" + ) + + dnp_drift = openmc.DNPDrift( + velocity_field=velocity_field, + boundary_map={"inlet": [1], "outlet": [2], "wall": [3, 4]}, + physical_group_map={ + "face_ids": [ + 0, 24, 12, 36, 7, 31, 19, 43, 15, 39, 21, 45, 2, + 26, 8,32, 4, 16, 10, 22, 29, 41, 35, 47 + ], + "physical_groups": [ + 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3 + ]}, + integrator="RK4", + integrator_dt=0.1, + external_travel_time=2.0, + ) + + # Export settings + settings = openmc.Settings() + settings.dnp_drift = dnp_drift + settings.export_to_xml() + + # Read settings from xml file + read_settings = openmc.Settings.from_xml() + + # Check consistency + assert read_settings.dnp_drift.velocity_field == velocity_field + assert read_settings.dnp_drift.boundary_map == { + "inlet": [1], + "outlet": [2], + "wall": [3, 4], + } + assert read_settings.dnp_drift.physical_group_map == { + "face_ids": [ + 0, 24, 12, 36, 7, 31, 19, 43, 15, 39, 21, 45, 2, + 26, 8,32, 4, 16, 10, 22, 29, 41, 35, 47 + ], + "physical_groups": [ + 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3 + ]} + assert read_settings.dnp_drift.integrator == "RK4" + assert read_settings.dnp_drift.integrator_dt == 0.1 + assert read_settings.dnp_drift.external_travel_time == 2.0 + + +def test_invalid_mesh_type(): + """Check one non-implemented mesh type.""" + with pytest.raises(TypeError): + velocity_field = openmc.VelocityField( + openmc.SphericalMesh(r_grid=[0.0, 1.0]), [] + ) + + +def test_inconsistent_values_and_cells_sizes(mesh): + """Check inconsistency in number of values / mesh cells.""" + with pytest.raises(ValueError): + velocity_field = openmc.VelocityField( + mesh=mesh, values=np.tile([0, 0, 10], (9, 1)) + ) + + +def test_missing_boundary_map_keys(mesh): + """Check missing required boundary map keys.""" + velocity_field = openmc.VelocityField(mesh=mesh, values=np.tile([0, 0, 10], (8, 1))) + + with pytest.raises(ValueError): + openmc.DNPDrift( + velocity_field=velocity_field, + boundary_map={"inlet": [1]}, + physical_group_map={"face_ids": [0], "physical_groups": [1]} + ) + + +def test_invalid_boundary_map_key(mesh): + """Check invalid boundary map keys.""" + velocity_field = openmc.VelocityField(mesh=mesh, values=np.tile([0, 0, 10], (8, 1))) + + with pytest.raises(ValueError): + openmc.DNPDrift( + velocity_field=velocity_field, + boundary_map={"inlet": [1], "outlet": [2], "dnp_drift": [3]}, + physical_group_map={"face_ids": [0], "physical_groups": [1]} + ) + + +def test_invalid_integrator_dt(mesh): + """Check negative integrator time step value.""" + velocity_field = openmc.VelocityField(mesh=mesh, values=np.tile([0, 0, 10], (8, 1))) + + with pytest.raises(ValueError): + openmc.DNPDrift( + velocity_field=velocity_field, + boundary_map={"inlet": [1], "outlet": [2]}, + physical_group_map={"face_ids": [0], "physical_groups": [1]}, + integrator_dt=-0.1 + ) + + +def test_invalid_external_travel_time(mesh): + """Check negative external travel time.""" + velocity_field = openmc.VelocityField(mesh=mesh, values=np.tile([0, 0, 10], (8, 1))) + + with pytest.raises(ValueError): + openmc.DNPDrift( + velocity_field=velocity_field, + boundary_map={"inlet": [1], "outlet": [2]}, + physical_group_map={"face_ids": [0], "physical_groups": [1]}, + external_travel_time=-1.0 + ) diff --git a/tests/unit_tests/test_temperature_field.py b/tests/unit_tests/test_temperature_field.py new file mode 100644 index 00000000000..a0208e629e8 --- /dev/null +++ b/tests/unit_tests/test_temperature_field.py @@ -0,0 +1,63 @@ +"""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 = [1., 2., 3., 4., 5., 6., 7., 8.] + + # 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 + + +def test_invalid_mesh_type(): + """Check one non-implemented mesh type.""" + mesh = openmc.SphericalMesh(r_grid=[0., 1.0]) + values = [] + + with pytest.raises(TypeError): + temperature_field = openmc.TemperatureField(mesh, values) + + +def test_inconsistent_values_and_cells_sizes(mesh): + """Check inconsistency in number of values / mesh cells.""" + values = [0.0] * 7 + + with pytest.raises(ValueError): + temperature_field = openmc.TemperatureField(mesh, values) + + +def test_invalid_values(mesh): + """Check negative values.""" + values = [-1.0] * 8 + + with pytest.raises(ValueError): + temperature_field = openmc.TemperatureField(mesh, values)