We extend the immersed boundary (IB) method to simulate the dynamics of a 2D dry foam by including the topological changes of the bubble network. In the article [Y. Kim, M.-C. Lai, and C. S. Peskin, J. Comput. Phys. 229:5194-5207,2010], we implemented an IB method for the foam problem in the two-dimensional case, and tested it by verifying the von Neumann relation which governs the coarsening of a two-dimensional dry foam. However, the method implemented in that article had an important limitation; we did not allow for the resolution of quadruple or higher order junctions into triple junctions. A total shrinkage of a bubble with more than four edges generates a quadruple or higher order junction. In reality, a higher order junction is unstable and resolves itself into triple junctions. We here extend the methodology previously introduced by allowing topological changes, and we illustrate the significance of such topological changes by comparing the behaviors of foams in which topological changes are allowed to those in which they are not.

We present an algorithm for the stochastic simulation of gene expression and heterogeneous population dynamics. The algorithm combines an exact method to simulate molecular-level fluctuations in single cells and a constant-number Monte Carlo method to simulate time-dependent statistical characteristics of growing cell populations. To benchmark performance, we compare simulation results with steady-state and time-dependent analytical solutions for several scenarios, including steady-state and time-dependent gene expression, and the effects on population heterogeneity of cell growth, division, and DNA replication. This comparison demonstrates that the algorithm provides an efficient and accurate approach to simulate how complex biological features influence gene expression. We also use the algorithm to model gene expression dynamics within “bet-hedging” cell populations during their adaption to environmental stress. These simulations indicate that the algorithm provides a framework suitable for simulating and analyzing realistic models of heterogeneous population dynamics combining molecular-level stochastic reaction kinetics, relevant physiological details and phenotypic variability.

A robust immersed boundary-lattice Boltzmann method (IB-LBM) is proposed to simulate fluid-structure interaction (FSI) problems in this work. Compared with the conventional IB-LBM, the current method employs the fractional step technique to solve the lattice Boltzmann equation (LBE) with a forcing term. Consequently, the non-physical oscillation of body force calculation, which is frequently encountered in the traditional IB-LBM, is suppressed greatly. It is of importance for the simulation of FSI problems. In the meanwhile, the no-slip boundary condition is strictly satisfied by using the velocity correction scheme. Moreover, based on the relationship between the velocity correction and forcing term, the boundary force can be calculated accurately and easily. A few test cases are first performed to validate the current method. Subsequently, a series of FSI problems, including the vortex-induced vibration of a circular cylinder, an elastic filament flapping in the wake of a fixed cylinder and sedimentation of particles, are simulated. Based on the good agreement between the current results and those in the literature, it is demonstrated that the proposed IB-LBM has the capability to handle various FSI problems effectively.

We describe a relativistic semi-Lagrangian scheme for the numerical solution of the relativistic Vlasov-Maxwell system. The implementation strategy on a modern non-unified memory access (NUMA) architecture using the OpenMP framework is discussed. We demonstrated that close to perfect scaling can be obtained on modern many-core, multi-socket systems. Application of this code to the problem of relativistic generation of high-harmonic laser radiation is demonstrated. The results are compared to particle-in-cell (PIC) simulations, indicating in particular that for warm plasma the Vlasov simulation is superior. We discuss the impact of plasma temperature on the radiation spectrum and show that the efficiency of harmonic generation depends on the plasma temperature.

We investigate the effectiveness of two extrapolation-based methods aiming to approximate the initial state required by an iterative solver in simulations of unsteady flow problems. The methods lead to about a ten-fold reduction in the iteration count while requiring only negligible computational overhead. They are particularly suitable for parallel computing since they are based almost exclusively on data stored locally on each processor. Performance has been evaluated in simulations of turbulent flow in a stenosed carotid artery and also in laminar flow in a very large domain containing the human intracranial arterial tree.

This paper is concerned with numerical methods for the Navier-Stokes-Nernst-Planck-Poisson equation system. The main goal is to construct and analyze some stable time stepping schemes for the time discretization and use a spectral method for the spatial discretization. The main contribution of the paper includes: 1) an useful stability inequality for the weak solution is derived; 2) a first order time stepping scheme is constructed, and the non-negativity of the concentration components of the discrete solution is proved. This is an important property since the exact solution shares the same property. Moreover, the stability of the scheme is established, together with a stability condition on the time step size; 3) a modified first order scheme is proposed in order to decouple the calculation of the velocity and pressure in the fluid field. This new scheme equally preserves the non-negativity of the discrete concentration solution, and is stable under a similar stability condition; 4) a stabilization technique is introduced to make the above mentioned schemes stable without restriction condition on the time step size; 5) finally we construct a second order finite difference scheme in time and spectral discretization in space. The numerical tests carried out in the paper show that all the proposed schemes possess some desirable properties, such as conditionally/unconditionally stability, first/second order convergence, non-negativity of the discrete concentrations, and so on.

