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()

Step 4: Element Balance Verification

Critical for validating Gibbs reactor correctness:

def verify_element_balance(feed_composition, product_composition,
                            element_matrix):
    """Verify element conservation.

    element_matrix: dict mapping element -> {species: count}
    Example: {"H": {"H2S": 2, "H2O": 2, "H2": 2}, "S": {"H2S": 1, "S8": 8}}
    """
    results = {}
    for element, species_counts in element_matrix.items():
        feed_total = sum(feed_composition.get(sp, 0) * count
                         for sp, count in species_counts.items())
        prod_total = sum(product_composition.get(sp, 0) * count
                         for sp, count in species_counts.items())

        if feed_total > 0:
            rel_error = abs(feed_total - prod_total) / feed_total
        else:
            rel_error = 0.0 if prod_total == 0 else float("inf")

        results[element] = {
            "feed": feed_total,
            "product": prod_total,
            "relative_error": rel_error
        }
    return results

def plot_element_balance_closure(df, save_path):
    """Plot element balance errors across all cases."""
    fig, ax = plt.subplots(figsize=(8, 6))

    elements = [col for col in df.columns if col.startswith("elem_err_")]

    for col in elements:
        elem_name = col.replace("elem_err_", "").upper()
        errors = df[col].replace(0, 1e-16)  # avoid log(0)
        ax.semilogy(range(len(errors)), sorted(errors),
                     label=elem_name, linewidth=2)

    ax.axhline(y=1e-10, color="red", linestyle="--",
               alpha=0.5, label="Target (1e-10)")
    ax.set_xlabel("Case index (sorted)")
    ax.set_ylabel("Relative element balance error")
    ax.set_title("Element Balance Closure")
    ax.legend()
    ax.grid(True, alpha=0.3)

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

Step 5: Jacobian Condition Number Analysis

def plot_jacobian_conditioning(df, save_path):
    """Analyze Jacobian conditioning across conditions."""
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # Condition number vs temperature
    ax = axes[0]
    ax.scatter(df["T_K"] - 273.15, df["jacobian_cond_number"],
               c=df["iterations"], cmap="viridis", alpha=0.5, s=20)
    ax.set_xlabel("Temperature (°C)")
    ax.set_ylabel("log₁₀(Condition Number)")
    ax.set_title("Jacobian Conditioning vs Temperature")
    ax.grid(True, alpha=0.3)

    # Condition number vs iterations
    ax = axes[1]
    ax.scatter(df["jacobian_cond_number"], df["iterations"],
               alpha=0.3, s=20)
    ax.set_xlabel("log₁₀(Condition Number)")
    ax.set_ylabel("Iterations to Convergence")
    ax.set_title("Conditioning vs Convergence Speed")
    ax.grid(True, alpha=0.3)

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

Step 6: Adiabatic vs Isothermal Comparison

def plot_adiabatic_vs_isothermal(df_iso, df_adi, save_path):
    """Compare equilibrium outcomes between modes."""
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # Temperature change in adiabatic mode
    ax = axes[0]
    dT = df_adi["outlet_T_K"] - df_adi["T_K"]
    ax.scatter(df_adi["T_K"] - 273.15, dT, alpha=0.5, s=20)
    ax.axhline(y=0, color="black", linestyle="-", alpha=0.3)
    ax.set_xlabel("Feed Temperature (°C)")
    ax.set_ylabel("ΔT (K)")
    ax.set_title("Adiabatic Temperature Change")
    ax.grid(True, alpha=0.3)

    # Iteration comparison
    ax = axes[1]
    ax.scatter(df_iso["iterations"], df_adi["iterations"], alpha=0.3, s=20)
    max_iter = max(df_iso["iterations"].max(), df_adi["iterations"].max())
    ax.plot([0, max_iter], [0, max_iter], "k--", alpha=0.5)
    ax.set_xlabel("Isothermal Iterations")
    ax.set_ylabel("Adiabatic Iterations")
    ax.set_title("Iteration Cost: Adiabatic vs Isothermal")
    ax.grid(True, alpha=0.3)
    ax.set_aspect("equal")

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

Step 7: Trace Species Behavior

def plot_trace_species(df, species_list, save_path):
    """Show how trace species evolve with temperature on log scale."""
    fig, ax = plt.subplots(figsize=(10, 7))

    for species in species_list:
        col = f"n_{species}"
        if col in df.columns:
            vals = df[col].replace(0, np.nan)
            ax.semilogy(df["T_K"] - 273.15, vals,
                         label=species, linewidth=2, marker="o", markersize=3)

    ax.set_xlabel("Temperature (°C)")
    ax.set_ylabel("Equilibrium moles")
    ax.set_title("Trace Species at Chemical Equilibrium")
    ax.legend()
    ax.grid(True, alpha=0.3)

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

Step 8: Validation Against Reference Data

def plot_validation_parity(neqsim_results, reference_results,
                            species_name, save_path):
    """Parity plot: NeqSim vs reference (NASA CEA, JANAF, etc.)."""
    fig, ax = plt.subplots(figsize=(7, 7))

    ax.scatter(reference_results, neqsim_results, s=40, alpha=0.7,
               edgecolors="black", linewidths=0.5)

    # Diagonal
    lims = [min(min(reference_results), min(neqsim_results)),
            max(max(reference_results), max(neqsim_results))]
    ax.plot(lims, lims, "k-", alpha=0.5, label="Perfect agreement")

    # ±10% bands
    ax.plot(lims, [l * 1.1 for l in lims], "r--", alpha=0.3, label="±10%")
    ax.plot(lims, [l * 0.9 for l in lims], "r--", alpha=0.3)

    ax.set_xlabel(f"Reference {species_name}")
    ax.set_ylabel(f"Calculated {species_name}")
    ax.set_title(f"Validation: {species_name}")
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_aspect("equal")

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

Figure Catalog for Gibbs Reactor Papers

Table Catalog

Key Validation Criteria

For a Gibbs reactor paper to be credible:

1. Element balance: Relative error < 1e-10 for ALL elements in ALL cases
2. Gibbs energy: Total G strictly decreases (or unchanged) each iteration
3. Reference agreement: AAD < 5% vs NASA CEA for major species at equilibrium
4. Trace species: Non-negative mole numbers (no unphysical negative values)
5. Mass balance: |mass_in - mass_out| / mass_in < 1e-10