Comput Phys Research Articles

Accurate discretization of diffusion in the LS-STAG cut-cell method using diamond cell techniques

In cut-cell methods for fluid flow in complex geometries such as the LS-STAG method (Nikfarjam et al. Comput Phys Commun 226:67-80, 2018), mesh non-orthogonality near the immersed boundary makes inaccurate 2-point formulas for computing the face-normal gradients in diffusive fluxes. A way to improve the accuracy is to compute the whole solution gradient at the cut-cell faces, thus decomposing the flux as an orthogonal contribution (using a standard 2-point formula) and non-orthogonal correction (using data at cell vertices). The solution at cut-cell vertices is then interpolated from cell-centered data and boundary conditions. This gradient reconstruction technique is commonly denominated “diamond cell method”.The diamond cell method is implemented in the LS-STAG code for cell-centered data (temperature equation in 2D geometries) and face-centered data (Navier-Stokes equations in 3D-extruded geometries) using various interpolation schemes for vertex reconstruction. The accuracy of the diamond cell discretization is firmly assessed on a series of 2D benchmark problems (Taylor-Couette solution, natural convection from a cylinder in an enclosure) by inspecting the formal order of accuracy and the heat flux distribution at the immersed boundary. Finally, the diamond cell technique is employed for enhancing the accuracy of the velocity gradients in 3D-extruded geometries which were previously computed with 2-point formulas. These computations illustrate the versatility of the diamond cell technique, which can be applied to other cut-cell methods.

Shift boundary material point method: an image-to-simulation workflow for solids of complex geometries undergoing large deformation

We introduce a mathematical framework designed to enable a simple image-to-simulation workflow for solids of complex geometries in the geometrically nonlinear regime. While the material point method is used to circumvent the mesh distortion issues commonly exhibited in Lagrangian meshes, a shifted domain technique originated from Main and Scovazzi (J Comput Phys 372:972–995, 2018) is used to represent the boundary conditions implicitly via a level set or signed distance function. Consequently, this method completely bypasses the need to generate high-quality conformal mesh to represent complex geometries and therefore allows modelers to select the space of the interpolation function without the constraints due to the geometric need. This important simplification enables us to simulate deformation of complex geometries inferred from voxel images. Verification examples on deformable body subjected to finite rotation have shown that the new shifted domain material point method is able to generate frame-indifferent results. Meanwhile, simulations using micro-CT images of a Hostun sand have demonstrated that this method is able to reproduce the quasi-brittle damage mechanisms of single grain without the excessively concentrated nodes commonly displayed in conformal meshes that represent 3D objects with local fine details.

Runtime optimization of a memory efficient CG solver for FFT-based homogenization: implementation details and scaling results for linear elasticity

The memory efficient CG algorithm of Kabel et al. (Comput Mech 54(6):1497–1514, 2014) reduces the memory requirements of a strain based implementation of the CG algorithm following Zeman et al. (J Comput Phys 229(21):8065–8071, 2010) for solving the equations of linear elasticity by 40%. But since the Fourier wave vectors have to be recalculated at several steps of the memory efficient algorithm, the runtime increases for a straightforward implementation.. We explain how to reduce the runtime overhead to a negligible size, and show that the memory efficient algorithm scales better than the standard algorithm with up to 256 MPI processes.

High-Order Low-Dissipation Targeted ENO Schemes for Ideal Magnetohydrodynamics

The recently proposed targeted ENO (TENO) schemes (Fu et al. J Comput Phys 305:333–359, 2016) are demonstrated to feature the controllable low numerical dissipation and sharp shock-capturing property in compressible gas dynamic simulations. However, the application of low-dissipation TENO schemes to ideal magnetohydrodynamics (MHD) is not straightforward. The complex interaction between fluid mechanics and electromagnetism induces extra numerical challenges, including simultaneously preserving the ENO-property, maintaining good numerical robustness and low dissipation as well as controlling divergence errors. In this paper, based on an unstaggered constrained transport framework to control the divergence error, we extend a set of high-order low-dissipation TENO schemes ranging from 5-point to 8-point stencils to solving the ideal MHD equations. A unique set of built-in parameters for each TENO scheme is determined. Different from the TENO schemes in Fu et al. (2016), a modified scale-separation formula is developed. The new formula can achieve stronger scale separation, and it is simpler and more efficient than the previous version as the computation cost of high-order global smoothness measure $${\tau _K}$$ is avoided. The performances of tailored schemes are systematically studied by several benchmark simulations. Numerical experiments demonstrate that the TENO schemes in the constrained transport framework are promising to simulate more complex MHD flows.

