import numpy as np import matplotlib.pyplot as plt from scipy import stats from scipy.interpolate import interp1d # Step 1: Parse data from file with headers and retrieve Reynolds number, radius def parse_file_with_header(filename): with open(filename, 'r') as f: lines = f.readlines() Re_tau = float(lines[0].split("=")[1].strip()) R = lines[3].split("=")[1].strip() R = float(R[0:len(R)-1]) data = np.loadtxt(filename, skiprows=6) y = data[:, 0] U = data[:, 1] return y, U, Re_tau, R # Step 2: Perform linear regression on log-transformed data def fit_power_law(y, U, R): transformed_y = y / R log_y = np.log(transformed_y) log_U = np.log(U) slope, intercept, r_value, p_value, std_err = stats.linregress(log_y, log_U) n = 1 / slope # n is the inverse of the slope U_max = np.exp(intercept) # U_max is exp of the intercept return n, U_max # Step 3: Read multiple datasets and perform regression def analyze_datasets(filenames): results = [] for filename in filenames: y, U, Re_tau, R = parse_file_with_header(filename) n, U_max = fit_power_law(y, U, R) results.append((Re_tau, n, U_max)) return np.array(results) # Example usage for filenames filenames = ["/work/Files_for_Flow/Data/Retau_5k_basic_stats.txt", "/work/Files_for_Flow/Data/Retau_10k_basic_stats.txt", "/work/Files_for_Flow/Data/Retau_20k_basic_stats.txt", "/work/Files_for_Flow/Data/Retau_40k_basic_stats.txt"] results = analyze_datasets(filenames) # Step 4: Interpolate n and U_max at Re_tau = 35061 def interpolate_n_Umax(results, target_Re): Re_taus = results[:, 0] ns = results[:, 1] U_maxs = results[:, 2] # Use linear interpolation for smoother results n_interpolator = interp1d(Re_taus, ns, kind='linear') Umax_interpolator = interp1d(Re_taus, U_maxs, kind='linear') n_35061 = n_interpolator(target_Re) Umax_35061 = Umax_interpolator(target_Re) return n_35061, Umax_35061 # Interpolate for Re_tau = 35061 target_Re = 35061 n_35061, Umax_35061 = interpolate_n_Umax(results, target_Re) print(f"Interpolated n for Re_tau={target_Re}: {n_35061}") print(f"Interpolated U_max for Re_tau={target_Re}: {Umax_35061}") # Step 5: Predict velocity profile for Re=35061 using the interpolated values def predict_velocity_profile(y, n, U_max, R): transformed_y = y / R U_predicted = U_max * (transformed_y ** (1 / n)) return U_predicted # Load the experimental data for validation (file for Re=35061) y_35k, U_35k, Re_tau_35k, R_35k = parse_file_with_header("/work/Files_for_Flow/Data/Retau_35k_basic_stats.txt") # Predict U for Re_tau=35061 using the interpolated n and U_max U_predicted_35k = predict_velocity_profile(y_35k, n_35061, Umax_35061, R=R_35k) # Step 6: Validate against experimental data average_error = np.mean(np.abs(U_predicted_35k - U_35k) / U_35k) * 100 # Error as percentage print(f"Average error in velocity prediction for Re_tau=35061: {average_error:.2f}%") # Enhanced plotting code plt.figure(figsize=(12, 10)) # Adjusted figure size for larger plots # Plot experimental data vs predicted data for Re_tau = 35061 plt.subplot(2, 2, 1) # Two rows, two columns, first subplot plt.plot(y_35k, U_35k, 'o', color="darkblue", markersize=6, label="Experimental Data") plt.plot(y_35k, U_predicted_35k, '-', color="orange", linewidth=2.5, label="Predicted Data") plt.xlabel("Distance from Wall (y/R)", fontsize=12) plt.ylabel("Velocity U (m/s)", fontsize=12) plt.title(f"Velocity Profile Comparison for $Re_\\tau$ = {target_Re}\nAverage Prediction Error = {average_error:.2f}%", fontsize=14) plt.grid(True) plt.legend(loc="upper left") # Plot n vs. Re_tau plt.subplot(2, 2, 2) # Second subplot plt.plot(results[:, 0], results[:, 1], 'bo-', label="n vs Re_tau", markersize=6) plt.xlabel("Re_tau", fontsize=12) plt.ylabel("Exponent n", fontsize=12) plt.title("Exponent n vs Re_tau", fontsize=14) plt.grid(True) # Plot U_max vs. Re_tau plt.subplot(2, 2, 3) # Third subplot plt.plot(results[:, 0], results[:, 2], 'ro-', label="U_max vs Re_tau", markersize=6) plt.xlabel("Re_tau", fontsize=12) plt.ylabel("U_max (m/s)", fontsize=12) plt.title("U_max vs Re_tau", fontsize=14) plt.grid(True) # Adjust layout to avoid overlapping plt.tight_layout() plt.show()