In this note, we extend the fast tensor-product algorithm for the simulation of time-dependent partial differential equations. We use the natural tensorization of the space-time domain to propose, after discretization, a set of independent problems, each one with the complexity of a single steady problem. This allows for an efficient parallel implementation that is already interesting on small architectures, but that can also be combined with standard domain-decomposition-based algorithms providing a further direction of parallelism on large computer platforms. Preliminary numerical simulations are presented for a one-dimensional unsteady heat equation

Maday, Yvon

Rønquist, Einar M

Yvon Maday

Einar M. Rønquist

Comptes Rendus Mathématique

C. R. Acad. Sci. Paris, Ser. I 346 (2008) 113–118http://france.elsevier.com/direct/CRASS1/Numerical AnalysisParallelization in time through tensor-product space–time solversYvon Maday a,b, Einar M. Rønquist ca Université Pierre et Marie Curie-Paris 6, UMR 7598, Laboratoire J.-L. Lions, 75252 Paris cedex 05, Franceb Division of Applied Mathematics, Brown University, 182 George Street, Providence, RI 02912, USAc Department of Mathematical Sciences, Norwegian University of Science and Technology, N-7491 Trondheim, NorwayReceived 20 August 2007; accepted 12 September 2007Available online 3 December 2007Presented by Olivier PironneauAbstractIn this Note, we extend the fast tensor-product algorithm for the simulation of time-dependent partial differential equations. Weuse the natural tensorization of the space–time domain to propose, after discretization, a set of independent problems, each one withthe complexity of a single steady problem. This allows for an efficient parallel implementation that is already interesting on smallarchitectures, but that can also be combined with standard domain-decomposition-based algorithms providing a further directionof parallelism on large computer platforms. Preliminary numerical simulations are presented for a one-dimensional unsteady heatequation. To cite this article: Y. Maday, E.M. Rønquist, C. R. Acad. Sci. Paris, Ser. I 346 (2008).© 2007 Published by Elsevier Masson SAS on behalf of Académie des sciences.RésuméAlgorithme de parallélisation en temps basé sur des techniques de produit tensoriel espace–temps. Dans cette Note ongénéralise à la simulation de phénomènes instationnaires l’algorithme de produit tensoriel. On utilise la tensorisation naturelle dudomaine espace–temps pour proposer, après discrétisation un ensemble de problèmes indépendants, chacun d’eux ayant la com-plexité d’un simple problème stationnaire. Ceci permet une mise en œuvre parallèle déjà interessante sur des petites architecturesmais elle peut être également combinée avec des techniques classiques de décomposition de domaine pour utiliser au mieux desarchitectures avec un nombre de processeurs plus important. Des premiers résultats sont présentés sur un problème de la chaleurinstationnaire monodimensionnel. Pour citer cet article : Y. Maday, E.M. Rønquist, C. R. Acad. Sci. Paris, Ser. I 346 (2008).© 2007 Published by Elsevier Masson SAS on behalf of Académie des sciences.Version française abrégéeLa résolution numérique de problèmes modélisés par des équations aux dérivées partielles est maintenant assez bienétablie pour des phénomènes stationnaires. Les techniques itératives comme le “multishooting” de [5], les approchesmultigrilles de [2], et plus récemment les approches pararéelles en temps de [3] permettent de traiter des phénomènesinstationnaires de façon efficace en tenant compte des architectures des calculateurs parallèles. La méthode non ité-rative, que nous présentons ici est une généralisation des techniques de tensorisation présentées par exemple dans [1]E-mail addresses: maday@ann.jussieu.fr (Y. Maday), ronquist@math.ntnu.no (E.M. Rønquist).1631-073X/$ – see front matter © 2007 Published by Elsevier Masson SAS on behalf of Académie des sciences.doi:10.1016/j.crma.2007.09.012114 Y. Maday, E.M. Rønquist / C. R. Acad. Sci. Paris, Ser. I 346 (2008) 113–118pour des discrétisations en espace seulement ; elle peut être combinée à des techniques classiques de décompositionde domaine en espace pour une plus grande efficacité.Considérons par exemple l’équation de la chaleur, Éqs. (1), (2). Après discrétisation en espace par une méthodeclassique, par exemple de différences finies de pas x = L/(N + 1) et une discrétisation en temps implicite, onobtient le système discret (3) où le vecteur um ∈ RN (resp. f m ∈ RN ) représente l’approximation de la solution (resp.la donnée de flux de chaleur) aux nœuds internes et au temps tm, et tm = tm − tm−1 est le pas de temps au niveau tm.La matrice A ∈ RN×N est classiquement symétrique définie positive.À chaque pas de temps on doit donc inverser l’opérateur de Helmholtz discret A + 1tmI et la détermination deum−1 est nécessaire pour calculer um. Pour briser la nature intrinsèquement séquentielle de cette résolution, on utilisel’algorithme de produit tensoriel rapide. Soit A ∈ RN×N et B ∈ RM×M ; on note tout d’abord C = B ⊗ A le produittensoriel de B et A : c’est la matrice C ∈ R(MN)×(MN) de coefficients cij = bklamn, 1  k, l  M , 1  m,n  N ,i = m + (k − 1)N et j = n + (l − 1)N . Avec ces notations, le système discret (3) s’écrit (6) où B est défini en (4), I test la matrice identité de dimension M et I x est la matrice identité de dimension N , de plus u ∈ RMN (resp. f ∈ RMN )représente le vecteur concaténé de la solution et du membre de droite selon (5) à tous les pas de temps. Si l’on supposetous les tm, m = 1, . . . ,M différents, la matrice B est diagonalisable, selon (7). Le système (6) peut donc prendrela forme (11) où l’inconnue w est définie en (10) et le membre de droite g est défini en (13). On remarque aisémentque ce système se déploie en M problèmes indépendants (12) qui peuvent donc être résolus en même temps sur desprocesseurs parallèles, et la solution {um,m = 1, . . . ,M} s’obtient par la transformation inverse (14).En notant par WH la complexité de résolution d’un système de type (3), la complexité de résolution de l’approchetensorisée (non parallélisée) est donnée en (15), en effet la diagonalisation de la matrice B compte pour O(M2) etle calcul des membres de droite (13) ainsi que la transformation inverse (14) compte pour O(M2N). Si l’on disposepar contre de P = M processeurs parallèles, la complexité par processeur devient (16), en ajoutant une complexitéde communication entre les processeurs Wcomm on obtient une efficacité de parallélisme S donnée par (17). Ainsi sil’on suppose d’une part que Wcomm  WH et en outre M  N et MN  WH on conclut à une efficité quasi optimaleS ≈ M .La tensorisation exposée ci dessus peut être effectuée pour des schémas en temps plus précis que le schémasimplicite d’ordre 1 de Euler, tant que la matrice B reste triangulaire inférieure. La simulation numérique suivantebasée sur un schémas implicite rétrograde d’ordre 3 (avec un choix de pas de temps tk = ρk−1t1, k > 1, avecρ = 1,2 se sorte que t15 ≈ 13t1, et T = ∑M=15k=1 tk = 0,592). L’erreur obtenue après un temps T égal à T =100T = 59,2 est 1,4 × 10−3 alors qu’elle est de 4,3 × 10−2 avec les schémas d’ordre 1.1. IntroductionThe simulation of three-dimensional steady problems governed by partial differential equations is nowadays quitecommon; the size of current computers allows one to get very good approximations, close to experimental data,for a large class of problems. The availability of parallel computers has resulted in many contributions related tothe parallelization of problems with respect to the spatial directions. A key ingredient has been to divide the globalproblem through domain decomposition, thus allowing smaller problems to be assigned to different processors.The next challenge, which has already started to be faced, is to solve time-dependent problems with the sameefficiency and ease. Roughly speaking, this amounts to adding an additional direction to the three spatial ones. Mostdiscretizations of time-dependent problems rely on using finite difference methods in time. This results in the need tosolve a sequence of stationary-like problems, for which domain decomposition is primarily the only way to exploitparallel computers. Indeed, the sequence of these stationary-like problems is recursive and is not obviously broken intoindependent subproblems. Among the emerging new methods, let us quote the multishooting technique introduced byNievergelt [5], and the multigrid approach introduced by Hackbusch [2]. More recently a new approach, named theparareal in time method has been introduced by Lions, Maday and Turinici [3] with various contributions showingthe interest of this predictor-corrector methodology.The approach we present here is of a completely different nature. First of all it is not iterative, at least in the versionwe propose here to tackle linear problems. Second, it may be of interest already with very few processors. Third, asa consequence of the diversity of available solvers, it can be combined with the parareal approach as we shall explainin the full paper associated with this Note. The method used here is purely algebraic. It is inspired by the fast tensorY. Maday, E.M. Rønquist / C. R. Acad. Sci. Paris, Ser. I 346 (2008) 113–118 115product solver that we introduced in the recent paper [1] for spatial discretizations, and has also some connection withthe classical way that leads to the harmonic version of a time dependent problem, and that looks for solutions of theform eiωtu(x). Finally, let us indicate that it can easily be combined with standard domain decomposition methods inspace, thus offering enhanced overall parallelism.2. The heat equationConsider the one-dimensional heat equation (with thermal diffusivity κ assumed to be constant)∂u∂t= κ ∂2u∂x2+ f (x, t) in ωx = (0,L), (1)subject to the homogeneous boundary conditions (resp. to the initial condition)u(x = 0, t) = u(x = L, t) = 0 ∀t ∈ ωt = [0, T ], (resp. u(x, t = 0) = u0(x)). (2)As an introductory example to the approach we will present, we consider the numerical solution of this problemusing finite difference methods in space, and an implicit method like the first order backward Euler scheme. The setof fully discrete equations can be expressed asum − um−1tm+ Aum = f m, m = 1, . . . ,M, (3)where the vector um ∈ RN denotes the numerical approximation of the solution at the N internal grid points at time tmand tm is the time step at time level tm; more precisely, tm = tm − tm−1. The vector f m ∈ RN represents the givendata at the internal grid points at time tm, and the matrix A ∈ RN×N , that represents the discrete diffusion operator, issymmetric positive definite.This simple approach to obtain the numerical solution results in the need to solve a system of algebraic equations ateach time step, with the discrete Helmholtz operator A + 1tmI as the system matrix. Note that this is indeed a purelysequential approach; we need to know um−1 before we can compute um by solving a discrete Helmholtz problem.3. Tensor-product representationWe will now present a reformulation of the model problem introduced in the previous section. This reformulationis based on a tensor product representation of time and space and will allow us to obtain a scheme with parallelprocessing possibilities in the time domain.First, let us briefly review the tensor product notation that we will be using. Let A ∈ RN×N and B ∈ RM×M . Thenthe tensor product of B and A, C = B ⊗ A , gives a matrix C ∈ R(MN)×(MN) with matrix entrees cij = bklamn, with1  k, l  M , 1  m,n  N , i = m + (k − 1)N and j = n + (l − 1)N .Associated with the backward Euler scheme (3), we define the matrix B to beB =⎛⎜⎜⎜⎜⎜⎜⎝1t1− 1t21t2− 1t31t3. . .− 1tM1tM⎞⎟⎟⎟⎟⎟⎟⎠. (4)We also define I t to be the identity matrix of dimension M (associated with the time domain) and I x to be theidentity matrix of dimension N (associated with the spatial domain). Finally, we define the concatenated solutionvector u ∈ RMN and the concatenated right-hand side f ∈ RMN asu =⎛⎜⎜⎜⎜⎝u1u2u3...uM⎞⎟⎟⎟⎟⎠ , f =⎛⎜⎜⎜⎜⎜⎝f 1 + 1t1u0f 2f 3...M⎞⎟⎟⎟⎟⎟⎠. (5)f116 Y. Maday, E.M. Rønquist / C. R. Acad. Sci. Paris, Ser. I 346 (2008) 113–118With this notation we can write (3) in tensor product notation as(B ⊗ I x + I t ⊗ A)u = ( I t ⊗ I x)f . (6)Note that this system of algebraic equations represents the numerical solution at M different time levels. Although thesystem matrix is sparse, it does not a priori seem like a good strategy to reduce the overall computational cost.In order to see the potential benefit of this approach, we proceed as follows. We first assume that it is possible todiagonalize the matrix B , i.e., we assume that we can write B asB = SΛS−1. (7)Here, the columns of S are the eigenvectors of B , while Λ is a diagonal matrix comprising the corresponding eigen-values of B . For our particular Euler backward scheme, the matrix B is lower triangular, and the eigenvalues are equalto the diagonal entries of B . Hence, the matrix B is diagonalizable only if all the time steps tm, m = 1, . . . ,M , aredifferent. Assuming that this is the case, we can write (6) as( SΛS−1 ⊗ I x + I t ⊗ A)u = ( I t ⊗ I x)f , (8)or( S ⊗ I x)(Λ ⊗ I x + I t ⊗ A)(S−1 ⊗ I x)u = ( I t ⊗ I x)f . (9)Introducingw = ( S−1 ⊗ I x)u (10)and premultiplying (9) by ( S−1 ⊗ I x) gives(Λ ⊗ I x + I t ⊗ A)w = ( S−1 ⊗ I x)f . (11)An equivalent formulation of (11) is thus(A + λmIx)wm = gm, m = 1, . . . ,M, (12)with a concatenated right-hand sideg = ( S−1 ⊗ I x)f . (13)We have here reformulated the solution in terms of M computational steps. However, in contrast to the systems (3)which represent M sequential steps, the systems (12) represent M independent steps. Hence, (12) can be solved inparallel in a completely decoupled fashion without any communication.Once the systems (12) have been solved for wm, m = 1, . . . ,M , we compute um, m = 1, . . . ,M , from (10), i.e.,u = ( S ⊗ I x)w. (14)4. Computational complexityWe now discuss the computational complexity associated with the approach presented in the previous section. Werecall that the model problem is a simple unsteady heat transfer problem discretized with a standard finite differencemethod in space and a backward Euler scheme in time. The extension to use finite element methods in space isstraightforward.The computational approach assumes that it is possible to diagonalize the matrix B defined in (4). The cost of thediagonalization (7) is O(M2) where M is the number of (standard) time steps following a backward Euler approach.The corresponding memory requirement is O(M2).Next, the cost of computing the transformed right-hand sides (13) is O(M2N), while the cost of solving the Msystems in (12) is the same as the standard approach for solving the M systems (3); the particular cost here dependson the chosen solver for these Helmholtz-type systems. Let us denote the cost per Helmholtz system as WH . The costof performing the final transformation (14) is O(M2N).The total (serial) cost W1 of the new approach is thereforeW1 = c1 · M2 + c2 · M2N + M · WH . (15)Y. Maday, E.M. Rønquist / C. R. Acad. Sci. Paris, Ser. I 346 (2008) 113–118 117The cost W̃1 of following a standard approach is just the last term on the right-hand side.However, assuming now the we have available P = M processors, the computational complexity per processor isWP = c1 · M2 + c2 · MN + WH . (16)In addition, we have a communication cost, Wcomm, associated with distributing the concatenated vector g properlyto all the processors, and collecting the concatenated vector w on each processor. These operations correspond totransposing data sets (all-to-all communication), and are needed in order to compute gm and um on each processor Pm,m = 1, . . . ,M ; see (12) and (14). The speedup S is thusS = W̃1WP + Wcomm =M · WHc1 · M2 + c2 · MN + WH + Wcomm . (17)If Wcomm  WH , M  N , and the work WH satisfies MN  WH , we see that S ≈ M , i.e., we get close to perfectspeedup. The assumption Wcomm  WH can be expected to be satisfied when considering the extension of this ap-proach to solve more realistic multi-dimensional problems. In this case, the assumption M  N can also easily befulfilled. The assumption MN  WH is also valid in many cases, in particular, if M is not too large, and in the contextof solving large multi-dimensional problems using unstructured discretizations. For more specialized situations wherefast tensor-product solvers can be used (e.g., FFT-based algorithms), the computational complexity per time step canbe expected to be a small multiple of N . In such cases, the assumption MN  WH may not be valid (at least forlarger values of M), and the solution approach proposed here may not give good speedup. For more realistic and moregeneral problems, the assumption MN  WH should be analyzed in the context of using iterative solvers for theindividual Helmholtz problems. Even if close of optimal iterative solvers are used, this assumption can be expected tobe valid if M is not too large.5. Generalization and numerical resultsThe method we have just proposed is first order in time, and since it requires that all the time steps are different,the accuracy will be related to the largest time step. In order to make the method more efficient, we propose to use ahigher order scheme in time with time steps tk = ρk−1t1, k > 1, with ρ larger but close to 1, e.g. ρ = 1.2 so thatt15 ≈ 13t1, and T = ∑M=15k=1 tk = 0.592. Note that, as can be expected, choosing ρ much closer to 1 may leadto instabilities due to numerical errors. The scheme we use here comprises a second order Crank–Nicolson schemefor the first time step, followed by a Backward Difference scheme of order 2 to get u2, then a Backward Differencescheme of order 3 to get uk for 2 < k  M = 15. The scheme is a little bit more involved but results in a matrixsimilar to B defined in (4) since B is still lower diagonal, with few terms (actually 3) under the diagonal, and the samediagonal terms as in (4). The complexity is identical as for the first order version we analyzed. We have performedthe simulation for a right-hand side that has an explicit exact solution u(x, t) = sin(πt) sin(πx). The error at timeT = 100T = 59.2 is 1.4 × 10−3 with the ‘high’ order scheme compared with 4.3 × 10−2 for the first order scheme.A pure spectral Galerkin method was here used for the spatial discretization.A further generalization, where we extend to the fullest the notion of a high order method, can be proposed. Itconsists of proposing a spectral method in time, and we refer to [4] for a first implementation that requires one towork in the complex field.AcknowledgementPart of this work has been done in the frame of the grant ANR-06-CIS6-007.References[1] T. Bjøntegaard, Y. Maday, E.M. Rønquist, Fast tensor-product solvers. Part I: Partially deformed three-dimensional domains, submitted forpublication, 2007 and http://www.ann.jussieu.fr/publications/2007/R07023.html.[2] W. Hackbusch, Parabolic multi-grid methods, in: R. Glowinski, J.-L. Lions (Eds.), Computing Methods in Applied Sciences and Engineering,vol. IV, North-Holland, 1984, pp. 189–197.[3] J.-L. Lions, Y. Maday, G. Turinici, A parareal in time discretization of pde’s, C. R. Acad. Sci. Paris, Ser. I 332 (2001) 661–668.118 Y. Maday, E.M. Rønquist / C. R. Acad. Sci. Paris, Ser. I 346 (2008) 113–118[4] Y. Maday, E.M. Rønquist, Fast tensor-product solvers: Part II: Spectral discretization in space and time, submitted for publication, 2007 andhttp://www.ann.jussieu.fr/publications/2007/R07038.html.[5] J. Nievergelt, Parallel methods for integration ordinary differential equations, Comm. ACM 7 (1964) 731–733.

