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