Skip to content

GitLab

Commit 0e6691d8 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

FETI-DP Dirichlet preconditioner.

parent a2ef2a34
Hide whitespace changes
Inline Side-by-side
Showing with 286 additions and 11 deletions
AMDiS/src/parallel/PetscSolverFeti.cc
View file @ 0e6691d8
......@@ -70,6 +70,62 @@ namespace AMDiS {
}
// y = PC * x
PetscErrorCode petscApplyFetiPrecon(PC pc, Vec x, Vec y)
{
void *ctx;
PCShellGetContext(pc, &ctx);
PetscFetiPreconData* data = static_cast<PetscFetiPreconData*>(ctx);
MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
int sizeB;
int sizeBound;
VecGetLocalSize(data->tmp_vec_b, &sizeB);
VecGetLocalSize(data->tmp_vec_bound0, &sizeBound);
PetscScalar *local_b;
VecGetArray(data->tmp_vec_b, &local_b);
PetscScalar *local_bound;
VecGetArray(data->tmp_vec_bound0, &local_bound);
for (int i = sizeB - sizeBound, j = 0; i < sizeB; i++, j++)
local_bound[j] = local_b[i];
VecRestoreArray(data->tmp_vec_b, &local_b);
VecRestoreArray(data->tmp_vec_bound0, &local_bound);
MatMult(*(data->mat_bound_bound), data->tmp_vec_bound0, data->tmp_vec_bound1);
MatMult(*(data->mat_interior_bound), data->tmp_vec_bound0, data->tmp_vec_interior);
KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
MatMult(*(data->mat_bound_interior), data->tmp_vec_interior, data->tmp_vec_bound0);
VecAXPBY(data->tmp_vec_bound0, 1.0, -1.0, data->tmp_vec_bound1);
VecGetArray(data->tmp_vec_b, &local_b);
VecGetArray(data->tmp_vec_bound0, &local_bound);
for (int i = sizeB - sizeBound, j = 0; i < sizeB; i++, j++)
local_b[i] = local_bound[j];
VecRestoreArray(data->tmp_vec_b, &local_b);
VecRestoreArray(data->tmp_vec_bound0, &local_bound);
MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);
return 0;
}
void PetscSolverFeti::updateDofData()
{
FUNCNAME("PetscSolverFeti::updateDofData()");
......@@ -379,14 +435,15 @@ namespace AMDiS {
// === Get correct index for all interior nodes. ===
int nInterior = 0;
nLocalInterior = 0;
for (DofMapping::iterator it = globalIndexB.begin();
it != globalIndexB.end(); ++it) {
it->second = nInterior + rStartB;
nInterior++;
it->second = nLocalInterior + rStartB;
nLocalInterior++;
}
nLocalBound = duals.size();
TEST_EXIT_DBG(nInterior + primals.size() + duals.size() ==
TEST_EXIT_DBG(nLocalInterior + primals.size() + duals.size() ==
static_cast<unsigned int>(admin->getUsedDofs()))
("Should not happen!\n");
......@@ -505,6 +562,8 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createFetiKsp()");
// === Create FETI-DP solver object. ===
petscFetiData.mat_primal_primal = &mat_primal_primal;
petscFetiData.mat_primal_b = &mat_primal_b;
petscFetiData.mat_b_primal = &mat_b_primal;
......@@ -529,6 +588,38 @@ namespace AMDiS {
KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
KSPSetOptionsPrefix(ksp_feti, "solver_feti_");
KSPSetFromOptions(ksp_feti);
// === Create FETI-DP Dirichlet preconditioner object. ===
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);
MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
MatScale(mat_lagrange_scaled, 0.5);
petscFetiPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
petscFetiPreconData.mat_interior_interior = &mat_interior_interior;
petscFetiPreconData.mat_bound_bound = &mat_bound_bound;
petscFetiPreconData.mat_interior_bound = &mat_interior_bound;
petscFetiPreconData.mat_bound_interior = &mat_bound_interior;
petscFetiPreconData.ksp_interior = &ksp_interior;
petscFetiPreconData.nInterior = nRankB - duals.size();
VecDuplicate(f_b, &(petscFetiPreconData.tmp_vec_b));
MatGetVecs(mat_bound_bound, PETSC_NULL, &(petscFetiPreconData.tmp_vec_bound0));
MatGetVecs(mat_bound_bound, PETSC_NULL, &(petscFetiPreconData.tmp_vec_bound1));
MatGetVecs(mat_interior_interior, PETSC_NULL, &(petscFetiPreconData.tmp_vec_interior));
KSPGetPC(ksp_feti, &precon_feti);
PCSetType(precon_feti, PCSHELL);
PCShellSetContext(precon_feti, static_cast<void*>(&petscFetiPreconData));
PCShellSetApply(precon_feti, petscApplyFetiPrecon);
}
......@@ -536,6 +627,8 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::destroyFetiKsp()");
// === Destroy FETI-DP solver object. ===
petscFetiData.mat_primal_primal = PETSC_NULL;
petscFetiData.mat_primal_b = PETSC_NULL;
petscFetiData.mat_b_primal = PETSC_NULL;
......@@ -549,6 +642,24 @@ namespace AMDiS {
MatDestroy(mat_feti);
KSPDestroy(ksp_feti);
// === Destroy FETI-DP Dirichlet preconditioner object. ===
KSPDestroy(ksp_interior);
petscFetiPreconData.mat_lagrange_scaled = NULL;
petscFetiPreconData.mat_interior_interior = NULL;
petscFetiPreconData.mat_bound_bound = NULL;
petscFetiPreconData.mat_interior_bound = NULL;
petscFetiPreconData.mat_bound_interior = NULL;
petscFetiPreconData.ksp_interior = NULL;
VecDestroy(petscFetiPreconData.tmp_vec_b);
VecDestroy(petscFetiPreconData.tmp_vec_bound0);
VecDestroy(petscFetiPreconData.tmp_vec_bound1);
VecDestroy(petscFetiPreconData.tmp_vec_interior);
MatDestroy(mat_lagrange_scaled);
}
......@@ -658,6 +769,8 @@ namespace AMDiS {
int nRowsOverallB = nOverallB * nComponents;
int nRowsRankPrimal = nRankPrimals * nComponents;
int nRowsOverallPrimal = nOverallPrimals * nComponents;
int nRowsInterior = nLocalInterior * nComponents;
int nRowsBound = nLocalBound * nComponents;
MatCreateMPIAIJ(PETSC_COMM_WORLD,
nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB,
......@@ -678,6 +791,25 @@ namespace AMDiS {
nRowsOverallPrimal, nRowsOverallB,
100, PETSC_NULL, 100, PETSC_NULL, &mat_primal_b);
// === Create matrices for Dirichlet FETI-DP preconditioner. ===
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsInterior, 100, PETSC_NULL,
&mat_interior_interior);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsBound, nRowsBound, 100, PETSC_NULL,
&mat_bound_bound);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsBound, 100, PETSC_NULL,
&mat_interior_bound);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsBound, nRowsInterior, 100, PETSC_NULL,
&mat_bound_interior);
// === Prepare traverse of sequentially created matrices. ===
......@@ -695,6 +827,13 @@ namespace AMDiS {
values.reserve(300);
valuesOther.reserve(300);
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. ===
......@@ -710,19 +849,27 @@ namespace AMDiS {
// 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 = primals.count(*cursor) != 0;
cols.clear();
values.clear();
colsOther.clear();
values.clear();
valuesOther.clear();
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 (primals.count(col(*icursor)) != 0) {
bool colPrimal = primals.count(col(*icursor)) != 0;
if (colPrimal) {
// Column is a primal variable.
TEST_EXIT_DBG(globalPrimalIndex.count(col(*icursor)))
......@@ -753,6 +900,47 @@ namespace AMDiS {
values.push_back(value(*icursor));
}
}
// === For preconditioner ===
if (!rowPrimal && !colPrimal) {
int rowIndex = globalIndexB[*cursor] - rStartB;
int colIndex = globalIndexB[col(*icursor)] - rStartB;
if (rowIndex < nLocalInterior) {
if (colIndex < nLocalInterior) {
int colIndex2 =
(globalIndexB[col(*icursor)] - rStartB) * nComponents + j;
colsLocal.push_back(colIndex2);
valuesLocal.push_back(value(*icursor));
} else {
int colIndex2 =
(globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;
colsLocalOther.push_back(colIndex2);
valuesLocalOther.push_back(value(*icursor));
}
} else {
if (colIndex < nLocalInterior) {
int colIndex2 =
(globalIndexB[col(*icursor)] - rStartB) * nComponents + j;
colsLocalOther.push_back(colIndex2);
valuesLocalOther.push_back(value(*icursor));
} else {
int colIndex2 =
(globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;
colsLocal.push_back(colIndex2);
valuesLocal.push_back(value(*icursor));
}
}
}
}
if (rowPrimal) {
......@@ -778,6 +966,37 @@ namespace AMDiS {
MatSetValues(mat_b_primal, 1, &rowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
}
// === For preconditioner ===
if (!rowPrimal) {
int rowIndex = globalIndexB[*cursor] - rStartB;
if (rowIndex < nLocalInterior) {
int rowIndex2 =
(globalIndexB[*cursor] - rStartB) * nComponents + i;
MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_interior_bound, 1, &rowIndex2, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
} else {
int rowIndex2 =
(globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
MatSetValues(mat_bound_bound, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_bound_interior, 1, &rowIndex2, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
}
}
}
}
}
......@@ -798,6 +1017,19 @@ namespace AMDiS {
MatAssemblyEnd(mat_primal_b, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_bound_bound, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_bound_bound, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_interior_bound, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_interior_bound, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_bound_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_bound_interior, MAT_FINAL_ASSEMBLY);
// === Create and fill PETSc's right hand side vectors. ===
VecCreate(PETSC_COMM_WORLD, &f_b);
......@@ -1100,7 +1332,15 @@ namespace AMDiS {
destroySchurPrimalKsp();
destroyFetiKsp();
destroyFetiKsp();
// === Destroy preconditioner data structures. ===
MatDestroy(mat_interior_interior);
MatDestroy(mat_bound_bound);
MatDestroy(mat_interior_bound);
MatDestroy(mat_bound_interior);
}
......@@ -1117,9 +1357,9 @@ namespace AMDiS {
solveFetiMatrix(vec);
} else {
solveReducedFetiMatrix(vec);
}
}
}
#endif
}
AMDiS/src/parallel/PetscSolverFeti.h
View file @ 0e6691d8
......@@ -87,6 +87,26 @@ namespace AMDiS {
};
struct PetscFetiPreconData {
/// Matrix of scaled Lagrange variables.
Mat *mat_lagrange_scaled;
Mat *mat_interior_interior, *mat_bound_bound, *mat_interior_bound, *mat_bound_interior;
/// Pointer to the solver for \ref PetscSolverFeti::mat_bb.
KSP *ksp_interior;
/// Temporal vector on the B variables.
Vec tmp_vec_b;
Vec tmp_vec_bound0, tmp_vec_bound1;
Vec tmp_vec_interior;
int nInterior;
};
/** \brief
* FETI-DP implementation based on PETSc.
......@@ -267,6 +287,21 @@ namespace AMDiS {
/// Data for MatMult operation in matrix \ref mat_feti
PetscFetiData petscFetiData;
Mat mat_lagrange_scaled;
PC precon_feti;
PetscFetiPreconData petscFetiPreconData;
Mat mat_interior_interior, mat_bound_bound, mat_interior_bound, mat_bound_interior;
KSP ksp_interior;
int nLocalInterior;
int nLocalBound;
};
#endif
......
Markdown is supported
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