PetscSolver.hh 6.48 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
23
 ******************************************************************************/


#ifdef HAVE_SEQ_PETSC

24
25
namespace AMDiS {

26
  template<typename M, typename V>
27
  PetscRunner<M,V>::PetscRunner(LinearSolverInterface* oem_)
28
29
30
31
32
33
34
35
    : oem(*oem_),
      kspPrefix(""),
      zeroStartVector(false),
      initialized(false),
      preconditioner(NULL)
  {
    PetscParameters params;
    matSolverPackage = false;
36

37
38
    // set the parameter prefix
    Parameters::get(oem.getName() + "->petsc prefix", kspPrefix); // important when more than one solver to configure
39

40
41
42
    // set the solver
    std::string solverName = "petsc";
    Parameters::get(oem.getName(), solverName);
43
    if (solverName == "petsc")
44
      Parameters::get(oem.getName() + "->ksp_type", solverName);
45

46
    std::string kspSolver = params.solverMap[solverName];
47

48
49
    if (params.matSolverPackage[kspSolver]) {
      // direct solvers
50
51
52
      petsc_options_insert_string(("-" + kspPrefix + "ksp_type preonly").c_str());
      petsc_options_insert_string(("-" + kspPrefix + "pc_type lu").c_str());
      petsc_options_insert_string(("-" + kspPrefix + "pc_factor_mat_solver_package " + (kspSolver != "direct" ? kspSolver : "umfpack")).c_str());
53
54
55
      oem.setMaxIterations(1);
      zeroStartVector = true;
      matSolverPackage = true;
56
    } else if (!params.emptyParam[kspSolver]) {
57
      // other solvers
58
      petsc_options_insert_string(("-" + kspPrefix + "ksp_type " + kspSolver).c_str());
59
    }
60

61
62
    // set the preconditioner
    setPrecon();
63

64
    if (oem.getInfo() >= 20)
65
      petsc_options_insert_string(("-" + kspPrefix + "ksp_monitor_true_residual").c_str());
66
    else if (oem.getInfo() >= 10)
67
      petsc_options_insert_string(("-" + kspPrefix + "ksp_monitor").c_str());
68

69
70
71
72
    // command line string
    std::string kspString = "";
    Parameters::get(oem.getName() + "->ksp", kspString);
    if (kspString != "")
73
      petsc_options_insert_string(kspString.c_str());
74
75
  }

76

77
78
79
80
81
  template<typename MatrixType, typename VectorType>
  void PetscRunner<MatrixType, VectorType>::init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix)
  {
    if (initialized)
      exit();
82

83
    createSubSolver(ksp, fullMatrix.matrix, kspPrefix);
84

85
    // set tolerance options from LinearSolverInterface-parameters
86
87
88
89
90
    KSPSetTolerances(ksp, oem.getRelative(), oem.getTolerance(), PETSC_DEFAULT, oem.getMaxIterations());

    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
91
92
93
94

    KSPGetPC(ksp, &pc);

    preconditioner->init(pc, A, fullMatrix);
95
96
97
    initialized = true;
  }

98

99
  template<typename MatrixType, typename VectorType>
100
  int PetscRunner<MatrixType, VectorType>::solve(const MatrixType& A, VectorType& x, const VectorType& b)
101
102
103
  {
    if (zeroStartVector)
      VecSet(x.vector, 0.0);
104

105
    PetscErrorCode solverError = KSPSolve(ksp, b.vector, x.vector);
106

107
108
109
110
    PetscInt nIter = 0;
    KSPGetIterationNumber(ksp, &nIter);
    PetscReal residual_norm = -1.0;
    KSPGetResidualNorm(ksp, &residual_norm);
111

112
113
114
115
116
117
    if (residual_norm <= 0.0) {
      Vec r;
      VecDuplicate(b.vector, &r);
      KSPBuildResidual(ksp, PETSC_NULL, PETSC_NULL, &r);
      VecNorm(r, NORM_2, &residual_norm);
    }
118

119
120
121
    oem.setErrorCode(solverError);
    oem.setIterations(nIter);
    oem.setResidual(residual_norm);
122

123
124
125
    return solverError;
  }

126

127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
  template<typename M, typename V>
  void PetscRunner<M,V>::createSubSolver(KSP &ksp_, Mat m, std::string kspPrefix_)
  {
    KSPCreate(PETSC_COMM_SELF, &ksp_);
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(ksp_, m, m);
#else
    KSPSetOperators(ksp_, m, m, SAME_NONZERO_PATTERN);
#endif
    KSPSetOptionsPrefix(ksp_, kspPrefix_.c_str());
    KSPSetFromOptions(ksp_);
  }


  template<typename M, typename V>
  void PetscRunner<M,V>::setSolver(KSP ksp_, std::string kspPrefix_,
143
				   KSPType kspType,  PCType pcType,
144
145
146
147
148
149
150
151
152
153
154
155
156
				   PetscReal rtol, PetscReal atol, PetscInt maxIt)
  {
    KSPSetType(ksp_, kspType);
    KSPSetTolerances(ksp_, rtol, atol, PETSC_DEFAULT, maxIt);
    if (kspPrefix_ != "")
      KSPSetOptionsPrefix(ksp_, kspPrefix_.c_str());
    KSPSetFromOptions(ksp_);

    PC pc_;
    KSPGetPC(ksp_, &pc_);
    PCSetType(pc_, pcType);
    PCSetFromOptions(pc_);
  }
157
158


159
160
161
162
163
164
165
166
167
  template<typename MatrixType, typename VectorType>
  void PetscRunner<MatrixType, VectorType>::setPrecon()
  {
    std::string precon = "";
    std::string initFileStr = oem.getName() + "->pc_type";
    Parameters::get(initFileStr, precon);
    if (!precon.size() || precon == "no" || precon == "0") {
      initFileStr = oem.getName() + "->left precon";
      Parameters::get(initFileStr, precon);
168
    }
169
170
171
172
173
174
175
    if (!precon.size() || precon == "no" || precon == "0") {
      initFileStr = oem.getName() + "->right precon";
      Parameters::get(initFileStr, precon);
    }
    // no preconditioner should be created
    if (!precon.size() || precon == "no" || precon == "0")
      return;
176

177
178
179
180
181
    PetscParameters params;
    if (preconditioner) {
      delete preconditioner;
      preconditioner = NULL;
    }
182

183
    // Creator for the left preconditioner
184
    CreatorInterfaceName< PetscPreconditionerInterface<MatrixType, VectorType> >* creator(NULL);
185
186

    try {
187
188
      creator = dynamic_cast<CreatorInterfaceName< PetscPreconditionerInterface<MatrixType, VectorType> >*>(
	CreatorMap<PetscPreconditionerInterface<MatrixType, VectorType> >::getCreator(precon, initFileStr) );
189
190
191
    } catch (...) {
      creator = NULL;
    }
192

193
194
195
196
197
    if (creator) {
      creator->setName(initFileStr);
      preconditioner = creator->create();
    } else if (!matSolverPackage && !params.emptyParam[precon]) {
      precon = (params.preconMap.find(precon) != params.preconMap.end() ? params.preconMap[precon] : precon);
198
      preconditioner = new PetscPreconditionerInterface<MatrixType, VectorType>(kspPrefix, precon);
199
200
201
202
203
204
205
    } else {
      ERROR_EXIT((std::string("There is no creator for the given preconditioner '") + precon + "'").c_str());
    }
  }
}

#endif // HAVE_SEQ_PETSC