# PetSc vs Sundials for serial numerical computations?

I am currently working on a physics problem that turns into a non-linear boundary value problem. I need an efficient numerical solver that I could run on my laptop with i5 dual core CPU. I am discretizing my 15 equations on an N x N x N cubic grid using fourth order finite difference derivatives. I am interested in using Newton-Krylov methods since the Jacobian that is generated is sparse. My ultimate goal is to solve this system on a 128 x 128 x 128 cube, which means I will be solving almost 15 million equations for the same number of variables. I am already in the process of learning SUNDIALS packages but I recently found the PETSc library, which many people seem to be praising for quality of the code and efficiency. PETSc also seems to have a large number of preconditioners needed for Krylov methods.

So I would like to ask how do both packages compare for serial computations since I would like to run the solver on my laptop (at least for testing purposes on a bit coarser grid). I might eventually move to my university's HPC cluster to solve the system on a finer grid. Also, my intention is to use these packages on a long term basis.

I am new to high performance computing so kindly excuse my ignorance. I would really appreciate any suggestions.

P.S- The equations I am working on come from Non-Abelian gauge theories and are in Tensor form. Individual components (15) equations are long and complicated...so I actually compute these equations using mathematica and then discretize them. These equations are similar to Maxwell's equations but more complicated and non-linear in the fields. I do not know how much this would help here but the equations are as follows

$(D_\mu F^{\mu \nu})^a = g \epsilon_{abc} (D_\nu \phi)^b \phi^c$, and

$D_\mu (D^\mu \phi)^a = - \lambda (\phi^b \phi^b -v^2) \phi^a$.

Here Einstein's summation convention is implied. The indices $\mu$ and $\nu$ run from 0 to 3. And index $a$ runs from 1 to 3. $D$ represents gauge covariant derivatives and $\phi$ and $F$ represent Higgs and gauge fields strength. You can check this article for more details....page 9. http://www-thphys.physics.ox.ac.uk/people/MaximeGabella/higgs.pdf

• Welcome to SciComp.SE. Can you be more specific about your problem? What are you solving for? Maybe, add the equations that you want to solve. Those things might help people to understand, and it would be easier for them to help. – nicoguaro Jun 12 '16 at 19:13
• Well, the equations I am working on come from Non-Abelian gauge theories and are in Tensor form. Individual components (15) equations are long and complicated...so I actually compute these equations using mathematica. These equations are similar to Maxwell's equations but more complicated and non-linear in the fields. You can check the equations on this article link below on page 9. The equations are highlighted in the box. But I do not know how useful this would be here as people here might not be familiar with this field... www-thphys.physics.ox.ac.uk/people/MaximeGabella/higgs.pdf – singularity Jun 12 '16 at 19:42
• You can write your equations in indicial notation, that way it won't be that messy. I don't know if they would help either, or if somebody is familiar with the topic. But it won't hurt, and you can let them decide. – nicoguaro Jun 12 '16 at 19:47

## 2 Answers

In general, you can do more with PETSc.

SUNDIALS is a collection of ODE solvers (in CVODE, Adams-Bashforth and BDF methods; in ARKODE, ARKIMEX methods) and DAE solvers (IDAS implements a BDF method) with sensitivity capabilities (the CVODES and IDAS variants), and a nonlinear solver (KINSOL). There are a few Krylov solvers in there (at least GMRES and CG) and a couple other fixed-point type methods, like Anderson iteration, plus interfaces to a few solvers, like KLU, and SuperLU.

PETSc has Krylov subspace methods, nonlinear solvers, ODE and DAE solvers, linear and nonlinear preconditioners (including a multigrid method), mesh management features, and interfaces to external software that provide more of these types of features, including a CVODE interface to SUNDIALS that uses GMRES.

The only things in SUNDIALS that aren't in PETSc right now (and may even be in newer versions of PETSc) are BDF methods, and BDF-based sensitivity calculations. PETSc has an implementation for sensitivity computations, may use different algorithms for computing sensitivities.

PETSc, on the other hand, has loads of stuff that SUNDIALS doesn't have, like more linear solvers, more preconditioners (as far as I remember, SUNDIALS doesn't have any native preconditioners), nonlinear preconditioning, mesh management, etc.

SUNDIALS is easier to learn, but PETSc will probably pay off more for you in the long run.

As mentioned in Geoff Oxberry's more complete answer, it should be noted that PETSc includes TSSUNDIALS, an interface to SUNDIALS.

If you configure PETSc with the --download-sundials option (see python2 ./configure --help | grep -A 2 sundials for other related options, such as using an existing SUNDIALS library), then you can use at least some of the SUNDIALS functionality from PETSc and thus might have an easier time making comparisons between SUNDIALS and PETSc's other TS implementations.

• I pointed out that PETSc has an interface to CVODE (this solver is what TSSUNDIALS interfaces to) in my answer, although I didn't mention it by the data type. – Geoff Oxberry Jun 14 '16 at 18:50
• Sorry that I missed that! – Patrick Sanan Jun 15 '16 at 6:48