PETSc Research Articles

HEMLAB Algorithm Applied to the High-Lift Japan Aerospace Exploration Agency Standard Model

The high-lift Japan Aerospace Exploration Agency (JAXA) standard model (HL-JSM) has been numerically analyzed in order to further validate the HEMLAB code for realistic aircraft configurations. The numerical algorithm is based on highly efficient edge-based data structure for a vertex-based finite volume algorithm on hybrid meshes. The data pattern is arranged to meet the requirements of a vertex-based finite volume algorithm by considering data access patterns and cache efficiency. A fully implicit version of the numerical algorithm has also been implemented based on the Portable, Extensible Toolkit for Scientific Computation (PETSc) library in order to improve the robustness of the algorithm. The resulting algebraic equations, including the one-equation Spalart–Allmaras, are solved in a monolithic manner using the restricted additive Schwarz preconditioner combined with the flexible generalized minimal residual method [FGMRES(m)] Krylov subspace algorithm. The numerical method is also combined with the metric-based anisotropic mesh refinement library Python Adaptive Mesh Geometry Suite–INRIA (pyAMG) from National Institute for Research in Digital Science and Technology in order to improve the numerical accuracy. The numerical algorithm is initially applied to the two-dimensional L1T2 (National High Lift Programme) high-lift system, and then the calculations around the HL-JSM are carried out at relatively high angles of attack. The numerical results with the anisotropic mesh refinement library pyAMG indicate significant improvements in numerical accuracy.

On cost-efficient parallel iterative solvers for 3D frequency-domain seismic multisource viscoelastic anisotropic wave modeling

Solving large sparse linear systems in 3D frequency-domain seismic wave modeling, especially in viscoelastic anisotropic media, poses significant challenges due to the increasing number of discrete moduli and nonzero elements in the linear system matrix. The computational load surpasses that of acoustic or viscoacoustic media, making it even more challenging when dealing with multisource problems. Popular scientific tools for solving a linear system, such as multifrontal massively parallel sparse direct solver (MUMPS), structured matrix package (STRUMPACK), and portable extensible toolkit for scientific computation (PETSc), can be used, but their applicability to our specific problem has not been comprehensively evaluated. Our study aims to tackle the challenges in solving large, sparse, complex-valued symmetric linear systems with multiple right-hand side vectors for 3D frequency-domain seismic wave modeling. We have adopted the preconditioned conjugate gradient iterative algorithms as the foundation for our research, introducing two highly cost-effective parallel iterative solvers: the parallel symmetric successive overrelaxation conjugate gradient (P-SSORCG) and the parallel incomplete Cholesky conjugate gradient (P-ICCG). These novel solvers are subjected to a comprehensive comparative analysis against well-established scientific tools, such as MUMPS, STRUMPACK, and PETSc, in the context of 3D frequency-domain seismic wave modeling. We show their promising performances in a practical 3D SEG/EAGE overthrust model and demonstrate that the grouped P-SSORCG offers an efficient alternative to parallel direct solvers, particularly in situations wherein computational resources are limited.

A parallel implementation of the confined–unconfined aquifer system model for subglacial hydrology: design, verification, and performance analysis (CUAS-MPI v0.1.0)

Abstract. The subglacial hydrological system affects (i) the motion of ice sheets through sliding, (ii) the location of lakes at the ice margin, and (iii) the ocean circulation by freshwater discharge directly at the grounding line or (iv) via rivers flowing over land. For modeling this hydrology system, a previously developed porous-media concept called the confined–unconfined aquifer system (CUAS) is used. To allow for realistic simulations at the ice sheet scale, we developed CUAS-MPI, an MPI-parallel C/C++ implementation of CUAS (MPI: Message Passing Interface), which employs the Portable, Extensible Toolkit for Scientific Computation (PETSc) infrastructure for handling grids and equation systems. We validate the accuracy of the numerical results by comparing them with a set of analytical solutions to the model equations, which involve two types of boundary conditions. We then investigate the scaling behavior of CUAS-MPI and show that CUAS-MPI scales up to 3840 MPI processes running a realistic Greenland setup on the Lichtenberg HPC system. Our measurements also show that CUAS-MPI reaches a throughput comparable to that of ice sheet simulations, e.g., the Ice-sheet and Sea-level System Model (ISSM). Lastly, we discuss opportunities for ice sheet modeling, explore future coupling possibilities of CUAS-MPI with other simulations, and consider throughput bottlenecks and limits of further scaling.

