Changeset 2538

Timestamp:

07/24/23 11:55:12 (10 months ago)

Author:

jderouillat

Message:

Backport 2506 : limitation of number of vertex max in remaper

Location:

XIOS3/trunk/extern/remap/src

Files:

• elt.hpp (modified) (6 diffs)
• intersection_ym.cpp (modified) (1 diff)
• libmapper.cpp (modified) (1 diff)
• mapper.cpp (modified) (19 diffs)
• mapper.hpp (modified) (1 diff)
• meshutil.cpp (modified) (7 diffs)
• polyg.cpp (modified) (1 diff)

XIOS3/trunk/extern/remap/src/elt.hpp

@@ -3 +3 @@
 #include <list>
 #include "triple.hpp"
+#include <vector>
 
-#define NMAX 10 /**< maximum number of vertices for polygons */
+#define NMAX 0 /**< maximum number of vertices for polygons */
 
 #define NOT_FOUND -1
@@ -58 +59 @@
         {
                 int k = 0;
+                vertex.resize(max_num_vert) ;
                 vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]);
                 for (int i = 1; i < max_num_vert; i++)
@@ -72 +74 @@
                 }
                 n = k;
-                x = barycentre(vertex, n);
+    vertex.resize(n) ;
+    vertex.shrink_to_fit();
+          allocate() ;
+
+                x = barycentre(vertex.data(), n);
         }
-
+  void allocate(void)
+  {
+    vertex.resize(n) ;
+    neighbour.resize(n) ;
+    d.resize(n) ;
+    edge.resize(n) ;
+    gradNeigh.resize(n) ;
+    neighId.resize(n) ;
+  }
         Elt& operator=(const Elt& rhs)
         {
@@ -87 +101 @@
                 is   = rhs.is;
 
-                for(int i = 0; i < NMAX; i++)
-                {
-                        neighbour[i] = rhs.neighbour[i];
-                        d[i]         = rhs.d[i];
-                        edge[i]      = rhs.edge[i];
-                        vertex[i]    = rhs.vertex[i];
-                        gradNeigh[i] = rhs.gradNeigh[i];
-                }
+                neighbour = rhs.neighbour;
+                d         = rhs.d;
+                edge      = rhs.edge;
+                vertex    = rhs.vertex;
+                gradNeigh = rhs.gradNeigh;
                 return *this;
         }
@@ -109 +120 @@
   void insert_vertex(int i, const Coord& v)
   {
+    vertex.resize(n+1) ;
+    edge.resize(n+1) ;
+    d.resize(n+1) ;
+    neighbour.resize(n+1) ;
+    gradNeigh.resize(n+1) ;
+    neighId.resize(n+1) ;
+
     for(int j=n; j > i ; j--)
     {
@@ -120 +138 @@
   }
   
-        int neighbour[NMAX];
-        double d[NMAX]; /**< distance of centre of small circle to origin, zero if great circle */
+        std::vector<int> neighbour;
+        std::vector<double> d; /**< distance of centre of small circle to origin, zero if great circle */
         double val;     /**< value (sample if src element, interpolated if dest element) */
-        Coord vertex[NMAX];
-        Coord edge[NMAX]; /**< edge normals: if great circle tangential to sphere, if small circle parallel to pole */
+        std::vector<Coord> vertex;
+        std::vector<Coord> edge; /**< edge normals: if great circle tangential to sphere, if small circle parallel to pole */
         Coord grad;    /**< gradient of the reconstructed linear function over this element */
-        Coord gradNeigh[NMAX]; /**< for weight computation: gradients for val=1 at individual neighbours */
-        struct GloId neighId[NMAX]; /**< weight computation needs to know global IDs for all sources with "link" */
+        std::vector<Coord> gradNeigh; /**< for weight computation: gradients for val=1 at individual neighbours */
+        std::vector<struct GloId> neighId; /**< weight computation needs to know global IDs for all sources with "link" */
         std::list<Polyg*> is; /**< intersections */
 };

XIOS3/trunk/extern/remap/src/intersection_ym.cpp

@@ -47 +47 @@
   Coord *b_gno   = new Coord[nb];
 
-  Coord OC=barycentre(a->vertex,a->n) ;
+  Coord OC=barycentre(a->vertex.data(),a->n) ;
   Coord Oz=OC ;
   Coord Ox=crossprod(Coord(0,0,1),Oz) ;

XIOS3/trunk/extern/remap/src/libmapper.cpp

@@ -138 +138 @@
 extern "C" void remap_get_weights(double* weights, int* src_indices, int* dst_indices)
 {
-        memcpy(weights, mapper->remapMatrix, mapper->nWeights*sizeof(double));
-        memcpy(src_indices, mapper->srcAddress, mapper->nWeights*sizeof(int));
-        memcpy(dst_indices, mapper->dstAddress, mapper->nWeights*sizeof(int));
+        memcpy(weights, mapper->remapMatrix.data(), mapper->nWeights*sizeof(double));
+        memcpy(src_indices, mapper->srcAddress.data(), mapper->nWeights*sizeof(int));
+        memcpy(dst_indices, mapper->dstAddress.data(), mapper->nWeights*sizeof(int));
         delete mapper;
 }

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

@@ -154 +154 @@
       nIntersections++;
   }
+
   /* overallocate for NMAX neighbours for each elements */
