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 "DisplacedProblem.h"
19 : #include "PenetrationLocator.h"
20 : #include "NearestNodeLocator.h"
21 : #include "MooseTypes.h"
22 : #include "MooseUtils.h"
23 : #include "CommandLine.h"
24 : #include "Console.h"
25 : #include "MultiMooseEnum.h"
26 : #include "Conversion.h"
27 : #include "Executioner.h"
28 : #include "MooseMesh.h"
29 : #include "ComputeLineSearchObjectWrapper.h"
30 : #include "Convergence.h"
31 : #include "ParallelParamObject.h"
32 :
33 : #include "libmesh/equation_systems.h"
34 : #include "libmesh/linear_implicit_system.h"
35 : #include "libmesh/nonlinear_implicit_system.h"
36 : #include "libmesh/petsc_linear_solver.h"
37 : #include "libmesh/petsc_matrix.h"
38 : #include "libmesh/petsc_nonlinear_solver.h"
39 : #include "libmesh/petsc_preconditioner.h"
40 : #include "libmesh/petsc_vector.h"
41 : #include "libmesh/sparse_matrix.h"
42 : #include "libmesh/petsc_solver_exception.h"
43 : #include "libmesh/simple_range.h"
44 :
45 : // PETSc includes
46 : #include <petsc.h>
47 : #include <petscsnes.h>
48 : #include <petscksp.h>
49 :
50 : // For graph coloring
51 : #include <petscmat.h>
52 : #include <petscis.h>
53 : #include <petscdm.h>
54 :
55 : // PetscDMMoose include
56 : #include "PetscDMMoose.h"
57 :
58 : // Standard includes
59 : #include <ostream>
60 : #include <fstream>
61 : #include <string>
62 :
63 : using namespace libMesh;
64 :
65 : void
66 0 : MooseVecView(NumericVector<Number> & vector)
67 : {
68 0 : PetscVector<Number> & petsc_vec = static_cast<PetscVector<Number> &>(vector);
69 0 : LibmeshPetscCallA(vector.comm().get(), VecView(petsc_vec.vec(), 0));
70 0 : }
71 :
72 : void
73 0 : MooseMatView(SparseMatrix<Number> & mat)
74 : {
75 0 : PetscMatrixBase<Number> & petsc_mat = static_cast<PetscMatrix<Number> &>(mat);
76 0 : LibmeshPetscCallA(mat.comm().get(), MatView(petsc_mat.mat(), 0));
77 0 : }
78 :
79 : void
80 0 : MooseVecView(const NumericVector<Number> & vector)
81 : {
82 0 : PetscVector<Number> & petsc_vec =
83 : static_cast<PetscVector<Number> &>(const_cast<NumericVector<Number> &>(vector));
84 0 : LibmeshPetscCallA(vector.comm().get(), VecView(petsc_vec.vec(), 0));
85 0 : }
86 :
87 : void
88 0 : MooseMatView(const SparseMatrix<Number> & mat)
89 : {
90 0 : PetscMatrixBase<Number> & petsc_mat =
91 : static_cast<PetscMatrix<Number> &>(const_cast<SparseMatrix<Number> &>(mat));
92 0 : LibmeshPetscCallA(mat.comm().get(), MatView(petsc_mat.mat(), 0));
93 0 : }
94 :
95 : namespace Moose
96 : {
97 : namespace PetscSupport
98 : {
99 :
100 : std::string
101 2695 : stringify(const LineSearchType & t)
102 : {
103 2695 : switch (t)
104 : {
105 2607 : case LS_BASIC:
106 2607 : return "basic";
107 0 : case LS_DEFAULT:
108 0 : return "default";
109 0 : case LS_NONE:
110 0 : return "none";
111 0 : case LS_SHELL:
112 0 : return "shell";
113 0 : case LS_L2:
114 0 : return "l2";
115 22 : case LS_BT:
116 22 : return "bt";
117 66 : case LS_CP:
118 66 : return "cp";
119 0 : case LS_CONTACT:
120 0 : return "contact";
121 0 : case LS_PROJECT:
122 0 : return "project";
123 0 : case LS_INVALID:
124 0 : mooseError("Invalid LineSearchType");
125 : }
126 0 : return "";
127 : }
128 :
129 : std::string
130 26896 : stringify(const MffdType & t)
131 : {
132 26896 : switch (t)
133 : {
134 26886 : case MFFD_WP:
135 26886 : return "wp";
136 10 : case MFFD_DS:
137 10 : return "ds";
138 0 : case MFFD_INVALID:
139 0 : mooseError("Invalid MffdType");
140 : }
141 0 : return "";
142 : }
143 :
144 : void
145 39627 : setSolverOptions(const SolverParams & solver_params, const MultiMooseEnum & dont_add_these_options)
146 : {
147 : // set PETSc options implied by a solve type
148 39627 : switch (solver_params._type)
149 : {
150 26779 : case Moose::ST_PJFNK:
151 80337 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
152 53558 : solver_params._prefix + "snes_mf_operator");
153 26779 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
154 53558 : solver_params._prefix + "mat_mffd_type",
155 53558 : stringify(solver_params._mffd_type));
156 26779 : break;
157 :
158 117 : case Moose::ST_JFNK:
159 117 : setSinglePetscOptionIfAppropriate(dont_add_these_options, solver_params._prefix + "snes_mf");
160 117 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
161 234 : solver_params._prefix + "mat_mffd_type",
162 234 : stringify(solver_params._mffd_type));
163 117 : break;
164 :
165 10687 : case Moose::ST_NEWTON:
166 10687 : break;
167 :
168 0 : case Moose::ST_FD:
169 0 : setSinglePetscOptionIfAppropriate(dont_add_these_options, solver_params._prefix + "snes_fd");
170 0 : break;
171 :
172 2044 : case Moose::ST_LINEAR:
173 6132 : setSinglePetscOptionIfAppropriate(
174 4088 : dont_add_these_options, solver_params._prefix + "snes_type", "ksponly");
175 6132 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
176 4088 : solver_params._prefix + "snes_monitor_cancel");
177 2044 : break;
178 : }
179 :
180 39627 : Moose::LineSearchType ls_type = solver_params._line_search;
181 39627 : if (ls_type == Moose::LS_NONE)
182 2453 : ls_type = Moose::LS_BASIC;
183 :
184 39627 : if (ls_type != Moose::LS_DEFAULT && ls_type != Moose::LS_CONTACT && ls_type != Moose::LS_PROJECT)
185 2695 : setSinglePetscOptionIfAppropriate(
186 5390 : dont_add_these_options, solver_params._prefix + "snes_linesearch_type", stringify(ls_type));
187 39627 : }
188 :
189 : void
190 95 : petscSetupDM(NonlinearSystemBase & nl, const std::string & dm_name)
191 : {
192 : PetscBool ismoose;
193 95 : DM dm = LIBMESH_PETSC_NULLPTR;
194 :
195 : // Initialize the part of the DM package that's packaged with Moose; in the PETSc source tree this
196 : // call would be in DMInitializePackage()
197 95 : LibmeshPetscCallA(nl.comm().get(), DMMooseRegisterAll());
198 : // Create and set up the DM that will consume the split options and deal with block matrices.
199 : PetscNonlinearSolver<Number> * petsc_solver =
200 95 : dynamic_cast<PetscNonlinearSolver<Number> *>(nl.nonlinearSolver());
201 95 : const char * snes_prefix = nullptr;
202 95 : std::string snes_prefix_str;
203 95 : if (nl.system().prefix_with_name())
204 : {
205 0 : snes_prefix_str = nl.system().prefix();
206 0 : snes_prefix = snes_prefix_str.c_str();
207 : }
208 95 : SNES snes = petsc_solver->snes(snes_prefix);
209 : // if there exists a DMMoose object, not to recreate a new one
210 95 : LibmeshPetscCallA(nl.comm().get(), SNESGetDM(snes, &dm));
211 95 : if (dm)
212 : {
213 95 : LibmeshPetscCallA(nl.comm().get(), PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose));
214 95 : if (ismoose)
215 0 : return;
216 : }
217 95 : LibmeshPetscCallA(nl.comm().get(), DMCreateMoose(nl.comm().get(), nl, dm_name, &dm));
218 95 : LibmeshPetscCallA(nl.comm().get(), DMSetFromOptions(dm));
219 95 : LibmeshPetscCallA(nl.comm().get(), DMSetUp(dm));
220 95 : LibmeshPetscCallA(nl.comm().get(), SNESSetDM(snes, dm));
221 95 : LibmeshPetscCallA(nl.comm().get(), DMDestroy(&dm));
222 : // We temporarily comment out this updating function because
223 : // we lack an approach to check if the problem
224 : // structure has been changed from the last iteration.
225 : // The indices will be rebuilt for every timestep.
226 : // TODO: figure out a way to check the structure changes of the
227 : // matrix
228 : // ierr = SNESSetUpdate(snes,SNESUpdateDMMoose);
229 : // CHKERRABORT(nl.comm().get(),ierr);
230 95 : }
231 :
232 : void
233 39975 : addPetscOptionsFromCommandline()
234 : {
235 : // commandline options always win
236 : // the options from a user commandline will overwrite the existing ones if any conflicts
237 : { // Get any options specified on the command-line
238 : int argc;
239 : char ** args;
240 :
241 39975 : LibmeshPetscCallA(PETSC_COMM_WORLD, PetscGetArgs(&argc, &args));
242 : #if PETSC_VERSION_LESS_THAN(3, 7, 0)
243 : LibmeshPetscCallA(PETSC_COMM_WORLD, PetscOptionsInsert(&argc, &args, NULL));
244 : #else
245 39975 : LibmeshPetscCallA(PETSC_COMM_WORLD,
246 : PetscOptionsInsert(LIBMESH_PETSC_NULLPTR, &argc, &args, NULL));
247 : #endif
248 : }
249 39975 : }
250 :
251 : void
252 39443 : petscSetOptionsHelper(const PetscOptions & po, FEProblemBase * const problem)
253 : {
254 : // Add any additional options specified in the input file
255 41915 : for (const auto & flag : po.flags)
256 : // Need to use name method here to pass a str instead of an EnumItem because
257 : // we don't care if the id attributes match
258 2472 : if (!po.dont_add_these_options.contains(flag.name()) ||
259 0 : po.user_set_options.contains(flag.name()))
260 2472 : setSinglePetscOption(flag.rawName().c_str());
261 :
262 : // Add option pairs
263 166662 : for (auto & option : po.pairs)
264 128666 : if (!po.dont_add_these_options.contains(option.first) ||
265 1447 : po.user_set_options.contains(option.first))
266 125782 : setSinglePetscOption(option.first, option.second, problem);
267 :
268 39443 : addPetscOptionsFromCommandline();
269 39443 : }
270 :
271 : void
272 1610 : petscSetOptions(const PetscOptions & po,
273 : const SolverParams & solver_params,
274 : FEProblemBase * const problem)
275 : {
276 1610 : PetscCallAbort(PETSC_COMM_WORLD, PetscOptionsClear(LIBMESH_PETSC_NULLPTR));
277 1610 : setSolverOptions(solver_params, po.dont_add_these_options);
278 1610 : petscSetOptionsHelper(po, problem);
279 1610 : }
280 :
281 : void
282 37833 : petscSetOptions(const PetscOptions & po,
283 : const std::vector<SolverParams> & solver_params_vec,
284 : FEProblemBase * const problem)
285 : {
286 37833 : PetscCallAbort(PETSC_COMM_WORLD, PetscOptionsClear(LIBMESH_PETSC_NULLPTR));
287 75850 : for (const auto & solver_params : solver_params_vec)
288 38017 : setSolverOptions(solver_params, po.dont_add_these_options);
289 37833 : petscSetOptionsHelper(po, problem);
290 37833 : }
291 :
292 : PetscErrorCode
293 44573 : petscSetupOutput(CommandLine * cmd_line)
294 : {
295 : PetscFunctionBegin;
296 44573 : char code[10] = {45, 45, 109, 111, 111, 115, 101};
297 44573 : const std::vector<std::string> argv = cmd_line->getArguments();
298 510217 : for (const auto & arg : argv)
299 : {
300 465644 : if (arg.compare(code) == 0)
301 : {
302 0 : Console::petscSetupOutput();
303 0 : break;
304 : }
305 : }
306 44573 : PetscFunctionReturn(PETSC_SUCCESS);
307 44573 : }
308 :
309 : PetscErrorCode
310 747226 : petscNonlinearConverged(SNES /*snes*/,
311 : PetscInt it,
312 : PetscReal /*xnorm*/,
313 : PetscReal /*snorm*/,
314 : PetscReal /*fnorm*/,
315 : SNESConvergedReason * reason,
316 : void * ctx)
317 : {
318 : PetscFunctionBegin;
319 747226 : FEProblemBase & problem = *static_cast<FEProblemBase *>(ctx);
320 :
321 : // execute objects that may be used in convergence check
322 747226 : problem.execute(EXEC_NONLINEAR_CONVERGENCE);
323 :
324 : // perform the convergence check
325 : Convergence::MooseConvergenceStatus status;
326 747226 : if (problem.getFailNextNonlinearConvergenceCheck())
327 : {
328 203 : status = Convergence::MooseConvergenceStatus::DIVERGED;
329 203 : problem.resetFailNextNonlinearConvergenceCheck();
330 : }
331 : else
332 : {
333 747023 : auto & convergence = problem.getConvergence(
334 747023 : problem.getNonlinearConvergenceNames()[problem.currentNonlinearSystem().number()]);
335 747023 : status = convergence.checkConvergence(it);
336 : }
337 :
338 : // convert convergence status to PETSc converged reason
339 747222 : switch (status)
340 : {
341 466681 : case Convergence::MooseConvergenceStatus::ITERATING:
342 466681 : *reason = SNES_CONVERGED_ITERATING;
343 466681 : break;
344 :
345 280286 : case Convergence::MooseConvergenceStatus::CONVERGED:
346 280286 : *reason = SNES_CONVERGED_FNORM_ABS;
347 280286 : break;
348 :
349 255 : case Convergence::MooseConvergenceStatus::DIVERGED:
350 255 : *reason = SNES_DIVERGED_DTOL;
351 255 : break;
352 : }
353 :
354 747222 : PetscFunctionReturn(PETSC_SUCCESS);
355 : }
356 :
357 : PetscErrorCode
358 836 : petscLinearConverged(
359 : KSP /*ksp*/, PetscInt it, PetscReal /*norm*/, KSPConvergedReason * reason, void * ctx)
360 : {
361 : PetscFunctionBegin;
362 836 : FEProblemBase & problem = *static_cast<FEProblemBase *>(ctx);
363 :
364 : // execute objects that may be used in convergence check
365 : // Right now, setting objects to execute on this flag would be ignored except in the
366 : // linear-system-only use case.
367 836 : problem.execute(EXEC_LINEAR_CONVERGENCE);
368 :
369 : // perform the convergence check
370 : Convergence::MooseConvergenceStatus status;
371 836 : if (problem.getFailNextSystemConvergenceCheck())
372 : {
373 0 : status = Convergence::MooseConvergenceStatus::DIVERGED;
374 0 : problem.resetFailNextSystemConvergenceCheck();
375 : }
376 : else
377 : {
378 836 : auto & convergence = problem.getConvergence(
379 836 : problem.getLinearConvergenceNames()[problem.currentLinearSystem().number()]);
380 836 : status = convergence.checkConvergence(it);
381 : }
382 :
383 : // convert convergence status to PETSc converged reason
384 836 : switch (status)
385 : {
386 748 : case Convergence::MooseConvergenceStatus::ITERATING:
387 748 : *reason = KSP_CONVERGED_ITERATING;
388 748 : break;
389 :
390 : // TODO: find a KSP code that works better for this case
391 88 : case Convergence::MooseConvergenceStatus::CONVERGED:
392 88 : *reason = KSP_CONVERGED_RTOL_NORMAL;
393 88 : break;
394 :
395 0 : case Convergence::MooseConvergenceStatus::DIVERGED:
396 0 : *reason = KSP_DIVERGED_DTOL;
397 0 : break;
398 : }
399 :
400 836 : PetscFunctionReturn(PETSC_SUCCESS);
401 : }
402 :
403 : PCSide
404 0 : getPetscPCSide(Moose::PCSideType pcs)
405 : {
406 0 : switch (pcs)
407 : {
408 0 : case Moose::PCS_LEFT:
409 0 : return PC_LEFT;
410 0 : case Moose::PCS_RIGHT:
411 0 : return PC_RIGHT;
412 0 : case Moose::PCS_SYMMETRIC:
413 0 : return PC_SYMMETRIC;
414 0 : default:
415 0 : mooseError("Unknown PC side requested.");
416 : break;
417 : }
418 : }
419 :
420 : KSPNormType
421 370226 : getPetscKSPNormType(Moose::MooseKSPNormType kspnorm)
422 : {
423 370226 : switch (kspnorm)
424 : {
425 0 : case Moose::KSPN_NONE:
426 0 : return KSP_NORM_NONE;
427 0 : case Moose::KSPN_PRECONDITIONED:
428 0 : return KSP_NORM_PRECONDITIONED;
429 370226 : case Moose::KSPN_UNPRECONDITIONED:
430 370226 : return KSP_NORM_UNPRECONDITIONED;
431 0 : case Moose::KSPN_NATURAL:
432 0 : return KSP_NORM_NATURAL;
433 0 : case Moose::KSPN_DEFAULT:
434 0 : return KSP_NORM_DEFAULT;
435 0 : default:
436 0 : mooseError("Unknown KSP norm type requested.");
437 : break;
438 : }
439 : }
440 :
441 : void
442 346426 : petscSetDefaultKSPNormType(FEProblemBase & problem, KSP ksp)
443 : {
444 716652 : for (const auto i : make_range(problem.numSolverSystems()))
445 : {
446 370226 : SolverSystem & sys = problem.getSolverSystem(i);
447 370226 : LibmeshPetscCallA(problem.comm().get(),
448 : KSPSetNormType(ksp, getPetscKSPNormType(sys.getMooseKSPNormType())));
449 : }
450 346426 : }
451 :
452 : void
453 346426 : petscSetDefaultPCSide(FEProblemBase & problem, KSP ksp)
454 : {
455 716652 : for (const auto i : make_range(problem.numSolverSystems()))
456 : {
457 370226 : SolverSystem & sys = problem.getSolverSystem(i);
458 :
459 : // PETSc 3.2.x+
460 370226 : if (sys.getPCSide() != Moose::PCS_DEFAULT)
461 0 : LibmeshPetscCallA(problem.comm().get(), KSPSetPCSide(ksp, getPetscPCSide(sys.getPCSide())));
462 : }
463 346426 : }
464 :
465 : void
466 345248 : petscSetKSPDefaults(FEProblemBase & problem, KSP ksp)
467 : {
468 345248 : auto & es = problem.es();
469 :
470 345248 : PetscReal rtol = es.parameters.get<Real>("linear solver tolerance");
471 345248 : PetscReal atol = es.parameters.get<Real>("linear solver absolute tolerance");
472 :
473 : // MOOSE defaults this to -1 for some dumb reason
474 345248 : if (atol < 0)
475 0 : atol = 1e-50;
476 :
477 345248 : PetscReal maxits = es.parameters.get<unsigned int>("linear solver maximum iterations");
478 :
479 : // 1e100 is because we don't use divtol currently
480 345248 : LibmeshPetscCallA(problem.comm().get(), KSPSetTolerances(ksp, rtol, atol, 1e100, maxits));
481 :
482 345248 : petscSetDefaultPCSide(problem, ksp);
483 :
484 345248 : petscSetDefaultKSPNormType(problem, ksp);
485 345248 : }
486 :
487 : void
488 337956 : petscSetDefaults(FEProblemBase & problem)
489 : {
490 : // We care about both nonlinear and linear systems when setting the SNES prefix because
491 : // SNESSetOptionsPrefix will also set its KSP prefix which could compete with linear system KSPs
492 682904 : for (const auto nl_index : make_range(problem.numNonlinearSystems()))
493 : {
494 344948 : NonlinearSystemBase & nl = problem.getNonlinearSystemBase(nl_index);
495 :
496 : //
497 : // prefix system matrices
498 : //
499 :
500 344948 : auto & lm_sys = nl.system();
501 : // This check is necessary because we still have at least the system matrix lying around even
502 : // when doing matrix-free, but critically even though the libmesh object(s) exist, the wrapped
503 : // PETSc Mat(s) do not
504 344948 : if (problem.solverParams(nl_index)._type != Moose::ST_JFNK)
505 690548 : for (auto & [mat_name, mat] : as_range(lm_sys.matrices_begin(), lm_sys.matrices_end()))
506 : {
507 345721 : libmesh_ignore(mat_name);
508 345721 : if (auto * const petsc_mat = dynamic_cast<PetscMatrixBase<Number> *>(mat.get()); petsc_mat)
509 : {
510 345721 : LibmeshPetscCallA(nl.comm().get(),
511 : MatSetOptionsPrefix(petsc_mat->mat(), (nl.name() + "_").c_str()));
512 : // We should call this here to ensure that options from the command line are properly
513 : // applied
514 345721 : LibmeshPetscCallA(nl.comm().get(), MatSetFromOptions(petsc_mat->mat()));
515 : }
516 : }
517 :
518 : //
519 : // prefix SNES/KSP
520 : //
521 :
522 : // dig out PETSc solver
523 344948 : auto * const petsc_solver = cast_ptr<PetscNonlinearSolver<Number> *>(nl.nonlinearSolver());
524 344948 : const char * snes_prefix = nullptr;
525 344948 : std::string snes_prefix_str;
526 344948 : if (nl.system().prefix_with_name())
527 : {
528 23800 : snes_prefix_str = nl.system().prefix();
529 23800 : snes_prefix = snes_prefix_str.c_str();
530 : }
531 344948 : SNES snes = petsc_solver->snes(snes_prefix);
532 : KSP ksp;
533 344948 : LibmeshPetscCallA(nl.comm().get(), SNESGetKSP(snes, &ksp));
534 :
535 344948 : LibmeshPetscCallA(nl.comm().get(), SNESSetMaxLinearSolveFailures(snes, 1000000));
536 :
537 344948 : LibmeshPetscCallA(nl.comm().get(), SNESSetCheckJacobianDomainError(snes, PETSC_TRUE));
538 :
539 : // In 3.0.0, the context pointer must actually be used, and the
540 : // final argument to KSPSetConvergenceTest() is a pointer to a
541 : // routine for destroying said private data context. In this case,
542 : // we use the default context provided by PETSc in addition to
543 : // a few other tests.
544 : {
545 344948 : LibmeshPetscCallA(
546 : nl.comm().get(),
547 : SNESSetConvergenceTest(snes, petscNonlinearConverged, &problem, LIBMESH_PETSC_NULLPTR));
548 : }
549 :
550 344948 : petscSetKSPDefaults(problem, ksp);
551 344948 : }
552 :
553 342956 : for (auto sys_index : make_range(problem.numLinearSystems()))
554 : {
555 : // dig out PETSc solver
556 5000 : LinearSystem & lin_sys = problem.getLinearSystem(sys_index);
557 5000 : PetscLinearSolver<Number> * petsc_solver = dynamic_cast<PetscLinearSolver<Number> *>(
558 5000 : lin_sys.linearImplicitSystem().get_linear_solver());
559 5000 : KSP ksp = petsc_solver->ksp();
560 :
561 5000 : if (problem.hasLinearConvergenceObjects())
562 88 : LibmeshPetscCallA(
563 : lin_sys.comm().get(),
564 : KSPSetConvergenceTest(ksp, petscLinearConverged, &problem, LIBMESH_PETSC_NULLPTR));
565 :
566 : // We dont set the KSP defaults here because they seem to clash with the linear solve parameters
567 : // set in FEProblemBase::solveLinearSystem
568 : }
569 337956 : }
570 :
571 : void
572 71260 : processSingletonMooseWrappedOptions(FEProblemBase & fe_problem, const InputParameters & params)
573 : {
574 71260 : setSolveTypeFromParams(fe_problem, params);
575 71260 : setLineSearchFromParams(fe_problem, params);
576 71260 : setMFFDTypeFromParams(fe_problem, params);
577 71260 : }
578 :
579 : #define checkPrefix(prefix) \
580 : mooseAssert(prefix[0] == '-', \
581 : "Leading prefix character must be a '-'. Current prefix is '" << prefix << "'"); \
582 : mooseAssert((prefix.size() == 1) || (prefix.back() == '_'), \
583 : "Terminating prefix character must be a '_'. Current prefix is '" << prefix << "'"); \
584 : mooseAssert(MooseUtils::isAllLowercase(prefix), \
585 : "PETSc prefixes should be all lower-case. What are you, a crazy person?")
586 :
587 : void
588 71260 : storePetscOptions(FEProblemBase & fe_problem,
589 : const std::string & prefix,
590 : const ParallelParamObject & param_object)
591 : {
592 71260 : const auto & params = param_object.parameters();
593 71260 : processSingletonMooseWrappedOptions(fe_problem, params);
594 :
595 : // The parameters contained in the Action
596 71260 : const auto & petsc_options = params.get<MultiMooseEnum>("petsc_options");
597 : const auto & petsc_pair_options =
598 71260 : params.get<MooseEnumItem, std::string>("petsc_options_iname", "petsc_options_value");
599 :
600 : // A reference to the PetscOptions object that contains the settings that will be used in the
601 : // solve
602 71260 : auto & po = fe_problem.getPetscOptions();
603 :
604 : // First process the single petsc options/flags
605 71260 : addPetscFlagsToPetscOptions(petsc_options, prefix, param_object, po);
606 :
607 : // Then process the option-value pairs
608 71260 : addPetscPairsToPetscOptions(
609 71260 : petsc_pair_options, fe_problem.mesh().dimension(), prefix, param_object, po);
610 71260 : }
611 :
612 : void
613 84752 : setSolveTypeFromParams(FEProblemBase & fe_problem, const InputParameters & params)
614 : {
615 : // Note: Options set in the Preconditioner block will override those set in the Executioner block
616 84752 : if (params.isParamValid("solve_type") && !params.isParamValid("_use_eigen_value"))
617 : {
618 : // Extract the solve type
619 34917 : const std::string & solve_type = params.get<MooseEnum>("solve_type");
620 70254 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
621 35337 : fe_problem.solverParams(i)._type = Moose::stringToEnum<Moose::SolveType>(solve_type);
622 34917 : }
623 84752 : }
624 :
625 : void
626 71260 : setLineSearchFromParams(FEProblemBase & fe_problem, const InputParameters & params)
627 : {
628 : // Note: Options set in the Preconditioner block will override those set in the Executioner block
629 71260 : if (params.isParamValid("line_search"))
630 : {
631 57735 : const auto & line_search = params.get<MooseEnum>("line_search");
632 114805 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
633 57070 : if (fe_problem.solverParams(i)._line_search == Moose::LS_INVALID || line_search != "default")
634 : {
635 : Moose::LineSearchType enum_line_search =
636 56744 : Moose::stringToEnum<Moose::LineSearchType>(line_search);
637 56744 : fe_problem.solverParams(i)._line_search = enum_line_search;
638 56744 : if (enum_line_search == LS_CONTACT || enum_line_search == LS_PROJECT)
639 : {
640 0 : NonlinearImplicitSystem * nl_system = dynamic_cast<NonlinearImplicitSystem *>(
641 0 : &fe_problem.getNonlinearSystemBase(i).system());
642 0 : if (!nl_system)
643 0 : mooseError("You've requested a line search but you must be solving an EigenProblem. "
644 : "These two things are not consistent.");
645 : PetscNonlinearSolver<Real> * petsc_nonlinear_solver =
646 0 : dynamic_cast<PetscNonlinearSolver<Real> *>(nl_system->nonlinear_solver.get());
647 0 : if (!petsc_nonlinear_solver)
648 0 : mooseError("Currently the MOOSE line searches all use Petsc, so you "
649 : "must use Petsc as your non-linear solver.");
650 : petsc_nonlinear_solver->linesearch_object =
651 0 : std::make_unique<ComputeLineSearchObjectWrapper>(fe_problem);
652 : }
653 : }
654 : }
655 71260 : }
656 :
657 : void
658 71260 : setMFFDTypeFromParams(FEProblemBase & fe_problem, const InputParameters & params)
659 : {
660 71260 : if (params.isParamValid("mffd_type"))
661 : {
662 70830 : const auto & mffd_type = params.get<MooseEnum>("mffd_type");
663 141247 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
664 70417 : fe_problem.solverParams(i)._mffd_type = Moose::stringToEnum<Moose::MffdType>(mffd_type);
665 : }
666 71260 : }
667 :
668 : template <typename T>
669 : void
670 54756 : checkUserProvidedPetscOption(const T & option, const ParallelParamObject & param_object)
671 : {
672 54756 : const auto & string_option = static_cast<const std::string &>(option);
673 54756 : if (string_option[0] != '-')
674 0 : param_object.mooseError("PETSc option '", string_option, "' does not begin with '-'");
675 54756 : }
676 :
677 : void
678 71260 : addPetscFlagsToPetscOptions(const MultiMooseEnum & petsc_flags,
679 : const std::string & prefix,
680 : const ParallelParamObject & param_object,
681 : PetscOptions & po)
682 : {
683 : checkPrefix(prefix);
684 :
685 : // Update the PETSc single flags
686 73845 : for (const auto & option : petsc_flags)
687 : {
688 2585 : checkUserProvidedPetscOption(option, param_object);
689 :
690 2585 : const std::string & string_option = option.name();
691 :
692 : /**
693 : * "-log_summary" cannot be used in the input file. This option needs to be set when PETSc is
694 : * initialized
695 : * which happens before the parser is even created. We'll throw an error if somebody attempts
696 : * to add this option later.
697 : */
698 2585 : if (option == "-log_summary" || option == "-log_view")
699 0 : mooseError("The PETSc option \"-log_summary\" or \"-log_view\" can only be used on the "
700 : "command line. Please "
701 : "remove it from the input file");
702 :
703 : // Update the stored items, but do not create duplicates
704 2585 : const std::string prefixed_option = prefix + string_option.substr(1);
705 2585 : if (!po.flags.isValueSet(prefixed_option))
706 : {
707 2585 : po.flags.setAdditionalValue(prefixed_option);
708 2585 : po.user_set_options.setAdditionalValue(prefixed_option);
709 : }
710 2585 : }
711 71260 : }
712 :
713 : void
714 57735 : setConvergedReasonFlags(FEProblemBase & fe_problem, const std::string & prefix)
715 : {
716 : checkPrefix(prefix);
717 57735 : libmesh_ignore(fe_problem); // avoid unused warnings for old PETSc
718 :
719 : #if !PETSC_VERSION_LESS_THAN(3, 14, 0)
720 : // the boolean in these pairs denote whether the user has specified any of the reason flags in the
721 : // input file
722 57735 : std::array<std::string, 2> reason_flags = {{"snes_converged_reason", "ksp_converged_reason"}};
723 :
724 57735 : auto & po = fe_problem.getPetscOptions();
725 :
726 173205 : for (const auto & reason_flag : reason_flags)
727 345590 : if (!po.flags.isValueSet(prefix + reason_flag) &&
728 114650 : (std::find_if(po.pairs.begin(),
729 : po.pairs.end(),
730 162125 : [&reason_flag, &prefix](auto & pair)
731 506895 : { return pair.first == (prefix + reason_flag); }) == po.pairs.end()))
732 114650 : po.pairs.emplace_back(prefix + reason_flag, "::failed");
733 : #endif
734 57735 : }
735 :
736 : void
737 71482 : addPetscPairsToPetscOptions(
738 : const std::vector<std::pair<MooseEnumItem, std::string>> & petsc_pair_options,
739 : const unsigned int mesh_dimension,
740 : const std::string & prefix,
741 : const ParallelParamObject & param_object,
742 : PetscOptions & po)
743 : {
744 : checkPrefix(prefix);
745 :
746 : // Setup the name value pairs
747 71482 : bool boomeramg_found = false;
748 71482 : bool strong_threshold_found = false;
749 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
750 71482 : bool superlu_dist_found = false;
751 71482 : bool fact_pattern_found = false;
752 71482 : bool tiny_pivot_found = false;
753 : #endif
754 71482 : std::string pc_description = "";
755 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
756 : // If users use HMG, we would like to set
757 71482 : bool hmg_found = false;
758 71482 : bool matptap_found = false;
759 71482 : bool hmg_strong_threshold_found = false;
760 : #endif
761 71482 : std::vector<std::pair<std::string, std::string>> new_options;
762 :
763 123653 : for (const auto & [option_name, option_value] : petsc_pair_options)
764 : {
765 52171 : checkUserProvidedPetscOption(option_name, param_object);
766 :
767 52171 : new_options.clear();
768 : const std::string prefixed_option_name =
769 52171 : prefix + static_cast<const std::string &>(option_name).substr(1);
770 :
771 : // Do not add duplicate settings
772 52171 : if (auto it =
773 52171 : MooseUtils::findPair(po.pairs, po.pairs.begin(), prefixed_option_name, MooseUtils::Any);
774 52171 : it == po.pairs.end())
775 : {
776 : #if !PETSC_VERSION_LESS_THAN(3, 9, 0)
777 52059 : if (option_name == "-pc_factor_mat_solver_package")
778 150 : new_options.emplace_back(prefix + "pc_factor_mat_solver_type", option_value);
779 : #else
780 : if (option_name == "-pc_factor_mat_solver_type")
781 : new_options.push_back(prefix + "pc_factor_mat_solver_package", option_value);
782 : #endif
783 :
784 : // Look for a pc description
785 78330 : if (option_name == "-pc_type" || option_name == "-sub_pc_type" ||
786 26271 : option_name == "-pc_hypre_type")
787 45893 : pc_description += option_value + ' ';
788 :
789 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
790 52059 : if (option_name == "-pc_type" && option_value == "hmg")
791 42 : hmg_found = true;
792 :
793 : // MPIAIJ for PETSc 3.12.0: -matptap_via
794 : // MAIJ for PETSc 3.12.0: -matmaijptap_via
795 : // MPIAIJ for PETSc 3.13 to 3.16: -matptap_via, -matproduct_ptap_via
796 : // MAIJ for PETSc 3.13 to 3.16: -matproduct_ptap_via
797 : // MPIAIJ for PETSc 3.17 and higher: -matptap_via, -mat_product_algorithm
798 : // MAIJ for PETSc 3.17 and higher: -mat_product_algorithm
799 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
800 52185 : if (hmg_found && (option_name == "-matptap_via" || option_name == "-matmaijptap_via" ||
801 126 : option_name == "-matproduct_ptap_via"))
802 0 : new_options.emplace_back(prefix + "mat_product_algorithm", option_value);
803 : #elif !PETSC_VERSION_LESS_THAN(3, 13, 0)
804 : if (hmg_found && (option_name == "-matptap_via" || option_name == "-matmaijptap_via"))
805 : new_options.emplace_back(prefix + "matproduct_ptap_via", option_value);
806 : #else
807 : if (hmg_found && (option_name == "-matproduct_ptap_via"))
808 : {
809 : new_options.emplace_back(prefix + "matptap_via", option_value);
810 : new_options.emplace_back(prefix + "matmaijptap_via", option_value);
811 : }
812 : #endif
813 :
814 104118 : if (option_name == "-matptap_via" || option_name == "-matmaijptap_via" ||
815 104118 : option_name == "-matproduct_ptap_via" || option_name == "-mat_product_algorithm")
816 0 : matptap_found = true;
817 :
818 : // For 3D problems, we need to set this 0.7
819 52059 : if (option_name == "-hmg_inner_pc_hypre_boomeramg_strong_threshold")
820 0 : hmg_strong_threshold_found = true;
821 : #endif
822 : // This special case is common enough that we'd like to handle it for the user.
823 52059 : if (option_name == "-pc_hypre_type" && option_value == "boomeramg")
824 20105 : boomeramg_found = true;
825 52059 : if (option_name == "-pc_hypre_boomeramg_strong_threshold")
826 0 : strong_threshold_found = true;
827 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
828 52059 : if ((option_name == "-pc_factor_mat_solver_package" ||
829 52978 : option_name == "-pc_factor_mat_solver_type") &&
830 919 : option_value == "superlu_dist")
831 303 : superlu_dist_found = true;
832 52059 : if (option_name == "-mat_superlu_dist_fact")
833 0 : fact_pattern_found = true;
834 52059 : if (option_name == "-mat_superlu_dist_replacetinypivot")
835 0 : tiny_pivot_found = true;
836 : #endif
837 :
838 52059 : if (!new_options.empty())
839 : {
840 150 : std::copy(new_options.begin(), new_options.end(), std::back_inserter(po.pairs));
841 300 : for (const auto & option : new_options)
842 150 : po.user_set_options.setAdditionalValue(option.first);
843 : }
844 : else
845 : {
846 51909 : po.pairs.push_back(std::make_pair(prefixed_option_name, option_value));
847 51909 : po.user_set_options.setAdditionalValue(prefixed_option_name);
848 : }
849 : }
850 : else
851 : {
852 : do
853 : {
854 112 : it->second = option_value;
855 112 : it = MooseUtils::findPair(po.pairs, std::next(it), prefixed_option_name, MooseUtils::Any);
856 112 : } while (it != po.pairs.end());
857 : }
858 52171 : }
859 :
860 : // When running a 3D mesh with boomeramg, it is almost always best to supply a strong threshold
861 : // value. We will provide that for the user here if they haven't supplied it themselves.
862 71482 : if (boomeramg_found && !strong_threshold_found && mesh_dimension == 3)
863 : {
864 894 : po.pairs.emplace_back(prefix + "pc_hypre_boomeramg_strong_threshold", "0.7");
865 894 : pc_description += "strong_threshold: 0.7 (auto)";
866 : }
867 :
868 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
869 71482 : if (hmg_found && !hmg_strong_threshold_found && mesh_dimension == 3)
870 : {
871 28 : po.pairs.emplace_back(prefix + "hmg_inner_pc_hypre_boomeramg_strong_threshold", "0.7");
872 28 : pc_description += "strong_threshold: 0.7 (auto)";
873 : }
874 :
875 : // Default PETSc PtAP takes too much memory, and it is not quite useful
876 : // Let us switch to use new algorithm
877 71482 : if (hmg_found && !matptap_found)
878 : {
879 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
880 42 : po.pairs.emplace_back(prefix + "mat_product_algorithm", "allatonce");
881 : #elif !PETSC_VERSION_LESS_THAN(3, 13, 0)
882 : po.pairs.emplace_back(prefix + "matproduct_ptap_via", "allatonce");
883 : #else
884 : po.pairs.emplace_back(prefix + "matptap_via", "allatonce");
885 : po.pairs.emplace_back(prefix + "matmaijptap_via", "allatonce");
886 : #endif
887 : }
888 : #endif
889 :
890 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
891 : // In PETSc-3.7.{0--4}, there is a bug when using superlu_dist, and we have to use
892 : // SamePattern_SameRowPerm, otherwise we use whatever we have in PETSc
893 71482 : if (superlu_dist_found && !fact_pattern_found)
894 : {
895 303 : po.pairs.emplace_back(prefix + "mat_superlu_dist_fact",
896 : #if PETSC_VERSION_LESS_THAN(3, 7, 5)
897 : "SamePattern_SameRowPerm");
898 : pc_description += "mat_superlu_dist_fact: SamePattern_SameRowPerm ";
899 : #else
900 : "SamePattern");
901 303 : pc_description += "mat_superlu_dist_fact: SamePattern ";
902 : #endif
903 : }
904 :
905 : // restore this superlu option
906 71482 : if (superlu_dist_found && !tiny_pivot_found)
907 : {
908 303 : po.pairs.emplace_back(prefix + "mat_superlu_dist_replacetinypivot", "1");
909 303 : pc_description += " mat_superlu_dist_replacetinypivot: true ";
910 : }
911 : #endif
912 : // Set Preconditioner description
913 71482 : if (!pc_description.empty() && prefix.size() > 1)
914 634 : po.pc_description += "[" + prefix.substr(1, prefix.size() - 2) + "]: ";
915 71482 : po.pc_description += pc_description;
916 71482 : }
917 :
918 : std::set<std::string>
919 286959 : getPetscValidLineSearches()
920 : {
921 2582631 : return {"default", "shell", "none", "basic", "l2", "bt", "cp"};
922 573918 : }
923 :
924 : InputParameters
925 399569 : getPetscValidParams()
926 : {
927 399569 : InputParameters params = emptyInputParameters();
928 :
929 399569 : MooseEnum solve_type("PJFNK JFNK NEWTON FD LINEAR");
930 399569 : params.addParam<MooseEnum>("solve_type",
931 : solve_type,
932 : "PJFNK: Preconditioned Jacobian-Free Newton Krylov "
933 : "JFNK: Jacobian-Free Newton Krylov "
934 : "NEWTON: Full Newton Solve "
935 : "FD: Use finite differences to compute Jacobian "
936 : "LINEAR: Solving a linear problem");
937 :
938 399569 : MooseEnum mffd_type("wp ds", "wp");
939 399569 : params.addParam<MooseEnum>("mffd_type",
940 : mffd_type,
941 : "Specifies the finite differencing type for "
942 : "Jacobian-free solve types. Note that the "
943 : "default is wp (for Walker and Pernice).");
944 :
945 1198707 : params.addParam<MultiMooseEnum>(
946 799138 : "petsc_options", getCommonPetscFlags(), "Singleton PETSc options");
947 1198707 : params.addParam<MultiMooseEnum>(
948 799138 : "petsc_options_iname", getCommonPetscKeys(), "Names of PETSc name/value pairs");
949 399569 : params.addParam<std::vector<std::string>>(
950 : "petsc_options_value",
951 : "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\"");
952 399569 : params.addParamNamesToGroup("solve_type petsc_options petsc_options_iname petsc_options_value "
953 : "mffd_type",
954 : "PETSc");
955 :
956 799138 : return params;
957 399569 : }
958 :
959 : MultiMooseEnum
960 415791 : getCommonSNESFlags()
961 : {
962 : return MultiMooseEnum(
963 : "-ksp_monitor_snes_lg -snes_ksp_ew -ksp_snes_ew -snes_converged_reason "
964 : "-snes_ksp -snes_linesearch_monitor -snes_mf -snes_mf_operator -snes_monitor "
965 : "-snes_test_display -snes_view -snes_monitor_cancel",
966 : "",
967 415791 : true);
968 : }
969 :
970 : MultiMooseEnum
971 415165 : getCommonKSPFlags()
972 : {
973 : return MultiMooseEnum(
974 415165 : "-ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_monitor", "", true);
975 : }
976 :
977 : MultiMooseEnum
978 415165 : getCommonPetscFlags()
979 : {
980 415165 : auto options = MultiMooseEnum("-dm_moose_print_embedding -dm_view", "", true);
981 415165 : options.addValidName(getCommonKSPFlags());
982 415165 : options.addValidName(getCommonSNESFlags());
983 415165 : return options;
984 0 : }
985 :
986 : MultiMooseEnum
987 415791 : getCommonSNESKeys()
988 : {
989 : return MultiMooseEnum("-snes_atol -snes_linesearch_type -snes_ls -snes_max_it -snes_rtol "
990 : "-snes_divergence_tolerance -snes_type",
991 : "",
992 415791 : true);
993 : }
994 :
995 : MultiMooseEnum
996 415165 : getCommonKSPKeys()
997 : {
998 : return MultiMooseEnum("-ksp_atol -ksp_gmres_restart -ksp_max_it -ksp_pc_side -ksp_rtol "
999 : "-ksp_type -sub_ksp_type",
1000 : "",
1001 415165 : true);
1002 : }
1003 : MultiMooseEnum
1004 415165 : getCommonPetscKeys()
1005 : {
1006 : auto options = MultiMooseEnum("-mat_fd_coloring_err -mat_fd_type -mat_mffd_type "
1007 : "-pc_asm_overlap -pc_factor_levels "
1008 : "-pc_factor_mat_ordering_type -pc_hypre_boomeramg_grid_sweeps_all "
1009 : "-pc_hypre_boomeramg_max_iter "
1010 : "-pc_hypre_boomeramg_strong_threshold -pc_hypre_type -pc_type "
1011 : "-sub_pc_type",
1012 : "",
1013 415165 : true);
1014 415165 : options.addValidName(getCommonKSPKeys());
1015 415165 : options.addValidName(getCommonSNESKeys());
1016 415165 : return options;
1017 0 : }
1018 :
1019 : bool
1020 564 : isSNESVI(FEProblemBase & fe_problem)
1021 : {
1022 564 : const PetscOptions & petsc = fe_problem.getPetscOptions();
1023 :
1024 : int argc;
1025 : char ** args;
1026 564 : LibmeshPetscCallA(fe_problem.comm().get(), PetscGetArgs(&argc, &args));
1027 :
1028 564 : std::vector<std::string> cml_arg;
1029 6826 : for (int i = 0; i < argc; i++)
1030 6262 : cml_arg.push_back(args[i]);
1031 :
1032 564 : if (MooseUtils::findPair(petsc.pairs, petsc.pairs.begin(), MooseUtils::Any, "vinewtonssls") ==
1033 1128 : petsc.pairs.end() &&
1034 564 : MooseUtils::findPair(petsc.pairs, petsc.pairs.begin(), MooseUtils::Any, "vinewtonrsls") ==
1035 608 : petsc.pairs.end() &&
1036 1172 : std::find(cml_arg.begin(), cml_arg.end(), "vinewtonssls") == cml_arg.end() &&
1037 608 : std::find(cml_arg.begin(), cml_arg.end(), "vinewtonrsls") == cml_arg.end())
1038 44 : return false;
1039 :
1040 520 : return true;
1041 564 : }
1042 :
1043 : void
1044 215447 : setSinglePetscOption(const std::string & name,
1045 : const std::string & value /*=""*/,
1046 : FEProblemBase * const problem /*=nullptr*/)
1047 : {
1048 215447 : static const TIMPI::Communicator comm_world(PETSC_COMM_WORLD);
1049 215447 : const TIMPI::Communicator & comm = problem ? problem->comm() : comm_world;
1050 215447 : LibmeshPetscCallA(comm.get(),
1051 : PetscOptionsSetValue(LIBMESH_PETSC_NULLPTR,
1052 : name.c_str(),
1053 : value == "" ? LIBMESH_PETSC_NULLPTR : value.c_str()));
1054 215447 : const auto lower_case_name = MooseUtils::toLower(name);
1055 20 : auto check_problem = [problem, &lower_case_name]()
1056 : {
1057 20 : if (!problem)
1058 0 : mooseError(
1059 : "Setting the option '",
1060 : lower_case_name,
1061 : "' requires passing a 'problem' parameter. Contact a developer of your application "
1062 : "to have them update their code. If in doubt, reach out to the MOOSE team on Github "
1063 : "discussions");
1064 215467 : };
1065 :
1066 : // Select vector type from user-passed PETSc options
1067 215447 : if (lower_case_name.find("-vec_type") != std::string::npos)
1068 : {
1069 10 : check_problem();
1070 20 : for (auto & solver_system : problem->_solver_systems)
1071 : {
1072 10 : auto & lm_sys = solver_system->system();
1073 30 : for (auto & [vec_name, vec] : as_range(lm_sys.vectors_begin(), lm_sys.vectors_end()))
1074 : {
1075 20 : libmesh_ignore(vec_name);
1076 20 : auto * const petsc_vec = cast_ptr<PetscVector<Number> *>(vec.get());
1077 20 : LibmeshPetscCallA(comm.get(), VecSetFromOptions(petsc_vec->vec()));
1078 : }
1079 : // The solution vectors aren't included in the system vectors storage
1080 10 : auto * petsc_vec = cast_ptr<PetscVector<Number> *>(lm_sys.solution.get());
1081 10 : LibmeshPetscCallA(comm.get(), VecSetFromOptions(petsc_vec->vec()));
1082 10 : petsc_vec = cast_ptr<PetscVector<Number> *>(lm_sys.current_local_solution.get());
1083 10 : LibmeshPetscCallA(comm.get(), VecSetFromOptions(petsc_vec->vec()));
1084 : }
1085 : }
1086 : // Select matrix type from user-passed PETSc options
1087 215437 : else if (lower_case_name.find("mat_type") != std::string::npos)
1088 : {
1089 10 : check_problem();
1090 :
1091 10 : bool found_matching_prefix = false;
1092 10 : for (const auto i : index_range(problem->_solver_systems))
1093 : {
1094 10 : const auto & solver_sys_name = problem->_solver_sys_names[i];
1095 10 : if (lower_case_name.find("-" + MooseUtils::toLower(solver_sys_name) + "_mat_type") ==
1096 : std::string::npos)
1097 0 : continue;
1098 :
1099 10 : if (problem->solverParams(i)._type == Moose::ST_JFNK)
1100 0 : mooseError(
1101 : "Setting option '", lower_case_name, "' is incompatible with a JFNK 'solve_type'");
1102 :
1103 10 : auto & lm_sys = problem->_solver_systems[i]->system();
1104 20 : for (auto & [mat_name, mat] : as_range(lm_sys.matrices_begin(), lm_sys.matrices_end()))
1105 : {
1106 10 : libmesh_ignore(mat_name);
1107 10 : if (auto * const petsc_mat = dynamic_cast<PetscMatrixBase<Number> *>(mat.get()); petsc_mat)
1108 : {
1109 : #ifdef DEBUG
1110 : const char * mat_prefix;
1111 : LibmeshPetscCallA(comm.get(), MatGetOptionsPrefix(petsc_mat->mat(), &mat_prefix));
1112 : mooseAssert(strcmp(mat_prefix, (solver_sys_name + "_").c_str()) == 0,
1113 : "We should have prefixed these matrices previously");
1114 : #endif
1115 10 : LibmeshPetscCallA(comm.get(), MatSetFromOptions(petsc_mat->mat()));
1116 : }
1117 : }
1118 10 : found_matching_prefix = true;
1119 10 : break;
1120 : }
1121 :
1122 10 : if (!found_matching_prefix)
1123 0 : mooseError("We did not find a matching solver system name for the petsc option '",
1124 : lower_case_name,
1125 : "'");
1126 : }
1127 215447 : }
1128 :
1129 : void
1130 66219 : setSinglePetscOptionIfAppropriate(const MultiMooseEnum & dont_add_these_options,
1131 : const std::string & name,
1132 : const std::string & value /*=""*/,
1133 : FEProblemBase * const problem /*=nullptr*/)
1134 : {
1135 66219 : if (!dont_add_these_options.contains(name))
1136 64865 : setSinglePetscOption(name, value, problem);
1137 66219 : }
1138 :
1139 : void
1140 0 : colorAdjacencyMatrix(PetscScalar * adjacency_matrix,
1141 : unsigned int size,
1142 : unsigned int colors,
1143 : std::vector<unsigned int> & vertex_colors,
1144 : const char * coloring_algorithm)
1145 : {
1146 : // Mat A will be a dense matrix from the incoming data structure
1147 : Mat A;
1148 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatCreate(PETSC_COMM_SELF, &A));
1149 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSetSizes(A, size, size, size, size));
1150 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSetType(A, MATSEQDENSE));
1151 : // PETSc requires a non-const data array to populate the matrix
1152 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSeqDenseSetPreallocation(A, adjacency_matrix));
1153 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1154 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1155 :
1156 : // Convert A to a sparse matrix
1157 : #if PETSC_VERSION_LESS_THAN(3, 7, 0)
1158 : LibmeshPetscCallA(PETSC_COMM_SELF, MatConvert(A, MATAIJ, MAT_REUSE_MATRIX, &A));
1159 : #else
1160 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A));
1161 : #endif
1162 :
1163 : ISColoring iscoloring;
1164 : MatColoring mc;
1165 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringCreate(A, &mc));
1166 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetType(mc, coloring_algorithm));
1167 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetMaxColors(mc, static_cast<PetscInt>(colors)));
1168 :
1169 : // Petsc normally colors by distance two (neighbors of neighbors), we just want one
1170 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetDistance(mc, 1));
1171 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetFromOptions(mc));
1172 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringApply(mc, &iscoloring));
1173 :
1174 : PetscInt nn;
1175 : IS * is;
1176 : #if PETSC_RELEASE_LESS_THAN(3, 12, 0)
1177 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringGetIS(iscoloring, &nn, &is));
1178 : #else
1179 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &nn, &is));
1180 : #endif
1181 :
1182 0 : if (nn > static_cast<PetscInt>(colors))
1183 0 : throw std::runtime_error("Not able to color with designated number of colors");
1184 :
1185 0 : for (int i = 0; i < nn; i++)
1186 : {
1187 : PetscInt isize;
1188 : const PetscInt * indices;
1189 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISGetLocalSize(is[i], &isize));
1190 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISGetIndices(is[i], &indices));
1191 0 : for (int j = 0; j < isize; j++)
1192 : {
1193 : mooseAssert(indices[j] < static_cast<PetscInt>(vertex_colors.size()), "Index out of bounds");
1194 0 : vertex_colors[indices[j]] = i;
1195 : }
1196 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISRestoreIndices(is[i], &indices));
1197 : }
1198 :
1199 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatDestroy(&A));
1200 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringDestroy(&mc));
1201 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringDestroy(&iscoloring));
1202 0 : }
1203 :
1204 : void
1205 3580 : dontAddPetscFlag(const std::string & flag, PetscOptions & petsc_options)
1206 : {
1207 3580 : if (!petsc_options.dont_add_these_options.contains(flag))
1208 3580 : petsc_options.dont_add_these_options.setAdditionalValue(flag);
1209 3580 : }
1210 :
1211 : void
1212 627 : dontAddNonlinearConvergedReason(FEProblemBase & fe_problem)
1213 : {
1214 627 : dontAddPetscFlag("-snes_converged_reason", fe_problem.getPetscOptions());
1215 627 : }
1216 :
1217 : void
1218 637 : dontAddLinearConvergedReason(FEProblemBase & fe_problem)
1219 : {
1220 637 : dontAddPetscFlag("-ksp_converged_reason", fe_problem.getPetscOptions());
1221 637 : }
1222 :
1223 : void
1224 0 : dontAddCommonKSPOptions(FEProblemBase & fe_problem)
1225 : {
1226 0 : auto & petsc_options = fe_problem.getPetscOptions();
1227 0 : for (const auto & flag : getCommonKSPFlags().getNames())
1228 0 : dontAddPetscFlag(flag, petsc_options);
1229 0 : for (const auto & key : getCommonKSPKeys().getNames())
1230 0 : dontAddPetscFlag(key, petsc_options);
1231 0 : }
1232 :
1233 : void
1234 626 : dontAddCommonSNESOptions(FEProblemBase & fe_problem)
1235 : {
1236 626 : auto & petsc_options = fe_problem.getPetscOptions();
1237 8138 : for (const auto & flag : getCommonSNESFlags().getNames())
1238 7512 : if (!petsc_options.dont_add_these_options.contains(flag))
1239 4382 : petsc_options.dont_add_these_options.setAdditionalValue(flag);
1240 5008 : for (const auto & key : getCommonSNESKeys().getNames())
1241 4382 : if (!petsc_options.dont_add_these_options.contains(key))
1242 2817 : petsc_options.dont_add_these_options.setAdditionalValue(key);
1243 626 : }
1244 :
1245 : std::unique_ptr<PetscMatrix<Number>>
1246 44 : createMatrixFromFile(const libMesh::Parallel::Communicator & comm,
1247 : Mat & mat,
1248 : const std::string & binary_mat_file,
1249 : const unsigned int mat_number_to_load)
1250 : {
1251 44 : LibmeshPetscCallA(comm.get(), MatCreate(comm.get(), &mat));
1252 : PetscViewer matviewer;
1253 44 : LibmeshPetscCallA(
1254 : comm.get(),
1255 : PetscViewerBinaryOpen(comm.get(), binary_mat_file.c_str(), FILE_MODE_READ, &matviewer));
1256 88 : for (unsigned int i = 0; i < mat_number_to_load; ++i)
1257 44 : LibmeshPetscCallA(comm.get(), MatLoad(mat, matviewer));
1258 44 : LibmeshPetscCallA(comm.get(), PetscViewerDestroy(&matviewer));
1259 :
1260 88 : return std::make_unique<PetscMatrix<Number>>(mat, comm);
1261 : }
1262 :
1263 : } // Namespace PetscSupport
1264 : } // Namespace MOOSE
|
