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