PetscSolver.cc 6.71 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
#include "parallel/PetscHelper.h"
30 31
#include "parallel/PetscSolver.h"
#include "parallel/StdMpi.h"
32

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

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

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

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

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


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


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

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

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

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

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

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

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

95
    fillPetscRhs(&b);
96

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

100
    solvePetscMatrix(x, NULL);
101

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

107
    return getErrorCode();
108
  }
109

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


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

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

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

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

140

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

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

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

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


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

167 168
    ISCreateGeneral(domainComm,
		    destIndex.size(),
169 170 171 172 173 174 175 176 177 178 179 180
		    &(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);
181
    ISDestroy(&destIs);
182 183 184
    VecScatterDestroy(&scatter);
  }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
185

186
  /// -> Helper-Funktion
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
187 188 189 190 191 192 193
  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
194 195
    bool isSym = true;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
    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]) {
216
		MSG("mat[%d][%d] does not exists: mat[%d][%d] = %e\n",
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
217 218 219
		    i, colsTrans[j],
		    colsTrans[j], i, valsTrans[j]);
	      } else {
220
		MSG("mat[%d][%d] does not exists: mat[%d][%d] = %e\n",
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
		    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);
      }
240 241
    }

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

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

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