commit 4148890596602c8f98cb11e65f22da2da0dfbefc Author: SchrodingerError Date: Fri Oct 18 20:03:55 2024 -0500 Homework 3 diff --git a/HW3/Bar.py b/HW3/Bar.py new file mode 100644 index 0000000..e578b98 --- /dev/null +++ b/HW3/Bar.py @@ -0,0 +1,24 @@ +# Bar.py + +from numpy import pi + +class Bar(): + def __init__(self, radius:float=0.1, k:float=0.5, h:float=0.0025): + self.radius:float = radius + self.k:float = k + self.h:float = h + + self.update_properties() + + def update_properties(self): + self.area:float = pi*self.radius**2 + self.p: float = 2*pi*self.radius + self.alpha = ((self.h*self.p)/(self.k*self.area))**0.5 + + def set_k(self, k:float): + self.k = k + self.update_properties() + +def set_k_ratio(ratio: float, bar1:'Bar', bar2:'Bar'): + '''Sets the value of bar2.k = ratio*bar1.k and updates internal properties''' + bar2.set_k(bar1.k*ratio) diff --git a/HW3/Derivations.pdf b/HW3/Derivations.pdf new file mode 100644 index 0000000..64602d4 Binary files /dev/null and b/HW3/Derivations.pdf differ diff --git a/HW3/Judson_Upchurch_HW3.pdf b/HW3/Judson_Upchurch_HW3.pdf new file mode 100644 index 0000000..1c633a0 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3.pdf differ diff --git a/HW3/Judson_Upchurch_HW3.zip b/HW3/Judson_Upchurch_HW3.zip new file mode 100644 index 0000000..4cdb4a8 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3.zip differ diff --git a/HW3/Judson_Upchurch_HW3/Bar.py b/HW3/Judson_Upchurch_HW3/Bar.py new file mode 100644 index 0000000..e578b98 --- /dev/null +++ b/HW3/Judson_Upchurch_HW3/Bar.py @@ -0,0 +1,24 @@ +# Bar.py + +from numpy import pi + +class Bar(): + def __init__(self, radius:float=0.1, k:float=0.5, h:float=0.0025): + self.radius:float = radius + self.k:float = k + self.h:float = h + + self.update_properties() + + def update_properties(self): + self.area:float = pi*self.radius**2 + self.p: float = 2*pi*self.radius + self.alpha = ((self.h*self.p)/(self.k*self.area))**0.5 + + def set_k(self, k:float): + self.k = k + self.update_properties() + +def set_k_ratio(ratio: float, bar1:'Bar', bar2:'Bar'): + '''Sets the value of bar2.k = ratio*bar1.k and updates internal properties''' + bar2.set_k(bar1.k*ratio) diff --git a/HW3/Judson_Upchurch_HW3/Judson_Upchurch_HW3.pdf b/HW3/Judson_Upchurch_HW3/Judson_Upchurch_HW3.pdf new file mode 100644 index 0000000..1c633a0 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/Judson_Upchurch_HW3.pdf differ diff --git a/HW3/Judson_Upchurch_HW3/case1.py b/HW3/Judson_Upchurch_HW3/case1.py new file mode 100644 index 0000000..4021a09 --- /dev/null +++ b/HW3/Judson_Upchurch_HW3/case1.py @@ -0,0 +1,244 @@ +# case1.py + +import numpy as np +import matplotlib.pyplot as plt + +from Bar import Bar, set_k_ratio + +def get_analytical_constants(bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Solve for C1, D1, C2, and D2 + Returns np.array([C1, D1, C2, D2])''' + epsilon = (bar1.k*bar1.area*bar1.alpha) / (bar2.k*bar2.area*bar2.alpha) + A = np.array([ + #C1 D1 C2 D2 + [0, 1, 0, 0], + [0, 0, np.sinh(bar2.alpha*l), np.cosh(bar2.alpha*l)], + [np.sinh(bar1.alpha*x_bar), 0, -1*np.sinh(bar2.alpha*x_bar), -1*np.cosh(bar2.alpha*x_bar)], + [epsilon*np.cosh(bar1.alpha*x_bar), 0, -1*np.cosh(bar2.alpha*x_bar), -1*np.sinh(bar2.alpha*x_bar)] + ]) + + B = np.array([ + 0, + 100, + 0, + 0 + ]) + + return np.linalg.solve(A, B) + +def get_analytical_datapoints(x_points, bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Generates the temperature distribution given x_points. + Returns np.array of tempature values''' + constants = get_analytical_constants(bar1=bar1,bar2=bar2,x_bar=x_bar,l=1) + C1, D1, C2, D2 = constants[0], constants[1], constants[2], constants[3] + + # Two temp functions for the left and right side of the interface + left_temp = lambda x: C1*np.sinh(bar1.alpha*x) + D1*np.cosh(bar1.alpha*x) + right_temp = lambda x: C2*np.sinh(bar2.alpha*x) + D2*np.cosh(bar2.alpha*x) + + temp_values = np.array([]) + for x in x_points: + if x <= x_bar: + temp = left_temp(x) + else: + temp = right_temp(x) + temp_values = np.append(temp_values, temp) + + return temp_values + +def get_analytical_heat_convection(x, bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Gets the analytical heat convection at x from heat flux conservation''' + constants = get_analytical_constants(bar1=bar1,bar2=bar2,x_bar=x_bar,l=1) + C1, D1, C2, D2 = constants[0], constants[1], constants[2], constants[3] + + t_prime = C2*bar2.alpha*np.cosh(bar2.alpha*x) + D2*bar2.alpha*np.sinh(bar2.alpha*x) + + q_dot = -bar2.k*bar2.area*t_prime + return q_dot + +def fdm_array(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1): + '''num_sections = total sections to cut the bar. Note that this is NOT evenly distributed over the entire + bar. Half of the sections are left of x_bar, half are to the right. + Returns the A matrix of the FDM equation Ax = B''' + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + k1 = 2 + bar1.alpha**2 * dx1**2 + k2 = 2 + bar2.alpha**2 * dx2**2 + + beta1 = (bar1.k * bar1.area) / (dx1) + beta2 = (bar2.k * bar2.area) / (dx2) + + #rho = beta1 - beta2 - (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + rho = beta1 + beta2 + (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + # Making the matrix + # There are num_sections - 1 rows and columns + # row 0, row n - 2, and row num_sections/2 - 1 are unique + # the rest are -1, k, -1 where k changes from k1 to k2 at num_sections/2 + + # Initialize matrix A with zeros + A = np.zeros((num_sections - 1, num_sections - 1)) + + # Fill the matrix A + for row in range(num_sections - 1): + if row == 0: + # First row for boundary at x = 0 (for example Dirichlet or Neumann condition) + A[row, row] = k1 + A[row, row+1] = -1 + elif row == num_sections // 2 - 1: + # Special row at the interface between bar1 and bar2 + A[row, row-1] = -beta1 + A[row, row] = rho # This is the rho value you computed + A[row, row+1] = -beta2 + elif row == num_sections - 2: + # Last row for boundary at x = l (again assuming Dirichlet or Neumann condition) + A[row, row-1] = -1 + A[row, row] = k2 + else: + # Interior rows for bar1 and bar2 (uniform grid) + if row < num_sections // 2 - 1: + # In bar1 region + A[row, row-1] = -1 + A[row, row] = k1 + A[row, row+1] = -1 + else: + # In bar2 region + A[row, row-1] = -1 + A[row, row] = k2 + A[row, row+1] = -1 + return A + +def fem_array(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1): + '''num_sections = total sections to cut the bar. Note that this is NOT evenly distributed over the entire + bar. Half of the sections are left of x_bar, half are to the right. + Returns the A matrix of the FEM equation Ax = B''' + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + k1 = (2/dx1 + 4*bar1.alpha**2*dx1/6) / (1/dx1 - bar1.alpha**2*dx1/6) + k2 = (2/dx2 + 4*bar2.alpha**2*dx2/6) / (1/dx2 - bar2.alpha**2*dx2/6) + + beta1 = (bar1.k * bar1.area) / (dx1) + beta2 = (bar2.k * bar2.area) / (dx2) + + rho = beta1 + beta2 + (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + # Making the matrix + # There are num_sections - 1 rows and columns + # row 0, row n - 2, and row num_sections/2 - 1 are unique + # the rest are -1, k, -1 where k changes from k1 to k2 at num_sections/2 + + # Initialize matrix A with zeros + A = np.zeros((num_sections - 1, num_sections - 1)) + + # Fill the matrix A + for row in range(num_sections - 1): + if row == 0: + # First row for boundary at x = 0 (for example Dirichlet or Neumann condition) + A[row, row] = k1 + A[row, row+1] = -1 + elif row == num_sections // 2 - 1: + # Special row at the interface between bar1 and bar2 + A[row, row-1] = -beta1 + A[row, row] = rho # This is the rho value you computed + A[row, row+1] = -beta2 + elif row == num_sections - 2: + # Last row for boundary at x = l (again assuming Dirichlet or Neumann condition) + A[row, row-1] = -1 + A[row, row] = k2 + else: + # Interior rows for bar1 and bar2 (uniform grid) + if row < num_sections // 2 - 1: + # In bar1 region + A[row, row-1] = -1 + A[row, row] = k1 + A[row, row+1] = -1 + else: + # In bar2 region + A[row, row-1] = -1 + A[row, row] = k2 + A[row, row+1] = -1 + return A + +def solve(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1, method="FDM"): + '''Solves using FDM or FEM as specified. Returns (temps, x_values)''' + if method == "FDM": + A = fdm_array(bar1, + bar2, + num_sections=num_sections, + x_bar=x_bar, + l=l) + elif method == "FEM": + A = fem_array(bar1, + bar2, + num_sections=num_sections, + x_bar=x_bar, + l=l) + + B = np.zeros((num_sections - 1)) + B[-1] = 100 + + # Add boundary conditions + interior_temps = np.linalg.solve(A, B) + all_temps = np.array([0]) + all_temps = np.append(all_temps, interior_temps) + all_temps = np.append(all_temps, 100) + + # Generate the x_values to return along with the temps + dx2 = (l - x_bar) / (num_sections / 2) + x_values = np.linspace(0, x_bar, num_sections // 2 + 1) # [0, x_bar] + x_values = np.append(x_values, np.linspace(x_bar+dx2, l, num_sections // 2)) + return all_temps, x_values + + +if __name__ == "__main__": + bar1 = Bar( + radius=0.1, + k=0.5, + h=0.25 + ) + bar2 = Bar( + radius=0.1, + k=2, + h=0.25 + ) + + set_k_ratio(0.5, bar1, bar2) + + x_bar = 0.5 + + x_values = np.linspace(0, 1, 1000) + temp_values = get_analytical_datapoints(x_points=x_values, + bar1=bar1, + bar2=bar2, + x_bar=x_bar) + + + # fdm tests + l = 1 + num_sections = 4 + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + fdm_temps, x_fdm = solve(bar1=bar1, + bar2=bar2, + num_sections=num_sections, + x_bar=x_bar, + l=1, + method="FDM") + fem_temps, x_fem = solve(bar1=bar1, + bar2=bar2, + num_sections=num_sections, + x_bar=x_bar, + l=1, + method="FEM") + + # Plot the data + plt.plot(x_values, temp_values, label='Temperature Distribution') + plt.axvline(x=x_bar, color='r', linestyle='--', label='x_bar') + plt.plot(x_fdm, fdm_temps, 'o-', label='Temperature Distribution (FDM)', color='g') + plt.plot(x_fem, fem_temps, 'o-', label='Temperature Distribution (FEM)', color='r') + plt.xlabel('x') + plt.ylabel('Temperature') + plt.title('Temperature Distribution Along the Bar') + plt.legend() + plt.grid(True) + + # Show plot + plt.show() \ No newline at end of file diff --git a/HW3/Judson_Upchurch_HW3/case2.py b/HW3/Judson_Upchurch_HW3/case2.py new file mode 100644 index 0000000..5231180 --- /dev/null +++ b/HW3/Judson_Upchurch_HW3/case2.py @@ -0,0 +1,230 @@ +# case2.py + +import numpy as np +import matplotlib.pyplot as plt + +from Bar import Bar + +def get_analytical_constants(bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Solve for C1, D1, C2, and D2 + Returns np.array([C1, D1, C2, D2])''' + epsilon = (bar1.k*bar1.area*bar1.alpha) / (bar2.k*bar2.area*bar2.alpha) + A = np.array([ + #C1 D1 C2 D2 + [1, 0, 0, 0], + [0, 0, np.sinh(bar2.alpha*l), np.cosh(bar2.alpha*l)], + [0, np.cosh(bar1.alpha*x_bar), -1*np.sinh(bar2.alpha*x_bar), -1*np.cosh(bar2.alpha*x_bar)], + [0, epsilon*np.sinh(bar1.alpha*x_bar), -1*np.cosh(bar2.alpha*x_bar), -1*np.sinh(bar2.alpha*x_bar)] + ]) + + B = np.array([ + 0, + 100, + 0, + 0 + ]) + + return np.linalg.solve(A, B) + +def get_analytical_datapoints(x_points, bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Generates the temperature distribution given x_points. + Returns np.array of tempature values''' + constants = get_analytical_constants(bar1=bar1,bar2=bar2,x_bar=x_bar,l=1) + C1, D1, C2, D2 = constants[0], constants[1], constants[2], constants[3] + + # Two temp functions for the left and right side of the interface + left_temp = lambda x: C1*np.sinh(bar1.alpha*x) + D1*np.cosh(bar1.alpha*x) + right_temp = lambda x: C2*np.sinh(bar2.alpha*x) + D2*np.cosh(bar2.alpha*x) + + temp_values = np.array([]) + for x in x_points: + if x <= x_bar: + temp = left_temp(x) + else: + temp = right_temp(x) + temp_values = np.append(temp_values, temp) + + return temp_values + +def fdm_array(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1): + '''num_sections = total sections to cut the bar. Note that this is NOT evenly distributed over the entire + bar. Half of the sections are left of x_bar, half are to the right. + Returns the A matrix of the FDM equation Ax = B''' + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + k1 = 2 + bar1.alpha**2 * dx1**2 + k2 = 2 + bar2.alpha**2 * dx2**2 + k1_prime = k1/2 + + beta1 = (bar1.k * bar1.area) / (dx1) + beta2 = (bar2.k * bar2.area) / (dx2) + + #rho = beta1 - beta2 - (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + rho = beta1 + beta2 + (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + # Making the matrix + # There are num_sections rows and columns + # row 0, row n - 2, and row num_sections/2 are unique + # the rest are -1, k, -1 where k changes from k1 to k2 at num_sections/2 + + # Initialize matrix A with zeros + A = np.zeros((num_sections, num_sections)) + + # Fill the matrix A + for row in range(num_sections): + if row == 0: + # First row for boundary at x = 0 ( Neumann condition) + A[row, row] = k1_prime + A[row, row+1] = -1 + elif row == num_sections // 2: + # Special row at the interface between bar1 and bar2 + A[row, row-1] = -beta1 + A[row, row] = rho # This is the rho value you computed + A[row, row+1] = -beta2 + elif row == num_sections - 1: + # Last row for boundary at x = l (again assuming Dirichlet) + A[row, row-1] = -1 + A[row, row] = k2 + else: + # Interior rows for bar1 and bar2 (uniform grid) + if row < num_sections // 2 - 1: + # In bar1 region + A[row, row-1] = -1 + A[row, row] = k1 + A[row, row+1] = -1 + else: + # In bar2 region + A[row, row-1] = -1 + A[row, row] = k2 + A[row, row+1] = -1 + return A + +def fem_array(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1): + '''num_sections = total sections to cut the bar. Note that this is NOT evenly distributed over the entire + bar. Half of the sections are left of x_bar, half are to the right. + Returns the A matrix of the FDM equation Ax = B''' + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + k1 = (2/dx1 + 4*bar1.alpha**2*dx1/6) / (1/dx1 - bar1.alpha**2*dx1/6) + k2 = (2/dx2 + 4*bar2.alpha**2*dx2/6) / (1/dx2 - bar2.alpha**2*dx2/6) + k1_prime = k1/2 + + beta1 = (bar1.k * bar1.area) / (dx1) + beta2 = (bar2.k * bar2.area) / (dx2) + + #rho = beta1 - beta2 - (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + rho = beta1 + beta2 + (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + # Making the matrix + # There are num_sections rows and columns + # row 0, row n - 2, and row num_sections/2 are unique + # the rest are -1, k, -1 where k changes from k1 to k2 at num_sections/2 + + # Initialize matrix A with zeros + A = np.zeros((num_sections, num_sections)) + + # Fill the matrix A + for row in range(num_sections): + if row == 0: + # First row for boundary at x = 0 ( Neumann condition) + A[row, row] = k1_prime + A[row, row+1] = -1 + elif row == num_sections // 2: + # Special row at the interface between bar1 and bar2 + A[row, row-1] = -beta1 + A[row, row] = rho # This is the rho value you computed + A[row, row+1] = -beta2 + elif row == num_sections - 1: + # Last row for boundary at x = l (again assuming Dirichlet) + A[row, row-1] = -1 + A[row, row] = k2 + else: + # Interior rows for bar1 and bar2 (uniform grid) + if row < num_sections // 2 - 1: + # In bar1 region + A[row, row-1] = -1 + A[row, row] = k1 + A[row, row+1] = -1 + else: + # In bar2 region + A[row, row-1] = -1 + A[row, row] = k2 + A[row, row+1] = -1 + return A + +def solve(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1, method="FDM"): + if method == "FDM": + A = fdm_array(bar1, + bar2, + num_sections=num_sections, + x_bar=x_bar, + l=l) + elif method == "FEM": + A = fem_array(bar1, + bar2, + num_sections=num_sections, + x_bar=x_bar, + l=l) + + B = np.zeros((num_sections)) + B[-1] = 100 + + # Append the boundary condition + interior_temps = np.linalg.solve(A, B) + all_temps = np.append(interior_temps, 100) + + # Generate the x_values to return along with the temps + dx2 = (l - x_bar) / (num_sections / 2) + x_values = np.linspace(0, x_bar, num_sections // 2 + 1) # [0, x_bar] + x_values = np.append(x_values, np.linspace(x_bar+dx2, l, num_sections // 2)) + return all_temps, x_values + + +if __name__ == "__main__": + bar1 = Bar( + radius=0.1, + k=0.5, + h=0.25 + ) + bar2 = Bar( + radius=0.1, + k=2, + h=0.25 + ) + + x_bar = 2/np.pi + + + x_values = np.linspace(0, 1, 1000) + temp_values = get_analytical_datapoints(x_points=x_values, + bar1=bar1, + bar2=bar2) + + + # fdm tests + l = 1 + num_sections = 6 + fdm_temps, x_fdm = solve(bar1=bar1, + bar2=bar2, + num_sections=num_sections, + x_bar=x_bar, + l=1, + method="FDM") + fem_temps, x_fem = solve(bar1=bar1, + bar2=bar2, + num_sections=num_sections, + x_bar=x_bar, + l=1, + method="FEM") + + # Plot the data + plt.plot(x_values, temp_values, label='Temperature Distribution') + plt.axvline(x=x_bar, color='r', linestyle='--', label='x_bar') + plt.plot(x_fdm, fdm_temps, 'o-', label='Temperature Distribution (FDM)', color='g') + plt.plot(x_fem, fem_temps, 'o-', label='Temperature Distribution (FEM)', color='r') + plt.xlabel('x') + plt.ylabel('Temperature') + plt.title('Temperature Distribution Along the Bar') + plt.legend() + plt.grid(True) + + # Show plot + plt.show() \ No newline at end of file diff --git a/HW3/Judson_Upchurch_HW3/common.py b/HW3/Judson_Upchurch_HW3/common.py new file mode 100644 index 0000000..bb54345 --- /dev/null +++ b/HW3/Judson_Upchurch_HW3/common.py @@ -0,0 +1,30 @@ +# common.py + +import numpy as np +from Bar import Bar + +def fdm_heat_extraction(t_0, t_1, dx, bar:'Bar', order=2): + '''Get the heat conduction at the point t_1 using Taylor series''' + if order == 1: + return -1 * bar.k * bar.area * (t_1 - t_0) / dx + elif order == 2: + return -1 * bar.k * bar.area * (((t_1 - t_0) / dx) + (bar.alpha**2 * dx * t_1 / 2)) + +def fem_heat_extraction(t_0, t_1, dx, bar:'Bar'): + '''Get the heat conduction at the point t_1 using FEM equation''' + term_1 = (-1/dx + bar.alpha**2*dx/6) * t_0 + term_2 = (1/dx + 2*bar.alpha**2*dx/6) * t_1 + return -1 * bar.k * bar.area * (term_1 + term_2) + +def calc_error(exact, q_1): + return np.abs((exact - q_1) / exact) + +def calc_beta(exact, q_1, q_2, dx_1, dx_2): + return np.log(np.abs((exact - q_1)/(exact - q_2))) / np.log(dx_1 / dx_2) + +def calc_extrapolated(q1, q2, q3): + '''Calcs the Richardson Extrapolation value based on 3 different meshes. + Assumes that q3 is 2x finer than q2 is 2x finer than q1''' + numerator = q1*q3 - q2**2 + denominator = q1 + q3 - 2*q2 + return numerator / denominator \ No newline at end of file diff --git a/HW3/Judson_Upchurch_HW3/constants.py b/HW3/Judson_Upchurch_HW3/constants.py new file mode 100644 index 0000000..5248508 --- /dev/null +++ b/HW3/Judson_Upchurch_HW3/constants.py @@ -0,0 +1,19 @@ +from numpy import pi + +h = 0.0025 # W / cm^2*K + +class Bar1(): + def __init__(self): + self.k = 0.5 # W / cm*K + self.r = 0.1 # cm + self.A = pi*self.r**2 # cm^2 + self.P = 2*pi*self.r # cm + +class Bar2(): + def __init__(self): + self.k = 0.5 # W / cm*K + self.r = 0.1 # cm + self.A = pi*self.r**2 # cm^2 + self.P = 2*pi*self.r # cm + + self.alpha = ((h*self.P)/(self.k*self.A))**0.5 \ No newline at end of file diff --git a/HW3/Judson_Upchurch_HW3/images/section1/x_bar_1/analytical_heat_convection_vs_k_ratio.png b/HW3/Judson_Upchurch_HW3/images/section1/x_bar_1/analytical_heat_convection_vs_k_ratio.png new file mode 100644 index 0000000..0f7c4b3 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section1/x_bar_1/analytical_heat_convection_vs_k_ratio.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section1/x_bar_1/analytical_temperature_vs_position.png b/HW3/Judson_Upchurch_HW3/images/section1/x_bar_1/analytical_temperature_vs_position.png new file mode 100644 index 0000000..bef0e63 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section1/x_bar_1/analytical_temperature_vs_position.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section1/x_bar_2/analytical_heat_convection_vs_k_ratio.png b/HW3/Judson_Upchurch_HW3/images/section1/x_bar_2/analytical_heat_convection_vs_k_ratio.png new file mode 100644 index 0000000..52186d4 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section1/x_bar_2/analytical_heat_convection_vs_k_ratio.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section1/x_bar_2/analytical_temperature_vs_position.png b/HW3/Judson_Upchurch_HW3/images/section1/x_bar_2/analytical_temperature_vs_position.png new file mode 100644 index 0000000..3e47b99 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section1/x_bar_2/analytical_temperature_vs_position.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/fdm_temperature_vs_position.png b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/fdm_temperature_vs_position.png new file mode 100644 index 0000000..af5c966 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/fdm_temperature_vs_position.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/fem_temperature_vs_position.png b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/fem_temperature_vs_position.png new file mode 100644 index 0000000..9411db2 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/fem_temperature_vs_position.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/heat_convection_vs_k_ratio.png b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/heat_convection_vs_k_ratio.png new file mode 100644 index 0000000..77c2a96 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_1/heat_convection_vs_k_ratio.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/fdm_temperature_vs_position.png b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/fdm_temperature_vs_position.png new file mode 100644 index 0000000..011d8ff Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/fdm_temperature_vs_position.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/fem_temperature_vs_position.png b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/fem_temperature_vs_position.png new file mode 100644 index 0000000..56630cc Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/fem_temperature_vs_position.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/heat_convection_vs_k_ratio.png b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/heat_convection_vs_k_ratio.png new file mode 100644 index 0000000..d040ec3 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section2/x_bar_2/heat_convection_vs_k_ratio.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section3/x_bar_1/heat_convergences.png b/HW3/Judson_Upchurch_HW3/images/section3/x_bar_1/heat_convergences.png new file mode 100644 index 0000000..fa8bcdd Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section3/x_bar_1/heat_convergences.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section3/x_bar_1/temp_convergences.png b/HW3/Judson_Upchurch_HW3/images/section3/x_bar_1/temp_convergences.png new file mode 100644 index 0000000..1e47f92 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section3/x_bar_1/temp_convergences.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section3/x_bar_2/heat_convergences.png b/HW3/Judson_Upchurch_HW3/images/section3/x_bar_2/heat_convergences.png new file mode 100644 index 0000000..78613b6 Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section3/x_bar_2/heat_convergences.png differ diff --git a/HW3/Judson_Upchurch_HW3/images/section3/x_bar_2/temp_convergences.png b/HW3/Judson_Upchurch_HW3/images/section3/x_bar_2/temp_convergences.png new file mode 100644 index 0000000..3aa08aa Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/images/section3/x_bar_2/temp_convergences.png differ diff --git a/HW3/Judson_Upchurch_HW3/main.py b/HW3/Judson_Upchurch_HW3/main.py new file mode 100644 index 0000000..e36efe4 --- /dev/null +++ b/HW3/Judson_Upchurch_HW3/main.py @@ -0,0 +1,722 @@ +# main.py + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import math + +import case1 +import case2 +import common +import Bar + + +bar1 = Bar.Bar( + radius=0.1, + k=0.5, + h=0.25 +) +bar2 = Bar.Bar( + radius=0.1, + k=2, + h=0.25 +) + +data_dict = {} +data_dict["k_ratios"] = np.array([1/16, 1/8, 1/4, 1/2, 1, 2, 4, 8, 16]) +data_dict["x_bar_values"] = np.array([1/2, 2/np.pi]) +data_dict["length"] = 1 + +def section_1(): + '''Analytical solutions with varying k_ratio''' + print("Section 1: Analytical") + + l = data_dict["length"] + + # Calculate the results for every k ratio for x_bar_1 + x_bar_1 = data_dict["x_bar_values"][0] + + data_dict["section1"] = {} + data_dict["section1"]["x_bar_1"] = {} + data_dict["section1"]["x_bar_1"]["x_values_fine"] = np.linspace(0, l, 50) + data_dict["section1"]["x_bar_1"]["x_values_fine"] = np.sort(np.append(data_dict["section1"]["x_bar_1"]["x_values_fine"], x_bar_1)) + data_dict["section1"]["x_bar_1"]["x_values_course"] = np.linspace(0, l, 9) + data_dict["section1"]["x_bar_1"]["x_values_course"] = np.sort(np.append(data_dict["section1"]["x_bar_1"]["x_values_course"], x_bar_1)) + + results = [] + results_course = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + result = case1.get_analytical_datapoints(data_dict["section1"]["x_bar_1"]["x_values_fine"], + bar1=bar1, + bar2=bar2, + x_bar=x_bar_1, + l=l) + result_course = case1.get_analytical_datapoints(data_dict["section1"]["x_bar_1"]["x_values_course"], + bar1=bar1, + bar2=bar2, + x_bar=x_bar_1, + l=l) + results.append(result) + results_course.append(result_course) + + data_dict["section1"]["x_bar_1"]["temp_results"] = np.array(results) + + df_x_bar_1 = pd.DataFrame(results_course, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in data_dict["section1"]["x_bar_1"]["x_values_course"]]) + + print("Analytical Temperature Results x_bar = 0.5:") + print(df_x_bar_1.to_string()) + print("\n" * 2) + + # Plotting x_bar_1 + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section1"]["x_bar_1"]["x_values_fine"], data_dict["section1"]["x_bar_1"]["temp_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_1, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section1/x_bar_1/analytical_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + # Calculate the temperature results for every k ratio for x_bar_1 + x_bar_2 = data_dict["x_bar_values"][1] + + data_dict["section1"]["x_bar_2"] = {} + data_dict["section1"]["x_bar_2"]["x_values_fine"] = np.linspace(0, l, 50) + data_dict["section1"]["x_bar_2"]["x_values_fine"] = np.sort(np.append(data_dict["section1"]["x_bar_2"]["x_values_fine"], x_bar_2)) + data_dict["section1"]["x_bar_2"]["x_values_course"] = np.linspace(0, l, 9) + data_dict["section1"]["x_bar_2"]["x_values_course"] = np.sort(np.append(data_dict["section1"]["x_bar_2"]["x_values_course"], x_bar_2)) + + results = [] + results_course = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) + result = case1.get_analytical_datapoints(data_dict["section1"]["x_bar_2"]["x_values_fine"], + bar1=bar1, + bar2=bar2, + x_bar=x_bar_2, + l=l) + result_course = case1.get_analytical_datapoints(data_dict["section1"]["x_bar_2"]["x_values_course"], + bar1=bar1, + bar2=bar2, + x_bar=x_bar_2, + l=l) + results.append(result) + results_course.append(result_course) + + data_dict["section1"]["x_bar_2"]["temp_results"] = np.array(results) + + df_x_bar_2 = pd.DataFrame(results_course, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in data_dict["section1"]["x_bar_2"]["x_values_course"]]) + + print("Analytical Temperature Results x_bar = 2/pi:") + print(df_x_bar_2.to_string()) + print("\n" * 2) + + # Plotting x_bar_2 + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section1"]["x_bar_2"]["x_values_fine"], data_dict["section1"]["x_bar_2"]["temp_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_2, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section1/x_bar_2/analytical_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + + + # Calculate the analytical heat convection leaving the rod at x=l for varying k ratio + + results_1 = [] + results_2 = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) + result_1 = case1.get_analytical_heat_convection(x=l, + bar1=bar1, + bar2=bar2, + x_bar=x_bar_1, + l=l) + result_2 = case1.get_analytical_heat_convection(x=l, + bar1=bar1, + bar2=bar2, + x_bar=x_bar_2, + l=l) + results_1.append([result_1]) + results_2.append([result_2]) + + data_dict["section1"]["x_bar_1"]["heat_results"] = np.array(results_1) + data_dict["section1"]["x_bar_2"]["heat_results"] = np.array(results_2) + + df_x_bar_1 = pd.DataFrame(results_1, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + df_x_bar_2 = pd.DataFrame(results_2, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + + print("Analytical Heat Convection Results x_bar = 0.5:") + print(df_x_bar_1.to_string()) + print("\n" * 2) + + print("Analytical Heat Convection Results x_bar = 2/pi:") + print(df_x_bar_2.to_string()) + print("\n" * 2) + + # Plotting x_bar_1 heat convection + plt.figure(figsize=(10, 6)) + plt.plot(data_dict["k_ratios"], data_dict["section1"]["x_bar_1"]["heat_results"], label=f'Heat Convection') + plt.xlabel('k2/k1') + plt.ylabel('Q_Dot (W)') + plt.legend() + plt.grid(True) + plt.savefig('images/section1/x_bar_1/analytical_heat_convection_vs_k_ratio.png', dpi=300) + #plt.show() + plt.clf() + + # Plotting x_bar_2 heat convection + plt.figure(figsize=(10, 6)) + plt.plot(data_dict["k_ratios"], data_dict["section1"]["x_bar_2"]["heat_results"], label=f'Heat Convection') + plt.xlabel('k2/k1') + plt.ylabel('Q_Dot (W)') + plt.legend() + plt.grid(True) + plt.savefig('images/section1/x_bar_2/analytical_heat_convection_vs_k_ratio.png', dpi=300) + #plt.show() + plt.clf() + + +def section_2(): + '''FDM and FEM with varying k_ratio''' + print("Section 2: FDM and FEM with Varying k_ratio") + + l = data_dict["length"] + + data_dict["section2"] = {} + data_dict["section2"]["num_sections"] = 8 + + + # Calculate the results for every k ratio for x_bar_1 + x_bar_1 = data_dict["x_bar_values"][0] + + # For every k_ratio, calculate the FDM and FEM temperature results + fdm_results = [] + fem_results = [] + x_fdm = [] + x_fem = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + fdm_temps, x_fdm = case1.solve(bar1=bar1, + bar2=bar2, + num_sections=data_dict["section2"]["num_sections"], + x_bar=x_bar_1, + l=l, + method="FDM") + fem_temps, x_fem = case1.solve(bar1=bar1, + bar2=bar2, + num_sections=data_dict["section2"]["num_sections"], + x_bar=x_bar_1, + l=l, + method="FEM") + + fdm_results.append(fdm_temps) + fem_results.append(fem_temps) + + data_dict["section2"]["x_bar_1"] = {} + data_dict["section2"]["x_bar_1"]["x_values"] = x_fdm # fdm and fem use same x_values + data_dict["section2"]["x_bar_1"]["fdm_results"] = np.array(fdm_results) + data_dict["section2"]["x_bar_1"]["fem_results"] = np.array(fem_results) + + df_fdm = pd.DataFrame(fdm_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in x_fdm]) + df_fem = pd.DataFrame(fem_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in x_fem]) + + print("FDM Temperature Results x_bar = 0.5:") + print(df_fdm.to_string()) + print("\n" * 2) + print("FEM Temperature Results x_bar = 0.5:") + print(df_fem.to_string()) + print("\n" * 2) + + # Now that the data is gathered for FDM and FEM, plot it + # Plotting x_bar_1, FDM + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section2"]["x_bar_1"]["x_values"], data_dict["section2"]["x_bar_1"]["fdm_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_1, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_1/fdm_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + # Plotting x_bar_1, FEM + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section2"]["x_bar_1"]["x_values"], data_dict["section2"]["x_bar_1"]["fem_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_1, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_1/fem_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + + # Calculate the results for every k ratio for x_bar_2 + x_bar_2 = data_dict["x_bar_values"][1] + + # For every k_ratio, calculate the FDM and FEM temperature results + fdm_results = [] + fem_results = [] + x_fdm = [] + x_fem = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + fdm_temps, x_fdm = case1.solve(bar1=bar1, + bar2=bar2, + num_sections=data_dict["section2"]["num_sections"], + x_bar=x_bar_2, + l=l, + method="FDM") + fem_temps, x_fem = case1.solve(bar1=bar1, + bar2=bar2, + num_sections=data_dict["section2"]["num_sections"], + x_bar=x_bar_2, + l=l, + method="FEM") + + fdm_results.append(fdm_temps) + fem_results.append(fem_temps) + + data_dict["section2"]["x_bar_2"] = {} + data_dict["section2"]["x_bar_2"]["x_values"] = x_fdm # fdm and fem use same x_values + data_dict["section2"]["x_bar_2"]["fdm_results"] = np.array(fdm_results) + data_dict["section2"]["x_bar_2"]["fem_results"] = np.array(fem_results) + + df_fdm = pd.DataFrame(fdm_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in x_fdm]) + df_fem = pd.DataFrame(fem_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in x_fem]) + + print("FDM Temperature Results x_bar = 2/pi:") + print(df_fdm.to_string()) + print("\n" * 2) + print("FEM Temperature Results x_bar = 2/pi:") + print(df_fem.to_string()) + print("\n" * 2) + + # Plotting x_bar_2, FDM + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section2"]["x_bar_2"]["x_values"], data_dict["section2"]["x_bar_2"]["fdm_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_2, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_2/fdm_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + # Plotting x_bar_2, FEM + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section2"]["x_bar_2"]["x_values"], data_dict["section2"]["x_bar_2"]["fem_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_2, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_2/fem_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + + + # After calculating temperature values, extract the heat convection at x=l + + # x_bar_1 + # for every k ratio, calculate the heat convection from fdm and fem temp results + fdm_heat_results = [] + fem_heat_results = [] + for i, k_ratio in enumerate(data_dict["k_ratios"]): + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + dx = data_dict["section2"]["x_bar_1"]["x_values"][-1] - data_dict["section2"]["x_bar_1"]["x_values"][-2] + (fdm_t0, fdm_t1) = data_dict["section2"]["x_bar_1"]["fdm_results"][i][-2:] + fdm_heat_result = common.fdm_heat_extraction(fdm_t0, fdm_t1, dx, bar2, order=2) + fdm_heat_results.append(fdm_heat_result) + + (fem_t0, fem_t1) = data_dict["section2"]["x_bar_1"]["fem_results"][i][-2:] + fem_heat_result = common.fem_heat_extraction(fem_t0, fem_t1, dx, bar2) + fem_heat_results.append(fem_heat_result) + + + data_dict["section2"]["x_bar_1"]["fdm_heat_results"] = np.array(fdm_heat_results) + data_dict["section2"]["x_bar_1"]["fem_heat_results"] = np.array(fem_heat_results) + + df_fdm = pd.DataFrame(fdm_heat_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + df_fem = pd.DataFrame(fem_heat_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + + print("FDM Heat Convection Results x_bar = 0.5:") + print(df_fdm.to_string()) + print("\n" * 2) + + print("FEM Heat Convection Results x_bar = 0.5:") + print(df_fem.to_string()) + print("\n" * 2) + + # x_bar_2 + # for every k ratio, calculate the heat convection from fdm and fem temp results + fdm_heat_results = [] + fem_heat_results = [] + for i, k_ratio in enumerate(data_dict["k_ratios"]): + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + dx = data_dict["section2"]["x_bar_2"]["x_values"][-1] - data_dict["section2"]["x_bar_2"]["x_values"][-2] + (fdm_t0, fdm_t1) = data_dict["section2"]["x_bar_2"]["fdm_results"][i][-2:] + fdm_heat_result = common.fdm_heat_extraction(fdm_t0, fdm_t1, dx, bar2, order=2) + fdm_heat_results.append(fdm_heat_result) + + (fem_t0, fem_t1) = data_dict["section2"]["x_bar_2"]["fem_results"][i][-2:] + fem_heat_result = common.fem_heat_extraction(fem_t0, fem_t1, dx, bar2) + fem_heat_results.append(fem_heat_result) + + + data_dict["section2"]["x_bar_2"]["fdm_heat_results"] = np.array(fdm_heat_results) + data_dict["section2"]["x_bar_2"]["fem_heat_results"] = np.array(fem_heat_results) + + df_fdm = pd.DataFrame(fdm_heat_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + df_fem = pd.DataFrame(fem_heat_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + + print("FDM Heat Convection Results x_bar = 2/pi:") + print(df_fdm.to_string()) + print("\n" * 2) + + print("FEM Heat Convection Results x_bar = 2/pi:") + print(df_fem.to_string()) + print("\n" * 2) + + # Plotting x_bar_1 heat convection + plt.figure(figsize=(10, 6)) + plt.plot(data_dict["k_ratios"], data_dict["section2"]["x_bar_1"]["fdm_heat_results"], label=f'FDM Heat Convection') + plt.plot(data_dict["k_ratios"], data_dict["section2"]["x_bar_1"]["fem_heat_results"], label=f'FEM Heat Convection') + #plt.plot(data_dict["k_ratios"], data_dict["section1"]["x_bar_1"]["heat_results"], label=f'Analytical Heat Convection') + plt.xlabel('k2/k1') + plt.ylabel('Q_Dot (W)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_1/heat_convection_vs_k_ratio.png', dpi=300) + #plt.show() + plt.clf() + + # Plotting x_bar_2 heat convection + plt.figure(figsize=(10, 6)) + plt.plot(data_dict["k_ratios"], data_dict["section2"]["x_bar_2"]["fdm_heat_results"], label=f'FDM Heat Convection') + plt.plot(data_dict["k_ratios"], data_dict["section2"]["x_bar_2"]["fem_heat_results"], label=f'FEM Heat Convection') + #plt.plot(data_dict["k_ratios"], data_dict["section1"]["x_bar_2"]["heat_results"], label=f'Analytical Heat Convection') + plt.xlabel('k2/k1') + plt.ylabel('Q_Dot (W)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_2/heat_convection_vs_k_ratio.png', dpi=300) + #plt.show() + plt.clf() + + + +def section_3(): + '''Convergence of FDM and FEM of interface temperature and heat convection''' + print("Section 3: Convergence of FDM and FEM with Varying k Ratio") + + l = data_dict["length"] + + data_dict["section3"] = {} + data_dict["section3"]["num_sections"] = [4, 8, 16, 32, 64, 128] + + for k, x_bar_1 in enumerate(data_dict["x_bar_values"]): + # Calculate the number of rows needed for two columns of subplots + num_k_ratios = len(data_dict["k_ratios"]) + num_rows = math.ceil(num_k_ratios / 2) + + # Define the figure for temperature convergence (two columns) + fig_temp, axs_temp = plt.subplots(num_rows, 2, figsize=(15, 5 * num_rows)) + + # Define the figure for heat convergence (two columns) + fig_heat, axs_heat = plt.subplots(num_rows, 2, figsize=(15, 5 * num_rows)) + + # Flatten the axs arrays for easier indexing + axs_temp = axs_temp.flatten() + axs_heat = axs_heat.flatten() + + # For every k ratio, iterate through every number of points, calculate convergence (compared to analytical and extrapolated), print and add a subplot + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + + # Get the analytical data + analy_interface_temp = case1.get_analytical_datapoints([x_bar_1], bar1, bar2, x_bar_1, l) + analy_heat = case1.get_analytical_heat_convection(l, bar1, bar2, x_bar_1, l) + + # Setup data arrays + fdm_interface_temps = [] + fem_interface_temps = [] + fdm_extrap_temps = [None, None] + fem_extrap_temps = [None, None] + fdm_heats = [] + fem_heats = [] + fdm_extrap_heats = [None, None] + fem_extrap_heats = [None, None] + fdm_interface_temp_errors = [] + fem_interface_temp_errors = [] + fdm_interface_temp_betas = [[None]] + fem_interface_temp_betas = [[None]] + fdm_heat_errors = [] + fem_heat_errors = [] + fdm_heat_betas = [[None]] + fem_heat_betas = [[None]] + fdm_interface_temp_extrap_errors = [None, None] + fem_interface_temp_extrap_errors = [None, None] + fdm_interface_temp_extrap_betas = [None, None] + fem_interface_temp_extrap_betas = [None, None] + fdm_heat_extrap_errors = [None, None] + fem_heat_extrap_errors = [None, None] + fdm_heat_extrap_betas = [None, None] + fem_heat_extrap_betas = [None, None] + + # Get the fdm and fem data + dx = 0 + for i, num_sections in enumerate(data_dict["section3"]["num_sections"]): + fdm_temps, x_fdm = case1.solve(bar1=bar1,bar2=bar2,num_sections=num_sections,x_bar=x_bar_1,l=l,method="FDM") + fem_temps, x_fem = case1.solve(bar1=bar1,bar2=bar2,num_sections=num_sections,x_bar=x_bar_1,l=l,method="FEM") + # Get the interface temperature + fdm_interface_temp = fdm_temps[num_sections // 2] + fem_interface_temp = fem_temps[num_sections // 2] + fdm_interface_temps.append(fdm_interface_temp) + fem_interface_temps.append(fem_interface_temp) + + # Get the 2 last temperature to calculate heat convection + dx_prev = dx + dx = x_fdm[-1] - x_fdm[-2] + (fdm_t0, fdm_t1) = fdm_temps[-2:] + fdm_heat = common.fdm_heat_extraction(fdm_t0, fdm_t1, dx, bar2, order=2) + (fem_t0, fem_t1) = fem_temps[-2:] + fem_heat = common.fem_heat_extraction(fem_t0, fem_t1, dx, bar2) + fdm_heats.append(fdm_heat) + fem_heats.append(fem_heat) + + # Calculated the Richardson extrpolation of interface temperature and heat convection + if i >= 2: + fdm_extrap_temp = common.calc_extrapolated(fdm_interface_temps[-3], fdm_interface_temps[-2], fdm_interface_temps[-1]) + fem_extrap_temp = common.calc_extrapolated(fem_interface_temps[-3], fem_interface_temps[-2], fem_interface_temps[-1]) + fdm_extrap_heat = common.calc_extrapolated(fdm_heats[-3], fdm_heats[-2], fdm_heats[-1]) + fem_extrap_heat = common.calc_extrapolated(fem_heats[-3], fem_heats[-2], fem_heats[-1]) + fdm_extrap_temps.append(fdm_extrap_temp) + fem_extrap_temps.append(fem_extrap_temp) + fdm_extrap_heats.append(fdm_extrap_heat) + fem_extrap_heats.append(fem_extrap_heat) + + # Calculate the errors in reference to extrapolated values + fdm_interface_temp_extrap_error = common.calc_error(fdm_extrap_temp, fdm_interface_temp) + fem_interface_temp_extrap_error = common.calc_error(fem_extrap_temp, fem_interface_temp) + fdm_heat_extrap_error = common.calc_error(fdm_extrap_heat, fdm_heat) + fem_heat_extrap_error = common.calc_error(fem_extrap_heat, fem_heat) + fdm_interface_temp_extrap_errors.append(fdm_interface_temp_extrap_error) + fem_interface_temp_extrap_errors.append(fem_interface_temp_extrap_error) + fdm_heat_extrap_errors.append(fdm_heat_extrap_error) + fem_heat_extrap_errors.append(fem_heat_extrap_error) + + # Calculate the betas in reference to extrapolated values + fdm_interface_temp_extrap_beta = common.calc_beta(fdm_extrap_temp, fdm_interface_temps[-1], fdm_interface_temps[-2], dx, dx_prev) + fem_interface_temp_extrap_beta = common.calc_beta(fem_extrap_temp, fem_interface_temps[-1], fem_interface_temps[-2], dx, dx_prev) + fdm_heat_extrap_beta = common.calc_beta(fdm_extrap_heat, fdm_heats[-1], fdm_heats[-2], dx, dx_prev) + fem_heat_extrap_beta = common.calc_beta(fem_extrap_heat, fem_heats[-1], fem_heats[-2], dx, dx_prev) + fdm_interface_temp_extrap_betas.append(fdm_interface_temp_extrap_beta) + fem_interface_temp_extrap_betas.append(fem_interface_temp_extrap_beta) + fdm_heat_extrap_betas.append(fdm_heat_extrap_beta) + fem_heat_extrap_betas.append(fem_heat_extrap_beta) + + + # Calculate the error and convergence rates for fdm temp, fem temp, fdm heat, and fem heat + fdm_interface_temp_error = common.calc_error(analy_interface_temp, fdm_interface_temp) + fem_interface_temp_error = common.calc_error(analy_interface_temp, fem_interface_temp) + fdm_heat_error = common.calc_error(analy_heat, fdm_heat) + fem_heat_error = common.calc_error(analy_heat, fem_heat) + fdm_interface_temp_errors.append(fdm_interface_temp_error) + fem_interface_temp_errors.append(fem_interface_temp_error) + fdm_heat_errors.append(fdm_heat_error) + fem_heat_errors.append(fem_heat_error) + + if i >= 1: # if i == 0 then we cannot calculate convergence + fdm_interface_temp_beta = common.calc_beta(analy_interface_temp, fdm_interface_temps[-1], fdm_interface_temps[-2], dx, dx_prev) + fem_interface_temp_beta = common.calc_beta(analy_interface_temp, fem_interface_temps[-1], fem_interface_temps[-2], dx, dx_prev) + fdm_heat_beta = common.calc_beta(analy_heat, fdm_heats[-1], fdm_heats[-2], dx, dx_prev) + fem_heat_beta = common.calc_beta(analy_heat, fem_heats[-1], fem_heats[-2], dx, dx_prev) + fdm_interface_temp_betas.append(fdm_interface_temp_beta) + fem_interface_temp_betas.append(fem_interface_temp_beta) + fdm_heat_betas.append(fdm_heat_beta) + fem_heat_betas.append(fem_heat_beta) + + + # Print the interface temps for this k ratio + table_data = [] + for i in range(len(data_dict["section3"]["num_sections"])): + table_data.append([ + analy_interface_temp[0], # True T + + fdm_interface_temps[i], # fdm result + fdm_extrap_temps[i], + fdm_interface_temp_errors[i][0], # fdm % error in reference to analytical + fdm_interface_temp_betas[i][0], # fdm beta + fdm_interface_temp_extrap_errors[i], # fdm % error in reference to extrapolated value + fdm_interface_temp_extrap_betas[i], # fdm beta + + fem_interface_temps[i], # fem result + fem_extrap_temps[i], + fem_interface_temp_errors[i][0], # fem % error + fem_interface_temp_betas[i][0], # fem beta + fem_interface_temp_extrap_errors[i], # fem extrap % error + fem_interface_temp_extrap_betas[i], # fem beta + ]) + + columns = ['True T', 'FDM T', 'FDM Extrap T', 'FDM Exact % Error', 'FDM Exact B', 'FDM Extrap % Error', 'FDM Extrap B', 'FEM T', 'FEM Extrap T', 'FEM Exact % Error', 'FEM Exact B', 'FEM Extrap % Error', 'FEM Extrap B'] + df = pd.DataFrame(table_data, index=[f'N = {i}' for i in data_dict["section3"]["num_sections"]], columns=columns) + + print(f"Convergence of FDM and FEM Temperature at x_bar = {x_bar_1:.3f} and k_ratio = {k_ratio}") + print(df.to_string()) + print("\n" * 2) + + # Print the heat convection results for this k ratio + table_data = [] + for i in range(len(data_dict["section3"]["num_sections"])): + table_data.append([ + analy_heat, # True Q + fdm_heats[i], # fdm result + fdm_extrap_heats[i], # fdm extrapolated heat + fdm_heat_errors[i], # fdm % error in reference to analytical + fdm_heat_betas[i], # fdm beta in reference to analytical + fdm_heat_extrap_errors[i], # fdm % error in reference to extrapolated value + fdm_heat_extrap_betas[i], # fdm beta + + fem_heats[i], # fem result + fem_extrap_heats[i], # fem extrapolated heat + fem_heat_errors[i], # fem % error + fem_heat_betas[i], # fem beta + fem_heat_extrap_errors[i], # fem % error + fem_heat_extrap_betas[i] # fem beta + ]) + + columns = ['True Q', 'FDM Q', 'FDM Extrap Q', 'FDM Exact % Error', 'FDM Exact B', 'FDM Extrap % Error', 'FDM Extrap B', 'FEM Q', 'FEM Extrap Q', 'FEM Exact % Error', 'FEM Exact B', 'FEM Extrap % Error', 'FEM Extrap B'] + df = pd.DataFrame(table_data, index=[f'N = {i}' for i in data_dict["section3"]["num_sections"]], columns=columns) + + print(f"Convergence of FDM and FEM Heat Convection with x_bar = {x_bar_1:.3f} and k_ratio = {k_ratio}") + print(df.to_string()) + print("\n" * 2) + + # Add a subplot of the interface temp convergence for this k ratio + ax_temp = axs_temp[idx] + + # Data to plot for temperature convergence + num_sections = data_dict["section3"]["num_sections"] + fdm_exact_errors = [e[0] for e in fdm_interface_temp_errors] + fem_exact_errors = [e[0] for e in fem_interface_temp_errors] + fdm_extrap_errors = fdm_interface_temp_extrap_errors + fem_extrap_errors = fem_interface_temp_extrap_errors + + # Plotting temperature lines + ax_temp.plot(num_sections, fdm_exact_errors, label="FDM Exact") + ax_temp.plot(num_sections, fem_exact_errors, label="FEM Exact") + ax_temp.plot(num_sections, fdm_extrap_errors, label="FDM Extrap") + ax_temp.plot(num_sections, fem_extrap_errors, label="FEM Extrap") + + ax_temp.set_xscale("log") + ax_temp.set_yscale("log") + ax_temp.set_xlabel("Number of Sections (N)") + ax_temp.set_ylabel("% Error") + ax_temp.set_title(f"k_ratio = {k_ratio}") + ax_temp.grid(True) + ax_temp.legend() + + # Add a subplot of the heat convergence for this k ratio + ax_heat = axs_heat[idx] + + # Data to plot for heat convergence + fdm_exact_heat_errors = fdm_heat_errors + fem_exact_heat_errors = fem_heat_errors + fdm_extrap_heat_errors = fdm_heat_extrap_errors + fem_extrap_heat_errors = fem_heat_extrap_errors + + # Plotting heat lines + ax_heat.plot(num_sections, fdm_exact_heat_errors, label="FDM Exact") + ax_heat.plot(num_sections, fem_exact_heat_errors, label="FEM Exact") + ax_heat.plot(num_sections, fdm_extrap_heat_errors, label="FDM Extrap") + ax_heat.plot(num_sections, fem_extrap_heat_errors, label="FEM Extrap") + + ax_heat.set_xscale("log") + ax_heat.set_yscale("log") + ax_heat.set_xlabel("Number of Sections (N)") + ax_heat.set_ylabel("% Error") + ax_heat.set_title(f"k_ratio = {k_ratio}") + ax_heat.grid(True) + ax_heat.legend() + + # Hide any unused subplots (in case of an odd number of k_ratios) + for i in range(num_k_ratios, len(axs_temp)): + fig_temp.delaxes(axs_temp[i]) + fig_heat.delaxes(axs_heat[i]) + + # Adjust layout for better spacing + fig_temp.tight_layout(rect=[0, 0.03, 1, 0.95]) + fig_heat.tight_layout(rect=[0, 0.03, 1, 0.95]) + + # Save the plots + if k == 0: + fig_temp.savefig("images/section3/x_bar_1/temp_convergences.png", dpi=300) + fig_heat.savefig("images/section3/x_bar_1/heat_convergences.png", dpi=300) + else: + fig_temp.savefig("images/section3/x_bar_2/temp_convergences.png", dpi=300) + fig_heat.savefig("images/section3/x_bar_2/heat_convergences.png", dpi=300) + + + +def main(): + # Make directories for plots + import os + import shutil + + base_dir = "images" + sub_dirs = ["section1", "section2", "section3"] + nested_dirs = ["x_bar_1", "x_bar_2"] + + # Create the base directory if it doesn't exist + if not os.path.exists(base_dir): + os.mkdir(base_dir) + + # Create or clear subdirectories and their nested directories + for sub_dir in sub_dirs: + section_path = os.path.join(base_dir, sub_dir) + if os.path.exists(section_path): + # Remove all contents of the directory + shutil.rmtree(section_path) + # Create the section directory + os.makedirs(section_path) + + # Create nested directories within each section + for nested_dir in nested_dirs: + nested_path = os.path.join(section_path, nested_dir) + os.makedirs(nested_path) + + # import sys + + # # Redirect all print statements to a file + # sys.stdout = open("output.txt", "w", encoding="utf-8") + + section_1() # Analytical solutions + section_2() # FDM and FEM temperature and heat convection + section_3() # Convergence of FDM and FEM temperature and heat convection + + # # Restore original stdout and close the file + # sys.stdout.close() + # sys.stdout = sys.__stdout__ + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/HW3/Judson_Upchurch_HW3/output.txt b/HW3/Judson_Upchurch_HW3/output.txt new file mode 100644 index 0000000..722312c Binary files /dev/null and b/HW3/Judson_Upchurch_HW3/output.txt differ diff --git a/HW3/__pycache__/Bar.cpython-312.pyc b/HW3/__pycache__/Bar.cpython-312.pyc new file mode 100644 index 0000000..5ed5bc9 Binary files /dev/null and b/HW3/__pycache__/Bar.cpython-312.pyc differ diff --git a/HW3/__pycache__/case1.cpython-312.pyc b/HW3/__pycache__/case1.cpython-312.pyc new file mode 100644 index 0000000..2901a60 Binary files /dev/null and b/HW3/__pycache__/case1.cpython-312.pyc differ diff --git a/HW3/__pycache__/case2.cpython-312.pyc b/HW3/__pycache__/case2.cpython-312.pyc new file mode 100644 index 0000000..830442f Binary files /dev/null and b/HW3/__pycache__/case2.cpython-312.pyc differ diff --git a/HW3/__pycache__/common.cpython-312.pyc b/HW3/__pycache__/common.cpython-312.pyc new file mode 100644 index 0000000..0548c34 Binary files /dev/null and b/HW3/__pycache__/common.cpython-312.pyc differ diff --git a/HW3/case1.py b/HW3/case1.py new file mode 100644 index 0000000..4021a09 --- /dev/null +++ b/HW3/case1.py @@ -0,0 +1,244 @@ +# case1.py + +import numpy as np +import matplotlib.pyplot as plt + +from Bar import Bar, set_k_ratio + +def get_analytical_constants(bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Solve for C1, D1, C2, and D2 + Returns np.array([C1, D1, C2, D2])''' + epsilon = (bar1.k*bar1.area*bar1.alpha) / (bar2.k*bar2.area*bar2.alpha) + A = np.array([ + #C1 D1 C2 D2 + [0, 1, 0, 0], + [0, 0, np.sinh(bar2.alpha*l), np.cosh(bar2.alpha*l)], + [np.sinh(bar1.alpha*x_bar), 0, -1*np.sinh(bar2.alpha*x_bar), -1*np.cosh(bar2.alpha*x_bar)], + [epsilon*np.cosh(bar1.alpha*x_bar), 0, -1*np.cosh(bar2.alpha*x_bar), -1*np.sinh(bar2.alpha*x_bar)] + ]) + + B = np.array([ + 0, + 100, + 0, + 0 + ]) + + return np.linalg.solve(A, B) + +def get_analytical_datapoints(x_points, bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Generates the temperature distribution given x_points. + Returns np.array of tempature values''' + constants = get_analytical_constants(bar1=bar1,bar2=bar2,x_bar=x_bar,l=1) + C1, D1, C2, D2 = constants[0], constants[1], constants[2], constants[3] + + # Two temp functions for the left and right side of the interface + left_temp = lambda x: C1*np.sinh(bar1.alpha*x) + D1*np.cosh(bar1.alpha*x) + right_temp = lambda x: C2*np.sinh(bar2.alpha*x) + D2*np.cosh(bar2.alpha*x) + + temp_values = np.array([]) + for x in x_points: + if x <= x_bar: + temp = left_temp(x) + else: + temp = right_temp(x) + temp_values = np.append(temp_values, temp) + + return temp_values + +def get_analytical_heat_convection(x, bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Gets the analytical heat convection at x from heat flux conservation''' + constants = get_analytical_constants(bar1=bar1,bar2=bar2,x_bar=x_bar,l=1) + C1, D1, C2, D2 = constants[0], constants[1], constants[2], constants[3] + + t_prime = C2*bar2.alpha*np.cosh(bar2.alpha*x) + D2*bar2.alpha*np.sinh(bar2.alpha*x) + + q_dot = -bar2.k*bar2.area*t_prime + return q_dot + +def fdm_array(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1): + '''num_sections = total sections to cut the bar. Note that this is NOT evenly distributed over the entire + bar. Half of the sections are left of x_bar, half are to the right. + Returns the A matrix of the FDM equation Ax = B''' + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + k1 = 2 + bar1.alpha**2 * dx1**2 + k2 = 2 + bar2.alpha**2 * dx2**2 + + beta1 = (bar1.k * bar1.area) / (dx1) + beta2 = (bar2.k * bar2.area) / (dx2) + + #rho = beta1 - beta2 - (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + rho = beta1 + beta2 + (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + # Making the matrix + # There are num_sections - 1 rows and columns + # row 0, row n - 2, and row num_sections/2 - 1 are unique + # the rest are -1, k, -1 where k changes from k1 to k2 at num_sections/2 + + # Initialize matrix A with zeros + A = np.zeros((num_sections - 1, num_sections - 1)) + + # Fill the matrix A + for row in range(num_sections - 1): + if row == 0: + # First row for boundary at x = 0 (for example Dirichlet or Neumann condition) + A[row, row] = k1 + A[row, row+1] = -1 + elif row == num_sections // 2 - 1: + # Special row at the interface between bar1 and bar2 + A[row, row-1] = -beta1 + A[row, row] = rho # This is the rho value you computed + A[row, row+1] = -beta2 + elif row == num_sections - 2: + # Last row for boundary at x = l (again assuming Dirichlet or Neumann condition) + A[row, row-1] = -1 + A[row, row] = k2 + else: + # Interior rows for bar1 and bar2 (uniform grid) + if row < num_sections // 2 - 1: + # In bar1 region + A[row, row-1] = -1 + A[row, row] = k1 + A[row, row+1] = -1 + else: + # In bar2 region + A[row, row-1] = -1 + A[row, row] = k2 + A[row, row+1] = -1 + return A + +def fem_array(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1): + '''num_sections = total sections to cut the bar. Note that this is NOT evenly distributed over the entire + bar. Half of the sections are left of x_bar, half are to the right. + Returns the A matrix of the FEM equation Ax = B''' + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + k1 = (2/dx1 + 4*bar1.alpha**2*dx1/6) / (1/dx1 - bar1.alpha**2*dx1/6) + k2 = (2/dx2 + 4*bar2.alpha**2*dx2/6) / (1/dx2 - bar2.alpha**2*dx2/6) + + beta1 = (bar1.k * bar1.area) / (dx1) + beta2 = (bar2.k * bar2.area) / (dx2) + + rho = beta1 + beta2 + (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + # Making the matrix + # There are num_sections - 1 rows and columns + # row 0, row n - 2, and row num_sections/2 - 1 are unique + # the rest are -1, k, -1 where k changes from k1 to k2 at num_sections/2 + + # Initialize matrix A with zeros + A = np.zeros((num_sections - 1, num_sections - 1)) + + # Fill the matrix A + for row in range(num_sections - 1): + if row == 0: + # First row for boundary at x = 0 (for example Dirichlet or Neumann condition) + A[row, row] = k1 + A[row, row+1] = -1 + elif row == num_sections // 2 - 1: + # Special row at the interface between bar1 and bar2 + A[row, row-1] = -beta1 + A[row, row] = rho # This is the rho value you computed + A[row, row+1] = -beta2 + elif row == num_sections - 2: + # Last row for boundary at x = l (again assuming Dirichlet or Neumann condition) + A[row, row-1] = -1 + A[row, row] = k2 + else: + # Interior rows for bar1 and bar2 (uniform grid) + if row < num_sections // 2 - 1: + # In bar1 region + A[row, row-1] = -1 + A[row, row] = k1 + A[row, row+1] = -1 + else: + # In bar2 region + A[row, row-1] = -1 + A[row, row] = k2 + A[row, row+1] = -1 + return A + +def solve(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1, method="FDM"): + '''Solves using FDM or FEM as specified. Returns (temps, x_values)''' + if method == "FDM": + A = fdm_array(bar1, + bar2, + num_sections=num_sections, + x_bar=x_bar, + l=l) + elif method == "FEM": + A = fem_array(bar1, + bar2, + num_sections=num_sections, + x_bar=x_bar, + l=l) + + B = np.zeros((num_sections - 1)) + B[-1] = 100 + + # Add boundary conditions + interior_temps = np.linalg.solve(A, B) + all_temps = np.array([0]) + all_temps = np.append(all_temps, interior_temps) + all_temps = np.append(all_temps, 100) + + # Generate the x_values to return along with the temps + dx2 = (l - x_bar) / (num_sections / 2) + x_values = np.linspace(0, x_bar, num_sections // 2 + 1) # [0, x_bar] + x_values = np.append(x_values, np.linspace(x_bar+dx2, l, num_sections // 2)) + return all_temps, x_values + + +if __name__ == "__main__": + bar1 = Bar( + radius=0.1, + k=0.5, + h=0.25 + ) + bar2 = Bar( + radius=0.1, + k=2, + h=0.25 + ) + + set_k_ratio(0.5, bar1, bar2) + + x_bar = 0.5 + + x_values = np.linspace(0, 1, 1000) + temp_values = get_analytical_datapoints(x_points=x_values, + bar1=bar1, + bar2=bar2, + x_bar=x_bar) + + + # fdm tests + l = 1 + num_sections = 4 + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + fdm_temps, x_fdm = solve(bar1=bar1, + bar2=bar2, + num_sections=num_sections, + x_bar=x_bar, + l=1, + method="FDM") + fem_temps, x_fem = solve(bar1=bar1, + bar2=bar2, + num_sections=num_sections, + x_bar=x_bar, + l=1, + method="FEM") + + # Plot the data + plt.plot(x_values, temp_values, label='Temperature Distribution') + plt.axvline(x=x_bar, color='r', linestyle='--', label='x_bar') + plt.plot(x_fdm, fdm_temps, 'o-', label='Temperature Distribution (FDM)', color='g') + plt.plot(x_fem, fem_temps, 'o-', label='Temperature Distribution (FEM)', color='r') + plt.xlabel('x') + plt.ylabel('Temperature') + plt.title('Temperature Distribution Along the Bar') + plt.legend() + plt.grid(True) + + # Show plot + plt.show() \ No newline at end of file diff --git a/HW3/case2.py b/HW3/case2.py new file mode 100644 index 0000000..5231180 --- /dev/null +++ b/HW3/case2.py @@ -0,0 +1,230 @@ +# case2.py + +import numpy as np +import matplotlib.pyplot as plt + +from Bar import Bar + +def get_analytical_constants(bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Solve for C1, D1, C2, and D2 + Returns np.array([C1, D1, C2, D2])''' + epsilon = (bar1.k*bar1.area*bar1.alpha) / (bar2.k*bar2.area*bar2.alpha) + A = np.array([ + #C1 D1 C2 D2 + [1, 0, 0, 0], + [0, 0, np.sinh(bar2.alpha*l), np.cosh(bar2.alpha*l)], + [0, np.cosh(bar1.alpha*x_bar), -1*np.sinh(bar2.alpha*x_bar), -1*np.cosh(bar2.alpha*x_bar)], + [0, epsilon*np.sinh(bar1.alpha*x_bar), -1*np.cosh(bar2.alpha*x_bar), -1*np.sinh(bar2.alpha*x_bar)] + ]) + + B = np.array([ + 0, + 100, + 0, + 0 + ]) + + return np.linalg.solve(A, B) + +def get_analytical_datapoints(x_points, bar1:'Bar', bar2:'Bar', x_bar=2/np.pi, l=1): + '''Generates the temperature distribution given x_points. + Returns np.array of tempature values''' + constants = get_analytical_constants(bar1=bar1,bar2=bar2,x_bar=x_bar,l=1) + C1, D1, C2, D2 = constants[0], constants[1], constants[2], constants[3] + + # Two temp functions for the left and right side of the interface + left_temp = lambda x: C1*np.sinh(bar1.alpha*x) + D1*np.cosh(bar1.alpha*x) + right_temp = lambda x: C2*np.sinh(bar2.alpha*x) + D2*np.cosh(bar2.alpha*x) + + temp_values = np.array([]) + for x in x_points: + if x <= x_bar: + temp = left_temp(x) + else: + temp = right_temp(x) + temp_values = np.append(temp_values, temp) + + return temp_values + +def fdm_array(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1): + '''num_sections = total sections to cut the bar. Note that this is NOT evenly distributed over the entire + bar. Half of the sections are left of x_bar, half are to the right. + Returns the A matrix of the FDM equation Ax = B''' + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + k1 = 2 + bar1.alpha**2 * dx1**2 + k2 = 2 + bar2.alpha**2 * dx2**2 + k1_prime = k1/2 + + beta1 = (bar1.k * bar1.area) / (dx1) + beta2 = (bar2.k * bar2.area) / (dx2) + + #rho = beta1 - beta2 - (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + rho = beta1 + beta2 + (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + # Making the matrix + # There are num_sections rows and columns + # row 0, row n - 2, and row num_sections/2 are unique + # the rest are -1, k, -1 where k changes from k1 to k2 at num_sections/2 + + # Initialize matrix A with zeros + A = np.zeros((num_sections, num_sections)) + + # Fill the matrix A + for row in range(num_sections): + if row == 0: + # First row for boundary at x = 0 ( Neumann condition) + A[row, row] = k1_prime + A[row, row+1] = -1 + elif row == num_sections // 2: + # Special row at the interface between bar1 and bar2 + A[row, row-1] = -beta1 + A[row, row] = rho # This is the rho value you computed + A[row, row+1] = -beta2 + elif row == num_sections - 1: + # Last row for boundary at x = l (again assuming Dirichlet) + A[row, row-1] = -1 + A[row, row] = k2 + else: + # Interior rows for bar1 and bar2 (uniform grid) + if row < num_sections // 2 - 1: + # In bar1 region + A[row, row-1] = -1 + A[row, row] = k1 + A[row, row+1] = -1 + else: + # In bar2 region + A[row, row-1] = -1 + A[row, row] = k2 + A[row, row+1] = -1 + return A + +def fem_array(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1): + '''num_sections = total sections to cut the bar. Note that this is NOT evenly distributed over the entire + bar. Half of the sections are left of x_bar, half are to the right. + Returns the A matrix of the FDM equation Ax = B''' + dx1 = x_bar / (num_sections / 2) + dx2 = (l - x_bar) / (num_sections / 2) + k1 = (2/dx1 + 4*bar1.alpha**2*dx1/6) / (1/dx1 - bar1.alpha**2*dx1/6) + k2 = (2/dx2 + 4*bar2.alpha**2*dx2/6) / (1/dx2 - bar2.alpha**2*dx2/6) + k1_prime = k1/2 + + beta1 = (bar1.k * bar1.area) / (dx1) + beta2 = (bar2.k * bar2.area) / (dx2) + + #rho = beta1 - beta2 - (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + rho = beta1 + beta2 + (bar1.k * bar1.area * dx1 * bar1.alpha**2) / 2 + (bar2.k * bar2.area * dx2 * bar2.alpha**2) / 2 + # Making the matrix + # There are num_sections rows and columns + # row 0, row n - 2, and row num_sections/2 are unique + # the rest are -1, k, -1 where k changes from k1 to k2 at num_sections/2 + + # Initialize matrix A with zeros + A = np.zeros((num_sections, num_sections)) + + # Fill the matrix A + for row in range(num_sections): + if row == 0: + # First row for boundary at x = 0 ( Neumann condition) + A[row, row] = k1_prime + A[row, row+1] = -1 + elif row == num_sections // 2: + # Special row at the interface between bar1 and bar2 + A[row, row-1] = -beta1 + A[row, row] = rho # This is the rho value you computed + A[row, row+1] = -beta2 + elif row == num_sections - 1: + # Last row for boundary at x = l (again assuming Dirichlet) + A[row, row-1] = -1 + A[row, row] = k2 + else: + # Interior rows for bar1 and bar2 (uniform grid) + if row < num_sections // 2 - 1: + # In bar1 region + A[row, row-1] = -1 + A[row, row] = k1 + A[row, row+1] = -1 + else: + # In bar2 region + A[row, row-1] = -1 + A[row, row] = k2 + A[row, row+1] = -1 + return A + +def solve(bar1:'Bar', bar2:'Bar', num_sections: int=6, x_bar=2/np.pi, l=1, method="FDM"): + if method == "FDM": + A = fdm_array(bar1, + bar2, + num_sections=num_sections, + x_bar=x_bar, + l=l) + elif method == "FEM": + A = fem_array(bar1, + bar2, + num_sections=num_sections, + x_bar=x_bar, + l=l) + + B = np.zeros((num_sections)) + B[-1] = 100 + + # Append the boundary condition + interior_temps = np.linalg.solve(A, B) + all_temps = np.append(interior_temps, 100) + + # Generate the x_values to return along with the temps + dx2 = (l - x_bar) / (num_sections / 2) + x_values = np.linspace(0, x_bar, num_sections // 2 + 1) # [0, x_bar] + x_values = np.append(x_values, np.linspace(x_bar+dx2, l, num_sections // 2)) + return all_temps, x_values + + +if __name__ == "__main__": + bar1 = Bar( + radius=0.1, + k=0.5, + h=0.25 + ) + bar2 = Bar( + radius=0.1, + k=2, + h=0.25 + ) + + x_bar = 2/np.pi + + + x_values = np.linspace(0, 1, 1000) + temp_values = get_analytical_datapoints(x_points=x_values, + bar1=bar1, + bar2=bar2) + + + # fdm tests + l = 1 + num_sections = 6 + fdm_temps, x_fdm = solve(bar1=bar1, + bar2=bar2, + num_sections=num_sections, + x_bar=x_bar, + l=1, + method="FDM") + fem_temps, x_fem = solve(bar1=bar1, + bar2=bar2, + num_sections=num_sections, + x_bar=x_bar, + l=1, + method="FEM") + + # Plot the data + plt.plot(x_values, temp_values, label='Temperature Distribution') + plt.axvline(x=x_bar, color='r', linestyle='--', label='x_bar') + plt.plot(x_fdm, fdm_temps, 'o-', label='Temperature Distribution (FDM)', color='g') + plt.plot(x_fem, fem_temps, 'o-', label='Temperature Distribution (FEM)', color='r') + plt.xlabel('x') + plt.ylabel('Temperature') + plt.title('Temperature Distribution Along the Bar') + plt.legend() + plt.grid(True) + + # Show plot + plt.show() \ No newline at end of file diff --git a/HW3/common.py b/HW3/common.py new file mode 100644 index 0000000..bb54345 --- /dev/null +++ b/HW3/common.py @@ -0,0 +1,30 @@ +# common.py + +import numpy as np +from Bar import Bar + +def fdm_heat_extraction(t_0, t_1, dx, bar:'Bar', order=2): + '''Get the heat conduction at the point t_1 using Taylor series''' + if order == 1: + return -1 * bar.k * bar.area * (t_1 - t_0) / dx + elif order == 2: + return -1 * bar.k * bar.area * (((t_1 - t_0) / dx) + (bar.alpha**2 * dx * t_1 / 2)) + +def fem_heat_extraction(t_0, t_1, dx, bar:'Bar'): + '''Get the heat conduction at the point t_1 using FEM equation''' + term_1 = (-1/dx + bar.alpha**2*dx/6) * t_0 + term_2 = (1/dx + 2*bar.alpha**2*dx/6) * t_1 + return -1 * bar.k * bar.area * (term_1 + term_2) + +def calc_error(exact, q_1): + return np.abs((exact - q_1) / exact) + +def calc_beta(exact, q_1, q_2, dx_1, dx_2): + return np.log(np.abs((exact - q_1)/(exact - q_2))) / np.log(dx_1 / dx_2) + +def calc_extrapolated(q1, q2, q3): + '''Calcs the Richardson Extrapolation value based on 3 different meshes. + Assumes that q3 is 2x finer than q2 is 2x finer than q1''' + numerator = q1*q3 - q2**2 + denominator = q1 + q3 - 2*q2 + return numerator / denominator \ No newline at end of file diff --git a/HW3/constants.py b/HW3/constants.py new file mode 100644 index 0000000..5248508 --- /dev/null +++ b/HW3/constants.py @@ -0,0 +1,19 @@ +from numpy import pi + +h = 0.0025 # W / cm^2*K + +class Bar1(): + def __init__(self): + self.k = 0.5 # W / cm*K + self.r = 0.1 # cm + self.A = pi*self.r**2 # cm^2 + self.P = 2*pi*self.r # cm + +class Bar2(): + def __init__(self): + self.k = 0.5 # W / cm*K + self.r = 0.1 # cm + self.A = pi*self.r**2 # cm^2 + self.P = 2*pi*self.r # cm + + self.alpha = ((h*self.P)/(self.k*self.A))**0.5 \ No newline at end of file diff --git a/HW3/images/section1/x_bar_1/analytical_heat_convection_vs_k_ratio.png b/HW3/images/section1/x_bar_1/analytical_heat_convection_vs_k_ratio.png new file mode 100644 index 0000000..0f7c4b3 Binary files /dev/null and b/HW3/images/section1/x_bar_1/analytical_heat_convection_vs_k_ratio.png differ diff --git a/HW3/images/section1/x_bar_1/analytical_temperature_vs_position.png b/HW3/images/section1/x_bar_1/analytical_temperature_vs_position.png new file mode 100644 index 0000000..bef0e63 Binary files /dev/null and b/HW3/images/section1/x_bar_1/analytical_temperature_vs_position.png differ diff --git a/HW3/images/section1/x_bar_2/analytical_heat_convection_vs_k_ratio.png b/HW3/images/section1/x_bar_2/analytical_heat_convection_vs_k_ratio.png new file mode 100644 index 0000000..52186d4 Binary files /dev/null and b/HW3/images/section1/x_bar_2/analytical_heat_convection_vs_k_ratio.png differ diff --git a/HW3/images/section1/x_bar_2/analytical_temperature_vs_position.png b/HW3/images/section1/x_bar_2/analytical_temperature_vs_position.png new file mode 100644 index 0000000..3e47b99 Binary files /dev/null and b/HW3/images/section1/x_bar_2/analytical_temperature_vs_position.png differ diff --git a/HW3/images/section2/x_bar_1/fdm_temperature_vs_position.png b/HW3/images/section2/x_bar_1/fdm_temperature_vs_position.png new file mode 100644 index 0000000..af5c966 Binary files /dev/null and b/HW3/images/section2/x_bar_1/fdm_temperature_vs_position.png differ diff --git a/HW3/images/section2/x_bar_1/fem_temperature_vs_position.png b/HW3/images/section2/x_bar_1/fem_temperature_vs_position.png new file mode 100644 index 0000000..9411db2 Binary files /dev/null and b/HW3/images/section2/x_bar_1/fem_temperature_vs_position.png differ diff --git a/HW3/images/section2/x_bar_1/heat_convection_vs_k_ratio.png b/HW3/images/section2/x_bar_1/heat_convection_vs_k_ratio.png new file mode 100644 index 0000000..77c2a96 Binary files /dev/null and b/HW3/images/section2/x_bar_1/heat_convection_vs_k_ratio.png differ diff --git a/HW3/images/section2/x_bar_2/fdm_temperature_vs_position.png b/HW3/images/section2/x_bar_2/fdm_temperature_vs_position.png new file mode 100644 index 0000000..011d8ff Binary files /dev/null and b/HW3/images/section2/x_bar_2/fdm_temperature_vs_position.png differ diff --git a/HW3/images/section2/x_bar_2/fem_temperature_vs_position.png b/HW3/images/section2/x_bar_2/fem_temperature_vs_position.png new file mode 100644 index 0000000..56630cc Binary files /dev/null and b/HW3/images/section2/x_bar_2/fem_temperature_vs_position.png differ diff --git a/HW3/images/section2/x_bar_2/heat_convection_vs_k_ratio.png b/HW3/images/section2/x_bar_2/heat_convection_vs_k_ratio.png new file mode 100644 index 0000000..d040ec3 Binary files /dev/null and b/HW3/images/section2/x_bar_2/heat_convection_vs_k_ratio.png differ diff --git a/HW3/images/section3/x_bar_1/heat_convergences.png b/HW3/images/section3/x_bar_1/heat_convergences.png new file mode 100644 index 0000000..fa8bcdd Binary files /dev/null and b/HW3/images/section3/x_bar_1/heat_convergences.png differ diff --git a/HW3/images/section3/x_bar_1/temp_convergences.png b/HW3/images/section3/x_bar_1/temp_convergences.png new file mode 100644 index 0000000..1e47f92 Binary files /dev/null and b/HW3/images/section3/x_bar_1/temp_convergences.png differ diff --git a/HW3/images/section3/x_bar_2/heat_convergences.png b/HW3/images/section3/x_bar_2/heat_convergences.png new file mode 100644 index 0000000..78613b6 Binary files /dev/null and b/HW3/images/section3/x_bar_2/heat_convergences.png differ diff --git a/HW3/images/section3/x_bar_2/temp_convergences.png b/HW3/images/section3/x_bar_2/temp_convergences.png new file mode 100644 index 0000000..3aa08aa Binary files /dev/null and b/HW3/images/section3/x_bar_2/temp_convergences.png differ diff --git a/HW3/main.py b/HW3/main.py new file mode 100644 index 0000000..e36efe4 --- /dev/null +++ b/HW3/main.py @@ -0,0 +1,722 @@ +# main.py + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import math + +import case1 +import case2 +import common +import Bar + + +bar1 = Bar.Bar( + radius=0.1, + k=0.5, + h=0.25 +) +bar2 = Bar.Bar( + radius=0.1, + k=2, + h=0.25 +) + +data_dict = {} +data_dict["k_ratios"] = np.array([1/16, 1/8, 1/4, 1/2, 1, 2, 4, 8, 16]) +data_dict["x_bar_values"] = np.array([1/2, 2/np.pi]) +data_dict["length"] = 1 + +def section_1(): + '''Analytical solutions with varying k_ratio''' + print("Section 1: Analytical") + + l = data_dict["length"] + + # Calculate the results for every k ratio for x_bar_1 + x_bar_1 = data_dict["x_bar_values"][0] + + data_dict["section1"] = {} + data_dict["section1"]["x_bar_1"] = {} + data_dict["section1"]["x_bar_1"]["x_values_fine"] = np.linspace(0, l, 50) + data_dict["section1"]["x_bar_1"]["x_values_fine"] = np.sort(np.append(data_dict["section1"]["x_bar_1"]["x_values_fine"], x_bar_1)) + data_dict["section1"]["x_bar_1"]["x_values_course"] = np.linspace(0, l, 9) + data_dict["section1"]["x_bar_1"]["x_values_course"] = np.sort(np.append(data_dict["section1"]["x_bar_1"]["x_values_course"], x_bar_1)) + + results = [] + results_course = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + result = case1.get_analytical_datapoints(data_dict["section1"]["x_bar_1"]["x_values_fine"], + bar1=bar1, + bar2=bar2, + x_bar=x_bar_1, + l=l) + result_course = case1.get_analytical_datapoints(data_dict["section1"]["x_bar_1"]["x_values_course"], + bar1=bar1, + bar2=bar2, + x_bar=x_bar_1, + l=l) + results.append(result) + results_course.append(result_course) + + data_dict["section1"]["x_bar_1"]["temp_results"] = np.array(results) + + df_x_bar_1 = pd.DataFrame(results_course, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in data_dict["section1"]["x_bar_1"]["x_values_course"]]) + + print("Analytical Temperature Results x_bar = 0.5:") + print(df_x_bar_1.to_string()) + print("\n" * 2) + + # Plotting x_bar_1 + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section1"]["x_bar_1"]["x_values_fine"], data_dict["section1"]["x_bar_1"]["temp_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_1, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section1/x_bar_1/analytical_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + # Calculate the temperature results for every k ratio for x_bar_1 + x_bar_2 = data_dict["x_bar_values"][1] + + data_dict["section1"]["x_bar_2"] = {} + data_dict["section1"]["x_bar_2"]["x_values_fine"] = np.linspace(0, l, 50) + data_dict["section1"]["x_bar_2"]["x_values_fine"] = np.sort(np.append(data_dict["section1"]["x_bar_2"]["x_values_fine"], x_bar_2)) + data_dict["section1"]["x_bar_2"]["x_values_course"] = np.linspace(0, l, 9) + data_dict["section1"]["x_bar_2"]["x_values_course"] = np.sort(np.append(data_dict["section1"]["x_bar_2"]["x_values_course"], x_bar_2)) + + results = [] + results_course = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) + result = case1.get_analytical_datapoints(data_dict["section1"]["x_bar_2"]["x_values_fine"], + bar1=bar1, + bar2=bar2, + x_bar=x_bar_2, + l=l) + result_course = case1.get_analytical_datapoints(data_dict["section1"]["x_bar_2"]["x_values_course"], + bar1=bar1, + bar2=bar2, + x_bar=x_bar_2, + l=l) + results.append(result) + results_course.append(result_course) + + data_dict["section1"]["x_bar_2"]["temp_results"] = np.array(results) + + df_x_bar_2 = pd.DataFrame(results_course, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in data_dict["section1"]["x_bar_2"]["x_values_course"]]) + + print("Analytical Temperature Results x_bar = 2/pi:") + print(df_x_bar_2.to_string()) + print("\n" * 2) + + # Plotting x_bar_2 + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section1"]["x_bar_2"]["x_values_fine"], data_dict["section1"]["x_bar_2"]["temp_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_2, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section1/x_bar_2/analytical_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + + + # Calculate the analytical heat convection leaving the rod at x=l for varying k ratio + + results_1 = [] + results_2 = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) + result_1 = case1.get_analytical_heat_convection(x=l, + bar1=bar1, + bar2=bar2, + x_bar=x_bar_1, + l=l) + result_2 = case1.get_analytical_heat_convection(x=l, + bar1=bar1, + bar2=bar2, + x_bar=x_bar_2, + l=l) + results_1.append([result_1]) + results_2.append([result_2]) + + data_dict["section1"]["x_bar_1"]["heat_results"] = np.array(results_1) + data_dict["section1"]["x_bar_2"]["heat_results"] = np.array(results_2) + + df_x_bar_1 = pd.DataFrame(results_1, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + df_x_bar_2 = pd.DataFrame(results_2, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + + print("Analytical Heat Convection Results x_bar = 0.5:") + print(df_x_bar_1.to_string()) + print("\n" * 2) + + print("Analytical Heat Convection Results x_bar = 2/pi:") + print(df_x_bar_2.to_string()) + print("\n" * 2) + + # Plotting x_bar_1 heat convection + plt.figure(figsize=(10, 6)) + plt.plot(data_dict["k_ratios"], data_dict["section1"]["x_bar_1"]["heat_results"], label=f'Heat Convection') + plt.xlabel('k2/k1') + plt.ylabel('Q_Dot (W)') + plt.legend() + plt.grid(True) + plt.savefig('images/section1/x_bar_1/analytical_heat_convection_vs_k_ratio.png', dpi=300) + #plt.show() + plt.clf() + + # Plotting x_bar_2 heat convection + plt.figure(figsize=(10, 6)) + plt.plot(data_dict["k_ratios"], data_dict["section1"]["x_bar_2"]["heat_results"], label=f'Heat Convection') + plt.xlabel('k2/k1') + plt.ylabel('Q_Dot (W)') + plt.legend() + plt.grid(True) + plt.savefig('images/section1/x_bar_2/analytical_heat_convection_vs_k_ratio.png', dpi=300) + #plt.show() + plt.clf() + + +def section_2(): + '''FDM and FEM with varying k_ratio''' + print("Section 2: FDM and FEM with Varying k_ratio") + + l = data_dict["length"] + + data_dict["section2"] = {} + data_dict["section2"]["num_sections"] = 8 + + + # Calculate the results for every k ratio for x_bar_1 + x_bar_1 = data_dict["x_bar_values"][0] + + # For every k_ratio, calculate the FDM and FEM temperature results + fdm_results = [] + fem_results = [] + x_fdm = [] + x_fem = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + fdm_temps, x_fdm = case1.solve(bar1=bar1, + bar2=bar2, + num_sections=data_dict["section2"]["num_sections"], + x_bar=x_bar_1, + l=l, + method="FDM") + fem_temps, x_fem = case1.solve(bar1=bar1, + bar2=bar2, + num_sections=data_dict["section2"]["num_sections"], + x_bar=x_bar_1, + l=l, + method="FEM") + + fdm_results.append(fdm_temps) + fem_results.append(fem_temps) + + data_dict["section2"]["x_bar_1"] = {} + data_dict["section2"]["x_bar_1"]["x_values"] = x_fdm # fdm and fem use same x_values + data_dict["section2"]["x_bar_1"]["fdm_results"] = np.array(fdm_results) + data_dict["section2"]["x_bar_1"]["fem_results"] = np.array(fem_results) + + df_fdm = pd.DataFrame(fdm_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in x_fdm]) + df_fem = pd.DataFrame(fem_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in x_fem]) + + print("FDM Temperature Results x_bar = 0.5:") + print(df_fdm.to_string()) + print("\n" * 2) + print("FEM Temperature Results x_bar = 0.5:") + print(df_fem.to_string()) + print("\n" * 2) + + # Now that the data is gathered for FDM and FEM, plot it + # Plotting x_bar_1, FDM + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section2"]["x_bar_1"]["x_values"], data_dict["section2"]["x_bar_1"]["fdm_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_1, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_1/fdm_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + # Plotting x_bar_1, FEM + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section2"]["x_bar_1"]["x_values"], data_dict["section2"]["x_bar_1"]["fem_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_1, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_1/fem_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + + # Calculate the results for every k ratio for x_bar_2 + x_bar_2 = data_dict["x_bar_values"][1] + + # For every k_ratio, calculate the FDM and FEM temperature results + fdm_results = [] + fem_results = [] + x_fdm = [] + x_fem = [] + for k_ratio in data_dict["k_ratios"]: + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + fdm_temps, x_fdm = case1.solve(bar1=bar1, + bar2=bar2, + num_sections=data_dict["section2"]["num_sections"], + x_bar=x_bar_2, + l=l, + method="FDM") + fem_temps, x_fem = case1.solve(bar1=bar1, + bar2=bar2, + num_sections=data_dict["section2"]["num_sections"], + x_bar=x_bar_2, + l=l, + method="FEM") + + fdm_results.append(fdm_temps) + fem_results.append(fem_temps) + + data_dict["section2"]["x_bar_2"] = {} + data_dict["section2"]["x_bar_2"]["x_values"] = x_fdm # fdm and fem use same x_values + data_dict["section2"]["x_bar_2"]["fdm_results"] = np.array(fdm_results) + data_dict["section2"]["x_bar_2"]["fem_results"] = np.array(fem_results) + + df_fdm = pd.DataFrame(fdm_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in x_fdm]) + df_fem = pd.DataFrame(fem_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=[f'x = {x:.3f}' for x in x_fem]) + + print("FDM Temperature Results x_bar = 2/pi:") + print(df_fdm.to_string()) + print("\n" * 2) + print("FEM Temperature Results x_bar = 2/pi:") + print(df_fem.to_string()) + print("\n" * 2) + + # Plotting x_bar_2, FDM + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section2"]["x_bar_2"]["x_values"], data_dict["section2"]["x_bar_2"]["fdm_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_2, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_2/fdm_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + # Plotting x_bar_2, FEM + plt.figure(figsize=(10, 6)) + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + plt.plot(data_dict["section2"]["x_bar_2"]["x_values"], data_dict["section2"]["x_bar_2"]["fem_results"][idx], label=f'k2/k1 = {k_ratio}') + plt.axvline(x=x_bar_2, color='r', linestyle='--', label='x_bar') + plt.xlabel('Position (cm)') + plt.ylabel('Temperature (°C)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_2/fem_temperature_vs_position.png', dpi=300) + #plt.show() + plt.clf() + + + + # After calculating temperature values, extract the heat convection at x=l + + # x_bar_1 + # for every k ratio, calculate the heat convection from fdm and fem temp results + fdm_heat_results = [] + fem_heat_results = [] + for i, k_ratio in enumerate(data_dict["k_ratios"]): + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + dx = data_dict["section2"]["x_bar_1"]["x_values"][-1] - data_dict["section2"]["x_bar_1"]["x_values"][-2] + (fdm_t0, fdm_t1) = data_dict["section2"]["x_bar_1"]["fdm_results"][i][-2:] + fdm_heat_result = common.fdm_heat_extraction(fdm_t0, fdm_t1, dx, bar2, order=2) + fdm_heat_results.append(fdm_heat_result) + + (fem_t0, fem_t1) = data_dict["section2"]["x_bar_1"]["fem_results"][i][-2:] + fem_heat_result = common.fem_heat_extraction(fem_t0, fem_t1, dx, bar2) + fem_heat_results.append(fem_heat_result) + + + data_dict["section2"]["x_bar_1"]["fdm_heat_results"] = np.array(fdm_heat_results) + data_dict["section2"]["x_bar_1"]["fem_heat_results"] = np.array(fem_heat_results) + + df_fdm = pd.DataFrame(fdm_heat_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + df_fem = pd.DataFrame(fem_heat_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + + print("FDM Heat Convection Results x_bar = 0.5:") + print(df_fdm.to_string()) + print("\n" * 2) + + print("FEM Heat Convection Results x_bar = 0.5:") + print(df_fem.to_string()) + print("\n" * 2) + + # x_bar_2 + # for every k ratio, calculate the heat convection from fdm and fem temp results + fdm_heat_results = [] + fem_heat_results = [] + for i, k_ratio in enumerate(data_dict["k_ratios"]): + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + dx = data_dict["section2"]["x_bar_2"]["x_values"][-1] - data_dict["section2"]["x_bar_2"]["x_values"][-2] + (fdm_t0, fdm_t1) = data_dict["section2"]["x_bar_2"]["fdm_results"][i][-2:] + fdm_heat_result = common.fdm_heat_extraction(fdm_t0, fdm_t1, dx, bar2, order=2) + fdm_heat_results.append(fdm_heat_result) + + (fem_t0, fem_t1) = data_dict["section2"]["x_bar_2"]["fem_results"][i][-2:] + fem_heat_result = common.fem_heat_extraction(fem_t0, fem_t1, dx, bar2) + fem_heat_results.append(fem_heat_result) + + + data_dict["section2"]["x_bar_2"]["fdm_heat_results"] = np.array(fdm_heat_results) + data_dict["section2"]["x_bar_2"]["fem_heat_results"] = np.array(fem_heat_results) + + df_fdm = pd.DataFrame(fdm_heat_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + df_fem = pd.DataFrame(fem_heat_results, index=[f'k2/k1 = {a}' for a in data_dict["k_ratios"]], columns=["Heat Convection"]) + + print("FDM Heat Convection Results x_bar = 2/pi:") + print(df_fdm.to_string()) + print("\n" * 2) + + print("FEM Heat Convection Results x_bar = 2/pi:") + print(df_fem.to_string()) + print("\n" * 2) + + # Plotting x_bar_1 heat convection + plt.figure(figsize=(10, 6)) + plt.plot(data_dict["k_ratios"], data_dict["section2"]["x_bar_1"]["fdm_heat_results"], label=f'FDM Heat Convection') + plt.plot(data_dict["k_ratios"], data_dict["section2"]["x_bar_1"]["fem_heat_results"], label=f'FEM Heat Convection') + #plt.plot(data_dict["k_ratios"], data_dict["section1"]["x_bar_1"]["heat_results"], label=f'Analytical Heat Convection') + plt.xlabel('k2/k1') + plt.ylabel('Q_Dot (W)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_1/heat_convection_vs_k_ratio.png', dpi=300) + #plt.show() + plt.clf() + + # Plotting x_bar_2 heat convection + plt.figure(figsize=(10, 6)) + plt.plot(data_dict["k_ratios"], data_dict["section2"]["x_bar_2"]["fdm_heat_results"], label=f'FDM Heat Convection') + plt.plot(data_dict["k_ratios"], data_dict["section2"]["x_bar_2"]["fem_heat_results"], label=f'FEM Heat Convection') + #plt.plot(data_dict["k_ratios"], data_dict["section1"]["x_bar_2"]["heat_results"], label=f'Analytical Heat Convection') + plt.xlabel('k2/k1') + plt.ylabel('Q_Dot (W)') + plt.legend() + plt.grid(True) + plt.savefig('images/section2/x_bar_2/heat_convection_vs_k_ratio.png', dpi=300) + #plt.show() + plt.clf() + + + +def section_3(): + '''Convergence of FDM and FEM of interface temperature and heat convection''' + print("Section 3: Convergence of FDM and FEM with Varying k Ratio") + + l = data_dict["length"] + + data_dict["section3"] = {} + data_dict["section3"]["num_sections"] = [4, 8, 16, 32, 64, 128] + + for k, x_bar_1 in enumerate(data_dict["x_bar_values"]): + # Calculate the number of rows needed for two columns of subplots + num_k_ratios = len(data_dict["k_ratios"]) + num_rows = math.ceil(num_k_ratios / 2) + + # Define the figure for temperature convergence (two columns) + fig_temp, axs_temp = plt.subplots(num_rows, 2, figsize=(15, 5 * num_rows)) + + # Define the figure for heat convergence (two columns) + fig_heat, axs_heat = plt.subplots(num_rows, 2, figsize=(15, 5 * num_rows)) + + # Flatten the axs arrays for easier indexing + axs_temp = axs_temp.flatten() + axs_heat = axs_heat.flatten() + + # For every k ratio, iterate through every number of points, calculate convergence (compared to analytical and extrapolated), print and add a subplot + for idx, k_ratio in enumerate(data_dict["k_ratios"]): + Bar.set_k_ratio(k_ratio, bar1, bar2) # Set the k value of bar2 = k*bar1 + + # Get the analytical data + analy_interface_temp = case1.get_analytical_datapoints([x_bar_1], bar1, bar2, x_bar_1, l) + analy_heat = case1.get_analytical_heat_convection(l, bar1, bar2, x_bar_1, l) + + # Setup data arrays + fdm_interface_temps = [] + fem_interface_temps = [] + fdm_extrap_temps = [None, None] + fem_extrap_temps = [None, None] + fdm_heats = [] + fem_heats = [] + fdm_extrap_heats = [None, None] + fem_extrap_heats = [None, None] + fdm_interface_temp_errors = [] + fem_interface_temp_errors = [] + fdm_interface_temp_betas = [[None]] + fem_interface_temp_betas = [[None]] + fdm_heat_errors = [] + fem_heat_errors = [] + fdm_heat_betas = [[None]] + fem_heat_betas = [[None]] + fdm_interface_temp_extrap_errors = [None, None] + fem_interface_temp_extrap_errors = [None, None] + fdm_interface_temp_extrap_betas = [None, None] + fem_interface_temp_extrap_betas = [None, None] + fdm_heat_extrap_errors = [None, None] + fem_heat_extrap_errors = [None, None] + fdm_heat_extrap_betas = [None, None] + fem_heat_extrap_betas = [None, None] + + # Get the fdm and fem data + dx = 0 + for i, num_sections in enumerate(data_dict["section3"]["num_sections"]): + fdm_temps, x_fdm = case1.solve(bar1=bar1,bar2=bar2,num_sections=num_sections,x_bar=x_bar_1,l=l,method="FDM") + fem_temps, x_fem = case1.solve(bar1=bar1,bar2=bar2,num_sections=num_sections,x_bar=x_bar_1,l=l,method="FEM") + # Get the interface temperature + fdm_interface_temp = fdm_temps[num_sections // 2] + fem_interface_temp = fem_temps[num_sections // 2] + fdm_interface_temps.append(fdm_interface_temp) + fem_interface_temps.append(fem_interface_temp) + + # Get the 2 last temperature to calculate heat convection + dx_prev = dx + dx = x_fdm[-1] - x_fdm[-2] + (fdm_t0, fdm_t1) = fdm_temps[-2:] + fdm_heat = common.fdm_heat_extraction(fdm_t0, fdm_t1, dx, bar2, order=2) + (fem_t0, fem_t1) = fem_temps[-2:] + fem_heat = common.fem_heat_extraction(fem_t0, fem_t1, dx, bar2) + fdm_heats.append(fdm_heat) + fem_heats.append(fem_heat) + + # Calculated the Richardson extrpolation of interface temperature and heat convection + if i >= 2: + fdm_extrap_temp = common.calc_extrapolated(fdm_interface_temps[-3], fdm_interface_temps[-2], fdm_interface_temps[-1]) + fem_extrap_temp = common.calc_extrapolated(fem_interface_temps[-3], fem_interface_temps[-2], fem_interface_temps[-1]) + fdm_extrap_heat = common.calc_extrapolated(fdm_heats[-3], fdm_heats[-2], fdm_heats[-1]) + fem_extrap_heat = common.calc_extrapolated(fem_heats[-3], fem_heats[-2], fem_heats[-1]) + fdm_extrap_temps.append(fdm_extrap_temp) + fem_extrap_temps.append(fem_extrap_temp) + fdm_extrap_heats.append(fdm_extrap_heat) + fem_extrap_heats.append(fem_extrap_heat) + + # Calculate the errors in reference to extrapolated values + fdm_interface_temp_extrap_error = common.calc_error(fdm_extrap_temp, fdm_interface_temp) + fem_interface_temp_extrap_error = common.calc_error(fem_extrap_temp, fem_interface_temp) + fdm_heat_extrap_error = common.calc_error(fdm_extrap_heat, fdm_heat) + fem_heat_extrap_error = common.calc_error(fem_extrap_heat, fem_heat) + fdm_interface_temp_extrap_errors.append(fdm_interface_temp_extrap_error) + fem_interface_temp_extrap_errors.append(fem_interface_temp_extrap_error) + fdm_heat_extrap_errors.append(fdm_heat_extrap_error) + fem_heat_extrap_errors.append(fem_heat_extrap_error) + + # Calculate the betas in reference to extrapolated values + fdm_interface_temp_extrap_beta = common.calc_beta(fdm_extrap_temp, fdm_interface_temps[-1], fdm_interface_temps[-2], dx, dx_prev) + fem_interface_temp_extrap_beta = common.calc_beta(fem_extrap_temp, fem_interface_temps[-1], fem_interface_temps[-2], dx, dx_prev) + fdm_heat_extrap_beta = common.calc_beta(fdm_extrap_heat, fdm_heats[-1], fdm_heats[-2], dx, dx_prev) + fem_heat_extrap_beta = common.calc_beta(fem_extrap_heat, fem_heats[-1], fem_heats[-2], dx, dx_prev) + fdm_interface_temp_extrap_betas.append(fdm_interface_temp_extrap_beta) + fem_interface_temp_extrap_betas.append(fem_interface_temp_extrap_beta) + fdm_heat_extrap_betas.append(fdm_heat_extrap_beta) + fem_heat_extrap_betas.append(fem_heat_extrap_beta) + + + # Calculate the error and convergence rates for fdm temp, fem temp, fdm heat, and fem heat + fdm_interface_temp_error = common.calc_error(analy_interface_temp, fdm_interface_temp) + fem_interface_temp_error = common.calc_error(analy_interface_temp, fem_interface_temp) + fdm_heat_error = common.calc_error(analy_heat, fdm_heat) + fem_heat_error = common.calc_error(analy_heat, fem_heat) + fdm_interface_temp_errors.append(fdm_interface_temp_error) + fem_interface_temp_errors.append(fem_interface_temp_error) + fdm_heat_errors.append(fdm_heat_error) + fem_heat_errors.append(fem_heat_error) + + if i >= 1: # if i == 0 then we cannot calculate convergence + fdm_interface_temp_beta = common.calc_beta(analy_interface_temp, fdm_interface_temps[-1], fdm_interface_temps[-2], dx, dx_prev) + fem_interface_temp_beta = common.calc_beta(analy_interface_temp, fem_interface_temps[-1], fem_interface_temps[-2], dx, dx_prev) + fdm_heat_beta = common.calc_beta(analy_heat, fdm_heats[-1], fdm_heats[-2], dx, dx_prev) + fem_heat_beta = common.calc_beta(analy_heat, fem_heats[-1], fem_heats[-2], dx, dx_prev) + fdm_interface_temp_betas.append(fdm_interface_temp_beta) + fem_interface_temp_betas.append(fem_interface_temp_beta) + fdm_heat_betas.append(fdm_heat_beta) + fem_heat_betas.append(fem_heat_beta) + + + # Print the interface temps for this k ratio + table_data = [] + for i in range(len(data_dict["section3"]["num_sections"])): + table_data.append([ + analy_interface_temp[0], # True T + + fdm_interface_temps[i], # fdm result + fdm_extrap_temps[i], + fdm_interface_temp_errors[i][0], # fdm % error in reference to analytical + fdm_interface_temp_betas[i][0], # fdm beta + fdm_interface_temp_extrap_errors[i], # fdm % error in reference to extrapolated value + fdm_interface_temp_extrap_betas[i], # fdm beta + + fem_interface_temps[i], # fem result + fem_extrap_temps[i], + fem_interface_temp_errors[i][0], # fem % error + fem_interface_temp_betas[i][0], # fem beta + fem_interface_temp_extrap_errors[i], # fem extrap % error + fem_interface_temp_extrap_betas[i], # fem beta + ]) + + columns = ['True T', 'FDM T', 'FDM Extrap T', 'FDM Exact % Error', 'FDM Exact B', 'FDM Extrap % Error', 'FDM Extrap B', 'FEM T', 'FEM Extrap T', 'FEM Exact % Error', 'FEM Exact B', 'FEM Extrap % Error', 'FEM Extrap B'] + df = pd.DataFrame(table_data, index=[f'N = {i}' for i in data_dict["section3"]["num_sections"]], columns=columns) + + print(f"Convergence of FDM and FEM Temperature at x_bar = {x_bar_1:.3f} and k_ratio = {k_ratio}") + print(df.to_string()) + print("\n" * 2) + + # Print the heat convection results for this k ratio + table_data = [] + for i in range(len(data_dict["section3"]["num_sections"])): + table_data.append([ + analy_heat, # True Q + fdm_heats[i], # fdm result + fdm_extrap_heats[i], # fdm extrapolated heat + fdm_heat_errors[i], # fdm % error in reference to analytical + fdm_heat_betas[i], # fdm beta in reference to analytical + fdm_heat_extrap_errors[i], # fdm % error in reference to extrapolated value + fdm_heat_extrap_betas[i], # fdm beta + + fem_heats[i], # fem result + fem_extrap_heats[i], # fem extrapolated heat + fem_heat_errors[i], # fem % error + fem_heat_betas[i], # fem beta + fem_heat_extrap_errors[i], # fem % error + fem_heat_extrap_betas[i] # fem beta + ]) + + columns = ['True Q', 'FDM Q', 'FDM Extrap Q', 'FDM Exact % Error', 'FDM Exact B', 'FDM Extrap % Error', 'FDM Extrap B', 'FEM Q', 'FEM Extrap Q', 'FEM Exact % Error', 'FEM Exact B', 'FEM Extrap % Error', 'FEM Extrap B'] + df = pd.DataFrame(table_data, index=[f'N = {i}' for i in data_dict["section3"]["num_sections"]], columns=columns) + + print(f"Convergence of FDM and FEM Heat Convection with x_bar = {x_bar_1:.3f} and k_ratio = {k_ratio}") + print(df.to_string()) + print("\n" * 2) + + # Add a subplot of the interface temp convergence for this k ratio + ax_temp = axs_temp[idx] + + # Data to plot for temperature convergence + num_sections = data_dict["section3"]["num_sections"] + fdm_exact_errors = [e[0] for e in fdm_interface_temp_errors] + fem_exact_errors = [e[0] for e in fem_interface_temp_errors] + fdm_extrap_errors = fdm_interface_temp_extrap_errors + fem_extrap_errors = fem_interface_temp_extrap_errors + + # Plotting temperature lines + ax_temp.plot(num_sections, fdm_exact_errors, label="FDM Exact") + ax_temp.plot(num_sections, fem_exact_errors, label="FEM Exact") + ax_temp.plot(num_sections, fdm_extrap_errors, label="FDM Extrap") + ax_temp.plot(num_sections, fem_extrap_errors, label="FEM Extrap") + + ax_temp.set_xscale("log") + ax_temp.set_yscale("log") + ax_temp.set_xlabel("Number of Sections (N)") + ax_temp.set_ylabel("% Error") + ax_temp.set_title(f"k_ratio = {k_ratio}") + ax_temp.grid(True) + ax_temp.legend() + + # Add a subplot of the heat convergence for this k ratio + ax_heat = axs_heat[idx] + + # Data to plot for heat convergence + fdm_exact_heat_errors = fdm_heat_errors + fem_exact_heat_errors = fem_heat_errors + fdm_extrap_heat_errors = fdm_heat_extrap_errors + fem_extrap_heat_errors = fem_heat_extrap_errors + + # Plotting heat lines + ax_heat.plot(num_sections, fdm_exact_heat_errors, label="FDM Exact") + ax_heat.plot(num_sections, fem_exact_heat_errors, label="FEM Exact") + ax_heat.plot(num_sections, fdm_extrap_heat_errors, label="FDM Extrap") + ax_heat.plot(num_sections, fem_extrap_heat_errors, label="FEM Extrap") + + ax_heat.set_xscale("log") + ax_heat.set_yscale("log") + ax_heat.set_xlabel("Number of Sections (N)") + ax_heat.set_ylabel("% Error") + ax_heat.set_title(f"k_ratio = {k_ratio}") + ax_heat.grid(True) + ax_heat.legend() + + # Hide any unused subplots (in case of an odd number of k_ratios) + for i in range(num_k_ratios, len(axs_temp)): + fig_temp.delaxes(axs_temp[i]) + fig_heat.delaxes(axs_heat[i]) + + # Adjust layout for better spacing + fig_temp.tight_layout(rect=[0, 0.03, 1, 0.95]) + fig_heat.tight_layout(rect=[0, 0.03, 1, 0.95]) + + # Save the plots + if k == 0: + fig_temp.savefig("images/section3/x_bar_1/temp_convergences.png", dpi=300) + fig_heat.savefig("images/section3/x_bar_1/heat_convergences.png", dpi=300) + else: + fig_temp.savefig("images/section3/x_bar_2/temp_convergences.png", dpi=300) + fig_heat.savefig("images/section3/x_bar_2/heat_convergences.png", dpi=300) + + + +def main(): + # Make directories for plots + import os + import shutil + + base_dir = "images" + sub_dirs = ["section1", "section2", "section3"] + nested_dirs = ["x_bar_1", "x_bar_2"] + + # Create the base directory if it doesn't exist + if not os.path.exists(base_dir): + os.mkdir(base_dir) + + # Create or clear subdirectories and their nested directories + for sub_dir in sub_dirs: + section_path = os.path.join(base_dir, sub_dir) + if os.path.exists(section_path): + # Remove all contents of the directory + shutil.rmtree(section_path) + # Create the section directory + os.makedirs(section_path) + + # Create nested directories within each section + for nested_dir in nested_dirs: + nested_path = os.path.join(section_path, nested_dir) + os.makedirs(nested_path) + + # import sys + + # # Redirect all print statements to a file + # sys.stdout = open("output.txt", "w", encoding="utf-8") + + section_1() # Analytical solutions + section_2() # FDM and FEM temperature and heat convection + section_3() # Convergence of FDM and FEM temperature and heat convection + + # # Restore original stdout and close the file + # sys.stdout.close() + # sys.stdout = sys.__stdout__ + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/HW3/output.txt b/HW3/output.txt new file mode 100644 index 0000000..722312c Binary files /dev/null and b/HW3/output.txt differ diff --git a/README.md b/README.md new file mode 100644 index 0000000..e69de29