PETSc Library Research Articles

Parallel finite-element codes for the Bogoliubov-de Gennes stability analysis of Bose-Einstein condensates

We present and distribute a parallel finite-element toolbox written in the free software FreeFEM for computing the Bogoliubov-de Gennes (BdG) spectrum of stationary solutions to one- and two-component Gross-Pitaevskii (GP) equations, in two or three spatial dimensions. The parallelization of the toolbox relies exclusively upon the recent interfacing of FreeFEM with the PETSc library. The latter contains itself a wide palette of state-of-the-art linear algebra libraries, graph partitioners, mesh generation and domain decomposition tools, as well as a suite of eigenvalue solvers that are embodied in the SLEPc library. Within the present toolbox, stationary states of the GP equations are computed by a Newton method. Branches of solutions are constructed using an adaptive step-size continuation algorithm. The combination of mesh adaptivity tools from FreeFEM with the parallelization features from PETSc makes the toolbox efficient and reliable for the computation of stationary states. Their BdG spectrum is computed using the SLEPc eigenvalue solver. We perform extensive tests and validate our programs by comparing the toolbox's results with known theoretical and numerical findings that have been reported in the literature. Program summaryProgram Title: FFEM_BdG_ddm_toolbox.zipCPC Library link to program files:https://doi.org/10.17632/w9hg964wpb.1Licensing provisions: GPLv3Programming language:FreeFEM (v 4.12) free software (www.freefem.org)Nature of problem: Among the plethora of configurations that may exist in Gross-Pitaevskii (GP) equations modeling one or two-component Bose-Einstein condensates, only the ones that are deemed spectrally stable (or even, in some cases, weakly unstable) have high probability to be observed in realistic ultracold atoms experiments. To investigate the spectral stability of solutions requires the numerical study of the linearization of GP equations, the latter commonly known as the Bogoliubov-de Gennes (BdG) spectral problem. The present software offers an efficient and reliable tool for the computation of eigenvalues (or modes) of the BdG problem for a given two- or three-dimensional GP configuration. Then, the spectral stability (or instability) can be inferred from its spectrum, thus predicting (or not) its observability in experiments.Solution method: The present toolbox in FreeFEM consists of the following steps. At first, the GP equations in two (2D) and three (3D) spatial dimensions are discretized by using P2 (piece-wise quadratic) Galerkin triangular (in 2D) or tetrahedral (in 3D) finite elements. For a given configuration of interest, mesh adaptivity in FreeFEM is deployed in order to reduce the size of the problem, thus reducing the toolbox's execution time. Then, stationary states of the GP equations are obtained by a Newton method whose backbone involves the use of a reliable and efficient linear solver judiciously selected from the PETSc1 library. Upon identifying stationary configurations, to trace branches of such solutions a parameter continuation method over the chemical potential in the GP equations (effectively controlling the number of atoms in a BEC) is employed with step-size adaptivity of the continuation parameter. Finally, the computation of the stability of branches of solutions (i.e. the BdG spectrum), is carried out by accurately solving, at each point in the parameter space, the underlying eigenvalue problem by using the SLEPc2 library. Three-dimensional computations are made affordable in the present toolbox by using the domain decomposition method (DDM). In the course of the computation, the toolbox stores not only the solutions but also the eigenvalues and respective eigenvectors emanating from the solution to the BdG problem. We offer examples for computing stationary configurations and their BdG spectrum in one- and two-component GP equations.Running time: From minutes to hours depending on the mesh resolution and space dimension.

ReMKiT1D - A framework for building reactive multi-fluid models of the tokamak scrape-off layer with coupled electron kinetics in 1D

