Spatially adaptive stochastic numerical methods for intrinsic fluctuations in reaction–diffusion systems

Stochastic partial differential equations are introduced for the continuum concentration fields of reaction–diffusion systems. The stochastic partial differential equations account for fluctuations arising from the finite number of molecules which diffusively migrate and react. Spatially adaptive stochastic numerical methods are developed for approximation of the stochastic partial differential equations. The methods allow for adaptive meshes with multiple levels of resolution, Neumann and Dirichlet boundary conditions, and domains having geometries with curved boundaries. A key issue addressed by the methods is the formulation of consistent discretizations for the stochastic driving fields at coarse-refined interfaces of the mesh and at boundaries. Methods are also introduced for the efficient generation of the required stochastic driving fields on such meshes. As a demonstration of the methods, investigations are made of the role of fluctuations in a biological model for microorganism direction sensing based on concentration gradients. Also investigated, a mechanism for spatial pattern formation induced by fluctuations. The discretization approaches introduced for SPDEs have the potential to be widely applicable in the development of numerical methods for the study of spatially extended stochastic systems.     

Coupling p-multigrid to geometric multigrid for discontinuous Galerkin formulations of the convection–diffusion equation

An improved p-multigrid algorithm for discontinuous Galerkin (DG) discretizations of convection–diffusion problems is presented. The general p  -multigrid algorithm for DG discretizations involves a restriction from the p=1$p=1$ to p=0$p=0$ discontinuous polynomial solution spaces. This restriction is problematic and has limited the efficiency of the p  -multigrid method. For purely diffusive problems, Helenbrook and Atkins have demonstrated rapid convergence using a method that restricts from a discontinuous to continuous polynomial solution space at p=1$p=1$. It is shown that this method is not directly applicable to the convection–diffusion (CD) equation because it results in a central-difference discretization for the convective term. To remedy this, ideas from the streamwise upwind Petrov–Galerkin (SUPG) formulation are used to devise a transition from the discontinuous to continuous space at p=1$p=1$ that yields an upwind discretization. The results show that the new method converges rapidly for all Peclet numbers.     

A First-Passage Kinetic Monte Carlo algorithm for complex diffusion–reaction systems

We develop an asynchronous event-driven First-Passage Kinetic Monte Carlo (FPKMC) algorithm for continuous time and space systems involving multiple diffusing and reacting species of spherical particles in two and three dimensions. The FPKMC algorithm presented here is based on the method introduced in Oppelstrup et al. [10] and is implemented in a robust and flexible framework. Unlike standard KMC algorithms such as the n-fold algorithm, FPKMC is most efficient at low densities where it replaces the many small hops needed for reactants to find each other with large first-passage hops sampled from exact time-dependent Green’s functions, without sacrificing accuracy. We describe in detail the key components of the algorithm, including the event-loop and the sampling of first-passage probability distributions, and demonstrate the accuracy of the new method. We apply the FPKMC algorithm to the challenging problem of simulation of long-term irradiation of metals, relevant to the performance and aging of nuclear materials in current and future nuclear power plants. The problem of radiation damage spans many decades of time-scales, from picosecond spikes caused by primary cascades, to years of slow damage annealing and microstructure evolution. Our implementation of the FPKMC algorithm has been able to simulate the irradiation of a metal sample for durations that are orders of magnitude longer than any previous simulations using the standard Object KMC or more recent asynchronous algorithms.     

Stability of operator splitting methods for systems with indefinite operators: Advection–diffusion–reaction systems

This brief paper presents an A-stability result for operator splitting type time integration methods applied to advection–diffusion–reaction equations with possibly indefinite source terms. These results extend our earlier work on diffusion–reaction systems [D.L. Ropp, J.N. Shadid, Stability of operator splitting methods for systems with indefinite operators: reaction–diffusion systems, J. Comput. Phys. 203 (2) (2005) 449–466]. The A-stability result presents sufficient conditions that control both low and high wave number instabilities. A corollary shows that if L-stable methods are used for the diffusion term the high wave number instability will be controlled more easily. Numerical results are presented that verify second-order convergence for the operator splitting methods and demonstrate control of instabilities on a chemotaxis problem by use of an L-stable diffusion integrator.     

Quantitative study of diffusion jumps in atomistic simulations of model gas–polymer systems

Microscopic mechanisms underlying the diffusion of particles in polymeric and other systems include ‘jumps’ that are said to provide for a substantial contribution to the overall particle displacement. Such jumps have been observed in molecular simulations and experimentally, leading to important qualitative conclusions. An efficient method has been proposed for the identification and quantitative processing of jumps, and successfully employed in simulations of gas–liquid alkane systems. In the present work, the same method is applied in equilibrium Molecular Dynamics simulations of methane-like molecules dispersed in polymer-like alkanes, at atmospheric pressure and temperature well above the polymer melting point. The systems studied were prepared and equilibrated and a linear diffusion regime was confirmed by means of various criteria. The occurrence of distinct jump events was clearly revealed and their average length and frequency were calculated. In this way, a random-walk-type diffusion coefficient, D _s,?jumps, of the penetrants, was obtained and found to be lower than the overall diffusion coefficient D _s,?MSD calculated by the mean square displacement method. This is a strong indication that the overall diffusion is a combination of longer jumps with other microscopic mechanisms such as smoother and more gradual displacements effected upon the diffusing particle by its surroundings.     