Parallelization in time through tensor-product space–time solvers

Hal-Diderot

Parallelization in time through tensor-product
space-time solvers
Yvon Maday, Einar M Rønquist
To cite this version:
Yvon Maday, Einar M Rønquist. Parallelization in time through tensor-product space-time
solvers. 2007. <hal-00175047>
HAL Id: hal-00175047
https://hal.archives-ouvertes.fr/hal-00175047
Submitted on 29 Sep 2007
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destine´e au de´poˆt et a` la diffusion de documents
scientifiques de niveau recherche, publie´s ou non,
e´manant des e´tablissements d’enseignement et de
recherche franc¸ais ou e´trangers, des laboratoires
publics ou prive´s.
Parallelization in time through tensor-product space-time solvers
Yvon Maday1, Einar M. Rønquist2
1 - Universite´ Pierre et Marie Curie-Paris6, UMR 7598, Laboratoire J.-L. Lions, Paris, F-75005 France and Division of
Applied Mathematics, Brown University, 182 George Street, Providence, RI 02912, USA
2 - Department of Mathematical Sciences, Norwegian University of Science and Technology, N-7491 Trondheim, Norway
Abstract
In this note, we extend the fast tensor-product algorithm for the simulation of time-dependent partial differential
equations. We use the natural tensorization of the space-time domain to propose, after discretization, a set of
independent problems, each one with the complexity of a single steady problem. This allows for an efficient parallel
implementation that is already interesting on small architectures, but that can also be combined with standard
domain-decomposition-based algorithms providing a further direction of parallelism on large computer platforms.
Preliminary numerical simulations are presented for a one-dimensional unsteady heat equation.
Re´sume´
Dans cette note on ge´ne´ralise a` la simulation de phe´nome`nes instationnaires l’algorithme de produit tensoriel.
On utilise la tensorisation naturelle du domaine espace-temps pour proposer, apre`s discre´tisation un ensemble de
proble`mes inde´pendants, chacun d’eux ayant la complexite´ d’un simple proble`me stationnaire. Ceci permet une
mise en œuvre paralle`le de´ja` interessante sur des petites architectures mais elle peut eˆtre e´galement combine´e avec
des techniques classiques de de´composition de domaine pour utiliser au mieux des architectures avec un nombre
de processeurs plus important. Des premiers re´sultats sont pre´sente´s sur un proble`me de la chaleur instationnaire
monodimensionnel.
Version franc¸aise abre´ge´e
La re´solution nume´rique de proble`mes mode´lise´s par des e´quations aux de´rive´es partielles est maintenant
assez bien e´tablie pour des phe´nome`nes stationnaires. Les techniques ite´ratives comme le “multishooting”
de [5], les approches multigrilles de [2], et plus re´cemment les approches parare´elles en temps de [3]
permettent de traiter des phe´nome`nes instationnaires de fac¸on efficace en tenant compte des architectures
des calculateurs paralle`les. La me´thode non ite´rative, que nous pre´sentons ici est une ge´ne´ralisation des
techniques de tensorisation pre´sente´es par exemple dans [1] pour des discre´tisations en espace seulement ;
Email addresses: maday@ann.jussieu.fr (Yvon Maday1), ronquist@math.ntnu.no (Einar M. Rønquist2).
Preprint submitted to Elsevier Science
elle peut eˆtre combine´e a` des techniques classiques de de´composition de domaine en espace pour une plus
grande efficacite´.
Conside´rons par exemple l’e´quation de la chaleur (1, 2, 3). Apre`s discre´tisation en espace par une
me´thode classique, par exemple de diffe´rences finies de pas ∆x = L/(N + 1) et une discre´tisation en
temps implicite, on obtient le syste`me discret (5) ou` le vecteur um ∈ RN (resp. fm ∈ RN ) repre´sente
l’approximation de la solution (resp. la donne´e de flux de chaleur) aux nœuds internes et au temps tm, et
∆tm = t
m − tm−1 est le pas de temps au niveau tm. La matrice A ∈ RN×N est classiquement syme´trique
de´finie positive.
A` chaque pas de temps on doit donc inverser l’ope´rateur de Helmholtz discret A + 1∆tm I et la de´ter-
mination de um−1 est ne´cessaire pour calculer um. Pour briser la nature intrinse`quement se´quentielle de
cette re´solution, on utilise l’algorithme de produit tensoriel rapide. Soit A ∈ RN×N et B ∈ RM×M ; on
note tout d’abord C = B ⊗ A le produit tensoriel de B et A : c’est la matrice C ∈ R(MN)×(MN) de
coefficients cij = bklamn, 1 ≤ k, l ≤ M , 1 ≤ m,n ≤ N , i = m + (k − 1)N et j = n + (l − 1)N . Avec ces
notations, le syste`me discret (5) s’e´crit (8) ou` B est de´fini en (6), It est la matrice identite´ de dimension
M et Ix est la matrice identite´ de dimension N , de plus u ∈ R
MN (resp. f ∈ RMN ) repre´sente le vecteur
concate´ne´ de la solution et du membre de droite selon (7) a` tous les pas de temps. Si l’on suppose tous
les ∆tm, m = 1, . . . ,M diffe´rents, la matrice B est diagonalisable, selon (9). Le syste`me (8) peut donc
prendre la forme (13) ou` l’inconnue w est de´finie en (12) et le membre de droite g est de´fini en (15). On
remarque aise´ment que ce syste`me se de´ploie en M proble`mes inde´pendants (14) qui peuvent donc eˆtre
re´solus en meˆme temps sur des processeurs paralle`les, et la solution {um,m = 1, . . . ,M} s’obtient par la
transformation inverse (16).
En notant par WH la complexite´ de re´solution d’un syste`me de type (5), la complexite´ de re´solution
de l’approche tensorise´e (non paralle´lise´e) est donne´e en (17), en effet la diagonalisation de la matrice
B compte pour O(M2) et le calcul des membres de droite (15) ainsi que la transformation inverse (16)
compte pour O(M2N). Si l’on dispose par contre de P = M processeurs paralle`les, la complexite´ par
processeur devient (18), en ajoutant une complexite´ de communication entre les processeurs Wcomm on
obtient une efficacite´ de paralle´lisme S donne´e par (19). Ainsi si l’on suppose d’une part queWcomm ≪WH
et en outre M ≪ N et MN ≪WH on conclut a` une efficite´ quasi optimale S ≈M .
La tensorisation expose´e ci dessus peut eˆtre effectue´e pour des sche´mas en temps plus pre´cis que le
sche´mas implicite d’ordre 1 de Euler, tant que la matrice B reste triangulaire infe´rieure. La simulation
nume´rique suivante base´e sur un sche´mas implicite re´trograde d’ordre 3 (avec un choix de pas de temps
∆tk = ρ
k−1∆t1, k > 1, avec ρ = 1.2 se sorte que ∆t15 ≈ 13∆t1, et ∆T =
∑M=15
k=1 ∆tk = 0.592). L’erreur
obtenue apre`s un temps T e´gal a` T = 100∆T = 59.2 est 1.4 · 10−3 alors qu’elle est de 4.3 · 10−2 avec les
sche´mas d’ordre 1.
1. Introduction
The simulation of three-dimensional steady problems governed by partial differential equations is nowa-
days quite common; the size of current computers allows one to get very good approximations, close to
experimental data, for a large class of problems. The availability of parallel computers has resulted in
many contributions related to the parallelization of problems with respect to the spatial directions. A key
ingredient has been to divide the global problem through domain decomposition, thus allowing smaller
problems to be assigned to different processors.
The next challenge, which has already started to be faced, is to solve time-dependent problems with the
same efficiency and ease. Roughly speaking, this amounts to adding an additional direction to the three
spatial ones. Most discretizations of time-dependent problems rely on using finite difference methods
2
in time. This results in the need to solve a sequence of stationary-like problems, for which domain
decomposition is primarily the only way to exploit parallel computers. Indeed, the sequence of these
stationary-like problems is recursive and is not obviously broken into independent subproblems. Among
the emerging new methods, let us quote the multishooting technique introduced by Nievergelt [5], and the
multigrid approach introduced by Hackbusch [2]. More recently a new approach, named the parareal in
time method has been introduced by Lions, Maday and Turinici [3] with various contributions showing
the interest of this predictor-corrector methodology.
The approach we present here is of a completely different nature. First of all it is not iterative, at least
in the version we propose here to tackle linear problems. Second, it may be of interest already with very
few processors. Third, as a consequence of the diversity of available solvers, it can be combined with
the parareal approach as we shall explain in the full paper associated with this Note. The method used
here is purely algebraic. It is inspired by the fast tensor product solver that we introduced in the recent
paper [1] for spatial discretizations, and has also some connection with the classical way that leads to the
harmonic version of a time dependent problem, and that looks for solutions of the form eiωtu(x). Finally,
let us indicate that it can easily be combined with standard domain decomposition methods in space,
thus offering enhanced overall parallelism.
2. The heat equation
Consider the one-dimensional unsteady heat equation
∂u
∂t
= κ
∂2u
∂x2
+ f(x, t) in ωx = (0, L), (1)
subject to the homogeneous boundary conditions
u(x = 0, t) = u(x = L, t) = 0 ∀ t ∈ ωt = [0, T ], (2)
and to the initial condition
u(x, t = 0) = g(x). (3)
Here, u(x, t) is the solution (the temperature), κ is the thermal diffusivity, and f(x, t) is a thermal heat
source that depends on space and time. We assume that the thermal diffusivity, κ, is a constant.
As an introductory example to the approach we will present, we consider the numerical solution of this
problem using finite difference methods in space. The spatial domain ωx is discretized using a uniform
grid with grid spacing ∆x = L/(N + 1). Using a standard second order finite difference formula for the
diffusion term, we end up with a system of ordinary differential equations of the form
du(t)
dt
+Au(t) = f(t). (4)
Here, the vector u ∈ RN represents the numerical solution at the N internal grid points (at a particular
time t), the vector f ∈ RN represents the given data at the internal grid points (at a particular time t),
and the matrix A ∈ RN×N , that represents the discrete diffusion operator, is symmetric positive definite.
Assume now that we discretize the system of ordinary differential equations using an implicit method
like the first order backward Euler scheme. The set of fully discrete equations can then be expressed as
um − um−1
∆tm
+Aum = fm , m = 1, . . .M, (5)
where um denotes the numerical approximation of the solution at the internal grid points at time tm and
∆tm is the time step at time level t
m; more precisely, ∆tm = t
m − tm−1.
3
This simple approach to obtain the numerical solution results in the need to solve a system of algebraic
equations at each time step, with the discrete Helmholtz operator A+ 1∆tm I as the system matrix. Note
that this is indeed a purely sequential approach; we need to know um−1 before we can compute um by
solving a discrete Helmholtz problem.
3. Tensor-product representation
We will now present a reformulation of the model problem introduced in the previous section. This
reformulation is based on a tensor product representation of time and space and will allow us to obtain
a scheme with parallel processing possibilities in the time domain.
First, let us briefly review the tensor product notation that we will be using. Let A ∈ RN×N and
B ∈ RM×M . Then the tensor product of B and A, C = B ⊗ A , gives a matrix C ∈ R(MN)×(MN) with
matrix entrees cij = bklamn, with 1 ≤ k, l ≤M , 1 ≤ m,n ≤ N , i = m+ (k − 1)N and j = n+ (l − 1)N .
Associated with the backward Euler scheme (5), we define the matrix B to be
B =


