Pyrit: A finite element based field simulation software written in Python

Abstract

Purpose

The purpose of this paper is to present the freely available finite element simulation software Pyrit.

Design/methodology/approach

In a first step, the design principles and the objective of the software project are defined. Then, the software’s structure is established: The software is organized in packages for which an overview is given. The structure is based on the typical steps of a simulation workflow, i.e., problem definition, problem-solving and post-processing. State-of-the-art software engineering principles are applied to ensure a high code quality at all times. Finally, the modeling and simulation workflow of Pyrit is demonstrated by three examples.

Findings

Pyrit is a field simulation software based on the finite element method written in Python to solve coupled systems of partial differential equations. It is designed as a modular software that is easily modifiable and extendable. The framework can, therefore, be adapted to various activities, i.e., research, education and industry collaboration.

Research limitations/implications

The focus of Pyrit are static and quasistatic electromagnetic problems as well as (coupled) heat conduction problems. It allows for both time domain and frequency domain simulations.

Originality/value

In research, problem-specific modifications and direct access to the source code of simulation tools are essential. With Pyrit, the authors present a computationally efficient and platform-independent simulation software for various electromagnetic and thermal field problems.

Citation

Bundschuh, J., Ruppert, M.G. and Späck-Leigsnering, Y. (2023), "Pyrit: A finite element based field simulation software written in Python", COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. 42 No. 5, pp. 1007-1020. https://doi.org/10.1108/COMPEL-01-2023-0013

1. Introduction

Today, electromagnetic (EM) field simulation software is essential for the design of electric equipment, e.g., electric power converters and transmission systems, electronics and particle accelerator components. Many different tools for EM field simulation based on the finite element (FE) method have been developed over time. Commercial software tools, e.g., Ansys Maxwell^® (Ansys Maxwell, 2022), CST EM STUDIO^® (CST EM Studio, 2022), Flux2D/3D^® (Altair Flux, 2022), JMAG^® (JMAG, 2022) and COMSOL Multiphysics^® (Comsol Multiphysics, 2022), are designed to solve universal EM field problems with standard computational methods. The design principles of commercial packages are, first of all, a broad applicability. Second, the software is made easy to use by powerful graphical user interfaces and template-based workflows. However, users cannot customize and access the internal routines and system solvers.

In a research setting, however, problem-specific modifications and direct access to the FE matrices is indispensable. Various available freeware FE tools offer different levels of abstraction: General FE research tools with a focus on numerical mathematics are, for example, Deal (Arndt et al., 2022), FEniCS (Alnæs et al., 2015) and GetDP (Dular et al., 1998). However, these tools typically require an in-depth knowledge of the FE method. Besides, research tools motivated by electrical engineering problems are, for example, Agros (Karban et al., 2013), FEMM (Meeker, 2022) and openCFS (openCFS, 2022) and for electrical machine problems Pyleecan (Bonneel et al., 2018). These tools allow to study and further develop formulations, discretization techniques and solver strategies for EM field simulation, thereby covering a large fraction of research in a computational EM research group. Nevertheless, there is an additional need for a research code built from scratch, with access to all basic routines, with full flexibility and possibly adapted to a few specific lines of research.

Pyrit [1] is a FE solver developed from the perspective of the electric field simulation workflow. To allow a fast prototyping and due to its popularity in industry and academia, it is written in Python [2]. In general, it solves static and quasistatic EM and heat transfer problems. One of the central aims of Pyrit is to provide a user-friendly and computationally efficient FE code with a template-based structure for standard and coupled EM and thermal problems. Furthermore, a plain user interface via Python scripts supports the implementation of new modeling and simulation ideas. Its inherent structure supports students to learn how to implement and use an FE solver.

2. Design principles

Pyrit is developed collaboratively based on the following principles: it is robust, portable, scriptable, extensible and thoroughly documented. Users organize the simulation workflow with Python scripts or Jupyter notebooks (Perkel, 2018). Therefore, the full capabilities of Python can be exploited, and studies can be tailored to the needs of individual users.

The object-oriented structure of Pyrit allows, in contrast to other FE software packages, different levels of abstraction: On one hand, the user can access and customize the basic and physics-independent FE routines. On the other hand, Pyrit offers classes that serve as templates and predefined solvers for a set of selected EM and thermal problems. This allows to model and simulate field problems without detailed knowledge of the weak formulation or of the FE method. In addition, Pyrit provides standard post-processing routines as well as the possibility to export results to Paraview and LaTeX (for plotting with pgfplots). Pyrit is user-friendly, as all packages, modules, classes, functions and methods are thoroughly documented. The documentation of the whole software is stripped directly from the source files and is therefore always up-to-date. Supplementary information such as tutorials further extends the documentation. Furthermore, the geometry and mesh generation of the open source software Gmsh (Geuzaine and Remacle, 2009) is used. To make the geometry generation more convenient, the package geometry provides a unified interface to Gmsh.

