Schur Complement Method Research Articles

An explicit substructuring method for overlapping domain decomposition based on stochastic calculus

In a recent paper [7], a hybrid supercomputing algorithm for elliptic equations has been proposed. The idea is that the interfacial nodal solutions solve a linear system, whose coefficients are expectations of functionals of stochastic differential equations confined within patches of about subdomain size. Compared to standard substructuring techniques, such as the Schur complement method for the skeleton, the hybrid approach produces an explicit and sparse shrunken matrix—hence suitable for substructuring again. The ultimate goal is to push strong scalability beyond the state of the art by leveraging the potential for parallelisation of stochastic calculus. Here, we present a major revamping of that framework, based on the insight of embedding the domain in a cover of overlapping circles (in two dimensions). This allows for efficient Fourier interpolation along the interfaces (now circumferences) and—crucially—for the evaluation of most of the interfacial system entries as the solution of small boundary value problems on a circle. This is both extremely efficient (as they can be solved in parallel and by the pseudospectral method) and free of Monte Carlo error. Stochastic numerics are only needed on the relatively few circles intersecting the domain boundary. In sum, the new formulation is significantly faster, simpler, and more accurate while retaining all of the advantageous properties of PDDSparse. Numerical experiments are included for the purpose of illustration.

Exploiting NOMA and RIS in Integrated Sensing and Communication

A novel integrated sensing and communication (ISAC) system is proposed, where a dual-functional base station is utilized to transmit the superimposed non-orthogonal multiple access (NOMA) communication signal for serving communication users and sensing targets simultaneously. Furthermore, a new reconfigurable intelligent surface (RIS)-aided-sensing structure, where a dedicated RIS is deployed to provide virtual line-of-sight (LoS) links for radar targets, is also proposed to address the significant path loss or blockage of LoS links for the sensing task. Based on this setup, the beampattern gain at the RIS for the radar target is derived and adopted as a sensing metric. The objective of this paper is to maximize the minimum beampattern gain by jointly optimizing active beamforming at the base station (BS), power allocation coefficients among NOMA users and passive beamforming at the RIS. To tackle the non-convexity of the formulated optimization problem, the beampattern gain and constraints are first transformed into more tractable forms. Then, an iterative block coordinate descent (IBCD) algorithm is proposed by employing successive convex approximation (SCA), Schur complement, semidefinite relaxation (SDR) and sequential rank-one constraint relaxation (SRCR) methods. To reduce the complexity of the proposed IBCD algorithm, a low-complexity iterative alternating optimization (IAO) algorithm is proposed. Particularly, the active beamforming is optimized by solving a semidefinite programming (SDP) problem and the closed-form solutions of the power allocation coefficients are derived. Numerical results show that: i) the proposed RIS-NOMA-ISAC system always outperforms the RIS-ISAC system without NOMA in beampattern gain and illumination power; ii) the low-complexity IAO algorithm achieves a comparable performance to that achieved by the IBCD algorithm. iii) high beampattern gain can be achieved by the proposed joint optimization algorithms in underloaded and overloaded communication scenarios.

Efficient numerical solution of the Fokker–Planck equation using physics-conforming finite element methods

Abstract We consider the Fokker–Planck equation (FPE) for the orientation probability density of fiber suspensions. Using the continuous Galerkin method, we express the numerical solution in terms of Lagrange basis functions that are associated with N nodes of a computational mesh for a domain in the 3D physical space and M nodes of a mesh for the surface of a unit sphere representing the configuration space. The NM time-dependent unknowns of our finite element approximations are probabilities corresponding to discrete space locations and orientation angles. The framework of alternating-direction methods enables us to update the numerical solution in parallel by solving N evolution equations on the sphere and M three-dimensional advection equations in each (pseudo-)time step. To ensure positivity preservation as well as the normalization property of the probability density, we perform algebraic flux correction for each equation and synchronize the correction factors corresponding to different orientation angles. The velocity field for the spatial advection step is obtained using a Schur complement method to solve a generalized system of the incompressible Navier–Stokes equations (NSE). Fiber-induced subgrid-scale effects are taken into account using an effective stress tensor that depends on the second- and fourth-order moments of the orientation density function. Numerical studies are performed for individual subproblems and for the coupled FPE-NSE system.

Hybridised multigrid preconditioners for a compatible finite‐element dynamical core

