Search

We propose an iteratively adaptive Multi-scale Finite Element Method (MsFEM) for elliptic PDEs with rough coefficients. The choice of the local boundary conditions for the multi-sale basis functions determines the accuracy of the MsFEM numerical solution, and one needs to incorporate the global information of the elliptic equation into the local boundary conditions of the multi-scale basis functions to recover the underlying fine-mesh solution of the equation. In our proposed iteratively adaptive method, we achieve this global-to-local information transfer through the combination of coarse-mesh solving using adaptive multi-scale basis functions and fine-mesh smoothing operations. In each iteration step, we first update the multi-scale basis functions based on the approximate numerical solutions of the previous iteration steps, and obtain the coarse-mesh approximate solution using a Galerkin projection. Then we apply several steps of smoothing operations to the coarse-mesh approximate solution on the underlying fine mesh to get the updated approximate numerical solution. The proposed algorithm can be viewed as a nonlinear two-level multi-grid method with the restriction and prolongation operators adapted to the approximate numerical solutions of the previous iteration steps. Convergence analysis of the proposed algorithm is carried out under the framework of two-level multi-grid method, and the harmonic coordinates are employed to establish the approximation property of the adaptive multi-scale basis functions. We demonstrate the efficiency of our proposed multi-scale methods through several numerical examples including a multi-scale coefficient problem, a high-contrast interface problem, and a convection-dominated diffusion problem.

In this paper, we discuss a general multiscale model reduction framework based on multiscale finite element methods. We give a brief overview of related multiscale methods. Due to page limitations, the overview focuses on a few related methods and is not intended to be comprehensive. We present a general adaptive multiscale model reduction framework, the Generalized Multiscale Finite Element Method. Besides the method's basic outline, we discuss some important ingredients needed for the method's success. We also discuss several applications. The proposed method allows performing local model reduction in the presence of high contrast and no scale separation.

We propose two numerical methods for accelerating the convergence of the standard fixed point method associated with a nonlinear and/or degenerate elliptic partial differential equation. The first method is linearly stable, while the second is provably convergent in the viscosity solution sense. In practice, the methods converge at a nearly linear complexity in terms of the number of iterations required for convergence. The methods are easy to implement and do not require the construction or approximation of the Jacobian. Numerical examples are shown for Bellman's equation, Isaacs' equation, Pucci's equations, the Monge---Ampere equation, a variant of the infinity Laplacian, and a system of nonlinear equations.

For viscous conservation laws, solutions contain smooth but high-contrast features, which require the use of fine grids to properly resolve. On coarse grids, these high-contrast jumps resemble shocks rather than their true viscous profiles, which could lead to issues in the numerical approximation of their underlying dynamics. In many cases, the equations of motion emit traveling wave solutions which can be used to represent the viscous profiles analytically. The traveling wave solutions can be thought of as a lower dimensional representation of the motion, since they contain information from the evolution equation, but are constant along certain time-space curves. Using a parameterized basis involving the traveling waves, along with the sparse + low-energy decompositions found in imaging sciences, we propose an approximation to viscous conservation laws which separates the coarse smooth component from the sharp fine one. Our method provides an appropriate approximation to the solution on a coarse grid, thereby accurately under-resolving the viscous profile. This is similar to the philosophy of shock capturing methods, in the sense that we want to capture the viscous front without needing to resolve the profile. Theoretical results on the consistency of our method are shown in general. We provide several computational examples for convex and non-convex fluxes.

We introduce a new concept of sparsity for the stochastic elliptic operator $-{\rm div}\left(a(x,\omega)\nabla(\cdot)\right)$, which reflects the compactness of its inverse operator in the stochastic direction and allows for spatially heterogeneous stochastic structure. This new concept of sparsity motivates a heterogeneous stochastic finite element method ({\bf HSFEM}) framework for linear elliptic equations, which discretizes the equations using the heterogeneous coupling of spatial basis with local stochastic basis to exploit the local stochastic structure of the solution space. We also provide a sampling method to construct the local stochastic basis for this framework using the randomized range finding techniques. The resulting HSFEM involves two stages and suits the multi-query setting: in the offline stage, the local stochastic structure of the solution space is identified; in the online stage, the equation can be efficiently solved for multiple forcing functions. An online error estimation and correction procedure through Monte Carlo sampling is given. Numerical results for several problems with high dimensional stochastic input are presented to demonstrate the efficiency of the HSFEM in the online stage.

