The adaptation of Crouzeix-Raviart finite element in the context of multi-scale finite element method (MsFEM) is studied and implemented on diffusion and advection-diffusion problems in perforated media. It is known that the approximation of boundary condition on coarse element edges when computing the multiscale basis functions critically influences the eventual accuracy of any MsFEM approaches. The weakly enforced continuity of Crouzeix-Raviart function space across element edges leads to a natural boundary condition for the multiscale basis functions which relaxes the sensitivity of our method to complex patterns of perforations. Another ingredient to our method is the application of bubble functions which is shown to be instrumental in maintaining high accuracy amid dense perforations. Additionally, the application of penalization method makes it possible to avoid complex unstructured domain and allows extensive use of simpler Cartesian meshes.

This paper presents a new and better suited formulation to implement the limiting projection to high-order schemes that make use of high-order local reconstructions for hyperbolic conservation laws. The scheme, so-called MCV-WENO4 (multi-moment Constrained finite Volume with WENO limiter of 4th order) method, is an extension of the MCV method of Ii & Xiao (2009) by adding the 1st order derivative (gradient or slope) at the cell center as an additional constraint for the cell-wise local reconstruction. The gradient is computed from a limiting projection using the WENO (weighted essentially non-oscillatory) reconstruction that is built from the nodal values at 5 solution points within 3 neighboring cells. Different from other existing methods where only the cell-average value is used in the WENO reconstruction, the present method takes account of the solution structure within each mesh cell, and thus minimizes the stencil for reconstruction. The resulting scheme has 4th-order accuracy and is of significant advantage in algorithmic simplicity and computational efficiency. Numerical results of one and two dimensional benchmark tests for scalar and Euler conservation laws are shown to verify the accuracy and oscillation-less property of the scheme.

An unstructured nodal spectral-element method for the Navier-Stokes equations is developed in this paper. The method is based on a triangular and tetrahedral rational approximation and an easy-to-implement nodal basis which fully enjoys the tensorial product property. It allows arbitrary triangular and tetrahedral mesh, affording greater flexibility in handling complex domains while maintaining all essential features of the usual spectral-element method. The details of the implementation and some numerical examples are provided to validate the efficiency and flexibility of the proposed method.

We construct an efficient numerical scheme for the quantum Fokker-Planck- Landau (FPL) equation that works uniformly from kinetic to fluid regimes. Such a scheme inevitably needs an implicit discretization of the nonlinear collision operator, which is difficult to invert. Inspired by work [9] we seek a linear operator to penalize the quantum FPL collision term QqFPL in order to remove the stiffness induced by the small Knudsen number. However, there is no suitable simple quantum operator serving the purpose and for this kind of operators one has to solve the complicated quantum Maxwellians (Bose-Einstein or Fermi-Dirac distribution). In this paper, we propose to penalize QqFPL by the "classical" linear Fokker-Planck operator. It is based on the observation that the classical Maxwellian, with the temperature replaced by the internal energy, has the same first five moments as the quantum Maxwellian. Numerical results for Bose and Fermi gases are presented to illustrate the efficiency of the scheme in both fluid and kinetic regimes.

This paper deals with boundary optimal control problems for the heat and Navier-Stokes equations and addresses the issue of defining controls in function spaces which are naturally associated to the volume variables by trace restriction. For this reason we reformulate the boundary optimal control problem into a distributed problem through a lifting function approach. The stronger regularity requirements which are imposed by standard boundary control approaches can then be avoided. Furthermore, we propose a new numerical strategy that allows to solve the coupled optimality system in a robust way for a large number of unknowns. The optimality system resulting from a finite element discretization is solved by a local multigrid algorithm with domain decomposition Vanka-type smoothers. The purpose of these smoothers is to solve the optimality system implicitly over subdomains with a small number of degrees of freedom, in order to achieve robustness with respect to the regularization parameters in the cost functional. We present the results of some test cases where temperature is the observed quantity and the control quantity corresponds to the boundary values of the fluid temperature in a portion of the boundary. The control region for the observed quantity is a part of the domain where it is interesting to match a desired temperature value.

We study compact finite difference methods for the Schrodinger-Poisson equation in a bounded domain and establish their optimal error estimates under proper regularity assumptions on wave function and external potential V(x). The Crank-Nicolson compact finite difference method and the semi-implicit compact finite difference method are both of order Ҩ(h4 + τ2) in discrete l2,H1 and l∞ norms with mesh size h and time step t. For the errors ofcompact finite difference approximation to the second derivative andPoisson potential are nonlocal, thus besides the standard energy method and mathematical induction method, the key technique in analysisis to estimate the nonlocal approximation errors in discrete l∞ and H1 norm by discrete maximum principle of elliptic equation and properties of some related matrix. Also some useful inequalities are established in this paper. Finally, extensive numerical re-sults are reported to support our error estimates of the numerical methods.