AbstractCompatible finite‐element discretisations for the atmospheric equations of motion have recently attracted considerable interest. Semi‐implicit timestepping methods require the repeated solution of a large saddle‐point system of linear equations. Preconditioning this system is challenging, since the velocity mass matrix is nondiagonal, leading to a dense Schur complement. Hybridisable discretisations overcome this issue: weakly enforcing continuity of the velocity field with Lagrange multipliers leads to a sparse system of equations, which has a similar structure to the pressure Schur complement in traditional approaches. We describe how the hybridised sparse system can be preconditioned with a non‐nested two‐level preconditioner. To solve the coarse system, we use the multigrid pressure solver that is employed in the approximate Schur complement method previously proposed by the some of the authors. Our approach significantly reduces the number of solver iterations. The method shows excellent performance and scales to large numbers of cores in the Met Office next‐generation climate and weather prediction model LFRic.

Mesh-free high-resolution simulation of cerebrocortical oxygen supply with fast Fourier preconditioning.

Oxygen transfer from blood vessels to cortical brain tissue is representative of a class of problems with mixed-domain character. Large-scale efficient computation of tissue oxygen concentration is dependent on the manner in which the tubular network of blood vessels is coupled to the tissue mesh. Models which explicitly resolve the interface between the tissue and vasculature with a contiguous mesh are prohibitively expensive for very dense cerebral microvasculature. We propose a mixed-domain mesh-free technique whereby a vascular anatomical network (VAN) represented as a thin directed graph serves for convection of blood oxygen, and the surrounding extravascular tissue is represented as a Cartesian grid of 3D voxels throughout which oxygen is transported by diffusion. We split the network and tissue meshes by the Schur complement method of domain decomposition to obtain a reduced set of system equations for the tissue oxygen concentration at steady state. The use of a Cartesian grid allows the corresponding matrix equation to be solved approximately with a fast Fourier transform-based Poisson solver, which serves as an effective preconditioner for Krylov subspace iteration. The performance of this method enables the steady-state simulation of cortical oxygen perfusion for anatomically accurate vascular networks down to single micron resolution without the need for supercomputers.

Использование метода декомпозиции области для распараллеливания моделирования течения вязкой несжимаемой среды методом LS-STAG и дополнительного предобуславливания

As a rule, the main part of the computational costs in the numerical solution of problems in continuum mechanics consists in the solving of large sparse systems of linear algebraic equations. For this reason, efficient parallelization of this particular procedure can significantly speed up the simulation. To solve this problem, two main approaches can be used. The simplest approach consists in parallelizing of matrix-vector operations in a usual iterative solver. It requires several synchronization points and exchanges of coefficients at each iteration of the solver, which does not significantly speed up the simulation process as a whole. Therefore, domain decomposition methods are preferable. These methods involve dividing the computational domain into subdomains, constructing and solving separate problems in them, as well as some procedure to coordinate the solution between subdomains to ensure global convergence. Subdomains can overlap, as in the Schwartz method used in OpenFOAM, or they can be separated by interface sections, on which their own interface task is built, as in the Schur complement method. The latter method is used in this research to construct a parallel algorithm for viscous incompressible flow simulation by using the immersed boundary method LS-STAG with cut-cells and level-set functions. The resulting matrix of the interface system has a block tridiagonal structure. To speed up prototyping, OpenMP parallel programming technology is used in the software implementation of the developed algorithm, so computational experiments are carried out only on systems with shared memory, in particular on individual nodes of the educational and experimental cluster of the Applied Mathematics Department, Bauman Moscow State Technical University. To verify and evaluate the effectiveness of the developed parallel algorithm, a well-studied test problem about simulation of two-dimensional flow around a stationary circular airfoil is considered. Computations on a sequence of meshes with different numbers of subdomains show that the parallel algorithm allows one to obtain the same numerical solution as the original algorithm of the LS-STAG method, and the computed values of the Strouhal number and drag coefficient are in good agreement with the experimental and computational data known in the literature. Experiments demonstrate that the developed algorithm with domain decomposition allows to accelerate simulation even in sequential mode by reducing the number of solver iterations, i.e. the domain decomposition method acts as an additional preconditioner. Due to this property, the acceleration is superlinear when simulating in parallel mode with developed algorithm. This effect persists up to a certain number of subdomains, which depends on the size of the problem.

Planar potential flow on Cartesian grids

