diff --git a/Scripts/MaCFP-4/Cone_analysis.py b/Scripts/MaCFP-4/Cone_analysis.py index 40f0490..3b03f5c 100644 --- a/Scripts/MaCFP-4/Cone_analysis.py +++ b/Scripts/MaCFP-4/Cone_analysis.py @@ -76,7 +76,7 @@ def set_plot_style(): #region functions # ------------------------------------ def average_cone_series(series_name: str)->pd.DataFrame: - """Calculate average mass and MLR for a test series.""" + """Calculate average mass and HRR for a test series.""" paths = list(DATA_DIR.glob(f"*/*{series_name}_[rR]*.csv")) paths = [p for p in paths if "TEMPLATE" not in str(p)] paths = [p for p in paths if p in Cone_Data] @@ -135,6 +135,68 @@ def average_cone_series(series_name: str)->pd.DataFrame: return df_average + +def average_Gas_series(series_name: str)->pd.DataFrame: + """Calculate average mass and MLR for a test series.""" + paths = list(DATA_DIR.glob(f"*/*{series_name}_[rR]*.csv")) + paths = [p for p in paths if "TEMPLATE" not in str(p)] + paths = [p for p in paths if p in Gasification_Data] + + Dataframes = [] + if len(paths) == 0: + raise Exception((f"No files found for series {series_name}", "red")) + + # Read data + for i, path in enumerate(paths): + df_raw = pd.read_csv(path) + + t_floor = df_raw["Time (s)"].iloc[0] + t_floor = np.ceil(t_floor) + t_ceil = df_raw["Time (s)"].iloc[-1] + t_ceil = np.floor(t_ceil) + + InterpT = np.arange(t_floor, t_ceil+1, 1) + length = len(InterpT) + df_interp = pd.DataFrame(index=range(length)) + for columns in df_raw.columns[:]: + df_interp[columns] = np.interp( + InterpT, df_raw["Time (s)"], df_raw[columns] + ) + df_interp = Calculate_dm_dt(df_interp) + Dataframes.append(df_interp) + + + merged_df = Dataframes[0] + for df in Dataframes[1:]: + merged_df = pd.merge( + merged_df, + df, + on="Time (s)", + how="outer", + suffixes=("", f" {int(len(merged_df.columns)/2+0.5)}"), + ) + + merged_df.rename(columns={'Mass (g)': "Mass (g) 1"}, inplace=True) + merged_df.rename(columns={'dm/dt': "dm/dt 1"}, inplace=True) + + #average + time_cols = merged_df.filter(regex=r'^Time \(s\)').columns + mass_cols = merged_df.filter(regex=r'^Mass \(g\)').columns + dmdt_cols = merged_df.filter(regex=r'^dm/dt').columns + + df_average = pd.DataFrame({'Time (s)': merged_df['Time (s)']}) + n=2 + sum = merged_df[dmdt_cols].rolling(2*n+1, min_periods=1,center=True).sum().sum(axis=1) + cnt = merged_df[dmdt_cols].rolling(2*n+1, min_periods=1,center=True).count().sum(axis=1) + df_average['dm/dt'] = sum / cnt # Series: mean of all non-NaN values in rows i-2..i+2 across all columns + + diff = merged_df[dmdt_cols].sub(df_average['dm/dt'], axis=0)**2 + sum_diff = diff.rolling(2*n+1, min_periods=1,center=True).sum().sum(axis=1) + df_average['unc dm/dt'] = np.sqrt(sum_diff/(cnt*(cnt-1))) + + return df_average + + def calculate_int_HRR(df:pd.DataFrame)->pd.DataFrame: """Calculate integral HRR.""" total_hrr = np.zeros(len(df)) @@ -149,7 +211,6 @@ def Calculate_dm_dt(df:pd.DataFrame): dt = df['Time (s)'].shift(-2) - df['Time (s)'].shift(2) df['dm/dt'] = (df['Mass (g)'].shift(2) - df['Mass (g)'].shift(-2)) / dt df['dm/dt'] = df['dm/dt'].interpolate(method='linear', limit_direction='both') #avoid nan_values - print(df) return df @@ -211,6 +272,9 @@ def Calculate_dm_dt(df:pd.DataFrame): ax_rate = ax_HRR.twinx() df_average = average_cone_series(set) + Duck, color = label_def(set.split('_')[0]) + Conditions = '_'.join(set.split('_')[2:]) + # plot average # Plot mass (left y-axis) ax_HRR.plot(df_average['Time (s)'], df_average['HRR (kW/m2)'], @@ -268,8 +332,8 @@ def Calculate_dm_dt(df:pd.DataFrame): ax_HRR.set_xlabel('Time (s)') ax_HRR.set_ylabel('HRR (kW/m2)') - # Figure title - fig_title = set + # Figure title + plt.title(Duck+"\n"+Conditions) # Legend fig.legend() @@ -337,12 +401,41 @@ def Calculate_dm_dt(df:pd.DataFrame): ax1.set_ylabel('Temperature [K]') fig1.tight_layout() ax1.legend() - - if dev == 'Cone': - fig1.savefig(str(base_dir) + '/Cone/Cone_{}_{}_{}_BackT.{}'.format(material, flux,orient,ex)) + + fig1.savefig(str(base_dir) + '/Cone/Cone_{}_{}_{}_BackT.{}'.format(material, flux,orient,ex)) plt.close(fig1) +# Back side temperature plots for all unique institutions, atmospheres and heating rates (when available) +linestyle = ['-',':','-.'] +colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', "#a686c4", '#8c564b'] +for idx,set in enumerate(Cone_sets): + fig1, ax1 = plt.subplots(figsize=(6, 4)) + paths_CONE_set = list(DATA_DIR.glob(f"*/{set}_[rR]*.csv")) + color_counter = 0 + Duck,x = label_def(set.split('_')[0]) + Conditions = '_'.join(set.split('_')[2:]) + for path in paths_CONE_set: + label, color = label_def(path.stem.split('_')[0]) + df = pd.read_csv(path) + for i in range(1, 4): # Check for Temperature 1, 2, 3 + temp_col = f'TC back {i} (K)' + if temp_col in df.columns: + ax1.plot(df['Time (s)'], df[temp_col], label=label, color=colors[color_counter], linestyle = linestyle[i-1]) + if 'TC Top (K)' in df.columns: + ax1.plot(df['Time (s)'], df['TC Top (K)'], label=label, color=colors[color_counter], dashes=[5, 10]) + color_counter = color_counter+1 + + ax1.set_ylim(bottom=250, top=1200) + ax1.set_xlabel('Time [s]') + ax1.set_ylabel('Temperature [K]') + + # Figure title + ax1.set_title(Duck+"\n"+Conditions) + fig1.tight_layout() + fig1.savefig(str(base_dir) + '/Cone/Individual/Cone_{}_BackT.{}'.format(set,ex)) + + plt.close(fig1) # ------------------------------------ @@ -391,6 +484,43 @@ def Calculate_dm_dt(df:pd.DataFrame): +# Indivdual Mass and mass loss rate plots for all unique institutions atmospheres and heating rates (gasification) +for set in Gas_sets: + fig1, ax1 = plt.subplots(figsize=(6, 4)) + fig2, ax2 = plt.subplots(figsize=(6, 4)) + paths_GAS_set = list(DATA_DIR.glob(f"*/{set}_[rR]*.csv")) + for path in paths_GAS_set: + df_raw = pd.read_csv(path) + df=Calculate_dm_dt(df_raw) + label, color = label_def(path.stem.split('_')[0]) + if institute == 'TIFP+UCT': + ax1.plot(df['Time (s)'],savgol_filter(df['dm/dt']/0.01,41,3),'-', label = label, color=color) + elif institute == 'FSRI': + ax1.plot(df['Time (s)'],savgol_filter(df['dm/dt']/0.00385,41,3),'-', label = label, color=color) + ax2.plot(df['Time (s)'], df['Mass (g)'], '.', label = label, color=color) + + ax1.set_ylim(bottom=0) + ax1.set_xlabel('Time [s]') + ax1.set_ylabel('Mass loss rate [g s$^{-1}$ m$^{-2}$]') + fig1.tight_layout() + handles1, labels1 = ax1.get_legend_handles_labels() + by_label1 = dict(zip(labels1, handles1)) + ax1.legend(by_label1.values(), by_label1.keys()) + + ax2.set_ylim(bottom=0) + ax2.set_xlabel('Time [s]') + ax2.set_ylabel('Mass [g]') + fig2.tight_layout() + handles2, labels2 = ax2.get_legend_handles_labels() + by_label2 = dict(zip(labels2, handles2)) + ax2.legend(by_label2.values(), by_label2.keys()) + + fig1.savefig(str(base_dir) + '/Cone/Gasification_{}_{}_MLR.{}'.format(material, flux,ex)) + fig2.savefig(str(base_dir) + '/Cone/Gasification_{}_{}_Mass.{}'.format(material, flux,ex)) + + + plt.close(fig1) + plt.close(fig2) @@ -449,7 +579,7 @@ def Calculate_dm_dt(df:pd.DataFrame): label = label_def(path.stem.split('_')[0])[0] +' ' df_raw = pd.read_csv(path) df=Calculate_dm_dt(df_raw) - ax1.plot(df['Time (s)'],df['TC Back (K)'],'-', label = label, color='#aec7e8') + ax1.plot(df['Time (s)'],df['TC back 1 (K)'],'-', label = label, color='#aec7e8') ax1.plot(df['Time (s)'],df['TC Top (K)'],'.', label = label + 'Top', color="#bcbd22") @@ -464,4 +594,85 @@ def Calculate_dm_dt(df:pd.DataFrame): fig1.savefig(str(base_dir) + '/Cone/Gasification_{}_{}_{}_BackT.{}'.format(material, flux,orient,ex)) - plt.close(fig1) \ No newline at end of file + plt.close(fig1) + + +# Back side temperature plots for all unique institutions, atmospheres and heating rates (when available) +linestyle = ['-',':','-.'] +colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', "#a686c4", '#8c564b'] +for idx,set in enumerate(Gas_sets): + fig1, ax1 = plt.subplots(figsize=(6, 4)) + paths_Gas_set = list(DATA_DIR.glob(f"*/{set}_[rR]*.csv")) + color_counter = 0 + Duck,x = label_def(set.split('_')[0]) + Conditions = '_'.join(set.split('_')[2:]) + for path in paths_Gas_set: + label, color = label_def(path.stem.split('_')[0]) + df = pd.read_csv(path) + for i in range(1, 4): # Check for Temperature 1, 2, 3 + temp_col = f'TC back {i} (K)' + if temp_col in df.columns: + ax1.plot(df['Time (s)'], df[temp_col], label=label, color=colors[color_counter], linestyle = linestyle[i-1]) + if 'TC Top (K)' in df.columns: + ax1.plot(df['Time (s)'], df['TC Top (K)'], label=label, color=colors[color_counter], dashes=[1, 1]) + color_counter = color_counter+1 + + ax1.set_ylim(bottom=250, top=1200) + ax1.set_xlabel('Time [s]') + ax1.set_ylabel('Temperature [K]') + + # Figure title + ax1.set_title(Duck+"\n"+Conditions) + + fig1.tight_layout() + fig1.savefig(str(base_dir) + '/Cone/Individual/Gasification_{}_BackT.{}'.format(set,ex)) + + plt.close(fig1) + + +# plot average per Gas_set (unique institutions, unique material, unique conditions) +for idx,set in enumerate(Gas_sets): + fig, ax = plt.subplots(figsize=(6, 4)) + df_average = average_Gas_series(set) + + Duck, color = label_def(set.split('_')[0]) + Conditions = '_'.join(set.split('_')[2:]) + + # plot average + # Plot mass (left y-axis) + if institute == 'TIFP+UCT': + area = 0.01 + elif institute == 'FSRI': + area = 0.00385 + ax.plot(df_average['Time (s)'], savgol_filter(df_average['dm/dt']/area,41,3), + label='average MLR', color='limegreen') + ax.fill_between(df_average['Time (s)'], + savgol_filter((df_average['dm/dt']-2*df_average['unc dm/dt'])/area,41,3), + savgol_filter((df_average['dm/dt']+2*df_average['unc dm/dt'])/area,41,3), + color='limegreen', alpha = 0.3) + + + #plot individual + paths_Gas_set = list(DATA_DIR.glob(f"*/{set}_[rR]*.csv")) + for path in paths_Gas_set: + df_raw = pd.read_csv(path) + df=Calculate_dm_dt(df_raw) + ax.plot(df['Time (s)'],savgol_filter(df['dm/dt']/area,41,3), '.',color ='black',markersize=0.001) + + # Set lower limits of both y-axes to 0 + ax.set_ylim(bottom=0) + + # Axes labels + ax.set_xlabel('Time (s)') + ax.set_ylabel('Mass loss rate [g s$^{-1}$ m$^{-2}$]') + + # Figure title + plt.title(Duck+"\n"+Conditions) + + # Legend + fig.legend() + + fig.tight_layout() + plt.savefig(str(base_dir) + f'/Cone/Average/{set}.{ex}') + plt.close(fig) + diff --git a/Scripts/MaCFP-4/Elemental_composition.py b/Scripts/MaCFP-4/Elemental_composition.py new file mode 100644 index 0000000..4542689 --- /dev/null +++ b/Scripts/MaCFP-4/Elemental_composition.py @@ -0,0 +1,70 @@ +import matplotlib.pyplot as plt +import numpy as np +from pathlib import Path + +#region Save plots as pdf or png +ex = 'pdf' #options 'pdf' or 'png + +# TO DO: when prelim document pushed to main repo replace +'../../../matl-db-organizing-committee/' #with +'../../Documents/' + +#region create subdirectories to save plots. +base_dir = Path('../../../matl-db-organizing-committee/SCRIPT_FIGURES') +Composition_dir = base_dir / 'Composition' +Composition_dir.mkdir(parents=True, exist_ok=True) + +# ------------------------------------ +#region data +# ------------------------------------ +species = ['C', 'H', 'N', 'S', 'O'] +virgin = [48.9, 6.1, 0, 0, 40.7] +char = [91.2, 2.8, 0, 0, 3.7] +pyrolyzate = [22.3, 8.6, 0, 0, 63.9] + + + +# ------------------------------------ +#region figure +# ------------------------------------ +fig, ax = plt.subplots(figsize=(5, 3)) + +# Set width of bars and positions +bar_width = 0.25 +x = np.arange(len(species)) + +# Create bars +bars1 = ax.bar(x - bar_width, virgin, bar_width, label='Virgin', color='#FFA500') +bars2 = ax.bar(x, char, bar_width, label='Char', color='#000000') +bars3 = ax.bar(x + bar_width, pyrolyzate, bar_width, label='Pyrolyzate', color='#1E90FF') + +# Add value labels on top of bars +def add_labels(bars): + for bar in bars: + height = bar.get_height() + if height > 0: + ax.text(bar.get_x() + bar.get_width()/2., height, + f'{height:.1f}', + ha='center', va='bottom', fontsize=10) + +add_labels(bars1) +add_labels(bars2) +add_labels(bars3) + +# Customize the plot +ax.set_ylabel('wt (%)', fontsize=11) +ax.set_xlabel('Species', fontsize=11) +ax.set_xticks(x) +ax.set_xticklabels(species) +ax.set_ylim(0, 100) +ax.legend(frameon=True) + +# Add grid for better readability +ax.yaxis.grid(True, linestyle='-', linewidth=0.5, color='gray', alpha=0.3) +ax.set_axisbelow(True) + +# Adjust layout +plt.tight_layout() + +# Save as PDF +plt.savefig(str(base_dir) + '/Composition/Elemental_Composition.{}'.format(ex), format='pdf', bbox_inches='tight') diff --git a/Scripts/MaCFP-4/Utils.py b/Scripts/MaCFP-4/Utils.py index f17f953..640e038 100644 --- a/Scripts/MaCFP-4/Utils.py +++ b/Scripts/MaCFP-4/Utils.py @@ -140,7 +140,6 @@ def make_institution_table( # ---------- Table construction logic ---------- inst_codes = [code for code in CODES if any(counts.get((code, mat,atm, hr), 0) > 0 for mat in materials for atm in atmospheres for hr in heating_rates)] - print(inst_codes) # Case 1: single atmosphere → columns = heating rates if len(atmospheres) == 1 and len(heating_rates) > 1: atm = atmospheres[0] diff --git a/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R1.csv b/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R1.csv index 01f3452..49a4510 100644 --- a/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R1.csv +++ b/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R1.csv @@ -1,4 +1,4 @@ -Time (s),Mass (g),TC Back (K),TC Top (K) +Time (s),Mass (g),TC back 1 (K),TC Top (K) 0,12.377,293.1,416.5 1,12.346,293.2,443 2,12.365,293.2,483.8 diff --git a/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R2.csv b/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R2.csv index 82d1ed0..49342d4 100644 --- a/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R2.csv +++ b/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R2.csv @@ -1,4 +1,4 @@ -Time (s),Mass (g),TC Back (K),TC Top (K) +Time (s),Mass (g),TC back 1 (K),TC Top (K) 0,12.739,297.1,414.8 1,12.713,297.1,432 2,12.734,297.1,492.8 diff --git a/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R1.csv b/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R1.csv index e4e64f1..24f3871 100644 --- a/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R1.csv +++ b/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R1.csv @@ -1,4 +1,4 @@ -Time (s),Mass (g),TC Back (K),TC Top (K) +Time (s),Mass (g),TC back 1 (K),TC Top (K) 0,13.114,296.3,396.4 1,13.134,296.4,539.5 2,13.112,296.2,578.6 diff --git a/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R2.csv b/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R2.csv index 40b4882..2de4e9b 100644 --- a/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R2.csv +++ b/Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R2.csv @@ -1,4 +1,4 @@ -Time (s),Mass (g),TC Back (K),TC Top (K) +Time (s),Mass (g),TC back 1 (K),TC Top (K) 0,11.821,298.6,427.3 1,11.751,298.7,517.6 2,11.742,298.6,564.1