code: def _apply_BCs(self):
          r
          if 'pore.bc_rate' in self.keys():
              # Update b
              ind = np.isfinite(self['pore.bc_rate'])
              self.b[ind] = self['pore.bc_rate'][ind]
          if 'pore.bc_value' in self.keys():
              f = np.abs(self.A.data).mean()
              # Update b (impose bc values)
              ind = np.isfinite(self['pore.bc_value'])
              self.b[ind] = self['pore.bc_value'][ind] * f
              # Update b (substract quantities from b to keep A symmetric)
              x_BC = np.zeros(self.b.shape)
              x_BC[ind] = self['pore.bc_value'][ind]
              self.b[~ind] -= (self.A.tocsr() * x_BC)[~ind]
              # Update A
              P_bc = self.toindices(ind)
              indrow = np.isin(self.A.row, P_bc)
              indcol = np.isin(self.A.col, P_bc)
              self.A.data[indrow] = 0  # Remove entries from A for all BC rows
              self.A.data[indcol] = 0  # Remove entries from A for all BC cols
              datadiag = self.A.diagonal()  # Add diagonal entries back into A
              datadiag[P_bc] = np.ones_like(P_bc, dtype=np.float64) * f
              self.A.setdiag(datadiag)
              self.A.eliminate_zeros()
docs: r"""
          Applies all the boundary conditions that have been specified, by
          adding values to the *A* and *b* matrices.