PaperlayerView exact revision on arXiv

Source-linked AI summary

PetIGA: A Framework for High-Performance Isogeometric Analysis

Lisandro Dalcin, Nathan Collier, Philippe Vignal, Adriano M. A. Cortes, V. M. Calo

arXiv:1305.4452v3cs.MSmath.NA

TL;DR

PetIGA addresses the need for scalable, reusable software for isogeometric analysis of partial differential equations. It builds parallel Galerkin assembly and solver interfaces on PETSc, then demonstrates the framework on nonlinear applications. The framework scales to 4096 cores on a Navier–Stokes problem, while multipatch support remains future work.

  • Problem

    Isogeometric analysis needs scalable reusable software for solving linear, nonlinear, static, and transient PDEs with NURBS-based discretizations.

  • Method

    PetIGA reuses PETSc algorithms and data structures for parallel matrix and vector assembly and interfaces with PETSc solvers for NURBS-based Galerkin discretizations.

  • Results

    PetIGA scales well to 4096 cores on the Navier–Stokes problem and supports challenging nonlinear applications in solid and fluid mechanics.

  • Takeaways & Limitations

    The framework provides a robust, versatile platform suited to large-scale isogeometric PDE applications.

  • Takeaways & Limitations

    Multipatch geometries are not yet supported and extending PetIGA while retaining its usability and distributed-memory capabilities is nontrivial.

Abstract

from arXiv · show

We present PetIGA, a code framework to approximate the solution of partial differential equations using isogeometric analysis. PetIGA can be used to assemble matrices and vectors which come from a Galerkin weak form, discretized with Non-Uniform Rational B-spline basis functions. We base our framework on PETSc, a high-performance library for the scalable solution of partial differential equations, which simplifies the development of large-scale scientific codes, provides a rich environment for prototyping, and separates parallelism from algorithm choice. We describe the implementation of PetIGA, and exemplify its use by solving a model nonlinear problem. To illustrate the robustness and flexibility of PetIGA, we solve some challenging nonlinear partial differential equations that include problems in both solid and fluid mechanics. We show strong scaling results on up to 4096 cores, which confirm the suitability of PetIGA for large scale simulations.

1. Introduction

The introduction motivates isogeometric analysis as a way to use CAD representations directly and presents PetIGA as a scalable, reusable PETSc-based framework for large-scale PDE computation.

  • Motivation: Isogeometric analysis uses NURBS bases to preserve CAD geometry while discretizing partial differential equations.This avoids constructing an intermediate analysis-suitable geometrical description.
  • Motivation: Higher-order continuous NURBS spaces support higher-order PDEs and can provide superior approximability per degree of freedom compared with C0 finite element spaces.Their inter-element continuity can be chosen arbitrarily across spatial dimensions.
  • Framework: The paper argues that reusable scientific software preserves accumulated expertise, while acknowledging maintenance, funding, and interoperability costs.The authors conclude that the benefits of reuse outweigh these risks.
  • Framework: PETSc combines scalable PDE algorithms and data structures with MPI-based communication, high-level collective interfaces, and runtime customization.Its components include matrices, vectors, Krylov methods, preconditioners, distributed arrays, and nonlinear solvers.
  • Framework: PetIGA reuses PETSc algorithms and data structures to provide parallel matrix and vector assembly and access to PETSc solvers.The framework targets isogeometric analysis and is freely available under active development.

2. Implementation

