An experimental study of interface relaxation methods for composite ...
Available online at www.sciencedirect.com
Applied Mathematical Modelling 32 (2008) 1620–1641 www.elsevier.com/locate/apm
An experimental study of interface relaxation methods for composite elliptic differential equations P. Tsompanopoulou *, E. Vavalis University of Thessaly, Department of Computer and Communication Engineering, 37 Glavani, GR 382 21 Volos, Greece Center for Research and Technology – Thessaly, Technology Park of Thessaly, GR 385 00 Volos, Greece Received 1 January 2007; received in revised form 1 April 2007; accepted 29 April 2007 Available online 6 July 2007
Abstract The recently proposed simulation framework of interface relaxation for developing multi-domain multi-physics simulation engines is considered. An experimental study of the behavior of two representative interface relaxation methods is presented. Three linear and one non-linear elliptic two-dimensional PDE problems are considered and they are coupled with both cartesian and general decompositions. The characteristics and the effectiveness of the proposed collaborative PDE solving framework in general, and of the two interface relaxation methods in particular are shown. 2007 Elsevier Inc. All rights reserved. PACS: 15A15; 15A09; 15A23 Keywords: Domain decomposition methods; Differential equations; Interface relaxation methods
1. Introduction The simulation and modeling of complex physical systems usually involve many components of different natures and often require parallel computing strategies. Conventional modeling strategies that are typically based on single physical descriptions over a domain seem to be unnatural to some extent and ineffective in many respects existing simulation software applies only to simpler geometrical shapes and physical situations. New approaches for modeling and simulating complex physical systems in a more natural way and, at the same time, in a more convenient manner have been proposed during the past two decades. There are mainly three such approaches for complex problems that involve differential operators. They are based on (1) the domain decomposition (substructuring), (2) the Schwarz splitting and (3) the interface relaxation techniques. The interface relaxation methods do not consider a multi-domain and multi-physics (potentially different PDE operators are applied on different subdomains) simulation engine as a coherent all. They rather treat it as a loosely coupled system of subproblems consisting of much simpler (both as far as the geometry and the * Corresponding author. Present address: University of Thessaly, Department of Computer and Communication Engineering, 37 Glavani, GR 382 21 Volos, Greece. Tel.: +30 2421074976; fax: +30 2421074997. E-mail address: [email protected] (P. Tsompanopoulou).
0307-904X/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.04.014
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1621
differential operator concerns) PDE problems. Assume that given some boundary conditions we either have the analytic solution or we can ‘‘easily’’ compute an approximate solution of these simpler subproblems. The interface relaxation is an iterative method that uses a library of exact or ‘‘single, simple-domain, simple operator’’ PDE solvers and a collection of smoothing procedures to solve a composite PDE problem as follows: (1) (2) (3) (4)
Guess solution values (and derivatives if needed) on all subdomain interfaces. Solve all single PDEs on all the subdomains independently using these estimated values on the interfaces. Compare and improve the values on the interfaces using a relaxer (see below). Return to Step 2 until satisfactory accuracy is achieved.
Collaborating PDE solvers based on interface relaxation have been already proposed [1–5] and have been rather extensively considered for air pollution [6,7], underwater acoustics [8] and gas turbine engine [9,10] related simulations. A relatively large class of interface relaxation method have been derived [11] and theoretically analyzed (e.g., [1,5,12,13]) for model problems and decompositions. Numerical experiments have been also already presented (see for example [11,12]) but most of them are mainly for one-dimensional problems or simple two-dimensional ones but with one-dimensional decompositions. For this paper, we have selected from the existing class of interface relaxation methods two typical representatives that combine most of their characteristics. Namely the ROB and the AVE methods have been considered. The first one smooths the estimation of both the solution and its normal derivative on the interfaces via a single per iteration step based on interface conditions of Robin type, while the latter relaxes the solution and its normal derivatives separately via weighted average schemes in two distinct steps. The rate of the convergence of both methods is independent of the numerical scheme and the discretization details depending only on the physical characteristics of the PDE problem and the high definition of the interfaces. For the AVE and ROB methods we have performed an extensive set of numerical experiments for general two-dimensional PDEs with general two-dimensional decompositions. The objective of this study is to present selected performance data that contribute to the understanding of the characteristics and the idiosyncracies of interface relaxation, further confirm certain existing theoretical results and elucidate some applicability and extensibility issues. The rest of this paper is organized as follows. Section 2 briefly presents the interface relaxation methods considered in this study and their implementation on a distributed environment using an Agents system. In Section 3 we define the PDE problems we considered and in Sections 4–6 we present the numerical results from our experimentation. Our concluding remarks can be found in Section 7. 2. Interface relaxation methods and the SciAgents framework For the two interface relaxation methods in this study we consider the boundary value problem p [ Lu ¼ f ; x 2 X ¼ Xi
ð1Þ
i¼1
subject to boundary conditions on oX which, for simplicity, are taken to be homogeneous Dirichlet. We denote the restrictions of u, L and f in Xi by ui, Li and fi respectively and so we have the set of local PDEs Li ui ¼ fi ; x 2 Xi i ¼ 1; . . . ; p: ð2Þ We assume that the Xi’s do not overlap and that certain conditions are imposed on the interfaces. These conditions may range from simple continuity conditions of the solution and/or its derivatives to general ones (e.g., conservation of momentum) implicitly given as follows: ! \ oui ouj ; uj ; ; J 1 ; J 2 ¼ 0 on Ci;j Xi Xj ; ð3Þ Gi;j ui ; ogi;j ogj;i where g denotes the outward normal to the interface derivative and where J1 and J2 are the jump quantities of the uis and its derivatives. Gi,j might be a function mapping on the interface or even a functional. For simplicity in the presentation of the two interface relaxation methods considered in this study we assume no particular interface conditions (besides the ones naturally imposed by the required continuity),
1622
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
and that X is decomposed in an one-dimensional fashion into the p non-overlapping subdomains Xi [xi1, xi], i = 1, . . . , p with xi1 < xi 2 X for i = 1, . . . , p. 2.1. The ROB method The ROB scheme is defined for the above model problem by the following algorithm: 1. Define:
duiþ1 dx ðkÞ
gii ¼ giiþ1
¼
x¼xi
ðkÞ dui
dx
ðkÞ þ ki uiþ1
x¼xi
þ x¼xi
ðkÞ ki ui x¼xi
9 > > > = > > > ;
i ¼ 1; . . . ; p 1:
ð0Þ
2. Choose initial guesses ui for the solutions on each subdomain Xi, i = 1, 2, . . . , p. ðkÞ 3. Compute the sequence of subdomain solutions ui ; k ¼ 1; 2; . . . by solving the following PDE problems: Lp upðkþ1Þ ¼ fp in Xp ðkþ1Þ ðkþ1Þ dup ðkþ1Þ u1 ¼0 dx þ k u ¼ gpp1 p1 p x¼x0 x¼x p1 x¼x p1 ðkþ1Þ du1 ðkþ1Þ 1 ðkþ1Þ þ k u ¼ g 1 1 1 up ¼0 dx x¼x1 x¼x1 x¼xp 9 ðkþ1Þ > Li ui ¼ fi in Xi > > > > ðkþ1Þ > dui ðkþ1Þ i dx þ ki1 ui ¼ gi1 = x¼x ; i ¼ 2; . . . ; p 1: i1 x¼xi1 > > > ðkþ1Þ > dui ðkþ1Þ > > þ ki ui ¼ gii ; dx ðkþ1Þ
L1 u1
¼ f1
x¼xi
in X1
x¼xi
This scheme, first proposed in [1], is based on a simple relaxation technique that involves the Robin interface conditions shown above. The DE problem is solved in each subdomain where the boundary conditions are provided from the previously computed solution and its normal derivative from the adjacent subdomains. The relaxation parameter ki controls the influence of the function value and/or its normal derivative on the smoothing Robin interface conditions. 2.2. The two-step average AVE method The AVE [11,12] method is a two-step iterative scheme described by the following algorithm: ð0Þ
1. Choose initial guesses ui for the solution on each subdomain Xi, i = 1, 2, . . . , p. ð2kþ1Þ 2. Compute the odd terms of the sequence of subdomain solutions ui by solving the following PDE problems: ð2kÞ ð2kÞ du du gii ¼ bi i þ ð1 bi Þ iþ1 ; i ¼ 1; . . . ; p 1; dx dx x¼xi x¼xi for i ¼ 2; . . . ; p 1 ð2kþ1Þ ð2kþ1Þ ¼ fp in Xp L1 u1 ¼ f1 in X1 L uð2kþ1Þ ¼ f in X Lp up i i i i ð2kþ1Þ ð2kþ1Þ dup duð2kþ1Þ u1 ¼0 ¼ gp1 dx i p1 ¼ gi1 x¼x0 i1 dx x¼xp1 x¼xi1 ð2kþ1Þ du1 ð2kþ1Þ ð2kþ1Þ ¼ g11 up dui ¼0 dx i ¼ gi x¼x1 x¼xp dx x¼xi
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 ð2kÞ
3. Compute the even terms of the sequence of subdomain solution ui problems: ð2kþ1Þ ð2kþ1Þ hii ¼ ai ui þ ð1 ai Þuiþ1 ; i ¼ 1; . . . ; p 1; x¼xi x¼xi for i ¼ 2; . . . ; p 1 ð2kþ2Þ Lp upð2kþ2Þ ¼ fp in Xp ¼ f1 in X1 ð2kþ2Þ L 1 u1 L u ¼ f in X i i i ð2kþ2Þ p1 uð2kþ2Þ i u1 ¼0 ¼ hp1 ð2kþ2Þ i1 p x¼x0 x¼x u ¼ h i p1 i1 x¼x ð2kþ2Þ ð2kþ2Þ up ð2kþ2Þ i1 i ¼ h11 u1 j x¼xp ¼ 0 u x¼x1 j ¼h i
x¼xi
1623
by solving the following PDE
i
The relaxation parameters ai and bi are to smooth the function and its normal derivative respectively, and they both take values in (0, 1). In the first step (odd terms), the Dirichlet problem is solved for each subdomain. The boundary values are computed as a convex combination of the previously computed solutions on adjacent domains. Then a convex combination of the normal derivatives of the previously computed solutions in each subdomain are used to smooth the derivatives on each interface. Using these estimates of the normal derivatives, the Neumann problem is solved in the second step (even terms) for all subdomains.
2.3. The SciAgents experimental framework We have implemented (mainly using C and Java) a whole class of Interface Relaxation methods in an agent-based framework for general two-dimensional decompositions of linear and non-linear elliptic PDE problems. Details about this implementation, which is known as SciAgents can be found in [8,14]. The SciAgents exploit the inherent parallelism in the interface relaxation methods using the Agents computing paradigm over a network of heterogeneous workstations. Specifically, SciAgents transform the physical problem (for example, the one associated with domain XIV in Fig. 2) into a network of local PDE solvers and interface relaxers (marked as solvers and mediators respectively in Fig. 1). This network can then easily mapped onto a heterogeneous distributed computing platform while the coordination is taking place according to generic iterative workflows like the ones described above. In all experiments we have calculated an initial guess for the solution by using the appropriate interpolant on each interface segment. For initial guesses of the normal derivatives we simply impose the correct sign (direction), by setting them equal to the unit outward normal vector. To calculate the required by the mediators derivatives on the interfaces the associated build-in to ELLPACK procedure was used as it is described in [14]. All experiments presented in this paper were run in parallel using single precision arithmetic on heterogeneous SUN workstations connected through an ethernet line. Each subdomain was assigned to a different machine and all the interfaces were assigned to an additional machine.
Fig. 1. The network of solvers and mediators for the composite problem associated with domain XIV shown in Fig. 2.
1624
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
3. The PDE problems under investigation For our experimental analysis we have considered four domains. They are depicted in Fig. 2 (two on its left figure and the other two on its right) together with the associated boundary conditions, the decompositions of each one of the four domains into subdomains and the numbering of these subdomains. Specifically we have considered [XI]: The domain presented on the left figure, consisting of the eight rectangular subdomains ðXIi ; i ¼ 0; . . . ; 7Þ made distinct by dotted lines. II [X ]: The domain presented on the left figure, consisting of the four subdomains ðXII i ; i ¼ 0; . . . ; 3Þ made distinct by solid lines. [XIII]: The domain presented on the right figure, consisting of the four rectangular subdomains ðXIII i ; i ¼ 0; . . . ; 3Þ made distinct by dotted lines. IV [X ]: The domain presented on the right figure, consisting of the four subdomains ðXIV i ; i ¼ 0; . . . ; 3Þ made distinct by solid lines. Next we define four composite Partial Differential Equations by defining, for each one of them, the differential operators associated with its subdomain as follows: [PDE1] Linear Helmholtz equations with constant coefficients (on 5 subdomains) and variable coefficients (on 3 subdomains). Du þ c2i u ¼ fi c20 ¼
1 expðx þ yÞ; 2
on XIi ;
i ¼ 0; . . . ; 7 with
c21 ¼ 10;
c24 ¼ 2ðsinððx þ yÞpÞ þ 4Þ;
c22 ¼ 16;
c25 ¼ 15;
c23 ¼ 20;
c26 ¼ 25;
ð4Þ
x þ y c27 ¼ 5 exp : 8
[PDE2] Linear Helmholtz equations with constant coefficients (on 2 subdomains) and variable coefficients – slightly more complicated than in PDE1 (on 2 subdomains). i ¼ 0; 1; 2; 3 with Du þ c2i u ¼ fi on XII i ; 1 c0 ¼ expðx þ yÞ; c1 ¼ 16; c2 ¼ 25 and ð5Þ 2 x þ y 5 c3 ¼ sinððx þ yÞpÞ þ exp þ 4: 2 8 [PDE3] Linear Helmholtz equation with constant coefficient (on 2 subdomains) and variable coefficients (on one subdomain) and a general elliptic operator with constant coefficients on one subdomain. Du þ 0:2u þ 60ðx2 þ y 2 þ 2Þ ¼ 0;
on XIII 0 ;
III Du þ 0:4u ¼ 0; on XIII 1 and X3 ; ou ou þ Du 10 þ 0:3u ¼ 0 on XIII 2 : ox oy
ð6Þ
[PDE4] Linear Helmholtz equation with constant coefficient (on 2 subdomains) and general non-linear differential operators with variable coefficients (on two subdomains). u Du þ 0:2u 1 þ þ 60ðx2 þ y 2 þ 2Þ ¼ 0; on XIV 0 ; 1000 IV IV Du þ 0:4u ¼ 0; on X1 and X3 ; ð7Þ 2 2 ou u ou 1 ou ou þ3 þ 0:3u ¼ 0; on XIV þ 1þ þ 2 : ox2 1000 oy 2 500 ox oy
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1625
Fig. 2. Two domains, XI and XII, on the left and two XIII and XIV on the right, with the associated boundary conditions, their decompositions and the numbering of these subdomains and their interface segments.
The right hand sides (fi’s) in both PDE1 and PDE2 have been selected so that the true solution is x2 + y2. Therefore, for these problems we can calculate the errors involved in any numerical approximation of their solutions. For each one of these four decompositions we denote the interface segments by I ds1 ;s2 , where d denotes the corresponding decomposition, and s1, s2 the indices of the two subdomains which are adjacent to this particular interface. For example, by I I2;3 , we denote the interface that lays between the subdomains XI2 and XI3 of the XI (see Fig. 2). The above given PDE problems coupled with the associated boundary conditions shown in Fig. 2 lead as to a population of composite elliptic PDE problems. This set of problems is small enough for easy experimentation but at the same time it encapsulates many important features that might affect the performance of simulations based on interface relaxation techniques. Specifically, we let us investigate the affect of the following physical parameters of the PDE problems:
Table 1 The theoretically determined ‘‘optimal’’ values of the interface relaxation parameters used to obtain the data associated with the dotteddashed lines in Figs. 3–6 Interface segment I I0;1 I I1;2 I I2;3 I I0;4 I I1;5 I I2;6 I I4;5 I I5;6 I I4;7 I I5;7
ROB relaxation parameter
AVE relaxation parameters
ki
ai
bi
3.162
0.626
0.374
4.000
0.442
0.559
4.472
0.472
0.528
3.164
0.627
0.375
3.873
0.508
0.492
5.000
0.444
0.556
3.873
0.721
0.179
5.000
0.437
0.564
4.734
0.398
0.600
4.734
0.448
0.552
1626
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
0
log(||u(k+1) – u (k)||1)
0 I Interface I 0,1
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
0
2
4
6
8
10
12
14
16
18
20
log(||u(k+1) – u (k)||1)
Interface
–1
I I2,3
–2
–3
–3
–4
–4
–5
–5
2
4
6
8
10
12
14
16
18
20
8
10
12
14
16
18
20
Interface
I 0,4
16
18
20
16
18
20
16
18
20
16
18
20
I
0
2
4
6
8
10
12
14
1
0
log(||u(k+1) – u (k)||1)
6
–6 0
1
log(||u(k+1) – u (k)||1)
4
–1
–2
0 I
Interface I1,5
–1 –2
–2 –3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
–6
1
1
0
0
–1
Interface I
I
Interface I 2,6
–1
–3
–6
2
I 1,2
0
0
–6
0
I
1
1
I 4,5
0
2
4
6
8
10
–2
–2
–3
–3
–4
–4
–5
–5
12
14
Interface I I5,6
–1
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
1
1
0
0
log(||u(k+1) – u (k)||1)
–6
Interface
Interface I I4,7
–1
Interface I I5,7
–1
–2
–2
–3
–3
–4
–4
–5
–5 –6
–6 0
2
4
6
8
10 12 iterations
14
16
18
20
0
2
4
6
8
10 iterations
12
14
Fig. 3. Reduction of the L1 norm of the difference of two successive iterants on each interface using the ROB scheme for PDE1. Solid lines, dotted lines and circles denote data using ki = 1, ki = 2 and ki = 4 for i = 0, . . . , 8 respectively. The dash-dotted lines denote data using the theoretically determined optimum values for ki’s shown in Table 1.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 0
0
–1
–1
log(||u(k) –u||∞/||u||∞)
Subdomain ΩI0
–2
–3
–4
–4
–5
–5
–6
–6
log(||u(k) –u||∞/||u||∞)
0
2
4
6
8
10
12
14
16
18
20 0
–1
–1
I
Subdomain Ω2
–2
–3
–3
–4
–4
–5
–5
4
6
8
10
12
14
16
18
20
18
20
18
20
18
20
I
Subdomain Ω3
–6 0
2
4
6
8
10
12
14
16
18
20
0
0
0
–1
–1
I
Subdomain Ω4
–2
log(||u(k) –u||∞/||u||∞)
2
–2
–6
log(||u(k) –u||∞/||u||∞)
0
0
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
–6
0
0
–1
–1
Subdomain ΩI
–2
2
4
6
8
10
12
–3
–3
–4
–4
–5
–5
16
I
0
2
4
6
8
10
12
14
16
Subdomain ΩI
–2
6
14
Subdomain Ω5
–2
–3
–6
Subdomain ΩI1
–2
–3
1627
7
–6
–6 0
2
4
6
8
10
iterations
12
14
16
18
20
0
2
4
6
8
10
12
14
16
iterations
Fig. 4. Reduction of the L1 norm of the relative errors in each subdomain using the ROB scheme for PDE1 (legends as in Fig. 3).
1628
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
0
log(||u(k+1) – u(k)||1)
0
Interface
–1
II0,1
–2
–2
–3
–3
–4
–4
–5
–5
–6
–6
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
18
20
18
20
18
20
18
20
18
20
1
1
0
0
Interface
I I2,3
Interface I I
–1
log(||u(k+1) – u(k)||1)
–1 –2
–2
–3
–3
–4
–4
–5
–5
–6
0,4
–6
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
1
1
0
0
Interface
log(||u(k+1) – u(k)||1)
–1
I I1,5
I
Interface I 2,6
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
–6
0
2
4
6
8
10
12
14
16
18
20
1
0
2
4
6
8
10
12
14
16
1
0
0
–1
log(||u(k+1) – u(k)||1)
Interface
I I4,5
I
Interface I 5,6
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
–6
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
1
1
0
0
log(||u(k+1) – u(k)||1)
Interface I I1,2
–1
Interface I
–1
I 4,7
Interface I I
–1
–2
–2
–3
–3
–4
–4
–5
–5
5,7
–6
–6
0
2
4
6
8
10 iterations
12
14
16
18
20
0
2
4
6
8
10 iterations
12
14
16
Fig. 5. Reduction of the L1 norm of the difference of two successive iterants on each interface using the AVE scheme for PDE1. Solid lines, dotted lines and circles denote data using ai = bi = 0.5, ai = bi = 0.25, and ai = bi = 0.75, for i = 0, . . . , 7 respectively. The dash-dotted lines denote data using the theoretically determined optimum values for the ai’s and bi’s shown in Table 1.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 0
–1
–1
log(||u(k) – u||∞/||u||∞)
0
Subdomain ΩI
–2
Subdomain Ω I 1
–2
0
–3
–3
–4
–4
–5
–5
–6
1629
–6 0
2
4
6
8
10
12
14
16
18
0
20
0
0
–1
–1
log(||u(k) – u||∞/||u||∞)
Subdomain ΩI
–2
2
4
6
8
10
–3
–3
–4
–4
–5
–5
–6
14
16
18
20
18
20
18
20
18
20
Subdomain ΩI
–2
2
12
3
–6 0
2
4
6
8
10
12
14
16
18
20
0 0
–1
–1
log(||u(k) – u||∞/||u||∞)
0
I 4
SubdomainΩ
–2
2
4
6
8
10
–3
–4
–4
–5
–5
14
16
I 5
Subdomain Ω
–2
–3
12
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0 0
–1
–1
log(||u(k) – u||∞/||u||∞)
0
I
Subdomain Ω6
–2
2
4
6
8
10
–3
–4
–4
–5
–5
14
16
I
Subdomain Ω7
–2
–3
12
–6
–6 0
2
4
6
8
10
iterations
12
14
16
18
20
0
2
4
6
8
10
12
14
16
iterations
Fig. 6. Reduction of the L1 norm of the relative errors in each subdomain using the AVE scheme for PDE1 (legends as in Fig. 5).
1630
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
0
log(||u(k+1)– u (k)||1)
0 I Interface I 0,1
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
log(||u(k+1)– u (k)||1)
2
4
6
8
10
12
14
16
18
20
16
18
20
18
20
18
20
18
20
0
0 I Interface I 2,3
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
Interface I
I 0,4
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
1
1
0
0
log(||u(k+1)– u (k)||1)
0 1
1
I Interface I 1,5
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
Interface I
I 2,6
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
1
1
0
0
log(||u(k+1)– u (k)||1)
I 1,2
–6
–6
I InterfaceI 4,5
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
Interface I
I 5,6
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
1
1
0
0
log(||u(k+1)– u (k)||1)
Interface I
Interface I
–1
I 4,7
Interface I
–1
–2
–2
–3
–3
–4
–4
–5
–5
I 5,7
–6
–6 0
2
4
6
8
10 12 iterations
14
16
18
20
0
2
4
6
8
10 12 iterations
14
16
Fig. 7. Reduction of the L1 norm of the difference of two successive iterants on each interface using the ROB scheme with the optimum values for the ki’s for PDE1. Solid lines, dotted lines and dash-dotted denote data using the 5-point-star finite difference scheme, the collocation and the finite element respectively.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1631
• The size of the domains (as moving from domain XII to XIV). • The type of the decomposition, by comparing results for cartesian (rectangular) decompositions associated with rectangular domains XI, XIII to the ones for non-cartesian decompositions associated with non-rectangular domains XII, XIV. • The difficulty of the differential operator, as moving from linear Helmholtz the general non-linear operators (considered in PDE4). Furthermore, we note that one might consider the decompositions of subdomains XI and XIII as ‘‘cartesian approximations’’ to the ones of subdomains XII and XIV respectively and as such they might be used to calculate ‘‘good’’ relaxation parameters for the non-cartesian decompositions (see [11–13]). Finally, it is worth to mention that the above problems can be found in the heart of nuclear reactor simulation systems.
4. A linear PDE with a cartesian decomposition: PDE1 For the experiments in this section we discretize each of the subdomains using rectangular meshes with common discretization parameter h = 0.1. Note that for this discretization the local to each subdomain grids match on the interfaces. The local PDE operators were discretized using the ELLPACK’s [15] 5-point-star module, which is an O(h2) finite difference scheme, on all subdomains. This resulted into linear systems with N = 441, 441, 861, 441, 231, 231, 451, 2051 algebraic equations and unknowns on subdomains 0, 1, . . . , 7 respectively. All these systems were solved using the ELLPACK’s Gauss Elimination module for banded matrices. To calculate appropriate values for the relaxation parameters we tried to utilize our one-dimensional theoretical results presented in [12]. For that we collapse appropriate rows or columns of rectangles by considering them as lines and their interface lines as interface points. For example to obtain the relaxation parameters on the interface segments I I0;4 ; I I4;7 and I I1;5 ; I I5;7 we consider the one-dimensional decomposition [0, 2], [2, 3] and [3, 8], and on the segments I I0;1 ; I I1;2 and I I2;3 the [0, 2], [2, 4], [4, 8] and [8, 10] while for the segments I I4;5 ; I I5;6 the [0, 2], [2, 4] and [4, 8] and finally for the segment I I2;6 we consider the decomposition [0, 2], [2, 3]. Furthermore, the ci’s involved in the theoretical expressions of the one-dimensional analysis of [12] were set equal to the average of the coefficient of ui in each subdomain. The new one-dimensional PDE problems were first scaled to [0, 1] by multiplying the ci’s with the square of the length of the interval and then the optimum parameters are computed using the formulas obtained in [12]. The values estimated for PDE1 by using the above procedure are given in Table 1. To examine the basic convergence properties of the methods we denote by u(k) the computed at the kth iteration solution and we present a series of plots concerning its accuracy on each subdomain and on each interface segment by using the following two measures: S i ¼ log kuðkÞ jCi uðk1Þ jCi k1 ; ðkÞ
Ei ¼ log ku jXi ui k1 ;
i ¼ 1; 2; . . . ; number of interfaces;
ð8Þ
i ¼ 1; 2; . . . ; number of subdomains:
ð9Þ
Table 2 The theoretically determined ‘‘optimal’’ values of the interface relaxation parameters used to obtain the data associated with the dotteddashed lines in Figs. 8 and 9 Interface segment I II 0;1 I II 0;2 I II 0;3 I II 1;2 I II 2;3
ROB relaxation parameter
AVE parameters relaxation
ki
ai
bi
5.2248
0.5695
0.4305
4.2360
0.5557
0.4443
4.4721
0.5419
0.4581
5.0000
0.4721
0.5279
5.0003
0.4444
0.5556
1632
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
log(||u(k+1) – u (k)||1)
0
0
Interface
–1
III0,1
–2
–2
–3
–3
–4
–4
–5
–5
–6
0,2
–6 0
2
4
6
8
10
12
14
16
18
20
1
0
2
4
6
8
10
12
14
16
18
20
16
18
20
18
20
18
20
1
0
log(||u(k+1) – u (k)||1)
Interface I II
–1
0
–1
Interface
I II 0,3
Interface I II
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
1,2
–6 0
2
4
6
8
10 12 iterations 1
14
16
18
20
0
2
4
6
8
10 12 iterations
14
log(||u(k+1) – u (k)||1)
0 II
Interface I 2,3
–1 –2 –3 –4 –5 –6 0
2
4
6
8
0
10 iterations 0
14
16
18
20
–1
–1
log(||u(k) – u ||∞/||u||∞)
12
Subdomain
–2
II Ω0
II
Subdomain Ω 1
–2
–3
–3
–4
–4
–5
–5
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0 0
–1
–1
log(||u(k) – u ||∞/||u||∞)
0
II
Subdomain Ω 2
–2
2
4
6
8
–3
–4
–4
–5
–5
–6
12
14
Subdomain
–2
–3
10
16
II
Ω3
–6 0
2
4
6
8
10 12 iterations
14
16
18
20
0
2
4
6
8
10 12 iterations
14
16
Fig. 8. Reduction of the L1 norm of the difference of two successive iterants on each interface and the L1 norm of the relative errors in each subdomain using the ROB scheme for PDE2. Solid lines, dotted lines and circles denote data using ki = 1, ki = 2 and ki = 4 for i = 0, . . . , 4 respectively. The dash-dotted lines denote data using the theoretically determined optimum values for ki’s shown in Table 2.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
0
log(||u(k+1) – u(k)||1)
0 –1
Interface
II I0,1
–2
–3
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
1 0
log(||u(k+1) – u(k)||1)
–6
Interface I 0,2
0
2
4
6
8
10
12
14
16
18
20
16
18
20
18
20
18
20
1 0
–1
Interface
I II 0,3
–2
–3
–3
–4
–4
–5
–5
0
2
4
6
8
10 12 iterations 1
14
Interface I II
–1
–2
–6
II
–1
–2
–6
1633
16
18
20
–6
1,2
0
2
4
6
8
10 iterations
12
14
log(||u(k+1) – u(k)||1)
0
Interface I II
–1
2,3
–2 –3 –4 –5
–6
0
2
4
6
8
10
12
14
16
18
20
log(||u(k)–u||∞/||u||∞)
iterations
0
0
–1
–1
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
–6
0
0
–1
–1
Subdomain Ω
–2
II
Subdomain Ω 1
–2
–3
–6
log(||u(k)–u||∞/||u||∞)
II
Subdomain Ω 0
–2
II 2
0
2
4
6
8
10
–3
–4
–4
–5
–5
14
16
Subdomain Ω II 3
–2
–3
12
–6
–6 0
2
4
6
8
10 12 iterations
14
16
18
20
0
2
4
6
8
10 iterations
12
14
16
Fig. 9. Reduction of the L1 norm of the difference of two successive iterants on each interface and of the L1 norm of the relative errors in each subdomain using the AVE scheme for PDE2. The solid lines, denote data using ai = 0.5, for i = 0, . . . , 4 and the dash-dotted lines denote data using the theoretically determined optimum values for the ai’s and bi’s as these are shown in Table 2.
1634
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
Both Figs. 3 and 4 clearly show that the estimated, by the above procedure, values of the relaxation parameters for the ROB scheme are really close to ‘‘optimum’’. These values result in convergence that is significant faster than the one associated with ki = 2 or 4. In the AVE case though as seen from Figs. 5 and 6 the theoretical determined relaxation parameters do not show such optimality. Nevertheless, it is seen that the rate of convergence, depicted by the slopes of the convergence lines, does not differ significantly in the case of the experimentally ‘‘best’’ relaxation parameters (found by systematic search). This suggests that there exists only a ‘‘somehow’’ start-up cost difference on the convergence rates. By comparing Figs. 4 and 6 one sees that ‘‘optimum’’ ROB seems to be slightly faster than ‘‘optimum’’ AVE while both schemes seem to be rather effective, achieving 3 significant digits in less that 4 iterations and 5 digits in about 10 iterations. We conclude this section by investigating the effect of the particular PDE discretization scheme one might use to solve the individual PDE subproblems. For this we plot in Fig. 7 the history of convergence for the ROB scheme applied to PDE1 using three different discretization modules from ELLPACK. Specifically, we have performed three experiments, one (depicted by dash-dotted lines in the graphs), using the finite element module on all subdomains, another using the 5-point-star finite difference module (solid lines) and another using the collocation module (dotted lines). Note that the finite element and the finite difference schemes are of order O(h2) while the collocation is of order O(h4). For all methods we used h = 0.1. These lead to linear systems of significant difference on the size and of rather different mathematical properties. For example, the collocation matrix has 1764, 1764, 3444, 1764, 924, 924, 1804, 8364 equations and unknowns, the finite element has 304, 306, 601, 306, 167, 167, 319, 1463 and the finite difference has 399, 418, 818, 399, 218, 227, 450, 2000 in subproblems i = 0, . . . , 7 respectively. It is worth to mention that our experimental data for the AVE method for PDE1 not presented here exhibit the same behavior as the one found in Fig. 7. It is therefore safe to claim that the convergence behavior of the Interface Relaxation methods is almost identical to all three cases indicating the natural expectation one might have drawn from the formulation and the preliminary analysis already developed (see for example [11,12]) that the rate of convergence of interface relaxation methods does not depend on local discretization parameters and choices.
5. Linear PDEs with general decompositions: PDE2 and PDE3 We now move into PDE2 which was solved using the same configuration as the one described for PDE1 above with only the following two differences. First we discretize the domain XII using triangular elements with h = 0.1 on all subdomains. On this we used the ELLPACK’s Finite Element discretization module which is also an O(h2) scheme. Secondly to estimate the ‘‘optimal’’ values for the relaxation parameters we use the ‘‘approximate’’ cartesian decomposition associated with XI by utilizing our one-dimensional theoretical results in the following way: We first map the interface segments of the general decomposition that are not straight lines parallel to x- or y-axis to their ‘‘closest’’ interface lines on the ‘‘approximate’’ decomposition. This map can be done using geometry or computational geometry tools and procedures coupled with appropriate objective funcTable 3 The values of the relative error for PDE2 in the L1 and L2 norms of the computed solutions in each subdomain by using the ELLPACK’s Finite Element module on each single subdomain (in the second and third columns) and by using the ROB interface relaxation method with different relaxation parameters (in the subsequent columns) Single domain
ROB relaxation k=1
XII 0 XII 1 XII 2 XII 3
k=2
k = opt
L1
L2
L1
L2
L1
L2
L1
L2
.155E2
.747E3
.424E2
.163E2
.379E02
.163E2
.342E2
.166E2
.845E4
.732E4
.145E3
.721E4
.128E3
.704E4
.101E3
.679E4
.119E3
.11E3
.227E3
.125E3
.214E3
.123E3
.199E3
.122E3
.116E3
.109E3
.219E3
.112E3
.190E3
.108E3
.187E3
.109E3
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1 0
0 –1
log(||u(k+1) – u (k)||∞)
Interface
III 0,1
–2
–3
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
1
–6
0,2
0
2
4
6
8
10
12
14
16
18
20
16
18
20
18
20
18
20
1 0
log(||u(k+1) – u (k)||∞)
0 –1
Interface I
II 0,3
Interface I II
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
Interface I II
–1
–2
–6
1635
1,2
–6 0
2
4
6
8
10 12 iterations 1
14
16
18
20
0
2
4
6
8
10 12 iterations
14
log(||u(k+1) – u (k)||∞)
0 –1
Interface I
II 2,3
–2 –3 –4 –5 –6
0
2
4
6
8
0
10 12 iterations 0
log(||u(k)–u ||∞/||u||∞)
–1
Subdomain Ω
II 0
18
20
Subdomain Ω
–2
–3
–3
–4
–4
–5
–5
II 1
–6
–6 0
log(||u(k)– u ||∞/||u||∞)
16
–1
–2
2
4
6
8
10
12
14
16
18
20
0
0
0
–1
–1
Subdomain Ω II
–2
–3
–4
–4
–5
–5
0
2
4
6
8
10 12 iterations
14
16
2
4
6
8
18
20
–6
10
12
14
16
Subdomain Ω II3
–2
2
–3
–6
14
0
2
4
6
8
10 12 iterations
14
16
Fig. 10. Reduction of the L1 norm of the difference of two successive iterants on each interface and of the L1 norm of the relative errors in each subdomain using the ROB scheme and the Finite Element module on all subdomains for PDE2 with ki = 4 for i = 0, . . . , 4. Solid lines denote data using as space discretization parameter h = 0.1, the dotted lines denote data using h = 0.2.
1636
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
8
8
7
7
6
6
5
5 y axis
y axis
tions. Here we have done it using a straightforward naı¨ve way. Specifically, we use the corresponding PDE operators on each subdomain and for interface I II 0;1 the one-dimensional partitioning [0, 2][2, 10], for interfaces II II I II and I the decomposition [0, 2][2, 8], for interface I II 0;3 2;3 1;2 the decomposition [0, 2][2, 3], while for interface I 0;2 II II we just took the average of the values from interfaces I 0;1 and I 0;3 . The values determined by the above procedure are shown in Table 2. Similarly as for PDE1 we present in Figs. 8 and 9 the convergence histories of the norms of the successive iterants and the relative errors for PDE2 during the first 20 iterations for both the ROB and AVE schemes and different values for the relaxation parameters ki, ai and bi. We first note that by comparing Figs. 3–6 with 8 and 9 one sees that both the rate of convergence and the effectiveness of the parameters of the Interface Relaxation methods are very similar for both the PDE1 and PDE2 regardless of the fact that these two problems differ on many parameters (e.g., different decomposition, different type and number of subdomains). We also note that different local PDE discretizations were also used. Furthermore, the common horizontal leveling of the lines in Figs. 8 and 9 represent the PDE discretization error which, in contrast to the PDE1 case, is due to the geometry of the non-rectangular subdomains for PDE2. Both methods converge very fast achieving the discretization error level in about 5 iterations. To examine the effect of the ROB (similar results, not presented here, were obtained for AVE) Interface Relaxation procedure on the accuracy of the computed solution for PDE2 we give in Table 3 the two norms of the relative errors of the computed solutions obtained by two different ways. Specifically, in the second and third column we list the errors obtained by first imposing Dirichlet boundary conditions with exact (correct) right hand sides on all interfaces and then solving the individual uncoupled PDE subproblems defined on each subdomain by using the Finite Element discretization with h = 0.1 on each one of them. In the subsequent columns we give the computed values of the same errors obtained by the experiment with which we obtained
4 3
3
2
2
1
1 0
0 0
1
2
3
4
5 x axis
6
7
8
9
8
8
7
7
6
6
5
5 y axis
y axis
4
4
3
2
2
1
1
0
1
2
3
4
5
x axis
6
7
8
9
1
2
3
4
0
1
2
3
4
5 x axis
6
7
8
9
5
6
7
8
9
4
3
0
0
0
x axis
Fig. 11. Contour plots of the first 4 iterants (u(k), k = 1, 2, 3, 4) of PDE2 computed by the ROB scheme using the theoretically determined optimum relaxation parameters.
log(||u(k+1) – u(k)||1)
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 2
2
1
1
0
0
–1
–1
–2
–2
–3
Interface
–4 0
log(||u(k+1) – u(k)||1)
1637
5
IV I0,1
10
15
20
25
30
35
40
2
1
1
0
0
–1
–1
–2
–2 Interface I IV
Interface I 0,2
–4 0
2
–3
IV
–3
–4
10
15
20
25
30
35
40
15
20 iterations
25
30
35
40
Interface I IV
–3
0,3
5
1,2
–4 0
5
10
15
20 iterations
25
30
35
40
0
5
10
2
log(||u(k+1) – u(k)||1)
1 0
–1 –2 –3 –4
Interface
0
5
I IV
2,3
10
15
20 iterations
25
30
35
40
Fig. 12. Reduction of the L1 norm of the difference of two successive iterants on each interface using the ROB scheme for PDE3 (legends as in Fig. 10).
the data in Fig. 8. As is seen the relaxation slightly reduces the accuracy indicating that the constant involved in the convergence rate of ROB is significantly small.
5
4
y
3
2
1
0 0
1
2
3 x
4
5
6
IV IV Fig. 13. Space discretization for the PDE4 using cartesian grids on subdomains XIV 0 and X2 and triangular elements on subdomains X1 and XIV with discretization respectively with h = 0.1 in both cases. 3
1638
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
To investigate the effect of the space discretization parameter h on the rate of convergence we present in Fig. 10 the history of the norm of the successive iterants and the norms of the relative errors, for PDE2 using the ROB scheme with ki = 4 for i = 1, . . . , p 1. The dotted lines represent data with h = 0.2 and the solid lines data with h = 0.1. It is seen that even though we approximately double the number of elements this has virtually no effect on the convergence. To conclude the presentation of our experimental data for PDE2 we plot in Fig. 11 the contours of the computed solution after the first four iterations. The drastic smoothing effect on the interfaces is clearly observed. We switch now to PDE3 and present in Fig. 12 the reduction of the norm of the successive iterants of the ROB scheme for the different values of the ki’s. The configuration used to obtain these data was the same as the one for PDE2. The ‘‘optimum’’ values determined from the theoretical analysis were k0 = 0.640536, k1 = 0.6510, k2 = 0.66155, k3 = 0.921287 and k4 = 1.08033. As is clearly seen these values lead to the fastest convergence. In addition, by comparing the slopes of the lines in Figs. 10 and 12, we see that the rate of convergence for the two problems PDE2 and PDE3 that differ on both the size of the subdomains and the PDE operators applied on each one of them, is approximately the same. Finally, we also see the effectiveness of the scheme, in the sense that it offers full accuracy in less than 20 iterations. 6. A non-linear PDE with general decomposition: PDE4 IV For the experiments in this section we have discretized subdomains XIV 0 and X2 using rectangular elements, IV IV and subdomains X1 and X3 using triangular elements as it is shown in Fig. 13. In all four subdomain discretizations we have used h = 0.1. It is worth to observe that the local discretization points do not necessarily match on the interfaces. Therefore, agents performing local interpolation were used (see [14] for details) before relaxing the interface values.
1000
800
800
600
600
T
1200
1000
T
1200
400
400
200
200 0
0 6
1
5 3
1
5
2
4
0
0 6
1 x
2
4 3
3
2
3
2
4
4
1
0 5
0 5
y
y
x
1200 1200
1000
1000 800
600
600 T
T
800
400
400
200
200
0 6
0 5
1
4
2
3
3
2 0 5 x
1
5
2
4 3
3
2
4
1
0
0 6
1 y
4 0 5
y
x
Fig. 14. Three-dimensional plots of the 1st, 2nd, 3rd and 25th iterants (u(k), k = 1, 2, 3, 25) of PDE4 computed by the ROB scheme using ki = 4 for i = 0, . . . , 4.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1639
We have linearized the PDE operators in (7) using the Newton’s method. The resulting linearized PDEs IV were discretized using the Ellpack’s 5-point-star module for subdomains XIV 0 and X2 and the finite element IV IV module for subdomains X1 and X3 . This way subdomains 0, 1, 2 and 3 resulted in 400, 570, 400 and 433 linear equations (and unknowns) respectively. These linear systems were solved using the Gauss Elimination module for banded matrices. In Fig. 14 we plot in three dimensions the computed by the ROB scheme solution after iterations 1, 2, 3 and 25 and in Fig. 15 the convergence history in terms of the difference of successive iterants with 4 different set of
2
1
1
0
0
–1
–1
–2
–2
log(||u(k+1) – u(k)||1)
2
–3
Interface IIV
–4 0
5
15
20
25
30
35
40
0 2
1
1
0
0
–1
–1
–2
–2
log(||u(k+1) – u(k)||1)
Interface I0,2
–4 10
2
IV
–3
IV
–3
0,1
–4
10
15
20
25
30
35
40
15
20 iterations
25
30
35
40
IV
–3
Interface I0,3
5
Interface I1,2
–4 0
5
10
15
20 iterations
25
30
35
40
0
5
10
2
log(||u(k+1) – u(k)||1)
1
0
–1
–2
Interface IIV
–3
2,3
–4 0
5
10
15
20 iterations
25
30
35
40
Fig. 15. Reduction of the L1 norm of the difference of two successive iterants on each interface using the ROB scheme for PDE4. Solid lines, dotted lines and circles denote data using ki = 1, ki = 2 and ki = 4 for i = 0, . . . , 3 respectively. The dash-dotted lines denote data using the theoretically determined optimum values for ki’s.
1640
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
values of the relaxation parameters. The ones represented by the dash-dotted lines are for the ‘‘optimum’’ values determined by using the ‘‘approximate’’ cartesian partition (defined by XIII) and the associated procedure used for PDE2 and described in Section 5. These values were k0 = 0.5338, k1 = 0.5977, k2 = 0.6615, k3 = 1.0357 and k4 = 1.0474 and as it is seen in Fig. 15 they lead to divergence. This is explained from the fact that the differential operators in the subdomains are far away from the Helmholtz model problem we have analyzed theoretically. Other than that the scheme seems to converge fast enough, offering two correct decimal digits in about 5 iterations. Additional experiments, not presented here, show that, as for the previous problems, the rate convergence of both schemes for PDE4 does not significantly depend on either the local discretization parameter h nor the local discretization method used locally on each subdomain. 7. Synopsis We hope that the presented study numerically further establishes intuition about the characteristics and dynamics of multi-physics, multi-domain simulation systems that are based on interface relaxation methods and it adds to the further work needed before these methods put into practical considerations. Specifically, it is our believe that the numerical experiments presented above give us enough confidence to claim that (1) Both methods converge rather rapidly achieving numerical convergence in about 10 iterations. (2) Regardless their totally different formulation (e.g., one step vs. two steps, one parameter vs. two parameters) and motivation and smoothing techniques (via Robin vs. Dirichlet–Neumann conditions) the convergence behavior of both methods is quite similar. (3) Although careful tuning interface relaxation parameters certainly leads to faster convergence, the rate of convergence does not seem to strongly depend on the values of these parameters. Therefore, just rough estimations of the ‘‘optimum’’ values are needed in many cases. Furthermore, techniques for utilizing theoretical results derived for simple model decompositions to tune general decompositions have been proposed and proved their success. (4) The much needed independence on the local discretization details is further confirmed and the concept of using different discretizations, possible custom tailored, for different subproblems is clearly (re-)proved. (5) Moving from linear to non-linear differential operators seems ‘‘easy’’ exhibiting possibly similar convergence properties. (6) Software reuse is pushed to its limit. Our whole SciAgents implementation consists of more than 1.5 Million lines of code, the great majority of which is consisting of legacy Ellpack modules ‘‘wearing agent jackets’’.
References [1] P.L. Lions, On the Schwarz alternating method III: a variant for nonoverlapping subdomains, in: R. Glowinski, G.H. Golub, G.A. Meurant, J. Periaux (Eds.), Domain Decomposition Methods for Partial Differential Equations, SIAM, 1990, pp. 202–223. [2] H.S. McFaddin, J.R. Rice, Collaborating PDE solvers, Appl. Numer. Math. 10 (1992) 279–295. [3] M. Mu, J.R. Rice, Modeling with collaborating PDE solvers – theory and practice, Comput. Syst. Eng. 6 (1995) 87–95. [4] J.R. Rice, An agent-based architecture for solving partial differential equations, SIAM News 31 (6) (1998). [5] M. Mu, Solving composite problems with interface relaxation, SIAM J. Sci. Comput. 20 (1999) 1394–1416. [6] E. Vavalis, A collaborating framework for air pollution simulations, in: NATO-ASI Series, 1999, pp. 349–358. [7] E. Vavalis, Runtime support for collaborative air pollution agents, Syst. Anal. Modell. Simul. 42 (2002) 1575–1600. [8] J. Boloni, D.C. Marinescu, J.R. Rice, P. Tsompanopoulou, E. Vavalis, Agent based scientific simulation and modeling, Concurr.: Pract. Exp. 12 (2000) 845–861. [9] S. Markus, E. Houstis, A. Catlin, J. Rice, P. Tsompanopoulou, E. Vavalis, D. Gottfried, K. Su, G. Balakrishnan, An agent-based netcentric framework for multidisciplinary problem solving environments, Int. J. Comput. Eng. Sci. 1 (2000) 33–60. [10] E.N. Houstis, A.C. Catlin, P. Tsompanopoulou, D. Gottfried, G. Balakrishnan, K. Su, J.R. Rice, Gasturbnlab: a multidisciplinary problem solving environment for gas turbine engine design on a network of non-homogeneous machines, J. Comput. Appl. Math. 149 (1) (2002) 83–100.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1641
[11] J.R. Rice, P. Tsompanopoulou, E. Vavalis, Interface relaxation methods for elliptic differential equations, Appl. Numer. Methods 32 (2000) 219–245. [12] J.R. Rice, P. Tsompanopoulou, E. Vavalis, Fine tuning interface relaxation methods for elliptic differential equations, Appl. Numer. Methods 43 (2002) 459–481. [13] P. Tsompanopoulou, E. Vavalis, Analysis of a geometry based interface relaxation method for elliptic differential equations, submitted for publication. [14] P. Tsompanopoulou, Collaborative PDEs: theory and practice, Ph.D. thesis, Mathematics Department, University of Crete, Greece, 2000. [15] J.R. Rice, R.F. Boisvert, Solving Elliptic Problems Using ELLPACK, Springer-Verlag, New York, NY, 1985.
Applied Mathematical Modelling 32 (2008) 1620–1641 www.elsevier.com/locate/apm
An experimental study of interface relaxation methods for composite elliptic differential equations P. Tsompanopoulou *, E. Vavalis University of Thessaly, Department of Computer and Communication Engineering, 37 Glavani, GR 382 21 Volos, Greece Center for Research and Technology – Thessaly, Technology Park of Thessaly, GR 385 00 Volos, Greece Received 1 January 2007; received in revised form 1 April 2007; accepted 29 April 2007 Available online 6 July 2007
Abstract The recently proposed simulation framework of interface relaxation for developing multi-domain multi-physics simulation engines is considered. An experimental study of the behavior of two representative interface relaxation methods is presented. Three linear and one non-linear elliptic two-dimensional PDE problems are considered and they are coupled with both cartesian and general decompositions. The characteristics and the effectiveness of the proposed collaborative PDE solving framework in general, and of the two interface relaxation methods in particular are shown. 2007 Elsevier Inc. All rights reserved. PACS: 15A15; 15A09; 15A23 Keywords: Domain decomposition methods; Differential equations; Interface relaxation methods
1. Introduction The simulation and modeling of complex physical systems usually involve many components of different natures and often require parallel computing strategies. Conventional modeling strategies that are typically based on single physical descriptions over a domain seem to be unnatural to some extent and ineffective in many respects existing simulation software applies only to simpler geometrical shapes and physical situations. New approaches for modeling and simulating complex physical systems in a more natural way and, at the same time, in a more convenient manner have been proposed during the past two decades. There are mainly three such approaches for complex problems that involve differential operators. They are based on (1) the domain decomposition (substructuring), (2) the Schwarz splitting and (3) the interface relaxation techniques. The interface relaxation methods do not consider a multi-domain and multi-physics (potentially different PDE operators are applied on different subdomains) simulation engine as a coherent all. They rather treat it as a loosely coupled system of subproblems consisting of much simpler (both as far as the geometry and the * Corresponding author. Present address: University of Thessaly, Department of Computer and Communication Engineering, 37 Glavani, GR 382 21 Volos, Greece. Tel.: +30 2421074976; fax: +30 2421074997. E-mail address: [email protected] (P. Tsompanopoulou).
0307-904X/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.04.014
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1621
differential operator concerns) PDE problems. Assume that given some boundary conditions we either have the analytic solution or we can ‘‘easily’’ compute an approximate solution of these simpler subproblems. The interface relaxation is an iterative method that uses a library of exact or ‘‘single, simple-domain, simple operator’’ PDE solvers and a collection of smoothing procedures to solve a composite PDE problem as follows: (1) (2) (3) (4)
Guess solution values (and derivatives if needed) on all subdomain interfaces. Solve all single PDEs on all the subdomains independently using these estimated values on the interfaces. Compare and improve the values on the interfaces using a relaxer (see below). Return to Step 2 until satisfactory accuracy is achieved.
Collaborating PDE solvers based on interface relaxation have been already proposed [1–5] and have been rather extensively considered for air pollution [6,7], underwater acoustics [8] and gas turbine engine [9,10] related simulations. A relatively large class of interface relaxation method have been derived [11] and theoretically analyzed (e.g., [1,5,12,13]) for model problems and decompositions. Numerical experiments have been also already presented (see for example [11,12]) but most of them are mainly for one-dimensional problems or simple two-dimensional ones but with one-dimensional decompositions. For this paper, we have selected from the existing class of interface relaxation methods two typical representatives that combine most of their characteristics. Namely the ROB and the AVE methods have been considered. The first one smooths the estimation of both the solution and its normal derivative on the interfaces via a single per iteration step based on interface conditions of Robin type, while the latter relaxes the solution and its normal derivatives separately via weighted average schemes in two distinct steps. The rate of the convergence of both methods is independent of the numerical scheme and the discretization details depending only on the physical characteristics of the PDE problem and the high definition of the interfaces. For the AVE and ROB methods we have performed an extensive set of numerical experiments for general two-dimensional PDEs with general two-dimensional decompositions. The objective of this study is to present selected performance data that contribute to the understanding of the characteristics and the idiosyncracies of interface relaxation, further confirm certain existing theoretical results and elucidate some applicability and extensibility issues. The rest of this paper is organized as follows. Section 2 briefly presents the interface relaxation methods considered in this study and their implementation on a distributed environment using an Agents system. In Section 3 we define the PDE problems we considered and in Sections 4–6 we present the numerical results from our experimentation. Our concluding remarks can be found in Section 7. 2. Interface relaxation methods and the SciAgents framework For the two interface relaxation methods in this study we consider the boundary value problem p [ Lu ¼ f ; x 2 X ¼ Xi
ð1Þ
i¼1
subject to boundary conditions on oX which, for simplicity, are taken to be homogeneous Dirichlet. We denote the restrictions of u, L and f in Xi by ui, Li and fi respectively and so we have the set of local PDEs Li ui ¼ fi ; x 2 Xi i ¼ 1; . . . ; p: ð2Þ We assume that the Xi’s do not overlap and that certain conditions are imposed on the interfaces. These conditions may range from simple continuity conditions of the solution and/or its derivatives to general ones (e.g., conservation of momentum) implicitly given as follows: ! \ oui ouj ; uj ; ; J 1 ; J 2 ¼ 0 on Ci;j Xi Xj ; ð3Þ Gi;j ui ; ogi;j ogj;i where g denotes the outward normal to the interface derivative and where J1 and J2 are the jump quantities of the uis and its derivatives. Gi,j might be a function mapping on the interface or even a functional. For simplicity in the presentation of the two interface relaxation methods considered in this study we assume no particular interface conditions (besides the ones naturally imposed by the required continuity),
1622
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
and that X is decomposed in an one-dimensional fashion into the p non-overlapping subdomains Xi [xi1, xi], i = 1, . . . , p with xi1 < xi 2 X for i = 1, . . . , p. 2.1. The ROB method The ROB scheme is defined for the above model problem by the following algorithm: 1. Define:
duiþ1 dx ðkÞ
gii ¼ giiþ1
¼
x¼xi
ðkÞ dui
dx
ðkÞ þ ki uiþ1
x¼xi
þ x¼xi
ðkÞ ki ui x¼xi
9 > > > = > > > ;
i ¼ 1; . . . ; p 1:
ð0Þ
2. Choose initial guesses ui for the solutions on each subdomain Xi, i = 1, 2, . . . , p. ðkÞ 3. Compute the sequence of subdomain solutions ui ; k ¼ 1; 2; . . . by solving the following PDE problems: Lp upðkþ1Þ ¼ fp in Xp ðkþ1Þ ðkþ1Þ dup ðkþ1Þ u1 ¼0 dx þ k u ¼ gpp1 p1 p x¼x0 x¼x p1 x¼x p1 ðkþ1Þ du1 ðkþ1Þ 1 ðkþ1Þ þ k u ¼ g 1 1 1 up ¼0 dx x¼x1 x¼x1 x¼xp 9 ðkþ1Þ > Li ui ¼ fi in Xi > > > > ðkþ1Þ > dui ðkþ1Þ i dx þ ki1 ui ¼ gi1 = x¼x ; i ¼ 2; . . . ; p 1: i1 x¼xi1 > > > ðkþ1Þ > dui ðkþ1Þ > > þ ki ui ¼ gii ; dx ðkþ1Þ
L1 u1
¼ f1
x¼xi
in X1
x¼xi
This scheme, first proposed in [1], is based on a simple relaxation technique that involves the Robin interface conditions shown above. The DE problem is solved in each subdomain where the boundary conditions are provided from the previously computed solution and its normal derivative from the adjacent subdomains. The relaxation parameter ki controls the influence of the function value and/or its normal derivative on the smoothing Robin interface conditions. 2.2. The two-step average AVE method The AVE [11,12] method is a two-step iterative scheme described by the following algorithm: ð0Þ
1. Choose initial guesses ui for the solution on each subdomain Xi, i = 1, 2, . . . , p. ð2kþ1Þ 2. Compute the odd terms of the sequence of subdomain solutions ui by solving the following PDE problems: ð2kÞ ð2kÞ du du gii ¼ bi i þ ð1 bi Þ iþ1 ; i ¼ 1; . . . ; p 1; dx dx x¼xi x¼xi for i ¼ 2; . . . ; p 1 ð2kþ1Þ ð2kþ1Þ ¼ fp in Xp L1 u1 ¼ f1 in X1 L uð2kþ1Þ ¼ f in X Lp up i i i i ð2kþ1Þ ð2kþ1Þ dup duð2kþ1Þ u1 ¼0 ¼ gp1 dx i p1 ¼ gi1 x¼x0 i1 dx x¼xp1 x¼xi1 ð2kþ1Þ du1 ð2kþ1Þ ð2kþ1Þ ¼ g11 up dui ¼0 dx i ¼ gi x¼x1 x¼xp dx x¼xi
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 ð2kÞ
3. Compute the even terms of the sequence of subdomain solution ui problems: ð2kþ1Þ ð2kþ1Þ hii ¼ ai ui þ ð1 ai Þuiþ1 ; i ¼ 1; . . . ; p 1; x¼xi x¼xi for i ¼ 2; . . . ; p 1 ð2kþ2Þ Lp upð2kþ2Þ ¼ fp in Xp ¼ f1 in X1 ð2kþ2Þ L 1 u1 L u ¼ f in X i i i ð2kþ2Þ p1 uð2kþ2Þ i u1 ¼0 ¼ hp1 ð2kþ2Þ i1 p x¼x0 x¼x u ¼ h i p1 i1 x¼x ð2kþ2Þ ð2kþ2Þ up ð2kþ2Þ i1 i ¼ h11 u1 j x¼xp ¼ 0 u x¼x1 j ¼h i
x¼xi
1623
by solving the following PDE
i
The relaxation parameters ai and bi are to smooth the function and its normal derivative respectively, and they both take values in (0, 1). In the first step (odd terms), the Dirichlet problem is solved for each subdomain. The boundary values are computed as a convex combination of the previously computed solutions on adjacent domains. Then a convex combination of the normal derivatives of the previously computed solutions in each subdomain are used to smooth the derivatives on each interface. Using these estimates of the normal derivatives, the Neumann problem is solved in the second step (even terms) for all subdomains.
2.3. The SciAgents experimental framework We have implemented (mainly using C and Java) a whole class of Interface Relaxation methods in an agent-based framework for general two-dimensional decompositions of linear and non-linear elliptic PDE problems. Details about this implementation, which is known as SciAgents can be found in [8,14]. The SciAgents exploit the inherent parallelism in the interface relaxation methods using the Agents computing paradigm over a network of heterogeneous workstations. Specifically, SciAgents transform the physical problem (for example, the one associated with domain XIV in Fig. 2) into a network of local PDE solvers and interface relaxers (marked as solvers and mediators respectively in Fig. 1). This network can then easily mapped onto a heterogeneous distributed computing platform while the coordination is taking place according to generic iterative workflows like the ones described above. In all experiments we have calculated an initial guess for the solution by using the appropriate interpolant on each interface segment. For initial guesses of the normal derivatives we simply impose the correct sign (direction), by setting them equal to the unit outward normal vector. To calculate the required by the mediators derivatives on the interfaces the associated build-in to ELLPACK procedure was used as it is described in [14]. All experiments presented in this paper were run in parallel using single precision arithmetic on heterogeneous SUN workstations connected through an ethernet line. Each subdomain was assigned to a different machine and all the interfaces were assigned to an additional machine.
Fig. 1. The network of solvers and mediators for the composite problem associated with domain XIV shown in Fig. 2.
1624
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
3. The PDE problems under investigation For our experimental analysis we have considered four domains. They are depicted in Fig. 2 (two on its left figure and the other two on its right) together with the associated boundary conditions, the decompositions of each one of the four domains into subdomains and the numbering of these subdomains. Specifically we have considered [XI]: The domain presented on the left figure, consisting of the eight rectangular subdomains ðXIi ; i ¼ 0; . . . ; 7Þ made distinct by dotted lines. II [X ]: The domain presented on the left figure, consisting of the four subdomains ðXII i ; i ¼ 0; . . . ; 3Þ made distinct by solid lines. [XIII]: The domain presented on the right figure, consisting of the four rectangular subdomains ðXIII i ; i ¼ 0; . . . ; 3Þ made distinct by dotted lines. IV [X ]: The domain presented on the right figure, consisting of the four subdomains ðXIV i ; i ¼ 0; . . . ; 3Þ made distinct by solid lines. Next we define four composite Partial Differential Equations by defining, for each one of them, the differential operators associated with its subdomain as follows: [PDE1] Linear Helmholtz equations with constant coefficients (on 5 subdomains) and variable coefficients (on 3 subdomains). Du þ c2i u ¼ fi c20 ¼
1 expðx þ yÞ; 2
on XIi ;
i ¼ 0; . . . ; 7 with
c21 ¼ 10;
c24 ¼ 2ðsinððx þ yÞpÞ þ 4Þ;
c22 ¼ 16;
c25 ¼ 15;
c23 ¼ 20;
c26 ¼ 25;
ð4Þ
x þ y c27 ¼ 5 exp : 8
[PDE2] Linear Helmholtz equations with constant coefficients (on 2 subdomains) and variable coefficients – slightly more complicated than in PDE1 (on 2 subdomains). i ¼ 0; 1; 2; 3 with Du þ c2i u ¼ fi on XII i ; 1 c0 ¼ expðx þ yÞ; c1 ¼ 16; c2 ¼ 25 and ð5Þ 2 x þ y 5 c3 ¼ sinððx þ yÞpÞ þ exp þ 4: 2 8 [PDE3] Linear Helmholtz equation with constant coefficient (on 2 subdomains) and variable coefficients (on one subdomain) and a general elliptic operator with constant coefficients on one subdomain. Du þ 0:2u þ 60ðx2 þ y 2 þ 2Þ ¼ 0;
on XIII 0 ;
III Du þ 0:4u ¼ 0; on XIII 1 and X3 ; ou ou þ Du 10 þ 0:3u ¼ 0 on XIII 2 : ox oy
ð6Þ
[PDE4] Linear Helmholtz equation with constant coefficient (on 2 subdomains) and general non-linear differential operators with variable coefficients (on two subdomains). u Du þ 0:2u 1 þ þ 60ðx2 þ y 2 þ 2Þ ¼ 0; on XIV 0 ; 1000 IV IV Du þ 0:4u ¼ 0; on X1 and X3 ; ð7Þ 2 2 ou u ou 1 ou ou þ3 þ 0:3u ¼ 0; on XIV þ 1þ þ 2 : ox2 1000 oy 2 500 ox oy
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1625
Fig. 2. Two domains, XI and XII, on the left and two XIII and XIV on the right, with the associated boundary conditions, their decompositions and the numbering of these subdomains and their interface segments.
The right hand sides (fi’s) in both PDE1 and PDE2 have been selected so that the true solution is x2 + y2. Therefore, for these problems we can calculate the errors involved in any numerical approximation of their solutions. For each one of these four decompositions we denote the interface segments by I ds1 ;s2 , where d denotes the corresponding decomposition, and s1, s2 the indices of the two subdomains which are adjacent to this particular interface. For example, by I I2;3 , we denote the interface that lays between the subdomains XI2 and XI3 of the XI (see Fig. 2). The above given PDE problems coupled with the associated boundary conditions shown in Fig. 2 lead as to a population of composite elliptic PDE problems. This set of problems is small enough for easy experimentation but at the same time it encapsulates many important features that might affect the performance of simulations based on interface relaxation techniques. Specifically, we let us investigate the affect of the following physical parameters of the PDE problems:
Table 1 The theoretically determined ‘‘optimal’’ values of the interface relaxation parameters used to obtain the data associated with the dotteddashed lines in Figs. 3–6 Interface segment I I0;1 I I1;2 I I2;3 I I0;4 I I1;5 I I2;6 I I4;5 I I5;6 I I4;7 I I5;7
ROB relaxation parameter
AVE relaxation parameters
ki
ai
bi
3.162
0.626
0.374
4.000
0.442
0.559
4.472
0.472
0.528
3.164
0.627
0.375
3.873
0.508
0.492
5.000
0.444
0.556
3.873
0.721
0.179
5.000
0.437
0.564
4.734
0.398
0.600
4.734
0.448
0.552
1626
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
0
log(||u(k+1) – u (k)||1)
0 I Interface I 0,1
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
0
2
4
6
8
10
12
14
16
18
20
log(||u(k+1) – u (k)||1)
Interface
–1
I I2,3
–2
–3
–3
–4
–4
–5
–5
2
4
6
8
10
12
14
16
18
20
8
10
12
14
16
18
20
Interface
I 0,4
16
18
20
16
18
20
16
18
20
16
18
20
I
0
2
4
6
8
10
12
14
1
0
log(||u(k+1) – u (k)||1)
6
–6 0
1
log(||u(k+1) – u (k)||1)
4
–1
–2
0 I
Interface I1,5
–1 –2
–2 –3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
–6
1
1
0
0
–1
Interface I
I
Interface I 2,6
–1
–3
–6
2
I 1,2
0
0
–6
0
I
1
1
I 4,5
0
2
4
6
8
10
–2
–2
–3
–3
–4
–4
–5
–5
12
14
Interface I I5,6
–1
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
1
1
0
0
log(||u(k+1) – u (k)||1)
–6
Interface
Interface I I4,7
–1
Interface I I5,7
–1
–2
–2
–3
–3
–4
–4
–5
–5 –6
–6 0
2
4
6
8
10 12 iterations
14
16
18
20
0
2
4
6
8
10 iterations
12
14
Fig. 3. Reduction of the L1 norm of the difference of two successive iterants on each interface using the ROB scheme for PDE1. Solid lines, dotted lines and circles denote data using ki = 1, ki = 2 and ki = 4 for i = 0, . . . , 8 respectively. The dash-dotted lines denote data using the theoretically determined optimum values for ki’s shown in Table 1.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 0
0
–1
–1
log(||u(k) –u||∞/||u||∞)
Subdomain ΩI0
–2
–3
–4
–4
–5
–5
–6
–6
log(||u(k) –u||∞/||u||∞)
0
2
4
6
8
10
12
14
16
18
20 0
–1
–1
I
Subdomain Ω2
–2
–3
–3
–4
–4
–5
–5
4
6
8
10
12
14
16
18
20
18
20
18
20
18
20
I
Subdomain Ω3
–6 0
2
4
6
8
10
12
14
16
18
20
0
0
0
–1
–1
I
Subdomain Ω4
–2
log(||u(k) –u||∞/||u||∞)
2
–2
–6
log(||u(k) –u||∞/||u||∞)
0
0
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
–6
0
0
–1
–1
Subdomain ΩI
–2
2
4
6
8
10
12
–3
–3
–4
–4
–5
–5
16
I
0
2
4
6
8
10
12
14
16
Subdomain ΩI
–2
6
14
Subdomain Ω5
–2
–3
–6
Subdomain ΩI1
–2
–3
1627
7
–6
–6 0
2
4
6
8
10
iterations
12
14
16
18
20
0
2
4
6
8
10
12
14
16
iterations
Fig. 4. Reduction of the L1 norm of the relative errors in each subdomain using the ROB scheme for PDE1 (legends as in Fig. 3).
1628
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
0
log(||u(k+1) – u(k)||1)
0
Interface
–1
II0,1
–2
–2
–3
–3
–4
–4
–5
–5
–6
–6
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
18
20
18
20
18
20
18
20
18
20
1
1
0
0
Interface
I I2,3
Interface I I
–1
log(||u(k+1) – u(k)||1)
–1 –2
–2
–3
–3
–4
–4
–5
–5
–6
0,4
–6
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
1
1
0
0
Interface
log(||u(k+1) – u(k)||1)
–1
I I1,5
I
Interface I 2,6
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
–6
0
2
4
6
8
10
12
14
16
18
20
1
0
2
4
6
8
10
12
14
16
1
0
0
–1
log(||u(k+1) – u(k)||1)
Interface
I I4,5
I
Interface I 5,6
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
–6
0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
1
1
0
0
log(||u(k+1) – u(k)||1)
Interface I I1,2
–1
Interface I
–1
I 4,7
Interface I I
–1
–2
–2
–3
–3
–4
–4
–5
–5
5,7
–6
–6
0
2
4
6
8
10 iterations
12
14
16
18
20
0
2
4
6
8
10 iterations
12
14
16
Fig. 5. Reduction of the L1 norm of the difference of two successive iterants on each interface using the AVE scheme for PDE1. Solid lines, dotted lines and circles denote data using ai = bi = 0.5, ai = bi = 0.25, and ai = bi = 0.75, for i = 0, . . . , 7 respectively. The dash-dotted lines denote data using the theoretically determined optimum values for the ai’s and bi’s shown in Table 1.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 0
–1
–1
log(||u(k) – u||∞/||u||∞)
0
Subdomain ΩI
–2
Subdomain Ω I 1
–2
0
–3
–3
–4
–4
–5
–5
–6
1629
–6 0
2
4
6
8
10
12
14
16
18
0
20
0
0
–1
–1
log(||u(k) – u||∞/||u||∞)
Subdomain ΩI
–2
2
4
6
8
10
–3
–3
–4
–4
–5
–5
–6
14
16
18
20
18
20
18
20
18
20
Subdomain ΩI
–2
2
12
3
–6 0
2
4
6
8
10
12
14
16
18
20
0 0
–1
–1
log(||u(k) – u||∞/||u||∞)
0
I 4
SubdomainΩ
–2
2
4
6
8
10
–3
–4
–4
–5
–5
14
16
I 5
Subdomain Ω
–2
–3
12
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0 0
–1
–1
log(||u(k) – u||∞/||u||∞)
0
I
Subdomain Ω6
–2
2
4
6
8
10
–3
–4
–4
–5
–5
14
16
I
Subdomain Ω7
–2
–3
12
–6
–6 0
2
4
6
8
10
iterations
12
14
16
18
20
0
2
4
6
8
10
12
14
16
iterations
Fig. 6. Reduction of the L1 norm of the relative errors in each subdomain using the AVE scheme for PDE1 (legends as in Fig. 5).
1630
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
0
log(||u(k+1)– u (k)||1)
0 I Interface I 0,1
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
log(||u(k+1)– u (k)||1)
2
4
6
8
10
12
14
16
18
20
16
18
20
18
20
18
20
18
20
0
0 I Interface I 2,3
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
Interface I
I 0,4
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
1
1
0
0
log(||u(k+1)– u (k)||1)
0 1
1
I Interface I 1,5
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
Interface I
I 2,6
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
1
1
0
0
log(||u(k+1)– u (k)||1)
I 1,2
–6
–6
I InterfaceI 4,5
–1
–1
–2
–2
–3
–3
–4
–4
–5
–5
Interface I
I 5,6
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
1
1
0
0
log(||u(k+1)– u (k)||1)
Interface I
Interface I
–1
I 4,7
Interface I
–1
–2
–2
–3
–3
–4
–4
–5
–5
I 5,7
–6
–6 0
2
4
6
8
10 12 iterations
14
16
18
20
0
2
4
6
8
10 12 iterations
14
16
Fig. 7. Reduction of the L1 norm of the difference of two successive iterants on each interface using the ROB scheme with the optimum values for the ki’s for PDE1. Solid lines, dotted lines and dash-dotted denote data using the 5-point-star finite difference scheme, the collocation and the finite element respectively.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1631
• The size of the domains (as moving from domain XII to XIV). • The type of the decomposition, by comparing results for cartesian (rectangular) decompositions associated with rectangular domains XI, XIII to the ones for non-cartesian decompositions associated with non-rectangular domains XII, XIV. • The difficulty of the differential operator, as moving from linear Helmholtz the general non-linear operators (considered in PDE4). Furthermore, we note that one might consider the decompositions of subdomains XI and XIII as ‘‘cartesian approximations’’ to the ones of subdomains XII and XIV respectively and as such they might be used to calculate ‘‘good’’ relaxation parameters for the non-cartesian decompositions (see [11–13]). Finally, it is worth to mention that the above problems can be found in the heart of nuclear reactor simulation systems.
4. A linear PDE with a cartesian decomposition: PDE1 For the experiments in this section we discretize each of the subdomains using rectangular meshes with common discretization parameter h = 0.1. Note that for this discretization the local to each subdomain grids match on the interfaces. The local PDE operators were discretized using the ELLPACK’s [15] 5-point-star module, which is an O(h2) finite difference scheme, on all subdomains. This resulted into linear systems with N = 441, 441, 861, 441, 231, 231, 451, 2051 algebraic equations and unknowns on subdomains 0, 1, . . . , 7 respectively. All these systems were solved using the ELLPACK’s Gauss Elimination module for banded matrices. To calculate appropriate values for the relaxation parameters we tried to utilize our one-dimensional theoretical results presented in [12]. For that we collapse appropriate rows or columns of rectangles by considering them as lines and their interface lines as interface points. For example to obtain the relaxation parameters on the interface segments I I0;4 ; I I4;7 and I I1;5 ; I I5;7 we consider the one-dimensional decomposition [0, 2], [2, 3] and [3, 8], and on the segments I I0;1 ; I I1;2 and I I2;3 the [0, 2], [2, 4], [4, 8] and [8, 10] while for the segments I I4;5 ; I I5;6 the [0, 2], [2, 4] and [4, 8] and finally for the segment I I2;6 we consider the decomposition [0, 2], [2, 3]. Furthermore, the ci’s involved in the theoretical expressions of the one-dimensional analysis of [12] were set equal to the average of the coefficient of ui in each subdomain. The new one-dimensional PDE problems were first scaled to [0, 1] by multiplying the ci’s with the square of the length of the interval and then the optimum parameters are computed using the formulas obtained in [12]. The values estimated for PDE1 by using the above procedure are given in Table 1. To examine the basic convergence properties of the methods we denote by u(k) the computed at the kth iteration solution and we present a series of plots concerning its accuracy on each subdomain and on each interface segment by using the following two measures: S i ¼ log kuðkÞ jCi uðk1Þ jCi k1 ; ðkÞ
Ei ¼ log ku jXi ui k1 ;
i ¼ 1; 2; . . . ; number of interfaces;
ð8Þ
i ¼ 1; 2; . . . ; number of subdomains:
ð9Þ
Table 2 The theoretically determined ‘‘optimal’’ values of the interface relaxation parameters used to obtain the data associated with the dotteddashed lines in Figs. 8 and 9 Interface segment I II 0;1 I II 0;2 I II 0;3 I II 1;2 I II 2;3
ROB relaxation parameter
AVE parameters relaxation
ki
ai
bi
5.2248
0.5695
0.4305
4.2360
0.5557
0.4443
4.4721
0.5419
0.4581
5.0000
0.4721
0.5279
5.0003
0.4444
0.5556
1632
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
log(||u(k+1) – u (k)||1)
0
0
Interface
–1
III0,1
–2
–2
–3
–3
–4
–4
–5
–5
–6
0,2
–6 0
2
4
6
8
10
12
14
16
18
20
1
0
2
4
6
8
10
12
14
16
18
20
16
18
20
18
20
18
20
1
0
log(||u(k+1) – u (k)||1)
Interface I II
–1
0
–1
Interface
I II 0,3
Interface I II
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
1,2
–6 0
2
4
6
8
10 12 iterations 1
14
16
18
20
0
2
4
6
8
10 12 iterations
14
log(||u(k+1) – u (k)||1)
0 II
Interface I 2,3
–1 –2 –3 –4 –5 –6 0
2
4
6
8
0
10 iterations 0
14
16
18
20
–1
–1
log(||u(k) – u ||∞/||u||∞)
12
Subdomain
–2
II Ω0
II
Subdomain Ω 1
–2
–3
–3
–4
–4
–5
–5
–6
–6 0
2
4
6
8
10
12
14
16
18
20
0 0
–1
–1
log(||u(k) – u ||∞/||u||∞)
0
II
Subdomain Ω 2
–2
2
4
6
8
–3
–4
–4
–5
–5
–6
12
14
Subdomain
–2
–3
10
16
II
Ω3
–6 0
2
4
6
8
10 12 iterations
14
16
18
20
0
2
4
6
8
10 12 iterations
14
16
Fig. 8. Reduction of the L1 norm of the difference of two successive iterants on each interface and the L1 norm of the relative errors in each subdomain using the ROB scheme for PDE2. Solid lines, dotted lines and circles denote data using ki = 1, ki = 2 and ki = 4 for i = 0, . . . , 4 respectively. The dash-dotted lines denote data using the theoretically determined optimum values for ki’s shown in Table 2.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1
0
log(||u(k+1) – u(k)||1)
0 –1
Interface
II I0,1
–2
–3
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
1 0
log(||u(k+1) – u(k)||1)
–6
Interface I 0,2
0
2
4
6
8
10
12
14
16
18
20
16
18
20
18
20
18
20
1 0
–1
Interface
I II 0,3
–2
–3
–3
–4
–4
–5
–5
0
2
4
6
8
10 12 iterations 1
14
Interface I II
–1
–2
–6
II
–1
–2
–6
1633
16
18
20
–6
1,2
0
2
4
6
8
10 iterations
12
14
log(||u(k+1) – u(k)||1)
0
Interface I II
–1
2,3
–2 –3 –4 –5
–6
0
2
4
6
8
10
12
14
16
18
20
log(||u(k)–u||∞/||u||∞)
iterations
0
0
–1
–1
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
–6
0
0
–1
–1
Subdomain Ω
–2
II
Subdomain Ω 1
–2
–3
–6
log(||u(k)–u||∞/||u||∞)
II
Subdomain Ω 0
–2
II 2
0
2
4
6
8
10
–3
–4
–4
–5
–5
14
16
Subdomain Ω II 3
–2
–3
12
–6
–6 0
2
4
6
8
10 12 iterations
14
16
18
20
0
2
4
6
8
10 iterations
12
14
16
Fig. 9. Reduction of the L1 norm of the difference of two successive iterants on each interface and of the L1 norm of the relative errors in each subdomain using the AVE scheme for PDE2. The solid lines, denote data using ai = 0.5, for i = 0, . . . , 4 and the dash-dotted lines denote data using the theoretically determined optimum values for the ai’s and bi’s as these are shown in Table 2.
1634
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
Both Figs. 3 and 4 clearly show that the estimated, by the above procedure, values of the relaxation parameters for the ROB scheme are really close to ‘‘optimum’’. These values result in convergence that is significant faster than the one associated with ki = 2 or 4. In the AVE case though as seen from Figs. 5 and 6 the theoretical determined relaxation parameters do not show such optimality. Nevertheless, it is seen that the rate of convergence, depicted by the slopes of the convergence lines, does not differ significantly in the case of the experimentally ‘‘best’’ relaxation parameters (found by systematic search). This suggests that there exists only a ‘‘somehow’’ start-up cost difference on the convergence rates. By comparing Figs. 4 and 6 one sees that ‘‘optimum’’ ROB seems to be slightly faster than ‘‘optimum’’ AVE while both schemes seem to be rather effective, achieving 3 significant digits in less that 4 iterations and 5 digits in about 10 iterations. We conclude this section by investigating the effect of the particular PDE discretization scheme one might use to solve the individual PDE subproblems. For this we plot in Fig. 7 the history of convergence for the ROB scheme applied to PDE1 using three different discretization modules from ELLPACK. Specifically, we have performed three experiments, one (depicted by dash-dotted lines in the graphs), using the finite element module on all subdomains, another using the 5-point-star finite difference module (solid lines) and another using the collocation module (dotted lines). Note that the finite element and the finite difference schemes are of order O(h2) while the collocation is of order O(h4). For all methods we used h = 0.1. These lead to linear systems of significant difference on the size and of rather different mathematical properties. For example, the collocation matrix has 1764, 1764, 3444, 1764, 924, 924, 1804, 8364 equations and unknowns, the finite element has 304, 306, 601, 306, 167, 167, 319, 1463 and the finite difference has 399, 418, 818, 399, 218, 227, 450, 2000 in subproblems i = 0, . . . , 7 respectively. It is worth to mention that our experimental data for the AVE method for PDE1 not presented here exhibit the same behavior as the one found in Fig. 7. It is therefore safe to claim that the convergence behavior of the Interface Relaxation methods is almost identical to all three cases indicating the natural expectation one might have drawn from the formulation and the preliminary analysis already developed (see for example [11,12]) that the rate of convergence of interface relaxation methods does not depend on local discretization parameters and choices.
5. Linear PDEs with general decompositions: PDE2 and PDE3 We now move into PDE2 which was solved using the same configuration as the one described for PDE1 above with only the following two differences. First we discretize the domain XII using triangular elements with h = 0.1 on all subdomains. On this we used the ELLPACK’s Finite Element discretization module which is also an O(h2) scheme. Secondly to estimate the ‘‘optimal’’ values for the relaxation parameters we use the ‘‘approximate’’ cartesian decomposition associated with XI by utilizing our one-dimensional theoretical results in the following way: We first map the interface segments of the general decomposition that are not straight lines parallel to x- or y-axis to their ‘‘closest’’ interface lines on the ‘‘approximate’’ decomposition. This map can be done using geometry or computational geometry tools and procedures coupled with appropriate objective funcTable 3 The values of the relative error for PDE2 in the L1 and L2 norms of the computed solutions in each subdomain by using the ELLPACK’s Finite Element module on each single subdomain (in the second and third columns) and by using the ROB interface relaxation method with different relaxation parameters (in the subsequent columns) Single domain
ROB relaxation k=1
XII 0 XII 1 XII 2 XII 3
k=2
k = opt
L1
L2
L1
L2
L1
L2
L1
L2
.155E2
.747E3
.424E2
.163E2
.379E02
.163E2
.342E2
.166E2
.845E4
.732E4
.145E3
.721E4
.128E3
.704E4
.101E3
.679E4
.119E3
.11E3
.227E3
.125E3
.214E3
.123E3
.199E3
.122E3
.116E3
.109E3
.219E3
.112E3
.190E3
.108E3
.187E3
.109E3
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 1
1 0
0 –1
log(||u(k+1) – u (k)||∞)
Interface
III 0,1
–2
–3
–3
–4
–4
–5
–5
0
2
4
6
8
10
12
14
16
18
20
1
–6
0,2
0
2
4
6
8
10
12
14
16
18
20
16
18
20
18
20
18
20
1 0
log(||u(k+1) – u (k)||∞)
0 –1
Interface I
II 0,3
Interface I II
–1
–2
–2
–3
–3
–4
–4
–5
–5
–6
Interface I II
–1
–2
–6
1635
1,2
–6 0
2
4
6
8
10 12 iterations 1
14
16
18
20
0
2
4
6
8
10 12 iterations
14
log(||u(k+1) – u (k)||∞)
0 –1
Interface I
II 2,3
–2 –3 –4 –5 –6
0
2
4
6
8
0
10 12 iterations 0
log(||u(k)–u ||∞/||u||∞)
–1
Subdomain Ω
II 0
18
20
Subdomain Ω
–2
–3
–3
–4
–4
–5
–5
II 1
–6
–6 0
log(||u(k)– u ||∞/||u||∞)
16
–1
–2
2
4
6
8
10
12
14
16
18
20
0
0
0
–1
–1
Subdomain Ω II
–2
–3
–4
–4
–5
–5
0
2
4
6
8
10 12 iterations
14
16
2
4
6
8
18
20
–6
10
12
14
16
Subdomain Ω II3
–2
2
–3
–6
14
0
2
4
6
8
10 12 iterations
14
16
Fig. 10. Reduction of the L1 norm of the difference of two successive iterants on each interface and of the L1 norm of the relative errors in each subdomain using the ROB scheme and the Finite Element module on all subdomains for PDE2 with ki = 4 for i = 0, . . . , 4. Solid lines denote data using as space discretization parameter h = 0.1, the dotted lines denote data using h = 0.2.
1636
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
8
8
7
7
6
6
5
5 y axis
y axis
tions. Here we have done it using a straightforward naı¨ve way. Specifically, we use the corresponding PDE operators on each subdomain and for interface I II 0;1 the one-dimensional partitioning [0, 2][2, 10], for interfaces II II I II and I the decomposition [0, 2][2, 8], for interface I II 0;3 2;3 1;2 the decomposition [0, 2][2, 3], while for interface I 0;2 II II we just took the average of the values from interfaces I 0;1 and I 0;3 . The values determined by the above procedure are shown in Table 2. Similarly as for PDE1 we present in Figs. 8 and 9 the convergence histories of the norms of the successive iterants and the relative errors for PDE2 during the first 20 iterations for both the ROB and AVE schemes and different values for the relaxation parameters ki, ai and bi. We first note that by comparing Figs. 3–6 with 8 and 9 one sees that both the rate of convergence and the effectiveness of the parameters of the Interface Relaxation methods are very similar for both the PDE1 and PDE2 regardless of the fact that these two problems differ on many parameters (e.g., different decomposition, different type and number of subdomains). We also note that different local PDE discretizations were also used. Furthermore, the common horizontal leveling of the lines in Figs. 8 and 9 represent the PDE discretization error which, in contrast to the PDE1 case, is due to the geometry of the non-rectangular subdomains for PDE2. Both methods converge very fast achieving the discretization error level in about 5 iterations. To examine the effect of the ROB (similar results, not presented here, were obtained for AVE) Interface Relaxation procedure on the accuracy of the computed solution for PDE2 we give in Table 3 the two norms of the relative errors of the computed solutions obtained by two different ways. Specifically, in the second and third column we list the errors obtained by first imposing Dirichlet boundary conditions with exact (correct) right hand sides on all interfaces and then solving the individual uncoupled PDE subproblems defined on each subdomain by using the Finite Element discretization with h = 0.1 on each one of them. In the subsequent columns we give the computed values of the same errors obtained by the experiment with which we obtained
4 3
3
2
2
1
1 0
0 0
1
2
3
4
5 x axis
6
7
8
9
8
8
7
7
6
6
5
5 y axis
y axis
4
4
3
2
2
1
1
0
1
2
3
4
5
x axis
6
7
8
9
1
2
3
4
0
1
2
3
4
5 x axis
6
7
8
9
5
6
7
8
9
4
3
0
0
0
x axis
Fig. 11. Contour plots of the first 4 iterants (u(k), k = 1, 2, 3, 4) of PDE2 computed by the ROB scheme using the theoretically determined optimum relaxation parameters.
log(||u(k+1) – u(k)||1)
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641 2
2
1
1
0
0
–1
–1
–2
–2
–3
Interface
–4 0
log(||u(k+1) – u(k)||1)
1637
5
IV I0,1
10
15
20
25
30
35
40
2
1
1
0
0
–1
–1
–2
–2 Interface I IV
Interface I 0,2
–4 0
2
–3
IV
–3
–4
10
15
20
25
30
35
40
15
20 iterations
25
30
35
40
Interface I IV
–3
0,3
5
1,2
–4 0
5
10
15
20 iterations
25
30
35
40
0
5
10
2
log(||u(k+1) – u(k)||1)
1 0
–1 –2 –3 –4
Interface
0
5
I IV
2,3
10
15
20 iterations
25
30
35
40
Fig. 12. Reduction of the L1 norm of the difference of two successive iterants on each interface using the ROB scheme for PDE3 (legends as in Fig. 10).
the data in Fig. 8. As is seen the relaxation slightly reduces the accuracy indicating that the constant involved in the convergence rate of ROB is significantly small.
5
4
y
3
2
1
0 0
1
2
3 x
4
5
6
IV IV Fig. 13. Space discretization for the PDE4 using cartesian grids on subdomains XIV 0 and X2 and triangular elements on subdomains X1 and XIV with discretization respectively with h = 0.1 in both cases. 3
1638
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
To investigate the effect of the space discretization parameter h on the rate of convergence we present in Fig. 10 the history of the norm of the successive iterants and the norms of the relative errors, for PDE2 using the ROB scheme with ki = 4 for i = 1, . . . , p 1. The dotted lines represent data with h = 0.2 and the solid lines data with h = 0.1. It is seen that even though we approximately double the number of elements this has virtually no effect on the convergence. To conclude the presentation of our experimental data for PDE2 we plot in Fig. 11 the contours of the computed solution after the first four iterations. The drastic smoothing effect on the interfaces is clearly observed. We switch now to PDE3 and present in Fig. 12 the reduction of the norm of the successive iterants of the ROB scheme for the different values of the ki’s. The configuration used to obtain these data was the same as the one for PDE2. The ‘‘optimum’’ values determined from the theoretical analysis were k0 = 0.640536, k1 = 0.6510, k2 = 0.66155, k3 = 0.921287 and k4 = 1.08033. As is clearly seen these values lead to the fastest convergence. In addition, by comparing the slopes of the lines in Figs. 10 and 12, we see that the rate of convergence for the two problems PDE2 and PDE3 that differ on both the size of the subdomains and the PDE operators applied on each one of them, is approximately the same. Finally, we also see the effectiveness of the scheme, in the sense that it offers full accuracy in less than 20 iterations. 6. A non-linear PDE with general decomposition: PDE4 IV For the experiments in this section we have discretized subdomains XIV 0 and X2 using rectangular elements, IV IV and subdomains X1 and X3 using triangular elements as it is shown in Fig. 13. In all four subdomain discretizations we have used h = 0.1. It is worth to observe that the local discretization points do not necessarily match on the interfaces. Therefore, agents performing local interpolation were used (see [14] for details) before relaxing the interface values.
1000
800
800
600
600
T
1200
1000
T
1200
400
400
200
200 0
0 6
1
5 3
1
5
2
4
0
0 6
1 x
2
4 3
3
2
3
2
4
4
1
0 5
0 5
y
y
x
1200 1200
1000
1000 800
600
600 T
T
800
400
400
200
200
0 6
0 5
1
4
2
3
3
2 0 5 x
1
5
2
4 3
3
2
4
1
0
0 6
1 y
4 0 5
y
x
Fig. 14. Three-dimensional plots of the 1st, 2nd, 3rd and 25th iterants (u(k), k = 1, 2, 3, 25) of PDE4 computed by the ROB scheme using ki = 4 for i = 0, . . . , 4.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1639
We have linearized the PDE operators in (7) using the Newton’s method. The resulting linearized PDEs IV were discretized using the Ellpack’s 5-point-star module for subdomains XIV 0 and X2 and the finite element IV IV module for subdomains X1 and X3 . This way subdomains 0, 1, 2 and 3 resulted in 400, 570, 400 and 433 linear equations (and unknowns) respectively. These linear systems were solved using the Gauss Elimination module for banded matrices. In Fig. 14 we plot in three dimensions the computed by the ROB scheme solution after iterations 1, 2, 3 and 25 and in Fig. 15 the convergence history in terms of the difference of successive iterants with 4 different set of
2
1
1
0
0
–1
–1
–2
–2
log(||u(k+1) – u(k)||1)
2
–3
Interface IIV
–4 0
5
15
20
25
30
35
40
0 2
1
1
0
0
–1
–1
–2
–2
log(||u(k+1) – u(k)||1)
Interface I0,2
–4 10
2
IV
–3
IV
–3
0,1
–4
10
15
20
25
30
35
40
15
20 iterations
25
30
35
40
IV
–3
Interface I0,3
5
Interface I1,2
–4 0
5
10
15
20 iterations
25
30
35
40
0
5
10
2
log(||u(k+1) – u(k)||1)
1
0
–1
–2
Interface IIV
–3
2,3
–4 0
5
10
15
20 iterations
25
30
35
40
Fig. 15. Reduction of the L1 norm of the difference of two successive iterants on each interface using the ROB scheme for PDE4. Solid lines, dotted lines and circles denote data using ki = 1, ki = 2 and ki = 4 for i = 0, . . . , 3 respectively. The dash-dotted lines denote data using the theoretically determined optimum values for ki’s.
1640
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
values of the relaxation parameters. The ones represented by the dash-dotted lines are for the ‘‘optimum’’ values determined by using the ‘‘approximate’’ cartesian partition (defined by XIII) and the associated procedure used for PDE2 and described in Section 5. These values were k0 = 0.5338, k1 = 0.5977, k2 = 0.6615, k3 = 1.0357 and k4 = 1.0474 and as it is seen in Fig. 15 they lead to divergence. This is explained from the fact that the differential operators in the subdomains are far away from the Helmholtz model problem we have analyzed theoretically. Other than that the scheme seems to converge fast enough, offering two correct decimal digits in about 5 iterations. Additional experiments, not presented here, show that, as for the previous problems, the rate convergence of both schemes for PDE4 does not significantly depend on either the local discretization parameter h nor the local discretization method used locally on each subdomain. 7. Synopsis We hope that the presented study numerically further establishes intuition about the characteristics and dynamics of multi-physics, multi-domain simulation systems that are based on interface relaxation methods and it adds to the further work needed before these methods put into practical considerations. Specifically, it is our believe that the numerical experiments presented above give us enough confidence to claim that (1) Both methods converge rather rapidly achieving numerical convergence in about 10 iterations. (2) Regardless their totally different formulation (e.g., one step vs. two steps, one parameter vs. two parameters) and motivation and smoothing techniques (via Robin vs. Dirichlet–Neumann conditions) the convergence behavior of both methods is quite similar. (3) Although careful tuning interface relaxation parameters certainly leads to faster convergence, the rate of convergence does not seem to strongly depend on the values of these parameters. Therefore, just rough estimations of the ‘‘optimum’’ values are needed in many cases. Furthermore, techniques for utilizing theoretical results derived for simple model decompositions to tune general decompositions have been proposed and proved their success. (4) The much needed independence on the local discretization details is further confirmed and the concept of using different discretizations, possible custom tailored, for different subproblems is clearly (re-)proved. (5) Moving from linear to non-linear differential operators seems ‘‘easy’’ exhibiting possibly similar convergence properties. (6) Software reuse is pushed to its limit. Our whole SciAgents implementation consists of more than 1.5 Million lines of code, the great majority of which is consisting of legacy Ellpack modules ‘‘wearing agent jackets’’.
References [1] P.L. Lions, On the Schwarz alternating method III: a variant for nonoverlapping subdomains, in: R. Glowinski, G.H. Golub, G.A. Meurant, J. Periaux (Eds.), Domain Decomposition Methods for Partial Differential Equations, SIAM, 1990, pp. 202–223. [2] H.S. McFaddin, J.R. Rice, Collaborating PDE solvers, Appl. Numer. Math. 10 (1992) 279–295. [3] M. Mu, J.R. Rice, Modeling with collaborating PDE solvers – theory and practice, Comput. Syst. Eng. 6 (1995) 87–95. [4] J.R. Rice, An agent-based architecture for solving partial differential equations, SIAM News 31 (6) (1998). [5] M. Mu, Solving composite problems with interface relaxation, SIAM J. Sci. Comput. 20 (1999) 1394–1416. [6] E. Vavalis, A collaborating framework for air pollution simulations, in: NATO-ASI Series, 1999, pp. 349–358. [7] E. Vavalis, Runtime support for collaborative air pollution agents, Syst. Anal. Modell. Simul. 42 (2002) 1575–1600. [8] J. Boloni, D.C. Marinescu, J.R. Rice, P. Tsompanopoulou, E. Vavalis, Agent based scientific simulation and modeling, Concurr.: Pract. Exp. 12 (2000) 845–861. [9] S. Markus, E. Houstis, A. Catlin, J. Rice, P. Tsompanopoulou, E. Vavalis, D. Gottfried, K. Su, G. Balakrishnan, An agent-based netcentric framework for multidisciplinary problem solving environments, Int. J. Comput. Eng. Sci. 1 (2000) 33–60. [10] E.N. Houstis, A.C. Catlin, P. Tsompanopoulou, D. Gottfried, G. Balakrishnan, K. Su, J.R. Rice, Gasturbnlab: a multidisciplinary problem solving environment for gas turbine engine design on a network of non-homogeneous machines, J. Comput. Appl. Math. 149 (1) (2002) 83–100.
P. Tsompanopoulou, E. Vavalis / Applied Mathematical Modelling 32 (2008) 1620–1641
1641
[11] J.R. Rice, P. Tsompanopoulou, E. Vavalis, Interface relaxation methods for elliptic differential equations, Appl. Numer. Methods 32 (2000) 219–245. [12] J.R. Rice, P. Tsompanopoulou, E. Vavalis, Fine tuning interface relaxation methods for elliptic differential equations, Appl. Numer. Methods 43 (2002) 459–481. [13] P. Tsompanopoulou, E. Vavalis, Analysis of a geometry based interface relaxation method for elliptic differential equations, submitted for publication. [14] P. Tsompanopoulou, Collaborative PDEs: theory and practice, Ph.D. thesis, Mathematics Department, University of Crete, Greece, 2000. [15] J.R. Rice, R.F. Boisvert, Solving Elliptic Problems Using ELLPACK, Springer-Verlag, New York, NY, 1985.
