diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 68848014d1..e5c0ca1600 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -88,6 +88,14 @@ you care. This element has the following attributes/sub-elements: *Default*: 0.0 +-------------------------------- +```` Element +-------------------------------- + +When the DAGMC mode is enabled, the OpenMC geometry will be read from the file +``dagmc.h5m``. If a :ref:`geometry.xml ` file is present with +``dagmc`` set to ``true``, it will be ignored. + -------------------------------- ```` Element -------------------------------- diff --git a/docs/source/usersguide/geometry.rst b/docs/source/usersguide/geometry.rst index bb800699b5..70a665c53c 100644 --- a/docs/source/usersguide/geometry.rst +++ b/docs/source/usersguide/geometry.rst @@ -398,3 +398,32 @@ if needed, lattices, the last step is to create an instance of .. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry .. _quadratic surfaces: http://en.wikipedia.org/wiki/Quadric + +-------------------------- +Using CAD-based Geometry +-------------------------- + +OpenMC relies on the Direct Accelerated Geometry Monte Carlo toolkit (`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:: + + settings = openmc.Settings() + settings.dagmc = True + +or in the :ref:`settings.xml ` file:: + + true + +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. + + **Note:** DAGMC geometries used in OpenMC are currently required to be clean, + meaning that all surfaces have been `imprinted and merged + `_ + successfully and that the model is `watertight + `_. Future + implementations of DAGMC geometry will support small volume overlaps and + un-merged surfaces. diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 48dcd8c782..d3e59c06e1 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -159,14 +159,25 @@ Prerequisites sudo apt install mpich libmpich-dev sudo apt install openmpi-bin libopenmpi-dev + * DAGMC_ toolkit for simulation using CAD-based geometries + + OpenMC supports particle tracking in CAD-based geometries via the Direct + Accelerated Geometry Monte Carlo (DAGMC) toolkit (`installation + instructions + `_). For use in + OpenMC, only the ``MOAB_DIR`` and ``BUILD_TALLY`` variables need to be + specified in the CMake configuration step. + | * git_ version control software for obtaining source code + .. _gfortran: http://gcc.gnu.org/wiki/GFortran .. _gcc: https://gcc.gnu.org/ .. _CMake: http://www.cmake.org .. _OpenMPI: http://www.open-mpi.org .. _MPICH: http://www.mpich.org .. _HDF5: https://www.hdfgroup.org/solutions/hdf5/ +.. _DAGMC: https://svalinn.github.io/DAGMC/index.html Obtaining the Source -------------------- @@ -236,6 +247,12 @@ openmp Enables shared-memory parallelism using the OpenMP API. The Fortran compiler being used must support OpenMP. (Default: on) +dagmc + Enables use of CAD-based DAGMC_ geometries. Please see the note about DAGMC in + the optional dependencies list for more information on this feature. The + installation directory for DAGMC should also be defined as `DAGMC_ROOT` in the + CMake configuration command. (Default: off) + coverage Compile and link code instrumented for coverage analysis. This is typically used in conjunction with gcov_. diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 84e5f5318f..6c3338854c 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -7,6 +7,7 @@ #ifdef __cplusplus #include "openmc/bank.h" + extern "C" { int openmc_fission_bank(openmc::Bank** ptr, int64_t* n); int openmc_source_bank(openmc::Bank** ptr, int64_t* n); diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 51d97d7375..4a129ef095 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -2,6 +2,10 @@ #ifndef OPENMC_DAGMC_H #define OPENMC_DAGMC_H +namespace openmc { +extern "C" const bool dagmc_enabled; +} + #ifdef DAGMC #include "DagMC.hpp" @@ -10,10 +14,6 @@ namespace openmc { -//============================================================================== -// Global variables -//============================================================================== - namespace model { extern moab::DagMC* DAG; @@ -26,15 +26,11 @@ extern moab::DagMC* DAG; extern "C" void load_dagmc_geometry(); extern "C" void free_memory_dagmc(); +extern "C" pugi::xml_document* read_uwuw_materials(); +bool get_uwuw_materials_xml(std::string& s); } // namespace openmc #endif // DAGMC #endif // OPENMC_DAGMC_H - -#ifdef DAGMC -extern "C" constexpr bool dagmc_enabled = true; -#else -extern "C" constexpr bool dagmc_enabled = false; -#endif diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 87d5a2f50d..23af2f624e 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -2,6 +2,9 @@ #include "openmc/constants.h" #include "openmc/container_util.h" +#ifdef DAGMC +#include "openmc/dagmc.h" +#endif #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/hdf5_interface.h" @@ -90,15 +93,30 @@ extern "C" void read_mg_cross_sections_header(); void read_cross_sections_xml() { - // Check if materials.xml exists + 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 - pugi::xml_document doc; doc.load_file(filename.c_str()); +#ifdef DAGMC + } +#endif + auto root = doc.document_element(); // Find cross_sections.xml file -- the first place to look is the diff --git a/src/dagmc.cpp b/src/dagmc.cpp index e6ccee24d5..1acef12629 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -1,16 +1,38 @@ #include "openmc/dagmc.h" #include "openmc/cell.h" +#include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/file_utils.h" #include "openmc/string_utils.h" #include "openmc/settings.h" #include "openmc/geometry.h" +#ifdef DAGMC + +#include "uwuw.hpp" +#include "dagmcmetadata.hpp" + +#endif + #include #include #include +#include + +namespace openmc { #ifdef DAGMC +const bool dagmc_enabled = true; +#else +const bool dagmc_enabled = false; +#endif + +} + +#ifdef DAGMC + +const std::string DAGMC_FILENAME = "dagmc.h5m"; namespace openmc { @@ -20,28 +42,96 @@ moab::DagMC* DAG; } // namespace model + +bool get_uwuw_materials_xml(std::string& s) { + UWUW uwuw(DAGMC_FILENAME.c_str()); + + std::stringstream ss; + bool uwuw_mats_present = false; + if (uwuw.material_library.size() != 0) { + uwuw_mats_present = true; + // write header + ss << "\n"; + ss << "\n"; + const auto& mat_lib = uwuw.material_library; + // write materials + for (auto mat : mat_lib) { ss << mat.second.openmc("atom"); } + // write footer + ss << ""; + s = ss.str(); + } + + return uwuw_mats_present; +} + +pugi::xml_document* read_uwuw_materials() { + pugi::xml_document* doc = nullptr; + + std::string s; + bool found_uwuw_mats = get_uwuw_materials_xml(s); + if (found_uwuw_mats) { + doc = new pugi::xml_document(); + pugi::xml_parse_result result = doc->load_string(s.c_str()); + } + return doc; +} + +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(); + } + + return found_uwuw_mats; +} + void load_dagmc_geometry() { if (!model::DAG) { model::DAG = new moab::DagMC(); } - int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC + /// Materials \\\ - moab::ErrorCode rval = model::DAG->load_file("dagmc.h5m"); + // create uwuw instance + UWUW uwuw(DAGMC_FILENAME.c_str()); + + // check for uwuw material definitions + bool using_uwuw = !uwuw.material_library.empty(); + + // notify user if UWUW materials are going to be used + if (using_uwuw) { + std::cout << "Found UWUW Materials in the DAGMC geometry file.\n"; + } + + int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC runs + + // load the DAGMC geometry + moab::ErrorCode rval = model::DAG->load_file(DAGMC_FILENAME.c_str()); MB_CHK_ERR_CONT(rval); + // initialize acceleration data structures rval = model::DAG->init_OBBTree(); MB_CHK_ERR_CONT(rval); - std::vector prop_keywords; - prop_keywords.push_back("mat"); - prop_keywords.push_back("boundary"); - - std::map ph; - model::DAG->parse_properties(prop_keywords, ph, ":"); + // parse model metadata + dagmcMetaData DMD(model::DAG); + if (using_uwuw) { + DMD.load_property_data(); + } + std::vector keywords {"temp", "mat", "density", "boundary"}; + std::map dum; + std::string delimiters = ":/"; + rval = model::DAG->parse_properties(keywords, dum, delimiters.c_str()); MB_CHK_ERR_CONT(rval); + /// Cells (Volumes) \\\ + // initialize cell objects model::n_cells = model::DAG->num_entities(3); @@ -69,35 +159,80 @@ void load_dagmc_geometry() model::universes[it->second]->cells_.push_back(i); } + // check for temperature assignment + std::string temp_value; + if (model::DAG->has_prop(vol_handle, "temp")) { + rval = model::DAG->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 { + c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * settings::temperature_default)); + } + + // MATERIALS + if (model::DAG->is_implicit_complement(vol_handle)) { - // assuming implicit complement is void for now - c->material_.push_back(MATERIAL_VOID); + if (model::DAG->has_prop(vol_handle, "mat")) { + // if the implicit complement has been assigned a material, use it + std::string comp_mat = DMD.volume_material_property_data_eh[vol_handle]; + // Note: material numbers are set by UWUW + int mat_number = uwuw.material_library[comp_mat].metadata["mat_number"].asInt(); + c->material_.push_back(mat_number); + } else { + // if no material is found, the implicit complement is void + c->material_.push_back(MATERIAL_VOID); + } continue; } - if (model::DAG->has_prop(vol_handle, "mat")){ - std::string mat_value; + // determine volume material assignment + std::string mat_value; + if (model::DAG->has_prop(vol_handle, "mat")) { rval = model::DAG->prop_value(vol_handle, "mat", mat_value); MB_CHK_ERR_CONT(rval); - to_lower(mat_value); - - if (mat_value == "void" || mat_value == "vacuum") { - c->material_.push_back(MATERIAL_VOID); - } else { - c->material_.push_back(std::stoi(mat_value)); - } } else { std::stringstream err_msg; err_msg << "Volume " << c->id_ << " has no material assignment."; fatal_error(err_msg.str()); } + + std::string cmp_str = mat_value; + to_lower(cmp_str); + // material void checks + if (cmp_str.find("void") != std::string::npos || + cmp_str.find("vacuum") != std::string::npos || + cmp_str.find("graveyard") != std::string::npos) { + c->material_.push_back(MATERIAL_VOID); + } else { + if (using_uwuw) { + // lookup material in uwuw if the were present + std::string uwuw_mat = DMD.volume_material_property_data_eh[vol_handle]; + if (uwuw.material_library.count(uwuw_mat) != 0) { + // Note: material numbers are set by UWUW + int mat_number = uwuw.material_library[uwuw_mat].metadata["mat_number"].asInt(); + c->material_.push_back(mat_number); + } else { + std::stringstream err_msg; + err_msg << "Material with value " << mat_value << " not found "; + err_msg << "in the UWUW material library"; + fatal_error(err_msg); + } + } else { + // if not using UWUW materials, we'll find this material + // later in the materials.xml + c->material_.push_back(std::stoi(mat_value)); + } + } } - // Allocate the cell overlap count if necessary. + // allocate the cell overlap count if necessary if (settings::check_overlaps) { model::overlap_check_count.resize(model::cells.size(), 0); } + /// Surfaces \\\ + // initialize surface objects int n_surfaces = model::DAG->num_entities(2); model::surfaces.resize(n_surfaces); @@ -110,8 +245,9 @@ void load_dagmc_geometry() s->id_ = model::DAG->id_by_index(2, i+1); s->dagmc_ptr_ = model::DAG; + // set BCs + std::string bc_value; if (model::DAG->has_prop(surf_handle, "boundary")) { - std::string bc_value; rval = model::DAG->prop_value(surf_handle, "boundary", bc_value); MB_CHK_ERR_CONT(rval); to_lower(bc_value); @@ -130,7 +266,8 @@ void load_dagmc_geometry() << "\" specified on surface " << s->id_; fatal_error(err_msg); } - } else { // if no BC property is found, set to transmit + } else { + // if no condition is found, set to transmit s->bc_ = BC_TRANSMIT; } @@ -147,5 +284,6 @@ void free_memory_dagmc() delete model::DAG; } + } #endif diff --git a/src/dagmc_header.F90 b/src/dagmc_header.F90 index 9895791ac0..4a97f7f0ce 100644 --- a/src/dagmc_header.F90 +++ b/src/dagmc_header.F90 @@ -13,6 +13,11 @@ module dagmc_header subroutine free_memory_dagmc() bind(C) end subroutine free_memory_dagmc + function read_uwuw_materials() result(doc) bind(C) + import C_PTR + type(C_PTR) :: doc + end function read_uwuw_materials + end interface end module dagmc_header diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 102a95dfb4..d78135d77f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -404,18 +404,29 @@ contains ! Display output message call write_message("Reading materials XML file...", 5) - ! Check if materials.xml exists - filename = trim(path_input) // "materials.xml" - inquire(FILE=filename, EXIST=file_exists) - if (.not. file_exists) then - call fatal_error("Material XML file '" // trim(filename) // "' does not & - &exist!") + doc % ptr = C_NULL_PTR + +#ifdef DAGMC + if (dagmc) then + doc % ptr = read_uwuw_materials() + end if +#endif + + if (.not. c_associated(doc % ptr)) then + ! Check if materials.xml exists + filename = trim(path_input) // "materials.xml" + inquire(FILE=filename, EXIST=file_exists) + if (.not. file_exists) then + call fatal_error("Material XML file '" // trim(filename) // "' does not & + &exist!") end if ! Parse materials.xml file call doc % load_file(filename) - root = doc % document_element() + end if + + root = doc % document_element() call read_materials(root % ptr) ! Set total number of nuclides and elements diff --git a/src/settings.cpp b/src/settings.cpp index dbbe9c6dad..78f90082c7 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -405,13 +405,6 @@ void read_settings_xml() #endif } -#ifdef _OPENMP - if (dagmc && omp_get_max_threads() > 1) { - warning("Forcing number of threads to 1 for DAGMC simulation."); - omp_set_num_threads(1); - } -#endif - // ========================================================================== // EXTERNAL SOURCE diff --git a/src/tracking.F90 b/src/tracking.F90 index 43df277156..d6d6d5534c 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -499,7 +499,7 @@ contains p % coord(1) % cell = i_cell-1 ! decrement for C++ indexing p % cell_instance = 1 p % material = cells(i_cell) % material(1) - p % sqrtKT = cells(i_cell) % sqrtKT(1) + p % sqrtKT = cells(i_cell) % sqrtKT(0) return end if #endif diff --git a/tests/regression_tests/dagmc/dagmc.h5m b/tests/regression_tests/dagmc/dagmc.h5m index 5f84286de2..49c28c4f5e 100644 Binary files a/tests/regression_tests/dagmc/dagmc.h5m and b/tests/regression_tests/dagmc/dagmc.h5m differ diff --git a/tests/regression_tests/dagmc/results_true.dat b/tests/regression_tests/dagmc/results_true.dat index fc13a9cdaf..36475c5bf4 100644 --- a/tests/regression_tests/dagmc/results_true.dat +++ b/tests/regression_tests/dagmc/results_true.dat @@ -1,5 +1,5 @@ k-combined: -1.115067E+00 5.423808E-02 +1.028803E+00 3.340602E-02 tally 1: -8.543144E+00 -1.530584E+01 +8.346847E+00 +1.453407E+01 diff --git a/tests/regression_tests/dagmc/test.py b/tests/regression_tests/dagmc/test.py index 9766c9732f..6c9be63bf6 100644 --- a/tests/regression_tests/dagmc/test.py +++ b/tests/regression_tests/dagmc/test.py @@ -38,8 +38,8 @@ def test_dagmc(): water = openmc.Material() water.add_nuclide('H1', 2.0, 'ao') water.add_nuclide('O16', 1.0, 'ao') - water.add_s_alpha_beta('c_H_in_H2O') water.set_density('g/cc', 1.0) + water.add_s_alpha_beta('c_H_in_H2O') water.id = 41 mats = openmc.Materials([u235, water]) diff --git a/tests/regression_tests/uwuw/__init__.py b/tests/regression_tests/uwuw/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/uwuw/dagmc.h5m b/tests/regression_tests/uwuw/dagmc.h5m new file mode 100644 index 0000000000..8fca1fff14 Binary files /dev/null and b/tests/regression_tests/uwuw/dagmc.h5m differ diff --git a/tests/regression_tests/uwuw/inputs_true.dat b/tests/regression_tests/uwuw/inputs_true.dat new file mode 100644 index 0000000000..87ceb944d8 --- /dev/null +++ b/tests/regression_tests/uwuw/inputs_true.dat @@ -0,0 +1,23 @@ + + + eigenvalue + 100 + 5 + 0 + + + -4 -4 -4 4 4 4 + + + true + + + + + 1 + + + 1 + total + + diff --git a/tests/regression_tests/uwuw/results_true.dat b/tests/regression_tests/uwuw/results_true.dat new file mode 120000 index 0000000000..3711895aba --- /dev/null +++ b/tests/regression_tests/uwuw/results_true.dat @@ -0,0 +1 @@ +../dagmc/results_true.dat \ No newline at end of file diff --git a/tests/regression_tests/uwuw/test.py b/tests/regression_tests/uwuw/test.py new file mode 100644 index 0000000000..f6713c2e4f --- /dev/null +++ b/tests/regression_tests/uwuw/test.py @@ -0,0 +1,41 @@ +import openmc +import openmc.capi +from openmc.stats import Box +from openmc.material import Materials + +import pytest +from tests.testing_harness import PyAPITestHarness + +pytestmark = pytest.mark.skipif( + not openmc.capi._dagmc_enabled(), + reason="DAGMC CAD geometry is not enabled.") + +class UWUWTest(PyAPITestHarness): + + def _build_inputs(self): + model = openmc.model.Model() + + # settings + model.settings.batches = 5 + model.settings.inactive = 0 + model.settings.particles = 100 + + source = openmc.Source(space=Box([-4, -4, -4], + [ 4, 4, 4])) + model.settings.source = source + + model.settings.dagmc = True + + model.settings.export_to_xml() + + # tally + tally = openmc.Tally() + tally.scores = ['total'] + tally.filters = [openmc.CellFilter(1)] + model.tallies = [tally] + + model.tallies.export_to_xml() + +def test_uwuw(): + harness = UWUWTest('statepoint.5.h5') + harness.main()