The implementation develops B-spline and NURBS basis functions, their tensor-product extensions, and derivatives needed for mapped higher-order isogeometric analysis.

  • B-spline basis functions: B-spline basis functions are defined from a non-decreasing knot vector using the Cox–de Boor recursion formula.The knot vector contains the real-valued knots that determine the basis construction.
  • B-spline basis functions: Basis derivatives are obtained by differentiating the B-spline expressions, with computationally efficient algorithms used for evaluating values and derivatives.The general kth-derivative formula follows repeated differentiation.
  • Tensor product basis functions: Tensor products extend one-dimensional B-spline bases to multidimensional spaces using separate knot vectors and polynomial degrees for each coordinate direction.The three-dimensional construction combines basis functions associated with the ξ, η, and ζ directions.
  • Tensor product basis functions: Multidimensional B-spline functions are indexed globally through a dimension-independent notation for parametric coordinates.For three dimensions, the global index combines the indices from the three coordinate directions.
  • NURBS basis functions and derivatives: NURBS basis functions are constructed rationally from B-spline basis functions and projective weights, with derivatives computed through derivatives of the weighting function.The implementation expresses first-, second-, and third-order derivatives using the chain rule.
  • NURBS basis functions and derivatives: Spatial derivatives require transforming parametric derivatives through the isoparametric mapping, including higher-order mapping derivatives for higher-order PDEs.The paper provides the derivation of these higher-order derivatives in an appendix.

2.4. Periodic boundary conditions

PetIGA imposes periodicity directly in the B-spline function space by unclamping knot vectors and eliminating redundant boundary basis functions. The section generalizes unclamping to arbitrary continuity orders and mapped geometries.

  • Periodic construction: Periodic boundary conditions are embedded in the function space by unclamping open knot vectors instead of imposing coefficient constraints.This approach is presented as simpler and more general than building constraint equations.
  • Continuity control: 0 ≤k ≤p −1 determines the periodic continuity and requires redefining k + 1 knot values at both domain ends.The construction changes knot values at the left and right boundaries.
  • Examples: C2 cubic B-spline examples show varying continuity across the periodic boundary after unclamping the open knot vector with equations (14) and (15).Unique basis functions are globally numbered and colored to identify those crossing the periodic interface.
  • Periodic construction: Continuity order k reduces the number of basis functions to n −k, with out-of-range indices wrapped modulo n −k.The mapping handles redundant basis functions eliminated at one domain boundary.

2.5. Adjacency graph

The adjacency-graph procedure identifies overlapping B-spline supports and uses them to determine sparse matrix structure. This supports efficient assembly and extends naturally to multiple dimensions.

  • Purpose: The adjacency graph records interacting degrees of freedom for sparse-matrix preallocation and specialized Jacobian differentiation.Knowing the nonzero pattern is important for efficient CSR or CSC matrix assembly and colored finite differences.
  • Algorithm: Algorithm 2 finds the left-most and right-most basis indices whose supports overlap the support of basis function N_i,p.In one dimension, every basis between ℓ and r is adjacent to the i-th basis function.
  • Algorithm: In dimensions higher than one, adjacency is computed from the corresponding index ranges for each spatial dimension.The multidimensional graph is assembled from the one-dimensional overlap intervals.
  • Example: Figure 3 applies algorithm 2 to cubic basis functions with varying inter-element continuity and displays the resulting sparse-matrix nonzero pattern.The basis functions and matrix structure are shown as corresponding panels.

2.6. Partitioning

PetIGA partitions isogeometric elements and nodes across processes, using local grids, ghost values, and locality-aware numbering to support parallel assembly. Its abstraction hides parallel communication and sparse data-structure management while users provide problem-specific evaluations.

  • PetIGA partitions elements and nodes across processes using data structures tailored to higher-continuous, tensor-product isogeometric spaces.
  • Block-contiguous global numbering improves locality and parallel performance, while natural numbering remains preferable for persistence tasks such as checkpoint/restart.
  • Local grids combine strictly-local and ghost nodes so each process can evaluate elements whose supporting basis functions span neighboring processes.
  • Global vectors are assembled from local arrays, and row-distributed sparse matrices concentrate many entries in denser diagonal blocks, improving matrix-vector products.
  • Residual assembly gathers off-process values into local vectors, evaluates integrands at quadrature points, and returns local contributions for global assembly.
  • Users provide problem-specific integrands while PetIGA hides sparse matrix, vector assembly, and parallelism; colored finite differences exploit sparsity but require more colors for higher-order bases.

2.9. Additional features

