From 8f20263c39a46237f3fa004ce6ec97568184c64f Mon Sep 17 00:00:00 2001 From: Thomas Witkowski Date: Tue, 15 May 2012 12:56:40 +0000 Subject: [PATCH] Replaces subdomain solver in FETI-DP to be a general PetscSolver. --- AMDiS/src/parallel/PetscSolverFeti.cc | 3 +- AMDiS/src/parallel/PetscSolverFeti.h | 2 +- AMDiS/src/parallel/PetscSolverFetiStructs.h | 6 +- AMDiS/src/parallel/PetscSolverGlobalMatrix.cc | 284 ++++++++++++++++++ AMDiS/src/parallel/PetscSolverGlobalMatrix.h | 6 + 5 files changed, 296 insertions(+), 5 deletions(-) diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index 01faf258..87dff5fd 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -16,6 +16,7 @@ #include "parallel/StdMpi.h" #include "parallel/MpiHelper.h" #include "parallel/SubDomainSolver.h" +#include "parallel/PetscSolverGlobalMatrix.h" #include "io/VtkWriter.h" namespace AMDiS { @@ -233,7 +234,7 @@ namespace AMDiS { MeshLevelData& levelData = meshDistributor->getMeshLevelData(); if (subDomainSolver == NULL) { - subDomainSolver = new SubDomainSolver(); + subDomainSolver = new PetscSolverGlobalMatrix(); if (meshLevel == 0) { subDomainSolver->setMeshDistributor(meshDistributor, diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h index 947e0cec..de7cb328 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.h +++ b/AMDiS/src/parallel/PetscSolverFeti.h @@ -259,7 +259,7 @@ namespace AMDiS { bool multiLevelTest; - SubDomainSolver *subDomainSolver; + PetscSolver *subDomainSolver; int meshLevel; diff --git a/AMDiS/src/parallel/PetscSolverFetiStructs.h b/AMDiS/src/parallel/PetscSolverFetiStructs.h index 0458dc16..e65eb3a4 100644 --- a/AMDiS/src/parallel/PetscSolverFetiStructs.h +++ b/AMDiS/src/parallel/PetscSolverFetiStructs.h @@ -24,7 +24,7 @@ #define AMDIS_PETSC_SOLVER_FETI_STRUCTS_H #include -#include "parallel/SubDomainSolver.h" +#include "parallel/PetscSolver.h" namespace AMDiS { @@ -43,7 +43,7 @@ namespace AMDiS { /// Temporal vecor in the primal variables. Vec tmp_vec_primal; - SubDomainSolver* subSolver; + PetscSolver* subSolver; }; @@ -65,7 +65,7 @@ namespace AMDiS { /// Temporal vector on the lagrange variables. Vec tmp_vec_lagrange; - SubDomainSolver* subSolver; + PetscSolver* subSolver; /// Pointer to the solver of the schur complement on the primal variables. KSP *ksp_schur_primal; diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc index dd7b0e68..b7817352 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc @@ -21,6 +21,11 @@ namespace AMDiS { { FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()"); + if (coarseSpaceMap != NULL) { + fillPetscMatrixWithCoarseSpace(mat); + return; + } + TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n"); TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n"); TEST_EXIT_DBG(mat)("No DOF matrix defined!\n"); @@ -124,10 +129,210 @@ namespace AMDiS { } + void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix *mat) + { + FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()"); + + vector feSpaces = getFeSpaces(mat); + + int nRowsRankInterior = interiorMap->getRankDofs(); + int nRowsOverallInterior = interiorMap->getOverallDofs(); + + if (subdomainLevel == 0) { + nGlobalOverallInterior = nRowsOverallInterior; + + MatCreateSeqAIJ(mpiCommLocal, nRowsRankInterior, nRowsRankInterior, + 60, PETSC_NULL, &matIntInt); + } else { + MatCreateMPIAIJ(mpiCommLocal, + nRowsRankInterior, nRowsRankInterior, + nRowsOverallInterior, nRowsOverallInterior, + 60, PETSC_NULL, 60, PETSC_NULL, &matIntInt); + } + + if (coarseSpaceMap) { + int nRowsRankCoarse = coarseSpaceMap->getRankDofs(); + int nRowsOverallCoarse = coarseSpaceMap->getOverallDofs(); + + MatCreateMPIAIJ(mpiCommGlobal, + nRowsRankCoarse, nRowsRankCoarse, + nRowsOverallCoarse, nRowsOverallCoarse, + 60, PETSC_NULL, 60, PETSC_NULL, &matCoarseCoarse); + + MatCreateMPIAIJ(mpiCommGlobal, + nRowsRankCoarse, nRowsRankInterior, + nRowsOverallCoarse, nGlobalOverallInterior, + 60, PETSC_NULL, 60, PETSC_NULL, &matCoarseInt); + + MatCreateMPIAIJ(mpiCommGlobal, + nRowsRankInterior, nRowsRankCoarse, + nGlobalOverallInterior, nRowsOverallCoarse, + 60, PETSC_NULL, 60, PETSC_NULL, &matIntCoarse); + } + + // === 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. === + + 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 = isCoarseSpace(feSpaces[i], *cursor); + + cols.clear(); + colsOther.clear(); + values.clear(); + valuesOther.clear(); + + // Traverse all columns. + for (icursor_type icursor = begin(cursor), icend = end(cursor); + icursor != icend; ++icursor) { + + bool colPrimal = isCoarseSpace(feSpaces[j], col(*icursor)); + + if (colPrimal) { + if (rowPrimal) { + cols.push_back(col(*icursor)); + values.push_back(value(*icursor)); + } else { + colsOther.push_back(col(*icursor)); + valuesOther.push_back(value(*icursor)); + } + } else { + if (rowPrimal) { + colsOther.push_back(col(*icursor)); + valuesOther.push_back(value(*icursor)); + } else { + cols.push_back(col(*icursor)); + values.push_back(value(*icursor)); + } + } + } // for each nnz in row + + + // === Set matrix values. === + + if (rowPrimal) { + int rowIndex = coarseSpaceMap->getMatIndex(i, *cursor); + for (unsigned int k = 0; k < cols.size(); k++) + cols[k] = coarseSpaceMap->getMatIndex(j, cols[k]); + + MatSetValues(matCoarseCoarse, 1, &rowIndex, cols.size(), + &(cols[0]), &(values[0]), ADD_VALUES); + + if (colsOther.size()) { + if (subdomainLevel == 0) { + for (unsigned int k = 0; k < colsOther.size(); k++) + colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]); + } else { + for (unsigned int k = 0; k < colsOther.size(); k++) + colsOther[k] = + interiorMap->getMatIndex(j, colsOther[k]) + rStartInterior; + } + + MatSetValues(matCoarseInt, 1, &rowIndex, colsOther.size(), + &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); + } + } else { + int localRowIndex = + (subdomainLevel == 0 ? interiorMap->getLocalMatIndex(i, *cursor) : + interiorMap->getMatIndex(i, *cursor)); + + for (unsigned int k = 0; k < cols.size(); k++) { + if (subdomainLevel == 0) + cols[k] = interiorMap->getLocalMatIndex(j, cols[k]); + else + cols[k] = interiorMap->getMatIndex(j, cols[k]); + } + + MatSetValues(matIntInt, 1, &localRowIndex, cols.size(), + &(cols[0]), &(values[0]), ADD_VALUES); + + if (colsOther.size()) { + int globalRowIndex = interiorMap->getMatIndex(i, *cursor); + + if (subdomainLevel != 0) + globalRowIndex += rStartInterior; + + for (unsigned int k = 0; k < colsOther.size(); k++) + colsOther[k] = coarseSpaceMap->getMatIndex(j, colsOther[k]); + + MatSetValues(matIntCoarse, 1, &globalRowIndex, colsOther.size(), + &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); + } + } + } + } + } + + // === Start global assembly procedure. === + + MatAssemblyBegin(matIntInt, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(matIntInt, MAT_FINAL_ASSEMBLY); + + if (coarseSpaceMap) { + MatAssemblyBegin(matCoarseCoarse, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(matCoarseCoarse, MAT_FINAL_ASSEMBLY); + + MatAssemblyBegin(matIntCoarse, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(matIntCoarse, MAT_FINAL_ASSEMBLY); + + MatAssemblyBegin(matCoarseInt, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(matCoarseInt, MAT_FINAL_ASSEMBLY); + } + + + // === Create solver for the non primal (thus local) variables. === + + KSPCreate(mpiCommLocal, &kspInterior); + KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN); + KSPSetOptionsPrefix(kspInterior, "interior_"); + KSPSetType(kspInterior, KSPPREONLY); + PC pcInterior; + KSPGetPC(kspInterior, &pcInterior); + PCSetType(pcInterior, PCLU); + if (subdomainLevel == 0) + PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK); + else + PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS); + KSPSetFromOptions(kspInterior); + } + + void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec) { FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()"); + if (coarseSpaceMap) { + fillPetscRhsWithCoarseSpace(vec); + return; + } + TEST_EXIT_DBG(vec)("No DOF vector defined!\n"); TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n"); @@ -159,6 +364,47 @@ namespace AMDiS { } + void PetscSolverGlobalMatrix::fillPetscRhsWithCoarseSpace(SystemVector *vec) + { + FUNCNAME("SubDomainSolver::fillPetscRhs()"); + + VecCreateMPI(mpiCommGlobal, + interiorMap->getRankDofs(), + nGlobalOverallInterior, + &rhsInterior); + + if (coarseSpaceMap) + VecCreateMPI(mpiCommGlobal, + coarseSpaceMap->getRankDofs(), + coarseSpaceMap->getOverallDofs(), + &rhsCoarseSpace); + + + for (int i = 0; i < vec->getSize(); i++) { + const FiniteElemSpace *feSpace = vec->getDOFVector(i)->getFeSpace(); + DOFVector::Iterator dofIt(vec->getDOFVector(i), USED_DOFS); + for (dofIt.reset(); !dofIt.end(); ++dofIt) { + int index = dofIt.getDOFIndex(); + if (isCoarseSpace(feSpace, index)) { + index = coarseSpaceMap->getMatIndex(i, index); + VecSetValue(rhsCoarseSpace, index, *dofIt, ADD_VALUES); + } else { + index = interiorMap->getMatIndex(i, index) + rStartInterior; + VecSetValue(rhsInterior, index, *dofIt, ADD_VALUES); + } + } + } + + VecAssemblyBegin(rhsInterior); + VecAssemblyEnd(rhsInterior); + + if (coarseSpaceMap) { + VecAssemblyBegin(rhsCoarseSpace); + VecAssemblyEnd(rhsCoarseSpace); + } + } + + void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo) { @@ -202,6 +448,44 @@ namespace AMDiS { } + void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol) + { + FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()"); + + Vec tmp; + if (mpiCommLocal.Get_size() == 1) + VecCreateSeq(mpiCommLocal, interiorMap->getRankDofs(), &tmp); + else + VecCreateMPI(mpiCommLocal, + interiorMap->getRankDofs(), + interiorMap->getOverallDofs(), + &tmp); + + PetscScalar *tmpValues, *rhsValues; + VecGetArray(tmp, &tmpValues); + VecGetArray(rhs, &rhsValues); + + for (int i = 0; i < interiorMap->getRankDofs(); i++) + tmpValues[i] = rhsValues[i]; + + VecRestoreArray(rhs, &rhsValues); + VecRestoreArray(tmp, &tmpValues); + + KSPSolve(kspInterior, tmp, tmp); + + VecGetArray(tmp, &tmpValues); + VecGetArray(sol, &rhsValues); + + for (int i = 0; i < interiorMap->getRankDofs(); i++) + rhsValues[i] = tmpValues[i]; + + VecRestoreArray(sol, &rhsValues); + VecRestoreArray(tmp, &tmpValues); + + VecDestroy(&tmp); + } + + void PetscSolverGlobalMatrix::destroyMatrixData() { FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()"); diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h index 653e94fc..c04abe1e 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h @@ -49,10 +49,16 @@ namespace AMDiS { void fillPetscMatrix(Matrix *mat); + void fillPetscMatrixWithCoarseSpace(Matrix *mat); + void fillPetscRhs(SystemVector *vec); + void fillPetscRhsWithCoarseSpace(SystemVector *vec); + void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo); + void solveGlobal(Vec &rhs, Vec &sol); + void destroyMatrixData(); void destroyVectorData(); -- GitLab