Implicit Time Stepping Research Articles

Guidance on the Selection of Efficient Computational Methods for Multimedia Fate Models

Dynamic multimedia fate models (MFMs) have to deal with the temporal and spatial variation of physical-chemical properties, environmental scenarios, and chemical emissions. In such complex simulation tools, an analytical solution is not practically feasible, and even a numerical approach requires a suitable choice of the method in order to obtain satisfying speed and reliability, particularly when certain combinations of modeling scenarios and chemical properties occur. In this paper, considering some examples of a wide range of realistic chemical and scenario properties, some sources of stiffness in MFM equations are pinpointed. Next, a comparison of the performances of several numerical schemes (chosen as representatives of three wide classes) is performed. The accuracy and the computational effort required by each method is evaluated, illustrating the general effectiveness of automatically adapted timesteps in numerical algorithms and the pros and cons of implicit timestepping. The results show that automatic error control methods can significantly improve the quality of the computed solutions and most often lead to relevant savings in computing time. Additionally, explicit and implicit methods are compared, indicating that an implicit method of medium order (around 5) is the best choice as a general purpose MFM computing engine.

Dusty gas with smoothed particle hydrodynamics - II. Implicit timestepping and astrophysical drag regimes

In a companion paper (Laibe & Price 2011b), we have presented an algorithm for simulating two-fluid gas and dust mixtures in Smoothed Particle Hydrodynamics (SPH). In this paper, we develop an implicit timestepping method that preserves the exact conservation of the both linear and angular momentum in the underlying SPH algorithm, but unlike previous schemes, allows the iterations to converge to arbitrary accuracy and is suited to the treatment of non- linear drag regimes. The algorithm presented in Paper I is also extended to deal with realistic astrophysical drag regimes, including both linear and non-linear Epstein and Stokes drag. The scheme is benchmarked against the test suite presented in Paper I, including i) the analytic solutions of the dustybox problem and ii) solutions of the dustywave, dustyshock, dustysedov and dustydisc obtained with explicit timestepping. We find that the implicit method is 1- 10 times faster than the explicit temporal integration when the ratio r between the the timestep and the drag stopping time is 1 < r < 1000.

Dusty gas with smoothed particle hydrodynamics - I. Algorithm and test suite

We present a new algorithm for simulating two-fluid gas and dust mixtures in Smoothed Particle Hydrodynamics (SPH), systematically addressing a number of key issues including the generalised SPH density estimate in multi-fluid systems, the consistent treatment of variable smoothing length terms, finite particle size, time step stability, thermal coupling terms and the choice of kernel and smoothing length used in the drag operator. We find that using double-hump shaped kernels improves the accuracy of the drag interpolation by a factor of several hundred compared to the use of standard SPH bell-shaped kernels, at no additional computational expense. In order to benchmark our algorithm, we have developed a comprehensive suite of standardised, simple test problems for gas and dust mixtures: dustybox, dustywave, dustyshock, dustysedov and dustydisc, the first three of which have known analytic solutions. We present the validation of our algorithm against all of these tests. In doing so, we show that the spatial resolution criterion \Delta < cs ts is a necessary condition in all gas+dust codes that becomes critical at high drag (i.e. small stopping time ts) in order to correctly predict the dynamics. Implicit timestepping and the implementation of realistic astrophysical drag regimes are addressed in a companion paper.

Femoral bone strains during antegrade nailing: A comparison of two entry points with identical nails using finite element analysis

BackgroundAntegrade femoral nailing has become the standard treatment for diaphyseal femoral shaft fractures. Concerns linger that improper location of the nail entry point may lead to iatrogenic fracture and further complications. This study used finite element analysis to compare the strain magnitude and distribution resulting from each of two entry points in the proximal femur during antegrade nailing. MethodsA finite element model was created from a CT scan of a 37year old male femur and of a standard antegrade nail. Using implicit time-stepping, the nail was inserted through piriformis and trochanteric entry points and strain was computed at 9 anatomic locations. FindingsThe strain levels were higher overall when inserting a nail through the trochanteric starting point. The highest strain occurred immediately medial and lateral to the trochanteric entry point. The posterior greater trochanter also showed very high strain levels during nail insertion. All strain values for nail insertion through the piriformis entry point were less than 2000μm/m. InterpretationThe trochanteric entry will have a much greater potential of iatrogenic fracture of the proximal femur during insertion of a nail. Strains with this entry point exceed the yield level of bone and the repeated loading with the progression of the nail could cause fissures or fractures. Caution should be taken during insertion of an antegrade nail when utilizing a lateral trochanteric starting point secondary to an increased risk of trochanteric fracture and lateral cortex fracture.

