Parallel Domain Decomposition Research Articles

Domain decomposition methods and acceleration techniques for the phase field fracture staggered solver

AbstractThe phase field modeling of fracture is able to simulate the nucleation and the propagation of complex crack patterns. However, the relatively small internal lengths that are required usually lead to very fine meshes and high computational costs, especially for three‐dimensional applications. In the present work, additional cost also comes from the implicit dynamics regularization of unstable crack propagations which potentially leads to a large variation of time steps when switching from quasi‐static to dynamic regimes. To reduce the time to solution in this context, this study proposes a domain decomposition framework and acceleration techniques for the phase field fracture staggered solver. The displacement subproblem and the phase field one are solved with parallel domain decomposition solvers. Dual domain decomposition methods provide low cost preconditioner well adapted to the phase field subproblem. For displacement subproblems undergoing unstable crack propagations, primal domain decomposition methods are preferred to be less sensitive to the treatment of floating substructures. Preconditioners performances are assessed and scalability studies over academic test cases, up to 324 subdomains, are presented. Finally, the robustness of the approach is illustrated on two semi‐industrial simulations.

A computational framework for pharmaco‐mechanical interactions in arterial walls using parallel monolithic domain decomposition methods

AbstractA computational framework is presented to numerically simulate the effects of antihypertensive drugs, in particular calcium channel blockers, on the mechanical response of arterial walls. A stretch‐dependent smooth muscle model by Uhlmann and Balzani is modified to describe the interaction of pharmacological drugs and the inhibition of smooth muscle activation. The coupled deformation‐diffusion problem is then solved using the finite element software FEDDLib and overlapping Schwarz preconditioners from the Trilinos package FROSch. These preconditioners include highly scalable parallel GDSW (generalized Dryja–Smith–Widlund) and RGDSW (reduced GDSW) preconditioners. Simulation results show the expected increase in the lumen diameter of an idealized artery due to the drug‐induced reduction of smooth muscle contraction, as well as a decrease in the rate of arterial contraction in the presence of calcium channel blockers. Strong and weak parallel scalability of the resulting computational implementation are also analyzed.

A parallel domain decomposition method for identifying the space-time dependent diffusion coefficients of 3D parabolic problems

A parallel domain decomposition method for identifying the space-time dependent diffusion coefficients of 3D parabolic problems

Massively parallel numerical simulation of 3D shock wave propagation based on JASMIN framework

An efficient massively parallel computing program for shock wave propagation named Shock3d is developed for shock wave damage mechanism in complex scenes. With the finite volume method, the program is built upon JASMIN (J parallel Adaptive Structured Mesh applications INfrastructure) and employs a component-based hierarchical organization structure. Using the patch-based core algorithm, Shock3d enables high parallel efficiency and strong parallel scalability with a parallel domain decomposition method and an image domain-based communication mode. It not only considers gravity, real gas effects, and boundaries in complex scenes, but also supports dynamic load balancing (DLB) and three-dimensional structured meshes and allows the parallel solution of hundreds of millions of mesh problems. Typical examples were used to verify the validity, parallelism, and high-resolution simulation results of the program developed. In the scalability test, where the propagation of shock waves through hundreds of millions of mesh elements was simulated on 300,000 CPU cores, the program delivered parallel efficiency of 60.84%.

A Dimension-Oblivious Domain Decomposition Method Based on Space-Filling Curves

In this paper we present an algebraic dimension-oblivious two-level domain decomposition solver for discretizations of elliptic partial differential equations. The proposed parallel domain decomposition solver is based on a space-filling curve partitioning approach that is applicable to any discretization, i.e., it directly operates on the assembled matrix equations. Moreover, it allows for the effective use of arbitrary processor numbers independent of the dimension of the underlying partial differential equation while maintaining its optimal convergence behavior. This is the core property required to attain a combination technique solver with extreme scalability which can utilize exascale parallel systems efficiently. Moreover, this approach provides a basis for the development of a fault-tolerant solver for the numerical treatment of high-dimensional Problems. To achieve the required data redundancy we are therefore concerned with large overlaps of our domain decomposition, which we construct in any dimension via space-filling curves. In this paper, we propose a space-filling curve based domain decompositoin solver and present its convergence properties and scaling behavior. The results of our numerical experiments clearly show that our approach provides optimal convergence and scaling behavior in an arbitrary dimension utilizing arbitrary processor numbers.

