Search

Background. Smoke plumes from wildland fires are significant hazards, impacting public health and facilitating fire spotting. In crossflow conditions, typically characterized by vertical shear within the planetary boundary layer, smoke plumes bend downstream of the fire. This bending confines the plume within the boundary layer, delaying its dispersion and increasing fire spotting risks. Aims. This study aims to investigate the bending height of smoke plumes influenced by a sheared crossflow and stratified atmospheric conditions, providing insights into smoke dispersion and fire behavior. Methods. Radar observations from a prescribed fire experiment reveal a large-scale billow-like vorticity pair associated with plume bending. Large Eddy Simulations (LES) in the CM1 model are used to replicate these observations. Key Results. A scaling analysis estimates the plume bending height by incorporating a modified Byram's convective number to account for sheared crossflow. The bending height increases with higher fire intensity, more unstable stratification, and weaker crossflow shear. Simulated results align well with the scaling analysis, highlighting the significance of shear and stratification in plume dynamics. Conclusions. Plume bending height is strongly influenced by fire intensity, atmospheric stratification, and crossflow shear. Incorporating these factors improves our understanding of smoke plume behavior. Implications. The findings enhance smoke dispersion modeling and risk assessment for fire spotting, aiding in the development of strategies to mitigate wildfire impacts.

Data-driven techniques are being increasingly applied to complement physics-based models in fire science. However, the lack of sufficiently large datasets continues to hinder the application of certain machine learning techniques. In this paper, we use simulated data to investigate the ability of neural networks to parameterize dynamics in fire science. In particular, we investigate neural networks that map five key parameters in fire spread to the first arrival time, and the corresponding inverse problem. By using simulated data, we are able to characterize the error, the required dataset size, and the convergence properties of these neural networks. For the inverse problem, we quantify the network's sensitivity in estimating each of the key parameters. The findings demonstrate the potential of machine learning in fire science, highlight the challenges associated with limited dataset sizes, and quantify the sensitivity of neural networks to estimate key parameters governing fire spread dynamics.

The increasing frequency and severity of wildfires highlight the need for accurate fire and plume spread models. We introduce an approach that effectively isolates and tracks fire and plume behavior across various spatial and temporal scales and image types, identifying physical phenomena in the system and providing insights useful for developing and validating models. Our method combines image segmentation and graph theory to delineate fire fronts and plume boundaries. We demonstrate that the method effectively distinguishes fires and plumes from visually similar objects. Results demonstrate the successful isolation and tracking of fire and plume dynamics across various image sources, ranging from synoptic-scale ($10^4$-$10^5$ m) satellite images to sub-microscale ($10^0$-$10^1$ m) images captured close to the fire environment. Furthermore, the methodology leverages image inpainting and spatio-temporal dataset generation for use in statistical and machine learning models.

We present a numerical method for the solution of diffusion problems in unbounded planar regions with complex geometries of absorbing and reflecting bodies. Our numerical method applies the Laplace transform to the parabolic problem, yielding a modified Helmholtz equation which is solved with a boundary integral method. Returning to the time domain is achieved by quadrature of the inverse Laplace transform by deforming along the so-called Talbot contour. We demonstrate the method for various complex geometries formed by disjoint bodies of arbitrary shape on which either uniform Dirichlet or Neumann boundary conditions are applied. The use of the Laplace transform bypasses constraints with traditional time-stepping methods and allows for integration over the long equilibration timescales present in diffusion problems in unbounded domains. Using this method, we demonstrate shielding effects where the complex geometry modulates the dynamics of capture to absorbing sets. In particular, we show examples where geometry can guide diffusion processes to particular absorbing sites, obscure absorbing sites from diffusing particles, and even find the exits of confining geometries, such as mazes.

We numerically investigate the hydrodynamics and membrane dynamics of multicomponent vesicles in two strongly confined geometries. This serves as a simplified model for red blood cells undergoing large deformations while traversing narrow constrictions. We propose a new parameterization for the bending modulus that remains positive for all lipid phase parameter values. For a multicomponent vesicle passing through a stenosis, we establish connections between various properties: lipid phase coarsening, size and flow profile of the lubrication layers, excess pressure, and the tank-treading velocity of the membrane. For a multicomponent vesicle passing through a contracting channel, we find that the lipid always phase separates so that the vesicle is stiffer in the front as it passes through the constriction. For both cases of confinement, we find that lipid coarsening is arrested under strong confinement and proceeds at a high rate upon relief from extreme confinement. The results may be useful for efficient sorting lipid domains using microfluidic flows by controlled release of vesicles passing through strong confinement