The PetscSF Scalable Communication Layer

PetscSF, the communication component of the Portable, Extensible Toolkit for Scientific Computation (PETSc), is designed to provide PETSc's communication infrastructure suitable for exascale computers that utilize GPUs and other accelerators. PetscSF provides a simple application programming interface (API) for managing common communication patterns in scientific computations by using a star-forest graph representation. PetscSF supports several implementations based on MPI and NVSHMEM, whose selection is based on the characteristics of the application or the target architecture. An efficient and portable model for network and intra-node communication is essential for implementing large-scale applications. The Message Passing Interface, which has been the de facto standard for distributed memory systems, has developed into a large complex API that does not yet provide high performance on the emerging heterogeneous CPU-GPU-based exascale systems. In this article, we discuss the design of PetscSF, how it can overcome some difficulties of working directly with MPI on GPUs, and we demonstrate its performance, scalability, and novel features.

Hybrid Nodal Integral/Finite Element Method for Time-Dependent Convection Diffusion Equation

AbstractNodal integral methods (NIM) are a class of efficient coarse mesh methods that use transverse averaging to reduce the governing partial differential equation(s) (PDE) into a set of ordinary differential equations (ODE). These ODEs or their approximations are analytically solved. These analytical solutions are used to then develop the numerical scheme. Transverse averaging is an essential step in the development of NIM, and hence the standard application of this approach gets restricted to domains that have boundaries parallel to one of the coordinate axes (in 2D) or coordinate planes (in 3D). The hybrid nodal-integral/finite element method (NI-FEM) reported here has been developed to extend the application of NIM to arbitrary domains. NI-FEM is based on the idea that the interior region and the regions with boundaries parallel to the coordinate axes (2D) or coordinate planes (3D) can be solved using NIM, and the rest of the domain can be discretized and solved using FEM. The crux of the hybrid NI-FEM is in developing interfacial conditions at the common interfaces between the NIM regions and FEM regions. Since the discrete variables in the two numerical approaches are different, this requires special treatment of the discrete quantities on the interface between the two different types of discretized elements. We here report the development of hybrid NI-FEM for the time-dependent convection-diffusion equation (CDE) in arbitrary domains. The resulting hybrid numerical scheme is implemented in a parallel framework in fortran and solved using portable, extensible toolkit for scientific computation (PETSc). The preliminary approach to domain decomposition is also discussed. Numerical solutions are compared with exact solutions, and the scheme is shown to be second-order accurate in both space and time. The order of approximations used for the development of the scheme is also shown to be second order. The hybrid method is more efficient compared to standalone conventional numerical schemes like FEM.

Toward performance-portable PETSc for GPU-based exascale systems

The Portable Extensible Toolkit for Scientific computation (PETSc) library delivers scalable solvers for nonlinear time-dependent differential and algebraic equations and for numerical optimization. The PETSc design for performance portability addresses fundamental GPU accelerator challenges and stresses flexibility and extensibility by separating the programming model used by the application from that used by the library, and it enables application developers to use their preferred programming model, such as Kokkos, RAJA, SYCL, HIP, CUDA, or OpenCL, on upcoming exascale systems. A blueprint for using GPUs from PETSc-based codes is provided, and case studies emphasize the flexibility and high performance achieved on current GPU-based systems.

MERF v3.0, a highly computationally efficient non-hydrostatic ocean model with implicit parallelism: Algorithms and validation experiments

