# Multilingual Offensive Language Identification
with Cross-lingual Embeddings

Tharindu Ranasinghe
University of Wolverhampton
Wolverhampton, UK
<EMAIL_ADDRESS>
Marcos Zampieri
Rochester Institute of Technology
Rochester, NY, USA
<EMAIL_ADDRESS>

###### Abstract

Offensive content is pervasive in social media and a reason for concern to
companies and government organizations. Several studies have been recently
published investigating methods to detect the various forms of such content
(e.g. hate speech, cyberbulling, and cyberaggression). The clear majority of
these studies deal with English partially because most annotated datasets
available contain English data. In this paper, we take advantage of English
data available by applying cross-lingual contextual word embeddings and
transfer learning to make predictions in languages with less resources. We
project predictions on comparable data in Bengali, Hindi, and Spanish and we
report results of 0.8415 F1 macro for Bengali, 0.8568 F1 macro for Hindi, and
0.7513 F1 macro for Spanish. Finally, we show that our approach compares
favorably to the best systems submitted to recent shared tasks on these three
languages, confirming the robustness of cross-lingual contextual embeddings
and transfer learning for this task.

##  1 Introduction

Offensive posts on social media result in a number of undesired consequences
to users. They have been investigated as triggers of suicide attempts and
ideation, and mental health problems Bonanno and Hymel (2013); Bannink et al.
(2014). One of the most common ways to cope with offensive content online is
training systems to be capable of recognizing offensive messages or posts.
Once recognized, such offensive content can be set aside for human moderation
or deleted from the respective platform (e.g. Facebook, Twitter), preventing
harm to users and controlling the spread of abusive behavior in social media.

There have been several recent studies published on automatically identifying
various kinds of offensive content such as abuse Mubarak et al. (2017),
aggression Kumar et al. (2018, 2020), cyber-bullying Rosa et al. (2019), and
hate speech Malmasi and Zampieri (2018). While there are a few studies
published on languages such as Arabic Mubarak et al. (2020) and Greek Pitenis
et al. (2020), most studies and datasets created so far include English data.
Data augmentation Ghadery and Moens (2020) and multilingual word embeddings
Pamungkas and Patti (2019) have been applied to take advantage of existing
English resources to improve the performance in systems dealing with languages
other than English. To the best of our knowledge, however, state-of-the-art
cross-lingual contextual embeddings such as XLM-R Conneau et al. (2019) have
not yet been applied to offensive language identification. To address this
gap, we evaluate the performance of cross-lingual contextual embeddings and
transfer learning (TL) methods in projecting predictions from English to other
languages. We show that our methods compare favorably to state-of-the-art
approaches submitted to recent shared tasks on all datasets. The main
contributions of this paper are the following:

  1. 1.

We apply cross-lingual contextual word embeddings to offensive language
identification. We take advantage of existing English data to project
predictions in three other languages: Bengali, Hindi, and Spanish.

  2. 2.

We tackle both off-domain and off-task data for Bengali. We show that not only
can these methods project predictions for different languages but also for
different domains (e.g. Twitter vs. Facebook) and tasks (e.g. binary vs.
three-way classification).

  3. 3.

We provide important resources to the community: the code, and the English
model will be freely available to everyone interested in working on low-
resource languages using the same methodology.

##  2 Related Work

There is a growing interest in the development of computational models to
identify offensive content online. Early approaches relied heavily on feature
engineering combined with traditional machine learning classifiers such as
naive bayes and support vector machines Xu et al. (2012); Dadvar et al.
(2013). More recently, neural networks such as LSTMs, bidirectional LSTMs, and
GRUs combined with word embeddings have proved to outperform traditional
machine learning methods in this task Aroyehun and Gelbukh (2018); Majumder et
al. (2018). In the last couple of years, transformer models like ELMO Peters
et al. (2018) and BERT Devlin et al. (2019) have been applied to offensive
language identification achieving competitive scores and topping the
leaderboards in recent shared tasks Liu et al. (2019); Ranasinghe et al.
(2019). Most of these approaches use existing pre-trained transformer models
which can also be used as text classification models.

The clear majority of studies on this topic deal with English Malmasi and
Zampieri (2017); Yao et al. (2019); Ridenhour et al. (2020) partially
motivated by the availability English resources (e.g. corpora, lexicon, and
pre-trained models). In recent years, a number of studies have been published
on other languages such as Arabic Mubarak et al. (2020), Danish Sigurbergsson
and Derczynski (2020), Dutch Tulkens et al. (2016), French Chiril et al.
(2019), Greek Pitenis et al. (2020), Italian Poletto et al. (2017), Portuguese
Fortuna et al. (2019), Slovene Fišer et al. (2017), and Turkish Çöltekin
(2020) creating new datasets and resources for these languages.

Recent competitions organized in 2020 such as TRAC Kumar et al. (2020) and
OffensEval Zampieri et al. (2020) have included datasets in multiple languages
providing participants with the opportunity to explore cross-lingual learning
models opening exciting new avenues for research on languages other than
English and, in particular, on low-resource languages. The aforementioned deep
learning methods require large annotated datasets to perform well which is not
always available for low-resource languages. In this paper, we address the
problem of data scarcity in offensive language identification by using
transfer learning and cross-lingual transformers from a resource rich language
like English to three other languages: Bengali, Hindi, and Spanish.

##  3 Data

We acquired datasets in English and three other languages: Bengali, Hindi, and
Spanish (listed in Table 1). The four datasets have been used in shared tasks
in 2019 and 2020 allowing us to compare the performance of our methods to
other approaches. As our English dataset, we chose the Offensive Language
Identification Dataset (OLID) Zampieri et al. (2019a), used in the
SemEval-2019 Task 6 (OffensEval) Zampieri et al. (2019b). OLID is arguably one
of the most popular offensive language datasets. It contains manually
annotated tweets with the following three-level taxonomy and labels:

  * A:

Offensive language identification - offensive vs. non-offensive;

  * B:

Categorization of offensive language - targeted insult or thread vs.
untargeted profanity;

  * C:

Offensive language target identification - individual vs. group vs. other.

We chose OLID due to the flexibility provided by its hierarchical annotation
model that considers multiple types of offensive content in a single taxonomy
(e.g. targeted insults to a group are often hate speech whereas targeted
insults to an individual are often cyberbulling). This allows us to map OLID
level A (offensive vs. non-offensive) to labels in the other three datasets.
OLID’s annotation model is intended to serve as a general-purpose model for
multiple abusive language detection sub-tasks Waseem et al. (2017). The
transfer learning strategy used in this paper provides us with an interesting
opportunity to evaluate how closely the OLID labels relate to the classes in
datasets annotated using different guidelines and sub-task definitions (e.g.
aggression and hate speech).

The Hindi dataset Mandl et al. (2019) was used in the HASOC 2019 shared task,
while the Spanish dataset Basile et al. (2019) was used in SemEval-2019 Task 5
(HatEval). They both contain Twitter data and two labels. The Bengali dataset
Bhattacharya et al. (2020) was used in the TRAC-2 shared task Kumar et al.
(2020) on aggression identification. It is different than the other three
datasets in terms of domain (Facebook instead of Twitter) and set of labels
(three classes instead of binary), allowing us to compare the performance of
cross-lingual embeddings on off-domain data and off-task data.

Lang. | Inst. | S |  Labels
---|---|---|---
Bengali | 4,000 | F |  overtly aggressive, covertly aggressive, non aggressive
English | 14,100 | T |  offensive, non-offensive
Hindi | 8,000 | T |  hate offensive, non hate-offensive
Spanish | 6,600 | T |  hateful, non-hateful

Table 1: Instances (Inst.), source (S) and labels in all datasets. F stands
for Facebook and T for Twitter.

##  4 Methodology

Transformer models have been used successfully for various NLP tasks Devlin et
al. (2019). Most of the tasks were focused on English language due to the fact
the most of the pre-trained transformer models were trained on English data.
Even though, there were several multilingual models like BERT-m Devlin et al.
(2019) there was much speculations about its ability to represent all the
languages Pires et al. (2019) and although BERT-m model showed some cross-
lingual characteristics it has not been trained on crosslingual data
Karthikeyan et al. (2020). The motivation behind this methodology was the
recently released cross-lingual transformer models - XLM-R Conneau et al.
(2019) which has been trained on 104 languages. The interesting fact about
XLM-R is that it is very compatible in monolingual benchmarks while achieving
the best results in cross-lingual benchmarks at the same time Conneau et al.
(2019). The main idea of the methodology is that we train a classification
model on a resource rich, typically English, using a cross-lingual transformer
model, save the weights of the model and when we initialise the training
process for a lower resource language, start with the saved weights from
English. This process is also known as transfer learning and is illustrated in
Figure 1.

Figure 1: Transfer learning strategy.

There are two main parts of the methodology. Subsection 4.1 describes the
classification architecture we used for all the languages. In Subsection 4.2
we describe the transfer learning strategies used to take advantage of English
offensive language data in predicting offense in less-resourced languages.

###  4.1 XLM-R for Text Classification

Similar to other transformer architectures XLM-R transformer architecture can
also be used for text classification tasks Conneau et al. (2019). XLM-R-large
model contains approximately 125M parameters with 12-layers, 768 hidden-
states, 3072 feed-forward hidden-states and 8-heads Conneau et al. (2019). It
takes an input of a sequence of no more than 512 tokens and outputs the
representation of the sequence. The first token of the sequence is always
[CLS] which contains the special classification embedding Sun et al. (2019).

For text classification tasks, XLM-R takes the final hidden state h of the
first token [CLS] as the representation of the whole sequence. A simple
softmax classifier is added to the top of XLM-R to predict the probability of
label c: as shown in Equation 1 where W is the task-specific parameter matrix.

$p(c|\textbf{h})=softmax(W\textbf{h})$ (1)

We fine-tune all the parameters from XLM-R as well as W jointly by maximising
the log-probability of the correct label. The architecture diagram of the
classification is shown in Figure 2. We specifically used the XLM-R large
model.

Figure 2: Text classification architecture with XLM-R.

###  4.2 Transfer-learning strategies

When we adopt XLM-R for multilingual offensive language identification, we
perform transfer learning in two different ways.

#### Inter-language transfer learning

We first trained the XLM-R classification model on first level of English
offensive language identification dataset (OLID) Zampieri et al. (2019a). Then
we save the weights of the XLM-R model as well as the softmax layer. We use
these saved weights from English to initialise the weights for a new language.
To explore this transfer learning aspect we experimented on Hindi language
which was released for HASOC 2019 shared task Mandl et al. (2019) and on
Spanish data released for Hateval 2019 Basile et al. (2019).

#### Inter-task and inter-language transfer learning

Similar to the inter-language transfer learning strategy, we first trained the
XLM-R classification model on the first level of English offensive language
identification dataset (OLID) Zampieri et al. (2019a). Then we only save the
weights of the XLM-R model and use the saved weights to initialise the weights
for a new language. We did not use the weights of the last softmax layer since
we wanted to test this strategy on data that has a different number of
offensive classes to predict. We explored this transfer learning aspect with
Bengali dataset released with TRAC - 2 shared task Kumar et al. (2020). As
described in the Section 3 the classifier should make a 3-way classification
in between ‘Overtly Aggressive’, ‘Covertly Aggressive’ and ‘Non Aggressive’
text data.

