analyze-gibbs-convergence

Skill: Analyze Gibbs Convergence

Purpose

Interpret Gibbs energy minimization convergence metrics, analyze Jacobian conditioning, verify element balance closure, and produce publication-quality figures for chemical equilibrium papers.

When to Use

• After running Gibbs reactor benchmark experiments
• When analyzing convergence of Newton-Raphson Gibbs minimization
• When comparing solver variants (baseline vs optimized)
• When investigating failure cases in chemical equilibrium

Analysis Procedure

Step 1: Load and Parse Results

import json
import pandas as pd
import numpy as np

def load_reactor_results(results_dir, solver_name):
    """Load JSONL results into DataFrame."""
    records = []
    with open(f"{results_dir}/raw/{solver_name}_results.jsonl") as f:
        for line in f:
            records.append(json.loads(line))
    return pd.DataFrame(records)

Step 2: Equilibrium Composition vs Temperature

The most important figure for a chemical equilibrium paper:

import matplotlib.pyplot as plt

def plot_equilibrium_composition(df, system_name, save_path):
    """Plot equilibrium mole fractions vs temperature for all species."""
    fig, ax = plt.subplots(figsize=(10, 7))

    species = [col for col in df.columns if col.startswith("n_")]

    for species_col in species:
        name = species_col.replace("n_", "")
        ax.semilogy(df["T_K"] - 273.15, df[species_col],
                     label=name, linewidth=2)

    ax.set_xlabel("Temperature (°C)", fontsize=12)
    ax.set_ylabel("Equilibrium mole fraction", fontsize=12)
    ax.set_title(f"Chemical Equilibrium — {system_name}", fontsize=14)
    ax.legend(loc="best", fontsize=10)
    ax.grid(True, alpha=0.3)
    ax.set_ylim(bottom=1e-12)

    plt.tight_layout()
    plt.savefig(save_path, dpi=300, bbox_inches="tight")
    plt.close()

Step 3: Convergence Iteration Analysis

def plot_iteration_heatmap(df, save_path):
    """Heatmap of iteration count in T-P space."""
    fig, ax = plt.subplots(figsize=(10, 7))

    pivot = df.pivot_table(values="iterations", index="P_bara",
                           columns="T_K", aggfunc="mean")

    im = ax.pcolormesh(pivot.columns - 273.15, pivot.index,
                        pivot.values, cmap="YlOrRd", shading="auto")
    plt.colorbar(im, ax=ax, label="Iterations")

    ax.set_xlabel("Temperature (°C)")
    ax.set_ylabel("Pressure (bara)")
    ax.set_title("Newton Iteration Count")
    ax.set_yscale("log")

    plt.tight_layout()
    plt.savefig(save_path, dpi=300, bbox_inches="tight")
    plt.close()