Non-hydrostatic ocean models are widely required in coastal and estuarine studies but are limited in their application due to their very high computational expense. To alleviate this issue, this study developed a highly computationally efficient non-hydrostatic ocean model named MERF v3.0 (the third version of the Marine Environment Research and Forecasting model) in which a semi-implicit and variable layers (SIVL) scheme is proposed for dynamic pressure calculation. The SIVL scheme provides freedom to choose the ratio between the number of dynamic pressure layers and the number of velocity layers as required. Moreover, the model adopted numerical calculation framework OpenArray and PETSc (Portable, Extensible Toolkit for Scientific Computation) to achieve implicit parallelism and a concise code structure. Five classical numerical experiments were performed to validate the simulation capacity of this non-hydrostatic ocean model, including the standing wave in a closed basin, surface solitary wave, lock-exchange problem, periodic wave over a bar and tidally induced internal lee wave. The numerical simulation results are consistent with the analytical solutions and experimental data published in the literature. Benefitting from the advantages of OpenArray and PETSc, MERF v3.0 achieves an 87.8% parallel efficiency in the strong scaling test and an 84.1% parallel efficiency in the weak scaling test on 3584 processes, whereas the OpenArray limits the parallel efficiency of MERF v3.0 on more processes, which is beyond the scope of this article and will be further examined in our next study.

Implementation of Thermomechanical Multiphysics in a Large-Scale Three-Dimensional Topology Optimization Code

<div class="section abstract"><div class="htmlview paragraph">Due to the inherent computational cost of multiphysics topology optimization methods, it is a common practice to implement these methods in two-dimensions. However most real-world multiphysics problems are best optimized in three-dimensions, leading to the necessity for large-scale multiphysics topology optimization codes. To aid in the development of these codes, this paper presents a general thermomechanical topology optimization method and describes how to implement the method into a preexisting large-scale three-dimensional topology optimization code. The weak forms of the Galerkin finite element models are fully derived for mechanical, thermal, and coupled thermomechanical physics models. The objective function for the topology optimization method is defined as the weighted sum of the mechanical and thermal compliance. The corresponding sensitivity coefficients are derived using the direct differentiation method and are verified using the complex-step method. The design variables are updated using the method of moving asymptotes so that the objective function is minimized resulting in maximum performance of the structure. Solution of the linear systems of equation is scaled across several supercomputer nodes using the Portable, Extensible Toolkit for Scientific Computation (PETSc) developed by Argonne National Laboratory. The finite element models, sensitivity analysis, optimization algorithm, and solver are all implemented in a freeware, density-based, topology optimization code (written in C++) that has been modified for thermomechanical problems. As such, the implementation of the thermomechanical model is provided. The code is applied to solve an example problem under thermal, mechanical, and thermomechanical boundary conditions. The resulting topologies and field variables are shown with a corresponding discussion.</div></div>

Efficient solution techniques for two-phase flow in heterogeneous porous media using exact Jacobians

Two efficient and scalable numerical solution methods will be compared using exact Jacobians to solve the fully coupled Newton systems arising during fully implicit discretization of the equations for two-phase flow in porous media. These methods use algebraic multigrid (AMG) to solve the linear systems in every Newton step. The algebraic multigrid methods rely on (i) a Schur Complement Reduction (SCR-AMG) and (ii) a Constrained Pressure Residual method (CPR-AMG) to decouple elliptic and hyperbolic contributions. Both methods employ automatic differentiation (AD) to calculate exact Jacobians and are compared with finite difference (FD) approximations of the Jacobian. The superiority of AD is shown by several numerical test cases from the field of CO2 geo-sequestration comprising two- and three-dimensional examples. A weak scaling test on JUQUEEN, a BlueGene/Q supercomputer, demonstrates the efficiency and scalability of both methods. To achieve maximum comparability and reproducibility, the Portable Extensible Toolkit for Scientific Computation (PETSc) is used as framework for the comparison of all solvers.

A scalable matrix-free spectral element approach for unsteady PDE constrained optimization using PETSc/TAO