PetIGA supports standard tensor-product quadrature and additional basis options, while igakit addresses the difficulty of creating NURBS geometries. A bent pipe example illustrates igakit’s successive geometry-construction operations.

  • PetIGA uses tensor-product Gauss–Legendre quadrature and also supports Gauss–Lobatto rules and Lagrange-interpolant tensor-product bases.
  • NURBS volume representation is cumbersome and largely manual, motivating igakit as a package for manipulating knot vectors, control points, and hand-defined geometries.
  • The bent pipe is built by generating an annulus, extruding it into a pipe, revolving it into a 90° elbow, and joining the pieces parametrically.

3. Example

The example solves the nonlinear Bratu boundary-value problem with PetIGA’s Galerkin formulation and user-defined residual and Jacobian routines. The implementation configures discretization, boundary conditions, and parallel setup while the framework handles quadrature and assembly details.

  • 3.1. The Bratu equation: The Bratu equation is a nonlinear second-order boundary-value problem modeling steady-state reaction and heat conduction, with no solution above λmax = 6.80812.
  • PetIGA derives a finite-dimensional Galerkin problem by representing trial and weighting functions as combinations of basis functions and control variables.
  • Newton’s method requires a Jacobian J = ∂R/∂U, whose entries are implemented from the nonlinear residual formulation.
  • The residual routine receives quadrature context, local control variables, an output array for local residual contributions, and problem-specific parameters.
  • At each quadrature point, PetIGA computes the solution and gradient, evaluates the residual integrand, and internally handles quadrature weights and the geometry Jacobian.
  • The main program initializes PETSc and the IGA object, sets dimensions, degrees of freedom, options, boundary values, and parallel data structures before solving.

4. Applications

PetIGA is applied to nonlinear problems in solid mechanics, fluid mechanics, and phase-field modeling using isogeometric discretizations and PETSc-based solution strategies. The examples include large-deformation elasticity, time-dependent nonlinear systems, Cahn–Hilliard phase separation, and Navier–Stokes–Korteweg flows.

  • Applications: The applications cover nonlinear partial differential equations from solid mechanics, fluid mechanics, and phase-field modeling.The framework is demonstrated on challenging problems selected to highlight isogeometric analysis strengths.
  • Solid mechanics: A Neo-Hookean cylindrical tube is solved with Newton’s method in an updated-Lagrangian approach over 15 load steps.The tube uses 16 × 64 × 4 quadratic B-spline functions, with MUMPS solving the linearized systems.
  • Time-dependent problems: Time-dependent nonlinear problems are discretized in space with finite-dimensional B-spline subspaces and advanced using the generalized-α method.The method filters high-frequency modes and is second-order and unconditionally stable when its parameters are selected through ρ∞.
  • Phase-field modeling: The Cahn–Hilliard example models binary-mixture phase separation in dimensionless form with periodic boundary conditions and three-dimensional snapshots.Its mobility and chemical potential are nonlinear functions of concentration, and λ represents interface thickness.
  • Fluid mechanics: The Navier–Stokes–Korteweg three-bubble problem is solved in two and three dimensions, showing density and velocity evolution.The three-dimensional visualization uses density isosurfaces at ρ = {0.15, 0.55} and velocity vectors colored by magnitude.

5. Performance

Performance is evaluated through controlled strong-scaling tests for the incompressible Navier–Stokes problem on Stampede. The tests fix discretization, time-step, Newton, and GMRES choices to isolate code performance from changing solver behavior.

  • Test design: Strong scaling is measured for the incompressible Navier–Stokes code using controlled runs with fixed solver iteration counts.Each test uses ten time steps, two Newton iterations per step, and 30 GMRES iterations per Newton iteration.
  • Scaling considerations: Increasing problem size raises the linear-system condition number, while weaker domain-decomposition preconditioners at higher processor counts increase required iterations.
  • Test design: The test problem is flow between two plates with no-slip walls and an initial steady-state laminar flow profile.The domain uses C1 quadratic B-splines with full Gauss–Legendre quadrature and three points per direction per element.
  • Solver configuration: The preconditioner uses block Jacobi with one block per process and ILU(0) within each block.
  • Results: Parallel efficiency on multiple Stampede nodes surpasses 80% and exceeds 90% in all but one case, with local sizes as small as 512 elements per core.The runs use a distributed-memory system with 2 × 8-core Intel Xeon E5-2680 processors per node.

