6 #include <dune/common/unused.hh> 8 #include <amdis/Initfile.hpp> 9 #include <amdis/linearalgebra/RunnerInterface.hpp> 10 #include <amdis/linearalgebra/SolverInfo.hpp> 11 #include <amdis/linearalgebra/petsc/Interface.hpp> 14 #define KSPDGGMRES "dggmres" 45 template <
class Matrix,
class Vector>
49 using M =
typename Matrix::BaseMatrix;
50 using X =
typename Vector::BaseVector;
51 using Y =
typename Vector::BaseVector;
75 void init(M
const& A)
override 81 KSPCreate(PETSC_COMM_SELF, &ksp_);
84 PETSc::KSPSetOperators(ksp_, A, A);
85 initKSP(ksp_, prefix_);
98 int solve([[maybe_unused]] M
const& A, X& x, Y
const& b,
SolverInfo& solverInfo)
override 100 assert(initialized_);
102 KSPSolve(ksp_, b, x);
105 KSPReasonView(ksp_, PETSC_VIEWER_STDOUT_WORLD);
107 KSPConvergedReason reason;
108 KSPGetConvergedReason(ksp_, &reason);
109 return reason > 0 ? 0 :
int(reason);
112 KSP ksp() {
return ksp_; }
116 virtual void initKSP(KSP ksp, std::string prefix)
const 119 auto kspType = Parameters::get<std::string>(prefix +
"->ksp");
120 std::string kspTypeStr = kspType.value_or(
"default");
125 if (kspTypeStr ==
"direct")
126 initDirectSolver(ksp, prefix);
127 else if (kspTypeStr !=
"default") {
128 KSPSetType(ksp, kspTypeStr.c_str());
131 initKSPParameters(ksp, kspTypeStr.c_str(), prefix);
134 if (kspTypeStr !=
"preonly")
135 KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
142 KSPMonitorSet(ksp, PETSc::KSPMonitorCyclic,
nullptr,
nullptr);
143 }
else if (info > 1) {
144 KSPMonitorSet(ksp, PETSc::KSPMonitorNoisy,
nullptr,
nullptr);
148 auto maxIt = Parameters::get<PetscInt>(prefix +
"->max it");
149 auto rTol = Parameters::get<PetscReal>(prefix +
"->rtol");
150 auto aTol = Parameters::get<PetscReal>(prefix +
"->atol");
151 auto dTol = Parameters::get<PetscReal>(prefix +
"->divtol");
152 if (maxIt || rTol || aTol || dTol)
153 KSPSetTolerances(ksp,
154 rTol.value_or(PETSC_DEFAULT), aTol.value_or(PETSC_DEFAULT),
155 dTol.value_or(PETSC_DEFAULT), maxIt.value_or(PETSC_DEFAULT));
157 auto kspPrefixParam = Parameters::get<std::string>(prefix +
"->prefix");
158 if (kspPrefixParam) {
159 std::string kspPrefix = kspPrefixParam.value() +
"_";
160 KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
161 KSPSetFromOptions(ksp);
166 initPC(pc, prefix +
"->pc");
170 KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD);
176 virtual void initDirectSolver(KSP ksp, std::string prefix)
const 178 KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
179 KSPSetType(ksp, KSPRICHARDSON);
184 auto matSolverType = Parameters::get<std::string>(prefix +
"->mat solver");
186 PETSc::PCFactorSetMatSolverType(pc, matSolverType.value().c_str());
189 KSPGetOperators(ksp, &Amat, &Pmat);
191 MatGetType(Pmat, &type);
192 if (std::strcmp(type, MATSEQAIJ) == 0)
193 PETSc::PCFactorSetMatSolverType(pc, MATSOLVERUMFPACK);
194 else if (std::strcmp(type, MATAIJ) == 0)
195 PETSc::PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
201 virtual void initPC(PC pc, std::string prefix)
const 204 auto pcType = Parameters::get<std::string>(prefix);
205 std::string pcTypeStr =
"default";
207 pcTypeStr = pcType.value();
208 PCSetType(pc, pcTypeStr.c_str());
211 if (pcTypeStr ==
"lu") {
214 auto matSolverType = Parameters::get<std::string>(prefix +
"->mat solver");
216 PETSc::PCFactorSetMatSolverType(pc, matSolverType.value().c_str());
217 }
else if (pcTypeStr ==
"bjacobi") {
219 PetscInt n_local, first;
222 PCBJacobiGetSubKSP(pc, &n_local, &first, &ksps);
223 for (PetscInt i = 0; i < n_local; ++i)
224 initKSP(ksps[i], prefix +
"->sub ksp");
225 }
else if (pcTypeStr ==
"ksp") {
229 PCKSPGetKSP(pc, &ksp);
230 initKSP(ksp, prefix +
"->ksp");
233 auto pcPrefixParam = Parameters::get<std::string>(prefix +
"->prefix");
235 std::string pcPrefix = pcPrefixParam.value() +
"_";
236 PCSetOptionsPrefix(pc, pcPrefix.c_str());
237 PCSetFromOptions(pc);
243 virtual void initKSPParameters(KSP ksp,
char const* ksptype, std::string prefix)
const 246 if (std::strcmp(ksptype, KSPRICHARDSON) == 0)
248 auto scale = Parameters::get<PetscReal>(prefix +
"->scale");
250 KSPRichardsonSetScale(ksp, scale.value());
253 else if (std::strcmp(ksptype, KSPGMRES) == 0 ||
254 std::strcmp(ksptype, KSPFGMRES) == 0 ||
255 std::strcmp(ksptype, KSPLGMRES) == 0 ||
256 std::strcmp(ksptype, KSPDGGMRES) == 0 ||
257 std::strcmp(ksptype, KSPPGMRES) == 0)
259 auto restart = Parameters::get<PetscInt>(prefix +
"->restart");
261 KSPGMRESSetRestart(ksp, restart.value());
263 auto gramschmidt = Parameters::get<std::string>(prefix +
"->gramschmidt");
265 if (gramschmidt.value() ==
"classical") {
266 KSPGMRESSetOrthogonalization(ksp, KSPGMRESClassicalGramSchmidtOrthogonalization);
267 auto refinementType = Parameters::get<std::string>(prefix +
"->refinement type");
268 if (refinementType) {
269 if (refinementType.value() ==
"never")
270 KSPGMRESSetCGSRefinementType(ksp, KSP_GMRES_CGS_REFINE_NEVER);
271 else if (refinementType.value() ==
"ifneeded")
272 KSPGMRESSetCGSRefinementType(ksp, KSP_GMRES_CGS_REFINE_IFNEEDED);
273 else if (refinementType.value() ==
"always")
274 KSPGMRESSetCGSRefinementType(ksp, KSP_GMRES_CGS_REFINE_ALWAYS);
276 }
else if (gramschmidt.value() ==
"modified") {
277 KSPGMRESSetOrthogonalization(ksp, KSPGMRESModifiedGramSchmidtOrthogonalization);
282 else if (std::strcmp(ksptype, KSPBCGSL) == 0)
284 auto ell = Parameters::get<PetscInt>(prefix +
"->ell");
286 KSPBCGSLSetEll(ksp, ell.value());
289 else if (std::strcmp(ksptype, KSPGCR) == 0)
291 auto restart = Parameters::get<PetscInt>(prefix +
"->restart");
293 KSPGCRSetRestart(ksp, restart.value());
302 bool initialized_ =
false;
PetscRunner(std::string const &prefix)
Constructor.
Definition: PetscRunner.hpp:63
Wrapper around PETSc KSP and PC objects to solve a linear system.
Definition: PetscRunner.hpp:46
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
void init(M const &A) override
Implements RunnerInterface::init()
Definition: PetscRunner.hpp:75
int solve([[maybe_unused]] M const &A, X &x, Y const &b, SolverInfo &solverInfo) override
Implements RunnerInterface::solve()
Definition: PetscRunner.hpp:98
static std::optional< T > get(std::string const &key)
Get parameter-values from parameter-tree.
Definition: Initfile.hpp:25
void info(int noInfoLevel, std::string const &str, Args &&... args)
prints a message, if Environment::infoLevel() >= noInfoLevel
Definition: Output.hpp:105
Definition: SolverInfo.hpp:11
void exit() override
Implements RunnerInterface::exit()
Definition: PetscRunner.hpp:90
Interface for Runner / Worker types used in solver classes.
Definition: RunnerInterface.hpp:11
static Dune::MPIHelper::MPICommunicator comm()
Return the MPI_Comm object (or a fake communicator)
Definition: Environment.hpp:86