Scalable computational kernels for mortar finite element methods

Targeting simulations on parallel hardware architectures, this paper presents computational kernels for efficient computations in mortar finite element methods. Mortar methods enable a variationally consistent imposition of coupling conditions at high accuracy, but come with considerable numerical effort and cost for the evaluation of the mortar integrals to compute the coupling operators. In this paper, we identify bottlenecks in parallel data layout and domain decomposition that hinder an efficient evaluation of the mortar integrals. We then propose a set of computational strategies to restore optimal parallel communication and scalability for the core kernels devoted to the evaluation of mortar terms. We exemplarily study the proposed algorithmic components in the context of three-dimensional large-deformation contact mechanics, both for cases with fixed and dynamically varying interface topology, yet these concepts can naturally and easily be transferred to other mortar applications, e.g. classical meshtying problems. To restore parallel scalability, we employ overlapping domain decompositions of the interface discretization independent from the underlying volumes and then tackle parallel communication for the mortar evaluation by a geometrically motivated reduction of ghosting data. Using three-dimensional contact examples, we demonstrate strong and weak scalability of the proposed algorithms up to 480 parallel processes as well as study and discuss improvements in parallel communication related to mortar finite element methods. For the first time, dynamic load balancing is applied to mortar contact problems with evolving contact zones, such that the computational work is well balanced among all parallel processors independent of the current state of the simulation.

Acceleration of a parallel BDDC solver by using graphics processing units on subdomains

An approach to accelerating a parallel domain decomposition (DD) solver by graphics processing units (GPUs) is investigated. The solver is based on the Balancing Domain Decomposition Method by Constraints (BDDC), which is a nonoverlapping DD technique. Two kinds of local matrices are required by BDDC. First, dense matrices corresponding to local Schur complements of interior unknowns are constructed by the sparse direct solver. These are further used as part of the local saddle-point problems within BDDC. In the next step, the local matrices are copied to GPUs. Repeated multiplications of local vectors with the dense matrix of the Schur complement are performed for each subdomain. In addition, factorizations and backsubstitutions with the dense saddle-point subdomain matrices are also performed on GPUs. Detailed times of main components of the algorithm are measured on a benchmark Poisson problem. The method is also applied to an unsteady problem of incompressible flow, where the Krylov subspace iterations are performed repeatedly in each time step. The results demonstrate the potential of the approach to speed up realistic simulations up to 5 times with a preference towards large subdomains.

A closest point method library for PDEs on surfaces with parallel domain decomposition solvers and preconditioners

The DD-CPM software library provides a set of tools for the discretization and solution of problems arising from the closest point method (CPM) for partial differential equations on surfaces. The solvers are built on top of the well-known PETSc framework, and are supplemented by custom domain decomposition (DD) preconditioners specific to the CPM. These solvers are fully compatible with distributed memory parallelism through MPI. This library is particularly well suited to the solution of elliptic and parabolic equations, including many reaction-diffusion equations. The software is detailed herein, and a number of sample problems and benchmarks are demonstrated. Finally, the parallel scalability is measured.

A highly parallel fully implicit domain decomposition method for the simulation of the hemodynamics of a patient-specific artery at the full-body scale

