diff --git a/OpenPNM/Algorithms/__Drainage__.py b/OpenPNM/Algorithms/__Drainage__.py
index b537c1f331..de3276db43 100644
--- a/OpenPNM/Algorithms/__Drainage__.py
+++ b/OpenPNM/Algorithms/__Drainage__.py
@@ -574,7 +574,10 @@ def _calc_fractional_filling(self, element, pressure):
             temp_Pc = phys.models[key]['Pc']  # Store old Pc
             phys.models[key]['Pc'] = pressure
             # Regenerate Physics model and capture output locally
-            Snwp[phys.Pnet] = phys.models[key].run()
+            locs = phys._map(element=element,
+                             locations=phys._get_indices(element),
+                             target=self._net)
+            Snwp[locs] = phys.models[key].run()
             # Re-populate the residual element with the non-wetting phase
             if sp.any(self[element+'.residual']):
                 Snwp[self[element+'.residual']] = 1.0
diff --git a/OpenPNM/Algorithms/__GenericLinearTransport__.py b/OpenPNM/Algorithms/__GenericLinearTransport__.py
index db632b39fa..18053da137 100644
--- a/OpenPNM/Algorithms/__GenericLinearTransport__.py
+++ b/OpenPNM/Algorithms/__GenericLinearTransport__.py
@@ -641,10 +641,16 @@ def _do_inner_iteration_stage(self, guess, **kwargs):
                                          mode='update')
                     prop1 = 'pore.source_nonlinear_s1_' + source_name
                     mask1 = ~sp.isnan(self[prop1])
-                    s1[~sp.isnan(self[prop1])] = s1[mask1] + self[prop1][mask1]
+                    s1_temp_1 = s1[:self.Np][mask1]
+                    s1_temp_2 = self[prop1][mask1]
+                    mask_temp_1 = self.Ps[~sp.isnan(self[prop1])]
+                    s1[mask_temp_1] = s1_temp_1 + s1_temp_2
                     prop2 = 'pore.source_nonlinear_s2_' + source_name
                     mask2 = ~sp.isnan(self[prop2])
-                    s2[~sp.isnan(self[prop2])] = s2[mask2] + self[prop2][mask2]
+                    s2_temp_1 = s2[:self.Np][mask2]
+                    s2_temp_2 = self[prop2][mask2]
+                    mask_temp_2 = self.Ps[~sp.isnan(self[prop2])]
+                    s2[mask_temp_2] = s2_temp_1 + s2_temp_2
         self.s1 = s1
         self.s2 = s2
         # Modifying A and b
diff --git a/OpenPNM/Network/__Cubic__.py b/OpenPNM/Network/__Cubic__.py
index a94b4574df..f56fd72cf2 100644
--- a/OpenPNM/Network/__Cubic__.py
+++ b/OpenPNM/Network/__Cubic__.py
@@ -169,7 +169,7 @@ def _label_surfaces(self):
             if 'pore.'+label not in self.keys():
                 self['pore.'+label] = False
         if 'pore.boundary' in self.keys():
-            internal = -self['pore.boundary']
+            internal = ~self['pore.boundary']
         else:
             internal = self['pore.all']
         self['pore.internal'] = internal
diff --git a/OpenPNM/Physics/models/multiphase.py b/OpenPNM/Physics/models/multiphase.py
index 1ee1777358..16f9b322ab 100644
--- a/OpenPNM/Physics/models/multiphase.py
+++ b/OpenPNM/Physics/models/multiphase.py
@@ -108,7 +108,7 @@ def late_throat_filling(network, phase, physics,
     criteria such as the 'throat.inv_Pc' array on a *Drainage* algorithm.
 
     """
-    Swp = sp.ones(physics.Nt,)
+    Swp = sp.ones([physics.Nt, ])
     if Pc > 0:
         Swp = Swp_star*(physics[throat_entry_pressure]/Pc)**eta
     values = (1-Swp)*(physics[throat_entry_pressure] <= Pc)
@@ -157,7 +157,6 @@ def late_pore_filling(physics, phase, network,
     criteria such as the 'pore.inv_Pc' array on a *Drainage* algorithm.
 
     """
-
     # If pc_star has not yet been calculated, do so
     if pc_star not in physics.keys():
         pores = phase.Ps
@@ -166,7 +165,7 @@ def late_pore_filling(physics, phase, network,
         temp = sp.array([sp.amin(prop[row]) for row in neighborTs])
         physics[pc_star] = temp[physics.Pnet]
 
-    Swp = sp.ones(physics.Np,)
+    Swp = sp.ones([physics.Np, ])
     if Pc > 0:
         Swp = Swp_star*(physics[pc_star]/Pc)**eta
     else:
diff --git a/OpenPNM/Utilities/IO.py b/OpenPNM/Utilities/IO.py
index 83f2b48c63..924c935c42 100644
--- a/OpenPNM/Utilities/IO.py
+++ b/OpenPNM/Utilities/IO.py
@@ -1100,7 +1100,9 @@ def load(cls, path,
         network.extend(throat_conns=new_conns, labels='new_conns')
         for item in network.props('pore'):
             item = item.split('.')[1]
-            network['throat.'+item] = _sp.nan
+            arr = _sp.ones_like(network['pore.'+item])[0]
+            arr = _sp.tile(A=arr, reps=[network.Nt, 1])*_sp.nan
+            network['throat.'+item] = _sp.squeeze(arr)
             network['throat.'+item][network.throats('new_conns')] = \
                 network['pore.'+item][Ts]
         network.trim(pores=Ts)