Fast and robust flow simulations in discrete fracture networks with GPGPUs

A new approach for flow simulation in very complex discrete fracture networks based on PDE-constrained optimization has been recently proposed in Berrone et al. (SIAM J Sci Comput 35(2):B487–B510, 2013b; J Comput Phys 256:838–853, 2014) with the aim of improving robustness with respect to geometrical complexities. This is an essential issue, in particular for applications requiring simulations on geometries automatically generated like the ones used for uncertainty quantification analyses and hydro-mechanical simulations. In this paper, implementation of this approach in order to exploit Nvidia Compute Unified Device Architecture is discussed with the main focus to speed up the linear algebra operations required by the approach, being this task the most computational demanding part of the approach. Furthermore, two different approaches for linear algebra operations and two storage formats for sparse matrices are compared in terms of computational efficiency and memory constraints.

A multiscale approach to hybrid RANS/LES wall modeling within a high‐order discontinuous Galerkin scheme using function enrichment

SummaryWe present a novel approach to hybrid Reynolds‐averaged Navier‐Stokes (RANS)/ large eddy simulation (LES) wall modeling based on function enrichment, which overcomes the common problem of the RANS‐LES transition and enables coarse meshes near the boundary. While the concept of function enrichment as an efficient discretization technique for turbulent boundary layers has been proposed in an earlier article by Krank & Wall (A new approach to wall modeling in LES of incompressible flow via function enrichment. J Comput Phys. 2016;316:94‐116), the contribution of this work is a rigorous derivation of a new multiscale turbulence modeling approach and a corresponding discontinuous Galerkin discretization scheme. In the near‐wall area, the Navier‐Stokes equations are explicitly solved for an LES and a RANS component in one single equation. This is done by providing the Galerkin method with an independent set of shape functions for each of these two methods; the standard high‐order polynomial basis resolves turbulent eddies, where the mesh is sufficiently fine and the enrichment automatically computes the ensemble‐averaged flow if the LES mesh is too coarse. As a result of the derivation, the RANS model is applied solely to the RANS degrees of freedom, which effectively prevents the typical issue of a log‐layer mismatch in attached boundary layers. As the full Navier‐Stokes equations are solved in the boundary layer, spatial refinement gradually yields wall‐resolved LES with exact boundary conditions. Numerical tests show the outstanding characteristics of the wall model regarding grid independence, superiority compared to equilibrium wall models in separated flows, and achieve a speed‐up by two orders of magnitude compared to wall‐resolved LES.

Non-satisfiability of a positivity condition for commutator-free exponential integrators of order higher than four

We consider commutator-free exponential integrators as put forward in Alverman and Fehske (J Comput Phys 230:5930–5956, 2011). For parabolic problems, it is important for the well-definedness that such an integrator satisfies a positivity condition such that essentially it only proceeds forward in time. We prove that this requirement implies maximal convergence order of four for real coefficients, which has been conjectured earlier by other authors.

Numerical simulations of emulsions in shear flows

We present a modification of a recently developed volume of fluid method for multiphase problems (Ii et al. in J Comput Phys 231(5):2328–2358, 2012), so that it can be used in conjunction with a fractional-step method and fast Poisson solver, and validate it with standard benchmark problems. We then consider emulsions of two-fluid systems and study their rheology in a plane Couette flow in the limit of vanishing inertia. We examine the dependency of the effective viscosity mu on the volume fraction varPhi (from 10 to 30%) and the Capillary number Ca (from 0.1 to 0.4) for the case of density and viscosity ratio 1. We show that the effective viscosity decreases with the deformation and the applied shear (shear-thinning) while exhibiting a non-monotonic behavior with respect to the volume fraction. We report the appearance of a maximum in the effective viscosity curve and compare the results with those of suspensions of rigid and deformable particles and capsules. We show that the flow in the solvent is mostly a shear flow, while it is mostly rotational in the suspended phase; moreover, this behavior tends to reverse as the volume fraction increases. Finally, we evaluate the contributions to the total shear stress of the viscous stresses in the two fluids and of the interfacial force between them.

