Repository: PMEAL/OpenPNM

Issue:
Fix MultiPhase to work with v3 changes

Discussion:

Patch:
diff --git a/openpnm/contrib/_multiphase.py b/openpnm/contrib/_multiphase.py
index f07acc2ac..1b0d6c203 100644
--- a/openpnm/contrib/_multiphase.py
+++ b/openpnm/contrib/_multiphase.py
@@ -2,7 +2,7 @@ import logging
 import numpy as np
 import openpnm.models.misc as misc
 from openpnm.phase import GenericPhase as GenericPhase
-from openpnm.utils import SettingsAttr, Docorator, TypedList
+from openpnm.utils import SettingsAttr, Docorator, TypedSet
 
 
 logger = logging.getLogger(__name__)
@@ -11,7 +11,7 @@ docstr = Docorator()
 
 @docstr.dedent
 class MultiPhaseSettings:
-    r"""
+    """
     Parameters
     ----------
     %(PhaseSettings.parameters)s
@@ -32,51 +32,45 @@ class MultiPhaseSettings:
     partition_coef_prefix : str
         The throat property which contains the partition coefficient values
     """
-    phases = TypedList(types=[str])
+    phases = TypedSet(types=[str])
     throat_occupancy = 'manual'
     partition_coef_prefix = 'throat.partition_coef'
 
 
 @docstr.dedent
 class MultiPhase(GenericPhase):
-    r"""
+    """
     Creates Phase object that represents a multiphase system consisting of
     a given list of GenericPhases.
 
     Parameters
     ----------
-    phases : list of OpenPNM Phase objects
-        A list containing the phase objects which comprise the multiphase
-        mixture
+    phases : list[GenericPhase]
+        A list containing the phase objects that make up the multiphase system
     %(GenericPhase.parameters)s
 
+    Notes
+    -----
+    This class assumes that only a SINGLE phase exists in each pore/throat.
+
     """
 
     def __init__(self, phases=[], settings=None, **kwargs):
         self.settings = SettingsAttr(MultiPhaseSettings, settings)
         super().__init__(settings=self.settings, **kwargs)
 
-        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)
 
-        # 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)
+        # Initialize the partition coefficient, K
+        self._K = np.ones(self.Nt, dtype=float)
+        prefix = self.settings["partition_coef_prefix"]
+        self[f"{prefix}.global"] = self._K
 
         # 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!')
-
     def __getitem__(self, key):
         try:
             vals = super().__getitem__(key)
@@ -84,45 +78,44 @@ class MultiPhase(GenericPhase):
             vals = self.interleave_data(key)
         return vals
 
-    def _get_phases(self):
+    @property
+    def phases(self):
         phases = {self.project[item].name: self.project[item]
                   for item in self.settings['phases']}
         return phases
 
-    phases = property(fget=_get_phases)
+    @property
+    def K(self):
+        self._build_K()
+        return self._K
 
-    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)
+    @K.setter
+    def K(self, value):
+        self._K = value
 
     def add_phases(self, phases):
-        r"""
-        Adds supplied phases to the MultiPhase object and initializes
-        occupancy to 0. It does not return a value.
+        """
+        Adds supplied phases to MultiPhase object and sets occupancy to 0.
 
         Parameters
         ----------
-        phases : List[GenericPhase]
+        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.
-                self[f'throat.occupancy.{phase.name}'] = 0.
+            self.settings['phases'].add(phase.name)
+            self[f'pore.occupancy.{phase.name}'] = 0.0
+            self[f'throat.occupancy.{phase.name}'] = 0.0
 
     def set_binary_partition_coef(self, phases, model, **kwargs):
-        r"""
+        """
         Sets binary partition coefficient as defined by the interface
         concentration ratio of phase 1 to phase 2.
 
         Parameters
         ----------
-        phases : List[GenericPhase]
+        phases : list[GenericPhase]
             List of the two phases for which the binary partition
             coefficient model is being added.
         model : OpenPNM model
