Repository: pmeal/openpnm
Issue: Problem with 'throat.partition_coef.all' not being set as iterative_prop nor updating when regenerate_models is called manually

I am writing a model where the partition coefficient changes with concentration. So it's important that 'throat.partition_coef.all' updates as 'throat.diffusive_conductance' (based on multiphase_diffusion) depends on this property to know partition coefficient. 

Changed files:
scripts/example_mass_partitioning.py, openpnm/utils/misc.py, openpnm/utils/__init__.py, openpnm/core/ModelsMixin.py, openpnm/algorithms/ReactiveTransport.py, tests/unit/models/physics/DiffusiveConductanceTest.py, tests/unit/algorithms/GenericTransportTest.py, openpnm/phases/MultiPhase.py, tests/unit/phases/MultiPhaseTest.py, openpnm/models/physics/diffusive_conductance.py, scripts/example_multiphase_diffusion.py, openpnm/phases/Air.py

diff --git a/openpnm/algorithms/ReactiveTransport.py b/openpnm/algorithms/ReactiveTransport.py
index 4f50c065f..5b495bc99 100644
--- a/openpnm/algorithms/ReactiveTransport.py
+++ b/openpnm/algorithms/ReactiveTransport.py
@@ -368,14 +368,16 @@ class ReactiveTransport(GenericTransport):
 
     def _get_iterative_props(self):
         r"""
-        Find and return properties that need to be iterated while running the
-        algorithm
+        Find and return properties that need to be iterated while running
+        the algorithm.
 
         Notes
         -----
-        This method was moved from ReactiveTransport class to GenericTransport
-        because source terms are not necessarily the only properties that need
-        iteration during an algorithm (ex. concentration-dependent conductance)
+        This method was moved from ReactiveTransport class to
+        GenericTransport because source terms are not necessarily the only
+        properties that need iteration during an algorithm (ex.
+        concentration-dependent conductance)
+
         """
         import networkx as nx
         phase = self.project.phases(self.settings['phase'])
diff --git a/openpnm/core/ModelsMixin.py b/openpnm/core/ModelsMixin.py
index 37b18c7bd..649d70d60 100644
--- a/openpnm/core/ModelsMixin.py
+++ b/openpnm/core/ModelsMixin.py
@@ -1,6 +1,7 @@
 import inspect
 from openpnm.utils import PrintableDict, logging, Workspace
 from openpnm.utils.misc import is_valid_propname
+from openpnm.utils import prettify_logger_message
 logger = logging.getLogger(__name__)
 ws = Workspace()
 
@@ -395,9 +396,9 @@ class ModelsMixin:
         # Only regenerate model if regen_mode is correct
         if self.settings['freeze_models']:
             # Don't run ANY models if freeze_models is set to True
-            logger.warning(
-                f"{prop} was not run since freeze_models is set to True in obj. settings"
-            )
+            msg = (f"{prop} was not run since freeze_models is set to"
+                   " True in object settings.")
+            logger.warning(prettify_logger_message(msg))
         elif regen_mode == 'constant':
             # Only regenerate if data not already in dictionary
             if prop not in self.keys():
@@ -406,8 +407,9 @@ class ModelsMixin:
             try:
                 self[prop] = model(target=self, **kwargs)
             except KeyError as e:
-                logger.error(prop + ' was not run since the following '
-                             + 'property is missing: ' + e.__str__())
+                msg = (f"{prop} was not run since the following property"
+                       f" is missing: {e}")
+                logger.error(prettify_logger_message(msg))
                 self.models[prop]['regen_mode'] = 'deferred'
 
     def remove_model(self, propname=None, mode=['model', 'data']):
diff --git a/openpnm/models/physics/diffusive_conductance.py b/openpnm/models/physics/diffusive_conductance.py
index f2801ee7f..1e1fcb46f 100644
--- a/openpnm/models/physics/diffusive_conductance.py
+++ b/openpnm/models/physics/diffusive_conductance.py
@@ -498,6 +498,7 @@ def multiphase_diffusion(
     throat_diffusivity="throat.diffusivity",
     conduit_lengths="throat.conduit_lengths",
     conduit_shape_factors="throat.poisson_shape_factors",
+    partition_coef_global="throat.partition_coef.all"
 ):
     r"""
     Calculates the diffusive conductance of conduits in network.
@@ -532,18 +533,18 @@ def multiphase_diffusion(
     information.
 
     """
-    network = target.project.network
-    throats = network.map_throats(throats=target.Ts, origin=target)
+    network = target.network
+    throats = target.throats(target=network)
     phase = target.project.find_phase(target)
-    cn = network["throat.conns"][throats]
+    cn = network.conns[throats]
+
     # Getting equivalent areas
-    A1 = network[pore_area][cn[:, 0]]
+    A1, A2 = network[pore_area][cn].T
     At = network[throat_area][throats]
-    A2 = network[pore_area][cn[:, 1]]
     # Getting conduit lengths
-    L1 = network[conduit_lengths + ".pore1"][throats]
-    Lt = network[conduit_lengths + ".throat"][throats]
-    L2 = network[conduit_lengths + ".pore2"][throats]
+    L1 = network[f"{conduit_lengths}.pore1"][throats]
+    Lt = network[f"{conduit_lengths}.throat"][throats]
+    L2 = network[f"{conduit_lengths}.pore2"][throats]
     # Getting shape factors
     try:
         SF1 = phase[conduit_shape_factors + ".pore1"][throats]
