CPS713 Case Study II) CFD and Numerical Relativity NAS Part I -- Initial NAS Benchmark Setup Geoffrey Fox NPAC 111 College Place Syracuse NY 13244-4100 Abstract of CPS713-Case Study II Initial NAS Benchmark Setup This is first of three foilsets on CFD and NAS Benchmarks This describes the four basic NAS benchmarks and their relation to the Navier Stokes Equations in the the 5 component CFD equations We use opportunity to discuss time discretization and stepping in general based on Hirsch CFD book and Numerical Recipes Stability, Implicit and Explicit formulations are introduced Beam Warming Equations and their stability NAS Parallel Benchmarks Here we summarize three foilsets CPS713NAS-I,II,III This Section is Based on two publications The NAS Parallel Benchmarks David Bailey, John Barton, Thomas Lasinski, Horst Simon NASA Ames Sec. 3: A Methodology for Benchmarking some CFD Kernels on Highly Parallel Machines Sisira Weeratunga, Eric Barszcz, Rod Fatoohi, V. Venkatakrishnan Summary of Sec 3: NAS Parallel Benchmarks:The CFD Kernels A Collection of Iterative PDE Solvers embedded in a Pseudo CFD application program is proposed (and now used by community) for the performance evaluation of CFD Codes on highly parallel processors. The solvers are also now used as a general benchmark for different parallel machines. The pseudo application program is stripped of complexities associated with real CFD application programs, thereby enabling a simpler description of the algorithms. However it is capable of reproducing the essential computation and data access/movement characteristics of large scale state of the art CFD codes. These latter features make the CFD kernels suitable for CPS713 as a framework for discussing essential computer science issues in CFD. Overview of Four Computational (PDE Solver) Kernels in NAS Parallel Benchmarks NASA analysed their implicit CFD codes and identified four computational kernels which typified the majority of three-dimensional Euler (non Viscous) and Navier-Stokes (Viscous) CFD programs using finite difference or finite volume methods on a structured grid. First Two Computational Kernels 1)Solution of multiple independent systems of non diagonally dominant block tridiagonal equations with a 5 by 5 block size. (i.e. a matrix made out of 5 by 5 blocks with one nonzero block on diagonal and one nonzero block either side of diagonal) 2)Solution of multiple independent systems of scalar (i.e. not block) pentadiagonal equations. (i.e. matrix with nonzeros on diagonal and two nonzero elements either side of it.) These two kernels are representative of computations associated with implicit operator in the NASA ARC3D code. Second Two Computational Kernels 3) Regular Sparse 5 by 5 block matrix -vector multiplication i.e. Sparse matrix-vector multiplication where each operation is multiplication by a 5 by 5 submatrix into a five dimensional subvector This is typical of explicit part of almost all CFD solvers using structured grids. 4) Regular Sparse 5 by 5 block lower and upper triangular system solution i.e same basic 5 by 5 submatrix into 5 dimensional vector operation with the arithmetic that of using LU decomposition on a sparse matrix. This is typical of implicit solver used in newer codes such as INS3D-LU The First Two Simplifications of the NAS Parallel Benchmarks 1)Absence of Realistic Boundary Conditions: In a real CFD code,one would normally transform from a complicated physical domain into a computational domain where regular (structured) grid is appropriate. This transformation increases comptational complexity and storage requirements of real codes. 2)Simplified artificial dissipation The benchmark adds a simple linear fourth order "artificial dissipation" term. In real codes an adaptive (dependent on pressure changes) mix of nonlinear second and fourth-difference terms are used. This simplification does not change computational structure significantly. Artificial Dissipation is used to improve numerical stability of solver. The benchmark is designed to be very stable with more dissipation than normal. The Second Two Simplifications of the NAS Parallel Benchmarks 3) Time Differencing Schemes There are several clever time differencing schemes which have been developed to improve stability and accuracy of Navier Stokes Equations Upwind Differencing Flux-vector Splitting Total Variation Diminishing -- TVD Flux-difference Splitting The benchmarks only use the combination of a particular implicit differencing with artificial dissipation 4) Absence of Turbulence Models Forming the Reynolds Averaged Equations can introduce additional (typically one or two) differential equations and change values of parameters such as viscosity. This is absent in benchmark. Basic Formalism for 5 Component CFD We can write the five basic variables of CFD as a five component vector given in Hirsch Chapters 16 and 22 In detail the five components are given by: Note this uses Ideal Gas Equation of State P V = R T with V proportional to inverse of r What is relation of NAS benchmarks to Navier Stokes Equations -- Terms E F G? The terms E(U) F(U) and G(U) correspond to the non viscous terms in the Navier Stokes equations Note that the enthalpy H = u(5) + pressure p and this appears in 5th component of E F and G T is temperature u(5) = r times E Total energy E per unit mass = Internal Energy (cv T) plus external Energy (ux2 + uy2 + uz2)/2 Equation of State is for ideal gas pressure p = r . T . constant(r) . What is relation of NAS benchmarks to Navier Stokes Equations -- Terms T V W ? The terms T (U,Ux) V(U,Uh) and W(U,Uz) include A numerical viscosity term proportional to d(1..5)x,h,z which ends up if all the values of d are equal (which they almost are) as a equation like: An approximation to Viscosity term in Navier Stokes Equation where ALL cross derivatives as listed below, are neglected This approximation is implied by equation (3.1a) which doesn't allow such terms as those above. Note also that numerical coefficient (k3 . k4) seems to be missing in 5th component of T U and W The Ñ Temperature term for evolution of r E ( 5th componentof U) What is relation of NAS benchmarks to Navier Stokes Equations-- Term H and Boundary Conditions ? The term H is to viewed as an external force so arranged that a simple function exactly satisfies equations The simple function is a linear combination of powers ( 0 through 4 ) of x h and z. The boundary conditions and initial conditions are similarily set to have a very simple form Note H does not depend on the unknown U in spite of what memo says! H depends only on independent variables x h and z. What is relation of NAS benchmarks to Navier Stokes Equations -- Space and Time Variables Note computational domain is the unit cube in x h and z space We use these greek labels to indicate one has (notionally) transformed from original physical domain x h and z. (As noted the transformation matrix has been ignored) We use t not t for two reasons Firstly we probably have transformed time as well as space in map to computational domain Secondly we only want steady state solution U and so we can view t as labelling iteration and not necessarily as time in a diffusive equation. This implies all that counts is accuracy at infinite time t and not accuracy at intermediate finite time t. Treatment of Time Discretization for PDE's We will follow treatments in Volume 1 of Numerical Computation of Internal and External Flows Fundamentals of Numerical Discretization C. Hirsch , Wiley 1988 And Numerical Recipes by William Press, Saul Teukolsky, William Vetterling, and Brian Flannery Cambridge University Press 1992 Note two key approachs: Explicit: Values at next time explicitly calculable from those at previous time step(s) Implicit: Values at next time given implicitly by matrix equations involving previous time step(s). Treatment of Time Stepping for PDE's : Numerical Recipes -- Chapter 19 Equation (19.1.1) (19.1.6) are flux-conservative equations of general type of NAS benchmarks Equation (19.1.11) is FTCS (Forward Time Centered Space) for (19.1.6) Von Neumann Stability Analysis for constant coefficient equations Lax Method Equation (19.1.14) Courant Stability Condition Equation (19.1.17) Physical Interpretation in Fig 19.1.3 Numerical Viscosity implied in Lax's Method Equation (19.1.19) Comparison of Different Methods for Solving Linear Convection Equation This discusses figures 8.3.6 to 8.3.8 of Hirsch Convection Equation is (19.1.6) in Numerical Recipes Lax-Friedrich's Scheme is called Lax Method - Equations (19.1.14,15) in Numerical Recipes First Order Upwind Scheme is Equation (19.1.27) and Figure 19.1.4 in Numerical Recipes Note effect of numerical viscosity in Lax and Upwind methods which unphysically damps wave propagation in last two figures. The Leapfrog Method is Figure 19.1.5 and Equation 19.1.30 in Numerical Recipes The Lax-Wendroff Method is Equations (19.1.37,38) and Figure 19.1.7 of Numerical Recipes Both Leapfrog and Lax-Wendroff have formal errors which are second order in time step. Other two methods are only first order accurate. The smaller numerical viscosity is good in figures 8.3.7,8 but causes some unphysical oscillations in figure 8.3.6 for both second order methods Overview of Stability and Accuracy of Partial Differential Equation Solvers Note accuracy determines formally how fast the errors tend to zero as the time step d t tends to zero. Stability determines if the errors build up or reinforce each other as one iterates in time All time differencing schemes are valid representations of equations i.e. reduce to differential equation in limit of dx and d t tend to zero However not all differencing schemes are stable! We can investigate stability for simple linear equations with von Neumann's Analysis Let u(n)j be replaced by u(n)j + e(n)j where: n labels time iteration and j spatial discretization If the equations are linear in u, the e(n)j obeys same equations as u(n)j but with simple (zero) boundary conditions If equation not linear, then this analysis still governs local behavior and gives a stability condition must must be satisfied by linearized equations To linearize something like u ¶ (u) / ¶ t, replace varying u outside derivative by its local value. Solution of Homogenous Equations by Fourier Analysis If the Velocity v in convection equation is constant and we have the natural boundary conditions for e then Fourier Analysis diagonalizes difference equations Note that general equation: Solution of Linear Stability Equations for the Linear Convection Equation We expand error e in fashion suggested by Fourier solution General Stability Analysis of some two-step Methods -- Basic Iteration Equation See sec. 11.1 and table 11.1 of Hirsch, Volume 1. This analysis introduced by Beam and Warming (People not Hot Rods) This is implicit unless b = 0 when it becomes explicit This set of equations is defined in Equation (11.1.22) of Hirsch after he discusses a more general set involving K(n-1) Differential Operator and its Eigenvalues for Beam Warming Equations Note we follow same notation as NAS Benchmark Paper Equation (3.11) Hirsch calls our For Stability Analysis write the linear form: K = W U where W is for one dimension and second order differencing a 5 by 5 matrix When analysing the convection equation with K = -v ¶ / ¶x , we get W = -ivk is pure imaginary Stability Conditions for Beam Warming Equations Note that Explicit Euler was always unstable for convection equation but for our case, the artificial viscosity ensures that explicit Euler is stable for small enough time step D t Implicit Euler is always stable for dissipative systems because these by definition have Re(W) negative. The effect of Nonlinearity on Implicit Equations Taking the special case q=0 and b=1 used in the NAS benchmark Equation (3.12)