In this manuscript we present the recently developed flexible framework for building both fluid and electron kinetic models of the tokamak Scrape-Off Layer in 1D - ReMKiT1D (Reactive Multi-fluid and Kinetic Transport in 1D). The framework can handle systems of non-linear ODEs, various 1D PDEs arising in fluid modelling, as well as PDEs arising from the treatment of the electron kinetic equation. As such, the framework allows for flexibility in fluid models of the Scrape-Off Layer while allowing the easy addition of kinetic electron effects. We focus on presenting both the high-level design decisions that allow for model flexibility, as well as the most important implementation aspects. A significant number of verification and performance tests are presented, as well as a step-by-step walkthrough of a simple example for setting up models using the Python interface. Program summaryProgram title: ReMKiT1DCPC Library link to program files:https://doi.org/10.17632/j47ym66xzj.1Developer's repository link:https://github.com/ukaea/ReMKiT1D and https://github.com/ukaea/ReMKiT1D-PythonLicensing provisions: GPLv3Programming language: Fortran, PythonSupplementary material:https://doi.org/10.14468/fdq7-z869Nature of problem: The flexible generation and modification of 1D models pertaining to multi-fluid simulations of the tokamak Scrape-Off Layer (SOL) with electron kinetics and reaction support. This would then allow both for rapid iteration on reduced models as well as the evaluation of kinetic electron effects in equilibria and during transients, following the formalism previously developed for SOL-KiT [1]. The framework was not only envisioned as the successor to SOL-KiT, but a tool that would allow users to construct their own models coupled with electron kinetics capabilities.Solution method: The framework is written heavily utilizing Object-Oriented design principles, in particular using an extended version of the puppeteer pattern as presented by Rouson et al. [2], as well as the heavy use of the strategy pattern/dependency injection. The Fortran code is MPI parallel and utilizes the PETSc library for implicit time-stepping. MPI parallelization is extended to distribution function Legendre harmonics, allowing for improved strong scaling. Initialization of the Fortran framework is done using JSON configuration files generated by an accompanying Python interface, and data analysis is standardized using widely used data formats such as HDF5.Additional comments including restrictions and unusual features: The present manuscript focuses on the design and high-level implementation of the framework, as well as the demonstration of the workflow and various verification and performance benchmarking tests. Some details are avoided for the sake of brevity at various points, and these are meant to be available as part of the general code documentation or tutorials offered on the main repositories.

Scalable implicit solvers with dynamic mesh adaptation for a relativistic drift-kinetic Fokker–Planck–Boltzmann model

In this work we consider a relativistic drift-kinetic model for runaway electrons along with a Fokker–Planck operator for small-angle Coulomb collisions, a radiation damping operator, and a secondary knock-on (Boltzmann) collision source. We develop a new scalable fully implicit solver utilizing finite volume and conservative finite difference schemes and dynamic mesh adaptivity. A new data management framework in the PETSc library based on the p4est library is developed to enable simulations with dynamic adaptive mesh refinement (AMR), distributed memory parallelization, and dynamic load balancing of computational work. This framework and the runaway electron solver building on the framework are able to dynamically capture both bulk Maxwellian at the low-energy region and a runaway tail at the high-energy region. To effectively capture features via the AMR algorithm, a new AMR indicator prediction strategy is proposed that is performed alongside the implicit time evolution of the solution. This strategy is complemented by the introduction of computationally cheap feature-based AMR indicators that are analyzed theoretically. Numerical results quantify the advantages of the prediction strategy in better capturing features compared with nonpredictive strategies; and we demonstrate trade-offs regarding computational costs. The robustness with respect to model parameters, algorithmic scalability, and parallel scalability are demonstrated through several benchmark problems including manufactured solutions and solutions of different physics models. We focus on demonstrating the advantages of using implicit time stepping and AMR for runaway electron simulations.

A comparison of Algebraic Multigrid Bidomain solvers on hybrid CPU–GPU architectures

The numerical simulation of cardiac electrophysiology is a highly challenging problem in scientific computing. The Bidomain system is the most complete mathematical model of cardiac bioelectrical activity. It consists of an elliptic and a parabolic partial differential equation (PDE), of reaction–diffusion type, describing the spread of electrical excitation in the cardiac tissue. The two PDEs are coupled with a stiff system of ordinary differential equations (ODEs), representing ionic currents through the cardiac membrane. Developing efficient and scalable preconditioners for the linear systems arising from the discretization of such computationally challenging model is crucial in order to reduce the computational costs required by the numerical simulations of cardiac electrophysiology. In this work, focusing on the Bidomain system as a model problem, we have benchmarked two popular implementations of the Algebraic Multigrid (AMG) preconditioner embedded in the PETSc library and we have studied the performance on the calibration of specific parameters. We have conducted our analysis on modern HPC architectures, performing scalability tests on multi-core and multi-GPUs settings. The results have shown that, for our problem, although scalability is verified on CPUs, GPUs are the optimal choice, since they yield the best performance in terms of solution time.

Shape sensitivity of thermoacoustic oscillations in an annular combustor with a 3D adjoint Helmholtz solver

We model acoustic oscillations driven by velocity-coupled heat release rate fluctuations. We obtain an inhomogeneous wave equation and convert it to the frequency domain with a modal decomposition. We impose acoustic boundary conditions and, using a finite element discretization with Lagrange elements, express this as a nonlinear eigenvalue problem for the complex angular frequency, ω. We solve this using the open source platform FEniCS combined with the SLEPc, PETSc and OpenMPI libraries. In Hadamard form we write the derivative of the eigenvalue, ω, with respect to generic geometry changes. This requires the solution of the adjoint equation, which we obtain with the same method as the direct equation. The output is a thermoacoustic Helmholtz solver that cheaply calculates the effect of generic shape changes on the growth rate and frequency of thermoacoustic oscillations. We then consider symmetry-preserving and symmetry-breaking geometry modifications. For demonstration we model a 30kW laboratory-scale annular combustor (MICCA from EM2C). We parametrize the surfaces of the three-dimensional geometry with NURBS using control points. For the plenum and the combustion chamber, we find the eigenvalue shape derivatives with respect to the parameters of these control points. We apply two different strategies, perpendicular boundary movements and control point perturbations, to implement shape changes proportional to these shape derivatives, thereby reducing the growth rate of the unstable mode by increasing the phase shift between the pressure and the heat release rate oscillations. This computational method shows how to significantly alter thermoacoustic oscillations by making small geometry changes. The framework in this paper can handle arbitrarily complex three-dimensional geometries, which will be useful for the design of industrial combustion systems.

Multirate partitioned Runge–Kutta methods for coupled Navier–Stokes equations

Earth system models are complex integrated models of atmosphere, ocean, sea ice, and land surface. Coupling the components can be a significant challenge due to the difference in physics, temporal, and spatial scales. This study explores multirate partitioned Runge–Kutta methods for the fluid–fluid interaction problem and demonstrates its parallel performance by using the PETSc library. We consider compressible Navier–Stokes equations with gravity coupled through a rigid-lid interface. Our large-scale numerical experiments reveal that multirate partitioned Runge–Kutta coupling schemes (1) can conserve total mass; (2) have second-order accuracy in time; and (3) provide favorable strong- and weak-scaling performance on modern computing architectures. We also show that the speedup factors of multirate partitioned Runge–Kutta methods match theoretical expectations over their base (single-rate) method.

A Parallel Solver for FSI Problems with Fictitious Domain Approach

We present and analyze a parallel solver for the solution of fluid structure interaction problems described by a fictitious domain approach. In particular, the fluid is modeled by the non-stationary incompressible Navier–Stokes equations, while the solid evolution is represented by the elasticity equations. The parallel implementation is based on the PETSc library and the solver has been tested in terms of robustness with respect to mesh refinement and weak scalability by running simulations on a Linux cluster.

MPI Parallelization of Numerical Simulations for Pulsed Vacuum Arc Plasma Plumes Based on a Hybrid DSMC/PIC Algorithm

The transport characteristics of the unsteady flow field in rarefied plasma plumes is crucial for a pulsed vacuum arc in which the particle distribution varies from 1016 to 1022 m−3. The direct simulation Monte Carlo (DSMC) method and particle-in-cell (PIC) method are generally combined to study this kind of flow field. The DSMC method simulates the motion of neutral particles, while the PIC method simulates the motion of charged ions. A hybrid DSMC/PIC algorithm is investigated here to determine the unsteady axisymmetric flow characteristics of vacuum arc plasma plume expansion. Numerical simulations are found to be consistent with the experiments performed in the plasma mass and energy analyzer (EQP). The electric field is solved by Poisson’s equation, which is usually computationally expensive. The compressed sparse row (CSR) format is used to store the huge diluted matrix and PETSc library to solve Poisson’s equation through parallel calculations. Double weight factors and two timesteps under two grid sets are investigated using the hybrid DSMC/PIC algorithm. The fine PIC grid is nested in the coarse DSMC grid. Therefore, METIS is used to divide the much smaller coarse DSMC grid when dynamic load imbalances arise. Two parameters are employed to evaluate and distribute the computational load of each process. Due to the self-adaption of the dynamic-load-balancing parameters, millions of grids and more than 150 million particles are employed to predict the transport characteristics of the rarefied plasma plume. Atomic Ti and Ti2+ are injected into the small cylinders. The comparative analysis shows that the diffusion rate of Ti2+ is faster than that of atomic Ti under the electric field, especially in the z-direction. The fully diffuse reflection wall model is adopted, showing that neutral particles accumulate on the wall, while charged ions do not—due to their self-consistent electric field. The maximum acceleration ratio is about 17.94.

Conservative Projection Between Finite Element and Particle Bases

Particle-in-cell methods employ particle representations of unknown fields but also employ continuum fields for other parts of the problem. Thus projection between particle and continuum bases is required. Moreover, we often need to enforce conservation constraints on this projection. We derive a mechanism for enforcement based on weak equality and implement it in the PETSc libraries. Scalability is demonstrated to more than one billion particles.

Utilization of parallel computing for mathematical modeling of high-enthalpy flows

The complexity of problems in gas dynamics is increasing from year to year, which makes sequential computational codes more and more inadequate and brings the importance of creating their parallel implementations. The problem of building the parallel implementations is not simple, especially in case of parallelization of the existing sequential codes. In this paper we present a method of simulation of high-enthalpy flows and describe the peculiarities of its parallel implementation utilizing PETSc library, which gives a wide variety of tools to simplify creation of parallel codes. The results for calculation of OREX reentry module as well as analysis of scalability of the parallel code implementation are given.

Optimization of Sparse Distributed Computations

We address the problem of the optimization of sparse matrix-vector product (SpMV) on homogeneous distributed systems. For this purpose, we propose three approaches based on partitioning the matrix into row blocks. These blocks are defined by a set of a fixed number of rows and a set of contiguous (resp. non-contiguous) rows containing a fixed number of non-zero elements. These approaches lead to solve some specific NP-hard scheduling problems. Thus, adequate heuristics are designed. We analyse the theoretical performance of the proposed approaches and validate them by a series of experiments. This work represents an important step in an overall objective which is to determine the best-balanced distribution for the SpMV computation on a distributed system. In order to validate our approaches for sparse matrix distribution, we compare them to hypergraph model as well as to PETSc library for SpMV distribution on a homogenous multicore cluster. Experimentations show that our approaches provide performances 2 times better than hypergraph and 49 times better than PETSc.

A Robust Algebraic Domain Decomposition Preconditioner for Sparse Normal Equations

Solving the normal equations corresponding to large sparse linear least-squares problems is an important and challenging problem. For very large problems, an iterative solver is needed and, in general, a preconditioner is required to achieve good convergence. In recent years, a number of preconditioners have been proposed. These are largely serial, and reported results demonstrate that none of the commonly used preconditioners for the normal equations matrix is capable of solving all sparse least-squares problems. Our interest is thus in designing new preconditioners for the normal equations that are efficient and robust and can be implemented in parallel. Our proposed preconditioners can be constructed efficiently and algebraically without any knowledge of the problem and without any assumption on the least-squares matrix except that it is sparse. We exploit the structure of the symmetric positive definite normal equations matrix and use the concept of algebraic local symmetric positive semidefinite splittings to introduce two-level Schwarz preconditioners for least-squares problems. The condition number of the preconditioned normal equations matrix is shown to be theoretically bounded independently of the number of subdomains in the splitting. This upper bound can be adjusted using a single parameter $\tau$ that the user can specify. We discuss how the new preconditioners can be implemented on top of the PETSc library using only 150 lines of Fortran, C, or Python code. Problems arising from practical applications are used to compare the performance of the proposed new preconditioner with that of other preconditioners.

A mass conserving arbitrary Lagrangian–Eulerian formulation for three‐dimensional multiphase fluid flows

AbstractThe arbitrary Lagrangian–Eulerian (ALE) framework presented in [Sahin and Guventurk, An arbitrary Lagrangian–Eulerian framework with exact mass conservation for the numerical simulation of 2D rising bubble problem.Int J Numer Methods Eng. 2017;112:2110–2134] has been extended to simulate three‐dimensional immiscible incompressible multiphase fluid flows. The div‐stable side centered unstructured finite volume formulation is used for the discretization of the incompressible Navier–Stokes equations. A novel kinematic boundary condition that satisfies the local discrete geometric conservation law (DGCL) is implemented to ensure the mass conservation of each species at the machine precision. The pressure field is treated to be discontinuous across the interface with the discontinuous treatment of density and viscosity. Surface tension force is treated as a tangent force and discretized in a semi‐implicit form. Two different approaches for the computation of unit normal vector have been implemented: the least squares biquadratic surface fitting (LSBSF) and the mean weighted by sine and edge length reciprocals (MWSELR). The resulting large system of algebraic equations is solved in a fully coupled manner in order to improve the time step restrictions. As a preconditioner, an approximate matrix factorization similar to that of the projection method is employed and the parallel algebraic multigrid solver BoomerAMG provided by the HYPRE library, which is accessed through the PETSc library, has been utilized for the scaled discrete Laplacian of pressure and the diagonal blocks of mesh deformation equations. The accuracy of the proposed method is initially validated on the static bubble problem, since the surface tension force is highly sensitive to the accurate evaluation of the unit normal vector and the inaccuracies significantly contribute to unphysical velocities, called parasitic currents. The calculations indicate that the parasitic currents can be reduced to machine precision for the MWSELR method. Then, the proposed approach is applied to the single bubble rising in a viscous quiescent liquid for both low and high density ratios. The calculations produce accurate predictions of the bubble shape, center of mass, rise velocity, and so forth. Furthermore, the mass of each species is conserved at machine precision and discontinuous pressure field is obtained in order to avoid errors due to the incompressibility restriction in the vicinity of liquid‐liquid interfaces at large density and viscosity ratios.

Efficient Simulations of Large Scale Convective Heat Transfer Problems

We describe an approach for efficient solution of large scale convective heat transfer problems, formulated as coupled unsteady heat conduction and incompressible fluid flow equations. The original problem is discretized in time using classical implicit methods, while stabilized finite elements are used for space discretization. The algorithm employed for the discretization of the fluid flow problem uses Picard's iterations to solve the arising nonlinear equations. Both problems, heat transfer and Navier-Stokes quations, give rise to large sparse systems of linear equations. The systems are solved using iterative GMRES solver with suitable preconditioning. For the incompressible flow equations we employ a special preconditioner based on algebraic multigrid (AMG) technique. The paper presents algorithmic and implementation details of the solution procedure, which is suitably tuned, especially for ill conditioned systems arising from discretizations of incompressible Navier-Stokes equations. We describe parallel implementation of the solver using MPI and elements of PETSC library. The scalability of the solver is favourably compared with other methods such as direct solvers and standard GMRES method with ILU preconditioning.

