Changeset 2313

Timestamp:

03/21/22 15:27:07 (2 years ago)

Author:

ymipsl

Message:

Take into account detect_missing_value and renormalize attribute for interpolate class transformation.

YM

Location:

XIOS/dev/dev_ym/XIOS_COUPLING/src

Files:

• distribution/weight_transform_connector.cpp (modified) (1 diff)
• distribution/weight_transform_connector.hpp (modified) (3 diffs)
• transformation/algorithm_transformation_weight.cpp (modified) (1 diff)
• transformation/algorithm_transformation_weight.hpp (modified) (1 diff)
• transformation/axis_algorithm/axis_algorithm_interpolate.cpp (modified) (1 diff)
• transformation/domain_algorithm/domain_algorithm_interpolate.cpp (modified) (1 diff)

XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/weight_transform_connector.cpp

@@ -5 +5 @@
 
   CWeightTransformConnector::CWeightTransformConnector(shared_ptr<CLocalView> srcView, shared_ptr<CLocalView> dstView, unordered_map<int, std::vector<int>>& indexMap, 
-                                                   unordered_map<int, std::vector<double>>& weightMap) : srcView_(srcView), dstView_(dstView)
+                                                       unordered_map<int, std::vector<double>>& weightMap,  bool detectMissingValue, bool renormalize) : 
+                                                       srcView_(srcView), dstView_(dstView), detectMissingValue_(detectMissingValue), renormalize_(renormalize)
   {
     computeConnector(indexMap, weightMap) ;

XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/weight_transform_connector.hpp

@@ -23 +23 @@
       int srcSize_ ;
       int dstSize_ ;
+      bool detectMissingValue_ ;
+      bool renormalize_ ;
 
      void computeConnector(unordered_map<int, std::vector<int>>& indexMap, unordered_map<int, std::vector<double>>& weightMap) ;
@@ -29 +31 @@
 
     CWeightTransformConnector(shared_ptr<CLocalView> srcView, shared_ptr<CLocalView> dstView, unordered_map<int, std::vector<int>>& indexMap, 
-                              unordered_map<int, std::vector<double>>& weightMap) ;
+                              unordered_map<int, std::vector<double>>& weightMap, bool detectMissingValue, bool renormalize) ;
  
     void transfer(int repeat, int sizeT, const CArray<double,1>& dataIn, CArray<double,1>& dataOut)
@@ -36 +38 @@
       int srcSlice = srcSize_*sizeT ;
       dataOut.resize(repeat* dstSlice) ;
-      dataOut=0 ; // what to do about missing value => next step ?
-
+      double defaultValue = std::numeric_limits<double>::quiet_NaN();
+      dataOut = defaultValue ; // what to do about missing value => next step ?
       const double* input = dataIn.dataFirst()  ;
       double* output = dataOut.dataFirst()  ;
+      vector<bool> isFirst(dstSlice) ;
+      size_t pos ;
+      
+      if (renormalize_)
+      {
+        double inVal ;
+        vector<double> renormalizeFactor(repeat*dstSlice,1) ;
+        double* renorm = renormalizeFactor.data() ;
+        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, renorm+=dstSlice)
+        {
+          const double* in = input;
+          double* out = output ;
+          double* ren = renorm ;
+          isFirst.assign(dstSlice,true) ;
+          pos=0 ;
+          int k=0 ;
+          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT, ren+=sizeT)
+            for(int j=0 ; j<nWeights_[i] ; j++,k++)
+              for(int l=0; l<sizeT; l++)
+              {
+                inVal=in[connector_[k]*sizeT+l] ;
+                if (!std::isnan(inVal))
+                {
+                  if (isFirst[pos+l])  { out[l] = inVal*weights_[k] ; isFirst[pos+l]=false ;}
+                  else out[l] += inVal*weights_[k] ;
+                }
+                else  ren[l] -= weights_[k] ;
+              }
+          for(int i=0; i<dstSlice; i++) output[i] /= renorm[i] ; 
+        }
+      }
+      else if (detectMissingValue_)
+      {
+        double inVal ;
+        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice)
+        {
+          const double* in = input;
+          double* out = output ;
+          isFirst.assign(dstSlice,true) ;
+          pos=0 ;
+          int k=0 ;
+          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
+            for(int j=0 ; j<nWeights_[i] ; j++,k++)
+              for(int l=0; l<sizeT; l++)
+              {
+                inVal=in[connector_[k]*sizeT+l] ;
+                if (!std::isnan(inVal))
+                {
+                  if (isFirst[pos+l])  { out[l] = inVal*weights_[k] ; isFirst[pos+l]=false ;}
+                  else out[l] += inVal*weights_[k] ;
+                }
+              }
+        }
 
