<jupyter_start><jupyter_text>AIAA 2020 Figures for model comparison (Figure 5)<jupyter_code>import os
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib as mpl
import matplotlib.pyplot as plt
import pickle
from pandas import Timestamp

import matplotlib.dates as mdates 

from mhdpy import *

import mhdpy

finalanalysisfolder = mhdpy.fp.gen_path('mhdlab','Analysis', 'Lee', 'MHDGEN','Final')
dsst = mhdpy.load.loadprocesseddata(os.path.join(finalanalysisfolder, 'Data','dsst'))

ds_p = xr.load_dataset(os.path.join(finalanalysisfolder, 'Data','spectral_data', 'fit_params.cdf'))
ds_p_stderr = xr.load_dataset(os.path.join(finalanalysisfolder, 'Data','spectral_data', 'fit_params_stderr.cdf'))
ds_cantera = xr.load_dataset(os.path.join(finalanalysisfolder,'Data', 'ds_cantera.cdf'))
ds_cfd = xr.load_dataset(os.path.join(finalanalysisfolder,'Data', 'ds_cfd.cdf'))

#TODO: consolidate downselection of all datasets
ds_cfd_sel = ds_cfd.sel(tf=13, method='nearest')

with open(os.path.join(finalanalysisfolder,'Data', 'da_ct.pickle'), 'rb') as file:
    da_ct = pickle.load(file)<jupyter_output><empty_output><jupyter_text>Modeling results  Wall heat transfersin cantera ds because the were the ones used in the calculation...<jupyter_code>ds_calor = ds_cantera[['wallheattransfer', 'chan_ht', 'comb_ht']]

wht_mean = ds_calor['wallheattransfer'].mean('date', keep_attrs=True)
wht_std = ds_calor['wallheattransfer'].std('date', keep_attrs=True)
plot.common.xr_errorbar(wht_mean, wht_std, huedim='tf')
da_calor_sel_mean = ds_calor.to_array('var').mean('date', keep_attrs=True).sel(tf=12.9, method='nearest').drop('tf')
da_calor_sel_std = ds_calor.to_array('var').std('date', keep_attrs=True).sel(tf=12.9, method='nearest').drop('tf')
fig ,axes = plot.common.xr_errorbar(da_calor_sel_mean,da_calor_sel_std, huedim='var')

# axes.set_ylabel('Heat Transfer')

axes.set_ylim(0)

axes.get_legend().set_bbox_to_anchor([0,0.1,1.2,1])

plt.xscale('log')
fig ,axes = plot.common.xr_errorbar(da_calor_sel_mean.drop('wallheattransfer', 'var'),da_calor_sel_std.drop('wallheattransfer', 'var'), huedim='var')

# axes.set_ylabel('Heat Transfer')

axes.set_ylim(0)

handles, labels = axes.get_legend_handles_labels()
axes.legend(handles, ['Channel', 'Combustor'])

axes.get_legend().set_bbox_to_anchor([0,0.1,1.2,1])

plt.xscale('log')
plt.ylabel('$\dot{Q} \  [kW]$')
ds_cfd_sel[['chan_ht', 'comb_ht']].to_array('var').plot(hue='var', marker='s')
fig ,axes = plot.common.xr_errorbar(da_calor_sel_mean.drop('wallheattransfer', 'var'),da_calor_sel_std.drop('wallheattransfer', 'var'), huedim='var')

handles_1, labels_1 = axes.get_legend_handles_labels()
plt.gca().set_prop_cycle(None)

lns = ds_cfd_sel[['chan_ht', 'comb_ht']].to_array('var').plot(hue='var', marker='s', ax=axes,)

# axes.set_ylabel('Heat Transfer')

axes.set_ylim(0)

# 
axes.legend([*handles_1, *lns], ['$\dot{Q}_{Channel, expt}$', '$\dot{Q}_{Combustor, expt}$', '$\dot{Q}_{Channel, CFD}$', '$\dot{Q}_{Combustor, CFD}$'], ncol=2)

axes.get_legend().set_bbox_to_anchor([0,0.1,1.05,1.02])

