// // 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 "parallel/PetscSolverFeti.h" #include "parallel/StdMpi.h" #include "parallel/MpiHelper.h" namespace AMDiS { using namespace std; #ifdef HAVE_PETSC_DEV // y = mat * x int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y) { // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi void *ctx; MatShellGetContext(mat, &ctx); PetscSchurPrimalData* data = static_cast(ctx); MatMult(*(data->mat_b_primal), x, data->tmp_vec_b); KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b); MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal); MatMult(*(data->mat_primal_primal), x, y); VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal); return 0; } // y = mat * x int petscMultMatFeti(Mat mat, Vec x, Vec y) { // F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L) void *ctx; MatShellGetContext(mat, &ctx); PetscFetiData* data = static_cast(ctx); // y = L inv(K_BB) trans(L) x MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b); KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b); MatMult(*(data->mat_lagrange), data->tmp_vec_b, y); // tmp_vec_primal = inv(S_PiPi) K_PiB inv(K_BB) trans(L) MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal); KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal); // tmp_vec_lagrange = L inv(K_BB) K_BPi tmp_vec_primal // = L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L) MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b); KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b); MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange); VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange); return 0; } void PetscSolverFeti::updateDofData() { FUNCNAME("PetscSolverFeti::updateDofData()"); TEST_EXIT(meshDistributor->getMesh()->getDim() == 2) ("Works for 2D problems only!"); TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1) ("Works for linear basis functions only!\n"); createPrimals(); createDuals(); createLagrange(); createIndexB(); } void PetscSolverFeti::createPrimals() { FUNCNAME("PetscSolverFeti::createPrimals()"); // === Define all vertices on the interior boundaries of the macro mesh === // === to be primal variables. === primals.clear(); DofContainerSet& vertices = meshDistributor->getBoundaryDofInfo().geoDofs[VERTEX]; TEST_EXIT_DBG(vertices.size())("No primal vertices on this rank!\n"); for (DofContainerSet::iterator it = vertices.begin(); it != vertices.end(); ++it) primals.insert(**it); // === Calculate the number of primals that are owned by the rank and === // === create local indices of the primals starting at zero. === globalPrimalIndex.clear(); nRankPrimals = 0; for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) if (meshDistributor->getIsRankDof(*it)) { globalPrimalIndex[*it] = nRankPrimals; nRankPrimals++; } // === Get overall number of primals and rank's displacement in the === // === numbering of the primals. === nOverallPrimals = 0; rStartPrimals = 0; mpi::getDofNumbering(meshDistributor->getMpiComm(), nRankPrimals, rStartPrimals, nOverallPrimals); // === Create global primal index for all primals. === for (DofMapping::iterator it = globalPrimalIndex.begin(); it != globalPrimalIndex.end(); ++it) it->second += rStartPrimals; MSG_DBG("nRankPrimals = %d nOverallPrimals = %d\n", nRankPrimals, nOverallPrimals); // === Communicate primal's global index from ranks that own the === // === primals to ranks that contain this primals but are not owning === // === them. === StdMpi > stdMpi(meshDistributor->getMpiComm()); RankToDofContainer& sendDofs = meshDistributor->getSendDofs(); for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) if (globalPrimalIndex.count(**dofIt)) stdMpi.getSendData(it->first).push_back(globalPrimalIndex[**dofIt]); stdMpi.updateSendDataSize(); RankToDofContainer& recvDofs = meshDistributor->getRecvDofs(); for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { bool recvFromRank = false; for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) if (primals.count(**dofIt) && meshDistributor->getIsRankDof(**dofIt) == false) { recvFromRank = true; break; } if (recvFromRank) stdMpi.recv(it->first); } stdMpi.startCommunication(); for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { int i = 0; for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { if (primals.count(**dofIt) && meshDistributor->getIsRankDof(**dofIt) == false) globalPrimalIndex[**dofIt] = stdMpi.getRecvData(it->first)[i++]; } } TEST_EXIT_DBG(primals.size() == globalPrimalIndex.size()) ("Number of primals %d, but number of global primals on this rank is %d!\n", primals.size(), globalPrimalIndex.size()); TEST_EXIT_DBG(nOverallPrimals > 0) ("There are zero primal nodes in domain!\n"); } void PetscSolverFeti::createDuals() { FUNCNAME("PetscSolverFeti::createDuals()"); // === Create for each dual node that is owned by the rank, the set === // === of ranks that contain this node (denoted by W(x_j)). === boundaryDofRanks.clear(); RankToDofContainer& sendDofs = meshDistributor->getSendDofs(); for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) { for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { // If DOF is not primal, i.e., its a dual node if (primals.count(**dofIt) == 0) { boundaryDofRanks[**dofIt].insert(mpiRank); boundaryDofRanks[**dofIt].insert(it->first); } } } // === Communicate these sets for all rank owned dual nodes to other === // === ranks that also have this node. === StdMpi > > stdMpi(meshDistributor->getMpiComm()); for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) if (primals.count(**dofIt) == 0) stdMpi.getSendData(it->first).push_back(boundaryDofRanks[**dofIt]); stdMpi.updateSendDataSize(); RankToDofContainer& recvDofs = meshDistributor->getRecvDofs(); for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { bool recvFromRank = false; for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) if (primals.count(**dofIt) == 0) { recvFromRank = true; break; } if (recvFromRank) stdMpi.recv(it->first); } stdMpi.startCommunication(); for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { int i = 0; for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) if (primals.count(**dofIt) == 0) boundaryDofRanks[**dofIt] = stdMpi.getRecvData(it->first)[i++]; } // === Create global index of the dual nodes on each rank. === duals.clear(); globalDualIndex.clear(); int nRankAllDofs = meshDistributor->getFeSpace()->getAdmin()->getUsedDofs(); nRankB = nRankAllDofs - primals.size(); nOverallB = 0; rStartB = 0; mpi::getDofNumbering(meshDistributor->getMpiComm(), nRankB, rStartB, nOverallB); DofContainer allBoundaryDofs; meshDistributor->getAllBoundaryDofs(allBoundaryDofs); int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size(); int nRankDuals = 0; for (DofContainer::iterator it = allBoundaryDofs.begin(); it != allBoundaryDofs.end(); ++it) { if (primals.count(**it) == 0) { duals.insert(**it); globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals; nRankDuals++; } } int nOverallDuals = nRankDuals; mpi::globalAdd(nOverallDuals); MSG_DBG("nRankDuals = %d nOverallDuals = %d\n", nRankDuals, nOverallDuals); } void PetscSolverFeti::createLagrange() { FUNCNAME("PetscSolverFeti::createLagrange()"); // === Reserve for each dual node, on the rank that owns this node, the === // === appropriate number of Lagrange constraints. === nRankLagrange = 0; for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) { if (meshDistributor->getIsRankDof(*it)) { dofFirstLagrange[*it] = nRankLagrange; int degree = boundaryDofRanks[*it].size(); nRankLagrange += (degree * (degree - 1)) / 2; } } // === Get the overall number of Lagrange constraints and create the === // === mapping dofFirstLagrange, that defines for each dual boundary === // === node the first Lagrange constraint global index. === nOverallLagrange = 0; rStartLagrange = 0; mpi::getDofNumbering(meshDistributor->getMpiComm(), nRankLagrange, rStartLagrange, nOverallLagrange); for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) if (meshDistributor->getIsRankDof(*it)) dofFirstLagrange[*it] += rStartLagrange; MSG_DBG("nRankLagrange = %d nOverallLagrange = %d\n", nRankLagrange, nOverallLagrange); // === Communicate dofFirstLagrange to all other ranks. === StdMpi > stdMpi(meshDistributor->getMpiComm()); RankToDofContainer& sendDofs = meshDistributor->getSendDofs(); for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { if (primals.count(**dofIt) == 0) { TEST_EXIT_DBG(dofFirstLagrange.count(**dofIt))("Should not happen!\n"); stdMpi.getSendData(it->first).push_back(dofFirstLagrange[**dofIt]); } } stdMpi.updateSendDataSize(); RankToDofContainer& recvDofs = meshDistributor->getRecvDofs(); for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { bool recvData = false; for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) if (primals.count(**dofIt) == 0) { recvData = true; break; } if (recvData) stdMpi.recv(it->first); } stdMpi.startCommunication(); for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { int counter = 0; for (unsigned int i = 0; i < it->second.size(); i++) if (primals.count(*(it->second[i])) == 0) dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++]; } } void PetscSolverFeti::createIndexB() { FUNCNAME("PetscSolverFeti::createIndeB()"); globalIndexB.clear(); DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin(); // === To ensure that all interior node on each rank are listen first in === // === the global index of all B nodes, insert all interior nodes first, === // === without defining a correct index. === for (int i = 0; i < admin->getUsedSize(); i++) if (admin->isDofFree(i) == false && primals.count(i) == 0) if (duals.count(i) == 0 && primals.count(i) == 0) globalIndexB[i] = -1; // === Get correct index for all interior nodes. === int nInterior = 0; for (DofMapping::iterator it = globalIndexB.begin(); it != globalIndexB.end(); ++it) { it->second = nInterior + rStartB; nInterior++; } TEST_EXIT_DBG(nInterior + primals.size() + duals.size() == static_cast(admin->getUsedDofs())) ("Should not happen!\n"); // === And finally, add the global indicies of all dual nodes. === for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) globalIndexB[*it] = globalDualIndex[*it]; } void PetscSolverFeti::createMatLagrange() { FUNCNAME("PetscSolverFeti::createMatLagrange()"); // === Create distributed matrix for Lagrange constraints. === MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankLagrange * nComponents, nRankB * nComponents, nOverallLagrange * nComponents, nOverallB * nComponents, 2, PETSC_NULL, 2, PETSC_NULL, &mat_lagrange); // === Create for all duals the corresponding Lagrange constraints. On === // === each rank we traverse all pairs (n, m) of ranks, with n < m, === // === that contain this node. If the current rank number is r, and === // === n == r, the rank sets 1.0 for the corresponding constraint, if === // === m == r, than the rank sets -1.0 for the corresponding === // === constraint. === for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) { TEST_EXIT_DBG(dofFirstLagrange.count(*it))("Should not happen!\n"); TEST_EXIT_DBG(boundaryDofRanks.count(*it))("Should not happen!\n"); // Global index of the first Lagrange constriant for this node. int index = dofFirstLagrange[*it]; // Copy set of all ranks that contain this dual node. vector W(boundaryDofRanks[*it].begin(), boundaryDofRanks[*it].end()); // Number of ranks that contain this dual node. int degree = W.size(); TEST_EXIT_DBG(globalDualIndex.count(*it))("Should not happen!\n"); int dualCol = globalDualIndex[*it]; for (int i = 0; i < degree; i++) { for (int j = i + 1; j < degree; j++) { if (W[i] == mpiRank || W[j] == mpiRank) { // Set the constraint for all components of the system. for (int k = 0; k < nComponents; k++) { int rowIndex = index * nComponents + k; int colIndex = dualCol * nComponents + k; double value = (W[i] == mpiRank ? 1.0 : -1.0); MatSetValue(mat_lagrange, rowIndex, colIndex, value, INSERT_VALUES); } } index++; } } } MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY); } void PetscSolverFeti::createSchurPrimalKsp() { FUNCNAME("PetscSolverFeti::createSchurPrimal()"); petscSchurPrimalData.mat_primal_primal = &mat_primal_primal; petscSchurPrimalData.mat_primal_b = &mat_primal_b; petscSchurPrimalData.mat_b_primal = &mat_b_primal; petscSchurPrimalData.ksp_b = &ksp_b; VecDuplicate(f_b, &(petscSchurPrimalData.tmp_vec_b)); VecDuplicate(f_primal, &(petscSchurPrimalData.tmp_vec_primal)); MatCreateShell(PETSC_COMM_WORLD, nRankPrimals * nComponents, nRankPrimals * nComponents, nOverallPrimals * nComponents, nOverallPrimals * nComponents, &petscSchurPrimalData, &mat_schur_primal); MatShellSetOperation(mat_schur_primal, MATOP_MULT, (void(*)(void))petscMultMatSchurPrimal); KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal); KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_"); KSPSetFromOptions(ksp_schur_primal); } void PetscSolverFeti::destroySchurPrimalKsp() { FUNCNAME("PetscSolverFeti::destroySchurPrimal()"); petscSchurPrimalData.mat_primal_primal = PETSC_NULL; petscSchurPrimalData.mat_primal_b = PETSC_NULL; petscSchurPrimalData.mat_b_primal = PETSC_NULL; petscSchurPrimalData.ksp_b = PETSC_NULL; VecDestroy(petscSchurPrimalData.tmp_vec_b); VecDestroy(petscSchurPrimalData.tmp_vec_primal); MatDestroy(mat_schur_primal); KSPDestroy(ksp_schur_primal); } void PetscSolverFeti::createFetiKsp() { FUNCNAME("PetscSolverFeti::createFetiKsp()"); petscFetiData.mat_primal_primal = &mat_primal_primal; petscFetiData.mat_primal_b = &mat_primal_b; petscFetiData.mat_b_primal = &mat_b_primal; petscFetiData.mat_lagrange = &mat_lagrange; petscFetiData.ksp_b = &ksp_b; petscFetiData.ksp_schur_primal = &ksp_schur_primal; VecDuplicate(f_b, &(petscFetiData.tmp_vec_b)); VecDuplicate(f_primal, &(petscFetiData.tmp_vec_primal)); MatGetVecs(mat_lagrange, PETSC_NULL, &(petscFetiData.tmp_vec_lagrange)); MatCreateShell(PETSC_COMM_WORLD, nRankLagrange, nRankLagrange, nOverallLagrange, nOverallLagrange, &petscFetiData, &mat_feti); MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti); KSPCreate(PETSC_COMM_WORLD, &ksp_feti); KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_feti, "solver_feti_"); KSPSetFromOptions(ksp_feti); } void PetscSolverFeti::destroyFetiKsp() { FUNCNAME("PetscSolverFeti::destroyFetiKsp()"); petscFetiData.mat_primal_primal = PETSC_NULL; petscFetiData.mat_primal_b = PETSC_NULL; petscFetiData.mat_b_primal = PETSC_NULL; petscFetiData.mat_lagrange = PETSC_NULL; petscFetiData.ksp_b = PETSC_NULL; petscFetiData.ksp_schur_primal = PETSC_NULL; VecDestroy(petscFetiData.tmp_vec_b); VecDestroy(petscFetiData.tmp_vec_primal); VecDestroy(petscFetiData.tmp_vec_lagrange); MatDestroy(mat_feti); KSPDestroy(ksp_feti); } void PetscSolverFeti::recoverSolution(Vec &vec_sol_b, Vec &vec_sol_primal, SystemVector &vec) { FUNCNAME("PetscSolverFeti::recoverSolution()"); // === Get local part of the solution for B variables. === PetscScalar *localSolB; VecGetArray(vec_sol_b, &localSolB); // === Create scatter to get solutions of all primal nodes that are === // === contained in rank's domain. === vector globalIsIndex, localIsIndex; globalIsIndex.reserve(globalPrimalIndex.size() * nComponents); localIsIndex.reserve(globalPrimalIndex.size() * nComponents); { int counter = 0; for (DofMapping::iterator it = globalPrimalIndex.begin(); it != globalPrimalIndex.end(); ++it) { for (int i = 0; i < nComponents; i++) { globalIsIndex.push_back(it->second * nComponents + i); localIsIndex.push_back(counter++); } } } IS globalIs, localIs; ISCreateGeneral(PETSC_COMM_SELF, globalIsIndex.size(), &(globalIsIndex[0]), PETSC_USE_POINTER, &globalIs); ISCreateGeneral(PETSC_COMM_SELF, localIsIndex.size(), &(localIsIndex[0]), PETSC_USE_POINTER, &localIs); Vec local_sol_primal; VecCreateSeq(PETSC_COMM_SELF, localIsIndex.size(), &local_sol_primal); VecScatter primalScatter; VecScatterCreate(vec_sol_primal, globalIs, local_sol_primal, localIs, &primalScatter); VecScatterBegin(primalScatter, vec_sol_primal, local_sol_primal, INSERT_VALUES, SCATTER_FORWARD); VecScatterEnd(primalScatter, vec_sol_primal, local_sol_primal, INSERT_VALUES, SCATTER_FORWARD); ISDestroy(globalIs); ISDestroy(localIs); VecScatterDestroy(primalScatter); PetscScalar *localSolPrimal; VecGetArray(local_sol_primal, &localSolPrimal); // === And copy from PETSc local vectors to the DOF vectors. === for (int i = 0; i < nComponents; i++) { DOFVector& dofVec = *(vec.getDOFVector(i)); for (DofMapping::iterator it = globalIndexB.begin(); it != globalIndexB.end(); ++it) { int petscIndex = (it->second - rStartB) * nComponents + i; dofVec[it->first] = localSolB[petscIndex]; } int counter = 0; for (DofMapping::iterator it = globalPrimalIndex.begin(); it != globalPrimalIndex.end(); ++it) { dofVec[it->first] = localSolPrimal[counter * nComponents + i]; counter++; } } VecRestoreArray(vec_sol_b, &localSolB); VecRestoreArray(local_sol_primal, &localSolPrimal); VecDestroy(local_sol_primal); } void PetscSolverFeti::fillPetscMatrix(Matrix *mat, SystemVector *vec) { FUNCNAME("PetscSolverFeti::fillPetscMatrix()"); nComponents = vec->getSize(); // === Create all sets and indices. === updateDofData(); // === Create matrices for the FETI-DP method. === int nRowsRankB = nRankB * nComponents; int nRowsOverallB = nOverallB * nComponents; int nRowsRankPrimal = nRankPrimals * nComponents; int nRowsOverallPrimal = nOverallPrimals * nComponents; MatCreateMPIAIJ(PETSC_COMM_WORLD, nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB, 100, PETSC_NULL, 100, PETSC_NULL, &mat_b_b); MatCreateMPIAIJ(PETSC_COMM_WORLD, nRowsRankPrimal, nRowsRankPrimal, nRowsOverallPrimal, nRowsOverallPrimal, 10, PETSC_NULL, 10, PETSC_NULL, &mat_primal_primal); MatCreateMPIAIJ(PETSC_COMM_WORLD, nRowsRankB, nRowsRankPrimal, nRowsOverallB, nRowsOverallPrimal, 100, PETSC_NULL, 100, PETSC_NULL, &mat_b_primal); MatCreateMPIAIJ(PETSC_COMM_WORLD, nRowsRankPrimal, nRowsRankB, nRowsOverallPrimal, nRowsOverallB, 100, PETSC_NULL, 100, PETSC_NULL, &mat_primal_b); // === Prepare traverse of sequentially created matrices. === 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 traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; vector cols, colsOther; vector values, valuesOther; cols.reserve(300); colsOther.reserve(300); values.reserve(300); valuesOther.reserve(300); // === Traverse all sequentially created matrices and add the values to === // === the global PETSc matrices. === for (int i = 0; i < nComponents; i++) { for (int j = 0; j < nComponents; j++) { if (!(*mat)[i][j]) continue; traits::col::type col((*mat)[i][j]->getBaseMatrix()); traits::const_value::type value((*mat)[i][j]->getBaseMatrix()); // Traverse all rows. for (cursor_type cursor = begin((*mat)[i][j]->getBaseMatrix()), cend = end((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) { bool rowPrimal = primals.count(*cursor) != 0; cols.clear(); values.clear(); colsOther.clear(); valuesOther.clear(); // Traverse all columns. for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { if (primals.count(col(*icursor)) != 0) { // Column is a primal variable. TEST_EXIT_DBG(globalPrimalIndex.count(col(*icursor))) ("No global primal index for DOF %d!\n", col(*icursor)); int colIndex = globalPrimalIndex[col(*icursor)] * nComponents + j; if (rowPrimal) { cols.push_back(colIndex); values.push_back(value(*icursor)); } else { colsOther.push_back(colIndex); valuesOther.push_back(value(*icursor)); } } else { // Column is not a primal variable. TEST_EXIT_DBG(globalIndexB.count(col(*icursor))) ("No global B index for DOF %d!\n", col(*icursor)); int colIndex = globalIndexB[col(*icursor)] * nComponents + j; if (rowPrimal) { colsOther.push_back(colIndex); valuesOther.push_back(value(*icursor)); } else { cols.push_back(colIndex); values.push_back(value(*icursor)); } } } if (rowPrimal) { TEST_EXIT_DBG(globalPrimalIndex.count(*cursor)) ("Should not happen!\n"); int rowIndex = globalPrimalIndex[*cursor] * nComponents + i; MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); if (colsOther.size()) MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(), &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); } else { TEST_EXIT_DBG(globalIndexB.count(*cursor)) ("Should not happen!\n"); int rowIndex = globalIndexB[*cursor] * nComponents + i; MatSetValues(mat_b_b, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); if (colsOther.size()) MatSetValues(mat_b_primal, 1, &rowIndex, colsOther.size(), &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); } } } } // === Start global assembly procedure. === MatAssemblyBegin(mat_b_b, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_b_b, MAT_FINAL_ASSEMBLY); MatAssemblyBegin(mat_primal_primal, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_primal_primal, MAT_FINAL_ASSEMBLY); MatAssemblyBegin(mat_b_primal, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_b_primal, MAT_FINAL_ASSEMBLY); MatAssemblyBegin(mat_primal_b, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_primal_b, MAT_FINAL_ASSEMBLY); // === Create and fill PETSc's right hand side vectors. === VecCreate(PETSC_COMM_WORLD, &f_b); VecSetSizes(f_b, nRankB * nComponents, nOverallB * nComponents); VecSetType(f_b, VECMPI); VecCreate(PETSC_COMM_WORLD, &f_primal); VecSetSizes(f_primal, nRankPrimals * nComponents, nOverallPrimals * nComponents); VecSetType(f_primal, VECMPI); for (int i = 0; i < nComponents; i++) { DOFVector::Iterator dofIt(vec->getDOFVector(i), USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { int index = dofIt.getDOFIndex(); if (primals.count(index)) { TEST_EXIT_DBG(globalPrimalIndex.count(index)) ("Should not happen!\n"); index = globalPrimalIndex[index] * nComponents + i; double value = *dofIt; VecSetValues(f_primal, 1, &index, &value, ADD_VALUES); } else { TEST_EXIT_DBG(globalIndexB.count(index)) ("Should not happen!\n"); index = globalIndexB[index] * nComponents + i; double value = *dofIt; VecSetValues(f_b, 1, &index, &value, ADD_VALUES); } } } VecAssemblyBegin(f_b); VecAssemblyEnd(f_b); VecAssemblyBegin(f_primal); VecAssemblyEnd(f_primal); // === Create and fill PETSc matrix for Lagrange constraints. === createMatLagrange(); // === Create PETSc solver for the Schur complement on primal variables. === createSchurPrimalKsp(); // === Create PETSc solver for the FETI-DP operator. === createFetiKsp(); } void PetscSolverFeti::solveFetiMatrix(SystemVector &vec) { FUNCNAME("PetscSolverFeti::solveFetiMatrix()"); // Create transpose of Lagrange matrix. Mat mat_lagrange_transpose; MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose); // === Create nested matrix which will contain the overall FETI system. === Mat A; Mat nestedA[3][3]; nestedA[0][0] = mat_b_b; nestedA[0][1] = mat_b_primal; nestedA[0][2] = mat_lagrange_transpose; nestedA[1][0] = mat_primal_b; nestedA[1][1] = mat_primal_primal; nestedA[1][2] = PETSC_NULL; nestedA[2][0] = mat_lagrange; nestedA[2][1] = PETSC_NULL; nestedA[2][2] = PETSC_NULL; MatCreateNest(PETSC_COMM_WORLD, 3, PETSC_NULL, 3, PETSC_NULL, &(nestedA[0][0]), &A); MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); int nRankNest = (nRankB + nRankPrimals) * nComponents + nRankLagrange; int nOverallNest = (nOverallB + nOverallPrimals) * nComponents + nOverallLagrange; int rStartNest = (rStartB + rStartPrimals) * nComponents + rStartLagrange; { // === Test some matrix sizes. === int matRow, matCol; MatGetLocalSize(A, &matRow, &matCol); TEST_EXIT_DBG(matRow == nRankNest)("Should not happen!\n"); mpi::globalAdd(matRow); TEST_EXIT_DBG(matRow == nOverallNest)("Should not happen!\n"); MatGetOwnershipRange(A, &matRow, &matCol); TEST_EXIT_DBG(matRow == rStartNest)("Should not happen!\n"); } // === Create rhs and solution vectors for the overall FETI system. === Vec f; VecCreate(PETSC_COMM_WORLD, &f); VecSetSizes(f, nRankNest, nOverallNest); VecSetType(f, VECMPI); Vec b; VecDuplicate(f, &b); // === Fill rhs vector by coping the primal and non primal PETSc vectors. === PetscScalar *local_f_b; VecGetArray(f_b, &local_f_b); PetscScalar *local_f_primal; VecGetArray(f_primal, &local_f_primal); { int size; VecGetLocalSize(f_b, &size); TEST_EXIT_DBG(size == nRankB * nComponents)("Should not happen!\n"); VecGetLocalSize(f_primal, &size); TEST_EXIT_DBG(size == nRankPrimals * nComponents)("Should not happen!\n"); } PetscScalar *local_f; VecGetArray(f, &local_f); int index = 0; for (int i = 0; i < nRankB * nComponents; i++) local_f[index++] = local_f_b[i]; for (int i = 0; i < nRankPrimals * nComponents; i++) local_f[index++] = local_f_primal[i]; VecRestoreArray(f, &local_f); VecRestoreArray(f_b, &local_f_b); VecRestoreArray(f_primal, &local_f_primal); // === Create solver and solve the overall FETI system. === KSP ksp; KSPCreate(PETSC_COMM_WORLD, &ksp); KSPSetOperators(ksp, A, A, SAME_NONZERO_PATTERN); KSPSetFromOptions(ksp); KSPSolve(ksp, f, b); // === Reconstruct FETI solution vectors. === Vec u_b, u_primal; VecDuplicate(f_b, &u_b); VecDuplicate(f_primal, &u_primal); PetscScalar *local_b; VecGetArray(b, &local_b); PetscScalar *local_u_b; VecGetArray(u_b, &local_u_b); PetscScalar *local_u_primal; VecGetArray(u_primal, &local_u_primal); index = 0; for (int i = 0; i < nRankB * nComponents; i++) local_u_b[i] = local_b[index++]; for (int i = 0; i < nRankPrimals * nComponents; i++) local_u_primal[i] = local_b[index++]; VecRestoreArray(f, &local_b); VecRestoreArray(u_b, &local_u_b); VecRestoreArray(u_primal, &local_u_primal); recoverSolution(u_b, u_primal, vec); VecDestroy(u_b); VecDestroy(u_primal); VecDestroy(b); VecDestroy(f); KSPDestroy(ksp); } void PetscSolverFeti::solveReducedFetiMatrix(SystemVector &vec) { FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()"); // === Create solver for the non primal (thus local) variables. === KSPCreate(PETSC_COMM_WORLD, &ksp_b); KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_b, "solver_b_"); KSPSetFromOptions(ksp_b); // RHS and solution vector. Vec vec_rhs; // Some temporary vectors. Vec tmp_b0, tmp_b1, tmp_lagrange0, tmp_primal0, tmp_primal1; MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange0); MatGetVecs(mat_lagrange, PETSC_NULL, &vec_rhs); MatGetVecs(mat_b_b, PETSC_NULL, &tmp_b0); MatGetVecs(mat_b_b, PETSC_NULL, &tmp_b1); MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal0); MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal1); // === Create new rhs === // vec_rhs = L * inv(K_BB) * f_b KSPSolve(ksp_b, f_b, tmp_b0); MatMult(mat_lagrange, tmp_b0, vec_rhs); // tmp_primal0 = M_PiB * inv(K_BB) * f_b MatMult(mat_primal_b, tmp_b0, tmp_primal0); // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_b VecAXPBY(tmp_primal0, -1.0, 1.0, f_primal); // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_b) KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0); // MatMult(mat_b_primal, tmp_primal0, tmp_b0); KSPSolve(ksp_b, tmp_b0, tmp_b0); MatMult(mat_lagrange, tmp_b0, tmp_lagrange0); // VecAXPBY(vec_rhs, 1.0, 1.0, tmp_lagrange0); // === Solve with FETI-DP operator. === KSPSolve(ksp_feti, vec_rhs, vec_rhs); // === Solve for u_primals. === VecCopy(f_primal, tmp_primal0); KSPSolve(ksp_b, f_b, tmp_b0); MatMult(mat_primal_b, tmp_b0, tmp_primal1); VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1); MatMultTranspose(mat_lagrange, vec_rhs, tmp_b0); KSPSolve(ksp_b, tmp_b0, tmp_b0); MatMult(mat_primal_b, tmp_b0, tmp_primal1); VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1); KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0); // === Solve for u_b. === VecCopy(f_b, tmp_b0); MatMultTranspose(mat_lagrange, vec_rhs, tmp_b1); VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1); MatMult(mat_b_primal, tmp_primal0, tmp_b1); VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1); KSPSolve(ksp_b, tmp_b0, tmp_b0); // === And recover AMDiS solution vectors. === recoverSolution(tmp_b0, tmp_primal0, vec); // === Destroy all data structures. === VecDestroy(vec_rhs); VecDestroy(tmp_b0); VecDestroy(tmp_b1); VecDestroy(tmp_lagrange0); VecDestroy(tmp_primal0); VecDestroy(tmp_primal1); KSPDestroy(ksp_b); MatDestroy(mat_b_b); MatDestroy(mat_primal_primal); MatDestroy(mat_b_primal); MatDestroy(mat_primal_b); MatDestroy(mat_lagrange); VecDestroy(f_b); VecDestroy(f_primal); destroySchurPrimalKsp(); destroyFetiKsp(); } void PetscSolverFeti::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo) { FUNCNAME("PetscSolverFeti::solvePetscMatrix()"); int debug = 0; Parameters::get("parallel->debug feti", debug); if (debug) { WARNING("FETI matrix is solved globally, thus without reducing to the lagrange multipliers!\n"); solveFetiMatrix(vec); } else { solveReducedFetiMatrix(vec); } } #endif }