From e8081207d7df265f495779a5735cdfe5b34f704a Mon Sep 17 00:00:00 2001
From: agnelson <agnelson@anl.gov>
Date: Thu, 30 Sep 2021 16:34:58 -0500
Subject: [PATCH] Modified the C API to support setting material temperatures
 and to support setting the event based flag. Fixed failing MPI test.

---
 include/openmc/capi.h          |  1 +
 include/openmc/material.h      |  8 ++---
 include/openmc/settings.h      |  2 +-
 openmc/executor.py             | 17 ++++++----
 openmc/lib/__init__.py         |  2 ++
 openmc/lib/material.py         |  7 ++++
 openmc/lib/settings.py         |  1 +
 openmc/model/model.py          | 52 ++++++++++++++----------------
 src/material.cpp               | 12 +++++++
 tests/unit_tests/test_lib.py   |  6 ++++
 tests/unit_tests/test_model.py | 59 +++++++++++++++++-----------------
 11 files changed, 100 insertions(+), 67 deletions(-)

diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 0929a11f7..679919dfa 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -88,6 +88,7 @@ int openmc_material_set_densities(
 int openmc_material_set_id(int32_t index, int32_t id);
 int openmc_material_get_name(int32_t index, const char** name);
 int openmc_material_set_name(int32_t index, const char* name);
+int openmc_material_set_temperature(int32_t index, double temperature);
 int openmc_material_set_volume(int32_t index, double volume);
 int openmc_material_filter_get_bins(
   int32_t index, const int32_t** bins, size_t* n);
diff --git a/include/openmc/material.h b/include/openmc/material.h
index 709d20573..b13cfe330 100644
--- a/include/openmc/material.h
+++ b/include/openmc/material.h
@@ -158,6 +158,10 @@ public:
   double density_;                      //!< Total atom density in [atom/b-cm]
   double density_gpcc_;                 //!< Total atom density in [g/cm^3]
   double volume_ {-1.0};                //!< Volume in [cm^3]
+
+  double temperature_ {-1};             //!< Default temperature for material
+                                        //!< A negative indicates no default
+                                        //!< temperature was specified.
   bool fissionable_ {
     false};                 //!< Does this material contain fissionable nuclides
   bool depletable_ {false}; //!< Is the material depletable?
@@ -195,10 +199,6 @@ private:
   // Private data members
   gsl::index index_;
 
-  //! \brief Default temperature for cells containing this material.
-  //!
-  //! A negative value indicates no default temperature was specified.
-  double temperature_ {-1};
 };
 
 //==============================================================================
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index 2d0261960..b96ab91fd 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -32,7 +32,7 @@ extern "C" bool cmfd_run;  //!< is a CMFD run?
 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)
+extern "C" bool event_based;    //!< use event-based mode (instead of history-based)
 extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular?
 extern bool material_cell_offsets; //!< create material cells offsets?
 extern "C" bool output_summary;    //!< write summary.h5?
