source: XIOS/trunk/src/filter/spatial_transform_filter.cpp @ 790

Last change on this file since 790 was 790, checked in by rlacroix, 8 years ago

Add a new field attribute grid_path to allow chaining spatial transformations.

Intermediate grids must be separated by comma and the source/destination grids must not be listed.

File size: 6.7 KB
Line 
1#include "spatial_transform_filter.hpp"
2#include "grid_transformation.hpp"
3#include "context.hpp"
4#include "context_client.hpp"
5
6namespace xios
7{
8  CSpatialTransformFilter::CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine)
9    : CFilter(gc, 1, engine)
10  { /* Nothing to do */ }
11
12  std::pair<boost::shared_ptr<CSpatialTransformFilter>, boost::shared_ptr<CSpatialTransformFilter> >
13  CSpatialTransformFilter::buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid)
14  {
15    if (!srcGrid || !destGrid)
16      ERROR("std::pair<boost::shared_ptr<CSpatialTransformFilter>, boost::shared_ptr<CSpatialTransformFilter> >"
17            "buildFilterGraph(CGarbageCollector& gc, CGrid* srcGrid, CGrid* destGrid)",
18            "Impossible to build the filter graph if either the source or the destination grid are null.");
19
20    boost::shared_ptr<CSpatialTransformFilter> firstFilter, lastFilter;
21    // Note that this loop goes from the last transformation to the first transformation
22    do
23    {
24      CGridTransformation* gridTransformation = destGrid->getTransformations();
25      CSpatialTransformFilterEngine* engine = CSpatialTransformFilterEngine::get(destGrid->getTransformations());
26      boost::shared_ptr<CSpatialTransformFilter> filter(new CSpatialTransformFilter(gc, engine));
27
28      if (!lastFilter)
29        lastFilter = filter;
30      else
31        filter->connectOutput(firstFilter, 0);
32
33      firstFilter = filter;
34      destGrid = gridTransformation->getGridSource();
35    }
36    while (destGrid != srcGrid);
37
38    return std::make_pair(firstFilter, lastFilter);
39  }
40
41  CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)
42    : gridTransformation(gridTransformation)
43  {
44    if (!gridTransformation)
45      ERROR("CSpatialTransformFilterEngine::CSpatialTransformFilterEngine(CGridTransformation* gridTransformation)",
46            "Impossible to construct a spatial transform filter engine without a valid grid transformation.");
47  }
48
49  std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> > CSpatialTransformFilterEngine::engines;
50
51  CSpatialTransformFilterEngine* CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)
52  {
53    if (!gridTransformation)
54      ERROR("CSpatialTransformFilterEngine& CSpatialTransformFilterEngine::get(CGridTransformation* gridTransformation)",
55            "Impossible to get the requested engine, the grid transformation is invalid.");
56
57    std::map<CGridTransformation*, boost::shared_ptr<CSpatialTransformFilterEngine> >::iterator it = engines.find(gridTransformation);
58    if (it == engines.end())
59    {
60      boost::shared_ptr<CSpatialTransformFilterEngine> engine(new CSpatialTransformFilterEngine(gridTransformation));
61      it = engines.insert(std::make_pair(gridTransformation, engine)).first;
62    }
63
64    return it->second.get();
65  }
66
67  CDataPacketPtr CSpatialTransformFilterEngine::apply(std::vector<CDataPacketPtr> data)
68  {
69    CDataPacketPtr packet(new CDataPacket);
70    packet->date = data[0]->date;
71    packet->timestamp = data[0]->timestamp;
72    packet->status = data[0]->status;
73
74    if (packet->status == CDataPacket::NO_ERROR)
75    {
76      packet->data.resize(gridTransformation->getGridDestination()->storeIndex_client.numElements());
77      apply(data[0]->data, packet->data);
78    }
79
80    return packet;
81  }
82
83  void CSpatialTransformFilterEngine::apply(const CArray<double, 1>& dataSrc, CArray<double,1>& dataDest)
84  {
85    CContextClient* client = CContext::getCurrent()->client;
86
87    const std::map<int, CArray<int,1> >& localIndexToSend = gridTransformation->getLocalIndexToSendFromGridSource();
88    const std::map<int, std::vector<std::vector<std::pair<int,double> > > >& localIndexToReceive = gridTransformation->getLocalIndexToReceiveOnGridDest();
89
90    dataDest = 0.0;
91
92    // Sending data from field sources to do transformations
93    std::map<int, CArray<int,1> >::const_iterator itbSend = localIndexToSend.begin(), itSend,
94                                                  iteSend = localIndexToSend.end();
95    int idxSendBuff = 0;
96    std::vector<double*> sendBuff(localIndexToSend.size());
97    for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
98    {
99      if (0 != itSend->second.numElements())
100        sendBuff[idxSendBuff] = new double[itSend->second.numElements()];
101    }
102
103    idxSendBuff = 0;
104    std::vector<MPI_Request> sendRequest;
105    for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
106    {
107      int destRank = itSend->first;
108      const CArray<int,1>& localIndex_p = itSend->second;
109      int countSize = localIndex_p.numElements();
110      for (int idx = 0; idx < countSize; ++idx)
111      {
112        sendBuff[idxSendBuff][idx] = dataSrc(localIndex_p(idx));
113      }
114      sendRequest.push_back(MPI_Request());
115      MPI_Isend(sendBuff[idxSendBuff], countSize, MPI_DOUBLE, destRank, 12, client->intraComm, &sendRequest.back());
116    }
117
118    // Receiving data on destination fields
119    std::map<int,std::vector<std::vector<std::pair<int,double> > > >::const_iterator itbRecv = localIndexToReceive.begin(), itRecv,
120                                                                                     iteRecv = localIndexToReceive.end();
121    int recvBuffSize = 0;
122    for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize = (recvBuffSize < itRecv->second.size())
123                                                                     ? itRecv->second.size() : recvBuffSize;
124    double* recvBuff;
125    if (0 != recvBuffSize) recvBuff = new double[recvBuffSize];
126    for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
127    {
128      MPI_Status status;
129      int srcRank = itRecv->first;
130      int countSize = itRecv->second.size();
131      MPI_Recv(recvBuff, recvBuffSize, MPI_DOUBLE, srcRank, 12, client->intraComm, &status);
132      int countBuff = 0;
133      MPI_Get_count(&status, MPI_DOUBLE, &countBuff);
134      if (countBuff != countSize)
135        ERROR("CSpatialTransformFilterEngine::apply(const CArray<double, 1>& dataSrc, CArray<double,1>& dataDest)",
136              "Incoherent between the received size and expected size");
137      for (int idx = 0; idx < countSize; ++idx)
138      {
139        const std::vector<std::pair<int,double> >& localIndex_p = itRecv->second[idx];
140        int numIndex = localIndex_p.size();
141        for (int i = 0; i < numIndex; ++i)
142        {
143          dataDest(localIndex_p[i].first) += recvBuff[idx] * localIndex_p[i].second;
144        }
145      }
146    }
147
148
149    if (!sendRequest.empty()) MPI_Waitall(sendRequest.size(), &sendRequest[0], MPI_STATUSES_IGNORE);
150    idxSendBuff = 0;
151    for (itSend = itbSend; itSend != iteSend; ++itSend, ++idxSendBuff)
152    {
153      if (0 != itSend->second.numElements())
154        delete [] sendBuff[idxSendBuff];
155    }
156    if (0 != recvBuffSize) delete [] recvBuff;
157  }
158} // namespace xios
Note: See TracBrowser for help on using the repository browser.