A Fourth-Order Kernel-Free Boundary Integral Method for the Modified Helmholtz Equation

Based on the kernel-free boundary integral method proposed by Ying and Henriquez (J Comput Phys 227(2):1046---1074, 2007), which is a second-order accurate method for general elliptic partial differential equations, this work develops it to be a fourth-order accurate version for the modified Helmholtz equation. The updated method is in line with the original one. Unlike the traditional boundary integral method, it does not need to know any analytical expression of the fundamental solution or Green's function in evaluation of boundary or volume integrals. Boundary value problems under consideration are reformulated into Fredholm boundary integral equations of the second kind, whose corresponding discrete forms are solved with the simplest Krylov subspace iterative method, the Richardson iteration. During each iteration, a Cartesian grid based nine-point compact difference scheme is used to discretize the simple interface problem whose solution is the boundary or volume integral in the BIEs. The resulting linear system is solved by a fast Fourier transform based solver, whose computational work is roughly proportional to the number of grid nodes in the Cartesian grid used. As the discrete boundary integral equations are well-conditioned, the iteration converges within an essentially fixed number of steps, independent of the mesh parameter. Numerical results are presented to verify the solution accuracy and demonstrate the algorithm efficiency.

An Alternative Formulation of Discontinous Galerkin Schemes for Solving Hamilton–Jacobi Equations

The aim of this paper is to develop an alternative formulation of discontinuous Galerkin (DG) schemes for approximating the viscosity solutions to nonlinear Hamilton–Jacobi (HJ) equations. The main difficulty in designing DG schemes lies in the inherent non-divergence form of HJ equations. One effective approach is to explore the elegant relationship between HJ equations and hyperbolic conservation laws: the standard DG scheme is applied to solve a conservation law system satisfied by the derivatives of the solution of the HJ equation. In this paper, we consider an alternative approach to directly solving the HJ equations, motivated by a class of successful direct DG schemes by Cheng and Shu (J Comput Phys 223(1):398–415, 2007), Cheng and Wang (J Comput Phys 268:134–153, 2014). The proposed scheme is derived based on the idea from the central-upwind scheme by Kurganov et al. (SIAM J Sci Comput 23(3):707–740, 2001). In particular, we make use of precise information of the local speeds of propagation at discontinuous element interface with the goal of adding adequate numerical viscosity and hence naturally capturing the viscosity solutions. A collection of numerical experiments is presented to demonstrate the performance of the method for solving general HJ equations with linear, nonlinear, smooth, non-smooth, convex, or non-convex Hamiltonians.

Hybridized SLAU2–HLLI and hybridized AUSMPW$$+$$–HLLI Riemann solvers for accurate, robust, and efficient magnetohydrodynamics (MHD) simulations, part I: one-dimensional MHD

SLAU2 and AUSMPW $$+$$ , both categorized as AUSM-type Riemann solvers, have been extensively developed in gasdynamics. They are based on a splitting of the numerical flux into advected and pressure parts. In this paper, these two Riemann solvers have been extended to magnetohydrodynamics (MHD). The SLAU2 Riemann solver has the favorable attribute that its dissipation for low-speed flows scales as $$O(M^{2})$$ , where M is the Mach number. This is the physical scaling required for low-speed flows, and the dissipation in SLAU2 for MHD is engineered to have this low Mach number scaling. The AUSMPW $$+$$ , when its pressure flux is replaced with that of SLAU2, has the same low Mach number scaling. At higher Mach numbers, however, the pressure-split Riemann solvers were found not to function well for some MHD Riemann problems, despite the fact that they were engineered to have a dissipation that scales as $$O(\vert M\vert )$$ for high Mach number flows. The HLLI Riemann solver (Dumbser and Balsara in J Comput Phys 304:275–319, 2016) has a dissipation that scales as $$O(\vert M\vert )$$ , which makes it unsuitable for low Mach number flows. However, it has very favorable performance for higher Mach number MHD flows. Since the two families of Riemann solvers perform very well over a range of intermediate Mach numbers, the best way to benefit from the mutually complementary strengths of both these Riemann solvers is to hybridize between them. The result is an all-speed Riemann solver for MHD. We, therefore, document hybridized SLAU2–HLLI and AUSMPW $$+$$ –HLLI Riemann solvers. The hybrid Riemann solvers suppress the oscillations that appeared in single-solver solutions, and they also preserve contact discontinuities, as well as Alfven waves, very well. Furthermore, their better resolution at low speeds has been demonstrated. We also present several stringent one-dimensional test problems.