An efficient edge based data structure for the compressible Reynolds‐averaged Navier–Stokes equations on hybrid unstructured meshes

AbstractAn efficient edge based data structure has been developed in order to implement an unstructured vertex based finite volume algorithm for the Reynolds‐averaged Navier–Stokes equations on hybrid meshes. In the present approach, the data structure is tailored to meet the requirements of the vertex based algorithm by considering data access patterns and cache efficiency. The required data are packed and allocated in a way that they are close to each other in the physical memory. Therefore, the proposed data structure increases cache performance and improves computation time. As a result, the explicit flow solver indicates a significant speed up compared to other open‐source solvers in terms of CPU time. A fully implicit version has also been implemented based on the PETSc library in order to improve the robustness of the algorithm. The resulting algebraic equations due to the compressible Navier–Stokes and the one equation Spalart–Allmaras turbulence equations are solved in a monolithic manner using the restricted additive Schwarz preconditioner combined with the FGMRES Krylov subspace algorithm. In order to further improve the computational accuracy, the multiscale metric based anisotropic mesh refinement library PyAMG is used for mesh adaptation. The numerical algorithm is validated for the classical benchmark problems such as the transonic turbulent flow around a supercritical RAE2822 airfoil and DLR‐F6 wing‐body‐nacelle‐pylon configuration. The efficiency of the data structure is demonstrated by achieving up to an order of magnitude speed up in CPU times.

GPU Preconditioning for Block Linear Systems Using Block Incomplete Sparse Approximate Inverses

Solving sparse triangular systems is the building block for incomplete LU- (ILU-) based preconditioning, but parallel algorithms, such as the level-scheduling scheme, are sometimes limited by available parallelism extracted from the sparsity pattern. In this study, the block version of the incomplete sparse approximate inverses (ISAI) algorithm is studied, and the block-ISAI is considered for preconditioning by proposing an efficient algorithm and implementation on graphical processing unit (GPU) accelerators. Performance comparisons are carried out between the proposed algorithm and serial and parallel block triangular solvers from PETSc and cuSPARSE libraries. The experimental results show that GMRES (30) with the proposed block-ISAI preconditioning achieves accelerations 1.4 × –6.9 × speedups over that using the cuSPARSE library on NVIDIA Tesla V100 GPU.

A new block preconditioner and improved finite element solver of Poisson-Nernst-Planck equation

The Poisson-Nernst-Planck (PNP) equation is one important continuum model for studying charge transport in ion channels, which is a common phenomenon and plays a key role in molecular biosciences. In this paper, to improve the current PNP solvers, according to the equation structure, a new block preconditioner was proposed and proved that the preconditioned linear system has bounded eigenvalues independent of mesh sizes under some conditions, thus guaranteeing that the convergence rate of the preconditioned linear solver is independent of mesh sizes. Meanwhile, the commonly-used solution decomposition schemes in classic continuum models for isolating the singularities were presented and then one was chosen for solving PNP when zero initial guess was used for the Newton method. Furthermore, an efficient and improved finite element PNP solver was proposed by further combining the two-grid method as an acceleration technique. Then the new program package was fulfilled based on the state-of-the-art finite element library FEniCS and the efficient scientific library PETSc. Finally, numerical simulations on a test model with analytical solution as well as some tests on protein cases were carried out to validate the new program package and our theoretical results.

Fully Parallel Mesh I/O Using PETSc DMPlex with an Application to Waveform Modeling

