We discuss the parallel implementation of a two-level algebraic ILU(k)-based domain decomposition preconditioner using the PETSc library. We present strategies to improve performance and minimize communication among processes during setup and application phases. We compare our implementation with an off-the-shelf preconditioner in PETSc for solving linear systems arising in reservoir simulation problems, and show that for some cases our implementation performs better.The multilevel preconditioner has been an active research area for the last 30 years. One of the main representatives of this class is the algebraic multigrid (AMG) method [19,22,30,31], when used as a preconditioner rather than as a solver. It is widely used in oil reservoir simulation as part of the CPR (constrained pressure residual) preconditioner [6,33]. Despite its general acceptance there is room for new alternatives, as there are problems where AMG can perform poorly, see, for instance [13]. Among these alternatives we find the two-level preconditioners. There are many variations within this family of preconditioners, but basically we can discern at least two subfamilies: algebraic [1,4,12,14,18,20,23,28,32] and operator dependent [15,16,17]. Within the operator dependent preconditioners we should highlight the spectral methods [21,25].Notice that {Γ J } 1≤J≤P form a disjoint partition of Γ. See Figure 2.Finally, we define extended subdomains Ω J and extended local interfaces Γ J , which incorporate the portions of the interface connected to the subdomain interior Ω Int J . We define Ω J as(a jk = 0 or a k j = 0) .(2.6)