The Design of a Python Library for the Automatic Definition and Simulation of Transient Ionization Fronts

In recent years, the interest in nonthermal plasma dynamics has grown significantly, within both industry and research. This has been driven by the development of several novel cold plasma technologies across a wide range of different fields, for example, for plasma medicine, chemical processing, pollution control, and surface treatment. The optimization of these technologies relies heavily upon the understanding of gas discharge plasmas: their generation, electrical characteristics, and interaction with their surroundings. Moreover, the manifestation of nonthermal plasmas in the form of streamers is of high relevance and critical importance to high voltage insulation technology, and has further significance to geophysical research concerning atmospheric discharges. The present work describes the development of the StrAFE (Streamers on Adaptive Finite Elements) package, a dedicated Python library built atop the popular open-source FEniCS finite element software, designed with the objective to simplify and to automate the solution of ionization front models. The library features support for mesh adaptivity, distributed memory parallelism, and an intuitive programming interface, while providing an exceptionally high level of user configurability. This article presents the software implementation, describes its features, and presents several code verification studies performed within simple and complex domains. It is concluded that the numerical results gained from this open-source framework are comparable to other well-established software in terms of accuracy. Therefore, it further demonstrates the great potential for open-source software to make significant contributions to technologies involving nonthermal plasmas, ionization fronts, and gas discharges.

The associate editor coordinating the review of this manuscript and approving it for publication was Agustin Leobardo Herrera-May . Deviation parameter for gaussian seed.

I. INTRODUCTION
Nonthermal plasmas belong to a class of non-equilibrium (non-Maxwellian) gas discharge plasma, where the electron temperature is far higher than that of the ion temperature. Such plasmas may be produced under the application of sufficiently intense electric fields to gaseous media. Also known as cold-or low-temperature plasmas, their properties have led to their application to various novel technologies, found across many different areas of engineering and science. For example, for surface processing [1], [2], air cleaning [3], [4], [5], plasma medicine [6], or chemical processing [7]. However, many aspects of nonthermal plasmas remains to be fully understood, and a deeper understanding of the related processes is crucial for the successful advancement and optimization of these developing technologies. Equally, the study of fast ionization fronts is also critical to the prevention of electrical failure in high-voltage equipment, as the initial development of streamers in high-field regions may result in partial (or complete) discharges across insulation, compromising system integrity. Streamer discharges are similarly important for the understanding of geophysical phenomena occurring in the upper atmosphere [8], and for certain plasma flow and propulsion applications [9]. These transient, filamentary-type discharges are multiscale in nature [10], rendering them highly difficult to characterize solely through experimental means. In recent times, advances in computational power, and the growing accessibility to high performance computing (HPC), has facilitated complex multiphysics modelling, allowing it to become more widely available. For example, time-dependent simulations of streamer discharges have been successfully demonstrated on standard consumer grade hardware, using commercially available software [11], [12]. Specialized schemes using custom codes, employing techniques which are particularly effective for multiscale modelling have also been developed, and demonstrate major advances in computational efficiency and speed [13], [14]. The popularity of these methods has grown rapidly, with models that incorporate aspects such as complex plasma chemistry [15], nonuniform electric fields [16], and solid dielectric boundaries [17], [18], becoming increasingly prevalent. Currently, the software developed must typically strike a balance between issues such as computational speed, resource usage, accessibility, adaptability, accuracy, and model fidelity. Further difficulty comes with the sheer number of necessary parameters involved within plasma simulations, much of which must be manually entered into simulation software, a process which is highly susceptible to human error [19]. Generally, custom modifications to complex, often low-level, code are required when setting up new problems, which can be time consuming, and, in many cases, require some degree of computational expertise.
In this work, the development and design of StrAFE (Streamers on Adaptive Finite Elements), a convenient and user-friendly Python library, for the automation and solution of transient ionization front problems, is described. The StrAFE system was built atop the open-source FEniCS finite element (FE) framework [20], [21], [22], which provides an already convenient yet highly flexible set of software tools for general-purpose finite element modelling. The design philosophy behind StrAFE follows directly in the footsteps of FEniCS: the rapid translation of physical models to efficient FE code; reduction in the setup time necessary to generate new models, by reducing the need for computational expertise; and to allow maximum user control and configurability.
StrAFE was developed using the now legacy (2019.1.0) version of FEniCS. Development is now focused on a next-generation version known as FEniCSx, and documentation for both versions of this platform may be found at docs.fenicsproject.org. In this article, the technical implementation of StrAFE using FEniCS is described, while also discussing the many benefits which can be gained by employing open-source code. Results from several code verification and comparison studies are then presented. It should be noted that the StrAFE library itself has not been made openly available. Rather, the focus of this article is to demonstrate how software like StrAFE can be achieved using FEniCS, which itself is open-source and fully available. Consequently, this article acts as a guide for interested parties to develop similar software, and encourages the use of open-source frameworks such as FEniCS for simulating transient ionization fronts developed in gases.

