Stabilized Moving Finite Elements for Convection Dominated Problems

The Moving Finite Element method (MFE) when applied to purely hyperbolic problems tends to move its nodes with the flow (often a good thing). But for steady or near steady problems the nodes flow past stationary regions of critical interest and pile up at the outflow. We report on efforts to develop moveable node versions of the “stabilized” finite element methods which have so successfully improved upon Galerkin in the fixed node setting. One method in particular (Galerkin-Δx TLSMFE with a “fix”) yields very promising results on our simple 1-D model problem. Its nodes lock onto and resolve sharp stationary features but also lock onto and move with the moving features of the solution.     

A Generic Stabilization Approach for Higher Order Discontinuous Galerkin Methods for Convection Dominated Problems

In this paper we present a stabilized Discontinuous Galerkin (DG) method for hyperbolic and convection dominated problems. The presented scheme can be used in several space dimension and with a wide range of grid types. The stabilization method preserves the locality of the DG method and therefore allows to apply the same parallelization techniques used for the underlying DG method. As an example problem we consider the Euler equations of gas dynamics for an ideal gas. We demonstrate the stability and accuracy of our method through the detailed study of several test cases in two space dimension on both unstructured and cartesian grids. We show that our stabilization approach preserves the advantages of the DG method in regions where stabilization is not necessary. Furthermore, we give an outlook to adaptive and parallel calculations in 3d.     

Efficient Energy-Preserving Exponential Integrators for Multi-component Hamiltonian Systems

Graph-Laplacians and their spectral embeddings play an important role in multiple areas of machine learning. This paper is focused on graph-Laplacian dimension reduction for the spectral clustering of data as a primary application, however, it can also be applied in data mining, data manifold learning, etc. Spectral embedding provides a low-dimensional parametrization of the data manifold which makes the subsequent task (e.g., clustering with k-means or any of its approximations) much easier. However, despite reducing the dimensionality of data, the overall computational cost may still be prohibitive for large data sets due to two factors. First, computing the partial eigendecomposition of the graph-Laplacian typically requires a large Krylov subspace. Second, after the spectral embedding is complete, one still has to operate with the same number of data points, which may ruin the efficiency of the approach. For example, clustering of the embedded data is typically performed with various relaxations of k-means which computational cost scales poorly with respect to the size of data set. Also, they become prone to getting stuck in local minima, so their robustness depends on the choice of initial guess. In this work, we switch the focus from the entire data set to a subset of graph vertices (target subset). We develop two novel algorithms for such low-dimensional representation of the original graph that preserves important global distances between the nodes of the target subset. In particular, it allows to ensure that target subset clustering is consistent with the spectral clustering of the full data set if one would perform such. That is achieved by a properly parametrized reduced-order model (ROM) of the graph-Laplacian that approximates accurately the diffusion transfer function of the original graph for inputs and outputs restricted to the target subset. Working with a small target subset reduces greatly the required dimension of Krylov subspace and allows to exploit the conventional algorithms (like approximations of k-means) in the regimes when they are most robust and efficient. This was verified in the numerical clustering experiments with both synthetic and real data. We also note that our ROM approach can be applied in a purely transfer-function-data-driven way, so it becomes the only feasible option for extremely large graphs that are not directly accessible. There are several uses for our algorithms. First, they can be employed on their own for representative subset clustering in cases when handling the full graph is either infeasible or simply not required. Second, they may be used for quality control. Third, as they drastically reduce the problem size, they enable the application of more sophisticated algorithms for the task under consideration (like more powerful approximations of k-means based on semi-definite programming (SDP) instead of the conventional Lloyd’s algorithm). Finally, they can be used as building blocks of a multi-level divide-and-conquer type algorithm to handle the full graph. The latter will be reported in a separate article.

     

An Accurate Polynomial Approximation of Exponential Integrators

Numerical time propagation of semi-linear equations such as reaction–diffusion, non-linear Schrödinger or semi-linear wave equations can be performed by the use of exponential time differencing. However, the evaluation of exponential integrators poses a serious technical complexity, particularly in multiple dimensions. In this paper we approach this difficulty by deriving simple polynomial series approximations of exponential integrators. Several numerical examples are presented.     

二维非线性对流扩散方程求解程序的测试与优化

