A Hybrid Lagrangian/Eulerian Collocated Velocity Advection and Projection Method for Fluid Simulation

We present a hybrid particle/grid approach for simulating incompressible fluids on collocated velocity grids. Our approach supports both particle-based Lagrangian advection in very detailed regions of the flow and efficient Eulerian grid-based advection in other regions of the flow. A novel Backward Semi-Lagrangian method is derived to improve accuracy of grid based advection. Our approach utilizes the implicit formula associated with solutions of the inviscid Burgers’ equation. We solve this equation using Newton's method enabled by C¹ continuous grid interpolation. We enforce incompressibility over collocated, rather than staggered grids. Our projection technique is variational and designed for B-spline interpolation over regular grids where multiquadratic interpolation is used for velocity and multilinear interpolation for pressure. Despite our use of regular grids, we extend the variational technique to allow for cut-cell definition of irregular flow domains for both Dirichlet and free surface boundary conditions.     

A second-order curvilinear to Cartesian transformation of immersed interfaces and boundaries. Application to fictitious domains and multiphase flows

A global methodology dealing with fictitious domains of all kinds on curvilinear grids is presented. The main idea is to transform the curvilinear framework and its associated elements (velocity, immersed interfaces…) into a Cartesian grid. On such grids, many operations can be performed much faster than on curvilinear grids. The method is coupled with a Thread Ray-casting algorithm which works on Cartesian grids only. This algorithm computes quickly the Heaviside function related to the interior of an object on an Eulerian grid. The approach is also coupled with an immersed boundary method (L²-penalty) or with phase advection methods such as VOF–PLIC, VOF–TVD, Front-tracking or Level-set approaches. Applications, convergence and speed tests are performed for shape initializations, immersed boundary methods, and interface tracking.     

Guidelines for Poisson Solvers on Irregular Domains with Dirichlet Boundary Conditions Using the Ghost Fluid Method

We consider the variable coefficient Poisson equation with Dirichlet boundary conditions on irregular domains. We present numerical evidence for the accuracy of the solution and its gradients for different treatments at the interface using the Ghost Fluid Method for Poisson problems of Gibou et al. (J. Comput. Phys. 176:205–227, 2002; 202:577–601, 2005). This paper is therefore intended as a guide for those interested in using the GFM for Poisson-type problems (and by consequence diffusion-like problems and Stefan-type problems) by providing the pros and cons of the different choices for defining the ghost values and locating the interface. We found that in order to obtain second-order-accurate gradients, both a quadratic (or higher order) extrapolation for defining the ghost values and a quadratic (or higher order) interpolation for finding the interface location are required. In the case where the ghost values are defined by a linear extrapolation, the gradients of the solution converge slowly (at most first order in average) and the convergence rate oscillates, even when the interface location is defined by a quadratic interpolation. The same conclusions hold true for the combination of a quadratic extrapolation for the ghost cells and a linear interpolation. The solution is second-order accurate in all cases. Defining the ghost values with quadratic extrapolations leads to a non-symmetric linear system with a worse conditioning than that of the linear extrapolation case, for which the linear system is symmetric and better conditioned. We conclude that for problems where only the solution matters, the method described by Gibou, F., Fedkiw, R., Cheng, L.-T. and Kang, M. in (J. Comput. Phys. 176:205–227, 2002) is advantageous since the linear system that needs to be inverted is symmetric. In problems where the solution gradient is needed, such as in Stefan-type problems, higher order extrapolation schemes as described by Gibou, F. and Fedkiw, R. in (J. Comput. Phys. 202:577–601, 2005) are desirable.     

Boundary control,stabilization and zero-pole dynamics for a non-linear distributed parameter system

