Approximate Solver Research Articles

Modified ghost fluid method for three‐dimensional compressible multimaterial flows with interfaces exhibiting large curvature and topological change

SummaryIn order to capture the material interface dynamics, especially under the impact of strong shocks, the key feature of the modified ghost fluid method (MGFM) is to construct a multimaterial Riemann problem normal to the interface and use its solution to define interface conditions. However, such process sometimes may not be easily or accurately implemented when the multidimensional interfaces come into contact or undergo significant deformations. In this article, a three‐dimensional interface treating procedure is developed for a wide range of compressible multimaterial flows. It utilizes the MGFM with an explicit approximate Riemann problem solver to define interface conditions. More importantly, a weighted average technique is designed to optimize the treatment for interfaces exhibiting large curvature and topological change. This remedies two defects of the traditional approach in these extreme cases. One is that the normal directions of interfacial ghost nodes may not be easily calculated. The other is that the interface conditions may not be accurately defined. The numerical methodology is validated through several typical problems, including gas‐liquid Riemann problem and shock‐bubble/droplet interaction. These results indicate that the developed method is capable of dealing with interfacial evolutions in three dimensions, especially when interfaces undergo merger, fragmentation, and other complex changes.

A scalable computational platform for particulate Stokes suspensions

We describe a computational framework for simulating suspensions of rigid particles in Newtonian Stokes flow. One central building block is a collision-resolution algorithm that overcomes the numerical constraints arising from particle collisions. This algorithm extends the well-known complementarity method for non-smooth multi-body dynamics to resolve collisions in dense rigid body suspensions. This approach formulates the collision resolution problem as a linear complementarity problem with geometric ‘non-overlapping’ constraints imposed at each time-step. It is then reformulated as a constrained quadratic programming problem and the Barzilai-Borwein projected gradient descent method is applied for its solution. This framework is designed to be applicable for any convex particle shape, e.g., spheres and spherocylinders, and applicable to any Stokes mobility solver, including the Rotne-Prager-Yamakawa approximation, Stokesian Dynamics, and PDE solvers (e.g., boundary integral and immersed boundary methods). In particular, this method imposes Newton's Third Law and records the entire contact network. Further, we describe a fast, parallel, and spectrally-accurate boundary integral method tailored for spherical particles, capable of resolving lubrication effects. We show weak and strong parallel scalings up to 8×104 particles with approximately 4×107 degrees of freedom on 1792 cores. We demonstrate the versatility of this framework with several examples, including sedimentation of particle clusters, and active matter systems composed of ensembles of particles driven to rotate.

Supporting hard queries over probabilistic preferences

Preference analysis is widely applied in various domains such as social choice and e-commerce. A recently proposed frame- work augments the relational database with a preference re- lation that represents uncertain preferences in the form of statistical ranking models, and provides methods to evaluate Conjunctive Queries (CQs) that express preferences among item attributes. In this paper, we explore the evaluation of queries that are more general and harder to compute. The main focus of this paper is on a class of CQs that cannot be evaluated by previous work. These queries are provably hard since relate variables that represent items be- ing compared. To overcome this hardness, we instantiate these variables with their domain values, rewrite hard CQs as unions of such instantiated queries, and develop several exact and approximate solvers to evaluate these unions of queries. We demonstrate that exact solvers that target specific common kinds of queries are far more efficient than gen- eral solvers. Further, we demonstrate that sophisticated ap- proximate solvers making use of importance sampling can be orders of magnitude more efficient than exact solvers, while showing good accuracy. In addition to supporting provably hard CQs, we also present methods to evaluate an important family of count queries, and of top-k queries.

CASSCF with Extremely Large Active Spaces Using the Adaptive Sampling Configuration Interaction Method.