3. Software structure

Pyrit’s modeling and simulation workflow is organized in packages. They are assigned to the three steps in Figure 1 that are described in this section. Because the toolbox package can be associated with all three steps, it is described here:

• toolbox: This package is a collection of different toolboxes with additional functionalities that are not solely associated with one package. It contains, inter alia, a material and geometry library (for the problem definition step), a circuit simulation and a time integration toolbox (for the problem solving step) and a post-processing and export toolbox (for the post-processing step).

3.1 Problem definition

Before a simulation is executed, Pyrit defines a field problem. This first step is the choice of the problem type. This includes the differential equation to be solved, the dimension of the problem and the time dependency. Subsequently, the geometry is built, and the materials, boundary conditions and excitations are assigned. The assignment relies on the concept of physical groups of Gmsh, as indicated by an UML class diagram in Figure 2. The pendant in Pyrit are regions. A region has one or more geometrical entities that share a material, a boundary condition or an excitation. Consequently, a geometrical entity can have at most one region.

The following packages are associated with the problem definition step:

• bdrycond: In this package, boundary conditions are defined. Besides the standard Dirichlet and Neumann boundary conditions, also Robin, periodic, antiperiodic and floating boundary conditions are supported (De Gersem et al., 2004b). There is a container class that manages all boundary conditions of a problem and can apply them to the system of equations of the corresponding problem.

• excitation: In this package, excitations for problems are defined. Its structure is similar to bdrycond. Typical excitations are charges, charge densities, currents or current densities. In addition to that, the field-circuit coupling models solid and stranded conductor (De Gersem et al., 2004a) are also implemented in this package. Like in the previous package, there is also a container class that manages all excitations of a problem.

• geometry: In this package, the interaction between Pyrit and Gmsh is handled via the Gmsh Python API. This includes primarily the generation of geometries and meshes and the extraction of the relevant data for Pyrit. There are also several wrapper classes for geometrical entities to simplify the geometry generation.

• material: In this package, the management of material information is implemented. Materials can be defined by giving them a set of material properties, e.g., a conductivity or a permeability. The values of the material properties can be constants or functions that can depend, in particular, on space, time or other field values. They can be scalar or tensor valued. With this, a nonlinear, inhomogeneous and anisotropic material can be defined. A container class manages all the materials of a problem.

• region: In this package, classes for the organization of materials, boundary conditions and excitations on the mesh are provided. Regions can be seen as discrete counterpart to physical groups and are defined on the mesh. The package includes a class for defining a region and a container class that manages different regions.

3.2 Problem solving

Now that the problem is completely defined, it can be solved. For that, a mesh object and a shape function object are required. The former can be generated and extracted from Gmsh with the geometry package. The latter is fixed by the chosen problem. With the shape function object, the required matrices and vectors can be computed and assembled into the final system of equations. This can be solved using a predefined solve routine defined in the problem class or any other solve method. Currently, there exist wrappers for direct and iterative solvers for sparse matrices from SciPy [3] and Pardiso [4].

The relevant packages in this step are:

• mesh: In this package, classes for different kinds of meshes are defined, e.g., a triangular mesh for two-dimensional (2D) domains or a tetrahedral mesh for three-dimensional domains. These classes are responsible for the mesh data. They store the coordinates of all nodes and the definition of the higher dimensional entities (edges, triangles and tetrahedra).

• problem: This package collects classes that represent standard problem formulations. They provide a framework for organizing all the data needed for defining a problem. Furthermore, the classes in the problem package implement a solve method that provides a convenient way to solve a problem, i.e., to setup and solve the system of equations. There are separate classes for static, harmonic and transient problems per problem type, i.e., the underlying differential equation.

• shapefunction: This package provides the core of the FE method, i.e., the different shape functions. There is one class per FE shape function and dimension of the domain. For example, there is a class for nodal shape functions in cylindrical coordinates or a class for edge shape functions in 2D Cartesian coordinates. These classes implement the routines to compute the FE matrices. Currently, nodal shape functions are implemented for two- and three-dimensional Cartesian and cylindrical coordinates. Furthermore, edge functions are implemented for 2D Cartesian and cylindrical coordinates.

3.3 Post-processing

The result from the problem-solving step is a vector or several vectors, depending on the problem, respectively. Each problem class is equipped with a corresponding solution class. These implement standard post-processing routines and allow for a straightforward visualization of the results. This step also includes the export of the data to files of different formats.

The post-processing step comprises the following package:

