I am studying the time-integration of the following paper,

Young, L. C. (1981). A finite-element method for reservoir simulation. Society of Petroleum Engineers Journal, 21(01), 115-128.

A copy (PDF) can be found on the author's ResearchGate: https://www.researchgate.net/publication/236367847_A_finite_element_method_for_reservoir_simulation

After spatial discretization with a particular finite-element method, the author ends up with the continuous time approximation for the pressure $p$, $$ (\rho'D) \frac{dp}{dt}+Bp = F $$ where $D$ is a diagonal mass matrix, and $B$ is a flow-coefficient matrix (analogous to the stiffness matrix) that depends on the pressure $p$. The factor $\rho' = d\rho/dp$ originates from an equation of state $\rho = \rho(p)$.

The author uses the following two-stage implicit time discretization:

$$ \frac{2}{\Delta t}(\rho' D)^{n}(p^{n+\frac{1}{2}} - p^n) + B^n p^{n+\frac{1}{2}} = F \\ (\rho' D)^{n + \frac{1}{2}} \frac{p^{n+1} - p^n}{\Delta t} + B^{n+\frac{1}{2}} p^{n+\frac{1}{2}} = F $$

In the first stage, the values from the previous time-step are used to calculate the effective mass-matrix and flow coefficient matrix $B$ thereby avoiding expensive non-linear iterations. In the second stage the matrices are re-evaluated at the middle of the time interval, but the stage is explicit.

Q1: Does this integration scheme have a name? If not, what would be the Butcher tableau?

Q2: Can the scheme be combined with a backward Euler method for adaptive time-stepping?

  • 1
    $\begingroup$ If it were a Runge-Kutta method, it cannot have terms such as $B^n p^{n+\frac12}$. That's because that is a term that comes from the right hand side of $\dot p = f(p)$, so everything in the right hand side must be evaluated at the same time -- but the term indicated does not follow this rule. $\endgroup$ Sep 25, 2023 at 23:23
  • $\begingroup$ You could look at $B_n$ as the initial guess of a fixed-point iteration, but I understand your point. The purpose of this simplification is to avoid having to solve a non-linear problem. I've seen this done in other sources too and in limit of small timesteps I guess it is consistent. The question is up to which step-size is this type of "linearization in time" justifiable. $\endgroup$
    – IPribec
    Sep 26, 2023 at 9:35
  • $\begingroup$ The first equation is first-order accurate, the second second-order. In total one should get at least the stability of a generic second-order explicit method. Thus approximately stable for $L\Delta t\le 1$, where $L$ is a Lipschitz constant for $(\rho'D)^{-1}B$. $\endgroup$ Sep 26, 2023 at 11:48
  • $\begingroup$ I didn't mean to imply that the method can't be good. I was just pointing out that it is not an RK scheme. $\endgroup$ Sep 26, 2023 at 17:11

1 Answer 1


As Wolfgang stated in the comments, this is not a traditional RK due to the inconsistent time evaluations within a stage. At first it would seem it can't even be cast as an additive RK since terms with inconsistent time are multiplied and not simply added. One trick (see section 2.4) is to artificially duplicate the ODE $$ \begin{align} \frac{dp}{dt} &= (\rho'(q) D(q))^{-1} (F - B(q) p) \\ \frac{dq}{dt} &= (\rho'(q) D(q))^{-1} (F - B(q) p) \end{align} $$ then apply an additive RK method. The first equation will be treated with an implicit method, and the second equation will be treated with an explicit method. The careful choice of where $p$ vs $q$ appear facilitate the nonlinear splitting. With some coefficient matching, it turns out the additive RK tableau is

$$ \begin{array}{cc|cc} \frac{1}{2} & 0 & 0 & 0 \\ \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ \hline 0 & 1 & 0 & 1 \end{array} $$

and the stage computations are

$$ \begin{align} P_1 &= p^n + \frac{\Delta t}{2}(\rho'(Q_1) D(Q_1))^{-1} (F - B(Q_1) P_1) \\ Q_1 &= q^n = p^n \\ P_2 &= P_1 \\ Q_2 &= P_2 \\ p^{n+1} &= q^{n+1} = p^n + \Delta t (\rho'(Q_2) D(Q_2))^{-1} (F - B(Q_2) P_2) \end{align} $$

So this additive RK method is combining the implicit and explicit midpoint methods. While these methods have been paired extensively (see, for instance section 2.3), there are several ways to do this. I've not seen this exact tableau before, but I wouldn't be surprised if it can be found somewhere in the literature.

To answer your second question, you could use a backward Euler step to get a lower order solution for error estimation. From the tableau, we can see this would not be an embedding so it would take an extra stage. A better idea would be to apply a different second order (or higher!) IMEX method with an embedding to the duplicated ODE above.

  • $\begingroup$ I was not aware of additive RK. Thanks for this excellent answer! $\endgroup$
    – IPribec
    Sep 28, 2023 at 7:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.