PetscSolver.cc 6.67 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20 21


22
#include "AMDiS.h"
23 24 25 26
#include "DOFMatrix.h"
#include "Global.h"
#include "Initfile.h"
#include "parallel/MeshDistributor.h"
27
#include "parallel/MpiHelper.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
#include "parallel/ParallelDofMapping.h"
29 30
#include "parallel/PetscSolver.h"
#include "parallel/StdMpi.h"
31

32
namespace AMDiS { namespace Parallel {
Thomas Witkowski's avatar
Thomas Witkowski committed
33

34 35
  PetscSolver::PetscSolver(std::string name)
    : super(name),
36
      kspPrefix(""),
37
      removeRhsNullspace(false),
38
      hasConstantNullspace(false),
39
      isSymmetric(false),
40
      handleDirichletRows(true)
41
  {
42
    std::string tmp("");
43
    Parameters::get(name + "->petsc prefix", tmp);
44 45
    setKspPrefix(tmp);

46
    std::string kspStr = "";
47
    Parameters::get("parallel->solver->petsc->ksp", kspStr);
48 49
    if (!kspStr.size())
      Parameters::get(name + "->ksp", kspStr);
50

51
    if (kspStr != "")
52
      petsc::options_insert_string(kspStr.c_str());
53

54 55
    Parameters::get(name + "->remove rhs null space", removeRhsNullspace);
    Parameters::get(name + "->has constant null space", hasConstantNullspace);
56
    Parameters::get(name + "->nullspace->const in comp",
57
		    constNullspaceComponent);
58 59 60
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
61
  PetscSolver::~PetscSolver()
62
  {}
Thomas Witkowski's avatar
Thomas Witkowski committed
63 64


65 66
  void PetscSolver::init(std::vector<const FiniteElemSpace*> &fe0,
			 std::vector<const FiniteElemSpace*> &fe1,
Thomas Witkowski's avatar
bo eh  
Thomas Witkowski committed
67
			 bool createGlobalMapping)
68 69 70 71
  {
    FUNCNAME("PetscSolver::init()");

    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
72 73 74
    meshDistributor->setBoundaryDofRequirement(getBoundaryDofRequirement());

    super::init(fe0, fe1, createGlobalMapping);
75 76
  }

77
  /// see seq::LinearSolver::solveLinearSystem(...)
78
  int PetscSolver::solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
79 80 81 82
				    SystemVector& x,
				    SystemVector& b,
				    bool createMatrixData,
				    bool storeMatrixData)
83
  {
84
    FUNCNAME("PetscSolver::solveLinearSystem()");
85

86
    TEST_EXIT(meshDistributor)("No meshDistributor provided. Should not happen!\n");
87

88
    MPI::COMM_WORLD.Barrier();
89
    Timer t;
90

91
    if (createMatrixData)
92 93
      fillPetscMatrix(const_cast< Matrix< DOFMatrix* >* >(A.getOriginalMat()));

94
    fillPetscRhs(&b);
95

96
    INFO(info, 8)("creation of parallel data structures needed %.5f seconds\n",
97
		  t.elapsed());
98

99
    solvePetscMatrix(x, NULL);
100

101 102 103 104
    if (!storeMatrixData) {
      destroyVectorData();
      destroyMatrixData();
    }
105

106
    return getErrorCode();
107
  }
108

109 110 111 112 113 114 115 116
  void PetscSolver::fillPetscMatrix(DOFMatrix* mat)
  {
    Matrix<DOFMatrix*> sysMat(1, 1);
    sysMat[0][0] = mat;
    fillPetscMatrix(&sysMat);
  }


117
  /// see seq::PetscSolver::solve(const MatrixType& A, VectorType& x, const VectorType& b)
118 119 120
  void PetscSolver::solve(Vec &rhs, Vec &sol)
  {
    PetscErrorCode solverError = KSPSolve(kspInterior, rhs, sol);
121

122 123
    PetscInt nIter = 0;
    KSPGetIterationNumber(kspInterior, &nIter);
124
    PetscReal residual_norm = -1.0;
125
    KSPGetResidualNorm(kspInterior, &residual_norm);
126

127
    if (residual_norm <= 0.0) {
128 129 130 131
      Vec r;
      createVec(*interiorMap, r);
      KSPBuildResidual(kspInterior, PETSC_NULL, PETSC_NULL, &r);
      VecNorm(r, NORM_2, &residual_norm);
132
    }
133

134 135 136
    setErrorCode(solverError);
    setIterations(nIter);
    setResidual(residual_norm);
137 138
  }

139

140
  /// was soll das denn genau machen?
141 142 143 144
  void PetscSolver::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolver::solveGlobal()");

145 146 147 148 149 150
    int s, ls;
    VecGetSize(rhs, &s);
    VecGetLocalSize(rhs, &ls);

