Search

This paper deals with the implementation of high order semi-Lagrangian particle methods to handle high dimensional Vlasov–Poisson systems. It is based on recent developments in the numerical analysis of particle methods and the paper focuses on specific algorithmic features to handle large dimensions. The methods are tested with uniform particle distributions in particular against a recent multi-resolution wavelet based method on a 4D plasma instability case and a 6D gravitational case. Conservation properties, accuracy and computational costs are monitored. The excellent accuracy/cost trade-off shown by the method opens new perspective for accurate simulations of high dimensional kinetic equations by particle methods. [ABSTRACT FROM AUTHOR]

In this paper, a segregated algorithm is proposed to suppress high-frequency oscillations in the velocity field for incompressible flows. In this context, a new velocity formula based on a reconstruction of face fluxes is defined eliminating high-frequency errors. In analogy to the Rhie–Chow interpolation, this approach is equivalent to including a flux-based pressure gradient with a velocity diffusion in the momentum equation. In order to guarantee second-order accuracy of the numerical solver, a set of conditions are defined for the reconstruction operator. To arrive at the final formulation, an outlook over the state of the art regarding velocity reconstruction procedures is presented comparing them through an error analysis. A new operator is then obtained by means of a flux difference minimization satisfying the required spatial accuracy. The accuracy of the new algorithm is analyzed by performing mesh convergence studies for unsteady Navier–Stokes problems with analytical solutions. The stabilization properties of the solver are then tested in a problem where spurious numerical oscillations arise for the velocity field. The results show a remarkable performance of the proposed technique eliminating high-frequency errors without losing accuracy. [ABSTRACT FROM AUTHOR]

This paper presents a stable numerical algorithm for the Brinkman equations by using weak Galerkin (WG) finite element methods. The Brinkman equations can be viewed mathematically as a combination of the Stokes and Darcy equations which model fluid flow in a multi-physics environment, such as flow in complex porous media with a permeability coefficient highly varying in the simulation domain. In such applications, the flow is dominated by Darcy in some regions and by Stokes in others. It is well known that the usual Stokes stable elements do not work well for Darcy flow and vice versa. The challenge of this study is on the design of numerical schemes which are stable for both the Stokes and the Darcy equations. This paper shows that the WG finite element method is capable of meeting this challenge by providing a numerical scheme that is stable and accurate for both Darcy and the Stokes dominated flows. Error estimates of optimal order are established for the corresponding WG finite element solutions. The paper also presents some numerical experiments that demonstrate the robustness, reliability, flexibility and accuracy of the WG method for the Brinkman equations. [ABSTRACT FROM AUTHOR]

Abstract: When one solves PDEs modelling physical phenomena, it is of great importance to take physical constraints into account. More precisely, numerical schemes have to be designed such that discrete solutions satisfy the same constraints as exact solutions. For instance, the underlying physical assumptions for the Euler equations are the positivity of both density and pressure variables. We consider in this paper an unstructured vertex-based tesselation in . Given a MUSCL finite volume scheme and given a reconstruction method (including a limiting process), the point is to determine whether the overall scheme ensures the positivity. The present work is issued from seminal papers from Perthame and Shu (On positivity preserving finite volume schemes for Euler equations, Numer. Math. 73 (1996) 119–130) and Berthon (Robustness of MUSCL schemes for 2D unstructured meshes, J. Comput. Phys. 218 (2) (2006) 495–509). They proved in different frameworks that under assumptions on the corresponding one-dimensional numerical flux, a suitable CFL condition guarantees that density and pressure remain positive. We first analyse Berthon’s method by presenting the ins and outs. We then propose a more general approach adding non geometric degrees of freedom. This approach includes an optimization procedure in order to make the CFL condition explicit and as less restrictive as possible. The reconstruction method is handled independently by means of τ-limiters and of an additional damping parameter. An algorithm is provided in order to specify the adjustments to make in a preexisting code based on a certain numerical flux. Numerical simulations are carried out to prove the accuracy of the method and its ability to deal with low densities and pressures. [Copyright &y& Elsevier]

Abstract: A systematic approach is presented for constructing higher-order immersed boundary and ghost fluid methods for CFD in general, and fluid–structure interaction problems in particular. Such methods are gaining popularity because they simplify a number of computational issues. These range from gridding the fluid domain, to designing and implementing Eulerian-based algorithms for challenging fluid–structure applications characterized by large structural motions and deformations or topological changes. However, because they typically operate on non body-fitted grids, immersed boundary and ghost fluid methods also complicate other issues such as the treatment of wall boundary conditions in general, and fluid–structure transmission conditions in particular. These methods also tend to be at best first-order space-accurate at the immersed interfaces. In some cases, they are also provably inconsistent at these locations. A methodology is presented in this paper for addressing this issue. It is developed for inviscid flows and prescribed structural motions. For the sake of clarity, but without any loss of generality, this methodology is described in one and two dimensions. However, its extensions to flow-induced structural motions and three dimensions are straightforward. The proposed methodology leads to a departure from the current practice of populating ghost fluid values independently from the chosen spatial discretization scheme. Instead, it accounts for the pattern and properties of a preferred higher-order discretization scheme, and attributes ghost values as to preserve the formal order of spatial accuracy of this scheme. It is illustrated in this paper by its application to various finite difference and finite volume methods. Its impact is also demonstrated by one- and two-dimensional numerical experiments that confirm its theoretically proven ability to preserve higher-order spatial accuracy, including in the vicinity of the immersed interfaces. [Copyright &y& Elsevier]