A. THE FINITE ELEMENT METHOD
When using the hydrodynamic description of plasmas, solving the set of hyperbolic conservation laws which arise (see section III-A) generally benefits from numerical methods that guarantee exact local conservation, e.g., finite volume methods (FVM) as used in [13], [14], and [23]. However, this work is instead based on the continuous Galerkin finite element method (CG-FEM), to take advantage of some FEMspecific features. Despite only maintaining global conservation, FEM offers greater flexibility in terms of the order of spatial discretization, the type of finite element, and can handle complex boundaries with relative ease, without the need for any specialized numerical schemes. CG-FEM is mature, versatile, and has been frequently applied in many areas of science and engineering. Thus, a detailed review of FEM will not be conducted here, for which the reader is referred to [21], [24], and [25], and references therein. Instead, only a brief outline of necessary aspects is provided, to act as support for later discussion within this article.
To begin, let u denote the exact solution of an unknown function, and u h be its approximation. FEM is defined when u h is sought as a linear combination of basis functions, φ i , such that: where U i are coefficients to φ i which approximate u. Suppose that u satisfies the differential equation: in a bounded domain , where L is a differential operator; then the weighted residual reads: for weight function (or test function) v. In the standard Galerkin method, weight functions v are chosen to belong to the same function space as the approximating basis functions φ i , then the Galerkin FEM problem is solved by requiring that the weighted residual be driven to zero, i.e., where the integration over is over each finite element which spatially discretizes the domain. Application of (4) over each element results in the linear system: where A is a matrix, b is a vector, and U is the vector of unknown coefficients to be obtained. The system (5) can then be solved using any preferred linear algebra solution method. In physical problems, the differential operator L often has order two, which due to (4) imposes the condition that φ i must be at least twice-differentiable (i.e., elements with a linear basis cannot be used, which are preferred for their lower computational cost). To lift this constraint, Green's identity can be applied to the integral (4), which reduces second-order derivative terms to two derivatives of order one, yielding the weak or variational formulation of (4). This is illustrated further in section III-G for the specific case of coupled driftdiffusion-reaction-Poisson equations.

B. THE FEniCS PROJECT
The FEniCS Project [20], [21], [22] is an open-source collection of software components, designed with the aim to automate and simplify general FE problem generation and solution. Licensed under LPGL-v3, the original platform was developed using C++, and features several core components under the main user interface named DOLFIN. Fig. 1 outlines the corresponding system architecture; using corresponding labels: (i) FEniCS Form Compiler (FFC) [21], [26], [27] -Handles just-in-time (JIT) compilation and generation of  [20], [28] -A language developed for the discretization of partial-differential equations (PDE), in a form that closely resembles the original mathematical expressions. (iv) Finite Element Automatic Tabulator (FIAT) [29] -Handles the generation of finite element basis functions and elements of arbitrary order. The components listed above work alongside the system assembler (v) to form the system matrix, to then be passed to the solver. FEniCS can additionally be programmed using Python, which significantly lowers the entry barrier to performing efficient FEM simulations. Furthermore, there is native support for distributed memory parallelization through MPI (vi), with built-in mesh partitioning and load balancing, performed using the PT-SCOTCH [30] or ParMETIS [31] libraries (vii). FEniCS can also be compiled with a selection of popular open-source linear algebra backends (viii), which provides efficient computation, e.g., PETSc [32], [33], uBLAS [34], or Epetra [35]. Direct programming through Python also facilitates the automated processing of data outputs, since scripts for data analysis can be constructed using any available third-party Python library (ix), and can directly communicate with FEniCS. The platform has proved highly capable in numerous studies across a range of disciplines [36], [37], [38], and has shown to adapt well to HPC settings with over 24,000 processes [22], [39]. StrAFE is primarily an automation tool built on top of FEniCS through its Python interface, but additionally implements crucial features such as adaptive mesh refinement for efficient multiscale modelling.

III. MATHEMATICAL MODELLING A. HYDRODYNAMIC DESCRIPTION OF NONTHERMAL PLASMAS
The hydrodynamic (or fluid) approximation of plasma arises from the zero-order moment of the Boltzmann equation [40], and becomes valid when charge concentrations are sufficiently high, or if the Knudsen number is within the bounds of validity [40]. Consider a plasma comprised of a set, N , of different charged species, interacting under the presence of an electric field. According to the hydrodynamic approximation, the spatial and temporal evolution of their concentrations can be modelled following a set of advection-diffusion-reaction equations, given by: for i ∈ N , where the term S i describes the sources and sinks of species i, and the flux ⃗ i is given by the sum of advective and diffusive components, where q i , n i , µ i , and D i are the electric charge, concentration, (positive) mobility, and diffusion coefficient of species i, respectively. ϕ is the potential, coupled to the charge concentrations through the Poisson equation, where ε is the permittivity. The self-consistent solution of (6) and (8) with appropriate boundary conditions, initial conditions, and charge sources, forms the basic hydrodynamic description of plasmas. The types of initial and boundary conditions which are supported in StrAFE are described in the following sections, as are several additional functions, which have been developed to improve modelling fidelity. Sections III-A to III-G firstly focuses on the mathematical formulation, while details regarding the usage and syntactical implementation of the software can be found in section IV.

B. USING TOWNSEND COEFFICIENTS
A basic description of plasma can be gained using a minimum of two species -electrons, and generic positive ions. Using a minimal model, the Townsend ionisation coefficient can be used to define an ionisation source term of the form: whereᾱ is the effective ionisation coefficient, and the subscript e refers to electrons. Alternatively,ᾱ can be split intō α = α − η, where α, η are the ionization and attachment coefficients, respectively, and the set N includes negative ions, with the source term: For advanced models using individual reaction rate coefficients, and considers specific plasma chemistry, section III-C applies.

C. PLASMA CHEMISTRY
For detailed studies which concern plasma composition, and which involve numerous active ionic species, source terms are computed following: where k j is the reaction rate for reaction j, n is the total number of reactions which involve species i, and R is the set of reactants involved in reaction j. Symbol h j is either +1 or -1 for sources and sinks, respectively. Examples which employ the plasma chemistry module of StrAFE have been demonstrated in section V.

