diff --git a/pyvisco/prony.py b/pyvisco/prony.py index 2a62cf9..d832878 100644 --- a/pyvisco/prony.py +++ b/pyvisco/prony.py @@ -459,24 +459,41 @@ def fit_freq(df_dis, df_master, opt=False): #Estimate instantaneous (E_0) and equilibrium (E_inf) modulus E_0 = df_dis.E_0 + if opt: + #Get measurement data + E_freq_meas = np.concatenate((df_master[stor]/E_0, + df_master[loss]/E_0)) + omega_meas = df_master['omega'].values - #Get measurement data - E_freq_meas = np.concatenate((df_master[stor]/E_0, - df_master[loss]/E_0)) - omega_meas = df_master['omega'].values - - #Get Prony series - tau_i = df_dis['tau_i'].values - x0 = np.hstack((alpha_i, tau_i)) + #Get Prony series + tau_i = df_dis['tau_i'].values + x0 = np.hstack((alpha_i, tau_i)) - #Define bounds - bnd_a = ((0,1),)*alpha_i.shape[0] + #Define bounds + bnd_a = ((0,1),)*alpha_i.shape[0] - #Perform minimization to obtain alpha_i - res = minimize(ls_res(E_freq_norm), alpha_i, - args=(tau_i, E_freq_meas, omega_meas), method='L-BFGS-B', bounds=bnd_a) - alpha_i = res.x - err = res.fun + #Perform minimization to obtain alpha_i + res = minimize(ls_res(E_freq_norm), alpha_i, + args=(tau_i, E_freq_meas, omega_meas), method='L-BFGS-B', bounds=bnd_a) + alpha_i = res.x + err = res.fun + else: + E_inf = df_dis.E_inf + #Assembly 'K_global' matrix [Kraus 2017, Eq. 22] + N = df_dis.nprony + K_stor = np.tril(np.ones((N,N)), -1) + np.diag([0.5] * N) + K_loss = (np.diag([0.5] * N) + + np.diag([0.1] * (N-1), 1) + np.diag([0.1] * (N-1), -1) + + np.diag([0.01] * (N-2), 2) + np.diag([0.01] * (N-2), -2) + + np.diag([0.001] * (N-3), 3) + np.diag([0.001] * (N-3), -3)) + K_global = np.vstack([K_stor, K_loss, np.ones((1,N))]) + #Assembly right-hand vector + E = np.concatenate((df_dis[stor]/(E_0-E_inf), + df_dis[loss]/(E_0-E_inf), + np.array([1]))) + + #Solve equation system + alpha_i, err = nnls(K_global, E) ############################################################################### # DEPRECATED - Error gets large for certain viscoelastic materials diff --git a/pyvisco/shift.py b/pyvisco/shift.py index d9b1daf..27570ab 100644 --- a/pyvisco/shift.py +++ b/pyvisco/shift.py @@ -204,10 +204,10 @@ def fit_poly(df_aT): xdata = df_aT['T'].values+273.15 ydata = df_aT['log_aT'].values - df_K['D4'][['C0', 'C1', 'C2', 'C3', 'C4']], _pcov = curve_fit(poly4, xdata, ydata) - df_K['D3'][['C0', 'C1', 'C2', 'C3']], _pcov = curve_fit(poly3, xdata, ydata) - df_K['D2'][['C0', 'C1', 'C2']], _pcov = curve_fit(poly2, xdata, ydata) - df_K['D1'][['C0', 'C1']], _pcov = curve_fit(poly1, xdata, ydata) + df_K.loc[['C0', 'C1', 'C2', 'C3', 'C4'], 'D4'], _pcov = curve_fit(poly4, xdata, ydata) + df_K.loc[['C0', 'C1', 'C2', 'C3'], 'D3'], _pcov = curve_fit(poly3, xdata, ydata) + df_K.loc[['C0', 'C1', 'C2'], 'D2'], _pcov = curve_fit(poly2, xdata, ydata) + df_K.loc[['C0', 'C1'], 'D1'], _pcov = curve_fit(poly1, xdata, ydata) # Celsius df_C = pd.DataFrame( @@ -220,10 +220,10 @@ def fit_poly(df_aT): xdata = df_aT['T'].values ydata = df_aT['log_aT'].values - df_C['D4'][['C0', 'C1', 'C2', 'C3', 'C4']], _pcov = curve_fit(poly4, xdata, ydata) - df_C['D3'][['C0', 'C1', 'C2', 'C3']], _pcov = curve_fit(poly3, xdata, ydata) - df_C['D2'][['C0', 'C1', 'C2']], _pcov = curve_fit(poly2, xdata, ydata) - df_C['D1'][['C0', 'C1']], _pcov = curve_fit(poly1, xdata, ydata) + df_C.loc[['C0', 'C1', 'C2', 'C3', 'C4'], 'D4'], _pcov = curve_fit(poly4, xdata, ydata) + df_C.loc[['C0', 'C1', 'C2', 'C3'], 'D3'], _pcov = curve_fit(poly3, xdata, ydata) + df_C.loc[['C0', 'C1', 'C2'], 'D2'], _pcov = curve_fit(poly2, xdata, ydata) + df_C.loc[['C0', 'C1'], 'D1'], _pcov = curve_fit(poly1, xdata, ydata) # Prep dataframes and add units for output df_C = df_C.T diff --git a/verification/verify_freq_master.ipynb b/verification/verify_freq_master.ipynb index a057805..3128e28 100644 --- a/verification/verify_freq_master.ipynb +++ b/verification/verify_freq_master.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -29,7 +29,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -45,7 +45,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -55,7 +55,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -72,7 +72,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -111,7 +111,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -149,7 +149,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -186,7 +186,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -223,7 +223,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -259,7 +259,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -299,7 +299,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -341,7 +341,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -429,7 +429,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": null, "metadata": { "scrolled": false },