Potential flow has many applications, including the modelling of unsteady flows in aerodynamics. For these models to work efficiently, it is best to avoid Biot–Savart interactions. This work presents a grid-based treatment of potential flows in two dimensions and its use in a vortex model for simulating unsteady aerodynamic flows. For flows consisting of vortex elements, the treatment follows the vortex-in-cell approach and solves the streamfunction–vorticity Poisson equation on a Cartesian grid after transferring the circulation from the vortices onto the grid. For sources and sinks, an analogous approach can be followed using the scalar potential. The combined velocity field due to vortices, sinks and sources can then be obtained using the Helmholtz decomposition. In this work, we use several key tools that ensure the approach works on arbitrary geometries, with and without sharp edges. First, the immersed boundary projection method is used to account for bodies in the flow and the resulting body-forcing Lagrange multiplier is identified as the bound vortex sheet strength. Second, sharp edges are treated by decomposing the vortex sheet strength into a singular and non-singular part. To enforce the Kutta condition, the non-singular part can then be constrained to remove the singularity introduced by the sharp edge. These constraints and the Poisson equation are formulated as a saddle-point system and solved using the Schur complement method. The lattice Green's function is used to efficiently solve the discrete Poisson equation with unbounded boundary conditions. The method and its accuracy are demonstrated for several problems.

Legged Robot State Estimation With Dynamic Contact Event Information

This letter presents a state estimation algorithm for the legged robot by defining the problem as a Maximum A Posteriori (MAP) estimation problem and solving the problem with the Gauss-Newton algorithm. Moreover, marginalization by the Schur Complement method is adopted to make a fixed size problem. Each component of the cost function and its Jacobian are derived utilizing the SO(3) manifold structure, while we reparameterize the state with nominal state and variation to make linear algebra and vector calculus applied properly. Furthermore, a slip rejection method is proposed to reduce the erroneous effect of fault modeling of kinematics models. The proposed algorithm is verified by comparison with the Invariant Extended Kalman Filter (IEKF) in real robot experiments on various environments.

Modal analysis of slow varying non-stationary vibration by model updating with Schur complement

This paper proposes a new algorithm for the online monitoring of slow varying modal parameters in vibrating structures subjected to unknown excitations. In this method, an Auto-Regressive (AR) model will be applied in a Short-Time Sliding Window (STSW) on measured signals. The model parameters will be determined and updated through the order and time obtained from the previous computational window. The recursive multivariable least-squares method is utilized in conjunction with the Schur complement method to find solutions. Numerical simulations of the gradual changing system acted by the random excitations and sinusoidal force are presented. Finally, the modal parameters for the multichannel data measured on an experimental steel plate emerging from the water are identified and compared with those obtained with the analytical method and Short-Time Fourier Transform (STFT), providing a powerful means of validating the proposed method.

A scalable parallel preconditioned conjugate gradient method for bundle adjustment

Bundle adjustment is a fundamental problem in computer vision, with important applications such as 3D structure reconstruction from 2D images. This paper focuses on large-scale bundle adjustment tasks, e.g., city-wide 3D reconstruction, which require highly efficient solutions. For this purpose, it is common to apply the Levenberg-Marquardt algorithm, whose bottleneck lies in solving normal equations. The majority of recent methods focus on achieving scalability through modern hardware such as GPUs and distributed systems. On the other hand, the core of the solution, i.e., the math underlying the optimizer for the normal equations, remains largely unimproved since the proposal of the classic parallel bundle adjustment (PBA) algorithm, which increasingly becomes a major limiting factor for the scalability of bundle adjustment. This paper proposes parallel preconditioned conjugate gradient (PPCG) method, a novel parallel method for bundle adjustment based on preconditioned conjugate gradient, which achieves significantly higher efficiency and scalability than existing methods on the algorithmic level. The main idea is to exploit the sparsity of the Hessian matrix and reduce its structure parameters through an effective parallel Schur complement method; the result of this step is then fed into our carefully designed PPCG method which reduces matrix operations that are either expensive (e.g., large matrix reverse or multiplications) or scales poorly to multi-processors (e.g., parallel Reduce operators). Extensive experiments demonstrate that PPCG outperforms existing optimizers by large margins, on a wide range of datasets.

A direction-splitting Navier–Stokes solver on co-located grids

We introduce a new finite-difference solver for the incompressible Navier–Stokes equations that exploits the direction-splitting method proposed by Guermond and Minev in 2010, but is formulated on a co-located grid. The main ingredients of the new solver are: i) the direction-splitting approach adopted for both the momentum and the pressure equations; and ii) the co-located grid approach. The solver is parallelised by the Schur-complement method, and achieves very high performance levels on thousands of processors. Several test cases are proposed to assess the accuracy and efficiency of the method.

Symmetric near‐field Schur's complement preconditioner for hierarchal electric field integral equation solver