在IA-64架构Itanium2处理器上,应用gprof和pfmon对二维非线性对流扩散方程求解程序源代码进行了性能测试.在分析给定程序的数据结构,子过程调用关系,重点子程序中循环体的迭代空间、数据空间、访同轨迹,输入输出数据量大小和程序结构等的基础上,应用子过程合并、循环变换、分支消除、循环顺序逆转、数组一维结构化为二维结构、输入参数给定等方法,改善了数据访问的时空局部性,程序性能有15%的提高.     

Error Analysis for hp-FEM Semi-Lagrangian Second Order BDF Method for Convection-Dominated Diffusion Problems

We present in this paper an analysis of a semi-Lagrangian second order Backward Difference Formula combined with hp-finite element method to calculate the numerical solution of convection diffusion equations in ℝ². Using mesh dependent norms, we prove that the a priori error estimate has two components: one corresponds to the approximation of the exact solution along the characteristic curves, which is O(Dt²+h^m+1(1+\frac\mathopen|logh|Dt))O(\Delta t^{2}+h^{m+1}(1+\frac{\mathopen{|}\log h|}{\Delta t})); and the second, which is O(Dt^p+|| [(u)\vec]-[(u)\vec]_h||_L^￥)O(\Delta t^{p}+\| \vec{u}-\vec{u}_{h}\|_{L^{\infty}}), represents the error committed in the calculation of the characteristic curves. Here, m is the degree of the polynomials in the finite element space, [(u)\vec]\vec{u} is the velocity vector, [(u)\vec]_h\vec{u}_{h} is the finite element approximation of [(u)\vec]\vec{u} and p denotes the order of the method employed to calculate the characteristics curves. Numerical examples support the validity of our estimates.     

Field-Aligned Interpolation for Semi-Lagrangian Gyrokinetic Simulations

This work is devoted to the study of field-aligned interpolation in semi-Lagrangian codes. In the context of numerical simulations of magnetic fusion devices, this approach is motivated by the observation that gradients of the solution along the magnetic field lines are typically much smaller than along a perpendicular direction. In toroidal geometry, field-aligned interpolation consists of a 1D interpolation along the field line, combined with 2D interpolations on the poloidal planes (at the intersections with the field line). A theoretical justification of the method is provided in the simplified context of constant advection on a 2D periodic domain: unconditional stability is proven, and error estimates are given which highlight the advantages of field-aligned interpolation. The same methodology is successfully applied to the solution of the gyrokinetic Vlasov equation, for which we present the ion temperature gradient (ITG) instability as a classical test-case: first we solve this in cylindrical geometry (screw-pinch), and next in toroidal geometry (circular Tokamak). In the first case, the algorithm is implemented in Selalib (semi-Lagrangian library), and the numerical simulations provide linear growth rates that are in accordance with the linear dispersion analysis. In the second case, the algorithm is implemented in the Gysela code, and the numerical simulations are benchmarked with those employing the standard (not aligned) scheme. Numerical experiments show that field-aligned interpolation leads to considerable memory savings for the same level of accuracy; substantial savings are also expected in reactor-scale simulations.     

ℋ-matrices for Convection-diffusion Problems with Constant Convection

L ^∞-coefficients. This paper analyses the application of hierarchical matrices to the convection-dominant convection-diffusion equation with constant convection. In the case of increasing convection, the convergence of a standard ℋ-matrix approximant towards the original matrix will deteriorate. We derive a modified partitioning and admissibility condition that ensures good convergence also for the singularly perturbed case. Received January 1, 2003; revised March 4, 2003 Published online: May 2, 2003     

Error of Partitioned Runge-Kutta Methods for Multiple Stiff Singular Perturbation Problems

The main purpose of this paper is to deal with error behaviour of partitioned Runge-Kutta methods for one-parameter multiple stiff singular perturbation problems whose stiffness is caused by a small parameter ε and some other factors and to present some quantitative convergence results. Received November 30, 1998; revised August 10, 1999     

Stabilization of the Spectral Element Method in Convection Dominated Flows by Recovery of Skew-Symmetry

We investigate stability properties of the spectral element method for advection dominated incompressible flows. In particular, properties of the widely used convective form of the nonlinear term are studied. We remark that problems which are usually associated with the nonlinearity of the governing Navier–Stokes equations also arise in linear scalar transport problems, which implicates advection rather than nonlinearity as a source of difficulty. Thus, errors arising from insufficient quadrature of the convective term, commonly referred to as ‘aliasing errors’, destroy the skew-symmetric properties of the convection operator. Recovery of skew-symmetry can be efficiently achieved by the use of over-integration. Moreover, we demonstrate that the stability problems are not simply connected to underresolution. We combine theory with analysis of the linear advection-diffusion equation in 2D and simulations of the incompressible Navier–Stokes equations in 2D of thin shear layers at a very high Reynolds number and in 3D of turbulent and transitional channel flow at moderate Reynolds number. For the Navier–Stokes equations, where the divergence-free constraint needs to be enforced iteratively to a certain accuracy, small divergence errors can be detrimental to the stability of the method and it is therefore advised to use additional stabilization (e.g. so-called filter-based stabilization, spectral vanishing viscosity or entropy viscosity) in order to assure a stable spectral element method.     