@@ -554,18 +555,19 @@ def multiphase_diffusion(
     # Interpolate pore phase property values to throats
     D1, D2 = phase[pore_diffusivity][cn].T
     Dt = phase.interpolate_data(propname=pore_diffusivity)[throats]
+
     # Find g for half of pore 1, throat, and half of pore 2 + apply shape factors
     g1 = (D1 * A1) / L1 * SF1
     g2 = (D2 * A2) / L2 * SF2
     gt = (Dt * At) / Lt * SFt
+
     # Ensure infinite conductance for elements with zero length
-    g1[L1 == 0] = _np.inf
-    g2[L2 == 0] = _np.inf
-    gt[Lt == 0] = _np.inf
-    # Get partition coefficient dictionary key from phase settings
-    partition_coef = phase.settings["partition_coef"]
+    for gi, Li in zip([g1, gt, g2], [L1, Lt, L2]):
+        gi[Li == 0] = _np.inf
+
     # Apply Henry's partitioning coefficient
-    K12 = phase[partition_coef][throats]
+    K12 = phase[partition_coef_global][throats]
     G12 = K12 * (1.0 / g1 + 0.5 / gt + K12 * (1.0 / g2 + 0.5 / gt)) ** (-1)
     G21 = 1.0 / K12 * (1.0 / g2 + 0.5 / gt + 1.0 / K12 * (1.0 / g1 + 0.5 / gt)) ** (-1)
+
     return _np.vstack((G12, G21)).T
diff --git a/openpnm/phases/Air.py b/openpnm/phases/Air.py
index 03665e22c..6b7407524 100644
--- a/openpnm/phases/Air.py
+++ b/openpnm/phases/Air.py
@@ -100,6 +100,7 @@ class Air(GenericPhase):
     Vol. 25, No. 1, January 2004, pp. 21-69
 
     """
+
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
 
diff --git a/openpnm/phases/MultiPhase.py b/openpnm/phases/MultiPhase.py
index 6ea0f3ecd..ec0380208 100644
--- a/openpnm/phases/MultiPhase.py
+++ b/openpnm/phases/MultiPhase.py
@@ -12,62 +12,76 @@ class MultiPhase(GenericPhase):
 
     Parameters
     ----------
-    network : OpenPNM Network object
+    network : GenericNetwork
         The network to which this phase object will be attached.
-
-    project : OpenPNM Project object, optional
-        The Project with which this phase should be associted.  If a
-        ``network`` is given then this is ignored and the Network's project
-        is used.  If a ``network`` is not given then this is mandatory.
-
-    name : string, optional
-        The name of the phase.  This is useful to keep track of the objects
-        throughout the simulation.  The name must be unique to the project.
+    project : Project, optional
+        The Project with which this phase should be associted. If a
+        ``network`` is given then this is ignored and the Network's
+        project is used. If a ``network`` is not given then this is
+        mandatory.
+    name : str, optional
+        The name of the phase. This is useful to keep track of the objects
+        throughout the simulation. The name must be unique to the project.
         If no name is given, one is generated.
 
     Examples
     --------
     >>> import scipy as sp
     >>> import openpnm as op
-    >>> pn = op.network.Cubic(shape=[5, 5, 5])
-    >>> air = op.phases.Air(network=pn, name='air')  # Create two pure phases
-    >>> water = op.phases.Water(network=pn, name='water')
-    >>> m = op.phases.MultiPhase(network=pn, phases=[air, water], name='multi')
-    >>> Ps = pn['pore.coords'][:, 0] < 3  # Pick some pores to be air filled
-    >>> Ts = pn.find_neighbor_throats(pores=Ps)  # Find neighboring throats
-    >>> Ts = pn.tomask(throats=Ts)  # Convert throat indices to mask
-    >>> m.set_occupancy(phase=air, Pvals=Ps, Tvals=Ts)  # Assign occupancies
-    >>> m.set_occupancy(phase=water, Pvals=~Ps, Tvals=~Ts)
-
-    Air and water have uniform viscosity values throughout, so the peak-to-peak
-    distance is 0, while the mixture phase has the viscosity of air in some
-    locations and water in others, hence a heterogenous viscosity array:
-
-    >>> np.ptp(water['pore.viscosity']) == 0
+    >>> from openpnm.phases import Air, Water, MultiPhase
+
+    >>> net = op.network.Cubic(shape=[5, 5, 5])
+    >>> air = Air(network=net, name='air')  # Create two pure phases
+    >>> water = Water(network=net, name='water')
+    >>> mphase = MultiPhase(network=net, phases=[air, water], name='multi')
+    >>> Ps = net['pore.coords'][:, 0] < 3  # Pick some pores to be air filled
+    >>> Ts = net.find_neighbor_throats(pores=Ps)  # Find neighboring throats
+    >>> Ts = net.tomask(throats=Ts)  # Convert throat indices to mask
+    >>> mphase.set_occupancy(phase=air, Pvals=Ps, Tvals=Ts)  # Assign occupancies
+    >>> mphase.set_occupancy(phase=water, Pvals=~Ps, Tvals=~Ts)
+
+    Air and water have uniform viscosity values throughout, so the
+    peak-to-peak distance is 0, while the mixture phase has the viscosity
+    of air in some locations and water in others, hence a heterogenous
+    viscosity array:
+
+    >>> np.ptp(water['pore.viscosity']) == 0.
     True
-    >>> np.ptp(air['pore.viscosity']) == 0
+    >>> np.ptp(air['pore.viscosity']) == 0.
     True
-    >>> np.ptp(m['pore.viscosity']) == 0
+    >>> np.ptp(mphase['pore.viscosity']) == 0.
     False
 
     """
+
     def __init__(self, phases=[], settings={}, **kwargs):
         super().__init__(**kwargs)
-        self.settings.update({'phases': [],
-                              'throat_occupancy': 'manual',
-                              'partition_coef': 'throat.partition_coef.all'})
+        self.settings.update(
+            {
+                'phases': [],
+                'throat_occupancy': 'manual',
+                'partition_coef_prefix': 'throat.partition_coef'
+            }
+        )
         self.settings.update(settings)
 
         self['pore.occupancy.all'] = np.zeros(self.Np, dtype=float)
         self['throat.occupancy.all'] = np.zeros(self.Nt, dtype=float)
 
+        # Pressure/temperature must come from constituent phases
         self.pop('pore.temperature', None)
         self.pop('pore.pressure', None)
-        # Initialize throat partition coefficient values to 1.0
-        partition_coef = self.settings['partition_coef']
-        self[partition_coef] = np.ones(self.Nt)
 
-        # Add supplied phases to the phases list and initialize occupancy to 0
+        # Add a dummy model for assembling the global parition
+        # coefficients array. The dummy model only accepts 'target' and
+        # kwargs. When a binary partition coefficient is added to the
+        # MultiPhase, we manually append the binary partition coefficient
+        # propname as a keyword argument to this dummy model, so that our
+        # dependency handler is notified of this dependency!
+        partition_coef_global = self.settings['partition_coef_prefix'] + ".all"
+        self.add_model(propname=partition_coef_global, model=_dummy)
+
+        # Add supplied phases to phases dict and initialize occupancy to 0
         self.add_phases(phases)
 
         logger.warning('MultiPhase is a beta feature. Functionality may change!')
@@ -79,13 +93,6 @@ class MultiPhase(GenericPhase):
             vals = self.interleave_data(key)
         return vals
 
-    def _update_occupancy(self):
-        # Update "occupancy.all" by summing up all occupancies
-        for elem in ["pore", "throat"]:
-            dict_ = self[f"{elem}.occupancy"]
-            dict_.pop(f"{elem}.occupancy.all")
-            self[f"{elem}.occupancy.all"] = np.sum(list(dict_.values()), axis=0)
-
     def _get_phases(self):
         phases = {self.project[item].name: self.project[item]
                   for item in self.settings['phases']}
@@ -93,148 +100,305 @@ class MultiPhase(GenericPhase):
 
     phases = property(fget=_get_phases)
 
+    def _update_occupancy(self):
+        r"""Updates 'occupancy.all' by summing up all occupancies"""
+        for elem in ["pore", "throat"]:
+            dict_ = self[f"{elem}.occupancy"]
+            dict_.pop(f"{elem}.occupancy.all")
+            self[f"{elem}.occupancy.all"] = np.sum(list(dict_.values()), axis=0)
+
     def add_phases(self, phases):
         r"""
-        Add supplied phases to the MultiPhase object and initialize occupancy to 0.
+        Adds supplied phases to the MultiPhase object and initializes
+        occupancy to 0. It does not return a value.
+
+        Parameters
+        ----------
+        phases : List[GenericPhase]
+
         """
         phases = np.array(phases, ndmin=1)
         for phase in phases:
             if phase.name not in self.settings['phases']:
                 self.settings['phases'].append(phase.name)
-                self[f'pore.occupancy.{phase.name}'] = 0.0
-                self[f'throat.occupancy.{phase.name}'] = 0.0
+                self[f'pore.occupancy.{phase.name}'] = 0.
+                self[f'throat.occupancy.{phase.name}'] = 0.
 
-    def set_binary_partition_coef(self, phases, model,
-                                  propname="throat.partition_coef", **kwargs):
+    def set_binary_partition_coef(self, phases, model, **kwargs):
         r"""
-        Set binary partition coefficient that defines the concentration ratio
-        of phase 1 to phase 2.
+        Sets binary partition coefficient as defined by the interface
+        concentration ratio of phase 1 to phase 2.
+
+        Parameters
+        ----------
+        phases : List[GenericPhase]
+            List of the two phases for which the binary partition
+            coefficient model is being added.
+        model : OpenPNM model
+            Model for calculating the binary partition coefficient.
+        kwargs : dict
+            Keyword arguments to be passed to the ``model``.
+
+        Example
+        -------
+        >>> import openpnm as op
+        >>> from openpnm.phases import Air, Water, MultiPhase
+        >>> from openpnm.models.misc import constant
+
+        >>> net = op.network.Cubic(shape=[5, 5, 5])
+        >>> air = Air(network=net, name='air')  # Create two pure phases
+        >>> water = Water(network=net, name='water')
+        >>> mphase = MultiPhase(network=net, phases=[air, water], name='multi')
+
+        # Now, assign some pores to air and the rest to water
+        >>> Ps = net['pore.coords'][:, 0] < 3  # Pick some pores to be air filled
+        >>> Ts = net.find_neighbor_throats(pores=Ps)  # Find neighboring throats
+        >>> Ts = net.tomask(throats=Ts)  # Convert throat indices to mask
+        >>> mphase.set_occupancy(phase=air, Pvals=Ps, Tvals=Ts)  # Assign occupancies
+        >>> mphase.set_occupancy(phase=water, Pvals=~Ps, Tvals=~Ts)
+
+        # Now, add an interface model for binary partition coefficient
+        >>> mphase.set_binary_partition_coef(phases=[air, water],
+        ...                                  model=constant,
+        ...                                  value=0.5)
+
+        # Now, verify that K12 for interface throats is 0.5
+        >>> Ts_interface = net.find_neighbor_throats(Ps, mode="xor")
+        >>> K12_interface = mphase["throat.partition_coef.all"][Ts_interface]
+        >>> assert K12_interface.mean() == 0.5
+
+        # Finally, verify that K12 for non-interface throats is 1.
+        >>> Ts_rest = ~ np.isin(net.Ts, Ts_interface)
+        >>> K12_rest = mphase["throat.partition_coef.all"][Ts_rest]
+        >>> assert K12_rest.mean() == 1.
+
         """
-        # Raise error if phases is not exactly two-element long
         if np.size(phases) != 2:
             raise Exception("'phases' must contain exactly two elements!")
-        # Add interface model to the MultiPhase object
-        self._add_interface_prop(propname, phases, model, **kwargs)
-        # Get binary partition coefficient
-        interface_prop_key = propname + f".{phases[0].name}:{phases[1].name}"
-        K12 = self[interface_prop_key]
-        # Modify the global partition coefficient vector
-        network = self.project.network
-        cn = network["throat.conns"]
-        occ_phase1 = self[f"pore.occupancy.{phases[0].name}"][cn]
-        occ_phase2 = self[f"pore.occupancy.{phases[1].name}"][cn]
-        partition_coef = self.settings["partition_coef"]
-        K_global = self[partition_coef]
-        mask = (occ_phase1[:, 0] + occ_phase2[:, 1]) == 2.0
-        K_global[mask] = K12[mask]
-        mask = (occ_phase2[:, 0] + occ_phase1[:, 1]) == 2.0
-        K_global[mask] = 1.0 / K12[mask]
+
+        # Add partition coefficient interface model to the MultiPhase
+        propname_prefix = self.settings["partition_coef_prefix"]
+        self._add_interface_prop(propname_prefix, phases, model, **kwargs)
+
+        # Update global partition coef. model's args (for dependency handling)
+        partition_coef_global = self.settings['partition_coef_prefix'] + ".all"
+        K_global_model = self.models[partition_coef_global]
+        # Append binary partition coefficient propname to model args
+        propname = self._format_interface_prop(propname_prefix, phases)
+        key = f"K_{phases[0].name}_{phases[1].name}"
+        K_global_model[key] = propname
+
+        # Regenerate 'throat.parition_coef.all'
+        self.regenerate_models(partition_coef_global)
 
     def _add_interface_prop(self, propname, phases, model, **kwargs):
         r"""
-        Helper method used to add interface models to the MultiPhase object by
-        augmenting ``propname`` --> ``propname.phases[0].name:phases[1].name``
+        Helper method used to add interface models to the ``MultiPhase``
+        object by augmenting.
+
+        Notes
+        -----
+        For convenience and as an OpenPNM convention, the ``propname`` is
+        augmented as outlined by the following example.
+
+        Example
+        -------
+        Let's say the two phases corresponding to the interface model are
+        named: 'air' and 'water', and the interface propname to be added
+        is 'throat.foo'. After augmentation, 'throat.foo.air:water' will
+        be the propname that's stored on the ``MultiPhase`` object.
+
+        Note that because of this convention, the order of the phases that
+        are passed to this method is important.
+
         """
         # Add "throat" keyword to the begining of propname if no identifier is found
         if propname.split(".")[0] not in ["pore", "throat"]:
             propname = f"throat.{propname}"
         # Check propname is throat property
-        if not propname.startswith("throat"):
-            raise Exception("'propname' must be a throat property such as: "
-                            + "'throat.partition_coef'.")
-        propname += f".{phases[0].name}:{phases[1].name}"
+        if propname.startswith("pore"):
+            raise Exception("'propname' must be a throat property")
+        # Add model to Multiphase
+        propname = self._format_interface_prop(propname, phases)
         self.add_model(propname, model, **kwargs)
 
+    def _format_interface_prop(self, propname, phases):
+        r"""
+        Returns formatted interface propname.
+
+        Parameters
+        ----------
+        propname : str
+            Dictionary key of the interface property.
+        phases : List[GenericPhase]
+            List of two phases that together form the interface property.
+
+        Example
+        -------
+        >>> import openpnm as op
+        >>> net = op.network.Cubic(shape=[5, 5, 5])
+        >>> air = op.phases.Air(network=net, name='air')
+        >>> water = op.phases.Water(network=net, name='water')
+        >>> mphase = op.phases.MultiPhase(network=net, phases=[air, water])
+        >>> mphase._format_interface_prop(propname="throat.foo", phases=[air, water])
+        'throat.foo.air:water'
+
+        """
+        prefix = propname
+        suffix = ":".join(phase.name for phase in phases)
+        return f"{prefix}.{suffix}"
+
+    def _get_phases_names(self, formatted_propname):
+        r"""Retrieves phases' names from a formatted propname"""
+        return formatted_propname.split(".")[2].split(":")
+
+    def _assemble_partition_coef_global(self):
+        r"""Updates the global partition coefficient array"""
+        conns = self.network.conns
+        partition_coef_models = []
+        prefix = self.settings["partition_coef_prefix"]
+        K_global = np.ones(self.Nt)
+        partition_coef_global = prefix + ".all"
+
+        # Find all binary partition coefficient models
+        for model_key in self.models.keys():
+            if model_key.startswith(prefix):
+                if model_key != partition_coef_global:
+                    partition_coef_models.append(model_key)
+
+        # Modify the global partition coefficient for each phase pair
+        for propname in partition_coef_models:
+            K12 = self[propname]
+            phases_names = self._get_phases_names(propname)
+            S1 = self[f"pore.occupancy.{phases_names[0]}"][conns]
+            S2 = self[f"pore.occupancy.{phases_names[1]}"][conns]
+            mask = (S1[:, 0] + S2[:, 1]) == 2.
+            K_global[mask] = K12[mask]
+            mask = (S2[:, 0] + S1[:, 1]) == 2.
+            K_global[mask] = 1 / K12[mask]
+
+        return K_global
+
     def interleave_data(self, prop):
         r"""
-        Gathers property values from component phases to build a single array
+        Gathers property values from component phases to build a single
+        array.
 
-        If the requested ``prop`` is not on this MultiPhase, then a search is
-        conducted on all associated phase objects, and values from each
+        If the requested ``prop`` is not on this MultiPhase, then a search
+        is conducted on all associated phase objects, and values from each
         are assembled into a single array.
 
         Parameters
         ----------
-        prop : string
-            The property to be retrieved
+        prop : str
+            The property to be retrieved.
 
         Returns
         -------
-        array : ND-array
+        ndarray
             An array containing the specified property retrieved from each
-            component phase and assembled based on the specified mixing rule
+            component phase and assembled based on the specified mixing
+            rule.
 
         """
-        element = prop.split('.')[0]
+        element = self._parse_element(prop)[0]
         vals = np.zeros([self._count(element=element)], dtype=float)
         # Retrieve property from constituent phases (weight = occupancy)
         try:
             for phase in self.phases.values():
-                vals += phase[prop] * self[element + '.occupancy.' + phase.name]
+                vals += phase[prop] * self[f"{element}.occupancy.{phase.name}"]
         # Otherwise - if not found - retrieve from super class
         except KeyError:
             vals = super().interleave_data(prop)
 
-        # Check for consistency of occypancy values (i.e. add up to 1)
-        if np.any(self[element + '.occupancy.all'] != 1.0):
+        # Check for consistency of occupancy values (i.e. add up to 1)
+        if np.any(self[f"{element}.occupancy.all"] != 1.0):
             self._update_occupancy()
-            if np.any(self[element + '.occupancy.all'] != 1.0):
+            if np.any(self[f"{element}.occupancy.all"] != 1.0):
                 raise Exception(f"Occupancy doesn't add to unity in all {element}s")
+
         return vals
 
-    def regenerate_models(self, **kwargs):
+    def regenerate_models(self, propnames=None, exclude=[], deep=False):
         r"""
         Regenerate models associated with the Multiphase object
 
-        This method works by first regenerating the models associated with the
-        constituent phases, and then regenerating Multiphase models.
+        This method works by first regenerating the models associated with
+        the constituent phases, and then regenerating Multiphase models.
+
+        Parameters
+        ----------
+        propnames : List[str] or str
+            The list of property names to be regenerated. If None are
+            given then ALL models are re-run (except for those whose
+            ``regen_mode`` is 'constant').
+        exclude : List[str]
+            Since the default behavior is to run ALL models, this can be
+            used to exclude specific models. It may be more convenient to
+            supply as list of 2 models to exclude than to specify 8 models
+            to include.
+        deep : bool
+            Specifies whether or not to regenerate models on all
+            associated objects. For instance, if ``True``, then all
+            Physics models will be regenerated when method is called on
+            the corresponding Phase. The default is ``False``. The method
+            does not work in reverse, so regenerating models on a Physics
+            will not update a Phase.
 
         """
         # Regenerate models associated with phases within MultiPhase object
         for phase in self.phases.values():
-            phase.regenerate_models(**kwargs)
+            phase.regenerate_models(propnames=propnames, exclude=exclude, deep=deep)
         # Regenerate models specific to MultiPhase object
-        super().regenerate_models(**kwargs)
+        super().regenerate_models(propnames=propnames, exclude=exclude, deep=deep)
 
     def _set_automatic_throat_occupancy(self, mode="mean"):
         r"""
-        Automatically interpolates throat occupancy based on adjacent pores.
+        Automatically interpolates throat occupancy based on that in
+        adjacent pores. This method doesn't return any value.
+
+        Parameters
+        ----------
+        mode : str
+            Interpolation method, ex. 'mean' sets the throat occupancy as
+            the average of that in adjacent pores, while 'min' sets it to
+            the minimum of the two. Options are 'mean', 'min', and 'max'.
+
         """
         self.settings['throat_occupancy'] = 'automatic'
         for phase in self.phases.values():
-            self.add_model(propname=f"throat.occupancy.{phase.name}",
-                           model=misc.from_neighbor_pores,
-                           prop=f"pore.occupancy.{phase.name}",
-                           mode=mode)
+            self.add_model(
+                propname=f"throat.occupancy.{phase.name}",
+                model=misc.from_neighbor_pores,
+                prop=f"pore.occupancy.{phase.name}",
+                mode=mode
+            )
 
     def set_occupancy(self, phase, Pvals=[], Tvals=[], pores=[], throats=[]):
         r"""
-        Specify occupancy of a phase in each pore and/or throat
+        Specifies occupancy of a phase in each pore and/or throat. This
+        method doesn't return any value.
 
         Parameters
         ----------
-        phase : OpenPNM Phase object
-            The phase whose occupancy is being specified
-
-        Pvals : array_like, float
-            The volume fraction of ``phase`` in each pore. This array must be
-            *``Np``*-long, except when ``pores`` is also passed, where in that
-            case they must be of equal length, with one value between 0 and 1,
-            for each pore in the network. If a scalar is received it is applied
-            to all pores. If nothing is passed, ``Pvals=1.0`` is assumed.
-
-        Tvals : array_like, float
-            The volume fraction of ``phase`` in each throat. This array must be
-            *``Nt``*-long, except when ``throats`` is also passed, where in that
-            case they must be of equal length, with one value between 0 and 1,
-            for each throat in the network. If a scalar is received it is applied
-            to all throats. If nothing is passed, ``Tvals=1.0`` is assumed.
-
-        pores : array_like, int
-            The location of pores for which the phase occupancy is to be set.
-
-        throats : array_like, int
-            The location of throats for which the phase occupancy is to be set.
+        phase : GenericPhase
+            The phase whose occupancy is being specified.
+        Pvals : array_like
+            The volume fraction of ``phase`` in each pore. This array must
+            be Np-long, except when ``pores`` is also specified, where in
+            that case they must be of equal length. If a scalar is
+            received, it is applied to all the pores. If nothing is passed,
+            ``Pvals=1.0`` will be assumed.
+        Tvals : array_like
+            The volume fraction of ``phase`` in each throat. This array
+            must be Nt-long, except when ``throats`` is also specified,
+            where in that case they must be of equal length. If a scalar
+            is received it is applied to all the throats. If nothing is
+            passed, ``Tvals=1.0`` will be assumed.
+        pores : array_like
+            The location of pores whose occupancy is to be set.
+        throats : array_like
+            The location of throats whose occupancy is to be set.
 
         """
         # TODO: pores/throats could also be masks
@@ -284,4 +448,10 @@ class MultiPhase(GenericPhase):
 
         if self.settings["throat_occupancy"] == "automatic":
             self.regenerate_models(propnames=f"throat.occupancy.{phase.name}")
+
         self._update_occupancy()
+
+
+def _dummy(target, **kwargs):
+    r"""Dummy model to store propnames as kwargs for dependency handling"""
+    return target._assemble_partition_coef_global()
diff --git a/openpnm/utils/__init__.py b/openpnm/utils/__init__.py
index 4ddd8861a..64b31a419 100644
--- a/openpnm/utils/__init__.py
+++ b/openpnm/utils/__init__.py
@@ -24,6 +24,7 @@ from .misc import unique_list
 from .misc import tic, toc
 from .misc import is_symmetric
 from .misc import nbr_to_str
+from .misc import prettify_logger_message
 from .Workspace import Workspace
 from .Project import Project
 
diff --git a/openpnm/utils/misc.py b/openpnm/utils/misc.py
index 088ff257a..dc4ea0ba0 100644
--- a/openpnm/utils/misc.py
+++ b/openpnm/utils/misc.py
@@ -612,3 +612,12 @@ def nbr_to_str(nbr, t_precision):
             * (round(nbr, t_precision) != int(nbr)))
     nbr_str = (str(int(round(nbr, t_precision) * 10**n)) + (f'e-{n}') * (n != 0))
     return nbr_str
+
+
+def prettify_logger_message(msg):
+    r"""Prettifies logger messages by breaking them up into multi lines"""
+    from textwrap import wrap
+    linewidth = 75
+    indent = "\\n" + " " * 13
+    temp = wrap(msg, width=linewidth)
+    return indent.join(temp)
diff --git a/scripts/example_mass_partitioning.py b/scripts/example_mass_partitioning.py
index 9b94903eb..e2363a1fa 100644
--- a/scripts/example_mass_partitioning.py
+++ b/scripts/example_mass_partitioning.py
@@ -27,8 +27,7 @@ mphase._set_automatic_throat_occupancy()
 mphase.set_occupancy(phase=air, pores=[0, 1, 2, 3, 4])
 mphase.set_occupancy(phase=water, pores=[5, 6, 7, 8, 9])
 const = op.models.misc.constant
-mphase.set_binary_partition_coef(propname="throat.partition_coef",
-                                 phases=[water, air], model=const, value=0.5)
+mphase.set_binary_partition_coef(phases=[water, air], model=const, value=0.5)
 
 # Define physics object
 phys = op.physics.Standard(network=net, phase=mphase, geometry=geom)
diff --git a/scripts/example_multiphase_diffusion.py b/scripts/example_multiphase_diffusion.py
index 45eda7259..2678e075f 100755
--- a/scripts/example_multiphase_diffusion.py
+++ b/scripts/example_multiphase_diffusion.py
@@ -13,16 +13,22 @@ top and the bottom of the network.
 
 """
 import openpnm as op
-import numpy as np
+from openpnm.network import Cubic
+from openpnm.phases import Air, Water, MultiPhase
+from openpnm.geometry import StickAndBall
+from openpnm.physics import Standard
+from openpnm.models.misc import constant
+from openpnm.models.physics.generic_source_term import linear
 import matplotlib.pyplot as plt
+import numpy as np
 np.random.seed(10)
 export = False
 
 # Define network, geometry and constituent phases
-net = op.network.Cubic(shape=[100, 100, 1])
-geom = op.geometry.StickAndBall(network=net, pores=net.Ps, throats=net.Ts)
-air = op.phases.Air(network=net, name="air")
-water = op.phases.Water(network=net, name="water")
+net = Cubic(shape=[100, 100, 1])
+geom = StickAndBall(network=net, pores=net.Ps, throats=net.Ts)
+air = Air(network=net, name="air")
+water = Water(network=net, name="water")
 water["pore.diffusivity"] = air["pore.diffusivity"] * 0.05
 
 # Define the regions to be occupied by the two phases (air and water)
@@ -34,18 +40,16 @@ ts_air = net.find_neighbor_throats(pores=ps_air, mode="xnor")
 ts_interface = net.find_neighbor_throats(pores=ps_water, mode="xor")
 
 # Define multiphase and set phase occupancy
-mphase = op.phases.MultiPhase(network=net, phases=[air, water], name="mphase")
+mphase = MultiPhase(network=net, phases=[air, water], name="mphase")
 mphase._set_automatic_throat_occupancy()
 mphase.set_occupancy(air, pores=ps_air, throats=ts_air)
 mphase.set_occupancy(water, pores=ps_water, throats=ts_water)
 
 # Define physics
-phys = op.physics.Standard(network=net, phase=mphase, geometry=geom)
+phys = Standard(network=net, phase=mphase, geometry=geom)
 # Assign a partition coefficient (concentration ratio)
 K_water_air = 0.5   # c @ water / c @ air
-const = op.models.misc.constant
-mphase.set_binary_partition_coef(propname="throat.partition_coef",
-                                 phases=[water, air], model=const, value=K_water_air)
+mphase.set_binary_partition_coef(phases=[water, air], model=constant, value=K_water_air)
 # Replace the "default" ordinary_diffusion w/ multiphase_diffusion conductance model
 mdiff = op.models.physics.diffusive_conductance.multiphase_diffusion
 phys.add_model(propname="throat.diffusive_conductance", model=mdiff)
@@ -55,19 +59,23 @@ fd = op.algorithms.FickianDiffusion(network=net, phase=mphase)
 # Set source term
 phys["pore.A1"] = -1e-8 * geom["pore.area"]
 phys["pore.A2"] = 0.0
-linear = op.models.physics.generic_source_term.linear
-phys.add_model(propname="pore.rxn", model=linear, X="pore.concentration",
-                A1="pore.A1", A2="pore.A2", regen_mode="deferred")
-rxn_pores = net.pores(["left", "right"])
+phys.add_model(
+    propname="pore.rxn",
+    model=linear,
+    X="pore.concentration",
+    A1="pore.A1", A2="pore.A2",
+    regen_mode="deferred"
+)
+rxn_pores = net.pores(["front", "back"])
 net.set_label("rxn", pores=rxn_pores)
 fd.set_source(propname="pore.rxn", pores=rxn_pores)
 
 # Set BCs and run simulation
 net.set_label("air", pores=ps_air)
-front_air = net.pores(["front", "air"], mode="and")
-back_air = net.pores(["back", "air"], mode="and")
-fd.set_value_BC(pores=front_air, values=1.0)
-fd.set_value_BC(pores=back_air, values=0.1)
+left_air = net.pores(["left", "air"], mode="and")
+right_air = net.pores(["right", "air"], mode="and")
+fd.set_value_BC(pores=left_air, values=1.0)
+fd.set_value_BC(pores=right_air, values=0.1)
 fd.run()
 
 # Post-processing
diff --git a/tests/unit/algorithms/GenericTransportTest.py b/tests/unit/algorithms/GenericTransportTest.py
index a70ffb3c0..f8cf4858d 100644
--- a/tests/unit/algorithms/GenericTransportTest.py
+++ b/tests/unit/algorithms/GenericTransportTest.py
@@ -176,9 +176,9 @@ class GenericTransportTest:
         m.set_occupancy(phase=water, pores=[3, 4, 5])
         const = op.models.misc.constant
         K_water_air = 0.5
-        m.set_binary_partition_coef(propname="throat.partition_coef",
-                                    phases=[water, air],
-                                    model=const, value=K_water_air)
+        m.set_binary_partition_coef(
+            phases=[water, air], model=const, value=K_water_air
+        )
         m._set_automatic_throat_occupancy()
         _ = op.physics.Standard(network=net, phase=m, geometry=geom)
         alg = op.algorithms.GenericTransport(network=net, phase=m)
diff --git a/tests/unit/models/physics/DiffusiveConductanceTest.py b/tests/unit/models/physics/DiffusiveConductanceTest.py
index c7ec3912c..e3bc2939b 100644
--- a/tests/unit/models/physics/DiffusiveConductanceTest.py
+++ b/tests/unit/models/physics/DiffusiveConductanceTest.py
@@ -71,9 +71,9 @@ class DiffusiveConductanceTest:
         m.set_occupancy(phase=water, pores=[3, 4, 5])
         const = op.models.misc.constant
         K_water_air = 0.5
-        m.set_binary_partition_coef(propname="throat.partition_coef",
-                                    phases=[water, air],
-                                    model=const, value=K_water_air)
+        m.set_binary_partition_coef(
+            phases=[water, air], model=const, value=K_water_air
+        )
         m._set_automatic_throat_occupancy()
         phys = op.physics.GenericPhysics(network=net, phase=m, geometry=geom)
         mdiff = op.models.physics.diffusive_conductance.multiphase_diffusion
diff --git a/tests/unit/phases/MultiPhaseTest.py b/tests/unit/phases/MultiPhaseTest.py
index 83266d78c..872fdbca9 100644
--- a/tests/unit/phases/MultiPhaseTest.py
+++ b/tests/unit/phases/MultiPhaseTest.py
@@ -1,11 +1,12 @@
 import pytest
 import numpy as np
-import scipy as sp
-import openpnm as op
 from numpy.testing import assert_allclose
+import openpnm as op
+from openpnm.models.misc import constant
 
 
 class MultiPhaseTest:
+
     def setup_class(self):
         self.net = op.network.Cubic(shape=[10, 10, 10])
         self.water = op.phases.Water(network=self.net, name="water")
@@ -13,51 +14,63 @@ class MultiPhaseTest:
         self.oil = op.phases.Water(network=self.net, name="oil")
 
     def test_multiphase_init(self):
-        m = op.phases.MultiPhase(network=self.net, phases=[self.air,
-                                                           self.water])
+        phases = [self.air, self.water]
+        m = op.phases.MultiPhase(network=self.net, phases=phases)
+
         assert np.all(m['pore.occupancy.all'] == 0.0)
         assert np.all(m['throat.occupancy.all'] == 0.0)
+
         assert self.air.name in m.settings['phases']
         assert self.water.name in m.settings['phases']
 
     def test_multiphase_no_occupancy_yet(self):
-        m = op.phases.MultiPhase(network=self.net, phases=[self.air,
-                                                           self.water])
+        phases = [self.air, self.water]
+        m = op.phases.MultiPhase(network=self.net, phases=phases)
         self.water['pore.temperature'] = 300
+
         assert np.all(self.water['pore.temperature'] == 300)
+
         with pytest.raises(Exception):
             m['pore.temperature']
 
     def test_multiphase_set_occupancy_w_indices_only(self):
-        m = op.phases.MultiPhase(network=self.net, phases=[self.air,
-                                                           self.water])
+        phases = [self.air, self.water]
+        m = op.phases.MultiPhase(network=self.net, phases=phases)
+
         Ps_water = np.array([0, 1, 2])
         Ps_water_mask = self.net.tomask(pores=Ps_water)
         Ts_water = np.array([4, 12, 22])
         Ts_water_mask = self.net.tomask(throats=Ts_water)
+
         m.set_occupancy(self.water, pores=Ps_water, throats=Ts_water)
+
         assert m["pore.occupancy.water"][Ps_water_mask].mean() == 1
         assert m["pore.occupancy.water"][~Ps_water_mask].mean() == 0
         assert m["throat.occupancy.water"][Ts_water_mask].mean() == 1
         assert m["throat.occupancy.water"][~Ts_water_mask].mean() == 0
 
     def test_multiphase_set_occupancy_w_values_only(self):
-        m = op.phases.MultiPhase(network=self.net, phases=[self.air,
-                                                           self.water])
+        phases = [self.air, self.water]
+        m = op.phases.MultiPhase(network=self.net, phases=phases)
+
         Pvals = np.array([0.5, 0.9, 0.01])
         Tvals = np.array([0.9, 0.01])
+
         # Pvals must be Np-long if not accompanied by "pores" argument
         with pytest.raises(Exception):
             m.set_occupancy(self.water, Pvals=Pvals)
+
         # Tvals must be Nt-long if not accompanied by "throats" argument
         with pytest.raises(Exception):
             m.set_occupancy(self.water, Tvals=Tvals)
+
         # Set pore occupancy
         m.set_occupancy(self.water, Pvals=1)
         assert m["pore.occupancy.water"].mean() == 1
         Pvals = np.ones(self.net.Np) * 0.5
         m.set_occupancy(self.water, Pvals=Pvals)
         assert m["pore.occupancy.water"].mean() == 0.5
+
         # Setting throat occupancy
         m.set_occupancy(self.water, Tvals=1)
         assert m["throat.occupancy.water"].mean() == 1
@@ -66,30 +79,45 @@ class MultiPhaseTest:
         assert m["throat.occupancy.water"].mean() == 0.54
 
     def test_multiphase_set_occupancy_w_pore_indices_and_Pvals(self):
-        m = op.phases.MultiPhase(network=self.net, phases=[self.air,
-                                                           self.water])
+        phases = [self.air, self.water]
+        m = op.phases.MultiPhase(network=self.net, phases=phases)
+
         Ps_water = np.array([0, 1, 2])
         Pvals = np.array([0.5, 0.9, 0.01])
         Ps_water_mask = self.net.tomask(Ps_water)
         Ts_water = np.array([4, 12, 22])
         Ts_water_mask = self.net.tomask(throats=Ts_water)
         Tvals = np.array([0.3, 0.4, 0.1])
+
         # Pvals/Tvals and pores/throats; same array length
-        m.set_occupancy(self.water, pores=Ps_water, Pvals=Pvals,
-                        throats=Ts_water, Tvals=Tvals)
+        m.set_occupancy(
+            phase=self.water,
+            pores=Ps_water,
+            throats=Ts_water,
+            Pvals=Pvals,
+            Tvals=Tvals
+        )
         assert_allclose(m["pore.occupancy.water"][Ps_water], Pvals)
         assert_allclose(m["throat.occupancy.water"][Ts_water], Tvals)
         assert m["pore.occupancy.water"][~Ps_water_mask].mean() == 0
         assert m["throat.occupancy.water"][~Ts_water_mask].mean() == 0
+
         # Pvals and pores; inconsistent size
         with pytest.raises(Exception):
             m.set_occupancy(self.water, pores=[1, 5, 10], Pvals=[0.5, 0])
+
         # Tvals and throats; inconsistent size
         with pytest.raises(Exception):
             m.set_occupancy(self.water, throats=[10, 52, 0], Tvals=[0.5, 0.01])
+
         # Pvals/Tvals.size = 1 and pores/throats.size > 1
-        m.set_occupancy(self.water, pores=Ps_water, Pvals=0.25,
-                        throats=Ts_water, Tvals=0.23)
+        m.set_occupancy(
+            phase=self.water,
+            pores=Ps_water,
+            throats=Ts_water,
+            Pvals=0.25,
+            Tvals=0.23
+        )
         assert m["pore.occupancy.water"][Ps_water].mean() == 0.25
         assert m["throat.occupancy.water"][Ts_water].mean() == 0.23
 
@@ -122,68 +150,81 @@ class MultiPhaseTest:
         m = op.phases.MultiPhase(network=self.net)
         self.water['pore.temperature'] = 300
         self.air['pore.temperature'] = 200
+
         assert np.all(self.water['pore.temperature'] == 300)
         assert np.all(self.air['pore.temperature'] == 200)
+
         Ps = self.net['pore.coords'][:, 0] < 3
         Ts = self.net.tomask(throats=self.net.find_neighbor_throats(Ps))
+
         m.set_occupancy(phase=self.water, Pvals=Ps, Tvals=Ts)
         m.set_occupancy(phase=self.air, Pvals=~Ps, Tvals=~Ts)
+
         assert np.all(m['pore.temperature'] >= 200)
         assert np.all(m['pore.temperature'] <= 300)
 
     def test_mutliphase_partition_coef(self):
-        m = op.phases.MultiPhase(network=self.net,
-                                 phases=[self.water, self.air, self.oil])
+        phases = [self.water, self.air, self.oil]
+        m = op.phases.MultiPhase(network=self.net, phases=phases)
+
         x, y, z = self.net["pore.coords"].T
         ps_water = self.net.Ps[(y <= 3) + (y >= 8)]
         ps_air = self.net.Ps[(y > 3) * (y < 6)]
         ps_oil = self.net.Ps[(y >= 6) * (y < 8)]
+
         # Phase arrangement (y-axis): W | A | O | W
         m.set_occupancy(phase=self.water, pores=ps_water)
         m.set_occupancy(phase=self.air, pores=ps_air)
         m.set_occupancy(phase=self.oil, pores=ps_oil)
-        const = op.models.misc.constant
+
         K_air_water = 2.0
         K_air_oil = 1.8
         K_water_oil = 0.73
-        m.set_binary_partition_coef(propname="throat.partition_coef",
-                                    phases=[self.air, self.water],
-                                    model=const, value=K_air_water)
-        m.set_binary_partition_coef(propname="throat.partition_coef",
-                                    phases=[self.air, self.oil],
-                                    model=const, value=K_air_oil)
-        m.set_binary_partition_coef(propname="throat.partition_coef",
-                                    phases=[self.water, self.oil],
-                                    model=const, value=K_water_oil)
+
+        m.set_binary_partition_coef(
+            phases=[self.air, self.water], model=constant, value=K_air_water
+        )
+        m.set_binary_partition_coef(
+            phases=[self.air, self.oil], model=constant, value=K_air_oil
+        )
+        m.set_binary_partition_coef(
+            phases=[self.water, self.oil], model=constant, value=K_water_oil
+        )
+
         K_aw = m["throat.partition_coef.air:water"]
         K_ao = m["throat.partition_coef.air:oil"]
         K_wo = m["throat.partition_coef.water:oil"]
         K_global = m["throat.partition_coef.all"]
+
         assert np.isclose(K_aw.mean(), K_air_water)
         assert np.isclose(K_ao.mean(), K_air_oil)
         assert np.isclose(K_wo.mean(), K_water_oil)
+
         # Get water-air interface throats
         tmp1 = self.net.find_neighbor_throats(ps_water, mode="xor")
         tmp2 = self.net.find_neighbor_throats(ps_air, mode="xor")
-        Ts_water_air_interface = np.intersect1d(tmp1, tmp2)
+        Ts_water_air = np.intersect1d(tmp1, tmp2)
+
         # Get air-oil interface throats
         tmp1 = self.net.find_neighbor_throats(ps_air, mode="xor")
         tmp2 = self.net.find_neighbor_throats(ps_oil, mode="xor")
-        Ts_air_oil_interface = np.intersect1d(tmp1, tmp2)
+        Ts_air_oil = np.intersect1d(tmp1, tmp2)
+
         # Get oil-water interface throats
         tmp1 = self.net.find_neighbor_throats(ps_oil, mode="xor")
         tmp2 = self.net.find_neighbor_throats(ps_water, mode="xor")
-        Ts_oil_water_interface = np.intersect1d(tmp1, tmp2)
+        Ts_oil_water = np.intersect1d(tmp1, tmp2)
+
         # K_global for water-air interface must be 1/K_air_water
-        assert np.isclose(K_global[Ts_water_air_interface].mean(), 1/K_air_water)
+        assert np.isclose(K_global[Ts_water_air].mean(), 1/K_air_water)
+
         # K_global for air-oil interface must be K_air_oil (not 1/K_air_oil)
-        assert np.isclose(K_global[Ts_air_oil_interface].mean(), K_air_oil)
+        assert np.isclose(K_global[Ts_air_oil].mean(), K_air_oil)
         # K_global for oil-water interface must be 1/K_water_oil
-        assert np.isclose(K_global[Ts_oil_water_interface].mean(), 1/K_water_oil)
+        assert np.isclose(K_global[Ts_oil_water].mean(), 1/K_water_oil)
         # K_global for single-phase regions must be 1.0
-        interface_throats = np.hstack((Ts_water_air_interface,
-                                       Ts_air_oil_interface,
-                                       Ts_oil_water_interface))
+        temp = Ts_water_air, Ts_air_oil, Ts_oil_water
+        interface_throats = np.hstack(temp)
         Ts_single_phase = np.setdiff1d(self.net.Ts, interface_throats)
         assert np.isclose(K_global[Ts_single_phase].mean(), 1.0)
 
@@ -199,6 +240,44 @@ class MultiPhaseTest:
         # The next line ideally should throw an Exception, but warning for now
         m.set_occupancy(phase=self.water, Pvals=1.5, Tvals=2.5)
 
+    def test_format_interface_prop(self):
+        phases = [self.water, self.air]
+        m = op.phases.MultiPhase(network=self.net, phases=phases)
+        propname = m._format_interface_prop(propname="throat.foo", phases=phases)
+        assert propname == "throat.foo.water:air"
+
+    def test_get_phases_names(self):
+        phases = [self.water, self.air]
+        m = op.phases.MultiPhase(network=self.net, phases=phases)
+        names = m._get_phases_names("throat.foo.air:water")
+        assert names[0] == "air"
+        assert names[1] == "water"
+
+    def test_assemble_partition_coef_global(self):
+        # Create network and 3 pure phases
+        net = op.network.Cubic([6, 1, 1])
+        air = op.phases.Air(network=net, name="air")
+        water = op.phases.Water(network=net, name="water")
+        oil = op.phases.Water(network=net, name="oil")
+
+        # Create MultiPhase and set phase occupancies
+        m = op.phases.MultiPhase(network=net, phases=[water, air, oil])
+        m.set_occupancy(air, pores=[0, 1])
+        m.set_occupancy(water, pores=[2, 3])
+        m.set_occupancy(oil, pores=[4, 5])
+
+        # Check if K_global is initially ones
+        K_actual = m._assemble_partition_coef_global()
+        assert_allclose(K_actual, 1.)
+
+        # Add binary partition coefficient models
+        m.set_binary_partition_coef([water, air], model=constant, value=0.7)
+        m.set_binary_partition_coef([water, oil], model=constant, value=1.1)
+
+        # Check if K_global is correctly populated after models were added
+        K_actual = m._assemble_partition_coef_global()
+        assert_allclose(K_actual, [1., 1/0.7, 1., 1.1, 1])
+
 
 if __name__ == '__main__':
 
@@ -207,5 +286,5 @@ if __name__ == '__main__':
     t.setup_class()
     for item in t.__dir__():
         if item.startswith('test'):
-            print('running test: '+item)
+            print(f'Running test: {item}')
             t.__getattribute__(item)()