A Class of Low Dissipative Schemes for Solving Kinetic Equations

We introduce an extension of the fast semi-Lagrangian scheme developed in J Comput Phys 255:680–698 (2013) in an effort to increase the spatial accuracy of the method. The basic idea of this extension is to modify the first-order accurate transport step of the original semi-Lagrangian scheme to allow for a general piecewise polynomial reconstruction of the distribution function. For each discrete velocity, we update the solution not at cell centers, but rather at the extreme points of the spatial reconstruction, the locations of which are different for each discrete velocity and change with time. Several approaches are discussed for evaluating the collision operator at these extreme points using only cell center values by making special assumption on the spatial variation of the collision operator. The result is a class of schemes that preserves the structure of the solution over very long times when compared to existing schemes in the literature. As a proof of concept, the new method is implemented in a one-dimensional setting, using piecewise linear reconstructions of the distribution function together with a related reconstruction of the collision operator. The method is derived both for the relatively simple Bhatnagar–Gross–Krook (BGK) operator as well as for the classical Boltzmann operator. Several numerical tests are used to assess the performance of the implementation, including comparisons with the original method in J Comput Phys 255:680–698 (2013) and with classical semi-Lagrangian methods of first and second order. In convergence tests, we observe uniform second-order accuracy across all flow regimes for the BGK operator and nearly second-order accuracy for the Boltzmann operator. In addition, we observe that the method outperforms the classical semi-Lagrangian approach, in particular when resolving fine solution structures in space.

Numerical study on the convergence to steady-state solutions of a new class of finite volume WENO schemes: triangular meshes

In this paper, we continue our research on the numerical study of convergence to steady-state solutions for a new class of finite volume weighted essentially non-oscillatory (WENO) schemes in Zhu and Shu (J Comput Phys 349:80–96, 2017), from tensor product meshes to triangular meshes. For the case of triangular meshes, this new class of finite volume WENO schemes was designed for time-dependent conservation laws in Zhu and Qiu (SIAM J Sci Comput 40(2):A903–A928, 2018) for the third- and fourth-order versions. In this paper, we extend the design to a new fifth-order version in the same framework to keep the essentially non-oscillatory property near discontinuities. Similar to the case of tensor product meshes in Zhu and Shu (2017), by performing such spatial reconstruction procedures together with a TVD Runge–Kutta time discretization, these WENO schemes do not suffer from slight post-shock oscillations that are responsible for the phenomenon wherein the residues of classical WENO schemes hang at a truncation error level instead of converging to machine zero. The third-, fourth-, and fifth-order finite volume WENO schemes in this paper can suppress the slight post-shock oscillations and have their residues settling down to a tiny number close to machine zero in steady-state simulations in our extensive numerical experiments.

A New Robust Carbuncle-Free Roe Scheme for Strong Shock

The paper devises a new robust carbuncle-free Roe Riemann solver for strong shock, different from hybrid method, entropy fix, Liou’s conjecture (J Comput Phys 160:623–648, 2000) and artificial viscosity (Rodionov in J Comput Phys 345:308–329, 2017). Roe scheme encounters carbuncle phenomenon, violates entropy condition and lacks positivity property. The remedy integrates shear viscosity into momentum flux to damp undesirable perturbation and prevents unstable vorticity mode from triggering shock instability. Compared to HLL scheme, shear viscosity is properly established through dimensional analysis and analogy method. The non-linear wave speeds are slightly modified by incorporating the neighboring cell information for positive conservation. The pressure-based sensing function is applied to preserve shear layer while keeping shock robustness. The resulting scheme is very easily implemented by converting the existing Roe code. The matrix stability analysis confirms that this approach is shock-stable and more robust than entropy fix. A series of numerical results demonstrate its potential features: positivity-preserving property, entropy-satisfying property, accurate boundary-layer resolution and high robustness against shock instability. In addition, it signifies that the reason to cause shock instability for Roe scheme may be not the pressure difference term in mass flux, but the inadequate shear viscosity.

