Source author record

Stefano Markidis

Stefano Markidis appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

26works
12topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

26 published item(s)

preprint2022arXiv

FFTc: An MLIR Dialect for Developing HPC Fast Fourier Transform Libraries

Discrete Fourier Transform (DFT) libraries are one of the most critical software components for scientific computing. Inspired by FFTW, a widely used library for DFT HPC calculations, we apply compiler technologies for the development of HPC Fourier transform libraries. In this work, we introduce FFTc, a domain-specific language, based on Multi-Level Intermediate Representation (MLIR), for expressing Fourier Transform algorithms. We present the initial design, implementation, and preliminary results of FFTc.

preprint2022arXiv

Large-Scale Direct Numerical Simulations of Turbulence Using GPUs and Modern Fortran

We present our approach to making direct numerical simulations of turbulence with applications in sustainable shipping. We use modern Fortran and the spectral element method to leverage and scale on supercomputers powered by the Nvidia A100 and the recent AMD Instinct MI250X GPUs, while still providing support for user software developed in Fortran. We demonstrate the efficiency of our approach by performing the world's first direct numerical simulation of the flow around a Flettner rotor at Re=30'000 and its interaction with a turbulent boundary layer. We present one of the first performance comparisons between the AMD Instinct MI250X and Nvidia A100 GPUs for scalable computational fluid dynamics. Our results show that one MI250X offers performance on par with two A100 GPUs and has a similar power efficiency.

preprint2022arXiv

On Physics-Informed Neural Networks for Quantum Computers

Physics-Informed Neural Networks (PINN) emerged as a powerful tool for solving scientific computing problems, ranging from the solution of Partial Differential Equations to data assimilation tasks. One of the advantages of using PINN is to leverage the usage of Machine Learning computational frameworks relying on the combined usage of CPUs and co-processors, such as accelerators, to achieve maximum performance. This work investigates the design, implementation, and performance of PINNs, using the Quantum Processing Unit (QPU) co-processor. We design a simple Quantum PINN to solve the one-dimensional Poisson problem using a Continuous Variable (CV) quantum computing framework. We discuss the impact of different optimizers, PINN residual formulation, and quantum neural network depth on the quantum PINN accuracy. We show that the optimizer exploration of the training landscape in the case of quantum PINN is not as effective as in classical PINN, and basic Stochastic Gradient Descent (SGD) optimizers outperform adaptive and high-order optimizers. Finally, we highlight the difference in methods and algorithms between quantum and classical PINNs and outline future research challenges for quantum PINN development.

preprint2022arXiv

Workflows to driving high-performance interactive supercomputing for urgent decision making

Interactive urgent computing is a small but growing user of supercomputing resources. However there are numerous technical challenges that must be overcome to make supercomputers fully suited to the wide range of urgent workloads which could benefit from the computational power delivered by such instruments. An important question is how to connect the different components of an urgent workload; namely the users, the simulation codes, and external data sources, together in a structured and accessible manner. In this paper we explore the role of workflows from both the perspective of marshalling and control of urgent workloads, and at the individual HPC machine level. Ultimately requiring two workflow systems, by using a space weather prediction urgent use-cases, we explore the benefit that these two workflow systems provide especially when one exploits the flexibility enabled by them interoperating.

preprint2020arXiv

Magnetohydrodynamic with embedded particle-in-cell simulation of the Geospace Environment Modeling dayside kinetic processes challenge event

We use the MHD with embedded particle-in-cell model (MHD-EPIC) to study the Geospace Environment Modeling (GEM) dayside kinetic processes challenge event at 01:50-03:00 UT on 2015-11-18, when the magnetosphere was driven by a steady southward IMF. In the MHD-EPIC simulation, the dayside magnetopause is covered by a PIC code so that the dayside reconnection is properly handled. We compare the magnetic fields and the plasma profiles of the magnetopause crossing with the MMS3 spacecraft observations. Most variables match the observations well in the magnetosphere, in the magnetosheath, and also during the current sheet crossing. The MHD-EPIC simulation produces flux ropes, and we demonstrate that some magnetic field and plasma features observed by the MMS3 spacecraft can be reproduced by a flux rope crossing event. We use an algorithm to automatically identify the reconnection sites from the simulation results. It turns out that there are usually multiple X-lines at the magnetopause. By tracing the locations of the X-lines, we find the typical moving speed of the X-line endpoints is about 70~km/s, which is higher than but still comparable with the ground-based observations.

