Changeset 743

Timestamp:

10/19/15 17:41:35 (9 years ago)

Author:

mhnguyen

Message:

Processing pole of rectangular grid with latitude = 90

+) Use average value of all points with latitude 90 (-90)

Test
+) On Curie
+) test_remap pass

Location:

XIOS/trunk

Files:

• inputs/REMAP/iodef.xml (modified) (3 diffs)
• src/node/domain.cpp (modified) (2 diffs)
• src/node/domain.hpp (modified) (1 diff)
• src/transformation/domain_algorithm_interpolate.cpp (modified) (4 diffs)
• src/transformation/domain_algorithm_interpolate.hpp (modified) (1 diff)

XIOS/trunk/inputs/REMAP/iodef.xml

@@ -12 +12 @@
      <field id="dst_field_regular_grid"  operation="instant" field_ref="tmp_field" grid_ref="dst_grid_regular"  read_access="true"/>
      <field id="dst_field_regular"  operation="instant" field_ref="tmp_field" domain_ref="dst_domain_regular"  read_access="true"/>
+     <field id="dst_field_regular_pole"  operation="instant" field_ref="src_field" domain_ref="dst_domain_regular_pole"  read_access="true"/>
 
    </field_definition>
@@ -22 +23 @@
      <file id="output_dst" name="output_dst" type="one_file">
         <field field_ref="dst_field" name="field" />
+     </file>
+     <file id="out_dst_regular_pole" name="out_dst_regular_pole" type="one_file">
+        <field field_ref="dst_field_regular_pole" name="field" />
      </file>
      <file id="output_dst_regular" name="output_dst_regular" type="one_file">
@@ -42 +46 @@
      <domain id="src_domain" />
      <domain id="dst_domain" domain_src="src_domain">
+       <interpolate_domain/>
+     </domain>
+     <domain id="dst_domain_regular_pole" domain_src="src_domain" ni_glo="180" nj_glo="90" type="rectilinear">
+       <generate_rectilinear_domain bounds_lat_start="-90" bounds_lat_end="90" lon_start="2" lon_end="360"/>
        <interpolate_domain/>
      </domain>

XIOS/trunk/src/node/domain.cpp

@@ -375 +375 @@
    }
 
-   void CDomain::fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat)
+   void CDomain::fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat,
+                                              bool isNorthPole, bool isSouthPole)
    {
      int i,j,k;
@@ -397 +398 @@
          boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) != ni_glo) ? (ibegin + i +1) * lonStep + bounds_lon_start
                                                                         : bounds_lon_end;
-
+       }
+
+     double bounds_lat_start_pole = bounds_lat_start;
+     double bounds_lat_end_pole   = bounds_lat_end;
+     if (isNorthPole) bounds_lat_start_pole = lat_start;
+     if (isSouthPole) bounds_lat_end_pole   = lat_end;
+
+     for(j=0;j<nj;++j)
+       for(i=0;i<ni;++i)
+       {
+         k=j*ni+i;
          boundsLat(1,k) = boundsLat(2,k) = (0 != (jbegin + j)) ? (jbegin + j) * latStep + bounds_lat_start
-                                                               : bounds_lat_start;
+                                                               : bounds_lat_start_pole;
          boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) != nj_glo) ? (jbegin + j +1) * latStep + bounds_lat_start
-                                                               : bounds_lat_end;
+                                                               : bounds_lat_end_pole;
        }
    }

XIOS/trunk/src/node/domain.hpp

@@ -131 +131 @@
          void sendLonLatArea(void);
          void computeConnectedServer(void) ;
-         void fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat);
+         void fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat,
+                                           bool isNorthPole = false, bool isSouthPole = false);
 
          static bool dispatchEvent(CEventServer& event);

XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

@@ -39 +39 @@
         int orderInterp = interpDomain_->order.getValue();
 
