Skip to content

GitLab

Commit 82a26536 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Work on lumped preconditioner for FETI-DP Stokes.

parent 739cafa2
Hide whitespace changes
Inline Side-by-side
AMDiS/src/parallel/PetscSolverFeti.cc
View file @ 82a26536
......@@ -108,7 +108,7 @@ namespace AMDiS {
if (fullInterface != NULL)
interfaceDofMap.init(levelData, feSpaces, uniqueFe);
if (fetiPreconditioner != FETI_NONE) {
if (fetiPreconditioner == FETI_DIRICHLET) {
TEST_EXIT(meshLevel == 0)
("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");
......@@ -133,7 +133,7 @@ namespace AMDiS {
dualDofMap.clear();
lagrangeMap.clear();
localDofMap.clear();
if (fetiPreconditioner != FETI_NONE)
if (fetiPreconditioner == FETI_DIRICHLET)
interiorDofMap.clear();
primalDofMap.setDofComm(meshDistributor->getDofComm());
......@@ -143,7 +143,7 @@ namespace AMDiS {
dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
if (fetiPreconditioner != FETI_NONE)
if (fetiPreconditioner == FETI_DIRICHLET)
interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
if (meshLevel == 0)
......@@ -172,7 +172,7 @@ namespace AMDiS {
primalDofMap.update();
dualDofMap.update();
localDofMap.update();
if (fetiPreconditioner != FETI_NONE)
if (fetiPreconditioner == FETI_DIRICHLET)
interiorDofMap.update();
if (fullInterface != NULL)
......@@ -504,7 +504,7 @@ namespace AMDiS {
if (meshLevel == 0) {
localDofMap[feSpace].insertRankDof(i, nLocalInterior);
if (fetiPreconditioner != FETI_NONE)
if (fetiPreconditioner == FETI_DIRICHLET)
interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
nLocalInterior++;
......@@ -820,8 +820,8 @@ namespace AMDiS {
MatNullSpace matNullspace;
MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullspace);
// MatSetNullSpace(mat_feti, matNullspace);
// KSPSetNullSpace(ksp_feti, matNullspace);
MatSetNullSpace(mat_feti, matNullspace);
KSPSetNullSpace(ksp_feti, matNullspace);
}
......@@ -834,10 +834,13 @@ namespace AMDiS {
switch (fetiPreconditioner) {
case FETI_DIRICHLET:
TEST_EXIT_DBG(meshLevel == 0)
case FETI_DIRICHLET:
TEST_EXIT(meshLevel == 0)
("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");
TEST_EXIT(!enableStokesMode)
("Stokes mode does not yet support the Dirichlet precondition!\n");
KSPCreate(PETSC_COMM_SELF, &ksp_interior);
KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior,
SAME_NONZERO_PATTERN);
......@@ -898,10 +901,14 @@ namespace AMDiS {
break;
case FETI_LUMPED:
fetiLumpedPreconData.enableStokesMode = enableStokesMode;
fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;
for (unsigned int i = 0; i < feSpaces.size(); i++) {
if (enableStokesMode && feSpaces[i] == fullInterface)
continue;
DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
DegreeOfFreedom d = it->first;
......@@ -1347,182 +1354,9 @@ namespace AMDiS {
MPI::COMM_WORLD.Barrier();
MSG("FETI-DP timing 02: %.5f seconds (creation of interior matrices)\n",
MPI::Wtime() - wtime);
}
// === Create matrices for FETI-DP preconditioner. ===
if (fetiPreconditioner != FETI_NONE) {
wtime = MPI::Wtime();
int nRowsDual = dualDofMap.getRankDofs();
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsDual, nRowsDual, 100, PETSC_NULL,
&mat_duals_duals);
MatSetOption(mat_duals_duals, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
if (fetiPreconditioner == FETI_DIRICHLET) {
int nRowsInterior = interiorDofMap.getRankDofs();
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsInterior, 100, PETSC_NULL,
&mat_interior_interior);
MatSetOption(mat_interior_interior,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsDual, 100, PETSC_NULL,
&mat_interior_duals);
MatSetOption(mat_interior_duals,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsDual, nRowsInterior, 100, PETSC_NULL,
&mat_duals_interior);
MatSetOption(mat_duals_interior,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
}
// === 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<row, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
vector<int> colsLocal, colsLocalOther;
vector<double> 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<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
// Traverse all rows.
for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()),
cend = end<row>((*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<nz>(cursor), icend = end<nz>(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);
}
if (printTimings) {
MPI::COMM_WORLD.Barrier();
MSG("FETI-DP timing 03: %.5f seconds (creation of preconditioner matrices)\n",
MPI::Wtime() - wtime);
}
}
// For the case, that we want to print the timings, we force the LU
// factorization of the local problems to be done here explicitly.
if (printTimings) {
// For the case, that we want to print the timings, we force the LU
// factorization of the local problems to be done here explicitly.
wtime = MPI::Wtime();
KSPSetUp(subdomain->getSolver());
KSPSetUpOnBlocks(subdomain->getSolver());
......@@ -1530,25 +1364,21 @@ namespace AMDiS {
MSG("FETI-DP timing 04: %.5f seconds (factorization of subdomain matrices)\n",
MPI::Wtime() - wtime);
}
// === Create matrices for FETI-DP preconditioner. ===
createPreconditionerMatrix(mat);
// === Create and fill PETSc matrix for Lagrange constraints. ===
createMatLagrange(feSpaces);
// === If required, run debug tests. ===
dbgMatrix();
// === Create PETSc solver for the Schur complement on primal variables. ===
// === Create PETSc solver for the Schur complement on primal variables. ===
createSchurPrimalKsp(feSpaces);
// === Create PETSc solver for the FETI-DP operator. ===
createFetiKsp(feSpaces);
}
......@@ -1560,7 +1390,193 @@ namespace AMDiS {
subdomain->fillPetscRhs(vec);
}
void PetscSolverFeti::createPreconditionerMatrix(Matrix<DOFMatrix*> *mat)
{
FUNCNAME("PetscSolverFeti::createPreconditionerMatrix()");
if (fetiPreconditioner != FETI_NONE)
return;
double wtime = MPI::Wtime();
vector<const FiniteElemSpace*> feSpaces = AMDiS::getComponentFeSpaces(*mat);
int nRowsDual = dualDofMap.getRankDofs();
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsDual, nRowsDual, 100, PETSC_NULL,
&mat_duals_duals);
MatSetOption(mat_duals_duals, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
if (fetiPreconditioner == FETI_DIRICHLET) {
int nRowsInterior = interiorDofMap.getRankDofs();
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsInterior, 100, PETSC_NULL,
&mat_interior_interior);
MatSetOption(mat_interior_interior,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsDual, 100, PETSC_NULL,
&mat_interior_duals);
MatSetOption(mat_interior_duals,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsDual, nRowsInterior, 100, PETSC_NULL,
&mat_duals_interior);
MatSetOption(mat_duals_interior,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
}
// === 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<row, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
vector<int> colsLocal, colsLocalOther;
vector<double> 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 rowComponent = 0; rowComponent < nComponents; rowComponent++) {
for (int colComponent = 0; colComponent < nComponents; colComponent++) {
DOFMatrix* dofMat = (*mat)[rowComponent][colComponent];
if (!dofMat)
continue;
const FiniteElemSpace *rowFeSpace = feSpaces[rowComponent];
const FiniteElemSpace *colFeSpace = feSpaces[colComponent];
traits::col<Matrix>::type col(dofMat->getBaseMatrix());
traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
// Traverse all rows.
for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()),
cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
if (isPrimal(rowFeSpace, *cursor))
continue;
colsLocal.clear();
colsLocalOther.clear();
valuesLocal.clear();
valuesLocalOther.clear();
// Traverse all columns.
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
if (isPrimal(colFeSpace, col(*icursor)))
continue;
if (!isDual(rowFeSpace, *cursor)) {
if (!isDual(colFeSpace, col(*icursor))) {
int colIndex =
interiorDofMap.getLocalMatIndex(colComponent, col(*icursor));
colsLocal.push_back(colIndex);
valuesLocal.push_back(value(*icursor));
} else {
int colIndex =
dualDofMap.getLocalMatIndex(colComponent, col(*icursor));
colsLocalOther.push_back(colIndex);
valuesLocalOther.push_back(value(*icursor));
}
} else {
if (!isDual(colFeSpace, col(*icursor))) {
int colIndex =
interiorDofMap.getLocalMatIndex(colComponent, col(*icursor));
colsLocalOther.push_back(colIndex);
valuesLocalOther.push_back(value(*icursor));
} else {
int colIndex =
dualDofMap.getLocalMatIndex(colComponent, col(*icursor));
colsLocal.push_back(colIndex);
valuesLocal.push_back(value(*icursor));
}
}
} // for each nnz in row
// === Set matrix values for preconditioner ===
switch (fetiPreconditioner) {
case FETI_DIRICHLET:
if (!isDual(rowFeSpace, *cursor)) {
int rowIndex =
interiorDofMap.getLocalMatIndex(rowComponent, *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(rowComponent, *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(rowFeSpace, *cursor)) {
int rowIndex = dualDofMap.getLocalMatIndex(rowComponent, *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_LUMPED) {
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_duals_duals, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_duals_duals, 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);
}
if (printTimings) {
MPI::COMM_WORLD.Barrier();
MSG("FETI-DP timing 03: %.5f seconds (creation of preconditioner matrices)\n",
MPI::Wtime() - wtime);
}
}
void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
{
FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
......
AMDiS/src/parallel/PetscSolverFeti.h
View file @ 82a26536
......@@ -125,6 +125,9 @@ namespace AMDiS {
/// to \ref mat_lagrange.
void createMatLagrange(vector<const FiniteElemSpace*> &feSpaces);
///
void createPreconditionerMatrix(Matrix<DOFMatrix*> *mat);
/// Creates PETSc KSP solver object for solving the Schur complement
/// system on the primal variables, \ref ksp_schur_primal
void createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces);
......
AMDiS/src/parallel/PetscSolverFetiOperators.cc
View file @ 82a26536
......@@ -238,8 +238,25 @@ namespace AMDiS {
PCShellGetContext(pc, &ctx);
FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);
Vec xvec, yvec;
if (data->enableStokesMode) {
Vec x_interface, x_lagrange, y_interface, y_lagrange;
VecNestGetSubVec(x, 0, &x_interface);
VecNestGetSubVec(x, 1, &x_lagrange);
VecNestGetSubVec(y, 0, &y_interface);
VecNestGetSubVec(y, 1, &y_lagrange);
xvec = x_lagrange;
yvec = y_lagrange;
VecCopy(x_interface, y_interface);
} else {
xvec = x;
yvec = y;
}
// Multiply with scaled Lagrange constraint matrix.
MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
MatMultTranspose(*(data->mat_lagrange_scaled), xvec, data->tmp_vec_b);
// === Restriction of the B nodes to the boundary nodes. ===
......@@ -280,7 +297,7 @@ namespace AMDiS {
// Multiply with scaled Lagrange constraint matrix.
MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);
MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, yvec);
return 0;
}
......
AMDiS/src/parallel/PetscSolverFetiStructs.h
View file @ 82a26536
......@@ -99,6 +99,10 @@ namespace AMDiS {
struct FetiLumpedPreconData {
/// Defines whether the precondition is applied for a Stokes like FETI-DP
/// problem with interface variables which handled in a different way.
bool enableStokesMode;
/// Matrix of scaled Lagrange variables.
Mat *mat_lagrange_scaled;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment