extfieldsofstudy:
Computer Science
Physics
Mathematics

text:
Dynamical p-adaptivity for LES of compressible flows in a high order DG framework

We investigate the possibility of reducing the computational burden of LES models by employing locally and dynamically adaptive polynomial degrees in the framework of a high order DG method. A degree adaptation technique especially featured to be effective for LES applications, that was previously developed by the authors and tested in the statically adaptive case, is applied here in a dynamically adaptive fashion. Two significant benchmarks are considered, comparing the results of adaptive and non adaptive simulations. The proposed dynamically adaptive approach allows for a significant reduction of the computational cost of representative LES computation, while allowing to maintain the level of accuracy guaranteed by LES carried out with constant, maximum polynomial degree values.

Introduction
During the last decades, Large Eddy Simulation (LES) has been extensively studied and has allowed to achieve results comparable to those of Direct Numerical Simulations (DNS) with a significant reduction in computational cost, enabling the use of high accuracy models to simulate turbulent flows beyond the extremely idealized configurations to which DNS was typically limited. While the increase of the available computing power makes even more fluid flow configurations amenable to LES approaches, two interlinked difficulties prevent LES from being a widespread tool for industrial applications. The first major issue is that LES simulations, due to their space and time resolution requirements, are still very computationally expensive. The second is that selecting the right spatial resolution for an LES of a unsteady turbulent flow in a complex geometry is not a easy task. On one hand, the resolution must be sufficient to solve the equations with adequate accuracy, but limited to contain the computational costs. On the other hand, the implicit LES (ILES) commonly employed uses the discretization as the filter, so that the mesh resolution sets implicitly the filter width and the threshold between the resolved scales and the modelled ones. Moreover, in a turbulent flow in complex geometry, turbulence characteristics change in time and in space throughout the computational domain and cannot be as easily estimated as it was the case in the simple geometries employed for many years in LES studies. Therefore, an adaptive LES, where the resolution and the filter width are automatically adapted to the flow conditions, would be of great benefit for the solution of both problems, in order to have an LES with the correct filter sizing and an efficient spatial discretization. The problem of estimating the right resolution is further complicated by the fact that, as pointed out by the authors in [34], using standard error estimators to adapt the resolution is not the optimal approach for LES, if refinement to DNS resolutions is to be avoided, so that specific strategies must be devised in this context.
The need of an adaptive LES was first stated by [26], but not many attempts have been made in this direction since then. In the numerical solution of PDEs, however, in the last decades considerable efforts have been devoted to the development of refinement indicators. Most of these are based on local estimates of the discretization error, see e.g. [4,5,20,28,29,37,36,35,39]. In the LES context, however, simply increasing the resolution in order to decrease the error leads to a DNS solution [23,29], which is in contrast with the goal of adaptive LES, which consists in adjusting the resolution in order to directly resolve only a prescribed amount of the turbulent scales. An alternative approach, pursued e.g. in [17,19,18], employs indicators aimed at obtaining a good representation of some quantity of interest, rather than of the physical scales of the turbulence. Adjoint-based adaptation methods proposed by other authors [11] have also shown good performance, but entail significant computational costs.
In the present work, we will consider a physically based refinement indicator proposed by the authors in [34,31], which is especially suited for LES. This adaptivity approach is implemented in a numerical framework based on the Discontinuous Galerkin (DG) method, that was already presented and validated in [1], which allows to extend seamlessly LES concept to unstructured meshes. DG methods combine the high order of accuracy and low dissipation/dispersion properties with good parallel performances, since most of the computations are local to the element and only interelement boundary fluxes must be communicated. The adaptation technique we employ is a p−adaptive approach, in which the polynomial degree is varied locally, rather than the mesh size, as in more standard h−adaptive methods. DG methods provide an interesting environment for p−adaptivity, since they do not require to enforce continuity constraints at the interelement boundaries. Furthermore, p−adaptive techniques are appealing since they allow to correct possible shortcomings of the computational mesh, as well as to perform dynamically adaptive simulations without extensive remeshing. Static polynomial adaptivity was applied to LES in [10,34,24,31] to efficiently simulate statistically steady phenomena, while an example of dynamical adaptation is described in [11] in an ILES context.
In this paper, we will show how the adaptive method of [34,31] can be applied locally and dynamically in time, in order to simulate efficiently and accurately transient phenomena. A number of time dependent numerical simulations of two different configurations is considered, namely a temporally evolving mixing layer and the interaction of a vortex impinging on a square cylinder with the cylinder wake. The results obtained show that a significant reduction of the computational cost of LES can be achieved, reducing the number of degrees of freedom employed to values of up to 50% of those required by constant maximum polynomial degree simulations, while maintaining essentially the same level of accuracy.
In Section 2, the model equations and numerical method are reviewed. In Section 3, the dynamical adaptivity approach we propose is outlined. In Sections 4, 5, results of numerical simulations on two significant benchmarks are presented. Finally, in Section 6 some conclusions are drawn and some perspectives of future development of this approach are discussed.

Model equations and numerical method
In this Section, the model equations and the numerical method will be briefly illustrated. For a more detailed description, we refer to [1], where this methodology was first introduced and validated for LES applications. The LES filtered Navier-Stokes equations for compressible flows can be written in dimensionless form as: where U = [ρ, ρ u T , ρ e] T are the prognostic resolved variables and F = F ca − F v + F sgs are the fluxes, composed by the advective fluxes the viscous fluxes and finally the subgrid fluxes Here, · and · represent the grid filter and the Favre filter operators, so that p denotes the filtered pressure, ρ h = ρ e+p the resolved enthalpy and σ and q the momentum and heat diffusive fluxes, respectively. Equations (1) are completed by the dimensionless state equation for a perfect gas where ‹ T is the filtered temperature. To close the system the constitutive equations must also be specified: where the rate of resolved strain tensor is defined as and the dynamic viscosity, according to a power law, is with α = 0.7. In the Local Discontinuous Galerkin approach [8], the equations are rewritten as a first order system where ϕ = [ u T , ‹ T ] T . A tessellation T h composed of non overlapping tetrahedral elements is defined in the domain Ω over which a discontinuous finite element space V h is defined in the following way where P p K (K) denotes the space of polynomial functions of total degree p K . The weak discrete form of equation (9) is derived using test functions belonging to the same space as the numerical solution. A modal DG approach is here applied, by using a hierarchical orthonormal polynomial basis for each element K in the finite dimensional space V h to represent the numerical approximation of the generic variable a expressed as where φ K l are the basis functions on element K, a (l) are the modal coefficients of the basis functions and n φ (K) + 1 is the number of basis functions required to span the polynomial space P p K (K) of degree p K , defined in R 3 as: To compute the advective flux {F} a the Rusanov flux [21] is employed while a centred flux is used for the diffusive fluxes {F} v , {F} sgs and for the gradient variables flux {ϕ}.
In the present framework, as explained in detail in [1], the filtering operators are built in the DG spatial discretization, so that the LES filtering · is equivalent to the projection onto the employed finite dimensional solution subspace while Favre filter operator is defined as In the simulations presented in this paper, two different subgrid models have been employed, the classical Smagorinsky model, in which the filter width is determined as and the anisotropic dynamic model [1,3], in which the filter width is included in the coefficients to be dynamically determined [2] using the Germano identity [15].

A dynamically p−adaptive approach
A physically based refinement indicator especially suited for LES has been proposed by the authors in [34,31]. This indicator is based on the classical structure function where · represents the expected value operator. The structure function has been widely used to study turbulence statistics and it is known to be directly related to subgrid stresses [6,7,14]. Larger values of the structure function calculated inside the element denote a poorly correlated velocity field, and the need of a higher resolution, while lower values denote a highly correlated velocity field, which is an indication of a very well resolved turbulent region or laminar conditions, hence the possibility of employing a lower resolution. The degree adaptation indicator is then defined as: where D iso ij is the structure function in isotropic conditions [25]. In [34], this procedure was shown to be effective in a statically p−adaptive framework, producing accurate results with a significant reduction in the computational cost. However, static adaptivity presents some limitations, since the resolution is fixed at the beginning of the simulation and constant in time. In the simulation of a transient phenomenon, for which a time resolved solution is sought rather than a statistical average, employing a dynamic adaptivity framework is necessary. In our implementation of the dynamic degree adaptation, indicator (18) is computed at runtime at time intervals ∆t i and averaged over time. Once the indicator has been averaged on a time interval ∆t a which is sufficiently larger than ∆t i , but still small with respect to the time scale of the motion of the main turbulent structures, a new polynomial distribution is computed based on the indicator values and the flow field approximation is updated, by either reducing or increasing the number of degrees of freedom employed in each element.
In this work, the admissible polynomial degrees ranged from 2 to 4. Two indicator thresholds 1 < 2 were employed. Cells with a value of the indicator (18) lower than the smaller threshold were assigned the lowest allowed polynomial degree 2, those with value higher than the largest threshold were assigned the highest polynomial degree 4 and the intermediate ones were assigned degree 3. Polynomial degrees were only allowed to be increased or decreased by one degree at each adaptation time. Since the solution is represented by a hierarchical basis, when lowering the polynomial degree, the contribution associated to the removed modes is simply discarded, while when raising the polynomial degree the coefficients of the newly added modes are initialized to zero, to be updated by the subsequent time evolution.
It should be remarked that the procedure outlined above can potentially lead to unbalances between the computational load of different processors in a parallel run. The implementation of a full dynamic load balancing is beyond the scope of the present work, whose focus is mostly on the assessment of the accuracy of the above methodology. As a consequence, in all the numerical experiments presented below, the computational load among the processors is approximately balanced at the beginning of the simulation only, which is clearly suboptimal with respect to full parallel efficiency. Indeed, the results will not in general be assessed based of the CPU time required, but rather on the total number of the degrees of freedom employed in the adaptive versus non adaptive simulations. However, it must be noted that,t even with this suboptimal configuration, the adaptive simulation always led to a net reduction in computational effort. In order to assess the effective reduction in the CPU time required, further work on code optimization is needed, along with the development and application of load balancing approaches such as those employed e.g. in [11,40]. 8

Temporally evolving mixing layer
In a first assessment of the adaptive LES approach outlined in the previous sections, an isothermal time developing mixing layer was simulated. The temporal mixing layer represents an interesting test case for LES, due to its simplicity and advantages from the computational viewpoint, as well as the complexity of the physics related to the mixing, see e.g. the discussions in [13,30,38,41]. The commonly employed configuration to represent a mixing layer flow was employed, imposing periodic boundary conditions on the mean flow and in the span wise directions. Thanks to the periodicity in the stream wise direction, uncertainties related to the imposition of inlet and outlet boundary conditions are avoided. In the simulations presented here, the anisotropic dynamic subgrid model [1,3] has been used, in which the filter width is included in the coefficients to be dynamically determined. A further application of this approach is presented in [2].
A sketch of the flow configuration is shown in Figure 1. The mixing layer is characterized by two parallel flows with different velocities U 1 and U 2 . The convection velocity of the isothermal mixing layer, defined as U c = (U 1 + U 2 )/2, is the velocity that transports the cortical structures at the centreline. In the configuration we considered, the convection velocity is assumed to be zero, so that the eddies do not travel inside the domain by means of the flow, but only move by mutual interaction. As a consequence, U 1 = −U 2 and |U 1 | = |U 2 | = (U 1 − U 2 )/2. The reference frame is chosen so that the x axis is aligned with the main direction of the flow centred at the middle of the domain, while the y and z axis denote transverse and lateral directions, respectively. The initial vorticity thickness, defined as is chosen as reference length for the non-dimensionalization. The initial velocity jump ∆U r = (U 1 − U 2 ) is used as reference velocity, the initial uniform density ρ r and temperature T r are used as reference density and temperature, respectively. The fundamental dimensionless groups are then The temporally evolving mixing layer develops from a specified initial condition. In these computations, we use a hyperbolic tangent where U 1 = 0.5 and U 2 = −0.5 so that U c = 0. Similarly to [16,12], a 3D incompressible disturbance is added to the base velocity profile to initiate the transition process. It consists of harmonic disturbances expressed as where A = 0.025(U 1 −U 2 ) and the decay in y-direction is governed by σ = 0.05. Moreover, a 3D white noise perturbation is added to the initial velocity field, generated for each mesh node by a logistic map rescaled to take values on [−A, A] and multiplied by the attenuation factor e −σy 2 . Uniform pressure p = 1 and density ρ = 1 are used to initialize the simulation. The longitudinal dimension of the computational domain must be large enough to allow at least for the merging of two principal vortical structure, so it must be taken as a multiple of the wavelength λ a characterizing the most unstable perturbation, which according to the analysis in [30] is given in this context by λ a = 7.66. In the present simulation, the size in the streamwise direction has been chosen L x = 8×λ a = 60.8, so that up to three vortices are allowed to merge. In the normal direction, the computational domain extends for −30 ≤ y ≤ 30, while in spanwise direction L z = 30.4. Periodic boundary conditions are applied in the statistically homogeneous x and z directions, while a sponge layer is employed at the top and bottom boundaries, so that the effective size of the computational domain in the y direction is equal to 80.
The computational mesh is built first subdividing the domain in N hx × N hy × N hz = 20 × 36 × 10 hexahedra. Each hexahedron is then split into 6 tetrahedra, yielding a total of 43200 tetrahedral elements. The actual resolution depends on the local polynomial order and can be computed as For p = 4, we get ∆ x = 5.11, ∆ y = 0.37, ∆ z = 0.5.
A preliminary study of the mixing layer in the present configuration was presented in [27], where the sensitivity of the results to the choice of spatial resolution and subgrid scale model was checked. The objective of the present work is to test the p−adaptivity technique in a time evolving flow. For this purpose, the results obtained with polynomial degree dynamically adapting in space and time are compared with the results obtained with constant fourth and third degree polynomials.
To set the values of the adaptation thresholds 1 , 2 , the structure function indicator was evaluated on the velocity field obtained with constant p = 4 at the final dimensionless time t = 150. Thresholds were then chosen in order to obtain, for that field, an adapted number of degrees of freedom (dofs) lower than 80% of the corresponding number of dofs in the constant p = 3 case. The resulting threshold values are then 1 = 10 −6 and 2 = 2 × 10 −3 . The dynamical p−adaptation was carried out by computing the indicator value at intervals ∆t i = 0.0375 and averaging them over intervals ∆t a = 0.075, while the value ∆t = 0.0075 was used for time integration.
The polynomial degree distribution, the density and the streamwise component of the velocity, at different instants, in the plane z = 0, are represented in Figure 2. In the initial condition, the elements with highest degree p = 4 are concentrated along the middle of the domain, where high resolution is required. Then, the region with higher polynomial degree extends in space during time, adapting to the merging of the vortices, to the diffusion of turbulent structures and to the growth of the layer thickness. An important quantity for  the characterization of the mixing layer is the momentum thickness, defined as The time evolution of this quantity is represented in Figure 3. Notice that slope variations in the time evolution are associated to the merging of the vortices. A difference in the behaviour of the p = 3 simulation with respect to the p = 4 one is evident after the second merging, while the p−adaptive simulations well reproduces the p = 4 results. No relevant differences between the results obtained using adaptive or constant polynomial degree distribution are observed in the mean density and velocity profiles, shown in Figure 4, nor in the mean normal stresses profiles, shown in Figures 5(a)-5(c). Instead, looking at the shear stresses in Figure 5(d), although the peak is slightly overestimated, the p−adaptive solution is in good agreement with that of the p = 4 case. From the same picture, the more dissipative character of the p = 3 solution is also apparent.
Looking at the turbulent kinetic energy profile in Figure 6(a), we can affirm that the turbulent energy amount of the flow is reproduced by the p−adaptive simulation with the same accuracy as in the p = 4 case. In Figure 6(b), the subgrid scale (SGS) kinetic energy in the p = 3 case is clearly larger than in the other cases, due to the lower resolution that requires more intense contribution from the model. The SGS energy for the p−adaptive case is closer to that of the p = 4 case, demonstrating that the two simulations have an equivalent effective resolution.
All the simulations were carried out on the Marconi cluster at     Core hours p = 4 12600 p = 3 9000 p−adaptive 9300 Table 1: Core hours required for the mixing layer simulations using uniform polynomial degree distribution p = 4, p = 3 and p−adaptive, respectively.
CINECA, using 272 KNL processors. In Figure 7, the total number of dofs for the adaptive simulation is shown as a function of time and compared with the values corresponding to uniform polynomial distributions, that are about 1.51 × 10 6 and 8.64 × 10 5 for p = 4 and p = 3, respectively. The gain using p−adaptivity is evident, since it requires less than 80% of the dofs for p = 3 and less then half of the dofs for the p = 4 case. The core hours necessary to perform the simulations are reported in Table 1. In spite of the absence of a specific dynamic load balancing procedure, the p−adaptive technique also yields a reduction of the required CPU time, which has a value only marginally larger than that of the p = 3 simulation. On the other hand, the sub-optimal nature of the present implementation is highlighted by the fact that the p−adaptive simulation required more CPU time than the p = 3 case, even though it involved less dofs.
Concluding this analysis, we can affirm that the p−adaptive solutions display almost the same accuracy as those obtained with constant degree four, but require less then the half the number of dofs and about 29% less CPU time for this test case, in spite of a sub-optimal implementation. Moreover, the structure function indicator which was proven successful in static adaptation in [34] is shown to be suitable also for simulation of transitional flows using a dynamic adaptation in time.

Interaction of Vortex and Square Section Cylinder
A body-vortex interaction flow represents another interesting time dependent problem to test dynamic p−adaptivity. Here, a vortex interacting with a square section cylinder at Re = 22000 and M a = 0.3 is considered. This configuration was chosen because the analogous case without a vortex has already been studied with statically adaptive techniques in [34]. The classical Smagorinsky SGS model was applied for this test. The results obtained with polynomial adaptivity are compared with the solution obtained with uniform in space and constant in time degree distribution. The constant degree four solution is taken as reference.
The computational domain employed is box shaped and a 2D sketch of the geometry is represented in Figure 8. Denoting the cylinder side with H, which is used as reference length, the inflow length has been taken equal to L f = 10H, while the outflow length is equal to L r = 20H. The cylinder is vertically centred with a distance of L s = 10H from the upper and lower boundaries. The domain is extruded in the spanwise direction by a length equal to L z = 4H. This configuration was employed e.g. in [9] and has been used to assess statically adaptive LES simulations of the standard flow around a square section cylinder in [34,31,32]. The mesh employed consists of 23816 tetrahedra arranged in an outer unstructured area and a structured O-grid mesh around the cylinder, which is then extruded in the spanwise direction. Considering a polynomial basis of degree 4, the resolution in space around the cylinder is ∆ n = 0.008H in the wall normal direction, ∆ t = 0.024H along the cylinder sides and ∆ z = 0.135 along the spanwise direction. As in [34], the thresholds for the structure function (SF) indicator are taken to be 1 = 5×10 −4 , 2 = 1×10 −2 , so as to obtain a number of dofs just lower than the one obtained with uniform degree 3.
To simulate the cylinder-vortex interaction, a vortex is superimposed to an initial condition corresponding to a developed flow around the cylinder, obtained by a previous LES simulation. The introduced vortex is parallel to the z direction and defined by the velocity components where C v denotes the vortex strength, R v is the vortex radius, r 2 = x 2 + y 2 is the distance in the x − y plane from the vortex centre and U ∞ is the uniform flow over which the vortex is superimposed [22]. The maximum radial velocity is obtained at r = R v and is given by v θmax = Cv R √ e . Regarding the other variables, it is known from [9] that for a viscous compressible vortex the radial pressure distribution given by the solution of Assuming a constant temperature distribution T = T ∞ , one obtains from (24) the profile Finally the density distribution follows the equation of state ρ = P/(R gas T ).
In [31,33] simple advection of this vortex was tested. The maximum radial velocity chosen is v θmax = 0.5 and the vortex radius is R v = 0.41. A reference simulation without vortex was used to calibrate the introduction time in order for the vortex to reach the cylinder in the instant of maximum lift. The dynamical p−adaptation strategy consisted in this case of adapting the polynomial degrees approximately three times during the time needed by a fluid particle to pass through the smallest elements on the cylinder walls at the free stream velocity U ∞ . This choice led to computing the indicator value at intervals ∆t i = 0.0016 and averaging them over intervals ∆t a = 0.05, while the value ∆t = 0.0004 was used for time integration.
The dynamically adaptive procedure was able also in this case to effectively represent the structures of the flow, both in the advected vortex region and in the wake of the obstacle, as can be seen in Figures 10-13. It is possible to note that a higher polynomial degree is employed in the advected vortex region and in the shear layers around the cylinder. The dynamic adaptation procedure is able to effectively increase the polynomial degree around the vortex and to follow it as it is advected downstream, leaving all the elements with no vortex activity at the lowest resolution. Notice that the centre of the vortex is represented with degree 4 polynomials in the first instants, but only degree 3 polynomials are employed later, since the vortex enters an area with a more refined mesh in which a lower polynomial degree is sufficient.
Among the different effects caused by the vortex-cylinder interaction, the most interesting are those on the forces acting on the cylinder. The history of the load coefficients before, during and after the interaction is reported in Figure 14.
In the developed flow around the square cylinder, the periodic oscillation of the force coefficients are related to the unsteadiness of the recirculation bubble on the upper and lower sides and to the vortex shedding. Reaching the upper corner of the cylinder at t = 5, the vortex changes the pressure distribution on the front of the cylinder and strongly modifies the recirculating bubble on the upper side. In this condition, it is very important to use an adequate resolution on the upper side and in the separation bubble. In Figure  14 it can be observed that the peaks of the lift and drag coefficients at t = 5 with p−adaptivity equal those in the reference solution p = 4, demonstrating the capability of the proposed dynamically adaptive approach to well represent this unsteady phenomenon. The values of the differences of the peaks for the lift and drag coefficients with respect to the vortex free solution are reported in Table 2. Also in this case the values obtained with p−adaptivity are in very good agreement with the reference solution.      configuration Table 2: Differences, in percentage, of the force coefficients with respect to the simulation without vortex, during the vortex interaction time, t = 4 − 8.5. For the lift, only the change in the subsequent peak is presented, for the drag both the difference in maximum and minimum.
The simulations were carried out on the Marconi cluster at CI-NECA, using 272 Brodwell processors. The adaptive simulation allows a reduction of about 60% of the required CPU time with respect to the constant p = 4 degree case. In Figure 15, the time history of the dofs number is shown. It is worth noting that, even though the number of dofs is affected by the flow regime inside the domain, its variations are anyway fairly limited and the total number remains less than the one employed for the uniform p = 3 simulation. The core hours necessary to perform the simulations are reported in Table 3. Also in this case,the p−adaptive technique yields a reduction of the required CPU time, which has a value only marginally larger than that of the p = 3 simulation. Again, the sub-optimal nature of the present implementation is highlighted by the fact that the p−adaptive simulation required more CPU time than the p = 3 case, even though it involved less dofs.
Core hours p = 4 756 p = 3 388 p−adaptive 462 Table 3: Core hours required for the vortex-cylinder simulations using uniform polynomial degree distribution p = 4, p = 3 and p−adaptive, respectively.

Conclusion
The need for adaptive LES approaches was first stated in [26], but adaptive LES simulations are not very common yet. DG methods provide a appropriate environment for adaptive approaches, since they do not require to enforce continuity constraints at the interelement boundaries. Furthermore, p−adaptive techniques are appealing, since they allow to correct possible shortcomings of the computational mesh as well as to perform dynamically adaptive simulations without extensive remeshing. Static polynomial adaptivity has been applied to LES in a DG context in [10,34,24,31]. An essential tool for such simulations is an adaptation criterion that, rather than simply increasing the resolution in order to decrease the error, which is known to lead to a DNS solution [23,29], tries instead to adjusting the resolution in order to directly resolve only a prescribed amount of the turbulent scales.
In this work, the physically based refinement indicator proposed by the authors in [34,31] was shown to be applicable locally and dynamically in time, in order to simulate efficiently and accurately transient phenomena. Numerical simulations of a temporally evolving mixing layer and of a vortex impinging on a square cylinder were performed. The results obtained show that a significant reduction of the computational cost of LES is achieved, reducing the number of degrees of freedom employed to values down to 50% of those required by constant maximum polynomial degree simulations, while maintaining essentially the same level of accuracy.
A main limitation of the present results is that, in this preliminary implementation of a dynamically adaptive approach, no dynamic load balancing has been performed. In future work on code optimization, load balancing approaches such as those employed in [11,40] will be considered and tested in order to achieve maximum parallel efficiency.


s2fieldsofstudy:
Computer Science
Engineering

text:
Accelerated Sparse Matrix-Based Computation of Electromagnetic Transients

This paper is related to research on parallelization methods for the simulation of electromagnetic transients (EMTs). It presents an automatic parallelization approach based on the solution of sparse matrices resulting from the formulation of network equations. Modified-augmented-nodal analysis is used to formulate network equations. The selected sparse matrix solver is parallelized and adapted to improve performance by pivot validity testing and partial refactorization. Refactorization is needed when dealing with varying topology networks and nonlinear models. The EMT solver employs a fully iterative method for nonlinear functions. Conventional computer CPU-based parallelization is achieved and does not require any user intervention for given arbitrary network topologies. The presented approach is tested on real networks with complex models, including nonlinearities and power-electronics converters for wind generator applications.

I. INTRODUCTION
T HE computing time reduction for the simulation of electromagnetic transients [1], [2] (EMTs) is a crucial research topic. The EMT-type [2] simulation methods are circuit based and can use very accurate models for an extended frequency range of power system phenomena. This qualifies them as being of wideband type. In fact the EMT approach is applicable to both slower electromechanical transients and much faster electromagnetic transients. The computation of electromechanical transients can be achieved with EMT-type solvers for very large networks [3] and requires significant computing time when compared to phasor-domain approaches, but even for smaller networks, the computing time can become a key factor due to numerical integration time-step constraints or model complexity level. More and more challenging simulation cases are created for studying modern power systems, those include, for example, HVDC systems [4] and wind generation [4].
There are several techniques for improving computational performance in EMT-type solvers. Such techniques include improvements in model performance using, for example, average-value models [5] for power-electronics based systems or circuit reduction [6]. Network reduction can be also achieved using frequency domain fitting [7] or through dynamic equivalents [8]. Other approaches include usage of multiple time-steps [9], waveform relaxation [10] and combinations of different methods [11].
Although for some models the numerical performance improvements can be generalized, the above methods cannot be applied automatically and require user intervention and various other manual modifications. A more direct path towards computational speed improvement in EMT-type numerical methods is through efficient sparse matrix solvers and parallelization. The latter has been initiated in realtime simulation tools [12]- [14] and followed also in off-line [2], [15], [16] applications. The network tearing approach for parallelization is based on the exact separation created by distributed parameter transmission line/cable models. Other methods [17], [18] can be applied for cutting networks at arbitrary locations.
An important problem in network solution parallelization methods, such as [16], [17], is that user intervention is required for setting the network separation locations. In [18] the derivation of bordered-block-diagonal matrices (diakoptics) allows to cut networks without using transmission line delays. But such methods are demonstrated for linear and fixed topology networks that do not require time consuming matrix updating and refactorizations. Moreover, large network test cases are created by replicating a given small network a few hundred times. Such idealistic situations do not occur for actual real power grids.
This paper targets off-line simulation methods and presents CPU-based parallelization for conventional multi-core computers using a sparse matrix solver, named KLU [19], [20]. Such a solver has been already applied for circuit simulation problems [20]- [23], including the simulation of electromagnetic transients in [22], [23]. This paper contributes several significant improvements over [22], [23]. Fully automatic parallelization is again achieved through block-triangular factorization (BTF), but new modification of the KLU solver allow to significantly improve performance. The modifications include pivot validity testing and partial refactorization for a fully iterative sparse matrix solver. Refactorization is required to account for varying topologies (switching) and nonlinear models. The capability to efficiently account for iterations and refactorization in realistic networks is a distinctive aspect in the presented research.
Another improvement in this work is more efficient implementation of parallelization through distributed thread memory.
The proposed approach is implemented in EMTP [24] and tested on realistic large scale networks. The demonstration on such networks/benchmarks constitutes another important contribution of this paper. The solution process includes seamless transition from load-flow solution to timedomain parallel computations using modified-augmentednodal analysis.
The proposed approach can be implemented in any EMT-type simulation tool. Moreover, the proposed improvements in KLU are applicable to other types of solvers. This paper starts by reviewing the main structures of the EMT-type solver for which the sparse matrix code is being replaced. It follows by the presentation of the new parallel sparse matrix solution process and then its implementation. Test cases are presented to demonstrate the numerical performance advantages.

II. REVIEW OF EMTP STRUCTURE AND SOLVER
The new approach proposed in this paper is implemented using the modified-augmented-nodal analysis (MANA) formulation method [24]. This method offers several advantages [25] over classical nodal analysis. Its formulation is briefly recalled here to relate to material presented in the following sections. In MANA the system of equations is generic and can use different types of unknowns in addition to voltage. In MANA the system of network equations (NEs) is given by (1) This sparse system of equations is solved for the vector of unknowns x at each time-point in time-domain, after updating the matrix of coefficients A and the vector of known variables b. The vector b contains history terms and other independent functions, such as independent sources. The system (1) is non-symmetric and can also accommodate generic equations, such as where the terms on the left contribute coefficients (k j ) into the A matrix for voltage (v j ) and current (i j ) unknowns, and b z is a cell in the b vector. This equation allows to integrate directly source and switch equations. For an ideal switch equation, for example, In switch closed position, k = 1 and k s = 0. When the ideal switch is open, k = 0 and k s = 1. It is also possible to model non-ideal switches by setting k = 1 and replacing k s by high and low resistance values. Other models, such as ideal transformers with tap control can be easily accommodated [25]. The NEs (1) are solved at each time-point after updating A for switch position changes, transformer tap changes or any other modifications in model equations. For nonlinear models (NMs), the NEs must be solved iteratively to achieve an accurate simultaneous solution. This is done by linearizing each model [25] at each operating point change and solving iteratively (1) until convergence for all nonlinear models. Model linearization results into a Norton equivalent with the Norton resistance contributing changes into A and the Norton current contributing updates into the b vector. This is a Newton method with A being the Jacobian matrix.
For time-varying models (TVMs), such as switches or transformer tap positions, it is also possible to update A iteratively without advancing to the next time-point. This accuracy option, marked as iterative time-varying method (ITVM), allows to achieve a simultaneous solution for the determination of all changes and dependencies between models at the same time-point. This process also includes the sequential re-calculation of control system equations [25].
In EMTP the NE (1) are currently solved on a single CPUcore using a non-symmetric direct sparse solver method based on minimum degree ordering [26]. The permuted matrix A does not result into decoupled blocks and cannot be used for parallelization.
The fact that many models, such as the coefficients in (3), are time-varying and also for continuous changes that may be required by model updates in A during iterations, it is necessary to recalculate the LU factors of A. The refactorization process of the complete matrix A is time-consuming. The number of refactorizations is case dependent and has fluctuating impact on numerical performance. This is an important aspect. The upgrading of the sparse matrix solver for accelerating computations is described in the following sections.

