Merge branch 'develop' into properties-exim

This commit is contained in:
Paul Romano 2021-07-14 15:18:31 -05:00
commit c9d84310d3
71 changed files with 4811 additions and 3405 deletions

3
.gitignore vendored
View file

@ -116,3 +116,6 @@ CMakeSettings.json
# Visual Studio Code configuration files
.vscode/
# Python pickle files
*.pkl

View file

@ -19,7 +19,7 @@ Each ``<surface>`` element can have the following attributes or sub-elements:
:name:
An optional string name to identify the surface in summary output
files. This string is limited to 52 characters for formatting purposes.
files.
*Default*: ""
@ -122,7 +122,6 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
:name:
An optional string name to identify the cell in summary output files.
This string is limmited to 52 characters for formatting purposes.
*Default*: ""
@ -235,7 +234,7 @@ the following attributes or sub-elements:
:name:
An optional string name to identify the lattice in summary output
files. This string is limited to 52 characters for formatting purposes.
files.
*Default*: ""
@ -301,7 +300,7 @@ the following attributes or sub-elements:
:name:
An optional string name to identify the hex_lattice in summary output
files. This string is limited to 52 characters for formatting purposes.
files.
*Default*: ""
@ -370,3 +369,40 @@ Here is an example of a properly defined 2d hexagonal lattice:
202
</universes>
</hex_lattice>
.. _dagmc_element:
----------------------------
``<dagmc_universe>`` Element
----------------------------
Each ``<dagmc_universe>`` element can have the following attributes or sub-elements:
:id:
A unique integer used to identify the universe.
*Default*: None
:name:
An optional string name to identify the surface in summary output
files.
*Default*: None
:auto_geom_ids:
Boolean value indicating whether the existing geometry IDs will be used or appended
to the existing ID space of natively defined OpenMC geometry entities.
*Default*: false
:auto_mat_ids:
Boolean value indicating whether the existing material IDs will be used or appended
to the existing ID space of natively defined OpenMC materials.
*Default*: false
:filename:
A required string indicating the file to be loaded representing the DAGMC universe.
*Default*: None

View file

@ -31,7 +31,7 @@ Each ``material`` element can have the following attributes or sub-elements:
:name:
An optional string name to identify the material in summary output
files. This string is limited to 52 characters for formatting purposes.
files.
*Default*: ""

View file

@ -88,14 +88,6 @@ you care. This element has the following attributes/sub-elements:
*Default*: 0.0
--------------------------------
``<dagmc>`` Element
--------------------------------
When the DAGMC mode is enabled, the OpenMC geometry will be read from the file
``dagmc.h5m``. If a :ref:`geometry.xml <io_geometry>` file is present with
``dagmc`` set to ``true``, it will be ignored.
----------------------------
``<delayed_photon_scaling>``
----------------------------

View file

