V&V-Driven Development: The pyFrost-GM Workflow

Why Re-implement a Code That Already Works?

pyFrost-GM is a Python 3 reimplementation of LoKI-GM, a global model for low-temperature plasma chemistry developed at IST Lisbon and released as MATLAB code at ESCAMPIG 2026. It is licensed under GPLv3. The source code is provided to clients as part of Deep Why's consulting and simulation engineering services — get in touch to discuss your project.

LoKI-GM is current, well-documented, and actively maintained. The reimplementation exists to add capabilities that MATLAB cannot offer:

• No license lock — fully open-source (GPLv3), zero MATLAB dependency
• A second, independent implementation — forces every modelling choice to be understood explicitly, surfacing subtleties that are invisible from inside a single codebase
• Cloud and HPC readiness — containerised (devcontainer, GitHub Codespaces), ready for parametric sweeps on remote infrastructure
• Automated V&V pipelines — digitise reference data, run sweeps, generate comparison reports programmatically
• Full ownership and extensibility — add new chemistries, couple with ML, integrate into larger workflows

Verification: MATLAB Baselines as the Development Harness

The original MATLAB code served as the numerical ground truth. Each Python module was developed against parity tests that compared outputs to MATLAB baselines:

• Rate coefficients — verified bit-identical to MATLAB values
• Trajectory evolution — species densities, temperatures, and electric field tracked within bounded tolerance across the full time integration
• pytest harness — parity tests ran continuously during development, catching regressions immediately

A module was finished when it matched its baseline, and not before. Verification baselines were added incrementally as mismatches turned up during development, or as a new figure pulled in a part of the model that didn't exist yet.

Figure 1: Reduced electric field $E/N$ convergence — an example of a verification target. pyFrost-GM (markers) vs. LoKI-GM reference (dashed). Each such comparison served as a gate before proceeding to the next development phase.

What Parity Testing Surfaced

Parity testing surfaced a subtlety invisible from inside a single codebase: 3 of 4 O₂ vibrational rate functions (VT/VV processes) had onGasTemperature dependency flags that froze rates at the initial gas temperature $T_g$ instead of updating them as $T_g$ evolved self-consistently. Whether that was an oversight or a deliberate simplification for the standard test cases doesn't matter — a second, independent implementation made it visible immediately. The cause was traced, documented, and a compatibility flag now reproduces either behaviour.

Incremental Validation: Expanding Figure by Figure

Verification establishes code-to-code agreement. Validation asks a different question: does the physics model reproduce reality?

Validation wasn't a fixed, complete target set defined at the outset. Development was incremental: starting with a small number of published figures, building only enough of the model to attempt them, checking against the data, and then letting each new figure pull in whatever additional physics it needed. Every figure that already passed had to keep passing as the code expanded.

For the oxygen DC discharge studied by Dias et al. (2023) — a cylindrical chamber (radius $R = 1$ cm, discharge current $I = 30$ mA) swept across 12 pressures from 0.19 to 10 Torr — figures were added incrementally. This produced comparisons for:

• Reduced electric field $E/N$ and gas temperature $T_g$
• Electron density $n_e$ and atomic oxygen $O(^3P)$ density
• Singlet delta oxygen $O_2(a^1\Delta_g)$ density
• Vibrational distribution functions $O_2(X, v)$ for $v = 0$ to $41$
• Reaction pathway contributions for $O(^3P)$ creation and destruction

Figure 2a: Gas temperature $T_g$ vs. pressure — pyFrost-GM compared to experimental measurements from Dias et al. (2023).

Figure 2b: $O_2(X, v)$ vibrational distribution function at 0.2, 1, and 10 Torr — demonstrating VDF resolution across operating regimes.

Coverage: Input Files vs. Literature Corpus

In terms of chemistry coverage, the O₂ discharge represents 1 of 3 vibrational chemistry folders (O₂, N₂, and CO₂), with nitrogen and CO₂ next in the backlog. But counting folders understates the coverage, and it measures the wrong thing. The pressure sweep is a parametric study of 12 operating points, each producing multiple physical observables, reconstructed specifically to compare against published data.

Counting reproduced figures rather than input files gives the more meaningful measure. pyFrost-GM is currently validated against 2 publications, covering 7 reproduced figures (6 from Dias et al. 2023, plus 1 for electron density from Alves et al. 2026). Each new chemistry adds its own targets — and the gaps in the coverage map (like nitrogen, CO₂, and more oxygen cases) are the backlog for future work.

Development Timeline

The full reimplementation — from first commit to validated O₂ pressure sweep — was completed by a single developer in approximately 30 calendar days.

The C++ Bridge: Reuse, Don't Rewrite

The Boltzmann EEDF solver (LoKI-B-cpp, from a previous project in the LoKI-Suite) was integrated via an 85-line pybind11 bridge (lokib_bridge.cpp) rather than being reimplemented. SHA-256 input hashing eliminates redundant Boltzmann solves when input parameters haven't changed, yielding a ~28× speedup compared to the subprocess approach.

The AI-Agentic Workflow

AI coding assistants accelerated the translation — but this is not a story about AI writing physics code autonomously. It's about a physicist using AI as a force multiplier under rigorous verification constraints.

Where AI Helped

• Boilerplate translation: MATLAB-to-Python syntax conversion, class scaffolding, file I/O patterns
• Test scaffolding: Generating pytest structure, comparison plotting scripts, CSV data handling
• Documentation and packaging: Devcontainer configuration, README generation, license setup

Where Physics Expertise Was Essential

• Root-causing the $T_g$-freezing discrepancy: Understanding why vibrational rates diverged required knowledge of detailed balance and self-consistent gas temperature coupling
• Interpreting validation discrepancies: Knowing when a 10% deviation is physically acceptable versus a sign of a code error
• Architectural decisions: How to structure the solver coupling, when to cache EEDF results, how to handle chemical equilibrium convergence

Beyond Validation: Capabilities Unlocked

With the verified and validated Python architecture in place, new capabilities become immediately accessible — things that would be difficult or impossible in the original MATLAB framework:

• Massive parametric studies — automated pressure, current, and geometry sweeps with DVC-tracked outputs for full reproducibility
• New chemistries — N₂, CO₂, and mixed-gas systems can be added by extending the chemistry database, not rewriting solver code
• ML integration — Python-native architecture enables direct coupling with machine learning for surrogate modeling, sensitivity analysis, or reaction rate optimisation
• Automated reporting — from simulation to comparison plots to PDF carousel, the entire validation pipeline runs programmatically
• Containerised deployment — one-click setup via devcontainer or GitHub Codespaces; ready for cloud HPC

Note: While this article focuses on oxygen discharges, additional chemistries such as N₂, CO₂, and more are available upon consulting request.

Bibliography