Advection-dominated Problems Research Articles

A Simplified and Flexible Variant of GCROT for Solving Nonsymmetric Linear Systems

There is a need for flexible iterative solvers that can solve large-scale ($>10^6$ unknowns) nonsymmetric sparse linear systems to a small tolerance. Among flexible solvers, flexible GMRES (FGMRES) is attractive because it minimizes the residual norm over a particular subspace. In practice, FGMRES is often restarted periodically to keep memory and work requirements reasonable; however, like restarted GMRES, restarted FGMRES can suffer from stagnation. This has led us to develop a flexible variant of the Krylov linear solver GCROT (generalized conjugate residual with inner orthogonalization and outer truncation). Unlike the original GCROT algorithm, the proposed GCROT variant uses a simplified truncation strategy similar to loose GMRES (LGMRES). This modification is motivated by numerical experiments that suggest the specific subspace retained in the outer iteration of GCROT is less important than its size. The flexible GCROT variant appears to be well suited for advection-dominated problems. In particular, when applied to an adjoint problem from computational aerodynamics, the proposed GCROT variant is robust and efficient compared with several popular truncated Krylov subspace methods. Finally, a flexible version of LGMRES is easily constructed by recognizing algorithmic similarities to GCROT.

Computational Approaches to Solving Equations Arising from Wound Healing

In the wound healing process, the cell movement associated with chemotaxis generally outweighs the movement associated with random motion, leading to advection-dominated mathematical models of wound healing. The equations in these models must be solved with care, but often inappropriate approaches are adopted. Two one-dimensional test problems arising from advection-dominated models of wound healing are solved using four algorithms--MATLAB's inbuilt routine pdepe.m, the Numerical Algorithms Group routine d03pcf.f, and two finite volume methods. The first finite volume method is based on a first-order upwinding treatment of chemotaxis terms and the second on a flux limiting approach. The first test problem admits an analytic solution which can be used to validate the numerical results by analyzing two measures of the error for each method: the average absolute difference and a mass balance error. These criteria as well as the visual comparison between the numerical methods and the exact solution lead us to conclude that flux limiting is the best approach to solving advection-dominated wound healing problems numerically in one dimension. The second test problem is a coupled nonlinear three species model of wound healing angiogenesis. Measurement of the mass balance error for this test problem further confirms our hypothesis that flux limiting is the most appropriate method for solving advection-dominated governing equations in wound healing models. We also consider two two-dimensional test problems arising from wound healing, one that admits an analytic solution and a more complicated problem of blood vessels growth into a devascularized wound bed. The results from the two-dimensional test problems also demonstrate that the flux limiting treatment of advective terms is ideal for an advection-dominated problem.

A Subgrid-Scale Stabilized Finite Element Method for Multicomponent Reactive Transport through Porous Media

Standard Galerkin finite element methods (GFEM) lack stability in solving advection-dominated solute transport in porous media. They usually require prohibitively fine grids and extremely small time steps to solve for advection-dominated problems. The algebraic subgrid-scale stabilized (ASGS) finite element method has been proved to overcome such problems for single-species reactive transport. Its potential for dealing with multicomponent reactive transport has not yet been explored. Here we present a numerical formulation of ASGS for steady and transient multicomponent reactive transport. Subgrid-scale transport equations are solved first by using an ASGS approximation and their solutions are substituted back into the grid-scale equations. A sequential iteration approach (SIA) is used to solve for coupled transport and chemical equations. Coupling of ASGS and SIA, ASGS+SIA, has been implemented in a reactive transport code, CORE2D V4, and verified for conservative solute transport. ASGS+SIA has been tested for a wide range of 1-D transient multicomponent reactive transport problems involving different types of chemical reactions such as: (1) Kinetically controlled aqueous species degradation, (2) Kinetic mineral dissolution, (3) Serial-parallel decay networks, and (4) Cation exchange and pyrite oxidation at local equilibrium. ASGS+SIA always provides accurate solutions and therefore offers an efficient option to solve for advection-dominated multicomponent reactive transport problems.

Adaptive Time-Stepping for Incompressible Flow Part I: Scalar Advection-Diffusion

Even the simplest advection-diffusion problems can exhibit multiple time scales. This means that robust variable step time integrators are a prerequisite if such problems are to be efficiently solved computationally. The performance of the second order trapezoid rule using an explicit Adams-Bashforth method for error control is assessed in this work. This combination is particularly well suited to long time integration of advection-dominated problems. Herein it is shown that a stabilized implementation of the trapezoid rule leads to a very effective integrator in other situations: specifically diffusion problems with rough initial data; and general advection-diffusion problems with different physical time scales governing the system evolution.

