do you have a tip how to workaround the global variable?

[alexlib]

Hi,

thanks for the great code, it's very fast.

I'm struggling though with a need for a global variable, or maybe an argument that is updated by the ODE itself

@nb.cfunc(lsoda_sig)
def particle_ode(t, y_, dydt, p_):
    y = nb.carray(y_,(3,))
    p = nb.carray(p_,(11,))
    rho1, rho, nu , Cd, g, rhop, d, zu, zl, Vc0, trec = p
   
    global tzl # <---- here 

                           
    zp   = y[0]                # particle position      [m]
    V    = y[1]                # particle velocity      [m/s]
    tzl =  y[2] 
    Cam = 0.5                  # added mass coefficient [-]

    # determine Vc
    if (y[0]<= zl):
        Vc = Vc0
        tzl = t     # < ----- here I need to adjust this as we move through with y[0]
    else:
        Vc = Vc0 * np.exp(-1*((t-tzl)/trec))
    FS = (rho(zp) - rho1)*Vc*g

    # force balance on the particle 
    dzpdt = V
    dVdt  = ( (rhop-rho(zp)) * Vp(d) * g \
            - 0.5 *Cd(V*d/nu(zp))* rho(zp) * Ap(d) * np.abs(V) * V \
            - FS \
            ) / (rhop*Vp(d) + Cam*rho(zp)*Vp(d))

    dydt = [dzpdt, dVdt, 0.0 ]
funcptr = particle_ode.address

[Nicholaswogan]

This if else behavior will probably make the integrator unstable. It puts a discontinuity in the ODEs. These cases are best handled by doing two integrations. One integration for each branch of the if else logic.

You specific case is hard to deal with with numbalsoda because you need event finding. you need to integrate to where y[0] <= zl, then stop, and start a new integration for y[0] > zl. numbalsoda doesn't have event finding.

[alexlib]

Thanks. Did exactly this :)
Thank for the great package