In the present work, a new type of coupled compact difference scheme has been proposed for the solution of computational acoustics and flow problems. The proposed scheme evaluates the first, the second and the fourth derivative terms simultaneously. Derived compact difference scheme has a significant spectral resolution and a physical dispersion relation preserving (DRP) ability over a considerable wavenumber range when a fourth order four stage Runge-Kutta scheme is used for the time integration. Central stencil has been used for the present numerical scheme to evaluate spatial derivative terms. Derived scheme has the capability of adding numerical diffusion adaptively to attenuate spurious high wavenumber oscillations responsible for numerical instabilities. The DRP nature of the proposed scheme across a wider wavenumber range provides accurate results for the model wave equations as well as computational acoustic problems. In addition to the attractive feature of adaptive diffusion, present scheme also helps to control spurious reflections from the domain boundaries and is projected as an alternative to the perfectly matched layer (PML) technique.

We study the computation of ground states and time dependent solutions of the Schrödinger-Poisson system (SPS) on a bounded domain in 2D (i.e. in two space dimensions). On a disc-shaped domain, we derive exact artificial boundary conditions for the Poisson potential based on truncated Fourier series expansion in θ, and propose a second order finite difference scheme to solve the r-variable ODEs of the Fourier coefficients. The Poisson potential can be solved within (MNlogN) arithmetic operations where M,N are the number of grid points in r-direction and the Fourier bases. Combined with the Poisson solver, a backward Euler and a semi-implicit/leap-frog method are proposed to compute the ground state and dynamics respectively. Numerical results are shown to confirm the accuracy and efficiency. Also we make it clear that backward Euler sine pseudospectral (BESP) method in [33] can not be applied to 2D SPS simulation.

The immersed boundary method (IBM) has been popular in simulating fluid structure interaction (FSI) problems involving flexible structures, and the recent introduction of the lattice Boltzmann method (LBM) into the IBM makes the method more versatile. In order to test the coupling characteristics of the IBM with the multiple-relaxation-time LBM (MRT-LBM), the three-dimensional (3D) balloon dynamics, including inflation, release and breach processes, are simulated. In this paper, some key issues in the coupling scheme, including the discretization of 3D boundary surfaces, the calculation of boundary force density, and the introduction of external force into the LBM, are described. The good volume conservation and pressure retention properties are verified by two 3D cases. Finally, the three FSI processes of a 3D balloon dynamics are simulated. The large boundary deformation and oscillation, obvious elastic wave propagation, sudden stress release at free edge, and recoil phenomena are all observed. It is evident that the coupling scheme of the IBM and MRT-LBM can handle complicated 3D FSI problems involving large deformation and large pressure gradients with very good accuracy and stability.

Knowledge of helium diffusion kinetics is critical for materials in which helium measurements are made, particulary for thermochronology. In most cases the helium ages were younger than expected, an observation attributes to diffusive loss of helium and the ejection of high energy alpha particles. Therefore it is important to accurately calculate the distribution of the source term within a sample. In this paper, the prediction of the helium concentrations as function of a spatially variable source term are considered. Both the forward and inverse solutions are presented. Under the assumption of radially symmetric geometry, an analytical solution is deduced based on the eigenfunction expansion. Two regularization methods, the Tikhonov regularization and the spectral cutoff regularization, are considered to obtain the regularized solution. Error estimates with optimal convergence order are shown between the exact solution and the regularized solution. Numerical examples are presented to illustrate the validity and effectiveness of the proposed methods

