source: XIOS/trunk/src/transformation/domain_algorithm_expand.cpp @ 978

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

Correcting various bugs relating to transformation

+) Fix the order of transformation selection
+) Correct domain transformation selection
+) Reorganize test_remap

Test
+) On Curie
+) All tests pass

File size: 13.5 KB
Line 
1/*!
2   \file domain_algorithm_expand.cpp
3   \author Ha NGUYEN
4   \since 08 Aug 2016
5   \date 19 Sep 2016
6
7   \brief Algorithm for expanding an domain.
8 */
9#include "domain_algorithm_expand.hpp"
10#include "expand_domain.hpp"
11#include "mesh.hpp"
12#include "domain.hpp"
13#include "grid.hpp"
14#include "grid_transformation_factory_impl.hpp"
15#include "context.hpp"
16#include "context_client.hpp"
17
18namespace xios {
19CGenericAlgorithmTransformation* CDomainAlgorithmExpand::create(CGrid* gridDst, CGrid* gridSrc,
20                                                               CTransformation<CDomain>* transformation,
21                                                               int elementPositionInGrid,
22                                                               std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
23                                                               std::map<int, int>& elementPositionInGridSrc2AxisPosition,
24                                                               std::map<int, int>& elementPositionInGridSrc2DomainPosition,
25                                                               std::map<int, int>& elementPositionInGridDst2ScalarPosition,
26                                                               std::map<int, int>& elementPositionInGridDst2AxisPosition,
27                                                               std::map<int, int>& elementPositionInGridDst2DomainPosition)
28{
29  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
30  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
31
32  CExpandDomain* expandDomain = dynamic_cast<CExpandDomain*> (transformation);
33  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
34  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
35
36  return (new CDomainAlgorithmExpand(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], expandDomain));
37}
38
39bool CDomainAlgorithmExpand::registerTrans()
40{
41  CGridTransformationFactory<CDomain>::registerTransformation(TRANS_EXPAND_DOMAIN, create);
42}
43
44CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,
45                                               CDomain* domainSource,
46                                               CExpandDomain* expandDomain)
47: CDomainAlgorithmTransformation(domainDestination, domainSource)
48{
49  if (domainDestination == domainSource)
50  {
51    ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)",
52           << "Domain source and domain destination are the same. Please make sure domain destination refers to domain source" << std::endl
53           << "Domain source " <<domainSource->getId() << std::endl
54           << "Domain destination " <<domainDestination->getId() << std::endl);
55  }
56
57//  if (!domainDestination->hasRefTo(domainSource))
58//  {
59//    ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)",
60//           << "Domain domain destination must refer to domain source (directly or indirectly) by domain_ref " << std::endl
61//           << "Domain source " <<domainSource->getId() << std::endl
62//           << "Domain destination " <<domainDestination->getId() << std::endl);
63//  }
64
65  this->type_ = (ELEMENT_MODIFICATION_WITH_DATA);
66  expandDomain->checkValid(domainDestination);
67
68  switch (expandDomain->type)
69  {
70    case CExpandDomain::type_attr::node :
71      expandDomainNodeConnectivity(domainDestination,
72                                   domainSource);
73      break;
74    case CExpandDomain::type_attr::edge :
75      expandDomainEdgeConnectivity(domainDestination,
76                                   domainSource);
77      break;
78    default:
79      break;
80  }
81}
82
83/*!
84 *  Expand domain with edge-type neighbor
85 *  \param[in/out] domainDestination domain destination and will be modified
86 *  \param[in] domainSource domain source
87 */
88void CDomainAlgorithmExpand::expandDomainEdgeConnectivity(CDomain* domainDestination,
89                                                          CDomain* domainSource)
90{
91  CContext* context = CContext::getCurrent();
92  CContextClient* client=context->client;
93
94  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
95  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
96  CArray<int,2> neighborsSrc;
97
98  int type = 1; // For edge
99  CMesh mesh;
100  mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
101  updateDomainAttributes(domainDestination, domainSource, neighborsSrc);
102}
103
104/*!
105 *  Expand domain with node-type neighbor
106 *  \param[in/out] domainDestination domain destination and will be modified
107 *  \param[in] domainSource domain source
108 */
109void CDomainAlgorithmExpand::expandDomainNodeConnectivity(CDomain* domainDestination,
110                                                          CDomain* domainSource)
111{
112  CContext* context = CContext::getCurrent();
113  CContextClient* client=context->client;
114
115  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
116  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
117  CArray<int,2> neighborsSrc;
118
119  int type = 0; // For node
120  CMesh mesh;
121  mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
122  updateDomainAttributes(domainDestination, domainSource, neighborsSrc);
123}
124
125/*!
126 *  Extend domain destination and update its attributes
127 *  Suppose that domain destination and domain source have the same values for all attributes (by inheritance)
128 *  \param [in/out] domainDestination domain destination
129 *  \param [in] domainSource domain source
130 *  \param [in] neighborsDomainSrc domain extended part
131 */
132void CDomainAlgorithmExpand::updateDomainAttributes(CDomain* domainDestination,
133                                                    CDomain* domainSource,
134                                                    CArray<int,2>& neighborsDomainSrc)
135{
136  CContext* context = CContext::getCurrent();
137  CContextClient* client=context->client;
138
139  // First of all, "copy" all attributes of domain source to domain destination
140  StdString domainDstRef = (!domainDestination->domain_ref.isEmpty()) ? domainDestination->domain_ref.getValue()
141                                                                      : "";
142  if (domainDstRef != domainSource->getId())
143  {
144    domainDestination->domain_ref.setValue(domainSource->getId());
145    domainDestination->solveRefInheritance(true);
146  }
147
148  if (domainDstRef.empty()) domainDestination->domain_ref.reset();
149  else domainDestination->domain_ref.setValue(domainDstRef);
150
151  // Now extend domain destination
152  int niGlob = domainSource->ni_glo;
153  CArray<bool,1>& mask_1d_src = domainSource->mask_1d;
154  CArray<int,1>& i_index_src = domainSource->i_index;
155  CArray<double,1>& lon_src = domainSource->lonvalue_1d;
156  CArray<double,1>& lat_src = domainSource->latvalue_1d;
157  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
158  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
159  CArray<int,1>& data_i_index_src = domainSource->data_i_index;
160
161  int oldNbLocal = i_index_src.numElements(), index, globalIndex;
162  // Uncompress data_i_index
163  CArray<int,1> data_i_index_src_full(oldNbLocal);
164  int nbUnMaskedPointOnLocalDomain = 0;
165  data_i_index_src_full = -1; // Suppose all values are masked
166  for (int idx = 0; idx < data_i_index_src.numElements(); ++idx)
167  {
168    int dataIdx = data_i_index_src(idx);
169    if ((0 <= dataIdx) && (dataIdx < oldNbLocal))
170    {
171      data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIdx;
172      ++nbUnMaskedPointOnLocalDomain;
173    }
174  }
175
176  CArray<bool,1>& mask_1d_dst = domainDestination->mask_1d;
177  CArray<int,1>& i_index_dst = domainDestination->i_index;
178  CArray<int,1>& j_index_dst = domainDestination->j_index;
179  CArray<double,1>& lon_dst = domainDestination->lonvalue_1d;
180  CArray<double,1>& lat_dst = domainDestination->latvalue_1d;
181  CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d;
182  CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d;
183  CArray<int,1>& data_i_index_dst = domainDestination->data_i_index;
184  CArray<int,1>& data_j_index_dst = domainDestination->data_j_index;
185
186  // Resize all array-like attributes of domain destination
187  int nbNeighbor    = neighborsDomainSrc.shape()[1];
188  int newNbLocalDst = nbNeighbor + oldNbLocal;
189  int nVertex       = bounds_lon_dst.shape()[0];
190
191  mask_1d_dst.resizeAndPreserve(newNbLocalDst);
192  i_index_dst.resizeAndPreserve(newNbLocalDst);
193  j_index_dst.resizeAndPreserve(newNbLocalDst);
194  lon_dst.resizeAndPreserve(newNbLocalDst);
195  lat_dst.resizeAndPreserve(newNbLocalDst);
196  bounds_lon_dst.resizeAndPreserve(nVertex, newNbLocalDst);
197  bounds_lat_dst.resizeAndPreserve(nVertex, newNbLocalDst);
198  CArray<int,1> data_i_index_dst_full(newNbLocalDst);
199  data_i_index_dst_full(Range(0,oldNbLocal-1)) = data_i_index_src_full;
200  data_i_index_dst_full(Range(oldNbLocal,newNbLocalDst-1)) = -1;
201
202  // 1. Fill in array relating to global index (i_index, j_index, transmap, etc, ...)
203  // Global index mapping between destination and source
204  this->transformationMapping_.resize(1);
205  this->transformationWeight_.resize(1);
206  TransformationIndexMap& transMap = this->transformationMapping_[0];
207  TransformationWeightMap& transWeight = this->transformationWeight_[0];
208
209  transMap.rehash(std::ceil(newNbLocalDst/transMap.max_load_factor()));
210  transWeight.rehash(std::ceil(newNbLocalDst/transWeight.max_load_factor()));
211  // First, index mapping for local domain
212  for (int idx = 0; idx < oldNbLocal; ++idx)
213  {
214    index = i_index_dst(idx);
215    transMap[index].push_back(index);
216    transWeight[index].push_back(1.0);
217  }
218  // Then, index mapping for extended part
219  for (int idx = 0; idx < nbNeighbor; ++idx)
220  {
221    index = idx + oldNbLocal;
222    globalIndex = neighborsDomainSrc(0,idx);
223    i_index_dst(index) = globalIndex;
224    j_index_dst(index) = 0;
225    transMap[globalIndex].push_back(globalIndex);
226    transWeight[globalIndex].push_back(1.0);
227  }
228
229  // 2. Exchange local info among domains (lon,lat,bounds,mask,etc,...)
230  CClientClientDHTDouble::Index2VectorInfoTypeMap localData;
231  localData.rehash(std::ceil(oldNbLocal/localData.max_load_factor()));
232  // Information exchanged among domains (attention to their order), number in parentheses presents size of data
233  // lon(1) + lat(1) + bounds_lon(nVertex) + bounds_lat(nVertex) + mask(1) + data_i_index(1)
234  int dataPackageSize =  1 + 1 + // lon + lat
235                         nVertex + nVertex + //bounds_lon + bounds_lat
236                         1 + // mask_1d_dst;
237                         1; // data_i_index
238  // Initialize database
239  for (int idx = 0; idx < oldNbLocal; ++idx)
240  {
241    index = i_index_src(idx);
242    localData[index].resize(dataPackageSize);
243    std::vector<double>& data = localData[index];
244
245    //Pack data
246    int dataIdx = 0;
247    data[dataIdx] = lon_src(idx);++dataIdx;
248    data[dataIdx] = lat_src(idx);++dataIdx;
249    for (int i = 0; i < nVertex; ++i)
250    {
251      data[dataIdx] = bounds_lon_src(i,idx); ++dataIdx;
252    }
253    for (int i = 0; i < nVertex; ++i)
254    {
255      data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx;
256    }
257    data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1; ++dataIdx;
258    data[dataIdx] = data_i_index_src_full(idx);
259  }
260
261  CClientClientDHTDouble dhtData(localData,client->intraComm);
262  CArray<size_t,1> neighborInd(nbNeighbor);
263  for (int idx = 0; idx < nbNeighbor; ++idx)
264    neighborInd(idx) = neighborsDomainSrc(0,idx);
265
266  // Compute local data on other domains
267  dhtData.computeIndexInfoMapping(neighborInd);
268  CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap();
269  CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it;
270  // Ok get neighbor data
271  size_t nIdx;
272  int nbUnMaskedPointOnExtendedPart = 0;
273  for (int idx = 0; idx < nbNeighbor; ++idx)
274  {
275    nIdx  = neighborInd(idx);
276    it = neighborData.find(nIdx);
277    if (ite != it)
278    {
279      index = idx + oldNbLocal;
280      std::vector<double>& data = it->second;
281      // Unpack data
282      int dataIdx = 0;
283      lon_dst(index) = data[dataIdx]; ++dataIdx;
284      lat_dst(index) = data[dataIdx]; ++dataIdx;
285      for (int i = 0; i < nVertex; ++i)
286      {
287        bounds_lon_dst(i,index) = data[dataIdx]; ++dataIdx;
288      }
289      for (int i = 0; i < nVertex; ++i)
290      {
291        bounds_lat_dst(i,index) = data[dataIdx]; ++dataIdx;
292      }
293      mask_1d_dst(index) = (1.0 == data[dataIdx]) ? true : false; ++dataIdx;
294      data_i_index_dst_full(index) = (int)(data[dataIdx]);
295      if (0 <= data_i_index_dst_full(index))
296      {
297        data_i_index_dst_full(index) = index;
298        ++nbUnMaskedPointOnExtendedPart;
299      }
300    }
301  }
302
303
304  // Finally, update data_i_index
305  int nbUnMaskedPointOnNewDstDomain = (nbUnMaskedPointOnExtendedPart + nbUnMaskedPointOnLocalDomain);
306  int count = 0, dataIdx;
307  for (int idx = 0; idx < newNbLocalDst; ++idx)
308  {
309    dataIdx = data_i_index_dst_full(idx);
310    if ((0 <= dataIdx) && (dataIdx < newNbLocalDst))
311    {
312      ++count;
313    }
314  }
315
316  data_i_index_dst.resize(count);
317  data_j_index_dst.resize(count);
318  data_j_index_dst = 0;
319
320  count = 0;
321  for (int idx = 0; idx < newNbLocalDst; ++idx)
322  {
323    dataIdx = data_i_index_dst_full(idx);
324    if ((0 <= dataIdx) && (dataIdx < newNbLocalDst))
325    {
326      data_i_index_dst(count) = dataIdx;
327      ++count;
328    }
329  }
330
331  // Update ni
332  domainDestination->ni.setValue(newNbLocalDst);
333
334}
335
336
337/*!
338  Compute the index mapping between domain on grid source and one on grid destination
339*/
340void CDomainAlgorithmExpand::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
341{
342
343}
344
345}
Note: See TracBrowser for help on using the repository browser.