We used a Nvidia Tesla K80 GPU to train the models. We divided the dataset
into a training set and a validation set using 0.8:0.2 split on the dataset.
We predominantly fine tuned the learning rate and number of epochs of the
classification model manually to obtain the best results for the validation
set. We obtained $1e^{-}5$ as the best value for learning rate and 3 as the
best value for number of epochs for all the languages. Training for English
language took around 1 hour while training for other languages took around 30
minutes. The code and the pretrained English model is available on
GitHub.111The public GitHub repository are available on
https://github.com/tharindudr/DeepOffense

##  5 Results and Evaluation

We evaluate the results obtained by all models using the test sets provided by
the organizers of each competition. We compared our results to the best
systems in TRAC-2 for Bengali, HASOC for Hindi, HatEval for Spanish in terms
of weighted and macro F1 score according to the metrics reported by the task
organizers - TRAC-2 reported only macro F1, HatEval reported only weighted F1,
and HASOC reported both.

Finally, we evaluate the improvement of the transfer learning strategy in the
performance of both BERT and XLM-R. We present the results along with the
majority class baseline for each language in Table 2. TL indicates that the
model used the inter language transfer learning strategy described in
Subsection 4.2.

Language | Model | M F1 | W F1
---|---|---|---
| XLM-R (TL) | 0.8415 | 0.8423
| Risch and Krestel (2020) | 0.8219 |
Bengali | BERT-m (TL) | 0.8197 | 0.8231
| XLM-R | 0.8142 | 0.8188
| BERT-m | 0.8132 | 0.8157
| Baseline | 0.2498 | 0.4491
| XLM-R (TL) | 0.8568 | 0.8580
| BERT-m (TL) | 0.8211 | 0.8220
Hindi | Bashar and Nayak (2019) | 0.8149 | 0.8202
| XLM-R | 0.8061 | 0.8072
| BERT-m | 0.8025 | 0.8030
| Baseline | 0.3510 | 0.3798
| XLM-R (TL) | 0.7513 | 0.7591
| BERT-m (TL) | 0.7319 | 0.7385
| Vega et al. (2019) |  | 0.7300
Spanish | Pérez and Luque (2019) |  | 0.7300
| XLM-R | 0.7224 | 0.7265
| BERT-m | 0.7215 | 0.7234
| Baseline | 0.3700 | 0.4348

Table 2: Results ordered by macro (M) F1 for Bengali and weighted (W) F1 for
Hindi and Spanish.

For Hindi, transfer learning with XLM-R cross lingual embeddings provided the
best results achieving 0.8568 and 0.8580 weighted and macro F1 score
respectively. In HASOC 2019 Mandl et al. (2019), the best model by Bashar and
Nayak (2019) scored 0.8149 Macro F1 and 0.8202 Weighted F1 using convolutional
neural networks.

For Spanish transfer learning with XLM-R cross lingual embeddings also
provided the best results achieving 0.7513 and 0.7591 macro and weighted F1
score respectively. The best two models in HatEval Basile et al. (2019) for
Spanish scored 0.7300 macro F1 score. Both models applied SVM classifiers
trained on a variety of features like character and word n-grams, POS tags,
offensive word lexicons, and embeddings.

The results for Bengali deserve special attention because the Bengali data is
off-domain with respect to the English data (Facebook instead of Twitter), and
it contains three labels (covertly aggressive, overtly aggressive, and not
aggressive) instead of the two labels present in the English dataset
(offensive and non-offensive). TL indicates that the model used the inter-
task, inter-domain, and inter-language transfer learning strategy described in
Subsection 4.2. Similar to the Hindi and Spanish, transfer learning with XLM-R
cross lingual embeddings provided the best results for Bengali achieving
0.7513 and 0.7591 macro and weighted F1 respectively thus outperforming the
other models by a significant margin. The best model in the TRAC-2 shared task
Kumar et al. (2020) scored 0.821 weighted F1 score in Bengali using a BERT-
based system.

We look closer at the test set predictions by XLM-R (TL) for Bengali in Figure
3. We observe that the performance for the non-aggressive class is
substantially better than the performance for the overtly aggressive and
covertly aggressive classes following a trend observed by the TRAC-2
participants including Risch and Krestel (2020).

Figure 3: Heat map of the Bengali test set predictions by XLM-R (TL).

Finally, it is clear that in all the experimental settings, the cross-lingual
embedding models fine-tuned with transfer learning, outperforms the best
system available for the three languages. Furthermore, the results show that
the cross-lingual nature of the XLM-R model provided a boost over the
multilingual BERT model in all languages tested.

##  6 Conclusion

This paper is the first study to apply cross-lingual contextual word
embeddings in offensive language identification projecting predictions from
English to other languages using benchmarked datasets from shared tasks on
Bengali Kumar et al. (2020), Hindi Mandl et al. (2019), and Spanish Basile et
al. (2019). We have showed that XLM-R with transfer learning outperforms all
of the other methods we tested as well as the best results obtained by
participants of the three competitions.

The results obtained by our models confirm that OLID’s general hierarchical
annotation model encompasses multiple types of offensive content such as
aggression, included in the Bengali dataset, and hate speech included in the
Hindi and Spanish datasets, allowing us to model different sub-tasks jointly
using the methods described in this paper. Furthermore, the results we
obtained for Bengali show that it is possible to achieve high performance
using transfer learning on off-domain (Twitter vs. Facebook) and off-task data
when the labels do not have a direct correspondence in the projected dataset
(two in English and three in Bengali). This opens exciting new avenues for
future research considering the multitude of phenomena (e.g. hate speech,
aggression, cyberbulling), annotation schemes and guidelines used in offensive
language datasets.

In future work, we would like to further evaluate our models using SOLID, a
novel large English dataset with over 9 million tweets Rosenthal et al.
(2020), along with datasets in four other languages (Arabic, Danish, Greek,
and Turkish) that were made available for the second edition of OffensEval
Zampieri et al. (2020). These datasets were collected using the same
methodology and were annotated according to OLID’s guidelines. Finally, we
would also like to apply our models to languages with even less resources
available to help coping with the problem of offensive language in social
media.

## Acknowledgement

We would like to thank the HASOC, HatEval, and OffensEval shared task
organizers for making the datasets used in this paper available.

We further thank the anonymous EMNLP reviewers for providing us with valuable
feedback to improve this paper.

