High Order Accurate Solution of Flow Past a Circular Cylinder

A high order method is applied to time-dependent incompressible flow around a circular cylinder geometry. The space discretization employs compact fourth-order difference operators. In time we discretize with a second-order semi-implicit scheme. A large linear system of equations is solved in each time step by a combination of outer and inner iterations. An approximate block factorization of the system matrix is used for preconditioning. Well posed boundary conditions are obtained by an integral formulation of boundary data including a condition on the pressure. Two-dimensional flow around a circular cylinder is studied for Reynolds numbers in the range 7 ≤ R ≤ 180. The results agree very well with the data known from numerical and experimental studies in the literature.     

An efficient solver for the RANS equations and a one-equation turbulence model

A three-stage Runge-Kutta (RK) scheme with multigrid and an implicit preconditioner has been shown to be an effective solver for the fluid dynamic equations. Using the algebraic turbulence model of Baldwin and Lomax, this scheme has been used to solve the compressible Reynolds-averaged Navier–Stokes (RANS) equations for transonic and low-speed flows. In this paper we focus on the convergence of the RK/Implicit scheme when the effects of turbulence are represented by the one-equation model of Spalart and Allmaras. With the present scheme the RANS equations and the partial differential equation of the turbulence model are solved in a loosely coupled manner. This approach allows the convergence behavior of each system to be examined. Point symmetric Gauss-Seidel supplemented with local line relaxation is used to approximate the inverse of the implicit operator of the RANS solver. To solve the turbulence equation we consider three alternative methods: diagonally dominant alternating direction implicit (DDADI), symmetric line Gauss-Seidel (SLGS), and a two-stage RK scheme with implicit preconditioning. Computational results are presented for airfoil flows, and comparisons are made with experimental data. We demonstrate that the two-dimensional RANS equations and a transport-type equation for turbulence modeling can be efficiently solved with an indirectly coupled algorithm that uses RK/Implicit schemes.     

A parallel fully coupled implicit domain decomposition method for numerical simulation of microfluidic mixing in 3D

A parallel fully coupled implicit fluid solver based on a Newton–Krylov–Schwarz algorithm is developed on top of the Portable, Extensible Toolkit for Scientific computation for the simulation of microfluidic mixing described by the three-dimensional unsteady incompressible Navier–Stokes equations. The popularly used fractional step method, originally designed for high Reynolds number flows, requires some modification of the inviscid-type pressure boundary condition in order to reduce the divergence error near the wall. On the other hand, the fully coupled approach works well without any special treatment of the boundary condition for low Reynolds number microchannel flows. A key component of the algorithm is an additive Schwarz preconditioner, which is used to accelerate the convergence of a linear Krylov-type solver for the saddle-point-type Jacobian systems. As a test case, we carefully study a three-dimensional passive serpentine micromixer and report the parallel performance of the algorithm obtained on a parallel machine with more than one hundred processors.     

Efficient Nonlinear Solvers for Nodal High-Order Finite Elements in 3D

Conventional high-order finite element methods are rarely used for industrial problems because the Jacobian rapidly loses sparsity as the order is increased, leading to unaffordable solve times and memory requirements. This effect typically limits order to at most quadratic, despite the favorable accuracy and stability properties offered by quadratic and higher order discretizations. We present a method in which the action of the Jacobian is applied matrix-free exploiting a tensor product basis on hexahedral elements, while much sparser matrices based on Q ₁ sub-elements on the nodes of the high-order basis are assembled for preconditioning. With this “dual-order” scheme, storage is independent of spectral order and a natural taping scheme is available to update a full-accuracy matrix-free Jacobian during residual evaluation. Matrix-free Jacobian application circumvents the memory bandwidth bottleneck typical of sparse matrix operations, providing several times greater floating point performance and better use of multiple cores with shared memory bus. Computational results for the p-Laplacian and Stokes problem, using block preconditioners and AMG, demonstrate mesh-independent convergence rates and weak (bounded) dependence on order, even for highly deformed meshes and nonlinear systems with several orders of magnitude dynamic range in coefficients. For spectral orders around 5, the dual-order scheme requires half the memory and similar time to assembled quadratic (Q ₂) elements, making it very affordable for general use.     

Two-Level Space–Time Domain Decomposition Methods for Flow Control Problems