A positivity-preserving ALE finite element scheme for convection–diffusion equations in moving domains

A new high-resolution scheme is developed for convection–diffusion problems in domains with moving boundaries. A finite element approximation of the governing equation is designed within the framework of a conservative Arbitrary Lagrangian Eulerian (ALE) formulation. An implicit flux-corrected transport (FCT) algorithm is implemented to suppress spurious undershoots and overshoots appearing in convection-dominated problems. A detailed numerical study is performed for P₁ finite element discretizations on fixed and moving meshes. Simulation results for a Taylor dispersion problem (moderate Peclet numbers) and for a convection-dominated problem (large Peclet numbers) are presented to give a flavor of practical applications.     

Adaptive mesh refinement for stochastic reaction–diffusion processes

We present an algorithm for adaptive mesh refinement applied to mesoscopic stochastic simulations of spatially evolving reaction–diffusion processes. The transition rates for the diffusion process are derived on adaptive, locally refined structured meshes. Convergence of the diffusion process is presented and the fluctuations of the stochastic process are verified. Furthermore, a refinement criterion is proposed for the evolution of the adaptive mesh. The method is validated in simulations of reaction–diffusion processes as described by the Fisher–Kolmogorov and Gray–Scott equations.     

Oscillations vs. chaotic waves: Attractor selection in bistable stochastic reaction–diffusion systems

In bistable systems, the long-term behavior of solutions depends on the location of the initial conditions. In a deterministic setting, where the initial condition is kept fixed in one particular basin of attraction, repeated numerical simulations will always lead to the same long-term behavior. The other possible asymptotic solution type will never be observed. This clear distinction does not hold anymore if the system is forced by random fluctuations. In this case, both asymptotic solutions can be attained, and the relative frequency of different long-term behaviors observed in many repeated simulation runs will follow a certain probability distribution. We present a simple reaction–diffusion model of spatial predator–prey interaction, where depending on the initial spatial distribution of the two populations either spatially homogeneous or spatiotemporal irregular oscillations may be observed. We show by repeated stochastic simulations that, when starting in the basin of attraction of the spatiotemporal irregular solution, in the randomly forced system the probability to observe spatially homogeneous oscillations instead of spatiotemporally irregular oscillations follows a non-trivial bimodal distribution.     

A self-organizing Lagrangian particle method for adaptive-resolution advection–diffusion simulations

We present a novel adaptive-resolution particle method for continuous parabolic problems. In this method, particles self-organize in order to adapt to local resolution requirements. This is achieved by pseudo forces that are designed so as to guarantee that the solution is always well sampled and that no holes or clusters develop in the particle distribution. The particle sizes are locally adapted to the length scale of the solution. Differential operators are consistently evaluated on the evolving set of irregularly distributed particles of varying sizes using discretization-corrected operators. The method does not rely on any global transforms or mapping functions. After presenting the method and its error analysis, we demonstrate its capabilities and limitations on a set of two- and three-dimensional benchmark problems. These include advection–diffusion, the Burgers equation, the Buckley–Leverett five-spot problem, and curvature-driven level-set surface refinement.     

Semi-Lagrangian multistep exponential integrators for index 2 differential–algebraic systems

Implicit-explicit (IMEX) multistep methods are very useful for the time discretization of convection diffusion PDE problems such as the Burgers equations and the incompressible Navier–Stokes equations. In the latter as well as in PDE models of plasma physics and of electromechanical systems, semi-discretization in space gives rise to differential–algebraic (DAE) system of equations often of index higher than 1. In this paper we propose a new class of exponential integrators for index 2 DAEs arising from the semi-discretization of PDEs with a dominating and typically nonlinear convection term. This class of problems includes the incompressible Navier–Stokes equations. The integration methods are based on the backward differentiation formulae (BDF) and they can be applied without modifications in the semi-Lagrangian integration of convection diffusion problems. The approach gives improved performance at low viscosity regimes.     

A fully conservative Eulerian–Lagrangian method for a convection–diffusion problem in a solenoidal field

Tracer transport is governed by a convection–diffusion problem modeling mass conservation of both tracer and ambient fluids. Numerical methods should be fully conservative, enforcing both conservation principles on the discrete level. Locally conservative characteristics methods conserve the mass of tracer, but may not conserve the mass of the ambient fluid. In a recent paper by the authors [T. Arbogast, C. Huang, A fully mass and volume conserving implementation of a characteristic method for transport problems, SIAM J. Sci. Comput. 28 (2006) 2001–2022], a fully conservative characteristic method, the Volume Corrected Characteristics Mixed Method (VCCMM), was introduced for potential flows. Here we extend and apply the method to problems with a solenoidal (i.e., divergence-free) flow field. The modification is a computationally inexpensive simplification of the original VCCMM, requiring a simple adjustment of trace-back regions in an element-by-element traversal of the domain. Our numerical results show that the method works well in practice, is less numerically diffuse than uncorrected characteristic methods, and can use up to at least about eight times the CFL limited time step.     

