diff --git a/HISTORY.rst b/HISTORY.rst index 804c94c..4fc1aae 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -1,7 +1,11 @@ ======= History ======= - +2024.1.5 -- Bugfix: Thermal conductivity + * If the Helfand moments fit in the thermal conductivity step failed it stopped the + entire job. This is correct, as well as some of the undelying causes for + convergence issues. + 2023.5.29 -- Converged with general approach for trajectory analysis 2023.5.6 -- Bugfix diff --git a/thermal_conductivity_step/analysis.py b/thermal_conductivity_step/analysis.py index 90357fc..68687a3 100644 --- a/thermal_conductivity_step/analysis.py +++ b/thermal_conductivity_step/analysis.py @@ -1,6 +1,8 @@ # -*- coding: utf-8 -*- """Routines to help do Green-Kubo and Helfand moments analysis.""" +import logging +import pprint import warnings import numpy as np @@ -8,6 +10,8 @@ from scipy.optimize import curve_fit, OptimizeWarning import statsmodels.tsa.stattools as stattools +logger = logging.getLogger("thermal_conductivity") + tensor_labels = [ ("xx", "red", "rgba(255,0,0,0.1)"), ("yy", "green", "rgba(0,255,0,0.1)"), @@ -283,7 +287,7 @@ def fit_green_kubo_integral(y, xs, sigma=None): return kappa, kappa_err, a, a_err, tau, tau_err, xs[1:nf], fy -def fit_helfand_moment(y, xs, sigma=None, start=1): +def fit_helfand_moment(y, xs, sigma=None, start=1, logger=logger): """Find the best linear fit to longest possible segment. Parameters @@ -313,7 +317,7 @@ def fit_helfand_moment(y, xs, sigma=None, start=1): # We know the curves curve near the origin, so ignore the first part i = int(start / dx) if i > len(y): - i = int(1 / dx) + i = len(y) // 2 popt, pcov, infodict, msg, ierr = curve_fit( axb, @@ -323,6 +327,14 @@ def fit_helfand_moment(y, xs, sigma=None, start=1): sigma=sigma[i:], absolute_sigma=True, ) + if logger.isEnabledFor(logging.DEBUG): + logger.debug("") + logger.debug(f"{popt=}") + logger.debug(f"{pcov=}") + logger.debug(pprint.pformat(infodict, compact=True)) + logger.debug(f"{msg=}") + logger.debug(f"{ierr=}") + slope = float(popt[0]) b = float(popt[1]) err = float(np.sqrt(np.diag(pcov)[0])) diff --git a/thermal_conductivity_step/thermal_conductivity.py b/thermal_conductivity_step/thermal_conductivity.py index 883b256..aad55ea 100644 --- a/thermal_conductivity_step/thermal_conductivity.py +++ b/thermal_conductivity_step/thermal_conductivity.py @@ -250,7 +250,7 @@ def description_text(self, P=None, short=False): f"Error describing thermal_conductivity flowchart: {e} in " f"{node}" ) - logger.critical( + self.logger.critical( f"Error describing thermal_conductivity flowchart: {e} in " f"{node}" ) @@ -260,7 +260,7 @@ def description_text(self, P=None, short=False): "Unexpected error describing thermal_conductivity flowchart: " f"{sys.exc_info()[0]} in {str(node)}" ) - logger.critical( + self.logger.critical( "Unexpected error describing thermal_conductivity flowchart: " f"{sys.exc_info()[0]} in {str(node)}" ) @@ -521,29 +521,50 @@ def analyze(self, indent="", P=None, style="full", run=None, **kwargs): start = max(fit0[i]["tau"]) * 3 else: start = 1 - slope, err, xs, ys = fit_helfand_moment( - self.M[i], ts, sigma=self.M_err[i], start=start - ) - fit.append( - { - "kappa": slope, - "stderr": err, - "xs": xs, - "ys": ys, - } - ) - v, e = fmt_err(slope, 2 * err) - if style == "1-line": - alpha = self.tensor_labels[i][0] - table["K" + alpha].append(v) - table["e" + alpha].append(e) - elif style == "full": - table["Method"].append("Helfand Moments" if i == 0 else "") - table["Dir"].append(self.tensor_labels[i][0]) - table["Kappa"].append(v) - table["±"].append("±") - table["95%"].append(e) - table["Units"].append("W/m/K" if i == 0 else "") + if start < 1: + start = 1 + + if self.logger.isEnabledFor(logging.DEBUG): + self.logger.debug("\n\n\n**********************\n") + self.logger.debug(f"{i=}") + for v1, v2, v3 in zip(ts[0:9], self.M[i][0:9], self.M_err[i][0:9]): + self.logger.debug(f" {v1:.3f} {v2:12.4e} {v3:12.4e}") + self.logger.debug("...") + for v1, v2, v3 in zip( + ts[-9:-1], self.M[i][-9:-1], self.M_err[i][-9:-1] + ): + self.logger.debug(f" {v1:.3f} {v2:12.4e} {v3:12.4e}") + self.logger.debug(f"{start=}") + self.logger.debug("--------") + self.logger.debug("") + + try: + slope, err, xs, ys = fit_helfand_moment( + self.M[i], ts, sigma=self.M_err[i], start=start, logger=self.logger + ) + fit.append( + { + "kappa": slope, + "stderr": err, + "xs": xs, + "ys": ys, + } + ) + v, e = fmt_err(slope, 2 * err) + if style == "1-line": + alpha = self.tensor_labels[i][0] + table["K" + alpha].append(v) + table["e" + alpha].append(e) + elif style == "full": + table["Method"].append("Helfand Moments" if i == 0 else "") + table["Dir"].append(self.tensor_labels[i][0]) + table["Kappa"].append(v) + table["±"].append("±") + table["95%"].append(e) + table["Units"].append("W/m/K" if i == 0 else "") + except Exception as e: + logger.warning("The fit of the Helfand moments failed. Continuing...") + logger.warning(e) self.plot_helfand_moment(self.M, ts, M_err=self.M_err, fit=fit)