Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
TempestOnlineMap.cpp
Go to the documentation of this file.
1 /*
2  * =====================================================================================
3  *
4  * Filename: TempestOnlineMap.hpp
5  *
6  * Description: Interface to the TempestRemap library to compute the consistent,
7  * and accurate high-order conservative remapping weights for overlap
8  * grids on the sphere in climate simulations.
9  *
10  * Author: Vijay S. Mahadevan (vijaysm), [email protected]
11  *
12  * =====================================================================================
13  */
14 
15 #include "Announce.h"
16 #include "DataArray3D.h"
17 #include "FiniteVolumeTools.h"
18 #include "FiniteElementTools.h"
19 #include "TriangularQuadrature.h"
20 #include "GaussQuadrature.h"
21 #include "GaussLobattoQuadrature.h"
22 #include "SparseMatrix.h"
23 #include "STLStringHelper.h"
24 #include "LinearRemapFV.h"
25 
26 #include "LinearRemapSE0.h"
27 #include "LinearRemapFV.h"
28 
29 #include "moab/Remapping/TempestOnlineMap.hpp"
30 #include "moab/Remapping/IntegerReprosum.hpp"
31 #include "moab/IntxMesh/IntxUtils.hpp"
32 #include "DebugOutput.hpp"
33 #include "moab/TupleList.hpp"
34 #include "moab/MeshTopoUtil.hpp"
35 
36 #include <fstream>
37 #include <cmath>
38 #include <cstdlib>
39 #include <numeric>
40 #include <algorithm>
41 
42 /*
43 #ifdef MOAB_HAVE_NETCDF
44 #ifdef MOAB_HAVE_NETCDFPAR
45 #include "netcdfcpp_par.hpp"
46 #else
47 #include "netcdfcpp.h"
48 #endif
49 #endif
50 */
51 
52 // #define USE_NATIVE_TEMPESTREMAP_ROUTINES
53 
54 ///////////////////////////////////////////////////////////////////////////////
55 
56 // #define VERBOSE
57 // #define VVERBOSE
58 // #define CHECK_INCREASING_DOF
59 
60 ///////////////////////////////////////////////////////////////////////////////
61 
62 #define MPI_CHK_ERR( err ) \
63  if( err ) \
64  { \
65  std::cout << "MPI Failure. ErrorCode (" << ( err ) << ") "; \
66  std::cout << "\nMPI Aborting... \n"; \
67  return moab::MB_FAILURE; \
68  }
69 
70 moab::TempestOnlineMap::TempestOnlineMap( moab::TempestRemapper* remapper ) : OfflineMap(), m_remapper( remapper )
71 {
72  // Get the references for the MOAB core objects
73  m_interface = m_remapper->get_interface();
74 #ifdef MOAB_HAVE_MPI
75  m_pcomm = m_remapper->get_parallel_communicator();
76 #endif
77 
78  // now let us re-update the reference to the input source mesh
79  m_meshInput = m_remapper->GetMesh( moab::Remapper::SourceMesh );
80  // now let us re-update the reference to the covering mesh
81  m_meshInputCov = m_remapper->GetCoveringMesh();
82  // now let us re-update the reference to the output target mesh
83  m_meshOutput = m_remapper->GetMesh( moab::Remapper::TargetMesh );
84  // now let us re-update the reference to the output target mesh
85  m_meshOverlap = m_remapper->GetMesh( moab::Remapper::OverlapMesh );
86 
87  is_parallel = remapper->is_parallel;
88  is_root = remapper->is_root;
89  rank = remapper->rank;
90  size = remapper->size;
91 
92  // set default order
93  m_input_order = m_output_order = 1;
94 
95  // unknown until a map file is read (ReadParallelMap sets it to n_a)
96  m_nTotDofs_SrcGlobal = -1;
97 
98  // Initialize dimension information from file
99  this->setup_sizes_dimensions();
100 }
101 
102 void moab::TempestOnlineMap::setup_sizes_dimensions()
103 {
104  if( m_meshInputCov )
105  {
106  std::vector< std::string > dimNames;
107  std::vector< int > dimSizes;
108  dimNames.push_back( "num_elem" );
109  dimSizes.push_back( m_meshInputCov->faces.size() );
110 
111  this->InitializeSourceDimensions( dimNames, dimSizes );
112  }
113 
114  if( m_meshOutput )
115  {
116  std::vector< std::string > dimNames;
117  std::vector< int > dimSizes;
118  dimNames.push_back( "num_elem" );
119  dimSizes.push_back( m_meshOutput->faces.size() );
120 
121  this->InitializeTargetDimensions( dimNames, dimSizes );
122  }
123 }
124 
125 ///////////////////////////////////////////////////////////////////////////////
126 
127 moab::TempestOnlineMap::~TempestOnlineMap()
128 {
129  m_interface = nullptr;
130 #ifdef MOAB_HAVE_MPI
131  m_pcomm = nullptr;
132 #endif
133  m_meshInput = nullptr;
134  m_meshOutput = nullptr;
135  m_meshOverlap = nullptr;
136 }
137 
138 ///////////////////////////////////////////////////////////////////////////////
139 
140 moab::ErrorCode moab::TempestOnlineMap::SetDOFmapTags( const std::string srcDofTagName,
141  const std::string tgtDofTagName )
142 {
143  moab::ErrorCode rval;
144 
145  int tagSize = 0;
146  tagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
147  rval =
148  m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );
149 
150  if( rval == moab::MB_TAG_NOT_FOUND && m_eInputType != DiscretizationType_FV )
151  {
152  MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for source mesh." );
153  }
154  else
155  MB_CHK_ERR( rval );
156 
157  tagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
158  rval =
159  m_interface->tag_get_handle( tgtDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagDest, MB_TAG_ANY );
160  if( rval == moab::MB_TAG_NOT_FOUND && m_eOutputType != DiscretizationType_FV )
161  {
162  MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for target mesh." );
163  }
164  else
165  MB_CHK_ERR( rval );
166 
167  return moab::MB_SUCCESS;
168 }
169 
170 ///////////////////////////////////////////////////////////////////////////////
171 
172 moab::ErrorCode moab::TempestOnlineMap::SetDOFmapAssociation( DiscretizationType srcType,
173  int srcOrder,
174  bool isSrcContinuous,
175  DataArray3D< int >* srcdataGLLNodes,
176  DataArray3D< int >* srcdataGLLNodesSrc,
177  DiscretizationType destType,
178  int destOrder,
179  bool isTgtContinuous,
180  DataArray3D< int >* tgtdataGLLNodes )
181 {
182  std::vector< bool > dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
183  std::vector< int > src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;
184 
185  // We are assuming that these are element based tags that are sized: np * np
186  m_srcDiscType = srcType;
187  m_destDiscType = destType;
188  m_input_order = srcOrder;
189  m_output_order = destOrder;
190 
191  bool vprint = is_root && false;
192 
193  // Compute and store the total number of source and target DoFs corresponding
194  // to number of rows and columns in the mapping.
195  // Now compute the mapping and store it for the covering mesh
196  int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
197  if( m_remapper->point_cloud_source )
198  {
199  assert( m_nDofsPEl_Src == 1 );
200  col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
201  col_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), -1 );
202  src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), -1 );
203  MB_CHK_ERR(
204  m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] ) );
205  srcTagSize = 1;
206  }
207  else
208  {
209  col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
210  col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, -1 );
211  src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, -1 );
212  MB_CHK_ERR(
213  m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] ) );
214  }
215 
216  m_nTotDofs_SrcCov = 0;
217  if( srcdataGLLNodes == nullptr )
218  {
219  /* we only have a mapping for elements as DoFs */
220  for( unsigned i = 0; i < col_gdofmap.size(); ++i )
221  {
222  auto gdof = src_soln_gdofs[i];
223  assert( gdof > 0 );
224  col_gdofmap[i] = gdof - 1;
225  col_dtoc_dofmap[i] = i;
226  if( vprint ) std::cout << "Col: " << i << ", " << col_gdofmap[i] << "\n";
227  m_nTotDofs_SrcCov++;
228  }
229  }
230  else
231  {
232  if( isSrcContinuous )
233  dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, false );
234  // Put these remap coefficients into the SparseMatrix map
235  for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
236  {
237  for( int p = 0; p < m_nDofsPEl_Src; p++ )
238  {
239  for( int q = 0; q < m_nDofsPEl_Src; q++ )
240  {
241  const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
242  const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
243  if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
244  {
245  m_nTotDofs_SrcCov++;
246  dgll_cgll_covcol_ldofmap[localDOF] = true;
247  }
248  if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
249  assert( src_soln_gdofs[offsetDOF] > 0 );
250  // For CGLL: weight matrix uses localDOF (continuous shared node index)
251  // as column index → col_gdofmap must be indexed by localDOF
252  // For DGLL: weight matrix uses offsetDOF (= elem*nP*nP + p*nP + q)
253  // as column index → col_gdofmap must be indexed by offsetDOF
254  if( isSrcContinuous )
255  {
256  col_gdofmap[localDOF] = src_soln_gdofs[offsetDOF] - 1;
257  col_dtoc_dofmap[offsetDOF] = localDOF;
258  }
259  else
260  {
261  col_gdofmap[offsetDOF] = src_soln_gdofs[offsetDOF] - 1;
262  col_dtoc_dofmap[offsetDOF] = offsetDOF;
263  }
264  }
265  }
266  }
267  }
268 
269  if( m_remapper->point_cloud_source )
270  {
271  assert( m_nDofsPEl_Src == 1 );
272  srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
273  srccol_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), -1 );
274  locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), -1 );
275  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] ) );
276  }
277  else
278  {
279  srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
280  srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, -1 );
281  locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, -1 );
282  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] ) );
283  }
284 
285  // Now compute the mapping and store it for the original source mesh
286  m_nTotDofs_Src = 0;
287  if( srcdataGLLNodesSrc == nullptr )
288  {
289  /* we only have a mapping for elements as DoFs */
290  for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
291  {
292  auto gdof = locsrc_soln_gdofs[i];
293  assert( gdof > 0 );
294  srccol_gdofmap[i] = gdof - 1;
295  srccol_dtoc_dofmap[i] = i;
296  m_nTotDofs_Src++;
297  }
298  }
299  else
300  {
301  if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, false );
302  // Put these remap coefficients into the SparseMatrix map
303  for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
304  {
305  for( int p = 0; p < m_nDofsPEl_Src; p++ )
306  {
307  for( int q = 0; q < m_nDofsPEl_Src; q++ )
308  {
309  const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
310  const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
311  if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
312  {
313  m_nTotDofs_Src++;
314  dgll_cgll_col_ldofmap[localDOF] = true;
315  }
316  if( !isSrcContinuous ) m_nTotDofs_Src++;
317  assert( locsrc_soln_gdofs[offsetDOF] > 0 );
318  if( isSrcContinuous )
319  {
320  srccol_gdofmap[localDOF] = locsrc_soln_gdofs[offsetDOF] - 1;
321  srccol_dtoc_dofmap[offsetDOF] = localDOF;
322  }
323  else
324  {
325  srccol_gdofmap[offsetDOF] = locsrc_soln_gdofs[offsetDOF] - 1;
326  srccol_dtoc_dofmap[offsetDOF] = offsetDOF;
327  }
328  }
329  }
330  }
331  }
332 
333  int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
334  if( m_remapper->point_cloud_target )
335  {
336  assert( m_nDofsPEl_Dest == 1 );
337  row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
338  row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), -1 );
339  tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), -1 );
340  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] ) );
341  tgtTagSize = 1;
342  }
343  else
344  {
345  row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
346  row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, -1 );
347  tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, -1 );
348  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] ) );
349  }
350 
351  // Now compute the mapping and store it for the target mesh
352  // To access the GID for each row: row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
353  m_nTotDofs_Dest = 0;
354  if( tgtdataGLLNodes == nullptr )
355  {
356  /* we only have a mapping for elements as DoFs */
357  for( unsigned i = 0; i < row_gdofmap.size(); ++i )
358  {
359  auto gdof = tgt_soln_gdofs[i];
360  assert( gdof > 0 );
361  row_gdofmap[i] = gdof - 1;
362  row_dtoc_dofmap[i] = i;
363  if( vprint ) std::cout << "Row: " << i << ", " << row_gdofmap[i] << "\n";
364  m_nTotDofs_Dest++;
365  }
366  }
367  else
368  {
369  if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, false );
370  // Put these remap coefficients into the SparseMatrix map
371  for( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
372  {
373  for( int p = 0; p < m_nDofsPEl_Dest; p++ )
374  {
375  for( int q = 0; q < m_nDofsPEl_Dest; q++ )
376  {
377  const int localDOF = ( *tgtdataGLLNodes )[p][q][j] - 1;
378  const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
379  if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
380  {
381  m_nTotDofs_Dest++;
382  dgll_cgll_row_ldofmap[localDOF] = true;
383  }
384  if( !isTgtContinuous ) m_nTotDofs_Dest++;
385  assert( tgt_soln_gdofs[offsetDOF] > 0 );
386  if( isTgtContinuous )
387  {
388  row_gdofmap[localDOF] = tgt_soln_gdofs[offsetDOF] - 1;
389  row_dtoc_dofmap[offsetDOF] = localDOF;
390  }
391  else
392  {
393  row_gdofmap[offsetDOF] = tgt_soln_gdofs[offsetDOF] - 1;
394  row_dtoc_dofmap[offsetDOF] = offsetDOF;
395  }
396  if( vprint )
397  std::cout << "Row: " << offsetDOF << ", " << localDOF << ", " << row_gdofmap[offsetDOF] << ", "
398  << m_nTotDofs_Dest << "\n";
399  }
400  }
401  }
402  }
403 
404  // Let us also allocate the local representation of the sparse matrix
405 #if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
406  if( is_root )
407  {
408  std::cout << "[" << rank << "] DoFs: row = " << m_nTotDofs_Dest << " (gdofmap.size=" << row_gdofmap.size()
409  << "), col_src = " << m_nTotDofs_Src << ", col_cov = " << m_nTotDofs_SrcCov
410  << " (gdofmap.size=" << col_gdofmap.size() << ")\n";
411  }
412 #endif
413 
414  // check monotonicity of row_gdofmap and col_gdofmap
415 #ifdef CHECK_INCREASING_DOF
416  for( size_t i = 0; i < row_gdofmap.size() - 1; i++ )
417  {
418  if( row_gdofmap[i] > row_gdofmap[i + 1] )
419  std::cout << " on rank " << rank << " in row_gdofmap[" << i << "]=" << row_gdofmap[i] << " > row_gdofmap["
420  << i + 1 << "]=" << row_gdofmap[i + 1] << " \n";
421  }
422  for( size_t i = 0; i < col_gdofmap.size() - 1; i++ )
423  {
424  if( col_gdofmap[i] > col_gdofmap[i + 1] )
425  std::cout << " on rank " << rank << " in col_gdofmap[" << i << "]=" << col_gdofmap[i] << " > col_gdofmap["
426  << i + 1 << "]=" << col_gdofmap[i + 1] << " \n";
427  }
428 #endif
429 
430  return moab::MB_SUCCESS;
431 }
432 
433 moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs( std::vector< int >& values_entities )
434 {
435  // col_gdofmap has global dofs , that should be in the list of values, such that
436  // row_dtoc_dofmap[offsetDOF] = localDOF;
437  // we need to find col_dtoc_dofmap such that: col_gdofmap[ col_dtoc_dofmap[i] ] == values_entities [i];
438  // we know that col_gdofmap[0..(nbcols-1)] = global_col_dofs -> in values_entities
439  // form first inverse
440  //
441  // resize and initialize to -1 to signal that this value should not be used, if not set below
442  col_dtoc_dofmap.resize( values_entities.size(), -1 );
443  for( size_t j = 0; j < values_entities.size(); j++ )
444  {
445  // values are 1 based, but rowMap, colMap are not
446  const auto it = colMap.find( values_entities[j] - 1 );
447  if( it != colMap.end() ) col_dtoc_dofmap[j] = it->second;
448  }
449  return moab::MB_SUCCESS;
450 }
451 
452 moab::ErrorCode moab::TempestOnlineMap::set_row_dc_dofs( std::vector< int >& values_entities )
453 {
454  // we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
455  // resize and initialize to -1 to signal that this value should not be used, if not set below
456  row_dtoc_dofmap.resize( values_entities.size(), -1 );
457  for( size_t j = 0; j < values_entities.size(); j++ )
458  {
459  // values are 1 based, but rowMap, colMap are not
460  const auto it = rowMap.find( values_entities[j] - 1 );
461  if( it != rowMap.end() ) row_dtoc_dofmap[j] = it->second;
462  }
463  return moab::MB_SUCCESS;
464 }
465 
466 // Compute which weight-matrix columns the migrated coverage covers.
467 // delivered[mc] == true iff some covering cell maps to matrix column mc.
468 static void compute_delivered_columns( int ncols, const std::vector< int >& col_dtoc_dofmap,
469  std::vector< bool >& delivered )
470 {
471  delivered.assign( ncols, false );
472  for( size_t k = 0; k < col_dtoc_dofmap.size(); k++ )
473  {
474  const int mc = col_dtoc_dofmap[k];
475  if( mc >= 0 && mc < ncols ) delivered[mc] = true;
476  }
477 }
478 
479 int moab::TempestOnlineMap::CountAbsentColumns( int& first_absent_gid ) const
480 {
481  first_absent_gid = -1;
482  const int ncols = m_nTotDofs_SrcCov;
483  std::vector< bool > delivered;
484  compute_delivered_columns( ncols, col_dtoc_dofmap, delivered );
485  int cnt = 0;
486  for( int mc = 0; mc < ncols; mc++ )
487  {
488  if( !delivered[mc] )
489  {
490  cnt++;
491  if( first_absent_gid < 0 && mc < (int)col_gdofmap.size() )
492  first_absent_gid = (int)col_gdofmap[mc] + 1; // col_gdofmap is 0-based
493  }
494  }
495  return cnt;
496 }
497 
498 int moab::TempestOnlineMap::DropAbsentColumns()
499 {
500  const int ncols = m_nTotDofs_SrcCov;
501  std::vector< bool > delivered;
502  compute_delivered_columns( ncols, col_dtoc_dofmap, delivered );
503 
504  // Zero every stored coefficient whose column was not supplied by coverage.
505  // The projection already treats these columns as zero-source (ApplyWeights
506  // leaves m_colVector at 0 for them), so this changes no projected value; it
507  // only lets the dual-map CAAS bounds loop skip them via its |w|<1e-50 test.
508  int dropped = 0;
509  for( int r = 0; r < m_weightMatrix.outerSize(); r++ )
510  {
511  for( WeightMatrix::InnerIterator it( m_weightMatrix, r ); it; ++it )
512  {
513  const int mc = (int)it.col();
514  if( mc < 0 || mc >= ncols || !delivered[mc] )
515  {
516  if( it.value() != 0.0 ) dropped++;
517  it.valueRef() = 0.0;
518  }
519  }
520  }
521  // Remove the explicit zeros so iterators no longer visit them.
522  m_weightMatrix.prune( []( const Eigen::Index&, const Eigen::Index&, const double& v ) { return v != 0.0; } );
523  return dropped;
524 }
525 ///////////////////////////////////////////////////////////////////////////////
526 
527 moab::ErrorCode moab::TempestOnlineMap::GenerateRemappingWeights( std::string strInputType,
528  std::string strOutputType,
529  const GenerateOfflineMapAlgorithmOptions& mapOptions,
530  const std::string& srcDofTagName,
531  const std::string& tgtDofTagName )
532 {
533 #ifdef MOAB_HAVE_NETCDF
534  NcError error( NcError::silent_nonfatal );
535 #endif
536 
537  moab::DebugOutput dbgprint( std::cout, rank, 0 );
538  dbgprint.set_prefix( "[TempestOnlineMap]: " );
539  moab::ErrorCode rval;
540 
541  const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
542  const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
543  const bool m_bPointCloud = m_bPointCloudSource || m_bPointCloudTarget;
544 
545  // Build a matrix of source and target discretization so that we know how
546  // to assign the global DoFs in parallel for the mapping weights.
547  // For example,
548  // for FV->FV: the rows represented target DoFs and cols represent source DoFs
549  try
550  {
551  // Check command line parameters (data type arguments)
552  STLStringHelper::ToLower( strInputType );
553  STLStringHelper::ToLower( strOutputType );
554 
555  DiscretizationType eInputType;
556  DiscretizationType eOutputType;
557 
558  if( strInputType == "fv" )
559  {
560  eInputType = DiscretizationType_FV;
561  }
562  else if( strInputType == "cgll" )
563  {
564  eInputType = DiscretizationType_CGLL;
565  }
566  else if( strInputType == "dgll" )
567  {
568  eInputType = DiscretizationType_DGLL;
569  }
570  else if( strInputType == "pcloud" )
571  {
572  eInputType = DiscretizationType_PCLOUD;
573  }
574  else
575  {
576  _EXCEPTION1( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", strInputType.c_str() );
577  }
578 
579  if( strOutputType == "fv" )
580  {
581  eOutputType = DiscretizationType_FV;
582  }
583  else if( strOutputType == "cgll" )
584  {
585  eOutputType = DiscretizationType_CGLL;
586  }
587  else if( strOutputType == "dgll" )
588  {
589  eOutputType = DiscretizationType_DGLL;
590  }
591  else if( strOutputType == "pcloud" )
592  {
593  eOutputType = DiscretizationType_PCLOUD;
594  }
595  else
596  {
597  _EXCEPTION1( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", strOutputType.c_str() );
598  }
599 
600  // set all required input params
601  m_bConserved = !mapOptions.fNoConservation;
602  m_eInputType = eInputType;
603  m_eOutputType = eOutputType;
604 
605  // Method flags
606  std::string strMapAlgorithm( "" );
607  int nMonotoneType = ( mapOptions.fMonotone ) ? ( 1 ) : ( 0 );
608 
609  // Make an index of method arguments
610  std::set< std::string > setMethodStrings;
611  {
612  int iLast = 0;
613  for( size_t i = 0; i <= mapOptions.strMethod.length(); i++ )
614  {
615  if( ( i == mapOptions.strMethod.length() ) || ( mapOptions.strMethod[i] == ';' ) )
616  {
617  std::string strMethodString = mapOptions.strMethod.substr( iLast, i - iLast );
618  STLStringHelper::RemoveWhitespaceInPlace( strMethodString );
619  if( strMethodString.length() > 0 )
620  {
621  setMethodStrings.insert( strMethodString );
622  }
623  iLast = i + 1;
624  }
625  }
626  }
627 
628  for( const auto& it : setMethodStrings )
629  {
630  // Piecewise constant monotonicity
631  if( it == "mono2" )
632  {
633  if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
634  {
635  _EXCEPTIONT( "--method \"mono2\" is only used when remapping to/from CGLL or DGLL grids" );
636  }
637  nMonotoneType = 2;
638 
639  // Piecewise linear monotonicity
640  }
641  else if( it == "mono3" )
642  {
643  if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
644  {
645  _EXCEPTIONT( "--method \"mono3\" is only used when remapping to/from CGLL or DGLL grids" );
646  }
647  nMonotoneType = 3;
648 
649  // Volumetric remapping from FV to GLL
650  }
651  else if( it == "volumetric" )
652  {
653  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType == DiscretizationType_FV ) )
654  {
655  _EXCEPTIONT( "--method \"volumetric\" may only be used for FV->CGLL or FV->DGLL remapping" );
656  }
657  strMapAlgorithm = "volumetric";
658 
659  // Inverse distance mapping
660  }
661  else if( it == "invdist" )
662  {
663  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
664  {
665  _EXCEPTIONT( "--method \"invdist\" may only be used for FV->FV remapping" );
666  }
667  strMapAlgorithm = "invdist";
668 
669  // Delaunay triangulation mapping
670  }
671  else if( it == "delaunay" )
672  {
673  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
674  {
675  _EXCEPTIONT( "--method \"delaunay\" may only be used for FV->FV remapping" );
676  }
677  strMapAlgorithm = "delaunay";
678 
679  // Bilinear
680  }
681  else if( it == "bilin" )
682  {
683  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
684  {
685  _EXCEPTIONT( "--method \"bilin\" may only be used for FV->FV remapping" );
686  }
687  strMapAlgorithm = "fvbilin";
688 
689  // Integrated bilinear (same as mono3 when source grid is CGLL/DGLL)
690  }
691  else if( it == "intbilin" )
692  {
693  if( m_eOutputType != DiscretizationType_FV )
694  {
695  _EXCEPTIONT( "--method \"intbilin\" may only be used when mapping to FV." );
696  }
697  if( m_eInputType == DiscretizationType_FV )
698  {
699  strMapAlgorithm = "fvintbilin";
700  }
701  else
702  {
703  strMapAlgorithm = "mono3";
704  }
705 
706  // Integrated bilinear with generalized Barycentric coordinates
707  }
708  else if( it == "intbilingb" )
709  {
710  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
711  {
712  _EXCEPTIONT( "--method \"intbilingb\" may only be used for FV->FV remapping" );
713  }
714  strMapAlgorithm = "fvintbilingb";
715  }
716  else
717  {
718  _EXCEPTION1( "Invalid --method argument \"%s\"", it.c_str() );
719  }
720  }
721 
722  m_nDofsPEl_Src =
723  ( m_eInputType == DiscretizationType_FV || m_eInputType == DiscretizationType_PCLOUD ? 1
724  : mapOptions.nPin );
725  m_nDofsPEl_Dest =
726  ( m_eOutputType == DiscretizationType_FV || m_eOutputType == DiscretizationType_PCLOUD ? 1
727  : mapOptions.nPout );
728 
729  // Set the source and target mesh objects
730  MB_CHK_ERR( SetDOFmapTags( srcDofTagName, tgtDofTagName ) );
731 
732  /// the tag should be created already in the e3sm workflow; if not, create it here
733  Tag areaTag;
734  rval = m_interface->tag_get_handle( "aream", 1, MB_TYPE_DOUBLE, areaTag,
735  MB_TAG_DENSE | MB_TAG_EXCL | MB_TAG_CREAT );
736  if( MB_ALREADY_ALLOCATED == rval )
737  {
738  if( is_root ) dbgprint.printf( 0, "aream tag already defined \n" );
739  }
740 
741  double local_areas[3] = { 0.0, 0.0, 0.0 }, global_areas[3] = { 0.0, 0.0, 0.0 };
742  if( !m_bPointCloudSource )
743  {
744  // Calculate Input Mesh Face areas
745  if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
746  local_areas[0] = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
747  // Set source element areas as tag on the source mesh
748  MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea ) );
749 
750  // Update coverage source mesh areas as well.
751  m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
752  }
753 
754  if( !m_bPointCloudTarget )
755  {
756  // Calculate Output Mesh Face areas
757  if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
758  local_areas[1] = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
759  // Set target element areas as tag on the target mesh
760  MB_CHK_ERR(
761  m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea ) );
762  }
763 
764  if( !m_bPointCloud )
765  {
766  // Calculate Face areas
767  if (m_meshOverlap)
768  {
769  // Verify that overlap mesh is in the correct order (sanity check)
770  assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->vecTargetFaceIx.size() );
771 
772  if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
773  local_areas[2] =
774  m_meshOverlap->CalculateFaceAreas( mapOptions.fSourceConcave || mapOptions.fTargetConcave );
775 
776 #ifdef MOAB_HAVE_MPI
777  // CalculateFaceAreas() sums every face it holds; it has no notion of ownership.
778  // In parallel the overlap mesh also carries ghost elements -- intersections whose
779  // target cell is owned by another rank -- which ConvertOverlapMeshSourceOrdered()
780  // flags by setting vecTargetFaceIx to -1 so they are skipped when the weights are
781  // accumulated. Summing the raw total therefore counts those intersections twice
782  // (once here, once on the owning rank) and the reduced "Recovered Area" overshoots
783  // the sphere. Subtract the ghost contribution so the reported area matches the
784  // area the map is actually built from.
785  if( m_pcomm && is_parallel )
786  {
787  double ghost_area = 0.0;
788  for( size_t iover = 0; iover < m_meshOverlap->faces.size(); iover++ )
789  if( m_meshOverlap->vecTargetFaceIx[iover] < 0 ) ghost_area += m_meshOverlap->vecFaceArea[iover];
790  local_areas[2] -= ghost_area;
791  }
792 #endif
793  }
794 
795  // store it as global output for now - used later in reduction
796  std::copy( local_areas, local_areas + 3, global_areas );
797 #ifdef MOAB_HAVE_MPI
798  // reduce the local source, target and overlap mesh areas to global areas
799  if( m_pcomm && is_parallel )
800  MPI_Reduce( local_areas, global_areas, 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
801 #endif
802  if( is_root )
803  {
804  dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", global_areas[0] );
805  dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", global_areas[1] );
806  if (m_meshOverlap) dbgprint.printf( 0, "Overlap Mesh Recovered Area: %1.15e\n", global_areas[2] );
807  }
808 
809  // Correct areas to match the areas calculated in the overlap mesh
810  constexpr bool fCorrectAreas = true;
811  if( fCorrectAreas && m_meshOverlap ) // In MOAB-TempestRemap, we will always keep this to be true
812  {
813  if( is_root ) dbgprint.printf( 0, "Correcting source/target areas to overlap mesh areas\n" );
814  DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() );
815  DataArray1D< double > dTargetArea( m_meshOutput->faces.size() );
816 
817  assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() );
818  assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() );
819  assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() );
820 
821  assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() );
822  assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() );
823 
824  for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
825  {
826  if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 )
827  continue; // skip this cell since it is ghosted
828 
829  // let us recompute the source/target areas based on overlap mesh areas
830  assert( static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() );
831  dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
832  assert( static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() );
833  dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
834  }
835 
836  // How much disagreement between the accumulated overlap area and the
837  // geometric area we are willing to absorb.
838  //
839  // For meshes whose domains coincide, the two agree to roundoff and the
840  // 1e-10 test simply avoids replacing a good geometric area with a
841  // slightly noisier accumulated one.
842  //
843  // For a regionally refined mesh they legitimately differ: a cell on the
844  // edge of the other mesh's domain is only partly covered, so its overlap
845  // contributions sum to a fraction of its geometric area. The 1e-10 test
846  // then rejects the correction precisely where it is needed, leaving the
847  // cell with its full area while it receives only part of the overlap --
848  // which drives its map row sum below one. Accept the accumulated area
849  // for any cell that received a contribution, and leave cells with no
850  // coverage at their geometric area (they are outside the other domain;
851  // zeroing them would divide by zero downstream).
852  const bool regional = ( m_remapper != nullptr && m_remapper->IsRegionalMesh() );
853 
854  auto accept_area = [regional]( double accumulated, double geometric ) -> bool {
855  if( regional ) return accumulated > 0.0;
856  return fabs( accumulated - geometric ) < 1.0e-10;
857  };
858 
859  size_t nUncoveredTarget = 0;
860  for( size_t i = 0; i < m_meshInputCov->faces.size(); i++ )
861  {
862  if( accept_area( dSourceArea[i], m_meshInputCov->vecFaceArea[i] ) )
863  m_meshInputCov->vecFaceArea[i] = dSourceArea[i];
864  }
865  for( size_t i = 0; i < m_meshOutput->faces.size(); i++ )
866  {
867  if( accept_area( dTargetArea[i], m_meshOutput->vecFaceArea[i] ) )
868  m_meshOutput->vecFaceArea[i] = dTargetArea[i];
869  else if( regional && dTargetArea[i] <= 0.0 )
870  nUncoveredTarget++;
871  }
872 
873  if( regional )
874  {
875  // A regional remap legitimately leaves parts of the target uncovered.
876  // Report it: silently emitting empty rows is what made the previous
877  // behaviour so hard to diagnose.
878  //
879  // Only the target count is meaningful. m_meshInputCov is the
880  // *coverage* mesh, so in parallel it also holds cells owned by other
881  // ranks -- these have no overlap here by construction, and their
882  // number grows with the ghost-layer count rather than with the size
883  // of any hole. Reporting that as "uncovered" would be alarming and
884  // wrong, so it is left out.
885  size_t nUncovered = nUncoveredTarget;
886 #ifdef MOAB_HAVE_MPI
887  if( m_pcomm && is_parallel )
888  {
889  size_t local = nUncoveredTarget;
890  MPI_Reduce( &local, &nUncovered, 1, MPI_UNSIGNED_LONG, MPI_SUM, 0, m_pcomm->comm() );
891  }
892 #endif
893  if( is_root && nUncovered > 0 )
894  dbgprint.printf( 0,
895  "Regional mesh: %zu target cells have no overlap coverage; their areas "
896  "are left geometric and their map rows will be empty\n",
897  nUncovered );
898  }
899  }
900 
901  // Set source mesh areas in map
902  if( !m_bPointCloudSource && eInputType == DiscretizationType_FV )
903  {
904  this->SetSourceAreas( m_meshInputCov->vecFaceArea );
905  if( m_meshInputCov->vecMask.size() )
906  {
907  this->SetSourceMask( m_meshInputCov->vecMask );
908  }
909  }
910 
911  // Set target mesh areas in map
912  if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV )
913  {
914  this->SetTargetAreas( m_meshOutput->vecFaceArea );
915  if( m_meshOutput->vecMask.size() )
916  {
917  this->SetTargetMask( m_meshOutput->vecMask );
918  }
919  }
920 
921  // NOTE: Mesh::CalculateFaceAreasFromOverlap() looks like the natural helper
922  // for the sub-area case, but it is the wrong tool here on two counts: it
923  // indexes by vecSourceFaceIx, which refers to the *coverage* mesh rather
924  // than m_meshInput, and it zeroes every area before accumulating, which
925  // would discard the geometric area of cells the overlap never touches. The
926  // per-cell correction above does the same accumulation with the right
927  // indexing and only replaces areas it can justify.
928  }
929 
930  // Finite volume input / Finite volume output
931  if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
932  {
933  // Generate reverse node array and edge map
934  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
935  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
936 
937  // Initialize coordinates for map
938  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
939  this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
940 
941  this->m_pdataGLLNodesIn = nullptr;
942  this->m_pdataGLLNodesOut = nullptr;
943 
944  // Finite volume input / Finite element output
945  MB_CHK_ERR( this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
946  mapOptions.nPout, false, nullptr ) );
947 
948  // Construct remap for FV-FV
949  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
950 
951  // Construct OfflineMap
952  if( strMapAlgorithm == "invdist" )
953  {
954  if( m_meshInputCov->faces.size() )
955  {
956  if( is_root ) dbgprint.printf( 0, "Calculating map (invdist)\n" );
957  LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
958  }
959  }
960  else if( strMapAlgorithm == "delaunay" ) // does not need intersection mesh
961  {
962  if( m_meshInputCov->faces.size() )
963  {
964  if( is_root ) dbgprint.printf( 0, "Calculating map (delaunay)\n" );
965  if (m_meshOverlap) LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
966  else
967  {
968  Mesh dummy;
969  LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, dummy, *this );
970  }
971  }
972  }
973  else if( strMapAlgorithm == "fvintbilin" )
974  {
975  if( m_meshInputCov->faces.size() )
976  {
977  if( is_root ) dbgprint.printf( 0, "Calculating map (intbilin)\n" );
978  LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
979  }
980  }
981  else if( strMapAlgorithm == "fvintbilingb" )
982  {
983  if( m_meshInputCov->faces.size() )
984  {
985  if( is_root ) dbgprint.printf( 0, "Calculating map (intbilingb)\n" );
986  LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
987  *this );
988  }
989  }
990  else if( strMapAlgorithm == "fvbilin" ) // does not need intersection mesh
991  {
992 #ifdef VERBOSE
993  if( is_root )
994  {
995  m_meshInputCov->Write( "SourceMeshMBTR.g" );
996  m_meshOutput->Write( "TargetMeshMBTR.g" );
997  }
998  else
999  {
1000  m_meshInputCov->Write( "SourceMeshMBTR" + std::to_string( rank ) + ".g" );
1001  m_meshOutput->Write( "TargetMeshMBTR" + std::to_string( rank ) + ".g" );
1002  }
1003 #endif
1004 
1005  if( m_meshInputCov->faces.size() )
1006  {
1007  if( is_root ) dbgprint.printf( 0, "Calculating map (bilin)\n" );
1008  if (m_meshOverlap) LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1009  else
1010  {
1011  Mesh dummy;
1012  LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, dummy, *this );
1013  }
1014  }
1015  }
1016  else
1017  {
1018  if( is_root ) dbgprint.printf( 0, "Calculating conservative FV-FV map\n" );
1019  if( m_meshInputCov->faces.size() )
1020  {
1021 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1022  LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
1023  ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this );
1024 #else
1025  LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
1026 #endif
1027  }
1028  }
1029  }
1030  else if( eInputType == DiscretizationType_FV )
1031  {
1032  DataArray3D< double > dataGLLJacobian;
1033 
1034  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1035  double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
1036  dataGLLNodesDest, dataGLLJacobian );
1037 
1038  double dNumericalArea = dNumericalArea_loc;
1039 #ifdef MOAB_HAVE_MPI
1040  if( m_pcomm )
1041  MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
1042 #endif
1043  if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
1044 
1045  // Initialize coordinates for map
1046  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1047  this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1048 
1049  this->m_pdataGLLNodesIn = nullptr;
1050  this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1051 
1052  // Generate the continuous Jacobian
1053  bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
1054 
1055  if( eOutputType == DiscretizationType_CGLL )
1056  {
1057  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
1058  }
1059  else
1060  {
1061  GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
1062  }
1063 
1064  // Generate reverse node array and edge map
1065  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1066  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
1067 
1068  // Finite volume input / Finite element output
1069  MB_CHK_ERR( this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
1070  mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
1071  &dataGLLNodesDest ) );
1072 
1073  // Generate remap weights
1074  if( strMapAlgorithm == "volumetric" )
1075  {
1076  if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL (volumetric)\n" );
1077  LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
1078  dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *this,
1079  nMonotoneType, fContinuous, mapOptions.fNoConservation );
1080  }
1081  else
1082  {
1083  if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL\n" );
1084  LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
1085  this->GetTargetAreas(), mapOptions.nPin, *this, nMonotoneType, fContinuous,
1086  mapOptions.fNoConservation );
1087  }
1088  }
1089  else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) )
1090  {
1091  DataArray3D< double > dataGLLJacobian;
1092  if( !m_bPointCloudSource )
1093  {
1094  // Generate reverse node array and edge map
1095  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1096  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
1097 
1098  // Initialize coordinates for map
1099  if( eInputType == DiscretizationType_FV )
1100  {
1101  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1102  }
1103  else
1104  {
1105  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1106  DataArray3D< double > dataGLLJacobianSrc;
1107  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1108  dataGLLJacobian );
1109  GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1110  dataGLLJacobianSrc );
1111  }
1112  }
1113  // else { /* Source is a point cloud dataset */ }
1114 
1115  if( !m_bPointCloudTarget )
1116  {
1117  // Generate reverse node array and edge map
1118  if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray();
1119  if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap( false );
1120 
1121  // Initialize coordinates for map
1122  if( eOutputType == DiscretizationType_FV )
1123  {
1124  this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput );
1125  }
1126  else
1127  {
1128  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1129  GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
1130  dataGLLJacobian );
1131  }
1132  }
1133  // else { /* Target is a point cloud dataset */ }
1134 
1135  // Finite volume input / Finite element output
1136  MB_CHK_ERR( this->SetDOFmapAssociation(
1137  eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1138  ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrcCov ),
1139  ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrc ),
1140  eOutputType, mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
1141  ( m_bPointCloudTarget ? nullptr : &dataGLLNodesDest ) ) );
1142 
1143  // Construct remap
1144  if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" );
1145  MB_CHK_ERR( LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ ) );
1146  }
1147  else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
1148  {
1149  DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;
1150 
1151  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1152  // generate metadata for the input meshes (both source and covering source)
1153  GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1154  dataGLLJacobianSrc );
1155  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1156  dataGLLJacobian );
1157 
1158  if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
1159  {
1160  _EXCEPTIONT( "Number of element does not match between metadata and "
1161  "input mesh" );
1162  }
1163 
1164  // Initialize coordinates for map
1165  this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1166  this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
1167 
1168  // Generate the continuous Jacobian for input mesh
1169  bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1170 
1171  if( eInputType == DiscretizationType_CGLL )
1172  {
1173  GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
1174  }
1175  else
1176  {
1177  GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
1178  }
1179 
1180  // Finite element input / Finite volume output
1181  MB_CHK_ERR( this->SetDOFmapAssociation( eInputType, mapOptions.nPin,
1182  ( eInputType == DiscretizationType_CGLL ), &dataGLLNodesSrcCov,
1183  &dataGLLNodesSrc, eOutputType, mapOptions.nPout, false, nullptr ) );
1184 
1185  // Generate remap
1186  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1187 
1188  if( strMapAlgorithm == "volumetric" )
1189  {
1190  _EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for"
1191  "GLL input mesh" );
1192  }
1193 
1194  this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1195  this->m_pdataGLLNodesOut = nullptr;
1196 
1197 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1198  LinearRemapSE4( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov, dataGLLJacobian,
1199  nMonotoneType, fContinuousIn, mapOptions.fNoConservation, mapOptions.fSparseConstraints,
1200  *this );
1201 #else
1202  LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
1203  mapOptions.fNoConservation, mapOptions.fSparseConstraints );
1204 #endif
1205  }
1206  else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
1207  {
1208  DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
1209  DataArray3D< double > dataGLLJacobianOut;
1210 
1211  // Input metadata
1212  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1213  // generate metadata for the input meshes (both source and covering source)
1214  GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1215  dataGLLJacobianSrc );
1216  // now coverage
1217  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1218  dataGLLJacobianIn );
1219  // Output metadata
1220  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1221  GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
1222  dataGLLJacobianOut );
1223 
1224  // Initialize coordinates for map
1225  this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1226  this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1227 
1228  // Generate the continuous Jacobian for input mesh
1229  bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1230 
1231  if( eInputType == DiscretizationType_CGLL )
1232  {
1233  GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
1234  }
1235  else
1236  {
1237  GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
1238  }
1239 
1240  // Generate the continuous Jacobian for output mesh
1241  bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
1242 
1243  if( eOutputType == DiscretizationType_CGLL )
1244  {
1245  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
1246  }
1247  else
1248  {
1249  GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
1250  }
1251 
1252  // Input Finite Element to Output Finite Element
1253  MB_CHK_ERR( this->SetDOFmapAssociation( eInputType, mapOptions.nPin,
1254  ( eInputType == DiscretizationType_CGLL ), &dataGLLNodesSrcCov,
1255  &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1256  ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest ) );
1257 
1258  this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1259  this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1260 
1261  // Generate remap
1262  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1263 
1264 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1265  LinearRemapGLLtoGLL_Integrated( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov,
1266  dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1267  this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1268  fContinuousIn, fContinuousOut, mapOptions.fSparseConstraints, *this );
1269 #else
1270  LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1271  this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1272  fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
1273 #endif
1274  }
1275  else
1276  {
1277  _EXCEPTIONT( "Not implemented" );
1278  }
1279 
1280 #ifdef MOAB_HAVE_EIGEN3
1281  copy_tempest_sparsemat_to_eigen3();
1282 #endif
1283 
1284 #ifdef MOAB_HAVE_MPI
1285  if (m_meshOverlap)
1286  {
1287  // Remove ghosted entities from overlap set
1288  moab::Range ghostedEnts;
1289  MB_CHK_ERR( m_remapper->GetOverlapAugmentedEntities( ghostedEnts ) );
1290  moab::EntityHandle m_meshOverlapSet = m_remapper->GetMeshSet( moab::Remapper::OverlapMesh );
1291  MB_CHK_SET_ERR( m_interface->remove_entities( m_meshOverlapSet, ghostedEnts ),
1292  "Deleting ghosted entities failed" );
1293  }
1294 #endif
1295  // Verify consistency, conservation and monotonicity, globally
1296  if( !mapOptions.fNoCheck )
1297  {
1298  if( is_root ) dbgprint.printf( 0, "Verifying map" );
1299  this->IsConsistent( 1.0e-8 );
1300  if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 );
1301 
1302  if( nMonotoneType != 0 )
1303  {
1304  this->IsMonotone( 1.0e-12 );
1305  }
1306  }
1307  }
1308  catch( Exception& e )
1309  {
1310  dbgprint.printf( 0, "%s", e.ToString().c_str() );
1311  return ( moab::MB_FAILURE );
1312  }
1313  catch( ... )
1314  {
1315  return ( moab::MB_FAILURE );
1316  }
1317  return moab::MB_SUCCESS;
1318 }
1319 
1320 ///////////////////////////////////////////////////////////////////////////////
1321 
1322 int moab::TempestOnlineMap::IsConsistent( double dTolerance )
1323 {
1324 #ifndef MOAB_HAVE_MPI
1325 
1326  return OfflineMap::IsConsistent( dTolerance );
1327 
1328 #else
1329 
1330  // Get map entries
1331  DataArray1D< int > dataRows;
1332  DataArray1D< int > dataCols;
1333  DataArray1D< double > dataEntries;
1334 
1335  // Calculate row sums
1336  DataArray1D< double > dRowSums;
1337  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1338  dRowSums.Allocate( m_mapRemap.GetRows() );
1339 
1340  for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1341  {
1342  dRowSums[dataRows[i]] += dataEntries[i];
1343  }
1344 
1345  // Verify all row sums are equal to 1
1346  int fConsistent = 0;
1347  for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
1348  {
1349  if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
1350  {
1351  fConsistent++;
1352  int rowGID = row_gdofmap[i];
1353  Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
1354  }
1355  }
1356 
1357  int ierr;
1358  int fConsistentGlobal = 0;
1359  ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1360  if( ierr != MPI_SUCCESS ) return -1;
1361 
1362  return fConsistentGlobal;
1363 #endif
1364 }
1365 
1366 ///////////////////////////////////////////////////////////////////////////////
1367 
1368 int moab::TempestOnlineMap::IsConservative( double dTolerance )
1369 {
1370 #ifndef MOAB_HAVE_MPI
1371 
1372  return OfflineMap::IsConservative( dTolerance );
1373 
1374 #else
1375  // return OfflineMap::IsConservative(dTolerance);
1376 
1377  int ierr;
1378  // Get map entries
1379  DataArray1D< int > dataRows;
1380  DataArray1D< int > dataCols;
1381  DataArray1D< double > dataEntries;
1382  const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
1383  const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
1384 
1385  // Calculate column sums
1386  std::vector< int > dColumnsUnique;
1387  std::vector< double > dColumnSums;
1388 
1389  int nColumns = m_mapRemap.GetColumns();
1390  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1391  dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
1392  dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
1393 
1394  for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
1395  {
1396  dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
1397 
1398  assert( dataCols[i] < m_nTotDofs_SrcCov );
1399 
1400  // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
1401  int colGID = this->GetColGlobalDoF( dataCols[i] ); // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1402  // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1403  dColumnsUnique[dataCols[i]] = colGID;
1404 
1405  // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
1406  // std::endl;
1407  }
1408 
1409  int rootProc = 0;
1410  std::vector< int > nElementsInProc;
1411  const int nDATA = 3;
1412  nElementsInProc.resize( size * nDATA );
1413  int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
1414  ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
1415  if( ierr != MPI_SUCCESS ) return -1;
1416 
1417  int nTotVals = 0, nTotColumns = 0; // nTotColumnsUnq = 0;
1418  std::vector< int > dColumnIndices;
1419  std::vector< double > dColumnSumsTotal;
1420  std::vector< int > displs, rcount;
1421  if( rank == rootProc )
1422  {
1423  displs.resize( size + 1, 0 );
1424  rcount.resize( size, 0 );
1425  int gsum = 0;
1426  for( int ir = 0; ir < size; ++ir )
1427  {
1428  nTotVals += nElementsInProc[ir * nDATA];
1429  nTotColumns += nElementsInProc[ir * nDATA + 1];
1430  // nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
1431 
1432  displs[ir] = gsum;
1433  rcount[ir] = nElementsInProc[ir * nDATA + 1];
1434  gsum += rcount[ir];
1435 
1436  // printf( "%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns, displs[ir], rcount[ir], gsum );
1437  }
1438 
1439  printf( "Total nnz: %d, global source elements = %d\n", nTotVals, gsum );
1440 
1441  dColumnIndices.resize( nTotColumns, -1 );
1442  dColumnSumsTotal.resize( nTotColumns, 0.0 );
1443  // dColumnSourceAreas.resize ( nTotColumns, 0.0 );
1444  }
1445 
1446  // Gather all ColumnSums to root process and accumulate
1447  // We expect that the sums of all columns equate to 1.0 within user specified tolerance
1448  // Need to do a gatherv here since different processes have different number of elements
1449  // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
1450  // MPI_SUM, 0, m_pcomm->comm());
1451  // Use .data() rather than &vec[0] -- on non-root ranks dColumnIndices /
1452  // dColumnSumsTotal are empty (only resized on root, see ~10 lines above),
1453  // and &vec[0] indexing into an empty vector is undefined behavior. The
1454  // .data() form returns nullptr for an empty vector, which MPI_Gatherv
1455  // ignores since recvcount on non-root paths is effectively zero.
1456  ierr = MPI_Gatherv( dColumnsUnique.data(), m_nTotDofs_SrcCov, MPI_INT, dColumnIndices.data(), rcount.data(),
1457  displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
1458  if( ierr != MPI_SUCCESS ) return -1;
1459  ierr = MPI_Gatherv( dColumnSums.data(), m_nTotDofs_SrcCov, MPI_DOUBLE, dColumnSumsTotal.data(), rcount.data(),
1460  displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
1461  if( ierr != MPI_SUCCESS ) return -1;
1462  // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
1463  // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
1464  // MPI_SUCCESS ) return -1;
1465 
1466  // Clean out unwanted arrays now
1467  dColumnSums.clear();
1468  dColumnsUnique.clear();
1469 
1470  // Verify all column sums equal the input Jacobian
1471  int fConservative = 0;
1472  if( rank == rootProc )
1473  {
1474  displs[size] = ( nTotColumns );
1475  // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
1476  std::map< int, double > dColumnSumsOnRoot;
1477  // std::map<int, double> dColumnSourceAreasOnRoot;
1478  for( int ir = 0; ir < size; ir++ )
1479  {
1480  for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
1481  {
1482  if( dColumnIndices[ips] < 0 ) continue;
1483  // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
1484  // assert( dColumnIndices[ips] < nTotColumnsUnq );
1485  dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips]; // / dColumnSourceAreas[ips];
1486  // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
1487  // dColumnSourceAreas[ dColumnIndices[ips] ]
1488  }
1489  }
1490 
1491  for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
1492  {
1493  // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
1494  if( fabs( it->second - 1.0 ) > dTolerance )
1495  {
1496  fConservative++;
1497  Announce( "TempestOnlineMap is not conservative in column "
1498  // "%i (%1.15e)", it->first, it->second );
1499  "%i (%1.15e)",
1500  it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
1501  }
1502  }
1503  }
1504 
1505  // TODO: Just do a broadcast from root instead of a reduction
1506  ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
1507  if( ierr != MPI_SUCCESS ) return -1;
1508 
1509  return fConservative;
1510 #endif
1511 }
1512 
1513 ///////////////////////////////////////////////////////////////////////////////
1514 
1515 int moab::TempestOnlineMap::IsMonotone( double dTolerance )
1516 {
1517 #ifndef MOAB_HAVE_MPI
1518 
1519  return OfflineMap::IsMonotone( dTolerance );
1520 
1521 #else
1522 
1523  // Get map entries
1524  DataArray1D< int > dataRows;
1525  DataArray1D< int > dataCols;
1526  DataArray1D< double > dataEntries;
1527 
1528  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1529 
1530  // Verify all entries are in the range [0,1]
1531  int fMonotone = 0;
1532  for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1533  {
1534  if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
1535  {
1536  fMonotone++;
1537 
1538  Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
1539  }
1540  }
1541 
1542  int ierr;
1543  int fMonotoneGlobal = 0;
1544  ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1545  if( ierr != MPI_SUCCESS ) return -1;
1546 
1547  return fMonotoneGlobal;
1548 #endif
1549 }
1550 
1551 ///////////////////////////////////////////////////////////////////////////////
1552 
1553 void moab::TempestOnlineMap::ComputeAdjacencyRelations( std::vector< std::unordered_set< int > >& vecAdjFaces,
1554  int nrings,
1555  const Range& entities,
1556  bool useMOABAdjacencies,
1557  Mesh* trMesh )
1558 {
1559  assert( nrings > 0 );
1560  assert( useMOABAdjacencies || trMesh != nullptr );
1561 
1562  const size_t nrows = vecAdjFaces.size();
1563  moab::MeshTopoUtil mtu( m_interface );
1564  for( size_t index = 0; index < nrows; index++ )
1565  {
1566  vecAdjFaces[index].insert( index ); // add self target face first
1567  {
1568  // Compute the adjacent faces to the target face
1569  if( useMOABAdjacencies )
1570  {
1571  moab::Range ents;
1572  // ents.insert( entities.index( entities[index] ) );
1573  ents.insert( entities[index] );
1574  moab::Range adjEnts;
1575  moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, nrings );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
1576  for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
1577  {
1578  // int adjIndex = m_interface->id_from_handle(*it)-1;
1579  int adjIndex = entities.index( *it );
1580  // printf("rank: %d, Element %lu, entity: %lu, adjIndex %d\n", rank, index, *it, adjIndex);
1581  if( adjIndex >= 0 ) vecAdjFaces[index].insert( adjIndex );
1582  }
1583  }
1584  else
1585  {
1586  /// Vector storing adjacent Faces.
1587  typedef std::pair< int, int > FaceDistancePair;
1588  typedef std::vector< FaceDistancePair > AdjacentFaceVector;
1589  AdjacentFaceVector adjFaces;
1590  Face& face = trMesh->faces[index];
1591  GetAdjacentFaceVectorByEdge( *trMesh, index, nrings * face.edges.size(), adjFaces );
1592 
1593  // Add the adjacent faces to the target face list
1594  for( auto adjFace : adjFaces )
1595  if( adjFace.first >= 0 )
1596  vecAdjFaces[index].insert( adjFace.first ); // map target face to source face
1597  }
1598  }
1599  }
1600 }
1601 
1602 moab::ErrorCode moab::TempestOnlineMap::ApplyWeights( moab::Tag srcSolutionTag,
1603  moab::Tag tgtSolutionTag,
1604  bool transpose,
1605  CAASType caasType,
1606  double default_projection )
1607 {
1608  std::vector< double > solSTagVals;
1609  std::vector< double > solTTagVals;
1610 
1611  moab::Range sents, tents;
1612  if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
1613  {
1614  if( m_remapper->point_cloud_source )
1615  {
1616  moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh );
1617  solSTagVals.resize( covSrcEnts.size(), default_projection );
1618  sents = covSrcEnts;
1619  }
1620  else
1621  {
1622  moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
1623  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1624  default_projection );
1625  sents = covSrcEnts;
1626  }
1627  if( m_remapper->point_cloud_target )
1628  {
1629  moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
1630  solTTagVals.resize( tgtEnts.size(), default_projection );
1631  tents = tgtEnts;
1632  }
1633  else
1634  {
1635  moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
1636  solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
1637  this->GetDestinationNDofsPerElement(),
1638  default_projection );
1639  tents = tgtEnts;
1640  }
1641  }
1642  else
1643  {
1644  moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
1645  moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
1646  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1647  default_projection );
1648  solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
1649  this->GetDestinationNDofsPerElement(),
1650  default_projection );
1651 
1652  sents = covSrcEnts;
1653  tents = tgtEnts;
1654  }
1655 
1656  // The tag data is np*np*n_el_src
1657  MB_CHK_SET_ERR( m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] ),
1658  "Getting local tag data failed" );
1659 
1660  // Compute the application of weights on the suorce solution data and store it in the
1661  // destination solution vector data Optionally, can also perform the transpose application of
1662  // the weight matrix. Set the 3rd argument to true if this is needed
1663  MB_CHK_SET_ERR( this->ApplyWeights( solSTagVals, solTTagVals, transpose ),
1664  "Applying remap operator onto source vector data failed" );
1665 
1666  // The tag data is np*np*n_el_dest
1667  MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1668  "Setting target tag data failed" );
1669 
1670  if( caasType != CAAS_NONE )
1671  {
1672  std::string tgtSolutionTagName;
1673  MB_CHK_SET_ERR( m_interface->tag_get_name( tgtSolutionTag, tgtSolutionTagName ), "Getting tag name failed" );
1674 
1675  // Perform CAAS iterations iteratively until convergence
1676  constexpr int nmax_caas_iterations = 10;
1677  double mismatch = 1.0;
1678  int caasIteration = 0;
1679  double initialMismatch = 0.0;
1680  while( ( fabs( mismatch / initialMismatch ) > 1e-15 && fabs( mismatch ) > 1e-15 ) &&
1681  caasIteration++ < nmax_caas_iterations ) // iterate until convergence or a maximum of 5 iterations
1682  {
1683  double dMassDiffPostGlobal;
1684  std::pair< double, double > mDefect =
1685  this->ApplyBoundsLimiting( solSTagVals, solTTagVals, caasType, caasIteration, mismatch );
1686 #ifdef MOAB_HAVE_MPI
1687  double dMassDiffPost = mDefect.second;
1688  MPI_Allreduce( &dMassDiffPost, &dMassDiffPostGlobal, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
1689 #else
1690  dMassDiffPostGlobal = mDefect.second;
1691 #endif
1692  if( caasIteration == 1 ) initialMismatch = mDefect.first;
1693  if( m_remapper->verbose && is_root )
1694  {
1695  printf( "Field {%s} -> CAAS iteration: %d, mass defect: %3.4e, post-CAAS: %3.4e\n",
1696  tgtSolutionTagName.c_str(), caasIteration, mDefect.first, dMassDiffPostGlobal );
1697  }
1698  mismatch = dMassDiffPostGlobal;
1699 
1700  // The tag data is np*np*n_el_dest
1701  MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1702  "Setting local tag data failed" );
1703  }
1704  }
1705 
1706  return moab::MB_SUCCESS;
1707 }
1708 
1709 moab::ErrorCode moab::TempestOnlineMap::ApplyWeightsWithDualMap( moab::Tag srcSolutionTag,
1710  moab::Tag tgtSolutionTag,
1711  TempestOnlineMap* loWeightMap,
1712  CAASType caasType )
1713 {
1714  // Setup entity ranges (same pattern as ApplyWeights(Tag, Tag))
1715  std::vector< double > solSTagVals, solTTagVals;
1716  moab::Range sents, tents;
1717 
1718  if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
1719  {
1720  if( m_remapper->point_cloud_source )
1721  {
1722  moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh );
1723  solSTagVals.resize( covSrcEnts.size(), 0.0 );
1724  sents = covSrcEnts;
1725  }
1726  else
1727  {
1728  moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
1729  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() *
1730  this->GetSourceNDofsPerElement(),
1731  0.0 );
1732  sents = covSrcEnts;
1733  }
1734  if( m_remapper->point_cloud_target )
1735  {
1736  moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
1737  solTTagVals.resize( tgtEnts.size(), 0.0 );
1738  tents = tgtEnts;
1739  }
1740  else
1741  {
1742  moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
1743  solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
1744  this->GetDestinationNDofsPerElement(),
1745  0.0 );
1746  tents = tgtEnts;
1747  }
1748  }
1749  else
1750  {
1751  moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
1752  moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
1753  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1754  0.0 );
1755  solTTagVals.resize(
1756  tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), 0.0 );
1757  sents = covSrcEnts;
1758  tents = tgtEnts;
1759  }
1760 
1761  // Read source tag data from coverage mesh
1762  MB_CHK_SET_ERR( m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] ),
1763  "Getting source tag data failed" );
1764 
1765  // Apply high-order projection only (no CAAS — bounds come from the low-order map)
1766  MB_CHK_SET_ERR( this->ApplyWeights( solSTagVals, solTTagVals, false ),
1767  "High-order projection failed" );
1768 
1769  // Write initial projection to target tag
1770  MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1771  "Setting target tag data failed" );
1772 
1773  if( caasType == CAAS_NONE || loWeightMap == nullptr ) return moab::MB_SUCCESS;
1774 
1775  // =====================================================================
1776  // Dual-map CAAS (Clip-And-Assured-Sum) — bit-for-bit port of MCT's
1777  // seq_nlmap_avNormArr (driver-mct/main/seq_nlmap_mod.F90).
1778  //
1779  // PURPOSE
1780  // Conservative, bounds-preserving remap of a source field x onto a
1781  // target mesh using TWO weight matrices: a high-order
1782  // non-monotone map (A, = `this`) and a low-order monotone &
1783  // conservative map (Am, = `loWeightMap`). The high-order map gives
1784  // accuracy; the low-order map gives the conservation reference and
1785  // the bounds-preservation safety net. This routine wires them
1786  // together using the Clip-And-Assured-Sum scheme of
1787  // Bradley, Bosler & Guba, "Conservation with bounded variation
1788  // and limiters in semi-Lagrangian transport schemes",
1789  // SIAM J. Sci. Comput. 41(5), 2019, doi:10.1137/18M1165414.
1790  //
1791  // NOTATION (matching the reference)
1792  // x source field values on the coverage mesh (solSTagVals)
1793  // A high-order map (this->m_weightMatrix)
1794  // Am low-order map (loWeightMap)
1795  // y_hi = A * x high-order projection (in solTTagVals)
1796  // y_lo = Am * x low-order projection (mass reference)
1797  // [lo, hi] per-row source-value bounds taken over A's stencil
1798  // gmins/gmaxs unscaled global min/max of the per-row [lo, hi] —
1799  // used as a final safety clip
1800  // norm8wt fractional-coverage weight (one scalar per source cell)
1801  // propagated by the E3SM driver as a side-channel tag;
1802  // when present, all bounds & redistribution arithmetic is
1803  // rescaled to match MCT's lnorm=.true. branch exactly
1804  //
1805  // ALGORITHM (one pass, FP-order-preserved vs MCT)
1806  // 1) y_hi = A * x (done above, in solTTagVals)
1807  // 2) y_lo = Am * x [Step 2]
1808  // 2b) Pull source norm8wt side-channel tag if present [Step 2b]
1809  // 3) Per-row bounds [lo, hi] over A's stencil columns [Step 3]
1810  // Divide source value by srcNorm8wt before tracking
1811  // min/max so bounds are over RECOVERED x, not (frac*x).
1812  // 4) y_lo == 0 mask: where the low-order projection is zero,
1813  // force y_hi = lo = hi = 0 to drop the cell. [Step 4]
1814  // 4b) Snapshot UNSCALED global extrema gmins/gmaxs from the
1815  // masked, but not-yet-norm-scaled, per-row [lo, hi]. [Step 4b]
1816  // 4c) mappedNorm8wt = Am * srcNorm8wt (or Am * 1 if absent). [Step 4c]
1817  // 4d) Scale per-row bounds: lo *= mappedNorm8wt, hi *= ... [Step 4d]
1818  // 5) Per-cell CAAS quantities (clipping defect, room to lower/raise) [Step 5]
1819  // 6) Reproducible global reductions of the per-cell quantities [Step 7]
1820  // 7) dM_total = dM_clip + (M_low - M_hi_unclipped) [Step 8]
1821  // 8) Redistribute the deficit across cells with room. [Step 9]
1822  // 9) Final hard clip to gmins/gmaxs (skipping yLow==0 cells).
1823  //
1824  // REPRODUCIBILITY MODEL
1825  // "BfB with MCT" means: for the same inputs, this routine produces
1826  // the same target values MCT produces, BIT FOR BIT, regardless of
1827  // MPI rank count or mesh decomposition. This requires three things
1828  // that the code below enforces explicitly:
1829  //
1830  // (i) Same area values. MCT uses 'aream' = area_b from the netcdf
1831  // map file. We read the same MOAB 'aream' tag (loaded by
1832  // iMOAB_LoadMapFile). Recomputing spherical-polygon areas via
1833  // lHuiller from mesh geometry is FP-different and is only used
1834  // as a final fallback for online-computed maps with no aream.
1835  //
1836  // (ii) Same FP operation order in the per-cell arithmetic. The
1837  // redistribute step computes `(hi - yc)/cap_g * dM_total`, NOT
1838  // the algebraically-equivalent `(hi - yc) * (dM_total/cap_g)`.
1839  // See Step 9 below for why. The dM_total computation also uses
1840  // MCT's exact two-step form (subtract, then add), preserving
1841  // the catastrophic-cancellation residual MCT carries.
1842  //
1843  // (iii) Order-independent global reductions. We use Worley's
1844  // IntegerReprosum (the MOAB port of shr_reprosum_int), which
1845  // is MCT's default reprosum path. It is decomposition- and
1846  // order-independent by construction (integer-vector MPI sum).
1847  // A Kahan + sort-by-gid summation lambda is also defined
1848  // below as a reference alternative but is not the active
1849  // reducer — using two different algorithms would defeat BfB.
1850  //
1851  // GUARDRAILS
1852  // Bounds extraction (Step 3) hard-aborts the run if any owned
1853  // high-order row references a coverage column not present on this
1854  // rank. Silently dropping such columns would produce
1855  // decomposition-dependent bounds and break BfB. The error message
1856  // tells the caller exactly which row/column/weight failed and
1857  // recommends widening the ghost-layer count (nghlay_cov in the
1858  // E3SM coupler driver). See lines below the bounds loop.
1859  //
1860  // EARLY RETURN
1861  // If caasType == CAAS_NONE or loWeightMap is null, we keep the raw
1862  // high-order projection that was already written to tgtSolutionTag
1863  // above. The dual-map machinery only runs when the caller explicitly
1864  // activates it with a non-null low-order map and a non-CAAS_NONE
1865  // filter type.
1866 
1867  const size_t nTargetDofs = solTTagVals.size();
1868  const size_t nSourceDofs = solSTagVals.size();
1869 
1870  // Map from target tag index to matrix row index. Both A and Am must share
1871  // the same row layout (same target mesh, same partitioning); this is true
1872  // because both maps are loaded onto the same intersection application.
1873  if( row_dtoc_dofmap.size() < nTargetDofs )
1874  {
1875  MB_CHK_SET_ERR( moab::MB_FAILURE, "row_dtoc_dofmap smaller than target tag size" );
1876  }
1877 
1878  // ----- Step 2: low-order projection y_lo = Am * x ---------------------
1879  std::vector< double > yLow( nTargetDofs, 0.0 );
1880  MB_CHK_SET_ERR( loWeightMap->ApplyWeights( solSTagVals, yLow, false ),
1881  "Low-order projection failed" );
1882 
1883  // ----- Step 2b: pull source norm8wt side-channel (if present) ---------
1884  //
1885  // BACKGROUND
1886  // When the E3SM driver coupler asks for a normalized projection
1887  // (lnorm=.true.), it does NOT send raw source values x to the
1888  // remapper. Instead, in seq_map_avNormArr it pre-multiplies each
1889  // source data field by a fractional-coverage weight `frac`, and
1890  // sends the products (frac * x) over to the intersection app
1891  // together with a parallel single-component tag named "norm8wt"
1892  // that carries (frac) on each source coverage cell.
1893  //
1894  // The driver later UN-DOES this pre-norm on the target side by
1895  // dividing each mapped data field by the mapped norm8wt — so the
1896  // final value on the target is (Am*(frac*x)) / (Am*frac). That
1897  // per-cell weighted average is the conservative answer when source
1898  // cells are only partially covered (e.g. land/ocean coastlines).
1899  //
1900  // WHY THE CAAS KERNEL NEEDS TO SEE norm8wt
1901  // To match MCT's seq_nlmap_avNormArr bit-for-bit, two things have
1902  // to happen INSIDE the CAAS kernel — neither can be done by the
1903  // driver as a post-pass:
1904  //
1905  // (a) Per-row bounds [lo, hi] must be the min/max of RECOVERED x
1906  // over the high-order stencil, not the min/max of (frac*x).
1907  // MCT does
1908  // tmp = solSTagVals[srcIdx]
1909  // tmp = tmp / xPrimeAV(natt+1, col) ! divide by frac
1910  // in sMat_avMult_and_calc_bounds before extending bounds, and
1911  // skips columns where frac == 0 (the field can't say anything
1912  // meaningful at a cell with no source coverage). Without this
1913  // divide, bounds would be 0-suppressed in coastal regions and
1914  // the CAAS clip would lose accuracy.
1915  //
1916  // (b) The mapped-norm8wt scale factor used in Step 4d must be the
1917  // LOW-ORDER projection of the ACTUAL source `frac`, not the
1918  // low-order projection of constant-1. MCT computes this in
1919  // the same mct_sMat_avMult call that produces avp_o data — the
1920  // natt+1 column gets sum_l w_lo[j,l] * frac(l), and that's
1921  // what the bounds get scaled by.
1922  //
1923  // FALLBACK
1924  // If no "norm8wt" tag exists on the intersection-side mesh (callers
1925  // that never pre-normed), we set hasNorm8wt=false. In that branch
1926  // srcNorm8wt is treated as constant-1 for both (a) and (b), which is
1927  // mathematically correct: with no pre-norm, frac would have been
1928  // 1.0 everywhere and the divide / scale are no-ops.
1929  //
1930  // SHAPE CONSTRAINT
1931  // norm8wt is single-component (one double per source coverage cell).
1932  // For the FV-FV configuration on the active CAAS path,
1933  // sents.size() == nSourceDofs and the tag values map directly to
1934  // solSTagVals indices. If a future caller wires an SE source layout
1935  // where nSourceDofs > sents.size() (multi-DOF per cell), the
1936  // per-cell norm8wt cannot be unambiguously expanded to per-DOF
1937  // values here — we deliberately fall back to the constant-1 path
1938  // rather than guess an expansion that would silently break BfB.
1939  std::vector< double > srcNorm8wt;
1940  bool hasNorm8wt = false;
1941  {
1942  moab::Tag normTag = nullptr;
1943  moab::ErrorCode rvalN = m_interface->tag_get_handle( "norm8wt", normTag );
1944  if( MB_SUCCESS == rvalN && normTag != nullptr )
1945  {
1946  // Single-component tag (one double per source coverage entity).
1947  // For FV-FV (the only configuration on the active CAAS path)
1948  // sents.size() == nSourceDofs. If the source layout is multi-DOF
1949  // (e.g. SE) the per-cell norm8wt cannot be unambiguously expanded
1950  // to per-DOF values here; bail to the constant-1 fallback rather
1951  // than guess.
1952  srcNorm8wt.resize( sents.size(), 0.0 );
1953  moab::ErrorCode rvalD = m_interface->tag_get_data( normTag, sents, &srcNorm8wt[0] );
1954  if( MB_SUCCESS == rvalD && srcNorm8wt.size() == nSourceDofs )
1955  {
1956  hasNorm8wt = true;
1957  }
1958  else
1959  {
1960  srcNorm8wt.clear();
1961  }
1962  }
1963  }
1964 
1965  // ----- Step 3: per-row bounds from HIGH-ORDER stencil -----------------
1966  // bounds(A, x): for each target row r, [lo, hi] = [min, max] of x over
1967  // A(r,:)'s nonzero columns. The proof in the reference paper requires
1968  // bounds to come from the larger (high-order) stencil so the constraint
1969  // set is provably nonempty.
1970  std::vector< double > lcl_lo( nTargetDofs, 1e308 );
1971  std::vector< double > lcl_hi( nTargetDofs, -1e308 );
1972 
1973  WeightMatrix& hiW = this->m_weightMatrix;
1974  for( size_t i = 0; i < nTargetDofs; i++ )
1975  {
1976  int r = row_dtoc_dofmap[i];
1977  if( r < 0 || r >= hiW.outerSize() ) continue;
1978  for( WeightMatrix::InnerIterator it( hiW, r ); it; ++it )
1979  {
1980  // it.col() is a matrix column index; map it to source vector index
1981  // by inverting col_dtoc_dofmap. For FV-FV with cell-based DOFs
1982  // and one-to-one mapping, this is the identity for owned columns.
1983  int mc = (int)it.col();
1984  // Search col_dtoc_dofmap[k]==mc; for typical FV cases the mapping
1985  // is dense and contiguous, so a linear scan over solSTagVals is
1986  // avoided by precomputing an inverse (below). For correctness we
1987  // fall back to scanning if the inverse is not available.
1988  // Build inverse once outside the loop (see below).
1989  (void)mc;
1990  }
1991  }
1992 
1993  // Precompute matrix-col -> source-vector-index inverse (cached per call;
1994  // O(nSourceDofs) construction). col_dtoc_dofmap has size nSourceDofs and
1995  // maps source-vector-index -> matrix-col.
1996  int maxMatCol = -1;
1997  for( size_t k = 0; k < nSourceDofs && k < col_dtoc_dofmap.size(); k++ )
1998  if( col_dtoc_dofmap[k] > maxMatCol ) maxMatCol = col_dtoc_dofmap[k];
1999  std::vector< int > col_inv( maxMatCol + 1, -1 );
2000  for( size_t k = 0; k < nSourceDofs && k < col_dtoc_dofmap.size(); k++ )
2001  if( col_dtoc_dofmap[k] >= 0 ) col_inv[col_dtoc_dofmap[k]] = (int)k;
2002 
2003  // Track whether any owned row's high-order stencil column failed to
2004  // resolve into the local coverage source vector. If that happens the
2005  // [lcl_lo, lcl_hi] bounds are computed over an INCOMPLETE stencil and
2006  // the CAAS clip + redistribute will produce decomposition-dependent
2007  // values — exactly the symptom seen as 1-2 ULP cross-rank-count drift
2008  // on file-loaded maps. Fail loudly with the offending coordinates so
2009  // the coverage layout (nghlay_cov in the calling code) can be widened.
2010  int bndsLocalErr = 0;
2011  int bndsFirstRowG = -1; // global target row id where the first failure happened
2012  int bndsFirstMc = -1; // matrix col index that failed to resolve
2013  double bndsFirstWgt = 0.0; // the dropped (nonzero) weight value
2014  int bndsFirstKind = 0; // 1 = mc out of maxMatCol; 2 = col_inv -> -1; 3 = srcIdx OOB
2015 
2016  for( size_t i = 0; i < nTargetDofs; i++ )
2017  {
2018  int r = row_dtoc_dofmap[i];
2019  if( r < 0 || r >= hiW.outerSize() ) continue;
2020  for( WeightMatrix::InnerIterator it( hiW, r ); it; ++it )
2021  {
2022  // Skip explicit-zero entries. Eigen's InnerIterator visits any
2023  // stored coefficient regardless of value; TempestRemap offline
2024  // maps routinely emit explicit zeros. MCT's
2025  // sMat_avMult_and_calc_bounds explicitly does
2026  // if (wgt == 0) cycle
2027  // before extending bounds (seq_nlmap_mod.F90:855). We can't use
2028  // an exact-zero compare here (FP-fragile), but 1e-50 is below
2029  // any physically meaningful map weight while still robust to
2030  // sign and denormal noise — entries this small can't shift the
2031  // per-row [lo, hi] enough to cross a clip threshold either.
2032  if( fabs( it.value() ) < 1e-50 ) continue;
2033  const int mc = (int)it.col();
2034 
2035  // Hard checks: a nonzero high-order weight at column mc means
2036  // this owned row genuinely depends on source-coverage column mc.
2037  // If we cannot resolve mc to a local source-vector index, the
2038  // 3-ring coverage on this rank is too narrow for the high-order
2039  // stencil. Either the caller asked for too few ghost layers,
2040  // or the map file references columns not present in any rank's
2041  // coverage (a catastrophic mismatch). Either way, silently
2042  // skipping corrupts the bounds and breaks BFB.
2043  if( mc < 0 || mc > maxMatCol )
2044  {
2045  if( !bndsLocalErr )
2046  {
2047  bndsLocalErr = 1;
2048  bndsFirstRowG = (r >= 0 && r < (int)row_gdofmap.size()) ? (int)row_gdofmap[r] : -1;
2049  bndsFirstMc = mc;
2050  bndsFirstWgt = it.value();
2051  bndsFirstKind = 1;
2052  }
2053  continue;
2054  }
2055  const int srcIdx = col_inv[mc];
2056  if( srcIdx < 0 )
2057  {
2058  if( !bndsLocalErr )
2059  {
2060  bndsLocalErr = 1;
2061  bndsFirstRowG = (r >= 0 && r < (int)row_gdofmap.size()) ? (int)row_gdofmap[r] : -1;
2062  bndsFirstMc = mc;
2063  bndsFirstWgt = it.value();
2064  bndsFirstKind = 2;
2065  }
2066  continue;
2067  }
2068  if( srcIdx >= (int)nSourceDofs )
2069  {
2070  if( !bndsLocalErr )
2071  {
2072  bndsLocalErr = 1;
2073  bndsFirstRowG = (r >= 0 && r < (int)row_gdofmap.size()) ? (int)row_gdofmap[r] : -1;
2074  bndsFirstMc = mc;
2075  bndsFirstWgt = it.value();
2076  bndsFirstKind = 3;
2077  }
2078  continue;
2079  }
2080  // solSTagVals[srcIdx] holds (frac * x) when the driver pre-normed
2081  // (hasNorm8wt true); divide by frac to recover x for bounds, and
2082  // skip the source cell when frac == 0 (matches MCT
2083  // seq_nlmap_mod.F90:857 "if xPrimeAV(natt+1,col) == 0 cycle").
2084  // When hasNorm8wt is false, solSTagVals already holds raw x.
2085  double v = solSTagVals[srcIdx];
2086  if( hasNorm8wt )
2087  {
2088  const double n = srcNorm8wt[srcIdx];
2089  if( fabs(n) < 1E-20 ) continue;
2090  v /= n;
2091  }
2092  if( v < lcl_lo[i] ) lcl_lo[i] = v;
2093  if( v > lcl_hi[i] ) lcl_hi[i] = v;
2094  }
2095  // If row had no nonzero columns in the high-order stencil, set bounds
2096  // to 0 (matching MCT's sMat_avMult_and_calc_bounds: "lop = 0; hip = 0").
2097  // Together with the y_lo == 0 masking step below, this forces the cell
2098  // to 0 — a "rare, local reduction in order to one" per the reference.
2099  if( lcl_lo[i] > lcl_hi[i] )
2100  {
2101  lcl_lo[i] = 0.0;
2102  lcl_hi[i] = 0.0;
2103  }
2104  }
2105 
2106  // Globalize: any rank with bndsLocalErr triggers a collective failure.
2107 #ifdef MOAB_HAVE_MPI
2108  {
2109  MPI_Comm comm = m_pcomm ? m_pcomm->comm() : MPI_COMM_SELF;
2110  int bndsGlobalErr = 0;
2111  MPI_Allreduce( &bndsLocalErr, &bndsGlobalErr, 1, MPI_INT, MPI_MAX, comm );
2112  if( bndsGlobalErr )
2113  {
2114  int myRank = 0;
2115  MPI_Comm_rank( comm, &myRank );
2116  if( bndsLocalErr )
2117  {
2118  static const char* kindStr[4] = { "?", "mc>maxMatCol", "col_inv[mc]==-1", "srcIdx>=nSourceDofs" };
2119  fprintf( stderr,
2120  "FATAL: ApplyWeightsWithDualMap bounds extraction dropped a nonzero "
2121  "high-order stencil column on rank %d.\n"
2122  " global_target_row=%d matrix_col=%d weight=%.17e reason=%s\n"
2123  " This means the source coverage on this rank does NOT contain a "
2124  "column the owned high-order row references — the 3-ring (or whatever) "
2125  "ghost layer setting is too narrow, or the map file was generated against "
2126  "a different mesh. Bounds computed over an incomplete stencil break BFB; "
2127  "aborting rather than silently producing wrong CAAS output.\n",
2128  myRank, bndsFirstRowG, bndsFirstMc, bndsFirstWgt, kindStr[bndsFirstKind] );
2129  fflush( stderr );
2130  }
2131  MPI_Abort( comm, 1 );
2132  }
2133  }
2134 #else
2135  if( bndsLocalErr )
2136  {
2137  static const char* kindStr[4] = { "?", "mc>maxMatCol", "col_inv[mc]==-1", "srcIdx>=nSourceDofs" };
2138  fprintf( stderr,
2139  "FATAL: ApplyWeightsWithDualMap bounds extraction dropped a nonzero "
2140  "high-order stencil column.\n"
2141  " global_target_row=%d matrix_col=%d weight=%.17e reason=%s\n",
2142  bndsFirstRowG, bndsFirstMc, bndsFirstWgt, kindStr[bndsFirstKind] );
2143  fflush( stderr );
2144  return moab::MB_FAILURE;
2145  }
2146 #endif
2147 
2148  // ----- Step 4: mask -- where y_lo == 0, zero out y_hi and bounds ------
2149  // (Per reference: "An exact 0 in the low-order field will mask the
2150  // high-order field unnecessarily, but that's OK: it's a rare, local
2151  // reduction in order to one, not a wrong value.")
2152  for( size_t i = 0; i < nTargetDofs; i++ )
2153  {
2154  if( yLow[i] == 0.0 )
2155  {
2156  solTTagVals[i] = 0.0;
2157  lcl_lo[i] = 0.0;
2158  lcl_hi[i] = 0.0;
2159  }
2160  }
2161 
2162  // ----- Step 4b: compute UNSCALED global extrema for the final safety
2163  // clip (Item 4). MCT's seq_nlmap_avNormArr clips the redistributed result
2164  // against gmins/gmaxs = global min/max of the per-row (post-mask) bounds,
2165  // not against the per-row bounds themselves. Snapshot here, BEFORE the
2166  // bounds get scaled by the mapped norm in Step 4d below.
2167  double g_lo = 1e308, g_hi = -1e308;
2168  for( size_t i = 0; i < nTargetDofs; i++ )
2169  {
2170  int r = row_dtoc_dofmap[i];
2171  if( r < 0 || r >= (int)row_gdofmap.size() ) continue; // not owned
2172  if( lcl_lo[i] < g_lo ) g_lo = lcl_lo[i];
2173  if( lcl_hi[i] > g_hi ) g_hi = lcl_hi[i];
2174  }
2175 #ifdef MOAB_HAVE_MPI
2176  {
2177  MPI_Comm comm = m_pcomm ? m_pcomm->comm() : MPI_COMM_SELF;
2178  double tmp_min = g_lo, tmp_max = g_hi;
2179  MPI_Allreduce( &tmp_min, &g_lo, 1, MPI_DOUBLE, MPI_MIN, comm );
2180  MPI_Allreduce( &tmp_max, &g_hi, 1, MPI_DOUBLE, MPI_MAX, comm );
2181  }
2182 #endif
2183 
2184  // ----- Step 4c: compute mapped norm8wt = low-order map applied to a
2185  // source-norm8wt vector. For target row i this equals
2186  // sum_l w_lo[i,l] * srcNorm8wt(l)
2187  // — equivalently, the value MCT carries in the natt+1 column of avp_o
2188  // after mct_sMat_avMult is applied to avp_i (whose norm8wt slot holds
2189  // frac post-pre-norm). When no "norm8wt" tag is available on the intx
2190  // side, fall back to applying the low-order map to a constant-1 vector
2191  // (equivalent to srcNorm8wt(l) == 1 everywhere — consistent with the
2192  // bounds-extraction fallback above).
2193  std::vector< double > mappedNorm8wt( nTargetDofs, 0.0 );
2194  if( hasNorm8wt )
2195  {
2196  MB_CHK_SET_ERR( loWeightMap->ApplyWeights( srcNorm8wt, mappedNorm8wt, false ),
2197  "Mapped-norm8wt computation (low-order on source norm8wt) failed" );
2198  }
2199  else
2200  {
2201  std::vector< double > srcOnes( nSourceDofs, 1.0 );
2202  MB_CHK_SET_ERR( loWeightMap->ApplyWeights( srcOnes, mappedNorm8wt, false ),
2203  "Mapped-norm8wt computation (low-order on ones) failed" );
2204  }
2205 
2206  // ----- Step 4d: scale per-row bounds by mapped norm8wt (Item 2).
2207  // MCT does this inside the CAAS loop:
2208  // if (lnorm) then
2209  // lo = lo*avp_o%rAttr(natt+1,j)
2210  // hi = hi*avp_o%rAttr(natt+1,j)
2211  // end if
2212  // Doing it once here propagates correctly into Step 5 (clipping) and
2213  // Step 9 (redistribution) which both use lcl_lo/lcl_hi. Note: g_lo/g_hi
2214  // were already snapshotted above and remain UNSCALED (matching MCT's
2215  // gmins/gmaxs which are the global min/max of the unscaled per-row bounds).
2216  for( size_t i = 0; i < nTargetDofs; i++ )
2217  {
2218  const double w = mappedNorm8wt[i];
2219  lcl_lo[i] *= w;
2220  lcl_hi[i] *= w;
2221  }
2222 
2223  // ----- Get target areas (per matrix-row) ------------------------------
2224  // For BFB with MCT we MUST use the same area values MCT does. MCT uses
2225  // 'aream' (= area_b from the netcdf map file, loaded once when the map
2226  // is read). iMOAB_LoadMapFile populates the 'aream' tag on the target
2227  // mesh from area_b when arearead != 0 (e.g. arearead=3 for F-maps).
2228  //
2229  // The CAAS path MUST NOT recompute spherical-polygon areas from mesh
2230  // geometry via lHuiller (or any other re-derivation): doing so differs
2231  // from area_b at FP precision and silently breaks BfB with MCT. If the
2232  // caller has not loaded an area-bearing map and there are no online
2233  // areas (m_dTargetAreas) either, fail the run loudly so the caller can
2234  // fix their map-load configuration instead of getting silent non-BfB
2235  // results.
2236  std::vector< double > tgtAreas( nTargetDofs, 0.0 );
2237  {
2238  std::vector< moab::EntityHandle > tentVec;
2239  tentVec.reserve( tents.size() );
2240  for( moab::Range::iterator it = tents.begin(); it != tents.end(); ++it )
2241  tentVec.push_back( *it );
2242 
2243  bool got_areas = false;
2244 
2245  // Preferred: pull the 'aream' tag from the target MOAB mesh — this
2246  // is the area_b value loaded by iMOAB_LoadMapFile and is byte-identical
2247  // to the 'aream' field MCT uses in seq_nlmap_avNormArr.
2248  moab::Tag aream_tag = nullptr;
2249  moab::ErrorCode rval = m_interface->tag_get_handle( "aream", aream_tag );
2250  if( MB_SUCCESS == rval && aream_tag != nullptr && !tentVec.empty() )
2251  {
2252  const size_t nents = std::min< size_t >( tentVec.size(), nTargetDofs );
2253  std::vector< double > aream_vals( nents, 0.0 );
2254  rval = m_interface->tag_get_data( aream_tag, &tentVec[0], (int)nents, &aream_vals[0] );
2255  if( MB_SUCCESS == rval )
2256  {
2257  for( size_t i = 0; i < nents; i++ ) tgtAreas[i] = aream_vals[i];
2258  got_areas = true;
2259  }
2260  }
2261 
2262  // Fallback: areas were computed online and live in OfflineMap's
2263  // m_dTargetAreas (indexed by matrix row). This is BFB with MCT only
2264  // when the same online-area code path is used on both couplers; it
2265  // is acceptable for runs that build the map online.
2266  if( !got_areas )
2267  {
2268  const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
2269  const size_t nRows = dTargetAreas.GetRows();
2270  if( nRows >= nTargetDofs )
2271  {
2272  for( size_t i = 0; i < nTargetDofs; i++ )
2273  {
2274  int r = row_dtoc_dofmap[i];
2275  if( r >= 0 && (size_t)r < nRows )
2276  tgtAreas[i] = dTargetAreas[r];
2277  }
2278  got_areas = true;
2279  }
2280  }
2281 
2282  // No fallback to lHuiller. Recomputing areas from mesh geometry is
2283  // not BFB with MCT and there is no safe silent default — abort.
2284  if( !got_areas )
2285  {
2286  MB_SET_ERR( moab::MB_FAILURE,
2287  "ApplyWeightsWithDualMap: no target-cell areas available. "
2288  "Neither the 'aream' tag (from iMOAB_LoadMapFile with "
2289  "arearead != 0) nor OfflineMap::GetTargetAreas() (from an "
2290  "online map build) provided areas. Recomputing areas from "
2291  "mesh geometry is not bit-for-bit with MCT and is no longer "
2292  "permitted in the CAAS path. Re-load the map file with an "
2293  "area-bearing arearead setting (e.g. arearead=3 for F-maps), "
2294  "or build the online map so target areas are populated." );
2295  }
2296  }
2297 
2298  // ----- Step 5: build per-cell CAAS weights ----------------------------
2299  // For BFB summation, we accumulate per-row (gid, value) pairs and reduce
2300  // them deterministically.
2301  std::vector< int > rowGids( nTargetDofs, -1 );
2302  std::vector< double > massLowPerRow( nTargetDofs, 0.0 );
2303  std::vector< double > massHiUnclippedPerRow( nTargetDofs, 0.0 ); // y_hi BEFORE clipping
2304  std::vector< double > clipDefectPerRow( nTargetDofs, 0.0 );
2305  std::vector< double > capLowPerRow( nTargetDofs, 0.0 );
2306  std::vector< double > capHighPerRow( nTargetDofs, 0.0 );
2307 
2308  for( size_t i = 0; i < nTargetDofs; i++ )
2309  {
2310  int r = row_dtoc_dofmap[i];
2311  if( r < 0 || r >= (int)row_gdofmap.size() )
2312  rowGids[i] = -1; // not owned by this rank
2313  else
2314  rowGids[i] = (int)row_gdofmap[r];
2315 
2316  const double area = tgtAreas[i];
2317  const double y = solTTagVals[i]; // y_hi (pre-clip)
2318  const double lo = lcl_lo[i];
2319  const double hi = lcl_hi[i];
2320  double yc = y; // clipped value
2321  double dm = 0.0;
2322  if( y < lo )
2323  {
2324  yc = lo;
2325  dm = ( y - lo ) * area; // negative: cell exceeded below
2326  }
2327  else if( y > hi )
2328  {
2329  yc = hi;
2330  dm = ( y - hi ) * area; // positive: cell exceeded above
2331  }
2332  clipDefectPerRow[i] = dm;
2333  capLowPerRow[i] = ( yc - lo ) * area; // room to subtract
2334  capHighPerRow[i] = ( hi - yc ) * area; // room to add
2335  massLowPerRow[i] = yLow[i] * area;
2336  // Per-row mass of the UNCLIPPED high-order projection. MCT reduces
2337  // exactly this quantity to obtain glbl_masses(natt+k) (M_hi_unclipped),
2338  // and then computes dM_total = dM_clip + (M_low - M_hi_unclipped) in
2339  // that 2-step order. We store the unclipped y here (rather than the
2340  // clipped yc as before) so MOAB's reprosum byte-matches MCT's, which
2341  // in turn lets the MCT-matching dM_total formula below produce the
2342  // same last bits as MCT.
2343  massHiUnclippedPerRow[i] = y * area;
2344  // Update solTTagVals to the clipped value for the next stage
2345  solTTagVals[i] = yc;
2346  }
2347 
2348 
2349  // ----- Step 6: BFB-deterministic global reductions --------------------
2350 
2351  // Reproducible global reductions via Worley's integer-vector algorithm
2352  // (moab::IntegerReprosum) — bit-identical to MCT's shr_reprosum_int
2353  // regardless of MPI rank count, mesh decomposition, or local iteration
2354  // order. This is MCT's default reprosum path (the namelist default
2355  // repro_sum_use_ddpdd=.false. on the MCT side). The integer-vector
2356  // algorithm is order-independent by construction (MPI_Allreduce with
2357  // MPI_SUM on int64), which eliminates the cross-rank-count ULP drift
2358  // that an order-sensitive reducer (e.g. Kahan or DDPDD) would otherwise
2359  // leak into the CAAS bounds and mass totals.
2360 #ifdef MOAB_HAVE_MPI
2361  MPI_Comm reduce_comm = m_pcomm ? m_pcomm->comm() : MPI_COMM_SELF;
2362 #else
2363  int reduce_comm = 0; // serial build: comm unused but kept for API symmetry
2364 #endif
2365 #ifdef MOAB_HAVE_MPI
2366  moab::IntegerReprosum repro( reduce_comm );
2367 #else
2368  moab::IntegerReprosum repro;
2369 #endif
2370  // Build the ownership mask once (rowGids[i] >= 0 ↔ owned).
2371  const std::vector< int >& reduce_mask = rowGids;
2372  // Batched reduction: one MPI_Allreduce for the per-field metadata
2373  // (gmax_exp / gmin_exp / max_nsummands across the 5 fields) and one
2374  // MPI_Allreduce for the concatenated integer-vector encoding of all 5
2375  // fields. Bit-for-bit identical to calling sum_masked() five times
2376  // separately (each field still uses its own per-field metadata and
2377  // decode pass), but goes from 10 collective calls to 2.
2378  const std::vector< std::vector< double > > caasFields = {
2379  massLowPerRow, massHiUnclippedPerRow, clipDefectPerRow,
2380  capLowPerRow, capHighPerRow };
2381  std::vector< double > caasGsums;
2382  repro.sum_masked_batch( caasFields, reduce_mask, caasGsums );
2383  const double M_low = caasGsums[0];
2384  const double M_hi_unclipped = caasGsums[1];
2385  const double dM_clip = caasGsums[2];
2386  const double cap_low_g = caasGsums[3];
2387  const double cap_high_g = caasGsums[4];
2388 
2389  // ----- Step 8: total mass deficit between low-order and clipped high-order
2390  // The redistribution must drive the (clipped) high-order solution back to
2391  // the low-order mass. MCT's seq_nlmap_avNormArr (line 616) computes this
2392  // in EXACTLY the following 2-step form, and floating-point rounding makes
2393  // it FP-different from the algebraically-equivalent (M_low - M_hi_clip):
2394  //
2395  // ! MCT (Fortran array assignment, evaluated element-wise)
2396  // gwts(k) = gwts(k) ! gwts(k) holds dM_clip after reprosum
2397  // + (glbl_masses(k) ! M_low
2398  // - glbl_masses(natt+k)) ! M_hi_unclipped
2399  //
2400  // Reproducing MCT's exact bit pattern requires:
2401  // (a) reducing the UNCLIPPED per-cell high-order mass directly via
2402  // reprosum, NOT deriving it as M_hi_clip + dM_clip — that derivation
2403  // drops 1-2 ULP because reprosum is exact only on its inputs.
2404  // => see massHiUnclippedPerRow above.
2405  // (b) computing dM_total in MCT's order: subtract first, then add.
2406  //
2407  // Sign: dM_total > 0 means low-order carries more mass than the clipped
2408  // high-order, so we need to ADD mass; dM_total < 0 means we need to REMOVE.
2409  //
2410  const double diff = M_low - M_hi_unclipped; // step 1: subtraction
2411  const double dM_total = dM_clip + diff; // step 2: addition (MCT order)
2412 
2413  // ----- Step 9: redistribute -------------------------------------------
2414  // For BfB with MCT seq_nlmap_avNormArr we MUST match its FP operation
2415  // order exactly. MCT does, per cell:
2416  // y = max(lo, min(hi, nl_avp_o(k,j))) ! re-clip
2417  // nl_avp_o(k,j) = y + ((hi - y)/tmp)*gwts(k) ! line 648 / 662
2418  // i.e. divide-then-multiply, with the loop-invariant denominator
2419  // (cap_high_g or cap_low_g) and numerator (dM_total) NOT precomputed
2420  // into a single `scale = dM_total/cap_g`. Doing so introduces ULP-level
2421  // per-cell differences (a/b*c reorders to (c/b)*a). Similarly the
2422  // earlier MOAB pattern computed `room = (hi-yc)*area` and then
2423  // `room/area`, which doesn't algebraically cancel in FP and added two
2424  // extra roundings per cell. The straightforward `(hi - yc)/cap_g *
2425  // dM_total` form below matches MCT bit-for-bit.
2426  if( dM_total > 0.0 && cap_high_g > 0.0 )
2427  {
2428  for( size_t i = 0; i < nTargetDofs; i++ )
2429  {
2430  const double area = tgtAreas[i];
2431  const double yc = solTTagVals[i];
2432  if( area > 0.0 )
2433  solTTagVals[i] = yc + ( ( lcl_hi[i] - yc ) / cap_high_g ) * dM_total;
2434  }
2435  }
2436  else if( dM_total < 0.0 && cap_low_g > 0.0 )
2437  {
2438  for( size_t i = 0; i < nTargetDofs; i++ )
2439  {
2440  const double area = tgtAreas[i];
2441  const double yc = solTTagVals[i];
2442  if( area > 0.0 )
2443  solTTagVals[i] = yc + ( ( yc - lcl_lo[i] ) / cap_low_g ) * dM_total;
2444  }
2445  }
2446 
2447  // Final hard clip for floating-point safety, against UNSCALED global
2448  // extrema (Item 4). MCT's seq_nlmap_avNormArr does:
2449  // if (avp_o(k,j) == 0) cycle ! 0-mask skip
2450  // nl_avp_o%rAttr(k,j) = max(gmins(k), min(gmaxs(k), nl_avp_o%rAttr(k,j)))
2451  // Per-row bounds (lcl_lo/lcl_hi) are now SCALED by mapped_norm8wt and so
2452  // would be a tighter clip than MCT applies; using global g_lo/g_hi keeps
2453  // the safety net loose, as the reference algorithm intends. The 0-mask
2454  // skip is critical: without it, target cells that were zeroed in Step 4
2455  // (yLow == 0) get bumped from 0 up to g_lo when g_lo > 0 (e.g.
2456  // positive-only fields like temperature/pressure), and the post-norm
2457  // divide in the driver then amplifies that wrong value by 1/wghts at
2458  // coastal coverage cells where wghts is tiny but nonzero.
2459  for( size_t i = 0; i < nTargetDofs; i++ )
2460  {
2461  if( fabs(yLow[i]) < 1E-40 ) continue;
2462  if( solTTagVals[i] < g_lo ) solTTagVals[i] = g_lo;
2463  if( solTTagVals[i] > g_hi ) solTTagVals[i] = g_hi;
2464  }
2465 
2466  // Store result back to the target tag
2467  MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
2468  "Setting target tag data failed" );
2469 
2470  return moab::MB_SUCCESS;
2471 }
2472 
2473 moab::ErrorCode moab::TempestOnlineMap::DefineAnalyticalSolution( moab::Tag& solnTag,
2474  const std::string& solnName,
2475  moab::Remapper::IntersectionContext ctx,
2476  sample_function testFunction,
2477  moab::Tag* clonedSolnTag,
2478  std::string cloneSolnName )
2479 {
2480  const bool outputEnabled = ( is_root );
2481  int discOrder;
2482  DiscretizationType discMethod;
2483  // moab::EntityHandle meshset;
2484  moab::Range entities;
2485  Mesh* trmesh;
2486  switch( ctx )
2487  {
2488  case Remapper::SourceMesh:
2489  // meshset = m_remapper->m_covering_source_set;
2490  trmesh = m_remapper->m_covering_source;
2491  entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
2492  : m_remapper->m_covering_source_entities );
2493  discOrder = m_nDofsPEl_Src;
2494  discMethod = m_eInputType;
2495  break;
2496 
2497  case Remapper::TargetMesh:
2498  // meshset = m_remapper->m_target_set;
2499  trmesh = m_remapper->m_target;
2500  entities =
2501  ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
2502  discOrder = m_nDofsPEl_Dest;
2503  discMethod = m_eOutputType;
2504  break;
2505 
2506  default:
2507  if( outputEnabled )
2508  std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
2509  return moab::MB_FAILURE;
2510  }
2511 
2512  // Let us create teh solution tag with appropriate information for name, discretization order
2513  // (DoF space)
2514  MB_CHK_ERR( m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
2515  MB_TAG_DENSE | MB_TAG_CREAT ) );
2516  if( clonedSolnTag != nullptr )
2517  {
2518  if( cloneSolnName.size() == 0 )
2519  {
2520  cloneSolnName = solnName + std::string( "Cloned" );
2521  }
2522  MB_CHK_ERR( m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
2523  *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT ) );
2524  }
2525 
2526  // Triangular quadrature rule
2527  const int TriQuadratureOrder = 10;
2528 
2529  if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
2530 
2531  TriangularQuadratureRule triquadrule( TriQuadratureOrder );
2532 
2533  const int TriQuadraturePoints = triquadrule.GetPoints();
2534 
2535  const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
2536  const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
2537 
2538  // Output data
2539  DataArray1D< double > dVar;
2540  DataArray1D< double > dVarMB; // re-arranged local MOAB vector
2541 
2542  // Nodal geometric area
2543  DataArray1D< double > dNodeArea;
2544 
2545  // Calculate element areas
2546  // trmesh->CalculateFaceAreas(fContainsConcaveFaces);
2547 
2548  if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
2549  {
2550  /* Get the spectral points and sample the functionals accordingly */
2551  const bool fGLL = true;
2552  const bool fGLLIntegrate = false;
2553 
2554  // Generate grid metadata
2555  DataArray3D< int > dataGLLNodes;
2556  DataArray3D< double > dataGLLJacobian;
2557 
2558  GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );
2559 
2560  // Number of elements
2561  int nElements = trmesh->faces.size();
2562 
2563  // Verify all elements are quadrilaterals
2564  for( int k = 0; k < nElements; k++ )
2565  {
2566  const Face& face = trmesh->faces[k];
2567 
2568  if( face.edges.size() != 4 )
2569  {
2570  _EXCEPTIONT( "Non-quadrilateral face detected; "
2571  "incompatible with --gll" );
2572  }
2573  }
2574 
2575  // Number of unique nodes (CGLL) or total element-local DOFs (DGLL)
2576  const bool fDiscontinuous = ( discMethod == DiscretizationType_DGLL );
2577  int iMaxNode = 0;
2578  if( fDiscontinuous )
2579  {
2580  // DGLL: each element has independent DOFs
2581  iMaxNode = nElements * discOrder * discOrder;
2582  }
2583  else
2584  {
2585  // CGLL: shared nodes at element boundaries
2586  for( int i = 0; i < discOrder; i++ )
2587  for( int j = 0; j < discOrder; j++ )
2588  for( int k = 0; k < nElements; k++ )
2589  if( dataGLLNodes[i][j][k] > iMaxNode )
2590  iMaxNode = dataGLLNodes[i][j][k];
2591  }
2592 
2593  // Get Gauss-Lobatto quadrature nodes
2594  DataArray1D< double > dG;
2595  DataArray1D< double > dW;
2596 
2597  GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
2598 
2599  // Get Gauss quadrature nodes
2600  const int nGaussP = 10;
2601 
2602  DataArray1D< double > dGaussG;
2603  DataArray1D< double > dGaussW;
2604 
2605  GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
2606 
2607  // Allocate data
2608  dVar.Allocate( iMaxNode );
2609  dVarMB.Allocate( discOrder * discOrder * nElements );
2610  dNodeArea.Allocate( iMaxNode );
2611 
2612  // Sample data
2613  for( int k = 0; k < nElements; k++ )
2614  {
2615  const Face& face = trmesh->faces[k];
2616 
2617  // Sample data at GLL nodes
2618  if( fGLL )
2619  {
2620  for( int i = 0; i < discOrder; i++ )
2621  {
2622  for( int j = 0; j < discOrder; j++ )
2623  {
2624 
2625  // Apply local map
2626  Node node;
2627  Node dDx1G;
2628  Node dDx2G;
2629 
2630  ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
2631 
2632  // Sample data at this point
2633  double dNodeLon = atan2( node.y, node.x );
2634  if( dNodeLon < 0.0 )
2635  {
2636  dNodeLon += 2.0 * M_PI;
2637  }
2638  double dNodeLat = asin( node.z );
2639 
2640  double dSample = ( *testFunction )( dNodeLon, dNodeLat );
2641 
2642  if( fDiscontinuous )
2643  dVar[k * discOrder * discOrder + j * discOrder + i] = dSample;
2644  else
2645  dVar[dataGLLNodes[j][i][k] - 1] = dSample;
2646  }
2647  }
2648  // High-order Gaussian integration over basis function
2649  }
2650  else
2651  {
2652  DataArray2D< double > dCoeff( discOrder, discOrder );
2653 
2654  for( int p = 0; p < nGaussP; p++ )
2655  {
2656  for( int q = 0; q < nGaussP; q++ )
2657  {
2658 
2659  // Apply local map
2660  Node node;
2661  Node dDx1G;
2662  Node dDx2G;
2663 
2664  ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
2665 
2666  // Cross product gives local Jacobian
2667  Node nodeCross = CrossProduct( dDx1G, dDx2G );
2668 
2669  double dJacobian =
2670  sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
2671 
2672  // Find components of quadrature point in basis
2673  // of the first Face
2674  SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
2675 
2676  // Sample data at this point
2677  double dNodeLon = atan2( node.y, node.x );
2678  if( dNodeLon < 0.0 )
2679  {
2680  dNodeLon += 2.0 * M_PI;
2681  }
2682  double dNodeLat = asin( node.z );
2683 
2684  double dSample = ( *testFunction )( dNodeLon, dNodeLat );
2685 
2686  // Integrate
2687  for( int i = 0; i < discOrder; i++ )
2688  {
2689  for( int j = 0; j < discOrder; j++ )
2690  {
2691 
2692  double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
2693 
2694  dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
2695 
2696  dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
2697  }
2698  }
2699  }
2700  }
2701  }
2702  }
2703 
2704  // Divide by area
2705  if( fGLLIntegrate )
2706  {
2707  for( size_t i = 0; i < dVar.GetRows(); i++ )
2708  {
2709  dVar[i] /= dNodeArea[i];
2710  }
2711  }
2712 
2713  // Let us rearrange the data based on DoF ID specification
2714  if( ctx == Remapper::SourceMesh )
2715  {
2716  for( unsigned j = 0; j < entities.size(); j++ )
2717  for( int p = 0; p < discOrder; p++ )
2718  for( int q = 0; q < discOrder; q++ )
2719  {
2720  const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
2721  dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]];
2722  }
2723  }
2724  else
2725  {
2726  for( unsigned j = 0; j < entities.size(); j++ )
2727  for( int p = 0; p < discOrder; p++ )
2728  for( int q = 0; q < discOrder; q++ )
2729  {
2730  const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
2731  dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]];
2732  }
2733  }
2734 
2735  // Set the tag data
2736  MB_CHK_ERR( m_interface->tag_set_data( solnTag, entities, &dVarMB[0] ) );
2737  }
2738  else
2739  {
2740  // assert( discOrder == 1 );
2741  if( discMethod == DiscretizationType_FV )
2742  {
2743  /* Compute an element-wise integral to store the sampled solution based on Quadrature
2744  * rules */
2745  // Resize the array
2746  dVar.Allocate( trmesh->faces.size() );
2747 
2748  std::vector< Node >& nodes = trmesh->nodes;
2749 
2750  // Loop through all Faces
2751  for( size_t i = 0; i < trmesh->faces.size(); i++ )
2752  {
2753  const Face& face = trmesh->faces[i];
2754 
2755  // Loop through all sub-triangles
2756  for( size_t j = 0; j < face.edges.size() - 2; j++ )
2757  {
2758 
2759  const Node& node0 = nodes[face[0]];
2760  const Node& node1 = nodes[face[j + 1]];
2761  const Node& node2 = nodes[face[j + 2]];
2762 
2763  // Triangle area
2764  Face faceTri( 3 );
2765  faceTri.SetNode( 0, face[0] );
2766  faceTri.SetNode( 1, face[j + 1] );
2767  faceTri.SetNode( 2, face[j + 2] );
2768 
2769  double dTriangleArea = CalculateFaceArea( faceTri, nodes );
2770 
2771  // Calculate the element average
2772  double dTotalSample = 0.0;
2773 
2774  // Loop through all quadrature points
2775  for( int k = 0; k < TriQuadraturePoints; k++ )
2776  {
2777  Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
2778  TriQuadratureG[k][2] * node2.x,
2779  TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
2780  TriQuadratureG[k][2] * node2.y,
2781  TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
2782  TriQuadratureG[k][2] * node2.z );
2783 
2784  double dMagnitude = node.Magnitude();
2785  node.x /= dMagnitude;
2786  node.y /= dMagnitude;
2787  node.z /= dMagnitude;
2788 
2789  double dLon = atan2( node.y, node.x );
2790  if( dLon < 0.0 )
2791  {
2792  dLon += 2.0 * M_PI;
2793  }
2794  double dLat = asin( node.z );
2795 
2796  double dSample = ( *testFunction )( dLon, dLat );
2797 
2798  dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
2799  }
2800 
2801  dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
2802  }
2803  }
2804  MB_CHK_ERR( m_interface->tag_set_data( solnTag, entities, &dVar[0] ) );
2805  }
2806  else /* discMethod == DiscretizationType_PCLOUD */
2807  {
2808  /* Get the coordinates of the vertices and sample the functionals accordingly */
2809  std::vector< Node >& nodes = trmesh->nodes;
2810 
2811  // Resize the array
2812  dVar.Allocate( nodes.size() );
2813 
2814  for( size_t j = 0; j < nodes.size(); j++ )
2815  {
2816  Node& node = nodes[j];
2817  double dMagnitude = node.Magnitude();
2818  node.x /= dMagnitude;
2819  node.y /= dMagnitude;
2820  node.z /= dMagnitude;
2821  double dLon = atan2( node.y, node.x );
2822  if( dLon < 0.0 )
2823  {
2824  dLon += 2.0 * M_PI;
2825  }
2826  double dLat = asin( node.z );
2827 
2828  double dSample = ( *testFunction )( dLon, dLat );
2829  dVar[j] = dSample;
2830  }
2831 
2832  MB_CHK_ERR( m_interface->tag_set_data( solnTag, entities, &dVar[0] ) );
2833  }
2834  }
2835 
2836  return moab::MB_SUCCESS;
2837 }
2838 
2839 moab::ErrorCode moab::TempestOnlineMap::ComputeMetrics( moab::Remapper::IntersectionContext ctx,
2840  moab::Tag& exactTag,
2841  moab::Tag& approxTag,
2842  std::map< std::string, double >& metrics,
2843  bool verbose )
2844 {
2845  const bool outputEnabled = ( is_root );
2846  int discOrder;
2847  // DiscretizationType discMethod;
2848  // moab::EntityHandle meshset;
2849  moab::Range entities;
2850  // Mesh* trmesh;
2851  switch( ctx )
2852  {
2853  case Remapper::SourceMesh:
2854  // meshset = m_remapper->m_covering_source_set;
2855  // trmesh = m_remapper->m_covering_source;
2856  entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
2857  : m_remapper->m_covering_source_entities );
2858  discOrder = m_nDofsPEl_Src;
2859  // discMethod = m_eInputType;
2860  break;
2861 
2862  case Remapper::TargetMesh:
2863  // meshset = m_remapper->m_target_set;
2864  // trmesh = m_remapper->m_target;
2865  entities =
2866  ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
2867  discOrder = m_nDofsPEl_Dest;
2868  // discMethod = m_eOutputType;
2869  break;
2870 
2871  default:
2872  if( outputEnabled )
2873  std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
2874  return moab::MB_FAILURE;
2875  }
2876 
2877  // Let us create teh solution tag with appropriate information for name, discretization order
2878  // (DoF space)
2879  std::string exactTagName, projTagName;
2880  const int ntotsize = entities.size() * discOrder * discOrder;
2881  std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
2882  MB_CHK_ERR( m_interface->tag_get_name( exactTag, exactTagName ) );
2883  MB_CHK_ERR( m_interface->tag_get_data( exactTag, entities, &exactSolution[0] ) );
2884  MB_CHK_ERR( m_interface->tag_get_name( approxTag, projTagName ) );
2885  MB_CHK_ERR( m_interface->tag_get_data( approxTag, entities, &projSolution[0] ) );
2886 
2887  const auto& ovents = m_remapper->m_overlap_entities;
2888 
2889  std::vector< double > errnorms( 4, 0.0 ), globerrnorms( 4, 0.0 ); // L1Err, L2Err, LinfErr
2890  double sumarea = 0.0;
2891  for( size_t i = 0; i < ovents.size(); ++i )
2892  {
2893  const int srcidx = m_remapper->m_overlap->vecSourceFaceIx[i];
2894  if( srcidx < 0 ) continue; // Skip non-overlapping entities
2895  const int tgtidx = m_remapper->m_overlap->vecTargetFaceIx[i];
2896  if( tgtidx < 0 ) continue; // skip ghost target faces
2897  const double ovarea = m_remapper->m_overlap->vecFaceArea[i];
2898  const double error = fabs( exactSolution[tgtidx] - projSolution[tgtidx] );
2899  errnorms[0] += ovarea * error;
2900  errnorms[1] += ovarea * error * error;
2901  errnorms[3] = ( error > errnorms[3] ? error : errnorms[3] );
2902  sumarea += ovarea;
2903  }
2904  errnorms[2] = sumarea;
2905 #ifdef MOAB_HAVE_MPI
2906  if( m_pcomm )
2907  {
2908  MPI_Reduce( &errnorms[0], &globerrnorms[0], 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
2909  MPI_Reduce( &errnorms[3], &globerrnorms[3], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
2910  }
2911 #else
2912  for( int i = 0; i < 4; ++i )
2913  globerrnorms[i] = errnorms[i];
2914 #endif
2915 
2916  globerrnorms[0] = ( globerrnorms[0] / globerrnorms[2] );
2917  globerrnorms[1] = std::sqrt( globerrnorms[1] / globerrnorms[2] );
2918 
2919  metrics.clear();
2920  metrics["L1Error"] = globerrnorms[0];
2921  metrics["L2Error"] = globerrnorms[1];
2922  metrics["LinfError"] = globerrnorms[3];
2923 
2924  if( verbose && is_root )
2925  {
2926  std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
2927  std::cout << "\t Total Intersection area = " << globerrnorms[2] << std::endl;
2928  std::cout << "\t L_1 error = " << globerrnorms[0] << std::endl;
2929  std::cout << "\t L_2 error = " << globerrnorms[1] << std::endl;
2930  std::cout << "\t L_inf error = " << globerrnorms[3] << std::endl;
2931  }
2932 
2933  return moab::MB_SUCCESS;
2934 }