The complete active space self-consistent field (CASSCF) method is the principal approach employed for studying strongly correlated systems. However, exact CASSCF can only be performed on small active spaces of ∼20 electrons in ∼20 orbitals due to exponential growth in the computational cost. We show that employing the Adaptive Sampling Configuration Interaction (ASCI) method as an approximate Full CI solver in the active space allows CASSCF-like calculations within chemical accuracy (<1 kcal/mol for relative energies) in active spaces with more than ∼50 active electrons in ∼50 active orbitals, significantly increasing the sizes of systems amenable to accurate multiconfigurational treatment. The main challenge with using any selected CI-based approximate CASSCF is the orbital optimization problem; they tend to exhibit large numbers of local minima in orbital space due to their lack of invariance to active-active rotations (in addition to the local minima that exist in exact CASSCF). We highlight methods that can avoid spurious local extrema as a practical solution to the orbital optimization problem. We employ ASCI-SCF to demonstrate a lack of polyradical character in moderately sized periacenes with up to 52 correlated electrons and compare against heat-bath CI on an iron porphyrin system with more than 40 correlated electrons.

A Robust Riemann Solver for Multiple Hydro-Elastoplastic Solid Mediums

We propose a robust approximate solver for the hydro-elastoplastic solid material, a general constitutive law extensively applied in explosion and high speed impact dynamics, and provide a natural transformation between the fluid and solid in the case of phase transitions. The hydrostatic components of the solid is described by a family of general Mie-Gr\"uneisen equation of state (EOS), while the deviatoric component includes the elastic phase, linearly hardened plastic phase and fluid phase. The approximate solver provides the interface stress and normal velocity by an iterative method. The well-posedness and convergence of our solver are proved with mild assumptions on the equations of state. The proposed solver is applied in computing the numerical flux at the phase interface for our compressible multi-medium flow simulation on Eulerian girds. Several numerical examples, including Riemann problems, shock-bubble interactions, implosions and high speed impact applications, are presented to validate the approximate solver.

Fast Depth Estimation for Light Field Cameras

Fast depth estimation for light field images is an important task for multiple applications such as image-based rendering and refocusing. Most previous approaches to light field depth estimation involve high computational costs. Therefore, in this study, we propose a fast depth estimation method based on multi-view stereo matching for light field images. Similar to other conventional methods, our method consists of initial depth estimation and refinement. For the initial estimation, we use a one-bit feature for each pixel and calculate matching costs by summing all combinations of viewpoints with a fast algorithm. To reduce computational time, we introduce an offline viewpoint selection strategy and cost volume interpolation. Our refinement process solves the minimization problem in which the objective function consists of ℓ <sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">1</sub> data and smoothness terms. Although this problem can be solved via a graph cuts algorithm, it is computationally expensive; therefore, we propose an approximate solver based on a fast-weighted median filter. Experiments on synthetic and real-world data show that our method achieves competitive accuracy with the shortest computational time of all methods.

An HLLC-type Riemann solver and high-resolution Godunov method for a two-phase model of reactive flow with general equations of state

A high-resolution Godunov method is developed for a two-phase model of reactive flow. The model considers general equations of state of Mie-Grüneisen form and different choices for the reaction rate, both relevant for simulations of the dynamical behavior of detonations in PBX-type granular explosives. The numerical approach employs a Riemann solver, and various options ranging from an exact solver to an approximate solver of HLLC type are described. A principal aim of the paper is an assessment of the approximate Riemann solvers in terms of the accuracy and efficiency of numerical solutions of the two-phase model. In particular, this study considers the state-sensitive behavior in different regimes of the solutions in a reverse-impact configuration, including compaction, run to detonation, and the evolution to compaction-led and reaction-led steady detonations. For these regimes, the convergence properties of the numerical approach are examined for the different choices of the Riemann solver, as is the error in the solutions at lower grid resolution. Finally, the computational performance of the time-stepping scheme is assessed for the different solvers.

On the use of the envelope model for down-ramp injection in laser-plasma accelerators

The path towards high 6D-brightness electron beams coming from plasma accelerators relies on numerical simulations to investigate and explain novel physics and parameter regimes. Particle-in-cell codes are a reliable tool for this purpose, but full-scale simulations can be computationally challenging, specifically when performing large parameter scans. Approximate models and solvers are an alternative to quickly assess new regimes, but one needs to understand the validity of these models for each specific problem. In this work, we investigate the usefulness of the ponderomotive guiding center solver to model density down-ramp injection in laser wakefield accelerators. We also report results for density down-ramp injection in the context of the European Plasma Research Accelerator with eXcellence in Applications, where we find that this method can produce , beams that fulfill all the requirements for the injector stage, with <1% energy spread and <100 nm emittance.