Comparison of local weak and strong form meshless methods for 2-D diffusion equation

A comparison between weak form meshless local Petrov–Galerkin method (MLPG) and strong form meshless diffuse approximate method (DAM) is performed for the diffusion equation in two dimensions. The shape functions are in both methods obtained by moving least squares (MLS) approximation with the polynomial weight function of the fourth order on the local support domain with 13 closest nodes. The weak form test functions are similar to the MLS weight functions but defined over the square quadrature domain. Implicit timestepping is used. The methods are tested in terms of average and maximum error norms on uniform and non-uniform node arrangements on a square without and with a hole for a Dirichlet jump problem and involvement of Dirichlet and Neumann boundary conditions. The results are compared also to the results of the finite difference and finite element method. It has been found that both meshless methods provide a similar accuracy and the same convergence rate. The advantage of DAM is in simpler numerical implementation and lower computational cost.

Implicit-explicit evolution of single black holes

Numerical simulations of binary black holes---an important predictive tool for the detection of gravitational waves---are computationally expensive, especially for binaries with high mass ratios or with rapidly spinning constituent holes. Existing codes for evolving binary black holes rely on explicit timestepping methods, for which the timestep size is limited by the smallest spatial scale through the Courant-Friedrichs-Lewy condition. Binary inspiral typically involves spatial scales (the spatial resolution required by a small or rapidly spinning hole) which are orders of magnitude smaller than the relevant (orbital, precession, and radiation-reaction) timescales characterizing the inspiral. Therefore, in explicit evolutions of binary black holes, the timestep size is typically orders of magnitude smaller than the relevant physical timescales. Implicit timestepping methods allow for larger timesteps, and they often reduce the total computational cost (without significant loss of accuracy) for problems dominated by spatial rather than temporal error, such as for binary-black-hole inspiral in corotating coordinates. However, fully implicit methods can be difficult to implement for nonlinear evolution systems like the Einstein equations. Therefore, in this paper we explore implicit-explicit (IMEX) methods and use them for the first time to evolve black-hole spacetimes. Specifically, as a first step toward IMEX evolution of a full binary-black-hole spacetime, we develop an IMEX algorithm for the generalized harmonic formulation of the Einstein equations and use this algorithm to evolve stationary and perturbed single-black-hole spacetimes. Numerical experiments explore the stability and computational efficiency of our method.

Uniform Convergence of Interlaced Euler Method for Stiff Stochastic Differential Equations

In contrast to stiff deterministic systems of ordinary differential equations, in general, the implicit Euler method for stiff stochastic differential equations is not effective. This paper introduces a new numerical method for stiff differential equations which consists of interlacing large implicit Euler time steps with a sequence of small explicit Euler time steps. We emphasize that uniform convergence with respect to the time scale separation parameter e is a desirable property of a stiff solver. We prove that the means and variances of this interlaced method converge uniformly in e for a suitably chosen test problem. We also illustrate the effectiveness of this method via some numerical examples.

An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm

This paper discusses a novel fully implicit formulation for a one-dimensional electrostatic particle-in-cell (PIC) plasma simulation approach. Unlike earlier implicit electrostatic PIC approaches (which are based on a linearized Vlasov–Poisson formulation), ours is based on a nonlinearly converged Vlasov–Ampére (VA) model. By iterating particles and fields to a tight nonlinear convergence tolerance, the approach features superior stability and accuracy properties, avoiding most of the accuracy pitfalls in earlier implicit PIC implementations. In particular, the formulation is stable against temporal (Courant–Friedrichs–Lewy) and spatial (aliasing) instabilities. It is charge- and energy-conserving to numerical round-off for arbitrary implicit time steps (unlike the earlier “energy-conserving” explicit PIC formulation, which only conserves energy in the limit of arbitrarily small time steps). While momentum is not exactly conserved, errors are kept small by an adaptive particle sub-stepping orbit integrator, which is instrumental to prevent particle tunneling (a deleterious effect for long-term accuracy). The VA model is orbit-averaged along particle orbits to enforce an energy conservation theorem with particle sub-stepping. As a result, very large time steps, constrained only by the dynamical time scale of interest, are possible without accuracy loss. Algorithmically, the approach features a Jacobian-free Newton–Krylov solver. A main development in this study is the nonlinear elimination of the new-time particle variables (positions and velocities). Such nonlinear elimination, which we term particle enslavement, results in a nonlinear formulation with memory requirements comparable to those of a fluid computation, and affords us substantial freedom in regards to the particle orbit integrator. Numerical examples are presented that demonstrate the advertised properties of the scheme. In particular, long-time ion acoustic wave simulations show that numerical accuracy does not degrade even with very large implicit time steps, and that significant CPU gains are possible.

Numerical approximation with Nitscheʼs coupling of transient Stokesʼ/Darcyʼs flow problems applied to hemodynamics

We consider a time dependent coupled Stokes/Darcy flow problem and study an approximation method based on a unified finite element scheme complemented with implicit time stepping. Our finite element formulation relies on a weighing strategy in which the physical and discretization parameters are taken into account to robustly enforce interface and boundary conditions by means of the Nitsche method. We study absolute stability and convergence of the scheme, and discuss the algebraic properties of the associated discrete problem. Finally, we present numerical experiments confirming the predicted convergence behavior and algebraic properties, and report an application to the computational analysis of blood flow and plasma filtration in arteries after the implantation of a vascular graft.

On the effect of cluster resolution in riser flows on momentum and reaction kinetic interaction

Fine grid simulations of reactive gas–solid flows in a riser were carried out using an Eulerian multi-fluid kinetic theory of granular flow (KTGF) approach. A translationally periodic section of the riser was used to replicate experimental data collected in the fully developed region of a tall riser. The spatial and temporal resolution was varied in designed experiments to find an appropriate compromise between overall numerical accuracy and computational time. Results revealed that, when first order implicit timestepping is used, no timestep independence could be reached with timestep sizes that are practically feasible. Timestep independence could only be achieved by using second order implicit timestepping. Grid independence was studied in terms of cell width and cell aspect ratio. Solution independence could be reached at a cell width in the range of 10 particle diameters, but no complete grid aspect ratio independence could be achieved. Results suggested that grids with an aspect ratio smaller than one might be necessary to attain grid independent solutions. When sufficiently fine grids are used, however, the effect of a change in aspect ratio is sufficiently small to attain accurate solutions with an aspect ratio of two or lower. Certain important conversion measures were identified for scaling between simulation results collected in a 3D cylindrical domain and those predicted by a 2D planar simulation. System behavior predicted using these scaling rules agreed well with experimental results.

Three-phase numerical model of water migration in partially frozen geological media: model formulation, validation, and applications

Water in the subsurface of the Earth’s cold regions—and possibly the subsurface of Mars—resides in the liquid, vapor, and ice phases. However, relatively few simulations addressing full three-phase, nonisothermal water dynamics in below-freezing porous media have been undertaken. This paper presents a nonisothermal, three-phase approach to modeling water migration in partially frozen porous media. Conservation equations for water (as ice, liquid, and vapor) and a single gas species (in the gas phase and dissolved in water) are coupled to a heat transport equation and solved by a finite-volume method with fully implicit time stepping. Particular attention is given to the method of spatial differencing when the pore space is partially filled with ice. The numerical model is able to reproduce freezing-induced water redistribution observed in laboratory experiments. Simulations of Earth permafrost dynamics and of the formation and evolution of a planetary-scale cryosphere on Mars demonstrate the new capabilities.

The immersed boundary method for advection–electrodiffusion with implicit timestepping and local mesh refinement

We describe an immersed boundary method for problems of fluid–solute-structure interaction. The numerical scheme employs linearly implicit timestepping, allowing for the stable use of timesteps that are substantially larger than those permitted by an explicit method, and local mesh refinement, making it feasible to resolve the steep gradients associated with the space charge layers as well as the chemical potential, which is used in our formulation to control the permeability of the membrane to the (possibly charged) solute. Low Reynolds number fluid dynamics are described by the time-dependent incompressible Stokes equations, which are solved by a cell-centered approximate projection method. The dynamics of the chemical species are governed by the advection–electrodiffusion equations, and our semi-implicit treatment of these equations results in a linear system which we solve by GMRES preconditioned via a fast adaptive composite-grid (FAC) solver. Numerical examples demonstrate the capabilities of this methodology, as well as its convergence properties.

