Solving a memory leak in Communicator_mpi

lib/cartesian/Cartesian_base.h

@@ -50,6 +50,9 @@ public:
 
     GridBase(const std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
 
+    virtual ~GridBase() = default;
+
+
     // Physics Grid information.
     std::vector<int> _simd_layout;// Which dimensions get relayed out over simd lanes.
     std::vector<int> _fdimensions;// (full) Global dimensions of array prior to cb removal

lib/cartesian/Cartesian_full.h

@@ -93,6 +93,7 @@ public:
 
         // Use a reduced simd grid
         _ldimensions[d] = _gdimensions[d] / _processors[d]; //local dimensions
+        //std::cout << _ldimensions[d] << "  " << _gdimensions[d] << "  " << _processors[d] << std::endl;
         assert(_ldimensions[d] * _processors[d] == _gdimensions[d]);
 
         _rdimensions[d] = _ldimensions[d] / _simd_layout[d]; //overdecomposition
@@ -137,6 +138,8 @@ public:
         block = block * _rdimensions[d];
       }
     };
+
+    virtual ~GridCartesian() = default;
 };
 }
 #endif

lib/communicator/Communicator_base.h

@@ -152,6 +152,7 @@ class CartesianCommunicator {
   // Constructor of any given grid
   ////////////////////////////////////////////////
   CartesianCommunicator(const std::vector<int> &pdimensions_in);
+  virtual ~CartesianCommunicator();
   
   ////////////////////////////////////////////////////////////////////////////////////////
   // Wraps MPI_Cart routines, or implements equivalent on other impls

lib/communicator/Communicator_mpi.cc

@@ -75,6 +75,13 @@ CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
   
   assert(Size==_Nprocessors);
 }
+
+CartesianCommunicator::~CartesianCommunicator(){
+  if (communicator && !MPI::Is_finalized())
+    MPI_Comm_free(&communicator);
+}
+
+
 void CartesianCommunicator::GlobalSum(uint32_t &u){
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator);
   assert(ierr==0);

lib/qcd/action/scalar/ScalarInteractionAction.h

@@ -30,22 +30,23 @@ directory
 #ifndef SCALAR_INT_ACTION_H
 #define SCALAR_INT_ACTION_H
 
-
 // Note: this action can completely absorb the ScalarAction for real float fields
 // use the scalarObjs to generalise the structure
 
-namespace Grid {
+namespace Grid
+{
 // FIXME drop the QCD namespace everywhere here
 
 template <class Impl, int Ndim>
-  class ScalarInteractionAction : public QCD::Action<typename Impl::Field> {
+class ScalarInteractionAction : public QCD::Action<typename Impl::Field>
+{
 public:
   INHERIT_FIELD_TYPES(Impl);
+
 private:
   RealD mass_square;
   RealD lambda;
 
-
   typedef typename Field::vector_object vobj;
   typedef CartesianStencil<vobj, vobj> Stencil;
 
@@ -54,17 +55,20 @@ namespace Grid {
   std::vector<int> directions;    //    = {0,1,2,3,0,1,2,3};  // forcing 4 dimensions
   std::vector<int> displacements; //  = {1,1,1,1, -1,-1,-1,-1};
 
-
 public:
-
+  ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l), displacements(2 * Ndim, 0), directions(2 * Ndim, 0)
-    ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l), displacements(2*Ndim,0), directions(2*Ndim,0){
+  {
-      for (int mu = 0 ; mu < Ndim; mu++){
+    for (int mu = 0; mu < Ndim; mu++)
-		directions[mu]         = mu; directions[mu+Ndim]    = mu;
+    {
-		displacements[mu]      =  1; displacements[mu+Ndim] = -1;
+      directions[mu] = mu;
+      directions[mu + Ndim] = mu;
+      displacements[mu] = 1;
+      displacements[mu + Ndim] = -1;
     }
   }
 
-    virtual std::string LogParameters() {
+  virtual std::string LogParameters()
+  {
     std::stringstream sstream;
     sstream << GridLogMessage << "[ScalarAction] lambda      : " << lambda << std::endl;
     sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl;
@@ -75,16 +79,19 @@ namespace Grid {
 
   virtual void refresh(const Field &U, GridParallelRNG &pRNG) {}
 
-    virtual RealD S(const Field &p) {
+  virtual RealD S(const Field &p)
+  {
     assert(p._grid->Nd() == Ndim);
     static Stencil phiStencil(p._grid, npoint, 0, directions, displacements);
     phiStencil.HaloExchange(p, compressor);
     Field action(p._grid), pshift(p._grid), phisquared(p._grid);
     phisquared = p * p;
     action = (2.0 * Ndim + mass_square) * phisquared - lambda / 24. * phisquared * phisquared;
-      for (int mu = 0; mu < Ndim; mu++) {
+    for (int mu = 0; mu < Ndim; mu++)
+    {
       //  pshift = Cshift(p, mu, +1);  // not efficient, implement with stencils
-	parallel_for (int i = 0; i < p._grid->oSites(); i++) {
+      parallel_for(int i = 0; i < p._grid->oSites(); i++)
+      {
         int permute_type;
         StencilEntry *SE;
         vobj temp2;
@@ -92,15 +99,21 @@ namespace Grid {
 
         SE = phiStencil.GetEntry(permute_type, mu, i);
         t_p = &p._odata[i];
-	  if ( SE->_is_local ) {
+        if (SE->_is_local)
+        {
           temp = &p._odata[SE->_offset];
-	    if ( SE->_permute ) {
+          if (SE->_permute)
+          {
             permute(temp2, *temp, permute_type);
             action._odata[i] -= temp2 * (*t_p) + (*t_p) * temp2;
-	    } else {
+          }
+          else
+          {
             action._odata[i] -= (*temp) * (*t_p) + (*t_p) * (*temp);
           }
-	  } else {
+        }
+        else
+        {
           action._odata[i] -= phiStencil.CommBuf()[SE->_offset] * (*t_p) + (*t_p) * phiStencil.CommBuf()[SE->_offset];
         }
       }
@@ -111,7 +124,8 @@ namespace Grid {
     return -(TensorRemove(sum(trace(action)))).real();
   };
 
-    virtual void deriv(const Field &p, Field &force) {
+  virtual void deriv(const Field &p, Field &force)
+  {
     assert(p._grid->Nd() == Ndim);
     force = (2.0 * Ndim + mass_square) * p - lambda / 12. * p * p * p;
     // move this outside
@@ -119,23 +133,31 @@ namespace Grid {
     phiStencil.HaloExchange(p, compressor);
 
     //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
-      for (int point = 0; point < npoint; point++) {
+    for (int point = 0; point < npoint; point++)
-	parallel_for (int i = 0; i < p._grid->oSites(); i++) {
+    {
+      parallel_for(int i = 0; i < p._grid->oSites(); i++)
+      {
         const vobj *temp;
         vobj temp2;
         int permute_type;
         StencilEntry *SE;
         SE = phiStencil.GetEntry(permute_type, point, i);
 
-	  if ( SE->_is_local ) {
+        if (SE->_is_local)
+        {
           temp = &p._odata[SE->_offset];
-	    if ( SE->_permute ) {
+          if (SE->_permute)
+          {
             permute(temp2, *temp, permute_type);
             force._odata[i] -= temp2;
-	    } else {
+          }
+          else
+          {
             force._odata[i] -= *temp;
           }
-	  } else {
+        }
+        else
+        {
           force._odata[i] -= phiStencil.CommBuf()[SE->_offset];
         }
       }