In this paper, numerical simulations of solid particle erosion phenomena on H-type finned circular/elliptic tube surface, which is of great significance to the antiwear design of heat exchanger, are presented. The Eulerian-Lagrangian approach is applied to simulate the dilute gas-solid flow through H-type finned circular/elliptic tubes. A semi-empirical model is adopted to predict the erosion rate. The dynamics behavior of the entrained solid particles in the flow is presented. The geometry of eroded tube surface is changed with the predicted erosion which is taken into account by a UDF and the flow field is re-solved for the eroded tube surface at every time step. The influences of ten parameters (the tube bundle arrangement, particle size, particle concentration, fluid Reynolds number, fin thickness, fin pitch, fin length, fin width, slit width and the transverse tube pitch) on the maximum erosion depth of the H-type circular/elliptic finned tube surface are investigated. Using H-type finned elliptic tube surface can effectively reduce the erosion rate of tube surface comparedwith that using H-type finned circular tube surface. The erosion in in-line arrangement is less severe than that in staggered arrangement. With the increase of particle size, particle concentration and the fluid Reynolds number, the erosion rate of the tube surface rises. The numerically predicted effect of Reynolds number is in good agreement with previous test data. Among the six geometry parameters, the most influential parameter is the transverse tube pitch.

Numerical error caused by “ghost forces” in a quasicontinuum method is studied in the context of dynamic problems. The error in the discrete W1,∞ norm is analyzed for the time scale (ε) and the time scale (1) with ε being the lattice spacing.

Variational image segmentation based on the Mumford and Shah model [31], together with implementation by the piecewise constant level-set method (PCLSM) [26], leads to fully nonlinear Total Variation (TV)-Allen-Cahn equations. The commonly-used numerical approaches usually suffer from the difficulties not only with the non-differentiability of the TV-term, but also with directly evolving the discontinuous piecewise constant-structured solutions. In this paper, we propose efficient dual algorithms to overcome these drawbacks. The use of a splitting-penalty method results in TV-Allen-Cahn type models associated with different “double-well” potentials, which allow for the implementation of the dual algorithm of Chambolle [8]. Moreover, we present a new dual algorithm based on an edge-featured penalty of the dual variable, which only requires to solve a vectorial Allen-Cahn type equation with linear ∇(div)-diffusion rather than fully nonlinear diffusion in the Chambolle’s approach. Consequently, more efficient numerical algorithms such as time-splitting method and Fast Fourier Transform (FFT) can be implemented. Various numerical tests show that two dual algorithms are much faster and more stable than the primal gradient descent approach, and the new dual algorithm is at least as efficient as the Chambolle’s algorithm but is more accurate. We demonstrate that the new method also provides a viable alternative for image restoration.

The incompatibilities between the initial and boundary data will cause singularities at the time-space corners, which in turn adversely affect the accuracy of the numerical schemes used to compute the solutions. We study the corner singularity issue for nonlinear evolution equations in 1D, and propose two remedy procedures that effectively recover much of the accuracy of the numerical scheme in use. Applications of the remedy procedures to the 1D viscous Burgers equation, and to the 1D nonlinear reaction-diffusion equation are presented. The remedy procedures are applicable to other nonlinear diffusion equations as well.

We have introduced a fully second order IMplicit/EXplicit (IMEX) time in-tegration technique for solving the compressible Euler equations plus nonlinear heat conduction problems (also known as the radiation hydrodynamics problems) in Kadioglu et al., J. Comp. Physics [22,24]. In this paper, we study the implications when this method is applied to the incompressible Navier-Stokes (N-S) equations. The IMEX method is applied to the incompressible flow equations in the following manner. The hyperbolic terms of the flow equations are solved explicitly exploiting the well understood explicit schemes. On the other hand, an implicit strategy is employed for the non-hyperbolic terms. The explicit part is embedded in the implicit step in such a way that it is solved as part of the non-linear function evaluation within the framework of the Jacobian-Free Newton Krylov (JFNK) method [8,29,31]. This is done to obtain a self-consistent implementation of the IMEX method that eliminates the potential order reduction in time accuracy due to the specific operator separation. We employ a simple yet quite effective fractional step projection methodology (similar to those in [11,19,21,30]) as our preconditioner inside the JFNK solver. We present results from several test calculations. For each test, we show second order time convergence. Finally, we present a study for the algorithm performance of the JFNK solver with the new projection method based preconditioner.

Realistic and accurate computer simulations of the particle-hemodynamics in arterial systems can be a valuable tool for numerous biomedical applications. Examples include optimal by-pass grafting and optimal drug-delivery, as well as best medical management concerning the cardio-vascular system. However, such numerical analyses require large computer resources which may become prohibitive for extended sets of arterial bifurcations. A remedy is to develop a hybrid model where the first few generations of the bifurcating arteries of interest are simulated in full 3-D, while a 1-D model is then coupled for subsequent bifurcations. Alternatively, a 1-D computer model can be directly employed to simulate fluid-particle transport in complex bifurcating networks.

