Changeset 1114

Timestamp:

05/03/17 15:07:53 (8 years ago)

Author:

ymipsl

Message:

Add conservative remapping option for quantity, ie, the integrated value over a cell area.

YM

Location:

XIOS/trunk

Files:

• extern/remap/src/mapper.cpp (modified) (14 diffs)
• extern/remap/src/mapper.hpp (modified) (2 diffs)
• src/config/interpolate_domain_attribute.conf (modified) (1 diff)
• src/transformation/domain_algorithm_interpolate.cpp (modified) (2 diffs)

XIOS/trunk/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];
@@ -250 +253 @@
                         recvElement[rank] = new int[nbRecvElement[rank]];
                         sendValue[rank]   = new double[nbRecvElement[rank]];
+                        sendArea[rank]   = new double[nbRecvElement[rank]];
                         if (order == 2)
                         {
@@ -278 +282 @@
                         {
                                 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val;
+                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area;
                                 if (order == 2)
                                 {
@@ -297 +302 @@
                         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 +324 @@
                         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)
                         {
@@ -357 +366 @@
                         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 +394 @@
     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 +412 @@
                         delete[] sendElement[rank];
                         delete[] recvValue[rank];
+                        delete[] recvArea[rank];
                         if (order == 2)
                         {
@@ -410 +423 @@
                         delete[] recvElement[rank];
                         delete[] sendValue[rank];
+                        delete[] sendArea[rank];
                         if (order == 2)
                                 delete[] sendGrad[rank];

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

@@ -34 +34 @@
        /** @param trgElts are the elements of the unstructured target grid
            Returns the timings for substeps: */
-       vector<double> computeWeights(int interpOrder, bool renormalize=false);
+       vector<double> computeWeights(int interpOrder, bool renormalize=false, bool quantity=false);
        int getNbWeights(void) { return nWeights ; }
 /*
@@ -51 +51 @@
 private:
        /** @return number of weights (local to cpu) */
-       int remap(Elt* elements, int nbElements, int order, bool renormalize=false);
+       int remap(Elt* elements, int nbElements, int order, bool renormalize=false, bool quantity=false);
 
        void buildMeshTopology();

XIOS/trunk/src/config/interpolate_domain_attribute.conf

@@ -2 +2 @@
 DECLARE_ATTRIBUTE(int, order)
 DECLARE_ATTRIBUTE(bool, renormalize)
+DECLARE_ATTRIBUTE(bool, quantity)
 
 /* Write interpolation weights into file */

XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

@@ -100 +100 @@
   int orderInterp = interpDomain_->order.getValue();
   bool renormalize ;
+  bool quantity ;
 
   if (interpDomain_->renormalize.isEmpty()) renormalize=true;
   else renormalize = interpDomain_->renormalize;
+
+  if (interpDomain_->quantity.isEmpty()) quantity=false;
+  else quantity = interpDomain_->quantity;
 
   const double poleValue = 90.0;
@@ -342 +346 @@
   mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked);
 
-  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize);
+  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity);
 
   std::map<int,std::vector<std::pair<int,double> > > interpMapValue;