Abstract: In this paper, we tackle the modeling and numerical simulation of sprays and aerosols, that is dilute gas–droplet flows for which polydispersity description is of paramount importance. Starting from a kinetic description for point particles experiencing transport either at the carrier phase velocity for aerosols or at their own velocity for sprays as well as evaporation, we focus on an Eulerian high order moment method in size and consider a system of partial differential equations (PDEs) on a vector of successive integer size moments of order 0 to N, N >2, over a compact size interval. There exists a stumbling block for the usual approaches using high order moment methods resolved with high order finite volume methods: the transport algorithm does not preserve the moment space. Indeed, reconstruction of moments by polynomials inside computational cells coupled to the evolution algorithm can create N-dimensional vectors which fail to be moment vectors: it is impossible to find a size distribution for which there are the moments. We thus propose a new approach as well as an algorithm which is second order in space and time with very limited numerical diffusion and allows to accurately describe the advection process and naturally preserves the moment space. The algorithm also leads to a natural coupling with a recently designed algorithm for evaporation which also preserves the moment space; thus polydispersity is accounted for in the evaporation and advection process, very accurately and at a very reasonable computational cost. These modeling and algorithmic tools are referred to as the Eulerian Multi Size Moment (EMSM) model. We show that such an approach is very competitive compared to multi-fluid approaches, where the size phase space is discretized into several sections and low order moment methods are used in each section, as well as with other existing high order moment methods. An accuracy study assesses the order of the method as well as the low level of numerical diffusion on structured meshes. Whereas the extension to unstructured meshes is provided, we focus in this paper on cartesian meshes and two 2D test-cases are presented: Taylor–Green vortices and turbulent free jets, where the accuracy and efficiency of the approach are assessed. [Copyright &y& Elsevier]

Abstract: Modeling of laser-plasma wakefield accelerators in an optimal frame of reference has been shown to produce orders of magnitude speed-up of calculations from first principles. Obtaining these speedups required mitigation of a high-frequency instability that otherwise limits effectiveness. In this paper, methods are presented which mitigated the observed instability, including an electromagnetic solver with tunable coefficients, its extension to accommodate Perfectly Matched Layers and Friedman’s damping algorithms, as well as an efficient large bandwidth digital filter. It is observed that choosing the frame of the wake as the frame of reference allows for higher levels of filtering or damping than is possible in other frames for the same accuracy. Detailed testing also revealed the existence of a singular time step at which the instability level is minimized, independently of numerical dispersion. A combination of the techniques presented in this paper prove to be very efficient at controlling the instability, allowing for efficient direct modeling of 10GeV class laser plasma accelerator stages. The methods developed in this paper may have broader application, to other Lorentz-boosted simulations and Particle-In-Cell simulations in general. [Copyright &y& Elsevier]

Abstract: This paper describes the algorithms and functionality of a new module developed to support overset grid assembly associated with performing time-dependent and adaptive moving body calculations of external aerodynamic flows using a multi-solver paradigm (i.e. different CFD solvers in different parts of the computational domain). We use the term “domain connectivity” in this paper to denote all the procedures that are involved in an overset grid assembly, and the module developed is referred henceforth as the domain-connectivity module. The domain-connectivity module coordinates the data transfer between different solvers applied in different parts of the computational domain – body fitted structured or unstructured to capture viscous near-wall effects, and Cartesian adaptive mesh refinement to capture effects away from the wall. The execution of the CFD solvers and the domain-connectivity module are orchestrated by a Python-based computational infrastructure. The domain-connectivity module is fully parallel and performs all its operations (identification of grid overlaps and determination of data interpolation strategy) on the partitioned grid data. In addition, the domain connectivity procedures are completely automated such that no user intervention or manual input is necessary. The capabilities and performance of the package are presented for several test problems, including flow over a NACA 0015 wing and an AGARD A2 slotted airfoil, hover simulation of a scaled V-22 rotor, and dynamic simulation of a UH-60A rotor in forward flight. A modification to the algorithm for improved domain connectivity solutions in problems with tight tolerances as well as heterogeneous grid clustering is also presented. [Copyright &y& Elsevier]

Abstract: In this paper we analyze the equilibrium limit of the constitutive model for two-phase granular mixtures introduced in Papalexandris (2004) , and develop an algorithm for its numerical approximation. At, equilibrium, the constitutive model reduces to a strongly coupled, overdetermined system of quasilinear elliptic partial differential equations with respect to the pressure and the volume fraction of the solid granular phase. First we carry a perturbation analysis based on standard hydrostatic-type scaling arguments which reduces the complexity of the coupling of the equations. The perturbed system is then supplemented by an appropriate compatibility condition which arises from the properties of the gradient operator. Further, based on the Helmholtz decomposition and Ladyzhenskaya’s decomposition theorem, we develop a projection-type, Successive-Over-Relaxation numerical method. This method is general enough and can be applied to a variety of continuum models of complex mixtures and mixtures with micro-structure. We also prove that this method is both stable and consistent hence, under standard assumptions, convergent. The paper concludes with the presentation of representative numerical results. [Copyright &y& Elsevier]

Abstract: The present paper describes an iterative procedure for BEM–BEM coupling. The paper presents suitable interface conditions and algorithms for iteratively coupling sub-domains modeled by three different boundary element time-domain formulations, namely: acoustic and elastodynamic BEM formulations based on time-dependent Green’s functions and non-linear time-domain approach which employs elastostatic Green’s functions and therefore requires domain discretization. Two examples are analyzed and at the end of the paper conclusions of the study are presented. [Copyright &y& Elsevier]

The formation of shock waves in solutions of hyperbolic conservation laws calls for locally adaptive numerical solution algorithms and requires a practical tool for identifying where adaption is needed. In this paper, a new smoothness indicator (SI) is used to identify “rough” solution regions and is implemented in locally adaptive algorithms. The SI is based on the weak local truncation error of the approximate solution. It was recently reported in S. Karni and A. Kurganov, Local error analysis for approximate solutions of hyperbolic conservation laws, where error analysis and convergence properties were established. The present paper is concerned with its implementation in scheme adaption and mesh adaption algorithms. The SI provides a general framework for adaption and is not restricted to a particular discretization scheme. The implementation in this paper uses the central-upwind scheme of A. Kurganov, S. Noelle, and G. Petrova, SIAM J. Sci. Comput.23, 707 (2001). The extension of the SI to two space dimensions is given. Numerical results in one and two space dimensions demonstrate the robustness of the proposed SI and its potential in reducing computational costs and improving the resolution of the solution. [Copyright &y& Elsevier]