-  int constNVertex = 4; // Value by default number of vertex for rectangular domain
+  const double poleValue = 90.0;
+  const int constNVertex = 4; // Value by default number of vertex for rectangular domain
   int nVertexSrc, nVertexDest;
   nVertexSrc = nVertexDest = constNVertex;
@@ -79 +80 @@
   }
 
+  bool isNorthPole = false;
+  bool isSouthPole = false;
+  std::map<int,std::vector<std::pair<int,double> > > interpMapValueNorthPole;
+  std::map<int,std::vector<std::pair<int,double> > > interpMapValueSouthPole;
+
   int localDomainDestSize = domainDest_->i_index.numElements();
   int niDest = domainDest_->ni.getValue(), njDest = domainDest_->nj.getValue();
@@ -110 +116 @@
   else
   {
+    if (poleValue == std::abs(domainDest_->lat_start)) isNorthPole = true;
+    if (poleValue == std::abs(domainDest_->lat_end)) isSouthPole = true;
+    if (isNorthPole && (0 == domainDest_->jbegin.getValue()))
+    {
+      int ibegin = domainDest_->ibegin.getValue();
+      for (i = 0; i < niDest; ++i)
+      {
+        interpMapValueNorthPole[i+ibegin];
+      }
+    }
+
+    if (isSouthPole && (domainDest_->nj_glo.getValue() == (domainDest_->jbegin.getValue() + njDest)))
+    {
+      int ibegin = domainDest_->ibegin.getValue();
+      int njGlo = domainDest_->nj_glo.getValue();
+      int niGlo = domainDest_->ni_glo.getValue();
+      for (i = 0; i < niDest; ++i)
+      {
+        k = (njGlo - 1)*niGlo + i + ibegin;
+        interpMapValueSouthPole[k];
+      }
+    }
+
     // Ok, fill in boundary values for rectangular domain
-    domainDest_->fillInRectilinearBoundLonLat(boundsLonDest, boundsLatDest);
     nVertexDest = constNVertex;
+    domainDest_->fillInRectilinearBoundLonLat(boundsLonDest, boundsLatDest, isNorthPole, isSouthPole);
   }
 
@@ -152 +181 @@
 
   std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
+  std::map<int,std::vector<std::pair<int,double> > >::const_iterator iteNorthPole = interpMapValueNorthPole.end(),
+                                                                     iteSouthPole = interpMapValueSouthPole.end();
   for (int idx = 0;  idx < mapper.nWeights; ++idx)
   {
     interpMapValue[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
-  }
+    if (iteNorthPole != interpMapValueNorthPole.find(mapper.targetWeightId[idx]))
+    {
+      interpMapValueNorthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
+    }
+
+    if (iteSouthPole != interpMapValueSouthPole.find(mapper.targetWeightId[idx]))
+    {
+      interpMapValueSouthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
+    }
+  }
+  int niGloDst = domainDest_->ni_glo.getValue();
+  processPole(interpMapValueNorthPole, niGloDst);
+  processPole(interpMapValueSouthPole, niGloDst);
+
+  if (!interpMapValueNorthPole.empty())
+  {
+     std::map<int,std::vector<std::pair<int,double> > >::iterator itNorthPole = interpMapValueNorthPole.begin();
+     for (; itNorthPole != iteNorthPole; ++itNorthPole)
+     {
+       if (!(itNorthPole->second.empty()))
+        itNorthPole->second.swap(interpMapValue[itNorthPole->first]);
+     }
+  }
+
+  if (!interpMapValueSouthPole.empty())
+  {
+     std::map<int,std::vector<std::pair<int,double> > >::iterator itSouthPole = interpMapValueSouthPole.begin();
+     for (; itSouthPole != iteSouthPole; ++itSouthPole)
+     {
+       if (!(itSouthPole->second.empty()))
+        itSouthPole->second.swap(interpMapValue[itSouthPole->first]);
+     }
+  }
+
   exchangeRemapInfo(interpMapValue);
 
   delete [] globalSrc;
   delete [] globalDst;
