source: XIOS/trunk/src/transformation/axis_algorithm_inverse.cpp @ 827

Last change on this file since 827 was 827, checked in by mhnguyen, 5 years ago

Implementing dynamic interpolation on axis

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

Test
+) On Curie
+) All test pass

File size: 5.8 KB
Line 
1/*!
2   \file axis_algorithm_inverse.hpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 29 June 2015
6
7   \brief Algorithm for inversing an axis..
8 */
9#include "axis_algorithm_inverse.hpp"
10#include "transformation_mapping.hpp"
11#include "context.hpp"
12#include "context_client.hpp"
13
14namespace xios {
15
16CAxisAlgorithmInverse::CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource)
17 : CAxisAlgorithmTransformation(axisDestination, axisSource)
18{
19  if (axisDestination->n_glo.getValue() != axisSource->n_glo.getValue())
20  {
21    ERROR("CAxisAlgorithmInverse::CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource)",
22           << "Two axis have different global size"
23           << "Size of axis source " <<axisSource->getId() << " is " << axisSource->n_glo.getValue()  << std::endl
24           << "Size of axis destionation " <<axisDestination->getId() << " is " << axisDestination->n_glo.getValue());
25  }
26
27//  this->computeIndexSourceMapping();
28}
29
30void CAxisAlgorithmInverse::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
31{
32  this->transformationMapping_.resize(1);
33  this->transformationWeight_.resize(1);
34
35  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0];
36  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[0];
37
38  int globalIndexSize = axisDestGlobalIndex_.size();
39  for (int idx = 0; idx < globalIndexSize; ++idx)
40  {
41    transMap[axisDestGlobalIndex_[idx]].push_back(axisDestGlobalSize_-axisDestGlobalIndex_[idx]-1);
42    transWeight[axisDestGlobalIndex_[idx]].push_back(1.0);
43  }
44
45  int niSrc   = axisSrc_->n.getValue();
46  int sizeSrc = axisSrc_->n_glo.getValue();
47  if (niSrc != sizeSrc) updateAxisValue();
48  else
49  {
50    for (int idx = 0; idx < sizeSrc; ++idx)
51    {
52      axisDest_->value(idx) = axisSrc_->value(sizeSrc-idx-1);
53    }
54  }
55}
56
57/*!
58  Update value on axis after inversing
59  After an axis is inversed, not only the data on it must be inversed but also the value
60*/
61void CAxisAlgorithmInverse::updateAxisValue()
62{
63  CContext* context = CContext::getCurrent();
64  CContextClient* client=context->client;
65
66  int niSrc     = axisSrc_->n.getValue();
67  int ibeginSrc = axisSrc_->begin.getValue();
68
69  CTransformationMapping transformationMap(axisDest_, axisSrc_);
70
71  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0];
72  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[0];
73
74  std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexMapFromDestToSource;
75  std::map<int, std::vector<int> >::const_iterator it = transMap.begin(), ite = transMap.end();
76  for (; it != ite; ++it)
77  {
78    globaIndexMapFromDestToSource[it->first].push_back(make_pair((it->second)[0], (transWeight[it->first])[0]));
79  }
80
81  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource);
82
83  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
84  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
85
86 // Sending global index of original grid source
87  std::map<int,std::vector<size_t> >::const_iterator itbSend = globalIndexToSend.begin(), itSend,
88                                                     iteSend = globalIndexToSend.end();
89 int sendBuffSize = 0;
90 for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size();
91
92 typedef double Scalar;
93 Scalar* sendBuff, *currentSendBuff;
94 if (0 != sendBuffSize) sendBuff = new Scalar[sendBuffSize];
95 for (StdSize idx = 0; idx < sendBuffSize; ++idx) sendBuff[idx] = NumTraits<Scalar>::sfmax();
96
97 int currentBuffPosition = 0;
98 for (itSend = itbSend; itSend != iteSend; ++itSend)
99 {
100   int destRank = itSend->first;
101   const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second;
102   int countSize = globalIndexOfCurrentGridSourceToSend.size();
103   for (int idx = 0; idx < (countSize); ++idx)
104   {
105     int index = globalIndexOfCurrentGridSourceToSend[idx] - ibeginSrc;
106     sendBuff[idx+currentBuffPosition] = (axisSrc_->value)(index);
107   }
108   currentSendBuff = sendBuff + currentBuffPosition;
109   MPI_Send(currentSendBuff, countSize, MPI_DOUBLE, destRank, 14, client->intraComm);
110   currentBuffPosition += countSize;
111 }
112
113 // Receiving global index of grid source sending from current grid source
114 std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv,
115                                                                                     iteRecv = globalIndexToReceive.end();
116 int recvBuffSize = 0;
117 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size();
118
119 Scalar* recvBuff, *currentRecvBuff;
120 if (0 != recvBuffSize) recvBuff = new Scalar [recvBuffSize];
121 for (StdSize idx = 0; idx < recvBuffSize; ++idx) recvBuff[idx] = NumTraits<Scalar>::sfmax();
122
123 int currentRecvBuffPosition = 0;
124 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
125 {
126   MPI_Status status;
127   int srcRank = itRecv->first;
128   int countSize = (itRecv->second).size();
129   currentRecvBuff = recvBuff + currentRecvBuffPosition;
130   MPI_Recv(currentRecvBuff, countSize, MPI_DOUBLE, srcRank, 14, client->intraComm, &status);
131   currentRecvBuffPosition += countSize;
132 }
133
134 int ibeginDest = axisDest_->begin.getValue();
135 currentRecvBuff = recvBuff;
136 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
137 {
138   int countSize = (itRecv->second).size();
139   for (int idx = 0; idx < countSize; ++idx, ++currentRecvBuff)
140   {
141     int ssize = (itRecv->second)[idx].size();
142     for (int i = 0; i < ssize; ++i)
143     {
144       int index = ((itRecv->second)[idx][i]).first - ibeginDest;
145       (axisDest_->value)(index) = *currentRecvBuff;
146     }
147   }
148 }
149
150 if (0 != sendBuffSize) delete [] sendBuff;
151 if (0 != recvBuffSize) delete [] recvBuff;
152}
153
154}
Note: See TracBrowser for help on using the repository browser.