The numerical simulation of blood flows in the human body with a certain level of clinical accuracy is important for the understanding of the human physiology. The success of the modeling relies on a robust numerical method with the corresponding software that can handle the complex geometry, the complex fluid flows and run efficiently on a supercomputer. In this work, we introduce a highly parallel domain decomposition method to solve the three-dimensional incompressible Navier-Stokes equations on a patient-specific artery at the full-body scale from neck to feet with 222 outlets and a minimum diameter around 1.0 mm. A locally refined, unstructured mesh is used to resolve the complex fluid flow. Moreover, a two-level method is introduced to determine the model parameters in the Windkessel outlet boundary condition to guarantee clinically correct flow distributions to 14 major regions. A fully implicit Newton-Krylov-Schwarz method is used to solve the nonlinear algebraic system at each time step and numerical experiments show that the proposed method is robust with respect to the complex geometry, the graph-based partition of the complex mesh, the ill-conditioned sparse systems with locally dense blocks, and different model parameters and is scalable with up to 15,360 processor cores. With the proposed method, one simulation of the blood flow in a full-body arterial network can be obtained in about 8 hours per cardiac cycle, which enables its potential use in a wide range of clinical scenarios.

A parallel adaptive finite-element approach for 3-D realistic controlled-source electromagnetic problems using hierarchical tetrahedral grids

SUMMARY To effectively and efficiently interpret or invert controlled-source electromagnetic (CSEM) data which are recorded in areas with the kind of complex geological environments and arbitrary topography that are typical, 3-D CSEM forward modelling software that can quickly solve large-scale problems, provide accurate electromagnetic responses for complex geo-electrical models and can be easily incorporated into inversion algorithms are required. We have developed a parallel goal-oriented adaptive mesh refinement finite-element approach for frequency-domain 3-D CSEM forward modelling with hierarchical tetrahedral grids that can offer accurate electromagnetic responses for large-scale complex models and that can efficiently serve for inversion. The approach uses the goal-oriented adaptive vector finite element method to solve the total electric field vector equation. The geo-electrical model is discretized by unstructured tetrahedral grids which can deal with complex underground geological models with arbitrary surface topography. Different from previous adaptive finite element software working on unstructured tetrahedral grids, we have utilized a novel mesh refinement technique named the longest edge bisection method to generate hierarchically refined grids. As the refined grids are nested into the coarse grids, the refinement technique can precisely map the electrical parameters of inversion grids onto the forward modelling grids so that the extra numerical errors generated by the inconsistency of electrical parameters between inversion grids and forward modelling grids are eliminated. In addition, we use the parallel domain-decomposition technique to further accelerate the computations, and the flexible generalized minimum residual solver (FGMRES) with an auxiliary Maxwell solver pre-conditioner to solve the final large-scale system of linear equations. In the end, we validate the performance of the proposed scheme using two synthetic models and one realistic model. We demonstrate that accurate electromagnetic fields can be obtained by comparison with the analytic solutions and that the code is highly scalable for large-scale problems with millions or even hundreds of millions of unknowns. For the synthetic 3-D model and the realistic model with complex geometry, our solutions match well with the results calculated by an existing 3-D CSEM forward modelling code. Both synthetic and realistic examples demonstrate that our newly developed code is an effective, efficient forward modelling engine for interpreting CSEM field data acquired in areas of complex geology and topography.

Convergence of restricted additive Schwarz with impedance transmission conditions for discretised Helmholtz problems

The Restricted Additive Schwarz method with impedance transmission conditions, also known as the Optimised Restricted Additive Schwarz (ORAS) method, is a simple overlapping one-level parallel domain decomposition method, which has been successfully used as an iterative solver and as a preconditioner for discretised Helmholtz boundary-value problems. In this paper, we give, for the first time, a convergence analysis for ORAS as an iterative solver—and also as a preconditioner—for nodal finite element Helmholtz systems of any polynomial order. The analysis starts by showing (for general domain decompositions) that ORAS is an unconventional finite element approximation of a classical parallel iterative Schwarz method, formulated at the PDE (non-discrete) level. This non-discrete Schwarz method was recently analysed in [Gong, Gander, Graham, Lafontaine, and Spence, Convergence of parallel overlapping domain decomposition methods for the Helmholtz equation], and the present paper gives a corresponding discrete version of this analysis. In particular, for domain decompositions in strips in 2-d, we show that, when the mesh size is small enough, ORAS inherits the convergence properties of the Schwarz method, independent of polynomial order. The proof relies on characterising the ORAS iteration in terms of discrete ‘impedance-to-impedance maps’, which we prove (via a novel weighted finite-element error analysis) converge as h → 0 h\rightarrow 0 in the operator norm to their non-discrete counterparts.