III. BLOCK TRIANGULAR FACTORIZATION AND KLU
The Block Triangular Factorization (BTF) method [20], [27] can be used during symbolic analysis of A to transform it into a block-diagonal form A BTF . To do so, two matrices are calculated and applied, the row permutation matrix P R and the column permutation matrix P C , and result into This process is able to automatically find the subnetworks decoupled by distributed parameter transmission line/cable models. The blocks in A BTF are decoupled (independent) at each solution time-point and can be solved in parallel. The resulting block sizes vary according to subnetwork sizes.
The blocks in A BTF can be solved using any efficient sparse matrix solver. Let us assume that A i is the ith block in A BTF and there are N blocks in the network.
The sparse solver method selected here is the KLU method. It is currently one of the most efficient packages for circuitbased simulations [20]. This choice was made after testing other packages, including SuperLU [28] and [26].
The KLU solver, as other sparse matrix solvers, applies ordering to reduce the fillings in the LU-factorization. The AMD (approximate minimum degree ordering) [19] and COLAMD (column approximate minimum degree ordering) [19] ordering methods are available in KLU. AMD allowed to achieve the best performance for the test cases studied in this paper. The symbolic analysis of AMD is conducted only once at the startup of the simulation process.
The second stage in KLU is the numerical analysis stage which involves two main steps. First the Gilbert-Peierls algorithm is used to find the nonzero patterns of L and U matrices. Once the nonzero pattern is established, a left-looking numerical factorization with partial pivoting is used to calculate the L and U numerical values. This process can be expressed by where P R i and P C i are now permutation matrices for the ith block. The complete LU factorization process [19] of KLU is summarized (pseudocode) in Fig. 1 for an n × n matrix A i . It is noticed that the matrix L i is calculated from left to right and column by column.
In the above code, I is the identity matrix.
In KLU it is also possible to apply a scaling matrix for the above equation, but it was not retained here since its calculation process can create inefficiencies and there were no advantages for the applications studied in this paper.
In the standard KLU approach, when at least one element of A i is modified, it is needed to perform full re-factorization with updated pivoting. This process becomes inefficient for the refactorization needs expressed in the previous section. In KLU it is also possible to request refactorization reusing the previously calculated nonzero pattern, but in that case there is no pivot validity testing and the solution accuracy may deteriorate significantly.
To improve the performance of the existing KLU code for the EMT computation methods presented in this paper, the following enhancements are contributed in this paper: ''Pivot validity test'' and ''Partial re-factorization'' using a parallel implementation.
In the following text, the acronym SMPEMT (sparse matrix package for EMTs) refers to the enhanced KLU version presented in this paper.

A. PIVOT VALIDITY TEST
One of the challenging problems is the refactorization process for time-varying and iterative model coefficients in A i . Instead of performing full refactorization (standard KLU) as explained above, significant computational gains can be achieved by applying pivot validity testing. This algorithm is described by the following steps for A i : 1. First call to SMPEMT: a. Perform symbolic analysis with AMD b. Perform full factorization to find L i and U i c. Go to step 3 2. Not first call to SMPEMT: a) Perform re-factorization using existing nonzero pattern with new enabled invalid pivot testing i. Invalid pivot found: go to step 1.b b) Go to step 3 3. Backward and Forward substitution The step 2.a is applied when matrix A i has changed. The structure (non-zero elements) of the matrix is established at the beginning and cannot change, only matrix cell values can change.
For the step 2, the pseudocode used in Fig. 1 starts with a modified matrix, but initially keeps the previously calculated permutations due to partial pivoting (step 2.b in Fig. 1). The pivot validity test is now applied instead of step 2.b in Fig. 1. The invalid pivot (pivot validity) test is based on a definable tolerance ε p . When ε p a new > a old (6) the old pivot a old is discarded and a new pivot a new is selected.
The process restarts as shown in the step 2.a.i above. The pivot search is performed in each column of L i . By default ε p = 1%. It is noticed that for circuit equations, the diagonal VOLUME 7, 2020 of A i is typically dominant and consequently the number of pivot changes is small.

B. PARTIAL RE-FACTORIZATION
To reduce the computational cost of refactorization for A i even further, it is possible to apply partial refactorization. In a given matrix A i it is possible to determine the cells occupied by NMs and TVMs. Those dynamic cells may change between solution time-points and iterations at a given time-point. The changes require refactorization (recalculating L i and U i ).
The matrix A i is reordered using AMD and can be written as A i = P cc P cd P dc P dd (7) where the subscripts c and d mean constant and dynamic respectively. The c columns do not have any dynamic parts, but the d columns have dynamic cells.
In the left-looking algorithm of Fig. 1, the columns of L i are derived one-by-one by solving for each column of A i . If, for example, the lower matrix decomposition of A i is stopped at its first dynamic column then L p i = L cc 0 P dc I dd (8) where L p i is a partially computed lower-triangular matrix, L cc is a lower-triangular matrix, P dc and L cc contain the columns of the static part of L p i and I dd is the identity matrix. Once P dd is determined (status of time-varying devices or iterative Norton equivalent) at a given solution time-point, the calculation process is continued until the replacement of I dd to obtain L i from L p i . The upper-triangular matrix U i is calculated within the calculation process of L i . For (8), the partial upper matrix factorization is available for all constant columns In the above approach it is not necessary to restart the partial factorization process for the complete set of columns of P dd . Better efficiency can be gained, if partial re-factorization is applied by restarting from the first left modified column j m d in P dd . As before, since KLU is a leftlooking solver, all unchanged columns to the left of j m d can maintain the previous contributions to the LU factors.
It is also possible to apply a permutation technique that forces P dd to contain only NMs and TVMs (similar to [29]). But such an approach interferes with the AMD ordering and creates extra fill-ins which hinder the performance gains.

