source: XIOS/dev/branch_openmp/src/filter/spatial_transform_filter.cpp @ 1482

Last change on this file since 1482 was 1482, checked in by yushan, 3 years ago

Branch EP merged with Dev_cmip6 @r1481

File size: 14.4 KB
Line 
1#include "mpi.hpp"
2#include "spatial_transform_filter.hpp"
3#include "grid_transformation.hpp"
4#include "context.hpp"
5#include "context_client.hpp"
6#include "timer.hpp"
7using namespace ep_lib;
8
9namespace xios
10{
11  CSpatialTransformFilter::CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, double outputValue, size_t inputSlotsCount)
12    : CFilter(gc, inputSlotsCount, engine), outputDefaultValue(outputValue)
13  { /* Nothing to do */ }
14
15  std::pair<boost::shared_ptr<CSpatialTransformFilter>, boost::shared_ptr<CSpatialTransformFilter> >
16  CSpatialTransformFilter::buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid, bool hasMissingValue, double missingValue)
17  {
18    if (!srcGrid || !destGrid)
19      ERROR("std::pair<boost::shared_ptr<CSpatialTransformFilter>, boost::shared_ptr<CSpatialTransformFilter> >"
20            "buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid)",
21            "Impossible to build the filter graph if either the source or the destination grid are null.");
22
23    boost::shared_ptr<CSpatialTransformFilter> firstFilter, lastFilter;
24    // Note that this loop goes from the last transformation to the first transformation
25    do
26    {
27      CGridTransformation* gridTransformation = destGrid->getTransformations();
28      CSpatialTransformFilterEngine* engine = CSpatialTransformFilterEngine::get(destGrid->getTransformations());
29      const std::vector<StdString>& auxInputs = gridTransformation->getAuxInputs();
30      size_t inputCount = 1 + (auxInputs.empty() ? 0 : auxInputs.size());
31      double defaultValue  = (hasMissingValue) ? std::numeric_limits<double>::quiet_NaN() : 0.0;
32
33
34      const CGridTransformationSelector::ListAlgoType& algoList = gridTransformation->getAlgoList() ;
35      CGridTransformationSelector::ListAlgoType::const_iterator it  ;
36
37      bool isSpatialTemporal=false ;
38      for (it=algoList.begin();it!=algoList.end();++it)  if (it->second.first == TRANS_TEMPORAL_SPLITTING) isSpatialTemporal=true ;
39
40      boost::shared_ptr<CSpatialTransformFilter> filter ;
41      if( isSpatialTemporal) filter = boost::shared_ptr<CSpatialTransformFilter>(new CSpatialTemporalFilter(gc, engine, gridTransformation, defaultValue, inputCount));
42      else filter = boost::shared_ptr<CSpatialTransformFilter>(new CSpatialTransformFilter(gc, engine, defaultValue, inputCount));
43
44     
45      if (!lastFilter)
46        lastFilter = filter;
47      else
48        filter->connectOutput(firstFilter, 0);
49
50      firstFilter = filter;
51      for (size_t idx = 0; idx < auxInputs.size(); ++idx)
52      {
53        CField* fieldAuxInput = CField::get(auxInputs[idx]);
54        fieldAuxInput->buildFilterGraph(gc, false);
55        fieldAuxInput->getInstantDataFilter()->connectOutput(firstFilter,idx+1);
56      }
57
58      destGrid = gridTransformation->getGridSource();
59    }
60    while (destGrid != srcGrid);
61
62    return std::make_pair(firstFilter, lastFilter);
63  }
64
65  void CSpatialTransformFilter::onInputReady(std::vector<CDataPacketPtr> data)
66  {
67    CSpatialTransformFilterEngine* spaceFilter = static_cast<CSpatialTransformFilterEngine*>(engine);
68    CDataPacketPtr outputPacket = spaceFilter->applyFilter(data, outputDefaultValue);
69    if (outputPacket)
70      onOutputReady(outputPacket);
71  }
72
73
74
75
76
77  CSpatialTemporalFilter::CSpatialTemporalFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, CGridTransformation* gridTransformation, double outputValue, size_t inputSlotsCount)
78    : CSpatialTransformFilter(gc, engine, outputValue, inputSlotsCount), record(0)
79  {
80      const CGridTransformationSelector::ListAlgoType& algoList = gridTransformation->getAlgoList() ;
81      CGridTransformationSelector::ListAlgoType::const_iterator it  ;
82
83      int pos ;
84      for (it=algoList.begin();it!=algoList.end();++it) 
85        if (it->second.first == TRANS_TEMPORAL_SPLITTING)
86        {
87          pos=it->first ;
88          if (pos < algoList.size()-1)
89            ERROR("SpatialTemporalFilter::CSpatialTemporalFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, CGridTransformation* gridTransformation, double outputValue, size_t inputSlotsCount))",
90                  "temporal splitting operation must be the last of whole transformation on same grid") ;
91        }
92         
93      CGrid* grid=gridTransformation->getGridDestination() ;
94
95      CAxis* axis = grid->getAxis(gridTransformation->getElementPositionInGridDst2AxisPosition().find(pos)->second) ;
96
97      nrecords = axis->index.numElements() ;
98  }
99
100
101  void CSpatialTemporalFilter::onInputReady(std::vector<CDataPacketPtr> data)
102  {
103    CSpatialTransformFilterEngine* spaceFilter = static_cast<CSpatialTransformFilterEngine*>(engine);
104    CDataPacketPtr outputPacket = spaceFilter->applyFilter(data, outputDefaultValue);
105
106    if (outputPacket)
107    {
108      size_t nelements=outputPacket->data.numElements() ;
109      if (!tmpData.numElements())
110      {
111        tmpData.resize(nelements);
112        tmpData=outputDefaultValue ;
113      }
114
115      nelements/=nrecords ;
116      size_t offset=nelements*record ;
117      for(size_t i=0;i<nelements;++i)  tmpData(i+offset) = outputPacket->data(i) ;
118   
119      record ++ ;
120      if (record==nrecords)
121      {
122        record=0 ;
123        CDataPacketPtr packet = CDataPacketPtr(new CDataPacket);
124        packet->date = data[0]->date;
125        packet->timestamp = data[0]->timestamp;
126        packet->status = data[0]->status;
127        packet->data.resize(tmpData.numElements());
128        packet->data = tmpData;
129        onOutputReady(packet);
130        tmpData.resize(0) ;
131      }
132    }
133  }
134
135
136  CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)
137    : gridTransformation(gridTransformation)
138  {
139    if (!gridTransformation)
140      ERROR("CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)",
141            "Impossible to construct a spatial transform filter engine without a valid grid transformation.");
142  }
143
144  //std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> > CSpatialTransformFilterEngine::engines;
145  std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> > *CSpatialTransformFilterEngine::engines_ptr = 0;
146
147  CSpatialTransformFilterEngine* CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)
148  {
149    if (!gridTransformation)
150      ERROR("CSpatialTransformFilterEngine& CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)",
151            "Impossible to get the requested engine, the grid transformation is invalid.");
152   
153    if(engines_ptr == NULL) engines_ptr = new std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> >;
154
155    //std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> >::iterator it = engines.find(gridTransformation);
156    std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> >::iterator it = engines_ptr->find(gridTransformation);
157    //if (it == engines.end())
158    if (it == engines_ptr->end())
159    {
160      boost::shared_ptr<CSpatialTransformFilterEngine> engine(new CSpatialTransformFilterEngine(gridTransformation));
161      //it = engines.insert(std::make_pair(gridTransformation, engine)).first;
162      it = engines_ptr->insert(std::make_pair(gridTransformation, engine)).first;
163    }
164
165    return it->second.get();
166  }
167
168  CDataPacketPtr CSpatialTransformFilterEngine::apply(std::vector<CDataPacketPtr> data)
169  {
170    /* Nothing to do */
171  }
172
173  CDataPacketPtr CSpatialTransformFilterEngine::applyFilter(std::vector<CDataPacketPtr> data, double defaultValue)
174  {
175    CDataPacketPtr packet(new CDataPacket);
176    packet->date = data[0]->date;
177    packet->timestamp = data[0]->timestamp;
178    packet->status = data[0]->status;
179
180    if (packet->status == CDataPacket::NO_ERROR)
181    {
182      if (1 < data.size())  // Dynamical transformations
183      {
184        std::vector<CArray<double,1>* > dataAuxInputs(data.size()-1);
185        for (size_t idx = 0; idx < dataAuxInputs.size(); ++idx) dataAuxInputs[idx] = &(data[idx+1]->data);
186        gridTransformation->computeAll(dataAuxInputs, packet->timestamp);
187      }
188      packet->data.resize(gridTransformation->getGridDestination()->storeIndex_client.numElements());
189      if (0 != packet->data.numElements())
190        (packet->data)(0) = defaultValue;
191      apply(data[0]->data, packet->data);
192    }
193
194    return packet;
195  }
196
197  void CSpatialTransformFilterEngine::apply(const CArray<double, 1>& dataSrc, CArray<double,1>& dataDest)
198  {
199    CTimer::get("CSpatialTransformFilterEngine::apply").resume(); 
200   
201    CContextClient* client = CContext::getCurrent()->client;
202
203    // Get default value for output data
204    bool ignoreMissingValue = false; 
205    double defaultValue = std::numeric_limits<double>::quiet_NaN();
206    if (0 != dataDest.numElements()) ignoreMissingValue = NumTraits<double>::isNan(dataDest(0));
207
208    const std::list<CGridTransformation::SendingIndexGridSourceMap>& listLocalIndexSend = gridTransformation->getLocalIndexToSendFromGridSource();
209    const std::list<CGridTransformation::RecvIndexGridDestinationMap>& listLocalIndexToReceive = gridTransformation->getLocalIndexToReceiveOnGridDest();
210    const std::list<size_t>& listNbLocalIndexToReceive = gridTransformation->getNbLocalIndexToReceiveOnGridDest();
211    const std::list<std::vector<bool> >& listLocalIndexMaskOnDest = gridTransformation->getLocalMaskIndexOnGridDest();
212    const std::vector<CGenericAlgorithmTransformation*>& listAlgos = gridTransformation->getAlgos();
213
214    CArray<double,1> dataCurrentDest(dataSrc.copy());
215
216    std::list<CGridTransformation::SendingIndexGridSourceMap>::const_iterator itListSend  = listLocalIndexSend.begin(),
217                                                                              iteListSend = listLocalIndexSend.end();
218    std::list<CGridTransformation::RecvIndexGridDestinationMap>::const_iterator itListRecv = listLocalIndexToReceive.begin();
219    std::list<size_t>::const_iterator itNbListRecv = listNbLocalIndexToReceive.begin();
220    std::list<std::vector<bool> >::const_iterator itLocalMaskIndexOnDest = listLocalIndexMaskOnDest.begin();
221    std::vector<CGenericAlgorithmTransformation*>::const_iterator itAlgo = listAlgos.begin();
222
223    for (; itListSend != iteListSend; ++itListSend, ++itListRecv, ++itNbListRecv, ++itLocalMaskIndexOnDest, ++itAlgo)
224    {
225      CArray<double,1> dataCurrentSrc(dataCurrentDest);
226      const CGridTransformation::SendingIndexGridSourceMap& localIndexToSend = *itListSend;
227
228      // Sending data from field sources to do transformations
229      std::map<int, CArray<int,1> >::const_iterator itbSend = localIndexToSend.begin(), itSend,
230                                                    iteSend = localIndexToSend.end();
231      int idxSendBuff = 0;
232      std::vector<double*> sendBuff(localIndexToSend.size());
233      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
234      {
235        if (0 != itSend->second.numElements())
236          sendBuff[idxSendBuff] = new double[itSend->second.numElements()];
237      }
238
239      idxSendBuff = 0;
240      std::vector<MPI_Request> sendRecvRequest(localIndexToSend.size() + itListRecv->size());
241      int position = 0;
242      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
243      {
244        int destRank = itSend->first;
245        const CArray<int,1>& localIndex_p = itSend->second;
246        int countSize = localIndex_p.numElements();
247        for (int idx = 0; idx < countSize; ++idx)
248        {
249          sendBuff[idxSendBuff][idx] = dataCurrentSrc(localIndex_p(idx));
250        }
251        MPI_Isend(sendBuff[idxSendBuff], countSize, MPI_DOUBLE, destRank, 12, client->intraComm, &sendRecvRequest[position++]);
252      }
253
254      // Receiving data on destination fields
255      const CGridTransformation::RecvIndexGridDestinationMap& localIndexToReceive = *itListRecv;
256      CGridTransformation::RecvIndexGridDestinationMap::const_iterator itbRecv = localIndexToReceive.begin(), itRecv,
257                                                                       iteRecv = localIndexToReceive.end();
258      int recvBuffSize = 0;
259      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += itRecv->second.size(); //(recvBuffSize < itRecv->second.size())
260                                                                       //? itRecv->second.size() : recvBuffSize;
261      double* recvBuff;
262      if (0 != recvBuffSize) recvBuff = new double[recvBuffSize];
263      int currentBuff = 0;
264      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
265      {
266        int srcRank = itRecv->first;
267        int countSize = itRecv->second.size();
268        MPI_Irecv(recvBuff + currentBuff, countSize, MPI_DOUBLE, srcRank, 12, client->intraComm, &sendRecvRequest[position++]);
269        currentBuff += countSize;
270      }
271      std::vector<MPI_Status> status(sendRecvRequest.size());
272      MPI_Waitall(sendRecvRequest.size(), &sendRecvRequest[0], &status[0]);
273
274      dataCurrentDest.resize(*itNbListRecv);
275      const std::vector<bool>& localMaskDest = *itLocalMaskIndexOnDest;
276      for (int i = 0; i < localMaskDest.size(); ++i)
277        if (localMaskDest[i]) dataCurrentDest(i) = 0.0;
278        else dataCurrentDest(i) = defaultValue;
279
280      std::vector<bool> localInitFlag(dataCurrentDest.numElements(), true);
281      currentBuff = 0;
282      bool firstPass=true; 
283      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
284      {
285        int countSize = itRecv->second.size();
286        const std::vector<std::pair<int,double> >& localIndex_p = itRecv->second;
287        (*itAlgo)->apply(localIndex_p,
288                         recvBuff+currentBuff,
289                         dataCurrentDest,
290                         localInitFlag,
291                         ignoreMissingValue,firstPass);
292
293        currentBuff += countSize;
294        firstPass=false ;
295      }
296
297      (*itAlgo)->updateData(dataCurrentDest);
298
299      idxSendBuff = 0;
300      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
301      {
302        if (0 != itSend->second.numElements())
303          delete [] sendBuff[idxSendBuff];
304      }
305      if (0 != recvBuffSize) delete [] recvBuff;
306    }
307    if (dataCurrentDest.numElements() != dataDest.numElements())
308    ERROR("CSpatialTransformFilterEngine::apply(const CArray<double, 1>& dataSrc, CArray<double,1>& dataDest)",
309          "Incoherent between the received size and expected size. " << std::endl
310          << "Expected size: " << dataDest.numElements() << std::endl
311          << "Received size: " << dataCurrentDest.numElements());
312
313    dataDest = dataCurrentDest;
314
315    CTimer::get("CSpatialTransformFilterEngine::apply").suspend() ;
316  }
317} // namespace xios
Note: See TracBrowser for help on using the repository browser.