Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
229 changes: 220 additions & 9 deletions Scripts/MaCFP-4/Cone_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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))
Expand All @@ -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


Expand Down Expand Up @@ -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)'],
Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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)


# ------------------------------------
Expand Down Expand Up @@ -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)



Expand Down Expand Up @@ -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")


Expand All @@ -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)
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)

70 changes: 70 additions & 0 deletions Scripts/MaCFP-4/Elemental_composition.py
Original file line number Diff line number Diff line change
@@ -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')
1 change: 0 additions & 1 deletion Scripts/MaCFP-4/Utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
2 changes: 1 addition & 1 deletion Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R1.csv
Original file line number Diff line number Diff line change
@@ -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
Expand Down
2 changes: 1 addition & 1 deletion Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_40kW_R2.csv
Original file line number Diff line number Diff line change
@@ -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
Expand Down
2 changes: 1 addition & 1 deletion Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R1.csv
Original file line number Diff line number Diff line change
@@ -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
Expand Down
2 changes: 1 addition & 1 deletion Wood/Calibration_Data/FSRI/FSRI_Wood_CAPA_N2_60kW_R2.csv
Original file line number Diff line number Diff line change
@@ -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
Expand Down