Local submesh integral

Question

Dear all,

This might be a weird question, but interesting at least to me. Suppose I am computing an integral double integral = assemble(integral_form);, now, is it possible to know the integral on each process? This means, a process p with a certain domain partition assigned, could compute the integral only on its mesh_p subdomain.

This operation should, in theory, be without any communication, since I will sum or reduce the results on each process to get the correct answer, but I really don't know if this is even possible.

Thanks!

Answer 1

There are many options how to do this:

1. Use Scalar with MPI_COMM_SELF. Example at the bottom.

2. Subclass Scalar or GenericTensor class so that it does not collect the value across ranks. Use void assemble(GenericTensor& A, const Form& a);. To do the former, just take Scalar.cpp and strip off MPI communication.

3. Use

   void assemble_[cells|exterior_facets|interior_facets](
            GenericTensor& A,
            const Form& a, UFC& ufc,
            std::shared_ptr<const MeshFunction<std::size_t> > domains,
            std::vector<double>* values);
   

   Collect values on each rank separetely. Some extra memory will be needed.

4. Define form using measure with which is local to each rank and assemble separately on each rank. This will not get rid of communication and will require some extra memory to store a cell/facet function.

5. Prepare local submeshes on MPI_COMM_SELF and assemble there. Lot of extra memory needed.

EDIT Example implementation of 1: Functional.ufl

element = FiniteElement("Lagrange", triangle, 1)
f = Coefficient(element)
I = f*dx

forms = [I]

main.cpp

#include <dolfin.h>
#include "Functional.h"

using namespace dolfin;

class Source : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = MPI::rank(MPI_COMM_WORLD);
  }
};

int main()
{
  UnitSquareMesh mesh(32, 32);
  Source f;
  Functional::Form_I I(mesh, f);

  TensorLayout tl(MPI_COMM_SELF, {}, 0, 0, {}, false);
  Scalar local_value;
  local_value.init(tl);
  assemble(local_value, I);
  info("Rank %d, integral %f", MPI::rank(MPI_COMM_WORLD),
       local_value.get_scalar_value());

  return 0;
}

Comments

Of course. It works with Fenics MPI:

C++

#include <mpi.h>
#include <stdio.h>

int main(int argc, char *argv[])
{
  int rank, ll, gg;

  MPI_Init(&argc, &argv);

  MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  ll = rank;
  gg = 0;

  MPI_Reduce(&ll, &gg, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);

  if (rank == 0)
  {
      printf("MPI_Reduce gg = %d \n",gg);
  }

  MPI_Finalize();
}

OUTPUT

sensei@Senseis-MBP:Public$ which mpicc
/Applications/FEniCS.app/Contents/Resources/bin/mpicc

sensei@Senseis-MBP:Public$ mpicc a.c

sensei@Senseis-MBP:Public$ which mpirun
/Applications/FEniCS.app/Contents/Resources/bin/mpirun

sensei@Senseis-MBP:Public$ mpirun -np 6 ./a.out 
MPI_Reduce gg = 15 

sensei@Senseis-MBP:Public$ mpirun -np 10 ./a.out 
MPI_Reduce gg = 45 

---

Right! Try MPI_Allreduce. dolfin::MPI::sum is implemented using MPI_Allreduce.

---

It seems to work MPI_Allreduce in Fenics, however, when compiling in Xcode it doesn't, and that's very weird.

The compiler seems the same, but I might overlook something here:

sensei@Senseis-MBP:Public$ which mpicc
/Applications/FEniCS.app/Contents/Resources/bin/mpicc

sensei@Senseis-MBP:Public$ mpicc --showme
/usr/bin/gcc -I/Applications/FEniCS.app/Contents/Resources/include -L/Applications/FEniCS.app/Contents/Resources/lib -lmpi -lm

sensei@Senseis-MBP:Public$ /usr/bin/gcc --version
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 6.1.0 (clang-602.0.53) (based on LLVM 3.6.0svn)
Target: x86_64-apple-darwin14.4.0
Thread model: posix

I will make a CMake project just to be sure.

---

Well, I don't know what's the difference between

seems to work MPI_Allreduce in Fenics

and

in Xcode it doesn't

But yes, there could be linking problems causing this.

---

Solved. It seems I added libmpistubs.dylib to the linking phase, and this caused the problem. Weird enough, MPI_Reduce wasn't affected.

PS. "In Fenics" means from the shell, in Xcode means running inside Xcode.