Parallel fully coupled methods for bound-preserving solution of subsurface flow and transport in porous media

As more powerful supercomputer systems with lots of memory and a large number of computing cores become available, the family of fully coupled algorithms is drawing more attention in scientific and engineering applications, due to its impressive robustness and scalability in extreme-scale simulations. In this work, we introduce and study some parallel domain decomposition preconditioned generalized Newton algorithms for solving the fully coupled and bound-preserving formulation of the subsurface flow and transport problem in porous media. In the approach, we present the active–set reduced–space (ASRS) method to guarantee the nonlinear consistency of the fully coupled system in a monolithic way, and meanwhile ensure the boundedness requirement of the solution. Furthermore, we focus on the application of the overlapping additive Schwarz preconditioning technique to accelerate the linear convergence of Newton iterations and be beneficial to the scalability of the inner linear solver. We present some numerical experiments to demonstrate the parallel scalability of the proposed algorithm on the Shaheen-II supercomputer with thousands of processors. In particular, the numerical results also show that the fully coupled framework has the nature of bound preservation and is efficient for the proposed reservoir flow problems.

Overlapping Schwarz methods with GenEO coarse spaces for indefinite and nonself-adjoint problems

Abstract Generalized eigenvalue problems on the overlap(GenEO) is a method for computing an operator-dependent spectral coarse space to be combined with local solves on subdomains to form a robust parallel domain decomposition preconditioner for elliptic PDEs. It has previously been proved, in the self-adjoint and positive-definite case, that this method, when used as a preconditioner for conjugate gradients, yields iteration numbers that are completely independent of the heterogeneity of the coefficient field of the partial differential operator. We extend this theory to the case of convection–diffusion–reaction problems, which may be nonself-adjoint and indefinite, and whose discretizations are solved with preconditioned GMRES. The GenEO coarse space is defined here using a generalized eigenvalue problem based on a self-adjoint and positive-definite subproblem. We prove estimates on GMRES iteration counts that are independent of the variation of the coefficient of the diffusion term in the operator and depend only very mildly on variations of the other coefficients. These are proved under the assumption that the subdomain diameter is sufficiently small and the eigenvalue tolerance for building the coarse space is sufficiently large. While the iteration number estimates do grow as the nonself-adjointness and indefiniteness of the operator increases, practical tests indicate the deterioration is much milder. Thus, we obtain an iterative solver that is efficient in parallel and very effective for a wide range of convection–diffusion–reaction problems.

An unstructured adaptive mesh refinement approach for computational fluid dynamics of reacting flows

A new Direct Numerical Simulation (DNS) code HAMISH with Adaptive Mesh Refinement (AMR) has been developed to simulate compressible reacting flow in a computationally economical manner. The focus is on problems where high gradients of temperature, density and species mass fraction remain localised, for example within the interior structure of flames. The resolution requirements of such local high-gradient regions often determine the global mesh spacing when uniform meshes are used. The numerical framework in HAMISH is based on an unstructured finite-volume approach together with a Runge-Kutta algorithm for time-stepping. The unstructured Cartesian mesh used in HAMISH allows for adaptive refinement and de-refinement of the local mesh depending on the demands of the flow physics and numerical accuracy. The code utilises an octree data structure and a Morton space-filling curve, allowing efficient cell division as well as straightforward parallel domain decomposition. The spatial discretisation scheme uses fourth-order polynomial reconstruction to evaluate the inter-cell fluxes. A third-order three-step explicit Runge-Kutta time-marching algorithm is used together with adaptive time-stepping with proportional-integral-type error control using an embedded Runge-Kutta scheme. The code is fully parallelised through domain decomposition over an arbitrary number of processors and exhibits a high level of parallel efficiency. In this paper the key capabilities of HAMISH are demonstrated based on a number of test cases. These test cases include (a) simulation of 1-D planar laminar premixed flames to demonstrate the capability of AMR in capturing the sharp gradients within the flame; (b) fully developed laminar channel flow, which shows that the code is capable of refining the mesh in the boundary layer next to the wall; (c) 2-D expanding flame under quiescent laminar condition, which demonstrates that the code is capable of dynamic local mesh refinement based on the position of the flame front; (d) 3-D non-reacting Taylor-Green vortex which demonstrates HAMISH can deal with vortical motion typical of turbulent flows; (e) 3-D premixed turbulent flame propagation under isotropic homogeneous decaying turbulence to demonstrate that the code can be employed for production level DNS of turbulent reacting flows.