Circulant embedding with QMC: analysis for elliptic PDE with lognormal coefficients

In a previous paper (Graham et al. in J Comput Phys 230:3668–3694, 2011), the authors proposed a new practical method for computing expected values of functionals of solutions for certain classes of elliptic partial differential equations with random coefficients. This method was based on combining quasi-Monte Carlo (QMC) methods for computing the expected values with circulant embedding methods for sampling the random field on a regular grid. It was found capable of handling fluid flow problems in random heterogeneous media with high stochastic dimension, but no convergence theory was provided. This paper provides a convergence analysis for the method in the case when the QMC method is a specially designed randomly shifted lattice rule. The convergence result depends on the eigenvalues of the underlying nested block circulant matrix and can be independent of the number of stochastic variables under certain assumptions. In fact the QMC analysis applies to general factorisations of the covariance matrix to sample the random field. The error analysis for the underlying fully discrete finite element method allows for locally refined meshes (via interpolation from a regular sampling grid of the random field). Numerical results on a non-regular domain with corner singularities in two spatial dimensions and on a regular domain in three spatial dimensions are included.

Simulation of impacts on elastic\u2013viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes

The flux-difference splitting finite volume method (Leveque in J Comput Phys 131:327–353, 1997; Leveque in Finite volume methods for hyperbolic problems. Cambridge: Cambridge University Press, 2002) is here employed to perform numerical simulation of impacts on elastic–viscoplastic solids on bidimensional non-uniform quadrilateral meshes. The formulation is second order accurate in space through flux limiters, embeds the corner transport upwind method, and uses a fractional-step method to compute the relaxation operator. Elastic–viscoplastic constitutive models falling within the framework of generalized standard materials (Halphen and Nguyen in J Mech 14:667–688, 1975) in small strains are considered. Many test cases are proposed and two particular viscoplastic constitutive models are studied, on which comparisons with finite element solutions show a very good accuracy of the finite volume solutions, both on stresses and viscoplastic strains.

High Spatial Order Energy Stable FDTD Methods for Maxwell’s Equations in Nonlinear Optical Media in One Dimension

In this paper, we consider electromagnetic (EM) wave propagation in nonlinear optical media in one spatial dimension. We model the EM wave propagation by the time-dependent Maxwell’s equations coupled with a system of nonlinear ordinary differential equations (ODEs) for the response of the medium to the EM waves. The nonlinearity in the ODEs describes the instantaneous electronic Kerr response and the residual Raman molecular vibrational response. The ODEs also include the single resonance linear Lorentz dispersion. For such model, we will design and analyze fully discrete finite difference time domain (FDTD) methods that have arbitrary (even) order in space and second order in time. It is challenging to achieve provable stability for fully discrete methods, and this depends on the choices of temporal discretizations of the nonlinear terms. In Bokil et al. (J Comput Phys 350:420–452, 2017), we proposed novel modifications of second-order leap-frog and trapezoidal temporal schemes in the context of discontinuous Galerkin methods to discretize the nonlinear terms in this Maxwell model. Here, we continue this work by developing similar time discretizations within the framework of FDTD methods. More specifically, we design fully discrete modified leap-frog FDTD methods which are proved to be stable under appropriate CFL conditions. These method can be viewed as an extension of the Yee-FDTD scheme to this nonlinear Maxwell model. We also design fully discrete trapezoidal FDTD methods which are proved to be unconditionally stable. The performance of the fully discrete FDTD methods are demonstrated through numerical experiments involving kink, antikink waves and third harmonic generation in soliton propagation.

Conservative High Order Positivity-Preserving Discontinuous Galerkin Methods for Linear Hyperbolic and Radiative Transfer Equations