preprint2020arXiv

Optimization of Tensor-product Operations in Nekbone on GPUs

In the CFD solver Nek5000, the computation is dominated by the evaluation of small tensor operations. Nekbone is a proxy app for Nek5000 and has previously been ported to GPUs with a mixed OpenACC and CUDA approach. In this work, we continue this effort and optimize the main tensor-product operation in Nekbone further. Our optimization is done in CUDA and uses a different, 2D, thread structure to make the computations layer by layer. This enables us to use loop unrolling as well as utilize registers and shared memory efficiently. Our implementation is then compared on both the Pascal and Volta GPU architectures to previous GPU versions of Nekbone as well as a measured roofline. The results show that our implementation outperforms previous GPU Nekbone implementations by 6-10%. Compared to the measured roofline, we obtain 77 - 92% of the peak performance for both Nvidia P100 and V100 GPUs for inputs with 1024 - 4096 elements and polynomial degree 9.

preprint2019arXiv

Performance Evaluation of Advanced Features in CUDA Unified Memory

CUDA Unified Memory improves the GPU programmability and also enables GPU memory oversubscription. Recently, two advanced memory features, memory advises and asynchronous prefetch, have been introduced. In this work, we evaluate the new features on two platforms that feature different CPUs, GPUs, and interconnects. We derive a benchmark suite for the experiments and stress the memory system to evaluate both in-memory and oversubscription performance. The results show that memory advises on the Intel-Volta/Pascal-PCIe platform bring negligible improvement for in-memory executions. However, when GPU memory is oversubscribed by about 50%, using memory advises results in up to 25% performance improvement compared to the basic CUDA Unified Memory. In contrast, the Power9-Volta-NVLink platform can substantially benefit from memory advises, achieving up to 34% performance gain for in-memory executions. However, when GPU memory is oversubscribed on this platform, using memory advises increases GPU page faults and results in considerable performance loss. The CUDA prefetch also shows different performance impact on the two platforms. It improves performance by up to 50% on the Intel-Volta/Pascal-PCI-E platform but brings little benefit to the Power9-Volta-NVLink platform.

preprint2019arXiv

Posit NPB: Assessing the Precision Improvement in HPC Scientific Applications

Floating-point operations can significantly impact the accuracy and performance of scientific applications on large-scale parallel systems. Recently, an emerging floating-point format called Posit has attracted attention as an alternative to the standard IEEE floating-point formats because it could enable higher precision than IEEE formats using the same number of bits. In this work, we first explored the feasibility of Posit encoding in representative HPC applications by providing a 32-bit Posit NAS Parallel Benchmark (NPB) suite. Then, we evaluate the accuracy improvement in different HPC kernels compared to the IEEE 754 format. Our results indicate that using Posit encoding achieves optimized precision, ranging from 0.6 to 1.4 decimal digit, for all tested kernels and proxy-applications. Also, we quantified the overhead of the current software implementation of Posit encoding as 4x-19x that of IEEE 754 hardware implementation. Our study highlights the potential of hardware implementations of Posit to benefit a broad range of HPC applications.

preprint2019arXiv

Studying dawn-dusk asymmetries of Mercury's magnetotail using MHD-EPIC simulations

MESSENGER has observed a lot of dawn-dusk asymmetries in Mercury's magnetotail, such as the asymmetries of the cross-tail current sheet thickness and the occurrence of flux ropes, dipolarization events and energetic electron injections. In order to obtain a global pictures of Mercury's magnetotail dynamics and the relationship between these asymmetries, we perform global simulations with the magnetohydrodynamics with embedded particle-in-cell (MHD-EPIC) model, where Mercury's magnetotail region is covered by a PIC code. Our simulations show that the dawnside current sheet is thicker, the plasma density is larger, and the electron pressure is higher than the duskside. Under a strong IMF driver, the simulated reconnection sites prefer the dawnside. We also found the dipolarization events and the planetward electron jets are moving dawnward while they are moving towards the planet, so that almost all dipolarization events and high-speed plasma flows concentrate in the dawn sector. The simulation results are consistent with MESSENGER observations.

preprint2019arXiv

