Exponential Krylov time integration for modeling multi-frequency optical response with monochromatic sources

Light incident on a layer of scattering material such as a piece of sugar or white paper forms a characteristic speckle pattern in transmission and reflection. The information hidden in the correlations of the speckle pattern with varying frequency, polarization and angle of the incident light can be exploited for applications such as biomedical imaging and high-resolution microscopy. Conventional computational models for multi-frequency optical response involve multiple solution runs of Maxwell's equations with monochromatic sources. Exponential Krylov subspace time solvers are promising candidates for improving efficiency of such models, as single monochromatic solution can be reused for the other frequencies without performing full time-domain computations at each frequency. However, we show that the straightforward implementation appears to have serious limitations. We further propose alternative ways for efficient solution through Krylov subspace methods. Our methods are based on two different splittings of the unknown solution into different parts, each of which can be computed efficiently. Experiments demonstrate a significant gain in computation time with respect to the standard solvers.Comment: 22 pages, 4 figure

A Fast Algorithm for Parabolic PDE-based Inverse Problems Based on Laplace Transforms and Flexible Krylov Solvers

We consider the problem of estimating parameters in large-scale weakly nonlinear inverse problems for which the underlying governing equations is a linear, time-dependent, parabolic partial differential equation. A major challenge in solving these inverse problems using Newton-type methods is the computational cost associated with solving the forward problem and with repeated construction of the Jacobian, which represents the sensitivity of the measurements to the unknown parameters. Forming the Jacobian can be prohibitively expensive because it requires repeated solutions of the forward and adjoint time-dependent parabolic partial differential equations corresponding to multiple sources and receivers. We propose an efficient method based on a Laplace transform-based exponential time integrator combined with a flexible Krylov subspace approach to solve the resulting shifted systems of equations efficiently. Our proposed solver speeds up the computation of the forward and adjoint problems, thus yielding significant speedup in total inversion time. We consider an application from Transient Hydraulic Tomography (THT), which is an imaging technique to estimate hydraulic parameters related to the subsurface from pressure measurements obtained by a series of pumping tests. The algorithms discussed are applied to a synthetic example taken from THT to demonstrate the resulting computational gains of this proposed method

MATEX: A Distributed Framework for Transient Simulation of Power Distribution Networks

We proposed MATEX, a distributed framework for transient simulation of power distribution networks (PDNs). MATEX utilizes matrix exponential kernel with Krylov subspace approximations to solve differential equations of linear circuit. First, the whole simulation task is divided into subtasks based on decompositions of current sources, in order to reduce the computational overheads. Then these subtasks are distributed to different computing nodes and processed in parallel. Within each node, after the matrix factorization at the beginning of simulation, the adaptive time stepping solver is performed without extra matrix re-factorizations. MATEX overcomes the stiff-ness hinder of previous matrix exponential-based circuit simulator by rational Krylov subspace method, which leads to larger step sizes with smaller dimensions of Krylov subspace bases and highly accelerates the whole computation. MATEX outperforms both traditional fixed and adaptive time stepping methods, e.g., achieving around 13X over the trapezoidal framework with fixed time step for the IBM power grid benchmarks.Comment: ACM/IEEE DAC 2014. arXiv admin note: substantial text overlap with arXiv:1505.0669

Explicit schemes for time propagating many-body wavefunctions

Accurate theoretical data on many time-dependent processes in atomic and molecular physics and in chemistry require the direct numerical solution of the time-dependent Schr\"odinger equation, thereby motivating the development of very efficient time propagators. These usually involve the solution of very large systems of first order differential equations that are characterized by a high degree of stiffness. We analyze and compare the performance of the explicit one-step algorithms of Fatunla and Arnoldi. Both algorithms have exactly the same stability function, therefore sharing the same stability properties that turn out to be optimum. Their respective accuracy however differs significantly and depends on the physical situation involved. In order to test this accuracy, we use a predictor-corrector scheme in which the predictor is either Fatunla's or Arnoldi's algorithm and the corrector, a fully implicit four-stage Radau IIA method of order 7. We consider two physical processes. The first one is the ionization of an atomic system by a short and intense electromagnetic pulse; the atomic systems include a one-dimensional Gaussian model potential as well as atomic hydrogen and helium, both in full dimensionality. The second process is the decoherence of two-electron quantum states when a time independent perturbation is applied to a planar two-electron quantum dot where both electrons are confined in an anharmonic potential. Even though the Hamiltonian of this system is time independent the corresponding differential equation shows a striking stiffness. For the one-dimensional Gaussian potential we discuss in detail the possibility of monitoring the time step for both explicit algorithms. In the other physical situations that are much more demanding in term of computations, we show that the accuracy of both algorithms depends strongly on the degree of stiffness of the problem.Comment: 24 pages, 14 Figure

An Arnoldi-frontal approach for the stability analysis of flows in a collapsible channel

In this paper, we present a new approach based on a combination of the Arnoldi and frontal methods for solving large sparse asymmetric and generalized complex eigenvalue problems. The new eigensolver seeks the most unstable eigensolution in the Krylov subspace and makes use of the efficiency of the frontal solver developed for the finite element methods. The approach is used for a stability analysis of flows in a collapsible channel and is found to significantly improve the computational efficiency compared to the traditionally used QZ solver or a standard Arnoldi method. With the new approach, we are able to validate the previous results obtained either on a much coarser mesh or estimated from unsteady simulations. New neutral stability solutions of the system have been obtained which are beyond the limits of previously used methods

Wave packet propagation by the Faber polynomial approximation in electrodynamics of passive media

Maxwell's equations for propagation of electromagnetic waves in dispersive and absorptive (passive) media are represented in the form of the Schr\"odinger equation $i\partial \Psi/\partial t = {H}\Psi$, where ${H}$ is a linear differential operator (Hamiltonian) acting on a multi-dimensional vector $\Psi$ composed of the electromagnetic fields and auxiliary matter fields describing the medium response. In this representation, the initial value problem is solved by applying the fundamental solution $\exp(-itH)$ to the initial field configuration. The Faber polynomial approximation of the fundamental solution is used to develop a numerical algorithm for propagation of broad band wave packets in passive media. The action of the Hamiltonian on the wave function $\Psi$ is approximated by the Fourier grid pseudospectral method. The algorithm is global in time, meaning that the entire propagation can be carried out in just a few time steps. A typical time step is much larger than that in finite differencing schemes, $\Delta t_F \gg \|H\|^{-1}$. The accuracy and stability of the algorithm is analyzed. The Faber propagation method is compared with the Lanczos-Arnoldi propagation method with an example of scattering of broad band laser pulses on a periodic grating made of a dielectric whose dispersive properties are described by the Rocard-Powels-Debye model. The Faber algorithm is shown to be more efficient. The Courant limit for time stepping, $\Delta t_C \sim \|H\|^{-1}$, is exceeded at least in 3000 times in the Faber propagation scheme.Comment: Latex, 17 pages, 4 figures (separate png files); to appear in J. Comput. Phy