Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
TempestOnlineMapIO.cpp
Go to the documentation of this file.
1 /*
2  * =====================================================================================
3  *
4  * Filename: TempestOnlineMapIO.cpp
5  *
6  * Description: All I/O implementations related to TempestOnlineMap
7  *
8  * Version: 1.0
9  * Created: 02/06/2021 02:35:41
10  *
11  * Author: Vijay S. Mahadevan (vijaysm), [email protected]
12  * Company: Argonne National Lab
13  *
14  * =====================================================================================
15  */
16 
17 #include "FiniteElementTools.h"
18 #include "moab/Remapping/TempestOnlineMap.hpp"
19 #include "moab/TupleList.hpp"
20 
21 #ifdef MOAB_HAVE_NETCDFPAR
22 #include "netcdfcpp_par.hpp"
23 #else
24 #include "netcdfcpp.h"
25 #endif
26 
27 #ifdef MOAB_HAVE_PNETCDF
28 #include <pnetcdf.h>
29 
30 #define ERR_PARNC( err ) \
31  if( err != NC_NOERR ) \
32  { \
33  fprintf( stderr, "Error at line %d: %s\n", __LINE__, ncmpi_strerror( err ) ); \
34  MPI_Abort( MPI_COMM_WORLD, 1 ); \
35  }
36 
37 #endif
38 
39 #ifdef MOAB_HAVE_MPI
40 
41 #define MPI_CHK_ERR( err ) \
42  if( err ) \
43  { \
44  std::cout << "MPI Failure. ErrorCode (" << ( err ) << ") "; \
45  std::cout << "\nMPI Aborting... \n"; \
46  return moab::MB_FAILURE; \
47  }
48 
49 int moab::TempestOnlineMap::rearrange_arrays_by_dofs( const std::vector< unsigned int >& gdofmap,
50  DataArray1D< double >& vecFaceArea,
51  DataArray1D< double >& dCenterLon,
52  DataArray1D< double >& dCenterLat,
53  DataArray2D< double >& dVertexLon,
54  DataArray2D< double >& dVertexLat,
55  std::vector< int >& masks,
56  unsigned& N, // will have the local, after
57  int nv,
58  int& maxdof )
59 {
60  // first decide maxdof, for partitioning
61  unsigned int localmax = 0;
62  for( unsigned i = 0; i < N; i++ )
63  if( gdofmap[i] > localmax ) localmax = gdofmap[i];
64 
65  // decide partitioning based on maxdof/size
66  MPI_Allreduce( &localmax, &maxdof, 1, MPI_INT, MPI_MAX, m_pcomm->comm() );
67  // maxdof is 0 based, so actual number is +1
68  // maxdof
69  int size_per_task = ( maxdof + 1 ) / size; // based on this, processor to process dof x is x/size_per_task
70  // so we decide to reorder by actual dof, such that task 0 has dofs from [0 to size_per_task), etc
71  moab::TupleList tl;
72  unsigned numr = 2 * nv + 3; // doubles: area, centerlon, center lat, nv (vertex lon, vertex lat)
73  tl.initialize( 3, 0, 0, numr, N ); // to proc, dof, then
74  tl.enableWriteAccess();
75  // populate
76  for( unsigned i = 0; i < N; i++ )
77  {
78  int gdof = gdofmap[i];
79  int to_proc = gdof / size_per_task;
80  int mask = masks[i];
81  if( to_proc >= size ) to_proc = size - 1; // the last ones go to last proc
82  int n = tl.get_n();
83  tl.vi_wr[3 * n] = to_proc;
84  tl.vi_wr[3 * n + 1] = gdof;
85  tl.vi_wr[3 * n + 2] = mask;
86  tl.vr_wr[n * numr] = vecFaceArea[i];
87  tl.vr_wr[n * numr + 1] = dCenterLon[i];
88  tl.vr_wr[n * numr + 2] = dCenterLat[i];
89  for( int j = 0; j < nv; j++ )
90  {
91  tl.vr_wr[n * numr + 3 + j] = dVertexLon[i][j];
92  tl.vr_wr[n * numr + 3 + nv + j] = dVertexLat[i][j];
93  }
94  tl.inc_n();
95  }
96 
97  // now do the heavy communication
98  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl, 0 );
99 
100  // after communication, on each processor we should have tuples coming in
101  // still need to order by global dofs; then rearrange input vectors
102  moab::TupleList::buffer sort_buffer;
103  sort_buffer.buffer_init( tl.get_n() );
104  tl.sort( 1, &sort_buffer );
105  // count how many are unique, and collapse
106  int nb_unique = 1;
107  for( unsigned i = 0; i < tl.get_n() - 1; i++ )
108  {
109  if( tl.vi_wr[3 * i + 1] != tl.vi_wr[3 * i + 4] ) nb_unique++;
110  }
111  vecFaceArea.Allocate( nb_unique );
112  dCenterLon.Allocate( nb_unique );
113  dCenterLat.Allocate( nb_unique );
114  dVertexLon.Allocate( nb_unique, nv );
115  dVertexLat.Allocate( nb_unique, nv );
116  masks.resize( nb_unique );
117  int current_size = 1;
118  vecFaceArea[0] = tl.vr_wr[0];
119  dCenterLon[0] = tl.vr_wr[1];
120  dCenterLat[0] = tl.vr_wr[2];
121  masks[0] = tl.vi_wr[2];
122  for( int j = 0; j < nv; j++ )
123  {
124  dVertexLon[0][j] = tl.vr_wr[3 + j];
125  dVertexLat[0][j] = tl.vr_wr[3 + nv + j];
126  }
127  for( unsigned i = 0; i < tl.get_n() - 1; i++ )
128  {
129  int i1 = i + 1;
130  if( tl.vi_wr[3 * i + 1] != tl.vi_wr[3 * i + 4] )
131  {
132  vecFaceArea[current_size] = tl.vr_wr[i1 * numr];
133  dCenterLon[current_size] = tl.vr_wr[i1 * numr + 1];
134  dCenterLat[current_size] = tl.vr_wr[i1 * numr + 2];
135  for( int j = 0; j < nv; j++ )
136  {
137  dVertexLon[current_size][j] = tl.vr_wr[i1 * numr + 3 + j];
138  dVertexLat[current_size][j] = tl.vr_wr[i1 * numr + 3 + nv + j];
139  }
140  masks[current_size] = tl.vi_wr[3 * i1 + 2];
141  current_size++;
142  }
143  else
144  {
145  vecFaceArea[current_size - 1] += tl.vr_wr[i1 * numr]; // accumulate areas; will come here only for cgll ?
146  }
147  }
148 
149  N = current_size; // or nb_unique, should be the same
150  return 0;
151 }
152 #endif
153 
154 ///////////////////////////////////////////////////////////////////////////////
155 
156 moab::ErrorCode moab::TempestOnlineMap::WriteParallelMap( const std::string& strFilename,
157  const std::map< std::string, std::string >& attrMap )
158 {
159  moab::ErrorCode rval;
160 
161  size_t lastindex = strFilename.find_last_of( "." );
162  std::string extension = strFilename.substr( lastindex + 1, strFilename.size() );
163 
164  // Write the map file to disk in parallel
165  if( extension == "nc" )
166  {
167  /* Invoke the actual call to write the parallel map to disk in SCRIP format */
168  rval = this->WriteSCRIPMapFile( strFilename.c_str(), attrMap );MB_CHK_ERR( rval );
169  }
170  else
171  {
172  /* Write to the parallel H5M format */
173  rval = this->WriteHDF5MapFile( strFilename.c_str() );MB_CHK_ERR( rval );
174  }
175 
176  return rval;
177 }
178 
179 ///////////////////////////////////////////////////////////////////////////////
180 
181 moab::ErrorCode moab::TempestOnlineMap::WriteSCRIPMapFile( const std::string& strFilename,
182  const std::map< std::string, std::string >& attrMap )
183 {
184  NcError error( NcError::silent_nonfatal );
185 
186 #ifdef MOAB_HAVE_NETCDFPAR
187  bool is_independent = true;
188  ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcFile::Replace, NcFile::Netcdf4 );
189  // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
190 #else
191  NcFile ncMap( strFilename.c_str(), NcFile::Replace );
192 #endif
193 
194  if( !ncMap.is_valid() )
195  {
196  _EXCEPTION1( "Unable to open output map file \"%s\"", strFilename.c_str() );
197  }
198 
199  // Attributes
200  // ncMap.add_att( "Title", "MOAB-TempestRemap Online Regridding Weight Generator" );
201  auto it = attrMap.begin();
202  while( it != attrMap.end() )
203  {
204  // set the map attributes
205  ncMap.add_att( it->first.c_str(), it->second.c_str() );
206  // increment iterator
207  it++;
208  }
209 
210  /**
211  * Need to get the global maximum of number of vertices per element
212  * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
213  *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
214  *when writing this out, other processes may have a different size for the same array. This is
215  *hence a mess to consolidate in h5mtoscrip eventually.
216  **/
217 
218  /* Let us compute all relevant data for the current original source mesh on the process */
219  DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
220  DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
221  DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
222  if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
223  {
224  this->InitializeCoordinatesFromMeshFV(
225  *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
226  ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
227  m_remapper->max_source_edges );
228 
229  vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
230  for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
231  vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
232  }
233  else
234  {
235  DataArray3D< double > dataGLLJacobianSrc;
236  this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
237  dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
238 
239  // Generate the continuous Jacobian for input mesh
240  GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
241 
242  if( m_srcDiscType == DiscretizationType_CGLL )
243  {
244  GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
245  }
246  else
247  {
248  GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
249  }
250  }
251 
252  if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
253  {
254  this->InitializeCoordinatesFromMeshFV(
255  *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
256  ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
257  m_remapper->max_target_edges );
258 
259  vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
260  for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
261  {
262  vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
263  }
264  }
265  else
266  {
267  DataArray3D< double > dataGLLJacobianDest;
268  this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
269  dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
270 
271  // Generate the continuous Jacobian for input mesh
272  GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
273 
274  if( m_destDiscType == DiscretizationType_CGLL )
275  {
276  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
277  }
278  else
279  {
280  GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
281  }
282  }
283 
284  // Map dimensions
285  unsigned nA = ( vecSourceFaceArea.GetRows() );
286  unsigned nB = ( vecTargetFaceArea.GetRows() );
287 
288  std::vector< int > masksA, masksB;
289  ErrorCode rval = m_remapper->GetIMasks( moab::Remapper::SourceMesh, masksA );MB_CHK_SET_ERR( rval, "Trouble getting masks for source" );
290  rval = m_remapper->GetIMasks( moab::Remapper::TargetMesh, masksB );MB_CHK_SET_ERR( rval, "Trouble getting masks for target" );
291 
292  // Number of nodes per Face
293  int nSourceNodesPerFace = dSourceVertexLon.GetColumns();
294  int nTargetNodesPerFace = dTargetVertexLon.GetColumns();
295 
296  // if source or target cells have triangles at poles, center of those triangles need to come from
297  // the original quad, not from center in 3d, converted to 2d again
298  // start copy OnlineMap.cpp tempestremap
299  // right now, do this only for source mesh; copy the logic for target mesh
300  for( unsigned i = 0; i < nA; i++ )
301  {
302  const Face& face = m_meshInput->faces[i];
303 
304  int nNodes = face.edges.size();
305  int indexNodeAtPole = -1;
306  if( 3 == nNodes ) // check if one node at the poles
307  {
308  for( int j = 0; j < nNodes; j++ )
309  if( fabs( fabs( dSourceVertexLat[i][j] ) - 90.0 ) < 1.0e-12 )
310  {
311  indexNodeAtPole = j;
312  break;
313  }
314  }
315  if( indexNodeAtPole < 0 ) continue; // continue i loop, do nothing
316  // recompute center of cell, from 3d data; add one 2 nodes at pole, and average
317  int nodeAtPole = face[indexNodeAtPole]; // use the overloaded operator
318  Node nodePole = m_meshInput->nodes[nodeAtPole];
319  Node newCenter = nodePole * 2;
320  for( int j = 1; j < nNodes; j++ )
321  {
322  int indexi = ( indexNodeAtPole + j ) % nNodes; // nNodes is 3 !
323  const Node& node = m_meshInput->nodes[face[indexi]];
324  newCenter = newCenter + node;
325  }
326  newCenter = newCenter * 0.25;
327  newCenter = newCenter.Normalized();
328 
329 #ifdef VERBOSE
330  double iniLon = dSourceCenterLon[i], iniLat = dSourceCenterLat[i];
331 #endif
332  // dSourceCenterLon, dSourceCenterLat
333  XYZtoRLL_Deg( newCenter.x, newCenter.y, newCenter.z, dSourceCenterLon[i], dSourceCenterLat[i] );
334 #ifdef VERBOSE
335  std::cout << " modify center of triangle from " << iniLon << " " << iniLat << " to " << dSourceCenterLon[i]
336  << " " << dSourceCenterLat[i] << "\n";
337 #endif
338  }
339 
340  // first move data if in parallel
341 #if defined( MOAB_HAVE_MPI )
342  int max_row_dof, max_col_dof; // output; arrays will be re-distributed in chunks [maxdof/size]
343  // if (size > 1)
344  {
345  int ierr = rearrange_arrays_by_dofs( srccol_gdofmap, vecSourceFaceArea, dSourceCenterLon, dSourceCenterLat,
346  dSourceVertexLon, dSourceVertexLat, masksA, nA, nSourceNodesPerFace,
347  max_col_dof ); // now nA will be close to maxdof/size
348  if( ierr != 0 )
349  {
350  _EXCEPTION1( "Unable to arrange source data %d ", nA );
351  }
352  // rearrange target data: (nB)
353  //
354  ierr = rearrange_arrays_by_dofs( row_gdofmap, vecTargetFaceArea, dTargetCenterLon, dTargetCenterLat,
355  dTargetVertexLon, dTargetVertexLat, masksB, nB, nTargetNodesPerFace,
356  max_row_dof ); // now nA will be close to maxdof/size
357  if( ierr != 0 )
358  {
359  _EXCEPTION1( "Unable to arrange target data %d ", nB );
360  }
361  }
362 #endif
363 
364  // Number of non-zeros in the remap matrix operator
365  int nS = m_weightMatrix.nonZeros();
366 
367 #if defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_NETCDFPAR )
368  int locbuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
369  int offbuf[3] = { 0, 0, 0 };
370  int globuf[5] = { 0, 0, 0, 0, 0 };
371  MPI_Scan( locbuf, offbuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
372  MPI_Allreduce( locbuf, globuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
373  MPI_Allreduce( &locbuf[3], &globuf[3], 2, MPI_INT, MPI_MAX, m_pcomm->comm() );
374 
375  // MPI_Scan is inclusive of data in current rank; modify accordingly.
376  offbuf[0] -= nA;
377  offbuf[1] -= nB;
378  offbuf[2] -= nS;
379 
380 #else
381  int offbuf[3] = { 0, 0, 0 };
382  int globuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
383 #endif
384 
385  std::vector< std::string > srcdimNames, tgtdimNames;
386  std::vector< int > srcdimSizes, tgtdimSizes;
387  {
388  if( m_remapper->m_source_type == moab::TempestRemapper::RLL && m_remapper->m_source_metadata.size() )
389  {
390  srcdimNames.push_back( "lat" );
391  srcdimNames.push_back( "lon" );
392  srcdimSizes.resize( 2, 0 );
393  srcdimSizes[0] = m_remapper->m_source_metadata[0];
394  srcdimSizes[1] = m_remapper->m_source_metadata[1];
395  }
396  else
397  {
398  srcdimNames.push_back( "num_elem" );
399  srcdimSizes.push_back( globuf[0] );
400  }
401 
402  if( m_remapper->m_target_type == moab::TempestRemapper::RLL && m_remapper->m_target_metadata.size() )
403  {
404  tgtdimNames.push_back( "lat" );
405  tgtdimNames.push_back( "lon" );
406  tgtdimSizes.resize( 2, 0 );
407  tgtdimSizes[0] = m_remapper->m_target_metadata[0];
408  tgtdimSizes[1] = m_remapper->m_target_metadata[1];
409  }
410  else
411  {
412  tgtdimNames.push_back( "num_elem" );
413  tgtdimSizes.push_back( globuf[1] );
414  }
415  }
416 
417  // Write output dimensions entries
418  unsigned nSrcGridDims = ( srcdimSizes.size() );
419  unsigned nDstGridDims = ( tgtdimSizes.size() );
420 
421  NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
422  NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
423 
424  NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
425  NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
426 
427 #ifdef MOAB_HAVE_NETCDFPAR
428  ncMap.enable_var_par_access( varSrcGridDims, is_independent );
429  ncMap.enable_var_par_access( varDstGridDims, is_independent );
430 #endif
431 
432  // write dimension names
433  {
434  char szDim[64];
435  for( unsigned i = 0; i < srcdimSizes.size(); i++ )
436  {
437  varSrcGridDims->set_cur( nSrcGridDims - i - 1 );
438  varSrcGridDims->put( &( srcdimSizes[nSrcGridDims - i - 1] ), 1 );
439  }
440 
441  for( unsigned i = 0; i < srcdimSizes.size(); i++ )
442  {
443  snprintf( szDim, 64, "name%i", i );
444  varSrcGridDims->add_att( szDim, srcdimNames[nSrcGridDims - i - 1].c_str() );
445  }
446 
447  for( unsigned i = 0; i < tgtdimSizes.size(); i++ )
448  {
449  varDstGridDims->set_cur( nDstGridDims - i - 1 );
450  varDstGridDims->put( &( tgtdimSizes[nDstGridDims - i - 1] ), 1 );
451  }
452 
453  for( unsigned i = 0; i < tgtdimSizes.size(); i++ )
454  {
455  snprintf( szDim, 64, "name%i", i );
456  varDstGridDims->add_att( szDim, tgtdimNames[nDstGridDims - i - 1].c_str() );
457  }
458  }
459 
460  // Source and Target mesh resolutions
461  NcDim* dimNA = ncMap.add_dim( "n_a", globuf[0] );
462  NcDim* dimNB = ncMap.add_dim( "n_b", globuf[1] );
463 
464  // Number of nodes per Face
465  NcDim* dimNVA = ncMap.add_dim( "nv_a", globuf[3] );
466  NcDim* dimNVB = ncMap.add_dim( "nv_b", globuf[4] );
467 
468  // Write coordinates
469  NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA );
470  NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB );
471 
472  NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA );
473  NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB );
474 
475  NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA, dimNVA );
476  NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB, dimNVB );
477 
478  NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA, dimNVA );
479  NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB, dimNVB );
480 
481  // Write masks
482  NcVar* varMaskA = ncMap.add_var( "mask_a", ncInt, dimNA );
483  NcVar* varMaskB = ncMap.add_var( "mask_b", ncInt, dimNB );
484 
485 #ifdef MOAB_HAVE_NETCDFPAR
486  ncMap.enable_var_par_access( varYCA, is_independent );
487  ncMap.enable_var_par_access( varYCB, is_independent );
488  ncMap.enable_var_par_access( varXCA, is_independent );
489  ncMap.enable_var_par_access( varXCB, is_independent );
490  ncMap.enable_var_par_access( varYVA, is_independent );
491  ncMap.enable_var_par_access( varYVB, is_independent );
492  ncMap.enable_var_par_access( varXVA, is_independent );
493  ncMap.enable_var_par_access( varXVB, is_independent );
494  ncMap.enable_var_par_access( varMaskA, is_independent );
495  ncMap.enable_var_par_access( varMaskB, is_independent );
496 #endif
497 
498  varYCA->add_att( "units", "degrees" );
499  varYCB->add_att( "units", "degrees" );
500 
501  varXCA->add_att( "units", "degrees" );
502  varXCB->add_att( "units", "degrees" );
503 
504  varYVA->add_att( "units", "degrees" );
505  varYVB->add_att( "units", "degrees" );
506 
507  varXVA->add_att( "units", "degrees" );
508  varXVB->add_att( "units", "degrees" );
509 
510  // Verify dimensionality
511  if( dSourceCenterLon.GetRows() != nA )
512  {
513  _EXCEPTIONT( "Mismatch between dSourceCenterLon and nA" );
514  }
515  if( dSourceCenterLat.GetRows() != nA )
516  {
517  _EXCEPTIONT( "Mismatch between dSourceCenterLat and nA" );
518  }
519  if( dTargetCenterLon.GetRows() != nB )
520  {
521  _EXCEPTIONT( "Mismatch between dTargetCenterLon and nB" );
522  }
523  if( dTargetCenterLat.GetRows() != nB )
524  {
525  _EXCEPTIONT( "Mismatch between dTargetCenterLat and nB" );
526  }
527  if( dSourceVertexLon.GetRows() != nA )
528  {
529  _EXCEPTIONT( "Mismatch between dSourceVertexLon and nA" );
530  }
531  if( dSourceVertexLat.GetRows() != nA )
532  {
533  _EXCEPTIONT( "Mismatch between dSourceVertexLat and nA" );
534  }
535  if( dTargetVertexLon.GetRows() != nB )
536  {
537  _EXCEPTIONT( "Mismatch between dTargetVertexLon and nB" );
538  }
539  if( dTargetVertexLat.GetRows() != nB )
540  {
541  _EXCEPTIONT( "Mismatch between dTargetVertexLat and nB" );
542  }
543 
544  varYCA->set_cur( (long)offbuf[0] );
545  varYCA->put( &( dSourceCenterLat[0] ), nA );
546  varYCB->set_cur( (long)offbuf[1] );
547  varYCB->put( &( dTargetCenterLat[0] ), nB );
548 
549  varXCA->set_cur( (long)offbuf[0] );
550  varXCA->put( &( dSourceCenterLon[0] ), nA );
551  varXCB->set_cur( (long)offbuf[1] );
552  varXCB->put( &( dTargetCenterLon[0] ), nB );
553 
554  varYVA->set_cur( (long)offbuf[0] );
555  varYVA->put( &( dSourceVertexLat[0][0] ), nA, nSourceNodesPerFace );
556  varYVB->set_cur( (long)offbuf[1] );
557  varYVB->put( &( dTargetVertexLat[0][0] ), nB, nTargetNodesPerFace );
558 
559  varXVA->set_cur( (long)offbuf[0] );
560  varXVA->put( &( dSourceVertexLon[0][0] ), nA, nSourceNodesPerFace );
561  varXVB->set_cur( (long)offbuf[1] );
562  varXVB->put( &( dTargetVertexLon[0][0] ), nB, nTargetNodesPerFace );
563 
564  varMaskA->set_cur( (long)offbuf[0] );
565  varMaskA->put( &( masksA[0] ), nA );
566  varMaskB->set_cur( (long)offbuf[1] );
567  varMaskB->put( &( masksB[0] ), nB );
568 
569  // Write areas
570  NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
571 #ifdef MOAB_HAVE_NETCDFPAR
572  ncMap.enable_var_par_access( varAreaA, is_independent );
573 #endif
574  varAreaA->set_cur( (long)offbuf[0] );
575  varAreaA->put( &( vecSourceFaceArea[0] ), nA );
576 
577  NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
578 #ifdef MOAB_HAVE_NETCDFPAR
579  ncMap.enable_var_par_access( varAreaB, is_independent );
580 #endif
581  varAreaB->set_cur( (long)offbuf[1] );
582  varAreaB->put( &( vecTargetFaceArea[0] ), nB );
583 
584  // Write SparseMatrix entries
585  DataArray1D< int > vecRow( nS );
586  DataArray1D< int > vecCol( nS );
587  DataArray1D< double > vecS( nS );
588  DataArray1D< double > dFracA( nA );
589  DataArray1D< double > dFracB( nB );
590 
591  moab::TupleList tlValRow, tlValCol;
592  unsigned numr = 1; //
593  // value has to be sent to processor row/nB for for fracA and col/nA for fracB
594  // vecTargetArea (indexRow ) has to be sent for fracA (index col?)
595  // vecTargetFaceArea will have to be sent to col index, with its index !
596  tlValRow.initialize( 2, 0, 0, numr, nS ); // to proc(row), global row , value
597  tlValCol.initialize( 3, 0, 0, numr, nS ); // to proc(col), global row / col, value
598  tlValRow.enableWriteAccess();
599  tlValCol.enableWriteAccess();
600  /*
601  *
602  dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
603  dFracB[ row ] += val ;
604  */
605  int offset = 0;
606 #if defined( MOAB_HAVE_MPI )
607  int nAbase = ( max_col_dof + 1 ) / size; // it is nA, except last rank ( == size - 1 )
608  int nBbase = ( max_row_dof + 1 ) / size; // it is nB, except last rank ( == size - 1 )
609 #endif
610  for( int i = 0; i < m_weightMatrix.outerSize(); ++i )
611  {
612  for( WeightMatrix::InnerIterator it( m_weightMatrix, i ); it; ++it )
613  {
614  vecRow[offset] = 1 + this->GetRowGlobalDoF( it.row() ); // row index
615  vecCol[offset] = 1 + this->GetColGlobalDoF( it.col() ); // col index
616  vecS[offset] = it.value(); // value
617 
618 #if defined( MOAB_HAVE_MPI )
619  {
620  // value M(row, col) will contribute to procRow and procCol values for fracA and fracB
621  int procRow = ( vecRow[offset] - 1 ) / nBbase;
622  if( procRow >= size ) procRow = size - 1;
623  int procCol = ( vecCol[offset] - 1 ) / nAbase;
624  if( procCol >= size ) procCol = size - 1;
625  int nrInd = tlValRow.get_n();
626  tlValRow.vi_wr[2 * nrInd] = procRow;
627  tlValRow.vi_wr[2 * nrInd + 1] = vecRow[offset] - 1;
628  tlValRow.vr_wr[nrInd] = vecS[offset];
629  tlValRow.inc_n();
630  int ncInd = tlValCol.get_n();
631  tlValCol.vi_wr[3 * ncInd] = procCol;
632  tlValCol.vi_wr[3 * ncInd + 1] = vecRow[offset] - 1;
633  tlValCol.vi_wr[3 * ncInd + 2] = vecCol[offset] - 1; // this is column
634  tlValCol.vr_wr[ncInd] = vecS[offset];
635  tlValCol.inc_n();
636  }
637 
638 #endif
639  offset++;
640  }
641  }
642 #if defined( MOAB_HAVE_MPI )
643  // need to send values for their row and col processors, to compute fractions there
644  // now do the heavy communication
645  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValCol, 0 );
646  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValRow, 0 );
647 
648  // we have now, for example, dFracB[ row ] += val ;
649  // so we know that on current task, we received tlValRow
650  // reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
651  // dFracB[ row ] += val ;
652  for( unsigned i = 0; i < tlValRow.get_n(); i++ )
653  {
654  // int fromProc = tlValRow.vi_wr[2 * i];
655  int gRowInd = tlValRow.vi_wr[2 * i + 1];
656  int localIndexRow = gRowInd - nBbase * rank; // modulo nBbase rank is from 0 to size - 1;
657  double wgt = tlValRow.vr_wr[i];
658  assert( localIndexRow >= 0 );
659  assert( nB - localIndexRow > 0 );
660  dFracB[localIndexRow] += wgt;
661  }
662  // to compute dFracA we need vecTargetFaceArea[ row ]; we know the row, and we can get the proc we need it from
663 
664  std::set< int > neededRows;
665  for( unsigned i = 0; i < tlValCol.get_n(); i++ )
666  {
667  int rRowInd = tlValCol.vi_wr[3 * i + 1];
668  neededRows.insert( rRowInd );
669  // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
670  }
671  moab::TupleList tgtAreaReq;
672  tgtAreaReq.initialize( 2, 0, 0, 0, neededRows.size() );
673  tgtAreaReq.enableWriteAccess();
674  for( std::set< int >::iterator sit = neededRows.begin(); sit != neededRows.end(); sit++ )
675  {
676  int neededRow = *sit;
677  int procRow = neededRow / nBbase;
678  if( procRow >= size ) procRow = size - 1;
679  int nr = tgtAreaReq.get_n();
680  tgtAreaReq.vi_wr[2 * nr] = procRow;
681  tgtAreaReq.vi_wr[2 * nr + 1] = neededRow;
682  tgtAreaReq.inc_n();
683  }
684 
685  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaReq, 0 );
686  // we need to send back the tgtArea corresponding to row
687  moab::TupleList tgtAreaInfo; // load it with tgtArea at row
688  tgtAreaInfo.initialize( 2, 0, 0, 1, tgtAreaReq.get_n() );
689  tgtAreaInfo.enableWriteAccess();
690  for( unsigned i = 0; i < tgtAreaReq.get_n(); i++ )
691  {
692  int from_proc = tgtAreaReq.vi_wr[2 * i];
693  int row = tgtAreaReq.vi_wr[2 * i + 1];
694  int locaIndexRow = row - rank * nBbase;
695  double areaToSend = vecTargetFaceArea[locaIndexRow];
696  // int remoteIndex = tgtAreaReq.vi_wr[3*i + 2] ;
697 
698  tgtAreaInfo.vi_wr[2 * i] = from_proc; // send back requested info
699  tgtAreaInfo.vi_wr[2 * i + 1] = row;
700  tgtAreaInfo.vr_wr[i] = areaToSend; // this will be tgt area at row
701  tgtAreaInfo.inc_n();
702  }
703  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaInfo, 0 );
704 
705  std::map< int, double > areaAtRow;
706  for( unsigned i = 0; i < tgtAreaInfo.get_n(); i++ )
707  {
708  // we have received from proc, value for row !
709  int row = tgtAreaInfo.vi_wr[2 * i + 1];
710  areaAtRow[row] = tgtAreaInfo.vr_wr[i];
711  }
712 
713  // we have now for rows the
714  // it is ordered by index, so:
715  // now compute reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
716  // tgtAreaInfo will have at index i the area we need (from row)
717  // there should be an easier way :(
718  for( unsigned i = 0; i < tlValCol.get_n(); i++ )
719  {
720  int rRowInd = tlValCol.vi_wr[3 * i + 1];
721  int colInd = tlValCol.vi_wr[3 * i + 2];
722  double val = tlValCol.vr_wr[i];
723  int localColInd = colInd - rank * nAbase; // < local nA
724  // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
725  auto itMap = areaAtRow.find( rRowInd ); // it should be different from end
726  if( itMap != areaAtRow.end() )
727  {
728  double areaRow = itMap->second; // we fished a lot for this !
729  dFracA[localColInd] += val / vecSourceFaceArea[localColInd] * areaRow;
730  }
731  }
732 
733 #endif
734  // Load in data
735  NcDim* dimNS = ncMap.add_dim( "n_s", globuf[2] );
736 
737  NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
738  NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
739  NcVar* varS = ncMap.add_var( "S", ncDouble, dimNS );
740 #ifdef MOAB_HAVE_NETCDFPAR
741  ncMap.enable_var_par_access( varRow, is_independent );
742  ncMap.enable_var_par_access( varCol, is_independent );
743  ncMap.enable_var_par_access( varS, is_independent );
744 #endif
745 
746  varRow->set_cur( (long)offbuf[2] );
747  varRow->put( vecRow, nS );
748 
749  varCol->set_cur( (long)offbuf[2] );
750  varCol->put( vecCol, nS );
751 
752  varS->set_cur( (long)offbuf[2] );
753  varS->put( &( vecS[0] ), nS );
754 
755  // Calculate and write fractional coverage arrays
756  NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
757 #ifdef MOAB_HAVE_NETCDFPAR
758  ncMap.enable_var_par_access( varFracA, is_independent );
759 #endif
760  varFracA->add_att( "name", "fraction of target coverage of source dof" );
761  varFracA->add_att( "units", "unitless" );
762  varFracA->set_cur( (long)offbuf[0] );
763  varFracA->put( &( dFracA[0] ), nA );
764 
765  NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
766 #ifdef MOAB_HAVE_NETCDFPAR
767  ncMap.enable_var_par_access( varFracB, is_independent );
768 #endif
769  varFracB->add_att( "name", "fraction of source coverage of target dof" );
770  varFracB->add_att( "units", "unitless" );
771  varFracB->set_cur( (long)offbuf[1] );
772  varFracB->put( &( dFracB[0] ), nB );
773 
774  // Add global attributes
775  // std::map<std::string, std::string>::const_iterator iterAttributes =
776  // mapAttributes.begin();
777  // for (; iterAttributes != mapAttributes.end(); iterAttributes++) {
778  // ncMap.add_att(
779  // iterAttributes->first.c_str(),
780  // iterAttributes->second.c_str());
781  // }
782 
783  ncMap.close();
784 
785  return moab::MB_SUCCESS;
786 }
787 
788 ///////////////////////////////////////////////////////////////////////////////
789 
790 moab::ErrorCode moab::TempestOnlineMap::WriteHDF5MapFile( const std::string& strOutputFile )
791 {
792  moab::ErrorCode rval;
793 
794  /**
795  * Need to get the global maximum of number of vertices per element
796  * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
797  *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
798  *when writing this out, other processes may have a different size for the same array. This is
799  *hence a mess to consolidate in h5mtoscrip eventually.
800  **/
801 
802  /* Let us compute all relevant data for the current original source mesh on the process */
803  DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
804  DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
805  DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
806  if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
807  {
808  this->InitializeCoordinatesFromMeshFV(
809  *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
810  ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
811  m_remapper->max_source_edges );
812 
813  vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
814  for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
815  vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
816  }
817  else
818  {
819  DataArray3D< double > dataGLLJacobianSrc;
820  this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
821  dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
822 
823  // Generate the continuous Jacobian for input mesh
824  GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
825 
826  if( m_srcDiscType == DiscretizationType_CGLL )
827  {
828  GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
829  }
830  else
831  {
832  GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
833  }
834  }
835 
836  if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
837  {
838  this->InitializeCoordinatesFromMeshFV(
839  *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
840  ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
841  m_remapper->max_target_edges );
842 
843  vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
844  for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
845  vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
846  }
847  else
848  {
849  DataArray3D< double > dataGLLJacobianDest;
850  this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
851  dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
852 
853  // Generate the continuous Jacobian for input mesh
854  GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
855 
856  if( m_destDiscType == DiscretizationType_CGLL )
857  {
858  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
859  }
860  else
861  {
862  GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
863  }
864  }
865 
866  moab::EntityHandle& m_meshOverlapSet = m_remapper->m_overlap_set;
867  int tot_src_ents = m_remapper->m_source_entities.size();
868  int tot_tgt_ents = m_remapper->m_target_entities.size();
869  int tot_src_size = dSourceCenterLon.GetRows();
870  int tot_tgt_size = m_dTargetCenterLon.GetRows();
871  int tot_vsrc_size = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
872  int tot_vtgt_size = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
873 
874  const int weightMatNNZ = m_weightMatrix.nonZeros();
875  moab::Tag tagMapMetaData, tagMapIndexRow, tagMapIndexCol, tagMapValues, srcEleIDs, tgtEleIDs;
876  rval = m_interface->tag_get_handle( "SMAT_DATA", 13, moab::MB_TYPE_INTEGER, tagMapMetaData,
877  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
878  rval = m_interface->tag_get_handle( "SMAT_ROWS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexRow,
879  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
880  rval = m_interface->tag_get_handle( "SMAT_COLS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexCol,
881  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
882  rval = m_interface->tag_get_handle( "SMAT_VALS", weightMatNNZ, moab::MB_TYPE_DOUBLE, tagMapValues,
883  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
884  rval = m_interface->tag_get_handle( "SourceGIDS", tot_src_size, moab::MB_TYPE_INTEGER, srcEleIDs,
885  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
886  rval = m_interface->tag_get_handle( "TargetGIDS", tot_tgt_size, moab::MB_TYPE_INTEGER, tgtEleIDs,
887  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
888  moab::Tag srcAreaValues, tgtAreaValues;
889  rval = m_interface->tag_get_handle( "SourceAreas", tot_src_size, moab::MB_TYPE_DOUBLE, srcAreaValues,
890  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
891  rval = m_interface->tag_get_handle( "TargetAreas", tot_tgt_size, moab::MB_TYPE_DOUBLE, tgtAreaValues,
892  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
893  moab::Tag tagSrcCoordsCLon, tagSrcCoordsCLat, tagTgtCoordsCLon, tagTgtCoordsCLat;
894  rval = m_interface->tag_get_handle( "SourceCoordCenterLon", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLon,
895  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
896  rval = m_interface->tag_get_handle( "SourceCoordCenterLat", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLat,
897  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
898  rval = m_interface->tag_get_handle( "TargetCoordCenterLon", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLon,
899  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
900  rval = m_interface->tag_get_handle( "TargetCoordCenterLat", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLat,
901  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
902  moab::Tag tagSrcCoordsVLon, tagSrcCoordsVLat, tagTgtCoordsVLon, tagTgtCoordsVLat;
903  rval = m_interface->tag_get_handle( "SourceCoordVertexLon", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLon,
904  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
905  rval = m_interface->tag_get_handle( "SourceCoordVertexLat", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLat,
906  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
907  rval = m_interface->tag_get_handle( "TargetCoordVertexLon", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLon,
908  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
909  rval = m_interface->tag_get_handle( "TargetCoordVertexLat", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLat,
910  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
911  moab::Tag srcMaskValues, tgtMaskValues;
912  if( m_iSourceMask.IsAttached() )
913  {
914  rval = m_interface->tag_get_handle( "SourceMask", m_iSourceMask.GetRows(), moab::MB_TYPE_INTEGER, srcMaskValues,
915  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
916  }
917  if( m_iTargetMask.IsAttached() )
918  {
919  rval = m_interface->tag_get_handle( "TargetMask", m_iTargetMask.GetRows(), moab::MB_TYPE_INTEGER, tgtMaskValues,
920  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
921  }
922 
923  std::vector< int > smatrowvals( weightMatNNZ ), smatcolvals( weightMatNNZ );
924  std::vector< double > smatvals( weightMatNNZ );
925  // const double* smatvals = m_weightMatrix.valuePtr();
926  // Loop over the matrix entries and find the max global ID for rows and columns
927  for( int k = 0, offset = 0; k < m_weightMatrix.outerSize(); ++k )
928  {
929  for( moab::TempestOnlineMap::WeightMatrix::InnerIterator it( m_weightMatrix, k ); it; ++it, ++offset )
930  {
931  smatrowvals[offset] = this->GetRowGlobalDoF( it.row() );
932  smatcolvals[offset] = this->GetColGlobalDoF( it.col() );
933  smatvals[offset] = it.value();
934  }
935  }
936 
937  /* Set the global IDs for the DoFs */
938  ////
939  // col_gdofmap [ col_ldofmap [ 0 : local_ndofs ] ] = GDOF
940  // row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
941  ////
942  int maxrow = 0, maxcol = 0;
943  std::vector< int > src_global_dofs( tot_src_size ), tgt_global_dofs( tot_tgt_size );
944  for( int i = 0; i < tot_src_size; ++i )
945  {
946  src_global_dofs[i] = srccol_gdofmap[i];
947  maxcol = ( src_global_dofs[i] > maxcol ) ? src_global_dofs[i] : maxcol;
948  }
949 
950  for( int i = 0; i < tot_tgt_size; ++i )
951  {
952  tgt_global_dofs[i] = row_gdofmap[i];
953  maxrow = ( tgt_global_dofs[i] > maxrow ) ? tgt_global_dofs[i] : maxrow;
954  }
955 
956  ///////////////////////////////////////////////////////////////////////////
957  // The metadata in H5M file contains the following data:
958  //
959  // 1. n_a: Total source entities: (number of elements in source mesh)
960  // 2. n_b: Total target entities: (number of elements in target mesh)
961  // 3. nv_a: Max edge size of elements in source mesh
962  // 4. nv_b: Max edge size of elements in target mesh
963  // 5. maxrows: Number of rows in remap weight matrix
964  // 6. maxcols: Number of cols in remap weight matrix
965  // 7. nnz: Number of total nnz in sparse remap weight matrix
966  // 8. np_a: The order of the field description on the source mesh: >= 1
967  // 9. np_b: The order of the field description on the target mesh: >= 1
968  // 10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2 =
969  // dGLL]
970  // 11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2 =
971  // dGLL]
972  // 12. conserved: Flag to specify whether the remap operator has conservation constraints: [0,
973  // 1]
974  // 13. monotonicity: Flags to specify whether the remap operator has monotonicity constraints:
975  // [0, 1, 2]
976  //
977  ///////////////////////////////////////////////////////////////////////////
978  int map_disc_details[6];
979  map_disc_details[0] = m_nDofsPEl_Src;
980  map_disc_details[1] = m_nDofsPEl_Dest;
981  map_disc_details[2] = ( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD
982  ? 0
983  : ( m_srcDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
984  map_disc_details[3] = ( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD
985  ? 0
986  : ( m_destDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
987  map_disc_details[4] = ( m_bConserved ? 1 : 0 );
988  map_disc_details[5] = m_iMonotonicity;
989 
990 #ifdef MOAB_HAVE_MPI
991  int loc_smatmetadata[13] = { tot_src_ents,
992  tot_tgt_ents,
993  m_remapper->max_source_edges,
994  m_remapper->max_target_edges,
995  maxrow + 1,
996  maxcol + 1,
997  weightMatNNZ,
998  map_disc_details[0],
999  map_disc_details[1],
1000  map_disc_details[2],
1001  map_disc_details[3],
1002  map_disc_details[4],
1003  map_disc_details[5] };
1004  rval = m_interface->tag_set_data( tagMapMetaData, &m_meshOverlapSet, 1, &loc_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1005  int glb_smatmetadata[13] = { 0,
1006  0,
1007  0,
1008  0,
1009  0,
1010  0,
1011  0,
1012  map_disc_details[0],
1013  map_disc_details[1],
1014  map_disc_details[2],
1015  map_disc_details[3],
1016  map_disc_details[4],
1017  map_disc_details[5] };
1018  int loc_buf[7] = {
1019  tot_src_ents, tot_tgt_ents, weightMatNNZ, m_remapper->max_source_edges, m_remapper->max_target_edges,
1020  maxrow, maxcol };
1021  int glb_buf[4] = { 0, 0, 0, 0 };
1022  MPI_Reduce( &loc_buf[0], &glb_buf[0], 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
1023  glb_smatmetadata[0] = glb_buf[0];
1024  glb_smatmetadata[1] = glb_buf[1];
1025  glb_smatmetadata[6] = glb_buf[2];
1026  MPI_Reduce( &loc_buf[3], &glb_buf[0], 4, MPI_INT, MPI_MAX, 0, m_pcomm->comm() );
1027  glb_smatmetadata[2] = glb_buf[0];
1028  glb_smatmetadata[3] = glb_buf[1];
1029  glb_smatmetadata[4] = glb_buf[2];
1030  glb_smatmetadata[5] = glb_buf[3];
1031 #else
1032  int glb_smatmetadata[13] = { tot_src_ents,
1033  tot_tgt_ents,
1034  m_remapper->max_source_edges,
1035  m_remapper->max_target_edges,
1036  maxrow,
1037  maxcol,
1038  weightMatNNZ,
1039  map_disc_details[0],
1040  map_disc_details[1],
1041  map_disc_details[2],
1042  map_disc_details[3],
1043  map_disc_details[4],
1044  map_disc_details[5] };
1045 #endif
1046  // These values represent number of rows and columns. So should be 1-based.
1047  glb_smatmetadata[4]++;
1048  glb_smatmetadata[5]++;
1049 
1050  if( this->is_root )
1051  {
1052  std::cout << " " << this->rank << " Writing remap weights with size [" << glb_smatmetadata[4] << " X "
1053  << glb_smatmetadata[5] << "] and NNZ = " << glb_smatmetadata[6] << std::endl;
1054  EntityHandle root_set = 0;
1055  rval = m_interface->tag_set_data( tagMapMetaData, &root_set, 1, &glb_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1056  }
1057 
1058  int dsize;
1059  const int numval = weightMatNNZ;
1060  const void* smatrowvals_d = smatrowvals.data();
1061  const void* smatcolvals_d = smatcolvals.data();
1062  const void* smatvals_d = smatvals.data();
1063  rval = m_interface->tag_set_by_ptr( tagMapIndexRow, &m_meshOverlapSet, 1, &smatrowvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1064  rval = m_interface->tag_set_by_ptr( tagMapIndexCol, &m_meshOverlapSet, 1, &smatcolvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1065  rval = m_interface->tag_set_by_ptr( tagMapValues, &m_meshOverlapSet, 1, &smatvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1066 
1067  /* Set the global IDs for the DoFs */
1068  const void* srceleidvals_d = src_global_dofs.data();
1069  const void* tgteleidvals_d = tgt_global_dofs.data();
1070  dsize = src_global_dofs.size();
1071  rval = m_interface->tag_set_by_ptr( srcEleIDs, &m_meshOverlapSet, 1, &srceleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1072  dsize = tgt_global_dofs.size();
1073  rval = m_interface->tag_set_by_ptr( tgtEleIDs, &m_meshOverlapSet, 1, &tgteleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1074 
1075  /* Set the source and target areas */
1076  const void* srcareavals_d = vecSourceFaceArea;
1077  const void* tgtareavals_d = vecTargetFaceArea;
1078  dsize = tot_src_size;
1079  rval = m_interface->tag_set_by_ptr( srcAreaValues, &m_meshOverlapSet, 1, &srcareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1080  dsize = tot_tgt_size;
1081  rval = m_interface->tag_set_by_ptr( tgtAreaValues, &m_meshOverlapSet, 1, &tgtareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1082 
1083  /* Set the coordinates for source and target center vertices */
1084  const void* srccoordsclonvals_d = &dSourceCenterLon[0];
1085  const void* srccoordsclatvals_d = &dSourceCenterLat[0];
1086  dsize = dSourceCenterLon.GetRows();
1087  rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLon, &m_meshOverlapSet, 1, &srccoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1088  rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLat, &m_meshOverlapSet, 1, &srccoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1089  const void* tgtcoordsclonvals_d = &m_dTargetCenterLon[0];
1090  const void* tgtcoordsclatvals_d = &m_dTargetCenterLat[0];
1091  dsize = vecTargetFaceArea.GetRows();
1092  rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLon, &m_meshOverlapSet, 1, &tgtcoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1093  rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLat, &m_meshOverlapSet, 1, &tgtcoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1094 
1095  /* Set the coordinates for source and target element vertices */
1096  const void* srccoordsvlonvals_d = &( dSourceVertexLon[0][0] );
1097  const void* srccoordsvlatvals_d = &( dSourceVertexLat[0][0] );
1098  dsize = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
1099  rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLon, &m_meshOverlapSet, 1, &srccoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1100  rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLat, &m_meshOverlapSet, 1, &srccoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1101  const void* tgtcoordsvlonvals_d = &( m_dTargetVertexLon[0][0] );
1102  const void* tgtcoordsvlatvals_d = &( m_dTargetVertexLat[0][0] );
1103  dsize = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
1104  rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLon, &m_meshOverlapSet, 1, &tgtcoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1105  rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLat, &m_meshOverlapSet, 1, &tgtcoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1106 
1107  /* Set the masks for source and target meshes if available */
1108  if( m_iSourceMask.IsAttached() )
1109  {
1110  const void* srcmaskvals_d = m_iSourceMask;
1111  dsize = m_iSourceMask.GetRows();
1112  rval = m_interface->tag_set_by_ptr( srcMaskValues, &m_meshOverlapSet, 1, &srcmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1113  }
1114 
1115  if( m_iTargetMask.IsAttached() )
1116  {
1117  const void* tgtmaskvals_d = m_iTargetMask;
1118  dsize = m_iTargetMask.GetRows();
1119  rval = m_interface->tag_set_by_ptr( tgtMaskValues, &m_meshOverlapSet, 1, &tgtmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1120  }
1121 
1122 #ifdef MOAB_HAVE_MPI
1123  const char* writeOptions = ( this->size > 1 ? "PARALLEL=WRITE_PART" : "" );
1124 #else
1125  const char* writeOptions = "";
1126 #endif
1127 
1128  // EntityHandle sets[3] = {m_remapper->m_source_set, m_remapper->m_target_set, m_remapper->m_overlap_set};
1129  EntityHandle sets[1] = { m_remapper->m_overlap_set };
1130  rval = m_interface->write_file( strOutputFile.c_str(), NULL, writeOptions, sets, 1 );MB_CHK_ERR( rval );
1131 
1132 #ifdef WRITE_SCRIP_FILE
1133  sstr.str( "" );
1134  sstr << ctx.outFilename.substr( 0, lastindex ) << "_" << proc_id << ".nc";
1135  std::map< std::string, std::string > mapAttributes;
1136  mapAttributes["Creator"] = "MOAB mbtempest workflow";
1137  if( !ctx.proc_id ) std::cout << "Writing offline map to file: " << sstr.str() << std::endl;
1138  this->Write( strOutputFile.c_str(), mapAttributes, NcFile::Netcdf4 );
1139  sstr.str( "" );
1140 #endif
1141 
1142  return moab::MB_SUCCESS;
1143 }
1144 
1145 ///////////////////////////////////////////////////////////////////////////////
1146 
1147 void print_progress( const int barWidth, const float progress, const char* message )
1148 {
1149  std::cout << message << " [";
1150  int pos = barWidth * progress;
1151  for( int i = 0; i < barWidth; ++i )
1152  {
1153  if( i < pos )
1154  std::cout << "=";
1155  else if( i == pos )
1156  std::cout << ">";
1157  else
1158  std::cout << " ";
1159  }
1160  std::cout << "] " << int( progress * 100.0 ) << " %\r";
1161  std::cout.flush();
1162 }
1163 
1164 ///////////////////////////////////////////////////////////////////////////////
1165 moab::ErrorCode moab::TempestOnlineMap::ReadParallelMap( const char* strSource,
1166  const std::vector< int >& owned_dof_ids,
1167  bool row_partition )
1168 {
1169  NcError error( NcError::silent_nonfatal );
1170 
1171  NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
1172  int nS = 0, nA = 0, nB = 0;
1173 #ifdef MOAB_HAVE_PNETCDF
1174  // some variables will be used just in the case netcdfpar reader fails
1175  int ncfile = -1, ret = 0;
1176  int ndims, nvars, ngatts, unlimited;
1177 #endif
1178 #ifdef MOAB_HAVE_NETCDFPAR
1179  bool is_independent = true;
1180  ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
1181  // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
1182 #else
1183  NcFile ncMap( strSource, NcFile::ReadOnly );
1184 #endif
1185 
1186 #define CHECK_EXCEPTION( obj, type, varstr ) \
1187  { \
1188  if( obj == NULL ) \
1189  { \
1190  _EXCEPTION3( "Map file \"%s\" does not contain %s \"%s\"", strSource, type, varstr ); \
1191  } \
1192  }
1193 
1194  // Read SparseMatrix entries
1195 
1196  if( ncMap.is_valid() )
1197  {
1198  NcDim* dimNS = ncMap.get_dim( "n_s" );
1199  CHECK_EXCEPTION( dimNS, "dimension", "n_s" );
1200 
1201  NcDim* dimNA = ncMap.get_dim( "n_a" );
1202  CHECK_EXCEPTION( dimNA, "dimension", "n_a" );
1203 
1204  NcDim* dimNB = ncMap.get_dim( "n_b" );
1205  CHECK_EXCEPTION( dimNB, "dimension", "n_b" );
1206 
1207  // store total number of nonzeros
1208  nS = dimNS->size();
1209  nA = dimNA->size();
1210  nB = dimNB->size();
1211 
1212  varRow = ncMap.get_var( "row" );
1213  CHECK_EXCEPTION( varRow, "variable", "row" );
1214 
1215  varCol = ncMap.get_var( "col" );
1216  CHECK_EXCEPTION( varCol, "variable", "col" );
1217 
1218  varS = ncMap.get_var( "S" );
1219  CHECK_EXCEPTION( varS, "variable", "S" );
1220 
1221 #ifdef MOAB_HAVE_NETCDFPAR
1222  ncMap.enable_var_par_access( varRow, is_independent );
1223  ncMap.enable_var_par_access( varCol, is_independent );
1224  ncMap.enable_var_par_access( varS, is_independent );
1225 #endif
1226  }
1227  else
1228  {
1229 #ifdef MOAB_HAVE_PNETCDF
1230  // read the file using pnetcdf directly, in parallel; need to have MPI, we do not check that anymore
1231  // why build wth pnetcdf without MPI ?
1232  // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
1233  ret = ncmpi_open( m_pcomm->comm(), strSource, NC_NOWRITE, MPI_INFO_NULL, &ncfile );
1234  ERR_PARNC( ret ); // bail out completely
1235  ret = ncmpi_inq( ncfile, &ndims, &nvars, &ngatts, &unlimited );
1236  ERR_PARNC( ret );
1237  // find dimension ids for n_S
1238  int ins;
1239  ret = ncmpi_inq_dimid( ncfile, "n_s", &ins );
1240  ERR_PARNC( ret );
1241  MPI_Offset leng;
1242  ret = ncmpi_inq_dimlen( ncfile, ins, &leng );
1243  ERR_PARNC( ret );
1244  nS = (int)leng;
1245  ret = ncmpi_inq_dimid( ncfile, "n_b", &ins );
1246  ERR_PARNC( ret );
1247  ret = ncmpi_inq_dimlen( ncfile, ins, &leng );
1248  ERR_PARNC( ret );
1249  nB = (int)leng;
1250 #else
1251  _EXCEPTION1( "cannot read the file %s", strSource );
1252 #endif
1253  }
1254 
1255  // Let us declare the map object for every process
1256  SparseMatrix< double >& sparseMatrix = this->GetSparseMatrix();
1257 
1258  int localSize = nS / size;
1259  long offsetRead = rank * localSize;
1260  // leftovers on last rank
1261  if( rank == size - 1 )
1262  {
1263  localSize += nS % size;
1264  }
1265 
1266  std::vector< int > vecRow, vecCol;
1267  std::vector< double > vecS;
1268  vecRow.resize( localSize );
1269  vecCol.resize( localSize );
1270  vecS.resize( localSize );
1271 
1272  if( ncMap.is_valid() )
1273  {
1274  varRow->set_cur( (long)( offsetRead ) );
1275  varRow->get( &( vecRow[0] ), localSize );
1276 
1277  varCol->set_cur( (long)( offsetRead ) );
1278  varCol->get( &( vecCol[0] ), localSize );
1279 
1280  varS->set_cur( (long)( offsetRead ) );
1281  varS->get( &( vecS[0] ), localSize );
1282 
1283  ncMap.close();
1284  }
1285  else
1286  {
1287 #ifdef MOAB_HAVE_PNETCDF
1288  // fill the local vectors with the variables from pnetcdf file; first inquire, then fill
1289  MPI_Offset start = (MPI_Offset)offsetRead;
1290  MPI_Offset count = (MPI_Offset)localSize;
1291  int varid;
1292  ret = ncmpi_inq_varid( ncfile, "S", &varid );
1293  ERR_PARNC( ret );
1294  ret = ncmpi_get_vara_double_all( ncfile, varid, &start, &count, &vecS[0] );
1295  ERR_PARNC( ret );
1296  ret = ncmpi_inq_varid( ncfile, "row", &varid );
1297  ERR_PARNC( ret );
1298  ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecRow[0] );
1299  ERR_PARNC( ret );
1300  ret = ncmpi_inq_varid( ncfile, "col", &varid );
1301  ERR_PARNC( ret );
1302  ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecCol[0] );
1303  ERR_PARNC( ret );
1304  ret = ncmpi_close( ncfile );
1305  ERR_PARNC( ret );
1306 #endif
1307  }
1308 
1309  // Now let us set the necessary global-to-local ID maps so that A*x operations
1310  // can be performed cleanly as if map was computed online
1311  row_dtoc_dofmap.clear();
1312  // row_dtoc_dofmap.reserve( nB / size );
1313  col_dtoc_dofmap.clear();
1314  rowMap.clear();
1315  colMap.clear();
1316  // col_dtoc_dofmap.reserve( 2 * nA / size );
1317  // row_dtoc_dofmap.resize( m_nTotDofs_Dest, UINT_MAX );
1318  // col_dtoc_dofmap.resize( m_nTotDofs_SrcCov, UINT_MAX );
1319 
1320 #ifdef MOAB_HAVE_MPI
1321  // bother with tuple list only if size > 1
1322  // otherwise, just fill the sparse matrix
1323  if( size > 1 )
1324  {
1325  std::vector< int > ownership;
1326  // the default trivial partitioning scheme
1327  int nDofs = nB; // this is for row partitioning
1328  if( !row_partition ) nDofs = nA; // column partitioning
1329 
1330  // assert(row_major_ownership == true); // this block is valid only for row-based partitioning
1331  ownership.resize( size );
1332  int nPerPart = nDofs / size;
1333  int nRemainder = nDofs % size; // Keep the remainder in root
1334  ownership[0] = nPerPart + nRemainder;
1335  for( int ip = 1, roffset = ownership[0]; ip < size; ++ip )
1336  {
1337  roffset += nPerPart;
1338  ownership[ip] = roffset;
1339  }
1340  moab::TupleList* tl = new moab::TupleList;
1341  unsigned numr = 1; //
1342  tl->initialize( 3, 0, 0, numr, localSize ); // to proc, row, col, value
1343  tl->enableWriteAccess();
1344  // populate
1345  for( int i = 0; i < localSize; i++ )
1346  {
1347  int rowval = vecRow[i] - 1; // dofs are 1 based in the file
1348  int colval = vecCol[i] - 1;
1349  int to_proc = -1;
1350  int dof_val = colval;
1351  if( row_partition ) dof_val = rowval;
1352 
1353  if( ownership[0] > dof_val )
1354  to_proc = 0;
1355  else
1356  {
1357  for( int ip = 1; ip < size; ++ip )
1358  {
1359  if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
1360  {
1361  to_proc = ip;
1362  break;
1363  }
1364  }
1365  }
1366 
1367  int n = tl->get_n();
1368  tl->vi_wr[3 * n] = to_proc;
1369  tl->vi_wr[3 * n + 1] = rowval;
1370  tl->vi_wr[3 * n + 2] = colval;
1371  tl->vr_wr[n] = vecS[i];
1372  tl->inc_n();
1373  }
1374  // heavy communication
1375  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl, 0 );
1376 
1377  if( owned_dof_ids.size() > 0 )
1378  {
1379  // we need to send desired dof to the rendez_vous point
1380  moab::TupleList tl_re; //
1381  tl_re.initialize( 2, 0, 0, 0, owned_dof_ids.size() ); // to proc, value
1382  tl_re.enableWriteAccess();
1383  // send first to rendez_vous point, decided by trivial partitioning
1384 
1385  for( size_t i = 0; i < owned_dof_ids.size(); i++ )
1386  {
1387  int to_proc = -1;
1388  int dof_val = owned_dof_ids[i] - 1; // dofs are 1 based in the file, partition from 0 ?
1389 
1390  if( ownership[0] > dof_val )
1391  to_proc = 0;
1392  else
1393  {
1394  for( int ip = 1; ip < size; ++ip )
1395  {
1396  if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
1397  {
1398  to_proc = ip;
1399  break;
1400  }
1401  }
1402  }
1403 
1404  int n = tl_re.get_n();
1405  tl_re.vi_wr[2 * n] = to_proc;
1406  tl_re.vi_wr[2 * n + 1] = dof_val;
1407 
1408  tl_re.inc_n();
1409  }
1410  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl_re, 0 );
1411  // now we know in tl_re where do we need to send back dof_val
1412  moab::TupleList::buffer sort_buffer;
1413  sort_buffer.buffer_init( tl_re.get_n() );
1414  tl_re.sort( 1, &sort_buffer ); // so now we order by value
1415 
1416  sort_buffer.buffer_init( tl->get_n() );
1417  int indexOrder = 2; // colVal
1418  if( row_partition ) indexOrder = 1; // rowVal
1419  //tl->sort( indexOrder, &sort_buffer );
1420 
1421  std::map< int, int > startDofIndex, endDofIndex; // indices in tl_re for values we want
1422  int dofVal = -1;
1423  if( tl_re.get_n() > 0 ) dofVal = tl_re.vi_rd[1]; // first dof val on this rank
1424  startDofIndex[dofVal] = 0;
1425  endDofIndex[dofVal] = 0; // start and end
1426  for( unsigned k = 1; k < tl_re.get_n(); k++ )
1427  {
1428  int newDof = tl_re.vi_rd[2 * k + 1];
1429  if( dofVal == newDof )
1430  {
1431  endDofIndex[dofVal] = k; // increment by 1 actually
1432  }
1433  else
1434  {
1435  dofVal = newDof;
1436  startDofIndex[dofVal] = k;
1437  endDofIndex[dofVal] = k;
1438  }
1439  }
1440 
1441  // basically, for each value we are interested in, index in tl_re with those values are
1442  // tl_re.vi_rd[2*startDofIndex+1] == valDof == tl_re.vi_rd[2*endDofIndex+1]
1443  // so now we have ordered
1444  // tl_re shows to what proc do we need to send the tuple (row, col, val)
1445  moab::TupleList* tl_back = new moab::TupleList;
1446  unsigned numr = 1; //
1447  // localSize is a good guess, but maybe it should be bigger ?
1448  // this could be bigger for repeated dofs
1449  tl_back->initialize( 3, 0, 0, numr, tl->get_n() ); // to proc, row, col, value
1450  tl_back->enableWriteAccess();
1451  // now loop over tl and tl_re to see where to send
1452  // form the new tuple, which will contain the desired dofs per task, per row or column distribution
1453 
1454  for( unsigned k = 0; k < tl->get_n(); k++ )
1455  {
1456  int valDof = tl->vi_rd[3 * k + indexOrder]; // 1 for row, 2 for column // first value, it should be
1457  for( int ire = startDofIndex[valDof]; ire <= endDofIndex[valDof]; ire++ )
1458  {
1459  int to_proc = tl_re.vi_rd[2 * ire];
1460  int n = tl_back->get_n();
1461  tl_back->vi_wr[3 * n] = to_proc;
1462  tl_back->vi_wr[3 * n + 1] = tl->vi_rd[3 * k + 1]; // row
1463  tl_back->vi_wr[3 * n + 2] = tl->vi_rd[3 * k + 2]; // col
1464  tl_back->vr_wr[n] = tl->vr_rd[k];
1465  tl_back->inc_n();
1466  }
1467  }
1468 
1469  // now communicate to the desired tasks:
1470  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl_back, 0 );
1471 
1472  tl_re.reset(); // clear memory, although this will go out of scope
1473  tl->reset();
1474  tl = tl_back;
1475  }
1476 
1477  int rindexMax = 0, cindexMax = 0;
1478  // populate the sparsematrix, using rowMap and colMap
1479  int n = tl->get_n();
1480  for( int i = 0; i < n; i++ )
1481  {
1482  int rindex, cindex;
1483  const int& vecRowValue = tl->vi_wr[3 * i + 1];
1484  const int& vecColValue = tl->vi_wr[3 * i + 2];
1485 
1486  std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
1487  if( riter == rowMap.end() )
1488  {
1489  rowMap[vecRowValue] = rindexMax;
1490  rindex = rindexMax;
1491  row_gdofmap.push_back( vecRowValue );
1492  // row_dtoc_dofmap.push_back( vecRowValue );
1493  rindexMax++;
1494  }
1495  else
1496  rindex = riter->second;
1497 
1498  std::map< int, int >::iterator citer = colMap.find( vecColValue );
1499  if( citer == colMap.end() )
1500  {
1501  colMap[vecColValue] = cindexMax;
1502  cindex = cindexMax;
1503  col_gdofmap.push_back( vecColValue );
1504  // col_dtoc_dofmap.push_back( vecColValue );
1505  cindexMax++;
1506  }
1507  else
1508  cindex = citer->second;
1509 
1510  sparseMatrix( rindex, cindex ) = tl->vr_wr[i];
1511  }
1512  tl->reset();
1513  }
1514  else
1515 #endif
1516  {
1517  int rindexMax = 0, cindexMax = 0;
1518 
1519  for( int i = 0; i < nS; i++ )
1520  {
1521  int rindex, cindex;
1522  const int& vecRowValue = vecRow[i] - 1; // the rows, cols are 1 based in the file
1523  const int& vecColValue = vecCol[i] - 1;
1524 
1525  std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
1526  if( riter == rowMap.end() )
1527  {
1528  rowMap[vecRowValue] = rindexMax;
1529  rindex = rindexMax;
1530  row_gdofmap.push_back( vecRowValue );
1531  // row_dtoc_dofmap.push_back( vecRowValue );
1532  rindexMax++;
1533  }
1534  else
1535  rindex = riter->second;
1536 
1537  std::map< int, int >::iterator citer = colMap.find( vecColValue );
1538  if( citer == colMap.end() )
1539  {
1540  colMap[vecColValue] = cindexMax;
1541  cindex = cindexMax;
1542  col_gdofmap.push_back( vecColValue );
1543  // col_dtoc_dofmap.push_back( vecColValue );
1544  cindexMax++;
1545  }
1546  else
1547  cindex = citer->second;
1548 
1549  sparseMatrix( rindex, cindex ) = vecS[i];
1550  }
1551  }
1552 
1553  m_nTotDofs_SrcCov = sparseMatrix.GetColumns();
1554  m_nTotDofs_Dest = sparseMatrix.GetRows();
1555 
1556 #ifdef MOAB_HAVE_EIGEN3
1557  this->copy_tempest_sparsemat_to_eigen3();
1558 #endif
1559 
1560  // Reset the source and target data first
1561  m_rowVector.setZero();
1562  m_colVector.setZero();
1563 
1564  return moab::MB_SUCCESS;
1565 }
1566 
1567 ///////////////////////////////////////////////////////////////////////////////