For time-dependent control problems, the class of sub-optimal algorithms is popular and the parallelization is usually applied in the spatial dimension only. In the paper, we develop a class of fully-optimal methods based on space–time domain decomposition methods for some boundary and distributed control of fluid flow and heat transfer problems. In the fully-optimal approach, we focus on the use of an inexact Newton solver for the necessary optimality condition arising from the implicit discretization of the optimization problem and the use of one-level and two-level space–time overlapping Schwarz preconditioners for the Jacobian system. We show that the numerical solution from the fully-optimal approach is generally better than the solution from the sub-optimal approach in terms of meeting the objective of the optimization problem. To demonstrate the robustness and parallel scalability and efficiency of the proposed algorithm, we present some numerical results obtained on a parallel computer with a few thousand processors.     

Extension of the Complete Flux Scheme to Systems of Conservation Laws

We present the extension of the complete flux scheme to advection-diffusion-reaction systems. For stationary problems, the flux approximation is derived from a local system boundary value problem for the entire system, including the source term vector. Therefore, the numerical flux vector consists of a homogeneous and an inhomogeneous component, corresponding to the advection-diffusion operator and the source term, respectively. For time-dependent systems, the numerical flux is determined from a quasi-stationary boundary value problem containing the time-derivative in the source term. Consequently, the complete flux scheme results in an implicit semidiscretization. The complete flux scheme is validated for several test problems.     

Using exact Jacobians in an implicit Newton–Krylov method

In an implicit Newton–Krylov method for inviscid, compressible fluid flow, the derivation of the analytic flux Jacobian can become quite complicated depending on the complexity of the numerical flux calculation. Practically, the derivation of the exact Jacobian by hand is an unrealistic option because of the enormous man-hour investment needed. In this work, automatic differentiation is used to evaluate the exact Jacobian of upwind schemes implemented in the flow solver QUADFLOW. We compare the use of exact Jacobians and Jacobians numerically approximated by first-order forward differences. For a two-dimensional airfoil under three different flight conditions (quasi-incompressible flow, compressible subsonic flow, and transonic flow), we show that the robustness and performance of the present finite volume scheme is significantly improved by using exact Jacobians.     

A shape design system using volumetric implicit PDEs

Solid modeling based on partial differential equations (PDEs) can potentially unify both geometric constraints and functional requirements within a single design framework to model real-world objects via its explicit, direct integration with parametric geometry. In contrast, implicit functions indirectly define geometric objects as the level-set of underlying scalar fields. To maximize the modeling potential of PDE-based methodology, in this paper we tightly couple PDEs with volumetric implicit functions in order to achieve interactive, intuitive shape representation, manipulation, and deformation. In particular, the unified approach can reconstruct the PDE geometry of arbitrary topology from scattered data points or a set of sketch curves. We make use of elliptic PDEs for boundary value problems to define the volumetric implicit function. The proposed implicit PDE model has the capability to reconstruct a complete solid model from partial information and facilitates the direct manipulation of underlying volumetric datasets via sketch curves and iso-surface sculpting, deformation of arbitrary interior regions, as well as a set of CSG operations inside the working space. The prototype system that we have developed allows designers to interactively sketch the curve outlines of the object, define intensity values and gradient directions, and specify interpolatory points in the 3D working space. The governing implicit PDE treats these constraints as generalized boundary conditions to determine the unknown scalar intensity values over the entire working space. The implicit shape is reconstructed with specified intensity value accordingly and can be deformed using a set of sculpting toolkits. We use the finite-difference discretization and variational interpolating approach with the localized iterative solver for the numerical integration of our PDEs in order to accommodate the diversity of generalized boundary and additional constraints.     

Using Distributed Computers to Deterministically Approximate Higher Dimensional Convection Diffusion Equations

A new approach to solving D> 3 spatial dimensional convection-diffusion equation on clusters of workstations is derived by exploiting the stability and scalability of the combination of a generalized D dimensional high-order compact (HOC) implicit finite difference scheme and parallelized GMRES(m). We then consider its application to multifactor Option pricing using the Black–Scholes equation and further show that an isotropic fourth order compact difference scheme is numerically stable and determine conditions under which its coefficient matrix is positive definite. The performance of GMRES(m) on distributed computers is limited by the inter-processor communication required by the matrix-vector multiplication. It is shown that the compact scheme requires approximately half the number of communications as a non-compact difference scheme of the same order of truncation error. As the dimensionality is increased, the ratio of computation that can be overlapped with communication also increases. CPU times and parallel efficiency graphs for single time step approximation of up to a 7D HOC scheme on 16 processors confirm the numerical stability constraint and demonstrate improved parallel scalability over non-compact difference schemes.     

