Module conduction.toolbox.conduction_toolbox

Expand source code
# coding: utf-8

if __name__ == "__main__":


    # variables
    upperBC = 298.0 # Kelvin
    lowerBC = 0.03 # W/m2

    minX, maxX = 0.0, 50000.0
    minY, maxY = -35000.0, 0.0



    import numpy as np
    import csv
    import time
    from scipy.ndimage import imread
    from conductionNd_serial import ConductionND
    try: range=xrange
    except: pass



    matfile = 'geology.png'
    mat_img = imread(matfile, flatten=True)
    mat_img = np.flipud(mat_img)
    resY, resX = mat_img.shape

    mat_index = np.unique(mat_img)[::-1]


    lith_name = []
    index = []
    k_lith = []
    H_lith = []


    propfile = 'material_properties.csv'
    with open(propfile, 'r') as f:
        reader = csv.reader(f, delimiter=',')
        reader.next()
        for row in reader:
            lith_name.append(row[0])
            index.append(int(row[1]))
            k_lith.append(float(row[2]))
            H_lith.append(float(row[3]))


    if len(index) != len(mat_index):
        msg = "number of lithologies in {} and {} must be identical"
        raise ValueError(msg.format(matfile, propfile))



    # Map properties to mesh
    k_field = np.zeros_like(mat_img, dtype=np.float)
    H_field = np.zeros_like(mat_img, dtype=np.float)
    M_field = np.zeros_like(mat_img, dtype=np.int)

    fmt = "{:2} mapped {:15} with k = {:.2f} W/m/K,  H = {:.2f} uW/m3"

    for i in range(0, len(index)):
        mask = mat_img == mat_index[i]
        M_field[mask] = index[i]
        k_field[mask] = k_lith[i]
        H_field[mask] = H_lith[i]
        print(fmt.format(index[i], lith_name[i], k_lith[i], H_lith[i]*1e6))
        
    k_field = k_field.ravel()
    H_field = H_field.ravel()
    M_field = M_field.ravel()




    # Setup mesh
    mesh = ConductionND((minX,minY), (maxX,maxY), (resX, resY))





    # update conductivity and heat production fields
    mesh.update_properties(k_field, H_field)

    # set boundary conditions
    mesh.boundary_condition('maxY', upperBC, flux=False)
    mesh.boundary_condition('minY', lowerBC, flux=True)





    # solve
    walltime = time.time()
    T = mesh.solve()
    qy, qx = mesh.heatflux()
    print("\nTemperature solved in {:.1f} s\n".format(time.time()-walltime))


    # save to disk
    print("Saving table...")
    fmt = '%13.3f,%13.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f'
    header = 'x,y,T(C),k(W/m/K),H(uW/m3),HF_x(mW/m2),HF_y(mW/m2)'
    out = np.hstack([mesh.coords, T.reshape(-1,1)-273.14,\
                     k_field.reshape(-1,1), H_field.reshape(-1,1)*1e6,\
                     qx.reshape(-1,1)*1e3, qy.reshape(-1,1)*1e3])
    np.savetxt('out.csv', out, delimiter=',', fmt=fmt, header=header)


    import matplotlib.pyplot as plt
    def plot_figure(field, title, cmap):
        print("Saving {} figure...".format(title))
        fig = plt.figure()
        ax1 = fig.add_subplot(111)
        im1 = ax1.imshow(field.reshape(resY,resX), cmap=cmap,\
                         extent=(minX, maxX, minY, maxY), origin='lower')
        ax1.set_title(title)
        fig.colorbar(im1)
        fig.savefig(title+'.png', bbox_inches='tight')
        plt.clf()

    plot_figure(M_field, "material_map", 'Greys')
    plot_figure(T-273.14, "temperature", 'jet') # Celsius
    plot_figure(k_field, "conductivity", 'BuPu') # W/m/K
    plot_figure(H_field*1e6, "heat_production", 'Reds') # uW/m3
    plot_figure(qx*1e3, "heatflux_x", 'plasma') # mW/m2
    plot_figure(qy*1e3, "heatflux_y", 'plasma') # mW/m2