@@ -131,36 +124,18 @@ class MultiPhase(GenericPhase):
             Keyword arguments to be passed to the ``model``.
 
         """
-        if np.size(phases) != 2:
-            raise Exception("'phases' must contain exactly two elements!")
-
+        assert len(phases) == 2
         # 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)
+        self._build_K()
 
     def _add_interface_prop(self, propname, phases, model, **kwargs):
-        r"""
-        Helper method used to add interface models to the ``MultiPhase``
-        object by augmenting.
+        """
+        Adds an interface model to the MultiPhase object. See Notes.
 
         Notes
         -----
-        For convenience and as an OpenPNM convention, the ``propname`` is
-        augmented as outlined by the following example.
-
-        Examples
-        --------
         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
@@ -181,54 +156,42 @@ class MultiPhase(GenericPhase):
         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.
-
-        """
+        """Formats propname as {propname}.{phase[0].name}:{phase[1].name}"""
         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 _get_phase_labels(self, formatted_propname):
+        """Retrieves phases names from a formatted propname"""
+        assert ":" in formatted_propname
+        return formatted_propname.split(".")[-1].split(":")
 
-    def _assemble_partition_coef_global(self):
-        r"""Updates the global partition coefficient array"""
+    def _get_interface_throats(self, phase1, phase2):
         conns = self.network.conns
-        partition_coef_models = []
+        occ1 = self[f"pore.occupancy.{phase1}"][conns]
+        occ2 = self[f"pore.occupancy.{phase2}"][conns]
+        idx12, = np.where((occ1[:, 0] == 1) & (occ2[:, 1] == 1))
+        idx21, = np.where((occ2[:, 0] == 1) & (occ1[:, 1] == 1))
+        return idx12, idx21
+
+    def _build_K(self):
+        """Updates the global partition coefficient array"""
         prefix = self.settings["partition_coef_prefix"]
-        K_global = np.ones(self.Nt)
-        partition_coef_global = prefix + ".all"
-
+        self._K = np.ones(self.Nt, dtype=float)
         # 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)
-
+        models = [k for k in self.models.keys() if k.startswith(prefix)]
         # 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
+        for model in models:
+            K12 = self[model]
+            phase1, phase2 = self._get_phase_labels(model)
+            idx12, idx21 = self._get_interface_throats(phase1, phase2)
+            self._K[idx12] = K12[idx12]
+            self._K[idx21] = 1 / K12[idx21]
+        # Store a reference in self as a propname for convenience
+        self[f"{prefix}.global"][:] = self._K
 
     def interleave_data(self, prop):
-        r"""
+        """
         Gathers property values from component phases to build a single
         array.
 
@@ -255,16 +218,9 @@ class MultiPhase(GenericPhase):
         try:
             for phase in self.phases.values():
                 vals += phase[prop] * self[f"{element}.occupancy.{phase.name}"]
-        # Otherwise - if not found - retrieve from super class
+        # Otherwise (if not found) retrieve from super class
         except KeyError:
             vals = super().interleave_data(prop)
-
-        # 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[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, propnames=None, exclude=[], deep=False):
@@ -276,11 +232,11 @@ class MultiPhase(GenericPhase):
 
         Parameters
         ----------
-        propnames : List[str] or str
+        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]
+        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
@@ -324,90 +280,42 @@ class MultiPhase(GenericPhase):
         """
         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
-            )
-
-    def set_occupancy(self, phase, Pvals=[], Tvals=[], pores=[], throats=[]):
+            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, *, pores=[], throats=[], values=1):
         r"""
-        Specifies occupancy of a phase in each pore and/or throat. This
+        Specifies occupancy of a phase in each pore or throat. This
         method doesn't return any value.
 
         Parameters
         ----------
         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
+        pores : ndarray
             The location of pores whose occupancy is to be set.
-        throats : array_like
+        throats : ndarray
             The location of throats whose occupancy is to be set.
