Higher Order Convergence Rates Research Articles

Implicit-explicit time integration for the immersed wave equation

Efficient solution strategies for wave propagation problems with complex geometries are essential in many engineering fields. Immersed boundary methods simplify mesh generation by embedding the domain of interest into an extended domain that is easy to mesh, introducing the challenge of dealing with cells that intersect the domain boundary. We consider the finite cell method that extends the weak form from the physical into the extended domain, multiplied by a small value, to stabilize badly cut cells. Combined with explicit time integration methods, the presence of finite cells with very little support in the physical domain results in tiny critical time step sizes. While the finite cell stabilization limits how small the critical time step size can become, this limit is still restrictive for many applications. Explicit transient analyses commonly use the spectral element method due to its natural way of obtaining diagonal mass matrices through nodal lumping. The resulting method is very efficient since nodal lumping renders the solution of the equation systems trivial while not reducing the accuracy or the critical time step size. The combination of the spectral element and finite cell methods is called the spectral cell method. Unfortunately, a direct application of nodal lumping in the spectral cell method is impossible due to the special quadrature necessary to treat the discontinuous integrand inside the cut cells. The existing approaches to lump the mass matrices of cut cells significantly reduce the accuracy of the approximation.We analyze an implicit-explicit (IMEX) time integration method to exploit the advantages of the nodal lumping scheme for uncut cells on one side and the unconditional stability of implicit time integration schemes for cut cells on the other. In this hybrid, immersed Newmark IMEX approach, we use explicit second-order central differences to integrate the uncut degrees of freedom that lead to a diagonal block in the mass matrix and an implicit trapezoidal Newmark method to integrate the remaining degrees of freedom (those supported by at least one cut cell). The immersed Newmark IMEX approach preserves the high-order convergence rates and the geometric flexibility of the finite cell method while retaining the efficiency of the nodal lumping scheme for uncut cells without compromising the critical time step size. We analyze a simple system of spring-coupled masses to highlight some of the essential characteristics of Newmark IMEX time integration. We then solve the scalar wave equation with immersed boundaries on two- and three-dimensional examples with significant geometric complexity to show that our approach is more efficient than state-of-the-art time integration schemes when comparing accuracy and runtime. While we focus on the finite and spectral cell methods, we expect other immersed methods to benefit equally from Newmark IMEX time integration.

Efficient simulation of two-dimensional time-fractional Navier–Stokes equations using RBF-FD approach

This study is focused on the Navier–Stokes equations. The numerical method employed in this study for solving the Navier–Stokes equations is a meshless approach that combines the radial basis function (RBF) and finite difference (FD) schemes on two-dimensional arbitrary domains. This method, known as radial basis function-generated finite differences (RBF-FD), utilizes the finite difference technique to discretize the time domain and the RBF-FD approximation to discretize the space domain. The RBF-FD method offers several advantages, including high-order convergence rates and the ability to effectively handle sparse node layouts within arbitrary domains. The efficacy and accuracy of the suggested approach are evaluated through the analysis of four unique instances: comprising Taylor–Green vortices, lid-driven cavity flow, backward-facing step flow, and flow around a square cylinder. In the last example, we examined a simulation of air pollution around buildings.

A discontinuous Galerkin approach for atmospheric flows with implicit condensation

We present a discontinuous Galerkin method for moist atmospheric dynamics, with and without warm rain. By considering a combined density for water vapour and cloud water, we avoid the need to model and compute a source term for condensation. We recover the vapour and cloud densities by solving a pointwise non-linear problem each time step. Consequently, we enforce the requirement for the water vapour not to be supersaturated implicitly. Together with an explicit time-stepping scheme, the method is highly parallelisable and can utilise high-performance computing hardware. Furthermore, the discretisation works on structured and unstructured meshes in two and three spatial dimensions. We illustrate the performance of our approach using several test cases in two and three spatial dimensions. In the case of a smooth, exact solution, we illustrate the optimal higher-order convergence rates of the method.

Discretization Error Estimation Using the Unsteady Error Transport Equations

