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 3001 : stringify(const LineSearchType & t)
102 : {
103 3001 : switch (t)
104 : {
105 2906 : case LS_BASIC:
106 2906 : 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 23 : case LS_BT:
116 23 : return "bt";
117 72 : case LS_CP:
118 72 : 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 29067 : stringify(const MffdType & t)
131 : {
132 29067 : switch (t)
133 : {
134 29057 : case MFFD_WP:
135 29057 : 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 42863 : setSolverOptions(const SolverParams & solver_params, const MultiMooseEnum & dont_add_these_options)
146 : {
147 : // set PETSc options implied by a solve type
148 42863 : switch (solver_params._type)
149 : {
150 28945 : case Moose::ST_PJFNK:
151 86835 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
152 57890 : solver_params._prefix + "snes_mf_operator");
153 28945 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
154 57890 : solver_params._prefix + "mat_mffd_type",
155 57890 : stringify(solver_params._mffd_type));
156 28945 : break;
157 :
158 122 : case Moose::ST_JFNK:
159 122 : setSinglePetscOptionIfAppropriate(dont_add_these_options, solver_params._prefix + "snes_mf");
160 122 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
161 244 : solver_params._prefix + "mat_mffd_type",
162 244 : stringify(solver_params._mffd_type));
163 122 : break;
164 :
165 11504 : case Moose::ST_NEWTON:
166 11504 : 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 2292 : case Moose::ST_LINEAR:
173 6876 : setSinglePetscOptionIfAppropriate(
174 4584 : dont_add_these_options, solver_params._prefix + "snes_type", "ksponly");
175 6876 : setSinglePetscOptionIfAppropriate(dont_add_these_options,
176 4584 : solver_params._prefix + "snes_monitor_cancel");
177 2292 : break;
178 : }
179 :
180 42863 : Moose::LineSearchType ls_type = solver_params._line_search;
181 42863 : if (ls_type == Moose::LS_NONE)
182 2752 : ls_type = Moose::LS_BASIC;
183 :
184 42863 : if (ls_type != Moose::LS_DEFAULT && ls_type != Moose::LS_CONTACT && ls_type != Moose::LS_PROJECT)
185 3001 : setSinglePetscOptionIfAppropriate(
186 6002 : dont_add_these_options, solver_params._prefix + "snes_linesearch_type", stringify(ls_type));
187 42863 : }
188 :
189 : void
190 105 : petscSetupDM(NonlinearSystemBase & nl, const std::string & dm_name)
191 : {
192 : PetscBool ismoose;
193 105 : 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 105 : 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 105 : dynamic_cast<PetscNonlinearSolver<Number> *>(nl.nonlinearSolver());
201 105 : const char * snes_prefix = nullptr;
202 105 : std::string snes_prefix_str;
203 105 : 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 105 : SNES snes = petsc_solver->snes(snes_prefix);
209 : // if there exists a DMMoose object, not to recreate a new one
210 105 : LibmeshPetscCallA(nl.comm().get(), SNESGetDM(snes, &dm));
211 105 : if (dm)
212 : {
213 105 : LibmeshPetscCallA(nl.comm().get(), PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose));
214 105 : if (ismoose)
215 0 : return;
216 : }
217 105 : LibmeshPetscCallA(nl.comm().get(), DMCreateMoose(nl.comm().get(), nl, dm_name, &dm));
218 105 : LibmeshPetscCallA(nl.comm().get(), DMSetFromOptions(dm));
219 105 : LibmeshPetscCallA(nl.comm().get(), DMSetUp(dm));
220 105 : LibmeshPetscCallA(nl.comm().get(), SNESSetDM(snes, dm));
221 105 : 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 105 : }
231 :
232 : void
233 43211 : 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 43211 : 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 43211 : LibmeshPetscCallA(PETSC_COMM_WORLD,
246 : PetscOptionsInsert(LIBMESH_PETSC_NULLPTR, &argc, &args, NULL));
247 : #endif
248 : }
249 43211 : }
250 :
251 : void
252 42653 : petscSetOptionsHelper(const PetscOptions & po, FEProblemBase * const problem)
253 : {
254 : // Add any additional options specified in the input file
255 45251 : 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 2598 : if (!po.dont_add_these_options.contains(flag.name()) ||
259 0 : po.user_set_options.contains(flag.name()))
260 2598 : setSinglePetscOption(flag.rawName().c_str());
261 :
262 : // Add option pairs
263 180554 : for (auto & option : po.pairs)
264 139424 : if (!po.dont_add_these_options.contains(option.first) ||
265 1523 : po.user_set_options.contains(option.first))
266 136388 : setSinglePetscOption(option.first, option.second, problem);
267 :
268 42653 : addPetscOptionsFromCommandline();
269 42653 : }
270 :
271 : void
272 1815 : petscSetOptions(const PetscOptions & po,
273 : const SolverParams & solver_params,
274 : FEProblemBase * const problem)
275 : {
276 1815 : PetscCallAbort(PETSC_COMM_WORLD, PetscOptionsClear(LIBMESH_PETSC_NULLPTR));
277 1815 : setSolverOptions(solver_params, po.dont_add_these_options);
278 1815 : petscSetOptionsHelper(po, problem);
279 1815 : }
280 :
281 : void
282 40838 : petscSetOptions(const PetscOptions & po,
283 : const std::vector<SolverParams> & solver_params_vec,
284 : FEProblemBase * const problem)
285 : {
286 40838 : PetscCallAbort(PETSC_COMM_WORLD, PetscOptionsClear(LIBMESH_PETSC_NULLPTR));
287 81886 : for (const auto & solver_params : solver_params_vec)
288 41048 : setSolverOptions(solver_params, po.dont_add_these_options);
289 40838 : petscSetOptionsHelper(po, problem);
290 40838 : }
291 :
292 : PetscErrorCode
293 48118 : petscSetupOutput(CommandLine * cmd_line)
294 : {
295 : PetscFunctionBegin;
296 48118 : char code[10] = {45, 45, 109, 111, 111, 115, 101};
297 48118 : const std::vector<std::string> argv = cmd_line->getArguments();
298 506842 : for (const auto & arg : argv)
299 : {
300 458724 : if (arg.compare(code) == 0)
301 : {
302 0 : Console::petscSetupOutput();
303 0 : break;
304 : }
305 : }
306 48118 : PetscFunctionReturn(PETSC_SUCCESS);
307 48118 : }
308 :
309 : PetscErrorCode
310 812780 : 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 812780 : FEProblemBase & problem = *static_cast<FEProblemBase *>(ctx);
320 :
321 : // execute objects that may be used in convergence check
322 812780 : problem.execute(EXEC_NONLINEAR_CONVERGENCE);
323 :
324 : // perform the convergence check
325 : Convergence::MooseConvergenceStatus status;
326 812780 : if (problem.getFailNextNonlinearConvergenceCheck())
327 : {
328 212 : status = Convergence::MooseConvergenceStatus::DIVERGED;
329 212 : problem.resetFailNextNonlinearConvergenceCheck();
330 : }
331 : else
332 : {
333 812568 : auto & convergence = problem.getConvergence(
334 812568 : problem.getNonlinearConvergenceNames()[problem.currentNonlinearSystem().number()]);
335 812568 : status = convergence.checkConvergence(it);
336 : }
337 :
338 : // convert convergence status to PETSc converged reason
339 812776 : switch (status)
340 : {
341 506175 : case Convergence::MooseConvergenceStatus::ITERATING:
342 506175 : *reason = SNES_CONVERGED_ITERATING;
343 506175 : break;
344 :
345 306337 : case Convergence::MooseConvergenceStatus::CONVERGED:
346 306337 : *reason = SNES_CONVERGED_FNORM_ABS;
347 306337 : break;
348 :
349 264 : case Convergence::MooseConvergenceStatus::DIVERGED:
350 264 : *reason = SNES_DIVERGED_DTOL;
351 264 : break;
352 : }
353 :
354 812776 : PetscFunctionReturn(PETSC_SUCCESS);
355 : }
356 :
357 : PetscErrorCode
358 2634 : petscLinearConverged(
359 : KSP /*ksp*/, PetscInt it, PetscReal /*norm*/, KSPConvergedReason * reason, void * ctx)
360 : {
361 : PetscFunctionBegin;
362 2634 : 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 2634 : problem.execute(EXEC_LINEAR_CONVERGENCE);
368 :
369 : // perform the convergence check
370 : Convergence::MooseConvergenceStatus status;
371 2634 : if (problem.getFailNextSystemConvergenceCheck())
372 : {
373 0 : status = Convergence::MooseConvergenceStatus::DIVERGED;
374 0 : problem.resetFailNextSystemConvergenceCheck();
375 : }
376 : else
377 : {
378 2634 : auto & convergence = problem.getConvergence(
379 2634 : problem.getLinearConvergenceNames()[problem.currentLinearSystem().number()]);
380 2634 : status = convergence.checkConvergence(it);
381 : }
382 :
383 : // convert convergence status to PETSc converged reason
384 2634 : switch (status)
385 : {
386 2201 : case Convergence::MooseConvergenceStatus::ITERATING:
387 2201 : *reason = KSP_CONVERGED_ITERATING;
388 2201 : break;
389 :
390 : // TODO: find a KSP code that works better for this case
391 433 : case Convergence::MooseConvergenceStatus::CONVERGED:
392 433 : *reason = KSP_CONVERGED_RTOL_NORMAL;
393 433 : break;
394 :
395 0 : case Convergence::MooseConvergenceStatus::DIVERGED:
396 0 : *reason = KSP_DIVERGED_DTOL;
397 0 : break;
398 : }
399 :
400 2634 : 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 403395 : getPetscKSPNormType(Moose::MooseKSPNormType kspnorm)
422 : {
423 403395 : 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 403395 : case Moose::KSPN_UNPRECONDITIONED:
430 403395 : 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 377409 : petscSetDefaultKSPNormType(FEProblemBase & problem, KSP ksp)
443 : {
444 780804 : for (const auto i : make_range(problem.numSolverSystems()))
445 : {
446 403395 : SolverSystem & sys = problem.getSolverSystem(i);
447 403395 : LibmeshPetscCallA(problem.comm().get(),
448 : KSPSetNormType(ksp, getPetscKSPNormType(sys.getMooseKSPNormType())));
449 : }
450 377409 : }
451 :
452 : void
453 377409 : petscSetDefaultPCSide(FEProblemBase & problem, KSP ksp)
454 : {
455 780804 : for (const auto i : make_range(problem.numSolverSystems()))
456 : {
457 403395 : SolverSystem & sys = problem.getSolverSystem(i);
458 :
459 : // PETSc 3.2.x+
460 403395 : if (sys.getPCSide() != Moose::PCS_DEFAULT)
461 0 : LibmeshPetscCallA(problem.comm().get(), KSPSetPCSide(ksp, getPetscPCSide(sys.getPCSide())));
462 : }
463 377409 : }
464 :
465 : void
466 376163 : petscSetKSPDefaults(FEProblemBase & problem, KSP ksp)
467 : {
468 376163 : auto & es = problem.es();
469 :
470 376163 : PetscReal rtol = es.parameters.get<Real>("linear solver tolerance");
471 376163 : PetscReal atol = es.parameters.get<Real>("linear solver absolute tolerance");
472 :
473 : // MOOSE defaults this to -1 for some dumb reason
474 376163 : if (atol < 0)
475 0 : atol = 1e-50;
476 :
477 376163 : PetscReal maxits = es.parameters.get<unsigned int>("linear solver maximum iterations");
478 :
479 : // 1e100 is because we don't use divtol currently
480 376163 : LibmeshPetscCallA(problem.comm().get(), KSPSetTolerances(ksp, rtol, atol, 1e100, maxits));
481 :
482 376163 : petscSetDefaultPCSide(problem, ksp);
483 :
484 376163 : petscSetDefaultKSPNormType(problem, ksp);
485 376163 : }
486 :
487 : void
488 368162 : 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 744005 : for (const auto nl_index : make_range(problem.numNonlinearSystems()))
493 : {
494 375843 : NonlinearSystemBase & nl = problem.getNonlinearSystemBase(nl_index);
495 :
496 : //
497 : // prefix system matrices
498 : //
499 :
500 375843 : 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 375843 : if (problem.solverParams(nl_index)._type != Moose::ST_JFNK)
505 752416 : for (auto & [mat_name, mat] : as_range(lm_sys.matrices_begin(), lm_sys.matrices_end()))
506 : {
507 376706 : libmesh_ignore(mat_name);
508 376706 : if (auto * const petsc_mat = dynamic_cast<PetscMatrixBase<Number> *>(mat.get()); petsc_mat)
509 : {
510 376706 : 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 376706 : LibmeshPetscCallA(nl.comm().get(), MatSetFromOptions(petsc_mat->mat()));
515 : }
516 : }
517 :
518 : //
519 : // prefix SNES/KSP
520 : //
521 :
522 : // dig out PETSc solver
523 375843 : auto * const petsc_solver = cast_ptr<PetscNonlinearSolver<Number> *>(nl.nonlinearSolver());
524 375843 : const char * snes_prefix = nullptr;
525 375843 : std::string snes_prefix_str;
526 375843 : if (nl.system().prefix_with_name())
527 : {
528 25986 : snes_prefix_str = nl.system().prefix();
529 25986 : snes_prefix = snes_prefix_str.c_str();
530 : }
531 375843 : SNES snes = petsc_solver->snes(snes_prefix);
532 : KSP ksp;
533 375843 : LibmeshPetscCallA(nl.comm().get(), SNESGetKSP(snes, &ksp));
534 :
535 375843 : LibmeshPetscCallA(nl.comm().get(), SNESSetMaxLinearSolveFailures(snes, 1000000));
536 :
537 375843 : 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 375843 : LibmeshPetscCallA(
546 : nl.comm().get(),
547 : SNESSetConvergenceTest(snes, petscNonlinearConverged, &problem, LIBMESH_PETSC_NULLPTR));
548 : }
549 :
550 375843 : petscSetKSPDefaults(problem, ksp);
551 375843 : }
552 :
553 373587 : for (auto sys_index : make_range(problem.numLinearSystems()))
554 : {
555 : // dig out PETSc solver
556 5425 : LinearSystem & lin_sys = problem.getLinearSystem(sys_index);
557 5425 : PetscLinearSolver<Number> * petsc_solver = dynamic_cast<PetscLinearSolver<Number> *>(
558 5425 : lin_sys.linearImplicitSystem().get_linear_solver());
559 5425 : KSP ksp = petsc_solver->ksp();
560 :
561 5425 : if (problem.hasLinearConvergenceObjects())
562 447 : 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 368162 : }
570 :
571 : void
572 76672 : processSingletonMooseWrappedOptions(FEProblemBase & fe_problem, const InputParameters & params)
573 : {
574 76672 : setSolveTypeFromParams(fe_problem, params);
575 76672 : setLineSearchFromParams(fe_problem, params);
576 76672 : setMFFDTypeFromParams(fe_problem, params);
577 76672 : }
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), "PETSc prefixes should be all lower-case")
585 :
586 : void
587 76672 : storePetscOptions(FEProblemBase & fe_problem,
588 : const std::string & prefix,
589 : const ParallelParamObject & param_object)
590 : {
591 76672 : const auto & params = param_object.parameters();
592 76672 : processSingletonMooseWrappedOptions(fe_problem, params);
593 :
594 : // The parameters contained in the Action
595 76672 : const auto & petsc_options = params.get<MultiMooseEnum>("petsc_options");
596 : const auto & petsc_pair_options =
597 76672 : params.get<MooseEnumItem, std::string>("petsc_options_iname", "petsc_options_value");
598 :
599 : // A reference to the PetscOptions object that contains the settings that will be used in the
600 : // solve
601 76672 : auto & po = fe_problem.getPetscOptions();
602 :
603 : // First process the single petsc options/flags
604 76672 : addPetscFlagsToPetscOptions(petsc_options, prefix, param_object, po);
605 :
606 : // Then process the option-value pairs
607 76672 : addPetscPairsToPetscOptions(
608 76672 : petsc_pair_options, fe_problem.mesh().dimension(), prefix, param_object, po);
609 76672 : }
610 :
611 : void
612 91098 : setSolveTypeFromParams(FEProblemBase & fe_problem, const InputParameters & params)
613 : {
614 : // Note: Options set in the Preconditioner block will override those set in the Executioner block
615 91098 : if (params.isParamValid("solve_type") && !params.isParamValid("_use_eigen_value"))
616 : {
617 : // Extract the solve type
618 37492 : const std::string & solve_type = params.get<MooseEnum>("solve_type");
619 75430 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
620 37938 : fe_problem.solverParams(i)._type = Moose::stringToEnum<Moose::SolveType>(solve_type);
621 37492 : }
622 91098 : }
623 :
624 : void
625 76672 : setLineSearchFromParams(FEProblemBase & fe_problem, const InputParameters & params)
626 : {
627 : // Note: Options set in the Preconditioner block will override those set in the Executioner block
628 76672 : if (params.isParamValid("line_search"))
629 : {
630 62186 : const auto & line_search = params.get<MooseEnum>("line_search");
631 123554 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
632 61368 : if (fe_problem.solverParams(i)._line_search == Moose::LS_INVALID || line_search != "default")
633 : {
634 : Moose::LineSearchType enum_line_search =
635 61008 : Moose::stringToEnum<Moose::LineSearchType>(line_search);
636 61008 : fe_problem.solverParams(i)._line_search = enum_line_search;
637 61008 : if (enum_line_search == LS_CONTACT || enum_line_search == LS_PROJECT)
638 : {
639 0 : NonlinearImplicitSystem * nl_system = dynamic_cast<NonlinearImplicitSystem *>(
640 0 : &fe_problem.getNonlinearSystemBase(i).system());
641 0 : if (!nl_system)
642 0 : mooseError("You've requested a line search but you must be solving an EigenProblem. "
643 : "These two things are not consistent.");
644 : PetscNonlinearSolver<Real> * petsc_nonlinear_solver =
645 0 : dynamic_cast<PetscNonlinearSolver<Real> *>(nl_system->nonlinear_solver.get());
646 0 : if (!petsc_nonlinear_solver)
647 0 : mooseError("Currently the MOOSE line searches all use Petsc, so you "
648 : "must use Petsc as your non-linear solver.");
649 : petsc_nonlinear_solver->linesearch_object =
650 0 : std::make_unique<ComputeLineSearchObjectWrapper>(fe_problem);
651 : }
652 : }
653 : }
654 76672 : }
655 :
656 : void
657 76672 : setMFFDTypeFromParams(FEProblemBase & fe_problem, const InputParameters & params)
658 : {
659 76672 : if (params.isParamValid("mffd_type"))
660 : {
661 76209 : const auto & mffd_type = params.get<MooseEnum>("mffd_type");
662 151868 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
663 75659 : fe_problem.solverParams(i)._mffd_type = Moose::stringToEnum<Moose::MffdType>(mffd_type);
664 : }
665 76672 : }
666 :
667 : template <typename T>
668 : void
669 59079 : checkUserProvidedPetscOption(const T & option, const ParallelParamObject & param_object)
670 : {
671 59079 : const auto & string_option = static_cast<const std::string &>(option);
672 59079 : if (string_option[0] != '-')
673 0 : param_object.mooseError("PETSc option '", string_option, "' does not begin with '-'");
674 59079 : }
675 :
676 : void
677 76672 : addPetscFlagsToPetscOptions(const MultiMooseEnum & petsc_flags,
678 : const std::string & prefix,
679 : const ParallelParamObject & param_object,
680 : PetscOptions & po)
681 : {
682 : checkPrefix(prefix);
683 :
684 : // Update the PETSc single flags
685 79386 : for (const auto & option : petsc_flags)
686 : {
687 2714 : checkUserProvidedPetscOption(option, param_object);
688 :
689 2714 : const std::string & string_option = option.name();
690 :
691 : /**
692 : * "-log_summary" cannot be used in the input file. This option needs to be set when PETSc is
693 : * initialized
694 : * which happens before the parser is even created. We'll throw an error if somebody attempts
695 : * to add this option later.
696 : */
697 2714 : if (option == "-log_summary" || option == "-log_view")
698 0 : mooseError("The PETSc option \"-log_summary\" or \"-log_view\" can only be used on the "
699 : "command line. Please "
700 : "remove it from the input file");
701 :
702 : // Update the stored items, but do not create duplicates
703 2714 : const std::string prefixed_option = prefix + string_option.substr(1);
704 2714 : if (!po.flags.isValueSet(prefixed_option))
705 : {
706 2714 : po.flags.setAdditionalValue(prefixed_option);
707 2714 : po.user_set_options.setAdditionalValue(prefixed_option);
708 : }
709 2714 : }
710 76672 : }
711 :
712 : void
713 62186 : setConvergedReasonFlags(FEProblemBase & fe_problem, const std::string & prefix)
714 : {
715 : checkPrefix(prefix);
716 62186 : libmesh_ignore(fe_problem); // avoid unused warnings for old PETSc
717 :
718 : #if !PETSC_VERSION_LESS_THAN(3, 14, 0)
719 : // the boolean in these pairs denote whether the user has specified any of the reason flags in the
720 : // input file
721 62186 : std::array<std::string, 2> reason_flags = {{"snes_converged_reason", "ksp_converged_reason"}};
722 :
723 62186 : auto & po = fe_problem.getPetscOptions();
724 :
725 186558 : for (const auto & reason_flag : reason_flags)
726 372259 : if (!po.flags.isValueSet(prefix + reason_flag) &&
727 123515 : (std::find_if(po.pairs.begin(),
728 : po.pairs.end(),
729 175105 : [&reason_flag, &prefix](auto & pair)
730 546507 : { return pair.first == (prefix + reason_flag); }) == po.pairs.end()))
731 123515 : po.pairs.emplace_back(prefix + reason_flag, "::failed");
732 : #endif
733 62186 : }
734 :
735 : void
736 76910 : addPetscPairsToPetscOptions(
737 : const std::vector<std::pair<MooseEnumItem, std::string>> & petsc_pair_options,
738 : const unsigned int mesh_dimension,
739 : const std::string & prefix,
740 : const ParallelParamObject & param_object,
741 : PetscOptions & po)
742 : {
743 : checkPrefix(prefix);
744 :
745 : // Setup the name value pairs
746 76910 : bool boomeramg_found = false;
747 76910 : bool strong_threshold_found = false;
748 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
749 76910 : bool superlu_dist_found = false;
750 76910 : bool fact_pattern_found = false;
751 76910 : bool tiny_pivot_found = false;
752 : #endif
753 76910 : std::string pc_description = "";
754 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
755 : // If users use HMG, we would like to set
756 76910 : bool hmg_found = false;
757 76910 : bool matptap_found = false;
758 76910 : bool hmg_strong_threshold_found = false;
759 : #endif
760 76910 : std::vector<std::pair<std::string, std::string>> new_options;
761 :
762 133275 : for (const auto & [option_name, option_value] : petsc_pair_options)
763 : {
764 56365 : checkUserProvidedPetscOption(option_name, param_object);
765 :
766 56365 : new_options.clear();
767 : const std::string prefixed_option_name =
768 56365 : prefix + static_cast<const std::string &>(option_name).substr(1);
769 :
770 : // Do not add duplicate settings
771 56365 : if (auto it =
772 56365 : MooseUtils::findPair(po.pairs, po.pairs.begin(), prefixed_option_name, MooseUtils::Any);
773 56365 : it == po.pairs.end())
774 : {
775 : #if !PETSC_VERSION_LESS_THAN(3, 9, 0)
776 56245 : if (option_name == "-pc_factor_mat_solver_package")
777 156 : new_options.emplace_back(prefix + "pc_factor_mat_solver_type", option_value);
778 : #else
779 : if (option_name == "-pc_factor_mat_solver_type")
780 : new_options.push_back(prefix + "pc_factor_mat_solver_package", option_value);
781 : #endif
782 :
783 : // Look for a pc description
784 84649 : if (option_name == "-pc_type" || option_name == "-sub_pc_type" ||
785 28404 : option_name == "-pc_hypre_type")
786 49568 : pc_description += option_value + ' ';
787 :
788 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
789 56245 : if (option_name == "-pc_type" && option_value == "hmg")
790 42 : hmg_found = true;
791 :
792 : // MPIAIJ for PETSc 3.12.0: -matptap_via
793 : // MAIJ for PETSc 3.12.0: -matmaijptap_via
794 : // MPIAIJ for PETSc 3.13 to 3.16: -matptap_via, -matproduct_ptap_via
795 : // MAIJ for PETSc 3.13 to 3.16: -matproduct_ptap_via
796 : // MPIAIJ for PETSc 3.17 and higher: -matptap_via, -mat_product_algorithm
797 : // MAIJ for PETSc 3.17 and higher: -mat_product_algorithm
798 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
799 56371 : if (hmg_found && (option_name == "-matptap_via" || option_name == "-matmaijptap_via" ||
800 126 : option_name == "-matproduct_ptap_via"))
801 0 : new_options.emplace_back(prefix + "mat_product_algorithm", option_value);
802 : #elif !PETSC_VERSION_LESS_THAN(3, 13, 0)
803 : if (hmg_found && (option_name == "-matptap_via" || option_name == "-matmaijptap_via"))
804 : new_options.emplace_back(prefix + "matproduct_ptap_via", option_value);
805 : #else
806 : if (hmg_found && (option_name == "-matproduct_ptap_via"))
807 : {
808 : new_options.emplace_back(prefix + "matptap_via", option_value);
809 : new_options.emplace_back(prefix + "matmaijptap_via", option_value);
810 : }
811 : #endif
812 :
813 112490 : if (option_name == "-matptap_via" || option_name == "-matmaijptap_via" ||
814 112490 : option_name == "-matproduct_ptap_via" || option_name == "-mat_product_algorithm")
815 0 : matptap_found = true;
816 :
817 : // For 3D problems, we need to set this 0.7
818 56245 : if (option_name == "-hmg_inner_pc_hypre_boomeramg_strong_threshold")
819 0 : hmg_strong_threshold_found = true;
820 : #endif
821 : // This special case is common enough that we'd like to handle it for the user.
822 56245 : if (option_name == "-pc_hypre_type" && option_value == "boomeramg")
823 21727 : boomeramg_found = true;
824 56245 : if (option_name == "-pc_hypre_boomeramg_strong_threshold")
825 0 : strong_threshold_found = true;
826 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
827 56245 : if ((option_name == "-pc_factor_mat_solver_package" ||
828 57206 : option_name == "-pc_factor_mat_solver_type") &&
829 961 : option_value == "superlu_dist")
830 311 : superlu_dist_found = true;
831 56245 : if (option_name == "-mat_superlu_dist_fact")
832 0 : fact_pattern_found = true;
833 56245 : if (option_name == "-mat_superlu_dist_replacetinypivot")
834 0 : tiny_pivot_found = true;
835 : #endif
836 :
837 56245 : if (!new_options.empty())
838 : {
839 156 : std::copy(new_options.begin(), new_options.end(), std::back_inserter(po.pairs));
840 312 : for (const auto & option : new_options)
841 156 : po.user_set_options.setAdditionalValue(option.first);
842 : }
843 : else
844 : {
845 56089 : po.pairs.push_back(std::make_pair(prefixed_option_name, option_value));
846 56089 : po.user_set_options.setAdditionalValue(prefixed_option_name);
847 : }
848 : }
849 : else
850 : {
851 : do
852 : {
853 120 : it->second = option_value;
854 120 : it = MooseUtils::findPair(po.pairs, std::next(it), prefixed_option_name, MooseUtils::Any);
855 120 : } while (it != po.pairs.end());
856 : }
857 56365 : }
858 :
859 : // When running a 3D mesh with boomeramg, it is almost always best to supply a strong threshold
860 : // value. We will provide that for the user here if they haven't supplied it themselves.
861 76910 : if (boomeramg_found && !strong_threshold_found && mesh_dimension == 3)
862 : {
863 984 : po.pairs.emplace_back(prefix + "pc_hypre_boomeramg_strong_threshold", "0.7");
864 984 : pc_description += "strong_threshold: 0.7 (auto)";
865 : }
866 :
867 : #if !PETSC_VERSION_LESS_THAN(3, 12, 0)
868 76910 : if (hmg_found && !hmg_strong_threshold_found && mesh_dimension == 3)
869 : {
870 28 : po.pairs.emplace_back(prefix + "hmg_inner_pc_hypre_boomeramg_strong_threshold", "0.7");
871 28 : pc_description += "strong_threshold: 0.7 (auto)";
872 : }
873 :
874 : // Default PETSc PtAP takes too much memory, and it is not quite useful
875 : // Let us switch to use new algorithm
876 76910 : if (hmg_found && !matptap_found)
877 : {
878 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
879 42 : po.pairs.emplace_back(prefix + "mat_product_algorithm", "allatonce");
880 : #elif !PETSC_VERSION_LESS_THAN(3, 13, 0)
881 : po.pairs.emplace_back(prefix + "matproduct_ptap_via", "allatonce");
882 : #else
883 : po.pairs.emplace_back(prefix + "matptap_via", "allatonce");
884 : po.pairs.emplace_back(prefix + "matmaijptap_via", "allatonce");
885 : #endif
886 : }
887 : #endif
888 :
889 : #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
890 : // In PETSc-3.7.{0--4}, there is a bug when using superlu_dist, and we have to use
891 : // SamePattern_SameRowPerm, otherwise we use whatever we have in PETSc
892 76910 : if (superlu_dist_found && !fact_pattern_found)
893 : {
894 311 : po.pairs.emplace_back(prefix + "mat_superlu_dist_fact",
895 : #if PETSC_VERSION_LESS_THAN(3, 7, 5)
896 : "SamePattern_SameRowPerm");
897 : pc_description += "mat_superlu_dist_fact: SamePattern_SameRowPerm ";
898 : #else
899 : "SamePattern");
900 311 : pc_description += "mat_superlu_dist_fact: SamePattern ";
901 : #endif
902 : }
903 :
904 : // restore this superlu option
905 76910 : if (superlu_dist_found && !tiny_pivot_found)
906 : {
907 311 : po.pairs.emplace_back(prefix + "mat_superlu_dist_replacetinypivot", "1");
908 311 : pc_description += " mat_superlu_dist_replacetinypivot: true ";
909 : }
910 : #endif
911 : // Set Preconditioner description
912 76910 : if (!pc_description.empty() && prefix.size() > 1)
913 706 : po.pc_description += "[" + prefix.substr(1, prefix.size() - 2) + "]: ";
914 76910 : po.pc_description += pc_description;
915 76910 : }
916 :
917 : std::set<std::string>
918 295827 : getPetscValidLineSearches()
919 : {
920 2662443 : return {"default", "shell", "none", "basic", "l2", "bt", "cp"};
921 591654 : }
922 :
923 : InputParameters
924 410305 : getPetscValidParams()
925 : {
926 410305 : InputParameters params = emptyInputParameters();
927 :
928 410305 : MooseEnum solve_type("PJFNK JFNK NEWTON FD LINEAR");
929 410305 : params.addParam<MooseEnum>("solve_type",
930 : solve_type,
931 : "PJFNK: Preconditioned Jacobian-Free Newton Krylov "
932 : "JFNK: Jacobian-Free Newton Krylov "
933 : "NEWTON: Full Newton Solve "
934 : "FD: Use finite differences to compute Jacobian "
935 : "LINEAR: Solving a linear problem");
936 :
937 410305 : MooseEnum mffd_type("wp ds", "wp");
938 410305 : params.addParam<MooseEnum>("mffd_type",
939 : mffd_type,
940 : "Specifies the finite differencing type for "
941 : "Jacobian-free solve types. Note that the "
942 : "default is wp (for Walker and Pernice).");
943 :
944 1230915 : params.addParam<MultiMooseEnum>(
945 820610 : "petsc_options", getCommonPetscFlags(), "Singleton PETSc options");
946 1230915 : params.addParam<MultiMooseEnum>(
947 820610 : "petsc_options_iname", getCommonPetscKeys(), "Names of PETSc name/value pairs");
948 410305 : params.addParam<std::vector<std::string>>(
949 : "petsc_options_value",
950 : "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\"");
951 410305 : params.addParamNamesToGroup("solve_type petsc_options petsc_options_iname petsc_options_value "
952 : "mffd_type",
953 : "PETSc");
954 :
955 820610 : return params;
956 410305 : }
957 :
958 : MultiMooseEnum
959 426658 : getCommonSNESFlags()
960 : {
961 : return MultiMooseEnum(
962 : "-ksp_monitor_snes_lg -snes_ksp_ew -ksp_snes_ew -snes_converged_reason "
963 : "-snes_ksp -snes_linesearch_monitor -snes_mf -snes_mf_operator -snes_monitor "
964 : "-snes_test_display -snes_view -snes_monitor_cancel",
965 : "",
966 426658 : true);
967 : }
968 :
969 : MultiMooseEnum
970 425994 : getCommonKSPFlags()
971 : {
972 : return MultiMooseEnum(
973 425994 : "-ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_monitor", "", true);
974 : }
975 :
976 : MultiMooseEnum
977 425994 : getCommonPetscFlags()
978 : {
979 425994 : auto options = MultiMooseEnum("-dm_moose_print_embedding -dm_view", "", true);
980 425994 : options.addValidName(getCommonKSPFlags());
981 425994 : options.addValidName(getCommonSNESFlags());
982 425994 : return options;
983 0 : }
984 :
985 : MultiMooseEnum
986 426658 : getCommonSNESKeys()
987 : {
988 : return MultiMooseEnum("-snes_atol -snes_linesearch_type -snes_ls -snes_max_it -snes_rtol "
989 : "-snes_divergence_tolerance -snes_type",
990 : "",
991 426658 : true);
992 : }
993 :
994 : MultiMooseEnum
995 425994 : getCommonKSPKeys()
996 : {
997 : return MultiMooseEnum("-ksp_atol -ksp_gmres_restart -ksp_max_it -ksp_pc_side -ksp_rtol "
998 : "-ksp_type -sub_ksp_type",
999 : "",
1000 425994 : true);
1001 : }
1002 : MultiMooseEnum
1003 425994 : getCommonPetscKeys()
1004 : {
1005 : auto options = MultiMooseEnum("-mat_fd_coloring_err -mat_fd_type -mat_mffd_type "
1006 : "-pc_asm_overlap -pc_factor_levels "
1007 : "-pc_factor_mat_ordering_type -pc_hypre_boomeramg_grid_sweeps_all "
1008 : "-pc_hypre_boomeramg_max_iter "
1009 : "-pc_hypre_boomeramg_strong_threshold -pc_hypre_type -pc_type "
1010 : "-sub_pc_type",
1011 : "",
1012 425994 : true);
1013 425994 : options.addValidName(getCommonKSPKeys());
1014 425994 : options.addValidName(getCommonSNESKeys());
1015 425994 : return options;
1016 0 : }
1017 :
1018 : bool
1019 604 : isSNESVI(FEProblemBase & fe_problem)
1020 : {
1021 604 : const PetscOptions & petsc = fe_problem.getPetscOptions();
1022 :
1023 : int argc;
1024 : char ** args;
1025 604 : LibmeshPetscCallA(fe_problem.comm().get(), PetscGetArgs(&argc, &args));
1026 :
1027 604 : std::vector<std::string> cml_arg;
1028 6724 : for (int i = 0; i < argc; i++)
1029 6120 : cml_arg.push_back(args[i]);
1030 :
1031 604 : if (MooseUtils::findPair(petsc.pairs, petsc.pairs.begin(), MooseUtils::Any, "vinewtonssls") ==
1032 1208 : petsc.pairs.end() &&
1033 604 : MooseUtils::findPair(petsc.pairs, petsc.pairs.begin(), MooseUtils::Any, "vinewtonrsls") ==
1034 648 : petsc.pairs.end() &&
1035 1252 : std::find(cml_arg.begin(), cml_arg.end(), "vinewtonssls") == cml_arg.end() &&
1036 648 : std::find(cml_arg.begin(), cml_arg.end(), "vinewtonrsls") == cml_arg.end())
1037 44 : return false;
1038 :
1039 560 : return true;
1040 604 : }
1041 :
1042 : void
1043 233124 : setSinglePetscOption(const std::string & name,
1044 : const std::string & value /*=""*/,
1045 : FEProblemBase * const problem /*=nullptr*/)
1046 : {
1047 233124 : static const TIMPI::Communicator comm_world(PETSC_COMM_WORLD);
1048 233124 : const TIMPI::Communicator & comm = problem ? problem->comm() : comm_world;
1049 233124 : LibmeshPetscCallA(comm.get(),
1050 : PetscOptionsSetValue(LIBMESH_PETSC_NULLPTR,
1051 : name.c_str(),
1052 : value == "" ? LIBMESH_PETSC_NULLPTR : value.c_str()));
1053 233124 : const auto lower_case_name = MooseUtils::toLower(name);
1054 20 : auto check_problem = [problem, &lower_case_name]()
1055 : {
1056 20 : if (!problem)
1057 0 : mooseError(
1058 : "Setting the option '",
1059 : lower_case_name,
1060 : "' requires passing a 'problem' parameter. Contact a developer of your application "
1061 : "to have them update their code. If in doubt, reach out to the MOOSE team on Github "
1062 : "discussions");
1063 233144 : };
1064 :
1065 : // Select vector type from user-passed PETSc options
1066 233124 : if (lower_case_name.find("-vec_type") != std::string::npos)
1067 : {
1068 10 : check_problem();
1069 20 : for (auto & solver_system : problem->_solver_systems)
1070 : {
1071 10 : auto & lm_sys = solver_system->system();
1072 30 : for (auto & [vec_name, vec] : as_range(lm_sys.vectors_begin(), lm_sys.vectors_end()))
1073 : {
1074 20 : libmesh_ignore(vec_name);
1075 20 : auto * const petsc_vec = cast_ptr<PetscVector<Number> *>(vec.get());
1076 20 : LibmeshPetscCallA(comm.get(), VecSetFromOptions(petsc_vec->vec()));
1077 : }
1078 : // The solution vectors aren't included in the system vectors storage
1079 10 : auto * petsc_vec = cast_ptr<PetscVector<Number> *>(lm_sys.solution.get());
1080 10 : LibmeshPetscCallA(comm.get(), VecSetFromOptions(petsc_vec->vec()));
1081 10 : petsc_vec = cast_ptr<PetscVector<Number> *>(lm_sys.current_local_solution.get());
1082 10 : LibmeshPetscCallA(comm.get(), VecSetFromOptions(petsc_vec->vec()));
1083 : }
1084 : }
1085 : // Select matrix type from user-passed PETSc options
1086 233114 : else if (lower_case_name.find("mat_type") != std::string::npos)
1087 : {
1088 10 : check_problem();
1089 :
1090 10 : bool found_matching_prefix = false;
1091 10 : for (const auto i : index_range(problem->_solver_systems))
1092 : {
1093 10 : const auto & solver_sys_name = problem->_solver_sys_names[i];
1094 10 : if (lower_case_name.find("-" + MooseUtils::toLower(solver_sys_name) + "_mat_type") ==
1095 : std::string::npos)
1096 0 : continue;
1097 :
1098 10 : if (problem->solverParams(i)._type == Moose::ST_JFNK)
1099 0 : mooseError(
1100 : "Setting option '", lower_case_name, "' is incompatible with a JFNK 'solve_type'");
1101 :
1102 10 : auto & lm_sys = problem->_solver_systems[i]->system();
1103 20 : for (auto & [mat_name, mat] : as_range(lm_sys.matrices_begin(), lm_sys.matrices_end()))
1104 : {
1105 10 : libmesh_ignore(mat_name);
1106 10 : if (auto * const petsc_mat = dynamic_cast<PetscMatrixBase<Number> *>(mat.get()); petsc_mat)
1107 : {
1108 : #ifdef DEBUG
1109 : const char * mat_prefix;
1110 : LibmeshPetscCallA(comm.get(), MatGetOptionsPrefix(petsc_mat->mat(), &mat_prefix));
1111 : mooseAssert(strcmp(mat_prefix, (solver_sys_name + "_").c_str()) == 0,
1112 : "We should have prefixed these matrices previously");
1113 : #endif
1114 10 : LibmeshPetscCallA(comm.get(), MatSetFromOptions(petsc_mat->mat()));
1115 : }
1116 : }
1117 10 : found_matching_prefix = true;
1118 10 : break;
1119 : }
1120 :
1121 10 : if (!found_matching_prefix)
1122 0 : mooseError("We did not find a matching solver system name for the petsc option '",
1123 : lower_case_name,
1124 : "'");
1125 : }
1126 233124 : }
1127 :
1128 : void
1129 71630 : setSinglePetscOptionIfAppropriate(const MultiMooseEnum & dont_add_these_options,
1130 : const std::string & name,
1131 : const std::string & value /*=""*/,
1132 : FEProblemBase * const problem /*=nullptr*/)
1133 : {
1134 71630 : if (!dont_add_these_options.contains(name))
1135 70193 : setSinglePetscOption(name, value, problem);
1136 71630 : }
1137 :
1138 : void
1139 0 : colorAdjacencyMatrix(PetscScalar * adjacency_matrix,
1140 : unsigned int size,
1141 : unsigned int colors,
1142 : std::vector<unsigned int> & vertex_colors,
1143 : const char * coloring_algorithm)
1144 : {
1145 : // Mat A will be a dense matrix from the incoming data structure
1146 : Mat A;
1147 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatCreate(PETSC_COMM_SELF, &A));
1148 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSetSizes(A, size, size, size, size));
1149 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSetType(A, MATSEQDENSE));
1150 : // PETSc requires a non-const data array to populate the matrix
1151 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatSeqDenseSetPreallocation(A, adjacency_matrix));
1152 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1153 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1154 :
1155 : // Convert A to a sparse matrix
1156 : #if PETSC_VERSION_LESS_THAN(3, 7, 0)
1157 : LibmeshPetscCallA(PETSC_COMM_SELF, MatConvert(A, MATAIJ, MAT_REUSE_MATRIX, &A));
1158 : #else
1159 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A));
1160 : #endif
1161 :
1162 : ISColoring iscoloring;
1163 : MatColoring mc;
1164 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringCreate(A, &mc));
1165 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetType(mc, coloring_algorithm));
1166 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetMaxColors(mc, static_cast<PetscInt>(colors)));
1167 :
1168 : // Petsc normally colors by distance two (neighbors of neighbors), we just want one
1169 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetDistance(mc, 1));
1170 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetFromOptions(mc));
1171 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringApply(mc, &iscoloring));
1172 :
1173 : PetscInt nn;
1174 : IS * is;
1175 : #if PETSC_RELEASE_LESS_THAN(3, 12, 0)
1176 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringGetIS(iscoloring, &nn, &is));
1177 : #else
1178 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &nn, &is));
1179 : #endif
1180 :
1181 0 : if (nn > static_cast<PetscInt>(colors))
1182 0 : throw std::runtime_error("Not able to color with designated number of colors");
1183 :
1184 0 : for (int i = 0; i < nn; i++)
1185 : {
1186 : PetscInt isize;
1187 : const PetscInt * indices;
1188 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISGetLocalSize(is[i], &isize));
1189 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISGetIndices(is[i], &indices));
1190 0 : for (int j = 0; j < isize; j++)
1191 : {
1192 : mooseAssert(indices[j] < static_cast<PetscInt>(vertex_colors.size()), "Index out of bounds");
1193 0 : vertex_colors[indices[j]] = i;
1194 : }
1195 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISRestoreIndices(is[i], &indices));
1196 : }
1197 :
1198 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatDestroy(&A));
1199 0 : LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringDestroy(&mc));
1200 0 : LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringDestroy(&iscoloring));
1201 0 : }
1202 :
1203 : void
1204 3748 : dontAddPetscFlag(const std::string & flag, PetscOptions & petsc_options)
1205 : {
1206 3748 : if (!petsc_options.dont_add_these_options.contains(flag))
1207 3748 : petsc_options.dont_add_these_options.setAdditionalValue(flag);
1208 3748 : }
1209 :
1210 : void
1211 656 : dontAddNonlinearConvergedReason(FEProblemBase & fe_problem)
1212 : {
1213 656 : dontAddPetscFlag("-snes_converged_reason", fe_problem.getPetscOptions());
1214 656 : }
1215 :
1216 : void
1217 666 : dontAddLinearConvergedReason(FEProblemBase & fe_problem)
1218 : {
1219 666 : dontAddPetscFlag("-ksp_converged_reason", fe_problem.getPetscOptions());
1220 666 : }
1221 :
1222 : void
1223 0 : dontAddCommonKSPOptions(FEProblemBase & fe_problem)
1224 : {
1225 0 : auto & petsc_options = fe_problem.getPetscOptions();
1226 0 : for (const auto & flag : getCommonKSPFlags().getNames())
1227 0 : dontAddPetscFlag(flag, petsc_options);
1228 0 : for (const auto & key : getCommonKSPKeys().getNames())
1229 0 : dontAddPetscFlag(key, petsc_options);
1230 0 : }
1231 :
1232 : void
1233 664 : dontAddCommonSNESOptions(FEProblemBase & fe_problem)
1234 : {
1235 664 : auto & petsc_options = fe_problem.getPetscOptions();
1236 8632 : for (const auto & flag : getCommonSNESFlags().getNames())
1237 7968 : if (!petsc_options.dont_add_these_options.contains(flag))
1238 4648 : petsc_options.dont_add_these_options.setAdditionalValue(flag);
1239 5312 : for (const auto & key : getCommonSNESKeys().getNames())
1240 4648 : if (!petsc_options.dont_add_these_options.contains(key))
1241 2988 : petsc_options.dont_add_these_options.setAdditionalValue(key);
1242 664 : }
1243 :
1244 : std::unique_ptr<PetscMatrix<Number>>
1245 48 : createMatrixFromFile(const libMesh::Parallel::Communicator & comm,
1246 : Mat & mat,
1247 : const std::string & binary_mat_file,
1248 : const unsigned int mat_number_to_load)
1249 : {
1250 48 : LibmeshPetscCallA(comm.get(), MatCreate(comm.get(), &mat));
1251 : PetscViewer matviewer;
1252 48 : LibmeshPetscCallA(
1253 : comm.get(),
1254 : PetscViewerBinaryOpen(comm.get(), binary_mat_file.c_str(), FILE_MODE_READ, &matviewer));
1255 96 : for (unsigned int i = 0; i < mat_number_to_load; ++i)
1256 48 : LibmeshPetscCallA(comm.get(), MatLoad(mat, matviewer));
1257 48 : LibmeshPetscCallA(comm.get(), PetscViewerDestroy(&matviewer));
1258 :
1259 96 : return std::make_unique<PetscMatrix<Number>>(mat, comm);
1260 : }
1261 :
1262 : } // Namespace PetscSupport
1263 : } // Namespace MOOSE
|
