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