The evolution of precipitates in stressed solids is modeled by coupling a quasi-steady diffusion equation and a linear elasticity equation with dynamic boundary conditions. The governing equations are solved numerically using a boundary integral method (BIM). A critical step in applying BIM is to develop fast algorithms to reduce the arithmetic operation count of matrix-vector multiplications. In this paper, we develop a fast adaptive treecode algorithm for the diffusion and elasticity problems in two dimensions (2D). We present a novel source dividing strategy to parallelize the treecode. Numerical results show that the speedup factor is nearly perfect up to a moderate number of processors. This approach of parallelization can be readily implemented in other treecodes using either uniform or non-uniform point distribution. We demonstrate the effectiveness of the treecode by computing the long-time evolution of a complicated microstructure in elastic media, which would be extremely difficult with a direct summation method due to CPU time constraint. The treecode speeds up computations dramatically while fulfilling the stringent precision requirement dictated by the spectrally accurate BIM.

This paper presents an approach using the method of separation of variables applied to 2D Helmholtz equations in the Cartesian coordinate. The solution is then computed by a series solutions resulted from solving a sequence of 1D problems, in which the 1D solutions are computed using pollution free difference schemes. Moreover, non-polluted numerical integration formulae are constructed to handle the integration due to the forcing term in the inhomogeneous 1D problems. Consequently, the computed solution does not suffer the pollution effect. Another attractive feature of this approach is that a direct method can be effectively applied to solve the tridiagonal matrix resulted from numerical discretization of the 1D Helmholtz equation. The method has been tested to compute 2D Helmholtz solutions simulating electromagnetic scattering from an open large cavity and rectangular waveguide.

We propose here a class of numerical schemes for the approximation of weak solutions to nonlinear hyperbolic systems in nonconservative form—the notion of solution being understood in the sense of Dal Maso, LeFloch, and Murat (DLM). The proposed numerical method falls within LeFloch-Mishra's framework of schemes with well-controlled dissipation (WCD), recently introduced for dealing with small-scale dependent shocks. We design WCD schemes which are consistent with a given nonconservative system at arbitrarily high-order and then analyze their linear stability. We then investigate several nonconservative hyperbolic models arising in complex fluid dynamics, and we numerically demonstrate the convergence of our schemes toward physically meaningful weak solutions.

This paper presents a fast surface voxelization technique for the mapping of tessellated triangular surface meshes to uniform and structured grids that provide a basis for CFD simulations with the lattice Boltzmann method (LBM). The core algorithm is optimized for massively parallel execution on graphics processing units (GPUs) and is based on a unique dissection of the inner body shell. This unique definition necessitates a topology based neighbor search as a preprocessing step, but also enables parallel implementation. More specifically, normal vectors of adjacent triangular tessellations are used to construct half-angles that clearly separate the per-triangle regions. For each triangle, the grid nodes inside the axis-aligned bounding box (AABB) are tested for their distance to the triangle in question and for certain well-defined relative angles. The performance of the presented grid generation procedure is superior to the performance of the GPU-accelerated flow field computations per time step which allows efficient fluid-structure interaction simulations, without noticeable performance loss due to the dynamic grid update.

The coordinate transformation offers a remarkable way to design cloaks that can steer electromagnetic fields so as to prevent waves from penetrating into the cloaked region (denoted by Ω0, where the objects inside are invisible to observers outside). The ideal circular and elliptic cylindrical cloaked regions are blown up from a point and a line segment, respectively so the transformed material parameters and the corresponding coefficients of the resulted equations are highly singular at the cloaking boundary ∂Ω0. The electric field or magnetic field is not continuous across ∂Ω0. The imposition of appropriate cloaking boundary conditions (CBCs) to achieve perfect concealment is a crucial but challenging issue.

Based upon the principle that a well-behaved electromagnetic field in the original space must be well-behaved in the transformed space as well, we obtain CBCs that intrinsically relate to the essential “pole” conditions of a singular transformation. We also find that for the elliptic cylindrical cloak, the CBCs should be imposed differently for the cosine-elliptic and sine-elliptic components of the decomposed fields. With these at our disposal, we can rigorously show that the governing equation in Ω0 can be decoupled from the exterior region , and the total fields in the cloaked region vanish under mild conditions. We emphasize that our proposal of CBCs is different from any existing ones.