In this paper, some numerical algorithms (spectral collocation method, block spectral collocation method, boundary value method, block boundary value method, implicit Runge–Kutta method, diagonally implicit Runge–Kutta method and total variation diminishing Runge–Kutta method) are used to solve the highly oscillatory second-order initial value problems. We first derive these methods for the first-order initial value problems, and then extend these methods to the highly oscillatory nonlinear systems by matrix analysis methods. These new methods preserve the accuracy of the original methods and the main advantages of these new methods are low storage requirements and high efficiency. Extensive numerical results are presented to demonstrate the convergence properties of these methods. [ABSTRACT FROM AUTHOR]

Fourier extensions have been shown to be an effective means for the approximation of smooth, nonperiodic functions on bounded intervals given their values on an equispaced, or in general, scattered grid. Related to this method are two parameters. These are the extension parameter T (the ratio of the size of the extended domain to the physical domain) and the oversampling ratio η (the number of sampling nodes per Fourier mode). The purpose of this paper is to investigate how the choice of these parameters affects the accuracy and stability of the approximation. Our main contribution is to document the following interesting phenomenon: namely, if the desired condition number of the algorithm is fixed in advance, then the particular choice of such parameters makes little difference to the algorithm's accuracy. As a result, one is free to choose T without concern that it is suboptimal. In particular, one may use the value T=2 - which corresponds to the case where the extended domain is precisely twice the size of the physical domain - for which there is known to be a fast algorithm for computing the approximation. In addition, we determine the resolution power (points-per-wavelength) of the approximation to be equal to Tη, and address the trade-off between resolution power and stability. [ABSTRACT FROM AUTHOR]

Abstract: This is a research expository paper regarding superparameterization, a class of multi-scale numerical methods designed to cope with the intermittent multi-scale effects of inhomogeneous geophysical turbulence where energy often inverse-cascades from the unresolved scales to the large scales through the effects of waves, jets, vortices, and latent heat release from moist processes. Original as well as sparse space–time superparameterization algorithms are discussed for the important case of moist atmospheric convection including the role of multi-scale asymptotic methods in providing self-consistent constraints on superparameterization algorithms and related deterministic and stochastic multi-cloud parameterizations. Test models for the statistical numerical analysis of superparameterization algorithms are discussed both to elucidate the performance of the basic algorithms and to test their potential role in efficient multi-scale data assimilation. The very recent development of grid-free seamless stochastic superparameterization methods for geophysical turbulence appropriate for “eddy-permitting” mesoscale ocean turbulence is presented here including a general formulation and illustrative applications to two-layer quasigeostrophic turbulence, and another difficult test case involving one-dimensional models of dispersive wave turbulence. This last test case has randomly generated solitons as coherent structures which collapse and radiate wave energy back to the larger scales, resulting in strong direct and inverse turbulent energy cascades. [Copyright &y& Elsevier]

Abstract: Stable partitioned algorithms for fluid–structure interaction (FSI) problems are developed and analyzed in this two-part paper. Part I describes an algorithm for incompressible flow coupled with compressible elastic solids, while Part II discusses an algorithm for incompressible flow coupled with structural shells. The numerical approach described here for structural shells uses Robin (mixed) interface conditions for the pressure and velocity in the fluid which are derived directly from the governing equations. The resulting added-mass partitioned (AMP) algorithm is stable even for very light structures, requires no sub-iterations per time step, and is second-order accurate. The stability and accuracy of the AMP algorithm is evaluated for linearized FSI model problems. A normal mode analysis is performed to show that the new AMP algorithm is stable, even for the case of very light structures when added-mass effects are large. Exact traveling wave solutions are derived for the FSI model problems, and these solutions are used to verify the stability and accuracy of the corresponding numerical results obtained from the AMP algorithm for the cases of light, medium and heavy structures. A summary comparison of the AMP algorithm developed here and the one in Part I is provided. [Copyright &y& Elsevier]

Abstract: The paper describes an algorithm for PDE-constrained optimization that controls numerical errors using error estimates and grid adaptation during the optimization process. A key aspect of the algorithm is the use of adjoint variables to estimate errors in the first-order optimality conditions. Multilevel optimization is used to drive the optimality conditions and their estimated errors below a specified tolerance. The error estimate requires two additional adjoint solutions, but only at the beginning and end of each optimization cycle. Moreover, the adjoint systems can be formed and solved with limited additional infrastructure beyond that found in typical PDE-constrained optimization algorithms. The approach is general and can accommodate both reduced-space and full-space formulations of the optimization problem. The algorithm is illustrated using the inverse design of a nozzle constrained by the quasi-one-dimensional Euler equations. [Copyright &y& Elsevier]

Abstract: The pseudo-spectral analytical time-domain (PSATD) particle-in-cell (PIC) algorithm solves the vacuum Maxwellʼs equations exactly, has no Courant time-step limit (as conventionally defined), and offers substantial flexibility in plasma and particle beam simulations. It is, however, not free of the usual numerical instabilities, including the numerical Cherenkov instability, when applied to relativistic beam simulations. This paper derives and solves the numerical dispersion relation for the PSATD algorithm and compares the results with corresponding behavior of the more conventional pseudo-spectral time-domain (PSTD) and finite difference time-domain (FDTD) algorithms. In general, PSATD offers superior stability properties over a reasonable range of time steps. More importantly, one version of the PSATD algorithm, when combined with digital filtering, is almost completely free of the numerical Cherenkov instability for time steps (scaled to the speed of light) comparable to or smaller than the axial cell size. [Copyright &y& Elsevier]

Abstract: In this paper, a novel algorithm for incompressible flows of two-phase granular mixtures is proposed. The algorithm constitutes a generalization of projection-type methods to multi-phase flows and employs a predictor–corrector time-marching scheme. For the computation of the phasial pressures, a Poisson equation and a second-order elliptic PDE with variable coefficients are solved at both the prediction and the correction stage. Further, the algorithm is combined with a regularization scheme for the numerical treatment of material interfaces which are a common characteristic of the flows of interest. The efficiency and robustness of the proposed algorithm are illustrated via numerical simulations of unsteady two-dimensional shear flows for a reference mixture of beach sand and water. [Copyright &y& Elsevier]