A review and numerical assessment of the random walk particle tracking method

We review the basic mathematical concepts of random walk particle tracking (RWPT) and its advantages and limitations. Three different numerical approaches to overcome the local mass conservation problem of the random walk methodology are examined: (i) the interpolation method, (ii) the reflection principle, and (iii) the generalized stochastic differential equations (GSDE). Analytical solutions of the spatial moments for a two-layer system are compared to model predictions using the different techniques and results demonstrate that the interpolation method reproduces correctly average velocity, but fails to reproduce macrodispersion at higher hydraulic conductivity contrasts between the two layers. On the contrary, the reflection principle and the GSDE approach underestimate average velocity, but reproduce macrodispersion better for high contrasts. The different behavior is based on an artificial shift of mass for increasing heterogeneities for the GSDE approach and the reflection principle, whereas the interpolation method suffers from the smoothing of the dispersion tensor. The behavior of these approaches was furthermore analyzed in two-dimensional heterogeneous hydraulic conductivity fields, which are characterized by different random function models. Solute transport was simulated correctly by all three approaches for the reference fields having Gaussian structures or non-Gaussian structures with an isotropic spatial correlation, even for a variance of the natural log of hydraulic conductivity of σ ln K 2 = 4. However, for the non-Gaussian model with a strong anisotropic spatial correlation and a variance of σ ln K 2 = 2 and higher, the interpolation method was the only technique modelling solute transport correctly. Furthermore, we discuss the general applicability of random walk particle tracking in comparison to the standard transport models and conclude that in advection-dominated problems using a high spatial discretization or requiring the performance of many model runs, RWPT represents a good alternative for modelling contaminant transport.

Residual-based stabilized higher-order FEM for advection-dominated problems

We reconsider the numerical solution of linear(ized) advection–diffusion–reaction problems using higher-order finite elements together with stabilized Galerkin methods of streamline-diffusion type (SUPG) and with shock-capturing stabilization. The analysis improves the a priori analysis in our previous paper [T. Knopp, G. Lube, G. Rapin, Stabilized finite element methods with shock capturing for advection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 2997–3013]. The theoretical results are supported by some numerical experiments.

An adaptive fixed mesh method for modelling and simulating of unsteady two-phase flows with sharp physical interfaces

We present in this paper a new computational method for simulation of two-phase flow problems with moving boundaries and sharp physical interfaces. An adaptive interface-capturing technique (ICT) of the Eulerian type is developed for capturing the motion of the interfaces (free surfaces) in an unsteady flow state. The adaptive method is mainly based on the relative boundary conditions of the zero pressure head, at which the interface is corresponding to a free surface boundary. The definition of the free surface boundary condition is used as a marker for identifying the position of the interface (free surface) in the two-phase flow problems. An initial-value-problem (IVP) partial differential equation (PDE) is derived from the dynamic conditions of the interface, and it is designed to govern the motion of the interface in time. In this adaptive technique, the Navier–Stokes equations written for two incompressible fluids together with the IVP are solved numerically over the flow domain. An adaptive mass conservation algorithm is constructed to govern the continuum of the fluid. The finite element method (FEM) is used for the spatial discretization and a fully coupled implicit time integration method is applied for the advancement in time. FE-stabilization techniques are added to the standard formulation of the discretization, which possess good stability and accuracy properties for the numerical solution. The adaptive technique is tested in simulation of some numerical examples. With the test problems presented here, we demonstrated that the adaptive technique is a simple tool for modelling and computation of complex motion of sharp physical interfaces in convection–advection-dominated flow problems. We also demonstrated that the IVP and the evolution of the interface function are coupled explicitly and implicitly to the system of the computed unknowns in the flow domain. Copyright © 2005 John Wiley & Sons, Ltd.

Comparison of Galerkin and control volume finite element for advection-diffusion problems

The control volume finite element method (CVFEM) was developed to combine the local numerical conservation property of control volume methods with the unstructured grid and generality of finite element methods (FEMs). Most implementations of CVFEM include mass-lumping and upwinding techniques typical of control volume schemes. In this work we compare, via numerical error analysis, CVFEM and FEM utilizing consistent and lumped mass implementations, and stabilized Petrov–Galerkin streamline upwind schemes in the context of advection–diffusion processes. For this type of problem, we find no apparent advantage to the local numerical conservation aspect of CVFEM as compared to FEM. The stabilized schemes improve accuracy and degree of positivity on coarse grids, and also reduce iteration counts for advection-dominated problems. Published in 2005 by John Wiley & Sons, Ltd.

