Changeset 2506 for XIOS2/trunk/extern/remap/src/mapper.cpp

Timestamp:

05/31/23 12:20:29 (13 months ago)

Author:

ymipsl

Message:

Supress limitation of number of vertex max in remaper.
YM

File:

• XIOS2/trunk/extern/remap/src/mapper.cpp (modified) (19 diffs)

XIOS2/trunk/extern/remap/src/mapper.cpp

@@ -145 +145 @@
       nIntersections++;
   }
+
   /* overallocate for NMAX neighbours for each elements */
+/*
   remapMatrix = new double[nIntersections*NMAX];
   srcAddress = new int[nIntersections*NMAX];
@@ -152 +154 @@
   sourceWeightId =new long[nIntersections*NMAX];
   targetWeightId =new long[nIntersections*NMAX];
-
+*/
 
   if (mpiRank == 0 && verbose) cout << "Remapping..." << endl;
@@ -186 +188 @@
 
   int *nbSendElement = new int[mpiSize];
+  int *nbSendVertex = new int[mpiSize];
   int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */
   double **recvValue = new double*[mpiSize];
   double **recvArea = new double*[mpiSize];
   double **recvGivenArea = new double*[mpiSize];
+  int **recvVertexOffset = new int*[mpiSize];
   Coord **recvGrad = new Coord*[mpiSize];
   GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */
@@ -210 +214 @@
     {
       sendElement[rank] = new int[nbSendElement[rank]];
-      recvValue[rank]   = new double[nbSendElement[rank]];
-      recvArea[rank]    = new double[nbSendElement[rank]];
-      recvGivenArea[rank] = new double[nbSendElement[rank]];
-      if (order == 2)
-      {
-        recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)];
-        recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)];
-      }
-      else
-        recvNeighIds[rank] = new GloId[nbSendElement[rank]];
-
       last = -1;
       index = -1;
@@ -235 +228 @@
   /* communicate sizes of source elements to be sent (index lists and later values and gradients) */
   int *nbRecvElement = new int[mpiSize];
+  int *nbRecvVertex = new int[mpiSize];
   MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator);
 
@@ -244 +238 @@
   double **sendArea = new double*[mpiSize];
   double **sendGivenArea = new double*[mpiSize];
+  int    **sendVertexOffset = new int*[mpiSize];
   Coord **sendGrad = new Coord*[mpiSize];
   GloId **sendNeighIds = new GloId*[mpiSize];