diff --git a/openmc/executor.py b/openmc/executor.py
index fa55e1366..ef1bc25dc 100644
--- a/openmc/executor.py
+++ b/openmc/executor.py
@@ -7,8 +7,8 @@ import openmc
 
 def process_CLI_arguments(volume=False, geometry_debug=False, particles=None,
                           plot=False, restart_file=None, threads=None,
-                          tracks=False, event_based=False,
-                          openmc_exec='openmc', mpi_args=None):
+                          tracks=False, event_based=None, openmc_exec='openmc',
+                          mpi_args=None):
     """Run an OpenMC simulation.
 
     Parameters
@@ -30,8 +30,9 @@ def process_CLI_arguments(volume=False, geometry_debug=False, particles=None,
         :envvar:`OMP_NUM_THREADS` environment variable).
     tracks : bool, optional
         Write tracks for all particles. Defaults to False.
-    event_based : bool, optional
-        Turns on event-based parallelism, instead of default history-based
+    event_based : None or bool, optional
+        Turns on event-based parallelism if True. If None, the value in
+        the Settings will be used.
     openmc_exec : str, optional
         Path to OpenMC executable. Defaults to 'openmc'.
     mpi_args : list of str, optional
@@ -61,8 +62,9 @@ def process_CLI_arguments(volume=False, geometry_debug=False, particles=None,
     if geometry_debug:
         args.append('-g')
 
-    if event_based:
-        args.append('-e')
+    if event_based is not None:
+        if event_based:
+            args.append('-e')
 
     if isinstance(restart_file, str):
         args += ['-r', restart_file]
@@ -70,6 +72,9 @@ def process_CLI_arguments(volume=False, geometry_debug=False, particles=None,
     if tracks:
         args.append('-t')
 
+    if plot:
+        args.append('-p')
+
     if mpi_args is not None:
         args = mpi_args + args
 
diff --git a/openmc/lib/__init__.py b/openmc/lib/__init__.py
index f1b9c291f..c14b0d9c2 100644
--- a/openmc/lib/__init__.py
+++ b/openmc/lib/__init__.py
@@ -61,4 +61,6 @@ from .math import *
 from .plot import *
 
 # Flag to denote whether or not openmc.lib.init has been called
+# TODO: Establish and use a flag in the C++ code to represent the status of the
+# openmc_init and openmc_finalize methods
 is_initialized = False
diff --git a/openmc/lib/material.py b/openmc/lib/material.py
index fde197d3d..c1800cb46 100644
--- a/openmc/lib/material.py
+++ b/openmc/lib/material.py
@@ -57,6 +57,9 @@ _dll.openmc_material_get_name.errcheck = _error_handler
 _dll.openmc_material_set_name.argtypes = [c_int32, c_char_p]
 _dll.openmc_material_set_name.restype = c_int
 _dll.openmc_material_set_name.errcheck = _error_handler
+_dll.openmc_material_set_temperature.argtypes = [c_int32, c_double]
+_dll.openmc_material_set_temperature.restype = c_int
+_dll.openmc_material_set_temperature.errcheck = _error_handler
 _dll.openmc_material_set_volume.argtypes = [c_int32, c_double]
 _dll.openmc_material_set_volume.restype = c_int
 _dll.openmc_material_set_volume.errcheck = _error_handler
@@ -165,6 +168,10 @@ class Material(_FortranObjectWithID):
             return None
         return volume.value
 
+    @temperature.setter
+    def temperature(self, temperature):
+        _dll.openmc_material_set_temperature(self._index, temperature)
+
     @volume.setter
     def volume(self, volume):
         _dll.openmc_material_set_volume(self._index, volume)
diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py
index 1eab4aa8a..2e9dd18df 100644
--- a/openmc/lib/settings.py
+++ b/openmc/lib/settings.py
@@ -34,6 +34,7 @@ class _Settings:
     restart_run = _DLLGlobal(c_bool, 'restart_run')
     run_CE = _DLLGlobal(c_bool, 'run_CE')
     verbosity = _DLLGlobal(c_int, 'verbosity')
+    event_based = _DLLGlobal(c_bool, 'event_based')
 
     @property
     def run_mode(self):
diff --git a/openmc/model/model.py b/openmc/model/model.py
index 44479d7e7..465ad3064 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -73,7 +73,7 @@ class Model:
     """
 
     def __init__(self, geometry=None, materials=None, settings=None,
-                 tallies=None, plots=None, intracomm=None):
+                 tallies=None, plots=None):
         self.geometry = openmc.Geometry()
         self.materials = openmc.Materials()
         self.settings = openmc.Settings()
@@ -145,8 +145,6 @@ class Model:
 
     @property
     def is_initialized(self):
-        # TODO: Replace openmc.lib.is_initialized with a direct ctypes access
-        # to simulation::initialized
         return openmc.lib.is_initialized
 
     @geometry.setter
@@ -222,7 +220,7 @@ class Model:
         return cls(geometry, materials, settings)
 
     def init_lib(self, threads=None, geometry_debug=False, restart_file=None,
-                 tracks=False, output=True, event_based=False):
+                 tracks=False, output=True, event_based=None):
         """Initializes the model in memory via the C API
 
         .. versionadded:: 0.13.0
@@ -243,8 +241,9 @@ class Model:
             Write tracks for all particles. Defaults to False.
         output : bool
             Capture OpenMC output from standard out
-        event_based : bool, optional
-            Turns on event-based parallelism, instead of default history-based
+        event_based : None or bool, optional
+            Turns on event-based parallelism if True. If None, the value in
+            the Settings will be used.
         """
 
         # TODO: right now the only way to set most of the above parameters via