We provide a new approach for the efficient matrix-free application of the transpose of the Jacobian for the spectral element method for the adjoint-based solution of partial differential equation (PDE) constrained optimization. This results in optimizations of nonlinear PDEs using explicit integrators where the integration of the adjoint problem is not more expensive than the forward simulation. Solving PDE constrained optimization problems entails combining expertise from multiple areas, including simulation, computation of derivatives, and optimization. The Portable, Extensible Toolkit for Scientific computation (PETSc) together with its companion package, the Toolkit for Advanced Optimization (TAO), is an integrated numerical software library that contains an algorithmic/software stack for solving linear systems, nonlinear systems, ordinary differential equations, differential algebraic equations, and large-scale optimization problems and, as such, is an ideal tool for performing PDE-constrained optimization. This paper describes an efficient approach in which the software stack provided by PETSc/TAO can be used for large-scale nonlinear time-dependent problems. Time integration can involve a range of high-order methods, both implicit and explicit. The PDE-constrained optimization algorithm used is gradient-based and seamlessly integrated with the simulation of the physical problem.

Hierarchical finite element-based multi-scale modelling of composite laminates

This paper presents a hierarchic finite element-based computational framework for the multi-scale modelling of composite laminates. Hierarchic finite elements allow changing the approximation order locally or globally without changing the underlying finite element mesh. Both micro- and macro-level structures are discretised with these elements. The macro-level structures of composite laminates are divided into several blocks during the pre-processing stage, and approximation orders are assigned to each block independently. Due to a sharp increase in the interlaminar stresses, higher approximation orders are used in the vicinity of free edges as compared to the rest of the problem domain. This freedom of assigning approximation orders independently to each block provides an efficient and accurate way for modelling composite laminates. The computation framework can either accept the user-defined ply-level homogenised elastic material properties or calculates these directly from the underlying representative volume element consisting of matrix and fibre using the computational homogenisation. The model developed for the computational homogenisation has the flexibility of unified imposition of representative volume element boundary conditions, which allows convenient switching between linear displacement, uniform traction and periodic boundary conditions. The computational framework has additional flexibly of high-performance computing and makes use of state-of-the-art computational libraries including Portable, Extensible Toolkit for Scientific Computation (PETSc) and the Mesh-Oriented datABase (MOAB). Symmetric cross-ply, angle-ply and quasi-isotropic laminates subjected to uniaxial loading are used as test cases to demonstrate the correct implementation and computational efficiency of the developed computational framework.

Parallel solution of saddle point systems with nested iterative solvers based on the Golub‐Kahan Bidiagonalization

SummaryThe Golub‐Kahan bidiagonalization is widely used in the singular value decomposition of rectangular matrices and has been generalized to an iterative solver for symmetric indefinite linear systems with a two‐by‐two block structure. In this work, we present a scalability study of this generalized solver as implemented in a recent release of the parallel numerical library PETSc (Portable, Extensible Toolkit for Scientific Computation). We present an improved solver performance for the two‐dimensional (2D) Stokes equations as compared to previous work. Furthermore, we investigate the performance of different parallel inner solvers in the outer Golub‐Kahan iteration for a three‐dimensional Stokes problem. The study includes parallel sparse direct solvers and multigrid methods. When increasing the number of cores for a fixed total problem size, the solver exhibits good speedups of up to 50% at the 1024 core count. For the tests in which the total problem size grows while the workload in each core stays constant, the parallel performance of the solver scales almost linearly with the increase in the core counts. In particular, the computation time increases only by about 15% when the number of cores increases from 80 to 1024 for a 2D test case.

Slate: extending Firedrake's domain-specific abstraction to hybridized solvers for geoscience and beyond