A numerical method for solving a scalar advection-dominated transport equation with concentration-dependent sources

A numerical scheme, upstream biased Eulerian algorithm for transport equations with sources (UpBEATES), is developed for solving a scalar advective-dominated transport equation with concentration-independent and -dependent source terms. A control-volume method is used for spatial discretization. Time integration is invoked to yield a discrete system of integrated-flux-integrated-source form equations. The Bott's upstream-biased Eulerian advection scheme [Moneatry and Weakly Review 117 (1989a) 1006; Moneatry and Weakly Review 117 (1989b) 2633] is employed for approximating advective fluxes. A two-level time weighting scheme is employed for the dispersive fluxes. An upstream-biased Eulerian algorithm is proposed for the concentration-dependent source term. Flux and source limiters are developed to ensure non-negative evolution of the scalar concentration field. Numerical experiments were presented to illustrate its performance in comparison with theoretical solutions and those of conventional methods. The proposed scheme is mass-conservative, produces non-negative concentration values, exhibits low numerical dispersion, and is efficient for advection-dominated problems with concentration-dependent source terms. Like other Eulerian schemes, the Courant–Friedrich–Levy (CFL) stability criterion has to be met.

H -ADAPTIVITY AND "HONEST" GFEM FOR ADVECTION-DOMINATED TRANSPORT

A standard two-dimensional Galerkin finite-element method (GFEM) code for coupled Navier-Stokes and energy equations is used with h -adaptive meshing based on a posteriori error estimation using the superconvergent patch recovery technique for solving a range of advection-dominated transport problems. It is demonstrated that such a method provides a highly effective, simple, and efficient way of dealing with the perennial problems in numerical modeling of advection-dominated transport, such as oscillations or wiggles with central difference-type discretizations (such as GFEM) and numerical ("false") diffusion when wiggle-suppressant schemes are used. Additionally, the auto-adaptive finite-element method provides a powerful means of achieving optimal solutions without having to predefine a mesh, which may be either inadequate or too expensive. A number of benchmark problems are presented as application examples for this method before solving a problem of natural convection in an air-filled cavity with various orientations, for which experimental results are available.

A lattice BGK model for advection and anisotropic dispersion equation

This paper presents a lattice Boltzmann model (LBM) for 2-D advection and anisotropic dispersion equation (AADE) based on the Bhatnagar, Gross and Krook (BGK) model. In the proposed model, the particle speed space is discretized using a rectangular lattice that has four speeds in nine directions, and the single relaxation time is assumed to be directionally dependent. To ensure that the collision is mass-invariant when the relaxation time is directionally dependent, the concentration is calculated from a weighted summation of the particle distribution functions. The proposed model was verified against benchmark problems and the finite difference solution of solute transport with spatially variable dispersion coefficients and non-uniform velocity field. The significant results are that it conserves mass perfectly and offers accurate and efficient solutions for both dispersion-dominated and advection-dominated problems.

Arnoldi and Crank-Nicolson methods for integration in time of the transport equation

A comparison is made between the Arnoldi reduction method and the Crank–Nicolson method for the integration in time of the advection–diffusion equation. This equation is first discretized in space by the classic finite element (FE) approach, leading to an unsymmetric first-order differential system, which is then solved by the aforementioned methods. Arnoldi reduces the native FE equations to a much smaller set to be efficiently integrated in the Arnoldi vector space by the Crank–Nicolson scheme, with the solution recovered back by a standard Rayleigh–Ritz procedure. Crank–Nicolson implements a time marching scheme directly on the original first-order differential system. The computational performance of both methods is investigated in two- and three-dimensional sample problems with a size up 30 000. The results show that in advection-dominated problems less then 100 Arnoldi vectors generally suffice to give results with a 10−3–10−4 difference relative to the direct Crank–Nicolson solution. However, while the CPU time with the Crank–Nicolson starts from zero and increases linearly with the number of time steps used in the simulation, the Arnoldi requires a large initial cost to generate the Arnoldi vectors with subsequently much less expensive dynamics for the time integration. The break-even point is problem-dependent at a number of time steps which may be for some problems up to one order of magnitude larger than the number of Arnoldi vectors. A serious limitation of Arnoldi is the requirement of linearity and time independence of the flow field. It is concluded that Arnoldi can be cheaper than Crank–Nicolson in very few instances, i.e. when the solution is needed for a large number of time values, say several hundreds or even 1000, depending on the problem. Copyright © 2001 John Wiley & Sons, Ltd.