Active colloidal systems with non-equilibrium self-organization is a long-standing, challenging area in biology. To understand how hydrodynamic flow may be used to actively control self-assembly of Janus particles (JPs), we use a model recently developed for the many-body hydrodynamics of amphiphilic JPs suspended in a viscous background flow (JFM, 941, 2022). We investigate how various morphologies arise from tuning the hydrophobic distribution of the JP-solvent interface. We find JPs assembled into uni-lamella, multi-lamella and striated structures. To introduce dynamics, we include a linear shear flow and a steady Taylor-Green mixing flow, and measure the collective dynamics of JP particles in terms of their (a) free energy from the hydrophobic interactions between the JPs, (b) order parameter for the ordering of JPs in terms of alignment of their directors, and (c) strain parameter that captures the deformation in the assembly. We characterize the effective material properties of the JP structures and find that the uni-lamellar structures increases orientation order under shear flow, the multilamellar structure behaves as a shear thinning fluid, and the striated structure possesses a yield stress. These numerical results provide insights into dynamic control of non-equilibrium active biological systems with similar self-organization.

Using a Cauchy integral formulation of the boundary integral equations, we simulate the erosion a porous medium comprised of up to 100 solid bodies embedded in a Stokes flow. The grains of the medium are resolved individually and erode under the action of surface shear stress. Through nonlinear feedback with the surrounding flow fields, microscopic changes in grain morphology give way to larger-scale features in the medium such as channelization. The Cauchy-integral formulation and associated quadrature formulas enable us to resolve dense configurations of nearly contacting bodies. We observe substantial anisotropy to develop over the course of erosion; that is, the configurations that result from erosion generally permit flow in the longitudinal direction more easily than in the transverse direction by up to a factor of six. These results suggest that the erosion of solid material from groundwater flows may contribute to previously observed anisotropy of natural porous media.

We consider the problem of determining the arrival statistics of unbiased planar random walkers to complex target configurations. In contrast to problems posed in finite domains, simple moments of the distribution, such as the mean (MFPT) and variance, are not defined and it is necessary to obtain the full arrival statistics. We describe several methods to obtain these distributions and other associated quantities such as splitting probabilities. One approach combines a Laplace transform of the underlying parabolic equation with matched asymptotic analysis followed by numerical transform inversion. The second approach is similar, but uses a boundary integral equation method to solve for the Laplace transformed variable. To validate the results of this theory, and to obtain the arrival time statistics in very general configurations of absorbers, we introduce an efficient Kinetic Monte Carlo (KMC) method that describes trajectories as a combination of large but exactly solvable projection steps. The effectiveness of these methodologies is demonstrated on a variety of challenging examples highlighting the applicability of these methods to a variety of practical scenarios, such as source inference. A particularly useful finding arising from these results is that homogenization theories, in which complex configurations are replaced by equivalent simple ones, are remarkably effective at describing arrival time statistics.

We develop a new model, to our knowledge, for the many-body hydrodynamics of amphiphilic Janus particles suspended in a viscous background flow. The Janus particles interact through a hydrophobic attraction potential that leads to self-assembly into bilayer structures. We adopt an efficient integral equation method for solving the screened Laplace equation for hydrophobic attraction and for solving the mobility problem for hydrodynamic interactions. The integral equation formulation accurately captures both interactions for near touched boundaries. Under a linear shear flow, we observe the tank-treading deformation in a two-dimensional vesicle made of Janus particles. The results yield measurements of inter-monolayer friction, membrane permeability, and at large shear rates, membrane rupture. The simulations studies include a vesicle in parabolic flow and vesicle-vesicle interactions in shear and extensional flows. The hydrodynamics of the Janus particles vesicle replicate the behaviour of an inextensible elastic vesicle membrane.