Abstract. Within the finite element community, discontinuous Galerkin (DG) and mixed finite element methods have become increasingly popular in simulating geophysical flows. However, robust and efficient solvers for the resulting saddle point and elliptic systems arising from these discretizations continue to be an ongoing challenge. One possible approach for addressing this issue is to employ a method known as hybridization, where the discrete equations are transformed such that classic static condensation and local post-processing methods can be employed. However, it is challenging to implement hybridization as performant parallel code within complex models whilst maintaining a separation of concerns between applications scientists and software experts. In this paper, we introduce a domain-specific abstraction within the Firedrake finite element library that permits the rapid execution of these hybridization techniques within a code-generating framework. The resulting framework composes naturally with Firedrake's solver environment, allowing for the implementation of hybridization and static condensation as runtime-configurable preconditioners via the Python interface to the Portable, Extensible Toolkit for Scientific Computation (PETSc), petsc4py. We provide examples derived from second-order elliptic problems and geophysical fluid dynamics. In addition, we demonstrate that hybridization shows great promise for improving the performance of solvers for mixed finite element discretizations of equations related to large-scale geophysical flows.

A three-dimensional hierarchic finite element-based computational framework for the analysis of composite laminates

A three-dimensional hierarchic finite element-based computational framework is developed for the investigation of inter-laminar stresses and displacements in composite laminates of finite width. As compared to the standard finite elements, hierarchic finite elements allow to change the order of approximation both locally and globally without modifying the underlying finite element mesh leading to very accurate results for relatively coarse meshes. In this paper, both symmetric cross-ply and angle-ply laminates subjected to uniaxial tension are considered as test cases. Tetrahedral elements are used for the discretisation of laminates and uniform or global p-refinement is used to increase the order of approximation. Each ply within laminates is modelled as a linear-elastic, homogenous and orthotropic material. With increasing the order of approximation, the developed computational framework is able to capture the complex profiles of inter-laminar stresses and displacements very accurately. Results are compared with reference results from the literature and found in a very good agreement. The computational model is implemented in our in-house finite element software library Mesh-Oriented Finite Element Method (MoFEM). The computational framework has additional flexibly of high-performance computing and makes use of the state-of-the-art computational libraries including Portable, Extensible Toolkit for Scientific Computation (PETSc) and the Mesh-Oriented datABase (MOAB).

An implicit high‐order preconditioned flux reconstruction method for low‐Mach‐number flow simulation with dynamic meshes

SummaryA fully implicit high‐order preconditioned flux reconstruction/correction procedure via reconstruction (FR/CPR) method is developed to solve the compressible Navier‐Stokes equations at low Mach numbers. A dual‐time stepping approach with the second‐order backward differentiation formula (BDF2) is employed to ensure temporal accuracy for unsteady flow simulation. When dynamic meshes are used to handle moving/deforming domains, the geometric conservation law is implicitly enforced to eliminate errors due to the resolution discrepancy between BDF2 and the spatial FR/CPR discretization. The large linear system resulted from the spatial and temporal discretizations is tackled with the restarted generalized minimal residual solver in the PETSc (portable, extensible toolkit for scientific computation) library. Through several benchmark steady and unsteady numerical tests, the preconditioned FR/CPR methods have demonstrated good convergence and accuracy for simulating flows at low Mach numbers. The new flow solver is then used to study the effects of Mach number on unsteady force generation over a plunging airfoil when operating in low‐Mach‐number flows. It is observed that weak compressibility has a significant impact on thrust generation but has a negligible effect on lift generation of an oscillating airfoil.

Parallel Implementation of a PETSc-Based Framework for the General Curvilinear Coastal Ocean Model

The General Curvilinear Coastal Ocean Model (GCCOM) is a 3D curvilinear, structured-mesh, non-hydrostatic, large-eddy simulation model that is capable of running oceanic simulations. GCCOM is an inherently computationally expensive model: it uses an elliptic solver for the dynamic pressure; meter-scale simulations requiring memory footprints on the order of 10 12 cells and terabytes of output data. As a solution for parallel optimization, the Fortran-interfaced Portable–Extensible Toolkit for Scientific Computation (PETSc) library was chosen as a framework to help reduce the complexity of managing the 3D geometry, to improve parallel algorithm design, and to provide a parallelized linear system solver and preconditioner. GCCOM discretizations are based on an Arakawa-C staggered grid, and PETSc DMDA (Data Management for Distributed Arrays) objects were used to provide communication and domain ownership management of the resultant multi-dimensional arrays, while the fully curvilinear Laplacian system for pressure is solved by the PETSc linear solver routines. In this paper, the framework design and architecture are described in detail, and results are presented that demonstrate the multiscale capabilities of the model and the parallel framework to 240 cores over domains of order 10 7 total cells per variable, and the correctness and performance of the multiphysics aspects of the model for a baseline experiment stratified seamount.