Globalized matrix-explicit Newton-GMRES for the high-order accurate solution of the Euler equations

Implicit methods for finite-volume schemes on unstructured grids typically rely on a matrix-free implementation of GMRES and an explicit first-order accurate Jacobian for preconditioning. Globalization is typically achieved using a global timestep or a CFL based local timestep. We show that robustness of the globalization can be improved by supplementing the pseudo-timestepping method commonly used with a line search method. The number of timesteps required for convergence can be reduced by using a timestep that scales with the local residual. We also show that it is possible to form the high-order Jacobian explicitly at a reasonable computational cost. This is demonstrated for cases using both limited and unlimited reconstruction. This exact Jacobian can be used for preconditioning and directly in the GMRES method. The benefits of improvements in preconditioning and the elimination of residual evaluations in the inner iterations of the matrix-free GMRES method are substantial. Computational results focus on second- and fourth-order accurate schemes with some results for the third-order scheme. Overall computational cost for the matrix-explicit method is lower than the matrix-free method for all cases. The fourth-order matrix-explicit scheme is a factor of 1.6-3 faster than the matrix-free scheme while requiring about 50-100% more memory.     

An iterative matrix-free method in implicit immersed boundary/continuum methods

The objective of this paper is to present an iterative solution strategy for implicit immersed boundary/continuum methods. An overview of the newly proposed immersed continuum method in conjunction with the traditional immersed boundary method will also be presented. As a key ingredient of the fully implicit time integration, a matrix-free combination of Newton–Raphson iteration and GMRES iterative linear solver is proposed.     

Unconditionally stable ADI scheme of higher-order for linear hyperbolic equations

An unconditionally stable alternating direction implicit (ADI) method of higher-order in space is proposed for solving two- and three-dimensional linear hyperbolic equations. The method is fourth-order in space and second-order in time. The solution procedure consists of a multiple use of one-dimensional matrix solver which produces a computational cost effective solver. Numerical experiments are conducted to compare the new scheme with the existing scheme based on second-order spatial discretization. The effectiveness of the new scheme is exhibited from the numerical results.     

Low-rank tensor completion with fractional-Jacobian-extended tensor regularization for multi-component visual data inpainting

Several low-rank tensor completion methods have been integrated with total variation (TV) regularization to retain edge information and promote piecewise smoothness. In this paper, we first construct a fractional Jacobian matrix to nonlocally couple the structural correlations across components and propose a fractional-Jacobian-extended tensor regularization model, whose energy functional was designed proportional to the mixed norm of the fractional Jacobian matrix. Consistent regularization could thereby be performed on each component, avoiding band-by-band TV regularization and enabling effective handling of the contaminated fine-grained and complex details due to the introduction of a fractional differential. Since the proposed spatial regularization is linear convex, we further produced a novel fractional generalization of the classical primal-dual resolvent to develop its solver efficiently. We then combined the proposed tensor regularization model with low-rank constraints for tensor completion and addressed the problem by employing the augmented Lagrange multiplier method, which provides a splitting scheme. Several experiments were conducted to illustrate the performance of the proposed method for RGB and multispectral image restoration, especially its abilities to recover complex structures and the details of multi-component visual data effectively.     

Polynomial Time-Marching for Nonreflecting Boundary Problems

The newly developed polynomial time-marching technique has been successfully extended to nonperiodic boundary condition cases. In this paper, a special non-periodic boundary condition, nonreflecting or absorbing boundary condition, is incorporated into the pseudospectral polynomial time-marching scheme. Thus, this accurate and stable time-dependent PDE solver can be applied to some open domain or free space problems. The balanced overall spectral accuracy is illustrated by some numerical experiments in the one-dimensional case. The error goes to zero at a rate faster than many fixed orders of the finite-difference scheme. The order of the absorbing boundary approximation is addressed in one-dimension. Also, in the two-dimensional case, a 2nd-order absorbing approximation has been incorporated into the polynomial time-marching scheme with Chebyshev collocation in space. Comparison with the previous finite-difference implementation indicates that the high stability and efficiency of the polynomial time-marching remains. The overall accuracy is mainly limited by the 2nd-order absorbing approximation.     

Efficient exponential compact higher order difference scheme for convection dominated problems

