[SciPy-User] odeint - Excess work done on this call

[nicky van foreest]

Hi Farid,

Please see below.

On 11 October 2014 14:31, Farid Afandiyev <faridafandiyev at gmail.com> wrote:

> Hello!
>
> 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 and P_1, P_2 varying
> pressure of the nodes respectively.
>
> I end up writing differential questions below. Never mind physical meaning
> I just want to get math correct.
>
> dP_1/dt = dVa/dt * 1/(Vc_1 * B)
> dP_2/dt = dVa/dt * 1/(Vc_2 * B)
>
> Here B is compressibility.
>
> dVa/dt = A * K * sqrt(P_1 - P_2)
>
> dVa is amount of "flow" between nodes.
>
> K is some constant coefficient and A is link "weight". (P_1-P_2) can
> change sign so I've adjusted 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
>     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)
>     state = state_array[odeIterations-1]
>     print state_array
>     P_1 = state[0]
>     Va_1 = state[1]
>     P_2 = state[2]
>     Va_2 = state[3]
>
> 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.157000s
> [[  4.00000000e+03   2.64000000e-01   2.00000000e+03   3.30000000e-01]
>  [  3.49242034e+03   2.30499743e-01   2.20303186e+03   3.63500257e-01]
>  [  3.09580400e+03   2.04323064e-01   2.36167840e+03   3.89676936e-01]
>  [  2.81015098e+03   1.85469965e-01   2.47593961e+03   4.08530035e-01]
>  [  2.63546127e+03   1.73940444e-01   2.54581549e+03   4.20059556e-01]
>  [  2.57173487e+03   1.69734501e-01   2.57130605e+03   4.24265499e-01]
>  [  2.57142857e+03   1.69714286e-01   2.57142857e+03   4.24285714e-01]
>

It appears that your P1 - P2 = 0 at this line. It might mean that the
pressure in both vessels is the same at that moment in time, and then the
flow should stop, right?  So, what should the integrator tell you after
this point in time?

Nicky


>  [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
>  [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
>  [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]]
>  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. Also I don't need super accuracy from the odeint. Couple of digits
> after dot is more than enough for me. Any help is greatly appreciated!
>
> Kind regards,
>
> Farid.
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20141014/563eb899/attachment.html>