1
∆t1
−
1
∆t2
1
∆t2
−
1
∆t3
1
∆t3
. . .
−
1
∆tM
1
∆tM


. (6)
We also define It to be the identity matrix of dimension M (associated with the time domain) and Ix
to be the identity matrix of dimension N (associated with the spatial domain). Finally, we define the
concatenated solution vector u ∈ RMN and the concatenated right hand side f ∈ RMN as
u =


u1
u2
u3
.
uM


, f =


f1
f2
f3
.
fM


. (7)
With this notation we can write the (5) in tensor product notation as
(B ⊗ Ix + It ⊗A) u = (It ⊗ Ix) f . (8)
Note that this system of algebraic equations represents the numerical solution at M different time levels.
Although the system matrix is sparse, it does not a priori seem like a good strategy to reduce the overall
computational cost.
In order to see the potential benefit of this approach, we proceed as follows. We first assume that it is
possible to diagonalize the matrix B, i.e., we assume that we can write B as
B = SΛS−1. (9)
4
Here, the columns of S are the eigenvectors ofB, while Λ is a diagonal matrix comprising the corresponding
eigenvalues of B. For our particular Euler backward scheme, the matrix B is lower triangular, and the
eigenvalues are equal to the diagonal entries of B. Hence, the matrix B is diagonalizable only if all the
time steps ∆tm, m = 1, . . . ,M are different. Assuming that this is the case, we can write (8) as
(SΛS−1 ⊗ Ix + It ⊗A) u = (It ⊗ Ix) f , (10)
or
(S ⊗ Ix) (Λ⊗ Ix + It ⊗A) (S
−1 ⊗ Ix) u = (It ⊗ Ix) f. (11)
Introducing
w = (S−1 ⊗ Ix) u (12)
and premultiplying (11) by (S−1 ⊗ Ix) gives
(Λ⊗ Ix + It ⊗A) w = (S
−1 ⊗ Ix) f. (13)
An equivalent formulation of (13) is thus
(A+ λmIx) w
m = gm, m = 1, . . . ,M, (14)
with a concatenated right hand side
g = (S−1 ⊗ Ix) f. (15)
We have here reformulated the solution in terms of M computational steps. However, in contrast to the
systems (5) which represent M sequential steps, the systems (14) represent M independent steps. Hence,
(14) can be solved in parallel in a completely decoupled fashion without any communication.
Once the systems (14) have been solved for wm, m = 1, . . . ,M , we compute um, m = 1, . . . ,M from
(12), i.e.,
u = (S ⊗ Ix) w. (16)
4. Computational complexity
We now discuss the computational complexity associated with the approach presented in the previous
section. We recall that the model problem is a simple unsteady heat transfer problem discretized with a
standard finite difference method in space and a backward Euler scheme in time. The extension to use
finite element methods in space is straightforward.
The computational approach assumes that it is possible to diagonalize the matrix B defined in (6).
The cost of the diagonalization (9) is O(M2) where M is the number of (standard) time steps following
a backward Euler approach. The corresponding memory requirement is O(M2).
Next, the cost of computing the transformed right hand sides (15) is O(M2N), while the cost of solving
the M systems in (14) is the same as the standard approach for solving the M systems (5); the particular
cost here depends on the chosen solver for these Helmholtz-type systems. Let us denote the cost per
Helmholtz system as WH . The cost of performing the final transformation (16) is O(M
2N).
The total (serial) cost W1 of the new approach is therefore
W1 = c1 ·M
2 + c2 ·M
2N +M ·WH . (17)
The cost W˜1 of following a standard approach is just the last term on the right hand side.
However, assuming now the we have available P = M processors, the computational complexity per
processor is
WP = c1 ·M
2 + c2 ·M N +WH . (18)
5
In addition, we have a communication cost, Wcomm, associated with collecting the concatenated vector w
on each processor (all-to-all communication), which is needed in order to compute um on each processor
Pm, m = 1, . . . ,M ; see (16). The speedup S is thus
S =
W˜1
WP +Wcomm
=
M ·WH
c1 ·M2 + c2 ·M N +WH +Wcomm
. (19)
If Wcomm ≪WH , M ≪ N , and the work WH satisfies MN ≪WH , we see that S ≈M , i.e., we get close
to perfect speedup. The assumption Wcomm ≪WH can be expected to be satisfied when considering the
extension of this approach to solve more realistic multi-dimensional problems. In this case, the assumption
M ≪ N can also easily be fulfilled, and the assumption MN ≪WH follows in most cases.
5. Generalization and numerical results
The method we have just proposed is first order in time, and since it requires that all the time steps
are different, the accuracy will be related to the largest time step. In order to make the method more
efficient, we propose to use a higher order scheme in time with time steps ∆tk = ρ
k−1∆t1, k > 1, with
ρ larger but close to 1, e.g. ρ = 1.2 so that ∆t15 ≈ 13∆t1, and ∆T =
∑M=15
k=1 ∆tk = 0.592. Note that,
as can be expected, choosing ρ much closer to 1 may lead to instabilities due to numerical errors. The
scheme we use here comprises a second order Crank-Nicolson scheme for the first time step, followed by
a Backward Difference scheme of order 2 to get u2, then a Backward Difference scheme of order 3 to get
uk for 2 < k ≤ M = 15. The scheme is a little bit more involved but results in a matrix similar to B
defined in (6) since B is still lower diagonal, with few terms (actually 3) under the diagonal, and the same
diagonal terms as in (6). The complexity is identical as for the first order version we analyzed. We have
performed the simulation for a right hand side that has an explicit exact solution u(x, t) = sin(pit) sin(pix).
The error at time T = 100∆T = 59.2 is 1.4 · 10−3 with the ”high” order scheme compared with 4.3 · 10−2
for the first order scheme. A pure spectral Galerkin method was here used for the spatial discretization.
A further generalization, where we extend to the fullest the notion of a high order method, can be
proposed. It consists of proposing a spectral method in time, and we refer to [4] for a first implementation
that requires one to work in the complex field.
Acknowledgements. Part of this work has been done in the frame of the grant ANR-06-CIS6-007.
References
[1] T. Bjøntegaard, Y. Maday, E.M. Rønquist Fast tensor-product solvers. Part I: Partially deformed three-dimensional
domains, submitted 2007 and http://www.ann.jussieu.fr/publications/2007/R07023.html.
[2] W. Hackbusch, Parabolic multi-grid methods, in Computing Methods in Applied Sciences and Engineering, VI, R.
Glowinski and J.-L. Lions, eds., North-Holland, 1984, pp. 189–197.
[3] J.-L. Lions, Y. Maday, and G. Turinici, A parareal in time discretization of pde’s, C.R. Acad. Sci. Paris, Serie I, 332
(2001), pp. 661–668.
[4] Y. Maday, E.M. Rønquist Fast tensor-product solvers: Part II: Spectral discretization in space and time, submitted
2007 and http://www.ann.jussieu.fr/publications/2007.
[5] J. Nievergelt, Parallel methods for integration ordinary differential equations, Comm. ACM, 7 (1964), pp. 731–733.
6


Parallelization in time through tensor-product space-time solvers

Rønquist, Einar M.

Numérisation de Documents Anciens Mathématiques

https://hal.archives-ouvertes.fr/hal-00175047

Parallelization in time through tensor-product space-time solvers

Abstract

Similar works

Full text

Available Versions

Comptes Rendus Mathématique

Hal-Diderot

Numérisation de Documents Anciens Mathématiques