Features

FEDM is designed to make the implementation and solution procedure of the model in FEniCS easier. The code automates the process of setting up the equations and source terms, reducing the model implementation time and probability of introducing errors in the code. Some general FEniCSand newly introduced features used in the solution procedure are described as follows.

1. Reading the input data. In this step, the configuration files for the given problem are read to obtain particle properties and the reaction scheme, and the transport and the reaction rate coefficients. In this way, the number of particles for which the problem is solved is imported into the code. The reaction matrices containing the partial reaction orders and the stoichiometric coefficients for the given species, required for the rate and the source term definitions, are read from the reacscheme.cfg. The file contains symbolically written reactions, dependences, energy losses or gains, and names of the files containing the reaction rate coefficient values. The transport and rate coefficients are read from the .cfg files and stored in the lists. Depending on the problem, they can be imported in the form of a function (written as a string of Python code in the .cfg file) or in the form of look-up tables. The files start with a commented header containing information about the coefficients, such as a reference from which the value is taken and the dependence. The transport and rate coefficients may depend on the reduced electric field, the mean electron energy, the electron and/or gas temperature, or may be a constant. The coefficients are stored as functions (written as a string of the Python or C++ code in the .cfg file) or in the form of look-up tables. The files start with a header, containing information about the coefficients such as a reference from which the coefficient value is taken and its dependence, and the coefficient values. In cases when the dependence is the constant or a function, only one value or function is stored in a file. Otherwise, data are stored in two columns, one representing the reduced electric field, the mean electron energy, or the electron temperature and the other representing values of the coefficients. Note that a single blank space must be used for separating the columns, instead of a tab. Also, all files should use utf8 encoding. The following functions are used for reading the input data:

   • read_speclist() is used to read the particle species list from the input file stored in the file_input folder, and obtain the number of species, their names and corresponding file names

   • read_particle_properties()is used to read the particle properties from the input files stored in the file_input folder. The function returns the mass and the charge of the particles

   • reaction_matrices() is used to read the reaction matrices from the reaction scheme

   • rate_coefficient_file_names() are used to read the name of the files containing the rate coefficient from the reaction scheme file

   • energy_loss() is used to read the energy loss for the given reaction from the reaction scheme file

   • reading_transport_coefficients() is used to read the transport coefficients from the input files

   • reading_rate_coefficients() is used to read the rate coefficients from the input files

2. Mesh for FEM discretisation. The mesh is either generated using a built-in function (structured mesh) or imported from the external (.xml or .xdmf) file.

3. Function spaces. In the examples presented here, a fully coupled solution approach is used, so the native mixed element functions are used for the function space definition. The mixed element list is obtained using the custom function Mixed_element_list() which creates a list of the function elements. Since the number of functions varies from problem to problem, the procedure of function definition is automatised by calling the custom function Function_definition().

4. Domain boundaries. Depending on the given problem, one may need to mark specific boundaries. This can be done using the custom function Marking_boundaries().

5. Evaluation of transport and rate coefficients. The imported transport and rate coefficients are interpolated in the code when needed.

6. Source term definition. Particle source term definition is done using the function Source_term(), while the energy source term definition is done using Energy_Source_term(), using the imported reaction matrices (described in step 1).

7. Fluxes. Depending on whether the transformation of the variables is done or not (e.g. logarithmic representation), the fluxes are defined using two functions:

   • For the case when the transformation of variables is not done, function weak_form_balance_equation() is used

   • For the logarithmic representation of the equation (where the problem is solved for the logarithm of the variable) weak_form_balance_equation_log_representation() is designed.

8. Weak formulation. For the small examples, the weak formulation can be done manually for each equation. However, for the cases where the problem is solved for multiple particle species, this can become tedious. In this case, one may use a set of functions to define the variational form of equations. For the balance equations of particle species, two functions exist:

   • For the case when the transformation of the variables is not done, function weak_form_balance_equation() is used

   • For the logarithmic representation of the equation (where the problem is solved for the logarithm of the variable) weak_form_balance_equation_log_representation() is used

     Both functions implement the backward differentiation formula (BDF) of the second order with variable time-step. In practice, the change of order can be done without redefining the time-stepping scheme in the time-loop. Instead, by defining the previous time-step size dt_old as some large number, the adaptive BDF is reduced to the first order. Moreover, by defining that dt_old = dt, adaptive BDF can be reduced to BDF with a constant time-step size. Therefore, the seamless initiation of the formula and switching between two orders is possible.

   • For the electron energy balance equation, the same function may be used as for particle balance equations.

   • For Poisson’s equation, a variational form may be defined manually or using the weak_form_Poisson_equation() function.

9. Boundary conditions. The boundary integrals (if required) are defined separately using the function Boundary_flux() and added to the weak formulation.

10. Solver. FEniCS has the option to access PETSc solvers available in the PETSc library. The function PETScSNESSolver() is used to set up the nonlinear solver. The advantage of this approach is that it is possible to tune the nonlinear and linear solver parameters, such as relative or absolute tolerance, line search type in SNES, the maximum number of iterations, etc. It should be noted that it is possible to define your own custom solvers (see the Cahn-Hilliard equation demo in the FEniCS repository).

11. The file output. The linear interpolation of the calculated results is used to determine the output value for the particular time step, which are then saved into the file using the function file_output().

12. Adaptive time stepping. The time-step size control is based on the L2 norm difference between the previous and the current step. Depending on the used approximation, the relevant variable may be mean electron energy, electron number density or the whole solution vector. The new time step is obtained using the function adaptive_timestep() which utilises the PID controller, described in the article.