Abstract Computational Fluid Dynamics (CFD) has gained significant utility in the analysis of diverse fluid flow scenarios, thanks to advances in computational power. Accurate error estimation techniques are crucial to ensure the reliability of CFD simulations, as errors can lead to misleading conclusions. This study focuses on the estimation of discretization errors in time-dependent simulations, building upon prior work addressing steady-state problems Wang et al. (2020, “Error Transport Equation Implementation in the Sensei CFD Code,” AIAA Paper No. 2020-1047). In this research, we employ unsteady error transport equations (ETE) to generate localized discretization error estimates within the framework of the finite volume CFD code SENSEI. For steady-state problems, the ETE only need to be solved once after the solution has converged, whereas the unsteady ETE need to be co-advanced with the primal solve. To enhance efficiency, we adopt a one-sided temporal stencil and develop a modified iterative correction process tailored to the unsteady ETE. The time-marching schemes utilized encompass second-order accurate singly diagonally implicit Runge–Kutta (SDIRK) and second-order backward differentiation formula (BDF2), both being implicit. To rigorously assess the accuracy of our error estimates, all test cases feature known analytical solutions, facilitating order-of-accuracy evaluations. Two test cases are considered: the 2D convected vortex for inviscid flow and a cross-term sinusoidal (CTS) manufactured solution for viscous flow. Results indicate higher-order convergence rates for the 2D convected vortex test case even when iterative correction is not applied, with similar observations in the CTS case, albeit not at the finest grid levels. Although the current implementation of iterative correction exhibits lower stability compared to the primal solve, it generally enhances the discretization error estimate. Notably, after iterative correction, the discretization error estimate for the unsteady ETE achieves higher-order accuracy across all grid levels in the 2D CTS manufactured solution.

A regularized high-order diffusive smoothed particle hydrodynamics scheme without tensile instability

In the present work, we derive a novel high-order weakly compressible smoothed particle hydrodynamics scheme based on an accurate approximation of the pressure gradient and on the use of numerical Riemann fluxes. Specifically, a switch between non-conservative and conservative formulations of the pressure gradient is adopted close to the free surface, in order to fulfill the dynamic free-surface boundary condition and, at the same time, prevent the onset of the tensile instability in inner regions of the fluid domain. The numerical diffusion is obtained using Riemann solvers, with reconstruction/limitation of the left and right states derived from the Monotonic Upstream-centered Scheme for Conservation Laws technique. These allow for a high-order convergence rate of the diffusive terms that, for increasing spatial resolutions, results in a low numerical dissipation without tuning parameters. Regular particle distributions, which are crucial for the model accuracy, are obtained thanks to recent improvements in Particle Shifting Techniques. These are taken into account within the constitutive equations through a quasi-Lagrangian formalism. The energy balance of such a non-conservative formulation is derived, and an in-depth analysis of the term contributing to numerical dissipation is performed. The numerical investigation is carried out on several problems, illustrating the advantages of the present scheme with respect to conservative formulations. Since the proposed formulation does not intrinsically guarantee momenta conservation, the latter are monitored showing that the overall errors are generally small.

A new analytical approach for the local radial point interpolation discretisation in space and applications to high-order in time schemes for two-dimensional fractional PDEs

In weak-form formulations of the local radial point interpolation method (LRPIM) for the solution of partial differential equations, space discretisation matrices have most often been obtained entirely through numerical integration. This work introduces a novel approach which derives closed-form expressions for obtaining the entries of the discretisation matrix for the solution of two-dimensional time-fractional diffusion problems. This analytical approach also yields a closed-form formula for the approximation of the Laplacian. These techniques are then applied for developing LRPIM-based numerical algorithms. Since the exact solutions usually have unbounded first-order time derivatives at time zero, a graded mesh is employed for a high-order in time approximation of the Caputo derivative. The analytical shape functions are used to develop a weak-form algorithm and a strong-form algorithm is developed using the analytical approximation of the Laplacian. We demonstrate that computed solutions obtained using the weak-form and strong-form algorithms have the accuracy levels consistent with the theoretical accuracy in space. An appropriate choice of the mesh grading parameter yields a high-order convergence rate in time. The unconditional stability and convergence of a LRPIM strong form algorithm on a uniform temporal mesh is established under the assumption that the exact solution has sufficient regularity.