plt.xscale('log')
plt.xlabel('K wt%')
plt.ylabel('$\dot{Q} \  [kW]$')
axes.set_title('')
fig.set_size_inches(7,5)
plt.ylim(0,27)<jupyter_output><empty_output><jupyter_text>Outlet Temperature<jupyter_code>ds_cantera['T'].mean('date',keep_attrs=True).plot(hue='tf', marker ='o')
ds_cantera['T'].mean('date',keep_attrs=True).sel(tf=13, method='nearest').plot(marker ='o', label='Cantera')
plt.legend()
plt.xscale('log')
ds_cantera['T'].mean('date',keep_attrs=True).sel(tf=13, method='nearest').plot(marker ='o', label='Cantera', color='r')
ds_cfd_sel['T_outlet'].plot(label='CFD', marker ='s', color='g')
# plt.legend()
plt.xscale('log')
plt.xlabel('K wt%')
plt.ylabel('$T_{outlet} \ [K]$')
plt.gca().set_title('')
plt.gcf().set_size_inches(7,5)<jupyter_output>C:\Users\aspit\anaconda3\lib\site-packages\xarray\core\nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Potassium number denisty<jupyter_code>ds_cantera['nK_cant'].mean('date',keep_attrs=True).plot(hue='tf', marker ='o')
ds_cantera['sigma'].mean('date',keep_attrs=True).plot(hue='tf', marker ='o')<jupyter_output>C:\Users\aspit\anaconda3\lib\site-packages\xarray\core\nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Modeling and Experiment Comparison for TF = 13 g/s<jupyter_code>nK_expt = ds_p['nK_m3']
nK_expt.name ='nK_expt'
nK_expt_err = ds_p_stderr['nK_m3']

da_nK_out = ds_cantera['nK_cant'].sel(tf=12.9, method='nearest').mean('date', keep_attrs=True)
nK_expt_mean = nK_expt.mean('date', keep_attrs=True).sel(tf=12.9, method='nearest')
nK_expt_std = nK_expt.std('date', keep_attrs=True).sel(tf=12.9, method='nearest')

fig, axes = plot.common.xr_errorbar(nK_expt_mean, nK_expt_std, label = 'Experiment')
da_nK_out.plot(marker='o', ax= axes, label = 'Cantera', color='r')
plt.legend()

plt.gca().set_title('')
plt.yscale('log')
plt.xscale('log')

plt.ylabel('$n_K [\\#/m^3]$')
nK_expt_mean = nK_expt.mean('date', keep_attrs=True).sel(tf=12.9, method='nearest')
nK_expt_std = nK_expt.std('date', keep_attrs=True).sel(tf=12.9, method='nearest')

fig, axes = plot.common.xr_errorbar(nK_expt_mean, nK_expt_std, label = 'Experiment')
da_nK_out.plot(marker='o', ax= axes, label = 'Cantera', color='r')
ds_cfd_sel['nK'].plot(marker='s', ax= axes, label = 'CFD', color='g')
plt.legend()

plt.gca().set_title('')
plt.yscale('log')
plt.xscale('log')
plt.ylim(5e20)


plt.ylabel('$n_K [\\#/m^3]$')
plt.xlabel('K wt%')
plt.gcf().set_size_inches(7,5)<jupyter_output>C:\Users\aspit\anaconda3\lib\site-packages\xarray\core\nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>calculate experimental resistance vs experimental nK<jupyter_code>ds = analysis.ct.assign_tc_general(dsst['nhr'],da_ct)

da_mean, da_std = analysis.gen.bin_gen(ds, curname='I', voltname='V')

da_mean_highV = da_mean.sel(voltage=slice(50,100))#.drop(0.0, 'Kwt')

resist = da_mean_highV.coords['voltage']/da_mean_highV
resist = resist.mean('voltage', keep_attrs=True)
resist.attrs = dict(long_name = 'Resistance (50-100V) ', units = 'ohms')
resist = resist.where(resist>0)
resist = resist.where(resist<5000)
resist.name = 'resistance'
ds = xr.merge([resist,nK_expt])#.mean('date', keep_attrs=True).plot.scatter(y='resistance', x='nK')
ds.plot.scatter(x='nK_expt', y='resistance')<jupyter_output><empty_output><jupyter_text>Assign nK as a coordinate to resistance<jupyter_code>seldicts = analysis.xr.gen_seldicts(ds['resistance'].isel(Kwt=0))
seldicts

rs = []

for seldict in seldicts:
    nK_sel = ds['nK_expt'].sel(seldict)
    r_sel = ds['resistance'].sel(seldict)
    r_sel = r_sel.assign_coords(Kwt=nK_sel.values).dropna('Kwt','all').rename(Kwt='nK_expt')
    r_sel = r_sel.expand_dims('date').expand_dims('tf')
    rs.append(r_sel)
    
r_nK = xr.merge(rs)['resistance']
r_nK.coords['nK_expt'].attrs = ds['nK_expt'].attrs
g = r_nK.mean('date', keep_attrs=True).plot(hue='tf', marker='o')<jupyter_output>C:\Users\aspit\anaconda3\lib\site-packages\xarray\core\nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Bin the nK coordinate<jupyter_code>#Test out calculating coord bin standard deviations without having to return it yet in fix_coord_grid_bins and break code
from scipy.stats import binned_statistic

