Numerical experiments with nonlinear multigrid waveform relaxation on a parallel processor

Applied Numerical Mathematics 8 (1991) 149-161 North-Holland

149

Numerical experiments with nonlinear multigrid waveform relaxation on a parallel processor

Stefan Vandewalle * and Robert Piessens Katholieke Universiteit Leuven, Department of Computer Science, Celestijnenlaan 2OOA, B-3001 Leuven, Belgium

Abstract

Vandewalle, S. and R. Piessens, Numerical experiments with nonlinear multigrid waveform relaxation on a parallel processor, Applied Numerical Mathematics 8 (1991) 1499161.

Standard time-stepping techniques for solving parabolic partial differential equations cannot be parallelized or vectorized efficiently unless the problem to be solved is very large and the number of processors is small. To overcome these limitations several researchers have presented methods that increase the parallel performance by calculating the solution on several time levels simultaneously. In earlier papers we have discussed such a method, based on a combination of the waveform relaxation method and multigrid. In this paper we extend our approach to the case of nonlinear parabolic partial differential equations. We illustrate the application of the new algorithm with numerical examples of initial-boundary-value and time-periodic type and we discuss its implementation on an Intel iPSC/2 hypercube. It is shown that also in the nonlinear case the multigrid waveform relaxation method outperforms the standard sequential solvers.

1. Introduction

In a previous study we compared the parallel performance of several standard time-stepping techniques for solving linear parabolic partial differential equations, [16]. It was shown that the standard methods can only be parallelized efficiently for problems that are large enough. More precisely, the arithmetic complexity of the calculations per time step and per processor must outweigh the communication overheads. We argumented that new parallel algorithms are needed for solving relatively small parabolic problems, i.e., problems with few unknowns per time level or problems solved on large scale parallel machines.

A direction taken by several researchers is to improve the multiprocessor performance by operating on several or on all time levels at once, see e.g. [1,12,20]. A related algorithm was published in [6] by Lubich and Ostermann. Their approach is different from the previous ones in that it can be defined without explicit reference to any time discretization technique or to any time levels. The algorithm, which is described for solving linear parabolic problems, is based on the use of the multigrid method and waveform relaxation, a technique for solving large systems

* Senior research assistant of the National Fund for Scientific Research (Belgium)

0168-9274/91/$03.50 0 1991 - Elsevier Science Publishers B.V. All rights reserved

150 S. Vandewalle, R. Piessens / Nonlinear multigrid waveform relaxation

of ordinary differential equations. The same algorithm was used by the current authors in [15,17]. We compared the technique to standard time-stepping methods and presented numerical results and timings obtained with an implementation on a parallel machine. A straightforward extension of the method was found to lead to a very efficient technique for solving linear time-periodic equations, [17]. A theoretical analysis of the resulting method, which we called periodic multigrid waveform relaxation, was given in [18]. In this paper we will present the extension of the multigrid waveform relaxation method to nonlinear problems of initial-boundary-value and time-periodic type. A short discussion of the nonlinear initial-boundary-value case already appeared in an earlier paper, [14]. There, we gave a small example to illustrate the basic concepts. However, we did not give any parallel performance results, nor did we compare the method to standard techniques.

In Section 2, we briefly discuss the standard waveform relaxation method and we direct the interested reader to further references where the method is more completely covered. The nonlinear multigrid waveform relaxation algorithm is presented in Section 3. In Section 4 we consider several implementation aspects and we address the issues of parallelism and vectorization. In Section 5 we illustrate the method by two examples and we compare its performance to that of a parallel implementation of the standard solution methods.

2. The waveform relaxation method

Waveform relaxation, also called dynamic iteration or Picard-Lindelof iteration, is an iterative technique for solving systems of ordinary differential equations. It is conceptually very similar to the standard relaxation techniques for solving systems of linear and nonlinear algebraic equations. Consider the following general system of N first-order differential equations,

&=m .Yl,.--, YN) fortE[t,, tf], i=l,..., N, (2.1)

supplemented with either an initial condition at time t,, y,( to) = ylo, or a time-periodicity condition, y,( to) = y,(t,). We assume that the system is well-posed, i.e., a solution exists and is locally unique. In particular, in the case of a time-periodic problem we assume that the period of the solution is known a priori, in casu equal to t, - t,.