PERFORMANCE EVALUATION OF DIFFERENT LINEAR EQUATION SOLVERS FOR SOLVING NONLINEAR FE PROBLEMS ON MULTICORE ARCHITECTURES

<p>The aim of this paper is to evaluate the performance of existing parallel linear equation solvers to solving large-scale, nonlinear finite element analysis problems on systems with distributed memory. The parallel approach allows us to take an advantage of the distributed memory enabling forming large system matrices and of multiple processing units to achieve significant speedups. Our study is based on comparison of parallel direct solver and parallel iterative solver implemented in SuperLU DIST library from Portable, Extensible Toolkit for Scientific Computation (PETSc). Both considered solvers are designed for distributed system memory model and are based on a Massage Passing Interface (MPI).</p><p>The efficiency of individual solvers is evaluated on a selected benchmark problems, with different solution strategies by comparing computation times and obtained speedups.</p>

Distributing Load Flow Computations Across System Operators Boundaries Using the Newton–Krylov–Schwarz Algorithm Implemented in PETSc

The upward trends in renewable energy penetration, cross-border flow volatility and electricity actors’ proliferation pose new challenges in the power system management. Electricity and market operators need to increase collaboration, also in terms of more frequent and detailed system analyses, so as to ensure adequate levels of quality and security of supply. This work proposes a novel distributed load flow solver enabling for better cross border flow analysis and fulfilling possible data ownership and confidentiality arrangements in place among the actors. The model exploits an Inexact Newton Method, the Newton–Krylov–Schwarz method, available in the portable, extensible toolkit for scientific computation (PETSc) libraries. A case-study illustrates a real application of the model for the TSO–TSO (transmission system operator) cross-border operation, analyzing the specific policy context and proposing a test case for a coordinated power flow simulation. The results show the feasibility of performing the distributed calculation remotely, keeping the overall simulation times only a few times slower than locally.

A PETSc parallel‐in‐time solver based on MGRIT algorithm

SummaryWe address the development of a modular implementation of the MGRIT (MultiGrid‐In‐Time) algorithm to solve linear and nonlinear systems that arise from the discretization of evolutionary models with a parallel‐in‐time approach in the context of the PETSc (the Portable, Extensible Toolkit for Scientific computing) library. Our aim is to give the opportunity of predicting the performance gain achievable when using the MGRIT approach instead of the Time Stepping integrator (TS). To this end, we analyze the performance parameters of the algorithm that provide a‐priori the best number of processing elements and grid levels to use to address the scaling of MGRIT, regarded as a parallel iterative algorithm proceeding along the time dimension.

Aproximação de uma solução para a equação de Richards-2D pelo método de volumes finitos com o auxílio dos algoritmos de Picard-Krylov

Este artigo propõe-se solucionar a equação de Richards pelo método de volumes finitos em duas dimensões empregando o método de Picard com maior eficiência computacional. Para tanto foram empregadas técnicas iterativas de resolução de sistemas lineares baseadas do espaço de Krylov com matrizes pré-condicionadoras, concomitantemente com auxílio da biblioteca numérica Portable, Extensible Toolkit for Scientific Computation (PETSc). Os resultados indicam que quando se resolve a equação de Richards considerando-se o método de PICARD-KRYLOV, não importando o modelo de avaliação do solo, a melhor combinação para resolução dos sistemas lineares é KSPBCGS PCSOR. Por outro lado, quando se utiliza as equações de van Genuchten (1980) opta-se pela combinação KSPCG PCSOR. Por fim, o artigo traz contribuições na área da matemática, especificamente em métodos numéricos aplicados na resolução de problemas de fluxo em meios porosos não saturados. Em particular a modelagem proposta também poderá auxiliar no entendimento da recarga de aquíferos subterrâneos