The Constrained Pressure Residual (CPR) preconditioner is a powerful and widely used technique in the numerical solution of systems of equations arising from the discretization of multiphase flow in porous media, particularly in reservoir simulation. It is designed to efficiently handle the strong coupling and ill-conditioning present in the linear systems generated by fully implicit or sequential formulations of flow and transport equations.

---

🌐 Background: Why CPR?

In reservoir simulation, the governing equations (typically mass and momentum conservation) are discretized using finite volume or finite element methods, leading to large, sparse, nonlinear systems. After linearization (e.g., via Newton-Raphson), one must solve a sequence of large linear systems of the form:

[
J \Delta x = -r
]

where:

• ( J ) is the Jacobian matrix,
• ( \Delta x ) is the update vector (pressure, saturation, etc.),
• ( r ) is the residual.

These systems are often ill-conditioned due to:

• High contrast in rock and fluid properties,
• Strong coupling between pressure and saturation,
• Different scales in variables (e.g., pressure vs. saturation).

Standard preconditioners (like ILU or AMG) may struggle with convergence.

---

✅ What is CPR?

CPR is a two-stage preconditioning strategy that exploits the physics-based block structure of the system. It typically assumes that pressure is the primary driver of flow, and other variables (like saturations or concentrations) are secondary.

The key idea:

Decouple pressure from the other variables, solve for pressure accurately first, then correct the full system.

---

🔧 CPR Preconditioner: Two-Stage Approach

Stage 1: Pressure Solve (Elliptic Part)

• Extract or approximate the pressure-pressure block (( A_{pp} )) from the Jacobian.

• Apply an effective solver (e.g., Algebraic Multigrid – AMG) to the constrained pressure residual system:

  [
  A_{pp} \Delta p = -r_p
  ]

  where ( r_p ) is the pressure component of the residual.

• Use this pressure update to predict updates in other variables (e.g., via explicit or implicit decoupling).

Stage 2: Full System Correction (Parabolic/Hyperbolic Part)

• Apply a second preconditioner (e.g., ILU, block ILU, or smoother) to the entire system, now better conditioned due to the prior pressure correction.
• This stage handles the remaining coupling (e.g., saturation, temperature, composition).

---

🧩 Typical CPR Structure

1. Construct a pressure matrix (often from the diagonal block ( A_{pp} ) or a reduced elliptic operator).
2. Precondition with AMG on the pressure block.
3. Apply a post-smoother (e.g., ILU(0), block-Jacobi) on the full system.

Mathematically, the CPR preconditioner ( M^{-1} ) can be written as:

[
M^{-1} = M_2^{-1} \cdot M_1^{-1}
]

where:

• ( M_1^{-1} ): pressure preconditioner (e.g., AMG on ( A_{pp} )),
• ( M_2^{-1} ): second-stage preconditioner (e.g., ILU on full ( J )).

---

📈 Variants of CPR

• CPR-1: Classic two-stage CPR; widely used in black-oil models.
• CPR-2: Extends CPR-1 by including saturation or compositional effects in the second stage more rigorously.
• CPR-3: Further extensions for thermal or compositional simulations.
• Multiscale CPR: Incorporates multiscale basis functions for pressure.
• CPR with Adaptive Thresholds: Dynamically selects which variables to include in the pressure system.

---

✅ Advantages

• Fast convergence for strongly coupled systems.
• Robustness in heterogeneous media.
• Leverages physical insight (pressure dominates flow).
• Works well with Newton-based solvers.

---

⚠️ Challenges

• Requires accurate identification of the pressure block.
• Performance depends on quality of AMG for the pressure system.
• May degrade with strong gravity or capillarity effects.
• Not always effective for highly nonlinear or compositional systems without modifications.

---

🛠️ Implementation Tips

• Use consistent upwinding or transmissibility-based operators to build ( A_{pp} ).
• Apply scaling to balance variables (e.g., pressure in psi, saturation in [0,1]).
• Use deflation or adaptive dropping in ILU for the second stage.
• Consider block preconditioning for compositional or thermal models.

---

📚 References

1. Wallis, J. R. (1983). Incomplete Gaussian elimination as a preconditioning for generalized conjugate gradient methods. SPE Reservoir Simulation Symposium.
2. Lacroix, S., Vassilevski, Y., Wheeler, J., & Wheeler, M. (2003). Iterative solution methods for modeling multiphase flow in porous media fully implicitly. SIAM Journal on Scientific Computing.
3. Cusini, M., et al. (2015). A general framework for the simulation of multiphase flow in porous media. Computational Geosciences.
4. Jansen, J. D. (2013). Adjoint-based optimization of multi-phase flow through porous media. Lecture notes on reservoir simulation.

---

🔚 Summary

The CPR preconditioner is a physics-informed, two-stage method that:

1. Solves for pressure using a robust elliptic solver (like AMG),
2. Corrects the full system with a secondary preconditioner.

It is a cornerstone of modern reservoir simulators (e.g., in Eclipse, INTERSECT, AD-GPRS, OPM) due to its efficiency and scalability.

Let me know if you’d like a pseudocode implementation or a comparison with other preconditioners (e.g., Schur complement, block ILU).