Large-scale PDE simulations using high-order finite-element methods on unstructured meshes are an indispensable tool in science and engineering. The widely used open-source PETSc library offers an efficient representation of generic unstructured meshes within its DMPlex module. This paper details our recent implementation of parallel mesh reading and topological interpolation (computation of edges and faces from a cell-vertex mesh) into DMPlex. We apply these developments to seismic wave propagation scenarios on Mars as an example application. The principal motivation is to overcome single-node memory limits and reach mesh sizes which were impossible before. Moreover, we demonstrate that scalability of I/O and topological interpolation goes beyond 12,000 cores, and memory-imposed limits on mesh size vanish.

SOL-KiT—Fully implicit code for kinetic simulation of parallel electron transport in the tokamak Scrape-Off Layer

Here we present a new code for modelling electron kinetics in the tokamak Scrape-Off Layer (SOL). SOL-KiT (Scrape-Off Layer Kinetic Transport) is a fully implicit 1D code with kinetic (or fluid) electrons, fluid (or stationary) ions, and diffusive neutrals. The code is designed for fundamental exploration of non-local physics in the SOL and utilizes an arbitrary degree Legendre polynomial decomposition of the electron distribution function, treating both electron–ion and electron–atom collisions. We present a novel method for ensuring particle and energy conservation in inelastic and superelastic collisions, as well as the first full treatment of the logical boundary condition in the Legendre polynomial formalism. To our knowledge, SOL-KiT is the first fully implicit arbitrary degree harmonic kinetic code, offering a conservative and self-consistent approach to fluid–kinetic comparison with its integrated fluid electron mode. In this paper we give the model equations and their discretizations, as well as showing the results of a number of verification/benchmarking simulations. Program summaryProgram Title: SOL-KiTCPC Library link to program files:http://dx.doi.org/10.17632/7p44sdfnnp.1Licencing provisions: GNU GPLv3Programming language: Fortran 90Nature of problem:Fluid models of parallel transport in the Scrape-Off Layer (SOL) fail to account for the fact that the electron distribution function is often far from a Maxwellian, and kinetic effects have been linked to discrepancies between experiment and fluid modelling[1]. A kinetic treatment of electrons in the SOL requires detailed accounting of collisional processes, especially those with neutral particles, as well as a proper implementation of the logical boundary condition at the material surface[2]. Furthermore, the ability to identify differences between fluid and kinetic modelling using self-consistent comparison is desirable.Solution method: Electrons are modelled either as a fluid, or kinetically, maintaining self-consistency between models. All equations are solved using finite difference and the implicit Euler method, with fixed-point iteration. The kinetic approach is based on solving the Vlasov–Fokker–Planck–Boltzmann equation, decomposed in Legendre polynomials. Equations for the harmonics (or electron fluid equations) are solved alongside fluid equations for the ions, Ampère–Maxwell’s law for the electric field, as well as a diffusive–reactive Collisional-Radiative model for the evolution of hydrogenic atomic states. Each individual operator is built into a matrix, combined with all other operators, and the matrix equation arising from the Euler method is solved using the PETSc library, with MPI parallelization.Additional comments:This article presents the physical and numerical outline of the code, and is accompanied by html documentation, a small test suite based on benchmarking runs presented below, as well as instructions and means for compiling and executing SOL-KiT. Special focus in the article is given to the novel numerical and model aspects in the greater context of the developed software.

PETSc DMNetwork

We present DMNetwork, a high-level package included in the PETSc library for the simulation of multiphysics phenomena over large-scale networked systems. The library aims at applications that have networked structures such as those in electrical, gas, and water distribution systems. DMNetwork provides data and topology management, parallelization for multiphysics systems over a network, and hierarchical and composable solvers to exploit the problem structure. DMNetwork eases the simulation development cycle by providing the necessary infrastructure through simple abstractions to define and query the network components. This article presents the design of DMNetwork, illustrates its user interface, and demonstrates its ability to solve multiphysics systems, such as an electric circuit, a network of power grid and water subnetworks, and transient hydraulic systems over large networks with more than 2 billion variables on extreme-scale computers using up to 30,000 processors.