@@ -456,7 +455,7 @@ class Model:
 
     def run(self, particles=None, threads=None, geometry_debug=False,
             restart_file=None, tracks=False, output=True, cwd='.',
-            openmc_exec='openmc', mpi_args=None, event_based=False):
+            openmc_exec='openmc', mpi_args=None, event_based=None):
         """Runs OpenMC. If the C API has been initialized, then the C API is
         used, otherwise, this method creates the XML files and runs OpenMC via
         a system call. In both cases this method returns the path to the last
@@ -494,8 +493,9 @@ class Model:
         mpi_args : list of str, optional
             MPI execute command and any additional MPI arguments to pass,
             e.g. ['mpiexec', '-n', '8'].
-        event_based : bool, optional
-            Turns on event-based parallelism, instead of default history-based
+        event_based : None or bool, optional
+            Turns on event-based parallelism if True. If None, the value in
+            the Settings will be used.
 
         Returns
         -------
@@ -523,21 +523,21 @@ class Model:
                         msg = f"{arg_name} must be set via Model.is_initialized(...)"
                         raise ValueError(msg)
 
+                init_particles = openmc.lib.settings.particles
                 if particles is not None:
-                    init_particles = openmc.lib.settings.particles
                     if isinstance(particles, Integral) and particles > 0:
                         openmc.lib.settings.particles = particles
 
-                # Event-based can be set at init-time or on a case-basis.
-                # Handle the argument here.
-                # TODO This will be dealt with in a future change to the C API
+                init_event_based = openmc.lib.settings.event_based
+                if event_based is not None:
+                    openmc.lib.settings.event_based = event_based
 
                 # Then run using the C API
                 openmc.lib.run(output)
 
                 # Reset changes for the openmc.run kwargs handling
-                if particles is not None:
-                    openmc.lib.settings.particles = init_particles
+                openmc.lib.settings.particles = init_particles
+                openmc.lib.settings.event_based = init_event_based
 
             else:
                 # Then run via the command line
@@ -697,10 +697,7 @@ class Model:
         if obj_type == 'cell' and attrib_name == 'volume':
             raise NotImplementedError(
                 'Setting a Cell volume is not supported!')
-        # Same with setting temperatures, TODO: update C API for this
-        if obj_type == 'material' and attrib_name == 'temperature':
-            raise NotImplementedError(
-                'Setting a Material temperature is not yet supported!')
+
         # And some items just dont make sense
         if obj_type == 'cell' and attrib_name == 'density':
             raise ValueError('Cannot set a Cell density!')
@@ -751,10 +748,10 @@ class Model:
             # Next lets keep what is in C API memory up to date as well
             if self.is_initialized:
                 lib_obj = obj_by_id[id_]
-                if attrib_name == 'temperature':
-                    lib_obj.set_temperature(value)
-                elif attrib_name == 'density':
+                if attrib_name == 'density':
                     lib_obj.set_density(value, density_units)
+                elif attrib_name == 'temperature' and obj_type == 'cell':
+                    lib_obj.set_temperature(value)
                 else:
                     setattr(lib_obj, attrib_name, value)
 
@@ -800,7 +797,8 @@ class Model:
 
         """
 
-        self._change_py_lib_attribs(names_or_ids, vector, 'cell', 'translation')
+        self._change_py_lib_attribs(names_or_ids, vector, 'cell',
+                                    'translation')
 
     def update_densities(self, names_or_ids, density, density_units='atom/b-cm'):
         """Update the density of a given set of materials to a new value
@@ -822,8 +820,8 @@ class Model:
 
         """
 
