ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/devel_omp/src/parallel/ForceMatrixDecomposition.cpp
(Generate patch)

Comparing branches/devel_omp/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1598 by mciznick, Wed Jul 27 14:26:53 2011 UTC vs.
Revision 1608 by mciznick, Tue Aug 9 01:58:56 2011 UTC

# Line 50 | Line 50 | ForceMatrixDecomposition::ForceMatrixDecomposition(Sim
50  
51   ForceMatrixDecomposition::ForceMatrixDecomposition(SimInfo* info, InteractionManager* iMan) :
52          ForceDecomposition(info, iMan) {
53 –
53          // In a parallel computation, row and colum scans must visit all
54          // surrounding cells (not just the 14 upper triangular blocks that
55          // are used when the processor can see all pairs)
# Line 323 | Line 322 | void ForceMatrixDecomposition::distributeInitialData()
322                                  }
323                          }
324                  }
325 +        }
326 +
327 +        Globals* simParams_ = info_->getSimParams();
328 +        if (simParams_->haveNeighborListReorderFreq())
329 +        {
330 +                neighborListReorderFreq = simParams_->getNeighborListReorderFreq();
331 +        } else
332 +        {
333 +                neighborListReorderFreq = 0;
334          }
335 +        reorderFreqCounter = 1;
336  
337          createGtypeCutoffMap();
338  
# Line 796 | Line 805 | void ForceMatrixDecomposition::collectData() {
805          pairwisePot += pot_temp[ii];
806   #endif
807  
808 <        cerr << "pairwisePot = " << pairwisePot << "\n";
808 >        //      cerr << "pairwisePot = " << pairwisePot << "\n";
809   }
810  
811   int ForceMatrixDecomposition::getNAtomsInRow() {
# Line 847 | Line 856 | Vector3d ForceMatrixDecomposition::getIntergroupVector
856          Vector3d d;
857  
858          d = snap_->cgData.position[cg2->getLocalIndex()] - snap_->cgData.position[cg1->getLocalIndex()];
859 < /*      cerr << "cg1_gid = " << cg1->getGlobalIndex() << "\tcg1_lid = " << cg1->getLocalIndex() << "\tcg1p = "
860 <                        << snap_->cgData.position[cg1->getLocalIndex()] << "\n";
861 <        cerr << "cg2_gid = " << cg2->getGlobalIndex() << "\tcg2_lid = " << cg2->getLocalIndex() << "\tcg2p = "
862 <                        << snap_->cgData.position[cg2->getLocalIndex()] << "\n";*/
859 >        /*      cerr << "cg1_gid = " << cg1->getGlobalIndex() << "\tcg1_lid = " << cg1->getLocalIndex() << "\tcg1p = "
860 >         << snap_->cgData.position[cg1->getLocalIndex()] << "\n";
861 >         cerr << "cg2_gid = " << cg2->getGlobalIndex() << "\tcg2_lid = " << cg2->getLocalIndex() << "\tcg2p = "
862 >         << snap_->cgData.position[cg2->getLocalIndex()] << "\n";*/
863  
864          snap_->wrapVector(d);
865          return d;
# Line 924 | Line 933 | bool ForceMatrixDecomposition::skipAtomPair(int atom1,
933   bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
934          int unique_id_1, unique_id_2;
935  
936 < //      cerr << "sap with atom1, atom2 =\t" << atom1 << "\t" << atom2 << "\n";
936 >        //      cerr << "sap with atom1, atom2 =\t" << atom1 << "\t" << atom2 << "\n";
937   #ifdef IS_MPI
938          // in MPI, we have to look up the unique IDs for each atom
939          unique_id_1 = AtomRowToGlobal[atom1];
# Line 980 | Line 989 | void ForceMatrixDecomposition::addForceToAtomRow(int a
989   #else
990          snap_->atomData.force[atom1] += fg;
991   #endif
992 + }
993 +
994 + void ForceMatrixDecomposition::addForceToAtomRowOMP(int atom1, Vector3d fg) {
995 +        #pragma omp critical
996 +        {
997 +                snap_->atomData.force[atom1] += fg;
998 +        }
999   }
1000  
1001   void ForceMatrixDecomposition::addForceToAtomColumn(int atom2, Vector3d fg) {
# Line 990 | Line 1006 | void ForceMatrixDecomposition::addForceToAtomColumn(in
1006   #endif
1007   }
1008  
1009 + void ForceMatrixDecomposition::addForceToAtomColumnOMP(int atom2, Vector3d fg) {
1010 +        #pragma omp critical
1011 +        {
1012 +                snap_->atomData.force[atom2] += fg;
1013 +        }
1014 + }
1015 +
1016   // filling interaction blocks with pointers
1017   void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat, int atom1, int atom2) {
1018 +
1019 +        idat.excluded = excludeAtomPair(atom1, atom2);
1020 +
1021 + #ifdef IS_MPI
1022 +        idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]);
1023 +        //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
1024 +        //                         ff_->getAtomType(identsCol[atom2]) );
1025 +
1026 +        if (storageLayout_ & DataStorage::dslAmat)
1027 +        {
1028 +                idat.A1 = &(atomRowData.aMat[atom1]);
1029 +                idat.A2 = &(atomColData.aMat[atom2]);
1030 +        }
1031 +
1032 +        if (storageLayout_ & DataStorage::dslElectroFrame)
1033 +        {
1034 +                idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
1035 +                idat.eFrame2 = &(atomColData.electroFrame[atom2]);
1036 +        }
1037 +
1038 +        if (storageLayout_ & DataStorage::dslTorque)
1039 +        {
1040 +                idat.t1 = &(atomRowData.torque[atom1]);
1041 +                idat.t2 = &(atomColData.torque[atom2]);
1042 +        }
1043 +
1044 +        if (storageLayout_ & DataStorage::dslDensity)
1045 +        {
1046 +                idat.rho1 = &(atomRowData.density[atom1]);
1047 +                idat.rho2 = &(atomColData.density[atom2]);
1048 +        }
1049 +
1050 +        if (storageLayout_ & DataStorage::dslFunctional)
1051 +        {
1052 +                idat.frho1 = &(atomRowData.functional[atom1]);
1053 +                idat.frho2 = &(atomColData.functional[atom2]);
1054 +        }
1055 +
1056 +        if (storageLayout_ & DataStorage::dslFunctionalDerivative)
1057 +        {
1058 +                idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
1059 +                idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
1060 +        }
1061 +
1062 +        if (storageLayout_ & DataStorage::dslParticlePot)
1063 +        {
1064 +                idat.particlePot1 = &(atomRowData.particlePot[atom1]);
1065 +                idat.particlePot2 = &(atomColData.particlePot[atom2]);
1066 +        }
1067 +
1068 +        if (storageLayout_ & DataStorage::dslSkippedCharge)
1069 +        {
1070 +                idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]);
1071 +                idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]);
1072 +        }
1073 +
1074 + #else
1075 +
1076 +        idat.atypes = make_pair(atypesLocal[atom1], atypesLocal[atom2]);
1077 +        //idat.atypes = make_pair( ff_->getAtomType(idents[atom1]),
1078 +        //                         ff_->getAtomType(idents[atom2]) );
1079 +
1080 +        if (storageLayout_ & DataStorage::dslAmat)
1081 +        {
1082 +                idat.A1 = &(snap_->atomData.aMat[atom1]);
1083 +                idat.A2 = &(snap_->atomData.aMat[atom2]);
1084 +        }
1085 +
1086 +        if (storageLayout_ & DataStorage::dslElectroFrame)
1087 +        {
1088 +                idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
1089 +                idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
1090 +        }
1091 +
1092 +        if (storageLayout_ & DataStorage::dslTorque)
1093 +        {
1094 +                idat.t1 = &(snap_->atomData.torque[atom1]);
1095 +                idat.t2 = &(snap_->atomData.torque[atom2]);
1096 +        }
1097 +
1098 +        if (storageLayout_ & DataStorage::dslDensity)
1099 +        {
1100 +                idat.rho1 = &(snap_->atomData.density[atom1]);
1101 +                idat.rho2 = &(snap_->atomData.density[atom2]);
1102 +        }
1103  
1104 +        if (storageLayout_ & DataStorage::dslFunctional)
1105 +        {
1106 +                idat.frho1 = &(snap_->atomData.functional[atom1]);
1107 +                idat.frho2 = &(snap_->atomData.functional[atom2]);
1108 +        }
1109 +
1110 +        if (storageLayout_ & DataStorage::dslFunctionalDerivative)
1111 +        {
1112 +                idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
1113 +                idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
1114 +        }
1115 +
1116 +        if (storageLayout_ & DataStorage::dslParticlePot)
1117 +        {
1118 +                idat.particlePot1 = &(snap_->atomData.particlePot[atom1]);
1119 +                idat.particlePot2 = &(snap_->atomData.particlePot[atom2]);
1120 +        }
1121 +
1122 +        if (storageLayout_ & DataStorage::dslSkippedCharge)
1123 +        {
1124 +                idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]);
1125 +                idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]);
1126 +        }
1127 + #endif
1128 + }
1129 +
1130 + // filling interaction blocks with pointers
1131 + void ForceMatrixDecomposition::fillInteractionDataOMP(InteractionDataPrv &idat, int atom1, int atom2) {
1132 +
1133          idat.excluded = excludeAtomPair(atom1, atom2);
1134  
1135   #ifdef IS_MPI
# Line 1120 | Line 1257 | void ForceMatrixDecomposition::unpackInteractionData(I
1257  
1258   }
1259  
1260 + void ForceMatrixDecomposition::unpackInteractionDataOMP(InteractionDataPrv &idat, int atom1, int atom2) {
1261 + #pragma omp critical
1262 +        {
1263 +                pairwisePot += idat.pot;
1264 +
1265 +                snap_->atomData.force[atom1] += idat.f1;
1266 +                snap_->atomData.force[atom2] -= idat.f1;
1267 +        }
1268 + }
1269 +
1270   void ForceMatrixDecomposition::reorderGroupCutoffs(vector<int> &order) {
1271          vector<int> tmp = vector<int> (groupToGtype.size());
1272  
# Line 1144 | Line 1291 | void ForceMatrixDecomposition::reorderPosition(vector<
1291                  tmp[i] = snap_->cgData.position[i];
1292          }
1293  
1294 <        vector<int> mapPos = vector<int>(nGroups_);
1294 >        vector<int> mapPos = vector<int> (nGroups_);
1295          for (int i = 0; i < nGroups_; ++i)
1296          {
1297                  snap_->cgData.position[i] = tmp[order[i]];
# Line 1164 | Line 1311 | void ForceMatrixDecomposition::reorderPosition(vector<
1311                  }
1312          }
1313  
1314 < /*      if (info_->getNCutoffGroups() > 0)
1315 <        {
1316 <                for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1317 <                {
1318 <                        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
1319 <                        {
1320 <                                printf("gbI:%d locI:%d x:%f y:%f z:%f\n", cg->getGlobalIndex(), cg->getLocalIndex(),
1321 <                                                cgConfig->position[cg->getLocalIndex()].x(), cgConfig->position[cg->getLocalIndex()].y(),
1322 <                                                cgConfig->position[cg->getLocalIndex()].z());
1323 <                        }
1324 <                }
1325 <        } else
1326 <        {
1327 <                // center of mass of the group is the same as position of the atom
1328 <                // if cutoff group does not exist
1329 <                printf("ERROR!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1330 <                //                      cgConfig->position = config->position;
1331 <        }*/
1314 >        /*      if (info_->getNCutoffGroups() > 0)
1315 >         {
1316 >         for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1317 >         {
1318 >         for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
1319 >         {
1320 >         printf("gbI:%d locI:%d x:%f y:%f z:%f\n", cg->getGlobalIndex(), cg->getLocalIndex(),
1321 >         cgConfig->position[cg->getLocalIndex()].x(), cgConfig->position[cg->getLocalIndex()].y(),
1322 >         cgConfig->position[cg->getLocalIndex()].z());
1323 >         }
1324 >         }
1325 >         } else
1326 >         {
1327 >         // center of mass of the group is the same as position of the atom
1328 >         // if cutoff group does not exist
1329 >         printf("ERROR!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1330 >         //                     cgConfig->position = config->position;
1331 >         }*/
1332   }
1333  
1334   void ForceMatrixDecomposition::reorderGroupList(vector<int> &order) {
# Line 1200 | Line 1347 | void ForceMatrixDecomposition::reorderMemory(vector<ve
1347  
1348   void ForceMatrixDecomposition::reorderMemory(vector<vector<CutoffGroup *> > &H_c_l) {
1349          int n = 0;
1350 < //      printf("Reorder memory time:%f!!!!!!!!!!!!!!!!!!!!!!!!!!!\n",
1351 < //                      info_->getSnapshotManager()->getCurrentSnapshot()->getTime());
1350 >        //      printf("Reorder memory time:%f!!!!!!!!!!!!!!!!!!!!!!!!!!!\n",
1351 >        //                      info_->getSnapshotManager()->getCurrentSnapshot()->getTime());
1352  
1353          /* record the reordered atom indices */
1354          vector<int> k = vector<int> (nGroups_);
# Line 1222 | Line 1369 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1369   }
1370  
1371   vector<vector<CutoffGroup *> > ForceMatrixDecomposition::buildLayerBasedNeighborList() {
1372 < //      printf("buildLayerBasedNeighborList; nGroups:%d\n", nGroups_);
1372 >        //      printf("buildLayerBasedNeighborList; nGroups:%d\n", nGroups_);
1373          // Na = nGroups_
1374          /* cell occupancy counter */
1375 < //      vector<int> k_c;
1375 >        //      vector<int> k_c;
1376          /* c_i - has cell containing atom i (size Na) */
1377          vector<int> c = vector<int> (nGroups_);
1378          /* l_i - layer containing atom i (size Na) */
1379 < //      vector<int> l;
1379 >        //      vector<int> l;
1380  
1381          RealType rList_ = (largestRcut_ + skinThickness_);
1382          Snapshot* snap_ = sman_->getCurrentSnapshot();
# Line 1248 | Line 1395 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1395          int cellIndex;
1396          int nCtot = nCells_.x() * nCells_.y() * nCells_.z();
1397  
1398 < //      k_c = vector<int> (nCtot, 0);
1398 >        //      k_c = vector<int> (nCtot, 0);
1399  
1400          SimInfo::MoleculeIterator mi;
1401          Molecule* mol;
# Line 1277 | Line 1424 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1424                          whichCell.y() = nCells_.y() * scaled.y();
1425                          whichCell.z() = nCells_.z() * scaled.z();
1426  
1427 < //                      printf("pos x:%f y:%f z:%f cell x:%d y:%d z:%d\n", rs.x(), rs.y(), rs.z(), whichCell.x(), whichCell.y(),
1428 < //                                      whichCell.z());
1427 >                        //                      printf("pos x:%f y:%f z:%f cell x:%d y:%d z:%d\n", rs.x(), rs.y(), rs.z(), whichCell.x(), whichCell.y(),
1428 >                        //                                      whichCell.z());
1429  
1430                          // find single index of this cell:
1431                          cellIndex = Vlinear(whichCell, nCells_);
# Line 1287 | Line 1434 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1434                  }
1435          }
1436  
1437 < //      int k_c_curr;
1438 < //      int k_c_max = 0;
1437 >        //      int k_c_curr;
1438 >        //      int k_c_max = 0;
1439          /* the cell-layer occupancy matrix */
1440          vector<vector<CutoffGroup *> > H_c_l = vector<vector<CutoffGroup *> > (nCtot);
1441  
# Line 1305 | Line 1452 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1452                          //                      {
1453                          //                              k_c_max = k_c_curr;
1454                          //                      }
1308 –
1455                          H_c_l[c[cg->getGlobalIndex()]].push_back(/*l[*/cg/*]*/);
1456                  }
1457          }
1458  
1459 <        reorderMemory(H_c_l);
1459 >        /* Frequency of reordering the memory */
1460 >        if (neighborListReorderFreq != 0)
1461 >        {
1462 >                if (reorderFreqCounter == neighborListReorderFreq)
1463 >                {
1464 >                        //printf("neighborListReorderFreq:%d\n", neighborListReorderFreq);
1465 >                        reorderMemory(H_c_l);
1466 >                        reorderFreqCounter = 1;
1467 >                } else
1468 >                {
1469 >                        reorderFreqCounter++;
1470 >                }
1471 >        }
1472  
1473          int m;
1474          /* the neighbor matrix */

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines