 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
An offline map between two Meshes. More...
#include <TempestOnlineMap.hpp>
Inheritance diagram for moab::TempestOnlineMap: Collaboration diagram for moab::TempestOnlineMap:Public Types | |
| enum | DiscretizationType { DiscretizationType_FV = 0 , DiscretizationType_CGLL = 1 , DiscretizationType_DGLL = 2 , DiscretizationType_PCLOUD = 3 } |
| enum | CAASType { CAAS_NONE = 0 , CAAS_GLOBAL = 1 , CAAS_LOCAL = 2 , CAAS_LOCAL_ADJACENT = 3 , CAAS_QLT = 4 } |
| typedef double(* | sample_function) (double, double) |
Public Member Functions | |||||
| TempestOnlineMap (moab::TempestRemapper *remapper) | |||||
| Generate the metadata associated with the offline map. More... | |||||
| virtual | ~TempestOnlineMap () | ||||
| Define a virtual destructor. More... | |||||
| moab::ErrorCode | GenerateRemappingWeights (std::string strInputType, std::string strOutputType, const GenerateOfflineMapAlgorithmOptions &mapOptions, const std::string &srcDofTagName="GLOBAL_ID", const std::string &tgtDofTagName="GLOBAL_ID") | ||||
| Generate the offline map, given the source and target mesh and discretization details. This method generates the mapping between the two meshes based on the overlap and stores the result in the SparseMatrix. More... | |||||
| moab::ErrorCode | ReadParallelMap (const char *strSource, const std::vector< int > &tgt_dof_ids, std::vector< double > &areaA, int &nA, std::vector< double > &areaB, int &nB) | ||||
| Generate the metadata associated with the offline map. More... | |||||
| moab::ErrorCode | WriteParallelMap (const std::string &strTarget, const std::map< std::string, std::string > &attrMap) | ||||
| Write the TempestOnlineMap to a parallel NetCDF file. More... | |||||
| virtual int | IsConsistent (double dTolerance) | ||||
| Determine if the map is first-order accurate. More... | |||||
| virtual int | IsConservative (double dTolerance) | ||||
| Determine if the map is conservative. More... | |||||
| virtual int | IsMonotone (double dTolerance) | ||||
| Determine if the map is monotone. More... | |||||
| const DataArray1D< double > & | GetGlobalSourceAreas () const | ||||
| If we computed the reduction, get the vector representing the source areas for all entities in the mesh More... | |||||
| const DataArray1D< double > & | GetGlobalTargetAreas () const | ||||
| If we computed the reduction, get the vector representing the target areas for all entities in the mesh More... | |||||
| void | PrintMapStatistics () | ||||
| Print information and metadata about the remapping weights. More... | |||||
| moab::ErrorCode | SetDOFmapTags (const std::string srcDofTagName, const std::string tgtDofTagName) | ||||
Store the tag names associated with global DoF ids for source and target meshes to be used for mapping.
| |||||
| moab::ErrorCode | SetDOFmapAssociation (DiscretizationType srcType, int srcOrder, bool isSrcContinuous, DataArray3D< int > *srcdataGLLNodes, DataArray3D< int > *srcdataGLLNodesSrc, DiscretizationType destType, int destOrder, bool isTgtContinuous, DataArray3D< int > *tgtdataGLLNodes) | ||||
| std::pair< double, double > | ApplyBoundsLimiting (std::vector< double > &dataInDouble, std::vector< double > &dataOutDouble, CAASType caasType=CAAS_GLOBAL, int caasIteration=0, double mismatch=0.0) | ||||
| void | ComputeAdjacencyRelations (std::vector< std::unordered_set< int > > &vecAdjFaces, int nrings, const Range &entities, bool useMOABAdjacencies=true, Mesh *trMesh=nullptr) | ||||
| int | GetSourceGlobalNDofs () | ||||
| Get the number of total Degrees-Of-Freedom defined on the source mesh. More... | |||||
| int | GetDestinationGlobalNDofs () | ||||
| Get the number of total Degrees-Of-Freedom defined on the destination mesh. More... | |||||
| int | GetSourceLocalNDofs () | ||||
| Get the number of local Degrees-Of-Freedom defined on the source mesh. More... | |||||
| int | GetDestinationLocalNDofs () | ||||
| Get the number of local Degrees-Of-Freedom defined on the destination mesh. More... | |||||
| int | GetSourceNDofsPerElement () | ||||
| Get the number of Degrees-Of-Freedom per element on the source mesh. More... | |||||
| int | GetDestinationNDofsPerElement () | ||||
| Get the number of Degrees-Of-Freedom per element on the destination mesh. More... | |||||
| void | SetSourceNDofsPerElement (int ns) | ||||
| Set the number of Degrees-Of-Freedom per element on the source mesh. More... | |||||
| void | SetDestinationNDofsPerElement (int nt) | ||||
| Get the number of Degrees-Of-Freedom per element on the destination mesh. More... | |||||
| int | GetRowGlobalDoF (int localID) const | ||||
| Get the global Degrees-Of-Freedom ID on the destination mesh. More... | |||||
| int | GetIndexOfRowGlobalDoF (int globalRowDoF) const | ||||
| Get the index of globaRowDoF. More... | |||||
| int | GetColGlobalDoF (int localID) const | ||||
| Get the global Degrees-Of-Freedom ID on the source mesh. More... | |||||
| int | GetIndexOfColGlobalDoF (int globalColDoF) const | ||||
| Get the index of globaColDoF. More... | |||||
| moab::ErrorCode | ApplyWeights (moab::Tag srcSolutionTag, moab::Tag tgtSolutionTag, bool transpose=false, CAASType caasType=CAAS_NONE, double default_projection=0.0) | ||||
Apply the weight matrix onto the source vector (tag) provided as input, and return the column vector (solution projection) in a tag, after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals More... | |||||
| moab::ErrorCode | DefineAnalyticalSolution (moab::Tag &exactSolnTag, const std::string &solnName, Remapper::IntersectionContext ctx, sample_function testFunction, moab::Tag *clonedSolnTag=NULL, std::string cloneSolnName="") | ||||
| Define an analytical solution over the given (source or target) mesh, as specificed in the context. This routine will define a tag that is compatible with the specified discretization method type and order and sample the solution exactly using the analytical function provided by the user. More... | |||||
| moab::ErrorCode | ComputeMetrics (Remapper::IntersectionContext ctx, moab::Tag &exactTag, moab::Tag &approxTag, std::map< std::string, double > &metrics, bool verbose=true) | ||||
| Compute the error between a sampled (exact) solution and a projected solution in various error norms. More... | |||||
| moab::ErrorCode | fill_row_ids (std::vector< int > &ids_of_interest) | ||||
| moab::ErrorCode | fill_col_ids (std::vector< int > &ids_of_interest) | ||||
| moab::ErrorCode | set_col_dc_dofs (std::vector< int > &values_entities) | ||||
| moab::ErrorCode | set_row_dc_dofs (std::vector< int > &values_entities) | ||||
| void | SetMeshInput (Mesh *imesh) | ||||
Private Member Functions | |
| moab::ErrorCode | LinearRemapNN_MOAB (bool use_GID_matching=false, bool strict_check=false) |
| Compute the remapping weights as a permutation matrix that relates DoFs on the source mesh to DoFs on the target mesh. More... | |
| void | LinearRemapFVtoFV_Tempest_MOAB (int nOrder) |
| Compute the remapping weights for a FV field defined on the source to a FV field defined on the target mesh. More... | |
| void | LinearRemapSE0_Tempest_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian) |
| Generate the OfflineMap for linear conserative element-average spectral element to element average remapping. More... | |
| void | LinearRemapSE4_Tempest_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian, int nMonotoneType, bool fContinuousIn, bool fNoConservation) |
| Generate the OfflineMap for cubic conserative element-average spectral element to element average remapping. More... | |
| void | LinearRemapFVtoGLL_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian, const DataArray1D< double > &dataGLLNodalArea, int nOrder, int nMonotoneType, bool fContinuous, bool fNoConservation) |
| Generate the OfflineMap for remapping from finite volumes to finite elements. More... | |
| void | LinearRemapGLLtoGLL2_MOAB (const DataArray3D< int > &dataGLLNodesIn, const DataArray3D< double > &dataGLLJacobianIn, const DataArray3D< int > &dataGLLNodesOut, const DataArray3D< double > &dataGLLJacobianOut, const DataArray1D< double > &dataNodalAreaOut, int nPin, int nPout, int nMonotoneType, bool fContinuousIn, bool fContinuousOut, bool fNoConservation) |
| Generate the OfflineMap for remapping from finite elements to finite elements. More... | |
| void | LinearRemapGLLtoGLL2_Pointwise_MOAB (const DataArray3D< int > &dataGLLNodesIn, const DataArray3D< double > &dataGLLJacobianIn, const DataArray3D< int > &dataGLLNodesOut, const DataArray3D< double > &dataGLLJacobianOut, const DataArray1D< double > &dataNodalAreaOut, int nPin, int nPout, int nMonotoneType, bool fContinuousIn, bool fContinuousOut) |
| Generate the OfflineMap for remapping from finite elements to finite elements (pointwise interpolation). More... | |
| moab::ErrorCode | WriteSCRIPMapFile (const std::string &strOutputFile, const std::map< std::string, std::string > &attrMap) |
| Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections. More... | |
| moab::ErrorCode | WriteHDF5MapFile (const std::string &filename) |
| Parallel I/O with NetCDF to write out the SCRIP file from multiple processors. More... | |
| template<typename SparseMatrixType > | |
| void | serializeSparseMatrix (const SparseMatrixType &mat, const std::string &filename) |
| void | setup_sizes_dimensions () |
| void | CAASLimiter (std::vector< double > &dataCorrectedField, std::vector< double > &dataLowerBound, std::vector< double > &dataUpperBound, double &dMass) |
| double | QLTLimiter (int caasIteration, std::vector< double > &dataCorrectedField, std::vector< double > &dataLowerBound, std::vector< double > &dataUpperBound, std::vector< double > &dMassDefect) |
| moab::ErrorCode | ApplyWeights (std::vector< double > &srcVals, std::vector< double > &tgtVals, bool transpose=false) |
Apply the weight matrix onto the source vector provided as input, and return the column vector (solution projection) after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals More... | |
An offline map between two Meshes.
Definition at line 61 of file TempestOnlineMap.hpp.
| typedef double( * moab::TempestOnlineMap::sample_function) (double, double) |
Definition at line 414 of file TempestOnlineMap.hpp.
| Enumerator | |
|---|---|
| CAAS_NONE | |
| CAAS_GLOBAL | |
| CAAS_LOCAL | |
| CAAS_LOCAL_ADJACENT | |
| CAAS_QLT | |
Definition at line 86 of file TempestOnlineMap.hpp.
87 { 88 CAAS_NONE = 0, 89 CAAS_GLOBAL = 1, 90 CAAS_LOCAL = 2, 91 CAAS_LOCAL_ADJACENT = 3, 92 CAAS_QLT = 4 93 };
| Enumerator | |
|---|---|
| DiscretizationType_FV | |
| DiscretizationType_CGLL | |
| DiscretizationType_DGLL | |
| DiscretizationType_PCLOUD | |
Definition at line 77 of file TempestOnlineMap.hpp.
78 { 79 DiscretizationType_FV = 0, 80 DiscretizationType_CGLL = 1, 81 DiscretizationType_DGLL = 2, 82 DiscretizationType_PCLOUD = 3 83 };
| moab::TempestOnlineMap::TempestOnlineMap | ( | moab::TempestRemapper * | remapper | ) |
Generate the metadata associated with the offline map.
Definition at line 64 of file TempestOnlineMap.cpp.
64 : OfflineMap(), m_remapper( remapper )
65 {
66 // Get the references for the MOAB core objects
67 m_interface = m_remapper->get_interface();
68 #ifdef MOAB_HAVE_MPI
69 m_pcomm = m_remapper->get_parallel_communicator();
70 #endif
71
72 // now let us re-update the reference to the input source mesh
73 m_meshInput = m_remapper->GetMesh( moab::Remapper::SourceMesh );
74 // now let us re-update the reference to the covering mesh
75 m_meshInputCov = m_remapper->GetCoveringMesh();
76 // now let us re-update the reference to the output target mesh
77 m_meshOutput = m_remapper->GetMesh( moab::Remapper::TargetMesh );
78 // now let us re-update the reference to the output target mesh
79 m_meshOverlap = m_remapper->GetMesh( moab::Remapper::OverlapMesh );
80
81 is_parallel = remapper->is_parallel;
82 is_root = remapper->is_root;
83 rank = remapper->rank;
84 size = remapper->size;
85
86 // set default order
87 m_input_order = m_output_order = 1;
88
89 // Initialize dimension information from file
90 this->setup_sizes_dimensions();
91 }
References moab::Remapper::get_interface(), moab::TempestRemapper::GetCoveringMesh(), moab::TempestRemapper::GetMesh(), is_parallel, moab::TempestRemapper::is_parallel, is_root, moab::TempestRemapper::is_root, m_input_order, m_interface, m_meshInput, m_meshInputCov, m_meshOutput, m_meshOverlap, m_output_order, m_remapper, moab::Remapper::OverlapMesh, rank, moab::TempestRemapper::rank, setup_sizes_dimensions(), size, moab::TempestRemapper::size, moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.
|
virtual |
Define a virtual destructor.
Definition at line 118 of file TempestOnlineMap.cpp.
119 {
120 m_interface = nullptr;
121 #ifdef MOAB_HAVE_MPI
122 m_pcomm = nullptr;
123 #endif
124 m_meshInput = nullptr;
125 m_meshOutput = nullptr;
126 m_meshOverlap = nullptr;
127 }
| std::pair< double, double > moab::TempestOnlineMap::ApplyBoundsLimiting | ( | std::vector< double > & | dataInDouble, |
| std::vector< double > & | dataOutDouble, | ||
| CAASType | caasType = CAAS_GLOBAL, |
||
| int | caasIteration = 0, |
||
| double | mismatch = 0.0 |
||
| ) |
ApplyBoundsLimiting - Apply bounds limiting to the data field
| dataInDouble | - input data field |
| dataOutDouble | - output data field |
| caasType | - type of limiter |
| caasIteration | - iteration number of limiter |
| moab::ErrorCode moab::TempestOnlineMap::ApplyWeights | ( | moab::Tag | srcSolutionTag, |
| moab::Tag | tgtSolutionTag, | ||
| bool | transpose = false, |
||
| CAASType | caasType = CAAS_NONE, |
||
| double | default_projection = 0.0 |
||
| ) |
Apply the weight matrix onto the source vector (tag) provided as input, and return the column vector (solution projection) in a tag, after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals
Definition at line 1674 of file TempestOnlineMap.cpp.
1679 {
1680 std::vector< double > solSTagVals;
1681 std::vector< double > solTTagVals;
1682
1683 moab::Range sents, tents;
1684 if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
1685 {
1686 if( m_remapper->point_cloud_source )
1687 {
1688 moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh );
1689 solSTagVals.resize( covSrcEnts.size(), default_projection );
1690 sents = covSrcEnts;
1691 }
1692 else
1693 {
1694 moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
1695 solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1696 default_projection );
1697 sents = covSrcEnts;
1698 }
1699 if( m_remapper->point_cloud_target )
1700 {
1701 moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
1702 solTTagVals.resize( tgtEnts.size(), default_projection );
1703 tents = tgtEnts;
1704 }
1705 else
1706 {
1707 moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
1708 solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
1709 this->GetDestinationNDofsPerElement(),
1710 default_projection );
1711 tents = tgtEnts;
1712 }
1713 }
1714 else
1715 {
1716 moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
1717 moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
1718 solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1719 default_projection );
1720 solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
1721 this->GetDestinationNDofsPerElement(),
1722 default_projection );
1723
1724 sents = covSrcEnts;
1725 tents = tgtEnts;
1726 }
1727
1728 // The tag data is np*np*n_el_src
1729 MB_CHK_SET_ERR( m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] ),
1730 "Getting local tag data failed" );
1731
1732 // Compute the application of weights on the suorce solution data and store it in the
1733 // destination solution vector data Optionally, can also perform the transpose application of
1734 // the weight matrix. Set the 3rd argument to true if this is needed
1735 MB_CHK_SET_ERR( this->ApplyWeights( solSTagVals, solTTagVals, transpose ),
1736 "Applying remap operator onto source vector data failed" );
1737
1738 // The tag data is np*np*n_el_dest
1739 MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1740 "Setting target tag data failed" );
1741
1742 if( caasType != CAAS_NONE )
1743 {
1744 std::string tgtSolutionTagName;
1745 MB_CHK_SET_ERR( m_interface->tag_get_name( tgtSolutionTag, tgtSolutionTagName ), "Getting tag name failed" );
1746
1747 // Perform CAAS iterations iteratively until convergence
1748 constexpr int nmax_caas_iterations = 10;
1749 double mismatch = 1.0;
1750 int caasIteration = 0;
1751 double initialMismatch = 0.0;
1752 while( ( fabs( mismatch / initialMismatch ) > 1e-15 && fabs( mismatch ) > 1e-15 ) &&
1753 caasIteration++ < nmax_caas_iterations ) // iterate until convergence or a maximum of 5 iterations
1754 {
1755 double dMassDiffPostGlobal;
1756 std::pair< double, double > mDefect =
1757 this->ApplyBoundsLimiting( solSTagVals, solTTagVals, caasType, caasIteration, mismatch );
1758 #ifdef MOAB_HAVE_MPI
1759 double dMassDiffPost = mDefect.second;
1760 MPI_Allreduce( &dMassDiffPost, &dMassDiffPostGlobal, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
1761 #else
1762 dMassDiffPostGlobal = mDefect.second;
1763 #endif
1764 if( caasIteration == 1 ) initialMismatch = mDefect.first;
1765 if( m_remapper->verbose && is_root )
1766 {
1767 printf( "Field {%s} -> CAAS iteration: %d, mass defect: %3.4e, post-CAAS: %3.4e\n",
1768 tgtSolutionTagName.c_str(), caasIteration, mDefect.first, dMassDiffPostGlobal );
1769 }
1770 mismatch = dMassDiffPostGlobal;
1771
1772 // The tag data is np*np*n_el_dest
1773 MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1774 "Setting local tag data failed" );
1775 }
1776 }
1777
1778 return moab::MB_SUCCESS;
1779 }
References moab::Remapper::CoveringMesh, MB_CHK_SET_ERR, MB_SUCCESS, moab::Range::size(), and moab::Remapper::TargetMesh.
Referenced by main().
|
private |
Apply the weight matrix onto the source vector provided as input, and return the column vector (solution projection) after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals
|
private |
| void moab::TempestOnlineMap::ComputeAdjacencyRelations | ( | std::vector< std::unordered_set< int > > & | vecAdjFaces, |
| int | nrings, | ||
| const Range & | entities, | ||
| bool | useMOABAdjacencies = true, |
||
| Mesh * | trMesh = nullptr |
||
| ) |
| vecAdjFaces | |
| nrings | |
| entities | |
| useMOABAdjacencies | |
| trMesh |
Vector storing adjacent Faces.
Definition at line 1625 of file TempestOnlineMap.cpp.
1630 {
1631 assert( nrings > 0 );
1632 assert( useMOABAdjacencies || trMesh != nullptr );
1633
1634 const size_t nrows = vecAdjFaces.size();
1635 moab::MeshTopoUtil mtu( m_interface );
1636 for( size_t index = 0; index < nrows; index++ )
1637 {
1638 vecAdjFaces[index].insert( index ); // add self target face first
1639 {
1640 // Compute the adjacent faces to the target face
1641 if( useMOABAdjacencies )
1642 {
1643 moab::Range ents;
1644 // ents.insert( entities.index( entities[index] ) );
1645 ents.insert( entities[index] );
1646 moab::Range adjEnts;
1647 moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, nrings );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
1648 for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
1649 {
1650 // int adjIndex = m_interface->id_from_handle(*it)-1;
1651 int adjIndex = entities.index( *it );
1652 // printf("rank: %d, Element %lu, entity: %lu, adjIndex %d\n", rank, index, *it, adjIndex);
1653 if( adjIndex >= 0 ) vecAdjFaces[index].insert( adjIndex );
1654 }
1655 }
1656 else
1657 {
1658 /// Vector storing adjacent Faces.
1659 typedef std::pair< int, int > FaceDistancePair;
1660 typedef std::vector< FaceDistancePair > AdjacentFaceVector;
1661 AdjacentFaceVector adjFaces;
1662 Face& face = trMesh->faces[index];
1663 GetAdjacentFaceVectorByEdge( *trMesh, index, nrings * face.edges.size(), adjFaces );
1664
1665 // Add the adjacent faces to the target face list
1666 for( auto adjFace : adjFaces )
1667 if( adjFace.first >= 0 )
1668 vecAdjFaces[index].insert( adjFace.first ); // map target face to source face
1669 }
1670 }
1671 }
1672 }
References moab::Range::begin(), moab::Range::end(), entities, ErrorCode, moab::MeshTopoUtil::get_bridge_adjacencies(), moab::Range::insert(), and MB_CHK_SET_ERR_CONT.
| moab::ErrorCode moab::TempestOnlineMap::ComputeMetrics | ( | Remapper::IntersectionContext | ctx, |
| moab::Tag & | exactTag, | ||
| moab::Tag & | approxTag, | ||
| std::map< std::string, double > & | metrics, | ||
| bool | verbose = true |
||
| ) |
Compute the error between a sampled (exact) solution and a projected solution in various error norms.
Definition at line 2143 of file TempestOnlineMap.cpp.
2148 {
2149 moab::ErrorCode rval;
2150 const bool outputEnabled = ( is_root );
2151 int discOrder;
2152 // DiscretizationType discMethod;
2153 // moab::EntityHandle meshset;
2154 moab::Range entities;
2155 // Mesh* trmesh;
2156 switch( ctx )
2157 {
2158 case Remapper::SourceMesh:
2159 // meshset = m_remapper->m_covering_source_set;
2160 // trmesh = m_remapper->m_covering_source;
2161 entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
2162 : m_remapper->m_covering_source_entities );
2163 discOrder = m_nDofsPEl_Src;
2164 // discMethod = m_eInputType;
2165 break;
2166
2167 case Remapper::TargetMesh:
2168 // meshset = m_remapper->m_target_set;
2169 // trmesh = m_remapper->m_target;
2170 entities =
2171 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
2172 discOrder = m_nDofsPEl_Dest;
2173 // discMethod = m_eOutputType;
2174 break;
2175
2176 default:
2177 if( outputEnabled )
2178 std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
2179 return moab::MB_FAILURE;
2180 }
2181
2182 // Let us create teh solution tag with appropriate information for name, discretization order
2183 // (DoF space)
2184 std::string exactTagName, projTagName;
2185 const int ntotsize = entities.size() * discOrder * discOrder;
2186 std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
2187 rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval );
2188 rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval );
2189 rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval );
2190 rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval );
2191
2192 const auto& ovents = m_remapper->m_overlap_entities;
2193
2194 std::vector< double > errnorms( 4, 0.0 ), globerrnorms( 4, 0.0 ); // L1Err, L2Err, LinfErr
2195 double sumarea = 0.0;
2196 for( size_t i = 0; i < ovents.size(); ++i )
2197 {
2198 const int srcidx = m_remapper->m_overlap->vecSourceFaceIx[i];
2199 if( srcidx < 0 ) continue; // Skip non-overlapping entities
2200 const int tgtidx = m_remapper->m_overlap->vecTargetFaceIx[i];
2201 if( tgtidx < 0 ) continue; // skip ghost target faces
2202 const double ovarea = m_remapper->m_overlap->vecFaceArea[i];
2203 const double error = fabs( exactSolution[tgtidx] - projSolution[tgtidx] );
2204 errnorms[0] += ovarea * error;
2205 errnorms[1] += ovarea * error * error;
2206 errnorms[3] = ( error > errnorms[3] ? error : errnorms[3] );
2207 sumarea += ovarea;
2208 }
2209 errnorms[2] = sumarea;
2210 #ifdef MOAB_HAVE_MPI
2211 if( m_pcomm )
2212 {
2213 MPI_Reduce( &errnorms[0], &globerrnorms[0], 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
2214 MPI_Reduce( &errnorms[3], &globerrnorms[3], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
2215 }
2216 #else
2217 for( int i = 0; i < 4; ++i )
2218 globerrnorms[i] = errnorms[i];
2219 #endif
2220
2221 globerrnorms[0] = ( globerrnorms[0] / globerrnorms[2] );
2222 globerrnorms[1] = std::sqrt( globerrnorms[1] / globerrnorms[2] );
2223
2224 metrics.clear();
2225 metrics["L1Error"] = globerrnorms[0];
2226 metrics["L2Error"] = globerrnorms[1];
2227 metrics["LinfError"] = globerrnorms[3];
2228
2229 if( verbose && is_root )
2230 {
2231 std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
2232 std::cout << "\t Total Intersection area = " << globerrnorms[2] << std::endl;
2233 std::cout << "\t L_1 error = " << globerrnorms[0] << std::endl;
2234 std::cout << "\t L_2 error = " << globerrnorms[1] << std::endl;
2235 std::cout << "\t L_inf error = " << globerrnorms[3] << std::endl;
2236 }
2237
2238 return moab::MB_SUCCESS;
2239 }
References entities, moab::error(), ErrorCode, MB_CHK_ERR, MB_SUCCESS, moab::Remapper::SourceMesh, moab::Remapper::TargetMesh, and verbose.
Referenced by main().
| moab::ErrorCode moab::TempestOnlineMap::DefineAnalyticalSolution | ( | moab::Tag & | exactSolnTag, |
| const std::string & | solnName, | ||
| Remapper::IntersectionContext | ctx, | ||
| sample_function | testFunction, | ||
| moab::Tag * | clonedSolnTag = NULL, |
||
| std::string | cloneSolnName = "" |
||
| ) |
Define an analytical solution over the given (source or target) mesh, as specificed in the context. This routine will define a tag that is compatible with the specified discretization method type and order and sample the solution exactly using the analytical function provided by the user.
Definition at line 1781 of file TempestOnlineMap.cpp.
1787 {
1788 moab::ErrorCode rval;
1789 const bool outputEnabled = ( is_root );
1790 int discOrder;
1791 DiscretizationType discMethod;
1792 // moab::EntityHandle meshset;
1793 moab::Range entities;
1794 Mesh* trmesh;
1795 switch( ctx )
1796 {
1797 case Remapper::SourceMesh:
1798 // meshset = m_remapper->m_covering_source_set;
1799 trmesh = m_remapper->m_covering_source;
1800 entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
1801 : m_remapper->m_covering_source_entities );
1802 discOrder = m_nDofsPEl_Src;
1803 discMethod = m_eInputType;
1804 break;
1805
1806 case Remapper::TargetMesh:
1807 // meshset = m_remapper->m_target_set;
1808 trmesh = m_remapper->m_target;
1809 entities =
1810 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
1811 discOrder = m_nDofsPEl_Dest;
1812 discMethod = m_eOutputType;
1813 break;
1814
1815 default:
1816 if( outputEnabled )
1817 std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
1818 return moab::MB_FAILURE;
1819 }
1820
1821 // Let us create teh solution tag with appropriate information for name, discretization order
1822 // (DoF space)
1823 rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
1824 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
1825 if( clonedSolnTag != nullptr )
1826 {
1827 if( cloneSolnName.size() == 0 )
1828 {
1829 cloneSolnName = solnName + std::string( "Cloned" );
1830 }
1831 rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
1832 *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
1833 }
1834
1835 // Triangular quadrature rule
1836 const int TriQuadratureOrder = 10;
1837
1838 if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
1839
1840 TriangularQuadratureRule triquadrule( TriQuadratureOrder );
1841
1842 const int TriQuadraturePoints = triquadrule.GetPoints();
1843
1844 const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
1845 const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
1846
1847 // Output data
1848 DataArray1D< double > dVar;
1849 DataArray1D< double > dVarMB; // re-arranged local MOAB vector
1850
1851 // Nodal geometric area
1852 DataArray1D< double > dNodeArea;
1853
1854 // Calculate element areas
1855 // trmesh->CalculateFaceAreas(fContainsConcaveFaces);
1856
1857 if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
1858 {
1859 /* Get the spectral points and sample the functionals accordingly */
1860 const bool fGLL = true;
1861 const bool fGLLIntegrate = false;
1862
1863 // Generate grid metadata
1864 DataArray3D< int > dataGLLNodes;
1865 DataArray3D< double > dataGLLJacobian;
1866
1867 GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );
1868
1869 // Number of elements
1870 int nElements = trmesh->faces.size();
1871
1872 // Verify all elements are quadrilaterals
1873 for( int k = 0; k < nElements; k++ )
1874 {
1875 const Face& face = trmesh->faces[k];
1876
1877 if( face.edges.size() != 4 )
1878 {
1879 _EXCEPTIONT( "Non-quadrilateral face detected; "
1880 "incompatible with --gll" );
1881 }
1882 }
1883
1884 // Number of unique nodes
1885 int iMaxNode = 0;
1886 for( int i = 0; i < discOrder; i++ )
1887 {
1888 for( int j = 0; j < discOrder; j++ )
1889 {
1890 for( int k = 0; k < nElements; k++ )
1891 {
1892 if( dataGLLNodes[i][j][k] > iMaxNode )
1893 {
1894 iMaxNode = dataGLLNodes[i][j][k];
1895 }
1896 }
1897 }
1898 }
1899
1900 // Get Gauss-Lobatto quadrature nodes
1901 DataArray1D< double > dG;
1902 DataArray1D< double > dW;
1903
1904 GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
1905
1906 // Get Gauss quadrature nodes
1907 const int nGaussP = 10;
1908
1909 DataArray1D< double > dGaussG;
1910 DataArray1D< double > dGaussW;
1911
1912 GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
1913
1914 // Allocate data
1915 dVar.Allocate( iMaxNode );
1916 dVarMB.Allocate( discOrder * discOrder * nElements );
1917 dNodeArea.Allocate( iMaxNode );
1918
1919 // Sample data
1920 for( int k = 0; k < nElements; k++ )
1921 {
1922 const Face& face = trmesh->faces[k];
1923
1924 // Sample data at GLL nodes
1925 if( fGLL )
1926 {
1927 for( int i = 0; i < discOrder; i++ )
1928 {
1929 for( int j = 0; j < discOrder; j++ )
1930 {
1931
1932 // Apply local map
1933 Node node;
1934 Node dDx1G;
1935 Node dDx2G;
1936
1937 ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
1938
1939 // Sample data at this point
1940 double dNodeLon = atan2( node.y, node.x );
1941 if( dNodeLon < 0.0 )
1942 {
1943 dNodeLon += 2.0 * M_PI;
1944 }
1945 double dNodeLat = asin( node.z );
1946
1947 double dSample = ( *testFunction )( dNodeLon, dNodeLat );
1948
1949 dVar[dataGLLNodes[j][i][k] - 1] = dSample;
1950 }
1951 }
1952 // High-order Gaussian integration over basis function
1953 }
1954 else
1955 {
1956 DataArray2D< double > dCoeff( discOrder, discOrder );
1957
1958 for( int p = 0; p < nGaussP; p++ )
1959 {
1960 for( int q = 0; q < nGaussP; q++ )
1961 {
1962
1963 // Apply local map
1964 Node node;
1965 Node dDx1G;
1966 Node dDx2G;
1967
1968 ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
1969
1970 // Cross product gives local Jacobian
1971 Node nodeCross = CrossProduct( dDx1G, dDx2G );
1972
1973 double dJacobian =
1974 sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
1975
1976 // Find components of quadrature point in basis
1977 // of the first Face
1978 SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
1979
1980 // Sample data at this point
1981 double dNodeLon = atan2( node.y, node.x );
1982 if( dNodeLon < 0.0 )
1983 {
1984 dNodeLon += 2.0 * M_PI;
1985 }
1986 double dNodeLat = asin( node.z );
1987
1988 double dSample = ( *testFunction )( dNodeLon, dNodeLat );
1989
1990 // Integrate
1991 for( int i = 0; i < discOrder; i++ )
1992 {
1993 for( int j = 0; j < discOrder; j++ )
1994 {
1995
1996 double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
1997
1998 dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
1999
2000 dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
2001 }
2002 }
2003 }
2004 }
2005 }
2006 }
2007
2008 // Divide by area
2009 if( fGLLIntegrate )
2010 {
2011 for( size_t i = 0; i < dVar.GetRows(); i++ )
2012 {
2013 dVar[i] /= dNodeArea[i];
2014 }
2015 }
2016
2017 // Let us rearrange the data based on DoF ID specification
2018 if( ctx == Remapper::SourceMesh )
2019 {
2020 for( unsigned j = 0; j < entities.size(); j++ )
2021 for( int p = 0; p < discOrder; p++ )
2022 for( int q = 0; q < discOrder; q++ )
2023 {
2024 const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
2025 dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]];
2026 }
2027 }
2028 else
2029 {
2030 for( unsigned j = 0; j < entities.size(); j++ )
2031 for( int p = 0; p < discOrder; p++ )
2032 for( int q = 0; q < discOrder; q++ )
2033 {
2034 const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
2035 dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]];
2036 }
2037 }
2038
2039 // Set the tag data
2040 rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval );
2041 }
2042 else
2043 {
2044 // assert( discOrder == 1 );
2045 if( discMethod == DiscretizationType_FV )
2046 {
2047 /* Compute an element-wise integral to store the sampled solution based on Quadrature
2048 * rules */
2049 // Resize the array
2050 dVar.Allocate( trmesh->faces.size() );
2051
2052 std::vector< Node >& nodes = trmesh->nodes;
2053
2054 // Loop through all Faces
2055 for( size_t i = 0; i < trmesh->faces.size(); i++ )
2056 {
2057 const Face& face = trmesh->faces[i];
2058
2059 // Loop through all sub-triangles
2060 for( size_t j = 0; j < face.edges.size() - 2; j++ )
2061 {
2062
2063 const Node& node0 = nodes[face[0]];
2064 const Node& node1 = nodes[face[j + 1]];
2065 const Node& node2 = nodes[face[j + 2]];
2066
2067 // Triangle area
2068 Face faceTri( 3 );
2069 faceTri.SetNode( 0, face[0] );
2070 faceTri.SetNode( 1, face[j + 1] );
2071 faceTri.SetNode( 2, face[j + 2] );
2072
2073 double dTriangleArea = CalculateFaceArea( faceTri, nodes );
2074
2075 // Calculate the element average
2076 double dTotalSample = 0.0;
2077
2078 // Loop through all quadrature points
2079 for( int k = 0; k < TriQuadraturePoints; k++ )
2080 {
2081 Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
2082 TriQuadratureG[k][2] * node2.x,
2083 TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
2084 TriQuadratureG[k][2] * node2.y,
2085 TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
2086 TriQuadratureG[k][2] * node2.z );
2087
2088 double dMagnitude = node.Magnitude();
2089 node.x /= dMagnitude;
2090 node.y /= dMagnitude;
2091 node.z /= dMagnitude;
2092
2093 double dLon = atan2( node.y, node.x );
2094 if( dLon < 0.0 )
2095 {
2096 dLon += 2.0 * M_PI;
2097 }
2098 double dLat = asin( node.z );
2099
2100 double dSample = ( *testFunction )( dLon, dLat );
2101
2102 dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
2103 }
2104
2105 dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
2106 }
2107 }
2108 rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
2109 }
2110 else /* discMethod == DiscretizationType_PCLOUD */
2111 {
2112 /* Get the coordinates of the vertices and sample the functionals accordingly */
2113 std::vector< Node >& nodes = trmesh->nodes;
2114
2115 // Resize the array
2116 dVar.Allocate( nodes.size() );
2117
2118 for( size_t j = 0; j < nodes.size(); j++ )
2119 {
2120 Node& node = nodes[j];
2121 double dMagnitude = node.Magnitude();
2122 node.x /= dMagnitude;
2123 node.y /= dMagnitude;
2124 node.z /= dMagnitude;
2125 double dLon = atan2( node.y, node.x );
2126 if( dLon < 0.0 )
2127 {
2128 dLon += 2.0 * M_PI;
2129 }
2130 double dLat = asin( node.z );
2131
2132 double dSample = ( *testFunction )( dLon, dLat );
2133 dVar[j] = dSample;
2134 }
2135
2136 rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
2137 }
2138 }
2139
2140 return moab::MB_SUCCESS;
2141 }
References entities, ErrorCode, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.
Referenced by main().
|
inline |
Definition at line 449 of file TempestOnlineMap.hpp.
450 {
451 ids_of_interest.reserve( col_gdofmap.size() );
452 // need to add 1
453 for( auto it = col_gdofmap.begin(); it != col_gdofmap.end(); it++ )
454 ids_of_interest.push_back( *it + 1 );
455 return moab::MB_SUCCESS;
456 }
References col_gdofmap, and MB_SUCCESS.
|
inline |
Definition at line 439 of file TempestOnlineMap.hpp.
440 {
441 ids_of_interest.reserve( row_gdofmap.size() );
442 // need to add 1
443 for( auto it = row_gdofmap.begin(); it != row_gdofmap.end(); it++ )
444 ids_of_interest.push_back( *it + 1 );
445
446 return moab::MB_SUCCESS;
447 }
References MB_SUCCESS, and row_gdofmap.
| moab::ErrorCode moab::TempestOnlineMap::GenerateRemappingWeights | ( | std::string | strInputType, |
| std::string | strOutputType, | ||
| const GenerateOfflineMapAlgorithmOptions & | mapOptions, | ||
| const std::string & | srcDofTagName = "GLOBAL_ID", |
||
| const std::string & | tgtDofTagName = "GLOBAL_ID" |
||
| ) |
Generate the offline map, given the source and target mesh and discretization details. This method generates the mapping between the two meshes based on the overlap and stores the result in the SparseMatrix.
the tag should be created already in the e3sm workflow; if not, create it here
Definition at line 690 of file TempestOnlineMap.cpp.
695 {
696 NcError error( NcError::silent_nonfatal );
697
698 moab::DebugOutput dbgprint( std::cout, rank, 0 );
699 dbgprint.set_prefix( "[TempestOnlineMap]: " );
700 moab::ErrorCode rval;
701
702 const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
703 const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
704 const bool m_bPointCloud = m_bPointCloudSource || m_bPointCloudTarget;
705
706 // Build a matrix of source and target discretization so that we know how
707 // to assign the global DoFs in parallel for the mapping weights.
708 // For example,
709 // for FV->FV: the rows represented target DoFs and cols represent source DoFs
710 try
711 {
712 // Check command line parameters (data type arguments)
713 STLStringHelper::ToLower( strInputType );
714 STLStringHelper::ToLower( strOutputType );
715
716 DiscretizationType eInputType;
717 DiscretizationType eOutputType;
718
719 if( strInputType == "fv" )
720 {
721 eInputType = DiscretizationType_FV;
722 }
723 else if( strInputType == "cgll" )
724 {
725 eInputType = DiscretizationType_CGLL;
726 }
727 else if( strInputType == "dgll" )
728 {
729 eInputType = DiscretizationType_DGLL;
730 }
731 else if( strInputType == "pcloud" )
732 {
733 eInputType = DiscretizationType_PCLOUD;
734 }
735 else
736 {
737 _EXCEPTION1( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", strInputType.c_str() );
738 }
739
740 if( strOutputType == "fv" )
741 {
742 eOutputType = DiscretizationType_FV;
743 }
744 else if( strOutputType == "cgll" )
745 {
746 eOutputType = DiscretizationType_CGLL;
747 }
748 else if( strOutputType == "dgll" )
749 {
750 eOutputType = DiscretizationType_DGLL;
751 }
752 else if( strOutputType == "pcloud" )
753 {
754 eOutputType = DiscretizationType_PCLOUD;
755 }
756 else
757 {
758 _EXCEPTION1( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", strOutputType.c_str() );
759 }
760
761 // set all required input params
762 m_bConserved = !mapOptions.fNoConservation;
763 m_eInputType = eInputType;
764 m_eOutputType = eOutputType;
765
766 // Method flags
767 std::string strMapAlgorithm( "" );
768 int nMonotoneType = ( mapOptions.fMonotone ) ? ( 1 ) : ( 0 );
769
770 // Make an index of method arguments
771 std::set< std::string > setMethodStrings;
772 {
773 int iLast = 0;
774 for( size_t i = 0; i <= mapOptions.strMethod.length(); i++ )
775 {
776 if( ( i == mapOptions.strMethod.length() ) || ( mapOptions.strMethod[i] == ';' ) )
777 {
778 std::string strMethodString = mapOptions.strMethod.substr( iLast, i - iLast );
779 STLStringHelper::RemoveWhitespaceInPlace( strMethodString );
780 if( strMethodString.length() > 0 )
781 {
782 setMethodStrings.insert( strMethodString );
783 }
784 iLast = i + 1;
785 }
786 }
787 }
788
789 for( auto it : setMethodStrings )
790 {
791 // Piecewise constant monotonicity
792 if( it == "mono2" )
793 {
794 if( nMonotoneType != 0 )
795 {
796 _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
797 }
798 if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
799 {
800 _EXCEPTIONT( "--method \"mono2\" is only used when remapping to/from CGLL or DGLL grids" );
801 }
802 nMonotoneType = 2;
803
804 // Piecewise linear monotonicity
805 }
806 else if( it == "mono3" )
807 {
808 if( nMonotoneType != 0 )
809 {
810 _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
811 }
812 if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
813 {
814 _EXCEPTIONT( "--method \"mono3\" is only used when remapping to/from CGLL or DGLL grids" );
815 }
816 nMonotoneType = 3;
817
818 // Volumetric remapping from FV to GLL
819 }
820 else if( it == "volumetric" )
821 {
822 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType == DiscretizationType_FV ) )
823 {
824 _EXCEPTIONT( "--method \"volumetric\" may only be used for FV->CGLL or FV->DGLL remapping" );
825 }
826 strMapAlgorithm = "volumetric";
827
828 // Inverse distance mapping
829 }
830 else if( it == "invdist" )
831 {
832 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
833 {
834 _EXCEPTIONT( "--method \"invdist\" may only be used for FV->FV remapping" );
835 }
836 strMapAlgorithm = "invdist";
837
838 // Delaunay triangulation mapping
839 }
840 else if( it == "delaunay" )
841 {
842 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
843 {
844 _EXCEPTIONT( "--method \"delaunay\" may only be used for FV->FV remapping" );
845 }
846 strMapAlgorithm = "delaunay";
847
848 // Bilinear
849 }
850 else if( it == "bilin" )
851 {
852 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
853 {
854 _EXCEPTIONT( "--method \"bilin\" may only be used for FV->FV remapping" );
855 }
856 strMapAlgorithm = "fvbilin";
857
858 // Integrated bilinear (same as mono3 when source grid is CGLL/DGLL)
859 }
860 else if( it == "intbilin" )
861 {
862 if( m_eOutputType != DiscretizationType_FV )
863 {
864 _EXCEPTIONT( "--method \"intbilin\" may only be used when mapping to FV." );
865 }
866 if( m_eInputType == DiscretizationType_FV )
867 {
868 strMapAlgorithm = "fvintbilin";
869 }
870 else
871 {
872 strMapAlgorithm = "mono3";
873 }
874
875 // Integrated bilinear with generalized Barycentric coordinates
876 }
877 else if( it == "intbilingb" )
878 {
879 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
880 {
881 _EXCEPTIONT( "--method \"intbilingb\" may only be used for FV->FV remapping" );
882 }
883 strMapAlgorithm = "fvintbilingb";
884 }
885 else
886 {
887 _EXCEPTION1( "Invalid --method argument \"%s\"", it.c_str() );
888 }
889 }
890
891 m_nDofsPEl_Src =
892 ( m_eInputType == DiscretizationType_FV || m_eInputType == DiscretizationType_PCLOUD ? 1
893 : mapOptions.nPin );
894 m_nDofsPEl_Dest =
895 ( m_eOutputType == DiscretizationType_FV || m_eOutputType == DiscretizationType_PCLOUD ? 1
896 : mapOptions.nPout );
897
898 // Set the source and target mesh objects
899 MB_CHK_ERR( SetDOFmapTags( srcDofTagName, tgtDofTagName ) );
900
901 /// the tag should be created already in the e3sm workflow; if not, create it here
902 Tag areaTag;
903 rval = m_interface->tag_get_handle( "aream", 1, MB_TYPE_DOUBLE, areaTag,
904 MB_TAG_DENSE | MB_TAG_EXCL | MB_TAG_CREAT );
905 if( MB_ALREADY_ALLOCATED == rval )
906 {
907 if( is_root ) dbgprint.printf( 0, "aream tag already defined \n" );
908 }
909
910 double local_areas[3] = { 0.0, 0.0, 0.0 }, global_areas[3] = { 0.0, 0.0, 0.0 };
911 if( !m_bPointCloudSource )
912 {
913 // Calculate Input Mesh Face areas
914 if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
915 local_areas[0] = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
916 // Set source element areas as tag on the source mesh
917 MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea ) );
918
919 // Update coverage source mesh areas as well.
920 m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
921 }
922
923 if( !m_bPointCloudTarget )
924 {
925 // Calculate Output Mesh Face areas
926 if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
927 local_areas[1] = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
928 // Set target element areas as tag on the target mesh
929 MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea ) );
930 }
931
932 if( !m_bPointCloud )
933 {
934 // Verify that overlap mesh is in the correct order (sanity check)
935 assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->vecTargetFaceIx.size() );
936
937 // Calculate Face areas
938 if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
939 local_areas[2] =
940 m_meshOverlap->CalculateFaceAreas( mapOptions.fSourceConcave || mapOptions.fTargetConcave );
941
942 // store it as global output for now - used later in reduction
943 std::copy( local_areas, local_areas + 3, global_areas );
944 #ifdef MOAB_HAVE_MPI
945 // reduce the local source, target and overlap mesh areas to global areas
946 if( m_pcomm && is_parallel )
947 MPI_Reduce( local_areas, global_areas, 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
948 #endif
949 if( is_root )
950 {
951 dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", global_areas[0] );
952 dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", global_areas[1] );
953 dbgprint.printf( 0, "Overlap Mesh Recovered Area: %1.15e\n", global_areas[2] );
954 }
955
956 // Correct areas to match the areas calculated in the overlap mesh
957 constexpr bool fCorrectAreas = true;
958 if( fCorrectAreas ) // In MOAB-TempestRemap, we will always keep this to be true
959 {
960 if( is_root ) dbgprint.printf( 0, "Correcting source/target areas to overlap mesh areas\n" );
961 DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() );
962 DataArray1D< double > dTargetArea( m_meshOutput->faces.size() );
963
964 assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() );
965 assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() );
966 assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() );
967
968 assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() );
969 assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() );
970
971 for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
972 {
973 if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 )
974 continue; // skip this cell since it is ghosted
975
976 // let us recompute the source/target areas based on overlap mesh areas
977 assert( static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() );
978 dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
979 assert( static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() );
980 dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
981 }
982
983 for( size_t i = 0; i < m_meshInputCov->faces.size(); i++ )
984 {
985 if( fabs( dSourceArea[i] - m_meshInputCov->vecFaceArea[i] ) < 1.0e-10 )
986 {
987 m_meshInputCov->vecFaceArea[i] = dSourceArea[i];
988 }
989 }
990 for( size_t i = 0; i < m_meshOutput->faces.size(); i++ )
991 {
992 if( fabs( dTargetArea[i] - m_meshOutput->vecFaceArea[i] ) < 1.0e-10 )
993 {
994 m_meshOutput->vecFaceArea[i] = dTargetArea[i];
995 }
996 }
997 }
998
999 // Set source mesh areas in map
1000 if( !m_bPointCloudSource && eInputType == DiscretizationType_FV )
1001 {
1002 this->SetSourceAreas( m_meshInputCov->vecFaceArea );
1003 if( m_meshInputCov->vecMask.size() )
1004 {
1005 this->SetSourceMask( m_meshInputCov->vecMask );
1006 }
1007 }
1008
1009 // Set target mesh areas in map
1010 if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV )
1011 {
1012 this->SetTargetAreas( m_meshOutput->vecFaceArea );
1013 if( m_meshOutput->vecMask.size() )
1014 {
1015 this->SetTargetMask( m_meshOutput->vecMask );
1016 }
1017 }
1018
1019 /*
1020 // Recalculate input mesh area from overlap mesh
1021 if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) {
1022 dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n");
1023 dbgprint.printf(0, "Recalculating source mesh areas\n");
1024 dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap);
1025 dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput);
1026 }
1027 */
1028 }
1029
1030 // Finite volume input / Finite volume output
1031 if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
1032 {
1033 // Generate reverse node array and edge map
1034 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1035 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
1036
1037 // Initialize coordinates for map
1038 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1039 this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
1040
1041 this->m_pdataGLLNodesIn = nullptr;
1042 this->m_pdataGLLNodesOut = nullptr;
1043
1044 // Finite volume input / Finite element output
1045 rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
1046 mapOptions.nPout, false, nullptr );MB_CHK_ERR( rval );
1047
1048 // Construct remap for FV-FV
1049 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1050
1051 // Construct OfflineMap
1052 if( strMapAlgorithm == "invdist" )
1053 {
1054 if( is_root ) dbgprint.printf( 0, "Calculating map (invdist)\n" );
1055 if( m_meshInputCov->faces.size() )
1056 LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1057 }
1058 else if( strMapAlgorithm == "delaunay" )
1059 {
1060 if( is_root ) dbgprint.printf( 0, "Calculating map (delaunay)\n" );
1061 if( m_meshInputCov->faces.size() )
1062 LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1063 }
1064 else if( strMapAlgorithm == "fvintbilin" )
1065 {
1066 if( is_root ) dbgprint.printf( 0, "Calculating map (intbilin)\n" );
1067 if( m_meshInputCov->faces.size() )
1068 LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1069 }
1070 else if( strMapAlgorithm == "fvintbilingb" )
1071 {
1072 if( is_root ) dbgprint.printf( 0, "Calculating map (intbilingb)\n" );
1073 if( m_meshInputCov->faces.size() )
1074 LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
1075 *this );
1076 }
1077 else if( strMapAlgorithm == "fvbilin" )
1078 {
1079 #ifdef VERBOSE
1080 if( is_root )
1081 {
1082 m_meshInputCov->Write( "SourceMeshMBTR.g" );
1083 m_meshOutput->Write( "TargetMeshMBTR.g" );
1084 }
1085 else
1086 {
1087 m_meshInputCov->Write( "SourceMeshMBTR" + std::to_string( rank ) + ".g" );
1088 m_meshOutput->Write( "TargetMeshMBTR" + std::to_string( rank ) + ".g" );
1089 }
1090 #endif
1091 if( is_root ) dbgprint.printf( 0, "Calculating map (bilin)\n" );
1092 if( m_meshInputCov->faces.size() )
1093 LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1094 }
1095 else
1096 {
1097 if( is_root ) dbgprint.printf( 0, "Calculating conservative FV-FV map\n" );
1098 if( m_meshInputCov->faces.size() )
1099 {
1100 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1101 LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
1102 ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this );
1103 #else
1104 LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
1105 #endif
1106 }
1107 }
1108 }
1109 else if( eInputType == DiscretizationType_FV )
1110 {
1111 DataArray3D< double > dataGLLJacobian;
1112
1113 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1114 double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
1115 dataGLLNodesDest, dataGLLJacobian );
1116
1117 double dNumericalArea = dNumericalArea_loc;
1118 #ifdef MOAB_HAVE_MPI
1119 if( m_pcomm )
1120 MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
1121 #endif
1122 if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
1123
1124 // Initialize coordinates for map
1125 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1126 this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1127
1128 this->m_pdataGLLNodesIn = nullptr;
1129 this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1130
1131 // Generate the continuous Jacobian
1132 bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
1133
1134 if( eOutputType == DiscretizationType_CGLL )
1135 {
1136 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
1137 }
1138 else
1139 {
1140 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
1141 }
1142
1143 // Generate reverse node array and edge map
1144 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1145 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
1146
1147 // Finite volume input / Finite element output
1148 rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
1149 mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
1150 &dataGLLNodesDest );MB_CHK_ERR( rval );
1151
1152 // Generate remap weights
1153 if( strMapAlgorithm == "volumetric" )
1154 {
1155 if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL (volumetric)\n" );
1156 LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
1157 dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *this,
1158 nMonotoneType, fContinuous, mapOptions.fNoConservation );
1159 }
1160 else
1161 {
1162 if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL\n" );
1163 LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
1164 this->GetTargetAreas(), mapOptions.nPin, *this, nMonotoneType, fContinuous,
1165 mapOptions.fNoConservation );
1166 }
1167 }
1168 else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) )
1169 {
1170 DataArray3D< double > dataGLLJacobian;
1171 if( !m_bPointCloudSource )
1172 {
1173 // Generate reverse node array and edge map
1174 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1175 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
1176
1177 // Initialize coordinates for map
1178 if( eInputType == DiscretizationType_FV )
1179 {
1180 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1181 }
1182 else
1183 {
1184 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1185 DataArray3D< double > dataGLLJacobianSrc;
1186 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1187 dataGLLJacobian );
1188 GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1189 dataGLLJacobianSrc );
1190 }
1191 }
1192 // else { /* Source is a point cloud dataset */ }
1193
1194 if( !m_bPointCloudTarget )
1195 {
1196 // Generate reverse node array and edge map
1197 if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray();
1198 if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap( false );
1199
1200 // Initialize coordinates for map
1201 if( eOutputType == DiscretizationType_FV )
1202 {
1203 this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput );
1204 }
1205 else
1206 {
1207 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1208 GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
1209 dataGLLJacobian );
1210 }
1211 }
1212 // else { /* Target is a point cloud dataset */ }
1213
1214 // Finite volume input / Finite element output
1215 rval = this->SetDOFmapAssociation(
1216 eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1217 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrcCov ),
1218 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrc ),
1219 eOutputType, mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
1220 ( m_bPointCloudTarget ? nullptr : &dataGLLNodesDest ) );MB_CHK_ERR( rval );
1221
1222 // Construct remap
1223 if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" );
1224 rval = LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ );MB_CHK_ERR( rval );
1225 }
1226 else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
1227 {
1228 DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;
1229
1230 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1231 // generate metadata for the input meshes (both source and covering source)
1232 GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1233 dataGLLJacobianSrc );
1234 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1235 dataGLLJacobian );
1236
1237 if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
1238 {
1239 _EXCEPTIONT( "Number of element does not match between metadata and "
1240 "input mesh" );
1241 }
1242
1243 // Initialize coordinates for map
1244 this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1245 this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
1246
1247 // Generate the continuous Jacobian for input mesh
1248 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1249
1250 if( eInputType == DiscretizationType_CGLL )
1251 {
1252 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
1253 }
1254 else
1255 {
1256 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
1257 }
1258
1259 // Finite element input / Finite volume output
1260 rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1261 &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1262 false, nullptr );MB_CHK_ERR( rval );
1263
1264 // Generate remap
1265 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1266
1267 if( strMapAlgorithm == "volumetric" )
1268 {
1269 _EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for"
1270 "GLL input mesh" );
1271 }
1272
1273 this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1274 this->m_pdataGLLNodesOut = nullptr;
1275
1276 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1277 LinearRemapSE4( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov, dataGLLJacobian,
1278 nMonotoneType, fContinuousIn, mapOptions.fNoConservation, mapOptions.fSparseConstraints,
1279 *this );
1280 #else
1281 LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
1282 mapOptions.fNoConservation );
1283 #endif
1284 }
1285 else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
1286 {
1287 DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
1288 DataArray3D< double > dataGLLJacobianOut;
1289
1290 // Input metadata
1291 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1292 // generate metadata for the input meshes (both source and covering source)
1293 GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1294 dataGLLJacobianSrc );
1295 // now coverage
1296 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1297 dataGLLJacobianIn );
1298 // Output metadata
1299 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1300 GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
1301 dataGLLJacobianOut );
1302
1303 // Initialize coordinates for map
1304 this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1305 this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1306
1307 // Generate the continuous Jacobian for input mesh
1308 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1309
1310 if( eInputType == DiscretizationType_CGLL )
1311 {
1312 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
1313 }
1314 else
1315 {
1316 GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
1317 }
1318
1319 // Generate the continuous Jacobian for output mesh
1320 bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
1321
1322 if( eOutputType == DiscretizationType_CGLL )
1323 {
1324 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
1325 }
1326 else
1327 {
1328 GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
1329 }
1330
1331 // Input Finite Element to Output Finite Element
1332 rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1333 &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1334 ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
1335
1336 this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1337 this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1338
1339 // Generate remap
1340 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1341
1342 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1343 LinearRemapGLLtoGLL_Integrated( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov,
1344 dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1345 this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1346 fContinuousIn, fContinuousOut, mapOptions.fSparseConstraints, *this );
1347 #else
1348 LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1349 this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1350 fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
1351 #endif
1352 }
1353 else
1354 {
1355 _EXCEPTIONT( "Not implemented" );
1356 }
1357
1358 #ifdef MOAB_HAVE_EIGEN3
1359 copy_tempest_sparsemat_to_eigen3();
1360 #endif
1361
1362 #ifdef MOAB_HAVE_MPI
1363 {
1364 // Remove ghosted entities from overlap set
1365 moab::Range ghostedEnts;
1366 rval = m_remapper->GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
1367 moab::EntityHandle m_meshOverlapSet = m_remapper->GetMeshSet( moab::Remapper::OverlapMesh );
1368 rval = m_interface->remove_entities( m_meshOverlapSet, ghostedEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
1369 }
1370 #endif
1371 // Verify consistency, conservation and monotonicity, globally
1372 if( !mapOptions.fNoCheck )
1373 {
1374 if( is_root ) dbgprint.printf( 0, "Verifying map" );
1375 this->IsConsistent( 1.0e-8 );
1376 if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 );
1377
1378 if( nMonotoneType != 0 )
1379 {
1380 this->IsMonotone( 1.0e-12 );
1381 }
1382 }
1383 }
1384 catch( Exception& e )
1385 {
1386 dbgprint.printf( 0, "%s", e.ToString().c_str() );
1387 return ( moab::MB_FAILURE );
1388 }
1389 catch( ... )
1390 {
1391 return ( moab::MB_FAILURE );
1392 }
1393 return moab::MB_SUCCESS;
1394 }
References dbgprint, moab::error(), ErrorCode, MB_ALREADY_ALLOCATED, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_EXCL, MB_TYPE_DOUBLE, and moab::Remapper::OverlapMesh.
Referenced by main().
|
inline |
Get the global Degrees-Of-Freedom ID on the source mesh.
Definition at line 576 of file TempestOnlineMap.hpp.
577 {
578 return col_gdofmap[localColID];
579 }
| int moab::TempestOnlineMap::GetDestinationGlobalNDofs | ( | ) |
Get the number of total Degrees-Of-Freedom defined on the destination mesh.
| int moab::TempestOnlineMap::GetDestinationLocalNDofs | ( | ) |
Get the number of local Degrees-Of-Freedom defined on the destination mesh.
|
inline |
Get the number of Degrees-Of-Freedom per element on the destination mesh.
Definition at line 594 of file TempestOnlineMap.hpp.
595 {
596 return m_nDofsPEl_Dest;
597 }
| const DataArray1D< double >& moab::TempestOnlineMap::GetGlobalSourceAreas | ( | ) | const |
If we computed the reduction, get the vector representing the source areas for all entities in the mesh
| const DataArray1D< double >& moab::TempestOnlineMap::GetGlobalTargetAreas | ( | ) | const |
If we computed the reduction, get the vector representing the target areas for all entities in the mesh
|
inline |
Get the index of globaColDoF.
Definition at line 581 of file TempestOnlineMap.hpp.
582 {
583 return globalColDoF + 1; // temporary
584 }
|
inline |
Get the index of globaRowDoF.
Definition at line 570 of file TempestOnlineMap.hpp.
571 {
572 return globalRowDoF + 1;
573 }
|
inline |
Get the global Degrees-Of-Freedom ID on the destination mesh.
Definition at line 565 of file TempestOnlineMap.hpp.
566 {
567 return row_gdofmap[localRowID];
568 }
| int moab::TempestOnlineMap::GetSourceGlobalNDofs | ( | ) |
Get the number of total Degrees-Of-Freedom defined on the source mesh.
| int moab::TempestOnlineMap::GetSourceLocalNDofs | ( | ) |
Get the number of local Degrees-Of-Freedom defined on the source mesh.
|
inline |
Get the number of Degrees-Of-Freedom per element on the source mesh.
Definition at line 587 of file TempestOnlineMap.hpp.
588 {
589 return m_nDofsPEl_Src;
590 }
|
virtual |
Determine if the map is conservative.
Definition at line 1444 of file TempestOnlineMap.cpp.
1445 {
1446 #ifndef MOAB_HAVE_MPI
1447
1448 return OfflineMap::IsConservative( dTolerance );
1449
1450 #else
1451 // return OfflineMap::IsConservative(dTolerance);
1452
1453 int ierr;
1454 // Get map entries
1455 DataArray1D< int > dataRows;
1456 DataArray1D< int > dataCols;
1457 DataArray1D< double > dataEntries;
1458 const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
1459 const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
1460
1461 // Calculate column sums
1462 std::vector< int > dColumnsUnique;
1463 std::vector< double > dColumnSums;
1464
1465 int nColumns = m_mapRemap.GetColumns();
1466 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1467 dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
1468 dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
1469
1470 for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
1471 {
1472 dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
1473
1474 assert( dataCols[i] < m_nTotDofs_SrcCov );
1475
1476 // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
1477 int colGID = this->GetColGlobalDoF( dataCols[i] ); // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1478 // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1479 dColumnsUnique[dataCols[i]] = colGID;
1480
1481 // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
1482 // std::endl;
1483 }
1484
1485 int rootProc = 0;
1486 std::vector< int > nElementsInProc;
1487 const int nDATA = 3;
1488 nElementsInProc.resize( size * nDATA );
1489 int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
1490 ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
1491 if( ierr != MPI_SUCCESS ) return -1;
1492
1493 int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
1494 std::vector< int > dColumnIndices;
1495 std::vector< double > dColumnSourceAreas;
1496 std::vector< double > dColumnSumsTotal;
1497 std::vector< int > displs, rcount;
1498 if( rank == rootProc )
1499 {
1500 displs.resize( size + 1, 0 );
1501 rcount.resize( size, 0 );
1502 int gsum = 0;
1503 for( int ir = 0; ir < size; ++ir )
1504 {
1505 nTotVals += nElementsInProc[ir * nDATA];
1506 nTotColumns += nElementsInProc[ir * nDATA + 1];
1507 nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
1508
1509 displs[ir] = gsum;
1510 rcount[ir] = nElementsInProc[ir * nDATA + 1];
1511 gsum += rcount[ir];
1512
1513 // printf( "%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns, displs[ir], rcount[ir], gsum );
1514 }
1515
1516 printf( "Total nnz: %d, global source elements = %d\n", nTotVals, gsum );
1517
1518 dColumnIndices.resize( nTotColumns, -1 );
1519 dColumnSumsTotal.resize( nTotColumns, 0.0 );
1520 // dColumnSourceAreas.resize ( nTotColumns, 0.0 );
1521 }
1522
1523 // Gather all ColumnSums to root process and accumulate
1524 // We expect that the sums of all columns equate to 1.0 within user specified tolerance
1525 // Need to do a gatherv here since different processes have different number of elements
1526 // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
1527 // MPI_SUM, 0, m_pcomm->comm());
1528 ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
1529 displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
1530 if( ierr != MPI_SUCCESS ) return -1;
1531 ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
1532 displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
1533 if( ierr != MPI_SUCCESS ) return -1;
1534 // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
1535 // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
1536 // MPI_SUCCESS ) return -1;
1537
1538 // Clean out unwanted arrays now
1539 dColumnSums.clear();
1540 dColumnsUnique.clear();
1541
1542 // Verify all column sums equal the input Jacobian
1543 int fConservative = 0;
1544 if( rank == rootProc )
1545 {
1546 displs[size] = ( nTotColumns );
1547 // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
1548 std::map< int, double > dColumnSumsOnRoot;
1549 // std::map<int, double> dColumnSourceAreasOnRoot;
1550 for( int ir = 0; ir < size; ir++ )
1551 {
1552 for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
1553 {
1554 if( dColumnIndices[ips] < 0 ) continue;
1555 // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
1556 assert( dColumnIndices[ips] < nTotColumnsUnq );
1557 dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips]; // / dColumnSourceAreas[ips];
1558 // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
1559 // dColumnSourceAreas[ dColumnIndices[ips] ]
1560 }
1561 }
1562
1563 for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
1564 {
1565 // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
1566 if( fabs( it->second - 1.0 ) > dTolerance )
1567 {
1568 fConservative++;
1569 Announce( "TempestOnlineMap is not conservative in column "
1570 // "%i (%1.15e)", it->first, it->second );
1571 "%i (%1.15e)",
1572 it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
1573 }
1574 }
1575 }
1576
1577 // TODO: Just do a broadcast from root instead of a reduction
1578 ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
1579 if( ierr != MPI_SUCCESS ) return -1;
1580
1581 return fConservative;
1582 #endif
1583 }
References size.
|
virtual |
Determine if the map is first-order accurate.
Definition at line 1398 of file TempestOnlineMap.cpp.
1399 {
1400 #ifndef MOAB_HAVE_MPI
1401
1402 return OfflineMap::IsConsistent( dTolerance );
1403
1404 #else
1405
1406 // Get map entries
1407 DataArray1D< int > dataRows;
1408 DataArray1D< int > dataCols;
1409 DataArray1D< double > dataEntries;
1410
1411 // Calculate row sums
1412 DataArray1D< double > dRowSums;
1413 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1414 dRowSums.Allocate( m_mapRemap.GetRows() );
1415
1416 for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1417 {
1418 dRowSums[dataRows[i]] += dataEntries[i];
1419 }
1420
1421 // Verify all row sums are equal to 1
1422 int fConsistent = 0;
1423 for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
1424 {
1425 if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
1426 {
1427 fConsistent++;
1428 int rowGID = row_gdofmap[i];
1429 Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
1430 }
1431 }
1432
1433 int ierr;
1434 int fConsistentGlobal = 0;
1435 ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1436 if( ierr != MPI_SUCCESS ) return -1;
1437
1438 return fConsistentGlobal;
1439 #endif
1440 }
|
virtual |
Determine if the map is monotone.
Definition at line 1587 of file TempestOnlineMap.cpp.
1588 {
1589 #ifndef MOAB_HAVE_MPI
1590
1591 return OfflineMap::IsMonotone( dTolerance );
1592
1593 #else
1594
1595 // Get map entries
1596 DataArray1D< int > dataRows;
1597 DataArray1D< int > dataCols;
1598 DataArray1D< double > dataEntries;
1599
1600 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1601
1602 // Verify all entries are in the range [0,1]
1603 int fMonotone = 0;
1604 for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1605 {
1606 if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
1607 {
1608 fMonotone++;
1609
1610 Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
1611 }
1612 }
1613
1614 int ierr;
1615 int fMonotoneGlobal = 0;
1616 ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1617 if( ierr != MPI_SUCCESS ) return -1;
1618
1619 return fMonotoneGlobal;
1620 #endif
1621 }
|
private |
Compute the remapping weights for a FV field defined on the source to a FV field defined on the target mesh.
Definition at line 121 of file TempestLinearRemap.cpp.
122 {
123 // Order of triangular quadrature rule
124 const int TriQuadRuleOrder = 4;
125
126 // Verify ReverseNodeArray has been calculated
127 if( m_meshInputCov->faces.size() > 0 && m_meshInputCov->revnodearray.size() == 0 )
128 {
129 _EXCEPTIONT( "ReverseNodeArray has not been calculated for m_meshInputCov" );
130 }
131
132 // Triangular quadrature rule
133 TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
134
135 // Number of coefficients needed at this order
136 #ifdef RECTANGULAR_TRUNCATION
137 int nCoefficients = nOrder * nOrder;
138 #endif
139 #ifdef TRIANGULAR_TRUNCATION
140 int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
141 #endif
142
143 // Number of faces you need
144 const int nRequiredFaceSetSize = nCoefficients;
145
146 // Fit weight exponent
147 const int nFitWeightsExponent = nOrder + 2;
148
149 // Announcements
150 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
151 dbgprint.set_prefix( "[LinearRemapFVtoFV_Tempest_MOAB]: " );
152 if( is_root )
153 {
154 dbgprint.printf( 0, "Finite Volume to Finite Volume Projection\n" );
155 dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
156 dbgprint.printf( 0, "Number of coefficients: %i\n", nCoefficients );
157 dbgprint.printf( 0, "Required adjacency set size: %i\n", nRequiredFaceSetSize );
158 dbgprint.printf( 0, "Fit weights exponent: %i\n", nFitWeightsExponent );
159 }
160
161 // Current overlap face
162 int ixOverlap = 0;
163 #ifdef VERBOSE
164 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
165 #endif
166 DataArray2D< double > dIntArray;
167 DataArray1D< double > dConstraint( nCoefficients );
168
169 // Loop through all faces on m_meshInputCov
170 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
171 {
172 // Output every 1000 elements
173 #ifdef VERBOSE
174 if( ixFirst % outputFrequency == 0 && is_root )
175 {
176 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
177 }
178 #endif
179 // Find the set of Faces that overlap faceFirst
180 int ixOverlapBegin = ixOverlap;
181 unsigned ixOverlapEnd = ixOverlapBegin;
182
183 for( ; ixOverlapEnd < m_meshOverlap->faces.size(); ixOverlapEnd++ )
184 {
185 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapEnd] != 0 ) break;
186 }
187
188 unsigned nOverlapFaces = ixOverlapEnd - ixOverlapBegin;
189
190 if( nOverlapFaces == 0 ) continue;
191
192 // Build integration array
193 BuildIntegrationArray( *m_meshInputCov, *m_meshOverlap, triquadrule, ixFirst, ixOverlapBegin, ixOverlapEnd,
194 nOrder, dIntArray );
195
196 // Set of Faces to use in building the reconstruction and associated
197 // distance metric.
198 AdjacentFaceVector vecAdjFaces;
199
200 GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixFirst, nRequiredFaceSetSize, vecAdjFaces );
201
202 // Number of adjacent Faces
203 int nAdjFaces = vecAdjFaces.size();
204
205 // Determine the conservative constraint equation
206 double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
207 dConstraint.Zero();
208 for( int p = 0; p < nCoefficients; p++ )
209 {
210 for( unsigned j = 0; j < nOverlapFaces; j++ )
211 {
212 dConstraint[p] += dIntArray[p][j];
213 }
214 dConstraint[p] /= dFirstArea;
215 }
216
217 // Build the fit array from the integration operator
218 DataArray2D< double > dFitArray;
219 DataArray1D< double > dFitWeights;
220 DataArray2D< double > dFitArrayPlus;
221
222 BuildFitArray( *m_meshInputCov, triquadrule, ixFirst, vecAdjFaces, nOrder, nFitWeightsExponent, dConstraint,
223 dFitArray, dFitWeights );
224
225 // Compute the inverse fit array
226 bool fSuccess = InvertFitArray_Corrected( dConstraint, dFitArray, dFitWeights, dFitArrayPlus );
227
228 // Multiply integration array and fit array
229 DataArray2D< double > dComposedArray( nAdjFaces, nOverlapFaces );
230 if( fSuccess )
231 {
232 // Multiply integration array and inverse fit array
233 for( int i = 0; i < nAdjFaces; i++ )
234 {
235 for( size_t j = 0; j < nOverlapFaces; j++ )
236 {
237 for( int k = 0; k < nCoefficients; k++ )
238 {
239 dComposedArray( i, j ) += dIntArray( k, j ) * dFitArrayPlus( i, k );
240 }
241 }
242 }
243
244 // Unable to invert fit array, drop to 1st order. In this case
245 // dFitArrayPlus(0,0) = 1 and all other entries are zero.
246 }
247 else
248 {
249 dComposedArray.Zero();
250 for( size_t j = 0; j < nOverlapFaces; j++ )
251 {
252 dComposedArray( 0, j ) += dIntArray( 0, j );
253 }
254 }
255
256 // Put composed array into map
257 for( unsigned i = 0; i < vecAdjFaces.size(); i++ )
258 {
259 for( unsigned j = 0; j < nOverlapFaces; j++ )
260 {
261 int& ixFirstFaceLoc = vecAdjFaces[i].first;
262 int& ixSecondFaceLoc = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
263 // int ixFirstFaceGlob = m_remapper->GetGlobalID(moab::Remapper::SourceMesh,
264 // ixFirstFaceLoc); int ixSecondFaceGlob =
265 // m_remapper->GetGlobalID(moab::Remapper::TargetMesh, ixSecondFaceLoc);
266
267 // signal to not participate, because it is a ghost target
268 if( ixSecondFaceLoc < 0 ) continue; // do not do anything
269
270 m_mapRemap( ixSecondFaceLoc, ixFirstFaceLoc ) +=
271 dComposedArray[i][j] / m_meshOutput->vecFaceArea[ixSecondFaceLoc];
272 }
273 }
274
275 // Increment the current overlap index
276 ixOverlap += nOverlapFaces;
277 }
278
279 return;
280 }
References dbgprint.
|
private |
Generate the OfflineMap for remapping from finite volumes to finite elements.
|
private |
Generate the OfflineMap for remapping from finite elements to finite elements.
Definition at line 1588 of file TempestLinearRemap.cpp.
1599 {
1600 // Triangular quadrature rule
1601 TriangularQuadratureRule triquadrule( 8 );
1602
1603 const DataArray2D< double >& dG = triquadrule.GetG();
1604 const DataArray1D< double >& dW = triquadrule.GetW();
1605
1606 // Get SparseMatrix represntation of the OfflineMap
1607 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1608
1609 // Sample coefficients
1610 DataArray2D< double > dSampleCoeffIn( nPin, nPin );
1611 DataArray2D< double > dSampleCoeffOut( nPout, nPout );
1612
1613 // Announcemnets
1614 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1615 dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
1616 if( is_root )
1617 {
1618 dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
1619 dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
1620 dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
1621 }
1622
1623 // Build the integration array for each element on m_meshOverlap
1624 DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );
1625
1626 // Number of overlap Faces per source Face
1627 DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
1628
1629 int ixOverlap = 0;
1630 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1631 {
1632 // Determine how many overlap Faces and triangles are present
1633 int nOverlapFaces = 0;
1634 size_t ixOverlapTemp = ixOverlap;
1635 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1636 {
1637 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1638 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
1639 {
1640 break;
1641 }
1642
1643 nOverlapFaces++;
1644 }
1645
1646 nAllOverlapFaces[ixFirst] = nOverlapFaces;
1647
1648 // Increment the current overlap index
1649 ixOverlap += nAllOverlapFaces[ixFirst];
1650 }
1651
1652 // Geometric area of each output node
1653 DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );
1654
1655 // Area of each overlap element in the output basis
1656 DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );
1657
1658 // Loop through all faces on m_meshInputCov
1659 ixOverlap = 0;
1660 #ifdef VERBOSE
1661 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1662 #endif
1663 if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );
1664
1665 // generic triangle used for area computation, for triangles around the center of overlap face;
1666 // used for overlap faces with more than 4 edges;
1667 // nodes array will be set for each triangle;
1668 // these triangles are not part of the mesh structure, they are just temporary during
1669 // aforementioned decomposition.
1670 Face faceTri( 3 );
1671 NodeVector nodes( 3 );
1672 faceTri.SetNode( 0, 0 );
1673 faceTri.SetNode( 1, 1 );
1674 faceTri.SetNode( 2, 2 );
1675
1676 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1677 {
1678 #ifdef VERBOSE
1679 // Announce computation progress
1680 if( ixFirst % outputFrequency == 0 && is_root )
1681 {
1682 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1683 }
1684 #endif
1685 // Quantities from the First Mesh
1686 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1687
1688 const NodeVector& nodesFirst = m_meshInputCov->nodes;
1689
1690 // Number of overlapping Faces and triangles
1691 int nOverlapFaces = nAllOverlapFaces[ixFirst];
1692
1693 if( !nOverlapFaces ) continue;
1694
1695 // // Calculate total element Jacobian
1696 // double dTotalJacobian = 0.0;
1697 // for (int s = 0; s < nPin; s++) {
1698 // for (int t = 0; t < nPin; t++) {
1699 // dTotalJacobian += dataGLLJacobianIn[s][t][ixFirst];
1700 // }
1701 // }
1702
1703 // Loop through all Overlap Faces
1704 for( int i = 0; i < nOverlapFaces; i++ )
1705 {
1706 // Quantities from the overlap Mesh
1707 const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + i];
1708
1709 const NodeVector& nodesOverlap = m_meshOverlap->nodes;
1710
1711 // Quantities from the Second Mesh
1712 int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1713
1714 // signal to not participate, because it is a ghost target
1715 if( ixSecond < 0 ) continue; // do not do anything
1716
1717 const NodeVector& nodesSecond = m_meshOutput->nodes;
1718
1719 const Face& faceSecond = m_meshOutput->faces[ixSecond];
1720
1721 int nbEdges = faceOverlap.edges.size();
1722 int nOverlapTriangles = 1;
1723 Node center; // not used if nbEdges == 3
1724 if( nbEdges > 3 )
1725 { // decompose from center in this case
1726 nOverlapTriangles = nbEdges;
1727 for( int k = 0; k < nbEdges; k++ )
1728 {
1729 const Node& node = nodesOverlap[faceOverlap[k]];
1730 center = center + node;
1731 }
1732 center = center / nbEdges;
1733 center = center.Normalized(); // project back on sphere of radius 1
1734 }
1735
1736 Node node0, node1, node2;
1737 double dTriArea;
1738
1739 // Loop over all sub-triangles of this Overlap Face
1740 for( int j = 0; j < nOverlapTriangles; j++ )
1741 {
1742 if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1743 {
1744 node0 = nodesOverlap[faceOverlap[0]];
1745 node1 = nodesOverlap[faceOverlap[1]];
1746 node2 = nodesOverlap[faceOverlap[2]];
1747 dTriArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1748 }
1749 else // decompose polygon in triangles around the center
1750 {
1751 node0 = center;
1752 node1 = nodesOverlap[faceOverlap[j]];
1753 int j1 = ( j + 1 ) % nbEdges;
1754 node2 = nodesOverlap[faceOverlap[j1]];
1755 nodes[0] = center;
1756 nodes[1] = node1;
1757 nodes[2] = node2;
1758 dTriArea = CalculateFaceArea( faceTri, nodes );
1759 }
1760
1761 for( int k = 0; k < triquadrule.GetPoints(); k++ )
1762 {
1763 // Get the nodal location of this point
1764 double dX[3];
1765
1766 dX[0] = dG( k, 0 ) * node0.x + dG( k, 1 ) * node1.x + dG( k, 2 ) * node2.x;
1767 dX[1] = dG( k, 0 ) * node0.y + dG( k, 1 ) * node1.y + dG( k, 2 ) * node2.y;
1768 dX[2] = dG( k, 0 ) * node0.z + dG( k, 1 ) * node1.z + dG( k, 2 ) * node2.z;
1769
1770 double dMag = sqrt( dX[0] * dX[0] + dX[1] * dX[1] + dX[2] * dX[2] );
1771
1772 dX[0] /= dMag;
1773 dX[1] /= dMag;
1774 dX[2] /= dMag;
1775
1776 Node nodeQuadrature( dX[0], dX[1], dX[2] );
1777
1778 // Find the components of this quadrature point in the basis
1779 // of the first Face.
1780 double dAlphaIn;
1781 double dBetaIn;
1782
1783 ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlphaIn, dBetaIn );
1784
1785 // Find the components of this quadrature point in the basis
1786 // of the second Face.
1787 double dAlphaOut;
1788 double dBetaOut;
1789
1790 ApplyInverseMap( faceSecond, nodesSecond, nodeQuadrature, dAlphaOut, dBetaOut );
1791
1792 /*
1793 // Check inverse map value
1794 if ((dAlphaIn < 0.0) || (dAlphaIn > 1.0) ||
1795 (dBetaIn < 0.0) || (dBetaIn > 1.0)
1796 ) {
1797 _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1798 dAlphaIn, dBetaIn);
1799 }
1800
1801 // Check inverse map value
1802 if ((dAlphaOut < 0.0) || (dAlphaOut > 1.0) ||
1803 (dBetaOut < 0.0) || (dBetaOut > 1.0)
1804 ) {
1805 _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1806 dAlphaOut, dBetaOut);
1807 }
1808 */
1809 // Sample the First finite element at this point
1810 SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
1811
1812 // Sample the Second finite element at this point
1813 SampleGLLFiniteElement( nMonotoneType, nPout, dAlphaOut, dBetaOut, dSampleCoeffOut );
1814
1815 // Overlap output area
1816 for( int s = 0; s < nPout; s++ )
1817 {
1818 for( int t = 0; t < nPout; t++ )
1819 {
1820 double dNodeArea = dSampleCoeffOut[s][t] * dW[k] * dTriArea;
1821
1822 dOverlapOutputArea[ixOverlap + i][s * nPout + t] += dNodeArea;
1823
1824 dGeometricOutputArea[ixSecond][s * nPout + t] += dNodeArea;
1825 }
1826 }
1827
1828 // Compute overlap integral
1829 int ixp = 0;
1830 for( int p = 0; p < nPin; p++ )
1831 {
1832 for( int q = 0; q < nPin; q++ )
1833 {
1834 int ixs = 0;
1835 for( int s = 0; s < nPout; s++ )
1836 {
1837 for( int t = 0; t < nPout; t++ )
1838 {
1839 // Sample the Second finite element at this point
1840 dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
1841 dSampleCoeffOut[s][t] * dSampleCoeffIn[p][q] * dW[k] * dTriArea;
1842
1843 ixs++;
1844 }
1845 }
1846
1847 ixp++;
1848 }
1849 }
1850 }
1851 }
1852 }
1853
1854 // Coefficients
1855 DataArray2D< double > dCoeff( nOverlapFaces * nPout * nPout, nPin * nPin );
1856
1857 for( int i = 0; i < nOverlapFaces; i++ )
1858 {
1859 // int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1860
1861 int ixp = 0;
1862 for( int p = 0; p < nPin; p++ )
1863 {
1864 for( int q = 0; q < nPin; q++ )
1865 {
1866 int ixs = 0;
1867 for( int s = 0; s < nPout; s++ )
1868 {
1869 for( int t = 0; t < nPout; t++ )
1870 {
1871 dCoeff[i * nPout * nPout + ixs][ixp] = dGlobalIntArray[ixp][ixOverlap + i][ixs] /
1872 dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1873
1874 ixs++;
1875 }
1876 }
1877
1878 ixp++;
1879 }
1880 }
1881 }
1882
1883 // Source areas
1884 DataArray1D< double > vecSourceArea( nPin * nPin );
1885
1886 for( int p = 0; p < nPin; p++ )
1887 {
1888 for( int q = 0; q < nPin; q++ )
1889 {
1890 vecSourceArea[p * nPin + q] = dataGLLJacobianIn[p][q][ixFirst];
1891 }
1892 }
1893
1894 // Target areas
1895 DataArray1D< double > vecTargetArea( nOverlapFaces * nPout * nPout );
1896
1897 for( int i = 0; i < nOverlapFaces; i++ )
1898 {
1899 // int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1900 int ixs = 0;
1901 for( int s = 0; s < nPout; s++ )
1902 {
1903 for( int t = 0; t < nPout; t++ )
1904 {
1905 vecTargetArea[i * nPout * nPout + ixs] = dOverlapOutputArea[ixOverlap + i][nPout * s + t];
1906
1907 ixs++;
1908 }
1909 }
1910 }
1911
1912 // Force consistency and conservation
1913 if( !fNoConservation )
1914 {
1915 ForceIntArrayConsistencyConservation( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType != 0 ) );
1916 }
1917
1918 // Update global coefficients
1919 for( int i = 0; i < nOverlapFaces; i++ )
1920 {
1921 int ixp = 0;
1922 for( int p = 0; p < nPin; p++ )
1923 {
1924 for( int q = 0; q < nPin; q++ )
1925 {
1926 int ixs = 0;
1927 for( int s = 0; s < nPout; s++ )
1928 {
1929 for( int t = 0; t < nPout; t++ )
1930 {
1931 dGlobalIntArray[ixp][ixOverlap + i][ixs] =
1932 dCoeff[i * nPout * nPout + ixs][ixp] * dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1933
1934 ixs++;
1935 }
1936 }
1937
1938 ixp++;
1939 }
1940 }
1941 }
1942
1943 #ifdef VVERBOSE
1944 // Check column sums (conservation)
1945 for( int i = 0; i < nPin * nPin; i++ )
1946 {
1947 double dColSum = 0.0;
1948 for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1949 {
1950 dColSum += dCoeff[j][i] * vecTargetArea[j];
1951 }
1952 printf( "Col %i: %1.15e\n", i, dColSum / vecSourceArea[i] );
1953 }
1954
1955 // Check row sums (consistency)
1956 for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1957 {
1958 double dRowSum = 0.0;
1959 for( int i = 0; i < nPin * nPin; i++ )
1960 {
1961 dRowSum += dCoeff[j][i];
1962 }
1963 printf( "Row %i: %1.15e\n", j, dRowSum );
1964 }
1965 #endif
1966
1967 // Increment the current overlap index
1968 ixOverlap += nOverlapFaces;
1969 }
1970
1971 // Build redistribution map within target element
1972 if( is_root ) dbgprint.printf( 0, "Building redistribution maps on target mesh\n" );
1973 DataArray1D< double > dRedistSourceArea( nPout * nPout );
1974 DataArray1D< double > dRedistTargetArea( nPout * nPout );
1975 std::vector< DataArray2D< double > > dRedistributionMaps;
1976 dRedistributionMaps.resize( m_meshOutput->faces.size() );
1977
1978 for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
1979 {
1980 dRedistributionMaps[ixSecond].Allocate( nPout * nPout, nPout * nPout );
1981
1982 for( int i = 0; i < nPout * nPout; i++ )
1983 {
1984 dRedistributionMaps[ixSecond][i][i] = 1.0;
1985 }
1986
1987 for( int s = 0; s < nPout * nPout; s++ )
1988 {
1989 dRedistSourceArea[s] = dGeometricOutputArea[ixSecond][s];
1990 }
1991
1992 for( int s = 0; s < nPout * nPout; s++ )
1993 {
1994 dRedistTargetArea[s] = dataGLLJacobianOut[s / nPout][s % nPout][ixSecond];
1995 }
1996
1997 if( !fNoConservation )
1998 {
1999 ForceIntArrayConsistencyConservation( dRedistSourceArea, dRedistTargetArea, dRedistributionMaps[ixSecond],
2000 ( nMonotoneType != 0 ) );
2001
2002 for( int s = 0; s < nPout * nPout; s++ )
2003 {
2004 for( int t = 0; t < nPout * nPout; t++ )
2005 {
2006 dRedistributionMaps[ixSecond][s][t] *= dRedistTargetArea[s] / dRedistSourceArea[t];
2007 }
2008 }
2009 }
2010 }
2011
2012 // Construct the total geometric area
2013 DataArray1D< double > dTotalGeometricArea( dataNodalAreaOut.GetRows() );
2014 for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
2015 {
2016 for( int s = 0; s < nPout; s++ )
2017 {
2018 for( int t = 0; t < nPout; t++ )
2019 {
2020 dTotalGeometricArea[dataGLLNodesOut[s][t][ixSecond] - 1] +=
2021 dGeometricOutputArea[ixSecond][s * nPout + t];
2022 }
2023 }
2024 }
2025
2026 // Compose the integration operator with the output map
2027 ixOverlap = 0;
2028
2029 if( is_root ) dbgprint.printf( 0, "Assembling map\n" );
2030
2031 // Map from source DOFs to target DOFs with redistribution applied
2032 DataArray2D< double > dRedistributedOp( nPin * nPin, nPout * nPout );
2033
2034 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
2035 {
2036 #ifdef VERBOSE
2037 // Announce computation progress
2038 if( ixFirst % outputFrequency == 0 && is_root )
2039 {
2040 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
2041 }
2042 #endif
2043 // Number of overlapping Faces and triangles
2044 int nOverlapFaces = nAllOverlapFaces[ixFirst];
2045
2046 if( !nOverlapFaces ) continue;
2047
2048 // Put composed array into map
2049 for( int j = 0; j < nOverlapFaces; j++ )
2050 {
2051 int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
2052
2053 // signal to not participate, because it is a ghost target
2054 if( ixSecondFace < 0 ) continue; // do not do anything
2055
2056 dRedistributedOp.Zero();
2057 for( int p = 0; p < nPin * nPin; p++ )
2058 {
2059 for( int s = 0; s < nPout * nPout; s++ )
2060 {
2061 for( int t = 0; t < nPout * nPout; t++ )
2062 {
2063 dRedistributedOp[p][s] +=
2064 dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
2065 }
2066 }
2067 }
2068
2069 int ixp = 0;
2070 for( int p = 0; p < nPin; p++ )
2071 {
2072 for( int q = 0; q < nPin; q++ )
2073 {
2074 int ixFirstNode;
2075 if( fContinuousIn )
2076 {
2077 ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
2078 }
2079 else
2080 {
2081 ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
2082 }
2083
2084 int ixs = 0;
2085 for( int s = 0; s < nPout; s++ )
2086 {
2087 for( int t = 0; t < nPout; t++ )
2088 {
2089 int ixSecondNode;
2090 if( fContinuousOut )
2091 {
2092 ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;
2093
2094 if( !fNoConservation )
2095 {
2096 smatMap( ixSecondNode, ixFirstNode ) +=
2097 dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
2098 }
2099 else
2100 {
2101 smatMap( ixSecondNode, ixFirstNode ) +=
2102 dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
2103 }
2104 }
2105 else
2106 {
2107 ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;
2108
2109 if( !fNoConservation )
2110 {
2111 smatMap( ixSecondNode, ixFirstNode ) +=
2112 dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
2113 }
2114 else
2115 {
2116 smatMap( ixSecondNode, ixFirstNode ) +=
2117 dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
2118 }
2119 }
2120
2121 ixs++;
2122 }
2123 }
2124
2125 ixp++;
2126 }
2127 }
2128 }
2129
2130 // Increment the current overlap index
2131 ixOverlap += nOverlapFaces;
2132 }
2133
2134 return;
2135 }
References center(), dbgprint, and ForceIntArrayConsistencyConservation().
|
private |
Generate the OfflineMap for remapping from finite elements to finite elements (pointwise interpolation).
Definition at line 2139 of file TempestLinearRemap.cpp.
2149 {
2150 // Gauss-Lobatto quadrature within Faces
2151 DataArray1D< double > dGL;
2152 DataArray1D< double > dWL;
2153
2154 GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );
2155
2156 // Get SparseMatrix represntation of the OfflineMap
2157 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
2158
2159 // Sample coefficients
2160 DataArray2D< double > dSampleCoeffIn( nPin, nPin );
2161
2162 // Announcemnets
2163 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
2164 dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
2165 if( is_root )
2166 {
2167 dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
2168 dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
2169 dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
2170 }
2171
2172 // Number of overlap Faces per source Face
2173 DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
2174
2175 int ixOverlap = 0;
2176
2177 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
2178 {
2179 size_t ixOverlapTemp = ixOverlap;
2180 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
2181 {
2182 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
2183
2184 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
2185
2186 nAllOverlapFaces[ixFirst]++;
2187 }
2188
2189 // Increment the current overlap index
2190 ixOverlap += nAllOverlapFaces[ixFirst];
2191 }
2192
2193 // Number of times this point was found
2194 DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );
2195
2196 ixOverlap = 0;
2197 #ifdef VERBOSE
2198 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
2199 #endif
2200 // Loop through all faces on m_meshInputCov
2201 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
2202 {
2203 #ifdef VERBOSE
2204 // Announce computation progress
2205 if( ixFirst % outputFrequency == 0 && is_root )
2206 {
2207 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
2208 }
2209 #endif
2210 // Quantities from the First Mesh
2211 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
2212
2213 const NodeVector& nodesFirst = m_meshInputCov->nodes;
2214
2215 // Number of overlapping Faces and triangles
2216 int nOverlapFaces = nAllOverlapFaces[ixFirst];
2217
2218 // Loop through all Overlap Faces
2219 for( int i = 0; i < nOverlapFaces; i++ )
2220 {
2221 // Quantities from the Second Mesh
2222 int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
2223
2224 // signal to not participate, because it is a ghost target
2225 if( ixSecond < 0 ) continue; // do not do anything
2226
2227 const NodeVector& nodesSecond = m_meshOutput->nodes;
2228 const Face& faceSecond = m_meshOutput->faces[ixSecond];
2229
2230 // Loop through all nodes on the second face
2231 for( int s = 0; s < nPout; s++ )
2232 {
2233 for( int t = 0; t < nPout; t++ )
2234 {
2235 size_t ixSecondNode;
2236 if( fContinuousOut )
2237 {
2238 ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
2239 }
2240 else
2241 {
2242 ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
2243 }
2244
2245 if( ixSecondNode >= fSecondNodeFound.GetRows() ) _EXCEPTIONT( "Logic error" );
2246
2247 // Check if this node has been found already
2248 if( fSecondNodeFound[ixSecondNode] ) continue;
2249
2250 // Check this node
2251 Node node;
2252 Node dDx1G;
2253 Node dDx2G;
2254
2255 ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );
2256
2257 // Find the components of this quadrature point in the basis
2258 // of the first Face.
2259 double dAlphaIn;
2260 double dBetaIn;
2261
2262 ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );
2263
2264 // Check if this node is within the first Face
2265 if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
2266 ( dBetaIn > 1.0 + 1.0e-10 ) )
2267 continue;
2268
2269 // Node is within the overlap region, mark as found
2270 fSecondNodeFound[ixSecondNode] = true;
2271
2272 // Sample the First finite element at this point
2273 SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
2274
2275 // Add to map
2276 for( int p = 0; p < nPin; p++ )
2277 {
2278 for( int q = 0; q < nPin; q++ )
2279 {
2280 int ixFirstNode;
2281 if( fContinuousIn )
2282 {
2283 ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
2284 }
2285 else
2286 {
2287 ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
2288 }
2289
2290 smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
2291 }
2292 }
2293 }
2294 }
2295 }
2296
2297 // Increment the current overlap index
2298 ixOverlap += nOverlapFaces;
2299 }
2300
2301 // Check for missing samples
2302 for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
2303 {
2304 if( !fSecondNodeFound[i] )
2305 {
2306 _EXCEPTION1( "Can't sample point %i", i );
2307 }
2308 }
2309
2310 return;
2311 }
References dbgprint.
|
private |
Compute the remapping weights as a permutation matrix that relates DoFs on the source mesh to DoFs on the target mesh.
Definition at line 58 of file TempestLinearRemap.cpp.
59 {
60 /* m_mapRemap size = (m_nTotDofs_Dest X m_nTotDofs_SrcCov) */
61
62 #ifdef VVERBOSE
63 {
64 std::ofstream output_file( "rowcolindices.txt", std::ios::out );
65 output_file << m_nTotDofs_Dest << " " << m_nTotDofs_SrcCov << " " << row_gdofmap.size() << " "
66 << row_ldofmap.size() << " " << col_gdofmap.size() << " " << col_ldofmap.size() << "\n";
67 output_file << "Rows \n";
68 for( unsigned iv = 0; iv < row_gdofmap.size(); iv++ )
69 output_file << row_gdofmap[iv] << " " << row_dofmap[iv] << "\n";
70 output_file << "Cols \n";
71 for( unsigned iv = 0; iv < col_gdofmap.size(); iv++ )
72 output_file << col_gdofmap[iv] << " " << col_dofmap[iv] << "\n";
73 output_file.flush(); // required here
74 output_file.close();
75 }
76 #endif
77
78 if( use_GID_matching )
79 {
80 std::map< unsigned, unsigned > src_gl;
81 for( unsigned it = 0; it < col_gdofmap.size(); ++it )
82 src_gl[col_gdofmap[it]] = it;
83
84 std::map< unsigned, unsigned >::iterator iter;
85 for( unsigned it = 0; it < row_gdofmap.size(); ++it )
86 {
87 unsigned row = row_gdofmap[it];
88 iter = src_gl.find( row );
89 if( strict_check && iter == src_gl.end() )
90 {
91 std::cout << "Searching for global target DOF " << row
92 << " but could not find correspondence in source mesh.\n";
93 assert( false );
94 }
95 else if( iter == src_gl.end() )
96 {
97 continue;
98 }
99 else
100 {
101 unsigned icol = src_gl[row];
102 unsigned irow = it;
103
104 // Set the permutation matrix in local space
105 m_mapRemap( irow, icol ) = 1.0;
106 }
107 }
108
109 return moab::MB_SUCCESS;
110 }
111 else
112 {
113 /* Create a Kd-tree to perform local queries to find nearest neighbors */
114
115 return moab::MB_FAILURE;
116 }
117 }
References col_gdofmap, m_nTotDofs_Dest, m_nTotDofs_SrcCov, MB_SUCCESS, and row_gdofmap.
|
private |
Generate the OfflineMap for linear conserative element-average spectral element to element average remapping.
|
private |
Generate the OfflineMap for cubic conserative element-average spectral element to element average remapping.
Definition at line 1164 of file TempestLinearRemap.cpp.
1169 {
1170 // Order of the polynomial interpolant
1171 int nP = dataGLLNodes.GetRows();
1172
1173 // Order of triangular quadrature rule
1174 const int TriQuadRuleOrder = 4;
1175
1176 // Triangular quadrature rule
1177 TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
1178
1179 int TriQuadraturePoints = triquadrule.GetPoints();
1180
1181 const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
1182
1183 const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
1184
1185 // Sample coefficients
1186 DataArray2D< double > dSampleCoeff( nP, nP );
1187
1188 // GLL Quadrature nodes on quadrilateral elements
1189 DataArray1D< double > dG;
1190 DataArray1D< double > dW;
1191 GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
1192
1193 // Announcements
1194 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1195 dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
1196 if( is_root )
1197 {
1198 dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
1199 dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
1200 dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
1201 }
1202
1203 // Get SparseMatrix represntation of the OfflineMap
1204 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1205
1206 // NodeVector from m_meshOverlap
1207 const NodeVector& nodesOverlap = m_meshOverlap->nodes;
1208 const NodeVector& nodesFirst = m_meshInputCov->nodes;
1209
1210 // Vector of source areas
1211 DataArray1D< double > vecSourceArea( nP * nP );
1212
1213 DataArray1D< double > vecTargetArea;
1214 DataArray2D< double > dCoeff;
1215
1216 #ifdef VERBOSE
1217 std::stringstream sstr;
1218 sstr << "remapdata_" << rank << ".txt";
1219 std::ofstream output_file( sstr.str() );
1220 #endif
1221
1222 // Current Overlap Face
1223 int ixOverlap = 0;
1224 #ifdef VERBOSE
1225 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1226 #endif
1227 // generic triangle used for area computation, for triangles around the center of overlap face;
1228 // used for overlap faces with more than 4 edges;
1229 // nodes array will be set for each triangle;
1230 // these triangles are not part of the mesh structure, they are just temporary during
1231 // aforementioned decomposition.
1232 Face faceTri( 3 );
1233 NodeVector nodes( 3 );
1234 faceTri.SetNode( 0, 0 );
1235 faceTri.SetNode( 1, 1 );
1236 faceTri.SetNode( 2, 2 );
1237
1238 // Loop over all input Faces
1239 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1240 {
1241 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1242
1243 if( faceFirst.edges.size() != 4 )
1244 {
1245 _EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
1246 }
1247 #ifdef VERBOSE
1248 // Announce computation progress
1249 if( ixFirst % outputFrequency == 0 && is_root )
1250 {
1251 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1252 }
1253 #endif
1254 // Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
1255 // on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
1256 // However, the relation with MOAB and Tempest will go out of the roof
1257
1258 // Determine how many overlap Faces and triangles are present
1259 int nOverlapFaces = 0;
1260 size_t ixOverlapTemp = ixOverlap;
1261 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1262 {
1263 // if( m_meshOverlap->vecTargetFaceIx[ixOverlapTemp] < 0 ) continue; // skip ghost target faces
1264 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1265 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
1266
1267 nOverlapFaces++;
1268 }
1269
1270 // No overlaps
1271 if( nOverlapFaces == 0 )
1272 {
1273 continue;
1274 }
1275
1276 // Allocate remap coefficients array for meshFirst Face
1277 DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );
1278
1279 // Find the local remap coefficients
1280 for( int j = 0; j < nOverlapFaces; j++ )
1281 {
1282 const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
1283 if( m_meshOverlap->vecFaceArea[ixOverlap + j] < 1.e-16 ) // machine precision
1284 {
1285 Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
1286 m_meshOverlap->vecFaceArea[ixOverlap + j] );
1287 int n = faceOverlap.edges.size();
1288 Announce( "Number nodes: %d", n );
1289 for( int k = 0; k < n; k++ )
1290 {
1291 Node nd = nodesOverlap[faceOverlap[k]];
1292 Announce( "Node %d %d : %1.10e %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
1293 }
1294 continue;
1295 }
1296
1297 // #ifdef VERBOSE
1298 // if ( is_root )
1299 // Announce ( "\tLocal ID: %i/%i = %i, areas = %2.8e", j + ixOverlap, nOverlapFaces,
1300 // m_remapper->lid_to_gid_covsrc[m_meshOverlap->vecSourceFaceIx[ixOverlap + j]],
1301 // m_meshOverlap->vecFaceArea[ixOverlap + j] );
1302 // #endif
1303
1304 int nbEdges = faceOverlap.edges.size();
1305 int nOverlapTriangles = 1;
1306 Node center; // not used if nbEdges == 3
1307 if( nbEdges > 3 )
1308 { // decompose from center in this case
1309 nOverlapTriangles = nbEdges;
1310 for( int k = 0; k < nbEdges; k++ )
1311 {
1312 const Node& node = nodesOverlap[faceOverlap[k]];
1313 center = center + node;
1314 }
1315 center = center / nbEdges;
1316 center = center.Normalized(); // project back on sphere of radius 1
1317 }
1318
1319 Node node0, node1, node2;
1320 double dTriangleArea;
1321
1322 // Loop over all sub-triangles of this Overlap Face
1323 for( int k = 0; k < nOverlapTriangles; k++ )
1324 {
1325 if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1326 {
1327 node0 = nodesOverlap[faceOverlap[0]];
1328 node1 = nodesOverlap[faceOverlap[1]];
1329 node2 = nodesOverlap[faceOverlap[2]];
1330 dTriangleArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1331 }
1332 else // decompose polygon in triangles around the center
1333 {
1334 node0 = center;
1335 node1 = nodesOverlap[faceOverlap[k]];
1336 int k1 = ( k + 1 ) % nbEdges;
1337 node2 = nodesOverlap[faceOverlap[k1]];
1338 nodes[0] = center;
1339 nodes[1] = node1;
1340 nodes[2] = node2;
1341 dTriangleArea = CalculateFaceArea( faceTri, nodes );
1342 }
1343 // Coordinates of quadrature Node
1344 for( int l = 0; l < TriQuadraturePoints; l++ )
1345 {
1346 Node nodeQuadrature;
1347 nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
1348 TriQuadratureG[l][2] * node2.x;
1349
1350 nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
1351 TriQuadratureG[l][2] * node2.y;
1352
1353 nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
1354 TriQuadratureG[l][2] * node2.z;
1355
1356 nodeQuadrature = nodeQuadrature.Normalized();
1357
1358 // Find components of quadrature point in basis
1359 // of the first Face
1360 double dAlpha;
1361 double dBeta;
1362
1363 ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );
1364
1365 // Check inverse map value
1366 if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
1367 ( dBeta > 1.0 + 1.0e-13 ) )
1368 {
1369 _EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
1370 "(%1.5e %1.5e)",
1371 j, l, dAlpha, dBeta );
1372 }
1373
1374 // Sample the finite element at this point
1375 SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );
1376
1377 // Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
1378 for( int p = 0; p < nP; p++ )
1379 {
1380 for( int q = 0; q < nP; q++ )
1381 {
1382 dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
1383 m_meshOverlap->vecFaceArea[ixOverlap + j];
1384 }
1385 }
1386 }
1387 }
1388 }
1389
1390 #ifdef VERBOSE
1391 output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
1392 for( int j = 0; j < nOverlapFaces; j++ )
1393 {
1394 for( int p = 0; p < nP; p++ )
1395 {
1396 for( int q = 0; q < nP; q++ )
1397 {
1398 output_file << dRemapCoeff[p][q][j] << " ";
1399 }
1400 }
1401 }
1402 output_file << std::endl;
1403 #endif
1404
1405 // Force consistency and conservation
1406 if( !fNoConservation )
1407 {
1408 double dTargetArea = 0.0;
1409 for( int j = 0; j < nOverlapFaces; j++ )
1410 {
1411 dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
1412 }
1413
1414 for( int p = 0; p < nP; p++ )
1415 {
1416 for( int q = 0; q < nP; q++ )
1417 {
1418 vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
1419 }
1420 }
1421
1422 const double areaTolerance = 1e-10;
1423 // Source elements are completely covered by target volumes
1424 if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
1425 {
1426 vecTargetArea.Allocate( nOverlapFaces );
1427 for( int j = 0; j < nOverlapFaces; j++ )
1428 {
1429 vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
1430 }
1431
1432 dCoeff.Allocate( nOverlapFaces, nP * nP );
1433
1434 for( int j = 0; j < nOverlapFaces; j++ )
1435 {
1436 for( int p = 0; p < nP; p++ )
1437 {
1438 for( int q = 0; q < nP; q++ )
1439 {
1440 dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
1441 }
1442 }
1443 }
1444
1445 // Target volumes only partially cover source elements
1446 }
1447 else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
1448 {
1449 double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;
1450
1451 vecTargetArea.Allocate( nOverlapFaces + 1 );
1452 for( int j = 0; j < nOverlapFaces; j++ )
1453 {
1454 vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
1455 }
1456 vecTargetArea[nOverlapFaces] = dExtraneousArea;
1457
1458 #ifdef VERBOSE
1459 Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
1460 m_meshInputCov->vecFaceArea[ixFirst] );
1461 #endif
1462 if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
1463 {
1464 _EXCEPTIONT( "Partial element area exceeds total element area" );
1465 }
1466
1467 dCoeff.Allocate( nOverlapFaces + 1, nP * nP );
1468
1469 for( int j = 0; j < nOverlapFaces; j++ )
1470 {
1471 for( int p = 0; p < nP; p++ )
1472 {
1473 for( int q = 0; q < nP; q++ )
1474 {
1475 dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
1476 }
1477 }
1478 }
1479 for( int p = 0; p < nP; p++ )
1480 {
1481 for( int q = 0; q < nP; q++ )
1482 {
1483 dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
1484 }
1485 }
1486 for( int j = 0; j < nOverlapFaces; j++ )
1487 {
1488 for( int p = 0; p < nP; p++ )
1489 {
1490 for( int q = 0; q < nP; q++ )
1491 {
1492 dCoeff[nOverlapFaces][p * nP + q] -=
1493 dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
1494 }
1495 }
1496 }
1497 for( int p = 0; p < nP; p++ )
1498 {
1499 for( int q = 0; q < nP; q++ )
1500 {
1501 dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
1502 }
1503 }
1504
1505 // Source elements only partially cover target volumes
1506 }
1507 else
1508 {
1509 Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
1510 m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
1511 _EXCEPTIONT( "Target grid must be a subset of source grid" );
1512 }
1513
1514 ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 )
1515 /*, m_remapper->lid_to_gid_covsrc[ixFirst]*/ );
1516
1517 for( int j = 0; j < nOverlapFaces; j++ )
1518 {
1519 for( int p = 0; p < nP; p++ )
1520 {
1521 for( int q = 0; q < nP; q++ )
1522 {
1523 dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
1524 }
1525 }
1526 }
1527 }
1528
1529 #ifdef VERBOSE
1530 // output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
1531 // for ( int j = 0; j < nOverlapFaces; j++ )
1532 // {
1533 // for ( int p = 0; p < nP; p++ )
1534 // {
1535 // for ( int q = 0; q < nP; q++ )
1536 // {
1537 // output_file << dRemapCoeff[p][q][j] << " ";
1538 // }
1539 // }
1540 // }
1541 // output_file << std::endl;
1542 #endif
1543
1544 // Put these remap coefficients into the SparseMatrix map
1545 for( int j = 0; j < nOverlapFaces; j++ )
1546 {
1547 int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
1548
1549 // signal to not participate, because it is a ghost target
1550 if( ixSecondFace < 0 ) continue; // do not do anything
1551
1552 for( int p = 0; p < nP; p++ )
1553 {
1554 for( int q = 0; q < nP; q++ )
1555 {
1556 if( fContinuousIn )
1557 {
1558 int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;
1559
1560 smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
1561 m_meshOverlap->vecFaceArea[ixOverlap + j] /
1562 m_meshOutput->vecFaceArea[ixSecondFace];
1563 }
1564 else
1565 {
1566 int ixFirstNode = ixFirst * nP * nP + p * nP + q;
1567
1568 smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
1569 m_meshOverlap->vecFaceArea[ixOverlap + j] /
1570 m_meshOutput->vecFaceArea[ixSecondFace];
1571 }
1572 }
1573 }
1574 }
1575 // Increment the current overlap index
1576 ixOverlap += nOverlapFaces;
1577 }
1578 #ifdef VERBOSE
1579 output_file.flush(); // required here
1580 output_file.close();
1581 #endif
1582
1583 return;
1584 }
References center(), dbgprint, and ForceConsistencyConservation3().
| void moab::TempestOnlineMap::PrintMapStatistics | ( | ) |
Print information and metadata about the remapping weights.
Definition at line 284 of file TempestLinearRemap.cpp.
285 {
286 int nrows = m_weightMatrix.rows(); // Number of rows
287 int ncols = m_weightMatrix.cols(); // Number of columns
288 int NNZ = m_weightMatrix.nonZeros(); // Number of non zero values
289 #ifdef MOAB_HAVE_MPI
290 // find out min/max for NNZ, ncols, nrows
291 // should work on std c++ 11
292 int arr3[6] = { NNZ, nrows, ncols, -NNZ, -nrows, -ncols };
293 int rarr3[6];
294 MPI_Reduce( arr3, rarr3, 6, MPI_INT, MPI_MIN, 0, m_pcomm->comm() );
295
296 int total[3];
297 MPI_Reduce( arr3, total, 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
298 if( !rank )
299 std::cout << "-> Rows (min/max/sum): (" << rarr3[1] << " / " << -rarr3[4] << " / " << total[1] << "), "
300 << " Cols (min/max/sum): (" << rarr3[2] << " / " << -rarr3[5] << " / " << total[2] << "), "
301 << " NNZ (min/max/sum): (" << rarr3[0] << " / " << -rarr3[3] << " / " << total[0] << ")\n";
302 #else
303 std::cout << "-> Rows: " << nrows << ", Cols: " << ncols << ", NNZ: " << NNZ << "\n";
304 #endif
305 }
Referenced by main().
|
private |
| moab::ErrorCode moab::TempestOnlineMap::ReadParallelMap | ( | const char * | strSource, |
| const std::vector< int > & | tgt_dof_ids, | ||
| std::vector< double > & | areaA, | ||
| int & | nA, | ||
| std::vector< double > & | areaB, | ||
| int & | nB | ||
| ) |
Generate the metadata associated with the offline map.
Read the OfflineMap from a NetCDF file.
Definition at line 1205 of file TempestOnlineMapIO.cpp.
1211 {
1212 NcError error( NcError::silent_nonfatal );
1213
1214 NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
1215 NcVar *varAreaA = NULL, *varAreaB = NULL;
1216 int nS = 0;
1217 #ifdef MOAB_HAVE_PNETCDF
1218 // some variables will be used just in the case netcdfpar reader fails
1219 int ncfile = -1;
1220 int ndims, nvars, ngatts, unlimited;
1221 #endif
1222 #ifdef MOAB_HAVE_NETCDFPAR
1223 bool is_independent = true;
1224 ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
1225 // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
1226 #else
1227 NcFile ncMap( strSource, NcFile::ReadOnly );
1228 #endif
1229
1230 #define CHECK_EXCEPTION( obj, type, varstr ) \
1231 { \
1232 if( obj == nullptr ) \
1233 { \
1234 _EXCEPTION3( "Map file \"%s\" does not contain %s \"%s\"", strSource, type, varstr ); \
1235 } \
1236 }
1237
1238 // Read SparseMatrix entries
1239
1240 if( ncMap.is_valid() )
1241 {
1242 NcDim* dimNS = ncMap.get_dim( "n_s" );
1243 CHECK_EXCEPTION( dimNS, "dimension", "n_s" );
1244
1245 NcDim* dimNA = ncMap.get_dim( "n_a" );
1246 CHECK_EXCEPTION( dimNA, "dimension", "n_a" );
1247
1248 NcDim* dimNB = ncMap.get_dim( "n_b" );
1249 CHECK_EXCEPTION( dimNB, "dimension", "n_b" );
1250
1251 // store total number of nonzeros
1252 nS = dimNS->size();
1253 nA = dimNA->size();
1254 nB = dimNB->size();
1255
1256 varRow = ncMap.get_var( "row" );
1257 CHECK_EXCEPTION( varRow, "variable", "row" );
1258
1259 varCol = ncMap.get_var( "col" );
1260 CHECK_EXCEPTION( varCol, "variable", "col" );
1261
1262 varS = ncMap.get_var( "S" );
1263 CHECK_EXCEPTION( varS, "variable", "S" );
1264
1265 varAreaA = ncMap.get_var( "area_a" );
1266 CHECK_EXCEPTION( varAreaA, "variable", "area_a" );
1267
1268 varAreaB = ncMap.get_var( "area_b" );
1269 CHECK_EXCEPTION( varAreaB, "variable", "area_b" );
1270
1271 #ifdef MOAB_HAVE_NETCDFPAR
1272 ncMap.enable_var_par_access( varRow, is_independent );
1273 ncMap.enable_var_par_access( varCol, is_independent );
1274 ncMap.enable_var_par_access( varS, is_independent );
1275 ncMap.enable_var_par_access( varAreaA, is_independent );
1276 ncMap.enable_var_par_access( varAreaB, is_independent );
1277 #endif
1278 }
1279 else
1280 {
1281 #ifdef MOAB_HAVE_PNETCDF
1282 // read the file using pnetcdf directly, in parallel; need to have MPI, we do not check that anymore
1283 // why build wth pnetcdf without MPI ?
1284 // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
1285 ERR_PARNC(
1286 ncmpi_open( m_pcomm->comm(), strSource, NC_NOWRITE, MPI_INFO_NULL, &ncfile ) ); // bail out completely
1287 ERR_PARNC( ncmpi_inq( ncfile, &ndims, &nvars, &ngatts, &unlimited ) );
1288 // find dimension ids for n_S
1289 int ins;
1290 ERR_PARNC( ncmpi_inq_dimid( ncfile, "n_s", &ins ) );
1291 MPI_Offset leng;
1292 ERR_PARNC( ncmpi_inq_dimlen( ncfile, ins, &leng ) );
1293 nS = (int)leng;
1294 ERR_PARNC( ncmpi_inq_dimid( ncfile, "n_a", &ins ) );
1295 ERR_PARNC( ncmpi_inq_dimlen( ncfile, ins, &leng ) );
1296 nA = (int)leng;
1297 ERR_PARNC( ncmpi_inq_dimid( ncfile, "n_b", &ins ) );
1298 ERR_PARNC( ncmpi_inq_dimlen( ncfile, ins, &leng ) );
1299 nB = (int)leng;
1300 #else
1301 _EXCEPTION1( "cannot read the file %s", strSource );
1302 #endif
1303 }
1304
1305 // Let us declare the map object for every process
1306 SparseMatrix< double >& sparseMatrix = this->GetSparseMatrix();
1307
1308 int localSize = nS / size;
1309 long offsetRead = rank * localSize;
1310 // leftovers on last rank
1311 if( rank == size - 1 )
1312 {
1313 localSize += nS % size;
1314 }
1315
1316 int localSizeA = nA / size;
1317 long offsetReadA = rank * localSizeA;
1318 // leftovers on last rank
1319 if( rank == size - 1 )
1320 {
1321 localSizeA += nA % size;
1322 }
1323
1324 int localSizeB = nB / size;
1325 long offsetReadB = rank * localSizeB;
1326 // leftovers on last rank
1327 if( rank == size - 1 )
1328 {
1329 localSizeB += nB % size;
1330 }
1331
1332 std::vector< int > vecRow, vecCol;
1333 std::vector< double > vecS;
1334 vecRow.resize( localSize );
1335 vecCol.resize( localSize );
1336 vecS.resize( localSize );
1337 vecAreaA.resize( localSizeA );
1338 vecAreaB.resize( localSizeB );
1339
1340 if( ncMap.is_valid() )
1341 {
1342 varRow->set_cur( (long)( offsetRead ) );
1343 varRow->get( &( vecRow[0] ), localSize );
1344
1345 varCol->set_cur( (long)( offsetRead ) );
1346 varCol->get( &( vecCol[0] ), localSize );
1347
1348 varS->set_cur( (long)( offsetRead ) );
1349 varS->get( &( vecS[0] ), localSize );
1350
1351 varAreaA->set_cur( (long)( offsetReadA ) );
1352 varAreaA->get( &( vecAreaA[0] ), localSizeA );
1353
1354 varAreaB->set_cur( (long)( offsetReadB ) );
1355 varAreaB->get( &( vecAreaB[0] ), localSizeB );
1356
1357 ncMap.close();
1358 }
1359 else
1360 {
1361 #ifdef MOAB_HAVE_PNETCDF
1362 // fill the local vectors with the variables from pnetcdf file; first inquire, then fill
1363 MPI_Offset start = (MPI_Offset)offsetRead;
1364 MPI_Offset count = (MPI_Offset)localSize;
1365 int varid;
1366 // get the sparse matrix values, row and column indices
1367 ERR_PARNC( ncmpi_inq_varid( ncfile, "S", &varid ) );
1368 ERR_PARNC( ncmpi_get_vara_double_all( ncfile, varid, &start, &count, &vecS[0] ) );
1369 ERR_PARNC( ncmpi_inq_varid( ncfile, "row", &varid ) );
1370 ERR_PARNC( ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecRow[0] ) );
1371 ERR_PARNC( ncmpi_inq_varid( ncfile, "col", &varid ) );
1372 ERR_PARNC( ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecCol[0] ) );
1373
1374 ERR_PARNC( ncmpi_inq_varid( ncfile, "area_a", &varid ) );
1375 MPI_Offset startA = (MPI_Offset)offsetReadA;
1376 MPI_Offset countA = (MPI_Offset)localSizeA;
1377 ERR_PARNC( ncmpi_get_vara_double_all( ncfile, varid, &startA, &countA, &vecAreaA[0] ) );
1378
1379 ERR_PARNC( ncmpi_inq_varid( ncfile, "area_b", &varid ) );
1380 MPI_Offset startB = (MPI_Offset)offsetReadB;
1381 MPI_Offset countB = (MPI_Offset)localSizeB;
1382 ERR_PARNC( ncmpi_get_vara_double_all( ncfile, varid, &startB, &countB, &vecAreaB[0] ) );
1383
1384 ERR_PARNC( ncmpi_close( ncfile ) );
1385 #endif
1386 }
1387
1388 // NOTE: we can check if the following workflow would help here.
1389 // Initial experiments did not show same index map.
1390 // Might be useful to understand why that is the case.
1391 // First, ensure that all local and global indices are computed
1392 // this->m_remapper->ComputeGlobalLocalMaps();
1393 // Next, reuse the global to local map
1394 // const auto& rowMap = m_remapper->gid_to_lid_tgt;
1395 // const auto& colMap = m_remapper->gid_to_lid_covsrc;
1396
1397 #ifdef MOAB_HAVE_MPI
1398 // bother with tuple list only if size > 1
1399 // otherwise, just fill the sparse matrix
1400 if( size > 1 )
1401 {
1402 std::vector< int > ownership;
1403 // the default trivial partitioning scheme
1404 int nDofs = nB; // this is for row partitioning
1405
1406 // assert(row_major_ownership == true); // this block is valid only for row-based partitioning
1407 ownership.resize( size );
1408 int nPerPart = nDofs / size;
1409 int nRemainder = nDofs % size; // Keep the remainder in root
1410 ownership[0] = nPerPart + nRemainder;
1411 for( int ip = 1, roffset = ownership[0]; ip < size; ++ip )
1412 {
1413 roffset += nPerPart;
1414 ownership[ip] = roffset;
1415 }
1416 moab::TupleList* tl = new moab::TupleList;
1417 unsigned numr = 1; //
1418 tl->initialize( 3, 0, 0, numr, localSize ); // to proc, row, col, value
1419 tl->enableWriteAccess();
1420 // populate
1421 for( int i = 0; i < localSize; i++ )
1422 {
1423 int rowval = vecRow[i] - 1; // dofs are 1 based in the file
1424 int colval = vecCol[i] - 1;
1425 int to_proc = -1;
1426 int dof_val = rowval;
1427
1428 if( ownership[0] > dof_val )
1429 to_proc = 0;
1430 else
1431 {
1432 for( int ip = 1; ip < size; ++ip )
1433 {
1434 if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
1435 {
1436 to_proc = ip;
1437 break;
1438 }
1439 }
1440 }
1441
1442 int n = tl->get_n();
1443 tl->vi_wr[3 * n] = to_proc;
1444 tl->vi_wr[3 * n + 1] = rowval;
1445 tl->vi_wr[3 * n + 2] = colval;
1446 tl->vr_wr[n] = vecS[i];
1447 tl->inc_n();
1448 }
1449 // heavy communication
1450 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl, 0 );
1451
1452 if( owned_dof_ids.size() > 0 )
1453 {
1454 // we need to send desired dof to the rendezvous point
1455 moab::TupleList tl_re; //
1456 tl_re.initialize( 2, 0, 0, 0, owned_dof_ids.size() ); // to proc, value
1457 tl_re.enableWriteAccess();
1458 // send first to rendez_vous point, decided by trivial partitioning
1459
1460 for( size_t i = 0; i < owned_dof_ids.size(); i++ )
1461 {
1462 int to_proc = -1;
1463 int dof_val = owned_dof_ids[i] - 1; // dofs are 1 based in the file, partition from 0 ?
1464
1465 if( ownership[0] > dof_val )
1466 to_proc = 0;
1467 else
1468 {
1469 for( int ip = 1; ip < size; ++ip )
1470 {
1471 if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
1472 {
1473 to_proc = ip;
1474 break;
1475 }
1476 }
1477 }
1478
1479 int n = tl_re.get_n();
1480 tl_re.vi_wr[2 * n] = to_proc;
1481 tl_re.vi_wr[2 * n + 1] = dof_val;
1482
1483 tl_re.inc_n();
1484 }
1485 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl_re, 0 );
1486 // now we know in tl_re where do we need to send back dof_val
1487 moab::TupleList::buffer sort_buffer;
1488 sort_buffer.buffer_init( tl_re.get_n() );
1489 tl_re.sort( 1, &sort_buffer ); // so now we order by value
1490
1491 sort_buffer.buffer_init( tl->get_n() );
1492
1493 std::map< int, int > startDofIndex, endDofIndex; // indices in tl_re for values we want
1494 int dofVal = -1;
1495 if( tl_re.get_n() > 0 ) dofVal = tl_re.vi_rd[1]; // first dof val on this rank
1496 startDofIndex[dofVal] = 0;
1497 endDofIndex[dofVal] = 0; // start and end
1498 for( unsigned k = 1; k < tl_re.get_n(); k++ )
1499 {
1500 int newDof = tl_re.vi_rd[2 * k + 1];
1501 if( dofVal == newDof )
1502 {
1503 endDofIndex[dofVal] = k; // increment by 1 actually
1504 }
1505 else
1506 {
1507 dofVal = newDof;
1508 startDofIndex[dofVal] = k;
1509 endDofIndex[dofVal] = k;
1510 }
1511 }
1512
1513 // basically, for each value we are interested in, index in tl_re with those values are
1514 // tl_re.vi_rd[2*startDofIndex+1] == valDof == tl_re.vi_rd[2*endDofIndex+1]
1515 // so now we have ordered
1516 // tl_re shows to what proc do we need to send the tuple (row, col, val)
1517 moab::TupleList* tl_back = new moab::TupleList;
1518 unsigned numr = 1; //
1519 // localSize is a good guess, but maybe it should be bigger ?
1520 // this could be bigger for repeated dofs
1521 tl_back->initialize( 3, 0, 0, numr, tl->get_n() ); // to proc, row, col, value
1522 tl_back->enableWriteAccess();
1523 // now loop over tl and tl_re to see where to send
1524 // form the new tuple, which will contain the desired dofs per task, per row or column distribution
1525
1526 for( unsigned k = 0; k < tl->get_n(); k++ )
1527 {
1528 int valDof = tl->vi_rd[3 * k + 1]; // 1 for row, 2 for column // first value, it should be
1529 for( int ire = startDofIndex[valDof]; ire <= endDofIndex[valDof]; ire++ )
1530 {
1531 int to_proc = tl_re.vi_rd[2 * ire];
1532 int n = tl_back->get_n();
1533 tl_back->vi_wr[3 * n] = to_proc;
1534 tl_back->vi_wr[3 * n + 1] = tl->vi_rd[3 * k + 1]; // row
1535 tl_back->vi_wr[3 * n + 2] = tl->vi_rd[3 * k + 2]; // col
1536 tl_back->vr_wr[n] = tl->vr_rd[k];
1537 tl_back->inc_n();
1538 }
1539 }
1540
1541 // now communicate to the desired tasks:
1542 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl_back, 0 );
1543
1544 tl_re.reset(); // clear memory, although this will go out of scope
1545 tl->reset();
1546 tl = tl_back;
1547 }
1548
1549 int rindexMax = 0, cindexMax = 0;
1550 // populate the sparsematrix, using rowMap and colMap
1551 int n = tl->get_n();
1552 for( int i = 0; i < n; i++ )
1553 {
1554 int rindex, cindex;
1555 const int vecRowValue = tl->vi_wr[3 * i + 1];
1556 const int vecColValue = tl->vi_wr[3 * i + 2];
1557
1558 const auto riter = rowMap.find( vecRowValue );
1559 if( riter == rowMap.end() )
1560 {
1561 rowMap[vecRowValue] = rindexMax;
1562 rindex = rindexMax;
1563 row_gdofmap.push_back( vecRowValue );
1564 // row_dtoc_dofmap.push_back( vecRowValue );
1565 rindexMax++;
1566 }
1567 else
1568 rindex = riter->second;
1569
1570 const auto citer = colMap.find( vecColValue );
1571 if( citer == colMap.end() )
1572 {
1573 colMap[vecColValue] = cindexMax;
1574 cindex = cindexMax;
1575 col_gdofmap.push_back( vecColValue );
1576 // col_dtoc_dofmap.push_back( vecColValue );
1577 cindexMax++;
1578 }
1579 else
1580 cindex = citer->second;
1581
1582 sparseMatrix( rindex, cindex ) = tl->vr_wr[i];
1583 }
1584 tl->reset();
1585 }
1586 else
1587 #endif
1588 {
1589 int rindexMax = 0, cindexMax = 0;
1590 for( int i = 0; i < nS; i++ )
1591 {
1592 int rindex, cindex;
1593 const int vecRowValue = vecRow[i] - 1; // the rows, cols are 1 based in the file
1594 const int vecColValue = vecCol[i] - 1;
1595
1596 const auto riter = rowMap.find( vecRowValue );
1597 if( riter == rowMap.end() )
1598 {
1599 rowMap[vecRowValue] = rindexMax;
1600 rindex = rindexMax;
1601 row_gdofmap.push_back( vecRowValue );
1602 // row_dtoc_dofmap.push_back( vecRowValue );
1603 rindexMax++;
1604 }
1605 else
1606 rindex = riter->second;
1607
1608 const auto citer = colMap.find( vecColValue );
1609 if( citer == colMap.end() )
1610 {
1611 colMap[vecColValue] = cindexMax;
1612 cindex = cindexMax;
1613 col_gdofmap.push_back( vecColValue );
1614 // col_dtoc_dofmap.push_back( vecColValue );
1615 cindexMax++;
1616 }
1617 else
1618 cindex = citer->second;
1619
1620 sparseMatrix( rindex, cindex ) = vecS[i];
1621 }
1622 }
1623
1624 m_nTotDofs_Src = sparseMatrix.GetColumns();
1625 m_nTotDofs_SrcCov = m_nTotDofs_Src;
1626 m_nTotDofs_Dest = sparseMatrix.GetRows();
1627 // TODO: make this flexible and read the order from map with help of metadata
1628 m_nDofsPEl_Src = 1; // always assume FV-FV maps are read from file
1629 m_nDofsPEl_Dest = 1; // always assume FV-FV maps are read from file
1630
1631 #ifdef MOAB_HAVE_EIGEN3
1632 this->copy_tempest_sparsemat_to_eigen3();
1633 #endif
1634
1635 // Reset the source and target data first
1636 m_rowVector.setZero();
1637 m_colVector.setZero();
1638
1639 #ifdef VERBOSE
1640 serializeSparseMatrix( m_weightMatrix, "map_operator_" + std::to_string( rank ) + ".txt" );
1641 #endif
1642 return moab::MB_SUCCESS;
1643 }
References CHECK_EXCEPTION, moab::TupleList::enableWriteAccess(), moab::error(), moab::TupleList::get_n(), moab::TupleList::inc_n(), moab::TupleList::initialize(), MB_SUCCESS, moab::TupleList::reset(), size, moab::TupleList::sort(), moab::TupleList::vi_rd, moab::TupleList::vi_wr, moab::TupleList::vr_rd, and moab::TupleList::vr_wr.
|
private |
| moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs | ( | std::vector< int > & | values_entities | ) |
Definition at line 656 of file TempestOnlineMap.cpp.
657 {
658 // col_gdofmap has global dofs , that should be in the list of values, such that
659 // row_dtoc_dofmap[offsetDOF] = localDOF;
660 // we need to find col_dtoc_dofmap such that: col_gdofmap[ col_dtoc_dofmap[i] ] == values_entities [i];
661 // we know that col_gdofmap[0..(nbcols-1)] = global_col_dofs -> in values_entities
662 // form first inverse
663 //
664 // resize and initialize to -1 to signal that this value should not be used, if not set below
665 col_dtoc_dofmap.resize( values_entities.size(), -1 );
666 for( size_t j = 0; j < values_entities.size(); j++ )
667 {
668 // values are 1 based, but rowMap, colMap are not
669 const auto it = colMap.find( values_entities[j] - 1 );
670 if( it != colMap.end() ) col_dtoc_dofmap[j] = it->second;
671 }
672 return moab::MB_SUCCESS;
673 }
References MB_SUCCESS.
| moab::ErrorCode moab::TempestOnlineMap::set_row_dc_dofs | ( | std::vector< int > & | values_entities | ) |
Definition at line 675 of file TempestOnlineMap.cpp.
676 {
677 // we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
678 // resize and initialize to -1 to signal that this value should not be used, if not set below
679 row_dtoc_dofmap.resize( values_entities.size(), -1 );
680 for( size_t j = 0; j < values_entities.size(); j++ )
681 {
682 // values are 1 based, but rowMap, colMap are not
683 const auto it = rowMap.find( values_entities[j] - 1 );
684 if( it != rowMap.end() ) row_dtoc_dofmap[j] = it->second;
685 }
686 return moab::MB_SUCCESS;
687 }
References MB_SUCCESS.
|
inline |
Get the number of Degrees-Of-Freedom per element on the destination mesh.
Definition at line 605 of file TempestOnlineMap.hpp.
606 { 607 m_nDofsPEl_Dest = nt; 608 }
| moab::ErrorCode moab::TempestOnlineMap::SetDOFmapAssociation | ( | DiscretizationType | srcType, |
| int | srcOrder, | ||
| bool | isSrcContinuous, | ||
| DataArray3D< int > * | srcdataGLLNodes, | ||
| DataArray3D< int > * | srcdataGLLNodesSrc, | ||
| DiscretizationType | destType, | ||
| int | destOrder, | ||
| bool | isTgtContinuous, | ||
| DataArray3D< int > * | tgtdataGLLNodes | ||
| ) |
Compute the association between the solution tag global DoF numbering and the local matrix numbering so that matvec operations can be performed consistently.
| srcType | The discretization type of the source mesh |
| srcOrder | The order of the discretization on the source mesh |
| isSrcContinuous | The continuity of the discretization on the source mesh |
| srcdataGLLNodes | The GLL nodes on the source mesh |
| srcdataGLLNodesSrc | The GLL nodes on the source mesh |
| destType | The discretization type of the destination mesh |
| destOrder | The order of the discretization on the destination mesh |
| isTgtContinuous | The continuity of the discretization on the destination mesh |
| tgtdataGLLNodes | The GLL nodes on the destination mesh |
Definition at line 163 of file TempestOnlineMap.cpp.
172 {
173 std::vector< bool > dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
174 std::vector< int > src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;
175
176 // We are assuming that these are element based tags that are sized: np * np
177 m_srcDiscType = srcType;
178 m_destDiscType = destType;
179 m_input_order = srcOrder;
180 m_output_order = destOrder;
181
182 bool vprint = is_root && false;
183
184 // Compute and store the total number of source and target DoFs corresponding
185 // to number of rows and columns in the mapping.
186 #ifdef VVERBOSE
187 {
188 src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, -1 );
189 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
190 locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
191 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
192 tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
193 rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
194
195 if( is_root )
196 {
197 {
198 std::ofstream output_file( "sourcecov-gids-0.txt" );
199 output_file << "I, GDOF\n";
200 for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
201 output_file << i << ", " << src_soln_gdofs[i] << "\n";
202
203 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
204 m_nTotDofs_SrcCov = 0;
205 if( isSrcContinuous )
206 dgll_cgll_covcol_ldofmap.resize(
207 m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
208 for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
209 {
210 for( int p = 0; p < m_nDofsPEl_Src; p++ )
211 {
212 for( int q = 0; q < m_nDofsPEl_Src; q++ )
213 {
214 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
215 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
216 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
217 {
218 m_nTotDofs_SrcCov++;
219 dgll_cgll_covcol_ldofmap[localDOF] = true;
220 }
221 output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
222 << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
223 }
224 }
225 }
226 output_file.flush(); // required here
227 output_file.close();
228 dgll_cgll_covcol_ldofmap.clear();
229 }
230
231 {
232 std::ofstream output_file( "source-gids-0.txt" );
233 output_file << "I, GDOF\n";
234 for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
235 output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
236
237 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
238 m_nTotDofs_Src = 0;
239 if( isSrcContinuous )
240 dgll_cgll_col_ldofmap.resize(
241 m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
242 for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
243 {
244 for( int p = 0; p < m_nDofsPEl_Src; p++ )
245 {
246 for( int q = 0; q < m_nDofsPEl_Src; q++ )
247 {
248 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
249 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
250 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
251 {
252 m_nTotDofs_Src++;
253 dgll_cgll_col_ldofmap[localDOF] = true;
254 }
255 output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
256 << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
257 }
258 }
259 }
260 output_file.flush(); // required here
261 output_file.close();
262 dgll_cgll_col_ldofmap.clear();
263 }
264
265 {
266 std::ofstream output_file( "target-gids-0.txt" );
267 output_file << "I, GDOF\n";
268 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
269 output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
270
271 output_file << "ELEMID, IDOF, GDOF, NDOF\n";
272 m_nTotDofs_Dest = 0;
273
274 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
275 {
276 output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
277 << m_nTotDofs_Dest << "\n";
278 m_nTotDofs_Dest++;
279 }
280
281 output_file.flush(); // required here
282 output_file.close();
283 }
284 }
285 else
286 {
287 {
288 std::ofstream output_file( "sourcecov-gids-1.txt" );
289 output_file << "I, GDOF\n";
290 for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
291 output_file << i << ", " << src_soln_gdofs[i] << "\n";
292
293 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
294 m_nTotDofs_SrcCov = 0;
295 if( isSrcContinuous )
296 dgll_cgll_covcol_ldofmap.resize(
297 m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
298 for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
299 {
300 for( int p = 0; p < m_nDofsPEl_Src; p++ )
301 {
302 for( int q = 0; q < m_nDofsPEl_Src; q++ )
303 {
304 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
305 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
306 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
307 {
308 m_nTotDofs_SrcCov++;
309 dgll_cgll_covcol_ldofmap[localDOF] = true;
310 }
311 output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
312 << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
313 }
314 }
315 }
316 output_file.flush(); // required here
317 output_file.close();
318 dgll_cgll_covcol_ldofmap.clear();
319 }
320
321 {
322 std::ofstream output_file( "source-gids-1.txt" );
323 output_file << "I, GDOF\n";
324 for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
325 output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
326
327 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
328 m_nTotDofs_Src = 0;
329 if( isSrcContinuous )
330 dgll_cgll_col_ldofmap.resize(
331 m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
332 for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
333 {
334 for( int p = 0; p < m_nDofsPEl_Src; p++ )
335 {
336 for( int q = 0; q < m_nDofsPEl_Src; q++ )
337 {
338 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
339 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
340 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
341 {
342 m_nTotDofs_Src++;
343 dgll_cgll_col_ldofmap[localDOF] = true;
344 }
345 output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
346 << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
347 }
348 }
349 }
350 output_file.flush(); // required here
351 output_file.close();
352 dgll_cgll_col_ldofmap.clear();
353 }
354
355 {
356 std::ofstream output_file( "target-gids-1.txt" );
357 output_file << "I, GDOF\n";
358 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
359 output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
360
361 output_file << "ELEMID, IDOF, GDOF, NDOF\n";
362 m_nTotDofs_Dest = 0;
363
364 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
365 {
366 output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
367 << m_nTotDofs_Dest << "\n";
368 m_nTotDofs_Dest++;
369 }
370
371 output_file.flush(); // required here
372 output_file.close();
373 }
374 }
375 }
376 #endif
377
378 // Now compute the mapping and store it for the covering mesh
379 int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
380 if( m_remapper->point_cloud_source )
381 {
382 assert( m_nDofsPEl_Src == 1 );
383 col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
384 col_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
385 src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
386 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] ) );
387 srcTagSize = 1;
388 }
389 else
390 {
391 col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
392 col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
393 src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
394 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] ) );
395 }
396
397 #ifdef ALTERNATE_NUMBERING_IMPLEMENTATION
398 unsigned maxSrcIndx = 0;
399
400 // for ( unsigned j = 0; j < m_covering_source_entities.size(); j++ )
401 std::vector< int > locdofs( srcTagSize );
402 std::map< Node, moab::EntityHandle > mapLocalMBNodes;
403 double elcoords[3];
404 for( unsigned iel = 0; iel < m_remapper->m_covering_source_entities.size(); ++iel )
405 {
406 EntityHandle eh = m_remapper->m_covering_source_entities[iel];
407 rval = m_interface->get_coords( &eh, 1, elcoords );MB_CHK_ERR( rval );
408 Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
409 mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
410 }
411
412 const NodeVector& nodes = m_remapper->m_covering_source->nodes;
413 for( unsigned j = 0; j < m_remapper->m_covering_source->faces.size(); j++ )
414 {
415 const Face& face = m_remapper->m_covering_source->faces[j];
416
417 Node centroid;
418 centroid.x = centroid.y = centroid.z = 0.0;
419 for( unsigned l = 0; l < face.edges.size(); ++l )
420 {
421 centroid.x += nodes[face[l]].x;
422 centroid.y += nodes[face[l]].y;
423 centroid.z += nodes[face[l]].z;
424 }
425 const double factor = 1.0 / face.edges.size();
426 centroid.x *= factor;
427 centroid.y *= factor;
428 centroid.z *= factor;
429
430 EntityHandle current_eh;
431 if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
432 {
433 current_eh = mapLocalMBNodes[centroid];
434 }
435
436 rval = m_interface->tag_get_data( m_dofTagSrc, ¤t_eh, 1, &locdofs[0] );MB_CHK_ERR( rval );
437 for( int p = 0; p < m_nDofsPEl_Src; p++ )
438 {
439 for( int q = 0; q < m_nDofsPEl_Src; q++ )
440 {
441 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
442 const int offsetDOF = p * m_nDofsPEl_Src + q;
443 maxSrcIndx = ( localDOF > maxSrcIndx ? localDOF : maxSrcIndx );
444 std::cout << "Col: " << current_eh << ", " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF
445 << ", " << localDOF << ", " << locdofs[offsetDOF] - 1 << ", " << maxSrcIndx << "\n";
446 }
447 }
448 }
449 #endif
450
451 m_nTotDofs_SrcCov = 0;
452 if( srcdataGLLNodes == nullptr )
453 {
454 /* we only have a mapping for elements as DoFs */
455 std::vector<int> sorted_tmp_dofs( src_soln_gdofs.size() );
456 std::copy( src_soln_gdofs.begin(), src_soln_gdofs.end(), sorted_tmp_dofs.begin() );
457 std::sort( sorted_tmp_dofs.begin(), sorted_tmp_dofs.end() );
458
459 for( unsigned i = 0; i < col_gdofmap.size(); ++i )
460 {
461 auto gdof = sorted_tmp_dofs[i];
462 // printf("%d: Column -- Unsorted: %d, Sorted: %d\n", rank, src_soln_gdofs[i], sorted_tmp_dofs[i]);
463 assert( gdof > 0 );
464 col_gdofmap[i] = gdof - 1;
465 col_dtoc_dofmap[i] = i;
466 if( vprint ) std::cout << "Col: " << i << ", " << col_gdofmap[i] << "\n";
467 m_nTotDofs_SrcCov++;
468 }
469 }
470 else
471 {
472 if( isSrcContinuous )
473 dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, false );
474 // Put these remap coefficients into the SparseMatrix map
475 for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
476 {
477 for( int p = 0; p < m_nDofsPEl_Src; p++ )
478 {
479 for( int q = 0; q < m_nDofsPEl_Src; q++ )
480 {
481 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
482 const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
483 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
484 {
485 m_nTotDofs_SrcCov++;
486 dgll_cgll_covcol_ldofmap[localDOF] = true;
487 }
488 if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
489 assert( src_soln_gdofs[offsetDOF] > 0 );
490 col_gdofmap[localDOF] = src_soln_gdofs[offsetDOF] - 1;
491 col_dtoc_dofmap[offsetDOF] = localDOF;
492 if( vprint )
493 std::cout << "Col: " << offsetDOF << ", " << localDOF << ", " << col_gdofmap[offsetDOF] << ", "
494 << m_nTotDofs_SrcCov << "\n";
495 }
496 }
497 }
498 }
499
500 if( m_remapper->point_cloud_source )
501 {
502 assert( m_nDofsPEl_Src == 1 );
503 srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
504 srccol_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
505 locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
506 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] ) );
507 }
508 else
509 {
510 srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
511 srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
512 locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
513 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] ) );
514 }
515
516 // Now compute the mapping and store it for the original source mesh
517 m_nTotDofs_Src = 0;
518 if( srcdataGLLNodesSrc == nullptr )
519 {
520 /* we only have a mapping for elements as DoFs */
521 std::vector< int > sorted_tmp_dofs( locsrc_soln_gdofs.size() );
522 std::copy( locsrc_soln_gdofs.begin(), locsrc_soln_gdofs.end(), sorted_tmp_dofs.begin() );
523 std::sort( sorted_tmp_dofs.begin(), sorted_tmp_dofs.end() );
524
525 for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
526 {
527 auto gdof = sorted_tmp_dofs[i];
528 assert( gdof > 0 );
529 srccol_gdofmap[i] = gdof - 1;
530 srccol_dtoc_dofmap[i] = i;
531 m_nTotDofs_Src++;
532 }
533 }
534 else
535 {
536 if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, false );
537 // Put these remap coefficients into the SparseMatrix map
538 for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
539 {
540 for( int p = 0; p < m_nDofsPEl_Src; p++ )
541 {
542 for( int q = 0; q < m_nDofsPEl_Src; q++ )
543 {
544 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
545 const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
546 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
547 {
548 m_nTotDofs_Src++;
549 dgll_cgll_col_ldofmap[localDOF] = true;
550 }
551 if( !isSrcContinuous ) m_nTotDofs_Src++;
552 assert( locsrc_soln_gdofs[offsetDOF] > 0 );
553 srccol_gdofmap[localDOF] = locsrc_soln_gdofs[offsetDOF] - 1;
554 srccol_dtoc_dofmap[offsetDOF] = localDOF;
555 }
556 }
557 }
558 }
559
560 int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
561 if( m_remapper->point_cloud_target )
562 {
563 assert( m_nDofsPEl_Dest == 1 );
564 row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
565 row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
566 tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
567 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] ) );
568 tgtTagSize = 1;
569 }
570 else
571 {
572 row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
573 row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
574 tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
575 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] ) );
576 }
577
578 // Now compute the mapping and store it for the target mesh
579 // To access the GID for each row: row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
580 m_nTotDofs_Dest = 0;
581 if( tgtdataGLLNodes == nullptr )
582 {
583 /* we only have a mapping for elements as DoFs */
584 std::vector< int > sorted_tmp_dofs( tgt_soln_gdofs.size() );
585 std::copy( tgt_soln_gdofs.begin(), tgt_soln_gdofs.end(), sorted_tmp_dofs.begin() );
586 std::sort( sorted_tmp_dofs.begin(), sorted_tmp_dofs.end() );
587
588 for( unsigned i = 0; i < row_gdofmap.size(); ++i )
589 {
590 auto gdof = sorted_tmp_dofs[i];
591 assert( gdof > 0 );
592 row_gdofmap[i] = gdof - 1;
593 row_dtoc_dofmap[i] = i;
594 if( vprint ) std::cout << "Row: " << i << ", " << row_gdofmap[i] << "\n";
595 m_nTotDofs_Dest++;
596 }
597 }
598 else
599 {
600 if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, false );
601 // Put these remap coefficients into the SparseMatrix map
602 for( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
603 {
604 for( int p = 0; p < m_nDofsPEl_Dest; p++ )
605 {
606 for( int q = 0; q < m_nDofsPEl_Dest; q++ )
607 {
608 const int localDOF = ( *tgtdataGLLNodes )[p][q][j] - 1;
609 const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
610 if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
611 {
612 m_nTotDofs_Dest++;
613 dgll_cgll_row_ldofmap[localDOF] = true;
614 }
615 if( !isTgtContinuous ) m_nTotDofs_Dest++;
616 assert( tgt_soln_gdofs[offsetDOF] > 0 );
617 row_gdofmap[localDOF] = tgt_soln_gdofs[offsetDOF] - 1;
618 row_dtoc_dofmap[offsetDOF] = localDOF;
619 if( vprint )
620 std::cout << "Row: " << offsetDOF << ", " << localDOF << ", " << row_gdofmap[offsetDOF] << ", "
621 << m_nTotDofs_Dest << "\n";
622 }
623 }
624 }
625 }
626
627 // Let us also allocate the local representation of the sparse matrix
628 #if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
629 if( vprint )
630 {
631 std::cout << "[" << rank << "]" << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_gdofmap.size()
632 << ", col = " << m_nTotDofs_Src << ", " << m_nTotDofs_SrcCov << ", " << col_gdofmap.size() << "\n";
633 // std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
634 }
635 #endif
636
637 // check monotonicity of row_gdofmap and col_gdofmap
638 #ifdef CHECK_INCREASING_DOF
639 for( size_t i = 0; i < row_gdofmap.size() - 1; i++ )
640 {
641 if( row_gdofmap[i] > row_gdofmap[i + 1] )
642 std::cout << " on rank " << rank << " in row_gdofmap[" << i << "]=" << row_gdofmap[i] << " > row_gdofmap["
643 << i + 1 << "]=" << row_gdofmap[i + 1] << " \n";
644 }
645 for( size_t i = 0; i < col_gdofmap.size() - 1; i++ )
646 {
647 if( col_gdofmap[i] > col_gdofmap[i + 1] )
648 std::cout << " on rank " << rank << " in col_gdofmap[" << i << "]=" << col_gdofmap[i] << " > col_gdofmap["
649 << i + 1 << "]=" << col_gdofmap[i + 1] << " \n";
650 }
651 #endif
652
653 return moab::MB_SUCCESS;
654 }
References MB_CHK_ERR, and MB_SUCCESS.
| moab::ErrorCode moab::TempestOnlineMap::SetDOFmapTags | ( | const std::string | srcDofTagName, |
| const std::string | tgtDofTagName | ||
| ) |
Store the tag names associated with global DoF ids for source and target meshes to be used for mapping.
| srcDofTagName | The tag name associated with global DoF ids for the source mesh |
| tgtDofTagName | The tag name associated with global DoF ids for the target mesh |
Definition at line 131 of file TempestOnlineMap.cpp.
133 {
134 moab::ErrorCode rval;
135
136 int tagSize = 0;
137 tagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
138 rval =
139 m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );
140
141 if( rval == moab::MB_TAG_NOT_FOUND && m_eInputType != DiscretizationType_FV )
142 {
143 MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for source mesh." );
144 }
145 else
146 MB_CHK_ERR( rval );
147
148 tagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
149 rval =
150 m_interface->tag_get_handle( tgtDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagDest, MB_TAG_ANY );
151 if( rval == moab::MB_TAG_NOT_FOUND && m_eOutputType != DiscretizationType_FV )
152 {
153 MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for target mesh." );
154 }
155 else
156 MB_CHK_ERR( rval );
157
158 return moab::MB_SUCCESS;
159 }
References ErrorCode, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_ANY, MB_TAG_NOT_FOUND, and MB_TYPE_INTEGER.
|
inline |
Definition at line 463 of file TempestOnlineMap.hpp.
464 { 465 m_meshInput = imesh; 466 };
References m_meshInput.
|
inline |
Set the number of Degrees-Of-Freedom per element on the source mesh.
Definition at line 600 of file TempestOnlineMap.hpp.
601 { 602 m_nDofsPEl_Src = ns; 603 }
|
private |
Definition at line 93 of file TempestOnlineMap.cpp.
94 {
95 if( m_meshInputCov )
96 {
97 std::vector< std::string > dimNames;
98 std::vector< int > dimSizes;
99 dimNames.push_back( "num_elem" );
100 dimSizes.push_back( m_meshInputCov->faces.size() );
101
102 this->InitializeSourceDimensions( dimNames, dimSizes );
103 }
104
105 if( m_meshOutput )
106 {
107 std::vector< std::string > dimNames;
108 std::vector< int > dimSizes;
109 dimNames.push_back( "num_elem" );
110 dimSizes.push_back( m_meshOutput->faces.size() );
111
112 this->InitializeTargetDimensions( dimNames, dimSizes );
113 }
114 }
Referenced by TempestOnlineMap().
|
private |
Parallel I/O with NetCDF to write out the SCRIP file from multiple processors.
Need to get the global maximum of number of vertices per element Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However, when writing this out, other processes may have a different size for the same array. This is hence a mess to consolidate in h5mtoscrip eventually.
Definition at line 829 of file TempestOnlineMapIO.cpp.
830 {
831 moab::ErrorCode rval;
832
833 /**
834 * Need to get the global maximum of number of vertices per element
835 * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
836 *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
837 *when writing this out, other processes may have a different size for the same array. This is
838 *hence a mess to consolidate in h5mtoscrip eventually.
839 **/
840
841 /* Let us compute all relevant data for the current original source mesh on the process */
842 DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
843 DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
844 DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
845 if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
846 {
847 this->InitializeCoordinatesFromMeshFV(
848 *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
849 ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
850 m_remapper->max_source_edges );
851
852 vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
853 for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
854 vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
855 }
856 else
857 {
858 DataArray3D< double > dataGLLJacobianSrc;
859 this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
860 dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
861
862 // Generate the continuous Jacobian for input mesh
863 GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
864
865 if( m_srcDiscType == DiscretizationType_CGLL )
866 {
867 GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
868 }
869 else
870 {
871 GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
872 }
873 }
874
875 if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
876 {
877 this->InitializeCoordinatesFromMeshFV(
878 *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
879 ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
880 m_remapper->max_target_edges );
881
882 vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
883 for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
884 vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
885 }
886 else
887 {
888 DataArray3D< double > dataGLLJacobianDest;
889 this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
890 dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
891
892 // Generate the continuous Jacobian for input mesh
893 GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
894
895 if( m_destDiscType == DiscretizationType_CGLL )
896 {
897 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
898 }
899 else
900 {
901 GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
902 }
903 }
904
905 moab::EntityHandle& m_meshOverlapSet = m_remapper->m_overlap_set;
906 int tot_src_ents = m_remapper->m_source_entities.size();
907 int tot_tgt_ents = m_remapper->m_target_entities.size();
908 int tot_src_size = dSourceCenterLon.GetRows();
909 int tot_tgt_size = m_dTargetCenterLon.GetRows();
910 int tot_vsrc_size = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
911 int tot_vtgt_size = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
912
913 const int weightMatNNZ = m_weightMatrix.nonZeros();
914 moab::Tag tagMapMetaData, tagMapIndexRow, tagMapIndexCol, tagMapValues, srcEleIDs, tgtEleIDs;
915 rval = m_interface->tag_get_handle( "SMAT_DATA", 13, moab::MB_TYPE_INTEGER, tagMapMetaData,
916 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
917 rval = m_interface->tag_get_handle( "SMAT_ROWS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexRow,
918 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
919 rval = m_interface->tag_get_handle( "SMAT_COLS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexCol,
920 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
921 rval = m_interface->tag_get_handle( "SMAT_VALS", weightMatNNZ, moab::MB_TYPE_DOUBLE, tagMapValues,
922 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
923 rval = m_interface->tag_get_handle( "SourceGIDS", tot_src_size, moab::MB_TYPE_INTEGER, srcEleIDs,
924 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
925 rval = m_interface->tag_get_handle( "TargetGIDS", tot_tgt_size, moab::MB_TYPE_INTEGER, tgtEleIDs,
926 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
927 moab::Tag srcAreaValues, tgtAreaValues;
928 rval = m_interface->tag_get_handle( "SourceAreas", tot_src_size, moab::MB_TYPE_DOUBLE, srcAreaValues,
929 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
930 rval = m_interface->tag_get_handle( "TargetAreas", tot_tgt_size, moab::MB_TYPE_DOUBLE, tgtAreaValues,
931 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
932 moab::Tag tagSrcCoordsCLon, tagSrcCoordsCLat, tagTgtCoordsCLon, tagTgtCoordsCLat;
933 rval = m_interface->tag_get_handle( "SourceCoordCenterLon", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLon,
934 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
935 rval = m_interface->tag_get_handle( "SourceCoordCenterLat", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLat,
936 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
937 rval = m_interface->tag_get_handle( "TargetCoordCenterLon", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLon,
938 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
939 rval = m_interface->tag_get_handle( "TargetCoordCenterLat", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLat,
940 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
941 moab::Tag tagSrcCoordsVLon, tagSrcCoordsVLat, tagTgtCoordsVLon, tagTgtCoordsVLat;
942 rval = m_interface->tag_get_handle( "SourceCoordVertexLon", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLon,
943 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
944 rval = m_interface->tag_get_handle( "SourceCoordVertexLat", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLat,
945 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
946 rval = m_interface->tag_get_handle( "TargetCoordVertexLon", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLon,
947 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
948 rval = m_interface->tag_get_handle( "TargetCoordVertexLat", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLat,
949 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
950 moab::Tag srcMaskValues, tgtMaskValues;
951 if( m_iSourceMask.IsAttached() )
952 {
953 rval = m_interface->tag_get_handle( "SourceMask", m_iSourceMask.GetRows(), moab::MB_TYPE_INTEGER, srcMaskValues,
954 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
955 }
956 if( m_iTargetMask.IsAttached() )
957 {
958 rval = m_interface->tag_get_handle( "TargetMask", m_iTargetMask.GetRows(), moab::MB_TYPE_INTEGER, tgtMaskValues,
959 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
960 }
961
962 std::vector< int > smatrowvals( weightMatNNZ ), smatcolvals( weightMatNNZ );
963 std::vector< double > smatvals( weightMatNNZ );
964 // const double* smatvals = m_weightMatrix.valuePtr();
965 // Loop over the matrix entries and find the max global ID for rows and columns
966 for( int k = 0, offset = 0; k < m_weightMatrix.outerSize(); ++k )
967 {
968 for( moab::TempestOnlineMap::WeightMatrix::InnerIterator it( m_weightMatrix, k ); it; ++it, ++offset )
969 {
970 smatrowvals[offset] = this->GetRowGlobalDoF( it.row() );
971 smatcolvals[offset] = this->GetColGlobalDoF( it.col() );
972 smatvals[offset] = it.value();
973 }
974 }
975
976 /* Set the global IDs for the DoFs */
977 ////
978 // col_gdofmap [ col_ldofmap [ 0 : local_ndofs ] ] = GDOF
979 // row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
980 ////
981 int maxrow = 0, maxcol = 0;
982 std::vector< int > src_global_dofs( tot_src_size ), tgt_global_dofs( tot_tgt_size );
983 for( int i = 0; i < tot_src_size; ++i )
984 {
985 src_global_dofs[i] = srccol_gdofmap[i];
986 maxcol = ( src_global_dofs[i] > maxcol ) ? src_global_dofs[i] : maxcol;
987 }
988
989 for( int i = 0; i < tot_tgt_size; ++i )
990 {
991 tgt_global_dofs[i] = row_gdofmap[i];
992 maxrow = ( tgt_global_dofs[i] > maxrow ) ? tgt_global_dofs[i] : maxrow;
993 }
994
995 ///////////////////////////////////////////////////////////////////////////
996 // The metadata in H5M file contains the following data:
997 //
998 // 1. n_a: Total source entities: (number of elements in source mesh)
999 // 2. n_b: Total target entities: (number of elements in target mesh)
1000 // 3. nv_a: Max edge size of elements in source mesh
1001 // 4. nv_b: Max edge size of elements in target mesh
1002 // 5. maxrows: Number of rows in remap weight matrix
1003 // 6. maxcols: Number of cols in remap weight matrix
1004 // 7. nnz: Number of total nnz in sparse remap weight matrix
1005 // 8. np_a: The order of the field description on the source mesh: >= 1
1006 // 9. np_b: The order of the field description on the target mesh: >= 1
1007 // 10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2 =
1008 // dGLL]
1009 // 11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2 =
1010 // dGLL]
1011 // 12. conserved: Flag to specify whether the remap operator has conservation constraints: [0,
1012 // 1]
1013 // 13. monotonicity: Flags to specify whether the remap operator has monotonicity constraints:
1014 // [0, 1, 2]
1015 //
1016 ///////////////////////////////////////////////////////////////////////////
1017 int map_disc_details[6];
1018 map_disc_details[0] = m_nDofsPEl_Src;
1019 map_disc_details[1] = m_nDofsPEl_Dest;
1020 map_disc_details[2] = ( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD
1021 ? 0
1022 : ( m_srcDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
1023 map_disc_details[3] = ( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD
1024 ? 0
1025 : ( m_destDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
1026 map_disc_details[4] = ( m_bConserved ? 1 : 0 );
1027 map_disc_details[5] = m_iMonotonicity;
1028
1029 #ifdef MOAB_HAVE_MPI
1030 int loc_smatmetadata[13] = { tot_src_ents,
1031 tot_tgt_ents,
1032 m_remapper->max_source_edges,
1033 m_remapper->max_target_edges,
1034 maxrow + 1,
1035 maxcol + 1,
1036 weightMatNNZ,
1037 map_disc_details[0],
1038 map_disc_details[1],
1039 map_disc_details[2],
1040 map_disc_details[3],
1041 map_disc_details[4],
1042 map_disc_details[5] };
1043 rval = m_interface->tag_set_data( tagMapMetaData, &m_meshOverlapSet, 1, &loc_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1044 int glb_smatmetadata[13] = { 0,
1045 0,
1046 0,
1047 0,
1048 0,
1049 0,
1050 0,
1051 map_disc_details[0],
1052 map_disc_details[1],
1053 map_disc_details[2],
1054 map_disc_details[3],
1055 map_disc_details[4],
1056 map_disc_details[5] };
1057 int loc_buf[7] = {
1058 tot_src_ents, tot_tgt_ents, weightMatNNZ, m_remapper->max_source_edges, m_remapper->max_target_edges,
1059 maxrow, maxcol };
1060 int glb_buf[4] = { 0, 0, 0, 0 };
1061 MPI_Reduce( &loc_buf[0], &glb_buf[0], 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
1062 glb_smatmetadata[0] = glb_buf[0];
1063 glb_smatmetadata[1] = glb_buf[1];
1064 glb_smatmetadata[6] = glb_buf[2];
1065 MPI_Reduce( &loc_buf[3], &glb_buf[0], 4, MPI_INT, MPI_MAX, 0, m_pcomm->comm() );
1066 glb_smatmetadata[2] = glb_buf[0];
1067 glb_smatmetadata[3] = glb_buf[1];
1068 glb_smatmetadata[4] = glb_buf[2];
1069 glb_smatmetadata[5] = glb_buf[3];
1070 #else
1071 int glb_smatmetadata[13] = { tot_src_ents,
1072 tot_tgt_ents,
1073 m_remapper->max_source_edges,
1074 m_remapper->max_target_edges,
1075 maxrow,
1076 maxcol,
1077 weightMatNNZ,
1078 map_disc_details[0],
1079 map_disc_details[1],
1080 map_disc_details[2],
1081 map_disc_details[3],
1082 map_disc_details[4],
1083 map_disc_details[5] };
1084 #endif
1085 // These values represent number of rows and columns. So should be 1-based.
1086 glb_smatmetadata[4]++;
1087 glb_smatmetadata[5]++;
1088
1089 if( this->is_root )
1090 {
1091 std::cout << " " << this->rank << " Writing remap weights with size [" << glb_smatmetadata[4] << " X "
1092 << glb_smatmetadata[5] << "] and NNZ = " << glb_smatmetadata[6] << std::endl;
1093 EntityHandle root_set = 0;
1094 rval = m_interface->tag_set_data( tagMapMetaData, &root_set, 1, &glb_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1095 }
1096
1097 int dsize;
1098 const int numval = weightMatNNZ;
1099 const void* smatrowvals_d = smatrowvals.data();
1100 const void* smatcolvals_d = smatcolvals.data();
1101 const void* smatvals_d = smatvals.data();
1102 rval = m_interface->tag_set_by_ptr( tagMapIndexRow, &m_meshOverlapSet, 1, &smatrowvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1103 rval = m_interface->tag_set_by_ptr( tagMapIndexCol, &m_meshOverlapSet, 1, &smatcolvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1104 rval = m_interface->tag_set_by_ptr( tagMapValues, &m_meshOverlapSet, 1, &smatvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1105
1106 /* Set the global IDs for the DoFs */
1107 const void* srceleidvals_d = src_global_dofs.data();
1108 const void* tgteleidvals_d = tgt_global_dofs.data();
1109 dsize = src_global_dofs.size();
1110 rval = m_interface->tag_set_by_ptr( srcEleIDs, &m_meshOverlapSet, 1, &srceleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1111 dsize = tgt_global_dofs.size();
1112 rval = m_interface->tag_set_by_ptr( tgtEleIDs, &m_meshOverlapSet, 1, &tgteleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1113
1114 /* Set the source and target areas */
1115 const void* srcareavals_d = vecSourceFaceArea;
1116 const void* tgtareavals_d = vecTargetFaceArea;
1117 dsize = tot_src_size;
1118 rval = m_interface->tag_set_by_ptr( srcAreaValues, &m_meshOverlapSet, 1, &srcareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1119 dsize = tot_tgt_size;
1120 rval = m_interface->tag_set_by_ptr( tgtAreaValues, &m_meshOverlapSet, 1, &tgtareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1121
1122 /* Set the coordinates for source and target center vertices */
1123 const void* srccoordsclonvals_d = &dSourceCenterLon[0];
1124 const void* srccoordsclatvals_d = &dSourceCenterLat[0];
1125 dsize = dSourceCenterLon.GetRows();
1126 rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLon, &m_meshOverlapSet, 1, &srccoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1127 rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLat, &m_meshOverlapSet, 1, &srccoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1128 const void* tgtcoordsclonvals_d = &m_dTargetCenterLon[0];
1129 const void* tgtcoordsclatvals_d = &m_dTargetCenterLat[0];
1130 dsize = vecTargetFaceArea.GetRows();
1131 rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLon, &m_meshOverlapSet, 1, &tgtcoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1132 rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLat, &m_meshOverlapSet, 1, &tgtcoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1133
1134 /* Set the coordinates for source and target element vertices */
1135 const void* srccoordsvlonvals_d = &( dSourceVertexLon[0][0] );
1136 const void* srccoordsvlatvals_d = &( dSourceVertexLat[0][0] );
1137 dsize = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
1138 rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLon, &m_meshOverlapSet, 1, &srccoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1139 rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLat, &m_meshOverlapSet, 1, &srccoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1140 const void* tgtcoordsvlonvals_d = &( m_dTargetVertexLon[0][0] );
1141 const void* tgtcoordsvlatvals_d = &( m_dTargetVertexLat[0][0] );
1142 dsize = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
1143 rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLon, &m_meshOverlapSet, 1, &tgtcoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1144 rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLat, &m_meshOverlapSet, 1, &tgtcoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1145
1146 /* Set the masks for source and target meshes if available */
1147 if( m_iSourceMask.IsAttached() )
1148 {
1149 const void* srcmaskvals_d = m_iSourceMask;
1150 dsize = m_iSourceMask.GetRows();
1151 rval = m_interface->tag_set_by_ptr( srcMaskValues, &m_meshOverlapSet, 1, &srcmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1152 }
1153
1154 if( m_iTargetMask.IsAttached() )
1155 {
1156 const void* tgtmaskvals_d = m_iTargetMask;
1157 dsize = m_iTargetMask.GetRows();
1158 rval = m_interface->tag_set_by_ptr( tgtMaskValues, &m_meshOverlapSet, 1, &tgtmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1159 }
1160
1161 #ifdef MOAB_HAVE_MPI
1162 const char* writeOptions = ( this->size > 1 ? "PARALLEL=WRITE_PART" : "" );
1163 #else
1164 const char* writeOptions = "";
1165 #endif
1166
1167 // EntityHandle sets[3] = {m_remapper->m_source_set, m_remapper->m_target_set, m_remapper->m_overlap_set};
1168 EntityHandle sets[1] = { m_remapper->m_overlap_set };
1169 rval = m_interface->write_file( strOutputFile.c_str(), NULL, writeOptions, sets, 1 );MB_CHK_ERR( rval );
1170
1171 #ifdef WRITE_SCRIP_FILE
1172 sstr.str( "" );
1173 sstr << ctx.outFilename.substr( 0, lastindex ) << "_" << proc_id << ".nc";
1174 std::map< std::string, std::string > mapAttributes;
1175 mapAttributes["Creator"] = "MOAB mbtempest workflow";
1176 if( !ctx.proc_id ) std::cout << "Writing offline map to file: " << sstr.str() << std::endl;
1177 this->Write( strOutputFile.c_str(), mapAttributes, NcFile::Netcdf4 );
1178 sstr.str( "" );
1179 #endif
1180
1181 return moab::MB_SUCCESS;
1182 }
References ErrorCode, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TAG_VARLEN, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, moab::TempestRemapper::RLL, and size.
| moab::ErrorCode moab::TempestOnlineMap::WriteParallelMap | ( | const std::string & | strTarget, |
| const std::map< std::string, std::string > & | attrMap | ||
| ) |
Write the TempestOnlineMap to a parallel NetCDF file.
Definition at line 193 of file TempestOnlineMapIO.cpp.
195 {
196 moab::ErrorCode rval;
197
198 size_t lastindex = strFilename.find_last_of( "." );
199 std::string extension = strFilename.substr( lastindex + 1, strFilename.size() );
200
201 // Write the map file to disk in parallel
202 if( extension == "nc" )
203 {
204 /* Invoke the actual call to write the parallel map to disk in SCRIP format */
205 rval = this->WriteSCRIPMapFile( strFilename.c_str(), attrMap );MB_CHK_ERR( rval );
206 }
207 else
208 {
209 /* Write to the parallel H5M format */
210 rval = this->WriteHDF5MapFile( strFilename.c_str() );MB_CHK_ERR( rval );
211 }
212
213 return rval;
214 }
References ErrorCode, and MB_CHK_ERR.
Referenced by main().
|
private |
Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections.
Parallel I/O with HDF5 to write out the remapping weights from multiple processors.
Need to get the global maximum of number of vertices per element Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However, when writing this out, other processes may have a different size for the same array. This is hence a mess to consolidate in h5mtoscrip eventually.
Definition at line 218 of file TempestOnlineMapIO.cpp.
220 {
221 NcError error( NcError::silent_nonfatal );
222
223 #ifdef MOAB_HAVE_NETCDFPAR
224 bool is_independent = true;
225 ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcFile::Replace, NcFile::Netcdf4 );
226 // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
227 #else
228 NcFile ncMap( strFilename.c_str(), NcFile::Replace );
229 #endif
230
231 if( !ncMap.is_valid() )
232 {
233 _EXCEPTION1( "Unable to open output map file \"%s\"", strFilename.c_str() );
234 }
235
236 // Attributes
237 // ncMap.add_att( "Title", "MOAB-TempestRemap Online Regridding Weight Generator" );
238 auto it = attrMap.begin();
239 while( it != attrMap.end() )
240 {
241 // set the map attributes
242 ncMap.add_att( it->first.c_str(), it->second.c_str() );
243 // increment iterator
244 it++;
245 }
246
247 /**
248 * Need to get the global maximum of number of vertices per element
249 * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
250 *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
251 *when writing this out, other processes may have a different size for the same array. This is
252 *hence a mess to consolidate in h5mtoscrip eventually.
253 **/
254
255 /* Let us compute all relevant data for the current original source mesh on the process */
256 DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
257 DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
258 DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
259 if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
260 {
261 this->InitializeCoordinatesFromMeshFV(
262 *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
263 ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
264 m_remapper->max_source_edges );
265
266 vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
267 for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
268 vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
269 }
270 else
271 {
272 DataArray3D< double > dataGLLJacobianSrc;
273 this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
274 dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
275
276 // Generate the continuous Jacobian for input mesh
277 GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
278
279 if( m_srcDiscType == DiscretizationType_CGLL )
280 {
281 GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
282 }
283 else
284 {
285 GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
286 }
287 }
288
289 if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
290 {
291 this->InitializeCoordinatesFromMeshFV(
292 *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
293 ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
294 m_remapper->max_target_edges );
295
296 vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
297 for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
298 {
299 vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
300 }
301 }
302 else
303 {
304 DataArray3D< double > dataGLLJacobianDest;
305 this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
306 dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
307
308 // Generate the continuous Jacobian for input mesh
309 GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
310
311 if( m_destDiscType == DiscretizationType_CGLL )
312 {
313 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
314 }
315 else
316 {
317 GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
318 }
319 }
320
321 // Map dimensions
322 unsigned nA = ( vecSourceFaceArea.GetRows() );
323 unsigned nB = ( vecTargetFaceArea.GetRows() );
324
325 std::vector< int > masksA, masksB;
326 ErrorCode rval = m_remapper->GetIMasks( moab::Remapper::SourceMesh, masksA );MB_CHK_SET_ERR( rval, "Trouble getting masks for source" );
327 rval = m_remapper->GetIMasks( moab::Remapper::TargetMesh, masksB );MB_CHK_SET_ERR( rval, "Trouble getting masks for target" );
328
329 // Number of nodes per Face
330 int nSourceNodesPerFace = dSourceVertexLon.GetColumns();
331 int nTargetNodesPerFace = dTargetVertexLon.GetColumns();
332
333 // if source or target cells have triangles at poles, center of those triangles need to come from
334 // the original quad, not from center in 3d, converted to 2d again
335 // start copy OnlineMap.cpp tempestremap
336 // right now, do this only for source mesh; copy the logic for target mesh
337 for( unsigned i = 0; i < nA; i++ )
338 {
339 const Face& face = m_meshInput->faces[i];
340
341 int nNodes = face.edges.size();
342 int indexNodeAtPole = -1;
343 if( 3 == nNodes ) // check if one node at the poles
344 {
345 for( int j = 0; j < nNodes; j++ )
346 if( fabs( fabs( dSourceVertexLat[i][j] ) - 90.0 ) < 1.0e-12 )
347 {
348 indexNodeAtPole = j;
349 break;
350 }
351 }
352 if( indexNodeAtPole < 0 ) continue; // continue i loop, do nothing
353 // recompute center of cell, from 3d data; add one 2 nodes at pole, and average
354 int nodeAtPole = face[indexNodeAtPole]; // use the overloaded operator
355 Node nodePole = m_meshInput->nodes[nodeAtPole];
356 Node newCenter = nodePole * 2;
357 for( int j = 1; j < nNodes; j++ )
358 {
359 int indexi = ( indexNodeAtPole + j ) % nNodes; // nNodes is 3 !
360 const Node& node = m_meshInput->nodes[face[indexi]];
361 newCenter = newCenter + node;
362 }
363 newCenter = newCenter * 0.25;
364 newCenter = newCenter.Normalized();
365
366 #ifdef VERBOSE
367 double iniLon = dSourceCenterLon[i], iniLat = dSourceCenterLat[i];
368 #endif
369 // dSourceCenterLon, dSourceCenterLat
370 XYZtoRLL_Deg( newCenter.x, newCenter.y, newCenter.z, dSourceCenterLon[i], dSourceCenterLat[i] );
371 #ifdef VERBOSE
372 std::cout << " modify center of triangle from " << iniLon << " " << iniLat << " to " << dSourceCenterLon[i]
373 << " " << dSourceCenterLat[i] << "\n";
374 #endif
375 }
376
377 // first move data if in parallel
378 #if defined( MOAB_HAVE_MPI )
379 int max_row_dof, max_col_dof; // output; arrays will be re-distributed in chunks [maxdof/size]
380 // if (size > 1)
381 {
382 int ierr = rearrange_arrays_by_dofs( srccol_gdofmap, vecSourceFaceArea, dSourceCenterLon, dSourceCenterLat,
383 dSourceVertexLon, dSourceVertexLat, masksA, nA, nSourceNodesPerFace,
384 max_col_dof ); // now nA will be close to maxdof/size
385 if( ierr != 0 )
386 {
387 _EXCEPTION1( "Unable to arrange source data %d ", nA );
388 }
389 // rearrange target data: (nB)
390 //
391 ierr = rearrange_arrays_by_dofs( row_gdofmap, vecTargetFaceArea, dTargetCenterLon, dTargetCenterLat,
392 dTargetVertexLon, dTargetVertexLat, masksB, nB, nTargetNodesPerFace,
393 max_row_dof ); // now nA will be close to maxdof/size
394 if( ierr != 0 )
395 {
396 _EXCEPTION1( "Unable to arrange target data %d ", nB );
397 }
398 }
399 #endif
400
401 // Number of non-zeros in the remap matrix operator
402 int nS = m_weightMatrix.nonZeros();
403
404 #if defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_NETCDFPAR )
405 int locbuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
406 int offbuf[3] = { 0, 0, 0 };
407 int globuf[5] = { 0, 0, 0, 0, 0 };
408 MPI_Scan( locbuf, offbuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
409 MPI_Allreduce( locbuf, globuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
410 MPI_Allreduce( &locbuf[3], &globuf[3], 2, MPI_INT, MPI_MAX, m_pcomm->comm() );
411
412 // MPI_Scan is inclusive of data in current rank; modify accordingly.
413 offbuf[0] -= nA;
414 offbuf[1] -= nB;
415 offbuf[2] -= nS;
416
417 #else
418 int offbuf[3] = { 0, 0, 0 };
419 int globuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
420 #endif
421
422 std::vector< std::string > srcdimNames, tgtdimNames;
423 std::vector< int > srcdimSizes, tgtdimSizes;
424 {
425 if( m_remapper->m_source_type == moab::TempestRemapper::RLL && m_remapper->m_source_metadata.size() )
426 {
427 srcdimNames.push_back( "lat" );
428 srcdimNames.push_back( "lon" );
429 srcdimSizes.resize( 2, 0 );
430 srcdimSizes[0] = m_remapper->m_source_metadata[0];
431 srcdimSizes[1] = m_remapper->m_source_metadata[1];
432 }
433 else
434 {
435 srcdimNames.push_back( "num_elem" );
436 srcdimSizes.push_back( globuf[0] );
437 }
438
439 if( m_remapper->m_target_type == moab::TempestRemapper::RLL && m_remapper->m_target_metadata.size() )
440 {
441 tgtdimNames.push_back( "lat" );
442 tgtdimNames.push_back( "lon" );
443 tgtdimSizes.resize( 2, 0 );
444 tgtdimSizes[0] = m_remapper->m_target_metadata[0];
445 tgtdimSizes[1] = m_remapper->m_target_metadata[1];
446 }
447 else
448 {
449 tgtdimNames.push_back( "num_elem" );
450 tgtdimSizes.push_back( globuf[1] );
451 }
452 }
453
454 // Write output dimensions entries
455 unsigned nSrcGridDims = ( srcdimSizes.size() );
456 unsigned nDstGridDims = ( tgtdimSizes.size() );
457
458 NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
459 NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
460
461 NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
462 NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
463
464 #ifdef MOAB_HAVE_NETCDFPAR
465 ncMap.enable_var_par_access( varSrcGridDims, is_independent );
466 ncMap.enable_var_par_access( varDstGridDims, is_independent );
467 #endif
468
469 // write dimension names
470 {
471 char szDim[64];
472 for( unsigned i = 0; i < srcdimSizes.size(); i++ )
473 {
474 varSrcGridDims->set_cur( nSrcGridDims - i - 1 );
475 varSrcGridDims->put( &( srcdimSizes[nSrcGridDims - i - 1] ), 1 );
476 }
477
478 for( unsigned i = 0; i < srcdimSizes.size(); i++ )
479 {
480 snprintf( szDim, 64, "name%i", i );
481 varSrcGridDims->add_att( szDim, srcdimNames[nSrcGridDims - i - 1].c_str() );
482 }
483
484 for( unsigned i = 0; i < tgtdimSizes.size(); i++ )
485 {
486 varDstGridDims->set_cur( nDstGridDims - i - 1 );
487 varDstGridDims->put( &( tgtdimSizes[nDstGridDims - i - 1] ), 1 );
488 }
489
490 for( unsigned i = 0; i < tgtdimSizes.size(); i++ )
491 {
492 snprintf( szDim, 64, "name%i", i );
493 varDstGridDims->add_att( szDim, tgtdimNames[nDstGridDims - i - 1].c_str() );
494 }
495 }
496
497 // Source and Target mesh resolutions
498 NcDim* dimNA = ncMap.add_dim( "n_a", globuf[0] );
499 NcDim* dimNB = ncMap.add_dim( "n_b", globuf[1] );
500
501 // Number of nodes per Face
502 NcDim* dimNVA = ncMap.add_dim( "nv_a", globuf[3] );
503 NcDim* dimNVB = ncMap.add_dim( "nv_b", globuf[4] );
504
505 // Write coordinates
506 NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA );
507 NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB );
508
509 NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA );
510 NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB );
511
512 NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA, dimNVA );
513 NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB, dimNVB );
514
515 NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA, dimNVA );
516 NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB, dimNVB );
517
518 // Write masks
519 NcVar* varMaskA = ncMap.add_var( "mask_a", ncInt, dimNA );
520 NcVar* varMaskB = ncMap.add_var( "mask_b", ncInt, dimNB );
521
522 #ifdef MOAB_HAVE_NETCDFPAR
523 ncMap.enable_var_par_access( varYCA, is_independent );
524 ncMap.enable_var_par_access( varYCB, is_independent );
525 ncMap.enable_var_par_access( varXCA, is_independent );
526 ncMap.enable_var_par_access( varXCB, is_independent );
527 ncMap.enable_var_par_access( varYVA, is_independent );
528 ncMap.enable_var_par_access( varYVB, is_independent );
529 ncMap.enable_var_par_access( varXVA, is_independent );
530 ncMap.enable_var_par_access( varXVB, is_independent );
531 ncMap.enable_var_par_access( varMaskA, is_independent );
532 ncMap.enable_var_par_access( varMaskB, is_independent );
533 #endif
534
535 varYCA->add_att( "units", "degrees" );
536 varYCB->add_att( "units", "degrees" );
537
538 varXCA->add_att( "units", "degrees" );
539 varXCB->add_att( "units", "degrees" );
540
541 varYVA->add_att( "units", "degrees" );
542 varYVB->add_att( "units", "degrees" );
543
544 varXVA->add_att( "units", "degrees" );
545 varXVB->add_att( "units", "degrees" );
546
547 // Verify dimensionality
548 if( dSourceCenterLon.GetRows() != nA )
549 {
550 _EXCEPTIONT( "Mismatch between dSourceCenterLon and nA" );
551 }
552 if( dSourceCenterLat.GetRows() != nA )
553 {
554 _EXCEPTIONT( "Mismatch between dSourceCenterLat and nA" );
555 }
556 if( dTargetCenterLon.GetRows() != nB )
557 {
558 _EXCEPTIONT( "Mismatch between dTargetCenterLon and nB" );
559 }
560 if( dTargetCenterLat.GetRows() != nB )
561 {
562 _EXCEPTIONT( "Mismatch between dTargetCenterLat and nB" );
563 }
564 if( dSourceVertexLon.GetRows() != nA )
565 {
566 _EXCEPTIONT( "Mismatch between dSourceVertexLon and nA" );
567 }
568 if( dSourceVertexLat.GetRows() != nA )
569 {
570 _EXCEPTIONT( "Mismatch between dSourceVertexLat and nA" );
571 }
572 if( dTargetVertexLon.GetRows() != nB )
573 {
574 _EXCEPTIONT( "Mismatch between dTargetVertexLon and nB" );
575 }
576 if( dTargetVertexLat.GetRows() != nB )
577 {
578 _EXCEPTIONT( "Mismatch between dTargetVertexLat and nB" );
579 }
580
581 varYCA->set_cur( (long)offbuf[0] );
582 varYCA->put( &( dSourceCenterLat[0] ), nA );
583 varYCB->set_cur( (long)offbuf[1] );
584 varYCB->put( &( dTargetCenterLat[0] ), nB );
585
586 varXCA->set_cur( (long)offbuf[0] );
587 varXCA->put( &( dSourceCenterLon[0] ), nA );
588 varXCB->set_cur( (long)offbuf[1] );
589 varXCB->put( &( dTargetCenterLon[0] ), nB );
590
591 varYVA->set_cur( (long)offbuf[0] );
592 varYVA->put( &( dSourceVertexLat[0][0] ), nA, nSourceNodesPerFace );
593 varYVB->set_cur( (long)offbuf[1] );
594 varYVB->put( &( dTargetVertexLat[0][0] ), nB, nTargetNodesPerFace );
595
596 varXVA->set_cur( (long)offbuf[0] );
597 varXVA->put( &( dSourceVertexLon[0][0] ), nA, nSourceNodesPerFace );
598 varXVB->set_cur( (long)offbuf[1] );
599 varXVB->put( &( dTargetVertexLon[0][0] ), nB, nTargetNodesPerFace );
600
601 varMaskA->set_cur( (long)offbuf[0] );
602 varMaskA->put( &( masksA[0] ), nA );
603 varMaskB->set_cur( (long)offbuf[1] );
604 varMaskB->put( &( masksB[0] ), nB );
605
606 // Write areas
607 NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
608 #ifdef MOAB_HAVE_NETCDFPAR
609 ncMap.enable_var_par_access( varAreaA, is_independent );
610 #endif
611 varAreaA->set_cur( (long)offbuf[0] );
612 varAreaA->put( &( vecSourceFaceArea[0] ), nA );
613
614 NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
615 #ifdef MOAB_HAVE_NETCDFPAR
616 ncMap.enable_var_par_access( varAreaB, is_independent );
617 #endif
618 varAreaB->set_cur( (long)offbuf[1] );
619 varAreaB->put( &( vecTargetFaceArea[0] ), nB );
620
621 // Write SparseMatrix entries
622 DataArray1D< int > vecRow( nS );
623 DataArray1D< int > vecCol( nS );
624 DataArray1D< double > vecS( nS );
625 DataArray1D< double > dFracA( nA );
626 DataArray1D< double > dFracB( nB );
627
628 moab::TupleList tlValRow, tlValCol;
629 unsigned numr = 1; //
630 // value has to be sent to processor row/nB for for fracA and col/nA for fracB
631 // vecTargetArea (indexRow ) has to be sent for fracA (index col?)
632 // vecTargetFaceArea will have to be sent to col index, with its index !
633 tlValRow.initialize( 2, 0, 0, numr, nS ); // to proc(row), global row , value
634 tlValCol.initialize( 3, 0, 0, numr, nS ); // to proc(col), global row / col, value
635 tlValRow.enableWriteAccess();
636 tlValCol.enableWriteAccess();
637 /*
638 dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
639 dFracB[ row ] += val ;
640 */
641 int offset = 0;
642 #if defined( MOAB_HAVE_MPI )
643 int nAbase = ( max_col_dof + 1 ) / size; // it is nA, except last rank ( == size - 1 )
644 int nBbase = ( max_row_dof + 1 ) / size; // it is nB, except last rank ( == size - 1 )
645 #endif
646 for( int i = 0; i < m_weightMatrix.outerSize(); ++i )
647 {
648 for( WeightMatrix::InnerIterator it( m_weightMatrix, i ); it; ++it )
649 {
650 vecRow[offset] = 1 + this->GetRowGlobalDoF( it.row() ); // row index
651 vecCol[offset] = 1 + this->GetColGlobalDoF( it.col() ); // col index
652 vecS[offset] = it.value(); // value
653
654 #if defined( MOAB_HAVE_MPI )
655 {
656 // value M(row, col) will contribute to procRow and procCol values for fracA and fracB
657 int procRow = ( vecRow[offset] - 1 ) / nBbase;
658 if( procRow >= size ) procRow = size - 1;
659 int procCol = ( vecCol[offset] - 1 ) / nAbase;
660 if( procCol >= size ) procCol = size - 1;
661 int nrInd = tlValRow.get_n();
662 tlValRow.vi_wr[2 * nrInd] = procRow;
663 tlValRow.vi_wr[2 * nrInd + 1] = vecRow[offset] - 1;
664 tlValRow.vr_wr[nrInd] = vecS[offset];
665 tlValRow.inc_n();
666 int ncInd = tlValCol.get_n();
667 tlValCol.vi_wr[3 * ncInd] = procCol;
668 tlValCol.vi_wr[3 * ncInd + 1] = vecRow[offset] - 1;
669 tlValCol.vi_wr[3 * ncInd + 2] = vecCol[offset] - 1; // this is column
670 tlValCol.vr_wr[ncInd] = vecS[offset];
671 tlValCol.inc_n();
672 }
673
674 #endif
675 offset++;
676 }
677 }
678 #if defined( MOAB_HAVE_MPI )
679 // need to send values for their row and col processors, to compute fractions there
680 // now do the heavy communication
681 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValCol, 0 );
682 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValRow, 0 );
683
684 // we have now, for example, dFracB[ row ] += val ;
685 // so we know that on current task, we received tlValRow
686 // reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
687 // dFracB[ row ] += val ;
688 for( unsigned i = 0; i < tlValRow.get_n(); i++ )
689 {
690 // int fromProc = tlValRow.vi_wr[2 * i];
691 int gRowInd = tlValRow.vi_wr[2 * i + 1];
692 int localIndexRow = gRowInd - nBbase * rank; // modulo nBbase rank is from 0 to size - 1;
693 double wgt = tlValRow.vr_wr[i];
694 assert( localIndexRow >= 0 );
695 assert( nB - localIndexRow > 0 );
696 dFracB[localIndexRow] += wgt;
697 }
698 // to compute dFracA we need vecTargetFaceArea[ row ]; we know the row, and we can get the proc we need it from
699
700 std::set< int > neededRows;
701 for( unsigned i = 0; i < tlValCol.get_n(); i++ )
702 {
703 int rRowInd = tlValCol.vi_wr[3 * i + 1];
704 neededRows.insert( rRowInd );
705 // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
706 }
707 moab::TupleList tgtAreaReq;
708 tgtAreaReq.initialize( 2, 0, 0, 0, neededRows.size() );
709 tgtAreaReq.enableWriteAccess();
710 for( std::set< int >::iterator sit = neededRows.begin(); sit != neededRows.end(); sit++ )
711 {
712 int neededRow = *sit;
713 int procRow = neededRow / nBbase;
714 if( procRow >= size ) procRow = size - 1;
715 int nr = tgtAreaReq.get_n();
716 tgtAreaReq.vi_wr[2 * nr] = procRow;
717 tgtAreaReq.vi_wr[2 * nr + 1] = neededRow;
718 tgtAreaReq.inc_n();
719 }
720
721 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaReq, 0 );
722 // we need to send back the tgtArea corresponding to row
723 moab::TupleList tgtAreaInfo; // load it with tgtArea at row
724 tgtAreaInfo.initialize( 2, 0, 0, 1, tgtAreaReq.get_n() );
725 tgtAreaInfo.enableWriteAccess();
726 for( unsigned i = 0; i < tgtAreaReq.get_n(); i++ )
727 {
728 int from_proc = tgtAreaReq.vi_wr[2 * i];
729 int row = tgtAreaReq.vi_wr[2 * i + 1];
730 int locaIndexRow = row - rank * nBbase;
731 double areaToSend = vecTargetFaceArea[locaIndexRow];
732 // int remoteIndex = tgtAreaReq.vi_wr[3*i + 2] ;
733
734 tgtAreaInfo.vi_wr[2 * i] = from_proc; // send back requested info
735 tgtAreaInfo.vi_wr[2 * i + 1] = row;
736 tgtAreaInfo.vr_wr[i] = areaToSend; // this will be tgt area at row
737 tgtAreaInfo.inc_n();
738 }
739 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaInfo, 0 );
740
741 std::map< int, double > areaAtRow;
742 for( unsigned i = 0; i < tgtAreaInfo.get_n(); i++ )
743 {
744 // we have received from proc, value for row !
745 int row = tgtAreaInfo.vi_wr[2 * i + 1];
746 areaAtRow[row] = tgtAreaInfo.vr_wr[i];
747 }
748
749 // we have now for rows the
750 // it is ordered by index, so:
751 // now compute reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
752 // tgtAreaInfo will have at index i the area we need (from row)
753 // there should be an easier way :(
754 for( unsigned i = 0; i < tlValCol.get_n(); i++ )
755 {
756 int rRowInd = tlValCol.vi_wr[3 * i + 1];
757 int colInd = tlValCol.vi_wr[3 * i + 2];
758 double val = tlValCol.vr_wr[i];
759 int localColInd = colInd - rank * nAbase; // < local nA
760 // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
761 auto itMap = areaAtRow.find( rRowInd ); // it should be different from end
762 if( itMap != areaAtRow.end() )
763 {
764 double areaRow = itMap->second; // we fished a lot for this !
765 dFracA[localColInd] += val / vecSourceFaceArea[localColInd] * areaRow;
766 }
767 }
768
769 #endif
770 // Load in data
771 NcDim* dimNS = ncMap.add_dim( "n_s", globuf[2] );
772
773 NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
774 NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
775 NcVar* varS = ncMap.add_var( "S", ncDouble, dimNS );
776 #ifdef MOAB_HAVE_NETCDFPAR
777 ncMap.enable_var_par_access( varRow, is_independent );
778 ncMap.enable_var_par_access( varCol, is_independent );
779 ncMap.enable_var_par_access( varS, is_independent );
780 #endif
781
782 varRow->set_cur( (long)offbuf[2] );
783 varRow->put( vecRow, nS );
784
785 varCol->set_cur( (long)offbuf[2] );
786 varCol->put( vecCol, nS );
787
788 varS->set_cur( (long)offbuf[2] );
789 varS->put( &( vecS[0] ), nS );
790
791 // Calculate and write fractional coverage arrays
792 NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
793 #ifdef MOAB_HAVE_NETCDFPAR
794 ncMap.enable_var_par_access( varFracA, is_independent );
795 #endif
796 varFracA->add_att( "name", "fraction of target coverage of source dof" );
797 varFracA->add_att( "units", "unitless" );
798 varFracA->set_cur( (long)offbuf[0] );
799 varFracA->put( &( dFracA[0] ), nA );
800
801 NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
802 #ifdef MOAB_HAVE_NETCDFPAR
803 ncMap.enable_var_par_access( varFracB, is_independent );
804 #endif
805 varFracB->add_att( "name", "fraction of source coverage of target dof" );
806 varFracB->add_att( "units", "unitless" );
807 varFracB->set_cur( (long)offbuf[1] );
808 varFracB->put( &( dFracB[0] ), nB );
809
810 // Add global attributes
811 // std::map<std::string, std::string>::const_iterator iterAttributes =
812 // mapAttributes.begin();
813 // for (; iterAttributes != mapAttributes.end(); iterAttributes++) {
814 // ncMap.add_att(
815 // iterAttributes->first.c_str(),
816 // iterAttributes->second.c_str());
817 // }
818
819 ncMap.close();
820
821 #ifdef VERBOSE
822 serializeSparseMatrix( m_weightMatrix, "map_operator_" + std::to_string( rank ) + ".txt" );
823 #endif
824 return moab::MB_SUCCESS;
825 }
References moab::TupleList::enableWriteAccess(), moab::error(), ErrorCode, moab::TupleList::get_n(), moab::TupleList::inc_n(), moab::TupleList::initialize(), MB_CHK_SET_ERR, MB_SUCCESS, nr, moab::TempestRemapper::RLL, size, moab::Remapper::SourceMesh, moab::Remapper::TargetMesh, moab::TupleList::vi_wr, and moab::TupleList::vr_wr.
|
private |
Definition at line 539 of file TempestOnlineMap.hpp.
|
private |
Definition at line 536 of file TempestOnlineMap.hpp.
Referenced by fill_col_ids(), and LinearRemapNN_MOAB().
|
private |
Definition at line 541 of file TempestOnlineMap.hpp.
|
private |
Definition at line 544 of file TempestOnlineMap.hpp.
|
private |
Definition at line 544 of file TempestOnlineMap.hpp.
|
private |
Definition at line 544 of file TempestOnlineMap.hpp.
|
private |
Definition at line 559 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 559 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 551 of file TempestOnlineMap.hpp.
|
private |
Definition at line 545 of file TempestOnlineMap.hpp.
|
private |
Definition at line 535 of file TempestOnlineMap.hpp.
|
private |
The original tag data and local to global DoF mapping to associate matrix values to solution
Definition at line 535 of file TempestOnlineMap.hpp.
|
private |
Definition at line 550 of file TempestOnlineMap.hpp.
|
private |
Definition at line 550 of file TempestOnlineMap.hpp.
|
private |
Definition at line 552 of file TempestOnlineMap.hpp.
|
private |
Definition at line 542 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
The reference to the moab::Core object that contains source/target and overlap sets.
Definition at line 523 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 554 of file TempestOnlineMap.hpp.
Referenced by SetMeshInput(), and TempestOnlineMap().
|
private |
Definition at line 555 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 556 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 557 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 549 of file TempestOnlineMap.hpp.
|
private |
Definition at line 549 of file TempestOnlineMap.hpp.
|
private |
Definition at line 546 of file TempestOnlineMap.hpp.
Referenced by LinearRemapNN_MOAB().
|
private |
Definition at line 546 of file TempestOnlineMap.hpp.
|
private |
Definition at line 546 of file TempestOnlineMap.hpp.
Referenced by LinearRemapNN_MOAB().
|
private |
Definition at line 542 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
The fundamental remapping operator object.
Definition at line 510 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 545 of file TempestOnlineMap.hpp.
|
private |
Definition at line 560 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 539 of file TempestOnlineMap.hpp.
|
private |
Definition at line 536 of file TempestOnlineMap.hpp.
Referenced by fill_row_ids(), and LinearRemapNN_MOAB().
|
private |
Definition at line 541 of file TempestOnlineMap.hpp.
|
private |
Definition at line 560 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
|
private |
Definition at line 539 of file TempestOnlineMap.hpp.
|
private |
Definition at line 536 of file TempestOnlineMap.hpp.
1.9.1.