2714.sagews

Droplet growth equation speed comparison in between pure Python and Cythonized code.

import math

def pcompute():
    # Environmental conditions
    S_inf = 1.05
    T_inf = 288.15
    T_inf_C = 288.15 - 273.15
    p = 101325

    R_v = 461.5
    e_v_inf = 611.2 * math.exp(17.67*T_inf_C/(T_inf_C + 243.5))
    rho_v_inf = (e_v_inf) / (R_v * T_inf)

    rho_lw = (999.8396 * T_inf_C**0. + 18.224944 * T_inf_C**1. + -7.92221*10**-3. * T_inf_C**2. +\
             -55.44846*10**-6. * T_inf_C**3. + 149.7562*10**-9. * T_inf_C**4 + -393.2952*10**-12. * T_inf_C**5.)\
             / (1 + 18.159725*10**-3. * T_inf_C)

    D_v = 2.11*10**-5 * (101325 / p) * (T_inf / 273.15)**1.94

    Fd = (rho_lw * R_v * T_inf) / (D_v * e_v_inf)
    Lv = (2500301.87956 - 2348.34983 * T_inf_C)
    ka = (2.382*10**-2 + 7.032*10**-5 * T_inf_C)
    Fk = (Lv/(R_v * T_inf) - 1)* (Lv*rho_lw)/(ka*T_inf)

    mu_a = 1.8325*10**-5 * ((296.16 + 120)/(T_inf + 120)) * (T_inf/296.16)**(3./2)
    R_d = 287.05
    rho_dry = (p - e_v_inf) / (R_d * T_inf)
    rho_air = rho_dry + rho_v_inf

    N_Sc = mu_a / (D_v * rho_air)

    dt = 0.01
    r_next = 5 * 10**-6
    tsec = 0

    print "Part b-2 - Elapsed times to reach 10, 25, 50, 100, 250 um radii:"

    while True:
        v_t = 842 * (2*r_next)**0.8
        N_Re = v_t * 2*r_next * rho_air / mu_a
        X = N_Sc**(1./3) * N_Re**(1./2)
        
        if X > 1.4:
            f_v = 0.78 + 0.308 * X
        else:
            f_v = 1.0 + 0.108 * X**2
        
        # Method 1
        #r_diff = (f_v * (S_inf - 1)) / ((Fd+Fk) * r_next)
        #r_next = r_next + r_diff * dt

        # Method 2
        r_next = math.sqrt(r_next**2 + dt*((2 * f_v * (S_inf - 1)) / (Fd+Fk)))

        if 10**-5 < r_next < 10.05*10**-6: print tsec
        if 25*10**-6 < r_next < 25.02*10**-6: print tsec
        if 50*10**-6 < r_next < 50.01*10**-6: print tsec
        if 100*10**-6 < r_next < 100.005*10**-6: print tsec
        if 250*10**-6 < r_next < 250.02*10**-6:
            print tsec
            break
        
        tsec+=0.01

if __name__=="__main__":
    pass

timeit("pcompute()", number=1)

1 loops, best of 3: 54 s per loop

%cython

cimport cython

cdef extern from "math.h":
    double sqrt(double)

cdef extern from "math.h":
    double exp(double)

@cython.cdivision(True)
def ccompute():
    # Environmental conditions
    cdef double S_inf = 1.05
    cdef double T_inf = 288.15
    cdef double T_inf_C = 288.15 - 273.15
    cdef int p = 101325

    cdef double R_v = 461.5
    cdef double e_v_inf = 611.2 * exp(17.67*T_inf_C/(T_inf_C + 243.5))
    cdef double rho_v_inf = (e_v_inf) / (R_v * T_inf)

    cdef double rho_lw = (999.8396 * T_inf_C**0. + 18.224944 * T_inf_C**1. + -7.92221*10**-3. * T_inf_C**2. +\
                         -55.44846*10**-6. * T_inf_C**3. + 149.7562*10**-9. * T_inf_C**4 + -393.2952*10**-12. * T_inf_C**5.)\
                         / (1 + 18.159725*10**-3. * T_inf_C)

    cdef double D_v = 2.11*10**-5. * (101325 / p) * (T_inf / 273.15)**1.94

    cdef double Fd = (rho_lw * R_v * T_inf) / (D_v * e_v_inf)
    cdef double Lv = (2500301.87956 - 2348.34983 * T_inf_C)
    cdef double ka = (2.382*10**-2. + 7.032*10**-5. * T_inf_C)
    cdef double Fk = (Lv/(R_v * T_inf) - 1) * (Lv*rho_lw)/(ka*T_inf)

    cdef double mu_a = 1.8325*10**-5. * ((296.16 + 120)/(T_inf + 120)) * (T_inf/296.16)**(3./2)
    cdef double R_d = 287.05
    cdef double rho_dry = (p - e_v_inf) / (R_d * T_inf)
    cdef double rho_air = rho_dry + rho_v_inf

    cdef double N_Sc = mu_a / (D_v * rho_air)

    cdef float dt = 0.01
    cdef double r_next = 5*10**-6.
    cdef float tsec = 0

    print "Part b-2 - Elapsed times to reach 10, 25, 50, 100, 250 um radii:"

    while True:
        v_t = 842 * (2*r_next)**0.8
        N_Re = v_t * 2*r_next * rho_air / mu_a
        X = N_Sc**(1./3) * N_Re**(1./2)
        
        if X > 1.4:
            f_v = 0.78 + 0.308 * X
        else:
            f_v = 1.0 + 0.108 * X**2.
        
        # Method 1
        #r_diff = (f_v * (S_inf - 1)) / ((Fd+Fk) * r_next)
        #r_next = r_next + r_diff * dt

        # Method 2
        r_next = sqrt(r_next**2. + dt*((2 * f_v * (S_inf - 1)) / (Fd+Fk)))

        if 10**-5. < r_next < 10.05*10**-6.: print tsec
        if 25*10**-6. < r_next < 25.02*10**-6.: print tsec
        if 50*10**-6. < r_next < 50.01*10**-6.: print tsec
        if 100*10**-6. < r_next < 100.005*10**-6.: print tsec
        if 250*10**-6. < r_next < 250.02*10**-6.:
            print tsec
            break
        
        tsec+=0.01

if __name__=="__main__":
    pass

timeit("ccompute()", number=1)

1 loops, best of 3: 325 ms per loop

54/0.325

166.153846153846