## References

  * Aroyehun and Gelbukh (2018) Segun Taofeek Aroyehun and Alexander Gelbukh. 2018.  Aggression detection in social media: Using deep neural networks, data augmentation, and pseudo labeling.  In _Proceedings of TRAC_.
  * Bannink et al. (2014) Rienke Bannink, Suzanne Broeren, Petra M van de Looij-Jansen, Frouwkje G de Waart, and Hein Raat. 2014.  Cyber and Traditional Bullying Victimization as a Risk Factor for Mental Health Problems and Suicidal Ideation in Adolescents.  _PloS one_ , 9(4).
  * Bashar and Nayak (2019) Md Abul Bashar and Richi Nayak. 2019.  QutNocturnal@ HASOC’19: CNN for hate speech and offensive content identification in Hindi language.  In _Proceedings of FIRE_.
  * Basile et al. (2019) Valerio Basile, Cristina Bosco, Elisabetta Fersini, Debora Nozza, Viviana Patti, Francisco Manuel Rangel Pardo, Paolo Rosso, and Manuela Sanguinetti. 2019\.  SemEval-2019 task 5: Multilingual detection of hate speech against immigrants and women in twitter.  In _Proceedings of SemEval_.
  * Bhattacharya et al. (2020) Shiladitya Bhattacharya, Siddharth Singh, Ritesh Kumar, Akanksha Bansal, Akash Bhagat, Yogesh Dawer, Bornini Lahiri, and Atul Kr. Ojha. 2020.  Developing a multilingual annotated corpus of misogyny and aggression.  In _Proceedings of TRAC_.
  * Bonanno and Hymel (2013) Rina A Bonanno and Shelley Hymel. 2013.  Cyber bullying and internalizing difficulties: Above and beyond the impact of traditional forms of bullying.  _Journal of youth and adolescence_ , 42(5):685–697.
  * Çöltekin (2020) Çağrı Çöltekin. 2020.  A Corpus of Turkish Offensive Language on Social Media.  In _Proceedings of LREC_.
  * Chiril et al. (2019) Patricia Chiril, Farah Benamara Zitoune, Véronique Moriceau, Marlène Coulomb-Gully, and Abhishek Kumar. 2019.  Multilingual and multitarget hate speech detection in tweets.  In _Proceedings of TALN_.
  * Conneau et al. (2019) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Unsupervised cross-lingual representation learning at scale.  _arXiv preprint arXiv:1911.02116_.
  * Dadvar et al. (2013) Maral Dadvar, Dolf Trieschnigg, Roeland Ordelman, and Franciska de Jong. 2013.  Improving Dyberbullying Detection with User Context.  In _Proceedings of ECIR_.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  Bert: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of NAACL_.
  * Fišer et al. (2017) Darja Fišer, Tomaž Erjavec, and Nikola Ljubešić. 2017.  Legal Framework, Dataset and Annotation Schema for Socially Unacceptable On-line Discourse Practices in Slovene.  In _Proceedings ALW_.
  * Fortuna et al. (2019) Paula Fortuna, Joao Rocha da Silva, Leo Wanner, Sérgio Nunes, et al. 2019.  A Hierarchically-labeled Portuguese Hate Speech Dataset.  In _Proceedings of ALW_.
  * Ghadery and Moens (2020) Erfan Ghadery and Marie-Francine Moens. 2020.  Liir at semeval-2020 task 12: A cross-lingual augmentation approach for multilingual offensive language identification.  _arXiv preprint arXiv:2005.03695_.
  * Karthikeyan et al. (2020) K Karthikeyan, Zihan Wang, Stephen Mayhew, and Dan Roth. 2020.  Cross-lingual ability of multilingual bert: An empirical study.  In _Proceedings of ICLR_.
  * Kumar et al. (2018) Ritesh Kumar, Atul Kr Ojha, Shervin Malmasi, and Marcos Zampieri. 2018.  Benchmarking aggression identification in social media.  In _Proceedings of TRAC_.
  * Kumar et al. (2020) Ritesh Kumar, Atul Kr. Ojha, Shervin Malmasi, and Marcos Zampieri. 2020.  Evaluating Aggression Identification in Social Media.  In _Proceedings of TRAC_.
  * Liu et al. (2019) Ping Liu, Wen Li, and Liang Zou. 2019.  NULI at SemEval-2019 task 6: Transfer learning for offensive language detection using bidirectional transformers.  In _Proceedings of SemEval_.
  * Majumder et al. (2018) Prasenjit Majumder, Thomas Mandl, et al. 2018.  Filtering Aggression from the Multilingual Social Media Feed.  In _Proceedings TRAC_.
  * Malmasi and Zampieri (2017) Shervin Malmasi and Marcos Zampieri. 2017.  Detecting Hate Speech in Social Media.  In _Proceedings of RANLP_.
  * Malmasi and Zampieri (2018) Shervin Malmasi and Marcos Zampieri. 2018.  Challenges in Discriminating Profanity from Hate Speech.  _Journal of Experimental & Theoretical Artificial Intelligence_, 30:1 – 16.
  * Mandl et al. (2019) Thomas Mandl, Sandip Modha, Prasenjit Majumder, Daksh Patel, Mohana Dave, Chintak Mandlia, and Aditya Patel. 2019.  Overview of the hasoc track at fire 2019: Hate speech and offensive content identification in indo-european languages.  In _Proceedings of FIRE_.
  * Mubarak et al. (2017) Hamdy Mubarak, Darwish Kareem, and Magdy Walid. 2017.  Abusive language detection on Arabic social media.  In _Proceedings of ALW_.
  * Mubarak et al. (2020) Hamdy Mubarak, Ammar Rashed, Kareem Darwish, Younes Samih, and Ahmed Abdelali. 2020\.  Arabic offensive language on twitter: Analysis and experiments.  _arXiv preprint arXiv:2004.02192_.
  * Pamungkas and Patti (2019) Endang Wahyu Pamungkas and Viviana Patti. 2019.  Cross-domain and cross-lingual abusive language detection: A hybrid approach with deep learning and a multilingual lexicon.  In _Proceedings ACL:SRW_.
  * Pérez and Luque (2019) Juan Manuel Pérez and Franco M Luque. 2019.  Atalaya at semeval 2019 task 5: Robust embeddings for tweet classification.  In _Proceedings of SemEval_.
  * Peters et al. (2018) Matthew Peters, Mark Neumann, Mohit Iyyer, Matt Gardner, Christopher Clark, Kenton Lee, and Luke Zettlemoyer. 2018.  Deep Contextualized Word Representations.  In _Proceedings of NAACL_.
  * Pires et al. (2019) Telmo Pires, Eva Schlinger, and Dan Garrette. 2019.  How multilingual is multilingual BERT?  In _Proceedings of ACL_.
  * Pitenis et al. (2020) Zeses Pitenis, Marcos Zampieri, and Tharindu Ranasinghe. 2020.  Offensive Language Identification in Greek.  In _Proceedings of LREC_.
  * Poletto et al. (2017) Fabio Poletto, Marco Stranisci, Manuela Sanguinetti, Viviana Patti, and Cristina Bosco. 2017.  Hate Speech Annotation: Analysis of an Italian Twitter Corpus.  In _Proceedings of CLiC-it_.
  * Ranasinghe et al. (2019) Tharindu Ranasinghe, Marcos Zampieri, and Hansi Hettiarachchi. 2019.  Brums at hasoc 2019: Deep learning models for multilingual hate speech and offensive language identification.  In _Proceedings of FIRE_.
  * Ridenhour et al. (2020) Michael Ridenhour, Arunkumar Bagavathi, Elaheh Raisi, and Siddharth Krishnan. 2020\.  Detecting Online Hate Speech: Approaches Using Weak Supervision and Network Embedding Models.  _arXiv preprint arXiv:2007.12724_.
  * Risch and Krestel (2020) Julian Risch and Ralf Krestel. 2020.  Bagging bert models for robust aggression identification.  In _Proceedings of TRAC_.
  * Rosa et al. (2019) Hugo Rosa, N Pereira, Ricardo Ribeiro, Paula Costa Ferreira, Joao Paulo Carvalho, S Oliveira, Luísa Coheur, Paula Paulino, AM Veiga Simão, and Isabel Trancoso. 2019.  Automatic cyberbullying detection: A systematic review.  _Computers in Human Behavior_ , 93:333–345.
  * Rosenthal et al. (2020) Sara Rosenthal, Pepa Atanasova, Georgi Karadzhov, Marcos Zampieri, and Preslav Nakov. 2020.  A Large-Scale Weakly Supervised Dataset for Offensive Language Identification.  In _arXiv preprint arXiv:2004.14454_.
  * Sigurbergsson and Derczynski (2020) Gudbjartur Ingi Sigurbergsson and Leon Derczynski. 2020.  Offensive Language and Hate Speech Detection for Danish.  In _Proceedings of LREC_.
  * Sun et al. (2019) Chi Sun, Xipeng Qiu, Yige Xu, and Xuanjing Huang. 2019.  How to fine-tune bert for text classification?  In _Chinese Computational Linguistics_ , pages 194–206.
  * Tulkens et al. (2016) Stéphan Tulkens, Lisa Hilte, Elise Lodewyckx, Ben Verhoeven, and Walter Daelemans. 2016.  A Dictionary-based Approach to Racism Detection in Dutch Social Media.  In _Proceedings of TA-COS_.
  * Vega et al. (2019) Luis Enrique Argota Vega, Jorge Carlos Reyes-Magaña, Helena Gómez-Adorno, and Gemma Bel-Enguix. 2019.  Mineriaunam at semeval-2019 task 5: Detecting hate speech in twitter using multiple features in a combinatorial framework.  In _Proceedings of SemEval_.
  * Waseem et al. (2017) Zeerak Waseem, Thomas Davidson, Dana Warmsley, and Ingmar Weber. 2017.  Understanding Abuse: A Typology of Abusive Language Detection Subtasks.  In _Proceedings of ALW_.
  * Xu et al. (2012) Jun-Ming Xu, Kwang-Sung Jun, Xiaojin Zhu, and Amy Bellmore. 2012.  Learning from bullying traces in social media.  In _Proceedings of NAACL_.
  * Yao et al. (2019) Mengfan Yao, Charalampos Chelmis, and Daphney-Stavroula Zois. 2019.  Cyberbullying Ends Here: Towards Robust Detection of Cyberbullying in Social Media.  In _Proceedings of WWW_.
  * Zampieri et al. (2019a) Marcos Zampieri, Shervin Malmasi, Preslav Nakov, Sara Rosenthal, Noura Farra, and Ritesh Kumar. 2019a.  Predicting the type and target of offensive posts in social media.  In _Proceedings of NAACL_.
  * Zampieri et al. (2019b) Marcos Zampieri, Shervin Malmasi, Preslav Nakov, Sara Rosenthal, Noura Farra, and Ritesh Kumar. 2019b.  SemEval-2019 Task 6: Identifying and Categorizing Offensive Language in Social Media (OffensEval).  In _Proceedings of SemEval_.
  * Zampieri et al. (2020) Marcos Zampieri, Preslav Nakov, Sara Rosenthal, Pepa Atanasova, Georgi Karadzhov, Hamdy Mubarak, Leon Derczynski, Zeses Pitenis, and Cagrı Cöltekin. 2020.  Semeval-2020 task 12: Multilingual offensive language identification in social media (offenseval 2020).  _Proceedings of SemEval_.


Alexander Moriggl, University of Innsbruck Techniker. 13, A-6020 Innsbruck,
Austria

# Semi-Lagrangian 4d, 5d, and 6d kinetic plasma simulation on large scale GPU
equipped supercomputers

Lukas Einkemmer11affiliationmark:  and Alexander Moriggl11affiliationmark:
11affiliationmark:  University Of Innsbruck, Austria
<EMAIL_ADDRESS>
###### Abstract

Running kinetic plasma physics simulations using grid-based solvers is very
demanding both in terms of memory as well as computational cost. This is
primarily due to the up to six-dimensional phase space and the associated
unfavorable scaling of the computational cost as a function of grid spacing
(often termed the curse of dimensionality). In this paper, we present 4d, 5d,
and 6d simulations of the Vlasov–Poisson equation with a split-step semi-
Lagrangian discontinuous Galerkin scheme on graphic processing units (GPUs).
The local communication pattern of this method allows an efficient
implementation on large-scale GPU-based systems and emphasizes the importance
of considering algorithmic and high-performance computing aspects in unison.
We demonstrate a single node performance above 2 TB/s effective memory
bandwidth (on a node with 4 A100 GPUs) and show excellent scaling (parallel
efficiency between 30% and 67%) for up to 1536 A100 GPUs on JUWELS Booster.

###### keywords:

fully kinetic simulation, general purpose computing on graphic processing
units (GPGPU), Vlasov–Poisson equation, semi-Lagrangian discontinuous Galerkin
method

##  1 Introduction

Computational fluid dynamics has become an extremely successful framework to
understand a range of natural phenomena ranging from astrophysics to everyday
technological applications. Despite these successes, it is increasingly
realized that many phenomena in plasma physics require a kinetic description
Howes et al. (2008); Zweibel and Yamada (2009). The most challenging aspect of
kinetic models, from a computational point of view, is certainly the up to
six-dimensional phase space. The traditional approach to overcome this
difficulty has been to employ so-called particle methods Verboncoeur (2005);
Hariri et al. (2016). While these simulations are comparably cheap and have
been successfully employed in some applications, it is well known that they
miss or do not resolve accurately certain physical phenomena, such as Landau
damping. In addition, the inherent numerical noise introduced makes it
extremely difficult to obtain accurate results for many physical phenomena.
More recently, sparse grids Kormann and Sonnendrücker (2016); Guo and Cheng
(2016) and low-rank methods Einkemmer and Lubich (2018); Einkemmer and Joseph
(2021) have been used to reduce memory requirements and computational
complexity. Nevertheless, the need to run fully kinetic simulations in high
dimension persists for a range of applications.

To run such simulations large supercomputers are needed and codes to
facilitate this for CPU based systems are available. See, e.g., Bigot et al.
(2013) for a scaling study of the 5-dimensional GYSELA code. More recently,
also simulations for the full 6 dimensional problem have been considered
Kormann et al. (2019). In the latter work scaling to up to 65,536 processes
has been reported on the SuperMUC HPC system of the Leibnitz Supercomputing
Center. We also note that the amount of communication required usually
increases as we increase the dimensionality of the problem. The reason for
this is that the ratio of boundary points (that need to be transferred) to
interior grid points increases significantly.

To numerically discretize these models, semi-Lagrangian schemes have proved
popular. By using information on the characteristics of the underlying partial
differential equations (PDEs), these methods do not suffer from the usual step
size restriction of explicit methods. A particular emphasize in collisionless
models, as we will consider here, is to introduce as little numerical
diffusion as possible into the approximation. It is well known that semi-
Lagrangian schemes based on spline interpolation (see, e.g., Sonnendrücker et
al. (1999) and Filbet and Sonnendrücker (2003) for a comparison of different
methods) and spectral methods (see, e.g., Klimas and Farrell (1994); Klimas
and Viñas (2018); Camporeale et al. (2016)) perform well according to this
metric.

However, the results for scaling to large supercomputers obtained in Kormann
et al. (2019) use Lagrange interpolation instead. The reason for this is that
in a massively parallel environment, the global communication patterns induced
by spline interpolation or spectral methods provide a significant impediment
to scalability. This, in particular, is also a serious issue on GPU-based
systems (even on a single node; see, e.g., Einkemmer (2020)). Recently, it has
been demonstrated that a different type of semi-Lagrangian scheme, the local
semi-Lagrangian discontinuous Galerkin approach, can be implemented
efficiently on both Tesla and consumer GPUs, while maintaining or even
exceeding the performance of spline-based semi-Lagrangian schemes Einkemmer
(2020).

In this paper, we present results, using our SLDG code, that demonstrates that
large-scale kinetic simulations of the Vlasov–Poisson equation can be made to
scale efficiently on modern GPU-based supercomputers. This is particularly
relevant as most (if not all) pre-exascale and exascale systems will make use
of GPUs and kinetic simulations are prime candidates that could exploit such
systems.

##  2 Problem and algorithmic background

###  2.1 Vlasov–Poisson equation

To describe the kinetic dynamics of a plasma we use the Vlasov–Poisson
equations

