How to begin writing scientific codes in C++ in Trilinos or PETSC style?

Question

My background: I have taken some courses on numerical analysis during my PhD and read a few books on the topic as well. I mostly work on low Reynolds number fluid mechanics and use boundary element method for solving Stokes equations on unstructured meshes. I mostly read the mathematical equations and write my own codes in Fortran 90.
I seldom use large linear algebra packages. If I need something, then I write my own code. I have realised this is quite unsustainable and unscalable. I need to start writing codes on top of optimised linear algebra packages like Trilinos or PETSc. It appears foolish not to take advantage of these optimised packages.
However, since I don't have a background in C++ or C, I am finding it a bit difficult to understand the philosophy behind how these large codes are written. For examples, how classes and templates are made. I am finding it difficult to just read the codes in PETSc or Trilinos and understand why they are written/organised in a certain way. At the moment, I am reading the book, 'Parallel Scientific Computing in C++ and MPI', by Karniadakis and Kirby to get some ideas but I need to see examples that have good amount of description.
Please provide me with some suggestions or ideas that teach how to start writing or even thinking about writing such large software libraries. Any books that you particularly recommend or video lectures?
Considering my experience in FORTRAN, I am also not sure how I will benefit from using C++ classes. Is it just for better organisation of the codes? How do I start thinking of implementing mathematical objects in C++?

Victor Eijkhout · Answer

Ed Bueler just wrote the perfect book for you: "PETSc for Partial Differential Equations: Numerical Solutions in C and Python"
https://epubs.siam.org/isbn/978-1-61197-630-4

Nachiket · Answer

Yes, I think "better organisation of the codes" as you say is an important reason to use C++ and more generally Object Oriented Programming (OOP). Due to better organization of code, OOP can reuse code more effectively. Reusing code means less code overall means fewer bugs to fix.
Another reason, in my opinion, is OOP allows us to express things more clearly. If A,B,C are matrices what is more clearer:
for (ii = 0; ii < n; ii++){
    for (jj = 0; jj < n; jj++){
      C[ii][jj] = A[ii][jj] + B[ii][jj]
  }
}

or
C = A + B

Certainly, the latter is clearer. By supporting operator overloading OOP allows us to write the second code which is clearer. Another example would be assembling element rhs vectors into global rhs vectors.
OOP allows us to do
GlobalRhsVector += ElementRhsVector
All the baggage of converting local indices to global indices is taken care of when we define how the "+" operator behaves.

VoB · Answer

Just too long for a comment
I'm a beginner with MPI too, and I'm using that book also. But it seems to me just devoted to the MPI approach. It just give a brief introduction to classes with an example or two, and it does not even talk about templates, which is something that you find in most scientific codes.
This is what I'm doing now to learn:

I started with the classic "Principle and practice ..." by Stroustrup, which is great for a beginner as I still am, in order to understand the most important structures.

MPI for the parallel computing part

Numerical Recipes, which has a third edition in C++. I think that it is the most important C++ book  about numerical algorithms.

Take this with a grain of salt, as I'm just a beginner in the C++ scientific programming, and someone will surely have better suggestions than me.

How to begin writing scientific codes in C++ in Trilinos or PETSC style?

3 Answers

Add your own answers!

Ask a Question