One simulation to have them all: performance of the Bias Assignment Method against N-body simulations

Abstract In this paper we demonstrate that the information encoded in one single (sufficiently large) N-body simulation can be used to reproduce arbitrary numbers of halo catalogues, using approximated realisations of dark matter density fields with different initial conditions. To this end we use as a reference one realisation (from an ensemble of 300) of the Minerva N-body simulations and the recently published Bias Assignment Method to extract the local and non-local bias linking the halo to the dark matter distribution. We use an approximate (and fast) gravity solver to generate 300 dark matter density fields from the down-sampled initial conditions of the reference simulation and sample each of these fields using the halo-bias and a kernel, both calibrated from the arbitrarily chosen realisation of the reference simulation. We show that the power spectrum, its variance and the three-point statistics are reproduced within $\sim 2\%$ (up to k ∼ 1.0 h Mpc−1), $\sim 5-10\%$ and $\sim 10\%$, respectively. Using a model for the real space power spectrum (with three free bias parameters), we show that the covariance matrices obtained from our procedure lead to parameter uncertainties that are compatible within $\sim 10\%$ with respect to those derived from the reference covariance matrix, and motivate approaches that can help to reduce these differences to $\sim 1\%$. Our method has the potential to learn from one simulation with moderate volumes and high-mass resolution and extrapolate the information of the bias and the kernel to larger volumes, making it ideal for the construction of mock catalogues for present and forthcoming observational campaigns such as Euclid or DESI.

Parallelization of the inverse fast multipole method with an application to boundary element method

We present an algorithm to parallelize the inverse fast multipole method (IFMM), which is an approximate direct solver for dense linear systems. The parallel scheme is based on a greedy coloring algorithm, where two nodes in the hierarchy with the same color are separated by at least σ nodes. We proved that when σ≥6, the workload associated with one color is embarrassingly parallel. However, the number of nodes in a group (color) may be small when σ=6. Therefore, we also explored σ=3, where a small fraction of the algorithm needs to be serialized, and the overall parallel efficiency was improved. We implemented the parallel IFMM using OpenMP for shared-memory machines. Successively, we applied it to a fast-multipole accelerated boundary element method (FMBEM) as a preconditioner, and compared its efficiency with (a) the original IFMM parallelized by linking a multi-threaded linear algebra library and (b) the commonly used parallel block-diagonal preconditioner. Our results showed that our parallel IFMM achieved at most 4× and 11× speedups over the reference method (a) and (b), respectively, in realistic examples involving more than one million variables.

An Efficient MLACA-SVD Solver for Superconducting Integrated Circuit Analysis

Inductance extraction for superconducting integrated circuits requires the accurate solution of structural current distributions. FastHenry is a well known, magnetoquasistatic solver suitable for this task. It is based upon the partial element equivalent circuit, integral equation method, with the structure discretized into hexahedral filaments. It employs the multilevel fast multipole algorithm (MLFMA) for compressed storage of the mutual inductance matrix, which accounts for most of the required memory. This MLFMA implementation is especially memory efficient, given certain approximations and algorithmic parameter choices. However, errors are introduced into the matrix representation. Here, a multilevel adaptive cross approximation solver with singular value decomposition recompression (MLACA-SVD) is presented as an alternative to FastHenry's existing MLFMA solver. MLACA-SVD compresses off-diagonal matrix blocks to a specified error tolerance, based on evaluating selected entries. Quadrature recipes are presented for guaranteed accuracy of matrix entry evaluation. Numerical results for examples of practical interest show that the MLACA-SVD memory scaling versus b (number of filaments) is practically identical to that of FastHenry's MLFMA, and is close to O(blog b). The MLACA-SVD requires less memory for the same solution accuracy, and furthermore offers complete control over matrix approximation errors. For the examples considered, it is found to be a more efficient solver. It is well suited to parallelization.

Split-step Padé solver for three-dimensional Cartesian acoustic parabolic equation in stair-step representation of ocean environment.