D. PHOTOIONIZATION
In some gases (e.g., atmospheric air), photoionization induced by emitted photons from radiative de-excitation processes, are thought to be a major contributor of free electrons to sustained discharges such as positive streamers [41], [42], [43]. To account for this, Zheleznyak's model [41] has been implemented with a photoelectron source term given by: The photoelectron source can be efficiently obtained by using the Helmholtz approximation as described by Bourdon et al. [44]: for j = 1, 2, 3. Symbols p, p q , and p O 2 are the total gas pressure, collisional quenching pressure of nitrogen (accounting for non-radiative de-excitation), and the partial pressure of oxygen, respectively. ν u is the impact ionisation frequency for level u, ν i is the ionisation frequency, ξ is the photoionization efficiency [44], and S ion is the summed source term over all contributing ionizing reactions. Lastly, A j and λ j are fitting parameters for the j-th equation, and the values which are used throughout the examples within section V follow those computed in [44].

E. LOCAL FIELD AND LOCAL MEAN ENERGY APPROXIMATIONS
In many previous studies, the local field approximation (LFA) has been used for simplicity, which assumes that transport and rate coefficients are dependent only on the local magnitude of the electric field, which becomes immediately available from the solution of (8), and then applying ⃗ E = −∇ϕ. The range of validity of the LFA has been shown to be limited, and may become inaccurate for high-or spatially-inhomogeneous field conditions [45], [46]. In StrAFE, the LFA is the default option, however, there also exist the option to use the local mean energy approximation (LMEA), which attempts to include (to some extent) non-local electron transport processes, thereby extending the range of validity to beyond that of the LFA. Under the LMEA, transport and rate coefficients are dependent on the local value of the electron energy, which is dynamically computed by appending an additional balance equation to the model, given by: where the electron density flux ⃗ ε is Here, n ε is the electron energy density,ē is the elementary charge, ⃗ e is the electron flux, and E j is the energy loss or gain during reaction j. µ ε and D ε are the energy mobility and energy diffusion coefficient, respectively; these can either be user defined, or optionally, be approximated from the electron transport parameters [40] as:

F. BOUNDARY CONDITIONS
For electrode boundaries, Dirichlet boundary conditions are prescribed for the electrostatic potential, while Neumann-zero conditions are applied at axes of symmetry.
For studies involving subdomains (e.g., solid barriers, electrodes, walls), the accurate reflection of physical behavior requires appropriate boundary conditions for the normal charge fluxes to be prescribed at interfaces. In StrAFE, the option exists to have either a Neumann-zero condition, or wall conditions following Hagelaar et al. [47]. Using similar notation to Jovanovic et al. [19], these are given by: for electrons, for heavy species, and for the electron energy density. The parameter r is the reflection coefficient of the species at the boundary, γ is the secondary emission coefficient, andε γ is the mean energy of VOLUME 11, 2023 secondary electrons. v th is the thermal velocity of the species identified by the subscript, given by: where k B is the Boltzmann constant, T i is the species temperature, and m i is the species mass. For the electron temperature, T e is calculated from the energy following: and for the electron energy velocity, equation (22) holds: Across interfacial boundaries with differing permittivity, a discontinuity in the electric displacement is induced by the accumulation of surface charge, σ s , from incoming fluxes (volume conduction is currently not considered). To incorporate this, internal subdomain boundaries require that: where the surface charge is computed based on the sum of boundary-directed fluxes, following: G. WEAK FORMULATION As described in section II-A, the application of the Galerkin method requires the reformulation of equations (6)- (24) to hold weakly with respect to an appropriate set of test functions. The weak formulation of the governing equations is fundamental to the implementation of StrAFE, which uses the UFL syntax to directly input math-like form expressions, which are passed to the system assembler. See [20] and [21] for a detailed description of UFL and of the FEniCS form compiler. In the interests of repeatability, the remainder of this section focuses on providing the mathematical definition of the variational problem. To maintain consistency with indexing in Python, u is hereby defined as a vector function with components u = u 0 , u 1 , u 2 , . . . , u p , with index starting at zero, each of which represents a scalar field in a domain , and where p varies depending on the problem type (i.e., the number of species considered, whether LFA or LMEA is used, etc.). The rules for the construction of this vector are as follows: • Index 0 is always the potential field, ϕ. • Index 1 is always the electron density, n e . • Index 2 to N stores the densities of all other species under consideration, e.g., (n + , n − , . . .).
• Index N + 1 to N + 3 stores the three components of the photoionization source term, S ph,1 , S ph,2 , and S ph,3 , if enabled.
• Index N + 4 stores the surface charge σ s , defined only on element facets, and only if there are subdomains.
• The last index of the vector (in Python, index −1) stores the electron energy density, n ε , if LMEA is enabled. In summary, the mapping of vector u can be expressed as: . . , v p , is stored for the construction of variational forms, with each function v j acting as the test function for the corresponding index in u.
Consider firstly the Poisson equation (8), from which one may construct the weak form following section II-A to be: where ∂ denotes the external boundary of domain , arising from Green's identity. The last integral includes the term ε ∇u 0 ·n which may be chosen to enforce Neumann boundary conditions such as (23), or, if a Neumann-zero condition is present, this integral vanishes. Different Neumann conditions can also be prescribed to disjoint sections of ∂ , simply by splitting this integral into boundary segments such that: where ω is a boundary segment, and separate integrations are performed over each segment j, respectively. For advectiondiffusion-reaction equations (6), the weak form becomes: for i ∈ N , and similarly for the energy balance equation (14), but with the right-hand-side replacing S i in (28). The Helmholtz weak form for photoionization (13) is obtained in the same way, yielding: for i = N + 1, N + 2, N + 3.