+/*
   remapMatrix = new double[nIntersections*NMAX];
   srcAddress = new int[nIntersections*NMAX];
@@ -161 +163 @@
   sourceWeightId =new long[nIntersections*NMAX];
   targetWeightId =new long[nIntersections*NMAX];
-
+*/
 
   if (mpiRank == 0 && verbose) cout << "Remapping..." << endl;
@@ -195 +197 @@
 
   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 */
@@ -219 +223 @@
     {
       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;
@@ -244 +237 @@
   /* 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);
 
@@ -253 +247 @@
   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++)
   {
@@ -268 +263 @@
     {
       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++)
       {
@@ -302 +322 @@
         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)
         {
@@ -308 +340 @@
           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];
@@ -325 +357 @@
       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++;
@@ -342 +377 @@
     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++;
@@ -348 +396 @@
       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++;
@@ -366 +416 @@
   }
         
-        MPI_Waitall(nbSendRequest, sendRequest, status);
+  MPI_Waitall(nbSendRequest, sendRequest, status);
   MPI_Waitall(nbRecvRequest, recvRequest, status);
   
@@ -397 +447 @@
 
       /* 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);
@@ -423 +486 @@
     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++;
     }
@@ -443 +506 @@
       delete[] recvArea[rank];
       delete[] recvGivenArea[rank];
+      delete[] recvVertexOffset[rank];
       if (order == 2)
       {
@@ -455 +519 @@
       delete[] sendArea[rank];
       delete[] sendGivenArea[rank];
+      delete[] sendVertexOffset[rank];
       if (order == 2)
         delete[] sendGrad[rank];
@@ -466 +531 @@
   delete[] nbSendElement;
   delete[] nbRecvElement;
+  delete[] nbSendVertex;
+  delete[] nbRecvVertex;
   delete[] sendElement;
   delete[] recvElement;
@@ -956 +1023 @@
 Mapper::~Mapper()
 {
-  delete [] remapMatrix;
-  delete [] srcAddress;
-  delete [] srcRank;
-  delete [] dstAddress;
-  delete [] sourceWeightId ;
-  delete [] targetWeightId ;
-  if (neighbourElements) delete [] neighbourElements;
+  //delete [] remapMatrix;
+  //delete [] srcAddress;
+  //delete [] srcRank;
+  //delete [] dstAddress;
+  //delete [] sourceWeightId ;
+  //delete [] targetWeightId ;
+  //if (neighbourElements) delete [] neighbourElements;
   CTimer::release() ;
  }

XIOS3/trunk/extern/remap/src/mapper.hpp

@@ -42 +42 @@
 */
        /* where weights are returned after call to `computeWeights` */
-       double *remapMatrix;
-       int *srcAddress;
-       int *srcRank;
-       int *dstAddress;
+       std::vector<double> remapMatrix;
+       std::vector<int> srcAddress;
+       std::vector<int> srcRank;
+       std::vector<int> dstAddress;
        int nWeights;
-       long int* sourceWeightId ;
-       long int* targetWeightId ;
+       std::vector<long int> sourceWeightId ;
+       std::vector<long int> targetWeightId ;
 
 private:

XIOS3/trunk/extern/remap/src/meshutil.cpp

@@ -24 +24 @@
   Coord *a_gno   = new Coord[na];
   
-  Coord OC=barycentre(a.vertex,a.n) ;
+  Coord OC=barycentre(a.vertex.data(),a.n) ;
   Coord Oz=OC ;
   Coord Ox=crossprod(Coord(0,0,1),Oz) ;
@@ -67 +67 @@
   {
     bary = bary * (-1.) ;
-    switchOrientation(a.n, a.vertex,a.edge,a.d) ;
+    switchOrientation(a.n, a.vertex.data(),a.edge.data(),a.d.data()) ;
   }
   
@@ -79 +79 @@
 void cptEltGeom(Elt& elt, const Coord &pole)
 {
-   orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x);
+   orient(elt.n, elt.vertex.data(), elt.edge.data(), elt.d.data(), elt.x);
    normals(elt, pole);
 //  Coord gg;
@@ -205 +205 @@
     elts[j]->val = 0;
 
-  Elt *neighbours[NMAX];
+//  Elt *neighbours[NMAX];
   for (int j = 0; j < N; j++)
   {
+    vector<Elt*> neighbours(elts[j]->n) ;
+
     for (int i = 0; i < elts[j]->n; i++) 
       if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; // no neighbour
@@ -223 +225 @@
         /* for weight computation all values are always kept zero and only set to one when used .. */
         neighbours[i]->val = 1;
-        elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours);
+        elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours.data());
         /* .. and right after zeroed again */
         neighbours[i]->val = 0;
@@ -234 +236 @@
     }
 
+/*  not needed anymore
     for(int i = elts[j]->n ; i < NMAX; i++)
     {
@@ -239 +242 @@
       elts[j]->neighId[i].ind = -1; // mark end
     }
+*/
     /* For the most naive algorithm the case where the element itself is one must also be considered.
        Thomas says this can later be optimized out. */
     elts[j]->val = 1;
-    elts[j]->grad = gradient(*(elts[j]), neighbours);
+    elts[j]->grad = gradient(*(elts[j]), neighbours.data());
     elts[j]->val = 0;
   }

XIOS3/trunk/extern/remap/src/polyg.cpp

@@ -304 +304 @@
   e.n = *((int *) & (buffer[pos]));
   pos += sizeof(int);
-
+  
+  e.allocate() ;
+  
   for (int i = 0; i < e.n; i++)
   {