Commit b81cf528 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Go into weekend.

parent 6de23628
...@@ -768,9 +768,9 @@ namespace AMDiS { ...@@ -768,9 +768,9 @@ namespace AMDiS {
fetiData.subSolver = subdomain; fetiData.subSolver = subdomain;
fetiData.ksp_schur_primal = &ksp_schur_primal; 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); lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
primalDofMap.createVec(fetiData.tmp_vec_primal); primalDofMap.createVec(fetiData.tmp_vec_primal0);
if (enableStokesMode == false) { if (enableStokesMode == false) {
MatCreateShell(mpiCommGlobal, MatCreateShell(mpiCommGlobal,
...@@ -782,6 +782,10 @@ namespace AMDiS { ...@@ -782,6 +782,10 @@ namespace AMDiS {
MatShellSetOperation(mat_feti, MATOP_MULT, MatShellSetOperation(mat_feti, MATOP_MULT,
(void(*)(void))petscMultMatFeti); (void(*)(void))petscMultMatFeti);
} else { } else {
localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
primalDofMap.createVec(fetiData.tmp_vec_primal1);
interfaceDofMap.createVec(fetiData.tmp_vec_interface);
MatCreateShell(mpiCommGlobal, MatCreateShell(mpiCommGlobal,
interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
...@@ -799,6 +803,27 @@ namespace AMDiS { ...@@ -799,6 +803,27 @@ namespace AMDiS {
KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000); KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
KSPSetFromOptions(ksp_feti); 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. === // === Create FETI-DP preconditioner object. ===
...@@ -913,9 +938,16 @@ namespace AMDiS { ...@@ -913,9 +938,16 @@ namespace AMDiS {
fetiData.subSolver = NULL; fetiData.subSolver = NULL;
fetiData.ksp_schur_primal = PETSC_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_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); MatDestroy(&mat_feti);
KSPDestroy(&ksp_feti); KSPDestroy(&ksp_feti);
...@@ -1520,30 +1552,39 @@ namespace AMDiS { ...@@ -1520,30 +1552,39 @@ namespace AMDiS {
FetiTimings::fetiPreconditioner); FetiTimings::fetiPreconditioner);
} }
// === Solve for u_primals. === // === Solve for u_primals. ===
VecCopy(subdomain->getVecRhsCoarse(), tmp_primal0); if (!enableStokesMode) {
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); MatMultTranspose(mat_lagrange, vecSol, tmp_b0);
subdomain->solveGlobal(tmp_b0, tmp_b0); VecAYPX(tmp_b0, -1.0, subdomain->getVecRhsInterior());
MatMult(subdomain->getMatCoarseInterior(), tmp_b0, tmp_primal1); } 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); KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
// === Solve for u_b. === // === Solve for u_b. ===
VecCopy(subdomain->getVecRhsInterior(), tmp_b0);
MatMultTranspose(mat_lagrange, vecSol, tmp_b1);
VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
MatMult(subdomain->getMatInteriorCoarse(), tmp_primal0, tmp_b1); 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); subdomain->solveGlobal(tmp_b0, tmp_b0);
// === And recover AMDiS solution vectors. === // === And recover AMDiS solution vectors. ===
recoverSolution(tmp_b0, tmp_primal0, vec); recoverSolution(tmp_b0, tmp_primal0, vec);
......
...@@ -48,30 +48,30 @@ namespace AMDiS { ...@@ -48,30 +48,30 @@ namespace AMDiS {
MatShellGetContext(mat, &ctx); MatShellGetContext(mat, &ctx);
FetiData* data = static_cast<FetiData*>(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(); 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); 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(), MatMult(data->subSolver->getMatCoarseInterior(),
data->tmp_vec_b, data->tmp_vec_primal); data->tmp_vec_b0, data->tmp_vec_primal0);
wtime01 = MPI::Wtime(); 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); FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);
MatMult(data->subSolver->getMatInteriorCoarse(), MatMult(data->subSolver->getMatInteriorCoarse(),
data->tmp_vec_primal, data->tmp_vec_b); data->tmp_vec_primal0, data->tmp_vec_b0);
wtime01 = MPI::Wtime(); 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); 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); VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
...@@ -97,32 +97,68 @@ namespace AMDiS { ...@@ -97,32 +97,68 @@ namespace AMDiS {
MatShellGetContext(mat, &ctx); MatShellGetContext(mat, &ctx);
FetiData* data = static_cast<FetiData*>(ctx); FetiData* data = static_cast<FetiData*>(ctx);
MatMultTranspose(*(data->mat_lagrange), x_lagrange, data->tmp_vec_b);
double wtime01 = MPI::Wtime(); // === Calculation of v_{B} ===
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
// 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);
FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange); // === Calculate action of FETI-DP operator ===
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(), 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(); 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); 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(), MatMult(data->subSolver->getMatInteriorCoarse(),
data->tmp_vec_primal, data->tmp_vec_b); data->tmp_vec_primal0, data->tmp_vec_b1);
wtime01 = MPI::Wtime(); 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); 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); FetiTimings::fetiSolve += (MPI::Wtime() - wtime);
......
...@@ -57,14 +57,16 @@ namespace AMDiS { ...@@ -57,14 +57,16 @@ namespace AMDiS {
Mat *mat_lagrange; Mat *mat_lagrange;
/// Temporal vectors on the B variables. /// Temporal vectors on the B variables.
Vec tmp_vec_b; Vec tmp_vec_b0, tmp_vec_b1;
/// Temporal vector on the primal variables. /// Temporal vector on the primal variables.
Vec tmp_vec_primal; Vec tmp_vec_primal0, tmp_vec_primal1;
/// Temporal vector on the lagrange variables. /// Temporal vector on the lagrange variables.
Vec tmp_vec_lagrange; Vec tmp_vec_lagrange;
Vec tmp_vec_interface;
PetscSolver* subSolver; PetscSolver* subSolver;
/// Pointer to the solver of the schur complement on the primal variables. /// Pointer to the solver of the schur complement on the primal variables.
......
...@@ -161,23 +161,13 @@ namespace AMDiS { ...@@ -161,23 +161,13 @@ namespace AMDiS {
bool isColCoarse = bool isColCoarse =
isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor)); isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor));
if (isColCoarse) { if (isColCoarse == isRowCoarse) {
if (isRowCoarse) {
cols.push_back(col(*icursor)); cols.push_back(col(*icursor));
values.push_back(value(*icursor)); values.push_back(value(*icursor));
} else { } else {
colsOther.push_back(col(*icursor)); colsOther.push_back(col(*icursor));
valuesOther.push_back(value(*icursor)); valuesOther.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));
}
}
} // for each nnz in row } // for each nnz in row
...@@ -341,8 +331,9 @@ namespace AMDiS { ...@@ -341,8 +331,9 @@ namespace AMDiS {
MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace); MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
} }
MatSetNullSpace(getMatInterior(), matNullspace); MSG("NULLSPACE IS NOT REMOVED!\n");
KSPSetNullSpace(kspInterior, matNullspace); // MatSetNullSpace(getMatInterior(), matNullspace);
// KSPSetNullSpace(kspInterior, matNullspace);
// === Remove null space, if requested. === // === Remove null space, if requested. ===
......
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