resolve conflicts

openmc/deplete/__init__.py

@@ -17,6 +17,7 @@ from .stepresult import *
 from .results import *
 from .integrators import *
 from .transfer_rates import *
+from .batchwise import *
 from . import abc
 from . import cram
 from . import helpers

openmc/deplete/abc.py

@@ -18,14 +18,16 @@ from warnings import warn
 from numpy import nonzero, empty, asarray
 from uncertainties import ufloat
 
-from openmc.checkvalue import check_type, check_greater_than, PathLike
+from openmc.checkvalue import checkvalue, check_type, check_greater_than, PathLike
 from openmc.mpi import comm
 from .stepresult import StepResult
 from .chain import Chain
 from .results import Results
 from .pool import deplete
 from .transfer_rates import TransferRates
-
+from openmc import Material, Cell, Universe
+from .batchwise import (BatchwiseCellGeometrical, BatchwiseCellTemperature,
+    BatchwiseMaterialRefuel)
 
 __all__ = [
     "OperatorResult", "TransportOperator",
@@ -636,6 +638,7 @@ class Integrator(ABC):
         self.source_rates = asarray(source_rates)
 
         self.transfer_rates = None
+        self.batchwise = None
 
         if isinstance(solver, str):
             # Delay importing of cram module, which requires this file
@@ -763,6 +766,14 @@ class Integrator(ABC):
         return (self.operator.prev_res[-1].time[-1],
                 len(self.operator.prev_res) - 1)
 
+    def _get_bos_from_batchwise(self, step_index, bos_conc):
+        """Get BOS from criticality batch-wise control
+        """
+        x = deepcopy(bos_conc)
+        # Get new vector after keff criticality control
+        x, root = self.batchwise.search_for_keff(x, step_index)
+        return x, root
+
     def integrate(
             self,
             final_step: bool = True,
@@ -797,6 +808,11 @@ class Integrator(ABC):
 
                 # Solve transport equation (or obtain result from restart)
                 if i > 0 or self.operator.prev_res is None:
+                    # Update geometry/material according to batchwise definition
+                    if self.batchwise and source_rate != 0.0:
+                        n, root = self._get_bos_from_batchwise(i, n)
+                    else:
+                        root = None
                     n, res = self._get_bos_data_from_operator(i, source_rate, n)
                 else:
                     n, res = self._get_bos_data_from_restart(i, source_rate, n)
@@ -810,9 +826,8 @@ class Integrator(ABC):
 
                 # Remove actual EOS concentration for next step
                 n = n_list.pop()
-
                 StepResult.save(self.operator, n_list, res_list, [t, t + dt],
-                                source_rate, self._i_res + i, proc_time, path)
+                                source_rate, self._i_res + i, proc_time, root, path)
 
                 t += dt
 
@@ -822,9 +837,13 @@ class Integrator(ABC):
             # solve)
             if output and final_step and comm.rank == 0:
                 print(f"[openmc.deplete] t={t} (final operator evaluation)")
+            if self.batchwise and source_rate != 0.0:
+                n, root = self._get_bos_from_batchwise(i+1, n)
+            else:
+                root = None
             res_list = [self.operator(n, source_rate if final_step else 0.0)]
             StepResult.save(self.operator, [n], res_list, [t, t],
-                         source_rate, self._i_res + len(self), proc_time)
+                         source_rate, self._i_res + len(self), proc_time, root)
             self.operator.write_bos_data(len(self) + self._i_res)
 
         self.operator.finalize()
@@ -857,6 +876,30 @@ class Integrator(ABC):
         self.transfer_rates.set_transfer_rate(material, components, transfer_rate,
                                       transfer_rate_units, destination_material)
 
+    def add_batchwise(self, obj, attr, **kwargs):
+        """Add batchwise operation to integrator scheme.
+
+        Parameters
+        ----------
+        attr : str
+            Type of batchwise operation to add. `Trans` stands for geometrical
+            translation, `refuel` for material refueling and `dilute` for material
+            dilute.
+        **kwargs
+            keyword arguments that are passed to the batchwise class.
+
+        """
+        check_value('attribute', attr, ('translation', 'rotation', 'temperature', 'refuel'))
+        if attr in ('translation', 'rotation'):
+            batchwise = BatchwiseCellGeometrical
+        elif attr == 'temperature':
+            batchwise = BatchwiseCellTemperature
+        elif attr == 'refuel':
+            batchwise = BatchwiseMaterialRefuel
+
+        self.batchwise = batchwise.from_params(obj, attr, self.operator,
+                                           self.operator.model, **kwargs)
+
 @add_params
 class SIIntegrator(Integrator):
     r"""Abstract class for the Stochastic Implicit Euler integrators

openmc/deplete/stepresult.py

@@ -74,7 +74,7 @@ class StepResult:
         self.mat_to_hdf5_ind = None
 
         self.data = None
-
+        self.batchwise = None
     def __repr__(self):
         t = self.time[0]
         dt = self.time[1] - self.time[0]
@@ -349,6 +349,10 @@ class StepResult:
             "depletion time", (1,), maxshape=(None,),
             dtype="float64")
 
+        handle.create_dataset(
+            "batchwise_root", (1,), maxshape=(None,),
+            dtype="float64")
+
     def _to_hdf5(self, handle, index, parallel=False):
         """Converts results object into an hdf5 object.
 
@@ -379,6 +383,7 @@ class StepResult:
         time_dset = handle["/time"]
         source_rate_dset = handle["/source_rate"]
         proc_time_dset = handle["/depletion time"]
+        root_dset = handle["/batchwise_root"]
 
         # Get number of results stored
         number_shape = list(number_dset.shape)
@@ -412,6 +417,10 @@ class StepResult:
             proc_shape[0] = new_shape
             proc_time_dset.resize(proc_shape)
 
+            root_shape = list(root_dset.shape)
+            root_shape[0] = new_shape
+            root_dset.resize(root_shape)
+
         # If nothing to write, just return
         if len(self.index_mat) == 0:
             return
@@ -435,6 +444,7 @@ class StepResult:
                 proc_time_dset[index] = (
                     self.proc_time / (comm.size * self.n_hdf5_mats)
                 )
+            root_dset[index] = self.batchwise
 
     @classmethod
     def from_hdf5(cls, handle, step):
@@ -458,11 +468,13 @@ class StepResult:
         else:
             # Older versions used "power" instead of "source_rate"
             source_rate_dset = handle["/power"]
+        root_dset = handle["/batchwise_root"]
 
         results.data = number_dset[step, :, :, :]
         results.k = eigenvalues_dset[step, :]
         results.time = time_dset[step, :]
         results.source_rate = source_rate_dset[step, 0]
+        results.batchwise = root_dset[step]
 
         if "depletion time" in handle:
             proc_time_dset = handle["/depletion time"]
@@ -509,7 +521,7 @@ class StepResult:
 
     @staticmethod
     def save(op, x, op_results, t, source_rate, step_ind, proc_time=None,
-             path: PathLike = "depletion_results.h5"):
+             root=None, path: PathLike = "depletion_results.h5"):
         """Creates and writes depletion results to disk
 
         Parameters
@@ -564,6 +576,7 @@ class StepResult:
         results.proc_time = proc_time
         if results.proc_time is not None:
             results.proc_time = comm.reduce(proc_time, op=MPI.SUM)
+        results.batchwise = root
 
         if not Path(path).is_file():
             Path(path).parent.mkdir(parents=True, exist_ok=True)

openmc/search.py

@@ -12,7 +12,7 @@ _SCALAR_BRACKETED_METHODS = ['brentq', 'brenth', 'ridder', 'bisect']
 
 
 def _search_keff(guess, target, model_builder, model_args, print_iterations,
-                 run_args, guesses, results):
+                 run_args, guesses, results, run_in_memory):
     """Function which will actually create our model, run the calculation, and
     obtain the result. This function will be passed to the root finding
     algorithm
@@ -51,7 +51,11 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations,
     model = model_builder(guess, **model_args)
 
     # Run the model and obtain keff
-    sp_filepath = model.run(**run_args)
+    if run_in_memory:
+        openmc.lib.run(**run_args)
+        sp_filepath = f'statepoint.{model.settings.batches}.h5'
+    else:
+        sp_filepath = model.run(**run_args)
     with openmc.StatePoint(sp_filepath) as sp:
         keff = sp.keff
 
@@ -70,7 +74,7 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations,
 def search_for_keff(model_builder, initial_guess=None, target=1.0,
                     bracket=None, model_args=None, tol=None,
                     bracketed_method='bisect', print_iterations=False,
-                    run_args=None, **kwargs):
+                    run_args=None, run_in_memory=False, **kwargs):
     """Function to perform a keff search by modifying a model parametrized by a
     single independent variable.
 
@@ -193,12 +197,26 @@ def search_for_keff(model_builder, initial_guess=None, target=1.0,
 
     # Add information to be passed to the searching function
     args['args'] = (target, model_builder, model_args, print_iterations,
-                    run_args, guesses, results)
+                    run_args, guesses, results, run_in_memory)
 
     # Create a new dictionary with the arguments from args and kwargs
     args.update(kwargs)
 
     # Perform the search
-    zero_value = root_finder(**args)
+    if not run_in_memory:
+        zero_value = root_finder(**args)
+        return zero_value
 
-    return zero_value, guesses, results
+    else:
+        try:
+            zero_value, root_res = root_finder(**args, full_output=True, disp=False)
+            if root_res.converged:
+                return zero_value, guesses, results
+
+            else:
+                print(f'WARNING: {root_res.flag}')
+                return guesses, results
+        # In case the root finder is not successful
+        except Exception as e:
+            print(f'WARNING: {e}')
+            return guesses, results

tests/unit_tests/test_deplete_batchwise.py

@@ -0,0 +1,135 @@
+""" Tests for Batchwise class """
+
+from pathlib import Path
+from math import exp
+
+import pytest
+import numpy as np
+
+import openmc
+import openmc.lib
+from openmc.deplete import CoupledOperator
+from openmc.deplete import (BatchwiseCellGeometrical, BatchwiseCellTemperature,
+    BatchwiseMaterialRefuel)
+
+CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
+
+@pytest.fixture
+def model():
+    f = openmc.Material(name="f")
+    f.add_element("U", 1, percent_type="ao", enrichment=4.25)
+    f.add_element("O", 2)
+    f.set_density("g/cc", 10.4)
+    f.temperature = 293.15
+    w = openmc.Material(name="w")
+    w.add_element("O", 1)
+    w.add_element("H", 2)
+    w.set_density("g/cc", 1.0)
+    w.depletable = True
+    h = openmc.Material(name='h')
+    h.set_density('g/cm3', 0.001598)
+    h.add_element('He', 2.4044e-4)
+    radii = [0.42, 0.45]
+    f.volume = np.pi * radii[0] ** 2
+    w.volume = np.pi * (radii[1]**2 - radii[0]**2)
+    materials = openmc.Materials([f, w, h])
+    surf_wh = openmc.ZPlane(z0=0)
+    surf_f = openmc.Sphere(r=radii[0])
+    surf_w = openmc.Sphere(r=radii[1], boundary_type='vacuum')
+    cell_w = openmc.Cell(fill=w, region=-surf_wh)
+    cell_h = openmc.Cell(fill=h, region=+surf_wh)
+    uni_wh = openmc.Universe(cells=(cell_w, cell_h))
+    cell_f = openmc.Cell(name='f',fill=f, region=-surf_f)
+    cell_wh = openmc.Cell(name='wh',fill=uni_wh, region=+surf_f & -surf_w)
+    geometry = openmc.Geometry([cell_f, cell_wh])
+    settings = openmc.Settings()
+    settings.particles = 1000
+    settings.inactive = 10
+    settings.batches = 50
+    return openmc.Model(geometry, materials, settings)
+
+@pytest.mark.parametrize("object_name, attribute", [
+    ('wh', 'translation'),
+    ('wh', 'rotation'),
+    ('f', 'temperature'),
+    ('f', 'refuel'),
+    ])
+def test_attributes(model, object_name, attribute):
+    """
+    Test classes attributes are set correctly
+    """
+    op = CoupledOperator(model, CHAIN_PATH)
+
+    bracket = [-1,1]
+    bracket_limit = [-10,10]
+    axis = 2
+    mat_vector = {'U235':0.1, 'U238':0.9}
+    if attribute in ('translation','rotation'):
+        bw = BatchwiseCellGeometrical(object_name, attribute, op, model, axis, bracket, bracket_limit)
+        assert bw.cell_materials == [cell.fill for cell in  \
+                model.geometry.get_cells_by_name(object_name)[0].fill.cells.values() \
+                if cell.fill.depletable]
+        assert bw.attrib_name == attribute
+        assert bw.axis == axis
+
+    elif attribute == 'temperature':
+        bw = BatchwiseCellTemperature(object_name, op, model, bracket, bracket_limit)
+    elif attribute == 'refuel':
+        bw = BatchwiseMaterialRefuel(object_name, op, model, mat_vector, bracket, bracket_limit)
+        assert bw.mat_vector == mat_vector
+
+    assert bw.bracket == bracket
+    assert bw.bracket_limit == bracket_limit
+    assert bw.burn_mats == op.burnable_mats
+    assert bw.local_mats == op.local_mats
+
+@pytest.mark.parametrize("attribute", [
+    ('translation'),
+    ('rotation')
+    ])
+def test_cell_methods(run_in_tmpdir, model, attribute):
+    """
+    Test cell base class internal method
+    """
+    bracket = [-1,1]
+    bracket_limit = [-10,10]
+    axis = 2
+    op = CoupledOperator(model, CHAIN_PATH)
+    integrator = openmc.deplete.PredictorIntegrator(
+        op, [1,1], 0.0, timestep_units = 'd')
+    integrator.add_batchwise('wh', attribute, axis=axis, bracket=bracket,
+                             bracket_limit=bracket_limit)
+
+    model.export_to_xml()
+    openmc.lib.init()
+    integrator.batchwise._set_cell_attrib(0)
+    assert integrator.batchwise._get_cell_attrib() == 0
+    vol = integrator.batchwise._calculate_volumes()
+    for cell in integrator.batchwise.cell_materials:
+        assert vol[str(cell.id)] == pytest.approx([mat.volume for mat in model.materials if mat.id == cell.id])
+    openmc.lib.finalize()
+
+def test_update_materials_methods(run_in_tmpdir):
+    """
+    It's the same methods for all classe so one class is enough
+    """
+    bracket = [-1,1]
+    bracket_limit = [-10,10]
+    axis = 2
+    op = CoupledOperator(model, CHAIN_PATH)
+    integrator = openmc.deplete.PredictorIntegrator(
+        op, [1,1], 0.0, timestep_units = 'd')
+    integrator.add_batchwise('wh', 'translation', axis=axis, bracket=bracket,
+                             bracket_limit=bracket_limit)
+
+    model.export_to_xml()
+    openmc.lib.init()
+    x = op.number.number
+    #Double  number of atoms of U238 in fuel
+    mat_index = op.number.index_mat['1']
+    nuc_index = op.number.index_nuc['U238']
+    vol = op.number.get_mat_volume('1')
+    op.number.number[mat_index][nuc_idx] *= 2
+    integrator.batchwise._update_volumes(op.number.number)
+    
+    vol = integrator.batchwise._calculate_volumes()