Adaptive High-Order Methods for Fluid FlowsAuthors: Lott & Deane

Applied Mathematics and Scientific Computation

The goal of this research is to construct an Adaptive Spectral Element Method to solve non-linear PDE’s.

We aim to create an efficient method for solving large scale PDE's on complex domains. We use adaptivity to maintain a high working accuracy for long time periods without wasting computational resources. Spectral Element Methods are computationally efficient for higher dimensional problems due to their tensor product formulation of the system operators.

The Spectral Element Method approximates the solution to the PDE by partitioning the domain into elements, and on each element, solving an integral equation. In order to solve the equation, one constructs a weighted sum of polynomial basis functions. In particular, SEM uses the Gauss-Lobatto-Legendre quadrature rules in order to obtain an exact solution for integrands of order 2P+1 or less, when using Legendre polynomials of degree P. This gives SEM an exponential convergence property that methods such as Finite Element, or Finite Difference methods don’t have. As we increase the polynomial degree in the SEM, we exponentially approach the solution to the integral equation. In low order methods, where only the number of elements is increased, one obtains only algebraic convergence to the solution.

Solution to 1D Viscous Burgers’ Equation using P-Adaptivity

We have implemented error estimators that determine when to increase or decrease the polynomial degree on each element in order to obtain an accurate solution without wasting computational resources in well behaved areas of the solution.

To the right, we show the computational work required to integrate a linear advection equation 5 periods while maintaining a phase shift error <.1 (Karniadakis).

Gauss Legendre Lobatto Nodes

4 Elements Polynomial Degree 5

Slope of solution at x=0, within 99.85% of the analytical value at max slope.

After discretizing in space and time, we are left with a matrix equation. In 1D, the matrices are block diagonal, with each block corresponding to a single element. The block size is (p+1)^2, where p is the polynomial degree of the approximation on that element. In higher dimensions, the system matrices are formed by applying kronecker tensor products of the 1D operators. Thus

There is very little storage overhead for the higher dimensional systems. Moreover, the kronecker tensor formulation reduces the order of operations from O(n^(d+2)) to O(n^(d+1)) for a d-dimensional calculation with n grid points.

Linear advection of a 2D Gaussian wave for one time peroid.

Uncoupled Stiffness Matrix (Laplacian Operator)

We have implemented 2D Advection Diffusion Equation, including non-linear advection. This code allows for periodic, and Dirichlet Boundary conditions. We have designed our code in an object oriented fashion to allow us to add new capabilities. Future Research plans include implementing the Navier-Stokes Equations, curved elements, adaptivity in higher dimensions, and Preconditioned Iterative Schemes.

We determine a global ordering of each node of the discretization, this ordering determines the structure of the coupled system operators. However, because we are using iterative solvers, we never need to construct the global system. Instead, we compute the contribution of each element to the solution, and then use a weighted sum to obtain the global solution, taking into account where the same node is used on two elements. For example the global solution at node 7 in the global ordering is obtained by adding the solution at node 7 from element 1 and node 1 from element 2, and then dividing by 2. The name of this weighted summation operation is referred to as Direct-Stiffness-Summation (DSS). In order to parallelize the SEM, each processor is assigned a group of elements whose solution it will compute. Each time a solution on those elements is computed, a parallel DSS is used to determine the nodal values at (non-local) element Boundaries. For example if element 1 and element 2 are on separate processors, each time a solution is computed on both of them, the solution at node 7 on processor 1 is summed with the solution at node 1 on processor 2. The result is then divided by two and stored on both processors as the value of the solution at that node. In order to determine the dependencies between processors for complicated geometries, a parallel bin sort is used. This past summer we implemented this parallel DSS technique using MPI.

Element 2

Element 1

References

[1] Anil Deane. Spectral and spectral-element methods: Lecture notes in high performance computational physics. NASA Contractor Report 203877, 1997.

[2] P.F. Fischer. An overlapping schwarz method for spectral element solution of the incompressible navier-stokes equations. Journal of Computational Physics, 1997

[3] S.J. Sherwin G.E. Karniadakis. Spectral/hp Element Methods for CFD. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford, 1999.

[4] Rainald Löhner. An adaptive finite element scheme for transient problems in cfd. Computational Methods in Applied Mechanics and Engineering, 61:323338, 1987.

[5] P. Aaron Lott. Project website. http://www.lcv.umd.edu/ ~palott/research/graduate/663/.

[6] E.H. Mund M.O. Deville, P.F. Fischer. High-Order Methods for Incompressible Fluid Flows. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2002.

.

Global Ordering

Element 1

Element 2

Local Ordering