Relying on a representative axial velocity profile, a physiological 1-D model has been developed and validated, which is capable of predicting with reasonable accuracy arterial flow, pressure field and elastic wall interaction as well as particle transport. The usefulness of the novel 1-D simulation approach is demonstrated via a comparison to 3-D blood flow and microsphere transport in a hepatic artery system, featuring as outlets one major branch and four small daughter vessels. Compared to the 3-D simulation, the 1-D analysis requires only about 1% of computational time. The hybrid modeling approach would be also applicable to the human respiratory tract to evaluate the fate of inhaled aerosols.

A simple and cost-effective way to simulate particle-hemodynamics is using a 1-D model for simulating arterial pressures and flow rates as well as microsphere transport, based on assumptions involving the use of a simple algebraic pressure-area relation, an exponential elasticity model for the vessels, and considering only unidirectional flow with a representative skewed velocity profile. In summary, the novel contributions are:

• Particle tracking in arteries via 1-D fluid modeling and selection of an averaged, skewed velocity profile based on 3-D simulation results to provide more realistic friction and inertia term values for modeling a flow system with bifurcations.

• The 1-D model can be coupled to a 3-D model so that simulations can be run for larger regions of vascular or lung-airway systems.

Three high order shock-capturing schemes are compared for large eddy simulations (LES) of temporally evolving mixing layers for different convective Mach numbers ranging from the quasi-incompressible regime to highly compressible supersonic regime. The considered high order schemes are fifth-order WENO (WENO5), seventh-order WENO (WENO7) and the associated eighth-order central spatial base scheme with the dissipative portion of WENO7 as a nonlinear post-processing filter step (WENO7fi). This high order nonlinear filter method of Yee & Sjögreen is designed for accurate and efficient simulations of shock-free compressible turbulence, turbulence with shocklets and turbulence with strong shocks with minimum tuning of scheme parameters. The LES results by WENO7fi using the same scheme parameter agree well with experimental results compiled by Barone et al., and published direct numerical simulations (DNS) work of Rogers & Moser and Pantano & Sarkar, whereas results by WENO5 and WENO7 compare poorly with experimental data and DNS computations.

In this study, a numerical framework of the high order well-balanced weighted compact nonlinear (WCN) schemes is proposed for the shallow water equations based on the work in [S. Zhang, S. Jiang, C.-W Shu, J. Comput. Phys. 227 (2008) 7294-7321]. We employ a special splitting technique for the source term proposed in [Y. Xing, C.-W Shu, J. Comput. Phys. 208 (2005) 206-227] to maintain the exact C-property, which can be proved theoretically. In the meantime, the genuine high order accuracy of the numerical scheme can be observed successfully, and small perturbation of the stationary state can be resolved and evolved well. In order to capture the strong discontinuities and large gradients, the fifth-order upwind weighted nonlinear interpolations together with the fourth/sixth order cell-centered compact scheme are used to construct different WCN schemes. In addition, the local characteristic projections are considered to further restrain the potential numerical oscillations. A variety of representative one- and two-dimensional examples are tested to demonstrate the good performance of the proposed schemes.

Amethod for non-rigid image registration that is suitable for large deformations is presented. Conventional registration methods embed the image in a B-spline object, and the image is evolved by deforming the B-spline object. In this work, we represent the image using B-spline and deform the image using a composition approach. We also derive a computationally efficient algorithm for calculating the B-spline coefficients and gradients of the image by adopting ideas from signal processing using image filters. We demonstrate the application of our method on several different types of 2D and 3D images and compare it with existing methods.

This paper presents a GPU-accelerated implementation of the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) method with an inexact nullspace filtering approach to find eigenvalues in electromagnetics analysis with higher-order FEM. The performance of the proposed approach is verified using the Kepler (Tesla K40c) graphics accelerator, and is compared to the performance of the implementation based on functions from the Intel MKL on the Intel Xeon (E5-2680 v3, 12 threads) central processing unit (CPU) executed in parallel mode. Compared to the CPU reference implementation based on the Intel MKL functions, the proposed GPU-based LOBPCG method with inexact nullspace filtering allowed us to achieve up to 2.9-fold acceleration.