@ -24,8 +24,6 @@ The current version of the summary file format is 6.0.
- **n_universes** (*int*) -- Number of unique universes in the
problem.
- **n_lattices** (*int*) -- Number of lattices in the problem.
- **dagmc** (*int*) -- Indicates that a DAGMC geometry was used
if present.
**/geometry/cells/cell <uid>/**
@ -49,6 +47,8 @@ The current version of the summary file format is 6.0.
- **lattice** (*int*) -- Unique ID of the lattice which fills the
cell. Only present if fill_type is set to 'lattice'.
- **region** (*char[]*) -- Region specification for the cell.
- **geom_type** (*char[]*) -- Type of geometry used to create the cell.
Either 'csg' or 'dagmc'.
**/geometry/surfaces/surface <uid>/**
@ -62,12 +62,27 @@ The current version of the summary file format is 6.0.
- **boundary_condition** (*char[]*) -- Boundary condition applied to
the surface. Can be 'transmission', 'vacuum', 'reflective', or
'periodic'.
- **geom_type** (*char[]*) -- Type of geometry used to create the cell.
Either 'csg' or 'dagmc'.
**/geometry/universes/universe <uid>/**
:Datasets:
- **cells** (*int[]*) -- Array of unique IDs of cells that appear in
the universe.
- **geom_type** (*char[]*) -- Type of geometry used to create the cell.
Either 'csg' or 'dagmc'.
- **filename** (*char[]*) -- Name of the DAGMC file representing this universe.
Only present for DAGMC Universes.
:Attributes:
- **auto_geom_ids** (*int*) -- ``1`` if geometry IDs of the DAGMC
model will be appended to the ID space of the natively defined
CSG geometry, ``0`` if the existing DAGMC IDs will be used.
- **auto_mat_ids** (*int*) -- ``1`` if UWUW material IDs of the DAGMC
model will be appended to the ID space of the natively defined
OpenMC materials, ``0`` if the existing UWUW IDs will be used.
**/geometry/lattices/lattice <uid>/**
@ -119,8 +134,8 @@ The current version of the summary file format is 6.0.
:Attributes: - **volume** (*double[]*) -- Volume of this material [cm^3]. Only
present if ``volume`` supplied
- **temperature** (*double[]*) -- Temperature of this material [K].
Only present in ``temperature`` supplied
- **depletable** (*int[]*) -- ``1`` if the material can be depleted,
Only present if ``temperature`` is supplied
- **depletable** (*int*) -- ``1`` if the material can be depleted,
``0`` otherwise. Always present
**/nuclides/**

View file

@ -31,7 +31,7 @@ The ``<tally>`` element accepts the following sub-elements:
:name:
An optional string name to identify the tally in summary output
files. This string is limited to 52 characters for formatting purposes.
files.
*Default*: ""

View file

@ -433,30 +433,64 @@ will handle creating the unverse::
Using CAD-based Geometry
--------------------------
OpenMC relies on the Direct Accelerated Geometry Monte Carlo toolkit (`DAGMC
<https://svalinn.github.io/DAGMC/>`_) to represent CAD-based geometry in a
surface mesh format. A DAGMC run can be enabled in OpenMC by setting the
``dagmc`` property to ``True`` in the model Settings either via the Python
:class:`openmc.settings` Python class::
Defining Geometry
-----------------
settings = openmc.Settings()
settings.dagmc = True
OpenMC relies on the `Direct Accelerated Geometry Monte Carlo`_ (DAGMC)
to represent CAD-based geometry in a surface mesh format. DAGMC geometries are
applied as universes in the OpenMC geometry file. A geometry represented
entirely by a DAGMC geometry will contain only the DAGMC universe. Using a
:class:`openmc.DAGMCUniverse` looks like the following::
or in the :ref:`settings.xml <io_settings>` file::
dag_univ = openmc.DAGMCUniverse(filename='dagmc.h5m')
geometry = openmc.Geometry(dag_univ)
geometry.export_to_xml()
<dagmc>true</dagmc>
The resulting ``geometry.xml`` file will be:
With ``dagmc`` set to true, OpenMC will load the DAGMC model (from a local file
named ``dagmc.h5m``) when initializing a simulation. If a `geometry.xml
<../io_formats/geometry.html>`_ is present as well, it will be ignored.
.. code-block:: xml
**Note:** DAGMC geometries used in OpenMC are currently required to be clean,
meaning that all surfaces have been `imprinted and merged
<https://svalinn.github.io/DAGMC/usersguide/trelis_workflow.html>`_
successfully and that the model is `watertight
<https://svalinn.github.io/DAGMC/usersguide/tools.html#make-watertight>`_. Future
implementations of DAGMC geometry will support small volume overlaps and
un-merged surfaces.
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<dagmc auto_ids="false" filename="dagmc.h5m" id="1" name="" />
</geometry>
DAGMC universes can also be used to fill CSG cells or lattice cells in a geometry::
cell.fill = dagmc_univ
It is important in these cases to understand the DAGMC model's position
with respect to the CSG geometry. DAGMC geometries can be plotted with
OpenMC to verify that the model matches one's expectations.
**Note:** DAGMC geometries used in OpenMC are currently required to be clean,
meaning that all surfaces have been `imprinted and merged
<https://svalinn.github.io/DAGMC/usersguide/trelis_workflow.html>`_
successfully and that the model is `watertight
<https://svalinn.github.io/DAGMC/usersguide/tools.html#make-watertight>`_. Future
implementations of DAGMC geometry will support small volume overlaps and
un-merged surfaces.
Cell, Surface, and Material IDs
-------------------------------
By default, DAGMC applies cell and surface IDs defined by the CAD engine that
the model originated in. If these IDs overlap with IDs in the CSG ID space,
this will result in an error. However, the ``auto_ids`` property of a DAGMC
universe can be set to set DAGMC cell and surface IDs by appending to the
existing CSG cell ID space in the OpenMC model.
Similar options exist for the material IDs of DAGMC models. If DAGMC material
assignments are based on natively defined OpenMC materials, no further work is
required. If DAGMC materials are assigned using the `University of Wisconsin
Unified Workflow`_ (UWUW), however, material IDs in the UWUW material library
may overlap with those used in the CSG geometry. In this case, overlaps in the
UWUW and OpenMC material ID space will cause an error. To automatically resolve
these ID overlaps, ``auto_ids`` can be set to ``True`` to append the UWUW
material IDs to the OpenMC material ID space.
.. _Direct Accelerated Geometry Monte Carlo: https://svalinn.github.io/DAGMC/
.. _University of Wisconsin Unified Workflow: https://svalinn.github.io/DAGMC/usersguide/uw2.html
-------------------------
Calculating Atoms Content

View file

@ -45,6 +45,23 @@ with the :envvar:`OPENMC_CROSS_SECTIONS` environment variable. It is recommended
to add a line in your ``.profile`` or ``.bash_profile`` setting the
:envvar:`OPENMC_CROSS_SECTIONS` environment variable.
RuntimeError: Failed to open HDF5 file with mode 'w': summary.h5
****************************************************************
This often occurs when working with the Python API and executing multiple OpenMC
runs in a script. If an :class:`openmc.StatePoint` is open in the Python interpreter,
the file handle of the statepoint file as well as the linked `summary.h5` file will
be unavailable for writing, causing this error to appear. To avoid this situation,
it is recommended that data be extracted from statepoint files in a context manager:
.. code-block:: python
with openmc.StatePoint('statepoint.10.h5') as sp:
k_eff = sp.k_combined
or that the :meth:`StatePoint.close` method is called before executing a subsequent
OpenMC run.
Geometry Debugging
******************

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

@ -243,7 +243,7 @@
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEVyEhLpgJFNv8T///98iBL0AAAAAWJLR0QDEQxM8gAAAAd0SU1FB+MHEhEzAhO1TdMAAAKlSURBVGje7ZrBscIwDETxwSWkn5TAgXCgBPqhBA6kyj/fDhCIJa2zZAwz0pmHpZWdSazd7Tw8PDw8PDw8vinCMBzW03HIsRLvhnv0HL7qD+IwjzXKzaNaxeEt9kz21RWEBV5XQbfka3pQWL4qgdLyNQkUcbwFT/FP4zjWt+D+++OY4laZQJgtnqNOwe5l9XkGmIIL/PEHUAGTeuc5P15wBbu34ucSIAXkX77h4xUtIBSXnxIAOhCLy98TANNfLj8lYBYQCuLPWmAWEBe9f90DUPdKy08JWB0U1HsoaAgQxPSnAgwBopz+VABQvoDnAnqTP0r8zealzfPcQqqAQSs/C6AKGLX0pwKs8uX0cwGaAHr6uYC9wSt46qDCB4RXBDTkMwWMhnxZQF3+s8pf1AZY5VsCWuVnAUQ+YPxB43X5koAiH035soCa/AaeBOw34m359AaQPCK/1oAAyJ8aIPBI+7QGRkD+3IBt+A6QPzeg34SH2pcauN+Kt9uXGljkse0jb6BP8AD+vwGKPLZ95A0UofbnDbAFj20/eQN+gD8h/LgRD25/8QCA2088AD/Oo8dPOoDo8ZMOoPPNeej4pwdAgUcfX9IDzHnnf5lnz88XnH/nSf4M8cIL7I+/P3yCP0G88P7W+v2z9ft36+8P9vuJ/X5r/f3Jfj83//5vff/R+v6Hvb9i78/Y+7vW94/N71/Z+2P2/pq9P2fv7+n5ATu/YOcn7PyGnR+x8yt6ftYN3PzOENCcH7LzS3Z+Ss9vO62DV5uPmgAXSz5+fs7O72n/QBQLwPwLrH+C9W/Q/hHWv8L6Z2j/ThZgvX+I9S/R/inWv8X6x2j/Guufo/17rH+Q9S/S/knWv0n7R2n/Kuufpf27tH+Y9i/vWP+0h4eHh4eHh8cW8QcxLJDBvLKoigAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAxOS0wNy0xOFQyMjo1MTowMi0wNTowMAMmdtQAAAAldEVYdGRhdGU6bW9kaWZ5ADIwMTktMDctMThUMjI6NTE6MDItMDU6MDBye85oAAAAAElFTkSuQmCC\n",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEV3C9wCswC1wFT///+Se5pWAAAAAWJLR0QDEQxM8gAAAAd0SU1FB+UGGhQuFt5F1RsAAAKlSURBVGje7ZrBscIwDETxwSWkn5TAgXCgBPqhBA6kyj/fDhCIJa2zZAwz0pmHpZWdSazd7Tw8PDw8PDw8vinCMBzW03HIsRLvhnv0HL7qD+IwjzXKzaNaxeEt9kz21RWEBV5XQbfka3pQWL4qgdLyNQkUcbwFT/FP4zjWt+D+++OY4laZQJgtnqNOwe5l9XkGmIIL/PEHUAGTeuc5P15wBbu34ucSIAXkX77h4xUtIBSXnxIAOhCLy98TANNfLj8lYBYQCuLPWmAWEBe9f90DUPdKy08JWB0U1HsoaAgQxPSnAgwBopz+VABQvoDnAnqTP0r8zealzfPcQqqAQSs/C6AKGLX0pwKs8uX0cwGaAHr6uYC9wSt46qDCB4RXBDTkMwWMhnxZQF3+s8pf1AZY5VsCWuVnAUQ+YPxB43X5koAiH035soCa/AaeBOw34m359AaQPCK/1oAAyJ8aIPBI+7QGRkD+3IBt+A6QPzeg34SH2pcauN+Kt9uXGljkse0jb6BP8AD+vwGKPLZ95A0UofbnDbAFj20/eQN+gD8h/LgRD25/8QCA2088AD/Oo8dPOoDo8ZMOoPPNeej4pwdAgUcfX9IDzHnnf5lnz88XnH/nSf4M8cIL7I+/P3yCP0G88P7W+v2z9ft36+8P9vuJ/X5r/f3Jfj83//5vff/R+v6Hvb9i78/Y+7vW94/N71/Z+2P2/pq9P2fv7+n5ATu/YOcn7PyGnR+x8yt6ftYN3PzOENCcH7LzS3Z+Ss9vO62DV5uPmgAXSz5+fs7O72n/QBQLwPwLrH+C9W/Q/hHWv8L6Z2j/ThZgvX+I9S/R/inWv8X6x2j/Guufo/17rH+Q9S/S/knWv0n7R2n/Kuufpf27tH+Y9i/vWP+0h4eHh4eHh8cW8QcxLJDBvLKoigAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAyMS0wNi0yNlQyMDo0NjoyMiswMDowMIjWIT4AAAAldEVYdGRhdGU6bW9kaWZ5ADIwMjEtMDYtMjZUMjA6NDY6MjIrMDA6MDD5i5mCAAAAAElFTkSuQmCC\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
@ -421,11 +421,11 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=6.\n",
"/home/shriwise/opt/openmc/openmc/openmc/mixin.py:68: IDWarning: Another Filter instance already exists with id=6.\n",
" warn(msg, IDWarning)\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=3.\n",
"/home/shriwise/opt/openmc/openmc/openmc/mixin.py:68: IDWarning: Another Filter instance already exists with id=3.\n",
" warn(msg, IDWarning)\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=2.\n",
"/home/shriwise/opt/openmc/openmc/openmc/mixin.py:68: IDWarning: Another Filter instance already exists with id=2.\n",
" warn(msg, IDWarning)\n"
]
}
@ -478,78 +478,82 @@
" %%%%%%%%%%%\n",
"\n",
" | The OpenMC Monte Carlo Code\n",
" Copyright | 2011-2019 MIT and OpenMC contributors\n",
" License | http://openmc.readthedocs.io/en/latest/license.html\n",
" Version | 0.11.0-dev\n",
" Git SHA1 | 61c911cffdae2406f9f4bc667a9a6954748bb70c\n",
" Date/Time | 2019-07-18 22:51:02\n",
" OpenMP Threads | 4\n",
" Copyright | 2011-2021 MIT and OpenMC contributors\n",
" License | https://docs.openmc.org/en/latest/license.html\n",
" Version | 0.13.0-dev\n",
" Git SHA1 | 3dd81a1316ac3b5a0633e4b7a290f3bc97a066d9\n",
" Date/Time | 2021-06-26 15:46:22\n",
" OpenMP Threads | 2\n",
"\n",
" Reading settings XML file...\n",
" Reading cross sections XML file...\n",
" Reading materials XML file...\n",
" Reading geometry XML file...\n",
" Reading U235 from /opt/data/hdf5/nndc_hdf5_v15/U235.h5\n",
" Reading U238 from /opt/data/hdf5/nndc_hdf5_v15/U238.h5\n",
" Reading O16 from /opt/data/hdf5/nndc_hdf5_v15/O16.h5\n",
" Reading H1 from /opt/data/hdf5/nndc_hdf5_v15/H1.h5\n",
" Reading B10 from /opt/data/hdf5/nndc_hdf5_v15/B10.h5\n",
" Reading Zr90 from /opt/data/hdf5/nndc_hdf5_v15/Zr90.h5\n",
" Maximum neutron transport energy: 20000000.000000 eV for U235\n",
" Reading U235 from /home/shriwise/opt/openmc/xs/nndc_hdf5/U235.h5\n",
" Reading U238 from /home/shriwise/opt/openmc/xs/nndc_hdf5/U238.h5\n",
" Reading O16 from /home/shriwise/opt/openmc/xs/nndc_hdf5/O16.h5\n",
" Reading H1 from /home/shriwise/opt/openmc/xs/nndc_hdf5/H1.h5\n",
" Reading B10 from /home/shriwise/opt/openmc/xs/nndc_hdf5/B10.h5\n",
" Reading Zr90 from /home/shriwise/opt/openmc/xs/nndc_hdf5/Zr90.h5\n",
" Minimum neutron data temperature: 294.0 K\n",
" Maximum neutron data temperature: 294.0 K\n",
" Reading tallies XML file...\n",
" Preparing distributed cell instances...\n",
" Writing summary.h5 file...\n",
" Maximum neutron transport energy: 20000000.0 eV for U235\n",
" Initializing source particles...\n",
"\n",
" ====================> K EIGENVALUE SIMULATION <====================\n",
"\n",
" Bat./Gen. k Average k\n",
" ========= ======== ====================\n",
" 1/1 0.96168\n",
" 2/1 0.96651\n",
" 3/1 1.00678\n",
" 4/1 0.98773\n",
" 5/1 1.01883\n",
" 6/1 1.02959\n",
" 7/1 0.99859 1.01409 +/- 0.01550\n",
" 8/1 1.03441 1.02086 +/- 0.01123\n",
" 9/1 1.06097 1.03089 +/- 0.01279\n",
" 10/1 1.06132 1.03698 +/- 0.01163\n",
" 11/1 1.04687 1.03863 +/- 0.00964\n",
" 12/1 1.02982 1.03737 +/- 0.00824\n",
" 13/1 1.03520 1.03710 +/- 0.00714\n",
" 14/1 0.99508 1.03243 +/- 0.00784\n",
" 15/1 1.03987 1.03317 +/- 0.00705\n",
" 16/1 1.02743 1.03265 +/- 0.00640\n",
" 17/1 1.02975 1.03241 +/- 0.00585\n",
" 18/1 0.99671 1.02966 +/- 0.00604\n",
" 19/1 1.02040 1.02900 +/- 0.00563\n",
" 20/1 1.02024 1.02842 +/- 0.00527\n",
" 1/1 0.99327\n",
" 2/1 1.05535\n",
" 3/1 1.02165\n",
" 4/1 1.02898\n",
" 5/1 1.02521\n",
" 6/1 1.04782\n",
" 7/1 1.07014 1.05898 +/- 0.01116\n",
" 8/1 1.07724 1.06507 +/- 0.00886\n",
" 9/1 1.06268 1.06447 +/- 0.00630\n",
" 10/1 0.99628 1.05083 +/- 0.01448\n",
" 11/1 1.07257 1.05446 +/- 0.01237\n",
" 12/1 1.01603 1.04897 +/- 0.01181\n",
" 13/1 1.01047 1.04415 +/- 0.01130\n",
" 14/1 1.03639 1.04329 +/- 0.01000\n",
" 15/1 1.00699 1.03966 +/- 0.00966\n",
" 16/1 0.99098 1.03524 +/- 0.00979\n",
" 17/1 1.03990 1.03562 +/- 0.00895\n",
" 18/1 1.02076 1.03448 +/- 0.00831\n",
" 19/1 0.99685 1.03179 +/- 0.00815\n",
" 20/1 1.04833 1.03290 +/- 0.00767\n",
" Creating state point statepoint.20.h5...\n",
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 3.4427e-01 seconds\n",
" Reading cross sections = 3.1628e-01 seconds\n",
" Total time in simulation = 3.7319e+00 seconds\n",
" Time in transport only = 3.6302e+00 seconds\n",
" Time in inactive batches = 4.9601e-01 seconds\n",
" Time in active batches = 3.2359e+00 seconds\n",
" Time synchronizing fission bank = 2.8100e-03 seconds\n",
" Sampling source sites = 2.4682e-03 seconds\n",
" SEND/RECV source sites = 3.2484e-04 seconds\n",
" Time accumulating tallies = 4.4538e-05 seconds\n",
" Total time for finalization = 9.3656e-04 seconds\n",
" Total time elapsed = 4.0859e+00 seconds\n",
" Calculation Rate (inactive) = 25201.2 particles/second\n",
" Calculation Rate (active) = 11588.7 particles/second\n",
" Total time for initialization = 2.5195e-01 seconds\n",
" Reading cross sections = 2.3830e-01 seconds\n",
" Total time in simulation = 3.6172e+00 seconds\n",
" Time in transport only = 3.6041e+00 seconds\n",
" Time in inactive batches = 4.5282e-01 seconds\n",
" Time in active batches = 3.1644e+00 seconds\n",
" Time synchronizing fission bank = 2.4821e-03 seconds\n",
" Sampling source sites = 2.0826e-03 seconds\n",
" SEND/RECV source sites = 3.8781e-04 seconds\n",
" Time accumulating tallies = 2.8497e-04 seconds\n",
" Time writing statepoints = 8.5742e-03 seconds\n",
" Total time for finalization = 3.1545e-04 seconds\n",
" Total time elapsed = 3.8763e+00 seconds\n",
" Calculation Rate (inactive) = 27604.9 particles/second\n",
" Calculation Rate (active) = 11850.8 particles/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
" k-effective (Collision) = 1.02889 +/- 0.00492\n",
" k-effective (Track-length) = 1.02842 +/- 0.00527\n",
" k-effective (Absorption) = 1.02637 +/- 0.00349\n",
" Combined k-effective = 1.02700 +/- 0.00291\n",
" Leakage Fraction = 0.01717 +/- 0.00107\n",
" k-effective (Collision) = 1.03497 +/- 0.00681\n",
" k-effective (Track-length) = 1.03290 +/- 0.00767\n",
" k-effective (Absorption) = 1.02663 +/- 0.00590\n",
" Combined k-effective = 1.03085 +/- 0.00535\n",
" Leakage Fraction = 0.01728 +/- 0.00077\n",
"\n"
]
}
@ -631,8 +635,8 @@
" <th>0</th>\n",
" <td>total</td>\n",
" <td>(nu-fission / (absorption + current))</td>\n",
" <td>1.023002</td>\n",
" <td>0.006647</td>\n",
" <td>1.025847</td>\n",
" <td>0.009998</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -640,7 +644,7 @@
],
"text/plain": [
" nuclide score mean std. dev.\n",
"0 total (nu-fission / (absorption + current)) 1.02e+00 6.65e-03"
"0 total (nu-fission / (absorption + current)) 1.03e+00 1.00e-02"
]
},
"execution_count": 21,
@ -712,8 +716,8 @@
" <td>0.625</td>\n",
" <td>total</td>\n",
" <td>((absorption + current) / (absorption + current))</td>\n",
" <td>0.694368</td>\n",
" <td>0.004606</td>\n",
" <td>0.695924</td>\n",
" <td>0.007275</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -724,7 +728,7 @@
"0 0.00e+00 6.25e-01 total \n",
"\n",
" score mean std. dev. \n",
"0 ((absorption + current) / (absorption + current)) 6.94e-01 4.61e-03 "
"0 ((absorption + current) / (absorption + current)) 6.96e-01 7.27e-03 "
]
},
"execution_count": 22,
@ -790,8 +794,8 @@
" <td>0.625</td>\n",
" <td>total</td>\n",
" <td>(nu-fission / nu-fission)</td>\n",
" <td>1.203099</td>\n",
" <td>0.009615</td>\n",
" <td>1.203449</td>\n",
" <td>0.01381</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -802,7 +806,7 @@
"0 0.00e+00 6.25e-01 total (nu-fission / nu-fission) \n",
"\n",
" mean std. dev. \n",
"0 1.20e+00 9.61e-03 "
"0 1.20e+00 1.38e-02 "
]
},
"execution_count": 23,
@ -869,8 +873,8 @@
" <td>1</td>\n",
" <td>total</td>\n",
" <td>(absorption / absorption)</td>\n",
" <td>0.749423</td>\n",
" <td>0.006089</td>\n",
" <td>0.74975</td>\n",
" <td>0.009032</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -881,7 +885,7 @@
"0 0.00e+00 6.25e-01 1 total (absorption / absorption) \n",
"\n",
" mean std. dev. \n",
"0 7.49e-01 6.09e-03 "
"0 7.50e-01 9.03e-03 "
]
},
"execution_count": 24,
@ -946,8 +950,8 @@
" <td>1</td>\n",
" <td>total</td>\n",
" <td>(nu-fission / absorption)</td>\n",
" <td>1.663727</td>\n",
" <td>0.014403</td>\n",
" <td>1.663575</td>\n",
" <td>0.020557</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -958,7 +962,7 @@
"0 0.00e+00 6.25e-01 1 total (nu-fission / absorption) \n",
"\n",
" mean std. dev. \n",
"0 1.66e+00 1.44e-02 "
"0 1.66e+00 2.06e-02 "
]
},
"execution_count": 25,
@ -1020,8 +1024,8 @@
" <td>0.625</td>\n",
" <td>total</td>\n",
" <td>((absorption + current) / (absorption + current))</td>\n",
" <td>0.984668</td>\n",
" <td>0.005509</td>\n",
" <td>0.984639</td>\n",
" <td>0.00883</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1032,7 +1036,7 @@
"0 0.00e+00 6.25e-01 total \n",
"\n",
" score mean std. dev. \n",
"0 ((absorption + current) / (absorption + current)) 9.85e-01 5.51e-03 "
"0 ((absorption + current) / (absorption + current)) 9.85e-01 8.83e-03 "
]
},
"execution_count": 26,
@ -1093,8 +1097,8 @@
" <td>0.625</td>\n",
" <td>total</td>\n",
" <td>(absorption / (absorption + current))</td>\n",
" <td>0.997439</td>\n",
" <td>0.007548</td>\n",
" <td>0.997372</td>\n",
" <td>0.011707</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1105,7 +1109,7 @@
"0 0.00e+00 6.25e-01 total \n",
"\n",
" score mean std. dev. \n",
"0 (absorption / (absorption + current)) 9.97e-01 7.55e-03 "
"0 (absorption / (absorption + current)) 9.97e-01 1.17e-02 "
]
},
"execution_count": 27,
@ -1168,8 +1172,8 @@
" <td>1</td>\n",
" <td>total</td>\n",
" <td>(((((((absorption + current) / (absorption + c...</td>\n",
" <td>1.023002</td>\n",
" <td>0.018791</td>\n",
" <td>1.025847</td>\n",
" <td>0.028224</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1180,7 +1184,7 @@
"0 0.00e+00 6.25e-01 1 total \n",
"\n",
" score mean std. dev. \n",
"0 (((((((absorption + current) / (absorption + c... 1.02e+00 1.88e-02 "
"0 (((((((absorption + current) / (absorption + c... 1.03e+00 2.82e-02 "
]
},
"execution_count": 28,
@ -1260,8 +1264,8 @@
" <td>6.250000e-01</td>\n",
" <td>(U238 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>6.659486e-07</td>\n",
" <td>5.627975e-09</td>\n",
" <td>6.670761e-07</td>\n",
" <td>7.788171e-09</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
@ -1270,8 +1274,8 @@
" <td>6.250000e-01</td>\n",
" <td>(U238 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>2.099901e-01</td>\n",
" <td>1.748379e-03</td>\n",
" <td>2.099931e-01</td>\n",
" <td>2.336727e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
@ -1280,8 +1284,8 @@
" <td>6.250000e-01</td>\n",
" <td>(U235 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>3.566329e-01</td>\n",
" <td>3.030782e-03</td>\n",
" <td>3.571972e-01</td>\n",
" <td>4.203420e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
@ -1290,8 +1294,8 @@
" <td>6.250000e-01</td>\n",
" <td>(U235 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>5.555466e-03</td>\n",
" <td>4.635318e-05</td>\n",
" <td>5.555637e-03</td>\n",
" <td>6.200519e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
@ -1300,8 +1304,8 @@
" <td>2.000000e+07</td>\n",
" <td>(U238 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>7.251304e-03</td>\n",
" <td>5.161998e-05</td>\n",
" <td>7.269233e-03</td>\n",
" <td>6.814579e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
@ -1310,8 +1314,8 @@
" <td>2.000000e+07</td>\n",
" <td>(U238 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>2.272661e-01</td>\n",
" <td>9.576939e-04</td>\n",
" <td>2.276519e-01</td>\n",
" <td>7.674381e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
@ -1320,8 +1324,8 @@
" <td>2.000000e+07</td>\n",
" <td>(U235 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>7.920169e-03</td>\n",
" <td>5.751231e-05</td>\n",
" <td>8.069479e-03</td>\n",
" <td>5.324215e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
@ -1330,8 +1334,8 @@
" <td>2.000000e+07</td>\n",
" <td>(U235 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>3.358280e-03</td>\n",
" <td>1.341281e-05</td>\n",
" <td>3.363165e-03</td>\n",
" <td>1.147434e-05</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1349,14 +1353,14 @@
"7 1 6.25e-01 2.00e+07 (U235 / total) \n",
"\n",
" score mean std. dev. \n",
"0 (nu-fission / flux) 6.66e-07 5.63e-09 \n",
"1 (scatter / flux) 2.10e-01 1.75e-03 \n",
"2 (nu-fission / flux) 3.57e-01 3.03e-03 \n",
"3 (scatter / flux) 5.56e-03 4.64e-05 \n",
"4 (nu-fission / flux) 7.25e-03 5.16e-05 \n",
"5 (scatter / flux) 2.27e-01 9.58e-04 \n",
"6 (nu-fission / flux) 7.92e-03 5.75e-05 \n",
"7 (scatter / flux) 3.36e-03 1.34e-05 "
"0 (nu-fission / flux) 6.67e-07 7.79e-09 \n",
"1 (scatter / flux) 2.10e-01 2.34e-03 \n",
"2 (nu-fission / flux) 3.57e-01 4.20e-03 \n",
"3 (scatter / flux) 5.56e-03 6.20e-05 \n",
"4 (nu-fission / flux) 7.27e-03 6.81e-05 \n",
"5 (scatter / flux) 2.28e-01 7.67e-04 \n",
"6 (nu-fission / flux) 8.07e-03 5.32e-05 \n",
"7 (scatter / flux) 3.36e-03 1.15e-05 "
]
},
"execution_count": 30,
@ -1385,11 +1389,11 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[6.65948580e-07]\n",
" [3.56632881e-01]]\n",
"[[[6.67076123e-07]\n",
" [3.57197248e-01]]\n",
"\n",
" [[7.25130446e-03]\n",
" [7.92016892e-03]]]\n"
" [[7.26923324e-03]\n",
" [8.06947893e-03]]]\n"
]
}
],
@ -1415,9 +1419,9 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[0.00555547]]\n",
"[[[0.00555564]]\n",
"\n",
" [[0.00335828]]]\n"
" [[0.00336316]]]\n"
]
}
],
@ -1437,8 +1441,8 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[0.22726611]\n",
" [0.00335828]]]\n"
"[[[0.22765194]\n",
" [0.00336316]]]\n"
]
}
],
@ -1501,7 +1505,7 @@
" <td>U238</td>\n",
" <td>nu-fission</td>\n",
" <td>0.000002</td>\n",
" <td>9.679304e-09</td>\n",
" <td>1.382881e-08</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
@ -1510,8 +1514,8 @@
" <td>6.250000e-01</td>\n",
" <td>U235</td>\n",
" <td>nu-fission</td>\n",
" <td>0.854805</td>\n",
" <td>5.239673e-03</td>\n",
" <td>0.858278</td>\n",
" <td>7.512185e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
@ -1520,8 +1524,8 @@
" <td>2.000000e+07</td>\n",
" <td>U238</td>\n",
" <td>nu-fission</td>\n",
" <td>0.082978</td>\n",
" <td>5.346135e-04</td>\n",
" <td>0.082753</td>\n",
" <td>7.469156e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
@ -1530,8 +1534,8 @@
" <td>2.000000e+07</td>\n",
" <td>U235</td>\n",
" <td>nu-fission</td>\n",
" <td>0.090632</td>\n",
" <td>5.981942e-04</td>\n",
" <td>0.091863</td>\n",
" <td>5.596607e-04</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1540,15 +1544,15 @@
"text/plain": [
" cell energy low [eV] energy high [eV] nuclide score mean \\\n",
"0 1 0.00e+00 6.25e-01 U238 nu-fission 1.60e-06 \n",
"1 1 0.00e+00 6.25e-01 U235 nu-fission 8.55e-01 \n",
"2 1 6.25e-01 2.00e+07 U238 nu-fission 8.30e-02 \n",
"3 1 6.25e-01 2.00e+07 U235 nu-fission 9.06e-02 \n",
"1 1 0.00e+00 6.25e-01 U235 nu-fission 8.58e-01 \n",
"2 1 6.25e-01 2.00e+07 U238 nu-fission 8.28e-02 \n",
"3 1 6.25e-01 2.00e+07 U235 nu-fission 9.19e-02 \n",
"\n",
" std. dev. \n",
"0 9.68e-09 \n",
"1 5.24e-03 \n",
"2 5.35e-04 \n",
"3 5.98e-04 "
"0 1.38e-08 \n",
"1 7.51e-03 \n",
"2 7.47e-04 \n",
"3 5.60e-04 "
]
},
"execution_count": 34,
@ -1605,8 +1609,8 @@
" <td>1.080060e-01</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>4.541188</td>\n",
" <td>0.025230</td>\n",
" <td>4.542677</td>\n",
" <td>0.037555</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
@ -1615,8 +1619,8 @@
" <td>1.166529e+00</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>2.001332</td>\n",
" <td>0.006754</td>\n",
" <td>2.019002</td>\n",
" <td>0.011992</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
@ -1625,8 +1629,8 @@
" <td>1.259921e+01</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>1.639292</td>\n",
" <td>0.011374</td>\n",
" <td>1.622688</td>\n",
" <td>0.011889</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
@ -1635,8 +1639,8 @@
" <td>1.360790e+02</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>1.821633</td>\n",
" <td>0.009590</td>\n",
" <td>1.834317</td>\n",
" <td>0.012322</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
@ -1645,8 +1649,8 @@
" <td>1.469734e+03</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>2.032395</td>\n",
" <td>0.009953</td>\n",
" <td>2.040064</td>\n",
" <td>0.012967</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
@ -1655,8 +1659,8 @@
" <td>1.587401e+04</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>2.120745</td>\n",
" <td>0.011090</td>\n",
" <td>2.107758</td>\n",
" <td>0.012781</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
@ -1665,8 +1669,8 @@
" <td>1.714488e+05</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>2.181709</td>\n",
" <td>0.013602</td>\n",
" <td>2.175480</td>\n",
" <td>0.014165</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
@ -1675,8 +1679,8 @@
" <td>1.851749e+06</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>2.013644</td>\n",
" <td>0.009219</td>\n",
" <td>1.983382</td>\n",
" <td>0.012931</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
@ -1685,8 +1689,8 @@
" <td>2.000000e+07</td>\n",
" <td>H1</td>\n",
" <td>scatter</td>\n",
" <td>0.372640</td>\n",
" <td>0.002903</td>\n",
" <td>0.372359</td>\n",
" <td>0.003689</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1695,25 +1699,25 @@
"text/plain": [
" cell energy low [eV] energy high [eV] nuclide score mean \\\n",
"0 3 1.00e-02 1.08e-01 H1 scatter 4.54e+00 \n",
"1 3 1.08e-01 1.17e+00 H1 scatter 2.00e+00 \n",
"2 3 1.17e+00 1.26e+01 H1 scatter 1.64e+00 \n",
"3 3 1.26e+01 1.36e+02 H1 scatter 1.82e+00 \n",
"4 3 1.36e+02 1.47e+03 H1 scatter 2.03e+00 \n",
"5 3 1.47e+03 1.59e+04 H1 scatter 2.12e+00 \n",
"1 3 1.08e-01 1.17e+00 H1 scatter 2.02e+00 \n",
"2 3 1.17e+00 1.26e+01 H1 scatter 1.62e+00 \n",
"3 3 1.26e+01 1.36e+02 H1 scatter 1.83e+00 \n",
"4 3 1.36e+02 1.47e+03 H1 scatter 2.04e+00 \n",
"5 3 1.47e+03 1.59e+04 H1 scatter 2.11e+00 \n",
"6 3 1.59e+04 1.71e+05 H1 scatter 2.18e+00 \n",
"7 3 1.71e+05 1.85e+06 H1 scatter 2.01e+00 \n",
"8 3 1.85e+06 2.00e+07 H1 scatter 3.73e-01 \n",
"7 3 1.71e+05 1.85e+06 H1 scatter 1.98e+00 \n",
"8 3 1.85e+06 2.00e+07 H1 scatter 3.72e-01 \n",
"\n",
" std. dev. \n",
"0 2.52e-02 \n",
"1 6.75e-03 \n",
"2 1.14e-02 \n",
"3 9.59e-03 \n",
"4 9.95e-03 \n",
"5 1.11e-02 \n",
"6 1.36e-02 \n",
"7 9.22e-03 \n",
"8 2.90e-03 "
"0 3.76e-02 \n",
"1 1.20e-02 \n",
"2 1.19e-02 \n",
"3 1.23e-02 \n",
"4 1.30e-02 \n",
"5 1.28e-02 \n",
"6 1.42e-02 \n",
"7 1.29e-02 \n",
"8 3.69e-03 "
]
},
"execution_count": 35,
@ -1728,6 +1732,16 @@
" filters=[openmc.CellFilter], filter_bins=[(moderator_cell.id,)])\n",
"slice_test.get_pandas_dataframe()"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"# Close the statepoint file as a matter of best practice\n",
"sp.close()"
]
}
],
"metadata": {
@ -1746,7 +1760,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
"version": "3.7.3"
}
},
"nbformat": 4,

View file

@ -181,14 +181,15 @@
"\n",
"settings.run_mode = \"fixed source\"\n",
"settings.batches = 10\n",
"settings.particles = 100"
"settings.particles = 100\n",
"settings.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And we'll indicate that we're using a CAD-based geometry."
"Next we'll apply the DAGMC model as the root universe of the geometry."
]
},
{
@ -197,9 +198,9 @@
"metadata": {},
"outputs": [],
"source": [
"settings.dagmc = True\n",
"\n",
"settings.export_to_xml()"
"dagmc_univ = openmc.DAGMCUniverse(filename='dagmc.h5m')\n",
"geometry = openmc.Geometry(root=dagmc_univ)\n",
"geometry.export_to_xml()"
]
},
{
@ -246,15 +247,15 @@
" Copyright | 2011-2020 MIT and OpenMC contributors\n",
" License | https://docs.openmc.org/en/latest/license.html\n",
" Version | 0.12.1-dev\n",
" Git SHA1 | 8ae407c90e927af62bfdc8f150e96bfd5d7b2ec2\n",
" Date/Time | 2021-01-06 09:17:05\n",
" Git SHA1 | 9fb298b039bcd447c1a745cd9fea47f0d04eebdf\n",
" Date/Time | 2021-03-04 13:34:25\n",
" MPI Processes | 1\n",
" OpenMP Threads | 2\n",
" OpenMP Threads | 4\n",
"\n",
" Reading settings XML file...\n",
" Reading cross sections XML file...\n",
" Reading materials XML file...\n",
" Reading DAGMC geometry...\n",
" Reading geometry XML file...\n",
"Set overlap thickness = 0\n",
"Set numerical precision = 0.001\n",
"Loading file dagmc.h5m\n",
@ -317,16 +318,16 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 2.4196e+01 seconds\n",
" Reading cross sections = 2.1428e+00 seconds\n",
" Total time in simulation = 1.9512e-01 seconds\n",
" Time in transport only = 1.9269e-01 seconds\n",
" Time in active batches = 1.9512e-01 seconds\n",
" Time accumulating tallies = 2.0220e-06 seconds\n",
" Time writing statepoints = 2.2461e-03 seconds\n",
" Total time for finalization = 1.8940e-06 seconds\n",
" Total time elapsed = 2.4401e+01 seconds\n",
" Calculation Rate (active) = 5125.02 particles/second\n",
" Total time for initialization = 1.7970e+02 seconds\n",
" Reading cross sections = 1.8083e+00 seconds\n",
" Total time in simulation = 7.2061e-01 seconds\n",
" Time in transport only = 7.1887e-01 seconds\n",
" Time in active batches = 7.2061e-01 seconds\n",
" Time accumulating tallies = 2.7050e-06 seconds\n",
" Time writing statepoints = 1.5381e-03 seconds\n",
" Total time for finalization = 1.1780e-06 seconds\n",
" Total time elapsed = 1.8043e+02 seconds\n",
" Calculation Rate (active) = 1387.71 particles/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
@ -352,7 +353,7 @@
"metadata": {},
"outputs": [],
"source": [
"unstructured_mesh = openmc.UnstructuredMesh(\"manifold.h5m\", library='moab')\n",
"unstructured_mesh = openmc.UnstructuredMesh(\"manifold.h5m\")\n",
"\n",
"mesh_filter = openmc.MeshFilter(unstructured_mesh)\n",
"\n",
@ -496,7 +497,7 @@
"text": [
"<?xml version='1.0' encoding='utf-8'?>\r\n",
"<tallies>\r\n",
" <mesh id=\"1\" library=\"moab\" type=\"unstructured\">\r\n",
" <mesh id=\"1\" type=\"unstructured\">\r\n",
" <filename>manifold.h5m</filename>\r\n",
" </mesh>\r\n",
" <filter id=\"1\" type=\"mesh\">\r\n",

View file

@ -10,7 +10,6 @@
#include <gsl/gsl>
#include "hdf5.h"
#include "pugixml.hpp"
#include "dagmc.h"
#include "openmc/constants.h"
#include "openmc/memory.h" // for unique_ptr
@ -63,16 +62,25 @@ namespace model {
class Universe
{
public:
int32_t id_; //!< Unique ID
vector<int32_t> cells_; //!< Cells within this universe
//! \brief Write universe information to an HDF5 group.
//! \param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
virtual void to_hdf5(hid_t group_id) const;
virtual bool find_cell(Particle &p) const;
BoundingBox bounding_box() const;
const GeometryType& geom_type() const { return geom_type_; }
GeometryType& geom_type() { return geom_type_; }
unique_ptr<UniversePartitioner> partitioner_;
private:
GeometryType geom_type_ = GeometryType::CSG;
};
//==============================================================================
@ -118,7 +126,17 @@ public:
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const = 0;
void to_hdf5(hid_t group_id) const;
virtual void to_hdf5_inner(hid_t group_id) const = 0;
//! Export physical properties to HDF5
//! \param[in] group HDF5 group to read from
void export_properties_hdf5(hid_t group) const;
//! Import physical properties from HDF5
//! \param[in] group HDF5 group to write to
void import_properties_hdf5(hid_t group);
//! Get the BoundingBox for this cell.
virtual BoundingBox bounding_box() const = 0;
@ -153,14 +171,6 @@ public:
//! \return Map with cell indexes as keys and instances as values
std::unordered_map<int32_t, vector<int32_t>> get_contained_cells() const;
//! Export physical properties to HDF5
//! \param[in] group HDF5 group to read from
void export_properties_hdf5(hid_t group) const;
//! Import physical properties from HDF5
//! \param[in] group HDF5 group to write to
void import_properties_hdf5(hid_t group);
protected:
void get_contained_cells_inner(
std::unordered_map<int32_t, vector<int32_t>>& contained_cells,
@ -172,10 +182,11 @@ public:
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
Fill type_; //!< Material, universe, or lattice
Fill type_; //!< Material, universe, or lattice
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
int32_t n_instances_{0}; //!< Number of instances of this cell
GeometryType geom_type_; //!< Geometric representation type (CSG, DAGMC)
//! \brief Index corresponding to this cell in distribcell arrays
int distribcell_index_{C_NONE};
@ -233,7 +244,7 @@ public:
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface, Particle* p) const;
void to_hdf5(hid_t group_id) const;
void to_hdf5_inner(hid_t group_id) const override;
BoundingBox bounding_box() const;
@ -262,28 +273,6 @@ protected:
vector<int32_t>::iterator start, const vector<int32_t>& rpn);
};
//==============================================================================
#ifdef DAGMC
class DAGCell : public Cell
{
public:
DAGCell();
bool contains(Position r, Direction u, int32_t on_surface) const;
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface, Particle* p) const;
BoundingBox bounding_box() const;
void to_hdf5(hid_t group_id) const;
moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance
int32_t dag_index_; //!< DagMC index of cell
};
#endif
//==============================================================================
//! Speeds up geometry searches by grouping cells in a search tree.
//
@ -350,8 +339,9 @@ struct CellInstanceHash {
void read_cells(pugi::xml_node node);
#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed);
class DAGUniverse;
#endif
} // namespace openmc

View file

@ -380,6 +380,15 @@ enum class RunMode {
// For non-accelerated regions on coarse mesh overlay
constexpr int CMFD_NOACCEL {-1};
//==============================================================================
// Geometry Constants
enum class GeometryType {
CSG,
DAG
};
} // namespace openmc
#endif // OPENMC_CONSTANTS_H

View file

@ -1,4 +1,3 @@
#ifndef OPENMC_DAGMC_H
#define OPENMC_DAGMC_H
@ -6,27 +5,146 @@ namespace openmc {
extern "C" const bool DAGMC_ENABLED;
}
#ifdef DAGMC
#include "DagMC.hpp"
// always include the XML interface header
#include "openmc/xml_interface.h"
#include "openmc/position.h"
//==============================================================================
// Functions that are always defined
//==============================================================================
namespace openmc {
namespace model {
extern moab::DagMC* DAG;
void read_dagmc_universes(pugi::xml_node node);
void check_dagmc_root_univ();
}
#ifdef DAGMC
#include "DagMC.hpp"
#include "openmc/cell.h"
#include "openmc/particle.h"
#include "openmc/position.h"
#include "openmc/surface.h"
class UWUW;
namespace openmc {
class DAGSurface : public Surface {
public:
DAGSurface(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx);
double evaluate(Position r) const override;
double distance(Position r, Direction u, bool coincident) const override;
Direction normal(Position r) const override;
Direction reflect(Position r, Direction u, Particle* p) const override;
inline void to_hdf5_inner(hid_t group_id) const override {};
// Accessor methods
const std::shared_ptr<moab::DagMC>& dagmc_ptr() const { return dagmc_ptr_; }
int32_t dag_index() const { return dag_index_; }
private:
std::shared_ptr<moab::DagMC> dagmc_ptr_; //!< Pointer to DagMC instance
int32_t dag_index_; //!< DagMC index of surface
};
class DAGCell : public Cell {
public:
DAGCell(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx);
bool contains(Position r, Direction u, int32_t on_surface) const override;
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface, Particle* p) const override;
BoundingBox bounding_box() const override;
void to_hdf5_inner(hid_t group_id) const override;
// Accessor methods
const std::shared_ptr<moab::DagMC>& dagmc_ptr() const { return dagmc_ptr_; }
int32_t dag_index() const { return dag_index_; }
private:
std::shared_ptr<moab::DagMC> dagmc_ptr_; //!< Pointer to DagMC instance
int32_t dag_index_; //!< DagMC index of cell
};
class DAGUniverse : public Universe {
public:
explicit DAGUniverse(pugi::xml_node node);
//! Create a new DAGMC universe
//! \param[in] filename Name of the DAGMC file
//! \param[in] auto_geom_ids Whether or not to automatically assign cell and surface IDs
//! \param[in] auto_mat_ids Whether or not to automatically assign material IDs
explicit DAGUniverse(const std::string& filename,
bool auto_geom_ids = false,
bool auto_mat_ids = false);
//! Initialize the DAGMC accel. data structures, indices, material assignments, etc.
void initialize();
//! Reads UWUW materials and returns an ID map
void read_uwuw_materials();
//! Indicates whether or not UWUW materials are present
//! \return True if UWUW materials are present, False if not
bool uses_uwuw() const;
//! Returns the index to the implicit complement's index in OpenMC for this DAGMC universe
int32_t implicit_complement_idx() const;
//! Transform UWUW materials into an OpenMC-readable XML format
//! \return A string representing a materials.xml file of the UWUW materials in this universe
std::string get_uwuw_materials_xml() const;
//! Writes the UWUW material file to XML (for debugging purposes)
void write_uwuw_materials_xml(const std::string& outfile = "uwuw_materials.xml") const;
//! Assign a material to a cell based
//! \param[in] mat_string The DAGMC material assignment string
//! \param[in] c The OpenMC cell to which the material is assigned
void legacy_assign_material(std::string mat_string,
std::unique_ptr<DAGCell>& c) const;
//! Generate a string representing the ranges of IDs present in the DAGMC model.
//! Contiguous chunks of IDs are represented as a range (i.e. 1-10). If there is
//! a single ID a chunk, it will be represented as a single number (i.e. 2, 4, 6, 8).
//! \param[in] dim Dimension of the entities
//! \return A string of the ID ranges for entities of dimension \p dim
std::string dagmc_ids_for_dim(int dim) const;
bool find_cell(Particle &p) const override;
void to_hdf5(hid_t universes_group) const override;
// Data Members
std::shared_ptr<moab::DagMC> dagmc_instance_; //!< DAGMC Instance for this universe
int32_t cell_idx_offset_; //!< An offset to the start of the cells in this universe in OpenMC's cell vector
int32_t surf_idx_offset_; //!< An offset to the start of the surfaces in this universe in OpenMC's surface vector
// Accessors
bool has_graveyard() const { return has_graveyard_; }
private:
std::string filename_; //!< Name of the DAGMC file used to create this universe
std::shared_ptr<UWUW> uwuw_; //!< Pointer to the UWUW instance for this universe
bool adjust_geometry_ids_; //!< Indicates whether or not to automatically generate new cell and surface IDs for the universe
bool adjust_material_ids_; //!< Indicates whether or not to automatically generate new material IDs for the universe
bool has_graveyard_; //!< Indicates if the DAGMC geometry has a "graveyard" volume
};
//==============================================================================
// Non-member functions
//==============================================================================
void load_dagmc_geometry();
void free_memory_dagmc();
void read_geometry_dagmc();
bool read_uwuw_materials(pugi::xml_document& doc);
bool get_uwuw_materials_xml(std::string& s);
int32_t next_cell(DAGUniverse* dag_univ, DAGCell* cur_cell, DAGSurface* surf_xed);
} // namespace openmc

View file

@ -366,6 +366,7 @@ public:
MOABMesh() = default;
MOABMesh(pugi::xml_node);
MOABMesh(const std::string& filename);
MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
// Overridden Methods
@ -412,6 +413,9 @@ private:
// Methods
//! Create the MOAB interface pointer
void create_interface();
//! Find all intersections with faces of the mesh.
//
//! \param[in] start Staring location
@ -503,7 +507,7 @@ private:
moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh
moab::EntityHandle tetset_; //!< EntitySet containing all tetrahedra
moab::EntityHandle kdtree_root_; //!< Root of the MOAB KDTree
unique_ptr<moab::Interface> mbi_; //!< MOAB instance
std::shared_ptr<moab::Interface> mbi_; //!< MOAB instance
unique_ptr<moab::AdaptiveKDTree> kdtree_; //!< MOAB KDTree instance
vector<moab::Matrix3> baryc_data_; //!< Barycentric data for tetrahedra
vector<std::string> tag_names_; //!< Names of score tags added to the mesh

View file

@ -29,7 +29,6 @@ extern bool check_overlaps; //!< check overlaps in geometry?
extern bool confidence_intervals; //!< use confidence intervals for results?
extern bool create_fission_neutrons; //!< create fission neutrons (fixed source)?
extern "C" bool cmfd_run; //!< is a CMFD run?
extern "C" bool dagmc; //!< indicator of DAGMC geometry
extern bool delayed_photon_scaling; //!< Scale fission photon yield to include delayed
extern "C" bool entropy_on; //!< calculate Shannon entropy?
extern bool event_based; //!< use event-based mode (instead of history-based)

View file

@ -8,7 +8,6 @@
#include "hdf5.h"
#include "pugixml.hpp"
#include "dagmc.h"
#include "openmc/boundary_condition.h"
#include "openmc/constants.h"
#include "openmc/memory.h" // for unique_ptr
@ -30,7 +29,7 @@ namespace model {
} // namespace model
//==============================================================================
//! Coordinates for an axis-aligned cuboid bounding a geometric object.
//! Coordinates for an axis-aligned cuboid that bounds a geometric object.
//==============================================================================
struct BoundingBox
@ -88,6 +87,7 @@ public:
int id_; //!< Unique ID
std::string name_; //!< User-defined name
std::shared_ptr<BoundaryCondition> bc_ {nullptr}; //!< Boundary condition
GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC)
bool surf_source_ {false}; //!< Activate source banking for the surface?
explicit Surface(pugi::xml_node surf_node);
@ -134,10 +134,13 @@ public:
//! Write all information needed to reconstruct the surface to an HDF5 group.
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const = 0;
void to_hdf5(hid_t group_id) const;
//! Get the BoundingBox for this surface.
virtual BoundingBox bounding_box(bool /*pos_side*/) const { return {}; }
protected:
virtual void to_hdf5_inner(hid_t group_id) const = 0;
};
class CSGSurface : public Surface
@ -146,33 +149,10 @@ public:
explicit CSGSurface(pugi::xml_node surf_node);
CSGSurface();
void to_hdf5(hid_t group_id) const;
protected:
virtual void to_hdf5_inner(hid_t group_id) const = 0;
};
//==============================================================================
//! A `Surface` representing a DAGMC-based surface in DAGMC.
//==============================================================================
#ifdef DAGMC
class DAGSurface : public Surface
{
public:
DAGSurface();
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
Direction reflect(Position r, Direction u, Particle* p) const;
void to_hdf5(hid_t group_id) const;
moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance
int32_t dag_index_; //!< DagMC index of surface
};
#endif
//==============================================================================
//! A plane perpendicular to the x-axis.
//

