Search

Abstract Next generation, conservative finite element hydrodynamic and transport models are vital for accurate and efficient ocean, estuary and river simulation. Numerical models such as these have been developed for decades by the U.S. Army Corps of Engineers at the Engineering Research and Development Center (ERDC). This paper focuses on recently developed implicit, multi-dimensional finite element 3D shallow water and transport models included in the Adaptive Hydraulics (AdH) numerical suite. These AdH 3D models benefit from their adaptive meshing capabilities to resolve sharp solution gradients, such as those often encountered with baroclinic wedges traveling up an estuary channel. This paper presents the AdH 3D mathematical formulation and solution procedure used to solve the weak finite element equations for these models, along with results for several common verification cases and an AdH Galveston Bay validation study. A novel Streamline Upwind Petrov–Galerkin (SUPG) method for 3D shallow water models is described which reduces to the AdH 2D shallow water SUPG formulation under certain conditions. Careful attention is placed on ensuring discrete consistency in the equation set. [ABSTRACT FROM AUTHOR]

The construction of local reduced-order models via multiscale basis functions has been an area of active research. In this paper, we propose online multiscale basis functions which are constructed using the offline space and the current residual. Online multiscale basis functions are constructed adaptively in some selected regions based on our error indicators. We derive an error estimator which shows that one needs to have an offline space with certain properties to guarantee that additional online multiscale basis function will decrease the error. This error decrease is independent of physical parameters, such as the contrast and multiple scales in the problem. The offline spaces are constructed using Generalized Multiscale Finite Element Methods (GMsFEM). We show that if one chooses a sufficient number of offline basis functions, one can guarantee that additional online multiscale basis functions will reduce the error independent of contrast. We note that the construction of online basis functions is motivated by the fact that the offline space construction does not take into account distant effects. Using the residual information, we can incorporate the distant information provided the offline approximation satisfies certain properties. In the paper, theoretical and numerical results are presented. Our numerical results show that if the offline space is sufficiently large (in terms of the dimension) such that the coarse space contains all multiscale spectral basis functions that correspond to small eigenvalues, then the error reduction by adding online multiscale basis function is independent of the contrast. We discuss various ways computing online multiscale basis functions which include a use of small dimensional offline spaces. [ABSTRACT FROM AUTHOR]

