2
$\begingroup$

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++?

Improve this question
$\endgroup$
3
  • $\begingroup$ When I'm not mistaken PETSC should come with the possibility of use statements (see documentation and e.g. mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/…). $\endgroup$
    – albert
    Commented Jul 13, 2020 at 16:58
  • $\begingroup$ C++ offers a great deal of abstraction. You can translate mathematical objects into implementations and after that developing and prototyping scientific software becomes easier. I would suggest just going into PETSc documentation and checking their examples. That is practically how I learned about 6 years ago. I would also suggest joining their user mailing list, but it is highly encouraged to first run a search on their mail archives before asking a question because probably the answer is there already. Also it helps if someone around you knows about PETSc, too $\endgroup$
    – Abdullah Ali Sivas
    Commented Jul 13, 2020 at 20:41
  • 1
    $\begingroup$ I learned both C++ and finite element coding by reading the deal.II documentation, which is exhaustive. I believe you can also write boundary element methods in deal.II, see here dealii.org/developer/doxygen/deal.II/Tutorial.html#topic under advanced techniques. $\endgroup$
    – cfdlab
    Commented Jul 14, 2020 at 8:18

3 Answers 3

Reset to default
2
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you very much. I will read that book once I have finished reading yours "Intro to High Performance Scientific Computing". I have just started reading your book. $\endgroup$
    – Debasish Das
    Commented May 31, 2021 at 14:17
  • 1
    $\begingroup$ I also have a book "Parallel Programming" that has a PETSc section. Btw, I hope you find my books through my website pages.tacc.utexas.edu/~eijkhout/istc/istc.html There are lots of these "free pdf" sites that have completely ancient copies of my books. $\endgroup$
    – Victor Eijkhout
    Commented Jun 1, 2021 at 0:02
  • $\begingroup$ ..... and in the meantime my webserver has gone away. The link for my textbooks is now theartofhpc.com $\endgroup$
    – Victor Eijkhout
    Commented May 22, 2022 at 1:42
1
$\begingroup$

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.

Improve this answer
$\endgroup$
0
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Just as an aside in Fortran this is part of the core language $\endgroup$
    – Ian Bush
    Commented Jul 14, 2020 at 8:21
  • $\begingroup$ That is true, Ian. But with PETSc, you can do stuff like A.Mult(x,y) which is $y=Ax$ when A is sparse. You can further overload operators to get y = A * x which is more Fortran-like. I don't think that is possible in Fortran unless sparse matrices are native types in Fortran now. $\endgroup$
    – Abdullah Ali Sivas
    Commented Jul 14, 2020 at 13:23
  • $\begingroup$ Yes - you would have to implement a sparse data type, but then either approach is straightforward $\endgroup$
    – Ian Bush
    Commented Jul 14, 2020 at 16:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.