We further investigate the high order positivity-preserving discontinuous Galerkin (DG) methods for linear hyperbolic and radiative transfer equations developed in Yuan et al. (SIAM J Sci Comput 38:A2987–A3019, 2016). The DG methods in Yuan et al. (2016) can maintain positivity and high order accuracy, but they rely both on the scaling limiter in Zhang and Shu (J Comput Phys 229:8918–8934, 2010) and a rotational limiter, the latter may alter cell averages of the unmodulated DG scheme, thereby affecting conservation. Even though a Lax–Wendroff type theorem is proved in Yuan et al. (2016), guaranteeing convergence to weak solutions with correct shock speed when such rotational limiter is applied, it would still be desirable if a conservative DG method without changing the cell averages can be obtained which has both high order accuracy and positivity-preserving capability. In this paper, we develop and analyze such a DG method for both linear hyperbolic equations and radiative transfer equations. In the one-dimensional case, the method uses traditional DG space $$P^k$$ of piecewise polynomials of degree at most k. A key result is proved that the unmodulated DG solver in this case can maintain positivity of the cell average if the inflow boundary value and the source term are both positive, therefore the positivity-preserving framework in Zhang and Shu (2010) can be used to obtain a high order conservative positivity-preserving DG scheme. Unfortunately, in two-dimensions this is no longer the case. We show that the unmodulated DG solver based either on $$P^k$$ or $$Q^k$$ spaces (piecewise kth degree polynomials or piecewise tensor-product kth degree polynomials) could generate negative cell averages. We augment the DG space with additional functions so that the positivity of cell averages from the unmodulated DG solver can be restored, thereby leading to high order conservative positivity-preserving DG scheme based on these augmented DG spaces following the framework in Zhang and Shu (2010). Computational results are provided to demonstrate the good performance of our DG schemes.

The BR1 Scheme is Stable for the Compressible Navier–Stokes Equations

In this work we prove that the original (Bassi and Rebay in J Comput Phys 131:267–279, 1997) scheme (BR1) for the discretization of second order viscous terms within the discontinuous Galerkin collocation spectral element method (DGSEM) with Gauss Lobatto nodes is stable. More precisely, we prove in the first part that the BR1 scheme preserves energy stability of the skew-symmetric advection term DGSEM discretization for the linearized compressible Navier–Stokes equations (NSE). In the second part, we prove that the BR1 scheme preserves the entropy stability of the recently developed entropy stable compressible Euler DGSEM discretization of Carpenter et al. (SIAM J Sci Comput 36:B835–B867, 2014) for the non-linear compressible NSE, provided that the auxiliary gradient equations use the entropy variables. Both parts are presented for fully three-dimensional, unstructured curvilinear hexahedral grids. Although the focus of this work is on the BR1 scheme, we show that the proof naturally includes the Local DG scheme of Cockburn and Shu.

WSGD-OSC Scheme for Two-Dimensional Distributed Order Fractional Reaction–Diffusion Equation

In this paper, a new numerical approximation is discussed for the two-dimensional distributed-order time fractional reaction–diffusion equation. Combining with the idea of weighted and shifted Grunwald difference (WSGD) approximation (Tian et al. in Math Comput 84:1703–1727, 2015; Wang and Vong in J Comput Phys 277:1–15, 2014) in time, we establish orthogonal spline collocation (OSC) method in space. A detailed analysis shows that the proposed scheme is unconditionally stable and convergent with the convergence order $$\mathscr {O}(\tau ^2+\Delta \alpha ^2+h^{r+1})$$ , where $$\tau , \Delta \alpha , h$$ and r are, respectively the time step size, step size in distributed-order variable, space step size, and polynomial degree of space. Interestingly, we prove that the proposed WSGD-OSC scheme converges with the second-order in time, where OSC schemes proposed previously (Fairweather et al. in J Sci Comput 65:1217–1239, 2015; Yang et al. in J Comput Phys 256:824–837, 2014) can at most achieve temporal accuracy of order which depends on the order of fractional derivatives in the equations and is usually less than two. Some numerical results are also given to confirm our theoretical prediction.