IV. StrAFE
Fig . 2 shows the architecture of the StrAFE library, and describes its interconnections to the core FEniCS system as laid out in [48] and in Fig. 1. StrAFE was developed with a strong emphasis on usability and flexibility, in-line with the original design philosophy for FEniCS. Classes provided through StrAFE have been designed to be intuitive, and to require only basic knowledge of the Python programming language. Sections IV-A to IV-E gives a broad overview of the components of StrAFE, providing details regarding their implementation and usage.

A. BASIC USAGE AND CLASSES
In this section, classes which are central to StrAFE's functionality are described, to highlight how the combination of FEniCS and Python has been used to significantly simplify the process of configuring simulation models. At the core of any simulation is the computational mesh. The inheritance of StrAFE from FEniCS therefore permits that any mesh format supported by FEniCS as listed in [21], is also supported by StrAFE. This includes internally gen-erated meshes using native FEniCS functions, and external meshes from any other software. Both 2D and 3D unstructured meshes are supported, and are compatible with mesh routines as described in section IV-B. A simulation can then be configured using a number of dedicated StrAFE classes, the most important of which are described in brief below: • DriftDiffusionProblem() -The main problem class, where all simulation settings can be accessed and changed. Uses intuitive .set syntax, such as set_base_mesh(mesh) to provide a meshed domain object mesh to the solver. All settings pertaining to the problem can be set similarly, including the choice of LFA or LMEA, the species under consideration, AMR schemes, time-stepping, etc.
• ChargedSpecies() -Defines a single instance of a charged species to be tracked within the simulation. This class stores the species properties, such as its name, mass, charge, and transport parameters (can be provided as Python functions or tabulated data, such as those computed using BOLSIG+ [49]). Initial conditions are also stored in this class, where several popular initial charge distributions can be chosen from, e.g., a gaussian [50], capsule [17], or uniform distribution.
• Electrode() -Defines nodes on the mesh which are to be marked as electrodes. Doing so allows the application of potential boundary conditions. This is also a child class of the more general DirichletCondition(), which can be used for the prescription of arbitrary Dirichlet conditions, and is not limited to the potential field.
• Wall() -Defines nodes on the mesh which are to be marked as a wall. Doing so allows the application of the wall conditions according to (17)- (19). This class stores the attributes of the wall, such as reflection and secondary emission properties. This is also a child class of the more general NeumannCondition(), which can be used for the prescription of arbitrary Neumann conditions, and is not limited to the charge fluxes.
A subsequent call to the .solve() method launches the system assembler, which automatically constructs the UFL forms of equations (26)- (29), based on the provided settings. This is subsequently passed to the solver environment, and the time-stepping loop is invoked. In this way, simulation problems can be configured using minimal lines of highlyreadable code. This facilitates the rapid development and refinement of computational models, and has strong support for integration with external, Python-developed code. This is particularly useful for performing batch parametric studies, or to conduct automatic post-processing tasks on the simulation outputs using other Python-based programs.

B. ADAPTIVE MESH REFINEMENT (AMR)
A part of the motivation to develop StrAFE, was the want of a flexible and transparent framework to study fast-transient ionisation waves, such as streamer discharges, with emphasis VOLUME 11, 2023 on divergent field and time-varying field conditions. As mentioned, streamer discharges are multiscale phenomena [10], and often exhibit sharp features and steep spatial gradients, particularly within their electric field and ionic density distributions. This imparts considerable computational difficulty for mesh-based methods like FEM, where the accuracy is highly dependent on mesh resolution, but where computational resources are often limited. Therefore, adaptive meshing techniques have essentially become a necessity for simulating streamers in complex domains, or for longer time periods. These methods automatically adapt the mesh elements in time and space to suit the evolving dynamics of the simulation. Most commonly, this requires refinement of mesh cells in regions of fast-moving or high-gradient behavior, and subsequent de-refinement of cells where the solution is once again slow-varying (also known as h-type adaptivity, as opposed to p-type adaptivity which locally adapts the order of the approximating basis functions). The criterion for refinement is typically based upon an error estimate, or alternatively, can simply be based upon some combination of cell or nodal values of the most recently computed solution.
At the time of development of StrAFE, FEniCS did not include native support for time-dependent adaptive meshing schemes that allowed for per-cell de-refinement, and was limited to refinement only. To include this feature, StrAFE implements several refinement routines for the purpose of achieving dynamic meshes in time and in space. This is done through periodic domain remeshing at a user-defined frequency. Fig. 3 shows a flowchart of the AMR algorithm which has been employed, and which is also described in further detail by the following algorithm: 1) At the beginning of the simulation, store a base coarse mesh M b , and a fine mesh M 0 to solve for initial conditions. 2) If it is the first iteration, project initial conditions onto M 0 , solve once for u 0 , and store a temporary copy u r ← u 0 . If not, assign the current solution u r ← u k . In either case, set a temporary mesh M t ← M b . 3) Construct a set of functions F (or combination of functions) from those stored in u r , which are chosen by the user, to be evaluated against the refinement criteria.