Abstract: Explicit time stepping renders many high-resolution computational schemes to become less efficient when dealing with non-uniform grids typical of many aeroacoustic applications. Asynchronous time stepping, i.e., updating the solution in different cell sizes according to their local rates, is known to be a promising way to improve the efficiency of explicit time-stepping methods without compromise in accuracy. In the present paper, a new asynchronous time-stepping algorithm is developed for the Compact Accurately Boundary-Adjusting high-REsolution Technique (CABARET) Euler method. This allows to significantly speedup the original single-step CABARET method with non-uniform grids and improves its accuracy at the same time. Numerical examples are provided and issues associated with the method performance on various grid resolutions are discussed. [Copyright &y& Elsevier]

Abstract: In this paper we present a new algorithm based on Cartesian meshes for the numerical approximation of kinetic models set in an arbitrary geometry. Due to the high dimensional property of kinetic models, numerical algorithms based on unstructured meshes are not really appropriate since most of numerical methods (semi-Lagrangian, spectral methods) are particularly efficient on structured grids. Here we propose to adapt the inverse Lax–Wendroff procedure, which has been recently introduced for conservation laws [21], for kinetic equations. Numerical simulations in and based on this approach are proposed for Boltzmann type operators (BGK, ES-BGK models). [Copyright &y& Elsevier]

• A linearly operator splitting algorithm is proposed for the nonlocal VCH equation. • The energy stabilities for both subproblems are proved. • The stability and convergence of the operator splitting algorithm are studied. • A semi-implicit SDC method is further used to improve time accuracy. Recently, the viscous Cahn-Hilliard (VCH) equation has been proposed as a phenomenological continuum model for phase separation in glass and polymer systems where intermolecular friction forces become important. Compared with the classical local VCH model, the nonlocal VCH model equipped with nonlocal diffusion operator can describe more practical phenomena for modeling phase transitions of microstructures in materials. This paper presents a high order fast explicit method based on operator splitting and spectral deferred correction (SDC) for solving the nonlocal VCH equation. We start with a second-order operator splitting spectral scheme, which is based on the Fourier spectral method and the strong stability preserving Runge-Kutta (SSP-RK) method. The scheme takes advantage of applying the fast Fourier transform (FFT) and avoiding nonlinear iteration. The stability and convergence analysis of the obtained numerical scheme are analyzed. To improve the temporal accuracy, the semi-implicit SDC method is then introduced. Various numerical simulations are performed to validate the accuracy and efficiency of the proposed method. [ABSTRACT FROM AUTHOR]

Abstract: In this paper, we derive a series of local structure-preserving algorithms for the “good” Boussinesq equation, including multisymplectic geometric structure-preserving algorithms, local energy-preserving algorithms and local momentum-preserving algorithms. The outstanding advantage of the proposed algorithms is that they conserve these local structures in any time-space region exactly. For example, the proposed local energy-preserving algorithms preserve the local energy conservation law in any local domain. Therefore, the local structure-preserving algorithms overcome the shortage of global structure-preserving algorithms on the boundary conditions. Especially, with suitable boundary conditions such as periodic or homogeneous boundary conditions, the local structure-preserving algorithms will be global structure-preserving algorithms. Numerical results verify the theoretical analysis. [Copyright &y& Elsevier]

Abstract: This paper is concerned with spectral Galerkin algorithms for solving high even-order two point boundary value problems in one dimension subject to homogeneous and nonhomogeneous boundary conditions. The proposed algorithms are extended to solve two-dimensional high even-order differential equations. The key to the efficiency of these algorithms is to construct compact combinations of Chebyshev polynomials of the third and fourth kinds as basis functions. The algorithms lead to linear systems with specially structured matrices that can be efficiently inverted. Numerical examples are included to demonstrate the validity and applicability of the proposed algorithms, and some comparisons with some other methods are made. [Copyright &y& Elsevier]

Abstract: The paper presents a detailed numerical study of an iterative solution to 3-D sound-hard acoustic scattering problems at high frequency considering the Combined Field Integral Equation (CFIE). We propose a combination of an OSRC preconditioning technique and a Fast Multipole Method which leads to a fast and efficient algorithm independently of both a frequency increase and a mesh refinement. The OSRC-preconditioned CFIE exhibits very interesting spectral properties even for trapping domains. Moreover, this analytic preconditioner shows highly-desirable advantages: sparse structure, ease of implementation and low additional computational cost. We first investigate the numerical behavior of the eigenvalues of the related integral operators, CFIE and OSRC-preconditioned CFIE, in order to illustrate the influence of the proposed preconditioner. We then apply the resolution algorithm to various and significant test-cases using a GMRES solver. The OSRC-preconditioning technique is combined to a Fast Multipole Method in order to deal with high-frequency 3-D cases. This variety of tests validates the effectiveness of the method and fully justifies the interest of such a combination. [Copyright &y& Elsevier]

Abstract: Jacobian-free Newton–Krylov (JFNK) method is an effective algorithm for solving large scale nonlinear equations. One of the most important advantages of JFNK method is that there is no necessity to form and store the Jacobian matrix of the nonlinear system when JFNK method is employed. However, an approximation of the Jacobian is needed for the purpose of preconditioning. In this paper, JFNK method is employed to solve a class of non-equilibrium radiation diffusion coupled to material thermal conduction equations, and two preconditioners are designed by linearizing the equations in two methods. Numerical results show that the two preconditioning methods can improve the convergence behavior and efficiency of JFNK method. [Copyright &y& Elsevier]