The Gauss-Seidel variant of the waveform relaxation algorithm for solving (2.1) is stated in Fig. 1. The system of N differential equations is not solved time step after time step. Instead, the differential equations are repeatedly solved one after the other, as equations with a single unknown. In the case of a nonlinear problem these equations don’t need to be solved exactly. The solution may be approximated, e.g., by solving the linear differential equation that arises after application of a Newton linearization around a sufficiently close approximation. This technique is usually referred to as the Newton waveform relaxation method.

The theoretical foundations of the waveform relaxation method have been discussed in a number of papers. A convergence analysis for linear systems is given by Miekkala and Nevanlinna in [7]. They study the Gauss-Seidel scheme and the variants of Jacobi, SOR and SSOR type. In [19] convergence is proven for nonlinear systems. In [5] the relation is established between the number of iterations and the accuracy order, i.e., the number of correct terms in the

S. Vandewalle, R. Piessens / Nonlinear muhigrid waveform relaxation 151

choose y,‘“‘(t) for t E [t,, t,], i = 1,. . . , N

repeat foreachi=l,...,N

solve dy”‘+“/dt =!;(I, yl(“+r), . . . ) y,‘“,“‘, yp+y y,it”‘l,...

with either y,(“+‘)( to) = yiO or yi(‘+‘)( to) = Y:~+“( tf) n:=n+l

until convergence.

Fig. 1. The Gauss-Seidel waveform relaxation algorithm.

Taylor expansion of the partially converged solution. The convergence of the time-periodic iteration is analyzed in [18].

The main advantage over standard techniques is the applicability of waveform relaxation as a multirate time-integration method. Indeed, each equation can be solved by using the time step needed to reflect the time behaviour of the local variable. The method has proven to be very effective in solving the large systems of differential equations that describe the behaviour of complex VLSI devices, [8,10,19]. For this application the method results in considerable savings, as the variables change at very different rates and as many variables are inactive during most of the simulation time. Further applications are found in the literature on the simulation of chemical distillation processes, [13]. A second advantage is the ease with which the method can be parallelized on shared-memory and distributed-memory multiprocessors, [8,10].

3. Multigrid waveform relaxation

Application of the numerical method of lines to a parabolic partial differential equation leads to a large system of ordinary differential equations,

& + L,(u,) = Fh. (3.1)

In the above formula we used h to denote the spatial mesh size; u,, is the vector of unknown functions that approximate the solution at the grid points; L, is the discretized, possibly nonlinear spatial operator; Fh is a known right-hand side vector, which is derived from given functions in the differential equation and the boundary conditions. Equation (3.1) is completed with either an initial condition, uh( to) = u~,~, or a time-periodicity condition, am = u,(t,).

The waveform relaxation method as discussed in Section 2 may be applied to solve the system of ordinary differential equations (3.1). However, this does not lead to a satisfactory algorithm due to slow convergence. For the discretized heat equation, the convergence rates of the Jacobi and Gauss-Seidel schemes are of order 1 - 0( h2). This was shown in [7] for the initial-value case and in [18] for the time-periodic case. A similar convergence behaviour is observed for more difficult and nonlinear problems. The convergence rapidly deteriorates with an increasing number of grid lines. In contrast to the case of a system of equations arising from a discretized


elliptic equation, overrelaxation in a SOR fashion does not significantly improve convergence properties.

The convergence can be accelerated if waveform relaxation is combined with the multigrid idea. Each of the standard multigrid operations is then replaced by a similar operation defined to operate on functions. The linear variant, based on defect correction multigrid, was described and analyzed in [6]. The method that we consider here is derived from the well-known full approximation scheme, which is a standard multigrid method for solving nonlinear elliptic problems, [3]. It reduces to the defect correction scheme when it is applied to a linear problem. The technique is based on the interplay of fine-grid smoothing and coarse-grid correction. The sequence of operations of the two-grid algorithm is given below. In the description we will use subscript h to denote functions on the fine grid, and H for the coarse grid. We assume that an initial approximation, U,,, is given on the fine grid. It may be derived from a known solution of a related problem, or it may be derived by using the nested iteration idea, which is discussed further. The following cycle of three steps (pre-smoothing, coarse-grid correction, post-smoothing) is repeatedly applied in order to update the approximation U,.

