source: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_domain.cpp @ 2291

Last change on this file since 2291 was 2291, checked in by ymipsl, 2 years ago

Improve reduction transformation

  • make the difference between reduction over geometry or reduction between process.
  • geometrical reduction :

domain -> axis
axis -> scalar
domain -> scalar

  • reduction across processes for redondant geometrical cell :

axis -> axis
scalar -> scalar

Reduction can be local (only for the geometrical cell owned by current process) or global, using the "local" attribute (bool) over the reduction.

YM
  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 12.5 KB
Line 
1/*!
2   \file axis_algorithm_reduce_domain.cpp
3   \author Ha NGUYEN
4   \since 23 June 2016
5   \date 23 June 2016
6
7   \brief Algorithm for reduce a domain to an axis
8 */
9#include "axis_algorithm_reduce_domain.hpp"
10#include "reduce_domain_to_axis.hpp"
11#include "axis.hpp"
12#include "domain.hpp"
13#include "grid.hpp"
14#include "grid_transformation_factory_impl.hpp"
15#include "reduction.hpp"
16
17namespace xios {
18shared_ptr<CGenericAlgorithmTransformation> CAxisAlgorithmReduceDomain::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
19                                                                   CTransformation<CAxis>* transformation,
20                                                                   int elementPositionInGrid,
21                                                                   std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
22                                                                   std::map<int, int>& elementPositionInGridSrc2AxisPosition,
23                                                                   std::map<int, int>& elementPositionInGridSrc2DomainPosition,
24                                                                   std::map<int, int>& elementPositionInGridDst2ScalarPosition,
25                                                                   std::map<int, int>& elementPositionInGridDst2AxisPosition,
26                                                                   std::map<int, int>& elementPositionInGridDst2DomainPosition)
27TRY
28{
29  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
30  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
31
32  CReduceDomainToAxis* reduceDomain = dynamic_cast<CReduceDomainToAxis*> (transformation);
33  int axisDstIndex   = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
34  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
35
36  return make_shared<CAxisAlgorithmReduceDomain>(isSource, axisListDestP[axisDstIndex], domainListSrcP[domainSrcIndex], reduceDomain);
37}
38CATCH
39
40bool CAxisAlgorithmReduceDomain::dummyRegistered_ = CAxisAlgorithmReduceDomain::registerTrans();
41bool CAxisAlgorithmReduceDomain::registerTrans()
42TRY
43{
44  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_REDUCE_DOMAIN_TO_AXIS, create);
45}
46CATCH
47
48
49CAxisAlgorithmReduceDomain::CAxisAlgorithmReduceDomain(bool isSource, CAxis* axisDestination, CDomain* domainSource, CReduceDomainToAxis* algo)
50 : CAlgorithmTransformationReduce(isSource), domainSrc_(domainSource), axisDest_(axisDestination)
51TRY
52{
53  switch (algo->operation)
54  {
55    case CReduceDomainToAxis::operation_attr::sum:
56       operator_ = EReduction::sum;
57      break;
58    case CReduceDomainToAxis::operation_attr::min:
59       operator_ = EReduction::min;
60      break;
61    case CReduceDomainToAxis::operation_attr::max:
62       operator_ = EReduction::max;
63      break;
64    case CReduceDomainToAxis::operation_attr::average:
65       operator_ = EReduction::average;
66      break;
67    default:
68        ERROR("CAxisAlgorithmReduceDomain::CAxisAlgorithmReduceDomain(CAxis* axisDestination, CDomain* domainSource, CReduceDomainToAxis* algo)",
69         << "Operation is wrongly defined. Supported operations: sum, min, max, average." << std::endl
70         << "Domain source " <<domainSource->getId() << std::endl
71         << "Axis destination " << axisDestination->getId());
72
73  }
74 
75  algo->checkValid(axisDestination, domainSource);
76  dir_ = (CReduceDomainToAxis::direction_attr::iDir == algo->direction)  ? iDir : jDir;
77
78  bool local = false;
79  if (!algo->local.isEmpty()) local=algo->local ;
80 
81  size_t nj_glo = domainSource->nj_glo ;
82  size_t ni_glo = domainSource->ni_glo ;
83
84  bool validAxis = axisDestination->checkGeometricAttributes(false) ;
85  if (validAxis && !local) 
86  {
87   
88    axisDestination->checkAttributes() ; 
89    if (dir_==jDir)
90    {
91        if (axisDestination->n_glo != domainSource->nj_glo)
92          ERROR("CAxisAlgorithmReduceDomain::CAxisAlgorithmReduceDomain(bool isSource, CAxis* axisDestination, CDomain* domainSource, CReduceDomainToAxis* algo)",
93            << "Extract domain along j, axis destination should have n_glo equal to nj_glo of domain source"
94            << "Domain source " <<domainSource->getId() << " has nj_glo " << domainSource->nj_glo << std::endl
95            << "Axis destination " << axisDestination->getId() << " has n_glo " << axisDestination->n_glo)
96    }
97    else
98    {
99      if (axisDestination->n_glo != domainSource->ni_glo)
100          ERROR("CAxisAlgorithmReduceDomain::CAxisAlgorithmReduceDomain(bool isSource, CAxis* axisDestination, CDomain* domainSource, CReduceDomainToAxis* algo)",
101            << "Extract domain along j, axis destination should have n_glo equal to ni_glo of domain source"
102            << "Domain source " <<domainSource->getId() << " has ni_glo " << domainSource->ni_glo << std::endl
103            << "Axis destination " << axisDestination->getId() << " has n_glo " << axisDestination->n_glo);
104    }
105  }
106  else
107  {
108    // create axis destination from source domain
109
110    axisDestination->resetGeometricAttributes();
111    if (dir_== jDir)
112    {
113      axisDestination->n_glo = domainSource->nj_glo ;
114     
115      CArray<size_t, 1> srcGlobalIndex ;
116      set<size_t> indexFullView;
117      set<size_t> indexWorkflowView;
118      domainSource->getLocalView(CElementView::FULL)->getGlobalIndexView(srcGlobalIndex) ;
119      for(int i=0; i<srcGlobalIndex.numElements(); i++) indexFullView.insert(srcGlobalIndex(i)/ni_glo) ;
120      domainSource->getLocalView(CElementView::WORKFLOW)->getGlobalIndexView(srcGlobalIndex) ;
121      for(int i=0; i<srcGlobalIndex.numElements(); i++) indexWorkflowView.insert(srcGlobalIndex(i)/ni_glo) ;
122
123      axisDestination-> n = indexFullView.size() ;
124      axisDestination-> index.resize(axisDestination-> n) ;
125      axisDestination-> mask.resize(axisDestination-> n) ;
126     
127      map<size_t,int> globalToLocalIndex ;
128      auto it=indexFullView.begin();
129      for(int i=0; it!=indexFullView.end(); ++i,++it) 
130      {
131        axisDestination->index(i) = *it ;
132        if (indexWorkflowView.count(*it)==0) axisDestination->mask(i) = false ;
133        else axisDestination->mask(i) = true ;
134        globalToLocalIndex[*it] = i ;
135      }
136      if (domainSource->hasLonLat && domainSource->type == CDomain::type_attr::rectilinear)
137      {
138        axisDestination->value.resize(axisDestination-> n) ;
139        axisDestination->axis_type.setValue(CAxis::axis_type_attr::Y) ;
140        if (domainSource->hasBounds) axisDestination-> bounds.resize(axisDestination-> n, 2) ;
141
142        domainSource->getLocalView(CElementView::FULL)->getGlobalIndexView(srcGlobalIndex) ;
143        for(int i=0; i<srcGlobalIndex.numElements(); i++) 
144        {
145          axisDestination->value(globalToLocalIndex[srcGlobalIndex(i)/ni_glo]) = domainSource->latvalue(i) ;
146          if (domainSource->hasBounds) 
147          {
148            axisDestination->bounds(globalToLocalIndex[srcGlobalIndex(i)/ni_glo,0]) = domainSource->bounds_latvalue(i,0) ;
149            axisDestination->bounds(globalToLocalIndex[srcGlobalIndex(i)/ni_glo,1]) = domainSource->bounds_latvalue(i,1) ;
150          }
151        } 
152      }
153     
154    }
155    else // dir_== iDir
156    {
157      axisDestination->n_glo = domainSource->ni_glo ;
158     
159      CArray<size_t, 1> srcGlobalIndex ;
160      set<size_t> indexFullView;
161      set<size_t> indexWorkflowView;
162      domainSource->getLocalView(CElementView::FULL)->getGlobalIndexView(srcGlobalIndex) ;
163      for(int i=0; i<srcGlobalIndex.numElements(); i++) indexFullView.insert(srcGlobalIndex(i)%ni_glo) ;
164      domainSource->getLocalView(CElementView::WORKFLOW)->getGlobalIndexView(srcGlobalIndex) ;
165      for(int i=0; i<srcGlobalIndex.numElements(); i++) indexWorkflowView.insert(srcGlobalIndex(i)%ni_glo) ;
166
167      axisDestination-> n = indexFullView.size() ;
168      axisDestination-> index.resize(axisDestination-> n) ;
169      axisDestination-> mask.resize(axisDestination-> n) ;
170      map<size_t,int> globalToLocalIndex ;
171      auto it=indexFullView.begin();
172      for(int i=0; it!=indexFullView.end(); ++i,++it) 
173      {
174        axisDestination->index(i) = *it ;
175        if (indexWorkflowView.count(*it)==0) axisDestination->mask(i) = false ;
176        else axisDestination->mask(i) = true ;
177        globalToLocalIndex[*it] = i ;
178      }
179   
180      if (domainSource->hasLonLat && domainSource->type == CDomain::type_attr::rectilinear)
181      {
182        axisDestination-> value.resize(axisDestination-> n) ;
183        axisDestination-> axis_type.setValue(CAxis::axis_type_attr::X) ;
184        if (domainSource->hasBounds) axisDestination->bounds.resize(axisDestination-> n, 2) ;
185
186        domainSource->getLocalView(CElementView::FULL)->getGlobalIndexView(srcGlobalIndex) ;
187        for(int i=0; i<srcGlobalIndex.numElements(); i++) 
188        {
189          axisDestination->value(globalToLocalIndex[srcGlobalIndex(i)%ni_glo]) = domainSource->lonvalue(i) ;
190          if (domainSource->hasBounds) 
191          {
192            axisDestination->bounds(globalToLocalIndex[srcGlobalIndex(i)%ni_glo,0]) = domainSource->bounds_lonvalue(i,0) ;
193            axisDestination->bounds(globalToLocalIndex[srcGlobalIndex(i)%ni_glo,1]) = domainSource->bounds_lonvalue(i,1) ;
194          }
195        } 
196      }
197   
198    }
199   
200    axisDestination->checkAttributes() ;   
201  }
202 
203  // compute needed index for tranformation
204
205  TransformationIndexMap& transMap = transformationMapping_;
206
207  CArray<size_t, 1> srcGlobalIndex ;
208  domainSource->getLocalView(CElementView::WORKFLOW)->getGlobalIndexView(srcGlobalIndex) ;
209  CArray<size_t, 1> dstGlobalIndex ;
210  axisDestination->getLocalView(CElementView::WORKFLOW)->getGlobalIndexView(dstGlobalIndex) ;
211 
212  if (local)
213  {
214   
215    set<size_t> dstGlobalIndexSet ;
216    for(int i=0; i<dstGlobalIndex.numElements() ; i++) dstGlobalIndexSet.insert(dstGlobalIndex(i)) ; // for mask
217
218    size_t dstInd ; 
219    for(int i=0; i<srcGlobalIndex.numElements(); i++) 
220    {
221      if (dir_== jDir) dstInd=srcGlobalIndex(i)/ni_glo ;
222      else dstInd=srcGlobalIndex(i)%ni_glo ;
223      if (dstGlobalIndexSet.count(dstInd)!=0) transMap[dstInd].push_back(srcGlobalIndex(i)) ; 
224    }
225  }
226  else
227  {
228    for(int i=0; i<dstGlobalIndex.numElements() ;i++)
229    {
230      if (dir_== jDir) 
231        for(size_t j=0; j<ni_glo;j++) transMap[dstGlobalIndex(i)].push_back(dstGlobalIndex(i)*ni_glo+j) ;
232      else
233        for(size_t j=0; j<nj_glo;j++) transMap[dstGlobalIndex(i)].push_back(ni_glo*j+dstGlobalIndex(i)) ;
234    }
235  }
236
237/*
238
239
240  TransformationIndexMap& transMap = transformationMapping_;
241
242  CArray<int,1>& axisDstIndex = axisDest_->index;
243  int ni_glo = domainSrc_->ni_glo, nj_glo = domainSrc_->nj_glo;
244  if (iDir == dir_)
245  {
246    if (local)
247    {
248      const CArray<int, 1>& i_index = domainSrc_-> i_index.getValue() ;
249      const CArray<int, 1>& j_index = domainSrc_-> j_index.getValue() ;
250      const CArray<bool,1>& localMask = domainSrc_-> localMask ;
251      int nbDomainIdx = i_index.numElements();
252     
253      for (int idxDomain = 0; idxDomain < nbDomainIdx; ++idxDomain)
254      {
255        if (localMask(idxDomain))
256        {
257          transMap[j_index(idxDomain)].push_back(j_index(idxDomain)* ni_glo + i_index(idxDomain));
258        }
259      }
260    }
261    else
262    {
263      int nbAxisIdx = axisDstIndex.numElements();
264      for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
265      {
266        int globalAxisIdx = axisDstIndex(idxAxis);
267        transMap[globalAxisIdx].resize(ni_glo);
268        for (int idx = 0; idx < ni_glo; ++idx)
269        {
270          transMap[globalAxisIdx][idx] = globalAxisIdx * ni_glo + idx;
271        }
272      }
273    }
274  }
275  else if (jDir == dir_)
276  {
277    int nbAxisIdx = axisDstIndex.numElements();
278    if (local)
279    {
280      const CArray<int, 1>& i_index = domainSrc_-> i_index.getValue() ;
281      const CArray<int, 1>& j_index = domainSrc_-> j_index.getValue() ;
282      const CArray<bool,1>& localMask = domainSrc_-> localMask ;
283      int nbDomainIdx = i_index.numElements();
284     
285      for (int idxDomain = 0; idxDomain < nbDomainIdx; ++idxDomain)
286      {
287        if (localMask(idxDomain))
288        {
289          transMap[i_index(idxDomain)].push_back(j_index(idxDomain)* ni_glo + i_index(idxDomain));
290        }
291      }
292    }
293    else
294    {
295      for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
296      {
297        int globalAxisIdx = axisDstIndex(idxAxis);
298        transMap[globalAxisIdx].resize(nj_glo);
299        for (int idx = 0; idx < nj_glo; ++idx)
300        {
301          transMap[globalAxisIdx][idx] = globalAxisIdx + ni_glo*idx;
302        }
303      }
304    }
305  }
306  else
307  {}
308
309  axisDestination->checkAttributes() ;
310*/
311  this->computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), axisDestination->getLocalView(CElementView::WORKFLOW)) ;
312
313}
314CATCH
315
316
317CAxisAlgorithmReduceDomain::~CAxisAlgorithmReduceDomain()
318TRY
319{
320}
321CATCH
322
323
324}
Note: See TracBrowser for help on using the repository browser.