In this work we show that the now standard lumped non-linear enhancement of root-locus design still persists for a non-linear distributed parameter boundary control system governed by a scalar viscous Burgers' equation. Namely, we construct a proportional error boundary feedback control law and show that closed-loop trajectories tend to trajectories of the open-loop zero dynamics as the gain parameters are increased to infinity. We also prove a robust version of this result, valid for perturbations by an unknown disturbance with arbitrary L² norm. For the controlled Burgers' equation forced by a disturbance we prove that, for all initial data in L²(0, 1), the closed-loop trajectories converge in L²(0, 1), uniformly in t∈[0, T] and in H¹(0, 1), uniformly in t∈[t₀, T] for any t₀>0, to the trajectories of the corresponding perturbed zero dynamics. We have also extended these results to include the case when additional boundary controls are included in the closed-loop system. This provides a proof of convergence of trajectories in case the zero dynamics is replaced by a non-homogeneous Dirichlet boundary controlled Burgers' equation. As an application of our convergence of trajectories results, we demonstrate that our boundary feedback control scheme provides a semiglobal exponential stabilizing feedback law in L², H¹ and L^∞ for the open-loop system consisting of Burgers' equation with Neumann boundary conditions and zero forcing term. We also show that this result is robust in the sense that if the open-loop system is perturbed by a sufficiently small non-zero disturbance then the resulting closed-loop system is ‘practically semiglobally stabilizable’ in L²-norm. Copyright © 1999 John Wiley & Sons, Ltd.     

A pseudo-spectral scheme for the approximate solution of a time-fractional diffusion equation

We consider an initial-boundary-value problem for a time-fractional diffusion equation with initial condition u₀(x) and homogeneous Dirichlet boundary conditions in a bounded interval [0, L]. We study a semidiscrete approximation scheme based on the pseudo-spectral method on Chebyshev–Gauss–Lobatto nodes. In order to preserve the high accuracy of the spectral approximation we use an approach based on the evaluation of the Mittag-Leffler function on matrix arguments for the integration along the time variable. Some examples are presented and numerical experiments illustrate the effectiveness of the proposed approach.     

A Nonoverlapping Domain Decomposition Method for Legendre Spectral Collocation Problems

We consider the Dirichlet boundary value problem for Poisson’s equation in an L-shaped region or a rectangle with a cross-point. In both cases, we approximate the Dirichlet problem using Legendre spectral collocation, that is, polynomial collocation at the Legendre–Gauss nodes. The L-shaped region is partitioned into three nonoverlapping rectangular subregions with two interfaces and the rectangle with the cross-point is partitioned into four rectangular subregions with four interfaces. In each rectangular subregion, the approximate solution is a polynomial tensor product that satisfies Poisson’s equation at the collocation points. The approximate solution is continuous on the entire domain and its normal derivatives are continuous at the collocation points on the interfaces, but continuity of the normal derivatives across the interfaces is not guaranteed. At the cross point, we require continuity of the normal derivative in the vertical direction. The solution of the collocation problem is first reduced to finding the approximate solution on the interfaces. The discrete Steklov–Poincaré operator corresponding to the interfaces is self-adjoint and positive definite with respect to the discrete inner product associated with the collocation points on the interfaces. The approximate solution on the interfaces is computed using the preconditioned conjugate gradient method. A preconditioner is obtained from the discrete Steklov–Poincaré operators corresponding to pairs of the adjacent rectangular subregions. Once the solution of the discrete Steklov–Poincaré equation is obtained, the collocation solution in each rectangular subregion is computed using a matrix decomposition method. The total cost of the algorithm is O(N ³), where the number of unknowns is proportional to N ².      

On the accuracy assessment of Laplacian models in MPS

