def find_physics(self, geometry=None, phase=None):
        r"""
        Find the Physics object(s) associated with a given Geometry, Phase,
        or combination.

        Parameters
        ----------
        geometry : OpenPNM Geometry Object
            The Geometry object for which the Physics object(s) are sought

        phase : OpenPNM Phase Object
            The Phase object for which the Physics object(s) are sought

        Returns
        -------
        A list containing the Physics object(s).  If only a ``geometry`` is
        specified the the Physics for all Phases is returned.  If only a
        ``phase`` is specified, then the Physics for all Geometries is
        returned.  If both ``geometry`` and ``phase`` is specified then
        the list only contains a single Physics.  If no Physics is found, the
        the list will be empty.  See the Notes section for more information.

        See Also
        --------
        grid

        Notes
        -----
        The Project has an ``grid`` attribute that shows the association of
        all objects.  If each Geometry represents a row and each Phase is a
        column, then each row/col intersection represents a Physics. This
        method finds the PHysics' at each intersection

        """

        if geometry and phase:
            physics = self.find_physics(geometry=geometry)
            phases = list(self.phases().values())
            phys = physics[phases.index(phase)]
            return phys
        elif geometry:
            result = []
            net = self.network
            geoPs = net['pore.'+geometry.name]
            geoTs = net['throat.'+geometry.name]
            for phase in self.phases().values():
                physics = self.find_physics(phase=phase)
                temp = None
                for phys in physics:
                    Ps = phase.map_pores(pores=phys.Ps, origin=phys)
                    physPs = phase.tomask(pores=Ps)
                    Ts = phase.map_throats(throats=phys.Ts, origin=phys)
                    physTs = phase.tomask(throats=Ts)
                    if np.all(geoPs == physPs) and np.all(geoTs == physTs):
                        temp = phys
                result.append(temp)
            return result
        elif phase:
            names = set(self.physics().keys())
            keys = set([item.split('.')[-1] for item in phase.keys()])
            hits = names.intersection(keys)
            phys = [self.physics().get(i, None) for i in hits]
            return phys
        else:
            phys = list(self.physics().values())
            return phys