Setting initial conditions for inflation with reaction–diffusion equation

We discuss the issue of setting appropriate initial conditions for inflation. Specifically, we consider natural inflation model and discuss the fine tuning required for setting almost homogeneous initial conditions over a region of order several times the Hubble size which is orders of magnitude larger than any relevant correlation length for field fluctuations. We then propose to use the special propagating front solutions of reaction–diffusion equations for localized field domains of smaller sizes. Due to very small velocities of these propagating fronts we find that the inflaton field in such a field domain changes very slowly, contrary to naive expectation of rapid roll down to the true vacuum. Continued expansion leads to the energy density in the Hubble region being dominated by the vacuum energy, thereby beginning the inflationary phase. Our results show that inflation can occur even with a single localized field domain of size smaller than the Hubble size. We discuss possible extensions of our results for different inflationary models, as well as various limitations of our analysis (e.g. neglecting self gravity of the localized field domain).     

Moving mesh methods for blowup in reaction–diffusion equations with traveling heat source

In this paper moving mesh methods are used to simulate the blowup in a reaction–diffusion equation with traveling heat source. The finite-time blowup occurs if the speed of the movement of the heat source remains sufficiently low, and the blowup procedure is not fixed at one point not like that for stationary heat source. As time goes to the blowup time, the blowup profile converges to a stationary state. In the simulation a new moving mesh algorithm is designed to deal with the difficulty caused by the delta function in the traveling heat source. The convergence rates are verified and new blowup figures are generated from the numerical experiments.     

Local DG method using WENO type limiters for convection–diffusion problems

The local discontinuous Galerkin (LDG) method is a spatial discretization procedure for convection–diffusion equations, which employs useful features from high resolution finite volume schemes, such as the exact or approximate Riemann solvers serving as numerical fluxes and limiters, which is termed as Runge–Kutta LDG (RKLDG) when TVD Runge–Kutta method is applied for time discretization. It has the advantage of flexibility in handling complicated geometry, h-p adaptivity, and efficiency of parallel implementation and has been used successfully in many applications. However, the limiters used to control spurious oscillations in the presence of strong shocks are less robust than the strategies of essentially non-oscillatory (ENO) and weighted ENO (WENO) finite volume and finite difference methods. In this paper, we investigated RKLDG methods with WENO and Hermite WENO (HWENO) limiters for solving convection–diffusion equations on unstructured meshes, with the goal of obtaining a robust and high order limiting procedure to simultaneously obtain uniform high order accuracy and sharp, non-oscillatory shock transition. Numerical results are provided to illustrate the behavior of these procedures.     

Synchronized stability in a reaction–diffusion neural network model

The reaction–diffusion neural network consisting of a pair of identical tri-neuron loops is considered. We present detailed discussions about the synchronized stability and Hopf bifurcation, deducing the non-trivial role that delay plays in different locations. The corresponding numerical simulations are used to illustrate the effectiveness of the obtained results. In addition, the numerical results about the effects of diffusion reveal that diffusion may speed up the tendency to synchronization and induce the synchronized equilibrium point to be stable. Furthermore, if the parameters are located in appropriate regions, multiple unstability and bistability or unstability and bistability may coexist.     

Two new upwind difference schemes for a coupled system of convection–diffusion equations arising from the steady MHD duct flow problems

In this paper, we develop two new upwind difference schemes for solving a coupled system of convection–diffusion equations arising from the steady incompressible MHD duct flow problem with a transverse magnetic field at high Hartmann numbers. Such an MHD duct flow is convection-dominated and its solution may exhibit localized phenomena such as boundary layers, namely, narrow boundary regions where the solution changes rapidly. Most conventional numerical schemes cannot efficiently solve the layer problems because they are lacking in either stability or accuracy. In contrast, the newly proposed upwind difference schemes can achieve a reasonable accuracy with a high stability, and they are capable of resolving high gradients near the layer regions without refining the grid. The accuracy of the first new upwind scheme is O(h + k) and the second one improves the accuracy to O(ε²(h + k) + ε(h² + k²) + (h³ + k³)), where 0 < ε ? 1/M ? 1 and M is the high Hartmann number. Numerical examples are provided to illustrate the performance of the newly proposed upwind difference schemes.     

An adaptive algorithm for simulation of stochastic reaction–diffusion processes

We propose an adaptive hybrid method suitable for stochastic simulation of diffusion dominated reaction–diffusion processes. For such systems, simulation of the diffusion requires the predominant part of the computing time. In order to reduce the computational work, the diffusion in parts of the domain is treated macroscopically, in other parts with the tau-leap method and in the remaining parts with Gillespie’s stochastic simulation algorithm (SSA) as implemented in the next subvolume method (NSM). The chemical reactions are handled by SSA everywhere in the computational domain. A trajectory of the process is advanced in time by an operator splitting technique and the timesteps are chosen adaptively. The spatial adaptation is based on estimates of the errors in the tau-leap method and the macroscopic diffusion. The accuracy and efficiency of the method are demonstrated in examples from molecular biology where the domain is discretized by unstructured meshes.