In this work, we determine the full expression for the global truncation error of hyperbolic partial differential equations (PDEs). In particular, we use theoretical analysis and symbolic algebra to find exact expressions for the coefficients of the generic global truncation error. Our analysis is valid for any hyperbolic PDE, be it linear or non-linear, and employing finite difference, finite volume, or finite element discretization in space, and advanced in time with a predictor-corrector, multistep, or a deferred correction method, belonging to the Method of Lines. Furthermore, we discuss the practical implications of this analysis. If we employ a stable numerical scheme and the orders of accuracy of the global solution error and the global truncation error agree, we make the following asymptotic observations: (a) the order of convergence at constant ratio of $\Delta t$ to $\Delta x$ is governed by the minimum of the orders of the spatial and temporal discretizations, and (b) convergence cannot even be guaranteed under only spatial or temporal refinement. An implication of (a) is that it is impractical to invest in a time-stepping method of order higher than the spatial discretization. In addition to (b), we demonstrate that under certain circumstances, the error can even monotonically increase with refinement only in space or only in time, and explain why this phenomenon occurs. To verify our theoretical findings, we conduct convergence studies of linear and non-linear advection equations using finite difference and finite volume spatial discretizations, and predictor-corrector and multistep time-stepping methods. Finally, we study the effect of slope limiters and monotonicity-preserving strategies on the order of accuracy.

Lipid bilayer membranes have a native (albeit small) permeability for water molecules. Under an external load, provided that the bilayer structure stays intact and does not suffer from poration or rupture, a lipid membrane deforms and its water influx/efflux is often assumed negligible in the absence of osmolarity. In this work we use boundary integral simulations to investigate the effects of water permeability on the vesicle hydrodynamics due to a mechanical load, such as the viscous stress from an external flow deforming a vesicle membrane in free space or pushing it through a confinement. Incorporating the membrane permeability into the framework of Helfrich free energy for an inextensible, elastic membrane as a model for a semipermeable vesicle, we illustrate that, in the absence of an osmotic stress gradient, the semipermeable vesicle is affected by water influx/efflux over a sufficiently long time or under a strong confinement. Our simulations quantify the conditions for water permeation to be negligible in terms of the time scales, flow strength, and confinement. These results shed light on how microfluidic confinement can be utilized to estimate membrane permeability.

Transport of viscous fluid through porous media is a direct consequence of the pore structure. Here we investigate transport through a specific class of two-dimensional porous geometries, namely those formed by fluid-mechanical erosion. We investigate the tortuosity and dispersion by analyzing the first two statistical moments of tracer trajectories. For most initial configurations, tortuosity decreases in time as a result of erosion increasing the porosity. However, we find that tortuosity can also increase transiently in certain cases. The porosity-tortuosity relationships that result from our simulations are compared with models available in the literature. Asymptotic dispersion rates are also strongly affected by the erosion process, as well as by the number and distribution of the eroding bodies. Finally, we analyze the pore size distribution of an eroding geometry. The simulations are performed by combining a high-fidelity boundary integral equation solver for the fluid equations, a second-order stable time stepping method to simulate erosion, and new numerical methods to stably and accurately resolve nearly-touching eroded bodies and particle trajectories near the eroding bodies.

The dynamics of an adhesive two-dimensional vesicle doublet under various flow conditions is investigated numerically using a high-order, adaptive-in-time boundary integral method. In a quiescent flow, two nearby vesicles move slowly towards each other under the adhesive potential, pushing out fluid between them to form a vesicle doublet at equilibrium. A lubrication analysis on such draining of a thin film gives the dependencies of draining time on adhesion strength and separation distance that are in good agreement with numerical results. In a planar extensional flow we find a stable vesicle doublet forms only when two vesicles collide head-on around the stagnation point. In a microfluid trap where the stagnation of an extensional flow is dynamically placed in the middle of a vesicle doublet through an active control loop, novel dynamics of a vesicle doublet are observed. Numerical simulations show that there exists a critical extensional flow rate above which adhesive interaction is overcome by the diverging stream, thus providing a simple method to measure the adhesion strength between two vesicle membranes. In a planar shear flow, numerical simulations reveal that a vesicle doublet may form provided that the adhesion strength is sufficiently large at a given vesicle reduced area. Once a doublet is formed, its oscillatory dynamics is found to depend on the adhesion strength and their reduced area. Furthermore the effective shear viscosity of a dilute suspension of vesicle doublets is found to be a function of the reduced area. Results from these numerical studies and analysis shed light on the hydrodynamic and rheological consequences of adhesive interactions between vesicles in a viscous fluid.

