MPI vs OPENMP usage for boundary element method

Question

I have an implementation question with MPI and OPENMP.

I have a large three-dimensional surface mesh divided into elements. I want to loop over each element (outer loop) and compute an integral over the whole surface (inner loop). There are two loops in total. I wish to parallelise this operation.

This kind of operation occurs in boundary element methods. For example, such an operation is needed to solve the electrostatic potential on an arbitrary surface. The first loop goes through every evaluation points $\boldsymbol{x}_0$ at the centre of each element and the second loop performs the integration for this particular $\boldsymbol{x}_0$:

$\phi (\boldsymbol{x}_0) = -\int\limits_S G(\boldsymbol{x}_0, \boldsymbol{x}) [\boldsymbol{n}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla} \phi (\boldsymbol{x})] + \phi (\boldsymbol{x})[\boldsymbol{n}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla}G(\boldsymbol{x}_0, \boldsymbol{x}) ]~ \mathrm{d}S$

where $G(\boldsymbol{x}_0, \boldsymbol{x})=1/4\pi|\boldsymbol{x}_0-\boldsymbol{x}|$ is the free-space Green's function for Laplace equation and $\phi$ is the unknown potential that we want to solve for some given boundary conditions.

MPI: Here I have two options.

1. I can generate the mesh in one processor and broadcast it to all the others. Then I can divide the integration work amongst the processors. Let's say I have N_elm elements and n_procs processors, then each processor be responsible for integrating (N_elm/n_procs) elements. (Question) While this will not require any communication when doing the integration, will this make n_procs copies of the mesh in the memory?

2. I generate the mesh in one processor, divide it into (N_elm/n_procs) elements, so that each processor only has a part of the domain. Each processor will be able to loop over its own element and integrate locally over its part of the mesh. However, for integrating over parts of mesh it doesn't have communication will be required. So, I can create a temporary array for each processor of size [(N_elm/n_procs) elements * 3] for this purpose. In total, each processor will then have created 2*[(N_elm/n_procs) elements * 3], one for its own mesh and a second temporary one which gets replaced.

(Question) Is 1 or 2 preferable, or there is a 3rd option?

OPENMP: Here the implementation is easier. (Question) What happens under the hood when I create a mesh and parallelise the for-loop? Are the threads given multiple copies of the entire mesh similar to number (1) above or the threads share the mesh between each other, more like number (2) above.

Answer 1

MPI Q1: In a naive implementation yes, you will end up with a copy of the mesh on each process. You can reduce this somewhat by only duplicating the parts of the mesh which contain the surface, since you are only doing a surface integral and not a volume integral. This in theory will reduce how much memory is duplicated per process ($O(n^2)$ vs. $O(n^3)$). I suppose pseudo-code for this might look something like this:

on each process:
  allocate space for elements rank * nelems / nprocs to (rank + 1) * nelems / nprocs
  allocate space for the surface of the mesh
  initialize

  for each local element n:
    get x0_n
    for each face on the surface:
      compute the surface integral using x0_n

MPI Q2: If you have enough memory available and don't have too many faces on the surface, then I think method 1 would be better, especially if you decide to implement the slight modification to method 1 I proposed. This method has the least inter-process communcations (basically none), and usually the number of boundary surface faces is much smaller than the number of elements in the domain. However, this number may still be quite large, and doesn't allow any parallelization beyond the granularity of a total surface integral given x0_n.

You can add a second layer of parallelization where the surface integral itself is divided up (for example via OpenMP or an MPI sub-communicator) and computed given a single value of x0_n and then do something equivalent to MPI_Reduce to combine the results for that surface integral, though this may or may not actually help.

If you're worried about running out of memory, then you can try to allocate the surface using something like MPI_Win_allocate_shared or other MPI window function, though I think there will be some communication overhead with this approach (hopefully not too much).

OpenMP Q1: Nominally no. The mesh is allocated once per process, and all threads on that process have shared memory access to the mesh.

Note that there is a caveat that OpenMP is starting to support hardware accelerators, so there might be some situations where the mesh will be duplicated in main memory and on the accelerator's memory.