l &-e-smoothing. Apply a number of Gauss-Seidel waveform relaxation steps to the current approximation U,. When the problem is nonlinear, a Newton waveform relaxation technique can be used. As in the case of an elliptic problem, the ordering of the equations may be lexicographic or red-black. The latter is usually preferred for a parallel implementation.

l Coarse-grid correction. - Project the current approximation, U,, onto the coarse grid by using a restriction

operator, JhH,

u, == J/Uh. (3.2)

The waveform relaxation restriction operators are straightforward extensions of similar restriction operators used when solving elliptic problems. They are different in that they operate on functions instead of on scalar values. E.g., the one-dimensional full-weighting restriction operator reads like

u,(t) := +$+1(t) + w(t) + U;-1(t)), (3.3)

where u1 and ui are functions at corresponding grid points on the coarse and on the fine grid.

- Calculate the right-hand side of the problem on the coarse grid,

The restriction operator 1: may be different from the previously used operator jhH. However, in order to avoid the derivative calculation in (3.4) the use of the same restriction operator may be advantageous. Indeed, in that case the two derivatives cancel.

- Solve the following system of ordinary differential equations on the coarse grid,

&+L#H) =FH (3.5)

with either an initial condition, UH(tO) = i,“U,,( to) or a time-periodicity condition,

&(IO) = UH(tf).

S. Vandewalle, R. Piessens / Nonlinear multigrid waveform relaxation 153

- Calculate the correction. Interpolate the correction to the fine grid by using a waveform extension of a standard prolongation formula, IL, and correct the current fine grid approximation,

u,:= u,+rgu,- i&H). (3.6)

l Post-smoothing. Perform a number of smoothing operations, e.g., again by Gauss-Seidel waveform relaxation.

This two-grid cycle may be applied recursively to solve the coarse grid equation. This leads to the multigrid scheme. The algorithm is completely defined by specifying the sequence of nested grids, the discretized operators, the inter-grid transfer operations (restriction and prolongation), the nature of the smoothing relaxations, and by specifying the order and the frequency with which the different grids are visited. For a discussion of the various options and parameter choices in the case of an elliptic equation we refer to [3,4]. The discussion does not change qualitatively if we consider parabolic problems. Finally, the algorithm can be combined with the idea of nested iteration. The initial approximate solution to the problem on fine grid is then derived by interpolation from an approximate solution calculated on the coarse grid. This leads to the waveform equivalent of the so-called full multigrid method.

4. Implementation aspects

4.1. Time discretization and time integration

The solution to the ordinary differential equations in the waveform relaxation method will usually be calculated by a numerical approximation procedure. We will discuss the case of solving an initial-boundary-value problem first. Any stiff integrator can be used to solve the differential equations that arise in the Gauss-Seidel relaxation step. A variable-stepsize variable-order integrator is often used in engineering applications when systems are modeled with very disparate time characteristics. Each differential equation may then be solved using the time steps needed to accurately resolve the local variable. In our implementation we have not taken advantage of this possibility. The properties of waveform relaxation as a multirate integration method for solving parabolic partial differential equations remain to be analyzed in a further study. Instead, for the examples in the next section we used a straightforward implementation of the second-order accurate trapezoidal rule (Crank-Nicolson method), with fixed and a priori determined time step. Furthermore the nonlinearity is handled by a Newton waveform relaxation procedure, where the equation is linearized and solved once. This choice of time integration greatly simplifies implementation and, as will be explained further, it allows for vectorization.

The (linearized) two-point boundary-value problems that arise in the smoothing step of the periodic multigrid waveform relaxation method may be solved, e.g., by shooting or by a discretization technique. We opted for the latter and applied the trapezoidal rule with constant time-increment. This leads to an easily solvable system of linear equations with a cyclic bidiagonal coefficient matrix.

In the waveform method several functions are associated with each grid point. For these functions, an effective and compact computer representation should be chosen. It should allow


