// // Software License for AMDiS // // Copyright (c) 2010 Dresden University of Technology // All rights reserved. // Authors: Simon Vey, Thomas Witkowski et al. // // This file is part of AMDiS // // See also license.opensource.txt in the distribution. #include #include #include "parallel/PetscSolver.h" #include "parallel/StdMpi.h" #include "parallel/ParallelDebug.h" #include "DOFVector.h" #include "Debug.h" #include "SystemVector.h" #include "petscksp.h" namespace AMDiS { using namespace std; PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *) { if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0) cout << "[0] Petsc-Iteration " << iter << ": " << rnorm << endl; return 0; } void PetscSolver::solve(AdaptInfo *adaptInfo, bool fixedMatrix) { FUNCNAME("PetscSolver::solve()"); TEST_EXIT(meshDistributor)("Should not happen!\n"); double wtime = MPI::Wtime(); fillPetscMatrix(systemMatrix, rhs); solvePetscMatrix(*solution, adaptInfo); INFO(info, 8)("solution of discrete system needed %.5f seconds\n", MPI::Wtime() - wtime); } void PetscSolver::setDofMatrix(DOFMatrix* mat, int dispMult, int dispAddRow, int dispAddCol) { FUNCNAME("PetscSolver::setDofMatrix()"); TEST_EXIT(mat)("No DOFMatrix!\n"); using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits= mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; traits::col::type col(mat->getBaseMatrix()); traits::const_value::type value(mat->getBaseMatrix()); typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; vector cols; vector values; cols.reserve(300); values.reserve(300); vector globalCols; // === Traverse all rows of the dof matrix and insert row wise the values === // === to the PETSc matrix. === for (cursor_type cursor = begin(mat->getBaseMatrix()), cend = end(mat->getBaseMatrix()); cursor != cend; ++cursor) { // Global index of the current row dof. int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor); // Test if the current row dof is a periodic dof. bool periodicRow = meshDistributor->isPeriodicDof(globalRowDof); if (!periodicRow) { // === Row DOF index is not periodic. === // Calculate PETSc row index. int rowIndex = globalRowDof * dispMult + dispAddRow; cols.clear(); values.clear(); for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { // Global index of the current column index. int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor)); // Test if the current col dof is a periodic dof. bool periodicCol = meshDistributor->isPeriodicDof(globalColDof); // Calculate PETSc col index. int colIndex = globalColDof * dispMult + dispAddCol; // Ignore all zero entries, expect it is a diagonal entry. // if (value(*icursor) == 0.0 && globalRowDof != globalColDof) // if (value(*icursor) == 0.0 && rowIndex != colIndex) // continue; if (!periodicCol) { // Calculate the exact position of the column index in the PETSc matrix. cols.push_back(globalColDof * dispMult + dispAddCol); values.push_back(value(*icursor)); } else { // === Row index is not periodic, but column index is. === // Create set of all periodic associations of the column index. std::set perAsc; std::set& perColAsc = meshDistributor->getPerDofAssociations(globalColDof); for (std::set::iterator it = perColAsc.begin(); it != perColAsc.end(); ++it) if (*it >= -3) perAsc.insert(*it); // Scale value to the number of periodic associations of the column index. double scaledValue = value(*icursor) * pow(0.5, static_cast(perAsc.size())); // === Create set of all matrix column indices due to the periodic === // === associations of the column DOF index. === vector newCols; // First, add the original matrix index. newCols.push_back(globalColDof); // And add all periodic matrix indices. for (std::set::iterator it = perAsc.begin(); it != perAsc.end(); ++it) { int nCols = static_cast(newCols.size()); for (int i = 0; i < nCols; i++) { TEST_EXIT_DBG(meshDistributor->isPeriodicDof(newCols[i], *it)) ("Wrong periodic DOF associations at boundary %d with DOF %d!\n", *it, newCols[i]); newCols.push_back(meshDistributor->getPeriodicMapping(newCols[i], *it)); } } for (unsigned int i = 0; i < newCols.size(); i++) { cols.push_back(newCols[i] * dispMult + dispAddCol); values.push_back(scaledValue); } } } MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); } else { // === Row DOF index is periodic. === // Because this row is periodic, we will have to add the entries of this // matrix row to multiple rows. The following maps store to each row an // array of column indices and values of the entries that must be added to // the PETSc matrix. map > colsMap; map > valsMap; // Traverse all column entries. for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { // Global index of the current column index. int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor)); // Ignore all zero entries, expect it is a diagonal entry. // if (value(*icursor) == 0.0 && globalRowDof != globalColDof) // continue; // === Add all periodic associations of both, the row and the column === // === indices to the set perAsc. === std::set perAsc; if (meshDistributor->isPeriodicDof(globalColDof)) { std::set& perColAsc = meshDistributor->getPerDofAssociations(globalColDof); for (std::set::iterator it = perColAsc.begin(); it != perColAsc.end(); ++it) if (*it >= -3) perAsc.insert(*it); } std::set& perRowAsc = meshDistributor->getPerDofAssociations(globalRowDof); for (std::set::iterator it = perRowAsc.begin(); it != perRowAsc.end(); ++it) if (*it >= -3) perAsc.insert(*it); // Scale the value with respect to the number of periodic associations. double scaledValue = value(*icursor) * pow(0.5, static_cast(perAsc.size())); // === Create all matrix entries with respect to the periodic === // === associations of the row and column indices. === vector > entry; // First, add the original entry. entry.push_back(make_pair(globalRowDof, globalColDof)); // Then, traverse the periodic associations of the row and column indices // and create the corresponding entries. for (std::set::iterator it = perAsc.begin(); it != perAsc.end(); ++it) { int nEntry = static_cast(entry.size()); for (int i = 0; i < nEntry; i++) { int perRowDof = 0; if (meshDistributor->getPeriodicMapping()[*it].count(entry[i].first)) perRowDof = meshDistributor->getPeriodicMapping(entry[i].first, *it); else perRowDof = entry[i].first; int perColDof; if (meshDistributor->getPeriodicMapping()[*it].count(entry[i].second)) perColDof = meshDistributor->getPeriodicMapping(entry[i].second, *it); else perColDof = entry[i].second; entry.push_back(make_pair(perRowDof, perColDof)); } } // === Translate the matrix entries to PETSc's matrix. for (vector >::iterator eIt = entry.begin(); eIt != entry.end(); ++eIt) { // Calculate row and column indices of the PETSc matrix. int rowIndex = eIt->first * dispMult + dispAddRow; int colIndex = eIt->second * dispMult + dispAddCol; colsMap[rowIndex].push_back(colIndex); valsMap[rowIndex].push_back(scaledValue); } } // === Finally, add all periodic rows to the PETSc matrix. === for (map >::iterator rowIt = colsMap.begin(); rowIt != colsMap.end(); ++rowIt) { TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size()) ("Should not happen!\n"); int rowIndex = rowIt->first; MatSetValues(petscMatrix, 1, &rowIndex, rowIt->second.size(), &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES); } } } } void PetscSolver::setDofVector(Vec& petscVec, DOFVector* vec, int dispMult, int dispAdd) { FUNCNAME("PetscSolver::setDofVector()"); // Traverse all used DOFs in the dof vector. DOFVector::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { // Calculate global row index of the dof. DegreeOfFreedom globalRowDof = meshDistributor->mapLocalToGlobal(dofIt.getDOFIndex()); // Calculate PETSc index of the row dof. int index = globalRowDof * dispMult + dispAdd; if (meshDistributor->isPeriodicDof(globalRowDof)) { std::set& perAsc = meshDistributor->getPerDofAssociations(globalRowDof); double value = *dofIt / (perAsc.size() + 1.0); VecSetValues(petscVec, 1, &index, &value, ADD_VALUES); for (std::set::iterator perIt = perAsc.begin(); perIt != perAsc.end(); ++perIt) { int mappedDof = meshDistributor->getPeriodicMapping(globalRowDof, *perIt); int mappedIndex = mappedDof * dispMult + dispAdd; VecSetValues(petscVec, 1, &mappedIndex, &value, ADD_VALUES); } } else { // The dof index is not periodic. double value = *dofIt; VecSetValues(petscVec, 1, &index, &value, ADD_VALUES); } } } void PetscSolver::createPetscNnzStructure(Matrix *mat) { FUNCNAME("PetscSolver::createPetscNnzStructure()"); TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n"); TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n"); int nRankRows = meshDistributor->getNumberRankDofs() * nComponents; d_nnz = new int[nRankRows]; o_nnz = new int[nRankRows]; for (int i = 0; i < nRankRows; i++) { d_nnz[i] = 0; o_nnz[i] = 0; } using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits = mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; typedef vector > MatrixNnzEntry; typedef map RankToDofContainer; // Stores to each rank a list of nnz entries (i.e. pairs of row and column index) // that this rank will send to. These nnz entries will be assembled on this rank, // but because the row DOFs are not DOFs of this rank they will be send to the // owner of the row DOFs. map sendMatrixEntry; // Maps to each DOF that must be send to another rank the rank number of the // receiving rank. map sendDofToRank; // First, create for all ranks we send data to MatrixNnzEntry object with 0 entries. for (RankToDofContainer::iterator it = meshDistributor->getRecvDofs().begin(); it != meshDistributor->getRecvDofs().end(); ++it) { sendMatrixEntry[it->first].resize(0); for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) sendDofToRank[**dofIt] = it->first; } std::set recvFromRank; for (RankToDofContainer::iterator it = meshDistributor->getSendDofs().begin(); it != meshDistributor->getSendDofs().end(); ++it) recvFromRank.insert(it->first); for (int i = 0; i < nComponents; i++) { for (int j = 0; j < nComponents; j++) { if (!(*mat)[i][j]) continue; Matrix bmat = (*mat)[i][j]->getBaseMatrix(); traits::col::type col(bmat); traits::const_value::type value(bmat); typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; for (cursor_type cursor = begin(bmat), cend = end(bmat); cursor != cend; ++cursor) { int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor); vector rows; rows.push_back(globalRowDof); vector rowRank; if (meshDistributor->getIsRankDof(*cursor)) { rowRank.push_back(meshDistributor->getMpiRank()); } else { // Find out who is the member of this DOF. TEST_EXIT_DBG(sendDofToRank.count(*cursor))("Should not happen!\n"); rowRank.push_back(sendDofToRank[*cursor]); } // Map the local row number to the global DOF index and create from it // the global PETSc row index of this DOF. int petscRowIdx = globalRowDof * nComponents + i; if (meshDistributor->getIsRankDof(*cursor)) { // === The current row DOF is a rank dof, so create the corresponding === // === nnz values directly on rank's nnz data. === // This is the local row index of the local PETSc matrix. int localPetscRowIdx = petscRowIdx - meshDistributor->getRstart() * nComponents; TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows) ("Should not happen! \n Debug info: localRowIdx = %d globalRowIndx = %d petscRowIdx = %d localPetscRowIdx = %d rStart = %d nCompontens = %d nRankRows = %d\n", *cursor, meshDistributor->mapLocalToGlobal(*cursor), petscRowIdx, localPetscRowIdx, meshDistributor->getRstart(), nComponents, nRankRows); // Traverse all non zero entries in this row. for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { int petscColIdx = meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j; // if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) { // The row DOF is a rank DOF, if also the column is a rank DOF, // increment the d_nnz values for this row, otherwise the o_nnz value. if (petscColIdx >= meshDistributor->getRstart() * nComponents && petscColIdx < meshDistributor->getRstart() * nComponents + nRankRows) d_nnz[localPetscRowIdx]++; else o_nnz[localPetscRowIdx]++; // } } } else { // === The current row DOF is not a rank dof, i.e., it will be created === // === on this rank, but after this it will be send to another rank === // === matrix. So we need to send also the corresponding nnz structure === // === of this row to the corresponding rank. === // Send all non zero entries to the member of the row DOF. int sendToRank = sendDofToRank[*cursor]; for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { // if (value(*icursor) != 0.0) { int petscColIdx = meshDistributor->mapLocalToGlobal(col(*icursor)) * nComponents + j; sendMatrixEntry[sendToRank]. push_back(make_pair(petscRowIdx, petscColIdx)); // } } } // if (isRankDof[*cursor]) ... else ... } // for each row in mat[i][j] } } // === Send and recv the nnz row structure to/from other ranks. === StdMpi stdMpi(meshDistributor->getMpiComm(), true); stdMpi.send(sendMatrixEntry); for (std::set::iterator it = recvFromRank.begin(); it != recvFromRank.end(); ++it) stdMpi.recv(*it); stdMpi.startCommunication(); // === Evaluate the nnz structure this rank got from other ranks and add it to === // === the PETSc nnz data structure. === for (map::iterator it = stdMpi.getRecvData().begin(); it != stdMpi.getRecvData().end(); ++it) { if (it->second.size() > 0) { for (unsigned int i = 0; i < it->second.size(); i++) { int r = it->second[i].first; int c = it->second[i].second; int localRowIdx = r - meshDistributor->getRstart() * nComponents; TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows) ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n", r, localRowIdx, nRankRows, it->first); if (c < meshDistributor->getRstart() * nComponents || c >= meshDistributor->getRstart() * nComponents + nRankRows) o_nnz[localRowIdx]++; else d_nnz[localRowIdx]++; } } } // The above algorithm for calculating the number of nnz per row over- // approximates the value, i.e., the number is always equal or larger to // the real number of nnz values in the global parallel matrix. For small // matrices, the problem may arise, that the result is larger than the // number of elements in a row. This is fixed in the following. if (nRankRows < 100) for (int i = 0; i < nRankRows; i++) d_nnz[i] = std::min(d_nnz[i], nRankRows); } void PetscSolver::fillPetscMatrix(Matrix *mat, SystemVector *vec) { FUNCNAME("PetscSolver::fillPetscMatrix()"); double wtime = MPI::Wtime(); int nRankRows = meshDistributor->getNumberRankDofs() * nComponents; int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents; // === Create PETSc vector (rhs, solution and a temporary vector). === VecCreate(PETSC_COMM_WORLD, &petscRhsVec); VecSetSizes(petscRhsVec, nRankRows, nOverallRows); VecSetType(petscRhsVec, VECMPI); VecCreate(PETSC_COMM_WORLD, &petscSolVec); VecSetSizes(petscSolVec, nRankRows, nOverallRows); VecSetType(petscSolVec, VECMPI); VecCreate(PETSC_COMM_WORLD, &petscTmpVec); VecSetSizes(petscTmpVec, nRankRows, nOverallRows); VecSetType(petscTmpVec, VECMPI); int recvAllValues = 0; int sendValue = static_cast(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz); meshDistributor->getMpiComm().Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM); if (!d_nnz || recvAllValues != 0) { if (d_nnz) { delete [] d_nnz; d_nnz = NULL; delete [] o_nnz; o_nnz = NULL; } createPetscNnzStructure(mat); lastMeshNnz = meshDistributor->getLastMeshChangeIndex(); } // === Create PETSc matrix with the computed nnz data structure. === MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows, 0, d_nnz, 0, o_nnz, &petscMatrix); #if (DEBUG != 0) INFO(info, 8)("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime); #endif #if (DEBUG != 0) int a, b; MatGetOwnershipRange(petscMatrix, &a, &b); TEST_EXIT(a == meshDistributor->getRstart() * nComponents) ("Wrong matrix ownership range!\n"); TEST_EXIT(b == meshDistributor->getRstart() * nComponents + nRankRows) ("Wrong matrix ownership range!\n"); #endif // === Transfer values from DOF matrices to the PETSc matrix. === for (int i = 0; i < nComponents; i++) for (int j = 0; j < nComponents; j++) if ((*mat)[i][j]) setDofMatrix((*mat)[i][j], nComponents, i, j); #if (DEBUG != 0) INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime); #endif MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY); #if (DEBUG != 0) INFO(info, 8)("Fill petsc matrix 3 needed %.5f seconds\n", TIME_USED(MPI::Wtime(), wtime)); #endif // === Transfer values from DOF vector to the PETSc vector. === for (int i = 0; i < nComponents; i++) setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i); VecAssemblyBegin(petscRhsVec); VecAssemblyEnd(petscRhsVec); INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", MPI::Wtime() - wtime); } void PetscSolver::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo) { FUNCNAME("PetscSolver::solvePetscMatrix()"); #if 0 // Set old solution to be initiual guess for petsc solver. for (int i = 0; i < nComponents; i++) setDofVector(petscSolVec, vec->getDOFVector(i), nComponents, i); VecAssemblyBegin(petscSolVec); VecAssemblyEnd(petscSolVec); #endif // === Init PETSc solver. === KSP solver; PC pc; providePetscSolver(solver, pc); // Do not delete the solution vector, use it for the initial guess. // KSPSetInitialGuessNonzero(solver, PETSC_TRUE); KSPCreate(PETSC_COMM_WORLD, &solver); KSPGetPC(solver, &pc); PCSetType(pc, PCFIELDSPLIT); #if 0 typedef map RankToDofContainer; typedef map DofIndexToBool; std::set boundaryDofsSet; std::set boundaryLocalDofs; RankToDofContainer& sendDofs = meshDistributor->getSendDofs(); for (RankToDofContainer::iterator rankIt = sendDofs.begin(); rankIt != sendDofs.end(); ++rankIt) for (DofContainer::iterator dofIt = rankIt->second.begin(); dofIt != rankIt->second.end(); ++dofIt) { boundaryLocalDofs.insert(**dofIt); for (int i = 0; i < nComponents; i++) boundaryDofsSet.insert(meshDistributor->mapLocalToGlobal(**dofIt) * nComponents + i); } vector boundaryDofs(boundaryDofsSet.begin(), boundaryDofsSet.end()); std::set otherBoundaryLocalDofs; RankToDofContainer& recvDofs = meshDistributor->getRecvDofs(); for (RankToDofContainer::iterator rankIt = recvDofs.begin(); rankIt != recvDofs.end(); ++rankIt) for (DofContainer::iterator dofIt = rankIt->second.begin(); dofIt != rankIt->second.end(); ++dofIt) otherBoundaryLocalDofs.insert(**dofIt); std::set interiorDofsSet; DofIndexToBool& isRankDof = meshDistributor->getIsRankDof(); for (DofIndexToBool::iterator dofIt = isRankDof.begin(); dofIt != isRankDof.end(); ++dofIt) if (dofIt->second && boundaryLocalDofs.count(dofIt->first) == 0 && otherBoundaryLocalDofs.count(dofIt->first) == 0) for (int i = 0; i < nComponents; i++) interiorDofsSet.insert(meshDistributor->mapLocalToGlobal(dofIt->first) * nComponents + i); vector interiorDofs(interiorDofsSet.begin(), interiorDofsSet.end()); IS interiorIs; ISCreateGeneral(PETSC_COMM_WORLD, interiorDofs.size(), &(interiorDofs[0]), PETSC_COPY_VALUES, &interiorIs); PCFieldSplitSetIS(pc, "interior", interiorIs); IS boundaryIs; ISCreateGeneral(PETSC_COMM_WORLD, boundaryDofs.size(), &(boundaryDofs[0]), PETSC_COPY_VALUES, &boundaryIs); PCFieldSplitSetIS(pc, "boundary", boundaryIs); #endif // === Run PETSc. === KSPSolve(solver, petscRhsVec, petscSolVec); // === Transfere values from PETSc's solution vectors to the dof vectors. PetscScalar *vecPointer; VecGetArray(petscSolVec, &vecPointer); int nRankDofs = meshDistributor->getNumberRankDofs(); for (int i = 0; i < nComponents; i++) { DOFVector &dofvec = *(vec.getDOFVector(i)); for (int j = 0; j < nRankDofs; j++) dofvec[meshDistributor->mapLocalToDofIndex(j)] = vecPointer[j * nComponents + i]; } VecRestoreArray(petscSolVec, &vecPointer); // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to more === // === than one partition. === meshDistributor->synchVector(vec); // === Print information about solution process. === int iterations = 0; KSPGetIterationNumber(solver, &iterations); MSG(" Number of iterations: %d\n", iterations); adaptInfo->setSolverIterations(iterations); double norm = 0.0; MatMult(petscMatrix, petscSolVec, petscTmpVec); VecAXPY(petscTmpVec, -1.0, petscRhsVec); VecNorm(petscTmpVec, NORM_2, &norm); MSG(" Residual norm: %e\n", norm); // === Destroy Petsc's variables. === MatDestroy(petscMatrix); VecDestroy(petscRhsVec); VecDestroy(petscSolVec); VecDestroy(petscTmpVec); KSPDestroy(solver); } void PetscSolver::providePetscSolver(KSP &solver, PC &pc) { FUNCNAME("PetscSolver::providePetscSolver()"); KSPCreate(PETSC_COMM_WORLD, &solver); KSPGetPC(solver, &pc); KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT); KSPSetType(solver, KSPBCGS); KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0); KSPSetFromOptions(solver); PCSetFromOptions(pc); } }