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