NASA ContractorR eport 198273 ICASE Report No. 96-4 ICASE SCHWARZ-BASED ALGORITHMS FOR COMPRESSIBLE FLOWS M. D. Tidriri AppDnwed f pub& n"=e 0ftrlbudo Unimfted NASA ContractN o. NASI-19480 January 1996 Institutef or ComputerA pplications in Science and Engineering NASA Langley Research Center Hampton, VA 23681-0001 Operatedb y UniversitiesS pace Research Association 19960401 137 NationalA eronautics and Space Administration Langley Research Center Hampton, Virginia 23681-0001 i4tT r/l ~~OE Schwarz-based algorithms for compressible flows M. D. Tidriri * Institute for Computer Applications in Science and Engineering NASA Langley Research Center, Hampton, VA 23681-0001 Abstract We investigate in this paper the application of Schwarz-based algorithms to com- pressible flows. First, we study the combination of these methods with defect-correction procedures. We then study the effect on the Schwarz-based methods of replacing the explicit treatment of the boundary conditions by an implicit one. In the last part of this paper we study the combination of these methods with Newton-Krylov matrix- free methods. Numerical experiments that show the performance of our approaches are then presented. •ICASE, MS 132C, NASA-LaRC, Hampton, VA 23681-0001 email:[email protected]. This work was supported by the National Aeronautics and Space Administration under NASA contract NASI- 19480 while the author was in residence at the Institute for Computer Applications in Science and Engineering. 1 Introduction To compute steady compressible flows one often uses an implicit discretization ap- proach, which leads to a large sparse linear system that must be solved at each time step. In the derivation of this system one often uses a defect-correction procedure, in which the left-hand side of the system is discretized with a lower order approxima- tion than that used for the right-hand side. This is due to storage considerations and computational complexity, and also to the fact that the resulting lower order matrix is better conditioned than the higher order matrix. The resulting schemes are only moderately implicit. In the case of structured, body-fitted grids, the linear system can easily be solved using approximate factorization (AF), which is among the most widely used methods for such grids. However, for unstructured grids, such techniques are no longer valid, and the system is solved using direct or iterative techniques. Because of the prohibitive computational costs and large memory requirements for the solu- tion of compressible flows, iterative methods are preferred. In these defect-correction methods, which are implemented in most CFD computer codes, the mismatch in the right- and left-hand side operators, together with explicit treatment of the boundary conditions, lead to a severely limited CFL number, which results in a slow convergence to steady state aerodynamic solutions. Many authors have tried to replace explicit boundary conditions with implicit ones (see for instance [25], [21], and [13]). Although they clearly demonstrate that high CFL numbers are possible, the reduction in CPU time is not clear cut. The investigation of defect-correction procedures based on Krylov methods, to- gether with implicit treatment of the boundary conditions has been done by the author in [24]. In [24] the author has also studied Newton-Krylov matrix-free (see also [3], [22], [23], and [10]) methods combined with mixed discretization in the implicitly de- fined Jacobian Preconditioner. The preconditioner based on incomplete factorizations studied in [24] is difficult to parallelize efficiently. The focus in this work is on the developement of algorithms that are suitable for the parallel computing environment. In this case, domain decomposition methods that allow the reduction of the global solution of a given problem to the solutions of local subproblems are preferred. We propose, therefore, to combine these methods with the preconditioned Newton-Krylov matrix-free methods developed in [24]. One of the domain decomposition algorithms that has potential applications on parallel computers is the additive Schwarz algorithm [8]. The other Schwarz-based method; the multiplicative Schwarz method [8] can also be used in the parallel en- vironment by using a multi-coloring process. The proposed algorithm is, therefore, to combine the Newton-Krylov matrix-free methods with the Schwarz-based methods. The combination of Newton-Krylov matrix-free with domain decomposition methods was first introduced by the author in [22] and [23]. More precisely, the author has com- bined the Newton-Krylov matrix-free method with the Domain Decomposition Time 1 Marching Algorithm that was introduced by Le Tallec and Tidriri in [11] (see also [22] and [23]). In the next section, we describe the Euler solver. In section 3, we describe the methodology studied in this paper. In section 4, a comprehensive study of Schwarz- based methods combined with defect-correction procedures with explicit and implicit boundary conditions is performed. We then study the combination of the Schwarz- based methods with the Newton-Krylov matrix-free methods. The last section is de- voted to some conclusions and extensions. 2 Description of the Euler solver 2.1 Governing Equations The bidimensional Euler Equations in conservative form writes Wt + F(W), + G(W)y 0, (1) where W = (p, pu, pv, e)T, F = (pu, pu2 + p, puv, u(e + p))T, and G = (pu, puv, pv2 + p,v(e + p))T. Above p is the density, u, v are the velocity components, e is the internal energy, p is the pressure defined by p = (y- 1)(e- (p(u 2 + v2)/2)), and finally, -yi s a constant with - z 1.4 for air. After changing the variables into the curvilinear coordinate 'r= t, ý = ý(X, Y), 7 = IxY) we obtain the following set of equations W, + (F)ý + (G), = 0, (2) where W and the contravariant flux vectors, F and 0, are defined in terms of the Cartesian fluxes and the Jacobian determinant of the coordinate system transforma- tion, through F =J-W S=J-1 (qtW + r7.r -F- yG), and 2 a(x,y,t) det (_ /Y) = From now on, the tilde in the expressions of WIVF,, and G will be omitted. 2.2 Finite volume scheme An implicit finite volume discretization of equation (2) can be written as (W~~±'-W~% A (cid:127)Ay+ (F ji,.- ½?§.zA +(Gý(cid:127)+1 _ +1 AA =0, (3) where the values are taken at the center of either the cell (i, j) or the interfaces of the cell (i, J) and its neighbours. To compute the fluxes above, we shall use a flux splitting approach, which is defined for F by (see [20]) F = F+ +F-, with similar expressions for G. F+ is associated with the positive eigenvalues whereas F- is associated with the negative ones, and G+, G-are defined analogously. Let JW = - wthen the implicit split-flux discretization of (3) is given by 8W- + A(jF(+(r + F-)-+' + 5.(G+ + G-)-+') = 0, where Jý is defined by F= -[Fi+ 112,j - Fi- 1/2,j] (4) and 8. is defined similarly. This yields the following nonlinear system f(Wn+l) = 0. (5) This nonlinear system will be solved by using the proposed approach of this paper, which is based on a Newton-Krylov method (see next sections). Now, we shall de- scribe the more standard defect-correction method, which is based on the following linearization of first order in time of the nonlinear system above 3 I + +± .A. + J + !B-.)]jWn =- -Ar (J&F'~+ 8,G'). The superscripts i and e above indicate that the implicit and explicit operators are dis- cretized using different schemes. The dots indicate that the difference operators apply to the product of the Jacobian matrices with SWI. The matrices A+, A-, B+, and B- are defined by A+-OF+ A-= F- OW' ~Ow' B+ - OG+ B-- = G- Ow' aw The compact form of the above equation corresponds to the following defect-correction procedure A8Wn = b. (6) The different fluxes above are computed using the Roe's approximate Riemann solver [17]. Three limiters are employed: minmod, Superbee, and Van Leer. The Jacobians are evaluated using first-order Roe's scheme, or the first-order flux-vector split scheme [20], which corresponds to the true partials of the positive and negative flux vectors as described earlier. However, in the context of defect-correction method the flux-vector split scheme has been shown to give improved convergence rates over the Roe matrices. Therefore, for the defect-correction approach the Jacobian matrices corresponding to the flux-vector split scheme are used in the left-hand side. This results in an inconsistent left and right-hand side operators. Remark 2.1 For most CFD codes, the implicit spatial differences are only first-order accurate. The higher-order matrix representation is difficult to obtain, even if it is possible the resulting matrix is very large, requires a lot of storage, large operation count in its evaluation, and may be very difficult to invert. Following this remark, the implicit spatial differences (the left-hand side) in equa- tion (6) are approximated, only, through a first-order accurate scheme. The ex- plicit spatial differences (right-hand side) in equation (6) are approximated using the higher-order formulations of Roe's scheme, that are based on the work of Osher and Chakravarthy [16]. 2.3 Explicit boundary conditions The boundary conditions are derived using the locally one-dimensional characteristic variable boundary conditions, which yields (for the derivations see for example [15]): 4 2.3.1 Farfield-Subsonic Inflow Pb = (1/2)Pa + Pi + sign((A')poco[k.,(ua - uj) + ky(va - vi)] Pb = pN + [(Pb - PF)/Co] ub = ua + k,[(P. - Pb)/(poco)]sign(A') Vb = V, + ky[(P. - Pb)/(poco)]sign(A') Above, the point a is outside the computational domain, point b is on the compu- tational boundary, and i is inside the computational domain. 2.3.2 Farfield-Subsonic Outflow Pb = P. Pb = P. + [(Pb--P a)/C'o] Ub = u. + k.[(Pa - Pb)/(poco)]sign(A')k Vb = Va + ky[(Pa - Pb)/(p~oo)]sign(A')k 2.3.3 Impermeable Surface Pb = P, :F poco Ub = Ur - kx(kxur + kyv') Vb = V, - k,(k(cid:2)u, + kyv,) Where the point r is the center of the first cell from the boundary and the minus sign in equation (2) is used if r is in the positive k direction from the boundary, and the plus sign is used if r is in the negative direction from the boundary. 2.3.4 Farfield-Supersonic Inflow In this case all eigenvalues have the same sign. Since we have an in inflow case all variables are specified. 2.3.5 Farfield-Supersonic Outflow In this case also, all eigenvalues have the same sign. But now we have an outflow case, therefore, all variables must be obtained from the solution in the computational domain. All variables are extrapolated from inside the computational domain to the boundary. 5 2.4 Implicit boundary conditions In the implicit form the above boundary conditions can be written in the form of operators formulated as functions of the conservation vector W: fb(W) = 0 (7) and are implemented implicitly through: orfb 6w = -b(w). CW Using these implicit boundary conditions the author showed in [24], that starting from a small initial CFL number (10), CFL may be adaptively advanced according to: -*f(W)lP-1 = CFL?. CFLn+1 Iif(W)fll(cid:127) where the uperscripts denote the iteration in time. This is the key to the successful implementation of the preconditioned Newton-Krylov matrix-free method studied in [24], and which we combine here with the Schwarz-based methods. 3 Description of the methodology Newton-Krylov methods first proposed by Brown and Saad [3], have been investigated for compressible Euler and Navier-Stokes equations using unstructured grids in [22], [23], and [10], and for structured grids in [4], and [5], and [24]. In [22] and [23], the author has studied both transonic and supersonic compress- ible Navier-Stokes flows. In [4], [5], and [24] a study of a convection-diffusion model problem, the full potential flows and the transonic compressible Euler flows have been performed. implicitly defined Jacobian preconditioner. The most effective preconditioner, ILU, is difficult to parallelize efficiently. On the other hand domain decomposition methods appear to be effective for the parallel solution of large systems of linear or nonlinear algebraic equations resulting from the application of finite element methods or finite difference methods to fluid dynamics problems. The alternating method introduced by H. A. Schwarz in 1890 [19] appears to be the earliest domain decomposition method. For two subdomains this algorithm is intrinsically sequential. Its extension to include the case of many subdomains was done by P. L. Lions [12]. As a consequence of this work, the additive Schwarz methods were developed. Another method, which is a direct generalization of the original alternating method is the multiplicative algorithm. These methods reduce the solution of the global problem on the global domain to the solution of subproblems on local subdomains, obtained by considering an overlapping subdivision of the global domain. 6 Most of the theory and applications of the Schwarz-based methods have been pri- marily performed for elliptic and parabolic boundary value problems discretized using finite element methods. In this paper we shall focus on their applications to the hy- perbolic problems. We shall also study their combination with the Newton-Krylov matrix-free methods studied in [24]. 3.1 Newton's Method Consider the following nonlinear system of equations f(W) =0, (8) where f is a nonlinear function from JR2 to IR2. Newton's method applied to (8) results in the following iteration "* Define u , an initial guess 0 "* For k = 0,1,2, ... until convergence do Solve J(Wk)SWk -f(Wk), (9) Set Wk+l = Wk ±+ Wk, (10) where J(Wk) = -- (Wk) is the sytem Jacobian. For the compressible Euler case (see section 2) this Jacobian corresponds to a higher- order matrix-representation. Using direct-methods to solve the system (9), the memory requirements and the computational complexity are prohibitive. In this case iterative methods are preferred and the system (9) is solved only approximately. The resulting method is called the inexact newton method [6], and corresponds to the following iteration "* Define u , an initial guess 0 "* For k = 0, 1, 2,--- until convergence do Solve J(Wk)SWk = -f(Wk), (11) Set Wk+j = Wk + aJWk, (12) where J(Wk) (Wk) denotes the sytem Jacobian as before, and a is a parameter C9W selected using a line search or trust region method ([3] and [7]). 7 3.2 Krylov methods The iterative methods we will use to solve the linear system (11) which we rewrite as JSw = -f, (13) where f and its Jacobian J are evaluated at the current iterate, are the Krylov method. If w is an initial guess for the true solution of (13), then letting w = wo + Z, we have 0 the equivalent system JZ = r° where r' -f Jwo is the initial residual. Let Km be the Krylov subspace - Km := Span{r0, Jr°,..., Jn-'r°O. Arnoldi's method and GMRES both find an approximate solution w,=wo+Zm, with ZEKm, such that either (-f -Jw) I Km for Arnoldi's method or IIf + JWmH112 = minwEwo±+K f + Jw11 (= minzEKmI[r° - JZ ) 2 2 for GMRES. Here, I1.H[2 denotes the Euclidien norm on JR2 and orthogonality is meant in the usual Euclidien sense. In these Krylov methods only the action of the Jacobian J times a vector w, and not J explicitly is required. In the context of problem (8), this action can be approximated by difference quotient of the form J(u)w f(u +± w)- f(u) where u is the current approximation to a root of (8) and e is a scalar. Selecting an optimal parameter c in the difference formula for approximating J(u)w might be a difficult problem. If e is too small then the rounding errors made in the numerator are amplified by a factor of order which leads to an inaccurate result. If on the other hand 6 is too large then the approximation of J(u)w will be poor. Any reasonable choice of 6 should attempt to reach a compromise between these two difficulties. The technique for choosing the scalar 6 we use here is: _ 11V11'2 . max{j(u,v)j,typujvl}. 8