Multidimensional Riemann solvers that have internal sub-structure in the strongly-interacting state have been formulated recently (D.S. Balsara (2012, 2014) [5,16] ). Any multidimensional Riemann solver operates at the grid vertices and takes as its input all the states from its surrounding elements. It yields as its output an approximation of the strongly interacting state, as well as the numerical fluxes. The multidimensional Riemann problem produces a self-similar strongly-interacting state which is the result of several one-dimensional Riemann problems interacting with each other. To compute this strongly interacting state and its higher order moments we propose the use of a Galerkin-type formulation to compute the strongly interacting state and its higher order moments in terms of similarity variables. The use of substructure in the Riemann problem reduces numerical dissipation and, therefore, allows a better preservation of flow structures, like contact and shear waves. In this second part of a series of papers we describe how this technique is extended to unstructured triangular meshes. All necessary details for a practical computer code implementation are discussed. In particular, we explicitly present all the issues related to computational geometry. Because these Riemann solvers are Mu ltidimensional and have S elf-similar strongly- I nteracting states that are obtained by C onsistency with the conservation law, we call them MuSIC Riemann solvers. (A video introduction to multidimensional Riemann solvers is available on http://www.nd.edu/~dbalsara/Numerical-PDE-Course .) The MuSIC framework is sufficiently general to handle general nonlinear systems of hyperbolic conservation laws in multiple space dimensions. It can also accommodate all self-similar one-dimensional Riemann solvers and subsequently produces a multidimensional version of the same. In this paper we focus on unstructured triangular meshes. As examples of different systems of conservation laws we consider the Euler equations of compressible gas dynamics as well as the equations of ideal magnetohydrodynamics (MHD). Several stringent test problems are solved for both PDE systems. [ABSTRACT FROM AUTHOR]

This paper concerns the stability of a class of Helmholtz problems in rectangular domains. A well known application is the electromagnetic scattering from a rectangular cavity embedded in an infinite ground plane. Error analysis of numerical methods for cavity problems relies heavily on the stability estimates. However, it is extremely difficult to derive an optimal stability bound with the explicit dependency on wave numbers. In this paper a high-order finite element approximation is proposed for calculating the stability bound. Numerical experiments show that the stability depends strongly on wave numbers in extreme case and it is almost independent on the wave numbers in an average sense. Our numerical results also help to understand the stability of the multi-frequency inverse problems. [ABSTRACT FROM AUTHOR]

Abstract: In this paper, we approximate the semi-classical Schrödinger equation in the presence of electromagnetic field by the Hagedorn wave packets approach. By operator splitting, the Hamiltonian is divided into the modified part and the residual part. The modified Hamiltonian, which is the main new idea of this paper, is chosen by the fact that Hagedorn wave packets are localized both in space and momentum so that a crucial correction term is added to the truncated Hamiltonian, and is treated by evolving the parameters associated with the Hagedorn wave packets. The residual part is treated by a Galerkin approximation. We prove that, with the modified Hamiltonian only, the Hagedorn wave packets dynamics give the asymptotic solution with error , where ε is the scaled Planck constant. We also prove that, the Galerkin approximation for the residual Hamiltonian can reduce the approximation error to , where k depends on the number of Hagedorn wave packets added to the dynamics. This approach is easy to implement, and can be naturally extended to the multidimensional cases. Unlike the high order Gaussian beam method, in which the non-constant cut-off function is necessary and some extra error is introduced, the Hagedorn wave packets approach gives a practical way to improve accuracy even when ε is not very small. [Copyright &y& Elsevier]

Abstract: In this paper we present a new meshless method to solve well-known problems in physics and engineering with either constant or variable material properties in 2D space. Harmonic and bi-harmonic problems are considered in this paper. The method constructs a set of basis functions, called Equilibrated Basis Functions (EqBFs), through a weighted residual integration over a fictitious domain embedding the main one. The bases can satisfy the governing partial differential equations approximately. Chebyshev polynomials of the first kind are employed to construct the EqBFs and exponential functions are used as the weights in the integrals. The parameters of the solution are arranged so that all the integrals can be decomposed into much simpler 1D ones over a normalized intervals. This reduces the computational efforts significantly. Either the EqBFs or the results of the integration process may be stored for further use. The validity of the results is examined through an extended patch test. A set of physical sample problems with variable/non-variable material properties; as the potential flow over a cylinder, steady-state heat conduction problems in an anisotropic inhomogeneous functionally graded material, potential and stokes flow through a channel of finite width obstructed by periodic cylinders, and the bending of thin elastic plates having constant or variable thickness are solved to demonstrate the capabilities of the method. As a preliminary study, we show that the method may effectively be used in a domain decomposition approach. [Copyright &y& Elsevier]

Highlights • Hamilton–Jacobi equation. • Kernel based scheme. • Unconditionally stable. • High order accuracy. • Weighted essentially non-oscillatory methodology. Abstract In this paper, a class of high order numerical schemes is proposed for solving Hamilton–Jacobi (H–J) equations. This work is regarded as an extension of our previous work for nonlinear degenerate parabolic equations, see Christlieb et al. [14] , which relies on a special kernel-based formulation of the solutions and successive convolution. When applied to the H–J equations, the newly proposed scheme attains genuinely high order accuracy in both space and time, and more importantly, it is unconditionally stable, hence allowing for much larger time step evolution compared with other explicit schemes and saving computational cost. A high order weighted essentially non-oscillatory methodology and a novel nonlinear filter are further incorporated to capture the correct viscosity solution. Furthermore, by coupling the recently proposed inverse Lax–Wendroff boundary treatment technique, this method is very flexible in handing complex geometry as well as general boundary conditions. We perform numerical experiments on a collection of numerical examples, including H–J equations with linear, nonlinear, convex or non-convex Hamiltonians. The efficacy and efficiency of the proposed scheme in approximating the viscosity solution of general H–J equations is verified. [ABSTRACT FROM AUTHOR]

Abstract We propose a new tool for the study of the parabolic p -curl problem modelling the magnetic field in a high-temperature superconductor (HTS), a problem involving an operator ∇ × (| ∇ × H | p − 2 ∇ × H) analogous to the p -Laplacian. This so-called nonlinear p -curl problem is computationally expensive in three space dimensions because sharp fronts in the magnetic field develop near the surface of HTS. The proposed technique is a relaxation model of the p -curl problem which is shown to provide a consistent and stable approximation according to first, an inner/outer layer analysis and secondly, according to a time scale expansion. The relaxation model leads to monotone approximations of the front because of the presence of anisotropic diffusion. Most importantly, it can be used to develop new schemes that inherit the properties of the relaxation model. In order to demonstrate the relevance of the relaxation model, this paper applies the approximation to the construction of a second-order nonlinear finite-difference time-domain (FDTD) method, similar to Yee's scheme. It is shown that the natural discretization of the p -curl inspired by Yee's scheme leads to an unstable scheme, while a discretization based on the relaxation model is stable. We propose two numerical discretizations that are respectively first and second order in time, but both second-order in space. The linear and nonlinear stability of the schemes is studied and it is demonstrated that the schemes are both discrete divergence free. We verify the numerical schemes using manufactured solutions and Mayergoyz's moving-front solution in both one and two space dimensions, as well as on a magnetization of a typical HTS. [ABSTRACT FROM AUTHOR]

Abstract The Harten–Lax–van Leer with contact (HLLC) approximate Riemann solver is known to be plagued by various forms of numerical shock instability. In this paper, we propose a new framework for developing shock stable versions of the HLLC scheme for the Euler system of equations without compromising on its contact and shear preserving ability. The proposed framework, termed as S elective W ave M odification (SWM), identifies and enhances the magnitude of the inherent HLL-type diffusive component within the HLLC scheme in the vicinity of a normal shock front. This is achieved by suitably increasing the magnitudes of the nonlinear wave speed estimates appearing in them through a dissipation parameter ϵ. A linear perturbation analysis is performed to gauge the effectiveness of the proposed framework in damping unwanted perturbations in physical quantities and is also used to estimate a von-Neumann-type stability bound on the CFL number associated with its use. Two distinct methods to estimate ϵ , which results in HLLC-SWM-E (Eigenvalue based) and HLLC-SWM-P (Pressure based) schemes, are proposed and compared. Further, a matrix based stability analysis is used to show that the proposed schemes can be configured to remain shock stable over a wide range of free stream Mach numbers. Numerical results confirm that the proposed schemes are capable of computing shock stable solutions for a wide variety of problems. On viscous flows, the HLLC-SWM-P variant is found to be the most accurate. Performance of these schemes suggests that it is possible to construct shock stable upwind schemes that simultaneously retains accuracy on the linearly degenerate wavefields. Highlights • A simple strategy to cure shock instability in the HLLC Riemann solver proposed. • Cure achieved by enhancing the inherent HLL-type diffusive component. • Two shock stable variants namely HLLC-SWM-E and HLLC-SWM-P are developed. • Robustness and accuracy of the scheme are demonstrated on several numerical examples. [ABSTRACT FROM AUTHOR]

Highlights • Accurate and efficient numerical simulation of granular flow. • Two regularized models of the μ (I) -rheology. • Application of the PFEM to frictional material simulation. • Validation against 3D experimental tests. Abstract This paper proposes two regularized models of the μ (I) -rheology and shows their application to the numerical simulation of 3D dense granular flows. The proposed regularizations are inspired by the Papanastasiou and Bercovier–Engleman methods, typically used to approximate the Bingham law. The key idea is to keep limited the value of the apparent viscosity for low shear rates without introducing a fixed cutoff. The proposed techniques are introduced into the Particle Finite Element Method (PFEM) framework to deal with the large deformations expected in free-surface granular flows. After showing the numerical drawbacks associated to the standard μ (I) -rheology, the two regularization strategies are derived and discussed. The regularized μ (I) -rheology is then applied to the simulation of the collapse of 2D and 3D granular columns. The numerical results show that the regularization techniques improve substantially the conditioning of the linear system without affecting the solution accuracy. A good agreement with the experimental tests and other numerical methods is obtained in all the analyzed problems. [ABSTRACT FROM AUTHOR]

Abstract This paper focuses on the question of how unconditional stability can be achieved via multistep ImEx schemes, in practice problems where both the implicit and explicit terms are allowed to be stiff. For a class of new ImEx multistep schemes that involve a free parameter, strategies are presented on how to choose the ImEx splitting and the time stepping parameter, so that unconditional stability is achieved under the smallest approximation errors. These strategies are based on recently developed stability concepts, which also provide novel insights into the limitations of existing semi-implicit backward differentiation formulas (SBDF). For instance, the new strategies enable higher order time stepping that is not otherwise possible with SBDF. With specific applications in nonlinear diffusion problems and incompressible channel flows, it is demonstrated how the unconditional stability property can be leveraged to efficiently solve stiff nonlinear or nonlocal problems without the need to solve nonlinear or nonlocal problems implicitly. Highlights • Unconditional stability for equations with both implicit and explicit stiff terms. • Explanations of stability limitations in semi-implicit backward differentiation. • Strategies for the simultaneous choice of splitting and time stepping coefficients. • Examples avoiding nonlinear implicit solves in stiff nonlinear diffusion equations. [ABSTRACT FROM AUTHOR]

Abstract In this paper we analyze a multilevel quasidiffusion (QD) method for solving time-dependent multigroup nonlinear radiative transfer problems which describe interaction of photons with matter. The multilevel method is formulated by means of the high-order radiative transfer equation and a set of low-order moment equations. The fully implicit scheme is used to discretize equations in time. The stability analysis is applied to the method in semi-continuous and discretized forms. To perform Fourier analysis, the system of equations of the multilevel method is linearized about an equilibrium solution. The effects of discretization with respect to different independent variables are studied. The multilevel method is shown to be stable and fast converging. We also consider a version of the method in which time evolution in the radiative transfer equation is treated by means of the α -approximation. The Fleck–Cummings test problem is used to demonstrate performance of the multilevel QD method and study its iterative stability. Highlights • Fast multilevel iterative method for solving time-dependent nonlinear radiative transfer problems. • Fourier analysis of the multilevel quasidiffusion method in semi-continuous and discretized forms. • Performance and stability of the multilevel method are analyzed on a model multigroup thermal radiative transfer problem. [ABSTRACT FROM AUTHOR]

Abstract In the interpolation method, in some cases, one often has a number of data points and its derivatives, which are obtained by sampling or experimentation. In this case, the problem of finding an approximating function passing through these points and coinciding with given values of its derivatives at these points is generally known as "Hermite interpolation". The Hermite interpolation is mostly a method of interpolating data points as a polynomial function that is faced with some challenges in high dimensions and on irregular domains. Radial basis functions take advantage of being flexible with respect to geometry, easy to implement in higher dimensions, and can also provide high order convergence. So, they can be applied as a suitable tool to high dimensional Hermite interpolation problem on irregular domains. In many applied systems, commonly available derivatives information is presented using fractional order derivatives instead of integer ones. For this purpose, in this paper, we assume that the values of an unknown function and its fractional derivatives at some distinct points are presented. Therefore, we intend to apply a new approach, which we call it as "fractional Hermite interpolation" with radial basis functions in high dimensions. Optimal recovery conditions for the fractional Hermite interpolant are investigated. Then, the existence and uniqueness of the solution in this type of generalized interpolation are proved. In order to increase the accuracy and stability of the method, Hilbert Schmidt's theory has also been used. The main advantages of the used method are its simplicity and efficiency in high dimensions, and over irregular domains. Finally, numerical results in one, two and three dimensions and a real-world problem are presented to support our theoretical analysis. [ABSTRACT FROM AUTHOR]

Abstract While high order methods became very popular as they allow to perform very accurate solutions with low computational time and memory cost, there is a lack of tools to visualize and post-treat the solutions given by these methods. Originally, visualization softwares were developed to post-process results from methods such that finite differences or usual finite elements and therefore process linear primitives. In this paper, we present a methodology to visualize results of high order methods. Our approach is based on the construction of an optimized affine approximation of the high order solution which can therefore be handled by any visualization software. A representation mesh is constructed and the process is guided by an a posteriori estimate which control the error between the numerical solution and its representation pointwise. This point by point control is crucial as under their picture form, data correspond to values mapped on elements where anyone can pick up a pointwise information. A strategy is established to ensure that discontinuities are well represented. These discontinuities come either from the physical problem (material change) or the numerical method (Discontinuous Galerkin method) and are pictured accurately. Several numerical examples are presented to demonstrate the potential of the method. Highlights • Objectives to faithfully represent high order solutions are defined. • An algorithm to obtain the representation is given. • Under assumptions, the representation is proved to satisfy the objectives. • Specific developments for approximation spaces based on curved elements are led. • Relevancy and benefits are outlined with several numerical examples. [ABSTRACT FROM AUTHOR]

Abstract In this paper, we consider numerical methods for solving the two-dimensional Vlasov–Poisson equation in the finite Larmor radius approximation regime. The model describes the behaviour of charged particles under a strong external magnetic field and the finite Larmor radius scaling. We discretise the equation under Particle-in-Cell method, where the characteristics equations are highly oscillatory system in the limit regime. We apply popular numerical integrators including splitting methods, multi-revolution composition methods, two-scale formulation method and limit solver to integrate the characteristics. We then highlight the strengths and drawbacks of each method. Finally, numerical experiments are presented, and comparisons on the accuracy, efficiency and long-time behaviour of the methods are made, pointing to the method with the best performance for this problem. [ABSTRACT FROM AUTHOR]

Abstract Quadrature by Expansion (QBX) is a quadrature method for approximating the value of the singular integrals encountered in the evaluation of layer potentials. It exploits the smoothness of the layer potential by forming locally-valid expansions which are then evaluated to compute the near or on-surface value of the potential. Recent work towards coupling of a Fast Multipole Method (FMM) to QBX yielded a first step towards the rapid evaluation of such integrals (and the solution of related integral equations), albeit with only empirically understood error behavior. In this paper, we improve upon this approach with a modified algorithm for which we give a comprehensive analysis of error and cost in the case of the Laplace equation in two dimensions. For the same levels of (user-specified) accuracy, the new algorithm empirically has cost-per-accuracy comparable to prior approaches. We provide experimental results to demonstrate scalability and numerical accuracy. Highlights • A new fast multipole algorithm for Quadrature by Expansion (QBX). • Comprehensive error and cost analysis for the case of the 2D Laplace kernel. • Empirical results match modeled error and scaling behavior. [ABSTRACT FROM AUTHOR]

Abstract In this paper, an advanced numerical technique is proposed for a low-frequency electromagnetic field simulation. The finite difference time domain (FDTD) method is one of the best methods for analyzing electromagnetic field problems, which require high levels of precision, such as a heterogeneous whole-body voxel human model. However, directly using the standard FDTD method in the quasi-static frequency band requires an extra-large number of iterations. We overcome this drawback using the hybrid of surface equivalence theorem and quasi-static FDTD (QS-FDTD) method. The surface equivalence theorem is used to consider arbitrary shape of source and the QS-FDTD is applied to obtain electromagnetic response quickly in the low-frequency region. The results of the proposed method are in good agreement with those from a commercial electromagnetic simulator. The proposed method is valid where the quasi-static approximation is satisfied. [ABSTRACT FROM AUTHOR]

Finite Element Methods (FEM) are becoming increasingly popular in modeling seismic wave propagation. These methods provide higher order accuracy, geometrical flexibility and adaptive gridding capabilities that are not easy to incorporate in traditional finite difference methods employed for generation of synthetic seismograms. Moreover, several studies have shown that Discontinuous Galerkin FEM (DGM) is a promising approach for modeling wave propagation in fractured media. Here we propose an Enriched Galerkin FEM (EGM) for elastic wave propagation. EGM uses the same bilinear form as DGM and the continuous Galerkin finite element spaces enriched by discontinuous piecewise constants or bilinear functions. EGM satisfies local equilibrium while reducing the degrees of freedom in DGM formulations. In this paper, we consider elastic wave propagation and derive optimal a priori error estimates for DGM and EGM. We present numerical examples in two spatial dimensions that confirm these theoretical results. In addition, we provide numerical comparisons with the Spectral element method. In previous work, DGM has been shown to be effective in modeling elastic wave propagation in fractured media using the linear slip model. We now extend these results to EGM with reduced computational costs over DGM. [ABSTRACT FROM AUTHOR]

We are interested in the approximation of a steady hyperbolic problem. In some cases, the solution can satisfy an additional conservation relation, at least when it is smooth. This is the case of an entropy. In this paper, we show, starting from the discretization of the original PDE, how to construct a scheme that is consistent with the original PDE and the additional conservation relation. Since one interesting example is given by the systems endowed by an entropy, we provide one explicit solution, and show that the accuracy of the new scheme is at most degraded by one order. In the case of a discontinuous Galerkin scheme and a Residual distribution scheme, we show how not to degrade the accuracy. This improves the recent results obtained in [1–4] in the sense that no particular constraints are set on quadrature formula and that a priori maximum accuracy can still be achieved. We study the behavior of the method on a non linear scalar problem. However, the method is not restricted to scalar problems. [ABSTRACT FROM AUTHOR]

We develop new approximations of the Mori–Zwanzig equation based on operator series expansions of the orthogonal dynamics propagator. In particular, we study the Faber polynomial series, which yields asymptotically optimal approximations converging at least R -superlinearly with the polynomial order. Numerical applications are presented for random wave propagation and harmonic chains of oscillators interacting on the Bethe lattice and on graphs with arbitrary topology. Most of the analysis and numerical results presented in this paper are for linear dynamical systems evolving from random initial states. The extension to nonlinear systems is discussed as well the challenges in the implementation of the proposed approximation methods. [ABSTRACT FROM AUTHOR]

For a large class of orthogonal basis functions, there has been a recent identification of expansion methods for computing accurate, stable approximations of a quantity of interest. This paper presents, within the context of uncertainty quantification, a practical implementation using basis adaptation, and coherence motivated sampling, which under assumptions has satisfying guarantees. This implementation is referred to as Basis Adaptive Sample Efficient Polynomial Chaos (BASE-PC). A key component of this is the use of anisotropic polynomial order which admits evolving global bases for approximation in an efficient manner, leading to consistently stable approximation for a practical class of smooth functionals. This fully adaptive, non-intrusive method, requires no a priori information of the solution, and has satisfying theoretical guarantees of recovery. A key contribution to stability is the use of a presented correction sampling for coherence-optimal sampling in order to improve stability and accuracy within the adaptive basis scheme. Theoretically, the method may dramatically reduce the impact of dimensionality in function approximation, and numerically the method is demonstrated to perform well on problems with dimension up to 1000. [ABSTRACT FROM AUTHOR]

The presence of a high-dimensional stochastic input domain with discontinuities poses major computational challenges in analyzing and quantifying the effects of the uncertainties in a physical system. In this paper, we propose a stochastic collocation method with adaptive mesh refinement (SCAMR) to deal with high dimensional stochastic systems with discontinuities. Specifically, the proposed approach uses generalized polynomial chaos (gPC) expansion with Legendre polynomial basis and solves for the gPC coefficients using the least squares method. It also implements an adaptive mesh (element) refinement strategy which checks for abrupt variations in the output based on a low-order gPC approximation error to track discontinuities or non-smoothness. In addition, the proposed method involves a criterion for checking possible dimensionality reduction and consequently, the decomposition of the original high-dimensional problem to a number of lower-dimensional subproblems. Specifically, this criterion checks all the existing interactions between input parameters of a specific problem based on the high-dimensional model representation (HDMR) method, and therefore automatically provides the subproblems which only involve interacting input parameters. The efficiency of the approach is demonstrated using examples of both smooth and non-smooth problems with number of input parameters up to 500, and the approach is compared against other existing algorithms. [ABSTRACT FROM AUTHOR]

The aim of this paper is to compute the mean field and variance of solutions to three-dimensional Maxwell's equations with random interfaces via shape calculus and pivoted low-rank approximation. Based on the perturbation theory and shape calculus, we characterize the statistical moments of solutions to Maxwell's equations with random interfaces in terms of the perturbation magnitude via the first order shape-Taylor expansion. In order to capture oscillations with high resolution close to the interface, an adaptive finite element method using Nédélec's third order edge elements of the first kind is employed to solve the deterministic Maxwell's equations with the mean interface to approximate the expectation of solutions. For the second moment computation, an efficient low-rank approximation of the pivoted Cholesky decomposition is proposed to compute the two-point correlation function to approximate the variance of solutions. Numerical experiments are presented to demonstrate our theoretical results. [ABSTRACT FROM AUTHOR]

This paper serves to treat boundary conditions numerically with high order accuracy in order to suit the two-stage fourth-order finite volume schemes for hyperbolic problems developed in J. Li and Z. Du (2016) [17] . As such, it is significant when capturing small scale structures near physical boundaries. Different from previous contributions in literature, the current approach constructs a fourth order accurate approximation to boundary conditions by only using the Jacobian matrix of the flux function (characteristic information) instead of its successive differentiation of governing equations leading to tensors of high ranks in the inverse Lax–Wendroff method. Technically, data in several ghost cells are constructed with interpolation so that the interior scheme can be implemented over boundary cells, and theoretical boundary condition has to be modified properly at intermediate stages so as to make the two-stage scheme over boundary cells fully consistent with that over interior cells. This is nonintuitive and highlights the fact that theoretical boundary conditions are only prescribed for continuous partial differential equations (PDEs), while they must be approximated in a consistent way (even though they could be exactly valued) when the PDEs are discretized . Several numerical examples are provided to illustrate the performance of the current approach when dealing with general boundary conditions. [ABSTRACT FROM AUTHOR]

Reconstructing the interface within a cell, based on volume fraction and normal direction, is a key part of multiphase flow solvers which make use of piecewise linear interface calculation (PLIC) such as the Coupled Level Set Volume of Fluid (CLSVOF) method. In this paper, we present an iterative method for interface reconstruction (IR) in general convex cells based on tetrahedral decomposition. By splitting the cell into tetrahedra prior to IR the volume of the truncated polyhedron can be calculated much more rapidly than using existing clipping and capping methods. In addition the root finding algorithm is designed to take advantage of the nature of the relationship between volume fraction and interface position by using a combination of Newton's and Muller's methods. In stand-alone tests of the IR algorithm on single cells with up to 20 vertices the proposed method was found to be 2 times faster than an implementation of an existing analytical method, while being easy to implement. It was also found to be 3.4–11.8 times faster than existing iterative methods using clipping and capping and combined with Brent's root finding method. Tests were then carried out of the IR method as part of a CLSVOF solver. For a sphere deformed by a prescribed velocity field the proposed method was found to be up to 33% faster than existing iterative methods. For simulations including the solution of the velocity field the maximum speed up was found to be approximately 52% for a case where 12% of cells lie on the interface. Analysis of the full simulation CPU time budget also indicates that while the proposed method has produced a considerable speed-up, further gains due to increasing the efficiency of the IR method are likely to be small as the IR step now represents only a small proportion of the run time. [ABSTRACT FROM AUTHOR]

The asymptotic, physical optics (PO) approximation is applicable to electromagnetic scattering analysis of electrically large, conducting objects, which is important for various applications, such as radar cross-section (RCS) calculations. It relates the incident magnetic field to the induced surface current at points on a scatterer's surface with line-of-sight visibility to the source. Here, the focus is on mesh-based PO, using a triangle element mesh with standard Rao–Wilton–Glisson basis functions to represent the surface current density. The multiple-reflection PO (MRPO) method involves applying the PO approximation to the field radiated by the present current solution towards the scatterer itself, to account for successive internal reflections. The visibility status between every pair of basis functions is required, taking into account all geometry. In a conventional implementation, internal reflected field calculation and internal shadowing determination runtimes both scale as O ( N 2 ) ; N denotes the number of mesh elements. In this paper an accelerated version of the mesh-based MRPO method is presented: fast MRPO (FMRPO). The multi-level, fast multipole method (MLFMM) is used to accelerate internal reflected field calculation. The key aspect to achieving acceleration for general geometries, is to incorporate internal shadowing into the MLFMM interaction tree. This requires an alteration to the inter-group interaction criterion of the MLFMM, as well as efficient evaluation of inter-group shadowing status flags, to preserve the beneficial cost-scaling property of the MLFMM. Algorithmic parameters are introduced, which control the accuracy of internal shadowing determination. The FMRPO runtime scales as O(NlogN), in the high-frequency limit. Numerical results are presented, demonstrating accuracy and efficiency for general scattering objects. [ABSTRACT FROM AUTHOR]

In this paper we develop the a posteriori error analysis of an augmented mixed-primal finite element method for the 2D and 3D versions of a stationary flow and transport coupled system, typically encountered in sedimentation–consolidation processes. The governing equations consist in the Brinkman problem with concentration-dependent viscosity, written in terms of Cauchy pseudo-stresses and bulk velocity of the mixture; coupled with a nonlinear advection – nonlinear diffusion equation describing the transport of the solids volume fraction. We derive two efficient and reliable residual-based a posteriori error estimators for a finite element scheme using Raviart–Thomas spaces of order k for the stress approximation, and continuous piecewise polynomials of degree ≤ k + 1 for both velocity and concentration. For the first estimator we make use of suitable ellipticity and inf–sup conditions together with a Helmholtz decomposition and the local approximation properties of the Clément interpolant and Raviart–Thomas operator to show its reliability, whereas the efficiency follows from inverse inequalities and localisation arguments based on triangle-bubble and edge-bubble functions. Next, we analyse an alternative error estimator, whose reliability can be proved without resorting to Helmholtz decompositions. Finally, we provide some numerical results confirming the reliability and efficiency of the estimators and illustrating the good performance of the associated adaptive algorithm for the augmented mixed-primal finite element method. [ABSTRACT FROM AUTHOR]

The advection equation is the basis for mathematical models of continuum mechanics. In the approximate solution of nonstationary problems it is necessary to inherit main properties of the conservatism and monotonicity of the solution. In this paper, the advection equation is written in the symmetric form, where the advection operator is the half-sum of advection operators in conservative (divergent) and non-conservative (characteristic) forms. The advection operator is skew-symmetric. Standard finite element approximations in space are used. The standard explicit two-level scheme for the advection equation is absolutely unstable. New conditionally stable regularized schemes are constructed, on the basis of the general theory of stability (well-posedness) of operator-difference schemes, the stability conditions of the explicit Lax–Wendroff scheme are established. Unconditionally stable and conservative schemes are implicit schemes of the second (Crank–Nicolson scheme) and fourth order. The conditionally stable implicit Lax–Wendroff scheme is constructed. The accuracy of the investigated explicit and implicit two-level schemes for an approximate solution of the advection equation is illustrated by the numerical results of a model two-dimensional problem. [ABSTRACT FROM AUTHOR]

This paper presents two discrete maximum principles (DMP) for the numerical solution of second order elliptic equations arising from the weak Galerkin finite element method. The results are established by assuming an h -acute angle condition for the underlying finite element triangulations. The mathematical theory is based on the well-known De Giorgi technique adapted in the finite element context. Some numerical results are reported to validate the theory of DMP. [ABSTRACT FROM AUTHOR]

The Bloch mode synthesis (BMS) model-reduction method adapts component mode synthesis techniques to unit-cell problems in order to obtain a reduced-order model that quickly produces band-structure frequencies for any wave vector, or vice versa. Fundamental to BMS is a partitioning of the real-space model into interior and boundary components, and subsequent reduction of the interior via truncated normal mode expansion. In this paper, two enhancements are presented for the BMS method that reduce both computation time and error in band-structure calculations. The first enhancement improves the accuracy of the interior reduction by approximating the participation of the residual modes rather than simply truncating them. The original formulation of BMS includes a modal reduction of the boundary that must be recomputed for every wave vector. This limits computational benefits and prevents the reduced-order model from being useful for the inverse band-structure problem (i.e., the k ( ω ) calculation). The second enhancement is a local boundary reduction that is independent of wave vector and thus does not suffer from the aforementioned limitations. [ABSTRACT FROM AUTHOR]

A new algorithm for numerical evaluation of spherical Bessel transform is proposed in this paper. In this method, the spherical Bessel function is approximately represented as an exponential sum with complex parameters. This is obtained by expressing an integral representation of spherical Bessel function in complex plane, and discretizing contour integrals along steepest descent paths and a contour path parallel to real axis using numerical quadrature rule with the double-exponential transformation. The number of terms in the expression is reduced using the modified balanced truncation method. The residual part of integrand is also expanded by exponential functions using Prony-like method. The spherical Bessel transform can be evaluated analytically on arbitrary points in half-open interval. [ABSTRACT FROM AUTHOR]

This paper describes a high order accurate method to calculate integrals over curved surfaces with boundaries. Given data locations that are arbitrarily distributed over the surface, together with some functional description of the surface and its boundary, the algorithm produces matching quadrature weights. This extends on the authors' earlier methods for integrating over the surface of a sphere and over arbitrarily shaped smooth closed surfaces by also considering domain boundaries. The core approach consists again of combining RBF-FD (radial basis function-generated finite difference) approximations for curved surface triangles, which together make up the full surface. The provided examples include both curved and flat domains. In the highly special case of equi-spaced nodes over a regular interval in 1-D, the method provides a new opportunity for improving on the classical Gregory enhancements of the trapezoidal rule. [ABSTRACT FROM AUTHOR]

In this paper, we develop a new tensor-product based preconditioner for discontinuous Galerkin methods with polynomial degrees higher than those typically employed. This preconditioner uses an automatic, purely algebraic method to approximate the exact block Jacobi preconditioner by Kronecker products of several small, one-dimensional matrices. Traditional matrix-based preconditioners require O ( p 2 d ) storage and O ( p 3 d ) computational work, where p is the degree of basis polynomials used, and d is the spatial dimension. Our SVD-based tensor-product preconditioner requires O ( p d + 1 ) storage, O ( p d + 1 ) work in two spatial dimensions, and O ( p d + 2 ) work in three spatial dimensions. Combined with a matrix-free Newton–Krylov solver, these preconditioners allow for the solution of DG systems in linear time in p per degree of freedom in 2D, and reduce the computational complexity from O ( p 9 ) to O ( p 5 ) in 3D. Numerical results are shown in 2D and 3D for the advection, Euler, and Navier–Stokes equations, using polynomials of degree up to p = 30 . For many test cases, the preconditioner results in similar iteration counts when compared with the exact block Jacobi preconditioner, and performance is significantly improved for high polynomial degrees p . [ABSTRACT FROM AUTHOR]

We present a randomized iterative method for approximating unknown function sequentially on arbitrary point set. The method is based on a recently developed sequential approximation (SA) method, which approximates a target function using one data point at each step and avoids matrix operations. The focus of this paper is on data sets with highly irregular distribution of the points. We present a nearest neighbor replacement (NNR) algorithm, which allows one to sample the irregular data sets in a near optimal manner. We provide mathematical justification and error estimates for the NNR algorithm. Extensive numerical examples are also presented to demonstrate that the NNR algorithm can deliver satisfactory convergence for the SA method on data sets with high irregularity in their point distributions. [ABSTRACT FROM AUTHOR]

In this paper we propose a method for the strong imposition of random Dirichlet boundary conditions in the Dynamical Low Rank (DLR) approximation of parabolic PDEs and, in particular, incompressible Navier Stokes equations. We show that the DLR variational principle can be set in the constrained manifold of all S rank random fields with a prescribed value on the boundary, expressed in low rank format, with rank smaller then S . We characterize the tangent space to the constrained manifold by means of a Dual Dynamically Orthogonal (Dual DO) formulation, in which the stochastic modes are kept orthonormal and the deterministic modes satisfy suitable boundary conditions, consistent with the original problem. The Dual DO formulation is also convenient to include the incompressibility constraint, when dealing with incompressible Navier Stokes equations. We show the performance of the proposed Dual DO approximation on two numerical test cases: the classical benchmark of a laminar flow around a cylinder with random inflow velocity, and a biomedical application for simulating blood flow in realistic carotid artery reconstructed from MRI data with random inflow conditions coming from Doppler measurements. [ABSTRACT FROM AUTHOR]

In this paper, we propose a Newton–Krylov solver and a Picard–Krylov solver for finite element discrete problem of stationary incompressible magnetohydrodynamic equations in three dimensions. Using a mixed finite element method, we discretize the velocity and the pressure by H 1 ( Ω ) -conforming finite elements and discretize the magnetic field by H ( curl , Ω ) -conforming edge elements. An efficient preconditioner is proposed to accelerate the convergence of GMRES method for solving linearized discrete problems. By extensive numerical experiments, we demonstrate the robustness of the Newton–Krylov solver for relatively large physical parameters and the optimality with respect to the number of degrees of freedom. Moreover, the numerical experiments show that the Newton–Krylov solver is more robust than the Picard–Krylov solver for large Reynolds number. [ABSTRACT FROM AUTHOR]

Hybrid polynomial correlated function expansion (H-PCFE) is a novel metamodel formulated by coupling polynomial correlated function expansion (PCFE) and Kriging. Unlike commonly available metamodels, H-PCFE performs a bi-level approximation and hence, yields more accurate results. However, till date, it is only applicable to medium scaled problems. In order to address this apparent void, this paper presents an improved H-PCFE, referred to as locally refined hp – adaptive H-PCFE. The proposed framework computes the optimal polynomial order and important component functions of PCFE, which is an integral part of H-PCFE, by using global variance based sensitivity analysis. Optimal number of training points are selected by using distribution adaptive sequential experimental design. Additionally, the formulated model is locally refined by utilizing the prediction error, which is inherently obtained in H-PCFE. Applicability of the proposed approach has been illustrated with two academic and two industrial problems. To illustrate the superior performance of the proposed approach, results obtained have been compared with those obtained using hp – adaptive PCFE. It is observed that the proposed approach yields highly accurate results. Furthermore, as compared to hp – adaptive PCFE, significantly less number of actual function evaluations are required for obtaining results of similar accuracy. [ABSTRACT FROM AUTHOR]

The present paper deals with the numerical solution of the incompressible Navier–Stokes equations using high-order discontinuous Galerkin (DG) methods for discretization in space. For DG methods applied to the dual splitting projection method, instabilities have recently been reported that occur for small time step sizes. Since the critical time step size depends on the viscosity and the spatial resolution, these instabilities limit the robustness of the Navier–Stokes solver in case of complex engineering applications characterized by coarse spatial resolutions and small viscosities. By means of numerical investigation we give evidence that these instabilities are related to the discontinuous Galerkin formulation of the velocity divergence term and the pressure gradient term that couple velocity and pressure. Integration by parts of these terms with a suitable definition of boundary conditions is required in order to obtain a stable and robust method. Since the intermediate velocity field does not fulfill the boundary conditions prescribed for the velocity, a consistent boundary condition is derived from the convective step of the dual splitting scheme to ensure high-order accuracy with respect to the temporal discretization. This new formulation is stable in the limit of small time steps for both equal-order and mixed-order polynomial approximations. Although the dual splitting scheme itself includes inf–sup stabilizing contributions, we demonstrate that spurious pressure oscillations appear for equal-order polynomials and small time steps highlighting the necessity to consider inf–sup stability explicitly. [ABSTRACT FROM AUTHOR]

This paper shows how to obtain boundary formulas for discrete-adjoint-based sensitivities in shape optimization problems. The analysis is carried out for inviscid flows in two-dimensional, unstructured triangular grids. The new formulation agrees with known continuous adjoint boundary formulas for the sensitivities, and the derivation sheds light into the required approximations and the reasons for the equivalence (or lack thereof) between boundary and domain-based adjoint sensitivities. Alternative approximations are also discussed. [ABSTRACT FROM AUTHOR]

Multilevel Monte Carlo (MLMC) is a recently proposed variation of Monte Carlo (MC) simulation that achieves variance reduction by simulating the governing equations on a series of spatial (or temporal) grids with increasing resolution. Instead of directly employing the fine grid solutions, MLMC estimates the expectation of the quantity of interest from the coarsest grid solutions as well as differences between each two consecutive grid solutions. When the differences corresponding to finer grids become smaller, hence less variable, fewer MC realizations of finer grid solutions are needed to compute the difference expectations, thus leading to a reduction in the overall work. This paper presents an extension of MLMC, referred to as multilevel control variates (MLCV) , where a low-rank approximation to the solution on each grid, obtained primarily based on coarser grid solutions, is used as a control variate for estimating the expectations involved in MLMC. Cost estimates as well as numerical examples are presented to demonstrate the advantage of this new MLCV approach over the standard MLMC when the solution of interest admits a low-rank approximation and the cost of simulating finer grids grows fast. [ABSTRACT FROM AUTHOR]

The particle–particle random phase approximation (pp-RPA) has been shown to be capable of describing double, Rydberg, and charge transfer excitations, for which the conventional time-dependent density functional theory (TDDFT) might not be suitable. It is thus desirable to reduce the computational cost of pp-RPA so that it can be efficiently applied to larger molecules and even solids. This paper introduces an O ( N 3 ) algorithm, where N is the number of orbitals, based on an interpolative separable density fitting technique and the Jacobi–Davidson eigensolver to calculate a few low-lying excitations in the pp-RPA framework. The size of the pp-RPA matrix can also be reduced by keeping only a small portion of orbitals with orbital energy close to the Fermi energy. This reduced system leads to a smaller prefactor of the cubic scaling algorithm, while keeping the accuracy for the low-lying excitation energies. [ABSTRACT FROM AUTHOR]

Several physical phenomena such as transformation of pollutants, energy, particles and many others can be described by the well-known convection–diffusion equation which is a combination of the diffusion and advection equations. In this paper, this equation is generalized with the concept of variable-order fractional derivatives. The generalized equation is called variable-order time fractional advection–diffusion equation (V-OTFA–DE). An accurate and robust meshless method based on the moving least squares (MLS) approximation and the finite difference scheme is proposed for its numerical solution on two-dimensional (2-D) arbitrary domains. In the time domain, the finite difference technique with a θ -weighted scheme and in the space domain, the MLS approximation are employed to obtain appropriate semi-discrete solutions. Since the newly developed method is a meshless approach, it does not require any background mesh structure to obtain semi-discrete solutions of the problem under consideration, and the numerical solutions are constructed entirely based on a set of scattered nodes. The proposed method is validated in solving three different examples including two benchmark problems and an applied problem of pollutant distribution in the atmosphere. In all such cases, the obtained results show that the proposed method is very accurate and robust. Moreover, a remarkable property so-called positive scheme for the proposed method is observed in solving concentration transport phenomena. [ABSTRACT FROM AUTHOR]

In this paper, we present an accurate and efficient wavelet-based adaptive weighted essentially non-oscillatory (WENO) scheme for hydrodynamics and ideal magnetohydrodynamics (MHD) equations arising from the hyperbolic conservation systems. The proposed method works with the finite difference weighted essentially non-oscillatory (FD-WENO) method in space and the third order total variation diminishing (TVD) Runge–Kutta (RK) method in time. The philosophy of this work is to use the lifted interpolating wavelets as not only detector for singularities but also interpolator. Especially, flexible interpolations can be performed by an inverse wavelet transformation. When the divergence cleaning method introducing auxiliary scalar field ψ is applied to the base numerical schemes for imposing divergence-free condition to the magnetic field in a MHD equation, the approximations to derivatives of ψ require the neighboring points. Moreover, the fifth order WENO interpolation requires large stencil to reconstruct high order polynomial. In such cases, an efficient interpolation method is necessary. The adaptive spatial differentiation method is considered as well as the adaptation of grid resolutions. In order to avoid the heavy computation of FD-WENO, in the smooth regions fixed stencil approximation without computing the non-linear WENO weights is used, and the characteristic decomposition method is replaced by a component-wise approach. Numerical results demonstrate that with the adaptive method we are able to resolve the solutions that agree well with the solution of the corresponding fine grid. [ABSTRACT FROM AUTHOR]

A fractional reaction–diffusion model with a moving boundary is presented in this paper. An efficient numerical method is constructed to solve this moving boundary problem. Our method makes use of a finite difference approximation for the temporal discretization, and spectral approximation for the spatial discretization. The stability and convergence of the method is studied, and the errors of both the semi-discrete and fully-discrete schemes are derived. Numerical examples, motivated by problems from developmental biology, show a good agreement with the theoretical analysis and illustrate the efficiency of our method. [ABSTRACT FROM AUTHOR]

In this paper, we propose and study two approaches to approximate the solution of the Bethe–Salpeter equation (BSE) by using structured iterative eigenvalue solvers. Both approaches are based on the reduced basis method and low-rank factorizations of the generating matrices. We also propose to represent the static screen interaction part in the BSE matrix by a small active sub-block, with a size balancing the storage for rank-structured representations of other matrix blocks. We demonstrate by various numerical tests that the combination of the diagonal plus low-rank plus reduced-block approximation exhibits higher precision with low numerical cost, providing as well a distinct two-sided error estimate for the smallest eigenvalues of the Bethe–Salpeter operator. The complexity is reduced to O ( N b 2 ) in the size of the atomic orbitals basis set, N b , instead of the practically intractable O ( N b 6 ) scaling for the direct diagonalization. In the second approach, we apply the quantized-TT (QTT) tensor representation to both, the long eigenvectors and the column vectors in the rank-structured BSE matrix blocks, and combine this with the ALS-type iteration in block QTT format. The QTT-rank of the matrix entities possesses almost the same magnitude as the number of occupied orbitals in the molecular systems, N o < N b , hence the overall asymptotic complexity for solving the BSE problem by the QTT approximation is estimated by O ( log ⁡ ( N o ) N o 2 ) . We confirm numerically a considerable decrease in computational time for the presented iterative approaches applied to various compact and chain-type molecules, while supporting sufficient accuracy. [ABSTRACT FROM AUTHOR]

In this paper we provide a new approach for constructing non-reflecting boundary conditions. The boundary conditions are based on summation-by-parts operators and derived without Laplace transformation in time. We prove that the new non-reflecting boundary conditions yield a well-posed problem and that the corresponding numerical approximation is unconditionally stable. The analysis is demonstrated on a hyperbolic system in two space dimensions, and the theoretical results are confirmed by numerical experiments. [ABSTRACT FROM AUTHOR]

This paper discusses the necessity and strategy to unify the development of a dynamic solver and a minimum action method (MAM) for a spatially extended system when employing the large deviation principle (LDP) to study the effects of small random perturbations. A dynamic solver is used to approximate the unperturbed system, and a minimum action method is used to approximate the LDP, which corresponds to solving an Euler–Lagrange equation related to but more complicated than the unperturbed system. We will clarify possible inconsistencies induced by independent numerical approximations of the unperturbed system and the LDP, based on which we propose to define both the dynamic solver and the MAM on the same approximation space for spatial discretization. The semi-discrete LDP can then be regarded as the exact LDP of the semi-discrete unperturbed system, which is a finite-dimensional ODE system. We achieve this methodology for the two-dimensional Navier–Stokes equations using a divergence-free approximation space. The method developed can be used to study the nonlinear instability of wall-bounded parallel shear flows, and be generalized straightforwardly to three-dimensional cases. Numerical experiments are presented. [ABSTRACT FROM AUTHOR]

In this paper we propose a new acoustic-convective splitting-based numerical scheme for the Kapila five-equation two-phase flow model. The splitting operator decouples the acoustic waves and convective waves. The resulting two submodels are alternately numerically solved to approximate the solution of the entire model. The Lagrangian form of the acoustic submodel is numerically solved using an HLLC-type Riemann solver whereas the convective part is approximated with an upwind scheme. The result is a simple method which allows for a general equation of state. Numerical computations are performed for standard two-phase shock tube problems. A comparison is made with a non-splitting approach. The results are in good agreement with reference results and exact solutions. [ABSTRACT FROM AUTHOR]

The aim of this paper is to build and validate some explicit high-order schemes, both in space and time, for simulating the dynamics of systems of nonlinear Schrödinger/Gross–Pitaevskii equations. The method is based on the combination of high-order IMplicit–EXplicit (IMEX) schemes in time and Fourier pseudo-spectral approximations in space. The resulting IMEXSP schemes are highly accurate, efficient and easy to implement. They are also robust when used in conjunction with an adaptive time stepping strategy and appear as an interesting alternative to time-splitting pseudo-spectral (TSSP) schemes. Finally, a complete numerical study is developed to investigate the properties of the IMEXSP schemes, in comparison with TSSP schemes, for one- and two-components systems of Gross–Pitaevskii equations. [ABSTRACT FROM AUTHOR]

We describe a new theoretical and numerical framework for magnetohydrodynamic (MHD) simulations with an incorporated anisotropic pressure tensor, which can play an important role in collisionless plasmas. The classical approach to handle the anisotropy is based on application of the double adiabatic approximation, assuming that the pressure tensor is well described only by those components that are oriented parallel and perpendicular to the local magnetic field. This gyrotropic assumption, however, fails around magnetically neutral regions, where the cyclotron period may become comparable to or even longer than the system's dynamical time, which causes a singularity in the mathematical expression. In this paper, we demonstrate that this singularity can be completely removed by direct use of the 2nd-moment of the Vlasov equation, combined with an ingenious gyrotropization model. Numerical tests are used to verify that our model properly reduces to the standard MHD results or the double adiabatic formulation in an asymptotic manner under the limit of fast isotropization and fast gyrotropization, respectively. [ABSTRACT FROM AUTHOR]