From the basis of the Gauss divergence theorem applied on a circular control volume that was put forward by Isshiki (2011) in deriving the MPS-based differential operators, a more general Laplacian model is further deduced from the current work which involves the proposal of an altered kernel function. The Laplacians of several functions are evaluated and the accuracies of various MPS Laplacian models in solving the Poisson equation that is subjected to both Dirichlet and Neumann boundary conditions are assessed. For regular grids, the Laplacian model with smaller N$N$ is generally more accurate, owing to the reduction of leading errors due to those higher-order derivatives appearing in the modified equation. For irregular grids, an optimal N$N$ value does exist in ensuring better global accuracy, in which this optimal value of N$N$ will increase when cases employing highly irregular grids are computed. Finally, the accuracies of these MPS Laplacian models are assessed in an incompressible flow problem.     

Convergence of an upwind control-volume mixed finite element method for convection–diffusion problems

Summary We consider a upwind control volume mixed finite element method for convection–diffusion problem on rectangular grids. These methods use the lowest order Raviart–Thomas mixed finite element space as the trial functional space and associate control-volumes, or covolumes, with the vector variable as well as the scalar variable. Chou et al. [6] established a one-half order convergence in discrete L ²-norms. In this paper, we establish a first order convergence for both the vector variable as well as the scalar variable in discrete L ²-norms.      

Flux-based level-set method for two-phase flows on unstructured grids

In this paper we deal with the application of the flux-based level set method to moving interface computations on unstructured grids. The focus lies on the overcoming of the known difficulties of level set methods, e.g. accurate computations of important geometric properties, reliable and precise reinitialization of the level set function and the adaption of standard discretization methods to the moving boundary case. The basic building block of our approach is the high-resolution flux-based level set method for general advection equation (Frolkovi? and Mikula in SIAM J Sci Comput 29(2):579–597, 2007, Frolkovi? and Wehner in Comput Vis Sci 12(6):626–650, 2009). We extend this method for the problem of reinitialization of the level set function on unstructured grids by using quadratic interpolation to compute distances for nodes close to the interface. To realize numerical simulation for some applications with moving boundaries, we adapt the approach of ghost fluid method (Gibou and Fedkiw in J Comput Phys 202:577–601, 2005) for unstructured grids. The idea is to describe the development of the moving boundary with a level set formulation while the computational grid remains fixed and the boundary conditions are enforced using some extrapolation. Our main motivation is the numerical solution of two-phase incompressible flow problems. Additionally to previously mentioned steps, we introduce further numerical schemes in the framework of finite volume discretization for the flow. Possible jumps of the pressure and the directional derivative of velocity at the interface are modeled directly within the method using the approach of extended approximation spaces. Besides that, an algorithm for the computations of curvature is considered that exhibits the second order accuracy for some examples. Numerical experiments are provided for the presented methods.     

Wavelet solution of variable order pseudodifferential equations

Sobolev spaces H ^m(x)(I) of variable order 0<m(x)<1 on an interval I⊂ℝ arise as domains of Dirichlet forms for certain quadratic, pure jump Feller processes X _t ∈ℝ with unbounded, state-dependent intensity of small jumps. For spline wavelets with complementary boundary conditions, we establish multilevel norm equivalences in H ^m(x)(I) and prove preconditioning and wavelet matrix compression results for the variable order pseudodifferential generators A of X.     

Condition Numbers of Approximate Schur Complements in Two- and Three-Dimensional Discretizations on Hierarchically Ordered Grids

We investigate multilevel incomplete factorizations of M-matrices arising from finite difference discretizations. The nonzero patterns are based on special orderings of the grid points. Hence, the Schur complements that result from block elimination of unknowns refer to a sequence of hierarchical grids. Having reached the coarsest grid, Gaussian elimination yields a complete decomposition of the last Schur complement. The main focus of this paper is a generalization of the recursive five-point/nine-point factorization method (which can be applied in two-dimensional problems) to matrices that stem from discretizations on three-dimensional cartesian grids. Moreover, we present a local analysis that considers fundamental grid cells. Our analysis allows to derive sharp bounds for the condition number associated with one factorization level (two-grid estimates). A comparison in case of the Laplace operator with Dirichlet boundary conditions shows: Estimating the relative condition number of the multilevel preconditioner by multiplying corresponding two-grid values gives the asymptotic bound O(h ^−0.347) for the two- respectively O(h ^−4/5) for the three-dimensional model problem. Received October 19, 1998; revised September 27, 1999     