View file

@ -53,6 +53,9 @@ public:
void set_filters(gsl::span<Filter*> filters);
//! Given already-set filters, set the stride lengths
void set_strides();
int32_t strides(int i) const {return strides_[i];}
int32_t n_filter_bins() const {return n_filter_bins_;}

View file

@ -27,7 +27,7 @@ class Cell(IDManagerMixin):
automatically be assigned.
name : str, optional
Name of the cell. If not specified, the name is the empty string.
fill : openmc.Material or openmc.Universe or openmc.Lattice or None or iterable of openmc.Material, optional
fill : openmc.Material or openmc.UniverseBase or openmc.Lattice or None or iterable of openmc.Material, optional
Indicates what the region of space is filled with
region : openmc.Region, optional
Region of space that is assigned to the cell.
@ -38,7 +38,7 @@ class Cell(IDManagerMixin):
Unique identifier for the cell
name : str
Name of the cell
fill : openmc.Material or openmc.Universe or openmc.Lattice or None or iterable of openmc.Material
fill : openmc.Material or openmc.UniverseBase or openmc.Lattice or None or iterable of openmc.Material
Indicates what the region of space is filled with. If None, the cell is
treated as a void. An iterable of materials is used to fill repeated
instances of a cell with different materials.
@ -156,7 +156,7 @@ class Cell(IDManagerMixin):
def fill_type(self):
if isinstance(self.fill, openmc.Material):
return 'material'
elif isinstance(self.fill, openmc.Universe):
elif isinstance(self.fill, openmc.UniverseBase):
return 'universe'
elif isinstance(self.fill, openmc.Lattice):
return 'lattice'
@ -278,11 +278,10 @@ class Cell(IDManagerMixin):
cv.check_type('cell.fill[i]', f, openmc.Material)
elif not isinstance(fill, (openmc.Material, openmc.Lattice,
openmc.Universe)):
openmc.UniverseBase)):
msg = ('Unable to set Cell ID="{0}" to use a non-Material or '
'Universe fill "{1}"'.format(self._id, fill))
raise ValueError(msg)
self._fill = fill
# Info about atom content can now be invalid