Wave Attenuation by Vegetation: Model Implementation and Validation Study

Wave–vegetation interaction is implemented in the WAVEWATCH III (WW3) model. The vegetation sink term followed the early formulations of Dalrymple et al. (Journal of Waterway, Port, Coastal, and Ocean Engineering, 1984, 110, 67–79), which focused on monochromatic waves and vegetation approximated as an array of rigid, vertical cylinders, and was later expanded by Mendez and Losada (Coastal Engineering, 2004, 51, 103–118) for random wave transformations over mildly sloping vegetation fields under breaking and nonbreaking conditions assuming a Rayleigh distribution of wave heights. First, validation is carried out for 63 laboratory cases (Anderson and Smith, 2014) with homogeneous vegetation fields for single and double-peak wave spectra. Then, a field case application is conducted to assess the wave attenuation in a wetland environment with spatially variable vegetation fields during stormy conditions. The case study uses data collected at the Magothy Bay located in the Chesapeake Bay, United States, during Hurricanes Jose and Maria in 2017. The domain decomposition parallelization and the implicit scheme have been used for the simulations to efficiently resolve complex shorelines and high-gradient wave zones, incorporating dominant physics in the complicated coastal zone, including wave breaking, wave–current interaction, bottom friction and scattering, wave–vegetation interaction, and nonlinearity (Abdolali et al., 2020). The lab validation and field application demonstrate that WW3 is an effective tool for evaluating the capacity of wetland natural or nature-based features to attenuate wave energy to achieve coastal flood risk reduction.

ARCHER - a new Three-Dimensional method of characteristics neutron transport code for Pebble-bed HTR with coarse mesh finite difference acceleration

A new deterministic code ARCHER using Three-Dimensional (3-D) Method of Characteristics (MOC) has been developed by Institute of Nuclear and new Energy Technology (INET), Tsinghua University, aiming to obtain high-fidelity transport solution of pebble-bed High Temperature gas-cooled Reactor (HTR) with explicit pebble-bed geometry. However, 3-D MOC usually suffers from the slow convergence rate without efficient acceleration methods. In this work, the Coarse Mesh Finite Difference (CMFD) based on regular mesh is implemented to accelerate 3-D MOC calculation, which is available for the problems with cylinder and cuboid geometry. The complex geometry of reactor core with randomly stacked spherical fuel elements and irregular helium flow area, as well as the topological relationships between fine Flat Source Region (FSR) mesh for MOC and coarse CMFD mesh, are the key but challenging tasks for high-fidelity transport solution. Constructive Solid Geometry (CSG) and Boolean operation are constructed to describe problem geometry, which has the ability to deal with complex topological relationships between meshes. In addition, tree grid structure is utilized to reduce the time complexity of ray tracing and obtain the grid mapping relationship between coarse mesh and fine mesh. Spatial domain decomposition parallel based on MPI is utilized to alleviate the memory pressure and pursue high computational performance. Ray parallel is further achieved through OpenMP with shared memory. Several problems have been used to verify the geometric modeling capability and computational performance of this new 3-D MOC code for pebble-bed HTR.

Parallel domain decomposition schemes based on finite volume element discretization for nonsteady-state diffusion equations on distorted meshes

