// // 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" #include "io/VtkWriter.h" namespace AMDiS { using namespace std; // 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); SchurPrimalData* data = static_cast(ctx); MatMult(*(data->mat_b_primal), x, data->tmp_vec_b); data->fetiSolver->solveLocalProblem(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) // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L) void *ctx; MatShellGetContext(mat, &ctx); FetiData* data = static_cast(ctx); MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b); data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b); MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange); 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); MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b); data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b); MatMult(*(data->mat_lagrange), data->tmp_vec_b, y); VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange); return 0; } // y = PC * x PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y) { // Get data for the preconditioner void *ctx; PCShellGetContext(pc, &ctx); FetiDirichletPreconData* data = static_cast(ctx); // Multiply with scaled Lagrange constraint matrix. MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b); // === Restriction of the B nodes to the boundary nodes. === int nLocalB; int nLocalDuals; VecGetLocalSize(data->tmp_vec_b, &nLocalB); VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals); PetscScalar *local_b, *local_duals; VecGetArray(data->tmp_vec_b, &local_b); VecGetArray(data->tmp_vec_duals0, &local_duals); for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++) local_duals[j] = local_b[i]; VecRestoreArray(data->tmp_vec_b, &local_b); VecRestoreArray(data->tmp_vec_duals0, &local_duals); // === K_DD - K_DI inv(K_II) K_ID === MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1); MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior); KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior); MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0); VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1); // === Prolongation from local dual nodes to B nodes. VecGetArray(data->tmp_vec_b, &local_b); VecGetArray(data->tmp_vec_duals0, &local_duals); for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++) local_b[i] = local_duals[j]; VecRestoreArray(data->tmp_vec_b, &local_b); VecRestoreArray(data->tmp_vec_duals0, &local_duals); // Multiply with scaled Lagrange constraint matrix. MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y); return 0; } // y = PC * x PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y) { // Get data for the preconditioner void *ctx; PCShellGetContext(pc, &ctx); FetiLumpedPreconData* data = static_cast(ctx); // Multiply with scaled Lagrange constraint matrix. MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b); // === Restriction of the B nodes to the boundary nodes. === int nLocalB; int nLocalDuals; VecGetLocalSize(data->tmp_vec_b, &nLocalB); VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals); PetscScalar *local_b, *local_duals; VecGetArray(data->tmp_vec_b, &local_b); VecGetArray(data->tmp_vec_duals0, &local_duals); for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++) local_duals[j] = local_b[i]; VecRestoreArray(data->tmp_vec_b, &local_b); VecRestoreArray(data->tmp_vec_duals0, &local_duals); // === K_DD === MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1); // === Prolongation from local dual nodes to B nodes. VecGetArray(data->tmp_vec_b, &local_b); VecGetArray(data->tmp_vec_duals1, &local_duals); for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++) local_b[i] = local_duals[j]; VecRestoreArray(data->tmp_vec_b, &local_b); VecRestoreArray(data->tmp_vec_duals0, &local_duals); // Multiply with scaled Lagrange constraint matrix. MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y); return 0; } PetscSolverFeti::PetscSolverFeti() : PetscSolver(), nComponents(-1), schurPrimalSolver(0) { FUNCNAME("PetscSolverFeti::PetscSolverFeti()"); string preconditionerName = ""; Parameters::get("parallel->solver->precon", preconditionerName); if (preconditionerName == "" || preconditionerName == "none") { MSG("Create FETI-DP solver with no preconditioner!\n"); fetiPreconditioner = FETI_NONE; } else if (preconditionerName == "dirichlet") { MSG("Create FETI-DP solver with Dirichlet preconditioner!\n"); fetiPreconditioner = FETI_DIRICHLET; } else if (preconditionerName == "lumped") { MSG("Create FETI-DP solver with lumped preconditioner!\n"); fetiPreconditioner = FETI_LUMPED; } else { ERROR_EXIT("Preconditioner \"%s\" not available!\n", preconditionerName.c_str()); } Parameters::get("parallel->feti->schur primal solver", schurPrimalSolver); TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1) ("Wrong solver \"%d\"for the Schur primal complement!\n", schurPrimalSolver); } void PetscSolverFeti::updateDofData() { FUNCNAME("PetscSolverFeti::updateDofData()"); 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. === /// Set of DOF indices that are considered to be primal variables. DofIndexSet primals; 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("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()); } 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 (globalPrimalIndex.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 (globalPrimalIndex.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 (globalPrimalIndex.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 (globalPrimalIndex.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 - globalPrimalIndex.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 (globalPrimalIndex.count(**it) == 0) { duals.insert(**it); globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals; nRankDuals++; } } int nOverallDuals = nRankDuals; mpi::globalAdd(nOverallDuals); MSG("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("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 (globalPrimalIndex.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 (globalPrimalIndex.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 (globalPrimalIndex.count(*(it->second[i])) == 0) dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++]; } } void PetscSolverFeti::createIndexB() { FUNCNAME("PetscSolverFeti::createIndeB()"); localIndexB.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 && globalPrimalIndex.count(i) == 0) if (duals.count(i) == 0 && globalPrimalIndex.count(i) == 0) localIndexB[i] = -1; // === Get correct index for all interior nodes. === nLocalInterior = 0; for (DofMapping::iterator it = localIndexB.begin(); it != localIndexB.end(); ++it) { it->second = nLocalInterior; nLocalInterior++; } nLocalDuals = duals.size(); TEST_EXIT_DBG(nLocalInterior + globalPrimalIndex.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) localIndexB[*it] = globalDualIndex[*it] - rStartB; } 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()"); if (schurPrimalSolver == 0) { MSG("Create iterative schur primal solver!\n"); schurPrimalData.mat_primal_primal = &mat_primal_primal; schurPrimalData.mat_primal_b = &mat_primal_b; schurPrimalData.mat_b_primal = &mat_b_primal; schurPrimalData.fetiSolver = this; VecCreateMPI(PETSC_COMM_WORLD, nRankB * nComponents, nOverallB * nComponents, &(schurPrimalData.tmp_vec_b)); VecCreateMPI(PETSC_COMM_WORLD, nRankPrimals * nComponents, nOverallPrimals * nComponents, &(schurPrimalData.tmp_vec_primal)); MatCreateShell(PETSC_COMM_WORLD, nRankPrimals * nComponents, nRankPrimals * nComponents, nOverallPrimals * nComponents, nOverallPrimals * nComponents, &schurPrimalData, &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); } else { MSG("Create direct schur primal solver!\n"); TEST_EXIT(ksp_b)("No solver object for local problem!\n"); double wtime = MPI::Wtime(); int nRowsRankPrimal = nRankPrimals * nComponents; int nRowsOverallPrimal = nOverallPrimals * nComponents; int nRowsRankB = nRankB * nComponents; int nRowsOverallB = nOverallB * nComponents; Mat matBPi; MatCreateMPIAIJ(PETSC_COMM_WORLD, nRowsRankB, nRowsRankPrimal, nRowsOverallB, nRowsOverallPrimal, 30, PETSC_NULL, 30, PETSC_NULL, &matBPi); Mat matPrimal; PetscInt nCols; const PetscInt *cols; const PetscScalar *values; int nLocalPrimals = globalPrimalIndex.size() * nComponents; map > > mat_b_primal_cols; for (int i = 0; i < nRankB * nComponents; i++) { PetscInt row = rStartB * nComponents + i; MatGetRow(mat_b_primal, row, &nCols, &cols, &values); for (int j = 0; j < nCols; j++) if (values[j] != 0.0) mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j])); MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values); } int maxLocalPrimal = mat_b_primal_cols.size(); mpi::globalMax(maxLocalPrimal); TEST_EXIT(mat_b_primal_cols.size() == globalPrimalIndex.size() * nComponents)("Should not happen!\n"); for (map > >::iterator it = mat_b_primal_cols.begin(); it != mat_b_primal_cols.end(); ++it) { Vec tmpVec; VecCreateSeq(PETSC_COMM_SELF, nRankB * nComponents, &tmpVec); for (unsigned int i = 0; i < it->second.size(); i++) VecSetValue(tmpVec, it->second[i].first, it->second[i].second, INSERT_VALUES); VecAssemblyBegin(tmpVec); VecAssemblyEnd(tmpVec); KSPSolve(ksp_b, tmpVec, tmpVec); PetscScalar *tmpValues; VecGetArray(tmpVec, &tmpValues); for (int i = 0; i < nRankB * nComponents; i++) MatSetValue(matBPi, rStartB * nComponents + i, it->first, tmpValues[i], ADD_VALUES); VecRestoreArray(tmpVec, &tmpValues); VecDestroy(&tmpVec); } MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY); MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal); MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN); MatDestroy(&matPrimal); MatDestroy(&matBPi); MatInfo minfo; MatGetInfo(mat_primal_primal, MAT_GLOBAL_SUM, &minfo); MSG("Schur primal matrix nnz = %f\n", minfo.nz_used); KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal); KSPSetOperators(ksp_schur_primal, mat_primal_primal, mat_primal_primal, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_"); KSPSetFromOptions(ksp_schur_primal); MSG("Creating Schur primal matrix needed %.5f seconds.\n", MPI::Wtime() - wtime); } } void PetscSolverFeti::destroySchurPrimalKsp() { FUNCNAME("PetscSolverFeti::destroySchurPrimal()"); if (schurPrimalSolver == 0) { schurPrimalData.mat_primal_primal = PETSC_NULL; schurPrimalData.mat_primal_b = PETSC_NULL; schurPrimalData.mat_b_primal = PETSC_NULL; schurPrimalData.fetiSolver = NULL; VecDestroy(&schurPrimalData.tmp_vec_b); VecDestroy(&schurPrimalData.tmp_vec_primal); MatDestroy(&mat_schur_primal); KSPDestroy(&ksp_schur_primal); } else { KSPDestroy(&ksp_schur_primal); } } void PetscSolverFeti::createFetiKsp() { FUNCNAME("PetscSolverFeti::createFetiKsp()"); // === Create FETI-DP solver object. === fetiData.mat_primal_b = &mat_primal_b; fetiData.mat_b_primal = &mat_b_primal; fetiData.mat_lagrange = &mat_lagrange; fetiData.fetiSolver = this; fetiData.ksp_schur_primal = &ksp_schur_primal; VecCreateMPI(PETSC_COMM_WORLD, nRankB * nComponents, nOverallB * nComponents, &(fetiData.tmp_vec_b)); VecCreateMPI(PETSC_COMM_WORLD, nRankLagrange * nComponents, nOverallLagrange * nComponents, &(fetiData.tmp_vec_lagrange)); VecCreateMPI(PETSC_COMM_WORLD, nRankPrimals * nComponents, nOverallPrimals * nComponents, &(fetiData.tmp_vec_primal)); MatCreateShell(PETSC_COMM_WORLD, nRankLagrange * nComponents, nRankLagrange * nComponents, nOverallLagrange * nComponents, nOverallLagrange * nComponents, &fetiData, &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); // === Create FETI-DP preconditioner object. === if (fetiPreconditioner != FETI_NONE) { MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled); MatScale(mat_lagrange_scaled, 0.5); } switch (fetiPreconditioner) { case FETI_DIRICHLET: KSPCreate(PETSC_COMM_SELF, &ksp_interior); KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_interior, "solver_interior_"); KSPSetFromOptions(ksp_interior); fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled; fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior; fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals; fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals; fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior; fetiDirichletPreconData.ksp_interior = &ksp_interior; VecCreateMPI(PETSC_COMM_WORLD, nRankB * nComponents, nOverallB * nComponents, &(fetiDirichletPreconData.tmp_vec_b)); MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals0)); MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals1)); MatGetVecs(mat_interior_interior, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_interior)); KSPGetPC(ksp_feti, &precon_feti); PCSetType(precon_feti, PCSHELL); PCShellSetContext(precon_feti, static_cast(&fetiDirichletPreconData)); PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon); break; case FETI_LUMPED: fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled; fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals; VecCreateMPI(PETSC_COMM_WORLD, nRankB * nComponents, nOverallB * nComponents, &(fetiLumpedPreconData.tmp_vec_b)); MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals0)); MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals1)); KSPGetPC(ksp_feti, &precon_feti); PCSetType(precon_feti, PCSHELL); PCShellSetContext(precon_feti, static_cast(&fetiLumpedPreconData)); PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon); break; default: break; } } void PetscSolverFeti::destroyFetiKsp() { FUNCNAME("PetscSolverFeti::destroyFetiKsp()"); // === Destroy FETI-DP solver object. === fetiData.mat_primal_b = PETSC_NULL; fetiData.mat_b_primal = PETSC_NULL; fetiData.mat_lagrange = PETSC_NULL; fetiData.fetiSolver = NULL; fetiData.ksp_schur_primal = PETSC_NULL; VecDestroy(&fetiData.tmp_vec_b); VecDestroy(&fetiData.tmp_vec_lagrange); VecDestroy(&fetiData.tmp_vec_primal); MatDestroy(&mat_feti); KSPDestroy(&ksp_feti); // === Destroy FETI-DP preconditioner object. === switch (fetiPreconditioner) { case FETI_DIRICHLET: KSPDestroy(&ksp_interior); fetiDirichletPreconData.mat_lagrange_scaled = NULL; fetiDirichletPreconData.mat_interior_interior = NULL; fetiDirichletPreconData.mat_duals_duals = NULL; fetiDirichletPreconData.mat_interior_duals = NULL; fetiDirichletPreconData.mat_duals_interior = NULL; fetiDirichletPreconData.ksp_interior = NULL; VecDestroy(&fetiDirichletPreconData.tmp_vec_b); VecDestroy(&fetiDirichletPreconData.tmp_vec_duals0); VecDestroy(&fetiDirichletPreconData.tmp_vec_duals1); VecDestroy(&fetiDirichletPreconData.tmp_vec_interior); MatDestroy(&mat_lagrange_scaled); break; case FETI_LUMPED: fetiLumpedPreconData.mat_lagrange_scaled = NULL; fetiLumpedPreconData.mat_duals_duals = NULL; VecDestroy(&fetiLumpedPreconData.tmp_vec_b); VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0); VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1); break; default: break; } } 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 = localIndexB.begin(); it != localIndexB.end(); ++it) { int petscIndex = it->second * 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) { FUNCNAME("PetscSolverFeti::fillPetscMatrix()"); nComponents = mat->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; int nRowsInterior = nLocalInterior * nComponents; int nRowsDual = nLocalDuals * nComponents; MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL, &mat_b_b); MatCreateMPIAIJ(PETSC_COMM_WORLD, nRowsRankPrimal, nRowsRankPrimal, nRowsOverallPrimal, nRowsOverallPrimal, 30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal); MatCreateMPIAIJ(PETSC_COMM_WORLD, nRowsRankB, nRowsRankPrimal, nRowsOverallB, nRowsOverallPrimal, 30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal); MatCreateMPIAIJ(PETSC_COMM_WORLD, nRowsRankPrimal, nRowsRankB, nRowsOverallPrimal, nRowsOverallB, 15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b); // === Create matrices for FETI-DP preconditioner. === if (fetiPreconditioner != FETI_NONE) MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsDual, nRowsDual, 30, PETSC_NULL, &mat_duals_duals); if (fetiPreconditioner == FETI_DIRICHLET) { MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsInterior, nRowsInterior, 30, PETSC_NULL, &mat_interior_interior); MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsInterior, nRowsDual, 30, PETSC_NULL, &mat_interior_duals); MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsDual, nRowsInterior, 30, PETSC_NULL, &mat_duals_interior); } // === 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); vector colsLocal, colsLocalOther; vector valuesLocal, valuesLocalOther; colsLocal.reserve(300); colsLocalOther.reserve(300); valuesLocal.reserve(300); valuesLocalOther.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 = globalPrimalIndex.count(*cursor) != 0; cols.clear(); colsOther.clear(); values.clear(); valuesOther.clear(); colsLocal.clear(); colsLocalOther.clear(); valuesLocal.clear(); valuesLocalOther.clear(); // Traverse all columns. for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { bool colPrimal = globalPrimalIndex.count(col(*icursor)) != 0; if (colPrimal) { // 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(localIndexB.count(col(*icursor))) ("No global B index for DOF %d!\n", col(*icursor)); int colIndex = (localIndexB[col(*icursor)] + rStartB) * nComponents + j; if (rowPrimal) { colsOther.push_back(colIndex); valuesOther.push_back(value(*icursor)); } else { cols.push_back(colIndex); values.push_back(value(*icursor)); } } // === For preconditioner === if (!rowPrimal && !colPrimal) { int rowIndex = localIndexB[*cursor]; int colIndex = localIndexB[col(*icursor)]; if (rowIndex < nLocalInterior) { if (colIndex < nLocalInterior) { int colIndex2 = localIndexB[col(*icursor)] * nComponents + j; colsLocal.push_back(colIndex2); valuesLocal.push_back(value(*icursor)); } else { int colIndex2 = (localIndexB[col(*icursor)] - nLocalInterior) * nComponents + j; colsLocalOther.push_back(colIndex2); valuesLocalOther.push_back(value(*icursor)); } } else { if (colIndex < nLocalInterior) { int colIndex2 = localIndexB[col(*icursor)] * nComponents + j; colsLocalOther.push_back(colIndex2); valuesLocalOther.push_back(value(*icursor)); } else { int colIndex2 = (localIndexB[col(*icursor)] - nLocalInterior) * nComponents + j; colsLocal.push_back(colIndex2); valuesLocal.push_back(value(*icursor)); } } } } // for each nnz in row 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(localIndexB.count(*cursor)) ("Should not happen!\n"); int localRowIndex = localIndexB[*cursor] * nComponents + i; for (unsigned int k = 0; k < cols.size(); k++) cols[k] -= rStartB * nComponents; MatSetValues(mat_b_b, 1, &localRowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); if (colsOther.size()) { int globalRowIndex = (localIndexB[*cursor] + rStartB) * nComponents + i; MatSetValues(mat_b_primal, 1, &globalRowIndex, colsOther.size(), &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); } } // === For preconditioner === if (!rowPrimal) { switch (fetiPreconditioner) { case FETI_DIRICHLET: { int rowIndex = localIndexB[*cursor]; if (rowIndex < nLocalInterior) { int rowIndex2 = localIndexB[*cursor] * nComponents + i; MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); if (colsLocalOther.size()) MatSetValues(mat_interior_duals, 1, &rowIndex2, colsLocalOther.size(), &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES); } else { int rowIndex2 = (localIndexB[*cursor] - nLocalInterior) * nComponents + i; MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); if (colsLocalOther.size()) MatSetValues(mat_duals_interior, 1, &rowIndex2, colsLocalOther.size(), &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES); } } break; case FETI_LUMPED: { int rowIndex = localIndexB[*cursor]; if (rowIndex >= nLocalInterior) { int rowIndex2 = (localIndexB[*cursor] - nLocalInterior) * nComponents + i; MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); } } break; default: break; } } } } } // === 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); // === Start global assembly procedure for preconditioner matrices. === if (fetiPreconditioner != FETI_NONE) { MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); } if (fetiPreconditioner == FETI_DIRICHLET) { MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY); MatAssemblyBegin(mat_interior_duals, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_interior_duals, MAT_FINAL_ASSEMBLY); MatAssemblyBegin(mat_duals_interior, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_duals_interior, MAT_FINAL_ASSEMBLY); } // === Create and fill PETSc matrix for Lagrange constraints. === createMatLagrange(); // === Create solver for the non primal (thus local) variables. === KSPCreate(PETSC_COMM_SELF, &ksp_b); KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_b, "solver_b_"); KSPSetFromOptions(ksp_b); // === Create PETSc solver for the Schur complement on primal variables. === createSchurPrimalKsp(); // === Create PETSc solver for the FETI-DP operator. === createFetiKsp(); } void PetscSolverFeti::fillPetscRhs(SystemVector *vec) { FUNCNAME("PetscSolverFeti::fillPetscRhs()"); int nComponents = vec->getSize(); VecCreateMPI(PETSC_COMM_WORLD, nRankB * nComponents, nOverallB * nComponents, &f_b); VecCreateMPI(PETSC_COMM_WORLD, nRankPrimals * nComponents, nOverallPrimals * nComponents, &f_primal); 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 (globalPrimalIndex.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(localIndexB.count(index)) ("Should not happen!\n"); index = (localIndexB[index] + rStartB) * nComponents + i; VecSetValue(f_b, index, *dofIt, INSERT_VALUES); } } } VecAssemblyBegin(f_b); VecAssemblyEnd(f_b); VecAssemblyBegin(f_primal); VecAssemblyEnd(f_primal); } void PetscSolverFeti::solveLocalProblem(Vec &rhs, Vec &sol) { FUNCNAME("PetscSolverFeti::solveLocalProblem()"); Vec tmp; VecCreateSeq(PETSC_COMM_SELF, nRankB * nComponents, &tmp); PetscScalar *tmpValues; VecGetArray(tmp, &tmpValues); PetscScalar *rhsValues; VecGetArray(rhs, &rhsValues); for (int i = 0; i < nRankB * nComponents; i++) tmpValues[i] = rhsValues[i]; VecRestoreArray(rhs, &rhsValues); VecRestoreArray(tmp, &tmpValues); KSPSolve(ksp_b, tmp, tmp); VecGetArray(tmp, &tmpValues); PetscScalar *solValues; VecGetArray(sol, &solValues); for (int i = 0; i < nRankB * nComponents; i++) solValues[i] = tmpValues[i]; VecRestoreArray(sol, &solValues); VecRestoreArray(tmp, &tmpValues); VecDestroy(&tmp); } void PetscSolverFeti::solveFetiMatrix(SystemVector &vec) { FUNCNAME("PetscSolverFeti::solveFetiMatrix()"); int nOverallNest = (nOverallB + nOverallPrimals + nOverallLagrange) * nComponents; Mat mat_lagrange_transpose; MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose); Mat A_global; MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, nOverallNest, nOverallNest, 30, PETSC_NULL, 30, PETSC_NULL, &A_global); Vec A_global_rhs, A_global_sol; MatGetVecs(A_global, &A_global_rhs, &A_global_sol); PetscInt nCols; const PetscInt *cols; const PetscScalar *vals; // === mat_b_b === for (int i = 0; i < nRankB * nComponents; i++) { int rowIndex = rStartB * nComponents + i; MatGetRow(mat_b_b, i, &nCols, &cols, &vals); MatSetValues(A_global, 1, &rowIndex, nCols, cols, vals, ADD_VALUES); MatRestoreRow(mat_b_b, rowIndex, &nCols, &cols, &vals); } PetscScalar *g_f_b; VecGetArray(f_b, &g_f_b); for (int i = 0; i < nRankB * nComponents; i++) VecSetValue(A_global_rhs, rStartB * nComponents + i, g_f_b[i], INSERT_VALUES); VecRestoreArray(f_b, &g_f_b); // === mat_primal_primal === for (int i = 0; i < nRankPrimals * nComponents; i++) { int rowIndex = rStartPrimals * nComponents + i; MatGetRow(mat_primal_primal, rowIndex, &nCols, &cols, &vals); int rowIndexA = nOverallB * nComponents + rowIndex; vector colsA(nCols); for (int j = 0; j < nCols; j++) colsA[j] = nOverallB * nComponents + cols[j]; MatSetValues(A_global, 1, &rowIndexA, nCols, &(colsA[0]), vals, ADD_VALUES); MatRestoreRow(mat_primal_primal, rowIndex, &nCols, &cols, &vals); } PetscScalar *g_f_primal; VecGetArray(f_primal, &g_f_primal); for (int i = 0; i < nRankPrimals * nComponents; i++) VecSetValue(A_global_rhs, (nOverallB + rStartPrimals) * nComponents + i, g_f_primal[i], INSERT_VALUES); VecRestoreArray(f_primal, &g_f_primal); // === mat_b_primal === for (int i = 0; i < nRankB * nComponents; i++) { int rowIndex = rStartB * nComponents + i; MatGetRow(mat_b_primal, rowIndex, &nCols, &cols, &vals); vector colsA(nCols); for (int j = 0; j < nCols; j++) colsA[j] = nOverallB * nComponents + cols[j]; MatSetValues(A_global, 1, &rowIndex, nCols, &(colsA[0]), vals, ADD_VALUES); MatRestoreRow(mat_b_primal, rowIndex, &nCols, &cols, &vals); } // === mat_primal_b === for (int i = 0; i < nRankPrimals * nComponents; i++) { int rowIndex = rStartPrimals * nComponents + i; MatGetRow(mat_primal_b, rowIndex, &nCols, &cols, &vals); int rowIndexA = nOverallB * nComponents + rowIndex; MatSetValues(A_global, 1, &rowIndexA, nCols, cols, vals, ADD_VALUES); MatRestoreRow(mat_primal_b, rowIndex, &nCols, &cols, &vals); } // === mat_lagrange === for (int i = 0; i < nRankLagrange * nComponents; i++) { int rowIndex = rStartLagrange * nComponents + i; MatGetRow(mat_lagrange, rowIndex, &nCols, &cols, &vals); int rowIndexA = (nOverallB + nOverallPrimals) * nComponents + rowIndex; MatSetValues(A_global, 1, &rowIndexA, nCols, cols, vals, ADD_VALUES); MatRestoreRow(mat_lagrange, rowIndex, &nCols, &cols, &vals); } // === mat_lagrange_transpose === for (int i = 0; i < nRankB * nComponents; i++) { int rowIndex = rStartB * nComponents + i; MatGetRow(mat_lagrange_transpose, rowIndex, &nCols, &cols, &vals); vector colsA(nCols); for (int j = 0; j < nCols; j++) colsA[j] = (nOverallB + nOverallPrimals) * nComponents + cols[j]; MatSetValues(A_global, 1, &rowIndex, nCols, &(colsA[0]), vals, ADD_VALUES); MatRestoreRow(mat_lagrange_transpose, rowIndex, &nCols, &cols, &vals); } MatAssemblyBegin(A_global, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A_global, MAT_FINAL_ASSEMBLY); VecAssemblyBegin(A_global_rhs); VecAssemblyEnd(A_global_rhs); // === Create solver and solve the overall FETI system. === KSP ksp; KSPCreate(PETSC_COMM_WORLD, &ksp); KSPSetOperators(ksp, A_global, A_global, SAME_NONZERO_PATTERN); KSPSetFromOptions(ksp); KSPSolve(ksp, A_global_rhs, A_global_sol); Vec u_b, u_primal; VecDuplicate(f_b, &u_b); VecDuplicate(f_primal, &u_primal); vector localBIndex, globalBIndex; localBIndex.reserve(nRankB * nComponents); globalBIndex.reserve(nRankB * nComponents); for (int i = 0; i < nRankB * nComponents; i++) { localBIndex.push_back(rStartB * nComponents + i); globalBIndex.push_back(rStartB * nComponents + i); } IS localBIs, globalBIs; ISCreateGeneral(PETSC_COMM_SELF, localBIndex.size(), &(localBIndex[0]), PETSC_USE_POINTER, &localBIs); ISCreateGeneral(PETSC_COMM_SELF, globalBIndex.size(), &(globalBIndex[0]), PETSC_USE_POINTER, &globalBIs); VecScatter vecscat; VecScatterCreate(A_global_sol, globalBIs, u_b, localBIs, &vecscat); VecScatterBegin(vecscat, A_global_sol, u_b, INSERT_VALUES, SCATTER_FORWARD); VecScatterEnd(vecscat, A_global_sol, u_b, INSERT_VALUES, SCATTER_FORWARD); ISDestroy(&localBIs); ISDestroy(&globalBIs); VecScatterDestroy(&vecscat); localBIndex.resize(0); globalBIndex.resize(0); localBIndex.reserve(nRankPrimals * nComponents); globalBIndex.reserve(nRankPrimals * nComponents); for (int i = 0; i < nRankPrimals * nComponents; i++) { localBIndex.push_back(rStartPrimals * nComponents + i); globalBIndex.push_back(nOverallB * nComponents + rStartPrimals * nComponents + i); } ISCreateGeneral(PETSC_COMM_SELF, localBIndex.size(), &(localBIndex[0]), PETSC_USE_POINTER, &localBIs); ISCreateGeneral(PETSC_COMM_SELF, globalBIndex.size(), &(globalBIndex[0]), PETSC_USE_POINTER, &globalBIs); VecScatterCreate(A_global_sol, globalBIs, u_primal, localBIs, &vecscat); VecScatterBegin(vecscat, A_global_sol, u_primal, INSERT_VALUES, SCATTER_FORWARD); VecScatterEnd(vecscat, A_global_sol, u_primal, INSERT_VALUES, SCATTER_FORWARD); ISDestroy(&localBIs); ISDestroy(&globalBIs); VecScatterDestroy(&vecscat); recoverSolution(u_b, u_primal, vec); KSPDestroy(&ksp); } void PetscSolverFeti::solveReducedFetiMatrix(SystemVector &vec) { FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()"); // RHS 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); VecDuplicate(f_b, &tmp_b0); VecDuplicate(f_b, &tmp_b1); MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal0); MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal1); // === Create new rhs === // d = L inv(K_BB) f_B - L inv(K_BB) K_BPi inv(S_PiPi) [f_Pi - K_PiB inv(K_BB) f_B] // vec_rhs = L * inv(K_BB) * f_B solveLocalProblem(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); solveLocalProblem(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); solveLocalProblem(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); solveLocalProblem(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); solveLocalProblem(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); VecDestroy(&f_b); VecDestroy(&f_primal); } void PetscSolverFeti::destroyMatrixData() { FUNCNAME("PetscSolverFeti::destroyMatrixData()"); MatDestroy(&mat_b_b); MatDestroy(&mat_primal_primal); MatDestroy(&mat_b_primal); MatDestroy(&mat_primal_b); MatDestroy(&mat_lagrange); // === Destroy preconditioner data structures. === if (fetiPreconditioner != FETI_NONE) MatDestroy(&mat_duals_duals); if (fetiPreconditioner == FETI_DIRICHLET) { MatDestroy(&mat_interior_interior); MatDestroy(&mat_interior_duals); MatDestroy(&mat_duals_interior); } destroySchurPrimalKsp(); destroyFetiKsp(); KSPDestroy(&ksp_b); } void PetscSolverFeti::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo) { FUNCNAME("PetscSolverFeti::solvePetscMatrix()"); int debug = 0; Parameters::get("parallel->debug feti", debug); if (debug) { solveFetiMatrix(vec); } else { solveReducedFetiMatrix(vec); } MeshDistributor::globalMeshDistributor->synchVector(vec); } }