Stabilized isogeometric collocation methods for scalar transport and incompressible fluid flow

In this work we adapt classical residual-based stabilization techniques to the spline collocation setting. Inspired by the Streamline-Upwind-Petrov-Galerkin and Pressure-Stabilizing-Petrov-Galerkin methods, our stabilized collocation schemes address spurious oscillations that can arise from advection and pressure instabilities. Numerical examples for the advection–diffusion equation, Stokes equations, and incompressible Navier–Stokes equations show the effectiveness of the proposed stabilized schemes while maintaining the high-order convergence rates and accuracy of standard isogeometric collocation on smooth problems.

On the simultaneous solution of structural membranes on all level sets within a bulk domain

A mechanical model and numerical method for structural membranes implied by all isosurfaces of a level-set function in a three-dimensional bulk domain are proposed. The mechanical model covers large displacements in the context of the finite strain theory and is formulated based on the tangential differential calculus. Alongside curved two-dimensional membranes embedded in three dimensions, also the simpler case of curved ropes (cables) in two-dimensional bulk domains is covered. The implicit geometries (shapes) are implied by the level sets and the boundaries of the structures are given by the intersection of the level sets with the boundary of the bulk domain. For the numerical analysis, the bulk domain is discretized using a background mesh composed by (higher-order) elements with the dimensionality of the embedding space. The elements are by no means aligned to the level sets, i.e., the geometries of the structures, which resembles a fictitious domain method, most importantly the Trace FEM. The proposed numerical method is a hybrid of the classical FEM and fictitious domain methods which may be labelled as “Bulk Trace FEM”. Numerical studies confirm higher-order convergence rates and the potential for new material models with continuously embedded sub-structures in bulk domains.

A stable radial basis function partition of unity method for solving convection–diffusion equations on surfaces

The radial basis function partition of unity (RBF-PU) method is used to solve various types of partial differential equations, since it is based on scattered node sets that can be easily implemented. In this paper, a framework for the RBF-PU method is proposed to solve the convection–diffusion equations on closed surfaces. The surface gradient and Laplace–Beltrami operators are discretized based on the RBF-PU method using the polyharmonic spline radial basis functions augmented with polynomials. Since the RBF-PU method is a local meshless method, an advantage is that the differential matrix obtained by discretizing the surface differential operator is sparse. Additionally, solving convection-dominated problems or transport equations directly using RBF-PU method can produce numerical oscillations. For this reason, we add hyperviscosity term to improve the stability of the RBF-PU method for solving the convection–diffusion equations on surfaces, and numerical experiments show that such a combination has higher-order convergence rates. Finally, we solve the discontinuous source problem on the torus and the Turing system with convection term on different surfaces to demonstrate the effectiveness of the proposed method.

A stable RBF-FD method for solving two-dimensional variable-order time fractional advection-diffusion equation

The meshless method is used in this paper to provide a numerical solution for a two-dimensional variable-order time fractional advection-diffusion equation. The radial basis function and the finite difference scheme on two-dimensional (2-D) arbitrary domains are the foundations of this accurate and reliable meshless approach. The name of this approach is radial basis function-generated finite differences (RBF-FD). The finite difference technique and the RBF-FD approximation are both employed for discrete time and space domains, respectively. The stability and convergence related to the time discretization of this approach are then discussed and theoretically proven. The advantage of the RBF-FD method is that it has high-order convergence rates and its ability to manage sparse node layouts in arbitrary domains. The use of polyharmonic splines (PHS) as a RBF and the absence of the shape parameter are also a positive points for the mentioned method. The effectiveness and precision of the proposed method are examined with three distinct examples, incorporating two model examples and a real-world application of contamination transportation phenomena.

High-order weak Galerkin scheme for [formula omitted]-elliptic interface problems

