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