# Index reduction of a DAE from a PDE system

I have a system of 2 non-linear, coupled PDEs that I would like to transform to a stiff ODE system to solve them using the method of lines. The 2 equations:

\begin{align} \frac{\partial \phi}{\partial t} &= - \phi^m \frac{P}{R - H(P)(R-1)}\\ \phi^m \frac{P}{R - H(P)(R-1)} &= \nabla \cdot [\phi^n(\nabla P + u_z)] \end{align}

$$\phi$$ and $$P$$ are my two unknowns, $$R$$ is a constant, $$H(P)$$ is the heaviside function of $$P$$ and $$u_z$$ is an upward directed unit vector.

The problem is that my second equation is an elliptic equation, and hence I have no time derivative in there. My idea is to transform them to Differential Algebraic Equations (DAE) using finite differences for the space derivatives and then to do index reduction to obtain an ODE system.

I have no experience with DAE, so I am stuck there. If we work in 1D in the z direction, I obtain with finite differences:

$$\frac{\partial \phi_i}{\partial t} = - \phi_i^m \frac{P_i}{R - H(P_i)(R-1)} \\ \phi_i^m \frac{P_i}{R - H(P_i)(R-1)} = \frac{\bigg[\phi_{i+\frac{1}{2}}^n (\frac{Pc_{i+1}-Pc_{i}}{\Delta z} + 1) - \phi_{i-\frac{1}{2}}^n (\frac{Pc_{i}-Pc_{i-1}}{\Delta z} + 1)\bigg]}{\Delta z}$$

From what I read, it seems similar to a Hessenberg index-1:

\begin{align} y' &= f(t, y, z) \\ 0 &= g(t, y, z) \end{align}

This would only require 1 additional derivative.

May I have a bit of help on how to start from there? Pointers to do that using symbolic computation like Sympy could also work!

Update 1

On @Wolfgang Bangerth's advice, I've rewritten my second equation by isolating P, my unknown for this equation:

$$\frac{\partial \phi_i}{\partial t} = - \phi_i^m \frac{P_i}{R - H(P_i)(R-1)} \\ (\frac{\phi_i^m \Delta z}{R - H(P_i)(R-1)} + \frac{\phi^n_{i+\frac{1}{2}}}{\Delta z} + \frac{\phi^n_{i-\frac{1}{2}}}{\Delta z}) P_i - \frac{\phi^n_{i-\frac{1}{2}}}{\Delta z} P_{i-1} - \frac{\phi^n_{i+\frac{1}{2}}}{\Delta z} P_{i+1}= (\phi_{i+\frac{1}{2}}^n - \phi_{i-\frac{1}{2}}^n)$$

• This system doesn't quite seem right. The first equation describes the evolution of $\phi$, and it contains $P$ is one of the quantities. One would think that the second equation therefore tells us how we can compute $P$ from $\phi$, but at the least you don't write it that way. Are you thinking of the second equation as a nonlinear equation for $P$? Jan 18 at 2:48
• Sorry if it is confusing. I wrote it that way because it makes sense physically, but less mathematically. Yeah the second equation is used to compute P. What I did before to solve the system was to use two nested Picard Iteration: A first one solving P implicitly and updating H(P) at each iteration, and a second one for calculating $\phi$ from the new values of P but it is quite inefficient because I wrote it myself. My goal would be to solve my system using DifferentialEquations.jl in julia. But I need an ODE system for that. Jan 18 at 9:23
• Would it be more readable for you if I isolate P on 1 side? I actually don't know what to do with H(P) as it depends on P.. Jan 18 at 9:26
• @lddingsite did you try discretizing the PDE into a DAE and using ModelingToolkit.jl structural simplify for the index reduction and tearing? It may need dummy derivatives to continue imposing the implicit constraint. Navier-Stokes is a classic example where custom solvers are not necessary due to the DAE transformations. Jan 20 at 0:44
• MTK needs work scaling (and there's a lot of work ongoing in that front), but 200x200x2 PDE should be fine. Otherwise you could try doing some of the reduction by hand and find the pattern, though I find that's a bit tedious and error prone. Jan 20 at 15:11

Your equation is not so different. The simplest approach to solving is to alternate: Use the previous $$\phi$$ to solve for $$P$$, then use $$P$$ in the equation for $$\phi$$ to advance by one time step. Then repeat. This is akin to a first-order (Lie) operator splitting, and in the porous media community is called the IMPES (implicit pressure, explicit saturation) method if you solve the equation for $$\phi$$ explicitly. You can make this a bit more accurate if you extrapolate from the previous values of $$\phi$$ to the next time step, and then solve for $$P$$ with that, instead of taking the old $$\phi$$. You can then take the previous and the new $$P$$ to solve the equation for $$\phi$$ with something like the Crank-Nicolson method.