74 const char* sOutputFileName,
75 const char* sInputFileName,
76 bool bAbortAfterInput)
77 :
DataManager(HP, OF, pS, dInitialTime, sOutputFileName, sInputFileName, bAbortAfterInput)
79 silent_cerr(
"fatal error: you are building SchurDataManager, "
80 "but mbdyn was compiled without MPI. "
81 "Something weird is happening. "
82 "Anyway, please compile with -DUSE_MPI "
83 "to enable parallel solution" << std::endl);
170 const char* sOutputFileName,
171 const char* sInputFileName,
172 bool bAbortAfterInput)
173 :
DataManager(HP, OF, pS, dInitialTime, sOutputFileName, sInputFileName, bAbortAfterInput),
190 wgtflag(WEIGHT_VERTICES),
192 Partitioner(PARTITIONER_DEFAULT),
196 iTotalExpConnections(0)
198 DEBUGCOUT(
"Entering SchurDataManager" << std::endl);
201 DataComm = MBDynComm.Dup();
203 MyRank = DataComm.Get_rank();
211 const char* sKeyWords[] = {
215 "number" "of" "connections",
225 "aerodynamic" "element",
286 pedantic_cerr(
"no explicit connections declared "
287 "for this input file" << std::endl);
292 pedantic_cerr(
"no explicit connections declared "
293 "for this input file" << std::endl);
300 silent_cerr(
"Error: \"begin: parallel;\" expected at line "
309 silent_cerr(
"Error: Weight flag expected "
324 }
else if (HP.
IsKeyWord(
"communication")) {
331 silent_cerr(
"invalid weight "
343 silent_cerr(
"the partition "
344 "assignment is not complete; "
345 "only " << i <<
" vertices "
355 silent_cerr(
"illegal value "
358 "assignment[" << i <<
"] "
360 <<
"; must be between 0 and "
372 silent_cerr(
"METIS partitioner not available; aborting..." << std::endl);
380 silent_cerr(
"CHACO partitioner not available; aborting..." << std::endl);
387 silent_cerr(
"unknown partitioner "
396 silent_cerr(
"Error: \"end: parallel;\" expected "
398 <<
"; aborting..." << std::endl);
404 silent_cerr(
"Unknown input at line "
409 case NUMBEROFCONNECTIONS:
420 silent_cerr(
"Error: <Connection> "
428 for (
int k = 0; k < 2; k++) {
460 case AERODYNAMICELEMENT:
497 silent_cerr(
"Error: invalid element type "
499 <<
"; aborting..." << std::endl);
507 silent_cerr(
"Error: label expected "
509 <<
" element type at line "
511 <<
"; aborting..." << std::endl);
518 silent_cerr(
"Error at line "
521 <<
"(" << j <<
") undefined; "
522 "aborting..." << std::endl);
555 silent_cerr(
"Error: invalid node type "
557 <<
"; aborting..." << std::endl);
565 silent_cerr(
"Error: label expected "
567 <<
"; aborting..." << std::endl);
574 silent_cerr(
"Error: at line "
577 <<
"(" << j <<
") undefined; "
578 <<
"aborting..." << std::endl);
585 silent_cerr(
"Unknown input at line "
587 <<
"; aborting..." << std::endl);
593 ASSERT(iNumNodes == iNumElems);
594 if (iNumNodes + iNumElems != 2*iTotalExpConnections) {
595 silent_cerr(
"Error: total number of nodes and elements"
596 " in the parallel section at line "
598 <<
" is not consistent; aborting..." << std::endl);
663 silent_cerr(
"SchurDataManager::HowManyDofs: "
664 "illegal request (" <<
unsigned(who) <<
")"
684 silent_cerr(
"SchurDataManager::GetDofsList: "
685 "illegal request (" <<
unsigned(who) <<
")"
701 DEBUGCOUT(
"Entering SchurDataManager::CreatePartition()" << std::endl);
703 std::vector<std::vector<int> > adj(iTotVertices);
704 std::vector<int> CommWgts(iTotVertices);
707 std::vector<int> VertexWgts(iTotVertices);
720 memset(Vertices.
pXadj, 0, (iTotVertices + 1)*
sizeof(
int));
726 int iMaxConnectionsPerVertex =
729 int GravityPos = 0, AirPropPos = 0;
731 unsigned int* pRotLab = NULL;
736 memset(pRotPos, 0, iNum*
sizeof(
int));
737 memset(pRotLab, 0, iNum*
sizeof(
unsigned int));
742 std::vector<const Node *> connectedNodes;
748 InitList(Vertices.
pAdjncy, iTotVertices*iMaxConnectionsPerVertex, ADJ_UNDEFINED);
751 for (
unsigned int i = 0; i <
iTotNodes; i++) {
756 Node** ppCurrNode = NULL;
763 for (ElemVecType::const_iterator pTmpEl =
Elems.begin();
764 pTmpEl !=
Elems.end();
775 pRotLab[iNumRt] = (*pTmpEl)->GetLabel();
779 (*pTmpEl)->GetConnectedNodes(connectedNodes);
783 (*pTmpEl)->WorkSpaceDim(&dimA, &dimB);
784 VertexWgts[iCount] = dimA * dimB;
786 CommWgts[iCount] = (*pTmpEl)->iGetNumDof()*(*pTmpEl)->iGetNumDof();
788 for (std::vector<const Node *>::const_iterator i = connectedNodes.begin();
789 i != connectedNodes.end();
792 Vertices.
pXadj[iCount + 1]++;
796 unsigned label = (*i)->GetLabel();
799 position = ppCurrNode - ppNodes;
805 VertexWgts[iCount] += (*ppCurrNode)->
iGetNumDof();
809 if ((iCount*iMaxConnectionsPerVertex) + Vertices.
pXadj[iCount + 1] - 1 < iTotVertices*iMaxConnectionsPerVertex) {
810 Vertices.
pAdjncy[(iCount*iMaxConnectionsPerVertex)
811 + Vertices.
pXadj[iCount + 1] - 1] = position;
815 Vertices.
pXadj[position + 1]++;
816 if ((position*iMaxConnectionsPerVertex) + Vertices.
pXadj[position + 1] - 1 < iTotVertices*iMaxConnectionsPerVertex) {
817 Vertices.
pAdjncy[(position*iMaxConnectionsPerVertex)
818 + Vertices.
pXadj[position + 1] - 1] = iCount;
822 CommWgts[
position] = (*ppCurrNode)->iGetNumDof();
826 if (iTotalExpConnections != 0) {
841 Vertices.
pXadj[iNdPos + 1]++;
842 Vertices.
pXadj[iTotNodes + iElPos + 1]++;
843 Vertices.
pAdjncy[(iNdPos*iMaxConnectionsPerVertex)
844 + Vertices.
pXadj[iNdPos + 1] - 1] = iElPos + iTotNodes;
845 Vertices.
pAdjncy[((iTotNodes + iElPos)*iMaxConnectionsPerVertex)
846 + Vertices.
pXadj[iTotNodes + iElPos + 1] - 1] = iNdPos;
852 if (Vertices.
pXadj[i] > iMaxConnectionsPerVertex) {
853 iMax = Vertices.
pXadj[i];
857 if (iMax <= iMaxConnectionsPerVertex) {
861 iMaxConnectionsPerVertex = iMax;
871 Pack(Vertices.
pAdjncy, iTotVertices*iMaxConnectionsPerVertex);
880 std::ofstream ofPartition;
882 ofPartition.open(
"partition.debug");
883 ofPartition <<
"# METIS-like Input File" << std::endl
884 <<
"# Column 1: computational weight" << std::endl
885 <<
"# Column 2: communication weight" << std::endl
886 <<
"# Total Vertices: " << iTotVertices << std::endl
887 <<
"# Nodes" << std::endl;
888 for (
unsigned int i = 0; i <
iTotNodes; i++) {
889 ofPartition <<
"# " << i <<
" Node Type: "
890 <<
"(" <<
psNodeNames[ppNodes[i]->GetNodeType()] <<
")"
891 <<
" Label: " << ppNodes[i]->GetLabel() << std::endl
892 << VertexWgts[i] <<
" " << CommWgts[i];
893 for (
int j = Vertices.
pXadj[i]; j < Vertices.
pXadj[i + 1]; j++) {
894 ofPartition <<
" " << Vertices.
pAdjncy[j];
896 ofPartition << std::endl;
898 ofPartition <<
"# Elements" << std::endl;
899 for (
unsigned int i = 0; i < iTotElem; i++) {
900 ofPartition <<
"# " << i + iTotNodes <<
" Element Type: "
902 <<
" Label: " <<
Elems[i]->GetLabel() << std::endl
905 for (
int j = Vertices.
pXadj[i + iTotNodes];
906 j < Vertices.
pXadj[i + iTotNodes + 1];
908 ofPartition <<
" " << Vertices.
pAdjncy[j];
910 ofPartition << std::endl;
927 vwgt = &VertexWgts[0];
932 silent_cout(
"communication weights currently unsupported by CHACO; trying to emulate..." << std::endl);
938 silent_cout(
"edges weights currently unsupported by MBDyn" << std::endl);
956 if (
wgtflag & WEIGHT_VERTICES) {
957 vwgt = &VertexWgts[0];
962 silent_cout(
"edges weights currently unsupported by MBDyn" << std::endl);
981 silent_cerr(
"no partition library is available; "
982 "partition assignments must be provided "
983 "manually" << std::endl);
998 InterfNodes.
pXadj = NULL;
1002 memset(InterfNodes.
pXadj, 0, DataCommSize*
sizeof(
int));
1008 memset(InterfNodes.
pAdjncy, 0, DataCommSize*iMaxInterfNodes*8*
sizeof(
int));
1009 InitList(InterfNodes.
pAdjncy, DataCommSize*iMaxInterfNodes*8, ADJ_UNDEFINED);
1017 for (
unsigned int i = 0; i <
iTotNodes; i++) {
1023 for (
int j = Vertices.
pXadj[i]; j < Vertices.
pXadj[i + 1]; j++) {
1026 InterfNodes.
pXadj[TmpPrc] += 1;
1027 int k = TmpPrc * iMaxInterfNodes * 2 + InterfNodes.
pXadj[TmpPrc] - 1;
1028 if (k < DataCommSize * iMaxInterfNodes * 2) {
1046 if (InterfNodes.
pXadj[i] > iMaxInterfNodes) {
1047 iMax = InterfNodes.
pXadj[i];
1051 DataComm.Allreduce(&iMax, &iRMax, 1, MPI::INT,
MPI::MAX);
1053 if (iMax <= iMaxInterfNodes) {
1057 iMaxInterfNodes = iMax;
1063 MPI::Request* pRReq = NULL;
1064 MPI::Request* pSReq = NULL;
1069 pRReq[i] = MPI::REQUEST_NULL;
1070 pSReq[i] = MPI::REQUEST_NULL;
1073 const int DIM_TAG = 10;
1080 pRReq[i] = DataComm.Irecv(InterfNodes.
pAdjncy + iMaxInterfNodes + i*iMaxInterfNodes*2,
1081 iMaxInterfNodes, MPI::INT, i, DIM_TAG);
1082 pSReq[i] = DataComm.Isend(InterfNodes.
pAdjncy + i*iMaxInterfNodes*2,
1083 iMaxInterfNodes, MPI::INT, i, DIM_TAG);
1127 unsigned int pTmpLab = pIV->
GetLabel();
1129 for (
int k = 0; k < iNumRt; k++) {
1130 if (pTmpLab == pRotLab[k]) {
1135 for (
int j = 0; j < iMyTotRot; j++) {
1136 if (pos != pMyRot[j]) {
1137 pMyRot[iMyTotRot] = pos;
1143 pMyRot[iMyTotRot] = pos;
1153 for (
int i = 0; i < iNumRt; i++) {
1154 if (iMyTotRot == 0) {
1155 color = MPI::UNDEFINED;
1160 silent_cout(
"InducedVelocity(" <<
Elems[pRotPos[i]]->GetLabel() <<
")"
1161 <<
" assigned to process "
1169 pIndVelComm[i] = MBDynComm.Split(color, key);
1171 IndVelComm[i] = MPI::COMM_WORLD.Split(color, key);
1175 r->InitializeIndVelComm(pIndVelComm + i);
1178 for (
int i = 0; i < iMyTotRot; i++) {
1190 for (
unsigned int i = 0; i < iTotVertices -
iTotNodes; i++) {
1202 MPI::Request::Waitall(DataCommSize, pRReq);
1203 MPI::Request::Waitall(DataCommSize, pSReq);
1206 bool bRecvFlag =
false, bSentFlag =
false;
1209 bRecvFlag = MPI::Request::Testall(DataCommSize, pRReq);
1213 bSentFlag = MPI::Request::Testall(DataCommSize, pSReq);
1216 if (bRecvFlag && bSentFlag) {
1225 std::sort(InterfNodes.
pAdjncy,
1226 InterfNodes.
pAdjncy + iMaxInterfNodes * 2 * DataCommSize);
1229 int* p = std::unique(InterfNodes.
pAdjncy,
1230 InterfNodes.
pAdjncy + iMaxInterfNodes * 2 * DataCommSize);
1240 unsigned int* llabels = NULL;
1254 int* pPosIntElems = NULL;
1268 if (iIVIsMine == 1) {
1270 for (std::vector<const Node *>::const_iterator j = connectedNodes.begin();
1271 j != connectedNodes.end();
1274 unsigned int* p = std::find(llabels, llabels + iNumIntNodes, (*j)->GetLabel());
1275 if (p != llabels + iNumIntNodes) {
1295 for (std::vector<const Node *>::const_iterator j = connectedNodes.begin();
1296 j != connectedNodes.end();
1299 unsigned int* p = std::find(llabels, llabels + iNumIntNodes, (*j)->GetLabel());
1300 if (p != llabels + iNumIntNodes) {
1338 if (TmpDofNum != 0) {
1339 if (i != pPosIntElems[i2Count]) {
1341 integer First = (pWithDofs)->iGetFirstIndex();
1345 for (
int j = 1; j < TmpDofNum; j++) {
1357 for (
int i = 1; i < iNumIntNodes + 1; i++) {
1367 if (ppNodes[InterfNodes.
pAdjncy[i]]->iGetNumDof() != 0) {
1368 integer First = ppNodes[InterfNodes.
pAdjncy[i]]->iGetFirstIndex();
1377 for (
unsigned int j = 1;
1378 j < ppNodes[InterfNodes.
pAdjncy[i]]->iGetNumDof();
1396 integer First = (pWithDofs)->iGetFirstIndex();
1399 for (
int j = 1; j < TmpDofNum; j++) {
1407 if ( Vertices.
pXadj != NULL) {
1411 if ( Vertices.
pAdjncy != NULL) {
1415 if ( InterfNodes.
pXadj != NULL) {
1419 if ( InterfNodes.
pAdjncy != NULL) {
1431 if ( pPosIntElems != NULL) {
1435 if ( lTypes != NULL) {
1439 if ( llabels != NULL) {
1449 silent_cout(
"Making partition file...");
1454 time_t tCurrTime(time(NULL));
1456 <<
"# Partition file for MBDyn. Time: " << ctime(&tCurrTime) << std::endl
1457 <<
"# Partition produced ";
1484 <<
"# Control data to verify the partition "
1487 <<
"Total number of processes: " << DataCommSize <<
";" << std::endl
1488 <<
"Process #: " <<
MyRank <<
";" << std::endl
1489 <<
"Total number of Nodes: " << iTotNodes <<
";" << std::endl
1490 <<
"Total number of Elements: " <<
Elems.size() <<
";" << std::endl
1491 <<
"Total number of Dofs: " <<
iTotDofs <<
";" << std::endl;
1496 <<
"Local Dofs number: " <<
iNumLocDofs << std::endl;
1500 <<
"Local Nodes: " << iNumLocNodes << std::endl;
1514 <<
"Local Elements: "<< iNumLocElems << std::endl;
1528 <<
"Interface Dofs number: " <<
iNumIntDofs << std::endl
1530 <<
"Interface Nodes: " << iNumIntNodes << std::endl;
1544 <<
"Elements whose internal dofs are interface dofs: "
1545 << iNumIntElems << std::endl;
1559 <<
"Local Dofs List:" << std::endl;
1563 if (!(iNumLocDofs%10)) {
1567 if (iNumLocDofs%10) {
1575 <<
"Local Interface Dofs List:" <<std::endl;
1585 silent_cout(
"... partition file done" << std::endl);
1595 for (; pOld < pList + dim; pOld++) {
1596 if (*pOld != ADJ_UNDEFINED) {
1607 for (
int i = 0; i < dim; i++) {
1615 for (
int i = 0; i < dim; i++) {
1625 unsigned int* pbegin =
pLabelsList + (ppFirst - ppNodes);
1626 unsigned int* p = std::find(pbegin, pbegin + dim, label);
1635 DEBUGCOUT(
"Entering SchurDataManager::AssRes()" << std::endl);
1636 ASSERT(pWorkVec != NULL);
1641 }
catch (ChangedEquationStructure) {
1644 if (MyElemIter.bGetCurr(pEl) ==
true) {
1645 silent_cerr(
"Jacobian reassembly requested by "
1648 "currently unsupported by Schur data manager."
1652 silent_cerr(
"Jacobian reassembly requested "
1653 "by an element; currently unsupported "
1654 "by Schur data manager." << std::endl);
1664 DEBUGCOUT(
"Entering SchurDataManager::AssJac()" << std::endl);
1674 DEBUGCOUT(
"Entering SchurDataManager::Update()" << std::endl);
1699 DEBUGCOUT(
"Entering SchurDataManager::DerivativesUpdate()" << std::endl);
1734 DEBUGCOUT(
"Entering SchurDataManager::BeforePredict()" << std::endl);
1759 DEBUGCOUT(
"Entering SchurDataManager::AfterPredict()" << std::endl);
1788 DEBUGCOUT(
"Entering SchurDataManager::AfterConvergence()" << std::endl);
1840 DEBUGCOUT(
"Entering SchurDataManager::Output()" << std::endl);
integer * GetDofsList(DofType who) const
enum SchurDataManager::PartitionLibrary Partitioner
ElemContainerType ElemContainer
#define MBDYN_EXCEPT_ARGS
#define DEBUGCOUTFNAME(fname)
void DerivativesUpdate(void) const
virtual integer GetInt(integer iDefval=0)
doublereal dGetTime(void) const
const integer iDefaultMaxInterfNodes
static int position(const MathParser::MathArgs &args)
#define SAFEDELETEARR(pnt)
void BeforePredict(VectorHandler &X, VectorHandler &XP, VectorHandler &XPrev, VectorHandler &XPPrev) const
int mbdyn_METIS_PartGraph(int iTotVertices, int *pXadj, int *pAdjncy, int *pVertexWgts, int *pCommWgts, int *pEdgeWgts, int DataCommSize, int *pParAmgProcs)
virtual void AssRes(VectorHandler &ResHdl, doublereal dCoef)
Node ** SearchNode(Node **ppFirst, int dim, unsigned int &label)
virtual Elem::Type GetElemType(void) const =0
VecIter< Elem * > MyElemIter
void AssRes(VectorHandler &ResHdl, doublereal dCoef)
void AfterPredict(void) const
std::ostream & Partition(void) const
void AssJac(MatrixHandler &JacHdl, doublereal dCoef)
void Output(bool force=false) const
virtual const InducedVelocity * pGetInducedVelocity(void) const
DriveCaller * pOutputMeter
unsigned int * pLabelsList
virtual void MakeRestart(void)
virtual bool IsKeyWord(const char *sKeyWord)
doublereal dLastRestartTime
struct DataManager::NodeDataStructure NodeData[Node::LASTNODETYPE]
virtual void AssJac(MatrixHandler &JacHdl, doublereal dCoef)
VectorHandler * pXPrimeCurr
integer iRestartIterations
#define ASSERT(expression)
Elem ** ppFindElem(Elem::Type Typ, unsigned int uElem) const
virtual void BeforePredict(VectorHandler &, VectorHandler &, VectorHandler &, VectorHandler &) const
virtual unsigned int iGetNumDof(void) const
void OutputPartition(void)
virtual doublereal dGet(const doublereal &dVar) const =0
virtual void GetConnectedNodes(std::vector< const Node * > &connectedNodes) const
virtual unsigned int iGetNumDof(void) const =0
void Init(const T *p, unsigned i)
const char * psNodeNames[]
Dof * pGetDofsList(void) const
virtual Node::Type GetNodeType(void) const =0
#define SAFENEWARRNOFILL(pnt, item, sz)
virtual void Output(OutputHandler &OH) const
integer HowManyDofs(DofType who) const
const char * psElemNames[]
SchurDataManager(MBDynParser &HP, unsigned OF, Solver *pS, doublereal dInitialTime, const char *sOutputFileName, const char *sInputFileName, bool bAbortAfterInput)
virtual int GetWord(void)
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
void Pack(int *pList, int dim)
void chaco_interface(const int iTotVertices, int *start, int *adjacency, int *vertex_weights, int *comm_weights, int *edge_weights, const int num_processors, int *pParAmgProcs)
#define SAFENEWARR(pnt, item, sz)
virtual void AfterPredict(VectorHandler &X, VectorHandler &XP)
virtual HighParser::ErrOut GetLineData(void) const
void AfterConvergence(void) const
unsigned int GetLabel(void) const
struct DataManager::ElemDataStructure ElemData[Elem::LASTELEMTYPE]
Node * pFindNode(Node::Type Typ, unsigned int uNode) const
void CreatePartition(void)
const integer iDefaultMaxNodesPerElem
const integer iDefaultMaxConnectionsPerVertex
VariableSubMatrixHandler * pWorkMat
void InitList(int *list, int dim, int value)
virtual void Update(const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)