Abstract: Recently, an implicit, nonlinearly consistent, energy- and charge-conserving one-dimensional (1D) particle-in-cell method has been proposed for multi-scale, full-f kinetic simulations [G. Chen et al., J. Comput. Phys. 230 (18) (2011)]. The method employs a Jacobian-free Newton–Krylov (JFNK) solver, capable of using very large timesteps without loss of numerical stability or accuracy. A fundamental feature of the method is the segregation of particle-orbit computations from the field solver, while remaining fully self-consistent. This paper describes a very efficient, mixed-precision hybrid CPU–GPU implementation of the 1D implicit PIC algorithm exploiting this feature. The JFNK solver is kept on the CPU in double precision (DP), while the implicit, charge-conserving, and adaptive particle mover is implemented on a GPU (graphics processing unit) using CUDA in single-precision (SP). Performance-oriented optimizations are introduced with the aid of the roofline model. The implicit particle mover algorithm is shown to achieve up to 400GOp/s on a Nvidia GeForce GTX580. This corresponds to 25% absolute GPU efficiency against the peak theoretical performance, and is about 100 times faster than an equivalent single-core CPU (Intel Xeon X5460) compiler-optimized execution. For the test case chosen, the mixed-precision hybrid CPU–GPU solver is shown to over-perform the DP CPU-only serial version by a factor of ∼100, without apparent loss of robustness or accuracy in a challenging long-timescale ion acoustic wave simulation. [Copyright &y& Elsevier]

Abstract: This paper deals with a new solid–fluid coupling algorithm between a rigid body and an unsteady compressible fluid flow, using an Embedded Boundary method. The coupling with a rigid body is a first step towards the coupling with a Discrete Element method. The flow is computed using a finite volume approach on a Cartesian grid. The expression of numerical fluxes does not affect the general coupling algorithm and we use a one-step high-order scheme proposed by Daru and Tenaud [V. Daru, C. Tenaud, J. Comput. Phys. (2004)]. The Embedded Boundary method is used to integrate the presence of a solid boundary in the fluid. The coupling algorithm is totally explicit and ensures exact mass conservation and a balance of momentum and energy between the fluid and the solid. It is shown that the scheme preserves uniform movement of both fluid and solid and introduces no numerical boundary roughness. The efficiency of the method is demonstrated on challenging one- and two-dimensional benchmarks. [Copyright &y& Elsevier]

Abstract: In this paper, a new flux-based one-step hybrid remapping method for multi-material arbitrary Lagrangian–Eulerian (ALE) approach is introduced. In the vicinity of material interfaces, the swept region is intersected with pure material polygons in the Lagrangian mesh to construct the material fluxes. Far from interfaces, the fluxes are constructed in a standard swept-region manner without intersections. This method is conservative, second-order accurate and linearity-preserving (in case of straight material interfaces), and faster than method based on intersections, as shown on selected numerical examples. [Copyright &y& Elsevier]

Abstract: Hydrocarbon systems have been intensively studied via numerical methods, including electronic structure computations, molecular dynamics and Monte Carlo simulations. Typically, these methods require an initial structural model (atomic positions and types, topology, etc.) that may be produced using scripts and/or modeling tools. For many systems, however, these building methods may be ineffective, as the user may have to specify the positions of numerous atoms while maintaining structural plausibility. In this paper, we present an interactive physically-based modeling tool to construct structural models of hydrocarbon systems. As the user edits the geometry of the system, atomic positions are also influenced by the Brenner potential, a well-known bond-order reactive potential. In order to be able to interactively edit systems containing numerous atoms, we introduce a new adaptive simulation algorithm, as well as a novel algorithm to incrementally update the forces and the total potential energy based on the list of updated relative atomic positions. The computational cost of the adaptive simulation algorithm depends on user-defined error thresholds, and our potential update algorithm depends linearly with the number of updated bonds. This allows us to enable efficient physically-based editing, since the computational cost is decoupled from the number of atoms in the system. We show that our approach may be used to effectively build realistic models of hydrocarbon structures that would be difficult or impossible to produce using other tools. [Copyright &y& Elsevier]

Abstract: High-order compact finite difference scheme with operator splitting technique for solving two-dimensional time fractional diffusion equation is considered in this paper. A Grünwald–Letnikov approximation is used for the Riemann–Liouville time derivative, and the second order spatial derivatives are approximated by the compact finite differences to obtain a fully discrete implicit scheme. Alternating direction implicit (ADI) method is used to split the original problem into two separate one-dimensional problems. The local truncation error is analyzed and the stability is discussed by the Fourier method. The proposed scheme is suitable when the order of the time fractional derivative γ lies in the interval . A correction term is added to maintain high accuracy when . Numerical results are provided to verify the accuracy and efficiency of the proposed algorithm. [Copyright &y& Elsevier]

Abstract: The reconstruction of the equilibrium of a plasma in a Tokamak is a free boundary problem described by the Grad–Shafranov equation in axisymmetric configuration. The right-hand side of this equation is a nonlinear source, which represents the toroidal component of the plasma current density. This paper deals with the identification of this nonlinearity source from experimental measurements in real time. The proposed method is based on a fixed point algorithm, a finite element resolution, a reduced basis method and a least-square optimization formulation. This is implemented in a software called Equinox with which several numerical experiments are conducted to explore the identification problem. It is shown that the identification of the profile of the averaged current density and of the safety factor as a function of the poloidal flux is very robust. [Copyright &y& Elsevier]

Abstract: In this paper, we present an efficient numerical method for evaluating rare failure probability. The method is based on a recently developed surrogate-based method from Li and Xiu [J. Li, D. Xiu, Evaluation of failure probability via surrogate models, J. Comput. Phys. 229 (2010) 8966–8980] for failure probability computation. The method by Li and Xiu is of hybrid nature, in the sense that samples of both the surrogate model and the true physical model are used, and its efficiency gain relies on using only very few samples of the true model. Here we extend the capability of the method to rare probability computation by using the idea of importance sampling (IS). In particular, we employ cross-entropy (CE) method, which is an effective method to determine the biasing distribution in IS. We demonstrate that, by combining with the CE method, a surrogate-based IS algorithm can be constructed and is highly efficient for rare failure probability computation—it incurs much reduced simulation efforts compared to the traditional CE-IS method. In many cases, the new method is capable of capturing failure probability as small as 10−12 ∼10−6 with only several hundreds samples. [Copyright &y& Elsevier]

