diff --git a/AMDiS/src/parallel/BddcMlSolver.cc b/AMDiS/src/parallel/BddcMlSolver.cc index 682ba3daf383862dde3a2517cae545a89583e3e2..54d69bbcac3cf60ce8de32bb0453e3760e74a4de 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 6cfb6f6b9e2b07f57365b745b796a39c2d68e169..db9e5f58eaf9a1976761e3e3d1c6518d9a575843 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 f6702593b510aed0ddffd4c927e9494775142354..60279806170a309bb3da6dceab4f61edac9d76af 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 3d5983f25c1464363b5b7849b9f238258506a221..e9da7c6f99cb393bb00e98c5116e0458c8fa1aa1 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); }