IV. IMPLEMENTATION OF SMPEMT
Both the sequential and parallel versions of the SMPEMT described in the preceding section have been implemented in this paper. The parallel version has been programmed using the OpenMP [30] multithreading approach. It can be advantageously used in any EMT-type package as a highlevel programming method.
The parallel solver is developed based on the distributed thread memory concept that was found to be the most efficient application of OpenMP. This approach minimizes thread waiting time and memory access delays.
Parallelization is applied to the blocks found from the BTF, that is for each matrix A i . Each matrix A i has its own memory allocation that is part of designated thread memory. Each thread launches a call to the sparse matrix solver for its own A i . The call allocates (done once, at startup only) memory on the thread stack and uses thread's own cache to solve. This type of allocation helps all threads to fetch data faster and write results quicker compared to getting data and writing data to the main thread memory. All threads are launched and joined using OpenMP statements.
Since there is a limited number of CPU cores and the computing gains are limited by the largest network blocks, it is necessary to apply a balancing technique for the given number of cores. In this paper, the following approaches have been implemented.
In the first approach a pre-programmed method allows to estimate the number of floating-point operations for the solution of each block (NFPO i ). The following formula is used for a matrix block of n × n: where j and m are the indices of columns, Llen and Ulen are the counts of non-zeros in L i and U i , respectively. This formula accounts for the LU-factorization based on the initial LU nonzero patterns. It also accounts for the backwardforward substitution step.
In the second approach, the number of non-zeros in each block (NNZ i ) is available from its nonzero pattern. However, for all the test cases presented in this paper, using NNZ i was less efficient than using NFPO i .
The blocks are assigned to cores using the number of available cores (N C ) and the factor k d = NFPO/N C , where NFPO is the total number for the entire system of equations. Since the number of blocks N could be higher than N C , it is necessary to verify the limiting k d for each assignment. If a given core is assigned a block with NFPO i less than k d then it can contain additional blocks until k d is reached or exceeded. This is basically a packing procedure for populating available cores.
If a block's NFPO i falls below a minimum size, then it must be packed into an assigned core since threading for such a block can become inefficient.
The same is applicable using NNZ i instead of NFPO i .

V. TEST CASES
The objective of this section is to test SMPEMT for practical power system grids. The tests are conducted directly without any topological/data changes (user intervention) on existing networks. This aspect is fundamental, since it is not desirable to ask software users to perform manipulations for helping implicated numerical issues. The presented realistic cases have different levels of complexities and allow to test different aspects. In all cases, the simulation is initialized to reach steady-state from a load-flow solution and a perturbation (fault) is applied in time-domain. The computer used for testing was a Lenovo Intel(R) Xeon(R) CPU E5-2650 v4 @ 2.20GHz, with two clusters of 12 cores, with cache of 30MB and RAM of 32GB.
In the following tables, the tested solvers are: EMTP [24] for the existing sparse matrix package [26] programmed on single core; KLU represents the standard KLU singlecore code without the improvements presented in this paper, SMPEMT1 is the upgraded KLU solver with the above ''Pivot validity test'' and SMPEMT2 is the upgraded KLU solver with both ''Pivot validity test'' and ''Partial re-factorization''. Both SMPEMT1 and SMPEMT2 are implemented in parallel and contrary to the standard KLU, these two methods are set to refactor only the changed matrices A i .
In all cases in this paper, SMPEMT2 was faster than SMPEMT1, but the gains depend on the presence and location of dynamic columns, which is related to the simulated case equations.
The simulation results are identical with all solvers. The applied models are listed for each case. The synchronous generators are modeled using the dq0 approach (see [32]). Details about solution methods including the iterative process for nonlinear models, can be found in [25].
In the text below, computational gain or performance gain is actually acceleration ratio between different solution methods with different core numbers.

A. FIRST TEST CASE: HYDRO-QUEBEC GRID
This test case is an upgraded version of the one presented in [3]. It is based on the actual Hydro-Quebec grid (HQ-grid). The initial network matrix A is shown in Fig. 2. Its BTF version is presented in Fig. 3. The size of A is 41815 × 41815 with 168025 non-zeros. The network has 29803 nodes. The BTF matrix version has 217 blocks, with the largest block size being 2898 × 2898.
The HQ-grid performance timings are presented in Table 1    It is apparent that the KLU method alone does not have performance gains. It is even less efficient than the EMTP method. This is due to the fact that re-factorization is applied to all blocks without the improvements of SMPEMT1 and SMPEMT2.
For this case, SMPEMT2 and SMPEMT1 gave similar timings.
The computational gain against existing EMTP solver is 1048/31 = 33.8 with 12 cores. The gain over the standard single-core KLU solver is 1216/31 = 39 with 12 cores.
Performance plots are presented in Fig. 4. The maximum gain over the single core SMPEMT2 version is 9.2 and there are no significant gains after 12 cores. This is mainly due to the limitation imposed by the largest block, increased memory exchange and thread management costs.
Comparison of simulation results between EMTP and SMPEMT2 for a fault event are shown in Fig. 5 and Fig. 6 (zoomed version). The waveforms are identical. The studied scenario is a fault on phase-a to ground occurring at 400 ms.   The breaker SWm on the right of a transmission line opens its phase-a at current-zero crossing near 472 ms.