Using the exact circular (resp., elliptic) Dirichlet-to-Neumann (DtN) non-reflecting boundary conditions to reduce the unbounded domain to a bounded domain, we introduce an accurate and efficient Fourier-Legendre spectral-element method (FLSEM) (resp., Mathieu-Legendre spectral-element method (MLSEM)) to simulate the circular cylindrical cloak (resp., elliptic cylindrical cloak). We provide ample numerical results to demonstrate that the perfect concealment of waves can be achieved for the ideal circular/elliptic cylindrical cloaks under our proposed CBCs and accurate numerical solvers.

Determining the drag of a flowover a rough surface is a guiding example for the need to take geometric micro-scale effects into account when computing a macroscale quantity. A well-known strategy to avoid a prohibitively expensive numerical resolution of micro-scale structures is to capture the micro-scale effects through some effective boundary conditions posed for a problem on a (virtually) smooth domain. The central objective of this paper is to develop a numerical scheme for accurately capturing the micro-scale effects at essentially the cost of twice solving a problem on a (piecewise) smooth domain at affordable resolution. Here and throughout the paper “smooth” means the absence of any micro-scale roughness. Our derivation is based on a “conceptual recipe” formulated first in a simplified setting of boundary value problems under the assumption of sufficient local regularity to permit asymptotic expansions in terms of the micro-scale parameter.

The proposed multiscale model relies then on an upscaling strategy similar in spirit to previous works by Achdou et al. [1], Jäger and Mikelic [29, 31], Friedmann et al. [24, 25], for incompressible fluids. Extensions to compressible fluids, although with several noteworthy distinctions regarding e.g. the “micro-scale size” relative to boundary layer thickness or the systematic treatment of different boundary conditions, are discussed in Deolmi et al. [16,17]. For proof of concept the general strategy is applied to the compressible Navier-Stokes equations to investigate steady, laminar, subsonic flow over a flat plate with partially embedded isotropic and anisotropic periodic roughness imposing adiabatic and isothermal wall conditions, respectively. The results are compared with high resolution direct simulations on a fully resolved rough domain.

This paper presents a novel high-order space-time method for hyperbolic conservation laws. Two important concepts, the staggered space-time mesh of the space-time conservation element/solution element (CE/SE) method and the local discontinuous basis functions of the space-time discontinuous Galerkin (DG) finite element method, are the two key ingredients of the new scheme. The staggered space-time mesh is constructed using the cell-vertex structure of the underlying spatial mesh. The universal definitions of CEs and SEs are independent of the underlying spatial mesh and thus suitable for arbitrarily unstructured meshes. The solution within each physical time step is updated alternately at the cell level and the vertex level. For this solution updating strategy and the DG ingredient, the new scheme here is termed as the discontinuous Galerkin cell-vertex scheme (DG-CVS). The high order of accuracy is achieved by employing high-order Taylor polynomials as the basis functions inside each SE. The present DG-CVS exhibits many advantageous features such as Riemann-solver-free, high-order accuracy, point-implicitness, compactness, and ease of handling boundary conditions. Several numerical tests including the scalar advection equations and compressible Euler equations will demonstrate the performance of the new method.

We couple different flow models, i.e. a finite element solver for the Navier-Stokes equations and a Lattice Boltzmann automaton, using the framework Peano as a common base. The new coupling strategy between the meso- and macroscopic solver is presented and validated in a 2D channel flow scenario. The results are in good agreement with theory and results obtained in similar works by Latt et al. In addition, the test scenarios show an improved stability of the coupled method compared to pure Lattice Boltzmann simulations.

In the present contribution we evaluate the heat flux prediction capabilities of second-order accurate Residual Distribution (RD) methods in the context of atmospheric (re-)entry problems around blunt bodies. Our departing point is the computation of subsonic air flows (with air modeled either as an inert ideal gas or as chemically reacting and possibly out of thermal equilibrium gas mixture) around probe-like geometries, as those typically employed into high enthalpy wind tunnels. We confirm the agreement between the solutions obtained with the RD method and the solutions computed with other Finite Volume (FV) based codes.

However, a straightforward application of the same numerical technique to hypersonic cases involving strong shocks exhibits severe deficiencies even on a geometry as simple as a 2D cylinder. In an attempt to mitigate this problem, we derive new variants of RD schemes. A comparison of these alternative strategies against established ones allows us to derive a diagnose for the shortcomings observed in the traditional RD schemes.