# bins = [0,3e20,9e20,1e22]
raw_coords = r_nK.coords['nK_expt'].values
bins = [0, *np.logspace(np.log10(raw_coords.min())+0.5, np.log10(raw_coords.max())+0.1, 3)]
bins

bin_means, bin_edges, binnumber = binned_statistic(raw_coords, raw_coords, 'mean', bins)
bin_means = bin_means[~np.isnan(bin_means)]
bin_std, bin_edges, binnumber = binned_statistic(raw_coords, raw_coords, 'std', bins)
bin_std = bin_std[~np.isnan(bin_std)]


r_nK_bin = analysis.xr.fix_coord_grid(r_nK, 'nK_expt',  bin_means, keep_attrs=True)
r_nK_bin, r_nK_bin_stack = analysis.ct.reset_da_ct(r_nK_bin, keep_attrs=True)
r_nK_bin
r_nK_mean = r_nK_bin.mean('date', keep_attrs=True)
r_nK_std = r_nK_bin.std('date', keep_attrs=True)
plot.common.xr_errorbar(r_nK_mean, r_nK_std, huedim='tf')
plt.xscale('log')
conductance = 1/r_nK_bin
conductance.name = 'conductance'
conductance.attrs = dict(long_name='Conductance (50-100V)', units = 'S')

conductance_mean = conductance.mean('date', keep_attrs=True)
conductance_std = conductance.std('date', keep_attrs=True)
plot.common.xr_errorbar(conductance_mean, conductance_std, huedim='tf')
plt.xscale('log')<jupyter_output><empty_output><jupyter_text>Plot with x error have to downselect to one total flow case right now<jupyter_code>r_nK_mean = r_nK_bin.mean('date', keep_attrs=True).sel(tf=12.9,method='nearest')
r_nK_std = r_nK_bin.std('date').sel(tf=12.9,method='nearest')
plot.common.xr_errorbar(r_nK_mean, r_nK_std, xerr=bin_std)
plt.xscale('log')
r_nK_bin
#qucikly make consistent with other plots...hacky....

r_nK_bin_kohm = r_nK_bin/1000
r_nK_bin_kohm.attrs = dict(long_name='$R_{expt}$ (50-100V)', units ='Kohm')

r_nK_kohm_mean = r_nK_bin_kohm.mean('date', keep_attrs=True).sel(tf=12.9,method='nearest')
r_nK_kohm_std = r_nK_bin_kohm.std('date').sel(tf=12.9,method='nearest')
plot.common.xr_errorbar(r_nK_kohm_mean, r_nK_kohm_std, xerr=bin_std)
plt.xscale('log')
conductance = 1/r_nK_bin
conductance.name = 'conductance'
conductance.attrs = dict(long_name='Conductance (50-100V)', units = 'S')

conductance_mean = conductance.mean('date', keep_attrs=True).sel(tf=12.9,method='nearest')
conductance_std = conductance.std('date', keep_attrs=True).sel(tf=12.9,method='nearest')
plot.common.xr_errorbar(conductance_mean, conductance_std, xerr=bin_std)
plt.xscale('log')<jupyter_output><empty_output><jupyter_text>calculate cantera resistance vs cantera nKthis requires binning with cantera unlike before, because now the values have been calculated for each separate measured heat transfer<jupyter_code>seldicts = analysis.xr.gen_seldicts(ds_cantera.isel(Kwt=0))

dss = []

for seldict in seldicts:
    nK_sel = ds_cantera['nK_cant'].sel(seldict)
    ds_sel = ds_cantera.drop('nK_cant').sel(seldict)
    ds_sel = ds_sel.assign_coords(Kwt=nK_sel.values).dropna('Kwt','all').rename(Kwt='nK_cant')
    ds_sel = ds_sel.expand_dims('date').expand_dims('tf')
    dss.append(ds_sel)
    
ds_cantera_nK = xr.merge(dss)

ds_cantera_nK = analysis.xr.fix_coord_grid_bins(ds_cantera_nK, 'nK_cant', [0,1e21, 5e21, 1e23])

ds_cantera_nK, ds_cantera_nK_stack = analysis.ct.reset_da_ct(ds_cantera_nK, keep_attrs=True)
ds_cantera_nK.coords['nK_cant'].attrs = ds_cantera['nK_cant'].attrs
fig, axes = plot.common.xr_errorbar(conductance_mean, conductance_std, xerr=bin_std)
axes.legend(['Experiment\n (Statistics over dates)'])
axes.get_legend().set_bbox_to_anchor([0,0,1.9,1])
plt.xscale('log')
plt.ylim(0,4e-3)

