Solving a system of 4 coupled PDEs representing variable diffusivity

I have four partial differential equations representing mass conservation of two compressible fluid phases (marked by subscripts $p1$ and $p2$) in two different continuum media (marked by subscripts $c1$ and $c2$) and they are coupled together by the term $T_{c1c2}$ as shown below.

Continuum 1: \begin{align} \text{Phase 1: } \frac{\partial }{\partial t}\left(\frac{\phi_{c1}S_{p1c1}}{B_{p1c1}} \right) = \nabla\left(\lambda_{p1c1}\nabla{P_{p1c1}} \right) - \frac{T_{c1c2}}{V} - \frac{q_{p1c1}}{V}, \\ \text{Phase-2:} \frac{\partial }{\partial t}\left(\frac{\phi_{c1}S_{p2c1}}{B_{p2c1}} \right) = \nabla\left(\lambda_{p2c1}\nabla{P_{p2c1}} \right) - \frac{T_{c1c2}}{V} - \frac{q_{p2c1}}{V}, \end{align} Continuum 2: \begin{align} \text{Phase 1: }\frac{\partial }{\partial t}\left(\frac{\phi_{c2}S_{p1c2}}{B_{p1c2}} \right) = \nabla\left(\lambda_{p1c2}\nabla{P_{p1c2}} \right) - \frac{T_{c1c2}}{V} - \frac{q_{p1c2}}{V}, \\ \text{Phase 2: } \frac{\partial }{\partial t}\left(\frac{\phi_{c2}S_{p2c2}}{B_{p2c2}} \right) = \nabla\left(\lambda_{p2c2}\nabla{P_{p2c2}} \right) - \frac{T_{c1c2}}{V} - \frac{q_{p2c2}}{V}, \end{align}

Additionally, there are 4 more linear equations which are: $$S_{p1c1} + S_{p2c1} = 1$$ $$S_{p1c2} + S_{p2c2} = 1$$ $$P_{p2c1} - P_{p1c1} = P_{cap1}$$ $$P_{p2c2} - P_{p1c2} = P_{cap2}$$

Variables: In these 4 linear equations, $P_{cap1}$ and $P_{cap2}$ are known. Except for $\phi_{c1}$, $\phi_{c2}$ and $V$, which are constants, all other variables vary with time $t$ and distance $(x,y)$. Also, $$\frac{1}{B_{p1c1}} = 1+c_{p1}\left(P_{p1c1} - P^{STC}\right) \text{and} \frac{1}{B_{p2c1}} = 1+c_{p2}\left(P_{p2c1} - P^{STC}\right)$$ where $c_{p1}$, $c_{p2}$ and $P^{STC}$ are constants.

Initial Condition: Dirichlet constant pressure

Boundary Condition: Neumann no-flow boundary condition

Objective: To solve for saturations and pressures as a function of space and time, i.e.

$S_{p1c1}$, $S_{p2c1}$, $S_{p1c2}$, $S_{p2c2}$ and $P_{p1c1}$, $P_{p2c1}$, $P_{p1c2}$, $P_{p2c2}$,

This system of equations could be discretized into linear system of equations and then solved for each time step, however, that is a complex task and would take a lot of time. I would, however, like to use some form of "tool" which can solve this coupled system of equations in a mesh grid format. I saw an answer in other post suggested to use FiPy, but the feedback from the OP was that it's way too slow. Please suggest what could be a good tool to solve this problem.

• Please edit the equations next time in your posts. Copy and pasting from latex would not work in most of the cases. Jul 4 '15 at 2:55
• I did edit it. Can you cite where did I miss it?
– user5510
Jul 4 '15 at 3:54
• The two important equations were not rendering, I fixed it. But this is not a place for discussing that but for talking about your question. Jul 4 '15 at 3:57
• That's 4 equations and 8 unknowns. How do you intend to close the system? What form do the $B_{*}$ variables take? Are they known? Jul 4 '15 at 16:33
• @BillBarth: Thanks for pointing that out. I have edited the question to include the missing stuff.
– user5510
Jul 5 '15 at 1:07