-        self._change_py_lib_attribs(names_or_ids, density, 'material', 'density',
-                                  density_units)
+        self._change_py_lib_attribs(names_or_ids, density, 'material',
+                                    'density', density_units)
 
     def update_cell_temperatures(self, names_or_ids, temperature):
         """Update the temperature of a set of cells to the given value
@@ -844,7 +842,7 @@ class Model:
         """
 
         self._change_py_lib_attribs(names_or_ids, temperature, 'cell',
-                                  'temperature')
+                                    'temperature')
 
     def update_material_temperatures(self, names_or_ids, temperature):
         """Update the temperature of a set of materials to the given value
@@ -865,7 +863,7 @@ class Model:
         """
 
         self._change_py_lib_attribs(names_or_ids, temperature, 'material',
-                                  'temperature')
+                                   'temperature')
 
     def update_material_volumes(self, names_or_ids, volume):
         """Update the volume of a set of materials to the given value
diff --git a/src/material.cpp b/src/material.cpp
index ca29ce740..fb4b2675e 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -1480,6 +1480,18 @@ extern "C" int openmc_material_set_name(int32_t index, const char* name)
   return 0;
 }
 
+extern "C" int openmc_material_set_temperature(int32_t index, double temperature)
+{
+  if (index >= 0 && index < model::materials.size()) {
+    auto& m {model::materials[index]};
+    m->temperature_ = temperature;
+    return 0;
+  } else {
+    set_errmsg("Index in materials array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+}
+
 extern "C" int openmc_material_set_volume(int32_t index, double volume)
 {
   if (index >= 0 && index < model::materials.size()) {
diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index 079b4f9eb..c21f85633 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -207,6 +207,11 @@ def test_material(lib_init):
     m.name = "Not hot borated water"
     assert m.name == "Not hot borated water"
 
+    assert m.temperature == 293.6
+    m.temperature = 400.
+    assert m.temperature == 400.
+    m.temperature == 293.6
+
 
 def test_properties_density(lib_init):
     m = openmc.lib.materials[1]
@@ -263,6 +268,7 @@ def test_settings(lib_init):
     assert settings.generations_per_batch == 1
     assert settings.particles == 100
     assert settings.seed == 1
+    assert settings.event_based is False
     settings.seed = 11
 
 
diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py
index 2760eadb4..092ccfc3c 100644
--- a/tests/unit_tests/test_model.py
+++ b/tests/unit_tests/test_model.py
@@ -5,7 +5,6 @@ from pathlib import Path
 from shutil import which
 import openmc
 import openmc.lib
-from openmc.mpi import DummyCommunicator
 
 
 @pytest.fixture(scope='function')
@@ -128,12 +127,14 @@ def test_init(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     assert test_model._cells_by_id == {}
     assert test_model._cells_by_name == {}
     assert test_model.is_initialized is False
-    assert isinstance(test_model.intracomm, DummyCommunicator)
+    assert hasattr(test_model.intracomm, 'rank')
 
     # Now check proper init of an actual model. Assume no interference between
     # parameters and so we can apply them all at once instead of testing one
     # parameter initialization at a time
-    test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
+    test_model = openmc.Model(geom, mats, settings, tals, plots)
+    if mpi_intracomm is not None:
+        test_model.intracomm = mpi_intracomm
     assert test_model.geometry is geom
     assert test_model.materials is mats
     assert test_model.settings is settings
@@ -154,17 +155,10 @@ def test_init(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
         '': [geom.root_universe.cells[3], geom.root_universe.cells[4]],
         'inf fuel': [geom.root_universe.cells[2].fill.cells[1]]}
     assert test_model.is_initialized is False
-    if mpi_intracomm is None:
-        assert isinstance(test_model.intracomm, DummyCommunicator)
-    else:
-        assert test_model.intracomm == mpi_intracomm
 
     # Finally test the parameter type checking by passing bad types and
     # obtaining the right exception types
-    if mpi_intracomm is not None:
-        def_params = [geom, mats, settings, tals, plots, mpi_intracomm]
-    else:
-        def_params = [geom, mats, settings, tals, plots]
+    def_params = [geom, mats, settings, tals, plots]
     for i in range(len(def_params)):
         args = def_params.copy()
         # Try an integer, as that is a bad type for all arguments
@@ -223,13 +217,14 @@ def test_from_xml(run_in_tmpdir, pin_model_attributes):
              test_model.geometry.root_universe.cells[4]],
         'inf fuel': [test_model.geometry.root_universe.cells[2].fill.cells[1]]}
     assert test_model.is_initialized is False
-    assert isinstance(test_model.intracomm, DummyCommunicator)
 
 
 def test_init_finalize_lib(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     # We are going to init and then make sure data is loaded
     mats, geom, settings, tals, plots, _, _ = pin_model_attributes
-    test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
+    test_model = openmc.Model(geom, mats, settings, tals, plots)
+    if mpi_intracomm is not None:
+        test_model.intracomm = mpi_intracomm
     test_model.init_lib(output=False)
 
     # First check that the API is advertised as initialized
@@ -257,6 +252,8 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm):
     # Create PWR pin cell model and write XML files
     openmc.reset_auto_ids()
     model = openmc.examples.pwr_pin_cell()
+    if mpi_intracomm is not None:
+        model.intracomm = mpi_intracomm
     model.init_lib(output=False)
 
     # Change fuel temperature and density and export properties
@@ -296,7 +293,9 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm):
 
 def test_run(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     mats, geom, settings, tals, plots, _, _ = pin_model_attributes
-    test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
+    test_model = openmc.Model(geom, mats, settings, tals, plots)
+    if mpi_intracomm is not None:
+        test_model.intracomm = mpi_intracomm
 
     # This case will run by getting the k-eff and tallies for command-line and
     # C API execution modes and ensuring they give the same result.
@@ -334,7 +333,9 @@ def test_run(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
 
 def test_plots(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     mats, geom, settings, tals, plots, _, _ = pin_model_attributes
-    test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
+    test_model = openmc.Model(geom, mats, settings, tals, plots)
+    if mpi_intracomm is not None:
+        test_model.intracomm = mpi_intracomm
 
     # This test cannot check the correctness of the plot, but it can
     # check that a plot was made and that the expected ppm and png files are
@@ -367,7 +368,9 @@ def test_plots(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
 
 def test_py_lib_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     mats, geom, settings, tals, plots, _, _ = pin_model_attributes
-    test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
+    test_model = openmc.Model(geom, mats, settings, tals, plots)
+    if mpi_intracomm is not None:
+        test_model.intracomm = mpi_intracomm
 
     test_model.init_lib(output=False)
 
@@ -439,19 +442,13 @@ def test_py_lib_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
 
     # Now lets do the material temperature updates.
     # Check initial conditions
-    with pytest.raises(NotImplementedError):
-        test_model.update_material_temperatures(['UO2'], 600.)
-    # TODO: When C API material.set_temperature is implemented, uncomment below
-    # assert abs(openmc.lib.materials[1].temperature - 293.6) < 1e-13
-    # # The temperature on the material will be None because its just the
-    # # default assert test_model.materials[0].temperature is None
-    # # Change the temperature
-    # test_model.update_material_temperatures(['UO2'], 600.)
-    # assert abs(openmc.lib.materials[1].temperature - 600.) < 1e-13
-    # assert abs(test_model.materials[0].temperature - 600.) < 1e-13
+    assert abs(openmc.lib.materials[1].temperature - 293.6) < 1e-13
+    # Change the temperature
+    test_model.update_material_temperatures(['UO2'], 600.)
+    assert abs(openmc.lib.materials[1].temperature - 600.) < 1e-13
+    assert abs(test_model.materials[0].temperature - 600.) < 1e-13
 
     # And finally material volume
-    # import pdb; pdb.set_trace()
     assert abs(openmc.lib.materials[1].volume - 0.4831931368640985) < 1e-13
     # The temperature on the material will be None because its just the default
     assert abs(test_model.materials[0].volume - 0.4831931368640985) < 1e-13
@@ -468,7 +465,9 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
         pin_model_attributes
     with open('test_chain.xml', 'w') as f:
         f.write(chain_file_xml)
-    test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
+    test_model = openmc.Model(geom, mats, settings, tals, plots)
+    if mpi_intracomm is not None:
+        test_model.intracomm = mpi_intracomm
 
     initial_mat = mats[0].clone()
     initial_u = initial_mat.get_nuclide_atom_densities()['U235'][1]
@@ -518,7 +517,9 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
 def test_calc_volumes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     mats, geom, settings, tals, plots, _, _ = pin_model_attributes
 
-    test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
+    test_model = openmc.Model(geom, mats, settings, tals, plots)
+    if mpi_intracomm is not None:
+        test_model.intracomm = mpi_intracomm
 
     # With no vol calcs, it should fail
     with pytest.raises(ValueError):