Abstract: Common efficient schemes for the incompressible Navier–Stokes equations, such as projection or fractional step methods, have limited temporal accuracy as a result of matrix splitting errors, or introduce errors near the domain boundaries (which destroy uniform convergence to the solution). In this paper we recast the incompressible (constant density) Navier–Stokes equations (with the velocity prescribed at the boundary) as an equivalent system, for the primary variables velocity and pressure. equation for the pressure. The key difference from the usual approaches occurs at the boundaries, where we use boundary conditions that unequivocally allow the pressure to be recovered from knowledge of the velocity at any fixed time. This avoids the common difficulty of an, apparently, over-determined Poisson problem. Since in this alternative formulation the pressure can be accurately and efficiently recovered from the velocity, the recast equations are ideal for numerical marching methods. The new system can be discretized using a variety of methods, including semi-implicit treatments of viscosity, and in principle to any desired order of accuracy. In this work we illustrate the approach with a 2-D second order finite difference scheme on a Cartesian grid, and devise an algorithm to solve the equations on domains with curved (non-conforming) boundaries, including a case with a non-trivial topology (a circular obstruction inside the domain). This algorithm achieves second order accuracy in the L ∞ norm, for both the velocity and the pressure. The scheme has a natural extension to 3-D. [Copyright &y& Elsevier]

Abstract: The paper introduces a high-order, adaptive finite-element method for simulation of sonic measurements acquired with borehole-eccentered logging instruments. The resulting frequency-domain based algorithm combines a Fourier series expansion in one spatial dimension with a two-dimensional high-order adaptive finite-element method (FEM), and incorporates a perfectly matched layer (PML) for truncation of the computational domain. The simulation method was verified for various model problems, including a comparison to a semi-analytical solution developed specifically for this purpose. Numerical results indicate that for a wireline sonic tool operating in a fast formation, the main propagation modes are insensitive to the distance from the center of the tool to the center of the borehole (eccentricity distance). However, new flexural modes arise with an increase in eccentricity distance. In soft formations, we identify a new dipole tool mode which arises as a result of tool eccentricity. [Copyright &y& Elsevier]

Abstract: Vesicles are locally-inextensible fluid membranes that can sustain bending. In this paper, we extend the study of Veerapaneni et al. [S.K. Veerapaneni, D. Gueyffier, G. Biros, D. Zorin, A numerical method for simulating the dynamics of 3D axisymmetric vesicles suspended in viscous flows, Journal of Computational Physics 228 (19) (2009) 7233–7249] to general non-axisymmetric vesicle flows in three dimensions. Although the main components of the algorithm are similar in spirit to the axisymmetric case (spectral approximation in space, semi-implicit time-stepping scheme), important new elements need to be introduced for a full 3D method. In particular, spatial quantities are discretized using spherical harmonics, and quadrature rules for singular surface integrals need to be adapted to this case; an algorithm for surface reparameterization is needed to ensure stability of the time-stepping scheme, and spectral filtering is introduced to maintain reasonable accuracy while minimizing computational costs. To characterize the stability of the scheme and to construct preconditioners for the iterative linear system solvers used in the semi-implicit time-stepping scheme, we perform a spectral analysis of the evolution operator on the unit sphere. By introducing these algorithmic components, we obtain a time-stepping scheme that circumvents the stability constraint on the time-step and achieves spectral accuracy in space. We present results to analyze the cost and convergence rates of the overall scheme. To illustrate the applicability of the new method, we consider a few vesicle-flow interaction problems: a single vesicle in relaxation, sedimentation, shear flows, and many-vesicle flows. [Copyright &y& Elsevier]

Abstract: The Huygens subgridding (HSG) is a subgridding technique developed for the numerical solution of the Maxwell equations. It relies on the theoretical equivalence of any physical volume with two or more fictitious volumes connected by equivalent currents. The application of this concept to the finite-difference time-domain (FDTD) method has been previously published in the one dimensional and two dimensional cases. In this paper the HSG is extended to the general three dimensional case, the exchange of the electromagnetic energy between the two FDTD grids is investigated theoretically, and some modifications to the HSG algorithm are presented with the objective of simplifying its implementation. [Copyright &y& Elsevier]

Abstract: We develop and study the high-order conservative and monotone optimization-based remap (OBR) of a scalar conserved quantity (mass) between two close meshes with the same connectivity. The key idea is to phrase remap as a global inequality-constrained optimization problem for mass fluxes between neighboring cells. The objective is to minimize the discrepancy between these fluxes and the given high-order target mass fluxes, subject to constraints that enforce physically motivated bounds on the associated primitive variable (density). In so doing, we separate accuracy considerations, handled by the objective functional, from the enforcement of physical bounds, handled by the constraints. The resulting OBR formulation is applicable to general, unstructured, heterogeneous grids. Under some weak requirements on grid proximity, but not on the cell types, we prove that the OBR algorithm is linearity preserving in one, two and three dimensions. The paper also examines connections between the OBR and the recently proposed flux-corrected remap (FCR), Liska et al. . We show that the FCR solution coincides with the solution of a modified version of OBR (M-OBR), which has the same objective but a simpler set of box constraints derived by using a “worst-case” scenario. Because M-OBR (FCR) has a smaller feasible set, preservation of linearity may be lost and accuracy may suffer for some grid configurations. Our numerical studies confirm this, and show that OBR delivers significant increases in robustness and accuracy. Preliminary efficiency studies of OBR reveal that it is only a factor of 2.1 slower than FCR, but admits 1.5 times larger time steps. [Copyright &y& Elsevier]

Abstract: In this paper, we present a method to compute compressible flows in 2D. It uses two steps: a Lagrangian step and a metric-based triangular mesh adaptation step. Computational mesh is locally adapted according to some metric field that depends on physical or geometrical data. This mesh adaptation step embeds a conservative remapping procedure to satisfy consistency with Euler equations. The whole method is no more Lagrangian. After describing mesh adaptation patterns, we recall the metric formalism. Then, we detail an appropriate remapping procedure which is first-order and relies on exact intersections. We give some hints about the parallel implementation. Finally, we present various numerical experiments which demonstrate the good properties of the algorithm. [ABSTRACT FROM AUTHOR]