Error estimation and adaptation for functional outputs in time-dependent flow problems

The paper presents an adjoint-based approach for determining global error in the time domain that is relevant to functional outputs from unsteady flow simulations. The algorithm is derived for the unsteady Euler equations that are discretized for second-order accuracy in both space and time and takes into account the effect of dynamic meshes. In addition to error due to temporal resolution, the formulation also takes into account algebraic error arising from partial convergence of the governing equations at each implicit time-step. The resulting error distributions are then used to drive adaptation of the temporal resolution and the convergence tolerances for the governing equations at each time-step. The method is demonstrated in the context of both time-integrated and instantaneous functionals and the results are compared against traditional adaptation methods.

3D phase-field simulations of interfacial dynamics in Newtonian and viscoelastic fluids

This work presents a three-dimensional finite-element algorithm, based on the phase-field model, for computing interfacial flows of Newtonian and complex fluids. A 3D adaptive meshing scheme produces fine grid covering the interface and coarse mesh in the bulk. It is key to accurate resolution of the interface at manageable computational costs. The coupled Navier–Stokes and Cahn–Hilliard equations, plus the constitutive equation for non-Newtonian fluids, are solved using second-order implicit time stepping. Within each time step, Newton iteration is used to handle the nonlinearity, and the linear algebraic system is solved by preconditioned Krylov methods. The phase-field model, with a physically diffuse interface, affords the method several advantages in computing interfacial dynamics. One is the ease in simulating topological changes such as interfacial rupture and coalescence. Another is the capability of computing contact line motion without invoking ad hoc slip conditions. As validation of the 3D numerical scheme, we have computed drop deformation in an elongational flow, relaxation of a deformed drop to the spherical shape, and drop spreading on a partially wetting substrate. The results are compared with numerical and experimental results in the literature as well as our own axisymmetric computations where appropriate. Excellent agreement is achieved provided that the 3D interface is adequately resolved by using a sufficiently thin diffuse interface and refined grid. Since our model involves several coupled partial differential equations and we use a fully implicit scheme, the matrix inversion requires a large memory. This puts a limit on the scale of problems that can be simulated in 3D, especially for viscoelastic fluids.

Numerical solution of the Hamilton–Jacobi–Bellman formulation for continuous time mean variance asset allocation

We solve the optimal asset allocation problem using a mean variance approach. The original mean variance optimization problem can be embedded into a class of auxiliary stochastic linear-quadratic (LQ) problems using the method in Zhou and Li (2000) and Li and Ng (2000). We use a finite difference method with fully implicit timestepping to solve the resulting nonlinear Hamilton–Jacobi–Bellman (HJB) PDE, and present the solutions in terms of an efficient frontier and an optimal asset allocation strategy. The numerical scheme satisfies sufficient conditions to ensure convergence to the viscosity solution of the HJB PDE. We handle various constraints on the optimal policy. Numerical tests indicate that realistic constraints can have a dramatic effect on the optimal policy compared to the unconstrained solution.

A Meshless Solution Technique for the Solution of 3D Unsaturated Zone Problems, Based on Local Hermitian Interpolation with Radial Basis Functions

Recent developments in meshless numerical methods have led to algorithms that can be used to solve arbitrarily large problems without the support of a connected mesh, and without the computational cost and numerical ill-conditioning issues usually associated with such solution techniques. This work applies the Local Hermitian Interpolation (LHI) method, based on local interpolation with Radial Basis Functions (RBFs), to the solution of 3D unsaturated porous media problems. The proposed implementation is capable of handling real soil properties, provided either as an analytical function or as a series of pointwise measurements. The technique is implemented with implicit and explicit timestepping, and is validated against two transient Richards’ equation models, of which one has a known analytical solution. In addition, a real-world infiltration problem based on a saturated–unsaturated formulation is modelled, using a realistic variation of soil properties with water-pressure.