Implementation of quadratic upstream interpolation schemes for solute transport into HYDRUS-1D

Numerical solution of the advection–dispersion equation, used to evaluate transport of solutes in porous media, requires discretization schemes for space and time stepping. We examine use of quadratic upstream interpolation schemes QUICK, QUICKEST, and the total variation diminution scheme ULTIMATE, and compare these with UPSTREAM and CENTRAL schemes in the HYDRUS-1D model. Results for purely convective transport show that quadratic schemes can reduce the oscillations compared to the CENTRAL scheme and numerical dispersion compared to the UPSTREAM scheme. When dispersion is introduced all schemes give similar results for Peclet number Pe < 2. All schemes show similar behavior for non-uniform grids that become finer in the direction of flow. When grids become coarser in the direction of flow, some schemes produce considerable oscillations, with all schemes showing significant clipping of the peak, but quadratic schemes extending the range of stability tenfold to Pe < 20. Similar results were also obtained for transport of a non-linear retarded solute transport (except the QUICK scheme) and for reactive transport (except the UPSTREAM scheme). Analysis of transient solute transport show that all schemes produce similar results for the position of the infiltration front for Pe = 2. When Pe = 10, the CENTRAL scheme produced significant oscillations near the infiltration front, compared to only minor oscillations for QUICKEST and no oscillations for the ULTIMATE scheme. These comparisons show that quadratic schemes have promise for extending the range of stability in numerical solutions of solute transport in porous media and allowing coarser grids.     

Optimal Error Estimates of the Legendre Tau Method for Second-Order Differential Equations

In this paper, we prove that the Legendre tau method has the optimal rate of convergence in L ²-norm, H ¹-norm and H ²-norm for one-dimensional second-order steady differential equations with three kinds of boundary conditions and in C([0,T];L ²(I))-norm for the corresponding evolution equation with the Dirichlet boundary condition. For the generalized Burgers equation, we develop a Legendre tau-Chebyshev collocation method, which can also be optimally convergent in C([0,T];L ²(I))-norm. Finally, we give some numerical examples.     

Bounded Error Schemes for the Wave Equation on Complex Domains

This paper considers the application of the method of boundary penalty terms (SAT) to the numerical solution of the wave equation on complex shapes with Dirichlet boundary conditions. A theory is developed, in a semi-discrete setting, that allows the use of a Cartesian grid on complex geometries, yet maintains the order of accuracy with only a linear temporal error-bound. A numerical example, involving the solution of Maxwell’s equations inside a 2-D circular wave-guide demonstrates the efficacy of this method in comparison to others (e.g., the staggered Yee scheme)—we achieve a decrease of two orders of magnitude in the level of the L₂-error.     

Stability of a Hot Smoluchowski Fluid

We study coupled nonlinear parabolic equations for a fluid described by a material density and a temperature , both functions of space and time. In one dimension, we find some stationary solutions corresponding to fixing the temperature on the boundary, with no-escape boundary conditions for the material. For the special case, where the temperature on the boundary is the same at both ends, the linearized equations for small perturbations about a stationary solution at uniform temperature and density are derived; they are subject to boundary conditions, Dirichlet for and no-flow conditions for the material. The spectrum of the generator L of time evolution, regarded as an operator on L ², is shown to be real, discrete and non-positive, even though L is not self-adjoint. This result is necessary for the stability of the stationary state, but might not be sufficient. The problem lies in the fact that L is not a sectorial operator, since its numerical range is C.     

A Class of Sparse Spectral Operators for Inversion of Powers of the Laplacian in N Dimensions