Abstract: The deterministic transport code HZETRN was developed for research scientists and design engineers studying the effects of space radiation on astronauts and instrumentation protected by various shielding materials and structures. In this work, several aspects of code verification are examined. First, a detailed derivation of the light particle (A ⩽4) and heavy ion (A >4) numerical marching algorithms used in HZETRN is given. References are given for components of the derivation that already exist in the literature, and discussions are given for details that may have been absent in the past. The present paper provides a complete description of the numerical methods currently used in the code and is identified as a key component of the verification process. Next, a new numerical method for light particle transport is presented, and improvements to the heavy ion transport algorithm are discussed. A summary of round-off error is also given, and the impact of this error on previously predicted exposure quantities is shown. Finally, a coupled convergence study is conducted by refining the discretization parameters (step-size and energy grid-size). From this study, it is shown that past efforts in quantifying the numerical error in HZETRN were hindered by single precision calculations and computational resources. It is determined that almost all of the discretization error in HZETRN is caused by the use of discretization parameters that violate a numerical convergence criterion related to charged target fragments below 50AMeV. Total discretization errors are given for the old and new algorithms to 100g/cm2 in aluminum and water, and the improved accuracy of the new numerical methods is demonstrated. Run time comparisons between the old and new algorithms are given for one, two, and three layer slabs of 100g/cm2 of aluminum, polyethylene, and water. The new algorithms are found to be almost 100 times faster for solar particle event simulations and almost 10 times faster for galactic cosmic ray simulations. [Copyright &y& Elsevier]

Abstract: Diffusion, lossy wave, and Klein–Gordon equations find numerous applications in practical problems across a range of diverse disciplines. The temporal dependence of all three Green’s functions are characterized by an infinite tail. This implies that the cost complexity of the spatio-temporal convolutions, associated with evaluating the potentials, scales as , where N s and N t are the number of spatial and temporal degrees of freedom, respectively. In this paper, we discuss two new methods to rapidly evaluate these spatio-temporal convolutions by exploiting their block-Toeplitz nature within the framework of accelerated Cartesian expansions (ACE). The first scheme identifies a convolution relation in time amongst ACE harmonics and the fast Fourier transform (FFT) is used for efficient evaluation of these convolutions. The second method exploits the rank deficiency of the ACE translation operators with respect to time and develops a recursive numerical compression scheme for the efficient representation and evaluation of temporal convolutions. It is shown that the cost of both methods scales as . Several numerical results are presented for the diffusion equation to validate the accuracy and efficacy of the fast algorithms developed here. [ABSTRACT FROM AUTHOR]

Abstract: We propose the backward phase flow method to implement the Fourier–Bros–Iagolnitzer (FBI)-transform-based Eulerian Gaussian beam method for solving the Schrödinger equation in the semi-classical regime. The idea of Eulerian Gaussian beams has been first proposed in . In this paper we aim at two crucial computational issues of the Eulerian Gaussian beam method: how to carry out long-time beam propagation and how to compute beam ingredients rapidly in phase space. By virtue of the FBI transform, we address the first issue by introducing the reinitialization strategy into the Eulerian Gaussian beam framework. Essentially we reinitialize beam propagation by applying the FBI transform to wavefields at intermediate time steps when the beams become too wide. To address the second issue, inspired by the original phase flow method, we propose the backward phase flow method which allows us to compute beam ingredients rapidly. Numerical examples demonstrate the efficiency and accuracy of the proposed algorithms. [ABSTRACT FROM AUTHOR]

Abstract: This paper introduces a wavepacket-transform-based Gaussian beam method for solving the Schrödinger equation. We focus on addressing two computational issues of the Gaussian beam method: how to generate a Gaussian beam representation for general initial conditions and how to perform long time propagation for any finite period of time. To address the first question, we introduce fast Gaussian wavepacket transforms and develop on top of them an efficient initialization algorithm for general initial conditions. Based on this new initialization algorithm, we address the second question by reinitializing the beam representation when the beams become too wide. Numerical examples in one, two, and three dimensions demonstrate the efficiency and accuracy of the proposed algorithms. The methodology can be readily generalized to deal with other semi-classical quantum mechanical problems. [ABSTRACT FROM AUTHOR]

Abstract: The time-marching approach has clear physical meaning and strict mathematical nature and has been applied in computation of compressible flows widely and extended to many uniform algorithms for All-Speed flows. Remedy for its weakness in the problem of checkerboard decoupling of pressure field for incompressible flows is envisaged with the time-marching momentum interpolation method (MIM) taken into account in this paper. Existing preconditioning methods for suppressing decoupling and time-marching MIM are analyzed for this purpose, and algorithms of time-marching MIM are proposed for steady and unsteady flows and for All-Speed flows. Asymptotic analysis shows that the supposed time-marching MIM has at least a third-order accuracy, better than the existing time-marching coupling methods, which only have an accuracy of the same order as the adopted scheme has. Effects of the time step sizes on the ability of the time-marching MIM to suppress the checkerboard pressure decoupling are particularly discussed in terms of the dual-time stepping approach, and it is revealed how the decreased sizes of either the pseudo- or physical-time step increases the possibility of decoupling and how Choi’s modification, in which the history of the interface velocity is decided by itself instead of the arithmetic average of the velocities on its adjacent nodes, removes the unphysical pressure oscillation with small size of the physical time step but leads to divergence with the pseudo-time step as well. As a remedy for the pseudo-time step, such methods are recommended as implicit methods and the local-time step method with a proposed modification of the time-marching MIM preventing accuracy loss due to very large time step size. Numerical experiments support the theoretical analyses and show the validity of the time-marching MIM proposed. [ABSTRACT FROM AUTHOR]