-      for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice)
+      }
+      else
       {
-        const double* in = input;
-        double* out = output ;
-        int k=0 ;
-        for(int i=0; i<dstSize_; i++, out+=sizeT)
-          for(int j=0 ; j<nWeights_[i] ; j++,k++)
-            for(int l=0; l<sizeT; l++) out[l] += in[connector_[k]*sizeT+l]*weights_[k] ;
+        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice)
+        {
+          const double* in = input;
+          double* out = output ;
+          isFirst.assign(dstSlice,true) ;
+          pos=0 ;
+          int k=0 ;
+          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
+            for(int j=0 ; j<nWeights_[i] ; j++,k++)
+              for(int l=0; l<sizeT; l++) 
+                if (isFirst[pos+l])  { out[l] = in[connector_[k]*sizeT+l]*weights_[k] ; isFirst[pos+l]=false ;}
+                else out[l] += in[connector_[k]*sizeT+l]*weights_[k] ;
+        }
       }
     }

XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/algorithm_transformation_weight.cpp

@@ -6 +6 @@
 
 
-  void CAlgorithmTransformationWeight::computeAlgorithm(shared_ptr<CLocalView> srcView, shared_ptr<CLocalView> dstView)
+  void CAlgorithmTransformationWeight::computeAlgorithm(shared_ptr<CLocalView> srcView, shared_ptr<CLocalView> dstView, bool detectMissingValue, bool renormalize)
   {
     this->computeRecvElement(srcView, dstView) ;
-    weightTransformConnector_ = make_shared<CWeightTransformConnector>(recvElement_->getView(CElementView::FULL), dstView, transformationMapping_, transformationWeight_) ; 
+    weightTransformConnector_ = make_shared<CWeightTransformConnector>(recvElement_->getView(CElementView::FULL), dstView, transformationMapping_, transformationWeight_, detectMissingValue, renormalize) ; 
   }
  

XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/algorithm_transformation_weight.hpp

@@ -15 +15 @@
     public:
 
-      CAlgorithmTransformationWeight(bool isSource) : CGenericAlgorithmTransformation(isSource) {}
+      CAlgorithmTransformationWeight(bool isSource) : CGenericAlgorithmTransformation(isSource) {} 
+                                                                                                 
       virtual ~CAlgorithmTransformationWeight() {};
       virtual void apply(int dimBefore, int dimAfter, const CArray<double,1>& dataIn, CArray<double,1>& dataOut);
       virtual void computeRecvElement(shared_ptr<CLocalView> srcView, shared_ptr<CLocalView> dstView);
     protected:
-      virtual void computeAlgorithm(shared_ptr<CLocalView> srcView, shared_ptr<CLocalView> dstView) ;
+      virtual void computeAlgorithm(shared_ptr<CLocalView> srcView, shared_ptr<CLocalView> dstView, bool detectMissingValue, bool renormalize) ;
 
       //! Map between global index of destination element and source element

XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate.cpp

@@ -68 +68 @@
   std::vector<CArray<double,1>* > dataAuxInputs ;
   computeRemap(dataAuxInputs) ;
-  this->computeAlgorithm(axisSource->getLocalView(CElementView::WORKFLOW), axisDestination->getLocalView(CElementView::WORKFLOW)) ;
+  this->computeAlgorithm(axisSource->getLocalView(CElementView::WORKFLOW), axisDestination->getLocalView(CElementView::WORKFLOW), false, false) ;
 }
 CATCH

XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_interpolate.cpp

@@ -102 +102 @@
   else computeRemap(); 
 
-    this->computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), domainDestination->getLocalView(CElementView::WORKFLOW)) ;
+    this->computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), domainDestination->getLocalView(CElementView::WORKFLOW), detectMissingValue, renormalize) ;
     
 }