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