diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 1423195786b8bb2ca4dd72bb98ff926ac1649870..59ef2bb442bdfcab3c2570a3e35a3ebaaedd36eb 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -768,9 +768,9 @@ namespace AMDiS {
fetiData.subSolver = subdomain;
fetiData.ksp_schur_primal = &ksp_schur_primal;
- localDofMap.createVec(fetiData.tmp_vec_b, nGlobalOverallInterior);
+ localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
- primalDofMap.createVec(fetiData.tmp_vec_primal);
+ primalDofMap.createVec(fetiData.tmp_vec_primal0);
if (enableStokesMode == false) {
MatCreateShell(mpiCommGlobal,
@@ -782,6 +782,10 @@ namespace AMDiS {
MatShellSetOperation(mat_feti, MATOP_MULT,
(void(*)(void))petscMultMatFeti);
} else {
+ localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
+ primalDofMap.createVec(fetiData.tmp_vec_primal1);
+ interfaceDofMap.createVec(fetiData.tmp_vec_interface);
+
MatCreateShell(mpiCommGlobal,
interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
@@ -789,7 +793,7 @@ namespace AMDiS {
interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
&fetiData, &mat_feti);
MatShellSetOperation(mat_feti, MATOP_MULT,
- (void(*)(void))petscMultMatFetiInterface);
+ (void(*)(void))petscMultMatFetiInterface);
}
KSPCreate(mpiCommGlobal, &ksp_feti);
@@ -799,6 +803,27 @@ namespace AMDiS {
KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
KSPSetFromOptions(ksp_feti);
+ if (enableStokesMode) {
+ Vec nullSpaceInterface;
+ Vec nullSpaceLagrange;
+ Vec nullSpaceBasis;
+
+ interfaceDofMap.createVec(nullSpaceInterface);
+ VecSet(nullSpaceInterface, 1.0);
+
+ lagrangeMap.createVec(nullSpaceLagrange);
+ VecSet(nullSpaceLagrange, 0.0);
+
+ Vec vecArray[2] = {nullSpaceInterface, nullSpaceLagrange};
+ VecCreateNest(mpiCommGlobal, 2, PETSC_NULL, vecArray, &nullSpaceBasis);
+
+
+ MatNullSpace matNullspace;
+ MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullspace);
+ MatSetNullSpace(mat_feti, matNullspace);
+ KSPSetNullSpace(ksp_feti, matNullspace);
+ }
+
// === Create FETI-DP preconditioner object. ===
@@ -913,9 +938,16 @@ namespace AMDiS {
fetiData.subSolver = NULL;
fetiData.ksp_schur_primal = PETSC_NULL;
- VecDestroy(&fetiData.tmp_vec_b);
+ VecDestroy(&fetiData.tmp_vec_b0);
VecDestroy(&fetiData.tmp_vec_lagrange);
- VecDestroy(&fetiData.tmp_vec_primal);
+ VecDestroy(&fetiData.tmp_vec_primal0);
+
+ if (enableStokesMode) {
+ VecDestroy(&fetiData.tmp_vec_b1);
+ VecDestroy(&fetiData.tmp_vec_primal1);
+ VecDestroy(&fetiData.tmp_vec_interface);
+ }
+
MatDestroy(&mat_feti);
KSPDestroy(&ksp_feti);
@@ -1520,30 +1552,39 @@ namespace AMDiS {
FetiTimings::fetiPreconditioner);
}
+
// === Solve for u_primals. ===
- VecCopy(subdomain->getVecRhsCoarse(), tmp_primal0);
- subdomain->solveGlobal(subdomain->getVecRhsInterior(), tmp_b0);
- MatMult(subdomain->getMatCoarseInterior(), tmp_b0, tmp_primal1);
- VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1);
- MatMultTranspose(mat_lagrange, vecSol, tmp_b0);
- subdomain->solveGlobal(tmp_b0, tmp_b0);
- MatMult(subdomain->getMatCoarseInterior(), tmp_b0, tmp_primal1);
+ if (!enableStokesMode) {
+ MatMultTranspose(mat_lagrange, vecSol, tmp_b0);
+ VecAYPX(tmp_b0, -1.0, subdomain->getVecRhsInterior());
+ } else {
+ MatMultTranspose(mat_lagrange, vecSolLagrange, tmp_b0);
+ VecAYPX(tmp_b0, -1.0, subdomain->getVecRhsInterior());
+
+ MatMult(subdomain->getMatInteriorCoarse(1), vecSolInterface, tmp_b1);
+ VecAXPY(tmp_b0, -1.0, tmp_b1);
+ }
+
+ subdomain->solveGlobal(tmp_b0, tmp_b1);
+ MatMult(subdomain->getMatCoarseInterior(), tmp_b1, tmp_primal0);
+ VecAYPX(tmp_primal0, -1.0, subdomain->getVecRhsCoarse());
+
+ if (enableStokesMode) {
+ MatMult(subdomain->getMatCoarse(0, 1), vecSolInterface, tmp_primal1);
+ VecAXPY(tmp_primal0, -1.0, tmp_primal1);
+ }
- VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1);
KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
-
- // === Solve for u_b. ===
- VecCopy(subdomain->getVecRhsInterior(), tmp_b0);
- MatMultTranspose(mat_lagrange, vecSol, tmp_b1);
- VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
+ // === Solve for u_b. ===
MatMult(subdomain->getMatInteriorCoarse(), tmp_primal0, tmp_b1);
- VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
+ VecAXPY(tmp_b0, -1.0, tmp_b1);
subdomain->solveGlobal(tmp_b0, tmp_b0);
+
// === And recover AMDiS solution vectors. ===
recoverSolution(tmp_b0, tmp_primal0, vec);
diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.cc b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
index 83ec0837198c472953d945576f691e36d8e745b7..1cf7a349ca4f07763ad3812bc3ac005583b9f1dc 100644
--- a/AMDiS/src/parallel/PetscSolverFetiOperators.cc
+++ b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
@@ -48,30 +48,30 @@ namespace AMDiS {
MatShellGetContext(mat, &ctx);
FetiData* data = static_cast<FetiData*>(ctx);
- MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
+ MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
double wtime01 = MPI::Wtime();
- data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
+ data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
- MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
+ MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
MatMult(data->subSolver->getMatCoarseInterior(),
- data->tmp_vec_b, data->tmp_vec_primal);
+ data->tmp_vec_b0, data->tmp_vec_primal0);
wtime01 = MPI::Wtime();
- KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
+ KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal0, data->tmp_vec_primal0);
FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);
MatMult(data->subSolver->getMatInteriorCoarse(),
- data->tmp_vec_primal, data->tmp_vec_b);
+ data->tmp_vec_primal0, data->tmp_vec_b0);
wtime01 = MPI::Wtime();
- data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
+ data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
- MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
+ MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y);
VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
@@ -97,32 +97,68 @@ namespace AMDiS {
MatShellGetContext(mat, &ctx);
FetiData* data = static_cast<FetiData*>(ctx);
- MatMultTranspose(*(data->mat_lagrange), x_lagrange, data->tmp_vec_b);
- double wtime01 = MPI::Wtime();
- data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
+ // === Calculation of v_{B} ===
- FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
+ // tmp_vec_b0 = J^{T} \lambda
+ MatMultTranspose(*(data->mat_lagrange), x_lagrange, data->tmp_vec_b0);
+ // tmp_vec_b1 = A_{B\Gamma} u_{\Gamma}
+ MatMult(data->subSolver->getMatInteriorCoarse(1), x_interface, data->tmp_vec_b1);
+ // tmp_vec_b0 = A_{B\Gamma} u_{\Gamma} + J^{T} \lambda
+ VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1);
+
+
+ // == Calculation of v_{\Pi}
+
+ // tmp_vec_primal0 = A_{\Pi\Gamma} u_{\Gamma}
+ MatMult(data->subSolver->getMatCoarse(0, 1), x_interface, data->tmp_vec_primal0);
+
+
+ // === Calculate action of FETI-DP operator ===
- MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
+ double wtime01 = MPI::Wtime();
+ // tmp_vec_b0 = A_{BB}^{-1} v_{B}
+ data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
+ FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
+ // tmp_vec_primal1 = A_{\Pi B} A_{BB}^{-1} v_{B}
MatMult(data->subSolver->getMatCoarseInterior(),
- data->tmp_vec_b, data->tmp_vec_primal);
+ data->tmp_vec_b0, data->tmp_vec_primal1);
+
+ // tmp_vec_primal0 = v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B}
+ VecAXPY(data->tmp_vec_primal0, -1.0, data->tmp_vec_primal1);
wtime01 = MPI::Wtime();
- KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
+ // tmp_vec_primal0 = S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
+ KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal0, data->tmp_vec_primal0);
FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);
+ // tmp_vec_b1 = A_{B\Pi} S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
MatMult(data->subSolver->getMatInteriorCoarse(),
- data->tmp_vec_primal, data->tmp_vec_b);
+ data->tmp_vec_primal0, data->tmp_vec_b1);
wtime01 = MPI::Wtime();
- data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
+ // tmp_vec_b1 = A_{BB}^{-1} A_{B\Pi} S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
+ data->subSolver->solveGlobal(data->tmp_vec_b1, data->tmp_vec_b1);
FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
- MatMult(*(data->mat_lagrange), data->tmp_vec_b, y_lagrange);
+ // tmp_vec_b0 = A_{BB}^{-1} v_{B} - A_{BB}^{-1} A_{B\Pi} S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
+ VecAXPY(data->tmp_vec_b0, -1.0, data->tmp_vec_b1);
+
+
+ // === Calculate projection to interface and constraint variable ===
+
+ // y_interface = A_{\Gamma B} tmp_vec_b0
+ MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, y_interface);
+
+ // tmp_vec_primal0 = S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
+ // tmp_vec_interface = A_{\Gamma \Pi} tmp_vec_primal0
+ MatMult(data->subSolver->getMatCoarse(1, 0), data->tmp_vec_primal0, data->tmp_vec_interface);
+ // y_interface = A_{\Gamma B} tmp_vec_b0 + A_{\Gamma \Pi} tmp_vec_primal0
+ VecAXPY(y_interface, 1.0, data->tmp_vec_interface);
- VecAXPBY(y_lagrange, 1.0, 1.0, data->tmp_vec_lagrange);
+ // y_lagrange = J tmp_vec_b0
+ MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y_lagrange);
FetiTimings::fetiSolve += (MPI::Wtime() - wtime);
diff --git a/AMDiS/src/parallel/PetscSolverFetiStructs.h b/AMDiS/src/parallel/PetscSolverFetiStructs.h
index 340f7c13a0e33d86c3fa9e16c98e9c48611ab803..3da72354227174ee24e7827f8d34311cc73bcd77 100644
--- a/AMDiS/src/parallel/PetscSolverFetiStructs.h
+++ b/AMDiS/src/parallel/PetscSolverFetiStructs.h
@@ -57,14 +57,16 @@ namespace AMDiS {
Mat *mat_lagrange;
/// Temporal vectors on the B variables.
- Vec tmp_vec_b;
+ Vec tmp_vec_b0, tmp_vec_b1;
/// Temporal vector on the primal variables.
- Vec tmp_vec_primal;
+ Vec tmp_vec_primal0, tmp_vec_primal1;
/// Temporal vector on the lagrange variables.
Vec tmp_vec_lagrange;
+ Vec tmp_vec_interface;
+
PetscSolver* subSolver;
/// Pointer to the solver of the schur complement on the primal variables.
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 75379f2b8cc4be899858186a103a2d6684eeb86b..b9da1f97ceb0086880cefb0f707931281c119243 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -161,22 +161,12 @@ namespace AMDiS {
bool isColCoarse =
isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor));
- if (isColCoarse) {
- if (isRowCoarse) {
- cols.push_back(col(*icursor));
- values.push_back(value(*icursor));
- } else {
- colsOther.push_back(col(*icursor));
- valuesOther.push_back(value(*icursor));
- }
+ if (isColCoarse == isRowCoarse) {
+ cols.push_back(col(*icursor));
+ values.push_back(value(*icursor));
} else {
- if (isRowCoarse) {
- colsOther.push_back(col(*icursor));
- valuesOther.push_back(value(*icursor));
- } else {
- cols.push_back(col(*icursor));
- values.push_back(value(*icursor));
- }
+ colsOther.push_back(col(*icursor));
+ valuesOther.push_back(value(*icursor));
}
} // for each nnz in row
@@ -341,8 +331,9 @@ namespace AMDiS {
MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
}
- MatSetNullSpace(getMatInterior(), matNullspace);
- KSPSetNullSpace(kspInterior, matNullspace);
+ MSG("NULLSPACE IS NOT REMOVED!\n");
+ // MatSetNullSpace(getMatInterior(), matNullspace);
+ // KSPSetNullSpace(kspInterior, matNullspace);
// === Remove null space, if requested. ===