Changeset 1460 for XIOS/dev/branch_openmp/src/transformation/domain_algorithm_interpolate.cpp

Timestamp:

03/22/18 10:43:20 (6 years ago)

Author:

yushan

Message:

branch_openmp merged with XIOS_DEV_CMIP6@1459

File:

• XIOS/dev/branch_openmp/src/transformation/domain_algorithm_interpolate.cpp (modified) (10 diffs)

XIOS/dev/branch_openmp/src/transformation/domain_algorithm_interpolate.cpp

@@ -53 +53 @@
   CContext* context = CContext::getCurrent();
   interpDomain_->checkValid(domainSource);
+
+  detectMissingValue = interpDomain_->detect_missing_value ;
+  renormalize = interpDomain_->renormalize ;
+  quantity = interpDomain_->quantity ;
+
+  if (interpDomain_->read_write_convention == CInterpolateDomain::read_write_convention_attr::fortran) fortranConvention=true ;
+  else fortranConvention=false ;
+  
   fileToReadWrite_ = "xios_interpolation_weights_";
 
@@ -100 +108 @@
   std::vector<double> srcPole(3,0), dstPole(3,0);
   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;
@@ -456 +457 @@
         recvTemp[recvSourceIndexBuff[idx]] += recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
       else
-        recvTemp[recvSourceIndexBuff[idx]] = 0.0;
+        recvTemp[recvSourceIndexBuff[idx]] = recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
     }
 
@@ -476 +477 @@
 void CDomainAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
 {
-  if (readFromFile_ && !writeToFile_)  
+  if (readFromFile_)  
     readRemapInfo();
   else
@@ -566 +567 @@
   }
 
-  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp);
+  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp, client->serverSize);
   const CClientServerMapping::GlobalIndexMap& globalIndexInterpSendToClient = domainIndexClientClientMapping.getGlobalIndexOnServer();
 
@@ -775 +776 @@
 
   int index = 0;
+  int indexOffset=0 ;
+  if (fortranConvention) indexOffset=1 ;
   for (it = itb; it !=ite; ++it)
   {
@@ -780 +783 @@
     for (int idx = 0; idx < tmp.size(); ++idx)
     {
-      dst_idx(index) = it->first + 1;
-      src_idx(index) = tmp[idx].first + 1;
+      dst_idx(index) = it->first + indexOffset;
+      src_idx(index) = tmp[idx].first + indexOffset;
       weights(index) = tmp[idx].second;
       ++index;
@@ -843 +846 @@
   size_t nbWeightGlo ;
 
+
   CContext* context = CContext::getCurrent();
   CContextClient* client=context->client;
@@ -848 +852 @@
   int clientSize = client->clientSize;
 
+
+  {
+    ifstream f(filename.c_str());
+    if (!f.good()) ERROR("void CDomainAlgorithmInterpolate::readInterpolationInfo",
+                      << "Attempt to read file weight :"  << filename << " which doesn't seem to exist." << std::endl
+                      << "Please check this file ");
+  }
+                  
   nc_open(filename.c_str(),NC_NOWRITE, &ncid) ;
   nc_inq_dimid(ncid,"n_weight",&weightDimId) ;
@@ -882 +894 @@
   nc_get_vara_long(ncid, dstIndexId, &start, &nbWeight, dstIndex) ;
 
-  for(size_t ind=0; ind<nbWeight;++ind)
-    interpMapValue[dstIndex[ind]-1].push_back(make_pair(srcIndex[ind]-1,weight[ind]));
-}
-
-}
+  int indexOffset=0 ;
+  if (fortranConvention) indexOffset=1 ;
+    for(size_t ind=0; ind<nbWeight;++ind)
+      interpMapValue[dstIndex[ind]-indexOffset].push_back(make_pair(srcIndex[ind]-indexOffset,weight[ind]));
+ }
+
+void CDomainAlgorithmInterpolate::apply(const std::vector<std::pair<int,double> >& localIndex,
+                                            const double* dataInput,
+                                            CArray<double,1>& dataOut,
+                                            std::vector<bool>& flagInitial,
+                                            bool ignoreMissingValue, bool firstPass  )
+{
+  int nbLocalIndex = localIndex.size();   
+  double defaultValue = std::numeric_limits<double>::quiet_NaN();
+   
+  if (detectMissingValue)
+  {
+     if (firstPass && renormalize)
+     {
+       renormalizationFactor.resize(dataOut.numElements()) ;
+       renormalizationFactor=1 ;
+     }
+
+    for (int idx = 0; idx < nbLocalIndex; ++idx)
+    {
+      if (NumTraits<double>::isnan(*(dataInput + idx)))
+      {
+        flagInitial[localIndex[idx].first] = false;
+        if (renormalize) renormalizationFactor(localIndex[idx].first)-=localIndex[idx].second ;
+      }
+      else
+      {
+        dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
+        flagInitial[localIndex[idx].first] = true; // Reset flag to indicate not all data source are nan
+      }
+    }
+
+    // If all data source are nan then data destination must be nan
+    for (int idx = 0; idx < nbLocalIndex; ++idx)
+    {
+      if (!flagInitial[localIndex[idx].first])
+        dataOut(localIndex[idx].first) = defaultValue;
+    }
+  }
+  else
+  {
+    for (int idx = 0; idx < nbLocalIndex; ++idx)
+    {
+      dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
+    }
+  }
+}
+
+void CDomainAlgorithmInterpolate::updateData(CArray<double,1>& dataOut)
+{
+  if (detectMissingValue && renormalize)
+  {
+    if (renormalizationFactor.numElements()>0) dataOut/=renormalizationFactor ; // In some case, process doesn't received any data for interpolation (mask)
+                                                                                // so renormalizationFactor is not initialized
+  }
+}
+
+}