Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "PetscSupport.h"
11 :
12 : // MOOSE includes
13 : #include "MooseApp.h"
14 : #include "FEProblem.h"
15 : #include "DisplacedProblem.h"
16 : #include "NonlinearSystem.h"
17 : #include "LinearSystem.h"
18 : #include "AuxiliarySystem.h"
19 : #include "DisplacedProblem.h"
20 : #include "PenetrationLocator.h"
21 : #include "NearestNodeLocator.h"
22 : #include "MooseTypes.h"
23 : #include "MooseUtils.h"
24 : #include "CommandLine.h"
25 : #include "Console.h"
26 : #include "MultiMooseEnum.h"
27 : #include "Conversion.h"
28 : #include "Executioner.h"
29 : #include "MooseMesh.h"
30 : #include "ComputeLineSearchObjectWrapper.h"
31 : #include "Convergence.h"
32 : #include "ParallelParamObject.h"
33 :
34 : #include "libmesh/equation_systems.h"
35 : #include "libmesh/linear_implicit_system.h"
36 : #include "libmesh/nonlinear_implicit_system.h"
37 : #include "libmesh/petsc_linear_solver.h"
38 : #include "libmesh/petsc_matrix.h"
39 : #include "libmesh/petsc_nonlinear_solver.h"
40 : #include "libmesh/petsc_preconditioner.h"
41 : #include "libmesh/petsc_vector.h"
42 : #include "libmesh/sparse_matrix.h"
43 : #include "libmesh/petsc_solver_exception.h"
44 : #include "libmesh/simple_range.h"
45 :
46 : // PETSc includes
47 : #include <petsc.h>
48 : #include <petscsnes.h>
49 : #include <petscksp.h>
50 : #include <petscmat.h>
51 : #include <petscis.h>
52 : #include <petscdm.h>
53 : #include <petscoptions.h>
54 :
55 : // PetscDMMoose include
56 : #include "PetscDMMoose.h"
57 :
58 : // Standard includes
59 : #include <ostream>
60 : #include <cctype>
61 : #include <fstream>
62 : #include <optional>
63 : #include <string>
64 :
65 : using namespace libMesh;
66 :
67 : void
68 0 : MooseVecView(NumericVector<Number> & vector)
69 : {
70 0 : PetscVector<Number> & petsc_vec = static_cast<PetscVector<Number> &>(vector);
71 0 : LibmeshPetscCallA(vector.comm().get(), VecView(petsc_vec.vec(), 0));
72 0 : }
73 :
74 : void
75 0 : MooseMatView(SparseMatrix<Number> & mat)
76 : {
77 0 : PetscMatrixBase<Number> & petsc_mat = static_cast<PetscMatrix<Number> &>(mat);
78 0 : LibmeshPetscCallA(mat.comm().get(), MatView(petsc_mat.mat(), 0));
79 0 : }
80 :
81 : void
82 0 : MooseVecView(const NumericVector<Number> & vector)
83 : {
84 0 : PetscVector<Number> & petsc_vec =
85 : static_cast<PetscVector<Number> &>(const_cast<NumericVector<Number> &>(vector));
86 0 : LibmeshPetscCallA(vector.comm().get(), VecView(petsc_vec.vec(), 0));
87 0 : }
88 :
89 : void
90 0 : MooseMatView(const SparseMatrix<Number> & mat)
91 : {
92 0 : PetscMatrixBase<Number> & petsc_mat =
93 : static_cast<PetscMatrix<Number> &>(const_cast<SparseMatrix<Number> &>(mat));
94 0 : LibmeshPetscCallA(mat.comm().get(), MatView(petsc_mat.mat(), 0));
95 0 : }
96 :
97 : namespace Moose
98 : {
99 : namespace PetscSupport
100 : {
101 :
102 : namespace
103 : {
104 :
105 : void
106 64 : applySystemVectorTypeOptions(FEProblemBase & problem, libMesh::System & lm_sys)
107 : {
108 228 : for (auto & [_, vec] : as_range(lm_sys.vectors_begin(), lm_sys.vectors_end()))
109 : {
110 164 : auto * const petsc_vec = cast_ptr<PetscVector<Number> *>(vec.get());
111 164 : LibmeshPetscCallA(problem.comm().get(), VecSetFromOptions(petsc_vec->vec()));
112 : }
113 :
114 : // The solution vectors aren't included in the system vectors storage.
115 64 : auto * petsc_vec = cast_ptr<PetscVector<Number> *>(lm_sys.solution.get());
116 64 : LibmeshPetscCallA(problem.comm().get(), VecSetFromOptions(petsc_vec->vec()));
117 64 : petsc_vec = cast_ptr<PetscVector<Number> *>(lm_sys.current_local_solution.get());
118 64 : LibmeshPetscCallA(problem.comm().get(), VecSetFromOptions(petsc_vec->vec()));
119 64 : }
120 :
121 : void
122 32 : applyVectorTypeOptions(FEProblemBase & problem)
123 : {
124 64 : for (const auto sys_index : make_range(problem.numSolverSystems()))
125 32 : applySystemVectorTypeOptions(problem, problem.getSolverSystem(sys_index).system());
126 :
127 32 : applySystemVectorTypeOptions(problem, problem.getAuxiliarySystem().system());
128 32 : }
129 :
130 : bool
131 340443 : petscOptionsHasName(::PetscOptions options,
132 : const std::string & name,
133 : const std::string & prefix = "")
134 : {
135 340443 : PetscBool found = PETSC_FALSE;
136 340443 : const char * const prefix_ptr = prefix.empty() ? nullptr : prefix.c_str();
137 340443 : LibmeshPetscCallA(PETSC_COMM_WORLD,
138 : PetscOptionsHasName(options, prefix_ptr, name.c_str(), &found));
139 340443 : return found;
140 : }
141 :
142 : bool
143 3 : hasMatrixFreeSolveType(const FEProblemBase & problem)
144 : {
145 6 : for (const auto sys_index : make_range(problem.numSolverSystems()))
146 3 : if (const auto solve_type = problem.solverParams(sys_index)._type;
147 3 : solve_type == Moose::ST_JFNK || solve_type == Moose::ST_PJFNK)
148 0 : return true;
149 :
150 3 : return false;
151 : }
152 :
153 : bool
154 127862 : mightBeMatTypeOption(const std::string & name)
155 : {
156 : static constexpr std::string_view mat_type_suffix = "mat_type";
157 255724 : return name.size() >= mat_type_suffix.size() &&
158 127862 : std::equal(mat_type_suffix.rbegin(),
159 127862 : mat_type_suffix.rend(),
160 127862 : name.rbegin(),
161 353658 : [](const char left, const char right)
162 : // tolower requires representability by unsigned char
163 : {
164 353658 : return static_cast<int>(left) ==
165 353658 : std::tolower(libMesh::cast_int<unsigned char>(right));
166 127862 : });
167 : }
168 :
169 : void
170 40465 : errorOnUnprefixedMatTypeOption(::PetscOptions options, FEProblemBase & problem)
171 : {
172 161860 : if (!petscOptionsHasName(options, "-mat_type"))
173 40462 : return;
174 :
175 : std::string error_string =
176 : "Setting option '-mat_type' is not supported without a solver-system prefix. Use an option "
177 3 : "such as '-" +
178 3 : problem.getSolverSystem(0).name() + "_mat_type' for assembled libMesh matrices.";
179 3 : if (hasMatrixFreeSolveType(problem))
180 : error_string +=
181 : " Attempting to change the matrix "
182 : "type for the MFFD matrix type used to represent the Jacobian for (P)JFNK solve "
183 0 : "types is not supported.";
184 3 : mooseError(error_string);
185 0 : }
186 :
187 : // Allow iterating over all systems or allowing caller to specify a specific system for which to
188 : // apply matrix type options
189 : void
190 363384 : applyMatrixTypeOptions(FEProblemBase & problem,
191 : const std::optional<std::size_t> system_index = std::nullopt)
192 : {
193 363384 : const auto begin = system_index.value_or(0);
194 363384 : const auto end = system_index ? begin + 1 : problem.numSolverSystems();
195 :
196 739308 : for (const auto sys_index : make_range(begin, end))
197 : {
198 375924 : auto & solver_system = problem.getSolverSystem(sys_index);
199 375924 : auto & lm_sys = solver_system.system();
200 :
201 : // Even in matrix-free modes the libMesh matrix wrappers can exist before the PETSc Mat does.
202 375924 : if (problem.solverParams(sys_index)._type == Moose::ST_JFNK)
203 149 : continue;
204 :
205 752804 : for (auto & [_, mat] : as_range(lm_sys.matrices_begin(), lm_sys.matrices_end()))
206 377029 : if (auto * const petsc_mat = dynamic_cast<PetscMatrixBase<Number> *>(mat.get()); petsc_mat)
207 : {
208 377029 : LibmeshPetscCallA(
209 : problem.comm().get(),
210 : MatSetOptionsPrefix(petsc_mat->mat(), (solver_system.name() + "_").c_str()));
211 377029 : LibmeshPetscCallA(problem.comm().get(), MatSetFromOptions(petsc_mat->mat()));
212 : }
213 : }
214 363384 : }
215 :
216 : } // namespace
217 :
218 : std::string
219 2490 : stringify(const LineSearchType & t)
220 : {
221 2490 : switch (t)
222 : {
223 2407 : case LS_BASIC:
224 4814 : return "basic";
225 0 : case LS_DEFAULT:
226 0 : return "default";
227 0 : case LS_NONE:
228 0 : return "none";
229 0 : case LS_SHELL:
230 0 : return "shell";
231 0 : case LS_L2:
232 0 : return "l2";
233 21 : case LS_BT:
234 42 : return "bt";
235 62 : case LS_CP:
236 124 : return "cp";
237 0 : case LS_CONTACT:
238 0 : return "contact";
239 0 : case LS_PROJECT:
240 0 : return "project";
241 0 : case LS_INVALID:
242 0 : mooseError("Invalid LineSearchType");
243 : }
244 0 : return "";
245 : }
246 :
247 : std::string
248 28074 : stringify(const MffdType & t)
249 : {
250 28074 : switch (t)
251 : {
252 28065 : case MFFD_WP:
253 56130 : return "wp";
254 9 : case MFFD_DS:
255 18 : return "ds";
256 0 : case MFFD_INVALID:
257 0 : mooseError("Invalid MffdType");
258 : }
259 0 : return "";
260 : }
261 :
262 : void
263 40651 : setSolverOptions(const SolverParams & solver_params, const MultiMooseEnum & dont_add_these_options)
264 : {
265 40651 : const auto prefix_with_dash = '-' + solver_params._prefix;
266 : // set PETSc options implied by a solve type
267 40651 : switch (solver_params._type)
268 : {
269 27958 : case Moose::ST_PJFNK:
270 55916 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
271 55916 : prefix_with_dash + "snes_mf_operator");
272 27958 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
273 55916 : prefix_with_dash + "mat_mffd_type",
274 55916 : stringify(solver_params._mffd_type));
275 27958 : break;
276 :
277 116 : case Moose::ST_JFNK:
278 116 : setSinglePetscOptionIfAppropriate(dont_add_these_options, prefix_with_dash + "snes_mf");
279 116 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
280 232 : prefix_with_dash + "mat_mffd_type",
281 232 : stringify(solver_params._mffd_type));
282 116 : break;
283 :
284 10539 : case Moose::ST_NEWTON:
285 10539 : break;
286 :
287 0 : case Moose::ST_FD:
288 0 : setSinglePetscOptionIfAppropriate(dont_add_these_options, prefix_with_dash + "snes_fd");
289 0 : break;
290 :
291 2038 : case Moose::ST_LINEAR:
292 6114 : setSinglePetscOptionIfAppropriate(
293 4076 : dont_add_these_options, prefix_with_dash + "snes_type", "ksponly");
294 4076 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
295 4076 : prefix_with_dash + "snes_monitor_cancel");
296 2038 : break;
297 : }
298 :
299 40651 : Moose::LineSearchType ls_type = solver_params._line_search;
300 40651 : if (ls_type == Moose::LS_NONE)
301 2341 : ls_type = Moose::LS_BASIC;
302 :
303 40651 : if (ls_type != Moose::LS_DEFAULT && ls_type != Moose::LS_CONTACT && ls_type != Moose::LS_PROJECT)
304 2490 : setSinglePetscOptionIfAppropriate(
305 4980 : dont_add_these_options, prefix_with_dash + "snes_linesearch_type", stringify(ls_type));
306 40651 : }
307 :
308 : void
309 40983 : addPetscOptionsFromCommandline(FEProblemBase * const problem)
310 : {
311 : // commandline options always win
312 : // the options from a user commandline will overwrite the existing ones if any conflicts
313 : int argc;
314 : char ** args;
315 :
316 40983 : LibmeshPetscCallA(PETSC_COMM_WORLD, PetscGetArgs(&argc, &args));
317 40983 : std::vector<const char *> cl_args(args + 1, args + argc);
318 40983 : const auto cl_argc = libMesh::cast_int<int>(cl_args.size());
319 :
320 : ::PetscOptions command_line_options;
321 40983 : LibmeshPetscCallA(PETSC_COMM_WORLD, PetscOptionsCreate(&command_line_options));
322 40983 : LibmeshPetscCallA(PETSC_COMM_WORLD,
323 : PetscOptionsInsertArgs(command_line_options, cl_argc, cl_args.data()));
324 40983 : LibmeshPetscCallA(PETSC_COMM_WORLD,
325 : PetscOptionsInsertArgs(LIBMESH_PETSC_NULLPTR, cl_argc, cl_args.data()));
326 :
327 40983 : if (!problem)
328 : {
329 541 : LibmeshPetscCallA(PETSC_COMM_WORLD, PetscOptionsDestroy(&command_line_options));
330 541 : return;
331 : }
332 :
333 40442 : errorOnUnprefixedMatTypeOption(command_line_options, *problem);
334 :
335 : // Some vector/matrix-type options may have been consumed before the PETSc database rebuild.
336 : // Replay only the command-line-controlled applications so input-file options handled through
337 : // setSinglePetscOption() do not pay the cost twice.
338 121317 : const bool have_vec_type = petscOptionsHasName(command_line_options, "-vec_type");
339 40439 : bool have_mat_type = false;
340 :
341 81088 : for (const auto sys_index : make_range(problem->numSolverSystems()))
342 : {
343 40658 : have_mat_type = petscOptionsHasName(
344 81316 : command_line_options, "-mat_type", problem->getSolverSystem(sys_index).name() + "_");
345 40658 : if (have_mat_type)
346 9 : break;
347 : }
348 :
349 40439 : if (have_vec_type)
350 9 : applyVectorTypeOptions(*problem);
351 40439 : if (have_mat_type)
352 9 : applyMatrixTypeOptions(*problem);
353 :
354 40439 : LibmeshPetscCallA(PETSC_COMM_WORLD, PetscOptionsDestroy(&command_line_options));
355 40980 : }
356 :
357 : void
358 40442 : petscSetOptionsHelper(const PetscOptions & po, FEProblemBase * const problem)
359 : {
360 : // Add any additional options specified in the input file
361 42684 : for (const auto & flag : po.flags)
362 : // Need to use name method here to pass a str instead of an EnumItem because
363 : // we don't care if the id attributes match
364 2242 : if (!po.dont_add_these_options.contains(flag.name()) ||
365 0 : po.user_set_options.contains(flag.name()))
366 8968 : setSinglePetscOption(flag.rawName().c_str());
367 :
368 : // Add option pairs
369 169464 : for (auto & option : po.pairs)
370 130198 : if (!po.dont_add_these_options.contains(option.first) ||
371 1176 : po.user_set_options.contains(option.first))
372 127855 : setSinglePetscOption(option.first, option.second, problem);
373 :
374 40442 : addPetscOptionsFromCommandline(problem);
375 40439 : }
376 :
377 : void
378 1483 : petscSetOptions(const PetscOptions & po,
379 : const SolverParams & solver_params,
380 : FEProblemBase * const problem)
381 : {
382 1483 : PetscCallAbort(PETSC_COMM_WORLD, PetscOptionsClear(LIBMESH_PETSC_NULLPTR));
383 1483 : setSolverOptions(solver_params, po.dont_add_these_options);
384 1483 : petscSetOptionsHelper(po, problem);
385 1483 : }
386 :
387 : void
388 38959 : petscSetOptions(const PetscOptions & po,
389 : const std::vector<SolverParams> & solver_params_vec,
390 : FEProblemBase * const problem)
391 : {
392 38959 : PetscCallAbort(PETSC_COMM_WORLD, PetscOptionsClear(LIBMESH_PETSC_NULLPTR));
393 78127 : for (const auto & solver_params : solver_params_vec)
394 39168 : setSolverOptions(solver_params, po.dont_add_these_options);
395 38959 : petscSetOptionsHelper(po, problem);
396 38956 : }
397 :
398 : PetscErrorCode
399 47663 : petscSetupOutput(CommandLine * cmd_line)
400 : {
401 : PetscFunctionBegin;
402 47663 : char code[10] = {45, 45, 109, 111, 111, 115, 101};
403 47663 : const std::vector<std::string> argv = cmd_line->getArguments();
404 563064 : for (const auto & arg : argv)
405 : {
406 515401 : if (arg.compare(code) == 0)
407 : {
408 0 : Console::petscSetupOutput();
409 0 : break;
410 : }
411 : }
412 47663 : PetscFunctionReturn(PETSC_SUCCESS);
413 47663 : }
414 :
415 : PetscErrorCode
416 747622 : petscNonlinearConverged(SNES /*snes*/,
417 : PetscInt it,
418 : PetscReal /*xnorm*/,
419 : PetscReal /*snorm*/,
420 : PetscReal /*fnorm*/,
421 : SNESConvergedReason * reason,
422 : void * ctx)
423 : {
424 : PetscFunctionBegin;
425 747622 : FEProblemBase & problem = *static_cast<FEProblemBase *>(ctx);
426 :
427 : // execute objects that may be used in convergence check
428 747622 : problem.execute(EXEC_NONLINEAR_CONVERGENCE);
429 :
430 : // perform the convergence check
431 : Convergence::MooseConvergenceStatus status;
432 747622 : if (problem.getFailNextNonlinearConvergenceCheck())
433 : {
434 221 : status = Convergence::MooseConvergenceStatus::DIVERGED;
435 221 : problem.resetFailNextNonlinearConvergenceCheck();
436 : }
437 : else
438 : {
439 747401 : auto & convergence = problem.getConvergence(
440 747401 : problem.getNonlinearConvergenceNames()[problem.currentNonlinearSystem().number()]);
441 747401 : status = convergence.checkConvergence(it);
442 : }
443 :
444 : // convert convergence status to PETSc converged reason
445 747619 : switch (status)
446 : {
447 462783 : case Convergence::MooseConvergenceStatus::ITERATING:
448 462783 : *reason = SNES_CONVERGED_ITERATING;
449 462783 : break;
450 :
451 284568 : case Convergence::MooseConvergenceStatus::CONVERGED:
452 284568 : *reason = SNES_CONVERGED_FNORM_ABS;
453 284568 : break;
454 :
455 268 : case Convergence::MooseConvergenceStatus::DIVERGED:
456 268 : *reason = SNES_DIVERGED_DTOL;
457 268 : break;
458 : }
459 :
460 747619 : PetscFunctionReturn(PETSC_SUCCESS);
461 : }
462 :
463 : PetscErrorCode
464 4635 : petscLinearConverged(
465 : KSP /*ksp*/, PetscInt it, PetscReal /*norm*/, KSPConvergedReason * reason, void * ctx)
466 : {
467 : PetscFunctionBegin;
468 4635 : FEProblemBase & problem = *static_cast<FEProblemBase *>(ctx);
469 :
470 : // execute objects that may be used in convergence check
471 : // Right now, setting objects to execute on this flag would be ignored except in the
472 : // linear-system-only use case.
473 4635 : problem.execute(EXEC_LINEAR_CONVERGENCE);
474 :
475 : // perform the convergence check
476 : Convergence::MooseConvergenceStatus status;
477 4635 : if (problem.getFailNextSystemConvergenceCheck())
478 : {
479 0 : status = Convergence::MooseConvergenceStatus::DIVERGED;
480 0 : problem.resetFailNextSystemConvergenceCheck();
481 : }
482 : else
483 : {
484 4635 : auto & convergence = problem.getConvergence(
485 4635 : problem.getLinearConvergenceNames()[problem.currentLinearSystem().number()]);
486 4635 : status = convergence.checkConvergence(it);
487 : }
488 :
489 : // convert convergence status to PETSc converged reason
490 4635 : switch (status)
491 : {
492 4113 : case Convergence::MooseConvergenceStatus::ITERATING:
493 4113 : *reason = KSP_CONVERGED_ITERATING;
494 4113 : break;
495 :
496 : // TODO: find a KSP code that works better for this case
497 522 : case Convergence::MooseConvergenceStatus::CONVERGED:
498 : #if PETSC_VERSION_LESS_THAN(3, 24, 0)
499 : *reason = KSP_CONVERGED_RTOL_NORMAL;
500 : #else
501 522 : *reason = KSP_CONVERGED_RTOL_NORMAL_EQUATIONS;
502 : #endif
503 522 : break;
504 :
505 0 : case Convergence::MooseConvergenceStatus::DIVERGED:
506 0 : *reason = KSP_DIVERGED_DTOL;
507 0 : break;
508 : }
509 :
510 4635 : PetscFunctionReturn(PETSC_SUCCESS);
511 : }
512 :
513 : PCSide
514 0 : getPetscPCSide(Moose::PCSideType pcs)
515 : {
516 0 : switch (pcs)
517 : {
518 0 : case Moose::PCS_LEFT:
519 0 : return PC_LEFT;
520 0 : case Moose::PCS_RIGHT:
521 0 : return PC_RIGHT;
522 0 : case Moose::PCS_SYMMETRIC:
523 0 : return PC_SYMMETRIC;
524 0 : default:
525 0 : mooseError("Unknown PC side requested.");
526 : break;
527 : }
528 : }
529 :
530 : KSPNormType
531 375285 : getPetscKSPNormType(Moose::MooseKSPNormType kspnorm)
532 : {
533 375285 : switch (kspnorm)
534 : {
535 0 : case Moose::KSPN_NONE:
536 0 : return KSP_NORM_NONE;
537 0 : case Moose::KSPN_PRECONDITIONED:
538 0 : return KSP_NORM_PRECONDITIONED;
539 375285 : case Moose::KSPN_UNPRECONDITIONED:
540 375285 : return KSP_NORM_UNPRECONDITIONED;
541 0 : case Moose::KSPN_NATURAL:
542 0 : return KSP_NORM_NATURAL;
543 0 : case Moose::KSPN_DEFAULT:
544 0 : return KSP_NORM_DEFAULT;
545 0 : default:
546 0 : mooseError("Unknown KSP norm type requested.");
547 : break;
548 : }
549 : }
550 :
551 : void
552 350937 : petscSetDefaultKSPNormType(FEProblemBase & problem, KSP ksp)
553 : {
554 726222 : for (const auto i : make_range(problem.numSolverSystems()))
555 : {
556 375285 : SolverSystem & sys = problem.getSolverSystem(i);
557 375285 : LibmeshPetscCallA(problem.comm().get(),
558 : KSPSetNormType(ksp, getPetscKSPNormType(sys.getMooseKSPNormType())));
559 : }
560 350937 : }
561 :
562 : void
563 350937 : petscSetDefaultPCSide(FEProblemBase & problem, KSP ksp)
564 : {
565 726222 : for (const auto i : make_range(problem.numSolverSystems()))
566 : {
567 375285 : SolverSystem & sys = problem.getSolverSystem(i);
568 :
569 : // PETSc 3.2.x+
570 375285 : if (sys.getPCSide() != Moose::PCS_DEFAULT)
571 0 : LibmeshPetscCallA(problem.comm().get(), KSPSetPCSide(ksp, getPetscPCSide(sys.getPCSide())));
572 : }
573 350937 : }
574 :
575 : void
576 349748 : petscSetKSPDefaults(FEProblemBase & problem, KSP ksp)
577 : {
578 349748 : auto & es = problem.es();
579 :
580 349748 : PetscReal rtol = es.parameters.get<Real>("linear solver tolerance");
581 349748 : PetscReal atol = es.parameters.get<Real>("linear solver absolute tolerance");
582 :
583 : // MOOSE defaults this to -1 for some dumb reason
584 349748 : if (atol < 0)
585 0 : atol = 1e-50;
586 :
587 349748 : PetscReal maxits = es.parameters.get<unsigned int>("linear solver maximum iterations");
588 :
589 : // 1e100 is because we don't use divtol currently
590 349748 : LibmeshPetscCallA(problem.comm().get(), KSPSetTolerances(ksp, rtol, atol, 1e100, maxits));
591 :
592 349748 : petscSetDefaultPCSide(problem, ksp);
593 :
594 349748 : petscSetDefaultKSPNormType(problem, ksp);
595 349748 : }
596 :
597 : void
598 363352 : petscSetDefaults(FEProblemBase & problem)
599 : {
600 : // Apply matrix-type options once the per-system matrix prefixes are known. This is different
601 : // from vectors: libMesh/PETSc vector construction already sees a global '-vec_type' option,
602 : // but prefixed matrix options such as '-nl0_mat_type' cannot match anything until we set the
603 : // matrix prefix here. Without this, a matrix may be constructed with the default type and keep
604 : // that type for the rest of the solve, unless we not only set the options prefix but also apply
605 : // the options to the matrix in this function call.
606 363352 : applyMatrixTypeOptions(problem);
607 :
608 712792 : for (const auto nl_index : make_range(problem.numNonlinearSystems()))
609 : {
610 349440 : NonlinearSystemBase & nl = problem.getNonlinearSystemBase(nl_index);
611 :
612 : // dig out PETSc solver
613 349440 : auto * const petsc_solver = cast_ptr<PetscNonlinearSolver<Number> *>(nl.nonlinearSolver());
614 :
615 : // Ensure we properly prefix SNES which in turn prefixes its KSP
616 349440 : const char * snes_prefix = nullptr;
617 349440 : std::string snes_prefix_str;
618 349440 : if (nl.system().prefix_with_name())
619 : {
620 24348 : snes_prefix_str = nl.system().prefix();
621 24348 : snes_prefix = snes_prefix_str.c_str();
622 : }
623 349440 : SNES snes = petsc_solver->snes(snes_prefix);
624 : KSP ksp;
625 349440 : LibmeshPetscCallA(nl.comm().get(), SNESGetKSP(snes, &ksp));
626 349440 : LibmeshPetscCallA(nl.comm().get(), SNESSetMaxLinearSolveFailures(snes, 1000000));
627 349440 : LibmeshPetscCallA(nl.comm().get(), SNESSetCheckJacobianDomainError(snes, PETSC_TRUE));
628 349440 : LibmeshPetscCallA(
629 : nl.comm().get(),
630 : SNESSetConvergenceTest(snes, petscNonlinearConverged, &problem, LIBMESH_PETSC_NULLPTR));
631 :
632 349440 : petscSetKSPDefaults(problem, ksp);
633 349440 : }
634 :
635 389804 : for (auto sys_index : make_range(problem.numLinearSystems()))
636 : {
637 : // dig out PETSc solver
638 26452 : LinearSystem & lin_sys = problem.getLinearSystem(sys_index);
639 26452 : auto & lm_lin_sys = lin_sys.linearImplicitSystem();
640 : auto * const petsc_solver =
641 26452 : dynamic_cast<PetscLinearSolver<Number> *>(lm_lin_sys.get_linear_solver());
642 : // Ensure we properly prefix KSP
643 26452 : if (lm_lin_sys.prefix_with_name())
644 732 : petsc_solver->init(lm_lin_sys.prefix().c_str());
645 : else
646 25720 : petsc_solver->init();
647 : // The KSP call here would initialize without a prefix if we hadn't "manually" performed
648 : // initialization above
649 26452 : KSP ksp = petsc_solver->ksp();
650 :
651 26452 : if (problem.hasLinearConvergenceObjects())
652 531 : LibmeshPetscCallA(
653 : lin_sys.comm().get(),
654 : KSPSetConvergenceTest(ksp, petscLinearConverged, &problem, LIBMESH_PETSC_NULLPTR));
655 :
656 : // We dont set the KSP defaults here because they seem to clash with the linear solve parameters
657 : // set in FEProblemBase::solveLinearSystem
658 : }
659 363352 : }
660 :
661 : void
662 73278 : processSingletonMooseWrappedOptions(FEProblemBase & fe_problem, const InputParameters & params)
663 : {
664 73278 : setSolveTypeFromParams(fe_problem, params);
665 73278 : setLineSearchFromParams(fe_problem, params);
666 73278 : setMFFDTypeFromParams(fe_problem, params);
667 73278 : }
668 :
669 : #define checkPrefix(prefix) \
670 : mooseAssert(prefix[0] == '-', \
671 : "Leading prefix character must be a '-'. Current prefix is '" << prefix << "'"); \
672 : mooseAssert((prefix.size() == 1) || (prefix.back() == '_'), \
673 : "Terminating prefix character must be a '_'. Current prefix is '" << prefix << "'"); \
674 : mooseAssert(MooseUtils::isAllLowercase(prefix), "PETSc prefixes should be all lower-case")
675 :
676 : void
677 73278 : storePetscOptions(FEProblemBase & fe_problem,
678 : const std::string & prefix,
679 : const ParallelParamObject & param_object)
680 : {
681 73278 : const auto & params = param_object.parameters();
682 73278 : processSingletonMooseWrappedOptions(fe_problem, params);
683 :
684 : // The parameters contained in the Action
685 73278 : const auto & petsc_options = params.get<MultiMooseEnum>("petsc_options");
686 : const auto & petsc_pair_options =
687 219834 : params.get<MooseEnumItem, std::string>("petsc_options_iname", "petsc_options_value");
688 :
689 : // A reference to the PetscOptions object that contains the settings that will be used in the
690 : // solve
691 73278 : auto & po = fe_problem.getPetscOptions();
692 :
693 : // First process the single petsc options/flags
694 73278 : addPetscFlagsToPetscOptions(petsc_options, prefix, param_object, po);
695 :
696 : // Then process the option-value pairs
697 146556 : addPetscPairsToPetscOptions(
698 73278 : petsc_pair_options, fe_problem.mesh().dimension(), prefix, param_object, po);
699 73278 : }
700 :
701 : void
702 86419 : setSolveTypeFromParams(FEProblemBase & fe_problem, const InputParameters & params)
703 : {
704 : // Note: Options set in the Preconditioner block will override those set in the Executioner block
705 330163 : if (params.isParamValid("solve_type") && !params.isParamValid("_use_eigen_value"))
706 : {
707 : // Extract the solve type
708 34888 : const std::string & solve_type = params.get<MooseEnum>("solve_type");
709 70174 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
710 35286 : fe_problem.solverParams(i)._type = Moose::stringToEnum<Moose::SolveType>(solve_type);
711 34888 : }
712 86419 : }
713 :
714 : void
715 73278 : setLineSearchFromParams(FEProblemBase & fe_problem, const InputParameters & params)
716 : {
717 : // Note: Options set in the Preconditioner block will override those set in the Executioner block
718 146556 : if (params.isParamValid("line_search"))
719 : {
720 59977 : const auto & line_search = params.get<MooseEnum>("line_search");
721 119395 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
722 59418 : if (fe_problem.solverParams(i)._line_search == Moose::LS_INVALID || line_search != "default")
723 : {
724 : Moose::LineSearchType enum_line_search =
725 59065 : Moose::stringToEnum<Moose::LineSearchType>(line_search);
726 59065 : fe_problem.solverParams(i)._line_search = enum_line_search;
727 59065 : if (enum_line_search == LS_CONTACT || enum_line_search == LS_PROJECT)
728 : {
729 0 : NonlinearImplicitSystem * nl_system = dynamic_cast<NonlinearImplicitSystem *>(
730 0 : &fe_problem.getNonlinearSystemBase(i).system());
731 0 : if (!nl_system)
732 0 : mooseError("You've requested a line search but you must be solving an EigenProblem. "
733 : "These two things are not consistent.");
734 : PetscNonlinearSolver<Real> * petsc_nonlinear_solver =
735 0 : dynamic_cast<PetscNonlinearSolver<Real> *>(nl_system->nonlinear_solver.get());
736 0 : if (!petsc_nonlinear_solver)
737 0 : mooseError("Currently the MOOSE line searches all use Petsc, so you "
738 : "must use Petsc as your non-linear solver.");
739 : petsc_nonlinear_solver->linesearch_object =
740 0 : std::make_unique<ComputeLineSearchObjectWrapper>(fe_problem);
741 : }
742 : }
743 : }
744 73278 : }
745 :
746 : void
747 73278 : setMFFDTypeFromParams(FEProblemBase & fe_problem, const InputParameters & params)
748 : {
749 146556 : if (params.isParamValid("mffd_type"))
750 : {
751 72770 : const auto & mffd_type = params.get<MooseEnum>("mffd_type");
752 145227 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
753 72457 : fe_problem.solverParams(i)._mffd_type = Moose::stringToEnum<Moose::MffdType>(mffd_type);
754 : }
755 73278 : }
756 :
757 : template <typename T>
758 : void
759 54683 : checkUserProvidedPetscOption(const T & option, const ParallelParamObject & param_object)
760 : {
761 54683 : const auto & string_option = static_cast<const std::string &>(option);
762 54683 : if (string_option[0] != '-')
763 0 : param_object.mooseError("PETSc option '", string_option, "' does not begin with '-'");
764 54683 : }
765 :
766 : void
767 73278 : addPetscFlagsToPetscOptions(const MultiMooseEnum & petsc_flags,
768 : std::string prefix,
769 : const ParallelParamObject & param_object,
770 : PetscOptions & po)
771 : {
772 73278 : prefix.insert(prefix.begin(), '-');
773 : checkPrefix(prefix);
774 :
775 : // Update the PETSc single flags
776 75623 : for (const auto & option : petsc_flags)
777 : {
778 2345 : checkUserProvidedPetscOption(option, param_object);
779 :
780 2345 : const std::string & string_option = option.name();
781 :
782 : /**
783 : * "-log_summary" cannot be used in the input file. This option needs to be set when PETSc is
784 : * initialized
785 : * which happens before the parser is even created. We'll throw an error if somebody attempts
786 : * to add this option later.
787 : */
788 2345 : if (option == "-log_summary" || option == "-log_view")
789 0 : mooseError("The PETSc option \"-log_summary\" or \"-log_view\" can only be used on the "
790 : "command line. Please "
791 : "remove it from the input file");
792 :
793 : // Update the stored items, but do not create duplicates
794 2345 : const std::string prefixed_option = prefix + string_option.substr(1);
795 2345 : if (!po.flags.isValueSet(prefixed_option))
796 : {
797 2345 : po.flags.setAdditionalValue(prefixed_option);
798 2345 : po.user_set_options.setAdditionalValue(prefixed_option);
799 : }
800 2345 : }
801 73278 : }
802 :
803 : void
804 59977 : setConvergedReasonFlags(FEProblemBase & fe_problem, std::string prefix)
805 : {
806 59977 : prefix.insert(prefix.begin(), '-');
807 : checkPrefix(prefix);
808 59977 : libmesh_ignore(fe_problem); // avoid unused warnings for old PETSc
809 :
810 : #if !PETSC_VERSION_LESS_THAN(3, 14, 0)
811 : // the boolean in these pairs denote whether the user has specified any of the reason flags in the
812 : // input file
813 119954 : std::array<std::string, 2> reason_flags = {{"snes_converged_reason", "ksp_converged_reason"}};
814 :
815 59977 : auto & po = fe_problem.getPetscOptions();
816 :
817 179931 : for (const auto & reason_flag : reason_flags)
818 : {
819 119954 : const auto full_flag = prefix + reason_flag;
820 237509 : if (!po.flags.isValueSet(full_flag) && !po.dont_add_these_options.contains(full_flag) &&
821 117555 : (std::find_if(po.pairs.begin(),
822 : po.pairs.end(),
823 160139 : [&full_flag](auto & pair)
824 515203 : { return pair.first == (full_flag); }) == po.pairs.end()))
825 117555 : po.pairs.emplace_back(full_flag, "::failed");
826 119954 : }
827 : #endif
828 59977 : }
829 :
830 : void
831 73499 : addPetscPairsToPetscOptions(
832 : const std::vector<std::pair<MooseEnumItem, std::string>> & petsc_pair_options,
833 : const unsigned int mesh_dimension,
834 : std::string prefix,
835 : const ParallelParamObject & param_object,
836 : PetscOptions & po)
837 : {
838 73499 : prefix.insert(prefix.begin(), '-');
839 : checkPrefix(prefix);
840 :
841 : // Setup the name value pairs
842 73499 : bool boomeramg_found = false;
843 73499 : bool strong_threshold_found = false;
844 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
845 73499 : bool superlu_dist_found = false;
846 73499 : bool fact_pattern_found = false;
847 73499 : bool tiny_pivot_found = false;
848 : #endif
849 73499 : std::string pc_description = "";
850 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
851 : // If users use HMG, we would like to set
852 73499 : bool hmg_found = false;
853 73499 : bool matptap_found = false;
854 73499 : bool hmg_strong_threshold_found = false;
855 : #endif
856 73499 : std::vector<std::pair<std::string, std::string>> new_options;
857 :
858 125837 : for (const auto & [option_name, option_value] : petsc_pair_options)
859 : {
860 52338 : checkUserProvidedPetscOption(option_name, param_object);
861 :
862 52338 : new_options.clear();
863 : const std::string prefixed_option_name =
864 52338 : prefix + static_cast<const std::string &>(option_name).substr(1);
865 :
866 : // Do not add duplicate settings
867 52338 : if (auto it =
868 52338 : MooseUtils::findPair(po.pairs, po.pairs.begin(), prefixed_option_name, MooseUtils::Any);
869 52338 : it == po.pairs.end())
870 : {
871 : #if !PETSC_VERSION_LESS_THAN(3, 9, 0)
872 52230 : if (option_name == "-pc_factor_mat_solver_package")
873 142 : new_options.emplace_back(prefix + "pc_factor_mat_solver_type", option_value);
874 : #else
875 : if (option_name == "-pc_factor_mat_solver_type")
876 : new_options.push_back(prefix + "pc_factor_mat_solver_package", option_value);
877 : #endif
878 :
879 : // Look for a pc description
880 79083 : if (option_name == "-pc_type" || option_name == "-sub_pc_type" ||
881 26853 : option_name == "-pc_hypre_type")
882 44804 : pc_description += option_value + ' ';
883 :
884 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
885 52230 : if (option_name == "-pc_type" && option_value == "hmg")
886 36 : hmg_found = true;
887 :
888 : // MPIAIJ for PETSc 3.12.0: -matptap_via
889 : // MAIJ for PETSc 3.12.0: -matmaijptap_via
890 : // MPIAIJ for PETSc 3.13 to 3.16: -matptap_via, -matproduct_ptap_via
891 : // MAIJ for PETSc 3.13 to 3.16: -matproduct_ptap_via
892 : // MPIAIJ for PETSc 3.17 and higher: -matptap_via, -mat_product_algorithm
893 : // MAIJ for PETSc 3.17 and higher: -mat_product_algorithm
894 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
895 52338 : if (hmg_found && (option_name == "-matptap_via" || option_name == "-matmaijptap_via" ||
896 108 : option_name == "-matproduct_ptap_via"))
897 0 : new_options.emplace_back(prefix + "mat_product_algorithm", option_value);
898 : #elif !PETSC_VERSION_LESS_THAN(3, 13, 0)
899 : if (hmg_found && (option_name == "-matptap_via" || option_name == "-matmaijptap_via"))
900 : new_options.emplace_back(prefix + "matproduct_ptap_via", option_value);
901 : #else
902 : if (hmg_found && (option_name == "-matproduct_ptap_via"))
903 : {
904 : new_options.emplace_back(prefix + "matptap_via", option_value);
905 : new_options.emplace_back(prefix + "matmaijptap_via", option_value);
906 : }
907 : #endif
908 :
909 104460 : if (option_name == "-matptap_via" || option_name == "-matmaijptap_via" ||
910 104460 : option_name == "-matproduct_ptap_via" || option_name == "-mat_product_algorithm")
911 0 : matptap_found = true;
912 :
913 : // For 3D problems, we need to set this 0.7
914 52230 : if (option_name == "-hmg_inner_pc_hypre_boomeramg_strong_threshold")
915 0 : hmg_strong_threshold_found = true;
916 : #endif
917 : // This special case is common enough that we'd like to handle it for the user.
918 52230 : if (option_name == "-pc_hypre_type" && option_value == "boomeramg")
919 19331 : boomeramg_found = true;
920 52230 : if (option_name == "-pc_hypre_boomeramg_strong_threshold")
921 0 : strong_threshold_found = true;
922 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
923 52230 : if ((option_name == "-pc_factor_mat_solver_package" ||
924 53382 : option_name == "-pc_factor_mat_solver_type") &&
925 1152 : option_value == "superlu_dist")
926 285 : superlu_dist_found = true;
927 52230 : if (option_name == "-mat_superlu_dist_fact")
928 0 : fact_pattern_found = true;
929 52230 : if (option_name == "-mat_superlu_dist_replacetinypivot")
930 0 : tiny_pivot_found = true;
931 : #endif
932 :
933 52230 : if (!new_options.empty())
934 : {
935 142 : std::copy(new_options.begin(), new_options.end(), std::back_inserter(po.pairs));
936 284 : for (const auto & option : new_options)
937 142 : po.user_set_options.setAdditionalValue(option.first);
938 : }
939 : else
940 : {
941 52088 : po.pairs.push_back(std::make_pair(prefixed_option_name, option_value));
942 52088 : po.user_set_options.setAdditionalValue(prefixed_option_name);
943 : }
944 : }
945 : else
946 : {
947 : do
948 : {
949 108 : it->second = option_value;
950 108 : it = MooseUtils::findPair(po.pairs, std::next(it), prefixed_option_name, MooseUtils::Any);
951 108 : } while (it != po.pairs.end());
952 : }
953 52338 : }
954 :
955 : // When running a 3D mesh with boomeramg, it is almost always best to supply a strong threshold
956 : // value. We will provide that for the user here if they haven't supplied it themselves.
957 73499 : if (boomeramg_found && !strong_threshold_found && mesh_dimension == 3)
958 : {
959 961 : po.pairs.emplace_back(prefix + "pc_hypre_boomeramg_strong_threshold", "0.7");
960 961 : pc_description += "strong_threshold: 0.7 (auto)";
961 : }
962 :
963 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
964 73499 : if (hmg_found && !hmg_strong_threshold_found && mesh_dimension == 3)
965 : {
966 24 : po.pairs.emplace_back(prefix + "hmg_inner_pc_hypre_boomeramg_strong_threshold", "0.7");
967 24 : pc_description += "strong_threshold: 0.7 (auto)";
968 : }
969 :
970 : // Default PETSc PtAP takes too much memory, and it is not quite useful
971 : // Let us switch to use new algorithm
972 73499 : if (hmg_found && !matptap_found)
973 : {
974 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
975 36 : po.pairs.emplace_back(prefix + "mat_product_algorithm", "allatonce");
976 : #elif !PETSC_VERSION_LESS_THAN(3, 13, 0)
977 : po.pairs.emplace_back(prefix + "matproduct_ptap_via", "allatonce");
978 : #else
979 : po.pairs.emplace_back(prefix + "matptap_via", "allatonce");
980 : po.pairs.emplace_back(prefix + "matmaijptap_via", "allatonce");
981 : #endif
982 : }
983 : #endif
984 :
985 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
986 : // In PETSc-3.7.{0--4}, there is a bug when using superlu_dist, and we have to use
987 : // SamePattern_SameRowPerm, otherwise we use whatever we have in PETSc
988 73499 : if (superlu_dist_found && !fact_pattern_found)
989 : {
990 285 : po.pairs.emplace_back(prefix + "mat_superlu_dist_fact",
991 : #if PETSC_VERSION_LESS_THAN(3, 7, 5)
992 : "SamePattern_SameRowPerm");
993 : pc_description += "mat_superlu_dist_fact: SamePattern_SameRowPerm ";
994 : #else
995 : "SamePattern");
996 285 : pc_description += "mat_superlu_dist_fact: SamePattern ";
997 : #endif
998 : }
999 :
1000 : // restore this superlu option
1001 73499 : if (superlu_dist_found && !tiny_pivot_found)
1002 : {
1003 285 : po.pairs.emplace_back(prefix + "mat_superlu_dist_replacetinypivot", "1");
1004 285 : pc_description += " mat_superlu_dist_replacetinypivot: true ";
1005 : }
1006 : #endif
1007 : // Set Preconditioner description
1008 73499 : if (!pc_description.empty() && prefix.size() > 1)
1009 729 : po.pc_description += "[" + prefix.substr(1, prefix.size() - 2) + "]: ";
1010 73499 : po.pc_description += pc_description;
1011 73499 : }
1012 :
1013 : std::set<std::string>
1014 156854 : getPetscValidLineSearches()
1015 : {
1016 156854 : return {"default", "shell", "none", "basic", "l2", "bt", "cp"};
1017 : }
1018 :
1019 : InputParameters
1020 204584 : getPetscValidParams()
1021 : {
1022 204584 : InputParameters params = emptyInputParameters();
1023 :
1024 818336 : MooseEnum solve_type("PJFNK JFNK NEWTON FD LINEAR");
1025 818336 : params.addParam<MooseEnum>("solve_type",
1026 : solve_type,
1027 : "PJFNK: Preconditioned Jacobian-Free Newton Krylov "
1028 : "JFNK: Jacobian-Free Newton Krylov "
1029 : "NEWTON: Full Newton Solve "
1030 : "FD: Use finite differences to compute Jacobian "
1031 : "LINEAR: Solving a linear problem");
1032 :
1033 818336 : MooseEnum mffd_type("wp ds", "wp");
1034 818336 : params.addParam<MooseEnum>("mffd_type",
1035 : mffd_type,
1036 : "Specifies the finite differencing type for "
1037 : "Jacobian-free solve types. Note that the "
1038 : "default is wp (for Walker and Pernice).");
1039 :
1040 613752 : params.addParam<MultiMooseEnum>(
1041 409168 : "petsc_options", getCommonPetscFlags(), "Singleton PETSc options");
1042 613752 : params.addParam<MultiMooseEnum>(
1043 409168 : "petsc_options_iname", getCommonPetscKeys(), "Names of PETSc name/value pairs");
1044 818336 : params.addParam<std::vector<std::string>>(
1045 : "petsc_options_value",
1046 : "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\"");
1047 613752 : params.addParamNamesToGroup("solve_type petsc_options petsc_options_iname petsc_options_value "
1048 : "mffd_type",
1049 : "PETSc");
1050 :
1051 409168 : return params;
1052 204584 : }
1053 :
1054 : MultiMooseEnum
1055 210670 : getCommonSNESFlags()
1056 : {
1057 : return MultiMooseEnum(
1058 : "-ksp_monitor_snes_lg -snes_ksp_ew -snes_converged_reason "
1059 : "-snes_ksp -snes_linesearch_monitor -snes_mf -snes_mf_operator -snes_monitor "
1060 : "-snes_test_display -snes_view -snes_monitor_cancel",
1061 : "",
1062 421340 : true);
1063 : }
1064 :
1065 : MultiMooseEnum
1066 209638 : getCommonKSPFlags()
1067 : {
1068 : return MultiMooseEnum(
1069 419276 : "-ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_monitor", "", true);
1070 : }
1071 :
1072 : MultiMooseEnum
1073 208974 : getCommonPetscFlags()
1074 : {
1075 208974 : auto options = MultiMooseEnum("-dm_moose_print_embedding -dm_view", "", true);
1076 208974 : options.addValidName(getCommonKSPFlags());
1077 208974 : options.addValidName(getCommonSNESFlags());
1078 208974 : return options;
1079 0 : }
1080 :
1081 : MultiMooseEnum
1082 210670 : getCommonSNESKeys()
1083 : {
1084 : return MultiMooseEnum("-snes_atol -snes_linesearch_type -snes_ls -snes_max_it -snes_rtol "
1085 : "-snes_divergence_tolerance -snes_type",
1086 : "",
1087 421340 : true);
1088 : }
1089 :
1090 : MultiMooseEnum
1091 209638 : getCommonKSPKeys()
1092 : {
1093 : return MultiMooseEnum("-ksp_atol -ksp_gmres_restart -ksp_max_it -ksp_pc_side -ksp_rtol "
1094 : "-ksp_type -sub_ksp_type",
1095 : "",
1096 419276 : true);
1097 : }
1098 : MultiMooseEnum
1099 208974 : getCommonPetscKeys()
1100 : {
1101 : auto options = MultiMooseEnum("-mat_fd_coloring_err -mat_fd_type -mat_mffd_type "
1102 : "-pc_asm_overlap -pc_factor_levels "
1103 : "-pc_factor_mat_ordering_type -pc_hypre_boomeramg_grid_sweeps_all "
1104 : "-pc_hypre_boomeramg_max_iter "
1105 : "-pc_hypre_boomeramg_strong_threshold -pc_hypre_type -pc_type "
1106 : "-sub_pc_type",
1107 : "",
1108 208974 : true);
1109 208974 : options.addValidName(getCommonKSPKeys());
1110 208974 : options.addValidName(getCommonSNESKeys());
1111 208974 : return options;
1112 0 : }
1113 :
1114 : bool
1115 560 : isSNESVI(FEProblemBase & fe_problem)
1116 : {
1117 560 : const PetscOptions & petsc = fe_problem.getPetscOptions();
1118 :
1119 : int argc;
1120 : char ** args;
1121 560 : LibmeshPetscCallA(fe_problem.comm().get(), PetscGetArgs(&argc, &args));
1122 :
1123 560 : std::vector<std::string> cml_arg;
1124 6822 : for (int i = 0; i < argc; i++)
1125 12524 : cml_arg.push_back(args[i]);
1126 :
1127 560 : if (MooseUtils::findPair(petsc.pairs, petsc.pairs.begin(), MooseUtils::Any, "vinewtonssls") ==
1128 1120 : petsc.pairs.end() &&
1129 560 : MooseUtils::findPair(petsc.pairs, petsc.pairs.begin(), MooseUtils::Any, "vinewtonrsls") ==
1130 600 : petsc.pairs.end() &&
1131 1160 : std::find(cml_arg.begin(), cml_arg.end(), "vinewtonssls") == cml_arg.end() &&
1132 600 : std::find(cml_arg.begin(), cml_arg.end(), "vinewtonrsls") == cml_arg.end())
1133 40 : return false;
1134 :
1135 520 : return true;
1136 560 : }
1137 :
1138 : void
1139 218858 : setSinglePetscOption(const std::string & name,
1140 : const std::string & value /*=""*/,
1141 : FEProblemBase * const problem /*=nullptr*/)
1142 : {
1143 218858 : static const TIMPI::Communicator comm_world(PETSC_COMM_WORLD);
1144 218858 : const TIMPI::Communicator & comm = problem ? problem->comm() : comm_world;
1145 218858 : LibmeshPetscCallA(comm.get(),
1146 : PetscOptionsSetValue(LIBMESH_PETSC_NULLPTR,
1147 : name.c_str(),
1148 : value == "" ? LIBMESH_PETSC_NULLPTR : value.c_str()));
1149 : // Create a single option data base so that we can use PETSC's internal option checking which
1150 : // is case insensitive. This is better than re-implementing case-insensitive checks here in this
1151 : // TU
1152 : ::PetscOptions single_option;
1153 218858 : LibmeshPetscCallA(comm.get(), PetscOptionsCreate(&single_option));
1154 218858 : LibmeshPetscCallA(comm.get(),
1155 : PetscOptionsSetValue(single_option,
1156 : name.c_str(),
1157 : value == "" ? LIBMESH_PETSC_NULLPTR : value.c_str()));
1158 23 : auto check_problem = [problem, &name]()
1159 : {
1160 23 : if (!problem)
1161 0 : mooseError(
1162 : "Setting the option '",
1163 : name,
1164 : "' requires passing a 'problem' parameter. Contact a developer of your application "
1165 : "to have them update their code. If in doubt, reach out to the MOOSE team on Github "
1166 : "discussions");
1167 218858 : };
1168 :
1169 : // Select vector type from user-passed PETSc options
1170 875432 : if (petscOptionsHasName(single_option, "-vec_type"))
1171 : {
1172 23 : check_problem();
1173 23 : applyVectorTypeOptions(*problem);
1174 : }
1175 : // First do a cheap suffix check so unrelated PETSc options do not pay for looping over every
1176 : // solver system. Once we know the name looks like a matrix-type option, rely on PETSc's
1177 : // option lookup for the actual case-insensitive and prefix-aware matching.
1178 218835 : else if (problem && mightBeMatTypeOption(name))
1179 : {
1180 23 : errorOnUnprefixedMatTypeOption(single_option, *problem);
1181 :
1182 23 : for (const auto i : index_range(problem->_solver_systems))
1183 : {
1184 23 : const auto & solver_sys_name = problem->_solver_sys_names[i];
1185 69 : if (!petscOptionsHasName(single_option, "-mat_type", solver_sys_name + "_"))
1186 0 : continue;
1187 :
1188 23 : if (problem->solverParams(i)._type == Moose::ST_JFNK)
1189 0 : mooseError("Setting option '", name, "' is incompatible with a JFNK 'solve_type'");
1190 :
1191 23 : applyMatrixTypeOptions(*problem, i);
1192 23 : break;
1193 : }
1194 : }
1195 :
1196 218858 : LibmeshPetscCallA(comm.get(), PetscOptionsDestroy(&single_option));
1197 218858 : }
1198 :
1199 : void
1200 68455 : setSinglePetscOptionIfAppropriate(const MultiMooseEnum & dont_add_these_options,
1201 : const std::string & name,
1202 : const std::string & value /*=""*/,
1203 : FEProblemBase * const problem /*=nullptr*/)
1204 : {
1205 68455 : if (!dont_add_these_options.contains(name))
1206 64106 : setSinglePetscOption(name, value, problem);
1207 68455 : }
1208 :
1209 : void
1210 0 : colorAdjacencyMatrix(PetscScalar * adjacency_matrix,
1211 : unsigned int size,
1212 : unsigned int colors,
1213 : std::vector<unsigned int> & vertex_colors,
1214 : const char * coloring_algorithm)
1215 : {
1216 : // Mat A will be a dense matrix from the incoming data structure
1217 : Mat A;
1218 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatCreate(PETSC_COMM_SELF, &A));
1219 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSetSizes(A, size, size, size, size));
1220 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSetType(A, MATSEQDENSE));
1221 : // PETSc requires a non-const data array to populate the matrix
1222 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSeqDenseSetPreallocation(A, adjacency_matrix));
1223 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1224 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1225 :
1226 : // Convert A to a sparse matrix
1227 : #if PETSC_VERSION_LESS_THAN(3, 7, 0)
1228 : LibmeshPetscCallA(PETSC_COMM_SELF, MatConvert(A, MATAIJ, MAT_REUSE_MATRIX, &A));
1229 : #else
1230 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A));
1231 : #endif
1232 :
1233 : ISColoring iscoloring;
1234 : MatColoring mc;
1235 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringCreate(A, &mc));
1236 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetType(mc, coloring_algorithm));
1237 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetMaxColors(mc, static_cast<PetscInt>(colors)));
1238 :
1239 : // Petsc normally colors by distance two (neighbors of neighbors), we just want one
1240 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetDistance(mc, 1));
1241 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetFromOptions(mc));
1242 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringApply(mc, &iscoloring));
1243 :
1244 : PetscInt nn;
1245 : IS * is;
1246 : #if PETSC_RELEASE_LESS_THAN(3, 12, 0)
1247 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringGetIS(iscoloring, &nn, &is));
1248 : #else
1249 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &nn, &is));
1250 : #endif
1251 :
1252 0 : if (nn > static_cast<PetscInt>(colors))
1253 0 : throw std::runtime_error("Not able to color with designated number of colors");
1254 :
1255 0 : for (int i = 0; i < nn; i++)
1256 : {
1257 : PetscInt isize;
1258 : const PetscInt * indices;
1259 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISGetLocalSize(is[i], &isize));
1260 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISGetIndices(is[i], &indices));
1261 0 : for (int j = 0; j < isize; j++)
1262 : {
1263 : mooseAssert(indices[j] < static_cast<PetscInt>(vertex_colors.size()), "Index out of bounds");
1264 0 : vertex_colors[indices[j]] = i;
1265 : }
1266 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISRestoreIndices(is[i], &indices));
1267 : }
1268 :
1269 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatDestroy(&A));
1270 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringDestroy(&mc));
1271 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringDestroy(&iscoloring));
1272 0 : }
1273 :
1274 : void
1275 40756 : dontAddPetscFlag(const std::string & flag, PetscOptions & petsc_options)
1276 : {
1277 40756 : if (!petsc_options.dont_add_these_options.contains(flag))
1278 31460 : petsc_options.dont_add_these_options.setAdditionalValue(flag);
1279 40756 : }
1280 :
1281 : void
1282 625 : dontAddNonlinearConvergedReason(FEProblemBase & fe_problem)
1283 : {
1284 1250 : dontAddPetscFlag("-snes_converged_reason", fe_problem.getPetscOptions());
1285 625 : }
1286 :
1287 : void
1288 634 : dontAddLinearConvergedReason(FEProblemBase & fe_problem)
1289 : {
1290 1268 : dontAddPetscFlag("-ksp_converged_reason", fe_problem.getPetscOptions());
1291 634 : }
1292 :
1293 : void
1294 664 : dontAddCommonKSPOptions(FEProblemBase & fe_problem)
1295 : {
1296 664 : auto & petsc_options = fe_problem.getPetscOptions();
1297 2656 : for (const auto & flag : getCommonKSPFlags().getNames())
1298 2656 : dontAddPetscFlag(flag, petsc_options);
1299 5312 : for (const auto & key : getCommonKSPKeys().getNames())
1300 5312 : dontAddPetscFlag(key, petsc_options);
1301 664 : }
1302 :
1303 : void
1304 664 : dontAddCommonSNESOptions(FEProblemBase & fe_problem)
1305 : {
1306 664 : dontAddCommonSNESOptions(fe_problem, "");
1307 664 : }
1308 :
1309 : void
1310 1696 : dontAddCommonSNESOptions(FEProblemBase & fe_problem, const std::string & prefix)
1311 : {
1312 1696 : auto & petsc_options = fe_problem.getPetscOptions();
1313 20352 : for (const auto & flag : getCommonSNESFlags().getNames())
1314 20352 : dontAddPetscFlag("-" + prefix + flag.substr(1), petsc_options);
1315 13568 : for (const auto & key : getCommonSNESKeys().getNames())
1316 13568 : dontAddPetscFlag("-" + prefix + key.substr(1), petsc_options);
1317 1696 : }
1318 :
1319 : std::unique_ptr<PetscMatrix<Number>>
1320 40 : createMatrixFromFile(const libMesh::Parallel::Communicator & comm,
1321 : Mat & mat,
1322 : const std::string & binary_mat_file,
1323 : const unsigned int mat_number_to_load)
1324 : {
1325 40 : LibmeshPetscCallA(comm.get(), MatCreate(comm.get(), &mat));
1326 : PetscViewer matviewer;
1327 40 : LibmeshPetscCallA(
1328 : comm.get(),
1329 : PetscViewerBinaryOpen(comm.get(), binary_mat_file.c_str(), FILE_MODE_READ, &matviewer));
1330 80 : for (unsigned int i = 0; i < mat_number_to_load; ++i)
1331 40 : LibmeshPetscCallA(comm.get(), MatLoad(mat, matviewer));
1332 40 : LibmeshPetscCallA(comm.get(), PetscViewerDestroy(&matviewer));
1333 :
1334 80 : return std::make_unique<PetscMatrix<Number>>(mat, comm);
1335 : }
1336 :
1337 : } // Namespace PetscSupport
1338 : } // Namespace MOOSE
|