Abstract: We propose a level set method for the semiclassical limit of the Schrödinger equation with discontinuous potentials. The discontinuities in the potential corresponds to potential barriers, at which incoming waves can be partially transmitted and reflected. Previously such a problem was handled by Jin and Wen using the Liouville equation – which arises as the semiclassical limit of the Schrödinger equation – with an interface condition to account for partial transmissions and reflections (S. Jin, X. Wen, SIAM J. Num. Anal. 44 (2006) 1801–1828). However, the initial data are Dirac-delta functions which are difficult to approximate numerically with a high accuracy. In this paper, we extend the level set method introduced in (S. Jin, H. Liu, S. Osher, R. Tsai, J. Comp. Phys. 210 (2005) 497–518) for this problem. Instead of directly discretizing the Delta functions, our proposed method decomposes the initial data into finite sums of smooth functions that remain smooth in finite time along the phase flow, and hence can be solved much more easily using conventional high order discretization schemes. Two ideas are introduced here: (1) The solutions of the problems involving partial transmissions and partial reflections are decomposed into a finite sum of solutions solving problems involving only complete transmissions and those involving only complete reflections. For problems involving only complete transmission or complete reflection, the method of JLOT applies and is used in our simulations; (2) A reinitialization technique is introduced so that waves coming from multiple transmissions and reflections can be combined seamlessly as new initial value problems. This is implemented by rewriting the sum of several delta functions as one delta function with a suitable weight, which can be easily implemented numerically. We carry out numerical experiments in both one and two space dimensions to verify this new algorithm. [ABSTRACT FROM AUTHOR]

Abstract: We accelerate the computation of spherical harmonic transforms, using what is known as the butterfly scheme. This provides a convenient alternative to the approach taken in the second paper from this series on “Fast algorithms for spherical harmonic expansions”. The requisite precomputations become manageable when organized as a “depth-first traversal” of the program’s control-flow graph, rather than as the perhaps more natural “breadth-first traversal” that processes one-by-one each level of the multilevel procedure. We illustrate the results via several numerical examples. [Copyright &y& Elsevier]

Abstract: The numerical solution to the Biot equations of 3-D consolidation is still a challenging task because of the ill-conditioning of the resulting algebraic system and the instabilities that may affect the pore pressure solution. Recently new approaches have been advanced based on mixed formulations. In the present paper a fully coupled 3-D mixed finite element model is developed with the aim at alleviating the pore pressure numerical oscillations at the interface between materials with different permeabilities. A solution algorithm is implemented that takes advantage of the block structure of the discretized problem. The proposed model is verified against well-known analytical solutions and successfully experimented with in realistic applications of soil consolidation. [Copyright &y& Elsevier]

Abstract: In this paper, the solution of the Riemann Problem for the one-dimensional, free-surface Shallow Water Equations over a bed step is analyzed both from a theoretical and a numerical point of view. Particular attention has been paid to the wave that is generated at the location of the bed discontinuity. Starting from the classical Shallow Water Equations, considering the bed level as an additional variable, and adding to the system an equation imposing its time invariance, we show that this wave is a contact wave, across which one of the Riemann invariants, namely the energy, is not constant. This is due to the fact that the relevant problem is nonconservative. We demonstrate that, in this type of system, Riemann Invariants do not generally hold in contact waves. Furthermore, we show that in this case the equations that link the flow variables across the contact wave are the Generalized Rankine–Hugoniot relations and we obtain these for the specific problem. From the numerical point of view, we present an accurate and efficient solver for the step Riemann Problem to be used in a finite-volume Godunov-type framework. Through a two-step predictor–corrector procedure, the solver is able to provide solutions with any desired accuracy. The predictor step uses a well-balanced Generalized Roe solver while the corrector step solves the exact nonlinear system of equations that consitutes the problem by means of an iterative procedure that starts from the predictor solution. In order to show the effectiveness and the accuracy of the proposed approach, we consider several step Riemann Problems and compare the exact solutions with the numerical results obtained by using a standard Roe approach far from the step and the novel two-step algorithm for the fluxes over the step, achieving good results. [Copyright &y& Elsevier]

Abstract: Typical VOF algorithms rely on an implicit slip that scales with mesh refinement, to allow contact lines to move along no-slip boundaries. As a result, solutions of contact line phenomena vary continuously with mesh spacing; this paper presents examples of that variation. A mesh-dependent dynamic contact angle model is then presented, that is based on fundamental hydrodynamics and serves as a more appropriate boundary condition at a moving contact line. This new boundary condition eliminates the stress singularity at the contact line; the resulting problem is thus well-posed and yields solutions that converge with mesh refinement. Numerical results are presented of a solid plate withdrawing from a fluid pool, and of spontaneous droplet spread at small capillary and Reynolds numbers. [Copyright &y& Elsevier]

Abstract: The development of efficient algorithms to analyze complex electromagnetic structures is of topical interest. Application of these algorithms in commercial solvers requires rigorous error controllability. In this paper we focus on the perfectly matched layer based multilevel fast multipole algorithm (PML-MLFMA), a dedicated technique constructed to efficiently analyze large planar structures. More specifically the crux of the algorithm, viz. the pertinent layered medium Green functions, is under investigation. Therefore, particular attention is paid to the plane wave decomposition for 2-D homogeneous space Green functions in very lossy media, as needed in the PML-MLFMA. The result of the investigations is twofold. First, upper bounds expressing the required number of samples in the plane wave decomposition as a function of a preset accuracy are rigorously derived. These formulas can be used in 2-D homogeneous (lossy) media MLFMAs. Second, a more heuristic approach to control the error of the PML-MLFMA’s Green functions is presented. The theory is verified by means of several numerical experiments. [Copyright &y& Elsevier]

Abstract: In this paper a novel predictor–corrector algorithm is presented for the solution of coupled gas-phase – particulate systems. The emphasis of this work is the study of soot formation, but the concepts can be applied to other systems. This algorithm couples a stiff ODE solver to a Monte Carlo population balance solver. Such coupling has been achieved previously for similar systems using a Strang operator splitting algorithm, however, that algorithm demonstrated several numerical issues which resulted in a high computational cost to acquire adequate precision. In particular a source-sink instability was identified whereby a large-magnitude source term present in the ODE system was competing with a similarly sized sink term in the population balance. This instability required that the splitting step size was very small in order to keep numerical error sufficiently low. A predictor–corrector algorithm has been formulated to negate this instability. An additional efficiency is gained with this algorithm as a principal computational cost of the Strang splitting algorithm is removed: the requirement to re-initialise the ODE solver every splitting step. The numerical convergence of the new algorithm is demonstrated, and its efficiency is compared to that of the Strang splitting algorithm. Substantial computation time savings are demonstrated, which allow a fixed error in three studied system functionals to be achieved with an order-of-magnitude reduction in computation time. [Copyright &y& Elsevier]