Dangerous radioactive waste leftover from the Cold War era nuclear weapons production continues to contaminate sixteen sites around the United States. Although many challenges and obstacles exist in decontaminating these sites, two particularly difficult tasks associated with cleanup of this waste are extracting and separating actinide elements from the remainder of the solution, containing other actinide elements and non actinide elements. Developing effective methods for performing these separations is possible by designing new chelating agents that form stable complexes with actinide elements, and by investigating the interactions between the chelating agents and the actinide elements. In this work, new chelating agents (or ligands) with potential to facilitate the separation of radioactive waste are designed for Th, Pa, and U using relativistic Density Functional Theory (DFT) inconjunction with a high-throughput algorithm. We show that both methodologies can be combined efficiently to accelerate discovery and design of new ligands for separation of the radioactive actinides. The main hypothesis that we test with this approach is that the strength of secondary coordination sphere (SCS) can be tuned to increase the selectivity of binding with different actinides. More specifically, we show that links that connect two of the catecholamide ligands via covalent interactions are then added to increase the overall stability of the complex. The effect of increase in the selectivity is also observed when non-covalent interaction is used between ligands. We apply this approach for Th, Pa, and U, and discover linkers that can be used with other ligands., Comment: 22 pages, 13 figures, Supporting information (5 pages)

We consider suspensions of rigid bodies in a two-dimensional viscous fluid. Even with high-fidelity numerical methods, unphysical contact between particles occurs because of spatial and temporal discretization errors. We apply the method of Lu et al. [Journal of Computational Physics, 347:160-182, 2017] where overlap is avoided by imposing a minimum separation distance. In its original form, the method discretizes interactions between different particles explicitly. Therefore, to avoid stiffness, a large minimum separation distance is used. In this paper, we extend the method of Lu et al. by treating all interactions implicitly. This new time stepping method is able to simulate dense suspensions with large time step sizes and a small minimum separation distance. The method is tested on various unbounded and bounded flows, and rheological properties of the resulting suspensions are computed.

We develop numerical methods to simulate the fluid-mechanical erosion of many bodies in two-dimensional Stokes flow. The broad aim is to simulate the erosion of a porous medium (e.g. groundwater flow) with grain-scale resolution. Our fluid solver is based on a second-kind boundary integral formulation of the Stokes equations that is discretized with a spectrally-accurate Nystrom method and solved with fast-multipole-accelerated GMRES. The fluid solver provides the surface shear stress which is used to advance solid boundaries. We regularize interface evolution via curvature penalization using the $\theta$-$L$ formulation, which affords numerically stable treatment of stiff terms and therefore permits large time steps. The overall accuracy of our method is spectral in space and second-order in time. The method is computationally efficient, with the fluid solver requiring $\mathcal{O}(N)$ operations per GMRES iteration, a mesh-independent number of GMRES iterations, and a one-time $\mathcal{O}(N^2)$ computation to compute the shear stress. We benchmark single-body results against analytical predictions for the limiting morphology and vanishing rate. Multibody simulations reveal the spontaneous formation of channels between bodies of close initial proximity. The channelization is associated with a dramatic reduction in the resistance of the porous medium, much more than would be expected from the reduction in grain size alone., Comment: 10 figures, 3 tables, 28 pages

Wildland fire dynamics is a complex turbulent dimensional process. Cellular automata (CA) is an efficient tool to predict fire dynamics, but the main parameters of the method are challenging to estimate. To overcome this challenge, we compute statistical distributions of the key parameters of a CA model using infrared images from controlled burns. Moreover, we apply this analysis to different spatial scales and compare the experimental results to a simple statistical model. By performing this analysis and making this comparison, several capabilities and limitations of CA are revealed., Comment: 14 pages, 7 figures, 2 tables

We consider the bi-Laplacian eigenvalue problem for the modes of vibration of a thin elastic plate with a discrete set of clamped points. A high-order boundary integral equation method is developed for efficient numerical determination of these modes in the presence of multiple localized defects for a wide range of two-dimensional geometries. The defects result in eigenfunctions with a weak singularity that is resolved by decomposing the solution as a superposition of Green's functions plus a smooth regular part. This method is applied to a variety of regular and irregular domains and two key phenomena are observed. First, careful placement of clamping points can entirely eliminate particular eigenvalues and suggests a strategy for manipulating the vibrational characteristics of rigid bodies so that undesirable frequencies are removed. Second, clamping of the plate can result in partitioning of the domain so that vibrational modes are largely confined to certain spatial regions. This numerical method gives a precision tool for tuning the vibrational characteristics of thin elastic plates., Comment: 21 pages, 13 figures