In this paper, several domain decomposition algorithms are considered to solve parallelly non-steady diffusion equations. The contributions of this work are threefold. The first one is that we present effective parallel domain decomposition schemes underlying the classical and monotone finite volume element discretization. These schemes adopt a three-stage technique, which is also used in the difference method and the finite volume method. Positivity of the parallel scheme based on monotone finite volume element method is proved. The second one is that we propose a new domain decomposition computational strategy to improve the correction step in prediction-correction frame such that the whole parallel process is simplified. The third one is that we consider a parallel algorithm build on quadratic finite volume element method, which preserves local conservation and has higher convergence rates. These parallel schemes are of intrinsic parallelism since we apply the prediction-correction technique to the interface of subdomains coming from an artificial division of computational domain or multi-medium physical domain. The proposed parallel schemes need only local communication among neighboring processors. Numerical results are presented to illuminate the accuracy, stability and parallelism of the parallel schemes.

An efficient two-level overlapping domain decomposition method for recovering unsteady sources of 3D parabolic problems

We develop a parallel two-level domain decomposition method for the 3D unsteady source identification problem governed by a parabolic partial differential equation (PDE). The domain of the PDE is firstly decomposed into several overlapping subdomains and the original inverse source identification problem is then transformed into smaller independent subproblems defined on these subdomains. Each subproblem is formulated as a PDE-constrained optimization problem with appropriate conditions prescribed on the inner boundaries and discretized by finite element method. The resulting coupled algebraic systems are solved simultaneously by restarted GMRES method with a space-time restricted additive Schwarz preconditioner. When forming the preconditioner, a second level of domain decomposition is introduced for each subdomain. The solutions of these subproblems are combined together to form an approximated global solution to the original inverse problem by discarding the overlapping parts of the solution. Since all the subproblems are solved independently, the two-level domain decomposition method provides higher degree of parallelism and saves much computing time. Numerical experiments conducted on a supercomputer with thousands of processor cores validate the efficiency and robustness of the proposed approach.

Spectral, Tensor and Domain Decomposition Methods for Fractional PDEs

Abstract Fractional PDEs have recently found several geophysics and imaging science applications due to their nonlocal nature and their flexibility in capturing sharp transitions across interfaces. However, this nonlocality makes it challenging to design efficient solvers for such problems. In this paper, we introduce a spectral method based on an ultraspherical polynomial discretization of the Caffarelli–Silvestre extension to solve such PDEs on rectangular and disk domains. We solve the discretized problem using tensor equation solvers and thus can solve higher-dimensional PDEs. In addition, we introduce both serial and parallel domain decomposition solvers. We demonstrate the numerical performance of our methods on a 3D fractional elliptic PDE on a cube as well as an application to optimization problems with fractional PDE constraints.

Discontinuous boundary elements for steady-state fluid flow problems in discrete fracture networks

Modeling fluid flow in three-dimensional (3D) Discrete Fracture Networks (DFNs) is of relevance in many engineering applications, such as shale oil/gas production, geothermal energy extraction, nuclear waste disposal and CO2 sequestration. A new Boundary Element Method (BEM) technique with discontinuous quadratic elements, in conjunction with a parallel Domain Decomposition Method (DDM), is presented for the simulation of the steady-state fluid flow in DFNs, consisting of stochastically generated 3D planar fractures, arbitrarily oriented, and having differing hydraulic properties. Numerical examples characterized by DFNs of increasing complexity are proposed to show the accuracy and the efficiency of the presented technique, that provides good approximations of the fluid flow around domain interfaces, where the solution usually displays sharp gradients, like around intersections between traces (the segments originated by the intersection between two fractures), intersections between traces and fracture boundaries, or intersections between fractures and wellbores. The conjunction with a DDM approach is a promising strategy to speed up the computations, by also exploiting the advantages of parallel computing techniques. The technique is implemented in the code PyDFN3D, available at https://github.com/BinWang0213/PyDFN3D.