A conservative modification to the ghost fluid method (GFM) is developed for compressible multiphase flows. The motivation is to eliminate or reduce the conservation error of the GFM without affecting its performance. We track the conservative variables near the material interface and use this information to modify the numerical solution for an interfacing cell when the interface has passed the cell. The modification procedure can be used on the GFM with any base schemes. In this paper we use the fifth order finite difference WENO scheme for the spatial discretization and the third order TVD Runge-Kutta method for the time discretization. The level set method is used to capture the interface. Numerical experiments show that the method is at least mass and momentum conservative and is in general comparable in numerical resolution with the original GFM.

This paper describes an application of the recently developed sparse scheme of the method of fundamental solutions (MFS) for the simulation of three-dimensional modified Helmholtz problems. The solution to the given problems is approximated by a two-step strategy which consists of evaluating the particular solution and the homogeneous solution. The homogeneous solution is approximated by the traditional MFS. The original dense system of the MFS formulation is condensed into a sparse system based on the exponential decay of the fundamental solutions. Hence, the homogeneous solution can be efficiently obtained. The method of particular solutions with polyharmonic spline radial basis functions and the localized method of approximate particular solutions in combination with the Gaussian radial basis function are employed to approximate the particular solution. Three numerical examples including a near singular problem are presented to show the simplicity and effectiveness of this approach.

The differential capacitance of electric double-layer capacitors is studied by developing a generalized model of the self-consistent Gaussian field theory. This model includes many-body effects of particles near the interface such as ionic sizes, the order of water alignment and electrostatic correlations, and thus can present more accurate predictions of the electric double-layer structure and hence the capacitance than traditional continuum theories. Analytical simplification of the model and efficient numerical method are introduced, in particular, the approximation of the self-Green's function which describes the self energy of a mobile ion. We show that, when the applied voltage on interfaces is small the dielectric effect of the electrode materials plays an important role. For large voltage, this effect is screened, but the dielectric saturation due to the alignment of the nearby water is shown to be essential. For 2:1 electrolytes, abnormal enhancement on the capacitance due to the dielectric electrode is observed, which is due to the interplay of the image charge effect and Born solvation energy in the self energy of ions.

Orbital-free density functional theory (OFDFT) is a quantum mechanical method in which the energy of a material depends only on the electron density and ionic positions. We examine some popular algorithms for optimizing the electron density distribution in OFDFT, explaining their suitability, benchmarking their performance, and suggesting some improvements. We start by describing the constrained optimization problem that encompasses electron density optimization. Next, we discuss the line search (including Wolfe conditions) and the nonlinear conjugate gradient and truncated Newton algorithms, as implemented in our open source OFDFT code. We finally focus on preconditioners derived from OFDFT energy functionals. Newly-derived preconditioners are successful for simulation cells of all sizes without regions of low electron-density and for small simulation cells with such regions.

In this paper, we propose a strong stability-preserving predictor-corrector (SSPC) method based on an implicit Runge-Kutta method to solve the acoustic- and elastic-wave equations. We first transform the wave equations into a system of ordinary differential equations (ODEs) and apply the local extrapolation method to discretize the spatial high-order derivatives, resulting in a system of semi-discrete ODEs. Then we use the SSPC method based on an implicit Runge-Kutta method to solve the semi-discrete ODEs and introduce a weighting parameter into the SSPC method. On top of such a structure, we develop a robust numerical algorithm to effectively suppress the numerical dispersion, which is usually caused by the discretization of wave equations when coarse grids are used or geological models have large velocity contrasts between adjacent layers. Meanwhile, we investigate the performance of the SSPC method including numerical errors and convergence rate, numerical dispersion, and stability criteria with different choices of the weighting parameter to solve 1-D and 2-D acoustic- and elastic-wave equations. When the SSPC is applied to seismic simulations, the computational efficiency is also investigated by comparing the SSPC, the fourth-order Lax-Wendroff correction (LWC) method, and the staggered-grid (SG) finite difference method. Comparisons of synthetic waveforms computed by the SSPC and analytic solutions for acoustic and elastic models are given to illustrate the accuracy and the validity of the SSPC method. Furthermore, several numerical experiments are conducted for the geological models including a 2-D homogeneous transversely isotropic (TI) medium, a two-layer elastic model, and the 2-D SEG/EAGE salt model. The results show that the SSPC can be used as a practical tool for large-scale seismic simulation because of its effectiveness in suppressing numerical dispersion even in the situations such as coarse grids, strong interfaces, or high frequencies.