Vesicle suspensions appear in many biological and industrial applications. These suspensions are characterized by rich and complex dynamics of vesicles due to their interaction with the bulk fluid, and their large deformations and nonlinear elastic properties. Many existing state-of-the-art numerical schemes can resolve such complex vesicle flows. However, even when using provably optimal algorithms, these simulations can be computationally expensive, especially for suspensions with a large number of vesicles. In this paper, we investigate the effect of a number of algorithmic empirical fixes in an attempt to make low-resolution simulations more stable and more predictive. Based on our empirical studies for a number of flow configurations, we propose a scheme that attempts to integrate these fixes in a systematic way. This low-resolution scheme is an extension of our previous work Quaife and Biros (2014) and Quaife and Biros (2016). Our low-resolution correction algorithms (LRCA) include anti-aliasing and membrane reparametrization for avoiding spurious oscillations in vesicles' membranes, adaptive time stepping and a repulsion force for handling vesicle collisions and, correction of vesicles' area and arc-length for maintaining physical vesicle shapes. We perform a systematic error analysis by comparing the low-resolution simulations of dilute and dense suspensions with their high-fidelity, fully resolved, counterparts. We observe that the LRCA enables both efficient and statistically accurate low-resolution simulations of vesicle suspensions, while it can be 10x to 100x faster.

We consider an efficient preconditioner for boundary integral equation (BIE) formulations of the two-dimensional Stokes equations in porous media. While BIEs are well-suited for resolving the complex porous geometry, they lead to a dense linear system of equations that is computationally expensive to solve for large problems. This expense is further amplified when a significant number of iterations is required in an iterative Krylov solver such as GMRES. In this paper, we apply a fast inexact direct solver, the inverse fast multipole method (IFMM), as an efficient preconditioner for GMRES. This solver is based on the framework of $\mathcal{H}^{2}$-matrices and uses low-rank compressions to approximate certain matrix blocks. It has a tunable accuracy $\varepsilon$ and a computational cost that scales as $\mathcal{O} (N \log^2 1/\varepsilon)$. We discuss various numerical benchmarks that validate the accuracy and confirm the efficiency of the proposed method. We demonstrate with several types of boundary conditions that the preconditioner is capable of significantly accelerating the convergence of GMRES when compared to a simple block-diagonal preconditioner, especially for pipe flow problems involving many pores., Comment: 23 pages, 15 figures, 9 tables

We study mixing in Stokesian vesicle suspensions in two dimensions on a cylindrical Couette apparatus using numerical simulations. The vesicle flow simulation is done using a boundary integral method and the advection-diffusion equation for the mixing of the solute is solved using a pseudo-spectral scheme. We study the effect of the area fraction, the viscosity contrast between the inside (the vesicles) and the outside (the bulk) fluid, the initial condition of the solute, and the mixing metric. We compare mixing in the suspension with mixing in the Couette apparatus without vesicles. On the one hand, the presence of vesicles in most cases, slightly suppresses mixing. This is because the solute can be only diffused across the vesicle interface and not advected. On the other hand, there exist spatial distributions of the solute for which the unperturbed Couette flow completely fails to mix whereas the presence of vesicles enables mixing. We derive a simple condition that relates the velocity and solute and can be used to characterize the cases in which the presence of vesicles promotes mixing., Comment: 25 pages, 14 figures

We construct a high-order adaptive time stepping scheme for vesicle suspensions with viscosity contrast. The high-order accuracy is achieved using a spectral deferred correction (SDC) method, and adaptivity is achieved by estimating the local truncation error with the numerical error of physically constant values. Numerical examples demonstrate that our method can handle suspensions with vesicles that are tumbling, tank-treading, or both. Moreover, we demonstrate that a user-prescribed tolerance can be automatically achieved for simulations with long time horizons.

An integral equation method for solving the Yukawa-Beltrami equation on a multiply-connected sub-manifold of the unit sphere is presented. A fundamental solution for the Yukawa-Beltrami operator is constructed. This fundamental solution can be represented by conical functions. Using a suitable representation formula, a Fredholm equation of the second kind with a compact integral operator needs to be solved. The discretization of this integral equation leads to a linear system whose condition number is bounded independent of the size of the system. Several numerical examples exploring the properties of this integral equation are presented.

In this paper, a new two-parameter family of regularized kernels is introduced, suitable for applying high-order time stepping to N-body systems. These high-order kernels are derived by truncating a Taylor expansion of the non-regularized kernel about $(r^2+\epsilon^2)$, generating a sequence of increasingly more accurate kernels. This paper proves the validity of this two-parameter family of regularized kernels, constructs error estimates, and illustrates the benefits of using high-order kernels through numerical experiments., Comment: 27 pages, 15 figures