4) For all functions f ∈ F, mark cells of M t for refinement
if f ≥ κ ℓ , where κ ℓ is a tolerance defined by a set of tolerances in the simulation settings, κ ℓ ∈ K. 5) (Optional) grow the region of refinement markers by marking all cells within a distance R from all marked cells in M t . 6) Refine M t based on these markers g ℓ times, then reassign M t ← M t . 7) For the number of total required refinement levels, ℓ, increment ℓ and repeat from step 4. 8) Return the fully refined mesh M t . 9) If it is the first iteration, re-project initial conditions onto M t and re-solve before continuing. 10) Project the current solution onto M t and re-assemble the variational problem. The refinement functions F, tolerances K, refinement multipliers g ℓ , refinement levels ℓ, and radial growth distance R, can all be user-defined as Python functions and passed to the problem environment, using similar .set syntax as described in section IV-A. The AMR routine also supports time as an input argument, which allows the refinement criteria to be adjusted automatically over the course of the simulation. Moreover, marker functions can be provided, which allows the refinement region to be limited to a user-specified sub-domain of the mesh. Additional routines which permit the static refinement of the base mesh M b around a point, a line, or a user provided region have also been developed. These are particularly useful for simulations involving static features, like dielectric boundaries or particles, which may require a denser mesh near its interfacial region, to ensure good numerical convergence. AMR would then refine on top of the already refined base mesh. Using this scheme, granular control of the mesh is provided to the user, as there exists many different methods to refine the mesh precisely where required over the course of the simulation. It is remarked that the AMR algorithm operates in parallel with MPI, and has been designed to minimize inter-process communication where possible. An example of an adaptively refined mesh using this scheme is shown in Appendix A.

C. TIME INTEGRATION
For time-stepping, StrAFE currently employs the θ-scheme, such that the time derivative of the function u as featured in (25)-(29) is discretized as: Setting the parameter θ = 0 or θ = 1 results in the first-order explicit and implicit Euler schemes, respectively, while the default θ = 1/2 recovers the Crank-Nicolson scheme, giving second order accuracy in time. Currently, StrAFE also has some initial support for adaptive time stepping, which when enabled, returns an estimate of the local truncation error, from which the user may implement custom timestep controllers, to update t for the next step based on their own specified tolerances. In future, the implementation of higher order schemes would be highly beneficial, to increase solution accuracy and numerical convergence.

D. SUBDOMAIN SUPPORT
In practice, many engineering systems are composite in nature, such that interfaces between materials (e.g., between gas and solid dielectrics) may exist inside the domain of interest, and may not necessarily be positioned at the external boundaries. In such cases, equations (23), (17)-(19) must be applied for internal nodes of the computational mesh. StrAFE currently supports subdomains for modelling solid dielectric material with different relative permittivity values.
To create subdomains in StrAFE, the user must first ensure that the attached mesh has nodes that align with the internal boundaries, which can be achieved with ease using available mesh generation software, e.g., gmsh [51]. Boundary and subdomain markers should also be generated from external software, unless the boundary interface conforms to simple shapes which can be described analytically, in which case, they can be marked using native FEn-iCS tools [21]. A subdomain is defined by passing the cell marker function, the boundary marker function, and the value of permittivity, reflection, secondary emission, etc., to the .set_solid_subdomain() method. If this subdomain boundary function is used, conditions given by (17)- (19) are automatically applied, and do not need to be manually provided by setting Neumann conditions. If not, the boundary defaults to a Neumann-zero condition.
An arbitrary number of solid subdomains can be defined, which are internally stored as a list of markers. These are automatically regenerated on each new mesh if AMR is enabled. The markers instruct the assembler to split the meshed domain into two sets, = s ∪ g , for solid and gas, respectively, where the drift-diffusion equations (6) are solved in g only, while Poisson's equation (8) is defined over the full domain . It is remarked that, although StrAFE currently does not provide a convenient interface to define additional physics within subdomains, this could theoretically be done in future, by simply appending the corresponding variational forms before passing to the solver. This may be useful for studies involving simultaneous charge transport in solid media, or subdomains which contain liquid dielectrics and fluid flow.

E. PROBLEM GENERATION AND PLASMA CHEMISTRY INPUT
In section IV-A, the ChargedSpecies() class was described, which is used to define and include a single charged species. However, in many nonthermal plasma applications, e.g., [7] and [15], knowledge of the plasma composition, and of the exact combination of charged species developed during a discharge, is of critical importance. Manual definition of each individual species using ChargedSpecies() would be tedious and error-prone. For this, StrAFE features automated generation of the species list from a plaintext file, containing a list of ionic species, neutral species, table of chemical reactions, rate coefficients, and energy loss coefficients. Using the TabulatedSpecies() class, this user provided file can be read, and the list of species is automatically generated and attached to the problem environment. Parameters can be provided within the text file as constants, plaintext functions, or a file path pointing toward tabulated data. The chemistry data is independent of the mesh data, therefore allowing studies in different gases to be achieved with ease, simply by switching to a different chemistry file.

F. SOLVER AND OTHER MISCELLANEOUS OPTIONS
On default settings, StrAFE solves the linear system using the solvers provided by the PETSc [33] backend through FEniCS. Newton's method is used for the outer iteration, whilst the generalized minimum residual method (GMRES), preconditioned using the Hypre [52] algebraic multigrid (AMG) method is used as the inner linear solver. However, the solver interface is fully exposed, and can be configured to meet the user's requirements. Other options include biconjugate gradient methods, various relaxation methods, and popular direct solvers such as MUMPS [53]. For a full list of supported linear algebra solvers and settings, the reader is referred to [21].
The versatility which StrAFE has inherited from FEniCS further enables other options to be readily changed, for example, the use of higher-order spatial discretization, or different element types. Moreover, it is possible to override or modify the variational form, which is assembled only when .solve() is called. This allows the implementation of additional physics, stabilization schemes, or other, to suit the exact application.