• solution: This package collects solution classes for standard problem formulations. For each class in the problem package, there is an associated class in this package. The problems solve routine returns a solution object of correct type that contains the simulation result. In the solution class, standard post-processing and visualization routines, specific for this problem type, are implemented. In the case of an electrostatic problem, this includes, for example, the computation and visualization of the electric displacement field from the electric scalar potential.

4. Quality management

In the development of Pyrit, state-of-the-art software engineering principles were applied to ensure a good code quality and a high usability and reproducibility. This section addresses testing, linting and the documentation of the code. Continuous integration (CI) ensures that the code fulfills these requirements.

A dedicated focus lies on developing the software collaboratively. For that reason, the version control system Git [5] is used. It keeps track of the evolving software project and serves as a backup. It also helps distributing the code between the developers and users.

To ensure the correctness and robustness of the code in presence and in future, there are many tests for Pyrit. They range from testing the functionality of single functions and methods over the interaction of modules to integration tests, i.e., the test of whole simulation procedures. The tests of the shape functions cover, in particular, the consistency and convergence of the implemented FE method.

The coding style, which can vary depending on the collaborator, and the readability of the code are optimized through the use of a linter, i.e., a static code analyzer that finds discrepancies between the code and a predefined style.

A crucial part for every software is its documentation, as it guides the user through the software’s functionalities. In Pyrit, the documentation is automatically generated based on the provided docstrings of the source files. This implies that the documentation file is always up-to-date. Furthermore, the documentation is enriched by application examples and tutorials of EM and thermal field problems.

The tasks in the described aspects of quality management, i.e., testing, linting and documentation, have to be executed regularly, but at least after a commit has been pushed. This is automated on Git with CI. After new commits have been pushed, a so-called pipeline is started where first Pyrit is installed with the recent dependencies and, then, tested and linted. The single tasks are executed via Docker [6] in containers on a separate machine (Merkel, 2014).

5. Examples

Three examples show the modeling and simulation workflow in Pyrit. The first two examples use a predefined problem class of Pyrit, i.e., the user models the EM problem. Subsequently, the predefined template is used to setup, solve and postprocess the FE system. The third example shows how Pyrit can be customized. There, a system of equations from the FE formulation is built and extended by a field-circuit coupling.

5.1 Template-based simulation example

As a starting point, a simple example shows the basic procedure of the problem definition. It is an electrostatic simulation of a plate capacitor in 2D Cartesian coordinates filled with two materials [see Figure 3(a)]. The bottom plate is grounded, and the top plate is set to a potential of 1 V. On the left and right side, homogeneous Neumann boundary conditions are imposed.

We start at a point where the mesh and the regions already exist. The geometry was build with the Gmsh Python API and imported into Pyrit with the geometry package:

Next, the materials and boundary conditions of this example are defined:

With the region IDs from Figure 3(a), materials and boundary conditions can be assigned to the regions:

With that, the problem is completely defined, and all data structures necessary to instantiate an object of the problem class ElectricProblemCartStatic are available. The solve routine then solves the problem and returns an object of the associated solution class ElectricSolutionCartStatic:

The solution object can be used to plot, for instance, the electric field strength [see Figure 3(b)]:

5.2 Template-based nonlinear simulation example

The second example is a nonlinear electrothermal simulation of a cylindrical resistor (see Figure 4). The resistor consists of a nonlinear material, which features a strongly field- and temperature-dependent electric conductivity (see Figure 5). Such nonlinear materials are used in field grading layers of high-voltage cable systems (Hussain and Hinrichsen, 2017; Ruppert et al., 2023). The material is located between two electrodes to which a transient voltage V(t) is applied (see Figure 6). The resistor is surrounded by a layer of soil. The inner surface of the resistor and the outer surface of the soil are at a fixed temperature of 60°C and 20°C, respectively. For a more detailed description of the model configuration and material parameters, see Ruppert et al. (2023).

The first step in the simulation process is the definition of the geometry. Pyrit offers several ways to do this. The user can import a step file, create a geometry using Pyrit’s geometry package or import an already existing Gmsh file. For this example, an already existing Gmsh file is imported:

Second, the materials and boundary conditions are defined. The excitation voltage and fixed temperatures are implemented as Dirichlet boundary conditions of the electric potential and temperature, respectively:

The homogeneous Neumann boundary conditions of the lower and upper boundary in the z-direction are automatically satisfied by the FE ansatz.

Finally, the material properties and boundary conditions are assigned to the geometry. A region maps the physical properties to corresponding edges and elements in the FE mesh:

As stated in Section 3, Pyrit offers classes for EM and heat transfer problems of any time dependence. This allows to solve field problems based on a template-based workflow. Here, the class ElectrothermalProblemAxiTransient, which represents a coupled transient electrothermal problem in cylindrical coordinates, is used. It is composed of two subproblems containing the boundary conditions of the electric and thermal subsystem, respectively:

The ElectrothermalProblemAxiTransient class offers a solve method for running the electrothermal simulation. It automatically handles the coupling of the two subproblems either by a weak coupling scheme or a successive substitution iteration. The nonlinearity arising in the electric subproblem due to the nonlinearity of the nonlinear material is solved using a damped Newton method:

The solutions returned by the solve method come with a set of post-processing and plotting routines:

Figure 7 exemplarily shows the plot generated by the method plot_temperature_at_position.

5.3 Customized example

In the third example, no predefined problem class is used. A field-circuit coupled problem with a magnetoquasistatic field problem is considered (Figure 8). It consists of an electrical circuit [see Figure 8(a)] with a transformer that is represented by a field model [see Figure 8(b)]. The lumped elements in the circuit are R=1 Ω, RL=10 Ω and C = 0.1 µF. The voltage source is Vs=V^ssin⁡(ωt), with V^s=1 V, the angular frequency ω=2πf, the frequency f=3 kHz and the time variable t. The field model of the transformer uses a stranded conductor model for both coils (De Gersem et al., 2004a). The primary coil (teal) and the secondary coil (orange) have 300 and 500 turns, respectively.

In the following, one system for the magnetoquasistatics and one for the circuit is built in the frequency domain. These are then coupled with field-circuit coupling matrices.

The geometry and the mesh of this example were generated and exported to Pyrit with the Gmsh Python API. Because the focus of this example is on the field-circuit coupling, these steps are not shown here. We start with an object on which the shape functions, the mesh and the materials and boundary conditions already exist:

Now, the right-hand side vector and the FE matrices are retrieved from the shape function object and combined into the system of equations of the magnetoquasistatic formulation:

The system representing the circuit [see Figure 8(a)] is built with the values of the lumped elements:

The coupling matrices from the field-circuit coupling are generated by a class that represents stranded conductors:

Inserting the boundary conditions with the method yields the system matrix and the system vector. With these, the new system is given by:

Finally, we insert the circuit system:

After solving the system, the circuit solution consists of the voltages and currents at the transformer. The boundary conditions are integrated into the field solution that is then evaluated with a predefined solution class for magnetoquasistatic problems:

The object provides post-processing routines:

The plot of the magnetic flux density at a phase of 0, generated by the method, is shown in Figure 9.

6. Conclusion

Pyrit is a FE software environment for simulating EM and thermal field problems of any time dependence (static, transient or harmonic). Its software structure is based on an engineering workflow for solving a field problem numerically, i.e., problem definition, problem solution and post-processing. Pyrit allows to solve field problems in a template-based workflow. Furthermore, full access to the FE core supports research and development on numerical techniques. Three examples showed how users can define, solve and post-process a field problem in an easy and general way.

Figures

Figure 1.

Assignment of the packages to the steps “problem definition,” “problem solving” and “post-processing” of a field simulation workflow

Figure 2.

UML class diagram of the assignment of materials, boundary conditions and excitations to physical groups. They, in turn, are then assigned to geometrical entities as points, lines, triangles and tetrahedra. An object of a class can have at most one object of a class to its right but one or more objects of a class to its left

Figure 3.

Plate capacitor with two materials. In (a), the model is shown with the region IDs in blue. In (b), the electric field strength is plotted

Figure 4.

Schematic of the second example. The cylindrical resistor consists of a nonlinear material (blue) with an inner radius of 0.1 m and an outer radius of 0.3 m. The material is located between two electrodes to which a transient voltage V(t) is applied, where t denotes the time. The resistor is surrounded by a layer of soil (gray) with a thickness of 0.7 m. The inner surface of the resistor and the outer surface of the soil are at a fixed temperature of 60°C and 20°C, respectively. The surrounding space is nonconductive

Figure 5.

Field- and temperature-dependent conductivity σ(E,T) of the nonlinear material. The electric field strength and the temperature are denoted as E and T, respectively

Figure 6.

Transient voltage V(t) over the time t

Figure 7.

Plot generated by the method plot_temperature_at_position of the temperature evaluated at ρ = 0.2m inside the resistor depicted in Figure 4 over time

Figure 8.

Configuration of the third example. The circuit for the field-circuit coupling (a) with a transformer resolved in a field model (b). The transformer is described in axisymmetric coordinates and consists of a highly permeable yoke (gray) and a primary windings (teal) and secondary windings (orange)

Figure 9.

Absolute values of the magnetic flux density at phase 0 in the transformer

Notes

1. Pyrit’s public wiki: https://git.rwth-aachen.de/jonas.bundschuh/pyrit-wiki

2. Python website: www.python.org/

3. https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html

4. www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2023-0/onemkl-pardiso-parallel-direct-sparse-solver-iface.html

5. Git website: https://git-scm.com/

6. Docker website: www.docker.com/