View file

@ -15,7 +15,7 @@ from .cell import Cell
from .material import Material
from .mixin import IDManagerMixin
from .surface import Surface
from .universe import Universe
from .universe import UniverseBase
_FILTER_TYPES = (
@ -408,8 +408,8 @@ class UniverseFilter(WithIDFilter):
Parameters
----------
bins : openmc.Universe, int, or iterable thereof
The Universes to tally. Either openmc.Universe objects or their
bins : openmc.UniverseBase, int, or iterable thereof
The Universes to tally. Either openmc.UniverseBase objects or their
Integral ID numbers can be used.
filter_id : int
Unique identifier for the filter
@ -417,14 +417,14 @@ class UniverseFilter(WithIDFilter):
Attributes
----------
bins : Iterable of Integral
openmc.Universe IDs.
openmc.UniverseBase IDs.
id : int
Unique identifier for the filter
num_bins : Integral
The number of filter bins
"""
expected_type = Universe
expected_type = UniverseBase
class MaterialFilter(WithIDFilter):

View file

@ -14,13 +14,13 @@ class Geometry:
Parameters
----------
root : openmc.Universe or Iterable of openmc.Cell, optional
root : openmc.UniverseBase or Iterable of openmc.Cell, optional
Root universe which contains all others, or an iterable of cells that
should be used to create a root universe.
Attributes
----------
root_universe : openmc.Universe
root_universe : openmc.UniverseBase
Root universe which contains all others
bounding_box : 2-tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
@ -32,7 +32,7 @@ class Geometry:
self._root_universe = None
self._offsets = {}
if root is not None:
if isinstance(root, openmc.Universe):
if isinstance(root, openmc.UniverseBase):
self.root_universe = root
else:
univ = openmc.Universe()
@ -50,7 +50,7 @@ class Geometry:
@root_universe.setter
def root_universe(self, root_universe):
check_type('root universe', root_universe, openmc.Universe)
check_type('root universe', root_universe, openmc.UniverseBase)
self._root_universe = root_universe
def add_volume_information(self, volume_calc):
@ -159,6 +159,11 @@ class Geometry:
for s1, s2 in periodic.items():
surfaces[s1].periodic_surface = surfaces[s2]
# Add any DAGMC universes
for elem in root.findall('dagmc_universe'):
dag_univ = openmc.DAGMCUniverse.from_xml_element(elem)
universes[dag_univ.id] = dag_univ
# Dictionary that maps each universe to a list of cells/lattices that
# contain it (needed to determine which universe is the root)
child_of = defaultdict(list)

View file

@ -488,7 +488,7 @@ class RectLattice(Lattice):
@Lattice.universes.setter
def universes(self, universes):
cv.check_iterable_type('lattice universes', universes, openmc.Universe,
cv.check_iterable_type('lattice universes', universes, openmc.UniverseBase,
min_depth=2, max_depth=3)
self._universes = np.asarray(universes)

View file

@ -206,8 +206,7 @@ class Model:
d.mkdir(parents=True)
self.settings.export_to_xml(d)
if not self.settings.dagmc:
self.geometry.export_to_xml(d)
self.geometry.export_to_xml(d)
# If a materials collection was specified, export it. Otherwise, look
# for all materials in the geometry and use that to automatically build

View file

@ -44,8 +44,6 @@ class Settings:
weight assigned to particles that are not killed after Russian
roulette. Value of energy should be a float indicating energy in eV
below which particle type will be killed.
dagmc : bool
Indicate that a CAD-based DAGMC geometry will be used.
delayed_photon_scaling : bool
Indicate whether to scale the fission photon yield by (EGP + EGD)/EGP
where EGP is the energy release of prompt photons and EGD is the energy
@ -269,8 +267,6 @@ class Settings:
self._material_cell_offsets = None
self._log_grid_bins = None
self._dagmc = False
self._event_based = None
self._max_particles_in_flight = None
@ -434,10 +430,6 @@ class Settings:
def log_grid_bins(self):
return self._log_grid_bins
@property
def dagmc(self):
return self._dagmc
@property
def event_based(self):
return self._event_based
@ -633,11 +625,6 @@ class Settings:
cv.check_type('photon transport', photon_transport, bool)
self._photon_transport = photon_transport
@dagmc.setter
def dagmc(self, dagmc):
cv.check_type('dagmc geometry', dagmc, bool)
self._dagmc = dagmc
@ptables.setter
def ptables(self, ptables):
cv.check_type('probability tables', ptables, bool)
@ -1125,11 +1112,6 @@ class Settings:
elem = ET.SubElement(root, "log_grid_bins")
elem.text = str(self._log_grid_bins)
def _create_dagmc_subelement(self, root):
if self._dagmc:
elem = ET.SubElement(root, "dagmc")
elem.text = str(self._dagmc).lower()
def _eigenvalue_from_xml_element(self, root):
elem = root.find('eigenvalue')
if elem is not None:
@ -1407,11 +1389,6 @@ class Settings:
if text is not None:
self.log_grid_bins = int(text)
def _dagmc_from_xml_element(self, root):
text = get_text(root, 'dagmc')
if text is not None:
self.dagmc = text in ('true', '1')
def export_to_xml(self, path='settings.xml'):
"""Export simulation settings to an XML file.
@ -1465,7 +1442,6 @@ class Settings:
self._create_max_particles_in_flight_subelement(root_element)
self._create_material_cell_offsets_subelement(root_element)
self._create_log_grid_bins_subelement(root_element)
self._create_dagmc_subelement(root_element)
# Clean the indentation in the file to be user-readable
clean_indentation(root_element)
@ -1539,7 +1515,6 @@ class Settings:
settings._max_particles_in_flight_from_xml_element(root)
settings._material_cell_offsets_from_xml_element(root)
settings._log_grid_bins_from_xml_element(root)
settings._dagmc_from_xml_element(root)
# TODO: Get volume calculations

View file

@ -153,9 +153,7 @@ class StatePoint:
return self
def __exit__(self, *exc):
self._f.close()
if self._summary is not None:
self._summary._f.close()
self.close()
@property
def cmfd_on(self):
@ -490,6 +488,14 @@ class StatePoint:
for tally_id in self.tallies:
self.tallies[tally_id].sparse = self.sparse
def close(self):
"""Close the statepoint HDF5 file and the corresponding
summary HDF5 file if present.
"""
self._f.close()
if self._summary is not None:
self._summary._f.close()
def add_volume_information(self, volume_calc):
"""Add volume information to the geometry within the file

View file

@ -123,7 +123,9 @@ class Summary:
def _read_surfaces(self):
for group in self._f['geometry/surfaces'].values():
surface = openmc.Surface.from_hdf5(group)
self._fast_surfaces[surface.id] = surface
# surface may be None for DAGMC surfaces
if surface:
self._fast_surfaces[surface.id] = surface
def _read_cells(self):
@ -176,7 +178,11 @@ class Summary:
def _read_universes(self):
for group in self._f['geometry/universes'].values():
universe = openmc.Universe.from_hdf5(group, self._fast_cells)
geom_type = group.get('geom_type')
if geom_type and geom_type[()].decode() == 'dagmc':
universe = openmc.DAGMCUniverse.from_hdf5(group)
else:
universe = openmc.Universe.from_hdf5(group, self._fast_cells)
self._fast_universes[universe.id] = universe
def _read_lattices(self):

View file

@ -432,6 +432,7 @@ class Surface(IDManagerMixin, ABC):
kwargs = {}
kwargs['surface_id'] = int(elem.get('id'))
kwargs['boundary_type'] = elem.get('boundary', 'transmission')
kwargs['name'] = elem.get('name')
coeffs = [float(x) for x in elem.get('coeffs').split()]
kwargs.update(dict(zip(cls._coeff_keys, coeffs)))
@ -453,13 +454,19 @@ class Surface(IDManagerMixin, ABC):
"""
# If this is a DAGMC surface, do nothing for now
geom_type = group.get('geom_type')
if geom_type and geom_type[()].decode() == 'dagmc':
return
surface_id = int(group.name.split('/')[-1].lstrip('surface '))
name = group['name'][()].decode() if 'name' in group else ''
surf_type = group['type'][()].decode()
bc = group['boundary_type'][()].decode()
coeffs = group['coefficients'][...]
kwargs = {'boundary_type': bc, 'name': name, 'surface_id': surface_id}
surf_type = group['type'][()].decode()
cls = _SURFACE_CLASSES[surf_type]
return cls(*coeffs, **kwargs)

View file

@ -1,18 +1,154 @@
from abc import ABC, abstractmethod
from collections import OrderedDict
from collections.abc import Iterable
from copy import copy, deepcopy
from numbers import Real
import random
from xml.etree import ElementTree as ET
import numpy as np
import openmc
import openmc.checkvalue as cv
from ._xml import get_text
from .mixin import IDManagerMixin
from .plots import _SVG_COLORS
class Universe(IDManagerMixin):
class UniverseBase(ABC, IDManagerMixin):
"""A collection of cells that can be repeated.
Attributes
----------
id : int
Unique identifier of the universe
name : str
Name of the universe
"""
next_id = 1
used_ids = set()
def __init__(self, universe_id=None, name=''):
# Initialize Universe class attributes
self.id = universe_id
self.name = name
self._volume = None
self._atoms = {}
# Keys - Cell IDs
# Values - Cells
self._cells = OrderedDict()
def __repr__(self):
string = 'Universe\n'
string += '{: <16}=\t{}\n'.format('\tID', self._id)
string += '{: <16}=\t{}\n'.format('\tName', self._name)
return string
@property
def name(self):
return self._name
@property
def volume(self):
return self._volume
@name.setter
def name(self, name):
if name is not None:
cv.check_type('universe name', name, str)
self._name = name
else:
self._name = ''
@volume.setter
def volume(self, volume):
if volume is not None:
cv.check_type('universe volume', volume, Real)
self._volume = volume
def add_volume_information(self, volume_calc):
"""Add volume information to a universe.
Parameters
----------
volume_calc : openmc.VolumeCalculation
Results from a stochastic volume calculation
"""
if volume_calc.domain_type == 'universe':
if self.id in volume_calc.volumes:
self._volume = volume_calc.volumes[self.id].n
self._atoms = volume_calc.atoms[self.id]
else:
raise ValueError('No volume information found for this universe.')
else:
raise ValueError('No volume information found for this universe.')
@abstractmethod
def create_xml_subelement(self, xml_element, memo=None):
"""Add the universe xml representation to an incoming xml element
Parameters
----------
xml_element : xml.etree.ElementTree.Element
XML element to be added to
memo : set or None
A set of object id's representing geometry entities already
written to the xml_element. This parameter is used internally
and should not be specified by users.
Returns
-------
None
"""
def clone(self, clone_materials=True, clone_regions=True, memo=None):
"""Create a copy of this universe with a new unique ID, and clones
all cells within this universe.
Parameters
----------
clone_materials : bool
Whether to create separates copies of the materials filling cells
contained in this universe.
clone_regions : bool
Whether to create separates copies of the regions bounding cells
contained in this universe.
memo : dict or None
A nested dictionary of previously cloned objects. This parameter
is used internally and should not be specified by the user.
Returns
-------
clone : openmc.Universe
The clone of this universe
"""
if memo is None:
memo = {}
# If no memoize'd clone exists, instantiate one
if self not in memo:
clone = deepcopy(self)
clone.id = None
# Clone all cells for the universe clone
clone._cells = OrderedDict()
for cell in self._cells.values():
clone.add_cell(cell.clone(clone_materials, clone_regions,
memo))
# Memoize the clone
memo[self] = clone
return memo[self]
class Universe(UniverseBase):
"""A collection of cells that can be repeated.
Parameters
@ -44,42 +180,22 @@ class Universe(IDManagerMixin):
"""
next_id = 1
used_ids = set()
def __init__(self, universe_id=None, name='', cells=None):
# Initialize Cell class attributes
self.id = universe_id
self.name = name
self._volume = None
self._atoms = {}
# Keys - Cell IDs
# Values - Cells
self._cells = OrderedDict()
super().__init__(universe_id, name)
if cells is not None:
self.add_cells(cells)
def __repr__(self):
string = 'Universe\n'
string += '{: <16}=\t{}\n'.format('\tID', self._id)
string += '{: <16}=\t{}\n'.format('\tName', self._name)
string = super().__repr__()
string += '{: <16}=\t{}\n'.format('\tGeom', 'CSG')
string += '{: <16}=\t{}\n'.format('\tCells', list(self._cells.keys()))
return string
@property
def name(self):
return self._name
@property
def cells(self):
return self._cells
@property
def volume(self):
return self._volume
@property
def bounding_box(self):
regions = [c.region for c in self.cells.values()
@ -90,20 +206,6 @@ class Universe(IDManagerMixin):
# Infinite bounding box
return openmc.Intersection([]).bounding_box
@name.setter
def name(self, name):
if name is not None:
cv.check_type('universe name', name, str)
self._name = name
else:
self._name = ''
@volume.setter
def volume(self, volume):
if volume is not None:
cv.check_type('universe volume', volume, Real)
self._volume = volume
@classmethod
def from_hdf5(cls, group, cells):
"""Create universe from HDF5 group
@ -133,24 +235,6 @@ class Universe(IDManagerMixin):
return universe
def add_volume_information(self, volume_calc):
"""Add volume information to a universe.
Parameters
----------
volume_calc : openmc.VolumeCalculation
Results from a stochastic volume calculation
"""
if volume_calc.domain_type == 'universe':
if self.id in volume_calc.volumes:
self._volume = volume_calc.volumes[self.id].n
self._atoms = volume_calc.atoms[self.id]
else:
raise ValueError('No volume information found for this universe.')
else:
raise ValueError('No volume information found for this universe.')
def find(self, point):
"""Find cells/universes/lattices which contain a given point
@ -480,66 +564,7 @@ class Universe(IDManagerMixin):
return universes
def clone(self, clone_materials=True, clone_regions=True, memo=None):
"""Create a copy of this universe with a new unique ID, and clones
all cells within this universe.
Parameters
----------
clone_materials : bool
Whether to create separates copies of the materials filling cells
contained in this universe.
clone_regions : bool
Whether to create separates copies of the regions bounding cells
contained in this universe.
memo : dict or None
A nested dictionary of previously cloned objects. This parameter
is used internally and should not be specified by the user.
Returns
-------
clone : openmc.Universe
The clone of this universe
"""
if memo is None:
memo = {}
# If no nemoize'd clone exists, instantiate one
if self not in memo:
clone = deepcopy(self)
clone.id = None
# Clone all cells for the universe clone
clone._cells = OrderedDict()
for cell in self._cells.values():
clone.add_cell(cell.clone(clone_materials, clone_regions,
memo))
# Memoize the clone
memo[self] = clone
return memo[self]
def create_xml_subelement(self, xml_element, memo=None):
"""Add the universe xml representation to an incoming xml element
Parameters
----------
xml_element : xml.etree.ElementTree.Element
XML element to be added to
memo : set or None
A set of object id's representing geometry entities already
written to the xml_element. This parameter is used internally
and should not be specified by users.
Returns
-------
None
"""
# Iterate over all Cells
for cell in self._cells.values():
@ -600,3 +625,163 @@ class Universe(IDManagerMixin):
cell._num_instances += 1
if not instances_only:
cell._paths.append(cell_path)
class DAGMCUniverse(UniverseBase):
"""A reference to a DAGMC file to be used in the model.
Parameters
----------
filename : str
Path to the DAGMC file used to represent this universe.
universe_id : int, optional
Unique identifier of the universe. If not specified, an identifier will
automatically be assigned.
name : str, optional
Name of the universe. If not specified, the name is the empty string.
auto_geom_ids : bool
Set IDs automatically on initialization (True) or report overlaps
in ID space between CSG and DAGMC (False)
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps
in ID space between OpenMC and UWUW materials (False)
Attributes
----------
id : int
Unique identifier of the universe
name : str
Name of the universe
filename : str
Path to the DAGMC file used to represent this universe.
auto_geom_ids : bool
Set IDs automatically on initialization (True) or report overlaps
in ID space between CSG and DAGMC (False)
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps
in ID space between OpenMC and UWUW materials (False)
"""
def __init__(self,
filename,
universe_id=None,
name='',
auto_geom_ids=False,
auto_mat_ids=False):
super().__init__(universe_id, name)
# Initialize class attributes
self.filename = filename
self.auto_geom_ids = auto_geom_ids
self.auto_mat_ids = auto_mat_ids
def __repr__(self):
string = super().__repr__()
string += '{: <16}=\t{}\n'.format('\tGeom', 'DAGMC')
string += '{: <16}=\t{}\n'.format('\tFile', self.filename)
return string
@property
def filename(self):
return self._filename
@filename.setter
def filename(self, val):
cv.check_type('DAGMC filename', val, str)
self._filename = val
@property
def auto_geom_ids(self):
return self._auto_geom_ids
@auto_geom_ids.setter
def auto_geom_ids(self, val):
cv.check_type('DAGMC automatic geometry ids', val, bool)
self._auto_geom_ids = val
@property
def auto_mat_ids(self):
return self._auto_mat_ids
@auto_mat_ids.setter
def auto_mat_ids(self, val):
cv.check_type('DAGMC automatic material ids', val, bool)
self._auto_mat_ids = val
def get_all_cells(self, memo=None):
return OrderedDict()
def get_all_materials(self, memo=None):
return OrderedDict()
def create_xml_subelement(self, xml_element, memo=None):
if memo and self in memo:
return
if memo is not None:
memo.add(self)
# Set xml element values
dagmc_element = ET.Element('dagmc_universe')
dagmc_element.set('id', str(self.id))
if self.auto_geom_ids:
dagmc_element.set('auto_geom_ids', 'true')
if self.auto_mat_ids:
dagmc_element.set('auto_mat_ids', 'true')
dagmc_element.set('filename', self.filename)
xml_element.append(dagmc_element)
@classmethod
def from_hdf5(cls, group):
"""Create DAGMC universe from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
Returns
-------
openmc.DAGMCUniverse
DAGMCUniverse instance
"""
id = int(group.name.split('/')[-1].lstrip('universe '))
fname = group['filename'][()].decode()
name = group['name'][()].decode() if 'name' in group else None
out = cls(fname, universe_id=id, name=name)
out.auto_geom_ids = bool(group.attrs['auto_geom_ids'])
out.auto_mat_ids = bool(group.attrs['auto_mat_ids'])
return out
@classmethod
def from_xml_element(cls, elem):
"""Generate DAGMC universe from XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
`<dagmc_universe>` element
Returns
-------
openmc.DAGMCUniverse
DAGMCUniverse instance
"""
id = int(get_text(elem, 'id'))
fname = get_text(elem, 'filename')
out = cls(fname, universe_id=id)
name = get_text(elem, 'name')
if name is not None:
out.name = name
out.auto_geom_ids = bool(elem.get('auto_geom_ids'))
out.auto_mat_ids = bool(elem.get('auto_mat_ids'))
return out

View file

@ -96,6 +96,10 @@ void sort_fission_bank()
const auto& site = simulation::fission_bank[i];
int64_t offset = site.parent_id - 1 - simulation::work_index[mpi::rank];
int64_t idx = simulation::progeny_per_particle[offset] + site.progeny_id;
if (idx >= simulation::fission_bank.size()) {
fatal_error("Mismatch detected between sum of all particle progeny and "
"shared fission bank size.");
}
sorted_bank[idx] = site;
}

View file

@ -22,7 +22,6 @@
#include "openmc/material.h"
#include "openmc/nuclide.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -194,6 +193,9 @@ Universe::to_hdf5(hid_t universes_group) const
// Create a group for this universe.
auto group = create_group(universes_group, fmt::format("universe {}", id_));
// Write the geometry representation type.
write_string(group, "geom_type", "csg", false);
// Write the contained cells.
if (cells_.size() > 0) {
vector<int32_t> cell_ids;
@ -204,6 +206,31 @@ Universe::to_hdf5(hid_t universes_group) const
close_group(group);
}
bool
Universe::find_cell(Particle& p) const {
const auto& cells {
!partitioner_
? cells_
: partitioner_->get_cells(p.r_local(), p.u_local())
};
for (auto it = cells.begin(); it != cells.end(); it++) {
int32_t i_cell = *it;
int32_t i_univ = p.coord(p.n_coord()-1).universe;
if (model::cells[i_cell]->universe_ != i_univ) continue;
// Check if this cell contains the particle;
Position r {p.r_local()};
Direction u {p.u_local()};
auto surf = p.surface();
if (model::cells[i_cell]->contains(r, u, surf)) {
p.coord(p.n_coord()-1).cell = i_cell;
return true;
}
}
return false;
}
BoundingBox Universe::bounding_box() const {
BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
if (cells_.size() == 0) {
@ -321,14 +348,79 @@ void Cell::import_properties_hdf5(hid_t group)
close_group(cell_group);
}
void
Cell::to_hdf5(hid_t cell_group) const {
// Create a group for this cell.
auto group = create_group(cell_group, fmt::format("cell {}", id_));
if (!name_.empty()) {
write_string(group, "name", name_, false);
}
write_dataset(group, "universe", model::universes[universe_]->id_);
to_hdf5_inner(group);
// Write fill information.
if (type_ == Fill::MATERIAL) {
write_dataset(group, "fill_type", "material");
std::vector<int32_t> mat_ids;
for (auto i_mat : material_) {
if (i_mat != MATERIAL_VOID) {
mat_ids.push_back(model::materials[i_mat]->id_);
} else {
mat_ids.push_back(MATERIAL_VOID);
}
}
if (mat_ids.size() == 1) {
write_dataset(group, "material", mat_ids[0]);
} else {
write_dataset(group, "material", mat_ids);
}
std::vector<double> temps;
for (auto sqrtkT_val : sqrtkT_)
temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);
write_dataset(group, "temperature", temps);
} else if (type_ == Fill::UNIVERSE) {
write_dataset(group, "fill_type", "universe");
write_dataset(group, "fill", model::universes[fill_]->id_);
if (translation_ != Position(0, 0, 0)) {
write_dataset(group, "translation", translation_);
}
if (!rotation_.empty()) {
if (rotation_.size() == 12) {
std::array<double, 3> rot {rotation_[9], rotation_[10], rotation_[11]};
write_dataset(group, "rotation", rot);
} else {
write_dataset(group, "rotation", rotation_);
}
}
} else if (type_ == Fill::LATTICE) {
write_dataset(group, "fill_type", "lattice");
write_dataset(group, "lattice", model::lattices[fill_]->id_);
}
close_group(group);
}
//==============================================================================
// CSGCell implementation
//==============================================================================
CSGCell::CSGCell() {} // empty constructor
// default constructor
CSGCell::CSGCell() {
geom_type_ = GeometryType::CSG;
}
CSGCell::CSGCell(pugi::xml_node cell_node)
{
geom_type_ = GeometryType::CSG;
if (check_for_node(cell_node, "id")) {
id_ = std::stoi(get_node_value(cell_node, "id"));
} else {
@ -565,16 +657,10 @@ CSGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) cons
//==============================================================================
void
CSGCell::to_hdf5(hid_t cell_group) const
CSGCell::to_hdf5_inner(hid_t group_id) const
{
// Create a group for this cell.
auto group = create_group(cell_group, fmt::format("cell {}", id_));
if (!name_.empty()) {
write_string(group, "name", name_, false);
}
write_dataset(group, "universe", model::universes[universe_]->id_);
write_string(group_id, "geom_type", "csg", false);
// Write the region specification.
if (!region_.empty()) {
@ -595,52 +681,9 @@ CSGCell::to_hdf5(hid_t cell_group) const
region_spec << " " << ((token > 0) ? surf_id : -surf_id);
}
}
write_string(group, "region", region_spec.str(), false);
write_string(group_id, "region", region_spec.str(), false);
}
// Write fill information.
if (type_ == Fill::MATERIAL) {
write_dataset(group, "fill_type", "material");
vector<int32_t> mat_ids;
for (auto i_mat : material_) {
if (i_mat != MATERIAL_VOID) {
mat_ids.push_back(model::materials[i_mat]->id_);
} else {
mat_ids.push_back(MATERIAL_VOID);
}
}
if (mat_ids.size() == 1) {
write_dataset(group, "material", mat_ids[0]);
} else {
write_dataset(group, "material", mat_ids);
}
vector<double> temps;
for (auto sqrtkT_val : sqrtkT_)
temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);
write_dataset(group, "temperature", temps);
} else if (type_ == Fill::UNIVERSE) {
write_dataset(group, "fill_type", "universe");
write_dataset(group, "fill", model::universes[fill_]->id_);
if (translation_ != Position(0, 0, 0)) {
write_dataset(group, "translation", translation_);
}
if (!rotation_.empty()) {
if (rotation_.size() == 12) {
array<double, 3> rot {rotation_[9], rotation_[10], rotation_[11]};
write_dataset(group, "rotation", rot);
} else {
write_dataset(group, "rotation", rotation_);
}
}
} else if (type_ == Fill::LATTICE) {
write_dataset(group, "fill_type", "lattice");
write_dataset(group, "lattice", model::lattices[fill_]->id_);
}
close_group(group);
}
BoundingBox CSGCell::bounding_box_simple() const {
@ -815,70 +858,6 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
}
}
//==============================================================================
// DAGMC Cell implementation
//==============================================================================
#ifdef DAGMC
DAGCell::DAGCell() : Cell{} { simple_ = true; };
std::pair<double, int32_t>
DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
{
Expects(p);
// if we've changed direction or we're not on a surface,
// reset the history and update last direction
if (u != p->last_dir() || on_surface == 0) {
p->history().reset();
p->last_dir() = u;
}
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history());
MB_CHK_ERR_CONT(rval);
int surf_idx;
if (hit_surf != 0) {
surf_idx = dagmc_ptr_->index_by_handle(hit_surf);
} else {
// indicate that particle is lost
surf_idx = -1;
dist = INFINITY;
}
return {dist, surf_idx};
}
bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
int result = 0;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->point_in_volume(vol, pnt, result, dir);
MB_CHK_ERR_CONT(rval);
return result;
}
void DAGCell::to_hdf5(hid_t group_id) const { return; }
BoundingBox DAGCell::bounding_box() const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
double min[3], max[3];
rval = dagmc_ptr_->getobb(vol, min, max);
MB_CHK_ERR_CONT(rval);
return {min[0], max[0], min[1], max[1], min[2], max[2]};
}
#endif
//==============================================================================
// UniversePartitioner implementation
//==============================================================================
@ -1031,9 +1010,6 @@ void read_cells(pugi::xml_node node)
// Count the number of cells.
int n_cells = 0;
for (pugi::xml_node cell_node: node.children("cell")) {n_cells++;}
if (n_cells == 0) {
fatal_error("No cells found in geometry.xml!");
}
// Loop over XML cell elements and populate the array.
model::cells.reserve(n_cells);
@ -1052,6 +1028,8 @@ void read_cells(pugi::xml_node node)
}
}
read_dagmc_universes(node);
// Populate the Universe vector and map.
for (int i = 0; i < model::cells.size(); i++) {
int32_t uid = model::cells[i]->universe_;
@ -1071,6 +1049,10 @@ void read_cells(pugi::xml_node node)
if (settings::check_overlaps) {
model::overlap_check_count.resize(model::cells.size(), 0);
}
if (model::cells.size() == 0) {
fatal_error("No cells were found in the geometry.xml file");
}
}
//==============================================================================
@ -1340,21 +1322,6 @@ openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end)
return 0;
}
#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed)
{
moab::EntityHandle surf =
surf_xed->dagmc_ptr_->entity_by_index(2, surf_xed->dag_index_);
moab::EntityHandle vol =
cur_cell->dagmc_ptr_->entity_by_index(3, cur_cell->dag_index_);
moab::EntityHandle new_vol;
cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol);
return cur_cell->dagmc_ptr_->index_by_handle(new_vol);
}
#endif
extern "C" int cells_size() { return model::cells.size(); }
} // namespace openmc