In this study, a robust and effective preconditioner for the fast method of moments-based hierarchal electric field integral equation solver is proposed using symmetric near-field Schur's complement method. In this preconditioner, near-field blocks are scaled to a diagonal block matrix and these near-field blocks are replaced with the scaled diagonal block matrix which reduces the near-field storage memory and the overall matrix–vector product time. Scaled diagonal block matrix is further used as a preconditioner and due to the block diagonal form of the final preconditioner, no additional fill-ins are introduced in its inverse. The symmetric property of the near-field blocks is exploited to reduce the preconditioner setup time. Near linear complexity of preconditioner set up and solve times is achieved by near-field block ordering, using graph bandwidth reduction algorithms and compressing the fill-in blocks in preconditioner computation. Preconditioner set up time is reduced to half by using the symmetric property and near-field block ordering. It has been shown using a complexity analysis that the cost of preconditioner construction in terms of computation and memory is linear. Numerical experiments demonstrate an average of 1.5–2.3× speed-up in the iterative solution time over null-field-based preconditioners.

Inner product free iterative solution and elimination methods for linear systems of a three-by-three block matrix form

Large scale systems of algebraic equations are frequently solved by iterative solution methods, such as the conjugate gradient method for symmetric or a generalized conjugate gradient or generalized minimum residual method for nonsymmetric linear systems. In practice, to get an acceptable elapsed computing time when solving large scale problems, one shall use parallel computer platforms. However, such methods involve orthogonalization of search vectors which requires computation of many inner products and, hence, needs global communication of data, which will be costly in computer times. In this paper, we propose various inner product free methods, such as the Chebyshev acceleration method. We study the solution of linear systems arising from optimal control problems for PDEs, such as the edge element discretization of the time-periodic eddy current optimal control problem. Following a discretize-then-optimize scheme, the resulting linear system is of a three-by-three block matrix form. Various solution methods based on an approximate Schur complement and inner product free iterative solution methods for this linear system are analyzed and compared with an earlier used method for two-by-two block matrices with square blocks. The convergence properties and implementation details of the proposed methods are analyzed to show their effectiveness and practicality. Both serial and parallel numerical experiments are presented to further investigate the performance of the proposed methods compared with some other existing methods.

A Novel Linking-Domain Extraction Decomposition Method for Parallel Electromagnetic Transient Simulation of Large-Scale AC/DC Networks

Domain decomposition of the network conductance matrix is one of the efficient approaches to solve large-scale networks in parallel, wherein the most commonly-used non-iterative method is the Schur complement (SC) method. However, the SC method could not obtain the network conductance matrix inversion directly, and the computational cost will increase fast when the overlapping domain expands. In this work, a novel Linking-Domain Extraction (LDE) based decomposition method is proposed, in which the network matrix is expressed as the sum of a linking-domain matrix (LDM) and a diagonal block matrix (DBM) composed of multiple block matrices in diagonal. Through mathematical analysis over LDM, one lemma about the nature of LDM and its proof are proposed. Based on this lemma, the general formulation of the inverse matrix of the sum of LDM and DBM can be found using the Woodbury matrix identity, and based on the formulation the network matrix inversion can be directly computed in parallel to significantly accelerate the matrix inversion process. Test systems were implemented on both the FPGA and GPU parallel architectures, and the simulation results and speed-ups over the SC method and Gauss-Jordan elimination demonstrate the validity and efficiency of the proposed LDE method.

English

We present a parallel solution algorithm for the transient heat equation in one and two spatial dimensions. The problem is discretized in space by the lowest-order conforming finite element method. Further, a one-step time integration scheme is used for the numerical solution of the arising system of ordinary differential equations. For the latter, the parareal method decomposing the time interval into subintervals is employed. It leads to parallel solution of smaller time-dependent problems. At each time slice a pseudostationary elliptic heat equation is solved by means of a domain decomposition method (DDM). In the 2d, case we employ a nonoverlapping Schur complement method, while in the 1d case an overlapping Schwarz DDM is employed. We document computational efficiency, as well as theoretical convergence rates of FEM semi-discretization schemes on numerical examples.

On the design of two-stage multiprojection methods for distributed memory systems

Solving large sparse linear systems, efficiently, on supercomputing infrastructures is a time-consuming component for a wide variety of simulation processes. An effective parallel solver should meet the required specifications, concerning both convergence behavior and scalability. Herewith, a class of two-stage algebraic domain decomposition preconditioning schemes based on the upper Schur complement method is proposed, in order to exploit appropriately distributed memory systems with multicore processors. The design of the method has been focused on homogeneous hybrid parallel systems, i.e., distributed and shared memory systems. However, the proposed method can also be applied to heterogeneous systems, such as cloud infrastructures, or hybrid parallel systems with accelerators, by modifying the workload distribution algorithm and taking into account the different network latencies and bandwidths. The first stage of the proposed schemes is related to the assignment of the subdomains among the workstations of the distributed system, whereas the second stage concerns the further redistribution of the subdomains to each core of a processor. The proposed method utilizes multiprojection techniques, based on semi-aggregated subdomains, leading to improved convergence behavior as the number of subdomains increases. Moreover, a subspace compression technique is used, in order to improve the performance of the preprocessing phase and reduce the memory requirements of the proposed scheme. The preconditioning schemes were combined with a parallel Krylov subspace method, i.e., the parallel preconditioned GMRES(m) method. The convergence behavior, the performance and the scalability of the proposed preconditioning schemes are examined and compared to existing state-of-the-art methods, by conducting several numerical experiments on supercomputing infrastructures.