This is part II of our paper in which we propose and develop a dynamically bi-orthogonal method (DyBO) to study a class of time-dependent stochastic partial differential equations (SPDEs) whose solutions enjoy a low-dimensional structure. In part I of our paper [9], we derived the DyBO formulation and proposed numerical algorithms based on this formulation. Some important theoretical results regarding consistency and bi-orthogonality preservation were also established in the first part along with a range of numerical examples to illustrate the effectiveness of the DyBO method. In this paper, we focus on the computational complexity analysis and develop an effective adaptivity strategy to add or remove modes dynamically. Our complexity analysis shows that the ratio of computational complexities between the DyBO method and a generalized polynomial chaos method (gPC) is roughly of order O((m/N"p)^3) for a quadratic nonlinear SPDE, where m is the number of mode pairs used in the DyBO method and N"p is the number of elements in the polynomial basis in gPC. The effective dimensions of the stochastic solutions have been found to be small in many applications, so we can expect m is much smaller than N"p and computational savings of our DyBO method against gPC are dramatic. The adaptive strategy plays an essential role for the DyBO method to be effective in solving some challenging problems. Another important contribution of this paper is the generalization of the DyBO formulation for a system of time-dependent SPDEs. Several numerical examples are provided to demonstrate the effectiveness of our method, including the Navier-Stokes equations and the Boussinesq approximation with Brownian forcing.

We propose a dynamically bi-orthogonal method (DyBO) to solve time dependent stochastic partial differential equations (SPDEs). The objective of our method is to exploit some intrinsic sparse structure in the stochastic solution by constructing the sparsest representation of the stochastic solution via a bi-orthogonal basis. It is well-known that the Karhunen-Loeve expansion (KLE) minimizes the total mean squared error and gives the sparsest representation of stochastic solutions. However, the computation of the KL expansion could be quite expensive since we need to form a covariance matrix and solve a large-scale eigenvalue problem. The main contribution of this paper is that we derive an equivalent system that governs the evolution of the spatial and stochastic basis in the KL expansion. Unlike other reduced model methods, our method constructs the reduced basis on-the-fly without the need to form the covariance matrix or to compute its eigendecomposition. In the first part of our paper, we introduce the derivation of the dynamically bi-orthogonal formulation for SPDEs, discuss several theoretical issues, such as the dynamic bi-orthogonality preservation and some preliminary error analysis of the DyBO method. We also give some numerical implementation details of the DyBO methods, including the representation of stochastic basis and techniques to deal with eigenvalue crossing. In the second part of our paper [11], we will present an adaptive strategy to dynamically remove or add modes, perform a detailed complexity analysis, and discuss various generalizations of this approach. An extensive range of numerical experiments will be provided in both parts to demonstrate the effectiveness of the DyBO method.

Developing an effective turbulence model is important for engineering applications as well as for fundamental understanding of the flow physics. We present a mathematical derivation of a closure relating the Reynolds stress to the mean strain rate for incompressible flows. A systematic multiscale analysis expresses the Reynolds stress in terms of the solutions of local periodic cell problems. We reveal an asymptotic structure of the Reynolds stress by invoking the frame invariant property of the cell problems and an iterative dynamic homogenization of large- and small-scale solutions. The recovery of the Smagorinsky model for homogeneous turbulence validates our derivation. Another example is the channel flow, where we derive a simplified turbulence model using the asymptotic structure near the wall. Numerical simulations at two Reynolds numbers (Re's) using our model agrees well with both experiments and Direct Numerical Simulations of turbulent channel flow.

In this paper, we investigate the performance of pseudo-spectral methods in computing nearly singular solutions of fluid dynamics equations. We consider two different ways of removing the aliasing errors in a pseudo-spectral method. The first one is the traditional 2/3 dealiasing rule. The second one is a high (36th) order Fourier smoothing which keeps a significant portion of the Fourier modes beyond the 2/3 cut-off point in the Fourier spectrum for the 2/3 dealiasing method. Both the 1D Burgers equation and the 3D incompressible Euler equations are considered. We demonstrate that the pseudo-spectral method with the high order Fourier smoothing gives a much better performance than the pseudo-spectral method with the 2/3 dealiasing rule. Moreover, we show that the high order Fourier smoothing method captures about $12 \sim 15%$ more effective Fourier modes in each dimension than the 2/3 dealiasing method. For the 3D Euler equations, the gain in the effective Fourier codes for the high order Fourier smoothing method can be as large as 20% over the 2/3 dealiasing method. Another interesting observation is that the error produced by the high order Fourier smoothing method is highly localized near the region where the solution is most singular, while the 2/3 dealiasing method tends to produce oscillations in the entire domain. The high order Fourier smoothing method is also found be very stable dynamically. No high frequency instability has been observed., Comment: 26 pages, 23 figures

In this paper, we investigate the locomotion of fish and birds by applying a new unsteady, flexible wing theory that takes into account the strong nonlinear dynamics semi-analytically. We also make extensive comparative study between the new approach and the modified vortex blob method inspired from Chorin's and Krasny's work. We first implement the modified vortex blob method for two examples and then discuss the numerical implementation of the nonlinear analytical mathematical model of Wu. We will demonstrate that Wu's method can capture the nonlinear effects very well by applying it to some specific cases and by comparing with the experiments available. In particular, we apply Wu's method to analyze Wagner's result for a wing abruptly undergoing an increase in incidence angle. Moreover, we study the vorticity generated by a wing in heaving, pitching and bending motion. In both cases, we show that the new method can accurately represent the vortex structure behind a flying wing and its influence on the bound vortex sheet on the wing.

In this paper we propose a modified multiscale finite element method for two-phase flow simulations in heterogeneous porous media. The main idea of the method is to use the global fine-scale solution at initial time to determine the boundary conditions of the basis functions. This method provides a significant improvement in two-phase flow simulations in porous media where the long-range effects are important. This is typical for some recent benchmark tests, such as the SPE comparative solution project [M. Christie, M. Blunt, Tenth spe comparative solution project: a comparison of upscaling techniques, SPE Reser. Eval. Eng. 4 (2001) 308-317], where porous media have a channelized structure. The use of global information allows us to capture the long-range effects more accurately compared to the multiscale finite element methods that use only local information to construct the basis functions. We present some analysis of the proposed method to illustrate that the method can indeed capture the long-range effect in channelized media.

The main goal of this paper is to design an efficient sampling technique for dynamic data integration using the Langevin algorithms. Based on a coarse-scale model of the problem, we compute the proposals of the Langevin algorithms using the coarse-scale gradient of the target distribution. To guarantee a correct and efficient sampling, each proposal is first tested by a Metropolis acceptance-rejection step with a coarse-scale distribution. If the proposal is accepted in the first stage, then a fine-scale simulation is performed at the second stage to determine the acceptance probability. Comparing with the direct Langevin algorithm, the new method generates a modified Markov chain by incorporating the coarse-scale information of the problem. Under some mild technical conditions we prove that the modified Markov chain converges to the correct posterior distribution. We would like to note that the coarse-scale models used in the simulations need to be inexpensive, but not necessarily very accurate, as our analysis and numerical simulations demonstrate. We present numerical examples for sampling permeability fields using two-point geostatistics. Karhunen-Loeve expansion is used to represent the realizations of the permeability field conditioned to the dynamic data, such as the production data, as well as the static data. The numerical examples show that the coarse-gradient Langevin algorithms are much faster than the direct Langevin algorithms but have similar acceptance rates.

In this paper, we propose a numerical method based on Wiener Chaos expansion and apply it to solve the stochastic Burgers and Navier-Stokes equations driven by Brownian motion. The main advantage of the Wiener Chaos approach is that it allows for the separation of random and deterministic effects in a rigorous and effective manner. The separation principle effectively reduces a stochastic equation to its associated propagator, a system of deterministic equations for the coefficients of the Wiener Chaos expansion. Simple formulas for statistical moments of the stochastic solution are presented. These formulas only involve the solutions of the propagator. We demonstrate that for short time solutions the numerical methods based on the Wiener Chaos expansion are more efficient and accurate than those based on the Monte Carlo simulations.

In this paper, we propose a pseudo-spectral multiscale method for simulating complex systems with more than one spatial scale. Using a spectral decomposition, we split the displacement into its mean and fluctuation parts. Under the assumption of localized nonlinear fluctuations, we separate the domain into an MD (Molecular Dynamics) subdomain and an MC (MacrosCopic) subdomain. An interfacial condition is proposed across the two scales, in terms of a time history treatment. In the special case of a linear system, this is the first exact interfacial condition for multiscale computations. Meanwhile, we design coarse grid equations using a matching differential operator approach. The coarse grid discretization is of spectral accuracy. We do not use a handshaking region in this method. Instead, we define a coarse grid over the whole domain and reassign the coarse grid displacement in the MD subdomain with an average of the MD solution. To reduce the computational cost, we compute the dynamics of the coarse grid displacement and relate it to the mean displacement. Our method is therefore called a pseudo-spectral multiscale method. It allows one to reach high resolution by balancing the accuracy at the coarse grid with that at the interface. Numerical experiments in one- and two-space dimensions are presented to demonstrate the accuracy and the robustness of the method.

In this paper, we present a mathematical framework of the bridging scale method (BSM), recently proposed by Liu et al. Under certain conditions, it had been designed for accurately and efficiently simulating complex dynamics with different spatial scales. From a clear and consistent derivation, we identify two error sources in this method. First, we use a linear finite element interpolation, and derive the coarse grid equations directly from Newton's second law. Numerical error in this length scale exists mainly due to inadequate approximation for the effects of the fine scale fluctuations. An modified linear element (MLE) scheme is developed to improve the accuracy. Secondly, we derive an exact multiscale interfacial condition to treat the interfaces between the molecular dynamics region ΩD and the complementary domain ΩC, using a time history kernel technique. The interfacial condition proposed in the original BSM may be regarded as a leading order approximation to the exact one (with respect to the coarsening ratio). This approximation is responsible for minor reflections across the interfaces, with a dependency on the choice of ΩD. We further illustrate the framework and analysis with linear and non-linear lattices in one-dimensional space. Copyright © 2005 John Wiley & Sons, Ltd.

We present an efficient and robust approach in the finite element framework for numerical solutions that exhibit multiscale behavior, with applications to singularly perturbed convection-diffusion problems. The first type of equation we study is the convection-dominated convection-diffusion equation, with periodic or random coefficients; the second type of equation is an elliptic equation with singularities due to discontinuous coefficients and non-smooth boundaries. In both cases, standard methods for purely hyperbolic or elliptic problems perform poorly due to sharp boundary and internal layers in the solution. We propose a framework in which the finite element basis functions are designed to capture the local small-scale behavior correctly. When the structure of the layers can be determined locally, we apply the multiscale finite element method, in which we solve the corresponding homogeneous equation on each element to capture the small scale features of the differential operator. We demonstrate the effectiveness of this method by computing the enhanced diffusivity scaling for a passive scalar in the cellular flow. We also carry out the asymptotic error analysis for its convergence rate and perform numerical experiments for verification. For a random flow with nonlocal layer structure, we use a variational principle to gain additional information in our attempt to design asymptotic basis functions. We also apply the same framework for elliptic equations with discontinuous coefficients or non-smooth boundaries. In that case, we construct local basis function near singularities using infinite element method in order to resolve extreme singularity. Numerical results on problems with various singularities confirm the efficiency and accuracy of this approach.

We study singularity formation of three-dimensional (3-D) vortex sheets without surface tension using a new approach. First, we derive a leading order approximation to the boundary integral equation governing the 3-D vortex sheet. This leading order equation captures the most singular contributions of the integral equation. By introducing an appropriate change of variables, we show that the leading order vortex sheet equation degenerates to a two-dimensional vortex sheet equation in the direction of the tangential velocity jump. This change of variables is guided by a careful analysis based on properties of certain singular integral operators, and is crucial in identifying the leading order singular behavior. Our result confirms that the tangential velocity jump is the physical driving force of the vortex sheet singularities. We also show that the singularity type of the three-dimensional problem is similar to that of the two-dimensional problem. Moreover, we introduce a model equation for 3-D vortex sheets. This model equation captures the leading order singularity structure of the full 3-D vortex sheet equation, and it can be computed efficiently using fast Fourier transform. This enables us to perform well-resolved calculations to study the generic type of 3-D vortex sheet singularities. We will provide detailed numerical results to support the analytic prediction, and to reveal the generic form of the vortex sheet singularity.

We propose a multiscale model reduction method for partial differential equations. The main purpose of this method is to derive an effective equation for multiscale problems without scale separation. An essential ingredient of our method is to decompose the harmonic coordinates into a smooth part and a highly oscillatory part so that the smooth part is invertible and the highly oscillatory part is small. Such a decomposition plays a key role in our construction of the effective equation. We show that the solution to the effective equation is in H^2, and can be approximated by a regular coarse mesh. When the multiscale problem has scale separation and a periodic structure, our method recovers the traditional homogenized equation. Furthermore, we provide error analysis for our method and show that the solution to the effective equation is close to the original multiscale solution in the H^1 norm. Numerical results are presented to demonstrate the accuracy and robustness of the proposed method for several multiscale problems without scale separation, including a problem with a high contrast coefficient.

We develop an efficient dynamically adaptive mesh generator for time-dependent problems in two or more dimensions. The mesh generator is motivated by the variational approach and is based on solving a new set of nonlinear elliptic PDEs for the mesh map. When coupled to a physical problem, the mesh map evolves with the underlying solution and maintains high adaptivity as the solution develops complicated structures and even singular behavior. The overall mesh strategy is simple to implement, avoids interpolation, and can be easily incorporated into a broad range of applications. The efficacy of the mesh is first demonstrated by two examples of blowing-up solutions to the 2-D semilinear heat equation. These examples show that the mesh can follow with high adaptivity a finite-time singularity process. The focus of applications presented here is however the baroclinic generation of vorticity in a strongly layered 2-D Boussinesq fluid, a challenging problem. The moving mesh follows effectively the flow resolving both its global features and the almost singular shear layers developed dynamically. The numerical results show the fast collapse to small scales and an exponential vorticity growth.

We present a brief review of the application of boundary integral methods in two dimensions to multicomponent fluid flows and multiphase problems in materials science. We focus on the recent development and outcomes of methods which accurately and efficiently include surface tension. In fluid flows, we examine the effects of surface tension on the Kelvin–Helmholtz and Rayleigh–Taylor instabilities in inviscid fluids, the generation of capillary waves on the free surface, and problems in Hele-Shaw flows involving pattern formation through the Saffman–Taylor instability, pattern selection, and singularity formation. In materials science, we discuss microstructure evolution in diffusional phase transformations, and the effects of the competition between surface and elastic energies on microstructure morphology. A common link between these different physical phenomena is the utility of an analysis of the appropriate equations of motion at small spatial scales to develop accurate and efficient time-stepping methods.

Morphological instabilities are common to pattern formation problems such as the non-equilibrium growth of crystals and directional solidification. Very small perturbations caused by noise originate convoluted interfacial patterns when surface tension is small. The generic mechanisms in the formation of these complex patterns are present in the simpler problem of a Hele-Shaw interface. Amid this extreme noise sensitivity, what is then the role played by small surface tension in the dynamic formation and selection of these patterns? What is the asymptotic behaviour of the interface in the limit as surface tension tends to zero? The ill-posedness of the zero-surface-tension problem and the singular nature of surface tension pose challenging difficulties in the investigation of these questions. Here, we design a novel numerical method that greatly reduces the impact of noise, and allows us to accurately capture and identify the singular contributions of extremely small surface tensions. The numerical method combines the use of a compact interface parametrization, a rescaling of the governing equations, and very high precision. Our numerical results demonstrate clearly that the zero-surface-tension limit is indeed singular. The impact of a surface-tension-induced complex singularity is revealed in detail. The singular effects of surface tension are first felt at the tip of the interface and subsequently spread around it. The numerical simulations also indicate that surface tension defines a length scale in the fingers developing in a later stage of the interface evolution.

Standard numerical methods for the Birkhoff-Rott equation for a vortex sheet are unstable due to the amplification of roundoff error by the Kelvin-Helmholtz instability. A nonlinear filtering method was used by Krasny to eliminate this spurious growth of round-off error and accurately compute the Birkhoff-Rott solution essentially up to the time it becomes singular. In this paper convergence is proved for the discretized Birkhoff-Rott equation with Krasny filtering and simulated roundoff error. The convergence is proved for a time almost up to the singularity time of the continuous solution. The proof is in an analytic function class and uses a discrete form of the abstract Cauchy-Kowalewski theorem. In order for the proof to work almost up to the singularity time, the linear and nonlinear parts of the equation, as well as the effects of Krasny filtering, are precisely estimated. The technique of proof applies directly to other ill-posed problems such as Rayleigh-Taylor unstable interfaces in incompressible, inviscid, and irrotational fluids, as well as to Saffman-Taylor unstable interfaces in Hele-Shaw cells.

A computational method is proposed for the dynamics of solids capable of twinning and phase transitions. In a two-dimensional, sharp-interface model of twinning, the stored-energy function is a nonconvex potential with multiple wells. The evolution of twin interfaces is governed by field equations and jump conditions of momentum balance, and by a kinetic relation expressing the interface velocity as a function of the local driving traction and interfacial orientation. A regularized version of the model is constructed based on the level-set method. A level-set function which changes signs across the interface is introduced. The evolution of this function is described by a Hamilton?Jacobi equation, whose velocity coefficient is determined by the kinetic relation. Jump conditions are thereby eliminated, allowing finite-difference discretization. Numerical simulations exhibit complex evolution of the interface, including cusp formation, needle growth, spontaneous tip splitting, and topological changes that result in microstructure refinement. The results capture experimentally observed phenomena in martensitic crystals.

We propose a multiscale finite element method for solving second order elliptic equations with rapidly oscillating coefficients. The main purpose is to design a numerical method which is capable of correctly capturing the large scale components of the solution on a coarse grid without accurately resolving all the small scale features in the solution. This is accomplished by incorporating the local microstructures of the differential operator into the finite element base functions. As a consequence, the base functions are adapted to the local properties of the differential operator. In this paper, we provide a detailed convergence analysis of our method under the assumption that the oscillating coefficient is of two scales and is periodic in the fast scale. While such a simplifying assumption is not required by our method, it allows us to use homogenization theory to obtain a useful asymptotic solution structure. The issue of boundary conditions for the base functions will be discussed. Our numerical experiments demonstrate convincingly that our multiscale method indeed converges to the correct solution, independently of the small scale in the homogenization limit. Application of our method to problems with continuous scales is also considered.

In this paper, we present a new formulation for computing the motion of a curvature-driven 3-D filament. This new numerical method has none of the high order time step stability constraints that are usually associated with curvature regularization. This result generalizes the previous work of Houet al. (1994) for 2-D fluid interfaces with surface tension. Applications to 2-D vortex sheets, 3-D motion by curvature, the Kirchhoff rod model, and nearly anti-parallel vortex filaments will be presented to demonstrate the robustness of the method.

Boundary integral methods to simulate interfacial flows are very sensitive to numerical instabilities. In addition, surface tension introduces nonlinear terms with high order spatial derivatives into the interface dynamics. This makes the spatial discretization even more difficult and, at the same time, imposes a severe time step constraint for stable explicit time integration methods. A proof of the convergence of a reformulated boundary integral method for two-density fluid interfaces with surface tension is presented. The method is based on a scheme introduced by Hou, Lowengrub and Shelley [ J. Comp. Phys. 114 (1994), pp. 312–338] to remove the high order stability constraint or stiffness. Some numerical filtering is applied carefully at certain places in the discretization to guarantee stability. The key of the proof is to identify the most singular terms of the method and to show, through energy estimates, that these terms balance one another. The analysis is at a time continuous-space discrete level but a fully discrete case for a simple Hele-Shaw interface is also studied. The time discrete analysis shows that the high order stiffness is removed and also provides an estimate of how the CFL constraint depends on the curvature and regularity of the solution. The robustness of the method is illustrated with several numerical examples. A numerical simulation of an unstably stratified two-density interfacial flow shows the roll-up of the interface; the computations proceed up to a time where the interface is about to pinch off and trapped bubbles of fluid are formed. The method remains stable even in the full nonlinear regime of motion. Another application of the method shows the process of drop formation in a falling single fluid.

In this paper, a hybrid approach which combines theimmersed interface methodwith thelevel set approachis presented. The fast version of the immersed interface method is used to solve the differential equations whose solutions and their derivatives may be discontinuous across the interfaces due to the discontinuity of the coefficients or/and singular sources along the interfaces. The moving interfaces then are updated using the newly developed fast level set formulation which involves computation only inside some small tubes containing the interfaces. This method combines the advantage of the two approaches and gives a second-order Eulerian discretization for interface problems. Several key steps in the implementation are addressed in detail. This new approach is then applied to Hele?Shaw flow, an unstable flow involving two fluids with very different viscosity.

In this paper, we study a multiscale finite element method for solving a class of elliptic problems arising from composite materials and flows in porous media, which contain many spatial scales. The method is designed to efficiently capture the large scale behavior of the solution without resolving all the small scale features. This is accomplished by constructing the multiscale finite element base functions that are adaptive to the local property of the differential operator. Our method is applicable to general multiple-scale problems without restrictive assumptions. The construction of the base functions is fully decoupled from element to element; thus, the method is perfectly parallel and is naturally adapted to massively parallel computers. For the same reason, the method has the ability to handle extremely large degrees of freedom due to highly heterogeneous media, which are intractable by conventional finite element (difference) methods. In contrast to some empirical numerical upscaling methods, the multiscale method is systematic and self- consistent, which makes it easier to analyze. We give a brief analysis of the method, with emphasis on the “resonant sampling” effect. Then, we propose an oversampling technique to remove the resonance effect. We demonstrate the accuracy and efficiency of our method through extensive numerical experiments, which include problems with random coefficients and problems with continuous scales. Parallel implementation and performance of the method are also addressed.

We prove nonlinear stability and convergence of certain boundary integral methods for time-dependent water waves in a two-dimensional, inviscid, irrotational, incompressible fluid, with or without surface tension. The methods are convergent as long as the underlying solution remains fairly regular (and a sign condition holds in the case without surface tension). Thus, numerical instabilities are ruled out even in a fully nonlinear regime. The analysis is based on delicate energy estimates, following a framework previously developed in the continuous case [Beale, Hou, and Lowengrub, Comm. Pure Appl. Math., 46 (1993), pp. 1269--1301]. No analyticity assumption is made for the physical solution. Our study indicates that the numerical methods must satisfy certain compatibility conditions in order to be stable. Violation of these conditions will lead to numerical instabilities. A breaking wave is calculated as an illustration.

We use a boundary integral technique to study the two space dimensional Mullins?Sekerka free boundary problem which originates from a study of solidification and liquidation of materials of negligible specific heat. This is an area preserving and curve shortening motion. Evolution equations for the free boundaries are derived in terms of the tangent angle and total arclength, which makes a small scale decomposition possible and the Fourier transform a powerful tool in numerical calculations. With this formulation, implicit schemes can be implemented to avoid the difficult numerical stiffness associated with explicit schemes. We can compute solutions up to the time when there is a topological change, i.e., when particles touch or break up. Our numerical results for systems of a single particle or multi-particles provide some valuable information in the particle dynamics, such as the circularization of each individual particle, and the mass transfer between different particles during particle interactions.

This paper considers numerical methods for computing propagating phase boundaries in solids described by the physical model introduced by Abeyaratne and Knowles. The model under consideration consists of a set of conservation laws supplemented with a kinetic relation and a nucleation criterion. Discontinuities between two different phases are undercompressive crossing waves in the general terminology of nonstrictly hyperbolic systems of conservation laws. This paper studies numerical methods designed for the computation of such crossing waves. We propose a Godunov-type method combining front tracking with a capturing method; we also consider Glimm's random choice scheme. Both methods share the property that the phase boundaries are sharply computed in the sense that there are no numerical interior points for the description of a phase boundary. This property is well known for the Glimm's scheme; on the other hand, our front tracking algorithm is designed so that it tracks phase boundaries but captures shock waves. Phase boundaries are sensitive to numerical dissipation effects, so the above property is essential to ensure convergence toward the correct entropy weak solution. Convergence of the Godnuov-type method is demonstrated numerically. Extensive numerical experiments show the practical interest of both approaches for computations of undercompressive crossing waves.

A level set formulation is derived for incompressible, immiscible Navier?Stokes equations separated by a free surface. The interface is identified as the zero level set of a smooth function. Eulerian finite difference methods based on this level set formulation are proposed. These methods are robust and efficient and are capable of computing interface singularities such as merging and reconnection. Numerical experiments are presented to demonstrate the effectiveness of the methods.

In this paper, we propose a general approach called Generalized Multiscale Finite Element Method (GMsFEM) for performing multiscale simulations for problems without scale separation over a complex input space. As in multiscale finite element methods (MsFEMs), the main idea of the proposed approach is to construct a small dimensional local solution space that can be used to generate an efficient and accurate approximation to the multiscale solution with a potentially high dimensional input parameter space. In the proposed approach, we present a general procedure to construct the offline space that is used for a systematic enrichment of the coarse solution space in the online stage. The enrichment in the online stage is performed based on a spectral decomposition of the offline space. In the online stage, for any input parameter, a multiscale space is constructed to solve the global problem on a coarse grid. The online space is constructed via a spectral decomposition of the offline space and by choosing the eigenvectors corresponding to the largest eigenvalues. The computational saving is due to the fact that the construction of the online multiscale space for any input parameter is fast and this space can be re-used for solving the forward problem with any forcing and boundary condition. Compared with the other approaches where global snapshots are used, the local approach that we present in this paper allows us to eliminate unnecessary degrees of freedom on a coarse-grid level. We present various examples in the paper and some numerical results to demonstrate the effectiveness of our method., Comment: Revised version

Many physically interesting problems involve propagation of free surfaces. Vortex-sheet roll-up in hydrodynamic instability, wave interactions on the ocean's free surface, the solidification problem for crystal growth and Hele-Shaw cells for pattern formation are some of the significant examples. These problems present a great challenge to physicists and applied mathematicians because the underlying problem is very singular. The physical solution is sensitive to small perturbations. Naïve discretisations may lead to numerical instabilities. Other numerical difficulties include singularity formation and possible change of topology in the moving free surfaces, and the severe time-stepping stability constraint due to the stiffness of high-order regularisation effects, such as surface tension.This paper reviews some of the recent advances in developing stable and efficient numerical algorithms for solving free boundary-value problems arising from fluid dynamics and materials science. In particular, we will consider boundary integral methods and the level-set approach for water waves, general multi-fluid interfaces, Hele–Shaw cells, crystal growth and solidification. We will also consider the stabilising effect of surface tension and curvature regularisation. The issue of numerical stability and convergence will be discussed, and the related theoretical results for the continuum equations will be addressed. This paper is not intended to be a detailed survey and the discussion is limited by both the taste and expertise of the author.

The 91st London Mathematical Society Durham Symposium took place from July 5th to 15th 2010, with more than 100 international participants attending. The Symposium focused on Numerical Analysis of Multiscale Problems and this book contains 10 invited articles from some of the meeting's key speakers, covering a range of topics of contemporary interest in this area. Articles cover the analysis of forward and inverse PDE problems in heterogeneous media, high-frequency wave propagation, atomistic-continuum modeling and high-dimensional problems arising in modeling uncertainty. Novel upscaling and preconditioning techniques, as well as applications to turbulent multi-phase flow, and to problems of current interest in materials science are all addressed. As such this book presents the current state-of-the-art in the numerical analysis of multiscale problems and will be of interest to both practitioners and mathematicians working in those fields.

A new formulation and new methods are presented for computing the motion of fluid interfaces with surface tension in two-dimensional, irrotational, and incompressible fluids. Through the Laplace-Young condition at the interface, surface tension introduces high-order terms, both nonlinear and nonlocal, into the dynamics. This leads to severe stability constraints for explicit time integration methods and makes the application of implicit methods difficult. This new formulation has all the nice properties for time integration methods that are associated with having a linear, constant coefficient, highest order term. That is, using this formulation, we give implicit time integration methods that have no high order time step stability constraint associated with surface tension and are explicit in Fourier space. The approach is based on a boundary integral formulation and applies more generally, even to problems beyond the fluid mechanical context. Here they are applied to computing with high resolution the motion of interfaces in Hele-Shaw flows and the motion of free surfaces in inviscid flows governed by the Euler equations. One Hele-Shaw computation shows the behavior of an expanding gas bubble over long-time as the interface undergoes successive tip-splittings and finger competition. A second computation shows the formation of a very ramified interface through the interaction of surface tension with an unstable density stratification. In Euler flows, the computation of a vortex sheet shows its roll-up through the Kelvin-Helmholtz instability. This motion culminates in the late time self-intersection of the interface, creating trapped bubbles of fluid. This is, we believe, a type of singularity formation previously unobserved for such flows in 2D. Finally, computations of falling plumes in an unstably stratified Boussinesq fluid show a very similar behavior.

A rigorous convergence result is given for a projection scheme for the Navies–Stokes equations in the presence of boundaries. The numerical scheme is based on a finite-difference approximation, and the pressure is chosen so that the computed velocity satisfies a discrete divergence-free condition. This choice for the pressure and the particular way that the discrete divergence is calculated near the boundary permit the error in the pressure to be controlled and the second-order convergence in the velocity and the pressure to the exact solution to be shown. Some simplifications in the calculation of the pressure in the case without boundaries are also discussed.

We introduce a new multiscale finite element method which is able to accurately capture solutions of elliptic interface problems with high contrast coefficients by using only coarse quasiuniform meshes, and without resolving the interfaces. A typical application would be the modelling of flow in a porous medium containing a number of inclusions of low (or high) permeability embedded in a matrix of high (respectively low) permeability. Our method is H 1 H^1 - conforming, with degrees of freedom at the nodes of a triangular mesh and requiring the solution of subgrid problems for the basis functions on elements which straddle the coefficient interface but which use standard linear approximation otherwise. A key point is the introduction of novel coefficient-dependent boundary conditions for the subgrid problems. Under moderate assumptions, we prove that our methods have (optimal) convergence rate of O ( h ) O(h) in the energy norm and O ( h 2 ) O(h^2) in the L 2 L_2 norm where h h is the (coarse) mesh diameter and the hidden constants in these estimates are independent of the “contrast” (i.e. ratio of largest to smallest value) of the PDE coefficient. For standard elements the best estimate in the energy norm would be O ( h 1 / 2 − ε ) \mathcal {O}(h^{1/2-\varepsilon }) with a hidden constant which in general depends on the contrast. The new interior boundary conditions depend not only on the contrast of the coefficients, but also on the angles of intersection of the interface with the element edges.

A rigorous convergence result is presented for a finite difference scheme for the Navier-Stokes equations which uses vorticity boundary conditions. The approximating scheme is based on the vorticity-stream function formulation of the Navier-Stokes equations. The no-slip boundary condition is satisfied approximately by using a boundary condition of vorticity creation type. Convergence with second-order accuracy in vorticity and velocity is established for general domains in two space dimensions. Generalization to three space dimensions is also considered.

Whether the 3D incompressible Euler and Navier–Stokes equations can develop a finite-time singularity from smooth initial data with finite energy has been one of the most long-standing open questions. We review some recent theoretical and computational studies which show that there is a subtle dynamic depletion of nonlinear vortex stretching due to local geometric regularity of vortex filaments. We also investigate the dynamic stability of the 3D Navier–Stokes equations and the stabilizing effect of convection. A unique feature of our approach is the interplay between computation and analysis. Guided by our local non-blow-up theory, we have performed large-scale computations of the 3D Euler equations using a novel pseudo-spectral method on some of the most promising blow-up candidates. Our results show that there is tremendous dynamic depletion of vortex stretching. Moreover, we observe that the support of maximum vorticity becomes severely flattened as the maximum vorticity increases and the direction of the vortex filaments near the support of maximum vorticity is very regular. Our numerical observations in turn provide valuable insight, which leads to further theoretical breakthrough. Finally, we present a new class of solutions for the 3D Euler and Navier–Stokes equations, which exhibit very interesting dynamic growth properties. By exploiting the special nonlinear structure of the equations, we prove nonlinear stability and the global regularity of this class of solutions.

Convergence of the grid-free point vortex method is proved for three-dimensional Euler equations with smooth solutions. Two new techniques are used to obtain consistency and stability of the method. The first one is Strang’s trick [Numer. Math., 6 (1964), pp. 37–46] which allows smooth approximate solutions to be constructed to the vortex method equations with arbitrarily small errors. This result is used to obtain consistency and nonlinear stability. The second tool is the use of a very special discrete negative norm in $l_1$ space for vorticity which gives rise to the linear stability result. Combining these two techniques proves uniform convergence of the method with second-order accuracy.

For general second-order elliptic partial differential equations, the Schwarz alternating procedure is proved to converge at a rate independent of the aspect ratio for L-shaped, T-shaped, and C-shaped domains. These results cover both continuous and discrete versions of the Schwarz algorithm. Moreover, they apply to the nonoverlapping Schur complement algorithms with the preconditioner proposed in Chan [SIAM J. Numer. Anal., 24 (1987), pp. 382–390]. In particular, it is shown that the condition number of the preconditioned interface operator is bounded by 2 for all L-shaped and T-shaped domains. This improves similar geometry-independent convergence results for the Schur complement algorithms obtained previously by Chan and Resasco [“Advances in Computer Methods for Partial Differential Equations, VI,” R. Vichnevetsky and R. S. Stepleman, eds., pp. 317–322].

Based on the observation that a point vortex approximation can be made spectrally accurate by using Van de Vooren’s desingularization, short time convergence of the point vortex method for both vortex sheets and the Boussinesq approximation are proved using analytic data. The spectral accuracy of the method allows a very simple proof to be obtained without using the Cauchy–Kowalewski theorem.

In this paper, we perform a careful numerical study of nearly singular solutions of the 3D incompressible Euler equations with smooth initial data. We consider the interaction of two perturbed antiparallel vortex tubes which was previously investigated by Kerr in [16, 19]. In our numerical study, we use both the pseudo-spectral method with the 2/3 dealiasing rule and the pseudo-spectral method with a high order Fourier smoothing. Moreover, we perform a careful resolution study with grid points as large as 1,536 × 1,024 × 3,072 to demonstrate the convergence of both numerical methods. Our computational results show that the maximum vorticity does not grow faster than doubly exponential in time while the velocity field remains bounded up to T = 19, beyond the singularity time T = 18.7 reported by Kerr in [16, 19]. The local geometric regularity of vortex lines near the region of maximum vorticity seems to play an important role in depleting the nonlinear vortex stretching dynamically.

The immersed boundary method is one of the most useful computational methods in studying fluid structure interaction. On the other hand, the Immersed Boundary method is also known to require small time steps to maintain stability when solved with an explicit method. Many implicit or approximately implicit methods have been proposed in the literature to remove this severe time step stability constraint, but none of them give satisfactory performance. In this paper, we propose an efficient semi-implicit scheme to remove this stiffness from the immersed boundary method for the Navier-Stokes equations. The construction of our semi-implicit scheme consists of two steps. First, we obtain a semi-implicit discretization which is proved to be unconditionally stable. This unconditionally stable semi-implicit scheme is still quite expensive to implement in practice. Next, we apply the small scale decomposition to the unconditionally stable semi-implicit scheme to construct our efficient semi-implicit scheme. Unlike other implicit or semi-implicit schemes proposed in the literature, our semi-implicit scheme can be solved explicitly in the spectral space. Thus the computational cost of our semi-implicit schemes is comparable to that of an explicit scheme. Our extensive numerical experiments show that our semi-implicit scheme has much better stability property than an explicit scheme. This offers a substantial computational saving in using the immersed boundary method.

The Immersed Boundary method has evolved into one of the most useful computational methods in studying fluid structure interaction. On the other hand, the Immersed Boundary method is also known to suffer from a severe timestep stability restriction when using an explicit time discretization. In this paper, we propose several efficient semi-implicit schemes to remove this stiffness from the Immersed Boundary method for the two-dimensional Stokes flow. First, we obtain a novel unconditionally stable semi-implicit discretization for the immersed boundary problem. Using this unconditionally stable discretization as a building block, we derive several efficient semi-implicit schemes for the immersed boundary problem by applying the Small Scale Decomposition to this unconditionally stable discretization. Our stability analysis and extensive numerical experiments show that our semi-implicit schemes offer much better stability property than the explicit scheme. Unlike other implicit or semi-implicit schemes proposed in the literature, our semi-implicit schemes can be solved explicitly in the spectral space. Thus the computational cost of our semi-implicit schemes is comparable to that of an explicit scheme, but with a much better stability property., Comment: 40 pages with 8 figures

Convergence is proved of the variable blob vortex method proposed [T. Hou, D. Papageorgiou, G. Papanicolaou, and J. Zhu, An adaptive vortex method for two-dimensional, viscous, and incompressible flows, preprint] for the incompressible Euler and Navier–Stokes equations without boundaries. This version of the variable blob method is designed to approximate directly the Navier–Stokes equations without matching the interior flow to the boundary layer solution. By choosing vortex blobs adaptively, the numerical errors due to local shear are minimized.

We prove consistency, stability and convergence of the point vortex approximation to the 2-D incompressible Euler equations with smooth solutions. We first show that the discretization error is second-order accurate. Then we show that the method is stable in l p norm. Consequently the method converge in l p norm for all time. The convergence is also illustrated by a numerical experiment

In this paper, we discuss some applications of multiscale finite element methods to two-phase immiscible flow simulations in heterogeneous porous media. We discuss some extensions of multiscale finite element methods which take into account limited global information. These methods are well suited for channelized porous media, where the long-range effects are important. This is typical for some recent benchmark tests, such as the SPE comparative solution project [M. Christie, M. Blunt, Tenth SPE comparative solution project: A comparison of upscaling techniques, SPE Reser. Eval. Engrg. 4 (2001) 308-317], where porous media has a channelized structure. The applications of multiscale finite element methods to inverse problems arisen in subsurface characterization are also discussed in the paper.