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

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

dev_omp

File size: 11.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"
6using namespace ep_lib;
7
8namespace xios
9{
10  CSpatialTransformFilter::CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, double outputValue, size_t inputSlotsCount)
11    : CFilter(gc, inputSlotsCount, engine), outputDefaultValue(outputValue)
12  { /* Nothing to do */ }
13
14  std::pair<boost::shared_ptr<CSpatialTransformFilter>, boost::shared_ptr<CSpatialTransformFilter> >
15  CSpatialTransformFilter::buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid, bool hasMissingValue, double missingValue)
16  {
17    if (!srcGrid || !destGrid)
18      ERROR("std::pair<boost::shared_ptr<CSpatialTransformFilter>, boost::shared_ptr<CSpatialTransformFilter> >"
19            "buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid)",
20            "Impossible to build the filter graph if either the source or the destination grid are null.");
21
22    boost::shared_ptr<CSpatialTransformFilter> firstFilter, lastFilter;
23    // Note that this loop goes from the last transformation to the first transformation
24    do
25    {
26      CGridTransformation* gridTransformation = destGrid->getTransformations();
27      CSpatialTransformFilterEngine* engine = CSpatialTransformFilterEngine::get(destGrid->getTransformations());
28      const std::vector<StdString>& auxInputs = gridTransformation->getAuxInputs();
29      size_t inputCount = 1 + (auxInputs.empty() ? 0 : auxInputs.size());
30      double defaultValue  = (hasMissingValue) ? std::numeric_limits<double>::quiet_NaN() : 0.0;
31      boost::shared_ptr<CSpatialTransformFilter> filter(new CSpatialTransformFilter(gc, engine, defaultValue, inputCount));
32
33      if (!lastFilter)
34        lastFilter = filter;
35      else
36        filter->connectOutput(firstFilter, 0);
37
38      firstFilter = filter;
39      for (size_t idx = 0; idx < auxInputs.size(); ++idx)
40      {
41        CField* fieldAuxInput = CField::get(auxInputs[idx]);
42        fieldAuxInput->buildFilterGraph(gc, false);
43        fieldAuxInput->getInstantDataFilter()->connectOutput(firstFilter,idx+1);
44      }
45
46      destGrid = gridTransformation->getGridSource();
47    }
48    while (destGrid != srcGrid);
49
50    return std::make_pair(firstFilter, lastFilter);
51  }
52
53  void CSpatialTransformFilter::onInputReady(std::vector<CDataPacketPtr> data)
54  {
55    CSpatialTransformFilterEngine* spaceFilter = static_cast<CSpatialTransformFilterEngine*>(engine);
56    CDataPacketPtr outputPacket = spaceFilter->applyFilter(data, outputDefaultValue);
57    if (outputPacket)
58      onOutputReady(outputPacket);
59  }
60
61  CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)
62    : gridTransformation(gridTransformation)
63  {
64    if (!gridTransformation)
65      ERROR("CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)",
66            "Impossible to construct a spatial transform filter engine without a valid grid transformation.");
67  }
68
69  //std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> > CSpatialTransformFilterEngine::engines;
70  std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> > *CSpatialTransformFilterEngine::engines_ptr = 0;
71
72  CSpatialTransformFilterEngine* CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)
73  {
74    if (!gridTransformation)
75      ERROR("CSpatialTransformFilterEngine& CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)",
76            "Impossible to get the requested engine, the grid transformation is invalid.");
77   
78    if(engines_ptr == NULL) engines_ptr = new std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> >;
79
80    //std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> >::iterator it = engines.find(gridTransformation);
81    std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> >::iterator it = engines_ptr->find(gridTransformation);
82    //if (it == engines.end())
83    if (it == engines_ptr->end())
84    {
85      boost::shared_ptr<CSpatialTransformFilterEngine> engine(new CSpatialTransformFilterEngine(gridTransformation));
86      //it = engines.insert(std::make_pair(gridTransformation, engine)).first;
87      it = engines_ptr->insert(std::make_pair(gridTransformation, engine)).first;
88    }
89
90    return it->second.get();
91  }
92
93  CDataPacketPtr CSpatialTransformFilterEngine::apply(std::vector<CDataPacketPtr> data)
94  {
95    /* Nothing to do */
96  }
97
98  CDataPacketPtr CSpatialTransformFilterEngine::applyFilter(std::vector<CDataPacketPtr> data, double defaultValue)
99  {
100    CDataPacketPtr packet(new CDataPacket);
101    packet->date = data[0]->date;
102    packet->timestamp = data[0]->timestamp;
103    packet->status = data[0]->status;
104
105    if (packet->status == CDataPacket::NO_ERROR)
106    {
107      if (1 < data.size())  // Dynamical transformations
108      {
109        std::vector<CArray<double,1>* > dataAuxInputs(data.size()-1);
110        for (size_t idx = 0; idx < dataAuxInputs.size(); ++idx) dataAuxInputs[idx] = &(data[idx+1]->data);
111        gridTransformation->computeAll(dataAuxInputs, packet->timestamp);
112      }
113      packet->data.resize(gridTransformation->getGridDestination()->storeIndex_client.numElements());
114      if (0 != packet->data.numElements())
115        (packet->data)(0) = defaultValue;
116      apply(data[0]->data, packet->data);
117    }
118
119    return packet;
120  }
121
122  void CSpatialTransformFilterEngine::apply(const CArray<double, 1>& dataSrc, CArray<double,1>& dataDest)
123  {
124    CContextClient* client = CContext::getCurrent()->client;
125
126    // Get default value for output data
127    bool ignoreMissingValue = false; 
128    double defaultValue = std::numeric_limits<double>::quiet_NaN();
129    if (0 != dataDest.numElements()) ignoreMissingValue = NumTraits<double>::isnan(dataDest(0));
130
131    const std::list<CGridTransformation::SendingIndexGridSourceMap>& listLocalIndexSend = gridTransformation->getLocalIndexToSendFromGridSource();
132    const std::list<CGridTransformation::RecvIndexGridDestinationMap>& listLocalIndexToReceive = gridTransformation->getLocalIndexToReceiveOnGridDest();
133    const std::list<size_t>& listNbLocalIndexToReceive = gridTransformation->getNbLocalIndexToReceiveOnGridDest();
134    const std::list<std::vector<bool> >& listLocalIndexMaskOnDest = gridTransformation->getLocalMaskIndexOnGridDest();
135    const std::vector<CGenericAlgorithmTransformation*>& listAlgos = gridTransformation->getAlgos();
136
137    CArray<double,1> dataCurrentDest(dataSrc.copy());
138
139    std::list<CGridTransformation::SendingIndexGridSourceMap>::const_iterator itListSend  = listLocalIndexSend.begin(),
140                                                                              iteListSend = listLocalIndexSend.end();
141    std::list<CGridTransformation::RecvIndexGridDestinationMap>::const_iterator itListRecv = listLocalIndexToReceive.begin();
142    std::list<size_t>::const_iterator itNbListRecv = listNbLocalIndexToReceive.begin();
143    std::list<std::vector<bool> >::const_iterator itLocalMaskIndexOnDest = listLocalIndexMaskOnDest.begin();
144    std::vector<CGenericAlgorithmTransformation*>::const_iterator itAlgo = listAlgos.begin();
145
146    for (; itListSend != iteListSend; ++itListSend, ++itListRecv, ++itNbListRecv, ++itLocalMaskIndexOnDest, ++itAlgo)
147    {
148      CArray<double,1> dataCurrentSrc(dataCurrentDest);
149      const CGridTransformation::SendingIndexGridSourceMap& localIndexToSend = *itListSend;
150
151      // Sending data from field sources to do transformations
152      std::map<int, CArray<int,1> >::const_iterator itbSend = localIndexToSend.begin(), itSend,
153                                                    iteSend = localIndexToSend.end();
154      int idxSendBuff = 0;
155      std::vector<double*> sendBuff(localIndexToSend.size());
156      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
157      {
158        if (0 != itSend->second.numElements())
159          sendBuff[idxSendBuff] = new double[itSend->second.numElements()];
160      }
161
162      idxSendBuff = 0;
163      std::vector<MPI_Request> sendRecvRequest(localIndexToSend.size() + itListRecv->size());
164      int position = 0;
165      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
166      {
167        int destRank = itSend->first;
168        const CArray<int,1>& localIndex_p = itSend->second;
169        int countSize = localIndex_p.numElements();
170        for (int idx = 0; idx < countSize; ++idx)
171        {
172          sendBuff[idxSendBuff][idx] = dataCurrentSrc(localIndex_p(idx));
173        }
174        MPI_Isend(sendBuff[idxSendBuff], countSize, MPI_DOUBLE, destRank, 12, client->intraComm, &sendRecvRequest[position++]);
175      }
176
177      // Receiving data on destination fields
178      const CGridTransformation::RecvIndexGridDestinationMap& localIndexToReceive = *itListRecv;
179      CGridTransformation::RecvIndexGridDestinationMap::const_iterator itbRecv = localIndexToReceive.begin(), itRecv,
180                                                                       iteRecv = localIndexToReceive.end();
181      int recvBuffSize = 0;
182      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += itRecv->second.size(); //(recvBuffSize < itRecv->second.size())
183                                                                       //? itRecv->second.size() : recvBuffSize;
184      double* recvBuff;
185      if (0 != recvBuffSize) recvBuff = new double[recvBuffSize];
186      int currentBuff = 0;
187      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
188      {
189        int srcRank = itRecv->first;
190        int countSize = itRecv->second.size();
191        MPI_Irecv(recvBuff + currentBuff, countSize, MPI_DOUBLE, srcRank, 12, client->intraComm, &sendRecvRequest[position++]);
192        currentBuff += countSize;
193      }
194      std::vector<MPI_Status> status(sendRecvRequest.size());
195      MPI_Waitall(sendRecvRequest.size(), &sendRecvRequest[0], &status[0]);
196
197      dataCurrentDest.resize(*itNbListRecv);
198      const std::vector<bool>& localMaskDest = *itLocalMaskIndexOnDest;
199      for (int i = 0; i < localMaskDest.size(); ++i)
200        if (localMaskDest[i]) dataCurrentDest(i) = 0.0;
201        else dataCurrentDest(i) = defaultValue;
202
203      std::vector<bool> localInitFlag(dataCurrentDest.numElements(), true);
204      currentBuff = 0;
205      bool firstPass=true; 
206      for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
207      {
208        int countSize = itRecv->second.size();
209        const std::vector<std::pair<int,double> >& localIndex_p = itRecv->second;
210        (*itAlgo)->apply(localIndex_p,
211                         recvBuff+currentBuff,
212                         dataCurrentDest,
213                         localInitFlag,
214                         ignoreMissingValue,firstPass);
215
216        currentBuff += countSize;
217        firstPass=false ;
218      }
219
220      (*itAlgo)->updateData(dataCurrentDest);
221
222      idxSendBuff = 0;
223      for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
224      {
225        if (0 != itSend->second.numElements())
226          delete [] sendBuff[idxSendBuff];
227      }
228      if (0 != recvBuffSize) delete [] recvBuff;
229    }
230    if (dataCurrentDest.numElements() != dataDest.numElements())
231    ERROR("CSpatialTransformFilterEngine::apply(const CArray<double, 1>& dataSrc, CArray<double,1>& dataDest)",
232          "Incoherent between the received size and expected size. " << std::endl
233          << "Expected size: " << dataDest.numElements() << std::endl
234          << "Received size: " << dataCurrentDest.numElements());
235
236    dataDest = dataCurrentDest;
237  }
238} // namespace xios
Note: See TracBrowser for help on using the repository browser.