algebraic constraint to terminate ODE integration with scipy

Question

I'm using Scipy 14.0 to solve a system of ordinary differential equations describing the dynamics of a gas bubble rising vertically (in the z direction) in a standing still fluid because of buoyancy forces. In particular, I have an equation expressing the rising velocity U as a function of bubble radius R, i.e. U=dz/dt=f(R), and one expressing the radius variation as a function of R and U, i.e. dR/dT=f(R,U). All the rest appearing in the code below are material properties. I'd like to implement something to account for the physical constraint on z which, obviously, is limited by the liquid height H. I consequently implemented a sort of z<=H constraint in order to stop integration in advance if needed: I used set_solout in order to do so. The situation is that the code runs and gives good results, but set_solout is not working at all (it seems like z_constraint is never called actually...). Do you know why? Is there somebody with a more clever idea, may be also in order to interrupt exactly when z=H (i.e. a final value problem) ? is this the right way/tool or should I reformulate the problem?

thanks in advance

Emi

from scipy.integrate import ode 

Db0 = 0.001 # init bubble radius
y0, t0 = [ Db0/2 , 0. ], 0. #init conditions
H = 1

def y_(t,y,g,p0,rho_g,mi_g,sig_g,H):
    R = y[0] 
    z = y[1]
    z_ = ( R**2 * g * rho_g ) / ( 3*mi_g )  #velocity
    R_ = ( R/3 * g * rho_g * z_ ) / ( p0 + rho_g*g*(H-z) + 4/3*sig_g/R ) #R dynamics    
    return [R_, z_]

def z_constraint(t,y):
    H = 1  #should rather be a variable..
    z = y[1] 
    if z >= H:
        flag = -1 
    else:
        flag = 0
    return flag

r = ode( y_ )
r.set_integrator('dopri5')
r.set_initial_value(y0, t0)
r.set_f_params(g, 5*1e5, 2000, 40, 0.31, H)
r.set_solout(z_constraint)

t1 = 6
dt = 0.1

while r.successful() and r.t < t1:
    r.integrate(r.t+dt)

score 4 · Answer 1 · edited May 23 '17 at 11:58

You're running into this issue. For set_solout to work correctly, it must be called right after set_integrator, before set_initial_value. If you introduce this modification into your code (and set a value for g), integration will terminate when z >= H, as you want.

To find the exact time when the bubble reached the surface, you can make a change of variables after the integration is terminated by solout and integrate back with respect to z (rather than t) to z = H. A paper that describes the technique is M. Henon, Physica 5D, 412 (1982); you may also find this discussion helpful. Here's a very simple example in which the time t such that y(t) = 0.5 is found, given dy/dt = -y:

import numpy as np
from scipy.integrate import ode

def f(t, y):
    """Exponential decay: dy/dt = -y."""
    return -y

def solout(t, y):
    if y[0] < 0.5:
        return -1
    else:
        return 0

y_initial = 1
t_initial = 0

r = ode(f).set_integrator('dopri5')
r.set_solout(solout)
r.set_initial_value(y_initial, t_initial)

# Integrate until solout constraint violated
r.integrate(2)

# New system with t as independent variable: see Henon's paper for details.
def g(y, t):
    return -1.0/y

r2 = ode(g).set_integrator('dopri5')
r2.set_initial_value(r.t, r.y)

r2.integrate(0.5)

y_final = r2.t
t_final = r2.y

# Error: difference between found and analytical solution
print t_final - np.log(2)

algebraic constraint to terminate ODE integration with scipy

1 Answers1