Present work is the development of a finite difference scheme based on Richardson extrapolation technique. It gives an exponential compact higher order scheme (ECHOS) for two-dimensional linear convection-diffusion equations (CDE). It uses a compact nine point stencil, over which the governing equations are discretized for both fine and coarse grids. The resulting algebraic systems are solved using a line iterative approach with alternate direction implicit (ADI) procedure. Combining the solutions over fine and coarse grids, initially a sixth order solution over coarse grid points is obtained. The resultant solution is then extended to finer grid by interpolation derived from the difference operator. The convergence of the iterative procedure is guaranteed as the coefficient matrix of the developed scheme satisfies the conditions required to be monotone. The higher order accuracy and better rate of convergence of the developed algorithm have been demonstrated by solving numerous model problems.     

A new way for constructing high accuracy shock-capturing generalized compact difference schemes

The generalized compact (GC) schemes and some of their important properties are presented. And a new way for constructing high order accuracy and high-resolution GC schemes is presented. The schemes constructed by using this way could satisfy some principles and demands prescribed in advance to ensure some desired properties to the schemes, such as the principle about suppression of the oscillations, the principle of stability, the order of accuracy and number of scheme points, etc. As two examples, a three-point third-order compact scheme and a three-point fifth-order GC scheme satisfying the principle about suppression of the oscillations and the principle of stability are described in this paper. Numerical results show that these schemes are shock-capturing. The time-dependent boundary conditions proposed by Thompson are well employed when the algorithm is applied to the Euler equations of gas dynamics. Fourier analysis shows that the resolution characteristics are spectral-like.     

A parallel linear solver for multilevel Toeplitz systems with possibly several right-hand sides

A Toeplitz matrix has constant diagonals; a multilevel Toeplitz matrix is defined recursively with respect to the levels by replacing the matrix elements with Toeplitz blocks. Multilevel Toeplitz linear systems appear in a wide range of applications in science and engineering. This paper discusses an MPI implementation for solving such a linear system by using the conjugate gradient algorithm. The implementation techniques can be generalized to other iterative Krylov methods besides conjugate gradient. These techniques include the use of an arbitrary dimensional process grid for handling the multilevel Toeplitz structure, a communication-hiding approach for performing matrix–vector multiplications, the incorporation of multilevel circulant preconditioning for accelerating convergence, an efficient orthogonalization manager for detecting linear dependence in block iterations, and an algorithmic rearrangement to eliminate all-reduce synchronizations. The combined use of these techniques leads to a scalable solver for large multilevel Toeplitz systems, possibly with several right-hand sides. We show experimental results on matrices of size up to the order of one billion with nearly perfect scaling by using up to 1024 MPI processes. We also demonstrate an application of the solver in parameter estimation for analyzing large-scale climate data.     

A CCD–ADI method for unsteady convection–diffusion equations

With a combined compact difference scheme for the spatial discretization and the Crank–Nicolson scheme for the temporal discretization, respectively, a high-order alternating direction implicit method (ADI) is proposed for solving unsteady two dimensional convection–diffusion equations. The method is sixth-order accurate in space and second-order accurate in time. The resulting matrix at each ADI computation step corresponds to a triple-tridiagonal system which can be effectively solved with a considerable saving in computing time. In practice, Richardson extrapolation is exploited to increase the temporal accuracy. The unconditional stability is proved by means of Fourier analysis for two dimensional convection–diffusion problems with periodic boundary conditions. Numerical experiments are conducted to demonstrate the efficiency of the proposed method. Moreover, the present method preserves the higher order accuracy for convection-dominated problems.     

Algebraic multigrid methods based on element preconditioning

This paper presents a new algebraic multigrid (AMG) solution strategy for large linear systems with a sparse matrix arising from a finite element discretization of some self-adjoint, second order, scalar, elliptic partial differential equation. The AMG solver is based on Ruge/Stübens method. Ruge/Stübens algorithm is robust for M-matrices, but unfortunately the “region of robustness“ between symmetric positive definite M-matrices and general symmetric positive definite matrices is very fuzzy.

For this reason the so-called element preconditioning technique is introduced in this paper. This technique aims at the construction of an M-matrix that is spectrally equivalent to the original stiffness matrix. This is done by solving small restricted optimization problems. AMG applied to the spectrally equivalent M-matrix instead of the original stiffness matrix is then used as a preconditioner in the conjugate gradient method for solving the original problem.

The numerical experiments show the efficiency and the robustness of the new preconditioning method for a wide class of problems including problems with anisotropic elements.