I am currently developing a domain decomposition method for the solution of the scattering problem. Basically I am solving a system of Helmholtz BVPs iteratively. I discretize the equations using finite element method over triangular or tetrahedral meshes. I'm developing the code towards my Phd thesis. I am aware of some of the existing finite element libraries out there such as deal.ii or DUNE and though I think they are great, with inspirational design and API, for learning purposes I wanted to develop my own little application from scratch.

I am at a point where I have my serial versions running and now I want to parallelize them. After all, it is one of the strengths of the domain decomposition framework to formulate algorithms that are easy to parallelize, at least in principle. In practice however, there are many details one must consider. Mesh management is one of them. If the applications is to achieve high resolution while scaling well to many CPUs the replication of an entire mesh on every CPU is inefficient.

I wanted to ask those developers who work on similar applications in high performance computing environments how they deal with this issue.

There is p4est library for distributed mesh management. I do not need AMR so it might be an overkill since I'm only interested in using uniform meshes and I'm not sure if it can refine triangular meshes. I could also simply create a uniform mesh then feed it into one of the mesh partitioners and do some post processing of the output.

The simplest approach seems to create a separate file for each partition containing mesh information relevant only to that particular partition. This file would be read by a single CPU which would be responsible for assembly of the discrete system on that portion of the mesh. Of course, some global partition connectivity/neighborhood information would also need to be stored in a file read by all CPUs for inter process communication.

What other approaches are out there? If some of you could share, what are some of the commonly used methodologies in the industry, or government research institutions related to handling this issue? I am quite new to programming a parallel finite element solver and I wanted to get a feel for whether or not I'm thinking about this problem correctly and how others are approaching it. Any advice or pointers to relevant research articles would be greatly appreciated!

Thanks in advance!

  • $\begingroup$ If you are looking for mesh partitioner - METIS would be good choice. Check also ParMETIS. Managing meshes is different story, ITAPS iMesh can be solution for you. Please check also answers for my question here: scicomp.stackexchange.com/questions/4750/… $\endgroup$ Commented Jun 12, 2013 at 5:44
  • $\begingroup$ @KrzysztofBzowski: have you used Scotch library perhaps as well? I was wondering what is the difference between Scotch and Metis when it comes to finite elements. The iMesh project seems very interesting. I will read more about it in the next few days. I know about deal.II, and DUNE. I remember looking into openMesh some time ago but figured that it would be easier to implement the functionality I needed from scratch. For sequential meshes, basically I adapted the half edge/face data structure presented in this paper link Thanks! $\endgroup$
    – midurad
    Commented Jun 13, 2013 at 4:21

3 Answers 3


If you are not using AMR and do not want to scale beyond 1K-4K cores then simply do this.

  1. Rank 0 reads the entire mesh and partitions it using METIS/Scotch etc. (Note: This is a serial operation).

  2. Rank 0 broadcasts the element/node partitioning info to all other ranks and frees the memory (used to store the mesh)

  3. All ranks read the nodes/elements they own (including ghost nodes) from the same input file (Note: 2000 ranks accessing the same input file might sound slow but is not in practice, though it may be bad for the file system but then we are doing it only once).

  4. All ranks need to create the local to global node/element/dof mappings for application of BCs and assembling of matrices and renumber the nodes.

After everything is said and done all data on a rank will be local so you should be able to scale well (memory wise). I do all this in about 100 lines (see lines 35-132 here) in a small code of mine.

Now if your mesh is too large (e.g., >100-250 million elements) that you cannot partition it using METIS on a single node and need ParMETIS/PT-Scotch then you have the additional work of partitioning it in parallel before all cores/ranks can read it. In such a scenario it might be easier to keep the partitioning phase separate from the main code for logistical reasons.

Btw AMR libs usually dont do tets. Also PETSc is good choice for parallelization of your code.