Characteristics of distribution and transport of petroleum contaminants in fracture-karst water in Zibo Area, Shandong Province, China

Fracture-karst water is an important water resource for the water supply in North China. Petroleum contamination is one of the most problematic types of the groundwater pollution. The characteristics of distribution and transport of the petroleum contaminants in fracture-karst water are different from those in porous water. The flow velocity of fracture-karst water is much faster than the velocity of porous water on an average. Therefore, contaminant transport in fracture-karst water is an absolute advection-dominated problem. The plume of the petroleum contamination may extend to several kilometers from pollution sources. It was not caused by the oil pool floating on the water table but by the oil components dissolved and scattered in groundwater. The distribution of the petroleum contaminants over space are concentrated in the strong conductive zone on the plane. On the vertical section the highest concentration of the oil contaminants appeared in the strata where the contamination sources were located. The concentrations of the oil contaminants in wells changed greatly over time. Therefore, the curves of concentration versus time fluctuated greatly. The reasons are as follows. (a) Fracture-karst water has a very great velocity. (b) Local flow fields which were caused by pumping and stoppage in some wells changed frequently. (c) In fracture-karst aquifer the transport channels are complicated. (d) Residual oil in vadose zone was leached after rainfall. It is of great practical value for the control and remediation of petroleum contamination in fracture-karst aquifer to understand those characteristics.

On exponential variable transformation based boundary element formulation for advection–diffusion problems

This paper presents an alternative boundary element formulation based on an exponential variable transformation for the steady state advection–diffusion problems. Use of this transformation converts the advection–diffusion equation into modified Helmholtz equation. Standard boundary element formulation is used to obtain the integral equation for the resulting modified Helmholtz equation. Three different strategies for the solution of this boundary integral equation are discussed. These include (a) the solution of the modified Helmholtz equation followed by transformation at post-processing stage, (b) use of inverse transformation before boundary element discretization, and (c) use of inverse transformation after the boundary element discretization to obtain the discretized equations for the advection–diffusion problem. Computational features of each of these strategies have been discussed. Numerical results are presented for a set of standard test problems. It is found that strategy (a) is the most accurate and efficient choice for diffusion-dominated ( Pe≪1) as well as low Peclét number ( Pe≤20) situations. Strategy (c) is as accurate as strategy (a), but is of little practical use owing to its computational inefficiency. Further, strategy (b), which is same as the boundary element formulation using the fundamental solution of the adjoint of the advection–diffusion equation, is the most accurate and efficient choice for advection-dominated problems.

Global and Local Error Analysis for the Residual-Free Bubbles Method Applied to Advection-Dominated Problems

We prove the stability and a priori global and local error analysis for the residual-free bubbles (RFB) finite element method applied to advection-dominated advection-diffusion problems.

Implicit and semi-implicit schemes: Algorithms

This study formulates general guidelines to extend an explicit code with a great variety of implicit and semi-implicit time integration schemes. The discussion is based on their specific implementation in the Versatile Advection Code, which is a general purpose software package for solving systems of non-linear hyperbolic (and/or parabolic) partial differential equations, using standard high resolution shock capturing schemes. For all combinations of explicit high resolution schemes with implicit and semi-implicit treatments, it is shown how second-order spatial and temporal accuracy for the smooth part of the solutions can be maintained. Strategies to obtain steady state and time accurate solutions implicitly are discussed. The implicit and semi-implicit schemes require the solution of large linear systems containing the Jacobian matrix. The Jacobian matrix itself is calculated numerically to ensure the generality of this implementation. Three options are discussed in terms of applicability, storage requirements and computational efficiency. One option is the easily implemented matrix-free approach, but the Jacobian matrix can also be calculated by using a general grid masking algorithm, or by an efficient implementation for a specific Lax–Friedrich-type total variation diminishing (TVD) spatial discretization. The choice of the linear solver depends on the dimensionality of the problem. In one dimension, a direct block tridiagonal solver can be applied, while in more than one spatial dimension, a conjugate gradient (CG)-type iterative solver is used. For advection-dominated problems, preconditioning is needed to accelerate the convergence of the iterative schemes. The modified block incomplete LU-preconditioner is implemented, which performs very well. Examples from two-dimensional hydrodynamic and magnetohydrodynamic computations are given. They model transonic stellar outflow and recover the complex magnetohydrodynamic bow shock flow in the switch-on regime found in De Sterck et al. [Phys. Plasmas, 5, 4015 (1998)]. Copyright © 1999 John Wiley & Sons, Ltd.

Implicit and semi‐implicit schemes: Algorithms

This study formulates general guidelines to extend an explicit code with a great variety of implicit and semi-implicit time integration schemes. The discussion is based on their specific implementation in the Versatile Advection Code, which is a general purpose software package for solving systems of non-linear hyperbolic (and/or parabolic) partial differential equations, using standard high resolution shock capturing schemes. For all combinations of explicit high resolution schemes with implicit and semi-implicit treatments, it is shown how second-order spatial and temporal accuracy for the smooth part of the solutions can be maintained. Strategies to obtain steady state and time accurate solutions implicitly are discussed. The implicit and semi-implicit schemes require the solution of large linear systems containing the Jacobian matrix. The Jacobian matrix itself is calculated numerically to ensure the generality of this implementation. Three options are discussed in terms of applicability, storage requirements and computational efficiency. One option is the easily implemented matrix-free approach, but the Jacobian matrix can also be calculated by using a general grid masking algorithm, or by an efficient implementation for a specific Lax–Friedrich-type total variation diminishing (TVD) spatial discretization. The choice of the linear solver depends on the dimensionality of the problem. In one dimension, a direct block tridiagonal solver can be applied, while in more than one spatial dimension, a conjugate gradient (CG)-type iterative solver is used. For advection-dominated problems, preconditioning is needed to accelerate the convergence of the iterative schemes. The modified block incomplete LU-preconditioner is implemented, which performs very well. Examples from two-dimensional hydrodynamic and magnetohydrodynamic computations are given. They model transonic stellar outflow and recover the complex magnetohydrodynamic bow shock flow in the switch-on regime found in De Sterck et al. [Phys. Plasmas, 5, 4015 (1998)]. Copyright © 1999 John Wiley & Sons, Ltd.

A numerical eigenvalue study of preconditioned non-equilibrium transport equations

The finite element integration of non-equilibrium contaminant transport in porous media yields sparse, unsymmetric, real or complex equations, which may be solved by iterative projection methods, such as Bi-CGSTAB and TFQMR, on condition that they are effectively preconditioned. To ensure a fast convergence, the eigenspectrum of the preconditioned equations has to be very compact around unity. Compactness is generally measured by the spectral condition number. In difficult advection-dominated problems, however, the condition number may be large and nevertheless, convergence may be good. A numerical study of the preconditioned eigenspectrum of a representative test case is performed using the incomplete triangular factorization. Preconditioning eliminates most of the original complex eigenvalues, and compactness is not necessarily jeopardized by a large condition number

A non-overlapping domain decomposition method for parabolic initial-boundary value problems

A non-overlapping domain decomposition method with adaptive interface conditions is applied to parabolic initial-boundary value problems in the full range from diffusion- to advection-dominated problems. The basic discretizations are the discontinuous Galerkin method in time and a stabilized Galerkin method in space. A convergence proof is available in appropriate Sobolev norms for the continuous elliptic problems arising in each time step. The numerical convergence rate is independent of the mesh size. Finally we extend the approach to more complex problems.

Parallel algorithms for semi-lagrangian advection

SUMMARY Numerical time step limitations associated with the explicit treatment of advection-dominated problems in computational fluid dynamics are often relaxed by employing Eulerian‐Lagrangian methods. These are also known as semi-Lagrangian methods in the atmospheric sciences. Such methods involve backward time integration of a characteristic equation to find the departure point of a fluid particle arriving at a Eulerian grid point. The value of the advected field at the departure point is obtained by interpolation. Both the trajectory integration and repeated interpolation influence accuracy. We compare the accuracy and performance of interpolation schemes based on piecewise cubic polynomials and cubic B-splines in the context of a distributed memory, parallel computing environment. The computational cost and interprocessor communication requirements for both methods are reported. Spline interpolation has better conservation properties but requires the solution of a global linear system, initially appearing to hinder a distributed memory implementation. The proposed parallel algorithm for multidimensional spline interpolation has almost the same communication overhead as local piecewise polynomial interpolation. We also compare various techniques for tracking trajectories given different values for the Courant number. Large Courant numbers require a high-order ODE solver involving multiple interpolations of the velocity field. # 1997 by John Wiley & Sons, Ltd.