Fast computation of multiphase flow in porous media by implicit discontinuous Galerkin schemes with optimal ordering of elements

We present a family of implicit discontinuous Galerkin schemes for purely advective multiphase flow in porous media in the absence of gravity and capillary forces. To advance the solution one time step, one must solve a discrete system of nonlinear equations. By reordering the grid cells, the nonlinear system can be shown to have a lower triangular block structure, where each block corresponds to the degrees-of-freedom in a single or a small number of cells. To reorder the system, we view the grid cells and the fluxes over cell interfaces as vertices and edges in a directed graph and use a standard topological sorting algorithm. Then the global system can be computed by processing the blocks sequentially using a standard Newton–Raphson algorithm for the degrees-of-freedom in each block. Decoupling the system offers greater control over the nonlinear solution procedure and reduces the computational costs, memory requirements, and complexity of the scheme significantly. In particular, the first-order version of the method may be at least as efficient as modern streamline methods when accuracy requirements or the dynamics of the flow allow for large implicit time steps.

An optimal, parallel, fully implicit Newton–Krylov solver for three-dimensional viscoresistive magnetohydrodynamics

The conceptual development and implementation of a scalable nonlinear solver for the time-dependent three-dimensional (3D) compressible resistive magnetohydrodynamics model (MHD) is discussed. The approach is based on Jacobian-free Newton–Krylov technology, preconditioned with multigrid methods for algorithmic scalability. The key to the approach is the reformulation of the hyperbolic MHD system into a parabolic one, which is amenable to multigrid techniques. Such reformulation (parabolization) aims to render the modified system block diagonally dominant (unlike the original MHD system, which is diagonally submissive for implicit time steps Δt larger than the explicit Courant–Friedrichs–Lewy time step). The algorithm has been tested on a variety of 2D and 3D configurations which demonstrate its excellent algorithmic scalability properties. In particular, it is shown that, serially, the CPU time for a given simulation scales linearly with the number of unknowns, and that large CPU gains (∼30) are attainable with respect to an explicit approach even for moderate grids (256×256). In parallel, the algorithm features excellent scalability properties up to thousands of processors (4096) and millions of unknowns (∼1.3×108).

Electric and magnetic losses modeled by a stable hybrid with explicit–implicit time-stepping for Maxwell’s equations

A stable hybridization of the finite-element method (FEM) and the finite-difference time-domain (FDTD) scheme for Maxwell’s equations with electric and magnetic losses is presented for two-dimensional problems. The hybrid method combines the flexibility of the FEM with the efficiency of the FDTD scheme and it is based directly on Ampère’s and Faraday’s law. The electric and magnetic losses can be treated implicitly by the FEM on an unstructured mesh, which allows for local mesh refinement in order to resolve rapid variations in the material parameters and/or the electromagnetic field. It is also feasible to handle larger homogeneous regions with losses by the explicit FDTD scheme connected to an implicitly time-stepped and lossy FEM region. The hybrid method shows second-order convergence for smooth scatterers. The bistatic radar cross section (RCS) for a circular metal cylinder with a lossy coating converges to the analytical solution and an accuracy of 2% is achieved for about 20 points per wavelength. The monostatic RCS for an airfoil that features sharp corners yields a lower order of convergence and it is found to agree well with what can be expected for singular fields at the sharp corners. A careful convergence study with resolutions from 20 to 140 points per wavelength provides accurate extrapolated results for this non-trivial test case, which makes it possible to use as a reference problem for scattering codes that model both electric and magnetic losses.

Convergence of Implicit Monotone Schemes with Applications in Multiphase Flow in Porous Media

Phase-based upstreaming, which is a commonly used spatial discretization for multiphase flow in reservoir simulation, naturally gives rise to implicit monotone schemes when implicit time-stepping is used. The nonlinear Gauss-Seidel and Jacobi algorithms are shown to converge to a unique bounded solution when applied to the resulting system of equations. Thus, for one-dimensional problems, we obtain an alternate, constructive proof that such schemes are well-defined and converge to the entropy solution of the conservation law for arbitrary CFL numbers. The accuracy of phase-based upstream solutions is studied for various spatial and temporal grids, revealing the importance of unconditional stability when nonuniform grids and/or variable porosity is involved.