+}
+
+void CDomainAlgorithmInterpolate::processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole,
+                                              int nbGlobalPointOnPole)
+{
+  CContext* context = CContext::getCurrent();
+  CContextClient* client=context->client;
+
+  MPI_Comm poleComme(MPI_COMM_NULL);
+  MPI_Comm_split(client->intraComm, interMapValuePole.empty() ? MPI_UNDEFINED : 1, 0, &poleComme);
+  if (MPI_COMM_NULL != poleComme)
+  {
+    int nbClientPole;
+    MPI_Comm_size(poleComme, &nbClientPole);
+
+    std::map<int,std::vector<std::pair<int,double> > >::iterator itePole = interMapValuePole.end(), itPole,
+                                                                 itbPole = interMapValuePole.begin();
+
+    int nbWeight = 0;
+    for (itPole = itbPole; itPole != itePole; ++itPole)
+       nbWeight += itPole->second.size();
+
+    std::vector<int> recvCount(nbClientPole,0);
+    std::vector<int> displ(nbClientPole,0);
+    MPI_Allgather(&nbWeight,1,MPI_INT,&recvCount[0],1,MPI_INT,poleComme) ;
+
+    displ[0]=0;
+    for(int n=1;n<nbClientPole;++n) displ[n]=displ[n-1]+recvCount[n-1] ;
+    int recvSize=displ[nbClientPole-1]+recvCount[nbClientPole-1] ;
+
+    std::vector<int> sendSourceIndexBuff(nbWeight);
+    std::vector<double> sendSourceWeightBuff(nbWeight);
+    int k = 0;
+    for (itPole = itbPole; itPole != itePole; ++itPole)
+    {
+      for (int idx = 0; idx < itPole->second.size(); ++idx)
+      {
+        sendSourceIndexBuff[k] = (itPole->second)[idx].first;
+        sendSourceWeightBuff[k] = (itPole->second)[idx].second;
+        ++k;
+      }
+    }
+
+    std::vector<int> recvSourceIndexBuff(recvSize);
+    std::vector<double> recvSourceWeightBuff(recvSize);
+
+    // Gather all index and weight for pole
+    MPI_Allgatherv(&sendSourceIndexBuff[0],nbWeight,MPI_INT,&recvSourceIndexBuff[0],&recvCount[0],&displ[0],MPI_INT,poleComme);
+    MPI_Allgatherv(&sendSourceWeightBuff[0],nbWeight,MPI_DOUBLE,&recvSourceWeightBuff[0],&recvCount[0],&displ[0],MPI_DOUBLE,poleComme);
+
+    std::map<int,double> recvTemp;
+    for (int idx = 0; idx < recvSize; ++idx)
+    {
+      if (recvTemp.end() != recvTemp.find(recvSourceIndexBuff[idx]))
+        recvTemp[recvSourceIndexBuff[idx]] += recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
+      else
+        recvTemp[recvSourceIndexBuff[idx]] = 0.0;
+    }
+
+    std::map<int,double>::const_iterator itRecvTemp, itbRecvTemp = recvTemp.begin(), iteRecvTemp = recvTemp.end();
+
+    for (itPole = itbPole; itPole != itePole; ++itPole)
+    {
+      itPole->second.clear();
+      for (itRecvTemp = itbRecvTemp; itRecvTemp != iteRecvTemp; ++itRecvTemp)
+          itPole->second.push_back(make_pair(itRecvTemp->first, itRecvTemp->second));
+    }
+  }
+
 }
 

XIOS/trunk/src/transformation/domain_algorithm_interpolate.hpp

@@ -30 +30 @@
 private:
   void readInterpolationInfo(std::string& filename, std::map<int,std::vector<std::pair<int,double> > >& interpMapValue);
+  void processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole,
+                   int nbGlobalPointOnPole);
   void computeRemap();
   void readRemapInfo();