Implementation of generic solving capabilities in ParaFEM

Key Personnel

PI/Co-I: Pankaj Pankaj, Francesc Levrero Florencio, Lee Margetts

Technical: Mark Filipiak

Relevant documents

eCSE Technical Report: Implementation of generic solving capabilities in ParaFEM

Project summary

In this eCSE project, the PETSc solver library has been interfaced to the ParaFEM finite element analysis (FEA) library.

FEA is widely used in engineering to simulate the deformation of structures under load. In FEA, the structure is discretised into small elements (often millions of these), each connected to its neighbours. This leads to a large system of equations to solve and at the heart of every FEA program is an equation solver. The orthopaedic engineering research group at the University of Edinburgh uses the ParaFEM FEA library to simulate the mechanical behaviour of trabecular bone, using a detailed geometry generated through X-ray microtomography in order to assess the response of bone to mechanical loads at the microstructural level (Figure 1).

Figure 1- A finite element mesh of a trabecular bone specimen.

ParaFEM (www.parafem.org.uk) is a freely available, portable library of subroutines for parallel FEA and can be used to solve various types of problems, including stress analysis, heat flow, eigenvalue and forced vibrations. ParaFEM has two different equation solvers (CG and BiCGStab). These work well but there are a range of other solvers that could be used and may be faster for certain problems (including the bone example above). PETSc (https://www.mcs.anl.gov/petsc) is an open source linear solver library with a comprehensive set of scalable parallel solvers and preconditioners (including parallel direct solvers via the external packages MUMPS and SuperLU).

The bulk of the work in the eCSE project was developing an interface to PETSc in ParaFEM. With this interface, a user can easily convert a ParaFEM driver program to use PETSc and can choose different solvers at run time (and change the solver used during the run).

A comparison of similar solvers provided by ParaFEM and PETSc shows that the relative performance depends on the problem. For example, for a medium-sized bone simulation (~25 million elements, corresponding to a 1 cm³ bone sample), the PETSc solver is 10% slower than the corresponding ParaFEM CG solver and requires about three times as much memory (which is expected as the solvers in ParaFEM are based on an element by element approach, i.e. the global stiffness matrix is never assembled). The main benefit of adding the PETSc solvers will be seen in future applications which will be able to use any of the large range of solvers now available. For instance, trabecular bone suffers from damage and may undergo buckling. Such systems (with damage or non-associative plasticity) have an unsymmetric stiffness matrix; systems after buckling, or which are undergoing softening, have an indefinite stiffness matrix. With the range of available linear solvers in PETSc, these types of linear algebraic systems can now be solved.

Achievement of objectives

1 Integration of iterative/direct external solvers in ParaFEM. The success metric is that the solver is implemented and tested.

An interface to the PETSc solver library has been added to ParaFEM and the results and performance of the ParaFEM and PETSc conjugate gradient/Jacobi (CG/Jacobi) solvers have been compared.

2 Assess the performance and scaling of the solver. As a success metric, we expect a decrease in linear algebraic solving time of at least 10% for problems of size ranging from medium to large (2.5M to 100M degrees of freedom, DOF from now onwards). Since PETSc has direct solving capabilities, we expect that the decrease in solving time will be more than 25% for small problems (<2.5M DOF), since direct solvers are well known to deal more efficiently with small linear algebraic systems. Additionally, different preconditioners will be tested and their performance benchmarked.

For a medium linear elasticity problem (~12M degrees-of-freedom), the PETSc CG/Jacobi solver was 30% faster than ParaFEM, but for a large nonlinear elasto-plastic analysis of trabecular bone (~26M DOF), PETSc was 10% slower than ParaFEM at large process counts. In all cases the peak memory needed by PETSc was three times larger than that needed by ParaFEM. Direct solvers were slower than iterative solvers for the small case tested (~35K DOF). Other preconditioners were also slower overall, including setup time.

3 Addition of a general interface to easily swap between solvers. We expect to have a solver which is able to deal with all kinds of linear algebraic systems that arise in FE systems - positive-definite, indefinite and unsymmetric. The user will have the choice to use CG for positive-definite linear algebraic systems, BiCGSTAB(l) for unsymmetric linear algebraic systems (these two from ParaFEM) or PETSc solvers (direct or iterative) for any possible system. We would like to let the user choose which solver is the most appropriate for the problem at hand.

The user has a choice between ParaFEM solvers (CG and BiCGStab, both with Jacobi preconditioner) and PETSc solvers. The choice of PETSc solver is done using a control file in PETSc format - this allows the user to choose the solver and preconditioner at run time.

4 Updating the ParaFEM documentation. The ParaFEM documentation will be updated during the project.

Source code documentation for the application programmer's interface (API) was added in ROBODOC format. A User Guide was produced in the format of the "Programming the Finite Element Method" text book, so that the text can be easily integrated into the relevant chapters and appendices. No online user tutorial was produced.

Summary of the Software

It is recommended to first register on the ParaFEM website http://parafem.org.uk. The ParaFEM code is available on SourceForge in the Subversion repository https://sourceforge.net/projects/parafem. The ParaFEM-PETSc interface is in the petsc branch - this will eventually be merged into the trunk. The current stable revision of the petsc branch is 2237. The software is easy to install on ARCHER - an install script and several job scripts for ARCHER are included in the repository.

Scientific Benefits

Several research groups have shown interest in using the nonlinear elasto-plastic analysis program (xx15), which can now choose from PETSc solvers at the various stages of the simulations. The University of Sheffield Insigneo Institute for in Silico Medicine is interested in using ParaFEM in their workflow and they will collaborate with ParaFEM developers (this is within the European project CompBioMed). Groups at the universities of Bristol and Manchester are also interested. ParaFEM solid mechanics simulations are being coupled to OpenFOAM fluid simulations in a PhD project at the University of Manchester. A research group in Australia wants to use the xx15 program to simulate juice extraction from sugar cane to optimise the milling process.