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