Ignore:
Timestamp:
11/21/17 10:47:57 (6 years ago)
Author:
yushan
Message:

dev_omp

Location:
XIOS/dev/branch_openmp/extern
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp

    r1335 r1338  
    132132  timings.push_back(cputime() - tic); 
    133133 
    134         tic = cputime(); 
    135         if (interpOrder == 2) { 
    136                 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
    137                 buildMeshTopology(); 
    138                 computeGrads(); 
    139         } 
    140         timings.push_back(cputime() - tic); 
    141  
    142         /* Prepare computation of weights */ 
    143         /* compute number of intersections which for the first order case 
    144            corresponds to the number of edges in the remap matrix */ 
    145         int nIntersections = 0; 
    146         for (int j = 0; j < targetElements.size(); j++) 
    147         { 
    148                 Elt &elt = targetElements[j]; 
    149                 for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
    150                         nIntersections++; 
    151         } 
    152         /* overallocate for NMAX neighbours for each elements */ 
    153         remapMatrix = new double[nIntersections*NMAX]; 
    154         srcAddress = new int[nIntersections*NMAX]; 
    155         srcRank = new int[nIntersections*NMAX]; 
    156         dstAddress = new int[nIntersections*NMAX]; 
     134  tic = cputime(); 
     135  if (interpOrder == 2)  
     136  { 
     137    if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
     138    buildMeshTopology(); 
     139    computeGrads(); 
     140  } 
     141  timings.push_back(cputime() - tic); 
     142 
     143  /* Prepare computation of weights */ 
     144  /* compute number of intersections which for the first order case 
     145     corresponds to the number of edges in the remap matrix */ 
     146  int nIntersections = 0; 
     147  for (int j = 0; j < targetElements.size(); j++) 
     148  { 
     149    Elt &elt = targetElements[j]; 
     150    for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
     151      nIntersections++; 
     152  } 
     153  /* overallocate for NMAX neighbours for each elements */ 
     154  remapMatrix = new double[nIntersections*NMAX]; 
     155  srcAddress = new int[nIntersections*NMAX]; 
     156  srcRank = new int[nIntersections*NMAX]; 
     157  dstAddress = new int[nIntersections*NMAX]; 
    157158  sourceWeightId =new long[nIntersections*NMAX]; 
    158159  targetWeightId =new long[nIntersections*NMAX]; 
    159160 
    160161 
    161         if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
    162         tic = cputime(); 
    163         nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
    164         timings.push_back(cputime() - tic); 
     162  if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
     163  tic = cputime(); 
     164  nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
     165  timings.push_back(cputime() - tic); 
    165166 
    166167  for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 
    167168 
    168         return timings; 
     169  return timings; 
    169170} 
    170171 
    … …  
    177178int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 
    178179{ 
    179         int mpiSize, mpiRank; 
    180         MPI_Comm_size(communicator, &mpiSize); 
    181         MPI_Comm_rank(communicator, &mpiRank); 
    182  
    183         /* create list of intersections (super mesh elements) for each rank */ 
    184         multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
    185         for (int j = 0; j < nbElements; j++) 
    186         { 
    187                 Elt& e = elements[j]; 
    188                 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
    189                         elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
    190         } 
    191  
    192         int *nbSendElement = new int[mpiSize]; 
    193         int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
    194         double **recvValue = new double*[mpiSize]; 
    195         double **recvArea = new double*[mpiSize]; 
    196         Coord **recvGrad = new Coord*[mpiSize]; 
    197         GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
    198         for (int rank = 0; rank < mpiSize; rank++) 
    199         { 
    200                 /* get size for allocation */ 
    201                 int last = -1; /* compares unequal to any index */ 
    202                 int index = -1; /* increased to starting index 0 in first iteration */ 
    203                 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
    204                 { 
    205                         if (last != it->first) 
    206                                 index++; 
    207                         (it->second)->id.ind = index; 
    208                         last = it->first; 
    209                 } 
    210                 nbSendElement[rank] = index + 1; 
    211  
    212                 /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
    213                 if (nbSendElement[rank] > 0) 
    214                 { 
    215                         sendElement[rank] = new int[nbSendElement[rank]]; 
    216                         recvValue[rank]   = new double[nbSendElement[rank]]; 
    217                         recvArea[rank]    = new double[nbSendElement[rank]]; 
    218                         if (order == 2) 
    219                         { 
    220                                 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
    221                                 recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
    222                         } 
    223                         else 
    224                                 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
    225  
    226                         last = -1; 
    227                         index = -1; 
    228                         for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
    229                         { 
    230                                 if (last != it->first) 
    231                                         index++; 
    232                                 sendElement[rank][index] = it->first; 
    233                                 last = it->first; 
    234                         } 
    235                 } 
    236         } 
    237  
    238         /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
    239         int *nbRecvElement = new int[mpiSize]; 
    240         MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
    241  
    242         /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
    243         int nbSendRequest = 0; 
    244         int nbRecvRequest = 0; 
    245         int **recvElement = new int*[mpiSize]; 
    246         double **sendValue = new double*[mpiSize]; 
    247         double **sendArea = new double*[mpiSize]; 
    248         Coord **sendGrad = new Coord*[mpiSize]; 
    249         GloId **sendNeighIds = new GloId*[mpiSize]; 
    250         MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
    251         MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
    252         for (int rank = 0; rank < mpiSize; rank++) 
    253         { 
    254                 if (nbSendElement[rank] > 0) 
    255                 { 
    256                         MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    257                         nbSendRequest++; 
    258                 } 
    259  
    260                 if (nbRecvElement[rank] > 0) 
    261                 { 
    262                         recvElement[rank] = new int[nbRecvElement[rank]]; 
    263                         sendValue[rank]   = new double[nbRecvElement[rank]]; 
    264                         sendArea[rank]   = new double[nbRecvElement[rank]]; 
    265                         if (order == 2) 
    266                         { 
    267                                 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
    268                                 sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
    269                         } 
    270                         else 
    271                         { 
    272                                 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
    273                         } 
    274                         MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    275                         nbRecvRequest++; 
    276                 } 
    277         } 
     180  int mpiSize, mpiRank; 
     181  MPI_Comm_size(communicator, &mpiSize); 
     182  MPI_Comm_rank(communicator, &mpiRank); 
     183 
     184  /* create list of intersections (super mesh elements) for each rank */ 
     185  multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
     186  for (int j = 0; j < nbElements; j++) 
     187  { 
     188    Elt& e = elements[j]; 
     189    for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     190    { 
     191      elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
     192      //std::cout<<"elementList["<<(*it)->id.rank<<"].size = "<< elementList[(*it)->id.rank].size()<<std::endl; 
     193    } 
     194  } 
     195 
     196  int *nbSendElement = new int[mpiSize]; 
     197  int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
     198  double **recvValue = new double*[mpiSize]; 
     199  double **recvArea = new double*[mpiSize]; 
     200  Coord **recvGrad = new Coord*[mpiSize]; 
     201  GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
     202  for (int rank = 0; rank < mpiSize; rank++) 
     203  { 
     204    /* get size for allocation */ 
     205    int last = -1; /* compares unequal to any index */ 
     206    int index = -1; /* increased to starting index 0 in first iteration */ 
     207    for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     208    { 
     209      if (last != it->first) 
     210        index++; 
     211      (it->second)->id.ind = index; 
     212      last = it->first; 
     213    } 
     214    nbSendElement[rank] = index + 1; 
     215 
     216    /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
     217    if (nbSendElement[rank] > 0) 
     218    { 
     219      sendElement[rank] = new int[nbSendElement[rank]]; 
     220      recvValue[rank]   = new double[nbSendElement[rank]]; 
     221      recvArea[rank]    = new double[nbSendElement[rank]]; 
     222      if (order == 2) 
     223      { 
     224        recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
     225        recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
     226      } 
     227      else 
     228        recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
     229 
     230      last = -1; 
     231      index = -1; 
     232      for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     233      { 
     234        if (last != it->first) 
     235          index++; 
     236        sendElement[rank][index] = it->first; 
     237        last = it->first; 
     238      } 
     239    } 
     240  } 
     241 
     242  /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
     243  int *nbRecvElement = new int[mpiSize]; 
     244  MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
     245 
     246  /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
     247  int nbSendRequest = 0; 
     248  int nbRecvRequest = 0; 
     249  int **recvElement = new int*[mpiSize]; 
     250  double **sendValue = new double*[mpiSize]; 
     251  double **sendArea = new double*[mpiSize]; 
     252  Coord **sendGrad = new Coord*[mpiSize]; 
     253  GloId **sendNeighIds = new GloId*[mpiSize]; 
     254  MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
     255  MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
     256  for (int rank = 0; rank < mpiSize; rank++) 
     257  { 
     258    if (nbSendElement[rank] > 0) 
     259    { 
     260      MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     261      nbSendRequest++; 
     262    } 
     263 
     264    if (nbRecvElement[rank] > 0) 
     265    { 
     266      recvElement[rank] = new int[nbRecvElement[rank]]; 
     267      sendValue[rank]   = new double[nbRecvElement[rank]]; 
     268      sendArea[rank]   = new double[nbRecvElement[rank]]; 
     269      if (order == 2) 
     270      { 
     271        sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
     272        sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
     273      } 
     274      else 
     275      { 
     276        sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
     277      } 
     278      MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     279      nbRecvRequest++; 
     280    } 
     281  } 
    278282 
    279283  MPI_Status *status = new MPI_Status[4*mpiSize]; 
    … …  
    282286  MPI_Waitall(nbRecvRequest, recvRequest, status); 
    283287 
    284         /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
    285         nbSendRequest = 0; 
    286         nbRecvRequest = 0; 
    287         for (int rank = 0; rank < mpiSize; rank++) 
    288         { 
    289                 if (nbRecvElement[rank] > 0) 
    290                 { 
    291                         int jj = 0; // jj == j if no weight writing 
    292                         for (int j = 0; j < nbRecvElement[rank]; j++) 
    293                         { 
    294                                 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
    295                                 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
    296                                 if (order == 2) 
    297                                 { 
    298                                         sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
    299                                         sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
    300                                         jj++; 
    301                                         for (int i = 0; i < NMAX; i++) 
    302                                         { 
    303                                                 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
     288  /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
     289  nbSendRequest = 0; 
     290  nbRecvRequest = 0; 
     291  for (int rank = 0; rank < mpiSize; rank++) 
     292  { 
     293    if (nbRecvElement[rank] > 0) 
     294    { 
     295      int jj = 0; // jj == j if no weight writing 
     296      for (int j = 0; j < nbRecvElement[rank]; j++) 
     297      { 
     298        sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
     299        sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
     300        if (order == 2) 
     301        { 
     302          sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
     303          sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
     304          jj++; 
     305          for (int i = 0; i < NMAX; i++) 
     306          { 
     307            sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
    304308            sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 
    305                                                 jj++; 
    306                                         } 
    307                                 } 
    308                                 else 
    309                                         sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
    310                         } 
    311                         MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    312                         nbSendRequest++; 
    313                         MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    314                         nbSendRequest++; 
    315                         if (order == 2) 
    316                         { 
    317                                 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 
    318                                                                 MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    319                                 nbSendRequest++; 
    320                                 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     309            jj++; 
     310          } 
     311        } 
     312        else 
     313          sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
     314      } 
     315      MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     316      nbSendRequest++; 
     317      MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     318      nbSendRequest++; 
     319      if (order == 2) 
     320      { 
     321        MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     322        nbSendRequest++; 
     323        MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    321324//ym  --> attention taille GloId 
    322                                 nbSendRequest++; 
    323                         } 
    324                         else 
    325                         { 
    326                                 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     325        nbSendRequest++; 
     326      } 
     327      else 
     328      { 
     329        MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    327330//ym  --> attention taille GloId 
    328                                 nbSendRequest++; 
    329                         } 
    330                 } 
    331                 if (nbSendElement[rank] > 0) 
    332                 { 
    333                         MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    334                         nbRecvRequest++; 
    335                         MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    336                         nbRecvRequest++; 
    337                         if (order == 2) 
    338                         { 
    339                                 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 
    340                                                 MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    341                                 nbRecvRequest++; 
    342                                 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     331        nbSendRequest++; 
     332      } 
     333    } 
     334    if (nbSendElement[rank] > 0) 
     335    { 
     336      MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     337      nbRecvRequest++; 
     338      MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     339      nbRecvRequest++; 
     340      if (order == 2) 
     341      { 
     342        MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     343        nbRecvRequest++; 
     344        MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    343345//ym  --> attention taille GloId 
    344                                 nbRecvRequest++; 
    345                         } 
    346                         else 
    347                         { 
    348                                 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     346        nbRecvRequest++; 
     347      } 
     348      else 
     349      { 
     350        MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    349351//ym  --> attention taille GloId 
    350                                 nbRecvRequest++; 
    351                         } 
    352                 } 
    353         } 
     352        nbRecvRequest++; 
     353      } 
     354    } 
     355  } 
    354356         
    355         MPI_Waitall(nbSendRequest, sendRequest, status); 
    356         MPI_Waitall(nbRecvRequest, recvRequest, status); 
     357  MPI_Waitall(nbSendRequest, sendRequest, status); 
     358  MPI_Waitall(nbRecvRequest, recvRequest, status); 
    357359         
    358360 
    359         /* now that all values and gradients are available use them to computed interpolated values on target 
    360            and also to compute weights */ 
    361         int i = 0; 
    362         for (int j = 0; j < nbElements; j++) 
    363         { 
    364                 Elt& e = elements[j]; 
    365  
    366                 /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
    367                    (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
    368                    accumulate them so that there is only one final weight between two elements */ 
    369                 map<GloId,double> wgt_map; 
    370  
    371                 /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
    372                 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
    373                 { 
    374                         /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
    375                         but it->id is id of the source element that it intersects */ 
    376                         int n1 = (*it)->id.ind; 
    377                         int rank = (*it)->id.rank; 
    378                         double fk = recvValue[rank][n1]; 
    379                         double srcArea = recvArea[rank][n1]; 
    380                         double w = (*it)->area; 
     361  /* now that all values and gradients are available use them to computed interpolated values on target 
     362    and also to compute weights */ 
     363  int i = 0; 
     364  for (int j = 0; j < nbElements; j++) 
     365  { 
     366    Elt& e = elements[j]; 
     367 
     368    /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
     369    (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
     370    accumulate them so that there is only one final weight between two elements */ 
     371    map<GloId,double> wgt_map; 
     372 
     373    /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
     374    for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     375    { 
     376      /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
     377      but it->id is id of the source element that it intersects */ 
     378      int n1 = (*it)->id.ind; 
     379      int rank = (*it)->id.rank; 
     380      double fk = recvValue[rank][n1]; 
     381      double srcArea = recvArea[rank][n1]; 
     382      double w = (*it)->area; 
    381383      if (quantity) w/=srcArea ; 
    382384 
    383                         /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
    384                         int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
    385                         GloId neighID = recvNeighIds[rank][kk]; 
    386                         wgt_map[neighID] += w; 
    387  
    388                         if (order == 2) 
    389                         { 
    390                                 for (int k = 0; k < NMAX+1; k++) 
    391                                 { 
    392                                         int kk = n1 * (NMAX + 1) + k; 
    393                                         GloId neighID = recvNeighIds[rank][kk]; 
    394                                         if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
    395                                 } 
    396  
    397                         } 
    398                 } 
     385      /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
     386      int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
     387      GloId neighID = recvNeighIds[rank][kk]; 
     388      wgt_map[neighID] += w; 
     389 
     390      if (order == 2) 
     391      { 
     392        for (int k = 0; k < NMAX+1; k++) 
     393        { 
     394          int kk = n1 * (NMAX + 1) + k; 
     395          GloId neighID = recvNeighIds[rank][kk]; 
     396          if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
     397        } 
     398 
     399      } 
     400    } 
    399401 
    400402    double renorm=0; 
    … …  
    404406 
    405407    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
    406                 { 
     408    { 
    407409      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm; 
    408                         else this->remapMatrix[i] = (it->second / e.area) / renorm; 
    409                         this->srcAddress[i] = it->first.ind; 
    410                         this->srcRank[i] = it->first.rank; 
    411                         this->dstAddress[i] = j; 
     410      else this->remapMatrix[i] = (it->second / e.area) / renorm; 
     411      this->srcAddress[i] = it->first.ind; 
     412      this->srcRank[i] = it->first.rank; 
     413      this->dstAddress[i] = j; 
    412414      this->sourceWeightId[i]= it->first.globalId ; 
    413415      this->targetWeightId[i]= targetGlobalId[j] ; 
    414                         i++; 
    415                 } 
    416         } 
    417  
    418         /* free all memory allocated in this function */ 
    419         for (int rank = 0; rank < mpiSize; rank++) 
    420         { 
    421                 if (nbSendElement[rank] > 0) 
    422                 { 
    423                         delete[] sendElement[rank]; 
    424                         delete[] recvValue[rank]; 
    425                         delete[] recvArea[rank]; 
    426                         if (order == 2) 
    427                         { 
    428                                 delete[] recvGrad[rank]; 
    429                         } 
    430                         delete[] recvNeighIds[rank]; 
    431                 } 
    432                 if (nbRecvElement[rank] > 0) 
    433                 { 
    434                         delete[] recvElement[rank]; 
    435                         delete[] sendValue[rank]; 
    436                         delete[] sendArea[rank]; 
    437                         if (order == 2) 
    438                                 delete[] sendGrad[rank]; 
    439                         delete[] sendNeighIds[rank]; 
    440                 } 
    441         } 
    442         delete[] status; 
    443         delete[] sendRequest; 
    444         delete[] recvRequest; 
    445         delete[] elementList; 
    446         delete[] nbSendElement; 
    447         delete[] nbRecvElement; 
    448         delete[] sendElement; 
    449         delete[] recvElement; 
    450         delete[] sendValue; 
    451         delete[] recvValue; 
    452         delete[] sendGrad; 
    453         delete[] recvGrad; 
    454         delete[] sendNeighIds; 
    455         delete[] recvNeighIds; 
    456         return i; 
     416      i++; 
     417    } 
     418  } 
     419 
     420  /* free all memory allocated in this function */ 
     421 /* for (int rank = 0; rank < mpiSize; rank++) 
     422  { 
     423    if (nbSendElement[rank] > 0) 
     424    { 
     425      delete[] sendElement[rank]; 
     426      delete[] recvValue[rank]; 
     427      delete[] recvArea[rank]; 
     428      if (order == 2) 
     429      { 
     430        delete[] recvGrad[rank]; 
     431      } 
     432      delete[] recvNeighIds[rank]; 
     433    } 
     434    if (nbRecvElement[rank] > 0) 
     435    { 
     436      delete[] recvElement[rank]; 
     437      delete[] sendValue[rank]; 
     438      delete[] sendArea[rank]; 
     439      if (order == 2) 
     440        delete[] sendGrad[rank]; 
     441      delete[] sendNeighIds[rank]; 
     442    } 
     443  } 
     444  delete[] status; 
     445  delete[] sendRequest; 
     446  delete[] recvRequest; 
     447  delete[] elementList; 
     448  delete[] nbSendElement; 
     449  delete[] nbRecvElement; 
     450  delete[] sendElement; 
     451  delete[] recvElement; 
     452  delete[] sendValue; 
     453  delete[] recvValue; 
     454  delete[] sendGrad; 
     455  delete[] recvGrad; 
     456  delete[] sendNeighIds; 
     457  delete[] recvNeighIds;*/ 
     458  return i; 
    457459} 
    458460 

  • XIOS/dev/branch_openmp/extern/src_ep_dev/ep_split.cpp

    r1295 r1338  
    9999    MPI_Gather_local(&key, 1, MPI_INT, key_loc.data(), 0, comm); 
    100100    MPI_Bcast_local(key_loc.data(), num_ep, MPI_INT, 0, comm); 
     101     
     102    std::sort(key_loc.begin(), key_loc.end()); 
    101103 
    102104    for(int i=0; i<num_ep; i++) 
Note: See TracChangeset for help on using the changeset viewer.