scipy odeint - Excess work done on this call

Question

I'm newbie both in calculus and Python/Scipy so I apologize if this question is too dumb. I'm trying to model flow between two pressure vessels. Let's say we have two points and a link between them like this

[$Vc_1$, $P_1$]----($A$)----[$Vc_2$, $P_2$]

$Vc_1$, $Vc_2$ are constants volumes of the nodes(vessels) and $P_1$, $P_2$ are varying pressures of the nodes respectively.

I end up writing differential questions below. Never mind physical meaning I just want to get math correct.

$\frac{\mathrm{dP_1} }{\mathrm{d} t} = \frac{\mathrm{dVa} }{\mathrm{d} t} \cdot \frac{1}{Vc_1 \cdot B}$

$\frac{\mathrm{dP_2} }{\mathrm{d} t} = \frac{\mathrm{dVa} }{\mathrm{d} t} \cdot \frac{1}{Vc_2 \cdot B}$

Here $B$ is compressibility.

$\frac{\mathrm{dVa} }{\mathrm{d} t} = A \cdot K \cdot \sqrt {P_1 - P_2}$

$\frac{\mathrm{dVa} }{\mathrm{d} t}$ is amount of "flow" or change of additional volume between nodes. $K$ is some constant coefficient and $A$ is link "throughput". ($P_1-P_2$) can change sign so I've adjusted for this in the software.

Below is Python program that I wrote to evaluate this.

#!/usr/bin/env python

import math
from scipy.integrate import odeint
from time import time
import numpy

B_compressibility = 0.0000033 # water compressibility
K = 0.747871759938 # coefficient

Vc_1 = 20
Vc_2 = 50
A = 0.01
P_1 = 4000
P_2 = 2000

def deriv(state, t):    
    _P_1 = state[0]
    _P_2 = state[2]   
    diff_P = _P_1 - _P_2
    flow_direction = math.copysign(1, diff_P)
    dVa = flow_direction * A * K * math.sqrt(abs(diff_P))
    dP_1 = -(dVa/Vc_1)/B_compressibility
    dP_2 =  (dVa/Vc_2)/B_compressibility
    #print 'IN  ', state
    #print 'OUT ', [dP_1, -dVa, dP_2, dVa]
    return [dP_1, -dVa, dP_2, dVa]

if __name__ == '__main__':
    Va_1 = Vc_1 * P_1 * B_compressibility
    Va_2 = Vc_2 * P_2 * B_compressibility
    odeIterations = 10
    timeperiod = numpy.linspace(0.0,1.0, odeIterations)    
    initial_state = [P_1, Va_1, P_2, Va_2]
    t0 = time()
    state_array = odeint(deriv, initial_state, timeperiod)
    t1 = time()
    print 'runtime %fs' %(t1-t0)        
    print state_array
    P_1, Va_1, P_2, Va_2 = state_array[odeIterations-1]

Below is output from program

Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.
runtime 0.041000s
[[  4.00000000e+003   2.64000000e-001   2.00000000e+003   3.30000000e-001]
 [  3.49242034e+003   2.30499743e-001   2.20303186e+003   3.63500257e-001]
 [  3.09580400e+003   2.04323064e-001   2.36167840e+003   3.89676936e-001]
 [  2.81015098e+003   1.85469965e-001   2.47593961e+003   4.08530035e-001]
 [  2.63546127e+003   1.73940444e-001   2.54581549e+003   4.20059556e-001]
 [  2.57173487e+003   1.69734501e-001   2.57130605e+003   4.24265499e-001]
 [  2.57142857e+003   1.69714286e-001   2.57142857e+003   4.24285714e-001]
 [  1.83357137e-299   1.80790662e-299   1.83145695e-299   1.87152166e-299]
 [  1.83276935e-299   1.80681296e-299   1.83182150e-299   1.87141230e-299]
 [  1.83379011e-299   1.80746916e-299   1.83320682e-299   1.83229543e-299]]
 lsoda--  at current t (=r1), mxstep (=i1) steps   
       taken on this call before reaching tout     
      In above message,  I1 =       500
      In above message,  R1 =  0.6110315150411E+00