We present an adaptive arbitrary-order accurate time-stepping numerical scheme for the flow of vesicles suspended in Stokesian fluids. Our scheme can be summarized as an approximate implicit spectral deferred correction (SDC) method. Applying a textbook fully implicit SDC scheme to vesicle flows is prohibitively expensive. For this reason we introduce several approximations. Our scheme is based on a semi-implicit linearized low-order time stepping method. (Our discretization is spectrally accurate in space.) We also use invariant properties of vesicle flows, constant area and boundary length in two dimensions, to reduce the computational cost of error estimation for adaptive time stepping. We present results in two dimensions for single-vesicle flows, constricted geometry flows, converging flows, and flows in a Couette apparatus. We experimentally demonstrate that the proposed scheme enables automatic selection of the step size and high-order accuracy.

We consider numerical algorithms for the simulation of the rheology of two-dimensional vesicles suspended in a viscous Stokesian fluid. The vesicle evolution dynamics is governed by hydrodynamic and elastic forces. The elastic forces are due to local inextensibility of the vesicle membrane and resistance to bending. Numerically resolving vesicle flows poses several challenges. For example, we need to resolve moving interfaces, address stiffness due to bending, enforce the inextensibility constraint, and efficiently compute the (non-negligible) long-range hydrodynamic interactions. Our method is based on the work of {\em Rahimian, Veerapaneni, and Biros, "Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method", Journal of Computational Physics, 229 (18), 2010}. It is a boundary integral formulation of the Stokes equations coupled to the interface mass continuity and force balance. We extend the algorithms presented in that paper to increase the robustness of the method and enable simulations with concentrated suspensions. In particular, we propose a scheme in which both intra-vesicle and inter-vesicle interactions are treated semi-implicitly. In addition we use special integration for near-singular integrals and we introduce a spectrally accurate collision detection scheme. We test the proposed methodologies on both unconfined and confined flows for vesicles whose internal fluid may have a viscosity contrast with the bulk medium. Our experiments demonstrate the importance of treating both intra-vesicle and inter-vesicle interactions accurately.

The discretization of the double-layer potential integral equation for the interior Dirichlet Laplace problem in a domain with smooth boundary results in a linear system that has a bounded condition number. Thus, the number of iterations required for the convergence of a Krylov method is, asymptotically, independent of the discretization size $N$. Using the Fast Multipole Method (FMM) to accelerate the matrix-vector products, we obtain an optimal $\mathcal{O}(N)$ solver. In practice, however, when the geometry is complicated, the number of Krylov iterations can be quite large---to the extend that necessitates the use of preconditioning. We summarize the different methodologies that have appeared in the literature (single-grid, multigrid, approximate sparse inverses) and we propose a new class of preconditioners based on an FMM-based spatial decomposition of the double-layer operator. We present an experimental study in which we compare the different approaches and we discuss the merits and shortcomings of our approach. Our method can be easily extended to other second-kind integral equations with non-oscillatory kernels in two and three dimensions.

We present an efficient integral equation approach to solve the heat equation, $u_t (\x) - \Delta u(\x) = F(\x,t)$, in a two-dimensional, multiply connected domain, and with Dirichlet boundary conditions. Instead of using integral equations based on the heat kernel, we take the approach of discretizing in time, first. This leads to a non-homogeneous modified Helmholtz equation that is solved at each time step. The solution to this equation is formulated as a volume potential plus a double layer potential.The volume potential is evaluated using a fast multipole-accelerated solver. The boundary conditions are then satisfied by solving an integral equation for the homogeneous modified Helmholtz equation. The integral equation solver is also accelerated by the fast multipole method (FMM). For a total of $N$ points in the discretization of the boundary and the domain, the total computational cost per time step is $O(N)$ or $O(N\log N)$.

We present a collection of integral equation methods for the solution to the two-dimensional, modified Helmholtz equation, $u(\x) - \alpha^2 \Delta u(\x) = 0$, in bounded or unbounded multiply-connected domains. We consider both Dirichlet and Neumann problems. We derive well-conditioned Fredholm integral equations of the second kind, which are discretized using high-order, hybrid Gauss-trapezoid rules. Our fast multipole-based iterative solution procedure requires only O(N) or $O(N\log N)$ operations, where N is the number of nodes in the discretization of the boundary. We demonstrate the performance of the methods on several numerical examples., Comment: Published in Computers & Mathematics with Applications