A High-Order Local Projection Stabilization Method for Natural Convection Problems

In this paper, we propose a local projection stabilization (LPS) finite element method applied to numerically solve natural convection problems. This method replaces the projection-stabilized structure of standard LPS methods by an interpolation-stabilized structure, which only acts on the high frequencies components of the flow. This approach gives rise to a method which may be cast in the variational multi-scale framework, and constitutes a low-cost, accurate solver (of optimal error order) for incompressible flows, despite being only weakly consistent. Numerical simulations and results for the buoyancy-driven airflow in a square cavity with differentially heated side walls at high Rayleigh numbers (up to \(Ra=10^7\)) are given and compared with benchmark solutions. Good accuracy is obtained with relatively coarse grids.     

Some Modified Runge-Kutta Methods for the Numerical Solution of Initial-Value Problems with Oscillating Solutions

Two new modified Runge-Kutta methods with minimal phase-lag are developed for the numerical solution of initial-value problems with oscillating solutions which can be analyzed to a system of first order ordinary differential equations. These methods are based on the well known Runge-Kutta RK5(4)7FEq1 method of Higham and Hall (1990) of order five. Also, based on the property of the phase-lag a new error control procedure is introduced. Numerical and theoretical results show that this new approach is more efficient compared with the well known Runge-Kutta Dormand-Prince RK5(4)7S method [see Dormand and Prince (1980)] and the well known Runge-Kutta RK5(4)7FEq1 method of Higham and Hall (1990).     

High-Order Time-Accurate Parallel Schemes for Parabolic Singularly Perturbed Problems with Convection

The first boundary value problem for a singularly perturbed parabolic equation of convection-diffusion type on an interval is studied. For the approximation of the boundary value problem we use earlier developed finite difference schemes, ɛ-uniformly of a high order of accuracy with respect to time, based on defect correction. New in this paper is the introduction of a partitioning of the domain for these ɛ-uniform schemes. We determine the conditions under which the difference schemes, applied independently on subdomains may accelerate (ɛ-uniformly) the solution of the boundary value problem without losing the accuracy of the original schemes. Hence, the simultaneous solution on subdomains can in principle be used for parallelization of the computational method. Received December 3, 1999; revised April 20, 2000     

Some Computer Assisted Proofs for Solutions of the Heat Convection Problems

This is a continuation of our previous results (Y. Watanabe, N. Yamamoto, T. Nakao, and T. Nishida, A Numerical Verification of Nontrivial Solutions for the Heat Convection Problem, to appear in the Journal of Mathematical Fluid Mechanics). In that work, the authors considered two-dimensional Rayleigh-Bénard convection and proposed an approach to prove existence of steady-state solutions based on an infinite dimensional fixed-point theorem using a Newton-like operator with spectral approximation and constructive error estimates. We numerically verified several exact nontrivial solutions which correspond to solutions bifurcating from the trivial solution. This paper shows more detailed results of verification for given Prandtl and Rayleigh numbers. In particular, we found a new and interesting solution branch which was not obtained in the previous study, and it should enable us to present important information to clarify the global bifurcation structure. All numerical examples discussed are take into account of the effects of rounding errors in the floating point computations.     

一种基于半拉格朗日的液体实时仿真方法

近些年,在计算机图形学与虚拟现实技术领域中,自然现象的模拟得到了广泛的关注和研究.如何快速且逼真地模拟自然现象,是此类研究的目的.以液体表面作为研究对象,总结了关于液体模拟近年来的部分研究成果;针对三维液体的复杂流体状态,提出了一种基于半拉格朗日的液体实时仿真方法,并对仿真结果进行了表面构建.该方法首先将Navier-Stokes 方程离散化,并通过求解构造的Poisson 方程得到每一时间步长的数值解,进而精确驱动粒子运动以构建真实液体表面;之后,利用液体表面追踪及Marching Cubes 表面重建,生成了真实的液体表面模型.实验结果表明,该仿真方法不但在运算过程中遵循经典的流体力学方程,从而保证了结果的真实性,并且运算速度快且能取得较好的视觉效果.在计算机游戏、电影制作以及医学等领域的仿真,均有广泛的应用前景.     