-  MPI_Request *sendRequest = new MPI_Request[5*mpiSize];
-  MPI_Request *recvRequest = new MPI_Request[5*mpiSize];
+  MPI_Request *sendRequest = new MPI_Request[6*mpiSize];
+  MPI_Request *recvRequest = new MPI_Request[6*mpiSize];
   for (int rank = 0; rank < mpiSize; rank++)
   {
@@ -259 +254 @@
     {
       recvElement[rank] = new int[nbRecvElement[rank]];
+      MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+  MPI_Status *status = new MPI_Status[6*mpiSize];
+  
+  MPI_Waitall(nbSendRequest, sendRequest, status);
+  MPI_Waitall(nbRecvRequest, recvRequest, status);
+
+  nbSendRequest = 0;
+  nbRecvRequest = 0;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbRecvElement[rank] > 0)
+    {
+      nbRecvVertex[rank]=0 ;
+      for (int j = 0; j < nbRecvElement[rank]; j++) nbRecvVertex[rank]+=sstree.localElements[recvElement[rank][j]].n+1;
+      MPI_Issend(&nbRecvVertex[rank], 1, MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+  
+    if (nbSendElement[rank] > 0)
+    {
+      MPI_Irecv(&nbSendVertex[rank], 1, MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+  
+  MPI_Waitall(nbSendRequest, sendRequest, status);
+  MPI_Waitall(nbRecvRequest, recvRequest, status);
+
+
+  /* for all indices that have been received from requesting ranks: pack values and gradients, then send */
+  nbSendRequest = 0;
+  nbRecvRequest = 0;
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbRecvElement[rank] > 0)
+    {
       sendValue[rank]   = new double[nbRecvElement[rank]];
       sendArea[rank]   = new double[nbRecvElement[rank]];
       sendGivenArea[rank] = new double[nbRecvElement[rank]];
+      sendVertexOffset[rank] = new int[nbRecvElement[rank]];
+
       if (order == 2)
       {
-        sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)];
-        sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)];
-      }
-      else
-      {
-        sendNeighIds[rank] = new GloId[nbRecvElement[rank]];
-      }
-      MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-      nbRecvRequest++;
-    }
-  }
-  MPI_Status *status = new MPI_Status[5*mpiSize];
-  
-  MPI_Waitall(nbSendRequest, sendRequest, status);
-        MPI_Waitall(nbRecvRequest, recvRequest, status);
-
-  /* for all indices that have been received from requesting ranks: pack values and gradients, then send */
-  nbSendRequest = 0;
-  nbRecvRequest = 0;
-  for (int rank = 0; rank < mpiSize; rank++)
-  {
-    if (nbRecvElement[rank] > 0)
-    {
+        sendNeighIds[rank] = new GloId[nbRecvVertex[rank]];
+        sendGrad[rank]    = new Coord[nbRecvVertex[rank]];
+      }
+      else sendNeighIds[rank] = new GloId[nbRecvElement[rank]];
+
       int jj = 0; // jj == j if no weight writing
+      sendVertexOffset[rank][0]=0 ;
       for (int j = 0; j < nbRecvElement[rank]; j++)
       {
@@ -293 +313 @@
         sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area;
         sendGivenArea[rank][j] = sstree.localElements[recvElement[rank][j]].given_area;
+        
+        if (j==0) 
+        {
+          if (order==2) sendVertexOffset[rank][j]=sstree.localElements[recvElement[rank][j]].n+1 ;
+          else sendVertexOffset[rank][j]=1 ;
+        }
+        else 
+        {
+          if (order == 2) sendVertexOffset[rank][j] = sendVertexOffset[rank][j-1]+sstree.localElements[recvElement[rank][j]].n+1 ;
+          else sendVertexOffset[rank][j] = sendVertexOffset[rank][j-1]+1 ;
+        }
+
         if (order == 2)
         {
@@ -299 +331 @@
           sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id;
           jj++;
-          for (int i = 0; i < NMAX; i++)
+          for (int i = 0; i <  sstree.localElements[recvElement[rank][j]].n ; i++)
           {
             sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i];
@@ -316 +348 @@
       MPI_Issend(sendGivenArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
       nbSendRequest++;
+      MPI_Issend(sendVertexOffset[rank],  nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+
       if (order == 2)
       {
-        MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+        MPI_Issend(sendGrad[rank], 3*nbRecvVertex[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
         nbSendRequest++;
-        MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
+        MPI_Issend(sendNeighIds[rank], 4*nbRecvVertex[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
 //ym  --> attention taille GloId
         nbSendRequest++;
@@ -333 +368 @@
     if (nbSendElement[rank] > 0)
     {
+      recvValue[rank]   = new double[nbSendElement[rank]];
+      recvArea[rank]    = new double[nbSendElement[rank]];
+      recvGivenArea[rank] = new double[nbSendElement[rank]];
+      recvVertexOffset[rank] = new int[nbSendElement[rank]];
+
+      if (order == 2)
+      {
+        recvNeighIds[rank] = new GloId[nbSendVertex[rank]];
+        recvGrad[rank]     = new Coord[nbSendVertex[rank]];
+      }
+      else recvNeighIds[rank] = new GloId[nbSendElement[rank]];
+
+
       MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
       nbRecvRequest++;
@@ -339 +387 @@
       MPI_Irecv(recvGivenArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
       nbRecvRequest++;
+      MPI_Irecv(recvVertexOffset[rank],  nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+
       if (order == 2)
       {
-        MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1),
-            MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+        MPI_Irecv(recvGrad[rank], 3*nbSendVertex[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
         nbRecvRequest++;
-        MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+        MPI_Irecv(recvNeighIds[rank], 4*nbSendVertex[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
 //ym  --> attention taille GloId
         nbRecvRequest++;
@@ -357 +407 @@
   }
         
-        MPI_Waitall(nbSendRequest, sendRequest, status);
+  MPI_Waitall(nbSendRequest, sendRequest, status);
   MPI_Waitall(nbRecvRequest, recvRequest, status);
   
@@ -388 +438 @@
 
       /* first order: src value times weight (weight = supermesh area), later divide by target area */
-      int kk = (order == 2) ? n1 * (NMAX + 1) : n1;
-      GloId neighID = recvNeighIds[rank][kk];
+//      int kk = (order == 2) ? n1 * (NMAX + 1) : n1;
+      int offset ;
+      int nvertex ;
+      if (n1==0) 
+      {
+        offset=0 ;
+        nvertex=recvVertexOffset[rank][0] ;
+      }
+      else 
+      {
+        offset = recvVertexOffset[rank][n1-1];
+        nvertex = recvVertexOffset[rank][n1]-recvVertexOffset[rank][n1-1];
+      }
+
+      GloId neighID = recvNeighIds[rank][offset];
       wgt_map[neighID] += w;
 
       if (order == 2)
       {
-        for (int k = 0; k < NMAX+1; k++)
+        for (int k = 0; k < nvertex; k++)
         {
-          int kk = n1 * (NMAX + 1) + k;
+          int kk = offset + k;
           GloId neighID = recvNeighIds[rank][kk];
           if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x);
@@ -414 +477 @@
     for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++)
     {
-      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm;
-      else this->remapMatrix[i] = (it->second / e.area) / renorm;
-      this->srcAddress[i] = it->first.ind;
-      this->srcRank[i] = it->first.rank;
-      this->dstAddress[i] = j;
-      this->sourceWeightId[i]= it->first.globalId ;
-      this->targetWeightId[i]= targetGlobalId[j] ;
+      if (quantity)  this->remapMatrix.push_back((it->second ) / renorm);
+      else this->remapMatrix.push_back((it->second / e.area) / renorm);
+      this->srcAddress.push_back(it->first.ind);
+      this->srcRank.push_back(it->first.rank);
+      this->dstAddress.push_back(j);
+      this->sourceWeightId.push_back(it->first.globalId) ;
+      this->targetWeightId.push_back(targetGlobalId[j]) ;
       i++;
     }
@@ -434 +497 @@
       delete[] recvArea[rank];
       delete[] recvGivenArea[rank];
+      delete[] recvVertexOffset[rank];
       if (order == 2)
       {
@@ -446 +510 @@
       delete[] sendArea[rank];
       delete[] sendGivenArea[rank];
+      delete[] sendVertexOffset[rank];
       if (order == 2)
         delete[] sendGrad[rank];
@@ -457 +522 @@
   delete[] nbSendElement;
   delete[] nbRecvElement;
+  delete[] nbSendVertex;
+  delete[] nbRecvVertex;
   delete[] sendElement;
   delete[] recvElement;
@@ -939 +1006 @@
 Mapper::~Mapper()
 {
-  delete [] remapMatrix;
-  delete [] srcAddress;
-  delete [] srcRank;
-  delete [] dstAddress;
-  if (neighbourElements) delete [] neighbourElements;
-}
-
-}
+
+}
+
+}