6. Conclusions

The paper concludes that PetIGA provides a scalable PETSc-based implementation of isogeometric analysis for linear, nonlinear, static, and transient problems. It reports strong scaling to 4096 cores, while identifying multipatch geometries as an important future extension.

  • Conclusions: PetIGA implements scalable isogeometric analysis for linear and nonlinear, static and transient problems.The framework is built on PETSc to provide a robust and versatile platform for solving partial differential equations.
  • Conclusions: Strong scaling reaches up to 4096 cores on the Navier–Stokes problem, supporting large-scale applications.
  • Conclusions: PetIGA also achieves excellent performance on shared-memory multicore laptop and desktop computers despite being primarily designed for distributed-memory environments.
  • Future work: The current implementation focuses on single-patch geometries, while multipatch support is identified as necessary for more complex simulations.Maintaining user-friendly features and distributed-memory capabilities makes this extension nontrivial and leaves it for future work.

Appendix A. Higher-order spatial derivatives

The appendix describes the isoparametric mapping and derives spatial derivatives of its inverse through matrix inversion and repeated differentiation. These derivatives are then used to obtain first-, second-, and third-order basis-function derivatives.

  • Geometric mapping: The geometric mapping uses the isoparametric concept, with control-point locations defining the physical-space mapping.
  • Inverse mapping: First spatial derivatives of the inverse mapping are computed by matrix inversion through an identity involving the Kronecker delta.
  • Inverse mapping: Second and third spatial derivatives of the inverse mapping follow by repeated differentiation.
  • Basis derivatives: First-, second-, and third-order spatial derivatives of the basis functions are obtained from the mapping derivatives.

Appendix B. Tutorial for the Bratu example

The tutorial walks through building and running a parallel PetIGA Bratu example, from environment setup and compilation to discretization inspection, convergence monitoring, and solution visualization.

  • The tutorial assumes PETSc and PetIGA are configured with MPI in a POSIX environment for parallel execution.
  • The example uses a simple makefile that includes PetIGA's variables and build rules to generate the Bratu executable.
  • A four-process run solves a two-dimensional, single-component problem on a 128×128 grid using quadratic C1 B-spline basis functions and three-point Gauss–Legendre quadrature.
  • The SNES function norm decreases from 5.266384548611e-02 to 2.309601078808e-11 over eight monitored iterations.
  • Runtime visualization provides a quick solution check, while VTK output supports higher-quality post-processing in ParaView.

Appendix C. Numerical differentiation comparison

This section benchmarks three numerical differentiation approaches for assembling the global Jacobian of a three-dimensional Bratu problem. Colored finite differences is over an order of magnitude slower than the explicitly coded Jacobian baseline, while local numerical differentiation is about four times slower.

  • Benchmark setup: The benchmark computes the global Jacobian matrix for a computationally challenging three-dimensional Bratu problem.The tests modify a prior Bratu code example to compute the Jacobian rather than solve the nonlinear problem.
  • Compared approaches: Explicit uses the explicitly coded, user-provided Jacobian routine.
  • Compared approaches: Coloring uses PETSc's global colored finite differences, while Local ND uses PetIGA's local numerical differentiation at quadrature points.
  • Benchmark setup: The benchmarks ran on a workstation with two 10-core Intel Xeon E5-2680 v2 CPUs, 128 GB memory, Linux 4.0.4, and GCC 4.9.2 with -Ofast.
  • Results: Over an order of magnitude slower colored finite differences and about four times slower local numerical differentiation were measured against the explicitly coded Jacobian baseline.The reported results summarize minimum wall-clock times from five successive global-Jacobian computations; parameters include processes, elements, polynomial degree, and continuity order.
Loading 1305.4452v3…