+        values : ndarray or float
+            Pore/throat occupancy values.
 
         """
-        # TODO: pores/throats could also be masks
-
-        Pvals = np.array(Pvals, ndmin=1)
-        Tvals = np.array(Tvals, ndmin=1)
         pores = np.array(pores, ndmin=1)
         throats = np.array(throats, ndmin=1)
 
-        # Check for size consistency of the arguments
-        if Pvals.size and pores.size:
-            if Pvals.size != 1 and pores.size != 1:
-                if Pvals.size != pores.size:
-                    raise Exception("Pvals and pores must be the same size.")
-        if Tvals.size and throats.size:
-            if Tvals.size != 1 and throats.size != 1:
-                if Tvals.size != throats.size:
-                    raise Exception("Tvals and throats must be the same size.")
-
-        # Check if the passed phase is already part of MultiPhase object
+        if not(pores.size ^ throats.size):
+            raise Exception("Must either pass 'pores' or 'throats'")
         if phase not in self.project:
             raise Exception(f"{phase.name} doesn't belong to this project")
-        # Add the passed phase to MultiPhase object if not found
         if phase.name not in self.settings['phases']:
             self.add_phases(phase)
 
-        # Check for value consistency of the arguments
-        if np.any(Pvals > 1.0) or np.any(Pvals < 0.0):
-            logger.critical('Received Pvals contain volume fractions outside '
-                            + 'the range of 0 to 1')
-        if np.any(Tvals > 1.0) or np.any(Tvals < 0.0):
-            logger.critical('Received Tvals contain volume fractions outside '
-                            + 'the range of 0 to 1')
-
-        if Pvals.size and not pores.size:
-            pores = self.pores()
-        if Tvals.size and not throats.size:
-            throats = self.throats()
-
         if pores.size:
-            Pvals = Pvals if Pvals.size else 1.0
-            self['pore.occupancy.' + phase.name][pores] = Pvals
+            self[f'pore.occupancy.{phase.name}'][pores] = values
         if throats.size:
-            Tvals = Tvals if Tvals.size else 1.0
-            self['throat.occupancy.' + phase.name][throats] = Tvals
+            self[f'throat.occupancy.{phase.name}'][throats] = values
 
         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/models/physics/diffusive_conductance/_funcs.py b/openpnm/models/physics/diffusive_conductance/_funcs.py
index 324af77a3..2fecc08a4 100644
--- a/openpnm/models/physics/diffusive_conductance/_funcs.py
+++ b/openpnm/models/physics/diffusive_conductance/_funcs.py
@@ -206,7 +206,7 @@ def multiphase_diffusion(target,
                          pore_diffusivity="pore.diffusivity",
                          throat_diffusivity="throat.diffusivity",
                          size_factors="throat.diffusive_size_factors",
-                         partition_coef_global="throat.partition_coef.all"):
+                         partition_coef_global="throat.partition_coef.global"):
     r"""
     Calculates the diffusive conductance of conduits in network.
 
