As mentioned above, the function below works, however its very slow. I am very interested in using faster/optimised numpy (or other) vectorized alternatives. I have not posted the entire script here due to it being too large.

My specific question is - are there suitable numpy (or other) functions that I can use to 1) reduce run time and 2) reduce code volume of this function, specifically the for loop?

Edit: mass, temp, U and dpdh are functions that carry out simple algebraic calculations and return constants

def my_system(t, y, n, hIn, min, mAlumina, cpAlumina, sa, V):
    dydt = np.zeros(3 * n) #setting up zeros array for solution (solving for [H0,Ts0,m0,H1,Ts1,m1,H2,Ts2,m2,..Hn,Tsn,mn])
# y = [h_0, Ts_0, m_0, ... h_n, Ts_n, m_n]
# y[0] = hin
# y[1] = Ts0
# y[2] = minL

i=0

## Using thermo
T = temp(y[i],P) #initial T
m = mass(y[i],P) #initial m

#initial values
dydt[i] = (min * (hIn - y[i])   (U(hIn,P,min) * sa * (y[i   1] - T))) / m # dH/dt (eq. 2)
dydt[i   1] = -(U(hIn,P,min) * sa * (y[i   1] - T)) / (mAlumina * cpAlumina) # dTs/dt from eq.3
dmdt = dydt[i] * dpdh(y[i], P) * V # dm/dt (holdup variation) eq. 4b
dydt[i   2] = min - dmdt # mass flow out (eq.4a)

for i in range(3, 3 * n, 3): #starting at index 3, and incrementing by 3 because we are solving for 'triplets' [h,Ts,m] in each loop

    ## Using thermo
    T = temp(y[i],P)
    m = mass(y[i],P)

    # [h, TS, mdot]
    dydt[i] = (dydt[i-1] * (y[i - 3] - y[i])   (U(y[i-3], P, dydt[i-1]) * sa * (y[i   1] - T))) /  m # dH/dt (eq.2), dydt[i-1] is the mass of the previous tank
    dydt[i   1] = -(U(y[i-3], P, dydt[i-1]) * sa * (y[i   1] - T)) / (mAlumina * cpAlumina) # dTs/dt eq. (3)
    dmdt = dydt[i] * dpdh(y[i], P) * V # Equation 4b
    dydt[i   2] = dydt[i-1] - dmdt # Equation 4a

return dydt

The functions mass, temp, U, and dpdh used inside the my_system function all take numbers as input, perform some simple algebraic operation and return a number (no need to optimise these I am just providing them for further context)

def temp(H,P):
    # returns temperature given enthalpy (after processing function)
    T = flasher.flash(H=H, P=P, zs=zs, retry=True).T
    return T

def mass(H, P):
    # returns mass holdup in mol
    m = flasher.flash(H=H, P=P, zs=zs, retry=True).rho()*V
    return m

def dpdh(H, P):
    res = flasher.flash(H=H, P=P, zs=zs, retry=True)
    if res.phase_count == 1:
        if res.phase == 'L':
            drho_dTf = res.liquid0.drho_dT()
        else:
            drho_dTf = res.gas.drho_dT()
    else:
        drho_dTf = res.bulk._equilibrium_derivative(of='rho', wrt='T', const='P')
    dpdh = drho_dTf/res.dH_dT_P()
    return dpdh

def U(H,P,m):
    # Given T, P, m
    air = Mixture(['nitrogen', 'oxygen'], Vfgs=[0.79, 0.21], H=H, P=P)
    mu = air.mu*1000/mWAir #mol/m.s
    cp = air.Cpm #J/mol.K
    kg = air.k #W/m.K
    g0 = m/areaBed #mol/m2.s
    a = sa*n/vTotal #m^2/m^3 #QUESTIONABLE
    psi = 1
    beta = 10
    pr = (mu*cp)/kg
    re = (6*g0)/(a*mu*psi)
 hfs = ((2.19*(re**1/3))   (0.78*(re**0.619)))*(pr**1/3)*(kg)/diameterParticle

    h = 1/((1/hfs)   ((diameterParticle/beta)/kAlumina))

    return h

Reference Image: enter image description here

CodePudding user response：

For improving the speed, you can see Numba, which is useable if you use NumPy a lot but not every code can be used with Numba. Apart from that, the formulation of the equation system is confusing. You are solving 3 equations and adding the result to a single dydt list by 3 elements each. You can simply create three lists, solve each equation and add them to their respective list. For this, you need to re-write my_system as:

import numpy as np

def my_system(t, RHS, hIn, Ts0, minL, mAlumina, cpAlumina, sa, V):

    # get initial boundary condition values
    y1 = RHS[0]
    y2 = RHS[1]
    y3 = RHS[2]

    ## Using thermo
    T = # calculate T
    m = # calculate m

    # [h, TS, mdot] solve dy1dt for h, dy2dt for TS and dy3dt for mdot
    dy1dt = # dH/dt (eq.2), y1 corresponds to initial or previous value of dy1dt
    dy2dt = # dTs/dt eq. (3), y2 corresponds to initial or previous value of dy2dt
    dmdt =  # Equation 4b
    dy3dt = # Equation 4a, y3 corresponds to initial or previous value of dy3dt


    # Left-hand side of ODE
    LHS = np.zeros([3,])

    LHS[0] = dy1dt
    LHS[1] = dy2dt
    LHS[2] = dy3dt

    return LHS

In this function, you can pass RHS as a list with initial values ([dy1dt, dy2dt, dy3dt]) which will be unpacked as y1, y2, and y3 respectively and use them for respective differential equations. The solved equations (next values) will be saved to dy1dt, dy2dt, and dy3dt which will be returned as a list LHS.

Now you can solve this using scipy.integrate.odeint. Therefore, you can leave the for loop structure and solve the equations by using this method as follows:

hIn = #some val 
Ts0 = #some val 
minL = #some val
mAlumina = #some vaL
cpAlumina = #some val
sa = #some val
V = #some val
P = #some val

## Using thermo
T = temp(hIn,P) #initial T
m = mass(hIn,P) #initial m

#initial values
y01 = # calculate dH/dt (eq. 2)
y02 = # calculate dTs/dt from eq.3
dmdt = # calculate dm/dt (holdup variation) eq. 4b
y03 = # calculatemass flow out (eq.4a)

n = # time till where you want to solve the equation system
y0 = [y01, y02, y03]
step_size = 1
t = np.linspace(0, n, int(n/step_size)) # use that start time to which initial values corresponds
res = odeint(my_sytem, y0, t, args=(hIn, Ts0, minL, mAlumina, cpAlumina, sa, V,), tfirst=True)

print(res[:,0]) # print results for dH/dt
print(res[:,1]) # print results for dTs/dt
print(res[:,2]) # print results for Equation 4a

Here, I have passed all the initial values as y0 and chosen a step size of 1 which you can change as per your need.