A fitted weak Galerkin (WG) finite element scheme has been demonstrated in the present paper for solving H(div)-elliptic equation with discontinuous coefficients and interface. In both L2 norm and H1 norm, error estimates of optimal orders are discussed for the H(div)-elliptic interface problems. High-order convergence rates have been achieved by using suitable WG approximation spaces of higher degrees. We have performed some typical numerical tests to confirm the theoretical findings of the proposed WG algorithm. Moreover, numerically it is shown that the proposed WG algorithm is able to accommodate geometrically complicated and very irregular interfaces having sharp edges, cusps, and tips.

Interpolation-based immersed finite element and isogeometric analysis

We introduce a new paradigm for immersed finite element and isogeometric methods based on interpolating function spaces from an unfitted background mesh into Lagrange finite element spaces defined on a foreground mesh that captures the domain geometry but is otherwise subject to minimal constraints on element quality or connectivity. This is a generalization of the concept of Lagrange extraction from the isogeometric analysis literature and also related to certain variants of the finite cell and material point methods, as well as the concept of free-form deformation from design. Crucially, the interpolation may be approximate without sacrificing high-order convergence rates, which distinguishes the present method from existing high-order finite cell, CutFEM, and immersogeometric approaches. The interpolation paradigm also permits non-invasive reuse of existing finite element software for immersed analysis. We analyze the properties of the interpolation-based immersed paradigm for a model problem and implement it on top of the open-source FEniCS finite element software, to apply it to a variety of problems in fluid, solid, and structural mechanics where we demonstrate high-order accuracy and applicability to practical geometries like trimmed spline patches.

Extended reacting boundary modeling of porous materials with thin coverings for time-domain room acoustic simulations

Modeling of acoustic boundary conditions has a significant impact on the accuracy of room acoustic simulations, which play an important role in the design phase of indoor built environments in order to improve the acoustical comfort. In this work, a numerical framework based on the discontinuous Galerkin (DG) method is presented for modeling extended reacting boundaries of porous absorbers covered by thin materials. The domain decomposition methodology is applied by treating the porous material as a subdomain. Equivalent fluid models are used to depict the acoustic properties of porous materials, whose effective density and compressibility as irrational functions are approximated by multipole rational functions in the frequency domain. By employing the auxiliary differential equation approach to calculate the time convolution, the augmented time-domain governing equations of porous materials can be expressed in the same unified hyperbolic form as the linear acoustic equations, which further enables a consistent upwind numerical flux formulation throughout the whole domain. The numerical coupling across the interface between propagation media is handled by solving the underlying Riemann problem. Compared to existing approaches with the DG method for extended reacting boundaries modeling for room acoustics, the derived upwind numerical flux formulation does not involve the computation of auxiliary variables. The presented framework yield a well-posed linear hyperbolic system with admissible boundary conditions as guided by the “uniform Kreiss condition” (Kreiss, 1970). Acoustic properties of the covering materials are illustrated by considering a limp permeable membrane model. A local time-stepping approach is utilized to improve computational efficiency. Numerical validations against analytical solutions in 1D are performed to verify the desired high-order convergence rate. A 3D case study on modeling spherical wave fronts demonstrates the broadband accuracy of the formulation.

On the efficiency of 5(4) RK-embedded pairs with high order compact scheme and Robin boundary condition for options valuation

When solving the American options with or without dividends, numerical methods often obtain lower convergence rates if further treatment is not implemented even using high-order schemes. In this article, we present a fast and explicit fourth-order compact scheme for solving the free boundary options. In particular, the early exercise features with the asset option and option sensitivity are computed based on a coupled of nonlinear PDEs with fixed boundaries for which a high order analytical approximation is obtained. Furthermore, we implement a new treatment at the left boundary by introducing a third-order Robin boundary condition. Rather than computing the optimal exercise boundary from the analytical approximation, we simply obtain it from the asset option based on the linear relationship at the left boundary. As such, a high order convergence rate can be achieved. We validate by examples that the improvement at the left boundary yields a fourth-order convergence rate without further implementation of mesh refinement, Rannacher time-stepping, and/or smoothing of the initial condition. Furthermore, we extensively compare the performance of our present method with several 5(4) Runge–Kutta pairs and observe that Dormand and Prince and Bogacki and Shampine 5(4) pairs are faster and provide more accurate numerical solutions. Based on numerical results and comparison with other existing methods, we can validate that the present method is very fast and provides more accurate solutions with very coarse grids.

