]> https://gitweb.dealii.org/ - release-papers.git/commitdiff
Add sections: ca, pft, mf
authorpeterrum <peterrmuench@aol.com>
Sun, 10 May 2020 14:34:16 +0000 (16:34 +0200)
committerpeterrum <peterrmuench@aol.com>
Sun, 10 May 2020 14:34:16 +0000 (16:34 +0200)
9.2/paper.bib
9.2/paper.tex

index 567e8fcc4b94330f54360342b7278bb10a829d27..f855c13a7757de7b2c0f67687c5ef4432f06e683 100644 (file)
@@ -898,3 +898,22 @@ year = {2009}
  keywords = {GPU, adaptive refinement, finite element methods, hanging nodes, matrix free},
 } 
 
+@article{hoefler2010scalable,
+  title={Scalable communication protocols for dynamic sparse data exchange},
+  author={Hoefler, Torsten and Siebert, Christian and Lumsdaine, Andrew},
+  journal={ACM Sigplan Notices},
+  volume={45},
+  number={5},
+  pages={159--168},
+  year={2010},
+  publisher={ACM New York, NY, USA}
+}
+
+@misc{munch2020hyperdeal,
+    title={hyper.deal: An efficient, matrix-free finite-element library for high-dimensional partial differential equations},
+    author={Peter Munch and Katharina Kormann and Martin Kronbichler},
+    year={2020},
+    eprint={2002.08110},
+    archivePrefix={arXiv},
+    primaryClass={cs.MS}
+}
index 5f42446016cc3bde0c4bd36f51fefee47f3400b3..d21c1c446e40554740b86c8c31f1a17b94f932d2 100644 (file)
@@ -14,6 +14,8 @@
 \usepackage{graphicx}
 \usepackage{xspace}
 
+%\renewcommand{\baselinestretch}{2.0}
+
 \usepackage[normalem]{ulem}
 
 \pgfplotsset{compat=1.9}
@@ -45,6 +47,7 @@
     Timo Heister,
 %     Luca Heltai,
     Martin Kronbichler,
+    Peter Munch,
 %     Ross Maguire Kynch,
 %     Matthias Maier,
 %     Jean-Paul Pelteret,
  \affil[9]{Institute for Computational Mechanics,
    Technical University of Munich,
    Boltzmannstr.~15, 85748 Garching, Germany.
-   {\texttt{kronbichler@lnm.mw.tum.de}}}
+   {\texttt{kronbichler/munch@lnm.mw.tum.de}}}
 %
 % \author[10]{Ross~Maguire~Kynch}
 % \affil[10]{Zienkiewicz Centre for Computational Engineering,
   3368 TAMU,
   College Station, TX 77845, USA.
   {\texttt{maier@math.tamu.edu}}}
+  
+  \author[9,13]{Peter Munch}
+  
+
+ \affil[13]{Institute of Materials Research, Materials Mechanics, 
+ Helmholtz-Zentrum Geesthacht, 
+ Max-Planck-Str. 1, 21502 Geesthacht, Germany.
+   {\texttt{peter.muench@hzg.de}}}  
+  
 %
 % \author[4]{Jean-Paul~Pelteret}
 %