$\displaystyle\partial_{t}f+v\cdot\nabla_{x}f-E(f)\cdot\nabla_{v}f=0,$ (1)
$\displaystyle E=-\nabla_{x}\phi,\quad-\Delta\phi=1-\rho,\quad\rho=\int
f\,dv,$ (2)

where $f(t,x,v)$, with $(x,v)\in\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{v}}$,
is the (up to) six-dimensional particle density in phase space, $\rho(t,x)$
the (up to) three-dimensional charge density and $E(t,x)$ is the (up to)
three-dimensional electric field. The Vlasov–Poisson equation is a hyperbolic
model that couples a transport equation self-consistently to an equation for
the electric field. The electric field is determined via a Poisson problem
from the charge density $\rho$.

The particle density is up to six-dimensional. Depending on the physical model
numerical simulations can be run using different number of dimensions in space
(denoted by $d_{x}$) and in velocity (denoted by $d_{v}$). In this paper we
consider problems with 2x2v ($d_{x}=2$ and $d_{v}=2$, four dimensional), 2x3v
($d_{x}=2$ and $d_{v}=3$, five dimensional), and 3x3v ($d_{x}=3$ and
$d_{v}=3$, six dimensional).

While we restrict ourselves to the Vlasov–Poisson model in this work, we
emphasize that the performance results obtained provide a good indication for
other kinetic models as well. Both on the plasma physics side (such as in case
of the Vlasov–Maxwell equations) as well as for kinetic problems from
radiative transport.

###  2.2 Time splitting

The first step in our algorithm for solving the Vlasov–Poisson equation is to
perform a time splitting. The approach we use here has been introduced in the
seminal paper Cheng and Knorr (1976) and has subsequently been generalized to,
e.g., the Vlasov–Poisson equation Crouseilles et al. (2015).

The idea of the splitting approach is to decompose the nonlinear
Vlasov–Poisson equation into a sequence of simpler and linear steps. The
strategy to compute the numerical solution at time $t^{n+1}=t^{n}+\Delta t$,
where we use $f^{n}\approx f(t_{n},x,v)$, is as follows

  * 1.

Solve $\partial_{t}f(t,x,v)+v\cdot\nabla_{x}f(t,x,v)=0$ with initial value
$f(0,x,v)=f^{n}(x,v)$ to obtain $f^{\star}(x,v)=f(\Delta t/2,x,v)$.

  * 2.

Compute $E$ using $f^{\star}$ by solving the Poisson problem (2).

  * 3.

Solve $\partial_{t}f(t,x,v)-E(x)\cdot\nabla_{v}f(t,x,v)=0$ with initial value
$f(0,x,v)=f^{\star}(x,v)$ to obtain $f^{\star\star}(x,v)=f(\Delta t,x,v)$.

  * 4.

Solve $\partial_{t}f(t,x,v)+v\cdot\nabla_{x}f(t,x,v)=0$ with initial value
$f(0,x,v)=f^{\star\star}(t,x,v)$ to obtain $f^{n+1}(x,v)=f(\Delta t/2,x,v)$.

This results in a second-order scheme (i.e. Strang splitting). Since the
advection speed never depends on the direction of the transport, the
characteristics can be determined analytically. For example, for step 1. we
have

$f^{\star}(x,v)=f^{n}(x-v\Delta t/2,v).$

This advection in $d_{x}$ dimensions can be further reduced to one-dimensional
advections by the following splitting procedure

  * 1a.

$f^{\star a}(x_{1},x_{2},x_{3},v)=f^{n}(x_{1}-v_{1}\Delta t/2,x_{2},x_{3},v)$.

  * 1b.

$f^{\star b}(x_{1},x_{2},x_{3},v)=f^{\star a}(x_{1},x_{2}-v_{2}\Delta
t/2,x_{3},v)$.

  * 1c.

$f^{\star}(x_{1},x_{2},x_{3},v)=f^{\star b}(x_{1},x_{2},x_{3}-v_{3}\Delta
t/2,v)$.

We emphasize that no error is introduced by this splitting as the
corresponding operators commute. We apply a similar procedure to step 3, which
results in the following

  * 3a.

$f^{\star a}(x,v_{1},v_{2},v_{3})=f^{n}(x,v_{1}+E_{1}(x)\Delta
t,v_{2},v_{3})$.

  * 3b.

$f^{\star b}(x,v_{1},v_{2},v_{3})=f^{\star a}(x,v_{1},v_{2}+E_{2}(x)\Delta
t,v_{3})$.

  * 3c.

$f^{\star}(x,v_{1},v_{2},v_{3})=f^{\star b}(x,v_{1},v_{2},v_{3}+E_{3}(x)\Delta
t)$.

Since the advection in the $v$-direction (step 3a-3c) does not change the
electric field, we have reduced the up to six-dimensional Vlasov–Poisson
equation to a sequence of one-dimensional linear advections and a single
Poisson solve. The scheme outlined here is second-order accurate Einkemmer and
Ostermann (2014c). High-order methods have been constructed as well Casas et
al. (2017) and generalization of this basic approach to more complicated
problems also exist (see, e.g., Crouseilles et al. (2015) for the
Vlasov–Maxwell equations).

The advections are treated using a so-called semi-Lagrangian method, which
evaluates the function approximation at the foot of the characteristic curves.
Since these end points of the characteristics do not lie on the grid points,
an interpolation technique has to be applied. We will describe this procedure
in detail in the next section.

###  2.3 Semi-Lagrangian discontinuous Galerkin method

Semi-Lagrangian methods have been widely used to solve the Vlasov equation. In
particular, spline interpolation Sonnendrücker et al. (1999); Filbet et al.
(2001) has emerged as a popular method as it does not introduce too much
numerical diffusion, is conservative, and gives good accuracy at a reasonable
computational cost.

However, in the construction of a spline each degree of freedom couples with
each other degree of freedom. Due to this global data dependency, spline
interpolation is not suited to distributed memory parallelism. To remedy this
in Crouseilles et al. (2009) a method that performs spline interpolation on
patches has been proposed. However, on today’s massively parallel systems with
relatively few degrees of freedom per direction on each MPI process
(especially in the six-dimensional case), it is unclear how much advantage
this approach provides. In this context, it should be noted that the massively
parallel implementation found in Kormann et al. (2019) forgoes using spline
interpolation and focuses exclusively on Lagrange interpolation. The same is
true for the GPU based implementation in Mehrenberger et al. (2013). While
Lagrange interpolation scales well it is also known to be very diffusive and
usually large stencils are required in order to obtain accurate results.

In the present work, we will use a semi-Lagrangian discontinuous Galerkin
scheme (sldg). This method, see e.g. Crouseilles et al. (2011); Rossmanith and
Seal (2011), is conservative by construction Einkemmer (2017), has similar or
better properties with respect to limiting diffusion compared to spline
interpolation (see Einkemmer (2019) for a comparison), and is completely local
(at most two neighboring cells are accessed to update a cell). Thus, it
provides a method that in addition to being ideally suited for today’s
massively parallel CPU and GPU systems, is competitive from a numerical point
of view. In addition, a well-developed convergence analysis is available for
the Vlasov–Poisson equation Einkemmer and Ostermann (2014c, b).

The sldg scheme computes in an efficient way $u^{n+1}(y)=u^{n}(y-at)$. The
one-dimensional advections in the previous section are written precisely in
this form (all the remaining variables are considered as parameters here).
This scheme first divides the domain into cells. The restriction of the
solution to each such cell is a polynomial of degree $k$. Note that there is
no continuity constraint across cell interfaces and the approximant is thus
discontinuous across this interface (see the left picture in figure 1). The
basis to form these polynomials are Lagrange polynomials which interpolate at
the Gauss–Legendre points scaled to each cell. Thus, the degrees of freedom
are function values $u_{ij}^{n}$ where the index $i$ is the cell index and the
index $j$ refers to the $j^{th}$ Gaus–Legendre point of cell $i$. Since the
translated function $u^{n+1}$ in general does not lie in this approximation
space, an $L^{2}$ projection is applied, see figure 1. Only two adjacent cells
are required in order to compute the values at a given cell, we denote the
index of them with $i^{*}$ and $i^{*}+1$. The resulting numerical scheme has
order $k+1$ and can be written as follows

$u^{n+1}_{ij}=\sum_{l}A_{jl}u^{n}_{i^{*}l}+\sum_{l}B_{jl}u^{n}_{i^{*}+1;l}.$
(3)

The matrices $A$ and $B$ are of size $(k+1)\times(k+1)$ and only depend on the
advection speed, which can be precomputed. More details about the specific
form of these two matrices can be found in Einkemmer (2019).

Figure 1: Illustrates the semi-Lagrangian discontinuous Galerkin method. When
the advection is complete (central picture), a projection has to be performed
which results in the red line in the right picture. This projection lies in
the approximation space and is used for subsequent computations. Note that
only two adjacent cells are required to compute the values for a certain cell.

###  2.4 Poisson solver

To determine the electric field, the Poisson equation (2) has to be solved
once in each time step. However, since this is only an up to three-dimensional
problem, the computational cost of this part of the algorithm is small. Using
FFT based methods is a common approach that is also supported in our SLDG
code.

When the spatial variables $x$ are partitioned over multiple compute nodes,
however, the FFT method requires global communication. Thus, we propose an
alternative method which requires only local communication to solve this
problem. Since the advection is already treated using a discontinuous Galerkin
scheme, we will use the same approximation space for the Poisson solver. This
has the added benefit that no interpolation to and from the equidistant grid
on which FFT based methods operate is necessary. The details of the
discontinuous Galerkin Poisson solver are outlined in the Appendix.

##  3 Implementation and parallelization

Our code is written in C++/CUDA and uses templates to treat problems with
different number of dimensions in physical and velocity space. This allows us
to only maintain one codebase, while still providing excellent performance in
all these configurations. Our code SLDG is available at
https://bitbucket.org/leinkemmer/sldg under the terms of the MIT license.

The main computational effort of the code is due to the following three parts:

Advection: In the advection step equation (3) is applied to each cell.This is
done in a CUDA kernel and the operation is largely memory bound on modern
GPUs. This step is required whether we run the code on a single GPU or on a
large supercomputer and its cost per degree of freedom is expected to be
independent of the number of MPI processes taking part in the simulation.

Communication: In order to run the code on multiple GPUs that possible span
multiple nodes a domain decomposition approach is used. Depending on the
problem, either only the velocity direction or both the space and velocity
directions are divided along $d_{x}$ or $d_{x}+d_{v}$ directions. Each of
these blocks is assigned to a single GPU. Performing an advection step
necessitates the exchange of boundary data from GPU to GPU. On a single node
this can be done using a call to cudaMemcpy, which makes use of NVLink if
available. Communication between nodes is done via MPI. Our implementation
supports the use of CUDA-Aware MPI which can make use of GPU memory in MPI
calls. Thus, in this configuration MPI routines can directly read from and
place data in GPU memory.

Especially for high-dimensional problems, this type of boundary exchange can
be relatively expensive. The reason for this is that the number of degrees of
freedom in each spatial and velocity direction that can be stored on a single
GPU is relatively small (compared to a three-dimensional problem). Thus,
boundary cells constitute a larger fraction of the overall degrees of freedom,
which due to the discrepancy between GPU memory and InfiniBand bandwidth means
that communication can take up a significant fraction of the total run time.
The data transfer is illustrated in figure 2. Let us emphasize, however, that
not the entire boundary of each data block associated to a GPU has to be
transferred to each neighboring GPU. Since advection is directional, depending
on the sign of the velocity field only either a send or a receive is required.
This is illustrated in figure 3. By taking this fact into account, only half
the boundary data has to be transferred.