Edit: Also see here and here.

  • $\begingroup$ Thanks for sharing your code! I will most likely take the route you have outlined above. It seems the least complicated and I already have an idea of how to go about it. In addition it will be a good exercise in MPI programming. You mentioned that AMR libs usually don't handle tets. Would it be because a naive refinement on say a quad-tree of a triangular meshes could lead to bad quality mesh? Refining quads seems easy but splitting a tet into four seems difficult if one wants to preserve quality. Is there a C++ wrapper for PETSc perhaps? I can use C, but C++ would be better. $\endgroup$
    – midurad
    Commented Jun 13, 2013 at 4:28
  • $\begingroup$ @midurad If you prefer C++ over C, you should consider Trilinos, which is a C++ library comparable to PETSc. Moreover, Trilinos has a package (Zoltan) which you can use for performing mesh partitioning. $\endgroup$
    – Dr_Sam
    Commented Jun 13, 2013 at 7:34
  • $\begingroup$ @midurad You only need very few MPI calls if you use PETSc. Refining tets should be easy but dealing (efficiently) with the associated dynamic data structures might require some thought and work. You should be able to use PETSc with C++ but given your requirements libmesh might be a viable option (I think it supports AMR and tets). $\endgroup$
    – stali
    Commented Jun 13, 2013 at 14:47
  • $\begingroup$ Thank you all for the information. That was very helpful. $\endgroup$
    – midurad
    Commented Jun 14, 2013 at 4:16

This may not come as a surprise to you given that I develop deal.II, but here's my perspective: When I talk to students, I typically tell them to develop their own prototype in the beginning so they can see how it's done. But then, once they've got something small running, I make them use a library that allows them to go so much further because they don't have to reinvent the wheel with basically every step they take.

In your case, you've already seen how to implement a simple Helmholtz solver. But you'll spend the next 6 months writing the code necessary to do it in parallel, you'll spend another 3 months if you want to use more complicated geometries. You'll then spend 6 more months if you want an efficient solver. And all of this time you're writing code that's already been written by someone else and that, in a sense, doesn't get you any closer to what you actually need to do for your PhD: develop something new that hasn't been done before. If you go down this road, you'll spend 2-3 years of your PhD time re-doing what others have done, and maybe 1 year doing something new.

The alternative is that you now spend 6 months learning one of the existing libraries, but after that you'll have 2-3 years where you really do new stuff, things where every other week you can walk into your adviser's office and show him/her something that's truly new, that runs on massively large scales, or is just very cool in other regards. I think you probably see where I'm going with this by now.

  • 3
    $\begingroup$ Honest question since you're clearly an authority on this: who is going to write the next generation of frameworks like deal.ii if no one in the current crop of PhD students tackles problems like this? We're already seeing a problematic trend of incoming PhD students that have never even compiled a program. It's a little disturbing to me that the average code dev skills seems to be on a continuous decline in computational scientists. $\endgroup$
    – Aurelius
    Commented Mar 14, 2014 at 16:22
  • 1
    $\begingroup$ It's a fair question. You need grad students as bone-headed and stubborn as I was :-) But my answer is that just because we probably need a few people who do it, that doesn't mean that we should encourage everyone to spend years of their lives repeating what others have already implemented. $\endgroup$ Commented Mar 14, 2014 at 18:35
  • 2
    $\begingroup$ Yeah, fair enough. IMO, the single biggest thing that has held back the CFD research world for the past 20 years has been a lack of software engineering talent and a rejection of modern software practices by the greybeards. Frameworks aside, so many PhD students are held back by bad legacy code and an inability to quickly construct the complex pieces of software needed for modern numerical methods on modern hardware. $\endgroup$
    – Aurelius
    Commented Mar 14, 2014 at 19:27
  • $\begingroup$ I don't disagree with the statement about the greybeards (though mine is turning grey as well these days...). But they also see that you have to choose between crufty legacy codes or reinventing the wheel when you have a new grad student. Very few people get to enjoy having success with the software they write (present author not withstanding), and you don't want to send a promising grad student down that road if you don't know that they can make a career of it. $\endgroup$ Commented Mar 14, 2014 at 21:29

This is not a complete answer.

For implementation of parallel domain decomposition methods, I encountered some complications. First, one can use many processors for one subdomain or feed one processor with many subdomains and one might want to implement both paradigms. Second, substructured form of domain decomposition methods requires separating out subdomains' (not element's) faces, edges, vertices. I do not think these complications are readily included in the parallel mesh management. The situation becomes simpler if you consider one processor for one subdomain and using the overlapping RAS/RASHO method. Even in this case, you would better manage your parallel layout yourself, because when your subdomains have overlap the usual parallel mesh management may require communications for the entire overlapping regions while the necessary data for RAS/RASHO are only on the boundaries.


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.