View file

@ -3,9 +3,6 @@
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/container_util.h"
#ifdef DAGMC
#include "openmc/dagmc.h"
#endif
#include "openmc/error.h"
#include "openmc/geometry_aux.h"
#include "openmc/file_utils.h"
@ -96,27 +93,12 @@ void read_cross_sections_xml()
{
pugi::xml_document doc;
std::string filename = settings::path_input + "materials.xml";
#ifdef DAGMC
std::string s;
bool found_uwuw_mats = false;
if (settings::dagmc) {
found_uwuw_mats = get_uwuw_materials_xml(s);
}
if (found_uwuw_mats) {
// if we found uwuw materials, load those
doc.load_file(s.c_str());
} else {
#endif
// Check if materials.xml exists
if (!file_exists(filename)) {
fatal_error("Material XML file '" + filename + "' does not exist.");
}
// Parse materials.xml file
doc.load_file(filename.c_str());
#ifdef DAGMC
}
#endif
auto root = doc.document_element();

View file

@ -1,16 +1,15 @@
#include "openmc/dagmc.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/container_util.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/material.h"
#include "openmc/string_utils.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
#ifdef DAGMC
#include "uwuw.hpp"
@ -37,189 +36,135 @@ const bool DAGMC_ENABLED = false;
namespace openmc {
const std::string DAGMC_FILENAME = "dagmc.h5m";
//==============================================================================
// DAGMC Universe implementation
//==============================================================================
namespace model {
moab::DagMC* DAG;
} // namespace model
std::string dagmc_file() {
std::string filename = settings::path_input + DAGMC_FILENAME;
if (!file_exists(filename)) {
fatal_error("Geometry DAGMC file '" + filename + "' does not exist!");
DAGUniverse::DAGUniverse(pugi::xml_node node) {
if (check_for_node(node, "id")) {
id_ = std::stoi(get_node_value(node, "id"));
} else {
fatal_error("Must specify the id of the DAGMC universe");
}
return filename;
if (check_for_node(node, "filename")) {
filename_ = get_node_value(node, "filename");
} else {
fatal_error("Must specify a file for the DAGMC universe");
}
adjust_geometry_ids_ = false;
if (check_for_node(node, "auto_geom_ids")) {
adjust_geometry_ids_ = get_node_value_bool(node, "auto_geom_ids");
}
adjust_material_ids_ = false;
if (check_for_node(node, "auto_mat_ids")) {
adjust_material_ids_ = get_node_value_bool(node, "auto_mat_ids");
}
initialize();
}
bool get_uwuw_materials_xml(std::string& s) {
std::string filename = dagmc_file();
UWUW uwuw(filename.c_str());
std::stringstream ss;
bool uwuw_mats_present = false;
if (uwuw.material_library.size() != 0) {
uwuw_mats_present = true;
// write header
ss << "<?xml version=\"1.0\"?>\n";
ss << "<materials>\n";
const auto& mat_lib = uwuw.material_library;
// write materials
for (auto mat : mat_lib) { ss << mat.second->openmc("atom"); }
// write footer
ss << "</materials>";
s = ss.str();
DAGUniverse::DAGUniverse(const std::string& filename,
bool auto_geom_ids,
bool auto_mat_ids)
: filename_(filename), adjust_geometry_ids_(auto_geom_ids), adjust_material_ids_(auto_mat_ids) {
// determine the next universe id
int32_t next_univ_id = 0;
for (const auto& u : model::universes) {
if (u->id_ > next_univ_id) next_univ_id = u->id_;
}
next_univ_id++;
return uwuw_mats_present;
// set the universe id
id_ = next_univ_id;
initialize();
}
bool read_uwuw_materials(pugi::xml_document& doc) {
std::string s;
bool found_uwuw_mats = get_uwuw_materials_xml(s);
if (found_uwuw_mats) {
pugi::xml_parse_result result = doc.load_string(s.c_str());
if (!result) {
throw std::runtime_error{"Error reading UWUW materials"};
}
}
return found_uwuw_mats;
}
void
DAGUniverse::initialize() {
geom_type() = GeometryType::DAG;
bool write_uwuw_materials_xml() {
std::string s;
bool found_uwuw_mats = get_uwuw_materials_xml(s);
// if there is a material library in the file
if (found_uwuw_mats) {
// write a material.xml file
std::ofstream mats_xml("materials.xml");
mats_xml << s;
mats_xml.close();
// determine the next cell id
int32_t next_cell_id = 0;
for (const auto& c : model::cells) {
if (c->id_ > next_cell_id) next_cell_id = c->id_;
}
cell_idx_offset_ = model::cells.size();
next_cell_id++;
return found_uwuw_mats;
}
void legacy_assign_material(std::string mat_string, DAGCell* c)
{
bool mat_found_by_name = false;
// attempt to find a material with a matching name
to_lower(mat_string);
for (const auto& m : model::materials) {
std::string m_name = m->name();
to_lower(m_name);
if (mat_string == m_name) {
// assign the material with that name
if (!mat_found_by_name) {
mat_found_by_name = true;
c->material_.push_back(m->id_);
// report error if more than one material is found
} else {
fatal_error(fmt::format(
"More than one material found with name {}. Please ensure materials "
"have unique names if using this property to assign materials.",
mat_string));
}
}
// determine the next surface id
int32_t next_surf_id = 0;
for (const auto& s : model::surfaces) {
if (s->id_ > next_surf_id) next_surf_id = s->id_;
}
surf_idx_offset_ = model::surfaces.size();
next_surf_id++;
// if no material was set using a name, assign by id
if (!mat_found_by_name) {
try {
auto id = std::stoi(mat_string);
c->material_.emplace_back(id);
} catch (const std::invalid_argument&) {
fatal_error(fmt::format(
"No material {} found for volume (cell) {}", mat_string, c->id_));
}
}
if (settings::verbosity >= 10) {
const auto& m = model::materials[model::material_map.at(c->material_[0])];
std::stringstream msg;
msg << "DAGMC material " << mat_string << " was assigned";
if (mat_found_by_name) {
msg << " using material name: " << m->name_;
} else {
msg << " using material id: " << m->id_;
}
write_message(msg.str(), 10);
}
}
void load_dagmc_geometry()
{
if (!model::DAG) {
model::DAG = new moab::DagMC();
}
// create a new DAGMC instance
dagmc_instance_ = std::make_shared<moab::DagMC>();
// --- Materials ---
// create uwuw instance
auto filename = dagmc_file();
UWUW uwuw(filename.c_str());
// read any UWUW materials from the file
read_uwuw_materials();
// check for uwuw material definitions
bool using_uwuw = !uwuw.material_library.empty();
bool using_uwuw = uses_uwuw();
// notify user if UWUW materials are going to be used
if (using_uwuw) {
write_message("Found UWUW Materials in the DAGMC geometry file.", 6);
}
int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC runs
// load the DAGMC geometry
moab::ErrorCode rval = model::DAG->load_file(filename.c_str());
filename_ = settings::path_input + filename_;
if (!file_exists(filename_)) {
fatal_error("Geometry DAGMC file '" + filename_ + "' does not exist!");
}
moab::ErrorCode rval = dagmc_instance_->load_file(filename_.c_str());
MB_CHK_ERR_CONT(rval);
// initialize acceleration data structures
rval = model::DAG->init_OBBTree();
rval = dagmc_instance_->init_OBBTree();
MB_CHK_ERR_CONT(rval);
// parse model metadata
dagmcMetaData DMD(model::DAG, false, false);
dagmcMetaData DMD(dagmc_instance_.get(), false, false);
DMD.load_property_data();
vector<std::string> keywords {"temp"};
std::vector<std::string> keywords {"temp"};
std::map<std::string, std::string> dum;
std::string delimiters = ":/";
rval = model::DAG->parse_properties(keywords, dum, delimiters.c_str());
rval = dagmc_instance_->parse_properties(keywords, dum, delimiters.c_str());
MB_CHK_ERR_CONT(rval);
// --- Cells (Volumes) ---
// initialize cell objects
int n_cells = model::DAG->num_entities(3);
int n_cells = dagmc_instance_->num_entities(3);
moab::EntityHandle graveyard = 0;
for (int i = 0; i < n_cells; i++) {
moab::EntityHandle vol_handle = model::DAG->entity_by_index(3, i+1);
moab::EntityHandle vol_handle = dagmc_instance_->entity_by_index(3, i + 1);
// set cell ids using global IDs
DAGCell* c = new DAGCell();
c->dag_index_ = i+1;
c->id_ = model::DAG->id_by_index(3, c->dag_index_);
c->dagmc_ptr_ = model::DAG;
c->universe_ = dagmc_univ_id; // set to zero for now
auto c = std::make_unique<DAGCell>(dagmc_instance_, i + 1);
c->id_ = adjust_geometry_ids_ ? next_cell_id++ : dagmc_instance_->id_by_index(3, c->dag_index());
c->universe_ = this->id_;
c->fill_ = C_NONE; // no fill, single universe
model::cells.emplace_back(c);
model::cell_map[c->id_] = i;
// Populate the Universe vector and dict
auto it = model::universe_map.find(dagmc_univ_id);
if (it == model::universe_map.end()) {
model::universes.push_back(make_unique<Universe>());
model::universes.back()->id_ = dagmc_univ_id;
model::universes.back()->cells_.push_back(i);
model::universe_map[dagmc_univ_id] = model::universes.size() - 1;
auto in_map = model::cell_map.find(c->id_);
if (in_map == model::cell_map.end()) {
model::cell_map[c->id_] = model::cells.size();
} else {
model::universes[it->second]->cells_.push_back(i);
warning(fmt::format("DAGMC Cell IDs: {}", dagmc_ids_for_dim(3)));
fatal_error(fmt::format("Cell ID {} exists in both DAGMC Universe {} "
"and the CSG geometry.", c->id_, this->id_));
}
// MATERIALS
// --- Materials ---
// determine volume material assignment
std::string mat_str = DMD.get_volume_property("material", vol_handle);
@ -241,12 +186,12 @@ void load_dagmc_geometry()
if (using_uwuw) {
// lookup material in uwuw if present
std::string uwuw_mat = DMD.volume_material_property_data_eh[vol_handle];
if (uwuw.material_library.count(uwuw_mat) != 0) {
if (uwuw_->material_library.count(uwuw_mat) != 0) {
// Note: material numbers are set by UWUW
int mat_number = uwuw.material_library.get_material(uwuw_mat).metadata["mat_number"].asInt();
int mat_number = uwuw_->material_library.get_material(uwuw_mat).metadata["mat_number"].asInt();
c->material_.push_back(mat_number);
} else {
fatal_error(fmt::format("Material with value {} not found in the "
fatal_error(fmt::format("Material with value '{}' not found in the "
"UWUW material library", mat_str));
}
} else {
@ -258,19 +203,25 @@ void load_dagmc_geometry()
std::string temp_value;
// no temperature if void
if (c->material_[0] == MATERIAL_VOID) continue;
if (c->material_[0] == MATERIAL_VOID) {
model::cells.emplace_back(std::move(c));
continue;
}
// assign cell temperature
const auto& mat = model::materials[model::material_map.at(c->material_[0])];
if (model::DAG->has_prop(vol_handle, "temp")) {
rval = model::DAG->prop_value(vol_handle, "temp", temp_value);
if (dagmc_instance_->has_prop(vol_handle, "temp")) {
rval = dagmc_instance_->prop_value(vol_handle, "temp", temp_value);
MB_CHK_ERR_CONT(rval);
double temp = std::stod(temp_value);
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
} else {
} else if (mat->temperature() > 0.0) {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature()));
} else {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * settings::temperature_default));
}
model::cells.emplace_back(std::move(c));
}
// allocate the cell overlap count if necessary
@ -278,40 +229,28 @@ void load_dagmc_geometry()
model::overlap_check_count.resize(model::cells.size(), 0);
}
if (!graveyard) {
warning("No graveyard volume found in the DagMC model."
"This may result in lost particles and rapid simulation failure.");
}
has_graveyard_ = graveyard;
// --- Surfaces ---
// initialize surface objects
int n_surfaces = model::DAG->num_entities(2);
int n_surfaces = dagmc_instance_->num_entities(2);
for (int i = 0; i < n_surfaces; i++) {
moab::EntityHandle surf_handle = model::DAG->entity_by_index(2, i+1);
moab::EntityHandle surf_handle = dagmc_instance_->entity_by_index(2, i+1);
// set cell ids using global IDs
DAGSurface* s = new DAGSurface();
s->dag_index_ = i+1;
s->id_ = model::DAG->id_by_index(2, s->dag_index_);
s->dagmc_ptr_ = model::DAG;
if (contains(settings::source_write_surf_id, s->id_)) {
s->surf_source_ = true;
}
auto s = std::make_unique<DAGSurface>(dagmc_instance_, i+1);
s->id_ = adjust_geometry_ids_ ? next_surf_id++ : dagmc_instance_->id_by_index(2, i+1);
// set BCs
std::string bc_value = DMD.get_surface_property("boundary", surf_handle);
to_lower(bc_value);
if (bc_value.empty() || bc_value == "transmit" || bc_value == "transmission") {
// Leave the bc_ a nullptr
// set to transmission by default (nullptr)
} else if (bc_value == "vacuum") {
s->bc_ = std::make_shared<VacuumBC>();
} else if (bc_value == "reflective" || bc_value == "reflect" || bc_value == "reflecting") {
s->bc_ = std::make_shared<ReflectiveBC>();
} else if (bc_value == "white") {
fatal_error("White boundary condition not supported in DAGMC.");
} else if (bc_value == "periodic") {
fatal_error("Periodic boundary condition not supported in DAGMC.");
} else {
@ -321,7 +260,7 @@ void load_dagmc_geometry()
// graveyard check
moab::Range parent_vols;
rval = model::DAG->moab_instance()->get_parent_meshsets(surf_handle, parent_vols);
rval = dagmc_instance_->moab_instance()->get_parent_meshsets(surf_handle, parent_vols);
MB_CHK_ERR_CONT(rval);
// if this surface belongs to the graveyard
@ -331,26 +270,417 @@ void load_dagmc_geometry()
}
// add to global array and map
model::surfaces.emplace_back(s);
model::surface_map[s->id_] = i;
auto in_map = model::surface_map.find(s->id_);
if (in_map == model::surface_map.end()) {
model::surface_map[s->id_] = model::surfaces.size();
} else {
warning(fmt::format("DAGMC Surface IDs: {}", dagmc_ids_for_dim(2)));
fatal_error(fmt::format("Surface ID {} exists in both Universe {} "
"and the CSG geometry.", s->id_, this->id_));
}
model::surfaces.emplace_back(std::move(s));
} // end surface loop
}
std::string
DAGUniverse::dagmc_ids_for_dim(int dim) const
{
// generate a vector of ids
std::vector<int> id_vec;
int n_ents = dagmc_instance_->num_entities(dim);
for (int i = 1; i <= n_ents; i++) {
id_vec.push_back(dagmc_instance_->id_by_index(dim, i));
}
return;
// sort the vector of ids
std::sort(id_vec.begin(), id_vec.end());
// generate a string representation of the ID range(s)
std::stringstream out;
int i = 0;
int start_id = id_vec[0]; // initialize with first ID
int stop_id;
// loop over all cells in the universe
while (i < n_ents) {
stop_id = id_vec[i];
// if the next ID is not in this contiguous set of IDS,
// figure out how to write the string representing this set
if (id_vec[i + 1] > stop_id + 1) {
if (start_id != stop_id) {
// there are several IDs in a row, print condensed version (i.e. 1-10, 12-20)
out << start_id << "-" << stop_id;
} else {
// only one ID in this contiguous block (i.e. 3, 5, 7, 9)
out << start_id;
}
// insert a comma as long as we aren't in the last ID set
if (i < n_ents - 1) { out << ", "; }
// if we are at the end of a set, set the start ID to the first value
// in the next set.
start_id = id_vec[++i];
}
i++;
}
return out.str();
}
void read_geometry_dagmc()
int32_t
DAGUniverse::implicit_complement_idx() const {
moab::EntityHandle ic;
moab::ErrorCode rval = dagmc_instance_->geom_tool()->get_implicit_complement(ic);
MB_CHK_SET_ERR_CONT(rval, "Failed to get implicit complement");
// off-by-one: DAGMC indices start at one
return cell_idx_offset_ + dagmc_instance_->index_by_handle(ic) - 1;
}
bool
DAGUniverse::find_cell(Particle &p) const {
// if the particle isn't in any of the other DagMC
// cells, place it in the implicit complement
bool found = Universe::find_cell(p);
if (!found && model::universe_map[this->id_] != model::root_universe) {
p.coord(p.n_coord() - 1).cell = implicit_complement_idx();
found = true;
}
return found;
}
void
DAGUniverse::to_hdf5(hid_t universes_group) const {
// Create a group for this universe.
auto group = create_group(universes_group, fmt::format("universe {}", id_));
// Write the geometry representation type.
write_string(group, "geom_type", "dagmc", false);
// Write other properties of the DAGMC Universe
write_string(group, "filename", filename_, false);
write_attribute(group, "auto_geom_ids", static_cast<int>(adjust_geometry_ids_));
write_attribute(group, "auto_mat_ids", static_cast<int>(adjust_material_ids_));
close_group(group);
}
bool
DAGUniverse::uses_uwuw() const
{
write_message("Reading DAGMC geometry...", 5);
load_dagmc_geometry();
model::root_universe = find_root_universe();
return !uwuw_->material_library.empty();
}
void free_memory_dagmc()
std::string
DAGUniverse::get_uwuw_materials_xml() const
{
delete model::DAG;
if (!uses_uwuw()) {
throw std::runtime_error("This DAGMC Universe does not use UWUW materials");
}
std::stringstream ss;
// write header
ss << "<?xml version=\"1.0\"?>\n";
ss << "<materials>\n";
const auto& mat_lib = uwuw_->material_library;
// write materials
for (auto mat : mat_lib) { ss << mat.second->openmc("atom"); }
// write footer
ss << "</materials>";
return ss.str();
}
void
DAGUniverse::write_uwuw_materials_xml(const std::string& outfile) const
{
if (!uses_uwuw()) {
throw std::runtime_error("This DAGMC universe does not use UWUW materials.");
}
std::string xml_str = get_uwuw_materials_xml();
// if there is a material library in the file
std::ofstream mats_xml(outfile);
mats_xml << xml_str;
mats_xml.close();
}
void
DAGUniverse::legacy_assign_material(std::string mat_string,
std::unique_ptr<DAGCell>& c) const
{
bool mat_found_by_name = false;
// attempt to find a material with a matching name
to_lower(mat_string);
for (const auto& m : model::materials) {
std::string m_name = m->name();
to_lower(m_name);
if (mat_string == m_name) {
// assign the material with that name
if (!mat_found_by_name) {
mat_found_by_name = true;
c->material_.push_back(m->id_);
// report error if more than one material is found
} else {
fatal_error(fmt::format(
"More than one material found with name '{}'. Please ensure materials "
"have unique names if using this property to assign materials.",
mat_string));
}
}
}
// if no material was set using a name, assign by id
if (!mat_found_by_name) {
try {
auto id = std::stoi(mat_string);
c->material_.emplace_back(id);
} catch (const std::invalid_argument&) {
fatal_error(fmt::format(
"No material '{}' found for volume (cell) {}", mat_string, c->id_));
}
}
if (settings::verbosity >= 10) {
const auto& m = model::materials[model::material_map.at(c->material_[0])];
std::stringstream msg;
msg << "DAGMC material " << mat_string << " was assigned";
if (mat_found_by_name) {
msg << " using material name: " << m->name_;
} else {
msg << " using material id: " << m->id_;
}
write_message(msg.str(), 10);
}
}
void
DAGUniverse::read_uwuw_materials() {
int32_t next_material_id = 0;
for (const auto& m : model::materials) {
next_material_id = std::max(m->id_, next_material_id);
}
next_material_id++;
uwuw_ = std::make_shared<UWUW>(filename_.c_str());
const auto& mat_lib = uwuw_->material_library;
if (mat_lib.size() == 0) return;
// if we're using automatic IDs, update the UWUW material metadata
if (adjust_material_ids_) {
for (auto& mat : uwuw_->material_library) {
mat.second->metadata["mat_number"] = next_material_id++;
}
}
std::stringstream ss;
ss << "<?xml version=\"1.0\"?>\n";
ss << "<materials>\n";
for (auto mat : mat_lib) { ss << mat.second->openmc("atom"); }
ss << "</materials>";
std::string mat_xml_string = ss.str();
// create a pugi XML document from this string
pugi::xml_document doc;
auto result = doc.load_string(mat_xml_string.c_str());
if (!result) {
fatal_error("Error processing XML created using DAGMC UWUW materials.");
}
pugi::xml_node root = doc.document_element();
for (pugi::xml_node material_node : root.children("material")) {
model::materials.push_back(std::make_unique<Material>(material_node));
}
}
//==============================================================================
// DAGMC Cell implementation
//==============================================================================
DAGCell::DAGCell(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
: Cell{}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx) {
geom_type_ = GeometryType::DAG;
simple_ = true;
};
std::pair<double, int32_t>
DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
{
Expects(p);
// if we've changed direction or we're not on a surface,
// reset the history and update last direction
if (u != p->last_dir()) { p->last_dir() = u; p->history().reset(); }
if (on_surface == 0) { p->history().reset(); }
const auto& univ = model::universes[p->coord(p->n_coord() - 1).universe];
DAGUniverse* dag_univ = static_cast<DAGUniverse*>(univ.get());
if (!dag_univ) fatal_error("DAGMC call made for particle in a non-DAGMC universe");
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history());
MB_CHK_ERR_CONT(rval);
int surf_idx;
if (hit_surf != 0) {
surf_idx = dag_univ->surf_idx_offset_ + dagmc_ptr_->index_by_handle(hit_surf);
} else {
// indicate that particle is lost
surf_idx = -1;
dist = INFINITY;
if (!dagmc_ptr_->is_implicit_complement(vol) || model::universe_map[dag_univ->id_] == model::root_universe) {
p->mark_as_lost(fmt::format("No intersection found with DAGMC cell {}", id_));
}
}
return {dist, surf_idx};
}
bool
DAGCell::contains(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
int result = 0;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->point_in_volume(vol, pnt, result, dir);
MB_CHK_ERR_CONT(rval);
return result;
}
void
DAGCell::to_hdf5_inner(hid_t group_id) const {
write_string(group_id, "geom_type", "dagmc", false);
}
BoundingBox
DAGCell::bounding_box() const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
double min[3], max[3];
rval = dagmc_ptr_->getobb(vol, min, max);
MB_CHK_ERR_CONT(rval);
return {min[0], max[0], min[1], max[1], min[2], max[2]};
}
//==============================================================================
// DAGSurface implementation
//==============================================================================
DAGSurface::DAGSurface(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
: Surface{}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx)
{
geom_type_ = GeometryType::DAG;
} // empty constructor
double
DAGSurface::evaluate(Position r) const
{
return 0.0;
}
double
DAGSurface::distance(Position r, Direction u, bool coincident) const
{
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(surf, pnt, dir, hit_surf, dist, NULL, 0, 0);
MB_CHK_ERR_CONT(rval);
if (dist < 0.0) dist = INFTY;
return dist;
}
Direction
DAGSurface::normal(Position r) const
{
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
double pnt[3] = {r.x, r.y, r.z};
double dir[3];
rval = dagmc_ptr_->get_angle(surf, pnt, dir);
MB_CHK_ERR_CONT(rval);
return dir;
}
Direction
DAGSurface::reflect(Position r, Direction u, Particle* p) const
{
Expects(p);
p->history().reset_to_last_intersection();
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
double pnt[3] = {r.x, r.y, r.z};
double dir[3];
rval = dagmc_ptr_->get_angle(surf, pnt, dir, &p->history());
MB_CHK_ERR_CONT(rval);
p->last_dir() = u.reflect(dir);
return p->last_dir();
}
//==============================================================================
// Non-member functions
//==============================================================================
void read_dagmc_universes(pugi::xml_node node) {
for (pugi::xml_node dag_node : node.children("dagmc_universe")) {
model::universes.push_back(std::make_unique<DAGUniverse>(dag_node));
model::universe_map[model::universes.back()->id_] = model::universes.size() - 1;
}
}
void check_dagmc_root_univ() {
const auto& ru = model::universes[model::root_universe];
if (ru->geom_type() == GeometryType::DAG) {
// if the root universe contains DAGMC geometry, warn the user
// if it does not contain a graveyard volume
auto dag_univ = dynamic_cast<DAGUniverse*>(ru.get());
if (dag_univ && !dag_univ->has_graveyard()) {
warning("No graveyard volume found in the DagMC model. "
"This may result in lost particles and rapid simulation failure.");
}
}
}
int32_t next_cell(DAGUniverse* dag_univ, DAGCell* cur_cell, DAGSurface* surf_xed)
{
moab::EntityHandle surf =
surf_xed->dagmc_ptr()->entity_by_index(2, surf_xed->dag_index());
moab::EntityHandle vol =
cur_cell->dagmc_ptr()->entity_by_index(3, cur_cell->dag_index());
moab::EntityHandle new_vol;
cur_cell->dagmc_ptr()->next_vol(surf, vol, new_vol);
return cur_cell->dagmc_ptr()->index_by_handle(new_vol) + dag_univ->cell_idx_offset_;
}
}
#endif
} // namespace openmc
#else
namespace openmc {
void read_dagmc_universes(pugi::xml_node node) {};
void check_dagmc_root_univ() {};
} // namespace openmc
#endif // DAGMC

View file

@ -51,9 +51,6 @@ void free_memory()
if (settings::event_based) {
free_event_queues();
}
#ifdef DAGMC
free_memory_dagmc();
#endif
}
}
@ -74,7 +71,6 @@ int openmc_finalize()
settings::confidence_intervals = false;
settings::create_fission_neutrons = true;
settings::electron_treatment = ElectronTreatment::LED;
settings::dagmc = false;
settings::delayed_photon_scaling = true;
settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0};
settings::entropy_on = false;

View file

@ -112,35 +112,13 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list)
// that.
if (i_cell == C_NONE) {
int i_universe = p.coord(p.n_coord() - 1).universe;
const auto& univ {*model::universes[i_universe]};
const auto& cells {
!univ.partitioner_
? model::universes[i_universe]->cells_
: univ.partitioner_->get_cells(p.r_local(), p.u_local())
};
for (auto it = cells.cbegin(); it != cells.cend(); it++) {
i_cell = *it;
// Make sure the search cell is in the same universe.
int i_universe = p.coord(p.n_coord() - 1).universe;
if (model::cells[i_cell]->universe_ != i_universe) continue;
// Check if this cell contains the particle.
Position r {p.r_local()};
Direction u {p.u_local()};
auto surf = p.surface();
if (model::cells[i_cell]->contains(r, u, surf)) {
p.coord(p.n_coord() - 1).cell = i_cell;
found = true;
break;
}
}
}
if (!found) {
return found;
const auto& univ {model::universes[i_universe]};
found = univ->find_cell(p);
}
if (!found) { return found; }
i_cell = p.coord(p.n_coord() - 1).cell;
// Announce the cell that the particle is entering.
if (found && (settings::verbosity >= 10 || p.trace())) {
auto msg = fmt::format(" Entering cell {}", model::cells[i_cell]->id_);

View file

@ -42,13 +42,6 @@ void update_universe_cell_count(int32_t a, int32_t b) {
void read_geometry_xml()
{
#ifdef DAGMC
if (settings::dagmc) {
read_geometry_dagmc();
return;
}
#endif
// Display output message
write_message("Reading geometry XML file...", 5);
@ -73,8 +66,25 @@ void read_geometry_xml()
read_cells(root);
read_lattices(root);
// Check to make sure a boundary condition was applied to at least one
// surface
bool boundary_exists = false;
for (const auto& surf : model::surfaces) {
if (surf->bc_) {
boundary_exists = true;
break;
}
}
if (settings::run_mode != RunMode::PLOTTING && !boundary_exists) {
fatal_error("No boundary conditions were applied to any surfaces!");
}
// Allocate universes, universe cell arrays, and assign base universe
model::root_universe = find_root_universe();
// if the root universe is DAGMC geometry, make sure the model is well-formed
check_dagmc_root_univ();
}
//==============================================================================

View file

@ -14,7 +14,6 @@
#include "openmc/capi.h"
#include "openmc/cross_sections.h"
#include "openmc/container_util.h"
#include "openmc/dagmc.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/hdf5_interface.h"
@ -1240,24 +1239,14 @@ void read_materials_xml()
pugi::xml_document doc;
bool using_dagmc_mats = false;
#ifdef DAGMC
if (settings::dagmc) {
using_dagmc_mats = read_uwuw_materials(doc);
// Check if materials.xml exists
std::string filename = settings::path_input + "materials.xml";
if (!file_exists(filename)) {
fatal_error("Material XML file '" + filename + "' does not exist!");
}
#endif
if (!using_dagmc_mats) {
// Check if materials.xml exists
std::string filename = settings::path_input + "materials.xml";
if (!file_exists(filename)) {
fatal_error("Material XML file '" + filename + "' does not exist!");
}
// Parse materials.xml file and get root element
doc.load_file(filename.c_str());
}
// Parse materials.xml file and get root element
doc.load_file(filename.c_str());
// Loop over XML material elements and populate the array.
pugi::xml_node root = doc.document_element();

View file

@ -160,7 +160,6 @@ StructuredMesh::bin_label(int bin) const {
//==============================================================================
UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) {
n_dimension_ = 3;
// check the mesh type
if (check_for_node(node, "type")) {
@ -1582,14 +1581,22 @@ MOABMesh::MOABMesh(const std::string& filename) {
initialize();
}
MOABMesh::MOABMesh(std::shared_ptr<moab::Interface> external_mbi) {
mbi_ = external_mbi;
filename_ = "unknown (external file)";
this->initialize();
}
void MOABMesh::initialize() {
// create MOAB instance
mbi_ = make_unique<moab::Core>();
// load unstructured mesh file
moab::ErrorCode rval = mbi_->load_file(filename_.c_str());
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to load the unstructured mesh file: " + filename_);
}
// Create the MOAB interface and load data from file
this->create_interface();
// Initialise MOAB error code
moab::ErrorCode rval = moab::MB_SUCCESS;
// Set the dimension
n_dimension_ = 3;
// set member range of tetrahedral entities
rval = mbi_->get_entities_by_dimension(0, n_dimension_, ehs_);
@ -1619,6 +1626,22 @@ void MOABMesh::initialize() {
build_kdtree(ehs_);
}
void
MOABMesh::create_interface()
{
// Do not create a MOAB instance if one is already in memory
if (mbi_) return;
// create MOAB instance
mbi_ = std::make_shared<moab::Core>();
// load unstructured mesh file
moab::ErrorCode rval = mbi_->load_file(filename_.c_str());
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to load the unstructured mesh file: " + filename_);
}
}
void
MOABMesh::build_kdtree(const moab::Range& all_tets)
{

View file

@ -9,6 +9,7 @@
#include "openmc/capi.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/dagmc.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/hdf5_interface.h"
@ -285,6 +286,10 @@ Particle::event_collide()
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) score_collision_derivative(*this);
#ifdef DAGMC
history().reset();
#endif
}
void
@ -318,9 +323,8 @@ void
Particle::event_death()
{
#ifdef DAGMC
if (settings::dagmc)
history().reset();
#endif
history().reset();
#endif
// Finish particle track output.
if (write_track()) {
@ -377,6 +381,11 @@ Particle::cross_surface()
int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
}
// if we're crossing a CSG surface, make sure the DAG history is reset
#ifdef DAGMC
if (surf->geom_type_ == GeometryType::CSG) history().reset();
#endif
// Handle any applicable boundary conditions.
if (surf->bc_ && settings::run_mode != RunMode::PLOTTING) {
surf->bc_->handle_particle(*this, *surf);
@ -387,17 +396,18 @@ Particle::cross_surface()
// SEARCH NEIGHBOR LISTS FOR NEXT CELL
#ifdef DAGMC
if (settings::dagmc) {
auto cellp = dynamic_cast<DAGCell*>(model::cells[cell_last(0)].get());
// TODO: off-by-one
auto surfp =
dynamic_cast<DAGSurface*>(model::surfaces[std::abs(surface()) - 1].get());
int32_t i_cell = next_cell(cellp, surfp) - 1;
// in DAGMC, we know what the next cell should be
if (surf->geom_type_ == GeometryType::DAG) {
auto surfp = dynamic_cast<DAGSurface*>(surf);
auto cellp = dynamic_cast<DAGCell*>(model::cells[cell_last(n_coord() - 1)].get());
auto univp = static_cast<DAGUniverse*>(model::universes[coord(n_coord() - 1).universe].get());
// determine the next cell for this crossing
int32_t i_cell = next_cell(univp, cellp, surfp) - 1;
// save material and temp
material_last() = material();
sqrtkT_last() = sqrtkT();
// set new cell value
coord(0).cell = i_cell;
coord(n_coord() - 1).cell = i_cell;
cell_instance() = 0;
material() = model::cells[i_cell]->material_[0];
sqrtkT() = model::cells[i_cell]->sqrtkT_[0];
@ -503,13 +513,11 @@ Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
// boundary, it is necessary to redetermine the particle's coordinates in
// the lower universes.
// (unless we're using a dagmc model, which has exactly one universe)
if (!settings::dagmc) {
n_coord() = 1;
if (!neighbor_list_find_cell(*this)) {
mark_as_lost("Couldn't find particle after reflecting from surface " +
std::to_string(surf.id_) + ".");
return;
}
n_coord() = 1;
if (surf.geom_type_ != GeometryType::DAG && !neighbor_list_find_cell(*this)) {
this->mark_as_lost("Couldn't find particle after reflecting from surface "
+ std::to_string(surf.id_) + ".");
return;
}
// Set previous coordinate going slightly past surface crossing

View file

@ -167,8 +167,7 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
int nu = static_cast<int>(nu_t);
if (prn(p.current_seed()) <= (nu_t - nu)) ++nu;
// Begin banking the source neutrons
// First, if our bank is full then don't continue
// If no neutrons were produced then don't continue
if (nu == 0) return;
// Initialize the counter of delayed neutrons encountered for each delayed
@ -179,13 +178,16 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
p.nu_bank().clear();
p.fission() = true;
int skipped = 0;
// Determine whether to place fission sites into the shared fission bank
// or the secondary particle bank.
bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
for (int i = 0; i < nu; ++i) {
// Counter for the number of fission sites successfully stored to the shared
// fission bank or the secondary particle bank
int n_sites_stored;
for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {
// Initialize fission site object with particle data
SourceSite site;
site.r = p.r();
@ -203,8 +205,14 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
int64_t idx = simulation::fission_bank.thread_safe_append(site);
if (idx == -1) {
warning("The shared fission bank is full. Additional fission sites created "
"in this generation will not be banked.");
skipped++;
"in this generation will not be banked. Results may be non-deterministic.");
// Decrement number of particle progeny as storage was unsuccessful. This
// step is needed so that the sum of all progeny is equal to the size
// of the shared fission bank.
p.n_progeny()--;
// Break out of loop as no more sites can be added to fission bank
break;
}
} else {
@ -229,14 +237,14 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
// If shared fission bank was full, and no fissions could be added,
// set the particle fission flag to false.
if (nu == skipped) {
if (n_sites_stored == 0) {
p.fission() = false;
return;
}
// If shared fission bank was full, but some fissions could be added,
// reduce nu accordingly
nu -= skipped;
// Set nu to the number of fission sites successfully stored. If the fission
// bank was not found to be full then these values are already equivalent.
nu = n_sites_stored;
// Store the total weight banked for analog fission tallies
p.n_bank() = nu;

View file

@ -107,8 +107,7 @@ create_fission_sites(Particle& p)
nu++;
}
// Begin banking the source neutrons
// First, if our bank is full then don't continue
// If no neutrons were produced then don't continue
if (nu == 0) return;
// Initialize the counter of delayed neutrons encountered for each delayed
@ -119,13 +118,16 @@ create_fission_sites(Particle& p)
p.nu_bank().clear();
p.fission() = true;
int skipped = 0;
// Determine whether to place fission sites into the shared fission bank
// or the secondary particle bank.
bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
for (int i = 0; i < nu; ++i) {
// Counter for the number of fission sites successfully stored to the shared
// fission bank or the secondary particle bank
int n_sites_stored;
for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {
// Initialize fission site object with particle data
SourceSite site;
site.r = p.r();
@ -162,8 +164,14 @@ create_fission_sites(Particle& p)
int64_t idx = simulation::fission_bank.thread_safe_append(site);
if (idx == -1) {
warning("The shared fission bank is full. Additional fission sites created "
"in this generation will not be banked.");
skipped++;
"in this generation will not be banked. Results may be non-deterministic.");
// Decrement number of particle progeny as storage was unsuccessful. This
// step is needed so that the sum of all progeny is equal to the size
// of the shared fission bank.
p.n_progeny()--;
// Break out of loop as no more sites can be added to fission bank
break;
}
} else {
@ -188,14 +196,14 @@ create_fission_sites(Particle& p)
// If shared fission bank was full, and no fissions could be added,
// set the particle fission flag to false.
if (nu == skipped) {
if (n_sites_stored == 0) {
p.fission() = false;
return;
}
// If shared fission bank was full, but some fissions could be added,
// reduce nu accordingly
nu -= skipped;
// Set nu to the number of fission sites successfully stored. If the fission
// bank was not found to be full then these values are already equivalent.
nu = n_sites_stored;
// Store the total weight banked for analog fission tallies
p.n_bank() = nu;

View file

@ -93,14 +93,6 @@ void write_geometry(hid_t file)
{
auto geom_group = create_group(file, "geometry");
#ifdef DAGMC
if (settings::dagmc) {
write_attribute(geom_group, "dagmc", 1);
close_group(geom_group);
return;
}
#endif
write_attribute(geom_group, "n_cells", model::cells.size());
write_attribute(geom_group, "n_surfaces", model::surfaces.size());
write_attribute(geom_group, "n_universes", model::universes.size());

View file

@ -9,7 +9,6 @@
#include "openmc/array.h"
#include "openmc/container_util.h"
#include "openmc/dagmc.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
@ -198,21 +197,21 @@ Surface::diffuse_reflect(Position r, Direction u, uint64_t* seed) const
return u/u.norm();
}
CSGSurface::CSGSurface() : Surface{} {};
CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface{surf_node} {};
void
CSGSurface::to_hdf5(hid_t group_id) const
Surface::to_hdf5(hid_t group_id) const
{
std::string group_name {"surface "};
group_name += std::to_string(id_);
hid_t surf_group = create_group(group_id, fmt::format("surface {}", id_));
hid_t surf_group = create_group(group_id, group_name);
if (geom_type_ == GeometryType::DAG) {
write_string(surf_group, "geom_type", "dagmc", false);
} else if (geom_type_ == GeometryType::CSG) {
write_string(surf_group, "geom_type", "csg", false);
if (bc_) {
write_string(surf_group, "boundary_type", bc_->type(), false);
} else {
write_string(surf_group, "boundary_type", "transmission", false);
if (bc_) {
write_string(surf_group, "boundary_type", bc_->type(), false);
} else {
write_string(surf_group, "boundary_type", "transmission", false);
}
}
if (!name_.empty()) {
@ -224,60 +223,12 @@ CSGSurface::to_hdf5(hid_t group_id) const
close_group(surf_group);
}
//==============================================================================
// DAGSurface implementation
//==============================================================================
#ifdef DAGMC
DAGSurface::DAGSurface() : Surface{} {} // empty constructor
double DAGSurface::evaluate(Position r) const
{
return 0.0;
}
double
DAGSurface::distance(Position r, Direction u, bool coincident) const
{
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(surf, pnt, dir, hit_surf, dist, NULL, 0, 0);
MB_CHK_ERR_CONT(rval);
if (dist < 0.0) dist = INFTY;
return dist;
}
Direction DAGSurface::normal(Position r) const
{
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
double pnt[3] = {r.x, r.y, r.z};
double dir[3];
rval = dagmc_ptr_->get_angle(surf, pnt, dir);
MB_CHK_ERR_CONT(rval);
return dir;
}
Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const
{
Expects(p);
p->history().reset_to_last_intersection();
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
double pnt[3] = {r.x, r.y, r.z};
double dir[3];
rval = dagmc_ptr_->get_angle(surf, pnt, dir, &p->history());
MB_CHK_ERR_CONT(rval);
p->last_dir() = u.reflect(dir);
return p->last_dir();
}
void DAGSurface::to_hdf5(hid_t group_id) const {}
#endif
CSGSurface::CSGSurface() : Surface{} {
geom_type_ = GeometryType::CSG;
};
CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface{surf_node} {
geom_type_ = GeometryType::CSG;
};
//==============================================================================
// Generic functions for x-, y-, and z-, planes.
@ -1036,9 +987,6 @@ void read_surfaces(pugi::xml_node node)
// Count the number of surfaces
int n_surfaces = 0;
for (pugi::xml_node surf_node : node.children("surface")) {n_surfaces++;}
if (n_surfaces == 0) {
fatal_error("No surfaces found in geometry.xml!");
}
// Loop over XML surface elements and populate the array. Keep track of
// periodic surfaces.
@ -1183,19 +1131,6 @@ void read_surfaces(pugi::xml_node node)
surf2.bc_ = surf1.bc_;
}
}
// Check to make sure a boundary condition was applied to at least one
// surface
bool boundary_exists = false;
for (const auto& surf : model::surfaces) {
if (surf->bc_) {
boundary_exists = true;
break;
}
}
if (settings::run_mode != RunMode::PLOTTING && !boundary_exists) {
fatal_error("No boundary conditions were applied to any surfaces!");
}
}
void free_memory_surfaces()

View file

@ -366,9 +366,18 @@ Tally::set_filters(gsl::span<Filter*> filters)
}
}
// Set the strides.
set_strides();
}
void
Tally::set_strides()
{
// Set the strides. Filters are traversed in reverse so that the last filter
// has the shortest stride in memory and the first filter has the longest
// stride.
auto n = filters_.size();
strides_.resize(n, 0);
int stride = 1;
for (int i = n-1; i >= 0; --i) {

View file

@ -1,4 +1,8 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<dagmc_universe filename="dagmc.h5m" id="1" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="40" name="no-void fuel">
<density units="g/cc" value="11" />
@ -22,7 +26,6 @@
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
<dagmc>true</dagmc>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>

View file

@ -26,6 +26,10 @@ def model():
model.settings.dagmc = True
# geometry
dag_univ = openmc.DAGMCUniverse("dagmc.h5m")
model.geometry = openmc.Geometry(dag_univ)
# tally
tally = openmc.Tally()
tally.scores = ['total']

View file

@ -1,4 +1,11 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<dagmc_universe auto_geom_ids="true" filename="dagmc.h5m" id="9" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
@ -9,12 +16,11 @@
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
<dagmc>true</dagmc>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<filter id="1" type="cell">
<bins>1</bins>
<bins>2</bins>
</filter>
<tally id="1">
<filters>1</filters>

View file

@ -27,13 +27,18 @@ class UWUWTest(PyAPITestHarness):
model.settings.export_to_xml()
# geometry
dag_univ = openmc.DAGMCUniverse("dagmc.h5m", auto_geom_ids=True)
model.geometry = openmc.Geometry(dag_univ)
# tally
tally = openmc.Tally()
tally.scores = ['total']
tally.filters = [openmc.CellFilter(1)]
tally.filters = [openmc.CellFilter(2)]
model.tallies = [tally]
model.tallies.export_to_xml()
model.export_to_xml()
def test_refl():
harness = UWUWTest('statepoint.5.h5')

Binary file not shown.

View file

@ -0,0 +1,55 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell fill="10" id="13" region="9 -10 11 -12 13 -14" universe="11" />
<dagmc_universe auto_geom_ids="true" filename="dagmc.h5m" id="9" />
<lattice id="10">
<pitch>24.0 24.0</pitch>
<dimension>2 2</dimension>
<lower_left>-24.0 -24.0</lower_left>
<universes>
9 9
9 9 </universes>
</lattice>
<surface boundary="reflective" coeffs="-24.0" id="9" name="left" type="x-plane" />
<surface boundary="reflective" coeffs="24.0" id="10" name="right" type="x-plane" />
<surface boundary="reflective" coeffs="-24.0" id="11" name="front" type="y-plane" />
<surface boundary="reflective" coeffs="24.0" id="12" name="back" type="y-plane" />
<surface boundary="reflective" coeffs="-10.0" id="13" name="bottom" type="z-plane" />
<surface boundary="reflective" coeffs="10.0" id="14" name="top" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="13" name="no-void fuel">
<density units="g/cc" value="10.29769" />
<nuclide ao="0.93120485" name="U234" />
<nuclide ao="0.00055815" name="U235" />
<nuclide ao="0.022408" name="U238" />
<nuclide ao="0.045829" name="O16" />
</material>
<material id="14" name="clad">
<density units="g/cc" value="6.55" />
<nuclide ao="0.021827" name="Zr90" />
<nuclide ao="0.00476" name="Zr91" />
<nuclide ao="0.0072758" name="Zr92" />
<nuclide ao="0.0073734" name="Zr94" />
<nuclide ao="0.0011879" name="Zr96" />
</material>
<material id="15" name="water">
<density units="g/cc" value="0.740582" />
<nuclide ao="0.049457" name="H1" />
<nuclide ao="0.024672" name="O16" />
<nuclide ao="8.0042e-06" name="B10" />
<nuclide ao="3.2218e-05" name="B11" />
<sab name="c_H_in_H2O" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>2</inactive>
<output>
<summary>false</summary>
</output>
</settings>

View file

@ -0,0 +1,2 @@
k-combined:
9.436168E-01 2.905559E-02

View file

@ -0,0 +1,80 @@
import numpy as np
import openmc
import openmc.lib
import pytest
from tests.testing_harness import PyAPITestHarness
pytestmark = pytest.mark.skipif(
not openmc.lib._dagmc_enabled(),
reason="DAGMC CAD geometry is not enabled.")
class DAGMCUniverseTest(PyAPITestHarness):
def _build_inputs(self):
model = openmc.model.Model()
### MATERIALS ###
fuel = openmc.Material(name='no-void fuel')
fuel.set_density('g/cc', 10.29769)
fuel.add_nuclide('U234', 0.93120485)
fuel.add_nuclide('U235', 0.00055815)
fuel.add_nuclide('U238', 0.022408)
fuel.add_nuclide('O16', 0.045829)
cladding = openmc.Material(name='clad')
cladding.set_density('g/cc', 6.55)
cladding.add_nuclide('Zr90', 0.021827)
cladding.add_nuclide('Zr91', 0.00476)
cladding.add_nuclide('Zr92', 0.0072758)
cladding.add_nuclide('Zr94', 0.0073734)
cladding.add_nuclide('Zr96', 0.0011879)
water = openmc.Material(name='water')
water.set_density('g/cc', 0.740582)
water.add_nuclide('H1', 0.049457)
water.add_nuclide('O16', 0.024672)
water.add_nuclide('B10', 8.0042e-06)
water.add_nuclide('B11', 3.2218e-05)
water.add_s_alpha_beta('c_H_in_H2O')
model.materials = openmc.Materials([fuel, cladding, water])
### GEOMETRY ###
# create the DAGMC universe
pincell_univ = openmc.DAGMCUniverse(filename='dagmc.h5m', auto_geom_ids=True)
# create a 2 x 2 lattice using the DAGMC pincell
pitch = np.asarray((24.0, 24.0))
lattice = openmc.RectLattice()
lattice.pitch = pitch
lattice.universes = [[pincell_univ] * 2] * 2
lattice.lower_left = -pitch
left = openmc.XPlane(x0=-pitch[0], name='left', boundary_type='reflective')
right = openmc.XPlane(x0=pitch[0], name='right', boundary_type='reflective')
front = openmc.YPlane(y0=-pitch[1], name='front', boundary_type='reflective')
back = openmc.YPlane(y0=pitch[1], name='back', boundary_type='reflective')
# clip the DAGMC geometry at +/- 10 cm w/ CSG planes
bottom = openmc.ZPlane(z0=-10.0, name='bottom', boundary_type='reflective')
top = openmc.ZPlane(z0=10.0, name='top', boundary_type='reflective')
bounding_region = +left & -right & +front & -back & +bottom & -top
bounding_cell = openmc.Cell(fill=lattice, region=bounding_region)
model.geometry = openmc.Geometry([bounding_cell])
# settings
model.settings.particles = 100
model.settings.batches = 10
model.settings.inactive = 2
model.settings.output = {'summary' : False}
model.export_to_xml()
def test_univ():
harness = DAGMCUniverseTest('statepoint.10.h5')
harness.main()

View file

@ -1,4 +1,11 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<dagmc_universe filename="dagmc.h5m" id="9" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
@ -9,7 +16,6 @@
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
<dagmc>true</dagmc>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>

View file

@ -27,13 +27,17 @@ class UWUWTest(PyAPITestHarness):
model.settings.export_to_xml()
# geometry
dag_univ = openmc.DAGMCUniverse("dagmc.h5m")
model.geometry = openmc.Geometry(root=dag_univ)
# tally
tally = openmc.Tally()
tally.scores = ['total']
tally.filters = [openmc.CellFilter(1)]
model.tallies = [tally]
model.tallies.export_to_xml()
model.export_to_xml()
def test_uwuw():
harness = UWUWTest('statepoint.5.h5')

View file

@ -0,0 +1,73 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" name="fuel" region="1 -2 3 -4 5 -6" universe="1" />
<cell id="2" material="2" name="clad" region="(-1 | 2 | -3 | 4 | -5 | 6) (7 -8 9 -10 11 -12)" universe="1" />
<cell id="3" material="3" name="water" region="(-7 | 8 | -9 | 10 | -11 | 12) (13 -14 15 -16 17 -18)" universe="1" />
<surface coeffs="-5.0" id="1" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="2" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="3" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="4" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="5" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="6" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="7" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="8" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="9" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-10" id="13" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="10" id="14" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-10" id="15" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="10" id="16" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-10" id="17" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="10" id="18" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="1" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="2" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="3" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" library="moab" type="unstructured">
<filename>test_mesh_tets.h5m</filename>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<tally id="1" name="unstructured mesh tally">
<filters>1</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
</tallies>

View file

@ -0,0 +1,98 @@
#include "moab/Core.hpp"
#include "openmc/capi.h"
#include "openmc/error.h"
#include "openmc/mesh.h"
#include "openmc/tallies/filter_mesh.h"
#include "openmc/tallies/tally.h"
#include <iostream>
int main(int argc, char* argv[])
{
using namespace openmc;
int openmc_err;
// Initialise OpenMC
openmc_err = openmc_init(argc, argv, nullptr);
if (openmc_err == -1) {
// This happens for the -h and -v flags
return EXIT_SUCCESS;
} else if (openmc_err) {
fatal_error(openmc_err_msg);
}
// Create MOAB interface
std::shared_ptr<moab::Interface> moabPtrLocal =
std::make_shared<moab::Core>();
// Load unstructured mesh file
std::string filename = "test_mesh_tets.h5m";
moab::ErrorCode rval = moabPtrLocal->load_file(filename.c_str());
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to load the unstructured mesh file: " + filename);
} else {
std::cout << "Loaded external MOAB mesh from file " << filename
<< std::endl;
}
// Add a new unstructured mesh to openmc using new constructor
model::meshes.push_back(std::make_unique<MOABMesh>(moabPtrLocal));
// Check we now have 2 copies of the shared ptr
if (moabPtrLocal.use_count() != 2) {
fatal_error("Incorrect number of MOAB shared pointers");
}
// Auto-assign mesh ID
model::meshes.back()->set_id(C_NONE);
int mesh_id = model::meshes.back()->id_;
// Check we now have 2 meshes and id was correctly set
if (model::meshes.size() != 2)
fatal_error("Wrong number of meshes.");
else if (mesh_id != 2)
fatal_error("Mesh ID is incorrect");
// Add a new mesh filter with auto-assigned ID
Filter* filter_ptr = Filter::create("mesh", C_NONE);
// Upcast pointer type
MeshFilter* mesh_filter = dynamic_cast<MeshFilter*>(filter_ptr);
if (!mesh_filter) {
fatal_error("Failed to create mesh filter");
}
// Pass in the index of our mesh to the filter
int32_t mesh_idx = model::meshes.size() - 1;
mesh_filter->set_mesh(mesh_idx);
// Create a tally with auto-assigned ID
model::tallies.push_back(make_unique<Tally>(C_NONE));
// Set tally name - matches that in test.py
model::tallies.back()->name_ = "external mesh tally";
// Set tally filter to our mesh filter
std::vector<Filter*> filters(1, filter_ptr);
model::tallies.back()->set_filters(filters);
// Set tally estimator
model::tallies.back()->estimator_ = TallyEstimator::TRACKLENGTH;
// Set tally score
std::vector<std::string> score_names(1, "flux");
model::tallies.back()->set_scores(score_names);
// Run OpenMC
openmc_err = openmc_run();
if (openmc_err)
fatal_error(openmc_err_msg);
// Deallocate memory
openmc_err = openmc_finalize();
if (openmc_err)
fatal_error(openmc_err_msg);
return EXIT_SUCCESS;
}

View file

@ -0,0 +1,275 @@
from pathlib import Path
import os
import shutil
import subprocess
from subprocess import CalledProcessError
import textwrap
import glob
from itertools import product
import openmc
import openmc.lib
import numpy as np
import pytest
from tests.regression_tests import config
from tests.testing_harness import PyAPITestHarness
pytestmark = pytest.mark.skipif(
not openmc.lib._dagmc_enabled(),
reason="DAGMC is not enabled.")
TETS_PER_VOXEL = 12
# Test that an external moab instance can be passed in through the C API
@pytest.fixture
def cpp_driver(request):
"""Compile the external source"""
# Get build directory and write CMakeLists.txt file
openmc_dir = Path(str(request.config.rootdir)) / 'build'
with open('CMakeLists.txt', 'w') as f:
f.write(textwrap.dedent("""
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
project(openmc_cpp_driver CXX)
add_executable(main main.cpp)
find_package(OpenMC REQUIRED HINTS {})
target_link_libraries(main OpenMC::libopenmc)
set_target_properties(main PROPERTIES CXX_STANDARD
14 CXX_STANDARD_REQUIRED YES CXX_EXTENSIONS NO)
set(CMAKE_CXX_FLAGS "-pedantic-errors")
add_compile_definitions(DAGMC=1)
""".format(openmc_dir)))
# Create temporary build directory and change to there
local_builddir = Path('build')
local_builddir.mkdir(exist_ok=True)
os.chdir(str(local_builddir))
if config['mpi']:
os.environ['CXX'] = 'mpicxx'
try:
print("Building driver")
# Run cmake/make to build the shared libary
subprocess.run(['cmake', os.path.pardir], check=True)
subprocess.run(['make'], check=True)
os.chdir(os.path.pardir)
yield "./build/main"
finally:
# Remove local build directory when test is complete
shutil.rmtree('build')
os.remove('CMakeLists.txt')
class ExternalMoabTest(PyAPITestHarness):
def __init__(self, executable, statepoint_name, model):
super().__init__(statepoint_name, model)
self.executable = executable
def _run_openmc(self):
if config['mpi']:
mpi_args = [config['mpiexec'], '-n', config['mpi_np']]
openmc.run(openmc_exec=self.executable,
mpi_args=mpi_args,
event_based=config['event'])
else:
openmc.run(openmc_exec=self.executable,
event_based=config['event'])
# Override some methods to do nothing
def _get_results(self):
pass
def _write_results(self, results_string):
pass
def _overwrite_results(self):
pass
def _test_output_created(self):
pass
# Directly compare results of unstructured mesh with internal and
# external moab
def _compare_results(self):
with openmc.StatePoint(self._sp_name) as sp:
# loop over the tallies and get data
ext_data = []
unstr_data = []
for tally in sp.tallies.values():
# Safety check that mesh filter is correct
if tally.contains_filter(openmc.MeshFilter):
flt = tally.find_filter(openmc.MeshFilter)
if isinstance(flt.mesh, openmc.UnstructuredMesh):
if tally.name == "external mesh tally":
ext_data = tally.get_reshaped_data(value='mean')
elif tally.name == "unstructured mesh tally":
unstr_data = tally.get_reshaped_data(value='mean')
# we expect these results to be the same to within at 8
# decimal places
decimals = 8
np.testing.assert_array_almost_equal(unstr_data,
ext_data, decimals)
@staticmethod
def get_mesh_tally_data(tally):
data = tally.get_reshaped_data(value='mean')
std_dev = tally.get_reshaped_data(value='std_dev')
data.shape = (data.size, 1)
std_dev.shape = (std_dev.size, 1)
return np.sum(data, axis=1), np.sum(std_dev, axis=1)
def _cleanup(self):
super()._cleanup()
output = glob.glob('tally*.vtk')
for f in output:
if os.path.exists(f):
os.remove(f)
def test_external_mesh(cpp_driver):
# Materials
materials = openmc.Materials()
fuel_mat = openmc.Material(name="fuel")
fuel_mat.add_nuclide("U235", 1.0)
fuel_mat.set_density('g/cc', 4.5)
materials.append(fuel_mat)
zirc_mat = openmc.Material(name="zircaloy")
zirc_mat.add_element("Zr", 1.0)
zirc_mat.set_density("g/cc", 5.77)
materials.append(zirc_mat)
water_mat = openmc.Material(name="water")
water_mat.add_nuclide("H1", 2.0)
water_mat.add_nuclide("O16", 1.0)
water_mat.set_density("atom/b-cm", 0.07416)
materials.append(water_mat)
materials.export_to_xml()
# Geometry
fuel_min_x = openmc.XPlane(-5.0, name="minimum x")
fuel_max_x = openmc.XPlane(5.0, name="maximum x")
fuel_min_y = openmc.YPlane(-5.0, name="minimum y")
fuel_max_y = openmc.YPlane(5.0, name="maximum y")
fuel_min_z = openmc.ZPlane(-5.0, name="minimum z")
fuel_max_z = openmc.ZPlane(5.0, name="maximum z")
fuel_cell = openmc.Cell(name="fuel")
fuel_cell.region = +fuel_min_x & -fuel_max_x & \
+fuel_min_y & -fuel_max_y & \
+fuel_min_z & -fuel_max_z
fuel_cell.fill = fuel_mat
clad_min_x = openmc.XPlane(-6.0, name="minimum x")
clad_max_x = openmc.XPlane(6.0, name="maximum x")
clad_min_y = openmc.YPlane(-6.0, name="minimum y")
clad_max_y = openmc.YPlane(6.0, name="maximum y")
clad_min_z = openmc.ZPlane(-6.0, name="minimum z")
clad_max_z = openmc.ZPlane(6.0, name="maximum z")
clad_cell = openmc.Cell(name="clad")
clad_cell.region = (-fuel_min_x | +fuel_max_x |
-fuel_min_y | +fuel_max_y |
-fuel_min_z | +fuel_max_z) & \
(+clad_min_x & -clad_max_x &
+clad_min_y & -clad_max_y &
+clad_min_z & -clad_max_z)
clad_cell.fill = zirc_mat
bounds = (10, 10, 10)
water_min_x = openmc.XPlane(x0=-bounds[0],
name="minimum x",
boundary_type='vacuum')
water_max_x = openmc.XPlane(x0=bounds[0],
name="maximum x",
boundary_type='vacuum')
water_min_y = openmc.YPlane(y0=-bounds[1],
name="minimum y",
boundary_type='vacuum')
water_max_y = openmc.YPlane(y0=bounds[1],
name="maximum y",
boundary_type='vacuum')
water_min_z = openmc.ZPlane(z0=-bounds[2],
name="minimum z",
boundary_type='vacuum')
water_max_z = openmc.ZPlane(z0=bounds[2],
name="maximum z",
boundary_type='vacuum')
water_cell = openmc.Cell(name="water")
water_cell.region = (-clad_min_x | +clad_max_x |
-clad_min_y | +clad_max_y |
-clad_min_z | +clad_max_z) & \
(+water_min_x & -water_max_x &
+water_min_y & -water_max_y &
+water_min_z & -water_max_z)
water_cell.fill = water_mat
# create a containing universe
geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell])
# Meshes
mesh_filename = "test_mesh_tets.h5m"
# Create a normal unstructured mesh to compare to
uscd_mesh = openmc.UnstructuredMesh(mesh_filename, 'moab')
# Create filters
uscd_filter = openmc.MeshFilter(mesh=uscd_mesh)
# Tallies
tallies = openmc.Tallies()
uscd_tally = openmc.Tally(name="unstructured mesh tally")
uscd_tally.filters = [uscd_filter]
uscd_tally.scores = ['flux']
uscd_tally.estimator = 'tracklength'
tallies.append(uscd_tally)
# Settings
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.particles = 100
settings.batches = 10
# Source setup
space = openmc.stats.Point()
angle = openmc.stats.Monodirectional((-1.0, 0.0, 0.0))
energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0])
source = openmc.Source(space=space, energy=energy, angle=angle)
settings.source = source
model = openmc.model.Model(geometry=geometry,
materials=materials,
tallies=tallies,
settings=settings)
harness = ExternalMoabTest(cpp_driver,
'statepoint.10.h5',
model)
# Run open MC and check results
harness.main()

View file

@ -29,7 +29,9 @@ def dagmc_model(request):
source = openmc.Source(space=source_box)
model.settings.source = source
model.settings.dagmc = True
# geometry
dagmc_universe = openmc.DAGMCUniverse('dagmc.h5m')
model.geometry = openmc.Geometry(dagmc_universe)
# tally
tally = openmc.Tally()

View file

@ -54,7 +54,6 @@ def test_export_to_xml(run_in_tmpdir):
s.log_grid_bins = 2000
s.photon_transport = False
s.electron_treatment = 'led'
s.dagmc = False
# Make sure exporting XML works
s.export_to_xml()
@ -111,4 +110,3 @@ def test_export_to_xml(run_in_tmpdir):
assert s.log_grid_bins == 2000
assert not s.photon_transport
assert s.electron_treatment == 'led'
assert not s.dagmc