Python Program to Implement Steady State Shock Equations Assignment Solution.

import numpy as np

import matplotlib.pyplot as plt

M = 6.0

gam = 5./3.

# ==========================================================================

# ==========================================================================

def ddx(x, F):

    Ns = F.size - 1

    dFdx = np.zeros(Ns + 1, dtype='double')

    dFdx[0] = (F[1] - F[0]) / (x[1] - x[0])

    for i in range(1, Ns):

        dFdx[i] = (F[i + 1] - F[i - 1]) / (x[i + 1] - x[i - 1])

    dFdx[Ns] = (F[Ns] - F[Ns - 1]) / (x[Ns] - x[Ns - 1])

    return dFdx

# ==========================================================================

# ==========================================================================

if __name__ == '__main__':

    electron = np.loadtxt('./electron_quantites.txt')

    ion = np.loadtxt('./ion_quantites.txt')

    x = electron[:,0] # positions

    rho = ion[:,1] # this is rho/rho0.

    ke = electron[:,2]

    Te = electron[:,1]

    Ti = ion[:,2]

    ki = ion[:,3]

    eta = ion[:,4]

    u = ion[:,5]

    p = ion[:,6]

    E = ion[:,7]

    # First, calculate du/dx

    dudx = ddx(x, u)

    # Calculate dTi/dx

    dTidx = ddx(x, Ti)

    # Calculate dTe/dx

    dTedx = ddx(x, Te)

    # From equation (1b), calculate p0u0^2 + p0

    eq1b = p + gam*M**2.*rho*u**2. - eta*dudx

    eq1b_right = gam*M**2 + p[0]

    # Equation (1c)

    eq1c = E*u + p*u - ki*dTidx - ke*dTedx - eta*u*dudx

    eq1c_right = E[0]/rho[0] + p[0]

    # Eq. 1a

    eq1a = rho*u

    #==========================================================================

    #==========================================================================

    # create plot

    fig, axes = plt.subplots(nrows = 2, ncols = 3, figsize=(7,7))

    # Plot dTidx

    axes[0,0].plot(x, rho)

    axes[0,0].set_title(r'$\rho$')

    axes[0,0].set_xlabel('Position')

    axes[0,0].set_ylabel('Plasma mass density')

    axes[0,0].grid(True)

    axes[0,1].plot(x, ke, label = 'ke')

    axes[0,1].set_title(r'$k_{e}$')

    axes[0,1].set_xlabel('Position')

    axes[0,1].set_ylabel('Electron Conductivity')

    axes[0, 1].grid(True)

    axes[0,2].plot(x, Te, label = 'Te')

    axes[0,2].set_title(r'$T_{e}$')

    axes[0,2].set_xlabel('Position')

    axes[0,2].set_ylabel('Electron Temperature')

    axes[0, 2].grid(True)

    axes[1,0].plot(x, Ti, label = 'Ti')

    axes[1,0].set_title(r'$T_{i}$')

    axes[1,0].set_xlabel('Position')

    axes[1,0].set_ylabel('Ion Conductivity')

    axes[1, 0].grid(True)

    axes[1,1].plot(x, ki, label = 'ki')

    axes[1,1].set_title(r'$k_{i}$')

    axes[1, 1].set_xlabel('Position')

    axes[1, 1].set_ylabel('Ion conductivity')

    axes[1, 1].grid(True)

    axes[1,2].plot(x, eta, label = 'eta')

    axes[1,2].set_title(r'$\eta$')

    axes[1, 2].set_xlabel('Position')

    axes[1, 2].set_ylabel('Ion viscosity')

    axes[1, 2].grid(True)

    plt.show()

    fig, axes = plt.subplots(nrows = 2, ncols = 3, figsize=(7,7))

    axes[0,0].plot(x, u, label = 'u')

    axes[0,0].set_title(r'$u$')

    axes[0, 0].set_xlabel('Position')

    axes[0, 0].set_ylabel('Bulk fluid velocity')

    axes[0, 0].grid(True)

    axes[0,1].plot(x, p, label = 'p')

    axes[0,1].set_title(r'$p$')

    axes[0, 1].set_xlabel('Position')

    axes[0, 1].set_ylabel('Plasma pressure')

    axes[0, 1].grid(True)

    axes[0,2].plot(x, E, label = 'E')

    axes[0,2].set_title(r'$E$')

    axes[0, 2].set_xlabel('Position')

    axes[0, 2].set_ylabel('Total Plasma Energy')

    axes[0, 2].grid(True)

    axes[1,0].plot(x, dudx, label = 'dudx')

    axes[1,0].set_title(r'$\frac{\partial u}{\partial x}$')

    axes[1, 0].set_xlabel('Position')

    axes[1, 0].set_ylabel('Rate of change of bulk velocity')

    axes[1, 0].grid(True)

    axes[1,1].plot(x, dTidx, label = 'dTidx')

    axes[1,1].set_title(r'$\frac{\partial T_{i}}{\partial x}$')

    axes[1, 1].set_xlabel('Position')

    axes[1, 1].set_ylabel('Rate of change of Ion Temp.')

    axes[1, 1].grid(True)

    axes[1,2].plot(x, dTedx, label = 'dTedx')

    axes[1,2].set_title(r'$\frac{\partial T_{e}}{\partial x}$')

    axes[1, 2].set_xlabel('Position')

    axes[1, 2].set_ylabel('Rate of change of Electron Temp.')

    axes[1, 2].grid(True)

    plt.show()

    fig = plt.figure(figsize=(7, 5), dpi=120)

    plot = plt.plot(x, (eq1b-eq1b_right)/eq1b_right*100., 'royalblue', linestyle='--', label='Eq. (1b) percent error')

    plot = plt.plot(x, (eq1c-eq1c_right)/eq1c_right*100., 'tomato', linestyle='-', label='Eq. (1c) percent error')

    plot = plt.plot(x, (eq1a - 1.0)/1.0 *100, 'purple', linestyle='-', label = 'Eq. (1a) percent error')

    plt.ylabel(r'Percent error [%]')

    plt.xlabel(r'Position $\hat{x}$')

    #plt.xlim([25.0, 32.5])

    plt.xlim([min(x), max(x)])

    legend = plt.legend(loc='best', shadow=False, fontsize='small')

    #plt.savefig('equation_error.eps', format='eps', dpi=1000)

    plt.grid(True)

    plt.show()

    plt.close()

    #==========================================================================

    #==========================================================================