Improved newton-raphson with schur complement methods for load flow analysis

<p>The determination of power and voltage in the power load flow for the purpose of design and operation of the power system is very crucial in the assessment of actual or predicted generation and load conditions. The load flow studies are of the utmost importance and the analysis has been carried out by computer programming to obtain accurate results within a very short period through a simple and convenient way. In this paper, Newton-Raphson method which is the most common, widely-used and reliable algorithm of load flow analysis is further revised and modified to improve the speed and the simplicity of the algorithm. There are 4 Newton-Raphson algorithms carried out, namely Newton-Raphson, Newton-Raphson constant Jacobian, Newton-Raphson Schur Complement and Newton-Raphson Schur Complement constant Jacobian. All the methods are implemented on IEEE 14-, 30-, 57- and 118-bus system for comparative analysis using MATLAB programming. The simulation results are then compared for assessment using measurement parameter of computation time and convergence rate. Newton-Raphson Schur Complement constant Jacobian requires the shortest computational time.</p>

Schur Complement-based Substructuring of Stiff Multibody Systems with Contact

Substructuring permits parallelization of physics simulation on multi-core CPUs. We present a new substructuring approach for solving stiff multibody systems containing both bilateral and unilateral constraints. Our approach is based on non-overlapping domain decomposition with the Schur complement method, which we extend to systems involving contact formulated as a mixed bounds linear complementarity problem. At each time step, we alternate between solving the subsystem and interface constraint impulses, which leads to the identification of the active constraints. By using the active constraints to compute the effective mass of subsystems within the interface solve, we obtain an exact solution. We demonstrate that our simulations have preferable behavior compared to standard iterative solvers and substructuring techniques based on the exchange of forces at interface bodies. We observe considerable speedups for structured simulations where a user-defined partitioning can be applied, and moderate speedups for unstructured simulations, such as piles of bodies. In the latter case, we propose an automatic partitioning strategy based on the degree of bodies in the constraint graph. Because our method makes use of direct solvers, we are able to achieve interactive and real-time frame rates for a number of challenging scenarios involving large mass ratios, redundant constraints, and ill-conditioned systems.

Fast 2D full-waveform modeling and inversion using the Schur complement approach

In most full-waveform inversion (FWI) problems, sufficient prior information is available to constrain the velocity of certain parts of the model, e.g., the water column or, in some cases, near-surface velocities. We take advantage of this situation and develop a fast Schur-complement-based forward modeling and inversion approach by partitioning the velocity model into two parts. The first part consists of the constrained zone that does not change during the inversion, whereas the second part is the anomalous zone to be updated during the inversion. For this decomposition, we partially factorize the governing system of linear equations by computing a Schur complement for the anomalous zone. The Schur complement system is then solved to compute the fields in the anomalous zone, which are then back substituted to compute the fields in the constrained region. For each successive modeling steps with new anomalous zone velocities, the corresponding Schur complement is easily computed using simple algebra. Because the anomalous part of the model is comparatively smaller than the whole model, considerable computational savings can be achieved using our Schur approach. Additionally, we showed that the Schur complement method maintains the accuracy of standard frequency-domain finite difference formulations, but this comes at a slightly higher peak memory requirement. Our FWI workflow shows reduced runtime by 15%–57% depending upon the depth of the water column without losing any accuracy compared to the standard method.

Higher-order simulation of impressed current cathodic protection systems

In this paper, a surface integral equation formulation for the prediction of corrosion-related electrostatic fields of conducting structures in an electrolyte is presented. The integral equation is discretized using the Locally Corrected Nyström method. An iterative technique based on the Newton-Raphson method is described for use with nonlinear boundary conditions, and a Schur complement method is discussed to help improve the nonlinear solution efficiency. An image plane method is applied when an appropriate planar electrolyte-insulating bounding surface is present. The methods are validated by analysis of a canonical problem with analytic solution and by comparison to literature data. Higher-order solution convergence is observed for the canonical problem results.