Help
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Ground Controller Lvl 1
Message 1 of 4

How to make a python code run faster? - Said Bolluk

I am a structural engineer and trying to plot a moment-curvature interaction diagram via python. There, in the code, are so much iteration required. I am aware that Matlab is one of the best platforms in terms of scientific calculations and iterations. However, I want to learn the essentials of python while using my profession. I will briefly explain the code.

 

There are couple forces acting on the section; compression and tension. Compression forces are due to the concrete layer and the upper steel layer while the tensions are caused by the lower steel layers. The c is the neutral axis where these forces equal the sum of external forces which is also the axial force acting on the member. For the first case, it assumed that the total axial load (N) acting on the member is equal to zero. This will result that the sum of all compression and the tension forces must be zero. This could be achieved by manipulating the neutral axis depth, c. As the sum of forces closes to the intended value N=0, the related c value will determine all the values. Is there any problem that can be fixed in this code? Is it possible to make the code run faster? I will attach the code below, thank you.

def moment_curvature():
    import numpy as np
    import matplotlib.pyplot as plt
    c_list = []
    M_list = []
    K_list = []
    strain = np.linspace(0, 0.0038, 10)
    for i in range(len(strain)):
        ec_i = strain[i]
        c_list.append(axis_depth(ec_i)[0])
        M_list.append(axis_depth(ec_i)[1])
        K_list.append(axis_depth(ec_i)[2])
    plt.plot(K_list, M_list)
    plt.xlabel('K (1/mm)')
    plt.ylabel('M (kN.m)')
    plt.show()
#######
def axis_depth(ec_i):
    import numpy as np
    ###
    h = 1000
    bw = 500
    fc = 40
    #fcd = fc / 1.5
    fs = 420
    fyd = fs / 1.15
    Es = 200000
    ecu = 0.0038
    e0 = 0.002
    ######
    neutral_axis = np.linspace(1, h, h*100)
    for k in range(len(neutral_axis)):
        c = neutral_axis[k]
        ### concrete layers
        sections = 100
        delta = c / sections
        area = delta * bw
        ### first layer of concrete: y_1, ec_1, fc_1, Fc_1
        y_1 = delta / 2
        ec_1 = ec_i * (1 - (y_1 / c))
        if ec_1 <= e0:
            fc_1 = fc * ((2 * ec_1 / e0)-(ec_1 / e0) ** 2)
        elif e0 < ec_1 <= ecu:
            fc_1 = fc * (1 - 0.15 * ((ec_1 - e0) / (ecu - e0 ))) 
        Fc_1 = fc_1 * area * 1e-3
        moment_arm = (h / 2) - y_1
        moment_1 = Fc_1 * moment_arm * 1e-3
        ###
        concrete_force = Fc_1
        concrete_moment = moment_1
        y = delta / 2
        for i in range(sections - 1):
            y += delta
            ec_layer = ec_i * (1 - (y / c))
            fc_layer = concrete_strain(ec_layer, e0, ecu, fc)
            Fc_layer = fc_layer * area * 1e-3
            concrete_force += Fc_layer
            moment_arm = (h / 2) - y
            moment_layer = Fc_layer * moment_arm * 1e-3 
            concrete_moment += moment_layer
        ### steel reinforcement properties: area, h_s, d_s, y_s
        steel_area = [1256.64, 628.32, 1256.64]
        #h_s = [50, 70, 950]
        d_s = [950, 930, 50]
        y_s = [-450, -430, 450]
        es_1 = ec_i * (1 - (d_s[0] / c))
        es_2 = ec_i * (1 - (d_s[1] / c))
        es_3 = ec_i * (1 - (d_s[2] / c))
        steel_stress = [es_1 * Es, es_2 * Es, es_3 * Es]
        steel_force = 0
        steel_moment = 0
        for j in range(len(steel_stress)):
            if steel_stress[j] >= 0:
                if abs(steel_stress[j]) >= fyd:
                    steel_stress[j] = fyd
            else:
                if abs(steel_stress[j]) >= fyd:
                    steel_stress[j] = -fyd
            steel_f = (steel_stress[j] * steel_area[j] * 1e-3) 
            steel_force += steel_f
            steel_moment += (steel_f * y_s[j] * 1e-3)
        ###
        N = 0
        total_force = concrete_force + steel_force
        total_moment = steel_moment + concrete_moment
        curvature = ec_i / c
        if (N - 0.1) < total_force < (N + 0.1):
            break
    return [c, total_moment, curvature]
        
######
def concrete_strain(ec_layer, e0, ecu, fc):
    if ec_layer <= e0:
        fc_layer = fc * ((2 * ec_layer / e0) - (ec_layer / e0) ** 2)
    elif e0 < ec_layer <= ecu:
        fc_layer = fc * (1 - 0.15 * ((ec_layer - e0) / (ecu - e0 ))) 
    return fc_layer
######
moment_curvature()




concrete_strees_strain.png

 

 

 

3 Replies
Highlighted
Commander Lvl 1
Message 2 of 4

Re: How to make a python code run faster? - Said Bolluk

Hey! As a first step, I suggest doing a profile of your code to better understand where things are taking a long time: https://docs.python.org/2/library/profile.html

Highlighted
Commander Lvl 2
Message 3 of 4

Re: How to make a python code run faster? - Said Bolluk

First step is always to try to improve the algorithm. Can the algorithm / steps be simplified to reduce calculations, eliminate work, etc.

 

Once the algorithm is optimized, next is to write in a faster language, C or assembly.

 

Please follow-up to let us know how you made out. For good karma, mark a reply as the answer if it helped!

Commander Lvl 1
Message 4 of 4

Re: How to make a python code run faster? - Said Bolluk

If you’re looking for a language that’s close to python but faster, check out Julia!