TensorFlow Doing HPC

TensorFlow is a popular emerging open-source programming framework supporting the execution of distributed applications on heterogeneous hardware. While TensorFlow has been initially designed for developing Machine Learning (ML) applications, in fact TensorFlow aims at supporting the development of a much broader range of application kinds that are outside the ML domain and can possibly include HPC applications. However, very few experiments have been conducted to evaluate TensorFlow performance when running HPC workloads on supercomputers. This work addresses this lack by designing four traditional HPC benchmark applications: STREAM, matrix-matrix multiply, Conjugate Gradient (CG) solver and Fast Fourier Transform (FFT). We analyze their performance on two supercomputers with accelerators and evaluate the potential of TensorFlow for developing HPC applications. Our tests show that TensorFlow can fully take advantage of high performance networks and accelerators on supercomputers. Running our TensorFlow STREAM benchmark, we obtain over 50% of theoretical communication bandwidth on our testing platform. We find an approximately 2x, 1.7x and 1.8x performance improvement when increasing the number of GPUs from two to four in the matrix-matrix multiply, CG and FFT applications respectively. All our performance results demonstrate that TensorFlow has high potential of emerging also as HPC programming framework for heterogeneous supercomputers.

preprint2016arXiv

Energy exchanges in reconnection outflows

Reconnection outflows are highly energetic directed flows that interact with the ambient plasma or with flows from other reconnection regions. Under these conditions the flow becomes highly unstable and chaotic, as any flow jets interacting with a medium. We report here massively parallel simulations of the two cases of interaction between outflow jets and between a single outflow with an ambient plasma. We find in both case the development of a chaotic magnetic field, subject to secondary reconnection events that further complicate the topology of the field lines. The focus of the present analysis is on the energy balance. We compute each energy channel (electromagnetic, bulk, thermal, for each species) and find where the most energy is exchanged and in what form. The main finding is that the largest energy exchange is not at the reconnection site proper but in the regions where the outflowing jets are destabilized.

preprint2016arXiv

On the ions acceleration via collisionless magnetic reconnection in laboratory plasmas

This work presents an analysis of the ion outflow from magnetic reconnection throughout fully kinetic simulations with typical laboratory plasmas values. A symmetric initial configuration for the density and magnetic field is considered across the current sheet. After analyzing the behavior of a set of nine simulations with a reduced mass ratio and with a permuted value of three initial electron temperature and magnetic field intensity, the best ion acceleration scenario is further studied with a realistic mass ratio in terms of the ion dynamics and energy budget. Interestingly, a series of shock waves structures are observed in the outflow, resembling the shock discontinuities found in recent magnetohydrodynamic (MHD) simulations. An analysis of the ion outflow at several distances from the reconnection point is presented, in light of possible laboratory applications. The analysis suggests that magnetic reconnection could be used as a tool for plasma acceleration, with applications ranging from electric propulsion to production of ion thermal beams.

preprint2015arXiv

Energetics of kinetic reconnection in a three-dimensional null points cluster

We performed three-dimensional Particle-in-Cell simulations of magnetic reconnection with multiple magnetic null points. Magnetic field energy conversion into kinetic energy was about five times higher than in traditional Harris sheet configuration. More than 85% of initial magnetic field energy was transferred to particle energy during 25 reversed ion cyclofrequencies. Magnetic reconnection in the cluster of null points evolved in three phases. During the first phase, ion beams were excited, that then gave part of their energy back to magnetic field in the second phase. In the third phase, magnetic reconnection occurs in many small patches around the current channels formed along the stripes of low magnetic field. Magnetic reconnection in null points presents essentially three-dimensional features, with no two dimensional symmetries or current sheets.

preprint2015arXiv

Energy dissipation in magnetic null points at kinetic scales