In this study, the three-dimensional Cartesian acoustic parabolic equation is fully solved in the stair-step representation of an ocean environment based on the split-step Padé algorithm. The sum representation of the Padé approximation of a full exponential operator in the split-step marching solution is used for parallel computing. A stabilized self-starter solution is derived in the same way. This solver is benchmarked against reference solutions for a rectangular duct propagation, an ideal wedge problem, and an axisymmetric ocean. Additionally, a solution for the Gaussian canyon problem is provided and compared with that obtained from the approximate parabolic equation solver, considering a leading-order cross-term.

MLACA With Modified Grouping Strategy for Efficient Superconducting Circuit Analysis

Electromagnetic analysis of superconducting integrated circuits is routinely required for inductance extraction. FastHenry is a well-known, applicable magnetoquasistatic analysis tool. It is accelerated with the multilevel fast multipole algorithm (MLFMA). A multilevel adaptive cross approximation solver with singular value decomposition recompression (MLACA-SVD) was recently proposed for FastHenry, as an alternative to the MLFMA. It was found to be more efficient, with regards to accuracy versus required memory. It uses the same octree grouping as FastHenry's MLFMA. Here, two modified grouping strategies are proposed to further improve MLACA-SVD efficiency by compressing interactions between larger groups, while maintaining scaling performance consistent with valid admissibility of interactions. The strategies are denoted “shell merging” and “wall merging.” Results show that the MLACA-SVD solver with merging requires about four times less memory than FastHenry's MLFMA, for similar accuracy. Shell merging is found to be slightly superior.

Exact and Approximate Weighted Model Integration with Probability Density Functions Using Knowledge Compilation

Weighted model counting has recently been extended to weighted model integration, which can be used to solve hybrid probabilistic reasoning problems. Such problems involve both discrete and continuous probability distributions. We show how standard knowledge compilation techniques (to SDDs and d-DNNFs) apply to weighted model integration, and use it in two novel solvers, one exact and one approximate solver. Furthermore, we extend the class of employable weight functions to actual probability density functions instead of mere polynomial weight functions.

A Discrete Dipole Approximation Solver Based on the COCG-FFT Algorithm and Its Application to Microwave Breast Imaging.

We introduce the discrete dipole approximation (DDA) for efficiently calculating the two-dimensional electric field distribution for our microwave tomographic breast imaging system. For iterative inverse problems such as microwave tomography, the forward field computation is the time limiting step. In this paper, the two-dimensional algorithm is derived and formulated such that the iterative conjugate orthogonal conjugate gradient (COCG) method can be used for efficiently solving the forward problem. We have also optimized the matrix-vector multiplication step by formulating the problem such that the nondiagonal portion of the matrix used to compute the dipole moments is block-Toeplitz. The computation costs for multiplying the block matrices times a vector can be dramatically accelerated by expanding each Toeplitz matrix to a circulant matrix for which the convolution theorem is applied for fast computation utilizing the fast Fourier transform (FFT). The results demonstrate that this formulation is accurate and efficient. In this work, the computation times for the direct solvers, the iterative solver (COCG), and the iterative solver using the fast Fourier transform (COCG-FFT) are compared with the best performance achieved using the iterative solver (COCG-FFT) in C++. Utilizing this formulation provides a computationally efficient building block for developing a low cost and fast breast imaging system to serve under-resourced populations.

A Bayesian approach to improving the Born approximation for inverse scattering with high-contrast materials††This paper is dedicated to the memory of Armin Lechleiter, a gifted mathematician and friend.

Time harmonic inverse scattering using accurate forward models is often computationally expensive. On the other hand, the use of computationally efficient solvers, such as the Born approximation, may fail if the targets do not satisfy the assumptions of the simplified model. In the Bayesian framework for inverse problems, one can construct a statistical model for the errors that are induced when approximate solvers are used, and hence increase the domain of applicability of the approximate model. In this paper, we investigate the error structure that is induced by the Born approximation and show that the Bayesian approximation error approach can be used to partially recover from these errors. In particular, we study the model problem of reconstruction of the index of refraction of a penetrable medium from measurements of the far field pattern of the scattered wave.

Accelerating multiscale simulation of complex geomodels by use of dynamically adapted basis functions