High-order quadrature on multi-component domains implicitly defined by multivariate polynomials

A high-order quadrature algorithm is presented for computing integrals over curved surfaces and volumes whose geometry is implicitly defined by the level sets of (one or more) multivariate polynomials. The algorithm recasts the implicitly defined geometry as the graph of an implicitly defined, multi-valued height function, and applies a dimension reduction approach needing only one-dimensional quadrature. In particular, we explore the use of Gauss-Legendre and tanh-sinh methods and demonstrate that the quadrature algorithm inherits their high-order convergence rates. Under the action of h-refinement with q fixed, the quadrature schemes yield an order of accuracy of 2q, where q is the one-dimensional node count; numerical experiments demonstrate up to 22nd order. Under the action of q-refinement with the geometry fixed, the convergence is approximately exponential, i.e., doubling q approximately doubles the number of accurate digits of the computed integral. Complex geometry is automatically handled by the algorithm, including, e.g., multi-component domains, tunnels, and junctions arising from multiple polynomial level sets, as well as self-intersections, cusps, and other kinds of singularities. A variety of numerical experiments demonstrates the quadrature algorithm on two- and three-dimensional problems, including: randomly generated geometry involving multiple high-curvature pieces; challenging examples involving high degree singularities such as cusps; adaptation to simplex constraint cells in addition to hyperrectangular constraint cells; and boolean operations to compute integrals on overlapping domains.

A higher-order finite element method with unstructured anisotropic mesh adaption for two phase flows with surface tension

A novel finite element framework is proposed for the numerical simulation of two phase flows with surface tension. The Level-Set (LS) method with piece-wise quadratic (P2) interpolation for the liquid–gas interface is used in order to reach higher-order convergence rates in regions with smooth interface. A balanced-force implementation of the continuum surface force model is used to take into account the surface tension and to solve static problems as accurately as possible. Given that this requires a balance between the discretization used for the LS function, and that used for the pressure field, an equal-order P2/P2/P2 scheme is proposed for the Navier–Stokes and LS advection equations, which are strongly coupled with each other. This fully implicit formulation is stabilized using the residual-based variational multiscale framework. In order to improve the accuracy and obtain optimal convergence rates with a minimum number of elements, an anisotropic mesh adaption method is proposed where the unstructured mesh is kept as fine as possible close to the zero iso-value of the P2 LS function. Elements are automatically stretched in regions with flat interface in order to keep the complexity fixed during the simulation. The accuracy and efficiency of this approach are demonstrated for two and three dimensional simulations of a rising bubble.

An efficient high-order meshless method for advection-diffusion equations on time-varying irregular domains

We present a high-order radial basis function finite difference (RBF-FD) framework for the solution of advection-diffusion equations on time-varying domains. Our framework is based on a generalization of the recently developed Overlapped RBF-FD method that utilizes a novel automatic procedure for computing RBF-FD weights on stencils in variable-sized regions around stencil centers. This procedure eliminates the overlap parameter δ, thereby enabling tuning-free assembly of RBF-FD differentiation matrices on moving domains. In addition, our framework utilizes a simple and efficient procedure for updating differentiation matrices on moving domains tiled by node sets of time-varying cardinality. Finally, advection-diffusion in time-varying domains is handled through a combination of rapid node set modification, a new high-order semi-Lagrangian method that utilizes the new tuning-free overlapped RBF-FD method, and a high-order time-integration method. The resulting framework has no tuning parameters and has O(Nlog⁡N) time complexity. We demonstrate high-orders of convergence for advection-diffusion equations on time-varying 2D and 3D domains for both small and large Peclet numbers. We also present timings that verify our complexity estimates. Finally, we utilize our method to solve a coupled 3D problem motivated by models of platelet aggregation and coagulation, once again demonstrating high-order convergence rates on a moving domain.

MDFEM: Multivariate decomposition finite element method for elliptic PDEs with lognormal diffusion coefficients using higher-order QMC and FEM

