From f17a9b61b986633ecd675fee63a52564c0c34087 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Thu, 12 May 2022 14:37:33 -0600 Subject: [PATCH 1/7] TrilinosCouplings: Adding mesh statistics --- .../scaling/TrilinosCouplings_Statistics.hpp | 2 +- .../examples/scaling/example_Poisson_stk.cpp | 145 +++++++++++++++++- 2 files changed, 140 insertions(+), 7 deletions(-) diff --git a/packages/trilinoscouplings/examples/scaling/TrilinosCouplings_Statistics.hpp b/packages/trilinoscouplings/examples/scaling/TrilinosCouplings_Statistics.hpp index ba570d4f7c56..1fea0b302a10 100644 --- a/packages/trilinoscouplings/examples/scaling/TrilinosCouplings_Statistics.hpp +++ b/packages/trilinoscouplings/examples/scaling/TrilinosCouplings_Statistics.hpp @@ -144,7 +144,7 @@ class MachineLearningStatistics_Hex3D { int node3 = edgeToNode(edge, 0); int node4 = edgeToNode(edge, 1); - int numElems = elemToEdge.dimension(0); + int numElems = elemToNode.dimension(0); for(int i=0; i +#include +#include #include "TrilinosCouplings_config.h" @@ -167,6 +169,12 @@ #include "stk_mesh/base/FieldBase.hpp" // for field_data, etc #include "stk_mesh/base/Types.hpp" // for BucketVector, EntityId +#ifdef HAVE_TRILINOSCOUPLINGS_MUELU +#include "TrilinosCouplings_Statistics.hpp" +#else +#error "BORK" +#endif + /*********************************************************/ /* Typedefs */ /*********************************************************/ @@ -481,16 +489,20 @@ int main(int argc, char *argv[]) { for (int i=0; i CoordFieldType; // get coordinates field @@ -537,6 +549,17 @@ int main(int argc, char *argv[]) { } // end loop over mesh parts int numNodesPerElem = cellType.getNodeCount(); + int numEdgesPerElem = cellType.getEdgeCount(); + + if(MyPID==0) { + std::cout<<"Cell Topology: "< MLStatistics(numGlobalElements); + bool do_statistics = !strcmp(cellType.getName(),"Hexahedron_8"); + std::cout<<"do_statistics = "< elemToNode(numLocalElems,numNodesPerElem); + Intrepid::FieldContainer elemToEdge(numLocalElems,numEdgesPerElem); + Intrepid::FieldContainer nodeCoord (numLocalNodes, spaceDim); + Intrepid::FieldContainer sigmaVal(numLocalElems); + + int elem_ct=0; + std::map,int> local_edge_hash; + std::vector > edge_vector; + + for (size_t bucketIndex = 0; bucketIndex < localElementBuckets.size(); ++bucketIndex) { + stk::mesh::Bucket &elemBucket = *localElementBuckets[bucketIndex]; + for (size_t elemIndex = 0; elemIndex < elemBucket.size(); ++elemIndex) { + stk::mesh::Entity elem = elemBucket[elemIndex]; + //TODO (Optimization) It's assumed all elements are the same type, so this is constant. + //TODO Therefore there's no need to do this everytime. + unsigned numNodes = bulkData.num_nodes(elem); + stk::mesh::Entity const* nodes = bulkData.begin_nodes(elem); + for (unsigned inode = 0; inode < numNodes; ++inode) { + double *coord = stk::mesh::field_data(*coords, nodes[inode]); + int lid = globalMapG.LID((int)bulkData.identifier(nodes[inode]) -1); + elemToNode(elem_ct,inode) = lid; + if(lid != -1) { + nodeCoord(lid,0) = coord[0]; + nodeCoord(lid,1) = coord[1]; + if(spaceDim==3) + nodeCoord(lid,2) = coord[2]; + } + }//end node loop + + auto data = cellType.getCellTopologyData(); + for(unsigned iedge=0; iedgeedge[iedge].node[0]; + int n1 = data->edge[iedge].node[1]; + int lid0 = globalMapG.LID((int)bulkData.identifier(nodes[n0]) -1); + int lid1 = globalMapG.LID((int)bulkData.identifier(nodes[n1]) -1); + if(lid0 != -1 && lid1 != -1) { + int lo = std::min(lid0,lid1); + int hi = std::max(lid0,lid1); + std::pair key(lo,hi); + if (local_edge_hash.find(key) == local_edge_hash.end()) { + int new_edge_id = edge_vector.size(); + local_edge_hash[key] = new_edge_id; + edge_vector.push_back(key); + elemToEdge(elem_ct,iedge) = new_edge_id; + } + else { + elemToEdge(elem_ct,iedge) = local_edge_hash[key]; + } + } + }//end edge loop + + sigmaVal(elem_ct) = 1;// Not doing sigma here + + elem_ct++; + }//end element loop + }//end bucket loop + + Intrepid::FieldContainer edgeToNode(edge_vector.size(), 2); + for(int i=0; i<(int)edge_vector.size(); i++) { + edgeToNode(i,0) = edge_vector[i].first; + edgeToNode(i,1) = edge_vector[i].second; + } + + + MLStatistics.Phase1(elemToNode,elemToEdge,edgeToNode,nodeCoord,sigmaVal); + + } +#endif + /**********************************************************************************/ /******************** DEFINE WORKSETS AND LOOP OVER THEM **************************/ /**********************************************************************************/ @@ -873,6 +970,17 @@ int main(int argc, char *argv[]) { Time.ResetStartTime(); } + + + + /**********************************************************************************/ + /***************************** STATISTICS (Part II) ******************************/ + /**********************************************************************************/ +#ifdef HAVE_TRILINOSCOUPLINGS_MUELU + if(do_statistics) + MLStatistics.Phase2a(worksetJacobDet,worksetCubWeights); +#endif + /**********************************************************************************/ /* Assemble into Global Matrix */ /**********************************************************************************/ @@ -917,6 +1025,7 @@ int main(int argc, char *argv[]) { }// end workset loop + StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete(); rhsVector.GlobalAssemble(); @@ -926,6 +1035,26 @@ int main(int argc, char *argv[]) { << Time.ElapsedTime() << " seconds" << std::endl; Time.ResetStartTime(); } +/**********************************************************************************/ +/***************************** STATISTICS (Part IIb) ******************************/ +/**********************************************************************************/ +#ifdef HAVE_TRILINOSCOUPLINGS_MUELU + if(do_statistics){ + MLStatistics.Phase2b(Teuchos::rcp(&StiffMatrix.Graph(),false),Teuchos::rcp(&nCoord,false)); + } +#endif +/**********************************************************************************/ +/***************************** STATISTICS (Part III) ******************************/ +/**********************************************************************************/ +#ifdef HAVE_TRILINOSCOUPLINGS_MUELU + if(do_statistics){ + MLStatistics.Phase3(); + Teuchos::ParameterList problemStatistics = MLStatistics.GetStatistics(); + if(MyPID==0) std::cout<<"*** Problem Statistics ***"< const Scalar exactSolution(const Scalar& x, const Scalar& y, const Scalar& z) { + // Patch test: tri-linear function is in the FE space and should be recovered + return 1. + x + y + z + x*y + x*z + y*z + x*y*z; + + // Patch test - tet: function is in the FE space and should be recovered - return 1. + x + y + z ; + // return 1. + x + y + z ; // Patch test - hex: tri-linear function is in the FE space and should be recovered // return 1. + x + y + z + x*y + x*z + y*z + x*y*z; From 574ccf342ec37c5c8b5106cb5dad5e2779e26041 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Fri, 13 May 2022 15:08:22 -0600 Subject: [PATCH 2/7] TrilinosCouplings: Fixing bug in STK driver --- .../examples/scaling/example_Poisson_stk.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp index 9dced139059b..a4ab53ffb0d7 100644 --- a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp +++ b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp @@ -988,6 +988,7 @@ int main(int argc, char *argv[]) { //"WORKSET CELL" loop: local cell ordinal is relative to numLocalElems //JJH runs from 0 to (#local cells - 1) int worksetCellOrdinal = 0; + for (size_t bucketIndex = 0; bucketIndex < localElementBuckets.size(); ++bucketIndex) { stk::mesh::Bucket &elemBucket = *localElementBuckets[bucketIndex]; @@ -1018,8 +1019,9 @@ int main(int argc, char *argv[]) { }// end cell row loop + worksetCellOrdinal++; }// end workset cell loop - worksetCellOrdinal++; + } //for (size_t bucketIndex = 0; ... @@ -1087,6 +1089,8 @@ int main(int argc, char *argv[]) { ML_Epetra::Apply_BCsToMatrixRows(&(ownedBoundaryNodes[0]), ownedBoundaryNodes.size(), StiffMatrix); ML_Epetra::Remove_Zeroed_Rows(StiffMatrix,0.0); + + if(MyPID==0) { std::cout << "Adjust global matrix and rhs due to BCs " << Time.ElapsedTime() From 5ddce0776420ca0a4946be6c216bd6c54aa9bc6e Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Fri, 13 May 2022 15:54:08 -0600 Subject: [PATCH 3/7] TrilinosCouplings: Oops --- .../trilinoscouplings/examples/scaling/example_Poisson_stk.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp index a4ab53ffb0d7..6b68880968b5 100644 --- a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp +++ b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp @@ -171,8 +171,6 @@ #ifdef HAVE_TRILINOSCOUPLINGS_MUELU #include "TrilinosCouplings_Statistics.hpp" -#else -#error "BORK" #endif /*********************************************************/ From 350153cc39aac11c8bd3063bea528e2287900706 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Fri, 13 May 2022 15:55:11 -0600 Subject: [PATCH 4/7] TrilinosCouplings: Oops --- .../trilinoscouplings/examples/scaling/example_Poisson_stk.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp index 6b68880968b5..2b7b8412cd68 100644 --- a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp +++ b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp @@ -691,7 +691,7 @@ int main(int argc, char *argv[]) { // in all other cases EpetraExt::MultiVectorToMatrixMarketFile(coordsFilename.c_str(),nCoord,0,0,false); -#if 1 +#if 0 // Processor-by-processor output char fn[80]; sprintf(fn,"test-%d-%s",MyPID,coordsFilename.c_str()); From aa3d1da0e4ad91201062d5db644d6059f56fb146 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Mon, 16 May 2022 08:20:54 -0600 Subject: [PATCH 5/7] TrilinosCouplings: Fixing ifdef guards --- .../examples/scaling/example_Poisson_stk.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp index 2b7b8412cd68..09c8cf81e4fe 100644 --- a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp +++ b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp @@ -553,8 +553,10 @@ int main(int argc, char *argv[]) { std::cout<<"Cell Topology: "< MLStatistics(numGlobalElements); +#else + bool do_statistics = !strcmp(cellType.getName(),"Hexahedron_8"); std::cout<<"do_statistics = "< elemToNode(numLocalElems,numNodesPerElem); Intrepid::FieldContainer elemToEdge(numLocalElems,numEdgesPerElem); @@ -974,7 +976,7 @@ int main(int argc, char *argv[]) { /**********************************************************************************/ /***************************** STATISTICS (Part II) ******************************/ /**********************************************************************************/ -#ifdef HAVE_TRILINOSCOUPLINGS_MUELU +#if defined(HAVE_TRILINOSCOUPLINGS_MUELU) && defined(HAVE_MUELU_EPETRA) if(do_statistics) MLStatistics.Phase2a(worksetJacobDet,worksetCubWeights); #endif @@ -1038,7 +1040,7 @@ int main(int argc, char *argv[]) { /**********************************************************************************/ /***************************** STATISTICS (Part IIb) ******************************/ /**********************************************************************************/ -#ifdef HAVE_TRILINOSCOUPLINGS_MUELU +#if defined(HAVE_TRILINOSCOUPLINGS_MUELU) && defined(HAVE_MUELU_EPETRA) if(do_statistics){ MLStatistics.Phase2b(Teuchos::rcp(&StiffMatrix.Graph(),false),Teuchos::rcp(&nCoord,false)); } @@ -1046,7 +1048,7 @@ int main(int argc, char *argv[]) { /**********************************************************************************/ /***************************** STATISTICS (Part III) ******************************/ /**********************************************************************************/ -#ifdef HAVE_TRILINOSCOUPLINGS_MUELU +#if defined(HAVE_TRILINOSCOUPLINGS_MUELU) && defined(HAVE_MUELU_EPETRA) if(do_statistics){ MLStatistics.Phase3(); Teuchos::ParameterList problemStatistics = MLStatistics.GetStatistics(); From c39587a1cca97b87104a7db7ffc6a0e614c8b25f Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Mon, 16 May 2022 11:05:14 -0600 Subject: [PATCH 6/7] TrilinosCouplings: Fixing ifdef guards --- .../trilinoscouplings/examples/scaling/example_Poisson_stk.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp index 09c8cf81e4fe..aa47c0331ce4 100644 --- a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp +++ b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp @@ -555,8 +555,6 @@ int main(int argc, char *argv[]) { #if defined(HAVE_TRILINOSCOUPLINGS_MUELU) && defined(HAVE_MUELU_EPETRA) MachineLearningStatistics_Hex3D MLStatistics(numGlobalElements); -#else - bool do_statistics = !strcmp(cellType.getName(),"Hexahedron_8"); std::cout<<"do_statistics = "< Date: Mon, 16 May 2022 13:25:51 -0600 Subject: [PATCH 7/7] TrilinosCouplings: Fixing ifdef guards --- .../trilinoscouplings/examples/scaling/example_Poisson_stk.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp index aa47c0331ce4..20575a44a6fa 100644 --- a/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp +++ b/packages/trilinoscouplings/examples/scaling/example_Poisson_stk.cpp @@ -547,7 +547,9 @@ int main(int argc, char *argv[]) { } // end loop over mesh parts int numNodesPerElem = cellType.getNodeCount(); +#if defined(HAVE_TRILINOSCOUPLINGS_MUELU) && defined(HAVE_MUELU_EPETRA) int numEdgesPerElem = cellType.getEdgeCount(); +#endif if(MyPID==0) { std::cout<<"Cell Topology: "<