Changeset 827 for XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp

Timestamp:

03/23/16 16:10:45 (8 years ago)

Author:

mhnguyen

Message:

Implementing dynamic interpolation on axis

+) Change grid transformation to make it more flexible
+) Make some small improvements

Test
+) On Curie
+) All test pass

File:

• XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp (modified) (10 diffs)

XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp

@@ -12 +12 @@
 #include "context_client.hpp"
 #include "utils.hpp"
+#include "grid.hpp"
+#include "distribution_client.hpp"
 
 namespace xios {
 
 CAxisAlgorithmInterpolate::CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis)
-: CAxisAlgorithmTransformation(axisDestination, axisSource)
+: CAxisAlgorithmTransformation(axisDestination, axisSource), coordinate_(), transPosition_()
 {
   interpAxis->checkValid(axisSource);
   order_ = interpAxis->order.getValue();
-  if (order_ >= axisSource->n_glo.getValue())
-  {
-    ERROR("CAxisAlgorithmInterpolate::CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis)",
-           << "Order of interpolation is greater than global size of axis source"
-           << "Size of axis source " <<axisSource->getId() << " is " << axisSource->n_glo.getValue()  << std::endl
-           << "Order of interpolation is " << order_ );
-  }
-
-  computeIndexSourceMapping();
+  if (!interpAxis->coordinate.isEmpty())
+  {
+    coordinate_ = interpAxis->coordinate.getValue();
+    this->idAuxInputs_.resize(1);
+    this->idAuxInputs_[0] = coordinate_;
+  }
+
+//  computeIndexSourceMapping();
 }
 
@@ -34 +35 @@
   Compute the index mapping between axis on grid source and one on grid destination
 */
-void CAxisAlgorithmInterpolate::computeIndexSourceMapping()
-{
-  CArray<double,1>& axisValue = axisSrc_->value;
-  CArray<bool,1>& axisMask = axisSrc_->mask;
-
+void CAxisAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
   CContext* context = CContext::getCurrent();
   CContextClient* client=context->client;
   int nbClient = client->clientSize;
-
+  CArray<bool,1>& axisMask = axisSrc_->mask;
   int srcSize  = axisSrc_->n_glo.getValue();
-  int numValue = axisValue.numElements();
-
-  std::vector<double> recvBuff(srcSize);
-  std::vector<int> indexVec(srcSize);
-
-  retrieveAllAxisValue(recvBuff, indexVec);
-  XIOSAlgorithms::sortWithIndex<double, CVectorStorage>(recvBuff, indexVec);
-  computeInterpolantPoint(recvBuff, indexVec);
+  std::vector<CArray<double,1> > vecAxisValue;
+
+  // Fill in axis value from coordinate
+  fillInAxisValue(vecAxisValue, dataAuxInputs);
+
+  for (int idx = 0; idx < vecAxisValue.size(); ++idx)
+  {
+    CArray<double,1>& axisValue = vecAxisValue[idx];
+    std::vector<double> recvBuff(srcSize);
+    std::vector<int> indexVec(srcSize);
+    retrieveAllAxisValue(axisValue, axisMask, recvBuff, indexVec);
+    XIOSAlgorithms::sortWithIndex<double, CVectorStorage>(recvBuff, indexVec);
+    computeInterpolantPoint(recvBuff, indexVec, idx);
+  }
 }
 
@@ -60 +64 @@
   \param [in] indexVec permutation index of axisValue
 */
-void CAxisAlgorithmInterpolate::computeInterpolantPoint(const std::vector<double>& axisValue, const std::vector<int>& indexVec)
+void CAxisAlgorithmInterpolate::computeInterpolantPoint(const std::vector<double>& axisValue, const std::vector<int>& indexVec, int transPos)
 {
   std::vector<double>::const_iterator itb = axisValue.begin(), ite = axisValue.end();
@@ -78 +82 @@
     itUpperBound = std::upper_bound(itb, ite, destValue);
     if ((ite != itUpperBound) && (sfmax == *itUpperBound)) itUpperBound = ite;
-
 
     // If the value is not in the range, that means we'll do extra-polation
@@ -118 +121 @@
     }
   }
-  computeWeightedValueAndMapping(interpolatingIndexValues);
+  computeWeightedValueAndMapping(interpolatingIndexValues, transPos);
 }
 
@@ -125 +128 @@
   \param [in] interpolatingIndexValues the necessary axis value to calculate the coeffs
 */
