GM_Julia

Published: Jul 14, 2025 by Peter Hill

GM_Julia is a zero-dimensional chemistry reaction solver for low-temperature plasmas. It takes a list of reactions, some details on the plasma species, and evolves them in time. These systems can be computationally intensive due to the number of reactions and their disparate timescales: typical chemistry sets in GM_Julia are around 600 reaction equations, with rate coefficients spanning many orders of magnitude. This makes for very numerically stiff systems, requiring specialised solvers. The original author made the decision to use Julia in order to take advantage of OrdinaryDiffEq.jl’s Rosenbrock integrators, designed for stiff systems.

In most ways, this was a fairly typical project: take a PhD code, and add packaging and tests, do some general refactoring. This was the first (and so far, only) Julia project we’ve taken on, so it was a bit of learning curve for me – not necessarily in terms of the language, which is fairly straightforward, but the language environment and how the idiomatic ways of using it. For example, almost everything in Julia expects you to use it from the REPL, which I found quite awkward in many respects. And while the package ecosystem is very well developed for science, the testing experience was very disappointing, with absolutely nothing coming to close to the might of pytest for Python. I also found the very slow compilation times combined with runtime type checking to be very frustrating. After working in Rust for quite some time this year, Julia is a bit of a let down!

Language gripes aside, this was a pretty straightforward project. After getting the code properly packaged and adding tests, I made three major changes. Firstly, I switched the input files from a custom implementation based on the EPOCH “deck” format to TOML files for most inputs parameters, but moving the reaction sets to a separate file with a custom format. Secondly, I enabled runtime evaluation of the rate equations. Previously, these required a separate “prerun” step to parse the equations from the input file, and write them as Julia to a new file which had to be included in the main source. This was a bit awkward to use and didn’t play well with Julia packaging, but luckily it was trivial to replace with eval (and making use of the RuntimeGeneratedFunctions.jl package to smooth over some mind boggling language issues). Of course, there are some security concerns with using eval, so this should only be used with trusted inputs. A more secure implementation might be to use a scripting language or custom parser with limited abilities. The last major change was to refactor how species were stored and referenced from a simple integer index in a vector, to a string name in a Dict. This made the code a lot more flexible, and will enable the researchers to expand and customise the chemistry sets a lot more easily.

This last point actually is a nice example of how the choice of data structure impacts the actual code. The initial design choice of a fixed-sized vector meant a separate structure was required to keep track of the mapping between species and its index into this vector. This made adding new species much more difficult, with several places in the code needing to be updated. A good question to ask yourself about any design is “how do I add one more X?”, where “X” is usually some kind of user choice – species, initial condition, solver, widget, and so on. Ideally, one more X requires touching only one place in the source code. Often when we find we need to touch multiple places, a change in design or data structure can help us simplify things.