Commit f624d7ca authored by Thomas Witkowski's avatar Thomas Witkowski

Added nullspace support.

parent 705566b3
...@@ -323,7 +323,7 @@ namespace AMDiS { ...@@ -323,7 +323,7 @@ namespace AMDiS {
int use_defaults_int = 1; int use_defaults_int = 1;
int parallel_division_int = 1; int parallel_division_int = 1;
int use_arithmetic_int = 1; int use_arithmetic_int = 0;
int use_adaptive_int = 0; int use_adaptive_int = 0;
// MSG("call to \"bddcml_setup_preconditioner\" with the following arguments (each in one line):\n"); // MSG("call to \"bddcml_setup_preconditioner\" with the following arguments (each in one line):\n");
// MSG(" %d\n", matrixtype); // MSG(" %d\n", matrixtype);
......
...@@ -27,14 +27,16 @@ namespace AMDiS { ...@@ -27,14 +27,16 @@ namespace AMDiS {
coarseSpaceMap(NULL), coarseSpaceMap(NULL),
mpiRank(-1), mpiRank(-1),
kspPrefix(""), kspPrefix(""),
removeRhsNullSpace(false) removeRhsNullspace(false),
hasConstantNullspace(false)
{ {
string kspStr = ""; string kspStr = "";
Parameters::get("parallel->solver->petsc->ksp", kspStr); Parameters::get("parallel->solver->petsc->ksp", kspStr);
if (kspStr != "") if (kspStr != "")
PetscOptionsInsertString(kspStr.c_str()); PetscOptionsInsertString(kspStr.c_str());
Parameters::get("parallel->remove rhs null space", removeRhsNullSpace); Parameters::get("parallel->remove rhs null space", removeRhsNullspace);
Parameters::get("parallel->has constant null space", hasConstantNullspace);
} }
......
...@@ -116,9 +116,9 @@ namespace AMDiS { ...@@ -116,9 +116,9 @@ namespace AMDiS {
kspPrefix = s; kspPrefix = s;
} }
void setRemoveRhsNullSpace(bool b) void setRemoveRhsNullspace(bool b)
{ {
removeRhsNullSpace = b; removeRhsNullspace = b;
} }
/// Adds a new vector to the basis of the operator's nullspace. /// Adds a new vector to the basis of the operator's nullspace.
...@@ -253,7 +253,9 @@ namespace AMDiS { ...@@ -253,7 +253,9 @@ namespace AMDiS {
/// If true, the constant null space is projected out of the RHS vector. It /// If true, the constant null space is projected out of the RHS vector. It
/// depends on the specific PETSc solver if it considers this value. /// depends on the specific PETSc solver if it considers this value.
bool removeRhsNullSpace; bool removeRhsNullspace;
bool hasConstantNullspace;
}; };
......
...@@ -363,19 +363,7 @@ namespace AMDiS { ...@@ -363,19 +363,7 @@ namespace AMDiS {
// === Remove Dirichlet BC DOFs. === // === Remove Dirichlet BC DOFs. ===
// removeDirichletBcDofs(vec); // removeDirichletBcDofs(vec);
// === Remove null space, if requested. ===
if (removeRhsNullSpace) {
TEST_EXIT_DBG(coarseSpaceMap == NULL)("Not supported!\n");
MSG("Remove constant null space from the RHS!\n");
MatNullSpace sp;
MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &sp);
MatNullSpaceRemove(sp, rhsInterior, PETSC_NULL);
MatNullSpaceDestroy(&sp);
}
} }
...@@ -397,35 +385,53 @@ namespace AMDiS { ...@@ -397,35 +385,53 @@ namespace AMDiS {
VecAssemblyEnd(petscSolVec); VecAssemblyEnd(petscSolVec);
} }
MatNullSpace matNullSpace;
MatNullSpace matNullspace;
Vec nullspaceBasis; Vec nullspaceBasis;
if (nullspace.size() > 0) { if (nullspace.size() > 0 || hasConstantNullspace) {
TEST_EXIT_DBG(nullspace.size() == 1)("Not yet implemented!\n"); TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
VecDuplicate(petscSolVec, &nullspaceBasis); if (nullspace.size() > 0) {
VecDuplicate(petscSolVec, &nullspaceBasis);
for (int i = 0; i < nComponents; i++)
setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true);
VecAssemblyBegin(nullspaceBasis);
VecAssemblyEnd(nullspaceBasis);
MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE),
1, &nullspaceBasis, &matNullspace);
for (int i = 0; i < nComponents; i++) MatMult(matIntInt, nullspaceBasis, petscSolVec);
setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true); PetscReal n;
VecNorm(petscSolVec, NORM_2, &n);
MSG("NORM IS: %e\n", n);
} else {
MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
}
VecAssemblyBegin(nullspaceBasis); KSPSetNullSpace(kspInterior, matNullspace);
VecAssemblyEnd(nullspaceBasis);
MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullspaceBasis, &matNullSpace); // === Remove null space, if requested. ===
KSPSetNullSpace(kspInterior, matNullSpace);
MatMult(matIntInt, nullspaceBasis, petscSolVec); if (removeRhsNullspace) {
PetscReal n; TEST_EXIT_DBG(coarseSpaceMap == NULL)("Not supported!\n");
VecNorm(petscSolVec, NORM_2, &n);
MSG("NORM IS: %e\n", n); MatNullSpaceRemove(matNullspace, rhsInterior, PETSC_NULL);
}
} else {
TEST_EXIT(removeRhsNullspace == false)
("No nullspace provided that should be removed from rhs!\n");
} }
// PETSc. // PETSc.
solve(rhsInterior, petscSolVec); solve(rhsInterior, petscSolVec);
if (nullspace.size() > 0) { if (nullspace.size() > 0) {
MatNullSpaceDestroy(&matNullSpace); MatNullSpaceDestroy(&matNullspace);
VecDestroy(&nullspaceBasis); VecDestroy(&nullspaceBasis);
} }
......
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