We use kinetic particle-in-cell and magnetohydrodynamic simulations supported by an observational dataset to investigate magnetic reconnection in clusters of null points in space plasma. The magnetic configuration under investigation is driven by fast adiabatic flux rope compression that dissipates almost half of the initial magnetic field energy. In this phase powerful currents are excited producing secondary instabilities, and the system is brought into a state of `intermittent turbulence' within a few ion gyro-periods. Reconnection events are distributed all over the simulation domain and energy dissipation is rather volume-filling. Numerous spiral null points interconnected via their spines form null lines embedded into magnetic flux ropes; null point pairs demonstrate the signatures of torsional spine reconnection. However, energy dissipation mainly happens in the shear layers formed by adjacent flux ropes with oppositely directed currents. In these regions radial null pairs are spontaneously emerging and vanishing, associated with electron streams and small-scale current sheets. The number of spiral nulls in the simulation outweighs the number of radial nulls by a factor of 5\---10, in accordance with Cluster observations in the Earth's magnetosheath. Twisted magnetic fields with embedded spiral null points might indicate the regions of major energy dissipation for future space missions such as Magnetospheric Multiscale Mission (MMS).

preprint2015arXiv

Magnetic null points in kinetic simulations of space plasmas

We present a systematic attempt to study magnetic null points and the associated magnetic energy conversion in kinetic Particle-in-Cell simulations of various plasma configurations. We address three-dimensional simulations performed with the semi-implicit kinetic electromagnetic code iPic3D in different setups: variations of a Harris current sheet, dipolar and quadrupolar magnetospheres interacting with the solar wind; and a relaxing turbulent configuration with multiple null points. Spiral nulls are more likely created in space plasmas: in all our simulations except lunar magnetic anomaly and quadrupolar mini-magnetosphere the number of spiral nulls prevails over the number of radial nulls by a factor of 3-9. We show that often magnetic nulls do not indicate the regions of intensive energy dissipation. Energy dissipation events caused by topological bifurcations at radial nulls are rather rare and short-lived. The so-called X-lines formed by the radial nulls in the Harris current sheet and lunar magnetic anomaly simulations are rather stable and do not exhibit any energy dissipation. Energy dissipation is more powerful in the vicinity of spiral nulls enclosed by magnetic flux ropes with strong currents at their axes (their cross-sections resemble 2D magnetic islands). These null lines reminiscent of Z-pinches efficiently dissipate magnetic energy due to secondary instabilities such as the two-stream or kinking instability, accompanied by changes in magnetic topology. Current enhancements accompanied by spiral nulls may signal magnetic energy conversion sites in the observational data.

preprint2015arXiv

On the electron dynamics during island coalescence in asymmetric magnetic reconnection

We present an analysis of the electron dynamics during rapid island merging in asymmetric magnetic reconnection. We consider a doubly periodic system with two asymmetric transitions. The upper layer is an asymmetric Harris sheet initially perturbed to promote a single reconnection site. The lower layer is a tangential discontinuity that promotes the formation of many X-points, separated by rapidly merging islands. Across both layers the magnetic field and the density have a strong jump, but the pressure is held constant. Our analysis focuses on the consequences of electron energization during island coalescence. We focus first on the parallel and perpendicular components of the electron temperature to establish the presence of possible anisotropies and non-gyrotropies. Thanks to the direct comparison between the two different layers simulated, we can distinguish three main types of behavior characteristic of three different regions of interest. The first type represents the regions where traditional asymmetric reconnections take place without involving island merging. The second type of regions instead show reconnection events between two merging islands. Finally, the third regions identifies the regions between two diverging island and where typical signature of reconnection is not observed. Electrons in these latter regions additionally show a flat-top distribution resulting from the saturation of a two-stream instability generated by the two interacting electron beams from the two nearest reconnection points. Finally, the analysis of agyrotropy shows the presence of a distinct double structure laying all over the lower side facing the higher magnetic field region. The distinguishing features found for the three types of regions investigated provide clear indicators to the recently launched MMS NASA mission for investigating magnetopause reconnection involving multiple islands.

preprint2015arXiv

Role of Z-pinches in magnetic reconnection in space plasmas

A widely accepted scenario of magnetic reconnection in collisionless space plasmas is the breakage of magnetic field lines in X-points. In laboratory, reconnection is commonly studied in pinches, current channels embedded into twisted magnetic fields. No model of magnetic reconnection in space plasmas considers both null-points and pinches as peers. We have performed a particle-in-cell simulation of magnetic reconnection in a three-dimensional configuration where null-points are present initially, and Z-pinches are formed during the simulation along the lines of spiral null-points. The non-spiral null-points are more stable than spiral ones, and no substantial energy dissipation is associated with them. On the contrary, turbulent magnetic reconnection in the pinches causes the magnetic energy to decay at a rate of ~1.5% per ion gyro period. Dissipation in similar structures is a likely scenario in space plasmas with large fraction of spiral null-points.

preprint2014arXiv

Electromagnetic Energy Conversion in Downstream Fronts from 3D Kinetic Reconnection

The electromagnetic energy equation is analyzed term by term in a 3D simulation of kinetic reconnection previously reported by \citet{vapirev2013formation}. The evolution presents the usual 2D-like topological structures caused by an initial perturbation independent of the third dimension. However, downstream of the reconnection site, where the jetting plasma encounters the yet unperturbed pre-existing plasma, a downstream front (DF) is formed and made unstable by the strong density gradient and the unfavorable local acceleration field. The energy exchange between plasma and fields is most intense at the instability, reaching several $pW/m^3$, alternating between load (energy going from fields to particles) and generator (energy going from particles to fields) regions. Energy exchange is instead purely that of a load at the reconnection site itself in a region focused around the x-line and elongated along the separatrix surfaces. Poynting fluxes are generated at all energy exchange regions and travel away from the reconnection site transporting an energy signal of the order of about $\mathbf S \approx 10^{-3} W/m^2$.

preprint2014arXiv

Kinetic plasma turbulence during the nonlinear stage of the Kelvin-Helmholtz instability

Using a full kinetic, implicit particle-in-cell code, iPiC3D, we studied the properties of plasma kinetic turbulence, such as would be found at the interface between the solar wind and the Earth magnetosphere at low latitude during northwards periods. In this case, in the presence of a magnetic field B oriented mostly perpendicular to the velocity shear, turbulence is fed by the disruption of a Kelvin-Helmholtz vortex chain via secondary instabilities, vortex pairing and non-linear interactions. We found that the magnetic energy spectral cascade between ion and electron inertial scales, $d_i$ and $d_e$, is in agreement with satellite observations and other previous numerical simulations; however, in our case the spectrum ends with a peak beyond $d_e$ due to the occurrence of the lower hybrid drift instability. The electric energy spectrum is influenced by effects of secondary instabilities: anomalous resistivity, fed by the development of the lower hybrid drift instability, steepens the spectral decay and, depending on the alignment or anti-alignment of B and the shear vorticity, peaks due to ion-Bernstein waves may dominate the spectrum around $d_i$. These waves are generated by counter-streaming flow structures, through flux freezing also responsible for reconnection of the in-plane component of the magnetic field, which then generates electron pressure anisotropy and flattening of the field-aligned component of the electron distribution function.

preprint2013arXiv

Kinetic Simulations of Plasmoid Chain Dynamics

The dynamics of a plasmoid chain is studied with three dimensional Particle-in-Cell simulations. The evolution of the system with and without a uniform guide field, whose strength is 1/3 the asymptotic magnetic field, is investigated. The plasmoid chain forms by spontaneous magnetic reconnection: the tearing instability rapidly disrupts the initial current sheet generating several small-scale plasmoids, that rapidly grow in size coalescing and kinking. The plasmoid kink is mainly driven by the coalescence process. It is found that the presence of guide field strongly influences the evolution of the plasmoid chain. Without a guide field, a main reconnection site dominates and smaller reconnection regions are included in larger ones, leading to an hierarchical structure of the plasmoid-dominated current sheet. On the contrary in presence of a guide field, plasmoids have approximately the same size and the hierarchical structure does not emerge, a strong core magnetic field develops in the center of the plasmoid in the direction of the existing guide field, and bump-on-tail instability, leading to the formation of electron holes, is detected in proximity of the plasmoids.

preprint2012arXiv

Collisionless magnetic reconnection in a plasmoid chain

The kinetic features of plasmoid chain formation and evolution are investigated by two dimensional Particle-in-Cell simulations. Magnetic reconnection is initiated in multiple X points by the tearing instability. Plasmoids form and grow in size by continuously coalescing. Each chain plasmoid exhibits a strong out-of plane core magnetic field and an out-of-plane electron current that drives the coalescing process. The disappearance of the X points in the coalescence process are due to anti-reconnection, a magnetic reconnection where the plasma inflow and outflow are reversed with respect to the original reconnection flow pattern. Anti-reconnection is characterized by the Hall magnetic field quadrupole signature. Two new kinetic features, not reported by previous studies of plasmoid chain evolution, are here revealed. First, intense electric fields develop in-plane normally to the separatrices and drive the ion dynamics in the plasmoids. Second, several bipolar electric field structures are localized in proximity of the plasmoid chain. The analysis of the electron distribution function and phase space reveals the presence of counter-streaming electron beams, unstable to the two stream instability, and phase space electron holes along the reconnection separatrices.

preprint2012arXiv

Rethinking Electrostatic Solvers in Particle Simulations for the Exascale Era

In preparation to the exascale era, an alternative approach to calculate the electrostatic forces in Particle Mesh (PM) methods is proposed. While the traditional techniques are based on the calculation of the electrostatic potential by solving the Poisson equation, in the new approach the electric field is calculated by solving the Ampere's law. When the Ampere's law is discretized explicitly in time, the electric field values on the mesh are simply updated from the previous values. In this way, the electrostatic solver becomes an embarrassingly parallel problem, making the algorithm extremely scalable and suitable for exascale computing platforms. An implementation of a one dimensional PM code is presented to show that the proposed method produces correct results, and it is a very promising algorithm for exascale PM simulations.

preprint2012arXiv

Three dimensional density cavities in guide field collisionless magnetic reconnection

Particle-in-Cell simulations of collisionless magnetic reconnection with a guide field reveal for the first time the three dimensional features of the low density regions along the magnetic reconnection separatrices, the so-called "cavities". It is found that structures with further lower density develop within the cavities. Because their appearance is similar to the rib shape, these formations are here called "low density ribs". Their location remains approximately fixed in time and their density progressively decreases, as electron currents along the cavities evacuate them. They develop along the magnetic field lines and are supported by a strong perpendicular electric field that oscillates in space. In addition, bipolar parallel electric field structures form as isolated spheres between the cavities and the outflow plasma, along the direction of the low density ribs and of magnetic field lines.

preprint2011arXiv

Development and performance analysis of a UPC Particle-in-Cell code

The development and the implementation of a Particle-in-Cell code written in the Unified Parallel C (UPC) language for plasma simulations with application to astrophysics and fusion nuclear energy machines are presented. A simple one dimensional electrostatic Particle-in-Cell code has been developed first to investigate the implementation details in the UPC language, and second to study the UPC performance on parallel computers. The initial simulations of plasmas with the UPC Particle-in-Cell code and a study of parallel speed-up of the UPC code up to 128 cores are shown.

preprint2011arXiv

Kinetic simulations of magnetic reconnection in presence of a background O+ population

Particle-in-Cell simulations of magnetic reconnection with an H+ current sheet and a mixed background plasma of H+ and O+ ions are completed using physical mass ratios. Four main results are shown. First, the O+ presence slightly decreases the reconnection rate and the magnetic reconnection evolution depends mainly on the lighter H+ ion species in the presented simulations. Second, the Hall magnetic field is characterized by a two-scale structure in presence of O+ ions: it reaches sharp peak values in a small area in proximity of the neutral line, and then decreases slowly over a large region. Third, the two background species initially separate in the outflow region because H+ and O+ ions are accelerated by different mechanisms occurring on different time scales and with different strengths. Fourth, the effect of a guide field on the O+ dynamics is studied: the O+ presence does not change the reconnected flux and all the characteristic features of guide field magnetic reconnection are still present. Moreover, the guide field introduces an O+ circulation pattern between separatrices that enhances high O+ density areas and depletes low O+ density regions in proximity of the reconnection fronts. The importance and the validity of these results are finally discussed.

preprint2011arXiv

The Energy Conserving Particle-in-Cell Method

A new Particle-in-Cell (PIC) method, that conserves energy exactly, is presented. The particle equations of motion and the Maxwell's equations are differenced implicitly in time by the midpoint rule and solved concurrently by a Jacobian-free Newton Krylov (JFNK) solver. Several tests show that the finite grid instability is eliminated in energy conserving PIC simulations, and the method correctly describes the two-stream and Weibel instabilities, conserving exactly the total energy. The computational time of the energy conserving PIC method increases linearly with the number of particles, and it is rather insensitive to the number of grid points and time step. The kinetic enslavement technique can be effectively used to reduce the problem matrix size and the number of JFNK solver iterations.