Changeset 1158 for XIOS/dev/dev_olga/extern/remap/src/mapper.cpp

Timestamp:

06/06/17 17:58:16 (7 years ago)

Author:

oabramkina

Message:

Two server levels: merging with trunk r1137.
There are bugs.

File:

• XIOS/dev/dev_olga/extern/remap/src/mapper.cpp (modified) (16 diffs)

XIOS/dev/dev_olga/extern/remap/src/mapper.cpp

@@ -109 +109 @@
 }
 
-vector<double> Mapper::computeWeights(int interpOrder, bool renormalize)
+vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity)
 {
         vector<double> timings;
@@ -152 +152 @@
         if (mpiRank == 0 && verbose) cout << "Remapping..." << endl;
         tic = cputime();
-        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize);
+        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity);
         timings.push_back(cputime() - tic);
 
@@ -166 +166 @@
    @param order is the order of interpolaton (must be 1 or 2).
 */
-int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize)
+int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity)
 {
         int mpiSize, mpiRank;
@@ -184 +184 @@
         int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */
         double **recvValue = new double*[mpiSize];
+        double **recvArea = new double*[mpiSize];
         Coord **recvGrad = new Coord*[mpiSize];
         GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */
@@ -205 +206 @@
                         sendElement[rank] = new int[nbSendElement[rank]];
                         recvValue[rank]   = new double[nbSendElement[rank]];
+                        recvArea[rank]    = new double[nbSendElement[rank]];
                         if (order == 2)
                         {
@@ -234 +236 @@
         int **recvElement = new int*[mpiSize];
         double **sendValue = new double*[mpiSize];
+        double **sendArea = new double*[mpiSize];
         Coord **sendGrad = new Coord*[mpiSize];
         GloId **sendNeighIds = new GloId*[mpiSize];
-        MPI_Request *sendRequest = new MPI_Request[3*mpiSize];
-        MPI_Request *recvRequest = new MPI_Request[3*mpiSize];
+        MPI_Request *sendRequest = new MPI_Request[4*mpiSize];
+        MPI_Request *recvRequest = new MPI_Request[4*mpiSize];
         for (int rank = 0; rank < mpiSize; rank++)
         {
@@ -250 +253 @@
                         recvElement[rank] = new int[nbRecvElement[rank]];
                         sendValue[rank]   = new double[nbRecvElement[rank]];
+                        sendArea[rank]   = new double[nbRecvElement[rank]];
                         if (order == 2)
                         {
@@ -263 +267 @@
                 }
         }
-        MPI_Status *status = new MPI_Status[3*mpiSize];
-        MPI_Waitall(nbRecvRequest, recvRequest, status);
+        MPI_Status *status = new MPI_Status[4*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 */
@@ -278 +283 @@
                         {
                                 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val;
+                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area;
                                 if (order == 2)
                                 {
@@ -297 +303 @@
                         MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
                         nbSendRequest++;
+                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+                        nbSendRequest++;
                         if (order == 2)
                         {
@@ -317 +325 @@
                         MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
                         nbRecvRequest++;
+                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+                        nbRecvRequest++;
                         if (order == 2)
                         {
@@ -334 +344 @@
                 }
         }
+        
+        MPI_Waitall(nbSendRequest, sendRequest, status);
         MPI_Waitall(nbRecvRequest, recvRequest, status);
-        MPI_Waitall(nbSendRequest, sendRequest, status);
+        
 
         /* now that all values and gradients are available use them to computed interpolated values on target
@@ -357 +369 @@
                         int rank = (*it)->id.rank;
                         double fk = recvValue[rank][n1];
+                        double srcArea = recvArea[rank][n1];
                         double w = (*it)->area;
+      if (quantity) w/=srcArea ;
 
                         /* 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];
-                        wgt_map[neighID] += (*it)->area;
+                        wgt_map[neighID] += w;
 
                         if (order == 2)
@@ -383 +397 @@
     for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++)
                 {
-                        this->remapMatrix[i] = (it->second / e.area) / renorm;
+      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;
@@ -400 +415 @@
                         delete[] sendElement[rank];
                         delete[] recvValue[rank];
+                        delete[] recvArea[rank];
                         if (order == 2)
                         {
@@ -410 +426 @@
                         delete[] recvElement[rank];
                         delete[] sendValue[rank];
+                        delete[] sendArea[rank];
                         if (order == 2)
                                 delete[] sendGrad[rank];