@@ -291,6 +303,271 @@ in the file that lists all changes for this release}, see \cite{changes91}.
 \label{subsec:bla1}
 
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Improved large-scale performance}
+\label{subsec:performance}
+
+Large-scale simulations with 304,128 cores have revealed bottlenecks in release 
+9.1 during initialization due to the usage of expensive collective operations 
+like \texttt{MPI\_Allgather()} and \texttt{MPI\_\allowbreak Alltoall()}. These 
+operations are used to retrieve or even store information about all processes. 
+For example, all processes have stored in  an array the number of degrees of 
+freedom each process owns. This information is in particular needed to set up 
+the  \texttt{Utilities::MPI::Partitioner} class, which contains the 
+point-to-point communication pattern for vector ghost-value updates and 
+compressions. In release 9.2, we have removed such arrays and have replaced 
+the \texttt{MPI\_Allgather}/\allowbreak\texttt{MPI\_\allowbreak Alltoall} 
+function calls by consensus algorithms~\cite{hoefler2010scalable}, which can be 
+found in the namespace \texttt{Utilities::\allowbreak MPI::\allowbreak ConsensusAlgorithms}: now, only the locally relevant information is computed 
+(and recomputed) when needed, using these algorithms. 
+
+Consensus algorithms are algorithms dedicated to efficient dynamic-sparse 
+communication patterns. In this context, the term ``dynamic-sparse'' means 
+that by the time this function is called, the other processes do not know 
+yet that they have to answer requests and
+each process only has to communicate with a small subset of processes of the 
+MPI communicator. We provide two flavors of the consensus algorithm: the two-step 
+approach \texttt{ConsensusAlgorithms::PEX} and the \texttt{ConsensusAlgorithms::NBX}, 
+which uses only point-to-point communications and a single \texttt{MPI\_IBarrier()}. 
+The class \texttt{ConsensusAlgorithms::Selector} selects one of the two previous 
+algorithms, depending on the number of processes.
+
+Due to the excellent scalability of the consensus algorithms, users are encouraged 
+to use them for their own dynamic-sparse problems by providing a list of target 
+processes and pack/unpack routines either by implementing the interface 
+\texttt{ConsensusAlgorithms::Process} or by providing \texttt{std::function} 
+objects to \texttt{ConsensusAlgorithms::AnonymousProcess}.
+
+
+The \texttt{ConsensusAlgorithms} are used by now internally in many places. These 
+places are appropriate starting points for users for their own application of the 
+\texttt{ConsensusAlgorithms} infrastructure.
+For example, to set up the partitioners, the new function 
+\texttt{compute\_index\_owner() } is used: given an index set containing the 
+locally owned indices and an index set containing the ghost indices, it returns 
+the owner of the ghost indices. Consensus algorithms are used now also in the 
+functions \texttt{compute\_point\_to\_point\_communication\_pattern()} and 
+\texttt{compute\_\allowbreak n\_\allowbreak point\_\allowbreak to\_\allowbreak point\_\allowbreak communications()}. 
+Furthermore, it is utilized in the class \texttt{NoncontiguousPartitioner} to 
+efficiently permute distributed solution vectors globally in an arbitrary order, e.g., 
+to interface with external libraries that prescribe a certain partitioning and 
+padding of the data.
+
+By replacing the collective communications during setup and removing the arrays 
+that contain information for each process (enabled by the application of consensus 
+algorithms and other modifications---a full list of modifications leading to this 
+improvement can be found online), we were able to significantly improve the setup 
+time for large-scale simulations and to solve a Poisson problem with multigrid 
+with 12T unknowns. 
+Figure~\ref{} compares the timings of a simulation (incl. setup) with the 
+previous release 9.1 and with the current release 9.2. 
+{\color{red}TODO[Peter/Martin] description of the results}
+
+The new code has been also applied to solve problems with adaptively refined 
+meshes with more than 4B unknowns. {\color{red}TODO[Timo]}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{A new fully distributed triangulation class}
+\label{subsec:pft}
+
+By release 9.1, \texttt{deal.II} had three types of triangulation classes: the 
+serial triangulation class \texttt{Triangulation} as well as the parallel 
+triangulation classes \texttt{parallel::shared::Triangulation} and \texttt{parallel::distributed::Triangulation}. The latter builds around a serial 
+triangulation and uses \texttt{p4est} as an oracle during adaptive mesh refinement. 
+All these triangulation classes have in common that the coarse grid is shared by 
+all processes and the actual mesh used for computations is constructed by repeated 
+local and/or global refinement, which adapts nicely to curved boundaries described 
+by the \texttt{Manifold} class. However, this way to construct a computational 
+mesh has its limitations in industrial applications where, often, the mesh comes 
+from an external CAD program in the form of a file that already contains millions 
+or tens of millions of cells with a similar number of vertices. In such a case, 
+refining a mesh is not practical, since it would increase the computational effort 
+and new vertices would not be placed on curved boundaries. A problem that arises 
+for such large grids in the context how meshes have been treated in \texttt{deal.II} 
+until now is that the coarse grid, i.e., potentially the whole mesh is shared by 
+all processes. It might be a major difficulty in MPI-only parallelized applications 
+on modern multi-core processors, since these applications might already run out of 
+main memory during reading the mesh. Not even increasing the number of processes 
+might help in this situation. 
+
+The new class \texttt{parallel::fullydistributed::Triangulation} targets this issue 
+by distributing also the coarse grid, which is the reason for the name of the 
+chosen namespace: it distributes the coarse grid as well as the refinement levels. Such 
+a triangulation can be created by providing a \texttt{TriangulationDescription::Description} struct to each process, containing 
+1) the relevant data to construct the local part of the coarse grid, 2) the 
+translation of the local coarse-cell IDs to globally unique IDs, 3) the hierarchy 
+of mesh refinements, and 4) the owner of the cells on the active mesh level as well 
+as on the multigrid levels. Once the triangulation is set up with this struct, no 
+changes to the mesh are allowed at the moment.
+The \texttt{TriangulationDescription::Description} struct can be filled manually or 
+by the utility functions from the \texttt{TriangulationDescription::Utilities} 
+namespace. The function \texttt{create\_\allowbreak description\_\allowbreak 
+from\_ \allowbreak triangulation()} can convert a base triangulation (partitioned 
+serial \texttt{Tri\-angulation} and \texttt{parallel::distributed::Triangulation}) 
+to such a struct. The advantage of this approach is that all known utility 
+functions from the namespaces \texttt{GridIn} and \texttt{GridTools} can be used 
+on the base triangulations before converting them to the structs. Since this 
+function suffers from the same main memory problems as described above, we also 
+provide the function \texttt{create\_description\_from\_triangulation\_in\_groups()}, 
+which creates the structs only on the master process in a process group. These 
+structs are filled one by one and are sent to the relevant processes once they are 
+ready. A sensible process group size might contain all processes of one compute node.
+
+
+The new (fully) distributed triangulation class works---in contrast to  
+\texttt{parallel::distributed::\allowbreak Tri\-an\-gu\-la\-tion}---not only for 2D- and 3D- but also for 
+1D-problems. It can be used in the context of geometric multigrid methods and 
+supports periodic boundary conditions. Furthermore, hanging nodes are supported.
+
+We intend to extend the usability of the new triangulation class in regard of 
+different aspects, e.g., I/O. In addition, we would like to enable adaptive mesh 
+refinement, a feature of the other triangulation classes, which is very much 
+appreciated by many users. For repartitioning, we plan to use an oracle approach 
+known from \texttt{parallel::distributed::Triangulation}. Here, we would like to 
+rely on a user-provided partitioner, which might also be a graph partitioner.
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Advances of the SIMD capabilities and the matrix-free infrastructure}
+\label{subsec:mf}
+
+%\begin{itemize}
+%\item ECL
+%\item VectorizedArrayType
+%\end{itemize}
+
+The class \texttt{VectorizedArray<Number>} is a key ingredient for the high 
+node-level performance of the matrix-free algorithms in deal.II. It is a wrapper 
+class around $n$ vector entries of type \texttt{Number} and delegates relevant 
+function calls to appropriate Intrinsics instructions. Up to release 9.1, the 
+vector length $n$ has been set at compile time of the library to the highest 
+possible value supported by the given processor architecture.
+
+The class \texttt{VectorizedArray} has been made more user-friendly by making 
+it compatible with the STL algorithms found in the header \texttt{<algorithm>}. 
+Now, it has following features:
+\begin{itemize}
+\item \texttt{VectorizedArray::size()} returns the vector length. This function 
+replaces the public static attribute \texttt{VectorizedArray::n\_array\_elements}, 
+which has been deprecated.
+\item \texttt{VectorizedArray::value\_type} contains the underlying number type of 
+the array.
+\item \texttt{VectorizedArray} has an output operator 
+\texttt{std::ostream\& operator<<(\&out, \&p)}.
+\item \texttt{VectorizedArray::begin()} and \texttt{VectorizedArray::end()} allow 
+range-based iteration over all vector entries.
+\end{itemize}
+Furthermore, the \texttt{VectorizedArray} class supports the following (tested) 
+algorithms: \texttt{std::\allowbreak ad\-vance()}, \texttt{std::distance()}, and \texttt{std::max\_element()}.
+
+It has been also extended with the second optional template argument 
+\texttt{VectorizedArray<Number, size>} with \texttt{size} being related to the 
+vector length, i.e., the number of lanes to be used and the instruction set to be 
+used. By default, the number is set to the highest value supported by the given 
+hardware, i.e., \texttt{VectorizedArray<double>} is translated on Skylake-based 
+processors to \texttt{VectorizedArray<double, 8>}. A full list of supported 
+vector lengths are presented in Table~\ref{tab:vectorizedarray}.
+
+All matrix-free related classes (like \texttt{MatrixFree} and \texttt{FEEvaluation}) 
+have been templated with the floating-point number type \texttt{Number} (e.g. \texttt{double} or \texttt{float}); the computations were performed implicitly 
+on \texttt{VectorizedArray<Number>} structs with the highest 
+instruction-set-architecture extension, with each lane responsible for a separate 
+cell (vectoriziation over elements). In release 9.2, all matrix-free classes 
+have been extended with a new optional template argument specifying the 
+\texttt{VectorizedArrayType}. This allows users to select the vector length/ISA and, 
+as a consequence, the number of cells to be processed at once. 
+\begin{table}
+\caption{Supported vector lengths of the class \texttt{VectorizedArray} and 
+the corresponding instruction-set-architecture extensions. }\label{tab:vectorizedarray}
+\centering
+\begin{tabular}{ccc}
+\toprule
+\textbf{double} & \textbf{float} & \textbf{ISA}\\
+\midrule
+VectorizedArray<double, 1> & VectorizedArray<float, 1> & (auto-vectorization) \\
+VectorizedArray<double, 2> & VectorizedArray<float, 4> & SSE2 \\ 
+VectorizedArray<double, 4> & VectorizedArray<float, 8> & AVX/AVX2 \\ 
+VectorizedArray<double, 8> & VectorizedArray<float, 16> & AVX-512 \\ 
+\bottomrule
+\end{tabular}
+
+\caption{Comparison of relevant SIMD-related classes in deal.II and C++20.}\label{tab:simd}
+\centering
+\begin{tabular}{cc}
+\toprule
+\textbf{VectorizedArray (deal.II)} & \textbf{std::simd (C++20)} \\
+\midrule
+VectorizedArray<Number> & std::experimental::native\_simd<Number> \\
+VectorizedArray<Number, size> & std::experimental::fixed\_size\_simd<Number, size> \\ \bottomrule
+\end{tabular}
+\end{table}
+
+In standard (2D/3D) matrix-free applications with moderate polynomial degrees, 
+we found that there is no reason to modify the default vector length of 
+\texttt{VectorizedArray}, since it reaches the highest possible computational 
+throughput despite of increased memory footprint. However, in the deal.II-based 
+library \texttt{hyper.deal}~\cite{munch2020hyperdeal}, where the same matrix-free 
+infrastructure was used for solving the 6D Vlasov--Poisson equation with high-order 
+discontinuous Galerkin methods (with more than 1024 degrees of freedom per cell), the 
+benefit of decreasing the number of cells processed by a single SIMD instruction was 
+shown. That library works with \texttt{MatrixFree} objects of different SIMD-vector 
+length in the same application, which would not have been possible before this 
+release.
+
+In a next step, we intend to support that users could set 
+\texttt{VectorizedArrayType} to \texttt{Number}, which would lead to the usage not 
+of \texttt{VectorizedArray<Number, 1>} but of a specialized code-path exploiting 
+vectorization within an element~\cite{KronbichlerKormann2019}.
+
+A side effect of introducing the new template argument \texttt{VectorizedArrayType} 
+in the \texttt{MatrixFree} classes is that any data structures 
+\texttt{VectorizedArrayType} can be processed if they support required 
+functionalities like \texttt{size()} or \texttt{value\_type}. In this context, we 
+would like to highlight that the new \texttt{C++20} feature \texttt{std::simd} 
+can be processed by the matrix-free infrastructure with minor internal 
+adjustment as an open pull request shows 
+(see \url{https://github.com/dealii/dealii/pull/9994}).  
+Table~\ref{tab:simd} gives a comparison of the deal.II-specific SIMD classes and 
+the equivalent C++20 classes. We welcome the standardization of the SIMD 
+parallelization paradigm in C++ and intend to replace step by step our own 
+wrapper class, which has been continuously developed over the last decade. We 
+would like to emphasize that the work invested in this class was not in vain, 
+since many performance-relevant utility functions implemented with \texttt{VectorizedArray} in mind (e.g., \texttt{vectorized\_load\_and\_transpose} 
+and \texttt{vectorized\_transpose\_and\_store}) will be still used, since they 
+have not become part of the standard.
+
+Further additions to the \texttt{MatrixFree} infrastructure consist of:
+\begin{itemize}
+\item a new variant of \texttt{MatrixFree::cell\_loop()}: It takes two
+\texttt{std::function} objects with ranges on the locally owned degrees of freedom, one
+with work to be scheduled before the cell operation first touches some
+unknowns and another with work to be executed after the cell operation last 
+touches them. The goal of
+these functions is to bring vector operations close to the time when the
+vector entries are accessed by the cell operation, which increases the cache
+hit rate of modern processors by improved temporal locality. 
+\item a new form of loop \texttt{MatrixFree::loop\_cell\_centric()}: This 
+kind of loop can be used in the context of discontinuous Galerkin methods, 
+where both cell and face integrals have to be evaluated. While in the case of 
+the traditional \texttt{loop}, cell and face integrals have been performed 
+independently, the new loop performs all cell and face integrals of a cell in 
+one go. This includes that each face integral has to be evaluated twice, but 
+entries have to be written into the solution vector only once with improved 
+data locality. Previous publications based on \texttt{deal.II} have shown the 
+relevance of the latter aspect for reaching higher performance. 
+\end{itemize}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Advances in GPU support}
+\label{subsec:gpu}
+
+\begin{itemize}
+\item overlapping of computation and communication in the case of CUDA-aware MPI
+\end{itemize}
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsection{New and improved tutorial and code gallery programs}
@@ -306,8 +583,19 @@ In addition, there are seven new tutorial programs:
 \item \texttt{step-47}
 \item \texttt{step-50}
 \item \texttt{step-58}
-\item \texttt{step-65}
-\item \texttt{step-67}
+\item \texttt{step-65} presents \texttt{TransfiniteInterpolationManifold}, a
+manifold class that can propagate curved boundary information into the
+interior of a computational domain, and \texttt{MappingQCache} for fast operations for
+expensive manifolds.
+\item \texttt{step-67} presents an explicit time integrator for the
+compressible Euler equations discretized with a high-order discontinuous
+Galerkin scheme using the matrix-free infrastructure. Besides the use of
+matrix-free evaluators for systems of equations and over-integration, it also
+presents \texttt{MatrixFreeOperators::CellwiseInverseMassMatrix}, a fast implementation
+of the action of the inverse mass matrix in the DG setting using tensor
+products. Furthermore, this tutorial demonstrates i.a. the usage of the new 
+pre and post operations which can be passed to \texttt{cell\_loop()} 
+(see also Subsection~\ref{subsec:mf}) and discusses performance-related aspects.
 \item \texttt{step-69}
 \item \texttt{step-70}
 \end{itemize}
@@ -324,7 +612,13 @@ deprecated classes and functions; and the majority change internal
 interfaces that are not usually used in external
 applications. However, some are worth mentioning:
 \begin{itemize}
-\item
+\item The functions:
+\begin{itemize}
+\item \texttt{DoFHandler::loccaly\_owned\_dofs\_per\_processor()}
+\item \texttt{DoFHandler::loccaly\_owned\_mg\_dofs\_per\_processor()}
+\end{itemize} have been deprecated. As discussed in Subsection~\ref{subsec:performance}, deal.II does not store information for all processes on all processes processes, but only the local information or the locally-relevant information. Users are asked to construct the global information on their own, e.g. by calling \texttt{Utilities::MPI::Allgather(locally\_owned\_info(), comm)}.
+\item 
+
 % \item The \texttt{VectorView} class was removed. We recommend either copying the
 %       vector subset into a \texttt{Vector} or using a \texttt{BlockVector}.
 % \item The function \texttt{Subscriptor::subscribe()}, used through the

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.