diff --git a/AMDiS/src/AdaptInstationary.cc b/AMDiS/src/AdaptInstationary.cc
index de4bf478266e879555842163176775c4d94a709c..996bf0c8064be013d87f039445658283350428e1 100644
--- a/AMDiS/src/AdaptInstationary.cc
+++ b/AMDiS/src/AdaptInstationary.cc
@@ -17,6 +17,11 @@
#include "ProblemTimeInterface.h"
#include "Serializer.h"
+#if HAVE_PARALLEL_DOMAIN_AMDIS
+#include "parallel/MeshDistributor.h"
+#include <petsc.h>
+#endif
+
namespace AMDiS {
AdaptInstationary::AdaptInstationary(std::string name,
@@ -31,8 +36,8 @@ namespace AMDiS {
{
FUNCNAME("AdaptInstationary::AdaptInstationary()");
-// MSG("You make use of the obsolete constructor AdaptInstationary::AdaptInstationary(...)!\n");
-// MSG("Please use the constructor that uses references instead of pointers!\n");
+ MSG("You make use of the obsolete constructor AdaptInstationary::AdaptInstationary(...)!\n");
+ MSG("Please use the constructor that uses references instead of pointers!\n");
initConstructor(problemStat, info, initialInfo, initialTimestampSet);
}
@@ -262,6 +267,8 @@ namespace AMDiS {
{
FUNCNAME("AdaptInstationary::oneTimestep()");
+ MSG("ONE TIMESTEP!\n");
+
adaptInfo->setTimestepIteration(0);
switch (strategy) {
@@ -294,6 +301,10 @@ namespace AMDiS {
TEST_EXIT(adaptInfo->getTimestep() > 0)("timestep <= 0!\n");
+#if HAVE_PARALLEL_DOMAIN_AMDIS
+ MeshDistributor::globalMeshDistributor->initParallelization();
+#endif
+
if (adaptInfo->getTimestepNumber() == 0) {
adaptInfo->setTime(adaptInfo->getStartTime());
initialAdaptInfo->setStartTime(adaptInfo->getStartTime());
@@ -315,10 +326,7 @@ namespace AMDiS {
if (breakWhenStable && (adaptInfo->getSolverIterations() == 0))
break;
- /*
- if (adaptInfo->forceBreak)
- break;
- */
+
// Check if there is a runtime limitation. If there is a runtime limitation
// and there is no more time for a next adaption loop, than return the error
// code for rescheduling the problem and break the adaption loop.
@@ -328,6 +336,11 @@ namespace AMDiS {
}
}
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ MeshDistributor::globalMeshDistributor->exitParallelization();
+ PetscFinalize();
+#endif
+
return errorCode;
}
diff --git a/AMDiS/src/AdaptStationary.cc b/AMDiS/src/AdaptStationary.cc
index 75c94d089272dd6ef9bd99dd0675ead73d40bfa1..a95c7184bc04e302e04a1bb4ad12e20cb67e7bd2 100644
--- a/AMDiS/src/AdaptStationary.cc
+++ b/AMDiS/src/AdaptStationary.cc
@@ -16,6 +16,11 @@
#include "ProblemIterationInterface.h"
#include <math.h>
+#if HAVE_PARALLEL_DOMAIN_AMDIS
+#include "parallel/MeshDistributor.h"
+#include <petsc.h>
+#endif
+
namespace AMDiS {
AdaptStationary::AdaptStationary(std::string name,
@@ -25,8 +30,8 @@ namespace AMDiS {
{
FUNCNAME("AdaptStationary::AdaptStationary()");
-// MSG("You make use of the obsolete constructor AdaptStationary::AdaptStationary(...)!\n");
-// MSG("Please use the constructor that uses references instead of pointers!\n");
+ MSG("You make use of the obsolete constructor AdaptStationary::AdaptStationary(...)!\n");
+ MSG("Please use the constructor that uses references instead of pointers!\n");
initialize();
}
@@ -45,6 +50,10 @@ namespace AMDiS {
{
FUNCNAME("AdaptStationary::adapt()");
+#if HAVE_PARALLEL_DOMAIN_AMDIS
+ MeshDistributor::globalMeshDistributor->initParallelization();
+#endif
+
// initial iteration
if (adaptInfo->getSpaceIteration() == -1) {
problemIteration->beginIteration(adaptInfo);
@@ -68,6 +77,11 @@ namespace AMDiS {
adaptInfo->incSpaceIteration();
}
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ MeshDistributor::globalMeshDistributor->exitParallelization();
+ PetscFinalize();
+#endif
+
return 0;
}
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index b5ed15196c31f02e49dc58fbd4eeddb87c87dc9a..e0310499e918a142c15c9736604217f3036c3c19 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -143,6 +143,7 @@ namespace AMDiS {
TEST_EXIT(feSpace)("No FE space has been defined for the mesh distributor!\n");
TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
+ TEST_EXIT(!initialized)("MeshDistributor is already initialized!\n");
#ifdef HAVE_ZOLTAN
int a = 0;
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 232037765d86319710bdb75c82fdf527f225f483..9b076f923310f109f033c725ea66829871787567 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -47,7 +47,7 @@ namespace AMDiS {
void *ctx;
MatShellGetContext(mat, &ctx);
- PetscFetiData* data = static_cast<PetscFetiData*>(ctx);
+ FetiData* data = static_cast<FetiData*>(ctx);
// y = L inv(K_BB) trans(L) x
MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
@@ -71,61 +71,139 @@ namespace AMDiS {
// y = PC * x
- PetscErrorCode petscApplyFetiPrecon(PC pc, Vec x, Vec y)
+ PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
{
+ // Get data for the preconditioner
void *ctx;
PCShellGetContext(pc, &ctx);
- PetscFetiPreconData* data = static_cast<PetscFetiPreconData*>(ctx);
+ FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
+ // Multiply with scaled Lagrange constraint matrix.
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);
+ // === Restriction of the B nodes to the boundary nodes. ===
- PetscScalar *local_bound;
- VecGetArray(data->tmp_vec_bound0, &local_bound);
+ int nLocalB;
+ int nLocalDuals;
+ VecGetLocalSize(data->tmp_vec_b, &nLocalB);
+ VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
- for (int i = sizeB - sizeBound, j = 0; i < sizeB; i++, j++)
- local_bound[j] = local_b[i];
+ PetscScalar *local_b, *local_duals;
+ VecGetArray(data->tmp_vec_b, &local_b);
+ VecGetArray(data->tmp_vec_duals0, &local_duals);
+
+ for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
+ local_duals[j] = local_b[i];
VecRestoreArray(data->tmp_vec_b, &local_b);
- VecRestoreArray(data->tmp_vec_bound0, &local_bound);
+ VecRestoreArray(data->tmp_vec_duals0, &local_duals);
+ // === K_DD - K_DI inv(K_II) K_ID ===
- MatMult(*(data->mat_bound_bound), data->tmp_vec_bound0, data->tmp_vec_bound1);
+ MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
- MatMult(*(data->mat_interior_bound), data->tmp_vec_bound0, data->tmp_vec_interior);
+ MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, 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);
+ MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);
+
+ VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);
+
+
+ // === Prolongation from local dual nodes to B nodes.
+
+ VecGetArray(data->tmp_vec_b, &local_b);
+ VecGetArray(data->tmp_vec_duals0, &local_duals);
+
+ for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
+ local_b[i] = local_duals[j];
+
+ VecRestoreArray(data->tmp_vec_b, &local_b);
+ VecRestoreArray(data->tmp_vec_duals0, &local_duals);
+
+
+ // Multiply with scaled Lagrange constraint matrix.
+ MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);
+
+ return 0;
+ }
+
+
+ // y = PC * x
+ PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
+ {
+ // Get data for the preconditioner
+ void *ctx;
+ PCShellGetContext(pc, &ctx);
+ FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);
+
+ // Multiply with scaled Lagrange constraint matrix.
+ MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
- VecAXPBY(data->tmp_vec_bound0, 1.0, -1.0, data->tmp_vec_bound1);
+ // === Restriction of the B nodes to the boundary nodes. ===
+ int nLocalB;
+ int nLocalDuals;
+ VecGetLocalSize(data->tmp_vec_b, &nLocalB);
+ VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
+ PetscScalar *local_b, *local_duals;
VecGetArray(data->tmp_vec_b, &local_b);
- VecGetArray(data->tmp_vec_bound0, &local_bound);
+ VecGetArray(data->tmp_vec_duals0, &local_duals);
- for (int i = sizeB - sizeBound, j = 0; i < sizeB; i++, j++)
- local_b[i] = local_bound[j];
+ for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
+ local_duals[j] = local_b[i];
VecRestoreArray(data->tmp_vec_b, &local_b);
- VecRestoreArray(data->tmp_vec_bound0, &local_bound);
+ VecRestoreArray(data->tmp_vec_duals0, &local_duals);
+
+
+ // === K_DD ===
+
+ MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
+
+ // === Prolongation from local dual nodes to B nodes.
+ VecGetArray(data->tmp_vec_b, &local_b);
+ VecGetArray(data->tmp_vec_duals1, &local_duals);
+
+ for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
+ local_b[i] = local_duals[j];
+ VecRestoreArray(data->tmp_vec_b, &local_b);
+ VecRestoreArray(data->tmp_vec_duals0, &local_duals);
+
+ // Multiply with scaled Lagrange constraint matrix.
MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);
return 0;
}
+ PetscSolverFeti::PetscSolverFeti()
+ : PetscSolver(),
+ nComponents(-1)
+ {
+ FUNCNAME("PetscSolverFeti::PetscSolverFeti()");
+
+ string preconditionerName = "";
+ Parameters::get("parallel->solver->precon", preconditionerName);
+ if (preconditionerName == "" || preconditionerName == "none") {
+ fetiPreconditioner = FETI_NONE;
+ } else if (preconditionerName == "dirichlet") {
+ fetiPreconditioner = FETI_DIRICHLET;
+ } else if (preconditionerName == "lumped") {
+ fetiPreconditioner = FETI_LUMPED;
+ } else {
+ ERROR_EXIT("Preconditioner \"%s\" not available!\n", preconditionerName.c_str());
+ }
+ }
+
+
void PetscSolverFeti::updateDofData()
{
FUNCNAME("PetscSolverFeti::updateDofData()");
@@ -438,7 +516,7 @@ namespace AMDiS {
it->second = nLocalInterior + rStartB;
nLocalInterior++;
}
- nLocalBound = duals.size();
+ nLocalDuals = duals.size();
TEST_EXIT_DBG(nLocalInterior + primals.size() + duals.size() ==
static_cast<unsigned int>(admin->getUsedDofs()))
@@ -561,23 +639,23 @@ namespace AMDiS {
// === 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;
- petscFetiData.mat_lagrange = &mat_lagrange;
- petscFetiData.ksp_b = &ksp_b;
- petscFetiData.ksp_schur_primal = &ksp_schur_primal;
+ fetiData.mat_primal_primal = &mat_primal_primal;
+ fetiData.mat_primal_b = &mat_primal_b;
+ fetiData.mat_b_primal = &mat_b_primal;
+ fetiData.mat_lagrange = &mat_lagrange;
+ fetiData.ksp_b = &ksp_b;
+ fetiData.ksp_schur_primal = &ksp_schur_primal;
- VecDuplicate(f_b, &(petscFetiData.tmp_vec_b));
- VecDuplicate(f_primal, &(petscFetiData.tmp_vec_primal));
- MatGetVecs(mat_lagrange, PETSC_NULL, &(petscFetiData.tmp_vec_lagrange));
+ VecDuplicate(f_b, &(fetiData.tmp_vec_b));
+ VecDuplicate(f_primal, &(fetiData.tmp_vec_primal));
+ MatGetVecs(mat_lagrange, PETSC_NULL, &(fetiData.tmp_vec_lagrange));
MatCreateShell(PETSC_COMM_WORLD,
nRankLagrange * nComponents, nRankLagrange * nComponents,
nOverallLagrange * nComponents, nOverallLagrange * nComponents,
- &petscFetiData, &mat_feti);
+ &fetiData, &mat_feti);
MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);
@@ -587,36 +665,54 @@ namespace AMDiS {
KSPSetFromOptions(ksp_feti);
- // === Create FETI-DP Dirichlet preconditioner object. ===
+ // === Create FETI-DP 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);
+ if (fetiPreconditioner != FETI_NONE) {
+ MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
+ MatScale(mat_lagrange_scaled, 0.5);
+ }
-
- 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.nLocalInterior = nLocalInterior;
- petscFetiPreconData.nLocalBound = nLocalBound;
-
- 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);
+ switch (fetiPreconditioner) {
+ case FETI_DIRICHLET:
+ 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);
+
+ fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
+ fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
+ fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
+ fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
+ fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
+ fetiDirichletPreconData.ksp_interior = &ksp_interior;
+
+ VecDuplicate(f_b, &(fetiDirichletPreconData.tmp_vec_b));
+ MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals0));
+ MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals1));
+ MatGetVecs(mat_interior_interior, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_interior));
+
+ KSPGetPC(ksp_feti, &precon_feti);
+ PCSetType(precon_feti, PCSHELL);
+ PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
+ PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
+
+ break;
+
+ case FETI_LUMPED:
+ fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
+ fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;
+
+ VecDuplicate(f_b, &(fetiLumpedPreconData.tmp_vec_b));
+ MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals0));
+ MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals1));
+
+ KSPGetPC(ksp_feti, &precon_feti);
+ PCSetType(precon_feti, PCSHELL);
+ PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
+ PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
+
+ break;
+ }
}
@@ -626,37 +722,50 @@ namespace AMDiS {
// === Destroy FETI-DP solver object. ===
- petscFetiData.mat_primal_primal = PETSC_NULL;
- petscFetiData.mat_primal_b = PETSC_NULL;
- petscFetiData.mat_b_primal = PETSC_NULL;
- petscFetiData.mat_lagrange = PETSC_NULL;
- petscFetiData.ksp_b = PETSC_NULL;
- petscFetiData.ksp_schur_primal = PETSC_NULL;
+ fetiData.mat_primal_primal = PETSC_NULL;
+ fetiData.mat_primal_b = PETSC_NULL;
+ fetiData.mat_b_primal = PETSC_NULL;
+ fetiData.mat_lagrange = PETSC_NULL;
+ fetiData.ksp_b = PETSC_NULL;
+ fetiData.ksp_schur_primal = PETSC_NULL;
- VecDestroy(&petscFetiData.tmp_vec_b);
- VecDestroy(&petscFetiData.tmp_vec_primal);
- VecDestroy(&petscFetiData.tmp_vec_lagrange);
+ VecDestroy(&fetiData.tmp_vec_b);
+ VecDestroy(&fetiData.tmp_vec_primal);
+ VecDestroy(&fetiData.tmp_vec_lagrange);
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);
+ // === Destroy FETI-DP preconditioner object. ===
+
+ switch (fetiPreconditioner) {
+ case FETI_DIRICHLET:
+ KSPDestroy(&ksp_interior);
+
+ fetiDirichletPreconData.mat_lagrange_scaled = NULL;
+ fetiDirichletPreconData.mat_interior_interior = NULL;
+ fetiDirichletPreconData.mat_duals_duals = NULL;
+ fetiDirichletPreconData.mat_interior_duals = NULL;
+ fetiDirichletPreconData.mat_duals_interior = NULL;
+ fetiDirichletPreconData.ksp_interior = NULL;
+
+ VecDestroy(&fetiDirichletPreconData.tmp_vec_b);
+ VecDestroy(&fetiDirichletPreconData.tmp_vec_duals0);
+ VecDestroy(&fetiDirichletPreconData.tmp_vec_duals1);
+ VecDestroy(&fetiDirichletPreconData.tmp_vec_interior);
+ MatDestroy(&mat_lagrange_scaled);
+ break;
+
+ case FETI_LUMPED:
+ fetiLumpedPreconData.mat_lagrange_scaled = NULL;
+ fetiLumpedPreconData.mat_duals_duals = NULL;
+
+ VecDestroy(&fetiLumpedPreconData.tmp_vec_b);
+ VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0);
+ VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1);
+ break;
+ }
}
@@ -753,6 +862,10 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::fillPetscMatrix()");
+ PetscLogStage stageFetiCreate;
+ PetscLogStageRegister("Stage 0 CREATE FETI-DP", &stageFetiCreate);
+ PetscLogStagePush(stageFetiCreate);
+
nComponents = vec->getSize();
// === Create all sets and indices. ===
@@ -767,7 +880,7 @@ namespace AMDiS {
int nRowsRankPrimal = nRankPrimals * nComponents;
int nRowsOverallPrimal = nOverallPrimals * nComponents;
int nRowsInterior = nLocalInterior * nComponents;
- int nRowsBound = nLocalBound * nComponents;
+ int nRowsDual = nLocalDuals * nComponents;
MatCreateMPIAIJ(PETSC_COMM_WORLD,
nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB,
@@ -789,23 +902,26 @@ namespace AMDiS {
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);
+ // === Create matrices for FETI-DP preconditioner. ===
+
+ if (fetiPreconditioner != FETI_NONE)
+ MatCreateSeqAIJ(PETSC_COMM_SELF,
+ nRowsDual, nRowsDual, 100, PETSC_NULL,
+ &mat_duals_duals);
+
+ if (fetiPreconditioner == FETI_DIRICHLET) {
+ MatCreateSeqAIJ(PETSC_COMM_SELF,
+ nRowsInterior, nRowsInterior, 100, PETSC_NULL,
+ &mat_interior_interior);
+
+ MatCreateSeqAIJ(PETSC_COMM_SELF,
+ nRowsInterior, nRowsDual, 100, PETSC_NULL,
+ &mat_interior_duals);
+
+ MatCreateSeqAIJ(PETSC_COMM_SELF,
+ nRowsDual, nRowsInterior, 100, PETSC_NULL,
+ &mat_duals_interior);
+ }
// === Prepare traverse of sequentially created matrices. ===
@@ -967,33 +1083,51 @@ namespace AMDiS {
// === 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);
+ switch (fetiPreconditioner) {
+ case FETI_DIRICHLET:
+ {
+ 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_duals, 1, &rowIndex2, colsLocalOther.size(),
+ &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
+ } else {
+ int rowIndex2 =
+ (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
+
+ MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
+ &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
+
+ if (colsLocalOther.size())
+ MatSetValues(mat_duals_interior, 1, &rowIndex2, colsLocalOther.size(),
+ &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
+
+ }
+ }
+ break;
- }
+ case FETI_LUMPED:
+ {
+ int rowIndex = globalIndexB[*cursor] - rStartB;
+
+ if (rowIndex >= nLocalInterior) {
+ int rowIndex2 =
+ (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
+
+ MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
+ &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
+ }
+ }
+ break;
+ }
}
-
-
}
}
}
@@ -1012,19 +1146,25 @@ namespace AMDiS {
MatAssemblyBegin(mat_primal_b, MAT_FINAL_ASSEMBLY);
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);
+ // === Start global assembly procedure for preconditioner matrices. ===
- MatAssemblyBegin(mat_interior_bound, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(mat_interior_bound, MAT_FINAL_ASSEMBLY);
+ if (fetiPreconditioner != FETI_NONE) {
+ MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY);
+ }
- MatAssemblyBegin(mat_bound_interior, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(mat_bound_interior, 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);
+ }
// === Create and fill PETSc's right hand side vectors. ===
@@ -1080,6 +1220,8 @@ namespace AMDiS {
// === Create PETSc solver for the FETI-DP operator. ===
createFetiKsp();
+
+ PetscLogStagePop();
}
@@ -1225,6 +1367,10 @@ namespace AMDiS {
// === Create solver for the non primal (thus local) variables. ===
+ PetscLogStage stagePreFeti;
+ PetscLogStageRegister("Stage 0 PRE FETI-DP", &stagePreFeti);
+ PetscLogStagePush(stagePreFeti);
+
KSPCreate(PETSC_COMM_WORLD, &ksp_b);
KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN);
KSPSetOptionsPrefix(ksp_b, "solver_b_");
@@ -1269,9 +1415,17 @@ namespace AMDiS {
// === Solve with FETI-DP operator. ===
- MSG("START FETI SOLVE!\n");
+ PetscLogStagePop();
+ PetscLogStage stageFeti0, stageFeti1, stageFeti2;
+ PetscLogStageRegister("Stage 0 FETI-DP", &stageFeti0);
+ PetscLogStageRegister("Stage 1 FETI-DP", &stageFeti1);
+ PetscLogStageRegister("Stage 2 FETI-DP", &stageFeti2);
+ PetscLogStagePush(stageFeti0);
+
KSPSolve(ksp_feti, vec_rhs, vec_rhs);
+ PetscLogStagePop();
+ PetscLogStagePush(stageFeti1);
// === Solve for u_primals. ===
@@ -1301,6 +1455,9 @@ namespace AMDiS {
KSPSolve(ksp_b, tmp_b0, tmp_b0);
+ PetscLogStagePop();
+ PetscLogStagePush(stageFeti2);
+
// === And recover AMDiS solution vectors. ===
@@ -1335,10 +1492,16 @@ namespace AMDiS {
// === Destroy preconditioner data structures. ===
- MatDestroy(&mat_interior_interior);
- MatDestroy(&mat_bound_bound);
- MatDestroy(&mat_interior_bound);
- MatDestroy(&mat_bound_interior);
+ if (fetiPreconditioner != FETI_NONE)
+ MatDestroy(&mat_duals_duals);
+
+ if (fetiPreconditioner == FETI_DIRICHLET) {
+ MatDestroy(&mat_interior_interior);
+ MatDestroy(&mat_interior_duals);
+ MatDestroy(&mat_duals_interior);
+ }
+
+ PetscLogStagePop();
}
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 0a29a9414cfd63a37c65dcfa6ae514d44753b051..6a03db6c9496dcdb78df7218a83170c2181388fb 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -60,7 +60,7 @@ namespace AMDiS {
* the system matrix reduced to the Lagrange multipliers.
* \ref petscMultMatFeti
*/
- struct PetscFetiData {
+ struct FetiData {
/// Pointers to the matrix containing the primal variables.
Mat *mat_primal_primal;
@@ -87,11 +87,11 @@ namespace AMDiS {
};
- struct PetscFetiPreconData {
+ struct FetiDirichletPreconData {
/// Matrix of scaled Lagrange variables.
Mat *mat_lagrange_scaled;
- Mat *mat_interior_interior, *mat_bound_bound, *mat_interior_bound, *mat_bound_interior;
+ Mat *mat_interior_interior, *mat_duals_duals, *mat_interior_duals, *mat_duals_interior;
/// Pointer to the solver for \ref PetscSolverFeti::mat_bb.
KSP *ksp_interior;
@@ -99,27 +99,42 @@ namespace AMDiS {
/// Temporal vector on the B variables.
Vec tmp_vec_b;
- Vec tmp_vec_bound0, tmp_vec_bound1;
-
- Vec tmp_vec_interior;
+ /// Temporal vector on the dual variables.
+ Vec tmp_vec_duals0, tmp_vec_duals1;
+ /// Temporal vector on the interior variables.
+ Vec tmp_vec_interior;
+ };
- int nLocalInterior;
- int nLocalBound;
+ struct FetiLumpedPreconData {
+ /// Matrix of scaled Lagrange variables.
+ Mat *mat_lagrange_scaled;
+
+ Mat *mat_duals_duals;
+
+ /// Temporal vector on the B variables.
+ Vec tmp_vec_b;
+
+ /// Temporal vector on the dual variables.
+ Vec tmp_vec_duals0, tmp_vec_duals1;
};
+ typedef enum {
+ FETI_NONE = 0,
+ FETI_DIRICHLET = 1,
+ FETI_LUMPED = 2
+ } FetiPreconditioner;
+
+
/** \brief
* FETI-DP implementation based on PETSc.
*/
class PetscSolverFeti : public PetscSolver
{
public:
- PetscSolverFeti()
- : PetscSolver(),
- nComponents(-1)
- {}
+ PetscSolverFeti();
/// Assemble the sequentially created matrices and vectors to the
/// global matrices and vectors required by the FETI-DP method.
@@ -288,22 +303,29 @@ namespace AMDiS {
Mat mat_feti;
/// Data for MatMult operation in matrix \ref mat_feti
- PetscFetiData petscFetiData;
+ FetiData fetiData;
+ /// Defines which preconditioner should be used to solve the reduced
+ /// FETI-DP system.
+ FetiPreconditioner fetiPreconditioner;
+
+ /// Preconditioner object for the reduced FETI-DP system.
+ PC precon_feti;
Mat mat_lagrange_scaled;
- PC precon_feti;
+ FetiDirichletPreconData fetiDirichletPreconData;
- PetscFetiPreconData petscFetiPreconData;
+ FetiLumpedPreconData fetiLumpedPreconData;
- Mat mat_interior_interior, mat_bound_bound, mat_interior_bound, mat_bound_interior;
+ Mat mat_interior_interior, mat_duals_duals, mat_interior_duals, mat_duals_interior;
KSP ksp_interior;
int nLocalInterior;
- int nLocalBound;
+ // Number of local nodes that are duals.
+ int nLocalDuals;
};
#endif