@@ -232,21 +232,19 @@ def multiphase_diffusion(target,
 
     """
     network = target.network
-    throats = target.throats(to_global=True)
-    phase = target
-    cn = network.conns[throats]
-    F = network[size_factors]
+    cn = network.conns
+    F1, Ft, F2 = network[size_factors].values()
 
     # Fetch model parameters
-    D1, D2 = phase[pore_diffusivity][cn].T
-    Dt = phase[throat_diffusivity][throats]
-    g1 = D1 * F[f"{size_factors}.pore1"][throats]
-    gt = Dt * F[f"{size_factors}.throat"][throats]
-    g2 = D2 * F[f"{size_factors}.pore2"][throats]
+    D1, D2 = target[pore_diffusivity][cn].T
+    Dt = target[throat_diffusivity]
+    g1 = D1 * F1
+    gt = Dt * Ft
+    g2 = D2 * F2
 
     # Apply Henry's partitioning coefficient
-    # Note: m12 = G21*c1 - G12*c2 NOT G12*c1 - G21*c2
-    K12 = phase[partition_coef_global][throats]
+    # Note: m12 = (G21*c1 - G12*c2)  NOT  (G12*c1 - G21*c2)
+    K12 = target[partition_coef_global]
     G21 = (1/g1 + 0.5/gt + K12 * (1/g2 + 0.5/gt)) ** -1
     G12 = K12 * G21
 
diff --git a/openpnm/phase/_generic.py b/openpnm/phase/_generic.py
index 7d497aeeb..74728c336 100644
--- a/openpnm/phase/_generic.py
+++ b/openpnm/phase/_generic.py
@@ -3,7 +3,6 @@ import numpy as np
 from openpnm.core import Domain
 from openpnm.utils import Workspace
 from openpnm.utils import Docorator, SettingsAttr
-import openpnm.models as mods
 
 
 docstr = Docorator()


Test patch:
diff --git a/tests/unit/models/physics/DiffusiveConductanceTest.py b/tests/unit/models/physics/DiffusiveConductanceTest.py
index c7ad3865d..9c16499c8 100644
--- a/tests/unit/models/physics/DiffusiveConductanceTest.py
+++ b/tests/unit/models/physics/DiffusiveConductanceTest.py
@@ -81,29 +81,29 @@ class DiffusiveConductanceTest:
         assert_allclose(actual, desired=0.117568, rtol=1e-5)
         del self.net["throat.diffusive_size_factors"]
 
-    # def test_multiphase_diffusion(self):
-    #     net = op.network.Cubic(shape=[1, 6, 1])
-    #     net['throat.diffusive_size_factors'] = self.size_factors_dict
-    #     air = op.phase.Air(network=net)
-    #     water = op.phase.Water(network=net)
-    #     from openpnm import contrib
-    #     m = contrib.MultiPhase(network=net, phases=[air, water])
-    #     m.set_occupancy(phase=air, pores=[0, 1, 2])
-    #     m.set_occupancy(phase=water, pores=[3, 4, 5])
-    #     const = op.models.misc.constant
-    #     K_water_air = 0.5
-    #     m.set_binary_partition_coef(phases=[water, air], model=const, value=K_water_air)
-    #     m._set_automatic_throat_occupancy()
-    #     mdiff = op.models.physics.diffusive_conductance.multiphase_diffusion
-    #     m.add_model(propname="throat.diffusive_conductance", model=mdiff)
-    #     g = m["throat.diffusive_conductance"]
-    #     # Diffusive conductance for MultiPhase must be Nt by 2 (not Nt by 1)
-    #     assert g.shape == (net.Nt, 2)
-    #     # Columns 1, 2 of conductance must be equal except for interface throat
-    #     assert_allclose(g[:, 0][[0, 1, 3, 4]], g[:, 1][[0, 1, 3, 4]])
-    #     # G12 and G21 at interface must differ (ratio must be K_water_air)
-    #     assert_allclose(g[2, 0] / g[2, 1], 1 / K_water_air)
-    #     assert_allclose(g.mean(), 7.54400e-7, rtol=1e-5)
+    def test_multiphase_diffusion(self):
+        net = op.network.Cubic(shape=[1, 6, 1])
+        net['throat.diffusive_size_factors'] = self.size_factors_dict
+        air = op.phase.Air(network=net)
+        water = op.phase.Water(network=net)
+        from openpnm import contrib
+        m = contrib.MultiPhase(network=net, phases=[air, water])
+        m.set_occupancy(phase=air, pores=[0, 1, 2])
+        m.set_occupancy(phase=water, pores=[3, 4, 5])
+        const = op.models.misc.constant
+        K_water_air = 0.5
+        m.set_binary_partition_coef(phases=[water, air], model=const, value=K_water_air)
+        m._set_automatic_throat_occupancy()
+        mdiff = op.models.physics.diffusive_conductance.multiphase_diffusion
+        m.add_model(propname="throat.diffusive_conductance", model=mdiff)
+        g = m["throat.diffusive_conductance"]
+        # Diffusive conductance for MultiPhase must be Nt by 2 (not Nt by 1)
+        assert g.shape == (net.Nt, 2)
+        # Columns 1, 2 of conductance must be equal except for interface throat
+        assert_allclose(g[:, 0][[0, 1, 3, 4]], g[:, 1][[0, 1, 3, 4]])
+        # G12 and G21 at interface must differ (ratio must be K_water_air)
+        assert_allclose(g[2, 0] / g[2, 1], 1 / K_water_air)
+        assert_allclose(g.mean(), 7.54400e-7, rtol=1e-5)
 
     def test_taylor_aris_diffusion(self):
         self.net['throat.diffusive_size_factors'] = self.size_factors_dict
diff --git a/tests/unit/phase/MultiPhase b/tests/unit/phase/MultiPhaseTest.py
similarity index 72%
rename from tests/unit/phase/MultiPhase
rename to tests/unit/phase/MultiPhaseTest.py
index e3e687fcc..e09794d75 100644
--- a/tests/unit/phase/MultiPhase
+++ b/tests/unit/phase/MultiPhaseTest.py
@@ -17,8 +17,8 @@ class MultiPhaseTest:
         phases = [self.air, self.water]
         m = op.contrib.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_allclose(m['pore.occupancy.air'], 0)
+        assert_allclose(m['throat.occupancy.water'], 0)
 
         assert self.air.name in m.settings['phases']
         assert self.water.name in m.settings['phases']
@@ -26,12 +26,8 @@ class MultiPhaseTest:
     def test_multiphase_no_occupancy_yet(self):
         phases = [self.air, self.water]
         m = op.contrib.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']
+        self.water['pore.temperature'] = 300.0
+        assert_allclose(m['pore.temperature'], 0)
 
     def test_multiphase_set_occupancy_w_indices_only(self):
         phases = [self.air, self.water]
@@ -42,12 +38,13 @@ class MultiPhaseTest:
         Ts_water = np.array([4, 12, 22])
         Ts_water_mask = self.net.to_mask(throats=Ts_water)
 
-        m.set_occupancy(self.water, pores=Ps_water, throats=Ts_water)
+        m.set_occupancy(self.water, pores=Ps_water)
+        m.set_occupancy(self.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
+        assert_allclose(m["pore.occupancy.water"][Ps_water_mask], 1)
+        assert_allclose(m["pore.occupancy.water"][~Ps_water_mask], 0)
+        assert_allclose(m["throat.occupancy.water"][Ts_water_mask], 1)
+        assert_allclose(m["throat.occupancy.water"][~Ts_water_mask], 0)
 
     def test_multiphase_set_occupancy_w_values_only(self):
         phases = [self.air, self.water]
@@ -56,26 +53,18 @@ class MultiPhaseTest:
         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
+        m.set_occupancy(self.water, pores=m.Ps, values=1)
+        assert_allclose(m["pore.occupancy.water"], 1)
         Pvals = np.ones(self.net.Np) * 0.5
-        m.set_occupancy(self.water, Pvals=Pvals)
-        assert m["pore.occupancy.water"].mean() == 0.5
+        m.set_occupancy(self.water, pores=m.Ps, values=Pvals)
+        assert_allclose(m["pore.occupancy.water"], 0.5)
 
         # Setting throat occupancy
-        m.set_occupancy(self.water, Tvals=1)
+        m.set_occupancy(self.water, throats=m.Ts, values=1)
         assert m["throat.occupancy.water"].mean() == 1
         Tvals = np.ones(self.net.Nt) * 0.54
-        m.set_occupancy(self.water, Tvals=Tvals)
+        m.set_occupancy(self.water, throats=m.Ts, values=Tvals)
         assert m["throat.occupancy.water"].mean() == 0.54
 
     def test_multiphase_set_occupancy_w_pore_indices_and_Pvals(self):
@@ -90,48 +79,37 @@ class MultiPhaseTest:
         Tvals = np.array([0.3, 0.4, 0.1])
 
         # Pvals/Tvals and pores/throats; same array length
-        m.set_occupancy(
-            phase=self.water,
-            pores=Ps_water,
-            throats=Ts_water,
-            Pvals=Pvals,
-            Tvals=Tvals
-        )
+        m.set_occupancy(self.water, pores=Ps_water, values=Pvals)
+        m.set_occupancy(self.water, throats=Ts_water, values=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
+        assert_allclose(m["pore.occupancy.water"][~Ps_water_mask], 0)
+        assert_allclose(m["throat.occupancy.water"][~Ts_water_mask], 0)
 
         # Pvals and pores; inconsistent size
         with pytest.raises(Exception):
-            m.set_occupancy(self.water, pores=[1, 5, 10], Pvals=[0.5, 0])
+            m.set_occupancy(self.water, pores=[1, 5, 10], values=[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])
+            m.set_occupancy(self.water, throats=[10, 52, 0], values=[0.5, 0.01])
 
         # Pvals/Tvals.size = 1 and pores/throats.size > 1
-        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
+        m.set_occupancy(phase=self.water, pores=Ps_water, values=0.25)
+        m.set_occupancy(phase=self.water, throats=Ts_water, values=0.23)
+        assert_allclose(m["pore.occupancy.water"][Ps_water], 0.25)
+        assert_allclose(m["throat.occupancy.water"][Ts_water], 0.23)
 
     def test_multiphase_automatic_throat_occupancy(self):
-        m = op.contrib.MultiPhase(network=self.net, phases=[self.air,
-                                                           self.water])
+        m = op.contrib.MultiPhase(network=self.net, phases=[self.air, self.water])
         pores = np.random.choice(self.net.Ps, size=100, replace=False)
         Pvals = np.random.random(pores.size)
-        m.set_occupancy(self.water, pores=pores, Pvals=Pvals)
+        m.set_occupancy(self.water, pores=pores, values=Pvals)
         P1, P2 = self.net["throat.conns"].T
         oc1, oc2 = [m["pore.occupancy.water"][x] for x in (P1, P2)]
         # Throats take average occupancy of adjacent pores
         m._set_automatic_throat_occupancy(mode="mean")
-        assert_allclose(m["throat.occupancy.water"], (oc1+oc2)/2)
+        assert_allclose(m["throat.occupancy.water"], np.mean([oc1, oc2], axis=0))
         # Throats take maximum occupancy of adjacent pores
         m._set_automatic_throat_occupancy(mode="max")
         assert_allclose(m["throat.occupancy.water"], np.maximum(oc1, oc2))
@@ -143,7 +121,8 @@ class MultiPhaseTest:
         m = op.contrib.MultiPhase(network=self.net)
         self.water['pore.temperature'] = 300
         assert np.all(self.water['pore.temperature'] == 300)
-        m.set_occupancy(phase=self.water, Pvals=1, Tvals=1)
+        m.set_occupancy(phase=self.water, pores=m.Ps, values=1)
+        m.set_occupancy(phase=self.water, throats=m.Ts, values=1)
         assert np.all(m['pore.temperature'] == 300)
 
     def test_multiphase_occupancy_set_two_phase(self):
@@ -157,13 +136,15 @@ class MultiPhaseTest:
         Ps = self.net['pore.coords'][:, 0] < 3
         Ts = self.net.to_mask(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)
+        m.set_occupancy(phase=self.water, pores=Ps, values=1)
+        m.set_occupancy(phase=self.water, throats=Ts, values=1)
+        m.set_occupancy(phase=self.air, pores=~Ps, values=1)
+        m.set_occupancy(phase=self.air, throats=~Ts, values=1)
 
         assert np.all(m['pore.temperature'] >= 200)
         assert np.all(m['pore.temperature'] <= 300)
 
-    def test_mutliphase_partition_coef(self):
+    def test_multiphase_partition_coef(self):
         phases = [self.water, self.air, self.oil]
         m = op.contrib.MultiPhase(network=self.net, phases=phases)
 
@@ -181,20 +162,14 @@ class MultiPhaseTest:
         K_air_oil = 1.8
         K_water_oil = 0.73
 
-        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
-        )
+        m.set_binary_partition_coef([self.air, self.water], model=constant, value=K_air_water)
+        m.set_binary_partition_coef([self.air, self.oil], model=constant, value=K_air_oil)
+        m.set_binary_partition_coef([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"]
+        K_global = m["throat.partition_coef.global"]
 
         assert np.isclose(K_aw.mean(), K_air_water)
         assert np.isclose(K_ao.mean(), K_air_oil)
@@ -237,8 +212,8 @@ class MultiPhaseTest:
 
     def test_multiphase_invalid_occupancy(self):
         m = op.contrib.MultiPhase(network=self.net, phases=[self.water, self.air])
-        # The next line ideally should throw an Exception, but warning for now
-        m.set_occupancy(phase=self.water, Pvals=1.5, Tvals=2.5)
+        m.set_occupancy(phase=self.water, pores=m.Ps, values=1.5)
+        m.set_occupancy(phase=self.water, throats=m.Ts, values=2.5)
 
     def test_format_interface_prop(self):
         phases = [self.water, self.air]
@@ -249,7 +224,7 @@ class MultiPhaseTest:
     def test_get_phases_names(self):
         phases = [self.water, self.air]
         m = op.contrib.MultiPhase(network=self.net, phases=phases)
-        names = m._get_phases_names("throat.foo.air:water")
+        names = m._get_phase_labels("throat.foo.air:water")
         assert names[0] == "air"
         assert names[1] == "water"
 
@@ -267,16 +242,14 @@ class MultiPhaseTest:
         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.)
+        assert_allclose(m.K, 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])
+        assert_allclose(m.K, [1., 1/0.7, 1., 1.1, 1])
 
 
 if __name__ == '__main__':