E field: The computation of the electric field proceeds in the following
steps. First, the charge density $\rho(x)=\int f(x,v)\,dv$ has to be
determined. This is first done on each GPU (using the CUB library NVIDIA
(2018)) and then a reduction operation in the velocity direction is performed
by calling MPI_Reduce on processes which span the entire velocity domain but
represent the same part of the spatial domain (as illustrated in Figure 4).

At the completion of this step $\rho$ is distributed over certain MPI root
processes which span the entire space domain. These processes then participate
in computing the electric field $E$. This can be done either by an FFT based
method, which requires global communication among the processes, or an
iterative conjugate gradient discontinuous Galerkin (CGdG) solver which only
requires nearest neighbor exchange in each iteration (as illustrated in Figure
5).

Finally, the electric field has to be made available to all processes that
span the velocity domain. To do this MPI_Broadcast is called, see figure 6.

Although from a computational perspective the computation of the electric
field is relatively cheap (as it is only an up to three-dimensional problem),
the associated data transfer can become an impediment to scaling, especially
on large systems.

Figure 2: Illustration of the data transfer between GPUs in the 1+1
dimensional case (i.e. in the $x-v$ plane.The color orange represents the
ghost cell data that potentially (depending on the direction of the velocity)
has to be sent to neighboring GPUs. Black arrows indicate the corresponding
data transfer  Figure 3: Demonstration of the advection in the velocity
direction and the required boundary data (orange blocks) from neighboring
processes. The direction of data transfer depends on the sign of the velocity
(the blue function below the illustration). Figure 4: Illustration of the data
transfer necessary to compute the charge density $\rho$ (Step 1 in computing
the E field). Figure 5: Illustration of the data transfer required to compute
the electric field from the charge density $\rho$ (Step 2 in computing the E
field). Figure 6: Illustration of the data transfer required to distribute the
electric field to the appropriate processes (Step 3 in computing the E field).

In a multi-node context each MPI process is assigned one or multiple blocks.
Each block manages a single GPU. The configuration we employ on most clusters
with multiple GPUs on each node is to have one MPI process per node (i.e. one
MPI process per shared memory domain) and consequently as many blocks per MPI
process as GPUs are on a node. However, on Juwels Booster this significantly
degrades the bandwidth that is available to conduct MPI communication between
different nodes. Thus, the multi-node simulations in this paper are run using
a single MPI process per GPU. On Juwels Booster therefore 4 MPI processes are
launched on each node and CUDA-aware MPI is used both for inter-node as well
as intra-node communication.

###  3.1 Validation

To validate our implementation, two commonly used test scenarios are
considered. The linear Landau damping problem given by the initial value

$f(0,x,v)=\frac{\exp\left(-\|v\|^{2}/2\right)}{(2\pi)^{d_{v}/2}}\left(1+\alpha\sum_{i=1}^{d_{x}}\cos(kx_{i})\right)$

on $[0,4\pi]^{d_{x}}\times[-6,6]^{d_{v}}$, $\alpha=10^{-2}$, and $k=0.5$ and a
bump-on-tail instability

$\displaystyle f(0,x,v)=$
$\displaystyle\frac{1}{(2\pi)^{d_{v}/2}}\left(0.9\exp\left(-\frac{v_{1}^{2}}{2}\right)+0.2\exp\left(-2(v_{1}-4.5)^{2}\right)\right)$
$\displaystyle\exp\left(-\frac{\sum_{i=2}^{d_{v}}v_{i}^{2}}{2}\right)\left(1+0.03\sum_{i=1}^{d_{x}}\cos(0.3x_{i})\right).$

on $[0,20\pi/3]^{d_{x}}\times[-9,9]^{d_{v}}$.

In the simulations presented here we use the $4^{\text{th}}$ order semi-
Lagrangian discontinuous Galerkin scheme (as it represents a good compromise
between accuracy and computational cost Einkemmer (2017, 2020)) and a time
step size of $0.1$. In figures 7 and 8 the time evolution of the electric
energy for these initial conditions is shown for a five and six-dimensional
configuration, respectively. The results agree well with what has been
reported in the literature. In addition, the Landau damping results show
excellent agreement with the analytically obtained decay rate. There is no
significant difference between the results obtained with the FFT and
discontinuous Galerkin Poisson solver. The SLDG code also includes a set of
unit tests that check a range of results against known analytic and numerical
solutions and compares the results obtain on the GPU to the results on the
CPU.

Figure 7:  Time evolution of the electric energy for the linear Landau damping
(top) and bump-on-tail instability (bottom) for the 2x3v case. For both test
cases the $4^{\text{th}}$ order method with a time step of $0.1$ is used. For
the linear Landau test case $72^{2}144^{3}$ degrees of freedom are used (the
domain is distributed over 8 GPUs), while for the bump-on-tail test scenario
$144^{5}$ degrees of freedom are used (the domain is distributed over 32
GPUs).

Figure 8:  Time evolution of the electric energy of the linear Landau damping
(top) and bump-on-tail instability (bottom) for the 3x3v case. For both test
case the $4^{\text{th}}$ order method with a time step size of $0.1$ is used.
For the linear Landau test case $36^{3}72^{3}$ degrees of freedom are used
(the domain is distributed over 8 GPUs), while for the bump-on-tail test
scenario $72^{6}$ degrees of freedom are used (the domain is distributed over
64 GPUs).

We further remark that the results obtained using single and double precision
overlap in the plots. Thus if one is not interested in too stringent
tolerances running simulations with the semi-Lagrangian discontinuous Galerkin
scheme in single precision is a viable option. This reduces memory consumption
by a factor of two. In particular, for GPUs this can be attractive as memory,
in general, is more limited on GPU based systems. Thus, larger problem sizes
can be treated per GPU. In addition, for some problems consumer cards can
potentially be used, which do not have strong double-precision performance. We
will present and discuss both double and single precision results in the
remainder of this paper.

##  4 Single node performance

The results in this section are obtained simulating the linear Landau test
case. If not otherwise stated, a fixed time step size of $0.1$ and the
$4^{\text{th}}$ order semi-Lagrangian discontinuous Galerkin method is used.
To measure the performance of our code, the achieved memory bandwidth (BW) is
used as a metric. This allows us, for our primarily memory bound problem, to
easily compare the performance of our implementation with the peak performance
of the hardware we run the simulations on. In addition, the run time can be
easily obtained from this metric as well.

For the 3x3v (i.e. the six-dimensional) problem solved using the second-order
Strang splitting scheme nine advections have to be performed in each time
step. This implies that the entire data has to be accessed 18 times (one read
and one write per advection). For the 2x2v and 2x3v problem, 6 advections have
to be performed and the data is consequently accessed 12 times. Moreover, for
calculating the electric field once more the entire data has to be accessed
(to compute the density $\rho$). This is done twice (once at the end of the
time step to compute diagnostic values such as the electric energy and once at
the middle of the time step). Therefore, we have

$\text{achieved
bandwidth}=\frac{F_{d_{x},d_{v}}\cdot\text{sizeof}(\texttt{fp})\cdot\text{dof}\cdot
10^{-9}}{\text{time per step in seconds}}\text{GB/s}$

where fp is single or double and $F_{d_{x},d_{v}}$ is $20$ for the 3x3v
problem and 14 otherwise. For comparison, we give in Table 1 the theoretical
peak bandwidth for the different GPUs that we will use in this work. These are
the values specified by the vendors and the actual bandwidth obtained, even in
a simple benchmark such as copy, is somewhat less.

GPU | memory bandwidth | memory size
---|---|---
GTX1080 | 320 GB/s | 8 GB
Titan V | 653 GB/s | 12 GB
V100 | 900 GB/s | 16 GB
A100 | 1 555 GB/s | 40 GB
Table 1: The different GPUs used in this work and their specifications.

###  4.1 Single GPU performance

First, we compare the performance of the FFT and discontinuous Galerkin
Poisson solver. The results are presented in table 2. The run time for the FFT
and the discontinuous Galerkin method is very similar, particular for the 2x3v
and 3x3v problem. For the 2x2v problem the FFT based scheme is faster by
approximately 30%. This can be explained by the fact that the degrees of
freedom per direction are reduced for higher dimensional problems (due to
computational constraints) and thus not as many iterations are needed to
obtain convergence. Since the FFT Poisson solver is faster and there is no
disadvantage of using it for single node computations (we will discuss its
role in scaling to larger systems later in this work), all of the computations
performed in this section use the FFT based Poisson solver.

¸  FFT CGdG dof time/step time/step Average iter. 2x2v $220^{4}$ 0.372s 0.477s
173 2x3v $72^{5}$ 0.299s 0.306s 55 3x3v $36^{6}$ 0.619s 0.646s 38

Table 2: Performance using different Poisson solvers on our local GPU cluster.
The computations are done in double precision on an A100 GPU.

In table 3 and 4 we show the achieved bandwidth in the 2x3v and the 3x3v case
using four different types of GPUs. We achieve approximately half the
bandwidth of what is theoretically possible on each platform. We consider this
an excellent result as the memory access pattern is significantly less uniform
than in a stencil code and there is a significant amount of computation that
is not taken into account by our estimate of the achieved bandwidth (e.g.
computing the electric field from the density $\rho$; see Einkemmer (2020) for
a more detailed discussion).

| GTX | TitanV | V100 | A100
---|---|---|---|---
dof | single | single | single | single | double
$56^{5}$ | 104 | 356 | 412 | 583 | 649
$60^{5}$ | 100 | 377 | 430 | 607 | 670
$64^{5}$ |  | 401 | 482 | 666 | 717
$68^{5}$ |  |  | 454 | 634 | 696
$72^{5}$ |  |  |  | 680 | 727
$76^{5}$ |  |  |  | 660 |
$80^{5}$ |  |  |  | 708 |
$84^{5}$ |  |  |  | 681 |
Table 3: Achieved bandwidth for the 2x3v case on a single GPU using single and double precision. Missing values imply that there is not enough memory on the specific GPU for the specified resolution. | GTX | TitanV | V100 | A100
---|---|---|---|---
dof | single | single | single | single | double
$28^{3}28^{3}$ | 64 | 281 | 332 | 415 | 503
$28^{3}32^{3}$ | 66 | 309 | 366 | 462 | 555
$32^{3}32^{3}$ |  | 355 | 418 | 534 | 590
$32^{3}36^{3}$ |  |  | 442 | 554 | 626
$36^{3}36^{3}$ |  |  |  | 460 | 563
$36^{3}40^{3}$ |  |  |  | 482 |
$40^{3}40^{3}$ |  |  |  | 539 |
Table 4: Achieved bandwidth for the 3x3v case on a single GPU using single and
double precision. Missing values imply that there is not enough memory on the
specific GPU for the specified resolution.

Since single GPU results for the 2x2v and 2x3v case were already presented in
Einkemmer (2020), we will focus on the six-dimensional 3x3v case in the
following. In figure 9 a comparison of the achieved bandwidth between
different orders of the method in the 3x3v case is shown. The 4th and 5th
order method achieve better performance than the 2nd and 3rd order method,
which is similar to the behavior observed in Einkemmer (2020) for the 2x2v and
2x3v case.

Figure 9: Obtained bandwidth in the 3x3v case as a function of the degrees of
freedom where different spatial orders of the method are used. The simulations
are done in single precision on a A100 GPU.

In figure 10 we show the timings for the different parts that constitute the
algorithm in the 2x3v and the 3x3v case using different GPUs. The chosen
number of degrees of freedom is the maximum amount possible on the specific
GPU. In the 2x3v case it can be seen once more that the choice of the Poisson
solver has no impact on the overall performance. The main difference between
the consumer GPU (the GTX1080) and the enterprise grade GPUs is that the
GTX1080 needs more time for the computation of the charge density $\rho$. This
comes from the fact that the GTX1080 is able to perform significantly less
floating-point operations per second than the other GPUs. Those are important
to compute the diagnostic quantities (such as electric energy) in the
reduction step. Similar results can be observed for the 3x3v case.

Figure 10: Timings for the different parts of the algorithm. A single
precision computation with the largest problem that fits into memory for the
GPU used has been conducted.

To conclude this section, we will briefly compare the performance between the
GPU implementation (which is the focus of the present work) and the CPU
implementation. A comparison for the 3x3v problem between the A100 GPU and a
Intel Xeon Gold 6226R dual socket system (having $2\times 16$ CPUs) is given
in figure 11. It can clearly be observed that the GPU drastically outperforms
the CPU based system. In fact, every part of the code is computed by the GPU
in significantly less time. This, in addition to the four and five dimensional
results presented in Einkemmer (2020), underscores the advantage of running
semi-Lagrangian Vlasov simulations on GPU based systems. We also note that
this is enabled by the semi-Lagrangian discontinuous Galerkin method which, in
contrast to spline or FFT based interpolation, maps extremely well to the the
GPU architecture.

Figure 11: Performance comparison of a CPU based system and a A100 GPU. The
GPU is 24.5 times faster than the CPUs. This comparison is obtained using
double precision.

###  4.2 Multiple GPU performance

In this section, we analyze the performance of the code using 4 A100 GPUs on a
single node. We increase the number of grid points in two of the velocity
dimensions proportionally to the number of GPU used. That is we use weak
scaling, which is the only metric of interest for these high-dimensional
problems. For the computation we use the Juwels Booster supercomputer. The
four A100 GPUs on a single node are connected via NVLink3 to each other. The
corresponding performance results are presented in table 5, where we also
report results for a four A100 node without using NVLink. As we can see, using
NVLink results in approximately a factor of $2$ improvement yielding a total
single node achieved bandwidth of between $2.05$ and $2.55$ TB/s.

dim | dof | without NVLink | with NVLink
---|---|---|---
2x2v | $220^{2}440^{2}$ | 1821 GB/s | 2437 GB/s
2x3v | $72^{2}144^{2}72$ | 1202 GB/s | 2567 GB/s
3x3v | $36^{4}72^{2}$ | 859 GB/s | 2052 GB/s
Table 5: Comparison of the achieved bandwidth using four A100 GPUs on one
computation node with and without NVLink. The computations shown here are
performed in double precision.

##  5 Scaling to multiple nodes

As can be observed in the previous sections, the resulting grid is still
relatively coarse. This is especially true in the five and six-dimensional
case. In order to run the simulations on finer grids, more nodes are required.
In this section, we will report weak scaling results for our code on Juwels
Booster. This is the most important metric for practical applications, where a
finer grid is often required in order to obtain physically relevant results.
The reported measured wall clock time is the average over the computed time
steps. We analyze the total computation time and the time for the pure
computation of the advections (denoted as Advection in the figures), the time
which is required for the data transfer before the advection takes place
(denoted as Communication), and the time for computing the electric field
(denoted as E field), which includes the computation of the charge density,
the time for solving the Poisson problem (including all the required
communication) and the broadcast of the electric field. The time step size is
chosen such that the CFL-number is close to one. Moreover, the largest
possible local grid which fits into GPU memory is chosen.

#### 2x2v

Two test cases are considered in the 4 dimensional case. First, the full
domain is partitioned and distributed over multiple GPUs, and second, only the
velocity domain is refined. The latter is quite common as some problems
require a higher resolution in velocity space, for example, to shift back in
time the (purely numerical) recurrence effect in linear Landau damping
simulations Einkemmer and Ostermann (2014a). In table 6 the considered grids
are listed.

local grid 2x2v
---
dim | $x_{1}$ | $x_{2}$ | $v_{1}$ | $v_{2}$
dof | 220 | 220 | 220 | 220

#GPUs | decomp. in $xv$
---|---
1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 2 | 2
8 | 1 | 2 | 2 | 2
16 | 2 | 2 | 2 | 2
32 | 2 | 2 | 2 | 4
64 | 2 | 2 | 4 | 4
128 | 2 | 4 | 4 | 4
256 | 4 | 4 | 4 | 4
512 | 4 | 4 | 4 | 8
1024 | 4 | 4 | 8 | 8

#GPUs | decomp. in $v$
---|---
1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 2 | 2
16 | 1 | 1 | 4 | 4
64 | 1 | 1 | 8 | 8
144 | 1 | 1 | 12 | 12
256 | 1 | 1 | 16 | 16
484 | 1 | 1 | 22 | 22
1024 | 1 | 1 | 32 | 32

Table 6: For the 2x2v problem we use a local grid of $220^{4}$ degrees of
freedom on each GPU. In the multi-node setting the GPUs are either equally
divided among the different coordinate axis (middle table) or only the
velocity direction is refined (bottom table).

Let us first consider the case when the full domain is distributed over
multiple GPUs. The corresponding results are shown in figure 12. A time step
of 0.015 is chosen in order to keep the CFL number below 1, and the CGdG
method is used to solve the Poisson problem. The largest grid here has
$880^{2}1760^{2}$ degrees of freedom. It can be observed that the time to
perform the advection remains constant. This holds for each configuration in
this chapter. The communication time in this setting is relatively small as at
each time step only 2 GB of data has to be transferred to neighboring MPI
processes. The time for computing the electric field increases from 4 to 16
GPUs since a larger problem has to be solved, i.e., going from a grid size of
$220^{2}$ to $440^{2}$ and thus more iterations are required to solve the
linear system. The system size stays then constant until 64 GPUs are used.
Subsequently, the problem size in physical space increases again. That is, a
grid size in the space direction of $880^{2}$ is used when 256 or more GPUs
are used. This again requires more iterations and the wall clock time
increases again. The total time per step rises from about 0.5 seconds to
around 1.0 seconds. Thus running the simulation on $1024$ GPUs a parallel
efficiency of approximately 50% is achieved.

Figure 12: Weak scaling in the 2x2v dimensional setting using multiple GPUs
for each dimension. A time step size of 0.015 and the CGdG method to solve the
Poisson problem is used.

Since the impact of the Poisson solver on the overall performance, due to the
large number of required iterations in the 2x2v setting is quite high, we also
investigate the use of a FFT based method to solve the Poisson problem. The
results are shown in figure 13. The FFT based solver is faster on a single
node but scaling is a challenge. Thus, the overall parallel efficiency
decreases a little, while overall the run time on 1024 GPUs is still reduced
by approximately 10% compared to the CGdG method.

Figure 13: Weak scaling in the 2x2v dimensional setting using multiple GPUs
for each dimension. A time step case size of 0.015 and a FFT based solver for
the Poisson problem is used.

Next, we analyze the behavior when the degrees of freedom are increased in the
velocity directions only. Here, since in physical space no parallel numerical
method for the computations is required, the FFT based Poisson solver is used
since it is faster. In figure 14 the weak scaling results are given. It can be
observed that the required time for computing the electric field increases
when more than 144 GPUs are used. Here, in contrast to the results shown in
figure 12 and 13, not the time for solving the Poisson problem increases, but
the required time to compute the charge density, i.e., the call to MPI_Reduce.
The reason for this is that now a very fine resolution of $7040^{2}$ grid
points in the velocity directions is used. The total run time increases from
about $0.36$ to $0.6$, thus a parallel efficiency of 60% is achieved on $1024$
GPUs.

Figure 14: Weak scaling in the 2x2v dimensional setting using multiple GPUs to
increase the degrees of freedom in the velocity dimensions only.

The simulations presented previously were done with a relatively small time
step size of $0.015$ in order to keep the CFL number in the spatial directions
below one. However, the semi-Lagrangian discontinuous Galerkin scheme is
stable also for larger time steps. We will thus investigate the parallel
efficiency for a CFL number larger than one. The main additional issue here is
that the amount of boundary data that need to be communicated to neighboring
processes increases linearly with the CFL number. We use a local grid with a
resolution of $216^{4}$ points and in total $256$ GPUs, thus distributing the
domain in each direction over 4 GPUs. The local grid here is slightly smaller
than in the previous simulations since additional memory is required to store
the boundary data. The CFL number increases from 1 to 5 by increasing the time
step size by $0.02$ until we reach a step size of $0.095$ (where the CFL
number becomes 5). The results are reported in figure 15. It can be observed
that the time for communication increases continuously, which is expected
since more boundary cells have to be send to the neighboring processes.
Furthermore, the time for computing the advection remains almost constant, and
the time for computing the electric field increases. This comes from the fact
that when solving the Poisson problem iteratively, the starting point is the
electric potential obtained at the previous time step. Thus, when a larger
time step is used, the initial estimate used in the iterative Poisson solver
is farther away from the solution and thus more iterations are required for
convergence. Nevertheless, we note that the run for a CFL number smaller than
$1$ to a CFL number of $5$ increases by approximately 50%. Thus, running the
numerical simulation at a CFL number of $5$ still results in a speedup of
approximately a factor $3.3$ (despite the increased cost of communication per
time step).

Figure 15: Timings of the different parts of the algorithm using CFL numbers
larger than one. A local grid of $216^{4}$ points is used and the simulation
is run on 256 GPUs. We increase the time step from $0.015$ (CFL=1) up to
$0.095$ (CFL=5).

#### 2x3v

In figure 16 and 17 the scaling results for the 2x3v case are shown. Again, as
in the 2x2v case, two configurations are considered, namely the full grid
refinement and increasing the degrees of freedom in the velocity domain only.
The considered configurations are listed in table 7.

local grid $2x3v$
---
dim | $x_{1}$ | $x_{2}$ | $v_{1}$ | $v_{2}$ | $v_{3}$
dof | 72 | 72 | 72 | 72 | 72

#GPUs | decomp. in $xv$
---|---
1 | 1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 2 | 2 | 1
8 | 1 | 1 | 2 | 2 | 2
16 | 1 | 2 | 2 | 2 | 2
32 | 2 | 2 | 2 | 2 | 2
64 | 2 | 2 | 2 | 2 | 4
128 | 2 | 2 | 2 | 4 | 4
256 | 2 | 2 | 4 | 4 | 4
512 | 2 | 4 | 4 | 4 | 4
1024 | 4 | 4 | 4 | 4 | 4
1536 | 4 | 4 | 4 | 4 | 6

#GPUs | decomp. in $v$
---|---
1 | 1 | 1 | 1 | 1 | 1
8 | 1 | 1 | 2 | 2 | 2
64 | 1 | 1 | 4 | 4 | 4
216 | 1 | 1 | 6 | 6 | 6
512 | 1 | 1 | 8 | 8 | 8
1000 | 1 | 1 | 10 | 10 | 10
1440 | 1 | 1 | 12 | 12 | 10

Table 7: For the 2x3v problem we use a local grid of $72^{5}$ degrees of
freedom on each GPU. In the multi-node setting the GPUs are either equally
divided among the different coordinate axis (middle table) or only the
velocity direction is refined (bottom table).

When the full domain is distributed over multiple GPUs, a time step of 0.05 is
chosen which implies a CFL number below 1. The CGdG method is used to solve
the Poisson problem. The largest problem has a grid of $288^{4}432$ and is run
on 1536 GPUs. In figure 16 it can be observed that the computation of the
electric field has much less impact on the overall performance than in the
2x2v setting. This is mainly because solving the Poisson problem with the CGdG
requires less iterations due the smaller number of degrees of freedom per
direction. Nevertheless, also here, as expected, an increase of required wall
time can be observed if the problem size increases. The time for the
communication which is required to perform the advection has more impact on
the overall performance than in the 2x2v setting. The CGdG method allows us to
limit the impact to performance and is the preferred choice for this five-
dimensional problem. For 1536 GPUs we achieve a parallel efficiency of 43%.

Figure 16: Weak scaling in the 2x3v case using multiple GPUs for each
dimension. A time step size of $0.05$ and the CGdG method to solve the Poisson
problem is used.

In problems where only additional degrees of freedom are added in the velocity
directions, similar results as have been obtained in the 2x2v case can be
observed. In this case the FFT based solver is the preferred choice. The
corresponding weak scaling results are shown in figure 17, where we observe
almost ideal scaling for up to 1440 GPUs, The parallel efficiency in this case
is approximately 67%.

Figure 17: Weak scaling in the 2x3v dimensional setting using multiple GPUs
for velocity dimensions only.

#### 3x3v

In the six dimensional setting the only reasonable choice is to distribute the
problem in each direction (if a distribution is done only in the velocity
direction the number of degrees of freedom in physical space is too small to
be of practical interest). The used configurations are listed in table 8.
Using 1536 GPUs, results in a grid of size $72^{1}108^{1}144^{4}$ in double
precision, which is still relatively coarse, and $80^{1}120^{1}160^{4}$ in
single precision.

local grid 3x3v double
---
dim | $x_{1}$ | $x_{2}$ | $x_{3}$ | $v_{1}$ | $v_{2}$ | $v_{3}$
dof | 36 | 36 | 36 | 36 | 36 | 36
local grid 3x3v single
dim | $x_{1}$ | $x_{2}$ | $x_{3}$ | $v_{1}$ | $v_{2}$ | $v_{3}$
dof | 40 | 40 | 40 | 40 | 40 | 40
#GPUs | decomp. in $xv$
1 | 1 | 1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 1 | 1 | 2 | 2
8 | 1 | 1 | 1 | 2 | 2 | 2
16 | 1 | 1 | 2 | 2 | 2 | 2
32 | 1 | 2 | 2 | 2 | 2 | 2
64 | 2 | 2 | 2 | 2 | 2 | 2
128 | 2 | 2 | 2 | 2 | 2 | 4
256 | 2 | 2 | 2 | 2 | 4 | 4
512 | 2 | 2 | 2 | 4 | 4 | 4
1024 | 2 | 2 | 4 | 4 | 4 | 4
1536 | 2 | 3 | 4 | 4 | 4 | 4
Table 8: Configurations in the 3x3v case using multiple GPUs.

In figure 18 the scaling results for the 3x3v case are shown. Due to the
coarse grid, a time step of 0.1 is chosen which is enough to keep the CFL
number below 1. The CGdG solver is the preferred choice to solve the Poisson
problem. It can be observed that the communication time increases
significantly until 64 GPUs are used and remains almost constant afterwards.
The reason for this behavior is that for a six-dimensional problem 64 GPUs are
required so that in each direction we use more than a single GPU. Due to the
nearest neighbor exchange any further parallelization does, in principle, not
contribute to the overall communication time. The time per step increases from
0.6 seconds using one GPU to 1.9 seconds using 1536 GPUs, thus a parallel
efficiency of approximately 30% is achieved.

Figure 18: Weak scaling in the 3x3v case using multiple GPUs. The spatial
configuration is listed in table 8. Computations are done in double precision.

In figure 19 the results using single precision are shown. Since less memory
is required, approximately 1.9 times more degrees of freedom can be considered
compared to double precision. We observe a similar parallel efficiency of
approximately 33%.

Figure 19: Weak scaling in the $3x3v$ case using multiple GPUs. The
configuration is listed in table 8. Computations are performed in single
precision.

##  6 Conclusion

In this work we have demonstrated that semi-Lagrangian discontinuous Galerkin
schemes can be used to run five and six-dimensional Vlasov simulations on
large scale GPU equipped supercomputers. The local nature of this method gives
excellent performance on a single node using 4 GPUs (between $2$ and $2.5$
TB/s depending on the configuration) and scales well to up to 1500 GPUs on
Juwels Booster with a parallel efficiency of between 30% and 67% (depending on
the configuration). Since GPU based supercomputers will play an important role
in exascale systems and, most likely, in the future of high-performance
computing, the results demonstrated in this work pave the way towards running
large scale Vlasov simulations on such systems.

## Acknowledgment

We acknowledge PRACE for awarding us access to JUWELS at GCS@FZJ, Germany. The
computational results presented have been achieved in part using the Vienna
Scientific Cluster (VSC).

## Funding

The author(s) disclosed receipt of following financial support for the
research, authorship, and/or publication of this article: This project has
received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
847476. The views and opinions expressed herein do not necessarily reflect
those of the European Commission.

## Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the
research, authorship, and/or publication of this article.

## References

  * Bigot et al. (2013) Bigot J, Grandgirard V, Latu G, Passeron C, Rozar F and Thomine O (2013) Scaling Gysela code beyond 32K-cores on BlueGene/Q.  In: _ESAIM: Proceedings_ , volume 43. EDP Sciences, pp. 117–135.
  * Camporeale et al. (2016) Camporeale E, Delzanno G, Bergen B and Moulton J (2016) On the velocity space discretization for the Vlasov–Poisson system: Comparison between implicit Hermite spectral and Particle-in-Cell methods.  _Computer Physics Communications_ 198: 47–58.
  * Casas et al. (2017) Casas F, Crouseilles N, Faou E and Mehrenberger M (2017) High-order Hamiltonian splitting for the Vlasov–Poisson equations.  _Numerische Mathematik_ 135(3): 769–801.
  * Cheng and Knorr (1976) Cheng C and Knorr G (1976) The integration of the Vlasov equation in configuration space.  _Journal of Computational Physics_ 22(3): 330–351.
  * Crouseilles et al. (2015) Crouseilles N, Einkemmer L and Faou E (2015) Hamiltonian splitting for the Vlasov–Maxwell equations.  _Journal of Computational Physics_ 283: 224–240.
  * Crouseilles et al. (2009) Crouseilles N, Latu G and Sonnendrücker E (2009) A parallel Vlasov solver based on local cubic spline interpolation on patches.  _Journal of Computational Physics_ 228(5): 1429–1446.
  * Crouseilles et al. (2011) Crouseilles N, Mehrenberger M and Vecil F (2011) Discontinuous Galerkin semi-Lagrangian method for Vlasov-Poisson.  In: _ESAIM: Proceedings_ , volume 32. EDP Sciences, pp. 211–230.
  * Einkemmer (2017) Einkemmer L (2017) A study on conserving invariants of the Vlasov equation in semi-Lagrangian computer simulations.  _Journal of Plasma Physics_ 83(2).
  * Einkemmer (2019) Einkemmer L (2019) A performance comparison of semi-Lagrangian discontinuous Galerkin and spline based Vlasov solvers in four dimensions.  _Journal of Computational Physics_ 376: 937–951.
  * Einkemmer (2020) Einkemmer L (2020) Semi-Lagrangian Vlasov simulation on GPUs.  _Computer Physics Communications_ 254: 107351.
  * Einkemmer and Joseph (2021) Einkemmer L and Joseph I (2021) A mass, momentum, and energy conservative dynamical low-rank scheme for the Vlasov equation.  _Journal of Computational Physics_ : 110495.
  * Einkemmer and Lubich (2018) Einkemmer L and Lubich C (2018) A low-rank projector-splitting integrator for the Vlasov–Poisson equation.  _SIAM Journal on Scientific Computing_ 40(5): B1330–B1360.
  * Einkemmer and Ostermann (2014a) Einkemmer L and Ostermann A (2014a) A strategy to suppress recurrence in grid-based Vlasov solvers.  _The European Physical Journal D_ 68: 197.
  * Einkemmer and Ostermann (2014b) Einkemmer L and Ostermann A (2014b) Convergence analysis of a discontinuous Galerkin/Strang splitting approximation for the Vlasov–Poisson equations.  _SIAM Journal on Numerical Analysis_ 52(2): 757–778.
  * Einkemmer and Ostermann (2014c) Einkemmer L and Ostermann A (2014c) Convergence analysis of Strang splitting for Vlasov-type equations.  _SIAM Journal on Numerical Analysis_ 52(1): 140–155.
  * Epshteyn and Rivière (2007) Epshteyn Y and Rivière B (2007) Estimation of penalty parameters for symmetric interior penalty Galerkin methods.  _Journal of Computational and Applied Mathematics_ 206(2): 843–872.
  * Filbet et al. (2001) Filbet F, Sonnendrücker E and Bertrand P (2001) Conservative numerical schemes for the Vlasov equation.  _Journal of Computational Physics_ 172(1): 166–187.
  * Filbet and Sonnendrücker (2003) Filbet F and Sonnendrücker E (2003) Comparison of Eulerian Vlasov solvers.  _Computer Physics Communications_ 150(3): 247–266.
  * Guo and Cheng (2016) Guo W and Cheng Y (2016) A sparse grid discontinuous Galerkin method for high-dimensional transport equations and its application to kinetic simulations.  _SIAM Journal on Scientific Computing_ 38(6): A3381–A3409.
  * Hariri et al. (2016) Hariri F, Tran T, Jocksch A, Lanti E, Progsch J, Messmer P, Brunner S, Gheller C and Villard L (2016) A portable platform for accelerated PIC codes and its application to GPUs using OpenACC.  _Computer Physics Communications_ 207: 69–82.
  * Howes et al. (2008) Howes G, Dorland W, Cowley S, Hammett G, Quataert E, Schekochihin A and Tatsuno T (2008) Kinetic simulations of magnetized turbulence in astrophysical plasmas.  _Physical Review Letters_ 100: 065004.
  * Klimas and Farrell (1994) Klimas A and Farrell W (1994) A splitting algorithm for Vlasov simulation with filamentation filtration.  _Journal of Computational Physics_ 110(1): 150–163.
  * Klimas and Viñas (2018) Klimas A and Viñas A (2018) Absence of recurrence in Fourier–Fourier transformed Vlasov–Poisson simulations.  _Journal of Plasma Physics_ 84(4).
  * Kormann et al. (2019) Kormann K, Reuter K and Rampp M (2019) A massively parallel semi-Lagrangian solver for the six-dimensional Vlasov–Poisson equation.  _The International Journal of High Performance Computing Applications_ 33(5): 924–947.
  * Kormann and Sonnendrücker (2016) Kormann K and Sonnendrücker E (2016) Sparse Grids for the Vlasov–Poisson Equation.  In: _Sparse Grids and Applications - Stuttgart 2014_. Springer International Publishing, pp. 163–190.
  * Mehrenberger et al. (2013) Mehrenberger M, Steiner C, Marradi L, Crouseilles N, Sonnendrücker E and Afeyan B (2013) Vlasov on GPU (VOG project).  In: _ESAIM: Proceedings_ , volume 43. pp. 37–58.
  * NVIDIA (2018) NVIDIA (2018) Cub.  https://nvlabs.github.io/cub/.
  * Rivière (2008) Rivière B (2008) _Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation_.  Society for Industrial and Applied Mathematics.  10.1137/1.9780898717440.
  * Rossmanith and Seal (2011) Rossmanith J and Seal D (2011) A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov–Poisson equations.  _Journal of Computational Physics_ 230(16): 6203–6232.
  * Sonnendrücker et al. (1999) Sonnendrücker E, Roche J, Bertrand P and Ghizzo A (1999) The semi-Lagrangian method for the numerical resolution of the Vlasov equation.  _Journal of Computational Physics_ 149(2): 201–220.
  * Verboncoeur (2005) Verboncoeur JP (2005) Particle simulation of plasmas: review and advances.  _Plasma Physics and Controlled Fusion_ 47(5A): A231.
  * Zweibel and Yamada (2009) Zweibel E and Yamada M (2009) Magnetic reconnection in astrophysical and laboratory plasmas.  _Annual Review of Astronomy and Astrophysics_ 47: 291–332.

##  Appendix A Appendix: Discontinuous Galerkin Poisson solver

First, we describe the 1d case which can then be easily extended to arbitrary
dimensions, as we will show. The same polynomial space

$\displaystyle V^{k}_{h}=\\{$ $\displaystyle\varphi\in
L^{2}(\Omega):\varphi\in\mathbb{P}^{k}(I_{i}),$ $\displaystyle
I_{i}=[x_{i},x_{i+1}],\quad i=0,\ldots,N-1\\},$

that is used to perform the advections is chosen. As basis for $V_{h}^{k}$ the
Lagrange polynomials interpolating on the $k+1$ Gauss–Legendre points $x_{il}$
restricted to the cell $I_{i}$ are chosen:

$\varphi_{ij}(x)=\left\\{\begin{aligned} \prod_{l=0,\ldots,k;l\neq
j}\frac{x-x_{il}}{x_{ij}-x_{il}}&\quad x\in I_{i}\\\ 0\qquad\quad&\quad
x\in\Omega\setminus I_{i}\end{aligned}\right.$

The first step in deriving the discontinuous Galerkin Poisson solver consists
in multiplying the Poisson equation (2) by a test function $v\in V_{h}^{k}$,
integrating over each cell $I_{i}=[x_{i},x_{i+1}]$ and summing up all terms.
For more details, see Rivière (2008). This gives, in 1d, the following
equation

$\begin{split}\sum_{i=0}^{N-1}\bigg{(}\int_{I_{i}}&\phi^{\prime}(x)v^{\prime}(x)\,dx-[\phi^{\prime}(x_{i})]\\{v(x_{i})\\}-\\\
&\\{\phi^{\prime}(x_{i})\\}[v(x_{i})]+\frac{\sigma}{h}[v(x_{n})][\phi(x_{n})]\bigg{)}\\\
&=\int_{\Omega}(1-\rho(x))v(x)\,dx,\end{split}$ (4)

where $\\{v(x_{i})\\}=\tfrac{1}{2}(v(x_{i}^{-})+v(x_{i}^{+}))$ is the average
of the function $v$ at the point $x_{i}$ and
$[v(x_{i})]=v(x_{i}^{-})-v(x_{i}^{+})$ is the jump height. From the left hand
side of (4) a symmetric bilinear form $a:V_{h}^{k}\times
V_{h}^{k}\to\mathbb{R}$ can be derived. Therefore, we are looking for
$\phi^{\text{DG}}\in V_{h}^{k}$ such that

$a(\phi^{\text{DG}},v)=l(v)\quad\forall v\in V_{h}^{k}$

where the linear form $l$ is defined by the right-hand-side of (4). Plugging
$\phi^{\text{DG}}(x)=\sum_{i}\sum_{j=0}^{k}\phi(x_{ij})\varphi_{ij}(x)$ and
$v=\varphi_{lm}$ into $a$ results in a linear system of the form $Au=b$ where
we are looking for the unknown $u$ specified by the coefficients
$\phi(x_{ij})$. This method is called the symmetric interior penalty Galerkin
(SIPG) method. It can be easily verified that this matrix $A$ is block
tridiagonal, and that it can be constructed from two small matrices $M$ and
$B$ of size $k+1$,

$A=\frac{1}{h}\begin{pmatrix}M&B&&&&&B^{T}\\\ B^{T}&M&B&&&&\\\ &&&&&&\\\
&&&&&&\\\ &&&&B^{T}&M&B\\\ B&&&&&B^{T}&M\\\ \end{pmatrix}$

To ensure coercivity of the bilinear form $a$ and positive definiteness of the
matrix $A$, the condition $\sigma\geq k^{2}$ is required, see Epshteyn and
Rivière (2007). The matrices $M$ and $B$ are defined as follows

$\displaystyle(M)_{i,j}=$
$\displaystyle+\tfrac{1}{2}\varphi_{j}^{\prime}(0)\varphi_{i}(0)+\tfrac{1}{2}\varphi_{j}(0)\varphi_{i}^{\prime}(0)+\sigma\varphi_{j}(0)\varphi_{i}(0)$
$\displaystyle-\tfrac{1}{2}\varphi_{j}^{\prime}(1)\varphi_{i}(1)-\tfrac{1}{2}\varphi_{j}(1)\varphi_{i}^{\prime}(1)+\sigma\varphi_{j}(0)\varphi_{i}(1)$
$\displaystyle+\tfrac{1}{2}\sum_{l=0}^{k}\varphi_{i}^{\prime}(x_{l})\varphi_{j}^{\prime}(x_{l})w_{l},$
$\displaystyle(B)_{i,j}=$
$\displaystyle-\tfrac{1}{2}\varphi_{j}^{\prime}(0)\varphi_{i}(1)+\tfrac{1}{2}\varphi_{j}(0)\varphi_{i}^{\prime}(1)-\sigma\varphi_{j}(0)\varphi_{i}(1),$

where we denote with $\varphi_{l}(x)$ the Lagrange polynomials in the interval
$[0,1]$ that interpolate the Gauss–Legendre quadrature points $x_{l}$. The
$w_{l}$ are the Gaussian quadrature weights. In the higher dimensional
setting, the bilinear form is given as follows

$\displaystyle a(w,v)=$ $\displaystyle\sum_{E}\int_{E}\nabla v\cdot\nabla
w-\sum_{e}\int_{e}\\{\nabla v\cdot\textbf{n}_{e}\\}[w]-$
$\displaystyle\sum_{e}\int_{e}\\{\nabla w\cdot
n_{e}\\}[v]+\sum_{e}\frac{\sigma}{|E|/|e|}\int_{e}[v][w],$

where $E$ in 2d are rectangles and in 3d are squares, and $e$ are the edges of
the rectangles or the faces of the squares. Usually, the denominator in the
last term is defined as $|e|^{\beta}$, but setting $\beta=1$ in 2d and
$\beta=\tfrac{1}{2}$ in 3d is almost equivalent if the length of the edges do
not vary much. Since in higher dimensions the basis is the tensor product of
the 1d basis, the resulting matrix can be written in the following form in 2d

$(h_{y}W_{y}\otimes A_{x}+h_{x}A_{y}\otimes W_{x})u=b$

and in 3D

$\displaystyle(h_{y}h_{z}W_{z}\otimes W_{y}\otimes
A_{x}+h_{x}h_{z}W_{z}\otimes A_{y}\otimes W_{x}+$ $\displaystyle
h_{x}h_{y}A_{z}\otimes W_{y}\otimes W_{x})u=b,$

where $W_{x}$, $W_{y}$ and $W_{z}$ are diagonal matrices with the Gaussian
quadrature weights on the diagonal.

To solve this linear system at every time step, we are using the conjugate
gradient method. The reason why we prefer using this method over a direct
method is, that as starting vector of this iterative method the electric
potential $\phi$ computed at the previous time step can be used. Therefore,
the number of iterations is small, since, if short time steps are used $\phi$
does not change dramatically. Moreover, the whole matrix is never constructed,
a function is implemented which performs the action of this matrix applied on
a vector. Thus, we use a completely matrix-free implementation.

If multiple compute nodes are used, only the boundary cells have to be
transferred to the neighboring processes at each iteration. Therefore the
method can be parallelized easily. Unfortunately, the matrix is singular since
the system does not contain a constraint on the mean value of $\phi$. In
principle this is not a problem since the electric field is the gradient of
the potential function $\phi$ and therefore this constant shift vanishes.
Nevertheless, it is difficult to construct a computationally cheap termination
criteria in the CG method when the mean is growing continuously. Therefore, we
embedded this condition in the CG method itself with almost no additional
computational cost. The resulting CG method, which is mean preserving, is
described in algorithm 1.

Algorithm 1 Mean preserving conjugate gradient method

$r_{0}=b-Au_{0},\text{mean}(u_{0})=0$

$r_{0}=r_{0}-\text{mean}(r_{0})$, $p_{0}=r_{0}$

$\eta=(b+r_{0})^{T}u_{0}$

while $i\leq\text{MAXITER}$ do

$[w_{i},\text{mean\textsubscript{w}}]=Ap_{i}$

$\alpha_{i}=r_{i}^{T}r_{i}/p_{i}^{T}w_{i}$

$u_{i+1}=u_{i}+\alpha_{i}p_{i}$

$r_{i+1}=r_{i}-\alpha_{i}(w_{i}-\text{mean\textsubscript{w}})$

$\eta=\eta+\alpha_{i}r_{i}^{T}r_{i}$

if $\|r_{i+1}\|_{2}<\text{tol}\|b\|_{2}$ then

break

end if

$\beta_{i}=r_{i+1}^{T}r_{i+1}/r^{T}_{i}r_{i}$

$p_{i+1}=r_{i+1}+\beta_{i}p_{i}$

end while

To ensure that the dG Poisson solver is working properly, we test it first on
a simple example with periodic boundary conditions,

$\displaystyle-\Delta\phi(x)$ $\displaystyle=\cos(0.5x).$

In figure 20 we observe that the error is decreasing with order $k+1$ (note
that the dG Poisson solver converges with the same order as the dG advection
solver).

Figure 20: Relative $L^{\infty}$ error as a function of the number of cells
for the discontinuous Galerkin Poisson solver in 2d and 3d.

In the Vlasov-Poisson simulations, we set the relative tolerance in the CG
method below $10^{-3}$ which is sufficient as can be seen in figures 7 and 8.
Moreover, $\sigma$ is set to its minimum allowed value. We observed, that when
$\sigma$ is increasing, also the condition number of the matrix $A$ is
increasing and more iterations are required to reach convergence.
Additionally, we observe that convergence takes (almost) place within
$(\prod_{l}^{dx}N_{l}(k_{l}+1))^{1/dx}$ iterations. Therefore, we set the
maximum number of allowed iterations to this threshold times a safety factor.