A number of different multiscale methods have been developed as a robust alternative to upscaling and as a means for accelerated reservoir simulation of high-resolution geomodels. In their basic setup, multiscale methods use a restriction operator to construct a reduced system of flow equations on a coarser grid, and a prolongation operator to map pressure unknowns from the coarse grid back to the original simulation grid. The prolongation operator consists of basis functions computed numerically by solving localized flow problems. One can use the resulting multiscale solver both as a CPR-preconditioner in fully implicit simulators or as an efficient approximate iterative linear solver in a sequential setting. The latter approach has been successful implemented in a commercial simulator. Recently, we have shown that you can obtain significantly faster convergence if you instead of using a single pair of prolongation-restriction operators apply a sequence of such operators, where some of the operators adapt to faults, fractures, facies, or other geobodies. Herein, we present how you can accelerate the convergence even further, if you also include additional basis functions that capture local changes in the pressure.

MOSES: A Streaming Algorithm for Linear Dimensionality Reduction.

This paper introduces Memory-limited Online Subspace Estimation Scheme (MOSES) for both estimating the principal components of streaming data and reducing its dimension. More specifically, in various applications such as sensor networks, the data vectors are presented sequentially to a user who has limited storage and processing time available. Applied to such problems, MOSES can provide a running estimate of leading principal components of the data that has arrived so far and also reduce its dimension. MOSES generalises the popular incremental Singular Vale Decomposition (iSVD) to handle thin blocks of data, rather than just vectors. This minor generalisation in part allows us to complement MOSES with a comprehensive statistical analysis, thus providing the first theoretically-sound variant of iSVD, which has been lacking despite the empirical success of this method. This generalisation also enables us to concretely interpret MOSES as an approximate solver for the underlying non-convex optimisation program. We find that MOSES consistently surpasses the state of the art in our numerical experiments with both synthetic and real-world datasets, while being computationally inexpensive.

Formulation of exactly balanced solvers for blood flow in elastic vessels and their application to collapsed states

In this work, numerical solvers based on extensions of the Roe and HLL schemes are adapted to deal with test cases involving extreme collapsing conditions in elastic vessels. To achieve this goal, the system is transformed to provide a conservation–law form, allowing to define Rankine–Hugoniot conditions. The approximate solvers allow to describe the inner states of the solution. Therefore, source term fixes can be used to prevent unphysical values of vessel area and, at the same time, the eigenvalues of the system control stability. Numerical solvers of different order are tested using a wide variety of Riemann problems, including extreme vessel collapse and blockage. In all cases, the robustness of the approximate solvers presented here is checked using first and third order methods in time and space, using the WENO reconstruction scheme in combination with the TVDRK3 method.

A robust adaptive iterative ensemble smoother scheme for practical history matching applications

Much of the recent work on history matching reservoir models has focused on the iterative Ensemble Smoother (iES) method. This is well suited for practical applications because it avoids frequent simulation restarts which can be expensive for traditional Ensemble Kalman Filter (EnKF) approaches. In this paper, we derive a novel, adaptive iES method (ES-LM) based on the Levenberg-Marquardt algorithm for nonlinear least squares optimization. We demonstrate that the solution of the linearized least squares subproblems from each iteration of ES-LM have a similar structure to those of the standard ensemble smoother update equation. This allows us to use the ensemble smoother as an approximate linear least squares solver and thus avoid expensive adjoint calculations. The proposed algorithm may therefore be seen as a variant of iES whose regularization parameter can be updated using the usual trust region method. The resulting parameter update equation is similar to those in existing iES methods, but our derivation provides new insights into iES that enable us to extend the adaptive ES-LM to two new formulations: (1) Auto ES-LM, where the measurement error has unknown variance that will be estimated by the algorithm; and (2) Robust ES-LM, where the measurement error follows a Cauchy distribution. Robust ES-LM, in particular, reduces the effect of outliers in the data, which can be critical for real-world applications. Although estimating measurement error and adding robustness are common for data assimilation with EnKF, they have not been incorporated into iES within the reservoir engineering community. The three resulting algorithms are easy for reservoir engineers to use for practical history matching applications without the need for deep understanding of the computational details. Finally, we demonstrate the effectiveness and robustness of these algorithms on two synthetic reservoir models.