// // 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 "AMDiS.h" #include "parallel/PetscSolverFeti.h" #include "parallel/PetscSolverFetiStructs.h" #include "parallel/StdMpi.h" #include "parallel/MpiHelper.h" #include "parallel/SubDomainSolver.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->subSolver->getMatIntCoarse(), x, data->tmp_vec_b); data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b); MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal); MatMult(data->subSolver->getMatCoarseCoarse(), 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->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b); MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange); MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal); KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal); MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b); data->subSolver->solveGlobal(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 (map::iterator it = data->localToDualMap.begin(); it != data->localToDualMap.end(); ++it) local_duals[it->second] = local_b[it->first]; 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 (map::iterator it = data->localToDualMap.begin(); it != data->localToDualMap.end(); ++it) local_b[it->first] = local_duals[it->second]; 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 (map::iterator it = data->localToDualMap.begin(); it != data->localToDualMap.end(); ++it) local_duals[it->second] = local_b[it->first]; 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 (map::iterator it = data->localToDualMap.begin(); it != data->localToDualMap.end(); ++it) local_b[it->first] = local_duals[it->second]; 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(), schurPrimalSolver(0), multiLevelTest(false), subDomainSolver(NULL), meshLevel(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); Parameters::get("parallel->multi level test", multiLevelTest); if (multiLevelTest) meshLevel = 1; } void PetscSolverFeti::initialize(vector feSpaces) { FUNCNAME("PetscSolverFeti::initialize()"); TEST_EXIT_DBG(meshLevel + 1 == meshDistributor->getMeshLevelData().getLevelNumber()) ("Mesh hierarchy does not contain %d levels!\n", meshLevel + 1); MeshLevelData& levelData = meshDistributor->getMeshLevelData(); if (subDomainSolver == NULL) { if (meshLevel == 0) { subDomainSolver = new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm); } else { subDomainSolver = new SubDomainSolver(meshDistributor, levelData.getMpiComm(meshLevel - 1), levelData.getMpiComm(meshLevel)); subDomainSolver->setLevel(meshLevel); } } primalDofMap.init(levelData, feSpaces, meshDistributor->getFeSpaces(), true, true); primalDofMap.setComputeMatIndex(true); dualDofMap.init(levelData, feSpaces, meshDistributor->getFeSpaces(), false, false); dualDofMap.setComputeMatIndex(true); lagrangeMap.init(levelData, feSpaces, meshDistributor->getFeSpaces(), true, true); lagrangeMap.setComputeMatIndex(true); if (meshLevel == 0) { localDofMap.init(levelData, feSpaces, meshDistributor->getFeSpaces(), false, false); localDofMap.setComputeMatIndex(true); } else { localDofMap.init(levelData, feSpaces, meshDistributor->getFeSpaces(), true, true); localDofMap.setComputeMatIndex(true); } if (fetiPreconditioner == FETI_DIRICHLET) { TEST_EXIT(meshLevel == 0) ("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n"); interiorDofMap.init(levelData, feSpaces, meshDistributor->getFeSpaces(), false, false); interiorDofMap.setComputeMatIndex(true); } } void PetscSolverFeti::createFetiData() { FUNCNAME("PetscSolverFeti::createFetiData()"); TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n"); TEST_EXIT(meshDistributor->getFeSpaces().size() > 0) ("No FE space defined in mesh distributor!\n"); MeshLevelData& levelData = meshDistributor->getMeshLevelData(); primalDofMap.clear(); dualDofMap.clear(); lagrangeMap.clear(); localDofMap.clear(); if (fetiPreconditioner == FETI_DIRICHLET) interiorDofMap.clear(); primalDofMap.setDofComm(meshDistributor->getDofComm()); lagrangeMap.setDofComm(meshDistributor->getDofComm()); primalDofMap.setMpiComm(levelData.getMpiComm(0), 0); dualDofMap.setMpiComm(levelData.getMpiComm(0), 0); lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0); localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel); if (meshLevel > 0) localDofMap.setDofComm(meshDistributor->getDofCommSd()); for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); createPrimals(feSpace); createDuals(feSpace); createIndexB(feSpace); } primalDofMap.update(); dualDofMap.update(); localDofMap.update(); if (fetiPreconditioner != FETI_NONE) interiorDofMap.update(); for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); createLagrange(feSpace); } lagrangeMap.update(); MSG("FETI-DP data created on mesh level %d\n", meshLevel); for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); MSG("FETI-DP data for %d-ith FE space:\n", i); MSG(" nRankPrimals = %d nOverallPrimals = %d\n", primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs); MSG(" nRankDuals = %d nOverallDuals = %d\n", dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs); MSG(" nRankLagrange = %d nOverallLagrange = %d\n", lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs); TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == static_cast(feSpace->getAdmin()->getUsedDofs())) ("Should not happen!\n"); } // If multi level test, inform sub domain solver about coarse space. subDomainSolver->setDofMapping(&primalDofMap, &localDofMap); } void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace) { 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. DofContainerSet& vertices = meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX]; DofIndexSet primals; 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. === for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) if (meshDistributor->getDofMap()[feSpace].isRankDof(*it)) primalDofMap[feSpace].insertRankDof(*it); else primalDofMap[feSpace].insertNonRankDof(*it); } void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createDuals()"); // === Create global index of the dual nodes on each rank. === DofContainer allBoundaryDofs; meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs); for (DofContainer::iterator it = allBoundaryDofs.begin(); it != allBoundaryDofs.end(); ++it) if (!isPrimal(feSpace, **it)) dualDofMap[feSpace].insertRankDof(**it); } void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createLagrange()"); boundaryDofRanks[feSpace].clear(); if (dualDofMap[feSpace].nLocalDofs == 0) return; // === Create for each dual node that is owned by the rank, the set === // === of ranks that contain this node (denoted by W(x_j)). === for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace); !it.end(); it.nextRank()) for (; !it.endDofIter(); it.nextDof()) { if (!isPrimal(feSpace, it.getDofIndex())) { boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank); boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank()); } } // === Communicate these sets for all rank owned dual nodes to other === // === ranks that also have this node. === StdMpi > > stdMpi(meshDistributor->getMpiComm()); for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace); !it.end(); it.nextRank()) for (; !it.endDofIter(); it.nextDof()) if (!isPrimal(feSpace, it.getDofIndex())) stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]); stdMpi.updateSendDataSize(); for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); !it.end(); it.nextRank()) { bool recvFromRank = false; for (; !it.endDofIter(); it.nextDof()) { if (!isPrimal(feSpace, it.getDofIndex())) { recvFromRank = true; break; } } if (recvFromRank) stdMpi.recv(it.getRank()); } stdMpi.startCommunication(); for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); !it.end(); it.nextRank()) { int i = 0; for (; !it.endDofIter(); it.nextDof()) if (!isPrimal(feSpace, it.getDofIndex())) boundaryDofRanks[feSpace][it.getDofIndex()] = stdMpi.getRecvData(it.getRank())[i++]; } // === Reserve for each dual node, on the rank that owns this node, the === // === appropriate number of Lagrange constraints. === int nRankLagrange = 0; DofMap& dualMap = dualDofMap[feSpace].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) { lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange); int degree = boundaryDofRanks[feSpace][it->first].size(); nRankLagrange += (degree * (degree - 1)) / 2; } else { lagrangeMap[feSpace].insertNonRankDof(it->first); } } lagrangeMap[feSpace].nRankDofs = nRankLagrange; } void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace) { FUNCNAME("PetscSolverFeti::createIndexB()"); DOFAdmin* admin = feSpace->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. === int nLocalInterior = 0; for (int i = 0; i < admin->getUsedSize(); i++) { if (admin->isDofFree(i) == false && isPrimal(feSpace, i) == false && isDual(feSpace, i) == false) { if (meshLevel == 0) { localDofMap[feSpace].insertRankDof(i, nLocalInterior); if (fetiPreconditioner != FETI_NONE) interiorDofMap[feSpace].insertRankDof(i, nLocalInterior); nLocalInterior++; } else { if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i)) localDofMap[feSpace].insertRankDof(i); else localDofMap[feSpace].insertNonRankDof(i); } } } // === And finally, add the global indicies of all dual nodes. === for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin(); it != dualDofMap[feSpace].getMap().end(); ++it) if (meshLevel == 0) localDofMap[feSpace].insertRankDof(it->first); else { if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first)) localDofMap[feSpace].insertRankDof(it->first); else localDofMap[feSpace].insertNonRankDof(it->first); } } void PetscSolverFeti::createMatLagrange(vector &feSpaces) { FUNCNAME("PetscSolverFeti::createMatLagrange()"); // === Create distributed matrix for Lagrange constraints. === MatCreateMPIAIJ(mpiComm, lagrangeMap.getRankDofs(), localDofMap.getRankDofs(), lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(), 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 (unsigned int k = 0; k < feSpaces.size(); k++) { DofMap &dualMap = dualDofMap[feSpaces[k]].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first)) ("Should not happen!\n"); // Global index of the first Lagrange constriant for this node. int index = lagrangeMap.getMatIndex(k, it->first); // Copy set of all ranks that contain this dual node. vector W(boundaryDofRanks[feSpaces[k]][it->first].begin(), boundaryDofRanks[feSpaces[k]][it->first].end()); // Number of ranks that contain this dual node. int degree = W.size(); for (int i = 0; i < degree; i++) { for (int j = i + 1; j < degree; j++) { if (W[i] == mpiRank || W[j] == mpiRank) { int colIndex = localDofMap.getMatIndex(k, it->first); double value = (W[i] == mpiRank ? 1.0 : -1.0); MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES); } index++; } } } } MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY); } void PetscSolverFeti::createSchurPrimalKsp(vector &feSpaces) { FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()"); if (schurPrimalSolver == 0) { MSG("Create iterative schur primal solver!\n"); schurPrimalData.subSolver = subDomainSolver; VecCreateMPI(mpiComm, localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &(schurPrimalData.tmp_vec_b)); VecCreateMPI(mpiComm, primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(), &(schurPrimalData.tmp_vec_primal)); MatCreateShell(mpiComm, primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(), &schurPrimalData, &mat_schur_primal); MatShellSetOperation(mat_schur_primal, MATOP_MULT, (void(*)(void))petscMultMatSchurPrimal); KSPCreate(mpiComm, &ksp_schur_primal); KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_"); KSPSetType(ksp_schur_primal, KSPGMRES); KSPSetFromOptions(ksp_schur_primal); } else { MSG("Create direct schur primal solver!\n"); double wtime = MPI::Wtime(); int nRowsRankPrimal = primalDofMap.getRankDofs(); int nRowsOverallPrimal = primalDofMap.getOverallDofs(); int nRowsRankB = localDofMap.getRankDofs(); int nRowsOverallB = localDofMap.getOverallDofs(); Mat matBPi; MatCreateMPIAIJ(mpiComm, nRowsRankB, nRowsRankPrimal, nRowsOverallB, nRowsOverallPrimal, 30, PETSC_NULL, 30, PETSC_NULL, &matBPi); Mat matPrimal; PetscInt nCols; const PetscInt *cols; const PetscScalar *values; map > > mat_b_primal_cols; for (int i = 0; i < nRowsRankB; i++) { PetscInt row = localDofMap.getStartDofs() + i; MatGetRow(subDomainSolver->getMatIntCoarse(), 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(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values); } TEST_EXIT(static_cast(mat_b_primal_cols.size()) == primalDofMap.getLocalDofs()) ("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, nRowsRankB, &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); subDomainSolver->solve(tmpVec, tmpVec); PetscScalar *tmpValues; VecGetArray(tmpVec, &tmpValues); for (int i = 0; i < nRowsRankB; i++) MatSetValue(matBPi, localDofMap.getStartDofs() + i, it->first, tmpValues[i], ADD_VALUES); VecRestoreArray(tmpVec, &tmpValues); VecDestroy(&tmpVec); } MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY); MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal); MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal); MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN); MatDestroy(&matPrimal); MatDestroy(&matBPi); MatInfo minfo; MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo); MSG("Schur primal matrix nnz = %f\n", minfo.nz_used); KSPCreate(mpiComm, &ksp_schur_primal); KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_"); KSPSetType(ksp_schur_primal, KSPPREONLY); PC pc_schur_primal; KSPGetPC(ksp_schur_primal, &pc_schur_primal); PCSetType(pc_schur_primal, PCLU); PCFactorSetMatSolverPackage(pc_schur_primal, MATSOLVERMUMPS); 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.subSolver = NULL; VecDestroy(&schurPrimalData.tmp_vec_b); VecDestroy(&schurPrimalData.tmp_vec_primal); } MatDestroy(&mat_schur_primal); KSPDestroy(&ksp_schur_primal); } void PetscSolverFeti::createFetiKsp(vector &feSpaces) { FUNCNAME("PetscSolverFeti::createFetiKsp()"); // === Create FETI-DP solver object. === fetiData.mat_lagrange = &mat_lagrange; fetiData.subSolver = subDomainSolver; fetiData.ksp_schur_primal = &ksp_schur_primal; VecCreateMPI(mpiComm, localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &(fetiData.tmp_vec_b)); VecCreateMPI(mpiComm, lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(), &(fetiData.tmp_vec_lagrange)); VecCreateMPI(mpiComm, primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(), &(fetiData.tmp_vec_primal)); MatCreateShell(mpiComm, lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(), &fetiData, &mat_feti); MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti); KSPCreate(mpiComm, &ksp_feti); KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(ksp_feti, "feti_"); KSPSetType(ksp_feti, KSPGMRES); KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000); 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, "precon_interior_"); KSPSetType(ksp_interior, KSPPREONLY); PC pc_interior; KSPGetPC(ksp_interior, &pc_interior); PCSetType(pc_interior, PCLU); PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK); 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(mpiComm, localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &(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)); for (unsigned int i = 0; i < feSpaces.size(); i++) { DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { DegreeOfFreedom d = it->first; int matIndexLocal = localDofMap.getLocalMatIndex(i, d); int matIndexDual = dualDofMap.getLocalMatIndex(i, d); fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual; } } 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; for (unsigned int i = 0; i < feSpaces.size(); i++) { DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(); for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { DegreeOfFreedom d = it->first; int matIndexLocal = localDofMap.getLocalMatIndex(i, d); int matIndexDual = dualDofMap.getLocalMatIndex(i, d); fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual; } } VecCreateMPI(mpiComm, localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &(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_lagrange = PETSC_NULL; fetiData.subSolver = 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 feSpaces = getFeSpaces(&vec); vector globalIsIndex, localIsIndex; globalIsIndex.reserve(primalDofMap.getLocalDofs()); localIsIndex.reserve(primalDofMap.getLocalDofs()); { int cnt = 0; for (unsigned int i = 0; i < feSpaces.size(); i++) { DofMap& dofMap = primalDofMap[feSpaces[i]].getMap(); for (DofMap::iterator it = dofMap.begin();it != dofMap.end(); ++it) { globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first)); localIsIndex.push_back(cnt++); } } TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs()) ("Should not happen!\n"); } 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. === int cnt = 0; for (int i = 0; i < vec.getSize(); i++) { DOFVector& dofVec = *(vec.getDOFVector(i)); for (DofMap::iterator it = localDofMap[feSpaces[i]].getMap().begin(); it != localDofMap[feSpaces[i]].getMap().end(); ++it) { int petscIndex = localDofMap.getLocalMatIndex(i, it->first); dofVec[it->first] = localSolB[petscIndex]; } for (DofMap::iterator it = primalDofMap[feSpaces[i]].getMap().begin(); it != primalDofMap[feSpaces[i]].getMap().end(); ++it) dofVec[it->first] = localSolPrimal[cnt++]; } VecRestoreArray(vec_sol_b, &localSolB); VecRestoreArray(local_sol_primal, &localSolPrimal); VecDestroy(&local_sol_primal); } void PetscSolverFeti::fillPetscMatrix(Matrix *mat) { FUNCNAME("PetscSolverFeti::fillPetscMatrix()"); // === Create all sets and indices. === vector feSpaces = getFeSpaces(mat); initialize(feSpaces); createFetiData(); // === Create matrices for the FETI-DP method. === subDomainSolver->fillPetscMatrix(mat); // === Create matrices for FETI-DP preconditioner. === if (fetiPreconditioner != FETI_NONE) { int nRowsDual = dualDofMap.getRankDofs(); MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsDual, nRowsDual, 60, PETSC_NULL, &mat_duals_duals); if (fetiPreconditioner == FETI_DIRICHLET) { int nRowsInterior = interiorDofMap.getRankDofs(); MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsInterior, nRowsInterior, 60, PETSC_NULL, &mat_interior_interior); MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsInterior, nRowsDual, 60, PETSC_NULL, &mat_interior_duals); MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsDual, nRowsInterior, 60, 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 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. === int nComponents = mat->getSize(); 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 = isPrimal(feSpaces[i], *cursor); 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 = isPrimal(feSpaces[j], col(*icursor)); if (!rowPrimal && !colPrimal) { if (!isDual(feSpaces[i], *cursor)) { if (!isDual(feSpaces[j], col(*icursor))) { int colIndex = interiorDofMap.getLocalMatIndex(j, col(*icursor)); colsLocal.push_back(colIndex); valuesLocal.push_back(value(*icursor)); } else { int colIndex = dualDofMap.getLocalMatIndex(j, col(*icursor)); colsLocalOther.push_back(colIndex); valuesLocalOther.push_back(value(*icursor)); } } else { if (!isDual(feSpaces[j], col(*icursor))) { int colIndex = interiorDofMap.getLocalMatIndex(j, col(*icursor)); colsLocalOther.push_back(colIndex); valuesLocalOther.push_back(value(*icursor)); } else { int colIndex = dualDofMap.getLocalMatIndex(j, col(*icursor)); colsLocal.push_back(colIndex); valuesLocal.push_back(value(*icursor)); } } } } // for each nnz in row // === Set matrix values for preconditioner === if (!rowPrimal) { switch (fetiPreconditioner) { case FETI_DIRICHLET: if (!isDual(feSpaces[i], *cursor)) { int rowIndex = interiorDofMap.getLocalMatIndex(i, *cursor); MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); if (colsLocalOther.size()) MatSetValues(mat_interior_duals, 1, &rowIndex, colsLocalOther.size(), &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES); } else { int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor); MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); if (colsLocalOther.size()) MatSetValues(mat_duals_interior, 1, &rowIndex, colsLocalOther.size(), &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES); } break; case FETI_LUMPED: if (isDual(feSpaces[i], *cursor)) { int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor); MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); } break; default: break; } } } } } // === 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(feSpaces); // === Create PETSc solver for the Schur complement on primal variables. === createSchurPrimalKsp(feSpaces); // === Create PETSc solver for the FETI-DP operator. === createFetiKsp(feSpaces); } void PetscSolverFeti::fillPetscRhs(SystemVector *vec) { FUNCNAME("PetscSolverFeti::fillPetscRhs()"); subDomainSolver->fillPetscRhs(vec); } void PetscSolverFeti::solveFetiMatrix(SystemVector &vec) { FUNCNAME("PetscSolverFeti::solveFetiMatrix()"); // The function was removed when FETI-DP was reformulated for mixed finite // elements. To reconstruct this function, check for revision number 2024. ERROR_EXIT("This function was removed!\n"); } 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; VecCreateMPI(mpiComm, localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &tmp_b0); VecCreateMPI(mpiComm, localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &tmp_b1); VecCreateMPI(mpiComm, primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(), &tmp_primal0); VecCreateMPI(mpiComm, primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(), &tmp_primal1); MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange0); MatGetVecs(mat_lagrange, PETSC_NULL, &vec_rhs); // === 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 subDomainSolver->solveGlobal(subDomainSolver->getRhsInterior(), tmp_b0); MatMult(mat_lagrange, tmp_b0, vec_rhs); // tmp_primal0 = M_PiB * inv(K_BB) * f_B MatMult(subDomainSolver->getMatCoarseInt(), tmp_b0, tmp_primal0); // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_B VecAXPBY(tmp_primal0, 1.0, -1.0, subDomainSolver->getRhsCoarseSpace()); // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_B) KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0); // MatMult(subDomainSolver->getMatIntCoarse(), tmp_primal0, tmp_b0); subDomainSolver->solveGlobal(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(subDomainSolver->getRhsCoarseSpace(), tmp_primal0); subDomainSolver->solveGlobal(subDomainSolver->getRhsInterior(), tmp_b0); MatMult(subDomainSolver->getMatCoarseInt(), tmp_b0, tmp_primal1); VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1); MatMultTranspose(mat_lagrange, vec_rhs, tmp_b0); subDomainSolver->solveGlobal(tmp_b0, tmp_b0); MatMult(subDomainSolver->getMatCoarseInt(), 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(subDomainSolver->getRhsInterior(), tmp_b0); MatMultTranspose(mat_lagrange, vec_rhs, tmp_b1); VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1); MatMult(subDomainSolver->getMatIntCoarse(), tmp_primal0, tmp_b1); VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1); subDomainSolver->solveGlobal(tmp_b0, tmp_b0); // === And recover AMDiS solution vectors. === recoverSolution(tmp_b0, tmp_primal0, vec); VecDestroy(&vec_rhs); VecDestroy(&tmp_b0); VecDestroy(&tmp_b1); VecDestroy(&tmp_lagrange0); VecDestroy(&tmp_primal0); VecDestroy(&tmp_primal1); } void PetscSolverFeti::destroyMatrixData() { FUNCNAME("PetscSolverFeti::destroyMatrixData()"); 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(); } void PetscSolverFeti::destroyVectorData() { FUNCNAME("PetscSolverFeti::destroyVectorData()"); subDomainSolver->destroyVectorData(); } 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); } }