V. SIMULATION EXAMPLES
Since StrAFE relies upon FEniCS for system assembly and for its solver routines, it is unnecessary for this article to concern itself with the detailed benchmarking and validation of the numerical assembly and solver components. Instead, the objective is to demonstrate the accuracy and capabilities of the StrAFE interface, specifically for streamer and plasma modelling, which has been subsequently developed. At the time of writing, there exists no standardized 'benchmark' problems for nonthermal plasma simulations, though this may change as the interest in such numerical codes increase, with studies like [50]. Instead, the present section demonstrates several examples of streamer discharge simulations, in multiple configurations which have been reported in past literature, and were conducted by multiple other research groups. For each case, the simulation domain, boundary conditions, and relevant model settings have been VOLUME 11, 2023 FIGURE 4. Time evolution of an axisymmetric positive streamers' electric field and electron density, at times t = 3, 9, and 15 ns, generated using StrAFE. Original study from [50]. Equipotential lines are spaced by 2 kV. This simulation completed in approximately 3-4 hours.
briefly provided before the results are presented. Where possible, direct comparisons to original studies have been made, though this was not always possible, as it was dependent on data availability. Unless otherwise specified, all simulations were performed using triangular linear Lagrange elements with StrAFE running inside a Docker container [54], on either 16-core (AMD Ryzen 9 5950X) or 18-core (Intel Xeon W-2295) workstation computers with 64GB memory. It is remarked that while all the necessary components for 3D simulations have been implemented, only 2D (or 2Daxisymmetric) simulations have thus far been conducted. This is primarily due to the need for code verification and comparison, for which there is a far greater number of documented 2D studies. Since full 3D simulations have only become possible recently [13], [14], it would be difficult to evaluate the accuracy of the implementation based on limited data. However, it would be of high interest to conduct 3D simulations in the future, especially to evaluate StrAFE under HPC settings.

A. AXISYMMETRIC POSITIVE STREAMER
For the hydrodynamic approximation used in StrAFE, [50] provides the most recent and comprehensive comparison of different codes used by different groups for the simulation of streamer discharges. At the time of writing, it is essentially the only study which has been conducted purely for the purposes of verification and comparison. As such, evaluating results attained using StrAFE against the multiple available datasets from [50] was prioritized.
Focus was placed on the 'case 3' simulation of [50], to allow for the implementation of the Helmholtz photoionization model to also be evaluated. The domain consisted of a 2D square box with dimensions (r, z) = [1.25, 1.25] cm, and was rotationally symmetric around r = 0. Dirichlet conditions of U 0 = 18.75 kV and 0 kV were prescribed at z = 1.25 cm and z = 0 cm, respectively, and Neumann-zero conditions were present on all four sides for all charge fluxes and photoionization source terms. Only electrons and generic positive ions were considered, and transport parameters were supplied as empirically fitted expressions given in [50], using the local field approximation. Initial conditions consisted of a gaussian-distributed positive charge density placed at (r 0 , z 0 ) = (0, 1.0) cm following the expression: where N 0 = 5 × 10 18 m -3 and s = 0.4 mm. A uniform background density for both electrons and positive ions was set to 10 9 m -3 , and photoionization was enabled with the three-term exponential fitting parameters as provided in [50]. An adaptive timestep between 1 and 3 ps was used, with AMR enabled to perform remeshing every 30 iterations based on the value of the electric field normalized by the maximum field, and of the magnitude of the net space charge. The use of 30 iterations was determined from initial trial-and-error testing, and was found to provide a good balance between solution accuracy and total required computational time (frequent remeshing can be become detrimental to computational speed, but may not necessarily increase the solution accuracy). There currently exist no automated method to determine the optimal number of iterations between AMR, though this would be of interest to develop in the future. The panels of Fig. 4 are split down the axis of symmetry, where the left and right color plots correspond to the electric field magnitude, and electron density, respectively. Three timesteps of the positive streamer evolution have been shown, at 3, 9, and 15 ns. The streamer length over time, and the maximum field strength over the length of the streamer, has been additionally compared in Fig. 5 to the data from five other groups who participated in [50]. Results from StrAFE compare well to all other codes, and is well within the expected margin of error considering the many potential differences in implementation [50]. The results from StrAFE were found to resemble most closely those from group DE, whose computation was performed using the commercially available COMSOL Multiphysics [19] software. With a run-time of around three to four hours for the present code, and considering imperfect parallel scaling (see Appendix B), StrAFE compares well with this widely-used commercial option.

B. DOUBLE-HEADED COUNTERPROPAGATING STREAMERS
Demonstrated in several studies, e.g., [55], [56], and [57], is the initiation of streamers from an initially charged seed, which evolves into one positive and one negative streamer propagating in opposite directions, away from the location of the source. The morphology of negative streamers is known to differ significantly from their positive counterparts, as the differences in their sources of electrons have been suggested to result in vastly different characteristics, such as the maximum electric field, streamer radius, and propagation velocity [55]. When both are simultaneously present in a single simulation domain, it can be challenging to find a single set of AMR criteria which is equally effective for both streamers, because in general, the refinement criteria and tolerances require to be adjusted on a case-by-case basis. In these simulations, StrAFE is used with a partitioned-AMR scheme such that two different refinement criteria are used for the top half (z ≥ 7.5 mm), and bottom half (z < 7.5 mm), of the domain, ensuring that sufficiently fine meshes are provided for both streamers. This demonstrates a useful feature which can significantly aid convergence for some simulations, where there exists some prior knowledge of how the discharge will evolve in space. To avoid over-refinement in the early stages of the simulation, the AMR scheme was also configured to be time-dependent. The usage of E/E max is effective as a refinement threshold, but not during the initial stages where E/E max ≈ 1 everywhere in the domain, due to the presence of a neutral seed. Therefore, this criterion is only introduced once t ≥ 4 ns.  The domain in this case was once again axisymmetric, forming a cylinder of dimensions (r, z) = [4], [15] mm, with a constant voltage of 65 kV applied at z = 15 mm, and 0 kV at z = 0 mm. In this study, a simplified set of plasma chemical reactions for air was used, reaction rates following Pancheshnyi and Starikovskii [58], while electron transport parameters were obtained using BOLSIG+ [49] with Phelps' cross-sectional data [59], [60], [61]. These can be found tabulated together in Appendix C. Equal densities of N 0 = 5× 10 18 m -3 electrons and N + 2 ions were placed at the center of the domain, again following the gaussian form of (31), but with s = 0.2 mm. A background ionisation level of 10 9 m -3 was again used, as was photoionization. Fig. 6 shows the electric field distribution at 2, 4, and 6 ns, while Fig. 7 shows the corresponding electron density. The double-headed streamer was successfully resolved, and the propagating characteristics of both positive and negative streamers are in agreement with past studies [55].

C. COUNTERPROPAGATING STREAMERS TOWARD EACH OTHER FROM NEEDLE ELECTRODES UNDER FAST-RISING RAMP VOLTAGE
In practical electrical systems, such as HV power and pulsed power equipment, electrical pre-breakdown and breakdown VOLUME 11, 2023 phenomena are seldom initiated in regions of homogeneous electric fields. In most cases, highly nonuniform and divergent field conditions, induced by high aspect ratio geometry, field redistribution at triple-junctions, or accumulated spacecharge, typically form the most pressing problems in HV insulating system design.
Here, StrAFE's handling of curved geometries in the form of sharp needle electrodes is demonstrated, alongside a time-varying Dirichlet condition representing a fast-rising applied voltage. The simulation in question was initially performed in [62], between two needle electrodes with a curvature radius of 25 µm, and with an interelectrode gap distance of 1.2 mm. Once again, the domain is axisymmetric around r = 0. In the original study, a waveform generated using the CST studio suite [63] from their experimental configuration was also used for simulation, but was not made openly available. However, since the simulation was performed only over the initial stages of the discharge, and during the rising edge, an alternative waveform was used in this study, which closely approximates the rising-edge of the original signal, using a linearly increasing ramp voltage of the form: where the rate of voltage rise dU /dt was chosen to be 75 kV/ns, to roughly align with the slope of the original waveform from [62]. The plasma chemistry model of Appendix C and photoionization was once again used, with initial background densities of 10 9 m -3 for electrons and N + 2 ions. No initial seed was necessary in this simulation, as the streamers would initiate directly from the electrode tips, where the electric field was maximally enhanced. The LMEA was additionally used in this simulation, with energy-dependent electron transport data computed using BOLSIG+ [49] with Phelps' cross-sectional data [59], [60], [61]. Results are presented in the same format as in [62] to facilitate direct comparison. These include Fig. 8(a) and 8(b), which show streak plots of the electric field magnitude and electron density along the axis of symmetry, respectively; and Fig. 9, which plots the electric field strength along the same axis, at timesteps between t = 160 and 260 ps.
Overall, excellent agreement with the results from [62] was found, including the time of streamer inception, the delayed propagation of the negative streamer, the electron density distribution, and propagation velocities. Minor differences are observed in the maximum electric field strength at the negative streamer head and in the position where the streamers collide. As indicated by the difference between the black markers of Fig. 9 and the obtained data, the positive streamer field agrees closely with [62], but there exists a larger discrepancy for peak field values with the negative streamer. However, the difference is within expectation considering the many possible differences in numerical implementation. The discrepencies are believed to be attributed to the differences in the waveform, the numerical implementation, or possible dif-  . Electric field magnitude down the axis of symmetry for various timesteps, during the development of counterpropagating streamers, generated using StrAFE. This figure corresponds to Fig. 9(a) of [62]. Solid black markers indicate the position of the peaks from the simulations conducted in [62]. ferences in the boundary conditions applied at the electrodes, as these were not specified in [62].