t = plt.gca().twinx()
# ds_cantera_nK['sigma'].mean('date', keep_attrs=True).sel(tf=12.9,method='nearest').plot(ax=t, marker ='o')
ds_cantera_nK['sigma'].mean('date', keep_attrs=True).sel(tf=12.9, method='nearest').plot(ax=t, color='r', marker ='o', label = 'Cantera Conductivity')
# t.set_ylim(0,5.1)
t.arrow(1e21,1,1e20,0, color = 'red', width = 0.1, head_length = 5e19)
t.legend(bbox_to_anchor=[0,0,1.9,0.5])
t.set_title('')<jupyter_output>C:\Users\aspit\anaconda3\lib\site-packages\xarray\core\nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Interploate Attempt cantera conductance<jupyter_code>ds_cantera['G_cantera']
G_cantera = ds_cantera_nK['G_cantera']
R_cantera = ds_cantera_nK['R_cantera']
# G.coords['nK_cant'].attrs = dict(long_name='$n_K$', units = '$\\#/m^3$')<jupyter_output><empty_output><jupyter_text>This match was made on the first try (using the same parameters at the AIAA abstract)<jupyter_code>fig, axes = plot.common.xr_errorbar(conductance_mean, conductance_std, xerr=bin_std, label = 'Experiment TF=12.4')
G_cantera.mean('date', keep_attrs=True).sel(tf=12.9, method='nearest').plot(color='r', marker ='o', label = 'Cylindrical Model')
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.ylim(0,0.03)
axes.set_title('')
plt.xlabel('$n_K \ [\#/m^3]$')
fig, axes = plot.common.xr_errorbar(r_nK_mean, r_nK_std, xerr=bin_std, label='Experiment')
R_cantera.mean('date', keep_attrs=True).sel(tf=12.9, method='nearest').plot(color='r', marker ='o', label = 'Cantera')
plt.xscale('log')
plt.yscale('log')
# plt.legend(bbox_to_anchor=[0,0,1,0.9])
# plt.ylim(0,0.03)
axes.set_title('')
plt.xlabel('$n_K \ [\#/m^3]$')
plt.ylabel('Resistance [Ohm]')
plt.ylim(1e2,1.5e3)<jupyter_output>C:\Users\aspit\anaconda3\lib\site-packages\xarray\core\nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Calculate CFD resistance vs nK<jupyter_code>ds_cfd_sel['resistance'].plot(marker='s', label = 'CFD', color='g')
ds_cfd_sel['nK'].plot(marker='s', label = 'CFD', color='g')
seldicts = analysis.xr.gen_seldicts(ds_cfd.isel(Kwt=0))

dss = []

for seldict in seldicts:
    nK_sel = ds_cfd['nK'].sel(seldict)
    ds_sel = ds_cfd.drop('nK').sel(seldict)
    ds_sel = ds_sel.assign_coords(Kwt=nK_sel.values).dropna('Kwt','all').rename(Kwt='nK')
    ds_sel = ds_sel.expand_dims('date').expand_dims('tf')
    dss.append(ds_sel)
    
ds_cfd_nK = xr.merge(dss)

ds_cfd_nK = analysis.xr.fix_coord_grid_bins(ds_cfd_nK, 'nK', [0,2e21, 5e21, 1e23])

ds_cfd_nK, ds_cfd_nK_stack = analysis.ct.reset_da_ct(ds_cfd_nK, keep_attrs=True)
ds_cfd_nK.coords['nK'].attrs = ds_cfd['nK'].attrs
fig, axes = plot.common.xr_errorbar(r_nK_mean, r_nK_std, xerr=bin_std, label='Experiment')
R_cantera.mean('date', keep_attrs=True).sel(tf=12.9, method='nearest').plot(color='r', marker ='o', label = 'Cylindrical Model')
ds_cfd_nK.sel(tf=12.9, method='nearest')['resistance'].plot(marker='s', label = 'CFD', color='g', ax=axes)
plt.xscale('log')
plt.yscale('log')
# plt.legend(bbox_to_anchor=[0,0,1,0.9])
plt.ylim(1e2,1.5e3)
axes.set_title('')

plt.ylabel('Resistance [Ohm]')
plt.gcf().set_size_inches(7,5)
axes.set_xlabel('$n_K \ [\#/m^3]$')<jupyter_output>C:\Users\aspit\anaconda3\lib\site-packages\xarray\core\nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)
{'fraction_non_alphanumeric': 0.18193294992285033, 'fraction_numerical': 0.020339458549586197, 'mean_word_length': 7.856521739130435, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '383515', 'license': 'MIT', 'n_tokens_mistral': 6326, 'n_tokens_neox': 6063, 'n_words': 813, 'path': 'public_data/AIAA 2020/.ipynb_checkpoints/Figures_ModelingComparison-checkpoint.ipynb', 'repo_name': 'MHDLab/MHDLab.github.io', 'size': '14258'}