odeint gives correct results up to 7th line and then something goes seriously wrong. I have searched in google and it looks like I'm not the only one who struggles with scipy. Everybody suggests increasing mxstep but that doesn't solve my problem. In addition it slows down the method significantly. Somebody suggested reducing accuracy but I don't know how to do that. Decreasing accuracy is OK if that help as I don't need super accuracy from the odeint. Couple of digits after dot is more than enough for me. Also I just need final values so in fact I want to decrease number of steps. Any help is greatly appreciated!

Answer 1

Typically, this situation is handled by:

• Using an "event function" or something to that effect to stop the integration when $P_{1} = P_{2}$. LSODA (the backend for scipy.integrate.odeint) does not have this capability, but other integrators, such as CVODE (part of SUNDIALS) do.
• Reformulating the right-hand side so it's continuously differentiable when $P_{1} = P_{2}$. It's a bit of a hack -- I remember an engineer explaining it to me when talking about process simulators -- but it gets the physics approximately right.

Typical existence theorems for solutions of differential equations assume that the right-hand side is at least Lipschitz continuous. This particular solver fails because you need at least continuous differentiability in your right-hand side in order for the integration formulas to make sense, and at $P_{1} = P_{2}$, your right-hand side is not continuously differentiable (implied because it is not Lipschitz).

Physically, flow should stop when the pressures equilibrate. When modeling the flow as 1-D, there can be pressure waves that flow back and forth between the vessels, but since you're modeling this situation as 0-D, there's no driving force for pressure waves, so when you reach $P_{1} = P_{2}$, stopping integration at that point should give you the correct result.

Answer 2

If you need to terminate the ODE solver at some point (P1==P2) you can use odespy. The solve method of the solvers accepts a terminate function that accepts the state u, the time t and the integration step step_no and returns a boolean. I used it today to stop the integration when any of the state elements is too low:

u,t = solver.solve(linspace(0,100,1000), terminate=lambda u,t,step_no: (u < 1e-6).any())

Answer 3

I modified your code to use the more up-to-date solve_ivp and it no longer gets the problem:

#!/usr/bin/env python

import math
from scipy.integrate import solve_ivp
from time import time
import numpy
import matplotlib.pyplot as plt


def deriv(t, state, A, B_compressibility, K, Vc_1, Vc_2):    
    _P_1 = state[0]
    _P_2 = state[2]   
    diff_P = _P_1 - _P_2
    flow_direction = math.copysign(1, diff_P)
    dVa = flow_direction * A * K * math.sqrt(abs(diff_P))
    dP_1 = -(dVa/Vc_1)/B_compressibility
    dP_2 =  (dVa/Vc_2)/B_compressibility
    #print 'IN  ', state
    #print 'OUT ', [dP_1, -dVa, dP_2, dVa]
    return [dP_1, -dVa, dP_2, dVa]

if __name__ == '__main__':
    Vc_1 = 20
    Vc_2 = 50
    A = 0.01
    P_1 = 4000
    P_2 = 2000
    B_compressibility = 0.0000033 # water compressibility
    K = 0.747871759938 # coefficient
     
    Va_1 = Vc_1 * P_1 * B_compressibility
    Va_2 = Vc_2 * P_2 * B_compressibility
    n_points = 100
    t_start = 0.0
    t_end = 2.0
    t_eval = numpy.linspace(t_start, t_end, n_points)    
    initial_state = [P_1, Va_1, P_2, Va_2]
    t0 = time()
    result = solve_ivp(deriv, [t_start, t_end], initial_state, t_eval=t_eval, args=( A, B_compressibility, K, Vc_1, Vc_2), method='Radau')
    t1 = time()
    print ('runtime %fs' %(t1-t0))        
    print (result.y)
    P_1, Va_1, P_2, Va_2 = result.y[0], result.y[1], result.y[2], result.y[3]
    
    for i in range(4):
        plt.figure()
        plt.plot(result.t, result.y[i])

Notes:

1. solve_ivp requires t, state as arguments instead of odeint's state, t
2. In this particular case, we get a smoother plot with Radau solver method. With the other solvers there are eratic oscillations that are probably explained by some numerical instability.