-void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues)
-{
-  std::map<int, std::vector<int> >& transMap = this->transformationMapping_;
-  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_;
+void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos)
+{
+  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[transPos];
+  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[transPos];
   std::map<int, std::vector<std::pair<int,double> > >::const_iterator itb = interpolatingIndexValues.begin(), it,
                                                                       ite = interpolatingIndexValues.end();
@@ -138 +141 @@
     const std::vector<std::pair<int,double> >& interpVal = it->second;
     int interpSize = interpVal.size();
+    transMap[globalIndexDest].resize(interpSize);
+    transWeight[globalIndexDest].resize(interpSize);
     for (int idx = 0; idx < interpSize; ++idx)
     {
@@ -149 +154 @@
         weight /= (interpVal[idx].second - interpVal[k].second);
       }
-      transMap[globalIndexDest].push_back(index);
-      transWeight[globalIndexDest].push_back(weight);
+      transMap[globalIndexDest][idx] = index;
+      transWeight[globalIndexDest][idx] = weight;
+      if (!transPosition_.empty())
+      {
+        (this->transformationPosition_[transPos])[globalIndexDest] = transPosition_[transPos];
+      }
     }
   }
@@ -160 +169 @@
   \param [in/out] indexVec mapping between values and global index of axis
 */
-void CAxisAlgorithmInterpolate::retrieveAllAxisValue(std::vector<double>& recvBuff, std::vector<int>& indexVec)
-{
-  CArray<double,1>& axisValue = axisSrc_->value;
-  CArray<bool,1>& axisMask = axisSrc_->mask;
-
+void CAxisAlgorithmInterpolate::retrieveAllAxisValue(const CArray<double,1>& axisValue, const CArray<bool,1>& axisMask,
+                                                     std::vector<double>& recvBuff, std::vector<int>& indexVec)
+{
   CContext* context = CContext::getCurrent();
   CContextClient* client=context->client;
@@ -230 +237 @@
 }
 
-}
+/*!
+  Fill in axis value dynamically from a field whose grid is composed of a domain and an axis
+  \param [in/out] vecAxisValue vector axis value filled in from input field
+*/
+void CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue,
+                                                const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+  if (coordinate_.empty())
+  {
+    vecAxisValue.resize(1);
+    vecAxisValue[0].resize(axisSrc_->value.numElements());
+    vecAxisValue[0] = axisSrc_->value;
+    this->transformationMapping_.resize(1);
+    this->transformationWeight_.resize(1);
+  }
+  else
+  {
+    CField* field = CField::get(coordinate_);
+    CGrid* grid = field->grid;
+
+    std::vector<CDomain*> domListP = grid->getDomains();
+    std::vector<CAxis*> axisListP = grid->getAxis();
+    if (domListP.empty() || axisListP.empty() || (1 < domListP.size()) || (1 < axisListP.size()))
+      ERROR("CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue)",
+             << "XIOS only supports dynamic interpolation with coordinate (field) associated with grid composed of a domain and an axis"
+             << "Coordinate (field) id = " <<field->getId() << std::endl
+             << "Associated grid id = " << grid->getId());
+
+    CDomain* dom = domListP[0];
+    size_t vecAxisValueSize = dom->i_index.numElements();
+    vecAxisValue.resize(vecAxisValueSize);
+    if (transPosition_.empty())
+    {
+      transPosition_.resize(vecAxisValueSize);
+      for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
+      {
+        transPosition_[idx].resize(2);
+        transPosition_[idx][0] = (dom->i_index)(idx);
+        transPosition_[idx][1] = (dom->j_index)(idx);
+      }
+    }
+    this->transformationMapping_.resize(vecAxisValueSize);
+    this->transformationWeight_.resize(vecAxisValueSize);
+    this->transformationPosition_.resize(vecAxisValueSize);
+
+    const std::vector<size_t>& globalIndexSendToServer = grid->getDistributionClient()->getGlobalDataIndexSendToServer();
+    const std::vector<int>& localDataSendToServer = grid->getDistributionClient()->getLocalDataIndexSendToServer();
+    std::vector<int> permutIndex(globalIndexSendToServer.size());
+    std::vector<int>::iterator itVec;
+    XIOSAlgorithms::fillInIndex(globalIndexSendToServer.size(), permutIndex);
+    XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexSendToServer, permutIndex);
+    size_t axisSrcSize = axisSrc_->index.numElements();
+    std::vector<int> globalDimension = grid->getGlobalDimension();
+    XIOSBinarySearchWithIndex<size_t> binSearch(globalIndexSendToServer);
+    for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
+    {
+      size_t axisValueSize = 0;
+      for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
+      {
+        size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
+        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndex, itVec))
+        {
+          ++axisValueSize;
+        }
+      }
+
+      vecAxisValue[idx].resize(axisValueSize);
+      axisValueSize = 0;
+      for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
+      {
+        size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
+        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndex, itVec))
+        {
+          vecAxisValue[idx](axisValueSize) = (*dataAuxInputs[0])(localDataSendToServer[*itVec]);
+          ++axisValueSize;
+        }
+      }
+    }
+  }
+}
+
+}