High-Order Multiderivative Time Integrators for Hyperbolic Conservation Laws

Multiderivative time integrators have a long history of development for ordinary differential equations, and yet to date, only a small subset of these methods have been explored as a tool for solving partial differential equations (PDEs). This large class of time integrators include all popular (multistage) Runge–Kutta as well as single-step (multiderivative) Taylor methods. (The latter are commonly referred to as Lax–Wendroff methods when applied to PDEs). In this work, we offer explicit multistage multiderivative time integrators for hyperbolic conservation laws. Like Lax–Wendroff methods, multiderivative integrators permit the evaluation of higher derivatives of the unknown in order to decrease the memory footprint and communication overhead. Like traditional Runge–Kutta methods, multiderivative integrators admit the addition of extra stages, which introduce extra degrees of freedom that can be used to increase the order of accuracy or modify the region of absolute stability. We describe a general framework for how these methods can be applied to two separate spatial discretizations: the discontinuous Galerkin (DG) method and the finite difference essentially non-oscillatory (FD-WENO) method. The two proposed implementations are substantially different: for DG we leverage techniques that are closely related to generalized Riemann solvers; for FD-WENO we construct higher spatial derivatives with central differences. Among multiderivative time integrators, we argue that multistage two-derivative methods have the greatest potential for multidimensional applications, because they only require the flux function and its Jacobian, which is readily available. Numerical results indicate that multiderivative methods are indeed competitive with popular strong stability preserving time integrators.     

Superconvergence of Discontinuous Finite Element Solutions for Transient Convection–diffusion Problems

We present a study of the local discontinuous Galerkin method for transient convection–diffusion problems in one dimension. We show that p-degree piecewise polynomial discontinuous finite element solutions of convection-dominated problems are O(Δx ^p+2) superconvergent at Radau points. For diffusion- dominated problems, the solution’s derivative is O(Δx ^p+2) superconvergent at the roots of the derivative of Radau polynomial of degree p+1. Using these results, we construct several asymptotically exact a posteriori finite element error estimates. Computational results reveal that the error estimates are asymptotically exact.This revised version was published online in July 2005 with corrected volume and issue numbers.     

多体系统动力学的常用积分器算法

本文将以单步法中的广义 α族积分器和多步法中的BDF族积分器为主要讨论对象,详细介绍大型多体系统动力学软件中常见类型的积分器的算法细节.每族积分器都给出了不止一套计算公式,而且其对应求解微分代数方程组(DAE)的index可以为1、2或者3.除此以外,本文还着重介绍了微分代数方程组的误差估计、变阶变步长策略等关键技术;并讨论了大型DAE问题求解过程中的初始条件分析、Jacobian矩阵复用等重要环节的算法实现;对于BDF积分器族,文中还详细描述了高阶格式的非绝对稳定性、速度变量的误差估计等瓶颈问题的解决方案.全文以多体系统动力学软件的积分器程序实现为目标,强调在满足给定精度的条件下,如何提高计算效率和保证仿真运行的鲁棒性.另外,本文也简要介绍了在某些应用场合中有很大潜力的显式积分器族.通过分析和比较,文中还将指出各种算法的优缺点以及可能的改进方向,希望能够为研究人员和程序开发者提供一定的参考.由于篇幅限制,本文只列出了几个标准的算例比较,作为文中内容的补充;并给出了几种积分器性能比较的一般性结论.文中几乎所有方法都经由作者程序实现、测试和比较,并且相关算法的实现细节也都已尽量列出,可以很容易地编程实现并应用到实际问题的求解中去.     

A Semi-Lagrangian Method for a Fokker-Planck Equation Describing Fiber Dynamics

A simplified Fokker-Planck model for the lay-down of fibers on a conveyor belt in the production process of nonwovens is investigated. It takes into account the motion of the fiber under the influence of turbulence. The emphasis in this paper is on the development of a numerical procedure to solve the model. We present a semi-Lagrangian scheme that accurately captures the fiber dynamics and conserves the mass. The scheme allows large time steps to be taken in numerical simulations and requires moderate computing times to obtain steady state solutions. Numerical results and examples are presented and compared for several selection of fiber parameters. The obtained results show that the semi-Lagrangian method is able to reproduce accurately the time development of functionals of the process that are important for the quality assessment of industrial fibers.