D. POSITIVE STREAMER WITH SOLID DIELECTRIC
Equally prevalent in, and important to, practical HV systems, are the interactions between nonthermal plasma discharges and dielectric surfaces. Knowledge of these processes would be hugely beneficial to the development of novel plasma technologies, e.g., surface treatment technologies [2], or any other apparatus involving dielectric barrier discharges [7].
Using the described subdomain support in section IV-D, simulations from [17] have been recreated, showing the electrostatic attraction of streamers toward solid dielectric surfaces, and their subsequent propagation across these surfaces. Briefly, the domain consisted of a square geometry with (x, y) = [40,40] mm. A solid dielectric subdomain was defined where 0 ≤ x ≤ 10 mm, and with relative permittivity ε r = 2, while a neutral seed composed of electrons and positive ions initiated the discharge (plasma chemistry was not used in this simulation, following the original. Only electrons, generic positive ions, and generic negative ions were included). The seed followed the capsule definition as in [17], placed with its center 0.5 mm away from the solid surface, and near the anode that was set to 100 kV. Townsend coefficients were computed using BOLSIG+ [49] as before. The base mesh was refined along the gas-solid interface, while AMR was enabled to perform a remesh every 30 iterations -this value was once again determined from manual testing. No significant difference in the solution was observed when only 20 iterations per AMR was used, but signs of numerical nonconvergence near the fast-moving streamer head at 40 iterations per AMR were observed. 30 was therefore chosen to reduce computational time and to provide good convergence, yet have minimal effect on the computed solution. Fig. 10 shows the time evolution of the streamer electron density at t = 10, 12.5, 15, and 20 ns. The streamer initiates from the charged seed due to the field enhancement at the bottom of the seed, as electrons drift upward. The streamer then begins to propagate before 10 ns, and is drawn toward the solid dielectric surface, before evolving into a surface streamer. Here, the streamer is observed to gain velocity, and accelerates down the surface, as in [17] and [18]. StrAFE was also shown capable of resolving the thin sheath region between the positive streamer and the surface, as described in [17]. The maximum electric field strength at the streamer head was found to be very similar to that of [17], as was the propagation velocity during the gas-stage and surface streamer stage. The electron density was similarly comparable, with the channel density stabilizing at around 10 21 m -3 . Additional studies on the effects of increasing the solid permittivity were also conducted, results which recreated characteristics as observed in [17], such as the thinner surface streamer with increased ε r , and a far more rapid surface attachment speed.

VI. CONCLUSION
In summary, the present article has introduced and described the StrAFE (Streamers on Adaptive Finite Elements) Python library, a powerful automation tool for the study of transient ionization fronts, using the hydrodynamic approximation.
Based on the open-source FEniCS framework, this article has presented the design and development of the library, showing how FEniCS and Python can be used to develop an effective platform for the simulation of nonthermal gas discharge plasmas, with emphasis on configurability, easeof-use, and a high degree of flexibility.
StrAFE was developed in Python, known for its accessibility, due to its English-like syntax. The implementation of a number of dedicated classes allows the substantial simplification of generating transient ionization-FE problems, while also streamlining the process of passing the problem to FEniCS for efficient solving, using its highspeed C++ backend. StrAFE employs the hydrodynamic approach of coupled drift-diffusion-reaction-Poisson equa-FIGURE 10. Positive streamer attaching to a solid dielectric surface (solid red line), and forming a surface streamer at timesteps t = 10, 12.5, 15, and 20 ns, generated using StrAFE. Equipotential lines are spaced by 5 kV. Originally studied in [17].
tions for nonthermal plasmas, and further includes crucial features necessary for multiscale modelling. These include: user-configurable adaptive mesh refinement and adaptive timestep controllers, support for problem generation from a single plasma chemistry file; support for both local field approximations and local mean energy approximations, and support for subdomains, e.g., solid dielectric materials. Since StrAFE inherits from FEniCS, it additionally offers full support for distributed memory parallelism through the message passing interface (MPI), and the same desktop-developed prototype code can be executed in HPC settings with little to no change.
This article has importantly showed how it can be used to significantly simplify the process of simulating streamer discharges. Several simulation examples have been presented in various domains, and were of varying complexity, which demonstrated the capabilities of StrAFE through comparison with existing work. It was found that the use of FEniCS compares very well to a range of existing and well-established codes, both commercial and custom-made.
Yet, its open-source nature offers far greater flexibility and code transparency, and can be used to develop user-friendly software with a very low entry barrier. Such innovations are exceptionally well-suited for interested parties who may be new, or unfamiliar, with the field of gas discharge modelling, and thus possesses immense educational value. From this study, it can be concluded that open-source frameworks such VOLUME 11, 2023 FIGURE 11. (a) Load balanced mesh, each color demarcates nodes which are owned by a unique computational process, (b) adaptive mesh around an initially gaussianly-distributed charge seed, (c) example of the dynamic AMR scheme tracking a streamer head, overlaid on top of the electric field color map. as FEniCS, have great potential to make significant contributions to the study of nonthermal plasmas, ionization fronts, and gas discharges; and can help to further increase the accessibility of computational modelling of complex phenomena in engineering and science.
In future, StrAFE will be used to explore nonthermal plasma dynamics in practical topologies of interest. The flexibility offered by the library further facilitates experimentation with different numerical schemes, the exploration of higher-order discretization, the possibility to use discontinuous Galerkin methods (DG-FEM), additional physics, and more. It is also of high interest to conduct three-dimensional simulations in the future, and to further evaluate StrAFE in HPC settings -3D features have theoretically been implemented, but not yet thoroughly tested. Fig. 11(a) shows a visualization of the load-balanced mesh generated by FEniCS, around the initial seed of the example from section V-A. Fig. 11(b) shows the corresponding mesh around the initial seed region. Mesh cells belonging to the same color group are owned by the same MPI process. Fig. 11(c) additionally shows the adaptive mesh overlaid on half of the domain at t = 8 ns, showing the fine mesh tracking the streamer head. Fig. 12 shows the results of a basic parallel scaling test, up to 16 MPI processes, conducted by varying the number of processes the program ran with, when simulating example V-A. Note that only physical cores were used here, the use of virtual cores was disabled, since it did not appear to offer any additional computational speedup. Table 1 encloses the plasma chemical reactions for air, which were included in several of the simulations presented in section V.