We introduce the multivariate decomposition finite element method (MDFEM) for elliptic PDEs with lognormal diffusion coefficients, that is, when the diffusion coefficient has the form a = exp(Z) where Z is a Gaussian random field defined by an infinite series expansion Z(y) = ∑j≥1yjϕj with yj ~ N(0,1) and a given sequence of functions {ϕj}j≥1. We use the MDFEM to approximate the expected value of a linear functional of the solution of the PDE which is an infinite-dimensional integral over the parameter space. The proposed algorithm uses the multivariate decomposition method (MDM) to compute the infinite-dimensional integral by a decomposition into finite-dimensional integrals, which we resolve using quasi-Monte Carlo (QMC) methods, and for which we use the finite element method (FEM) to solve different instances of the PDE. We develop higher-order quasi-Monte Carlo rules for integration over the finite-dimensional Euclidean space with respect to the Gaussian distribution by use of a truncation strategy. By linear transformations of interlaced polynomial lattice rules from the unit cube to a multivariate box of the Euclidean space we achieve higher-order convergence rates for functions belonging to a class of anchored Gaussian Sobolev spaces while taking into account the truncation error. These cubature rules are then used in the MDFEM algorithm. Under appropriate conditions, the MDFEM achieves higher-order convergence rates in terms of error versus cost, i.e., to achieve an accuracy of O(ε) the computational cost is O(ε-1/λ-d′/λ)=O(ε-(p*+d′/τ)/(1-p*)) where ε−1/λ and ε−d′/λ) are respectively the cost of the quasi-Monte Carlo cubature and the finite element approximations, with d′= d (1+δ′) for some δ′ ≥ 0 and d the physical dimension, and 0<p*≤(2+d′/τ)-1 is a parameter representing the sparsity of {ϕj}j≥1.

On a generalized energy conservation/dissipation time finite element method for Hamiltonian mechanics

Energy conserving and dissipative algorithm designs in Hamilton’s canonical equations via Petrov–Galerkin time finite element methodology are proposed in this paper that provide new avenues with high-order convergence rate, improved solution accuracy, and controllable numerical dissipation. Lagrange quadratic shape function in time with flexible interpolation points are considered to approximate the solution over a time interval. Instead of specifying the weight functions, two algorithmic parameters, namely, a principal root (ρq∞) and a spurious root (ρp∞), are introduced to formulate a generalized weight function, which enables us to introduce controllable numerical dissipation with respect to displacement and momenta while preserving high-order convergence, and features with improved solution accuracy. A family of third-order accurate time finite element algorithms with controllable dissipation and improved solution accuracy is presented in both the homogeneous and non-homogeneous dynamic problems; and via setting (ρq∞,ρp∞)=(1,1), this third-order family of algorithms directly leads to a new family of fourth-order accurate non-dissipative algorithms in general homogeneous problems; and the fourth-order accuracy is also preserved in non-homogeneous problems when the third-order time derivative of the external excitation has the order of O(Δtn)(n≤1). Numerical examples are performed to demonstrate the pros/cons for the conserving properties of various schemes in the proposed Petrov–Galerkin time finite element of algorithms.

A unified finite strain theory for membranes and ropes

The finite strain theory is reformulated in the frame of the Tangential Differential Calculus (TDC) resulting in a unification in a threefold sense. Firstly, ropes, membranes and three-dimensional continua are treated with one set of governing equations. Secondly, the reformulated boundary value problem applies to parametrized and implicit geometries. Therefore, the formulation is more general than classical ones as it does not rely on parametrizations implying curvilinear coordinate systems and the concept of co- and contravariant base vectors. This leads to the third unification: TDC-based models are applicable to two fundamentally different numerical approaches. On the one hand, one may use the classical Surface FEM where the geometry is discretized by curved one-dimensional elements for ropes and two-dimensional surface elements for membranes. On the other hand, it also applies to recent Trace FEM approaches where the geometry is immersed in a higher-dimensional background mesh. Then, the shape functions of the background mesh are evaluated on the trace of the immersed geometry and used for the approximation. As such, the Trace FEM is a fictitious domain method for partial differential equations on manifolds. The numerical results show that the proposed finite strain theory yields higher-order convergence rates independent of the numerical methodology, the dimension of the manifold, and the geometric representation type.