B. SECOND TEST CASE: T0-GRID
This test case is designed to stress computational performance using a system with integration of renewables. The top level view of the simulation network is presented in Fig. 10. The complete data of this network is available in [31]. This is a realistic 400 kV, 50 Hz network. The yellow boxes include the sub-transmission networks at 154 kV with pi-section models for transmission lines. The green boxes include 154 kV sub-transmission and wind generation.
There is a total of 72 synchronous generators. All generators are modeled with exciters and governors. The 10.5 kV generators connected to the 400 kV grid include total (dq) saturation data.
There is a total of 10 wind parks with aggregated wind generators. The wind parks are operating on reactive power control mode. The wind generators are of DFIG-type and include all control and protection systems. The DFIG converters are modeled using the detailed modeling (see DFIG DM in [4]) approach which includes nonlinear IGBT models with diodes in parallel. Further details on the IGBT model also used for HVDC converters, can be found in [4], [5].
A second version of this test case with average-value models (AVMs) [4] for converters is also tested below.
In addition to the above, the network contains: The initial network matrix A is shown in Fig. 7. Its BTF version is presented in Fig. 8. The size of A is 4703 × 4703 with 25117 non-zeros. The network has 3407 nodes. The BTF matrix version has 28 blocks, with the largest block size being 573 × 573 and the smallest block size is 3 × 3. It is apparent from Fig. 8 that the presence of the sub-transmission grids creates a more balanced distribution of blocks. The subtransmission networks are using pi-section type line models (short) and disable further decoupling in matrix blocks.
The complete network is initialized from its load-flow solution. The studied event is a (phase-a-to-ground) fault on the transmission line ADAPA_to_GOKCE connected between the lines ADAPA and GOKCE (see fault symbol in Fig. 10). The fault occurs at 1 s, the phase-a breaker on the left of the line receives the opening signal at 1.08 s and the one on the right at 1.1 s. The phase-a breaker on the left recloses at 1.48 s and the one on the right at 1.5 s. The reclosing is unsuccessful and all breakers (all left and right phases) receive the opening signal at 1.56 s to isolate the line.
The T0-grid performance timings are presented in Table 2 (DM) and Table 3 (AVM) for a numerical integration timestep of t = 10 µs and for the simulation time of 2 s. The average number of iterations per time-point is 5.91 in the DM case and 1.7 in the AVM case. The higher number in the DM case is mainly due to the nonlinear IGBT models in the converters. This is a numerically challenging benchmark.
It is observed again that the standard KLU method is slower than EMTP due to numerous refactorizations. In the DM case SMPEMT2 is able to gain over SMPEMT1 due to numerous nonlinearities. There is no difference between both methods for the AVM version.
A gain of 2580/204 = 12.6 over EMTP is achieved with 12 cores for the DM case. In the case of AVM, the gain is 253/27 = 9.4 with 8 cores. The same gains over the singlecore KLU method are 4379/204 = 21.5 and 336/27 = 12.4, respectively.
Performance plots with SMPEMT2 for the DM version are presented in Fig. 9. The maximum gain over the single core SMPEMT2 solution is 6.9 with 12 cores.
It is apparent from Table 3 that the largest block in this network (see subnetwork AI in Fig. 10) imposes a limit after 8 cores with the AVM case. The AVM converters do not require refactorization. For the case of DM converters, it is observed in Table 2 that the SMPEMT2 method is able to gain after 8 cores. This is due to partial re-factorization option of  SMPEMT2. On average, the largest block is able to reduce its computing timings with repetitive re-factorizations and perform faster.
The computational performance with AVM even on a single core, is already significantly better than with DM and that is also why it becomes difficult to improve it with more cores.
Although the approach proposed in this paper is based on the complete avoidance of user intervention in simulation cases, it has been observed that inserting a one-time-step delay transmission line for decoupling the wind generator model from its sub-transmission network, allows to achieve gains up to 13 cores for both DM and AVM versions. In the case of DM the maximum gain over 1-core becomes 7.7. For AVM the same gain increases to 4.7.
The simulation results from EMTP and SMPEMT2 are compared in Fig. 11. The waveforms are perfectly overlapped and confirm the validity of SMPEMT2 computations.

VI. DISCUSSION
The above demonstrations have been made for practical networks and it is apparent that performances of parallel computations vary according to simulated network topology and contents. In this paper the parallelization process is based on transmission line delay decoupling of networks into subnetworks. A limitative factor is the size and/or contents of resulting subnetworks. Further parallelization may be achieved by tearing subnetworks using diakoptics. This means cutting through generic branches when transmission lines are not available. It comes back to the usage of a bordered-blockdiagonal matrix in the form of In (11), the matrix A i (ith subnetwork found through BTF) is decomposed, for example, into two sub-subnetworks (matrices A i 1 and A i 2 ) interconnected by the border matrices S. In the MANA formulation recalled in section II, it is possible to use interconnecting wires modeled by ideal switches, in which case S d = 0 and b i s = 0. The switch currents are given by i i s . Any number of sub-subnetworks can be created. Equation (11) can be solved in parallel using, for example, the approach presented in [18], [34], [35], or the compensation method [29], [33]. The difficulty in such methods is that in practical networks, such as those presented in this paper, the matrices A i j encounter switching events and/or contain nonlinearities which require repetitive refactorizations and supplementary duplications of time-consuming matrix operations (see, for example, determination of LU factors in [18]) required to account for the border matrices S j . Moreover, the automatic determination of the bordered-block-diagonal system requires further research for arbitrary topology networks used in this paper. The system used in [18] is based on duplication of a small system and does not account for realistic network topologies.
It was not possible to use the approach of (11) to achieve further gains for the test cases of this paper.

VII. CONCLUSION
This paper presents a new sparse matrix based method for acceleration of computations in the simulation of electromagnetic transients. The proposed method accounts for varying topologies and the accurate iterative solution of nonlinear models. The proposed method is generic and demonstrated on practical systems.
The network equations are solved using a sparse matrix solver (KLU). Parallelization is programmed after applying block-triangular factorization. This is a fully automatic process and does not require any user intervention for arbitrary small and large scale power systems.
In addition to parallelization, the KLU solver is improved using pivot validity testing and partial refactorization. This is an important contribution for achieving shown significant computational gains.
The proposed method is applicable to any software tool for the computation of electromagnetic transients. Moreover, the proposed enhancements to the KLU solver are applicable to other power system computation tools.
The computational gains are demonstrated for practical and large networks. The demonstration benchmarks and results constitute another contribution of this paper for research on practical problem solutions.