for efficient evaluation of the functions and, more importantly, for efficient algebraic manipula- tion (mainly the summation of functions and multiplication with a scalar). For the latter reason we have represented each function as a vector of values evaluated at equidistant time levels. We denote the vector length by n, (number of time levels). For simplicity, the time-increment is chosen equal to the time-increment used by the integration formula.

In the case of a standard time-stepping procedure the unknowns need to be stored in the computer memory on one or two time levels only. In the case of waveform relaxation the values at every time level need to be available. Storage requirements are therefore very high. However, if storage is a problem, the time interval of interest may be partitioned into several smaller windows, which can be treated in sequence (in the case of solving an initial-boundary-value problem).

4.2. Parallelitation

The parallel implementation of the linear waveform relaxation algorithm and the paralleliza- tion of the standard time-stepping methods is discussed in detail in our studies [15,16]. The implementation of the nonlinear algorithm is basically very similar. Therefore we will only briefly review some of the main issues. The multigrid operators, i.e., pointwise relaxation, defect calculation needed in the determination of the coarse-grid problem right-hand side, restriction, prolongation and correction, are local operators, i.e., they only affect neighbouring grid points. The latter four operators can be parallelized straightforwardly as the operations at each grid point may be executed independently from the operations on the other grid points. The Gauss-Seidel relaxation process is parallelized by ordering the equations based on the idea of grid point colouring. In the case of a five-point finite difference stencil the grid points are coloured alternatively red and black, thus dividing them into two disjunct sets. The equations associated with the grid points in one set may be relaxed concurrently, as they only depend on the values associated with grid points in the other set. The global relaxation step is thus divided into two partial relaxation steps which are executed in sequence. A similar four-colour scheme is usually applied in the case of a nine-point finite difference stencil.

A classical data decomposition is used to evenly distribute the computational workload. The processors are arranged in a rectangular array and are mapped onto the domain of the partial differential equation, which, in our implementation, is restricted to be a rectangle. Each processor is responsible for doing all computations on the grid points in its subdomain. During the computation, communication with neighbouring processors is needed mainly to update local boundary values. Some other communication strategies are further used to improve the parallel performance. We’d like to mention in particular the use of an agglomeration strategy to reduce the communication complexity of the coarse grid operations, [16]. For a quantitative discussion and comparison of the communication complexities of the waveform relaxation algorithm and the corresponding standard time-stepping method we refer to [17]. The discussion and formulae given there immediately extend to the nonlinear case.

4.3. Vectorization

The use of vector processors for executing the multigrid waveform algorithm may result in a substantial reduction of computing time. Indeed, most of the operators can be expressed as


simple arithmetic operations on functions, i.e., on uectors. In contrast to the standard approach uectorization is in the time direction and not in the spatial direction! The vector speedup of the arithmetic part of the computation will mainly depend on the value of the function length, n,. It will be virtually independent of the size of the spatial grid, the number of multigrid levels, the multigrid cycle used and the number of processors. This contrasts sharply with standard vectorization, which only leads to a performance improvement when the number of grid points per processor is very large.

The only operation which is not perfectly vectorizable is the core of the ODE integrator used in the smoothing step. It consists of one or two first-order recurrence formulae of length n, (one in the case of an initial-value problem and two in the case of a time-periodic problem solved by discretization with a one-step formula). This reduces the possible gain through vectorization. It can be shown, however, that the cost of the recurrence relation makes out only a small fraction of the total computation cost, typically 5 to 10%. This leads to a possible vector speedup of 10, or more.

5. Numerical examples

We illustrate in this section the numerical and parallel properties of the waveform relaxation method discussed in Section 3. We have implemented the algorithm on an Intel iPSC/2 hypercube (without vector processors). For a description of this multiprocessor, its hardware characteristics and various performance benchmarks, we refer to [2].

5.1. A nonlinear initial-boundary-value problem

A solid body with initial temperature u,,, is submerged into the environment with temperature

u, and starts losing heat through radiation along its sides. This heat transfer problem is governed by the following quasi-linear parabolic partial differential equation,

peg = div( k( u)grad( u)),

supplemented with the nonlinear radiation boundary condition,

k(u)& = ee( 2.42 - u4),

(5.1)

(5 4

where u is the temperature distribution [K], p is the density [kg/m3], c is the specific heat [J/kg K], k is the thermal conductivity [W/m K], e is the emissivity [ -1 and E is the Stefan-Boltz- mann constant (5.75 lop8 [W/m2 K4]), see e.g. [9]. With a/an we denote the outside normal derivation operator.

We will numerically determine the time-accurate temperature distribution for t E [0, T] of a square two-dimensional surface of size I, that is

(x, y) E [-:r, +1] x [-:I, :,I

which satisfies (5.1) and (5.2). Thanks to symmetry the computational domain may be restricted to [0, :L] X [0, $11 with homogeneous-flux boundary conditions, i.e., au/an = 0, imposed along

156 S. Vandewalle. R. Piessens / Nonlinear multigrid waveform relaxation

6)

Gauss-Seidel WR

k-0

0 5 10 15 20 (min) 0 5 10 15 20 (min) 0 5 10 15 20 (min)

Fig. 2. Waveform relaxation for time integration of (5.1) with boundary conditions (5.2). The figure shows the approximations of the temperature u(k) (t, l/2, l/2) after k iteration steps.

Multigrid WR Full Multigrid WR

,(K)

the sides x = 0 and y = 0. The following numerical values are chosen for the problem parame- ters:

2.40 = 473, u, = 3, p = 7000,

c = 400, e= 1, k(u) = 50(1- 0.0008U),

I= 0.5, T= 1545 s.

The problem is discretized by standard second-order central differences on an equidistant grid with 17 or 33 grid lines in x- and y-direction (n, = ny).

In Fig. 2 we illustrate the convergence behaviour of waveform relaxation (WR). The figure to the left shows the time profile of successive approximations U(~) evaluated at a corner point of the domain, (x, y) = (ir, :I), and generated by the WR Gauss-Seidel algorithm. For all grid points (n, = nY = 17 in this case) a constant time profile equal to the initial condition was used as the starting approximation. As can be seen from the figure the convergence of the algorithm is very slow. The time interval where the approximation is satisfactory only gradually extends as more relaxations are applied. Corresponding results obtained with multigrid WR are shown in the middle figure, which clearly illustrates a dramatic improvement in convergence. Multigrid V-cycles are used with one red-black WR pre-smoothing step and one similar post-smoothing step. Standard coarsening is performed down to a coarse level with one unknown. The differential equation at each grid point is integrated in the smoothing step by solving the linear equation that results after applying one Newton linearization. The computational cost of one such V-cycle is approximately equal to S/3 work-units, with one work-unit being the cost of one WR Gauss-Seidel iteration. The figure to the right shows the starting approximation on the fine grid as obtained by the full multigrid procedure. The cost of this initial step is only about g/9 work-units. The subsequent approximations cannot be distinguished graphically from the initial approximation.


residual residual

1 discretization: n,=ny=17, n,=51 discretization: n,=n,=17, II,=51

10-l

10-z

1o-3

10”

10-S

lo6

lo-’

lo-*

10-P

-1

- 10-l

- 10-Z

-lo-)

- lo4

- 10-S

+$ - 104

,kc - lo-’

- 10-8

-lo+

I I I I 1 lb I 20 40 60 80 160 1 0 i 10 15 20 $5 :o

execution time on 1 processor (in set) execution time on 16 processors (in set)

residual discrelization: n,=n,=33, n,=51 discretization: n,=n,=33, n,=51

I

e

so b 50 I lb0 180 260 250 I 360 350 I 460 4

+ 4x

, 1 I I I I I , , 5 10 1s 20 25 30 35 40 4s

residual

1

10-l

10-Z

10-3

lOA

1o-5

10d

10-7

10”

10-9

execution time on 1 processor (in SK) execution tinge on 16 proct~~rs (III SK)

Fig. 3. Comparison of execution times on 1 and on 16 processors for solving problem (5.1) with multigrid waveform relaxation using V-cycles (solid line) and with the Crank-Nicolson time stepping method with fixed number of

V-cycles ( + ) or with full multigrid and variable number of V-cycles ( X ).

In Fig. 3 the performance of the multigrid waveform relaxation (MWR) method and the well-known Crank-Nicolson (C-N) method are compared. The latter is a standard time-stepping technique, which transforms the parabolic problem into a sequence of elliptic partial differential equations to be solved at successive time levels. In our implementation the nonlinear problem at each time level is solved by the multigrid full approximation scheme for nonlinear elliptic equations, [3]. To make a fair comparison of both methods we applied the C-N rule to integrate the ordinary differential equations that arise in the MWR smoothing step. We also used the same time step as the one used in the C-N method. Under these conditions a careful analysis


shows that the discrete equations of both techniques are similar. Consequently both solutions are identical. In this example a variable-stepsize time-integration is used with At, ( = ti - ti_,)

determined a priori. It is chosen equal to 5 seconds for the first 10 time steps and it is doubled every 10 steps. A total of 50 time steps is applied.

The accuracy, in casu the largest residual with respect to the discrete equations, is plotted versus the execution time on 1 and on 16 processors, for two different spatial discretizations. Figure 3 shows a smooth curue for the waveform method. The residual of the initial (constant) approximation decreases as more and more multigrid cycles are applied. The C-N results show up as discrete points. The calculations are advanced time step per time step in a total of t

seconds. The largest residual of the resulting approximation is represented by a “ +” or a “ X" sign at position (t, residual) in the figure. We have applied C-N in two different (but both common) ways. In one series of experiments a fixed number of V-cycles is applied to the initial approximation on each time level, which, here, is chosen equal to the solution obtained in the previous step. The results are denoted by “a” (1 V-cycle), “b” (2 V-cycles), “c” (3 V-cycles), “d” (4 V-cycles) or “e” (5 V-cycles). In the second /series of experiments an adaptive scheme is implemented. The initial approximation on each’time level is determined by the full multigrid method and multigrid V-cycles are applied until ‘the largest residual is below a certain specified tolerance. The corresponding results are denoted by “1” (tol. = 10e4), “2” (tol. = 10-5), “3” (tol. = 10F6) and “4” (tol. = lo-‘). Note that a measure of the largest residual may be calculated at negligible cost, e.g., by calculating the residual on a subgrid of the fine grid.

On one processor MWR turns out be as efficient as C-N when the latter uses a fixed number of V-cycles. The C-N method with the adaptive scheme is somewhat faster. On 16 processors however MWR clearly outperforms C-N. In the case of the 17 X 17 discretization the speedup of MWR is about 7.5, whereas the speedup of C-N is only 3. A speedup of 11 is achieved by MWR for the 33 x 33 problem, while the C-N speedup is about 6.5. The speedup of the adaptive C-N scheme is even worse, which is mainly due to parallel overheads in the full multigrid phase. The latter consists of operations on coarse grids, which are difficult to parallelize efficiently.

These performance differences are due to very dissimilar communication to computation ratios, which we may qualitatively summarize as follows. The arithmetic complexity of WR methods increases linearly with the number of time levels, n,, used for time discretization. The total length of the messages exchanged during the computation is proportional to n,. The number of messages however, and the sequential overhead due to program control are independent of n,. From the high message start-up time on most parallel machines it is clear that the communication time to calculation time ratio of the waveform methods will decrease with increasing function length. The larger the number of time levels, the better the parallel efficiency will be. The corresponding ratio of the C-N method is independent of n,. It equals to first approximation the “worst-case” ratio of the waveform methods, namely the ratio for n, = 2 (i.e., one time step). We refer to [17] for a more quantitative analysis.

In addition we should note that the use of vectorization would make the differences in execution time even larger. On a sixteen processor iPSC/2 and for the above problem sizes an additional vector speedup of a factor of 3 or 4 could be obtained with our implementation of the waveform method, see [17]. However, the problem size is much too small to get any vector speedup with the CrankNicolson method.

S. VandewaIle, R. Piessens / Nonlinear multigrid waveform relaxation 159

5.2. A nonlinear time-periodic problem

We consider the numerical solution of the Brusselator, a well-known nonlinear system of parabolic partial differential equations which models a chemical reaction-diffusion process, [ll]. The two-dimensional Brusselator model is described by the following equations defined over the unit square,

(5 -3)

The functions X( t, x, y) and Y(t, x, y) denote chemical concentrations of intermediate reaction products. A(t, x, y) and B(t, x, y) are the concentrations of the input reagents. The diffusion coefficients are given by D, = 0.004 and D, = 0.008. The chemical reactor length L is chosen as 0.15. We consider Dirichlet boundary conditions

X(t, x, Y> =A, Y(t, x, y) =B/A, (x, y) E aa.

We will determine the stable periodic solution of the Brusselator model in the case of the following homogeneous concentrations of the input reagents,

A(t, x, y) =2+sin(4t) and B(t, x, y) =5.45. (5.4) To this end the problem is discretized with central differences on an equidistant grid with

mesh size h = A. The problem is solved with the periodic multigrid waveform relaxation method (PMWR), using five multigrid levels. Smoothing is performed by block-wise periodic red-black waveform relaxation. The system of two periodic differential equations at each grid point is solved with one Newton linearization step. The linear problem is then solved by discretization with the trapezoidal rule and by direct solution of the resulting cyclic block-bidiagonal matrix. The time interval of interest [0, $r], i.e., one period, is discretized with constant time-increment At = &IL

The convergence is illustrated in Fig. 4, in which the Euclidean norm of the residual (with respect to the set of discrete equations) is plotted for consecutive approximations. The results

a: WR V(I,l) b: WR V(2.1) c: WR W(1,l) d: WR W(2,l)

C 1 I I I

b a , , I , , ,

0 1 2 3 , ,

4 ,

5 6 7 ,

8 9 10 11 12 13 14

iteration number 5

Fig. 4. PMWR convergence for solving (5.3) with (5.4), L = 0.15, h = l/32, At = ~/200.

1

‘(

1:


6-

5-

4-

3-

Oa 0:5 i 1:5 2 215 3 315 4 4:5

Fig. 5. Periodic solution obtained by PMWR (time domain and state space).

Table 1 Averaged convergence factor, periodic Brusselator problem

PMWR V&l) PMWR V(2,l) PMWR W(l, 1) PMWR W(2,l) PMWR F(1, 1) PMWR F(2,l)

0.163 0.163 0.063 0.041 0.064 0.041

indicated by solid lines are obtained by starting the iterative algorithm with a constant time profile in each grid point. The dashed lines represent the results obtained with the full multigrid procedure. In addition we tabulate in Table 1 the experimentally observed convergence factors, averaged over 10 iterations. These convergence factors are of the same order as the convergence factors obtained with the multigrid waveform relaxation method for solving initial-boundary- value problems. The computational complexity of solving a time-periodic parabolic partial differential equation with the periodic multigrid waveform algorithm is therefore similar to the complexity of solving only one initial-boundary-value problem. This corresponds to the theoretical results obtained for linear problems in [18].

The profile of the stable periodic solution at the mid-point of the chemical reactor is plotted in Fig. 5.

The figure to the left shows the time behaviour, whereas the figure to the right illustrates the limit cycle in the state diagram. This solution can also be determined by a straightforward time integration of the equations starting from some initial condition. The resulting process is illustrated in Fig. 6. The periodic solution is obtained when the transient part of the solution has

Fig. 6.

6-

5-

4-

3-

2-

l- Y(t.%,%)

0 I o:5 1 I:5 1 2:5 3 3:5 4 4:5

Periodic solution obtained by brute-force

5 Y(t,%.%)

-f

3

4:%

2

,, 1

X(t.%,%,

0 -

3 4 5 b

time integration (time domain and state space).


decayed. Integration over a long time interval may be required when the periodic solution is not strongly attracting. Furthermore it is difficult to detect in an automatic way that the periodic behaviour is attained. Hence this brute-force time-integration technique is not competitive with the new waveform method.

References

111

PI

]31 141

[51

161 ]71

I81

[91 ]I01

[III [I21

t131

[I41

1151

[I61

[I71

P. Bastian, J. Burmeister and G. Horton, Implementation of a parallel multigrid method for parabolic partial differential equations, in: W. Hackbusch, ed., Parallel Algorithms jor PDEs (Vieweg, Braunschweig, 1991) 18-27. L. Bomans and D. Roose, Benchmarking the iPSC/2 hypercube multiprocessor, Concurrency: Practice and

Experience 1 (1) (1989) 3-18. A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comp. 31 (1977) 333-390. W. Hackbusch, Multigrid Methods and Applications, Springer Series in Computational Mathematics 4 (Springer,

Berlin, 1985). F. Juang, Accuracy increase in waveform relaxation, Report No. UIUCDCS-R-88-1466, Department of Com- puter Science, University of Illinois at Urbana-Champaign, Urbana, IL (1988). C. Lubich and A. Ostermann, Multi-grid dynamic iteration for parabolic equations, BIT 27 (1987) 216-234. U. Miekkala and 0. Nevanlinna, Convergence of dynamic iteration methods for initial value problems, SIAM J.

Sci. Statist. Comput. 8 (4) (1987) 459-482. P. Odent, Electrical-level simulation of VLSI MOS circuits using multi-processor systems, Ph.D. Thesis, Katholieke Universiteit Leuven, Belgium (1990). M. Ozisik, Heat Conduction (Wiley, New York, 1980). L. Peterson and S. Mattisson, Partitioning tradeoffs for waveform relaxation in transient analysis circuit simulation, in: D. Walker and Q. Stout, eds. Proceedings of the Fifth Distributed Memory Computing Conference

(IEEE Computer Society, Washington, 1990) 612-621. I. Prigogine and R. Lefever, Symmetry breaking instabilities, .Z. Chem. Phys. 48 (4) (1968) 1695-1700.

J. Saltz and V. Naik, Towards developing robust algorithms for solving partial differential equations on MIMD machines, Parallel Comput. 6 (1988) 19-44. A. Skjellum, Concurrent dynamic simulation: multicomputer algorithms research applied to ordinary differential-algebraic process systems in chemical engineering, Ph.D. Thesis, California Institute of Technology,

Pasadena, CA (1990). S. Vandewalle and D. Roose, The parallel waveform relaxation multigrid method, in: G. Rodrigue, ed., Parallel Processing for Scientific Computing, Proceedings of the Third SIAM Conference on Parallel Processing for Scientific Computing, Los Angeles, CA, 1987 (SIAM, Philadelphia, PA, 1989) 152-156. S. Vandewalle and R. Piessens, A comparison of the Crank-Nicolson and waveform relaxation multigrid methods on the Intel hypercube, in: J. Mandel, S. McCormick, J. Dendy, C. Farhat, G. Lonsdale, S. Parter, J. Ruge and K. Stliben, eds., Proceedings of the Fourth Copper Mountain Conference on Multigrid Methods (SIAM, Philadelphia, PA, 1990) 417-434.

S. Vandewalle, R. Van Driessche and R. Piessens, The parallel performance of standard parabolic marching schemes, Znternat. J. High Speed Comput. (to appear). S. Vandewalle and R. Piessens, Efficient parallel algorithms for solving initial-boundary value and time-periodic parabolic partial differential equations, Report TW 139, Katholieke Universiteit Leuven, Belgium (1990); SIAM

J. Sci. Statist. Comput. (submitted). [18] S. Vandewalle and R. Piessens, On dynamic iteration methods for solving time-periodic differential equations,

Report TW 148, Katholieke Universiteit Leuven, Belgium (1991); SIAM J. Numer. Anal. (submitted). [19] J. White, F. Odeh, A.S. Sangiovanni-Vincentelli and A. Ruehli, Waveform relaxation: theory and practice,

Memorandum No. UCB/ERL M85/65, Electronics Research Laboratory, College of Engineering, University of California, Berkeley, CA (1985).

[20] D. Womble, A time-stepping algorithm for parallel computers, SIAM J. Sci. Statist. Comput. 11 (5) (1990) 8244837.

Date post:	14-May-2023
Category:	Documents
Upload:	kuleuven
View:	1 times
Download:	0 times

Numerical experiments with nonlinear multigrid waveform relaxation on a parallel processor

Documents