    MSG("Solve global %d %d\n", ls, s);

151 152 153 154
    ERROR_EXIT("Not implemented!\n");
  }


155
  /// -> Helper-Funktion, benötigt nur domainComm
156
  void PetscSolver::copyVec(Vec& originVec, Vec& destVec,
157
			    std::vector<int>& originIndex, std::vector<int>& destIndex)
158 159
  {
    IS originIs, destIs;
160 161
    ISCreateGeneral(domainComm,
		    originIndex.size(),
162 163 164 165
		    &(originIndex[0]),
		    PETSC_USE_POINTER,
		    &originIs);

166 167
    ISCreateGeneral(domainComm,
		    destIndex.size(),
168 169 170 171 172 173 174 175 176 177 178 179
		    &(destIndex[0]),
		    PETSC_USE_POINTER,
		    &destIs);

    VecScatter scatter;
    VecScatterCreate(originVec, originIs, destVec, destIs, &scatter);
    VecScatterBegin(scatter, originVec, destVec,
		    INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(scatter, originVec, destVec,
		  INSERT_VALUES, SCATTER_FORWARD);

    ISDestroy(&originIs);
180
    ISDestroy(&destIs);
181 182 183
    VecScatterDestroy(&scatter);
  }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
184

185
  /// -> Helper-Funktion
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
186 187 188 189 190 191 192
  bool PetscSolver::testMatrixSymmetric(Mat mat, bool advancedTest)
  {
    FUNCNAME("PetscSolver::testMatrixSymmetric()");

    Mat matTrans;
    MatTranspose(mat, MAT_INITIAL_MATRIX, &matTrans);

Thomas Witkowski's avatar
Thomas Witkowski committed
193 194
    bool isSym = true;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214
    if (advancedTest) {
      int rowStart, rowEnd;
      MatGetOwnershipRange(mat, &rowStart, &rowEnd);

      PetscInt nCols, nColsTrans;
      const PetscInt *cols, *colsTrans;
      const PetscScalar *vals, *valsTrans;

      for (int i = rowStart; i < rowEnd; i++) {
	MatGetRow(mat, i, &nCols, &cols, &vals);
	MatGetRow(matTrans, i, &nColsTrans, &colsTrans, &valsTrans);

	if (nCols != nColsTrans) {
	  MSG("Symmetric test fails: mat row %d: %d nnz   mat col %d: %d nnz\n",
	      i, nCols, i, nColsTrans);
	  isSym = false;
	} else {
	  for (int j = 0; j < nCols; j++) {
	    if (cols[j] != colsTrans[j]) {
	      if (cols[j] > colsTrans[j]) {
215
		MSG("mat[%d][%d] does not exists: mat[%d][%d] = %e\n",
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
216 217 218
		    i, colsTrans[j],
		    colsTrans[j], i, valsTrans[j]);
	      } else {
219
		MSG("mat[%d][%d] does not exists: mat[%d][%d] = %e\n",
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238
		    cols[j], i,
		    i, cols[j], vals[j]);
	      }
	      isSym = false;
	    } else {
	      double n = fabs(vals[j] - valsTrans[j]);
	      if (n > 1e-10) {
		MSG("value diff:  mat[%d][%d] = %e   mat[%d][%d] = %e\n",
		    i, cols[j], vals[j], colsTrans[j], i, valsTrans[j]);

		isSym = false;
	      }
	    }
	  }
	}

	MatRestoreRow(mat, i, &nCols, &cols, &vals);
	MatRestoreRow(matTrans, i, &nColsTrans, &colsTrans, &valsTrans);
      }
239 240
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
241 242 243 244 245
    MatAXPY(matTrans, -1.0, mat, DIFFERENT_NONZERO_PATTERN);
    double norm1, norm2;
    MatNorm(matTrans, NORM_FROBENIUS, &norm1);
    MatNorm(matTrans, NORM_INFINITY, &norm2);
    MatDestroy(&matTrans);
246

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
247
    MSG("Matrix norm test: %e %e\n", norm1, norm2);
248

Thomas Witkowski's avatar
Thomas Witkowski committed
249
    return (isSym && norm1 < 1e-10 && norm2 < 1e-10);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
250
  }
251
} } // end namespaces