For inversion of the Laplacian subject to Dirichlet boundary conditions and, more generally, for the kth power of the Laplacian subject to boundary conditions on the function and its first k – 1 derivatives in the normal coordinate, there is a sparse symmetric, well-conditioned, projection of the operator that results from an expansion in associated Legendre polynomials.     

Perspective silhouette of a general swept volume

We present an efficient and robust algorithm for computing the perspective silhouette of the boundary of a general swept volume. We also construct the topology of connected components of the silhouette. At each instant t, a three-dimensional object moving along a trajectory touches the envelope surface of its swept volume along a characteristic curve K^t. The same instance of the moving object has a silhouette curve L^t on its own boundary. The intersection K^t∩L^t contributes to the silhouette of the general swept volume. We reformulate this problem as a system of two polynomial equations in three variables. The connected components of the resulting silhouette curves are constructed by detecting the instances where the two curves K^t and L^t intersect each other tangentially on the surface of the moving object. We also consider a general case where the eye position changes while moving along a predefined path. The problem is reformulated as a system of two polynomial equations in four variables, where the zero-set is a two-manifold. By analyzing the topology of the zero-set, we achieve an efficient algorithm for generating a continuous animation of perspective silhouettes of a general swept volume.     

Linear Time Algorithms for Exact Distance Transform

In 2003, Maurer et al. (IEEE Trans. Pattern Anal. Mach. Intell. 25:265–270, 2003) published a paper describing an algorithm that computes the exact distance transform in linear time (with respect to image size) for the rectangular binary images in the k-dimensional space ℝ ^k and distance measured with respect to L _p -metric for 1≤p≤∞, which includes Euclidean distance L ₂. In this paper we discuss this algorithm from theoretical and practical points of view. On the practical side, we concentrate on its Euclidean distance version, discuss the possible ways of implementing it as signed distance transform, and experimentally compare implemented algorithms. We also describe the parallelization of these algorithms and discuss the computational time savings associated with them. All these implementations will be made available as a part of the CAVASS software system developed and maintained in our group (Grevera et al. in J. Digit. Imaging 20:101–118, 2007). On the theoretical side, we prove that our version of the signed distance transform algorithm, GBDT, returns the exact value of the distance from the geometrically defined object boundary. We provide a complete proof (which was not given of Maurer et al. (IEEE Trans. Pattern Anal. Mach. Intell. 25:265–270, 2003) that all these algorithms work correctly for L _p -metric with 1<p<∞. We also point out that the precise form of the algorithm from Maurer et al. (IEEE Trans. Pattern Anal. Mach. Intell. 25:265–270, 2003) is not well defined for L ₁ and L _∞ metrics. In addition, we show that the algorithm can be used to find, in linear time, the exact value of the diameter of an object, that is, the largest possible distance between any two of its elements.     

Fast distance transformation on irregular two-dimensional grids

In this article, we propose a new fast algorithm to compute the squared Euclidean distance transform (E²DT) on every two-dimensional (2-D) irregular isothetic grids (regular square grids, quadtree based grids, etc.). Our new fast algorithm is an extension of the E²DT method proposed by Breu et al. [3]. It is based on the implicit order of the cells in the grid, and builds a partial Voronoi diagram of the centers of background cells thanks to a data structure of lists. We compare the execution time of our method with the ones of others approaches we developed in previous works. In those experiments, we consider various kinds of classical 2-D grids in imagery to show the interest of our methodology, and to point out its robustness. We also show that our method may be interesting regarding an application where we extract an adaptive medial axis construction based on a quadtree decomposition scheme.     

Null controllability of the Burgers system with distributed controls

This paper is devoted to present some positive and negative controllability results for the viscous Burgers equation. More precisely, in the context of null controllability with distributed controls, we present sharp estimates of the minimal time of controllability T(r) of initial data of L²-norm less or equal to r. In particular, we see that (global) null controllability does not hold in general (unless the control is exerted everywhere). The same results apply to similar boundary controllability systems with one boundary control.