From f624d7ca08452725fe66303f9815be3527080156 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 12 Jun 2012 05:47:39 +0000
Subject: [PATCH] Added nullspace support.
---
AMDiS/src/parallel/BddcMlSolver.cc | 2 +-
AMDiS/src/parallel/PetscSolver.cc | 6 +-
AMDiS/src/parallel/PetscSolver.h | 8 ++-
AMDiS/src/parallel/PetscSolverGlobalMatrix.cc | 62 ++++++++++---------
4 files changed, 44 insertions(+), 34 deletions(-)
diff --git a/AMDiS/src/parallel/BddcMlSolver.cc b/AMDiS/src/parallel/BddcMlSolver.cc
index 682ba3da..54d69bbc 100644
--- a/AMDiS/src/parallel/BddcMlSolver.cc
+++ b/AMDiS/src/parallel/BddcMlSolver.cc
@@ -323,7 +323,7 @@ namespace AMDiS {
int use_defaults_int = 1;
int parallel_division_int = 1;
- int use_arithmetic_int = 1;
+ int use_arithmetic_int = 0;
int use_adaptive_int = 0;
// MSG("call to \"bddcml_setup_preconditioner\" with the following arguments (each in one line):\n");
// MSG(" %d\n", matrixtype);
diff --git a/AMDiS/src/parallel/PetscSolver.cc b/AMDiS/src/parallel/PetscSolver.cc
index 6cfb6f6b..db9e5f58 100644
--- a/AMDiS/src/parallel/PetscSolver.cc
+++ b/AMDiS/src/parallel/PetscSolver.cc
@@ -27,14 +27,16 @@ namespace AMDiS {
coarseSpaceMap(NULL),
mpiRank(-1),
kspPrefix(""),
- removeRhsNullSpace(false)
+ removeRhsNullspace(false),
+ hasConstantNullspace(false)
{
string kspStr = "";
Parameters::get("parallel->solver->petsc->ksp", kspStr);
if (kspStr != "")
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);
}
diff --git a/AMDiS/src/parallel/PetscSolver.h b/AMDiS/src/parallel/PetscSolver.h
index f6702593..60279806 100644
--- a/AMDiS/src/parallel/PetscSolver.h
+++ b/AMDiS/src/parallel/PetscSolver.h
@@ -116,9 +116,9 @@ namespace AMDiS {
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.
@@ -253,7 +253,9 @@ namespace AMDiS {
/// 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.
- bool removeRhsNullSpace;
+ bool removeRhsNullspace;
+
+ bool hasConstantNullspace;
};
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 3d5983f2..e9da7c6f 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -363,19 +363,7 @@ namespace AMDiS {
// === Remove Dirichlet BC DOFs. ===
// 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 {
VecAssemblyEnd(petscSolVec);
}
- MatNullSpace matNullSpace;
+
+ MatNullSpace matNullspace;
Vec nullspaceBasis;
- if (nullspace.size() > 0) {
- TEST_EXIT_DBG(nullspace.size() == 1)("Not yet implemented!\n");
+ if (nullspace.size() > 0 || hasConstantNullspace) {
+ 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++)
- setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true);
+ MatMult(matIntInt, nullspaceBasis, petscSolVec);
+ PetscReal n;
+ VecNorm(petscSolVec, NORM_2, &n);
+ MSG("NORM IS: %e\n", n);
+ } else {
+ MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
+ }
- VecAssemblyBegin(nullspaceBasis);
- VecAssemblyEnd(nullspaceBasis);
+ KSPSetNullSpace(kspInterior, matNullspace);
- MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullspaceBasis, &matNullSpace);
- KSPSetNullSpace(kspInterior, matNullSpace);
+ // === Remove null space, if requested. ===
- MatMult(matIntInt, nullspaceBasis, petscSolVec);
- PetscReal n;
- VecNorm(petscSolVec, NORM_2, &n);
- MSG("NORM IS: %e\n", n);
+ if (removeRhsNullspace) {
+ TEST_EXIT_DBG(coarseSpaceMap == NULL)("Not supported!\n");
+
+ MatNullSpaceRemove(matNullspace, rhsInterior, PETSC_NULL);
+ }
+ } else {
+ TEST_EXIT(removeRhsNullspace == false)
+ ("No nullspace provided that should be removed from rhs!\n");
}
+
// PETSc.
solve(rhsInterior, petscSolVec);
if (nullspace.size() > 0) {
- MatNullSpaceDestroy(&matNullSpace);
+ MatNullSpaceDestroy(&matNullspace);
VecDestroy(&nullspaceBasis);
}
--
GitLab