run_flash_experiments

Skill: Run Flash Experiments

Purpose

Execute NeqSim flash calculations in batch mode, collect metrics, and produce structured result files for paper-quality benchmarking.

When to Use

• Running a benchmark suite designed by the design_flash_benchmark skill
• Comparing baseline vs candidate algorithms
• Generating raw data for the validation agent

Execution Procedure

Step 1: Load Benchmark Config

import json

with open("benchmark_config.json") as f:
    config = json.load(f)

Step 2: Generate All Cases

import numpy as np
from itertools import product

def generate_all_cases(config):
    """Generate all benchmark cases from config."""
    cases = []
    case_id = 0

    for family in config["families"]:
        # Base composition
        base = family["base_composition"]

        # Composition variants
        names = list(base.keys())
        alpha = np.array([base[n] for n in names]) * family["dirichlet_concentration"]
        np.random.seed(42)  # Reproducible

        compositions = [base]  # Include base
        for _ in range(family["n_composition_variants"] - 1):
            x = np.random.dirichlet(alpha)
            compositions.append(dict(zip(names, x.tolist())))

        # PT grid
        T_vals = np.linspace(family["T_range_K"][0], family["T_range_K"][1], family["n_T"])
        P_vals = np.logspace(
            np.log10(family["P_range_bara"][0]),
            np.log10(family["P_range_bara"][1]),
            family["n_P"]
        )

        for comp in compositions:
            for T, P in product(T_vals, P_vals):
                cases.append({
                    "case_id": f"{family['name'][:2].upper()}-{case_id:05d}",
                    "family": family["name"],
                    "components": comp,
                    "T_K": float(T),
                    "P_bara": float(P)
                })
                case_id += 1

    return cases

Step 3: Run Single Flash Case

import time
from tools.neqsim_bootstrap import get_jneqsim
jneqsim = get_jneqsim()

SystemSrkEos = jneqsim.thermo.system.SystemSrkEos
SystemPrEos = jneqsim.thermo.system.SystemPrEos
ThermodynamicOperations = jneqsim.thermodynamicoperations.ThermodynamicOperations

EOS_MAP = {
    "SRK": SystemSrkEos,
    "PR": SystemPrEos,
}

def run_flash_case(case, eos_name="SRK", timing_repeats=3):
    """Run a single TPflash and return metrics."""
    EosClass = EOS_MAP[eos_name]

    # Create fluid system
    fluid = EosClass(case["T_K"], case["P_bara"])
    for comp_name, frac in case["components"].items():
        fluid.addComponent(comp_name, frac)
    fluid.setMixingRule("classic")

    ops = ThermodynamicOperations(fluid)

    # Warmup run
    try:
        ops.TPflash()
    except Exception:
        pass

    # Timed runs
    times_ns = []
    for _ in range(timing_repeats):
        # Reset and re-flash
        fluid2 = fluid.clone()
        ops2 = ThermodynamicOperations(fluid2)

        t0 = time.perf_counter_ns()
        try:
            ops2.TPflash()
            elapsed = time.perf_counter_ns() - t0
            times_ns.append(elapsed)

            fluid2.initProperties()
            n_phases = int(fluid2.getNumberOfPhases())
            beta_vapor = float(fluid2.getBeta(0)) if n_phases > 0 else 0.0
            converged = True
            error = None
        except Exception as e:
            elapsed = time.perf_counter_ns() - t0
            times_ns.append(elapsed)
            n_phases = -1
            beta_vapor = -1.0
            converged = False
            error = str(e)

    median_time_ms = float(np.median(times_ns)) / 1e6

    return {
        "case_id": case["case_id"],
        "family": case["family"],
        "eos": eos_name,
        "T_K": case["T_K"],
        "P_bara": case["P_bara"],
        "converged": converged,
        "cpu_time_ms": round(median_time_ms, 4),
        "n_phases": n_phases,
        "beta_vapor": round(beta_vapor, 6) if beta_vapor >= 0 else None,
        "error": error
    }