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 : #pragma once
11 :
12 : #include "SolverSystem.h"
13 : #include "ConstraintWarehouse.h"
14 : #include "MooseObjectWarehouse.h"
15 : #include "MooseObjectTagWarehouse.h"
16 : #include "PerfGraphInterface.h"
17 : #include "ComputeMortarFunctor.h"
18 : #include "MooseHashing.h"
19 :
20 : #include "libmesh/transient_system.h"
21 : #include "libmesh/nonlinear_implicit_system.h"
22 : #include "libmesh/linear_solver.h"
23 :
24 : // Forward declarations
25 : class FEProblemBase;
26 : class MoosePreconditioner;
27 : class JacobianBlock;
28 : class TimeIntegrator;
29 : class Predictor;
30 : class ElementDamper;
31 : class NodalDamper;
32 : class GeneralDamper;
33 : class GeometricSearchData;
34 : class IntegratedBCBase;
35 : class NodalBCBase;
36 : class DirichletBCBase;
37 : class ADDirichletBCBase;
38 : class DGKernelBase;
39 : class InterfaceKernelBase;
40 : class ScalarKernelBase;
41 : class DiracKernelBase;
42 : class NodalKernelBase;
43 : class Split;
44 : class KernelBase;
45 : class HDGKernel;
46 : class BoundaryCondition;
47 : class ResidualObject;
48 : class PenetrationInfo;
49 : class FieldSplitPreconditionerBase;
50 :
51 : // libMesh forward declarations
52 : namespace libMesh
53 : {
54 : template <typename T>
55 : class NumericVector;
56 : template <typename T>
57 : class SparseMatrix;
58 : template <typename T>
59 : class DiagonalMatrix;
60 : class DofMapBase;
61 : } // namespace libMesh
62 :
63 : /**
64 : * Nonlinear system to be solved
65 : *
66 : * It is a part of FEProblemBase ;-)
67 : */
68 : class NonlinearSystemBase : public SolverSystem, public PerfGraphInterface
69 : {
70 : public:
71 : NonlinearSystemBase(FEProblemBase & problem, libMesh::System & sys, const std::string & name);
72 : virtual ~NonlinearSystemBase();
73 :
74 : virtual void preInit() override;
75 : /// Update the mortar functors if the mesh has changed
76 : void reinitMortarFunctors();
77 :
78 : bool computedScalingJacobian() const { return _computed_scaling; }
79 :
80 : /**
81 : * Turn off the Jacobian (must be called before equation system initialization)
82 : */
83 : virtual void turnOffJacobian();
84 :
85 : virtual void solve() override = 0;
86 :
87 : virtual libMesh::NonlinearSolver<Number> * nonlinearSolver() = 0;
88 :
89 : virtual SNES getSNES() = 0;
90 :
91 : virtual unsigned int getCurrentNonlinearIterationNumber() = 0;
92 :
93 : /**
94 : * Returns true if this system is currently computing the pre-SMO residual for a solve.
95 : * @return Whether or not we are currently computing the pre-SMO residual.
96 : */
97 845824 : bool computingPreSMOResidual() { return _computing_pre_smo_residual; }
98 :
99 : // Setup Functions ////
100 : virtual void initialSetup() override;
101 : virtual void timestepSetup() override;
102 : virtual void customSetup(const ExecFlagType & exec_type) override;
103 : virtual void residualSetup() override;
104 : virtual void jacobianSetup() override;
105 :
106 : virtual void setupFiniteDifferencedPreconditioner() = 0;
107 :
108 : bool haveFiniteDifferencedPreconditioner() const
109 : {
110 : return _use_finite_differenced_preconditioner;
111 : }
112 : bool haveFieldSplitPreconditioner() const { return _fsp; }
113 :
114 : /**
115 : * Adds a kernel
116 : * @param kernel_name The type of the kernel
117 : * @param name The name of the kernel
118 : * @param parameters Kernel parameters
119 : */
120 : virtual void addKernel(const std::string & kernel_name,
121 : const std::string & name,
122 : InputParameters & parameters);
123 :
124 : /**
125 : * Adds a hybridized discontinuous Galerkin (HDG) kernel
126 : * @param kernel_name The type of the hybridized kernel
127 : * @param name The name of the hybridized kernel
128 : * @param parameters HDG kernel parameters
129 : */
130 : virtual void addHDGKernel(const std::string & kernel_name,
131 : const std::string & name,
132 : InputParameters & parameters);
133 :
134 : /**
135 : * Adds a NodalKernel
136 : * @param kernel_name The type of the nodal kernel
137 : * @param name The name of the kernel
138 : * @param parameters Kernel parameters
139 : */
140 : virtual void addNodalKernel(const std::string & kernel_name,
141 : const std::string & name,
142 : InputParameters & parameters);
143 :
144 : /**
145 : * Adds a scalar kernel
146 : * @param kernel_name The type of the kernel
147 : * @param name The name of the kernel
148 : * @param parameters Kernel parameters
149 : */
150 : void addScalarKernel(const std::string & kernel_name,
151 : const std::string & name,
152 : InputParameters & parameters);
153 :
154 : /**
155 : * Adds a boundary condition
156 : * @param bc_name The type of the boundary condition
157 : * @param name The name of the boundary condition
158 : * @param parameters Boundary condition parameters
159 : */
160 : void addBoundaryCondition(const std::string & bc_name,
161 : const std::string & name,
162 : InputParameters & parameters);
163 :
164 : #ifdef MOOSE_KOKKOS_ENABLED
165 : /**
166 : * Adds a Kokkos kernel
167 : * @param kernel_name The type of the kernel
168 : * @param name The name of the kernel
169 : * @param parameters Kernel parameters
170 : */
171 : virtual void addKokkosKernel(const std::string & kernel_name,
172 : const std::string & name,
173 : InputParameters & parameters);
174 :
175 : /**
176 : * Adds a Kokkos nodal kernel
177 : * @param kernel_name The type of the nodal kernel
178 : * @param name The name of the kernel
179 : * @param parameters Kernel parameters
180 : */
181 : virtual void addKokkosNodalKernel(const std::string & kernel_name,
182 : const std::string & name,
183 : InputParameters & parameters);
184 :
185 : /**
186 : * Adds a Kokkos boundary condition
187 : * @param bc_name The type of the boundary condition
188 : * @param name The name of the boundary condition
189 : * @param parameters Boundary condition parameters
190 : */
191 : void addKokkosBoundaryCondition(const std::string & bc_name,
192 : const std::string & name,
193 : InputParameters & parameters);
194 : #endif
195 :
196 : /**
197 : * Adds a Constraint
198 : * @param c_name The type of the constraint
199 : * @param name The name of the constraint
200 : * @param parameters Constraint parameters
201 : */
202 : void
203 : addConstraint(const std::string & c_name, const std::string & name, InputParameters & parameters);
204 :
205 : /**
206 : * Adds a Dirac kernel
207 : * @param kernel_name The type of the dirac kernel
208 : * @param name The name of the Dirac kernel
209 : * @param parameters Dirac kernel parameters
210 : */
211 : void addDiracKernel(const std::string & kernel_name,
212 : const std::string & name,
213 : InputParameters & parameters);
214 :
215 : /**
216 : * Adds a DG kernel
217 : * @param dg_kernel_name The type of the DG kernel
218 : * @param name The name of the DG kernel
219 : * @param parameters DG kernel parameters
220 : */
221 : void
222 : addDGKernel(std::string dg_kernel_name, const std::string & name, InputParameters & parameters);
223 :
224 : /**
225 : * Adds an interface kernel
226 : * @param interface_kernel_name The type of the interface kernel
227 : * @param name The name of the interface kernel
228 : * @param parameters interface kernel parameters
229 : */
230 : void addInterfaceKernel(std::string interface_kernel_name,
231 : const std::string & name,
232 : InputParameters & parameters);
233 :
234 : /**
235 : * Adds a damper
236 : * @param damper_name The type of the damper
237 : * @param name The name of the damper
238 : * @param parameters Damper parameters
239 : */
240 : void addDamper(const std::string & damper_name,
241 : const std::string & name,
242 : InputParameters & parameters);
243 :
244 : /**
245 : * Adds a split
246 : * @param split_name The type of the split
247 : * @param name The name of the split
248 : * @param parameters Split parameters
249 : */
250 : void
251 : addSplit(const std::string & split_name, const std::string & name, InputParameters & parameters);
252 :
253 : /**
254 : * Retrieves a split by name
255 : * @param name The name of the split
256 : */
257 : std::shared_ptr<Split> getSplit(const std::string & name);
258 :
259 : /**
260 : * Retrieves all splits
261 : */
262 55469 : MooseObjectWarehouseBase<Split> & getSplits() { return _splits; }
263 :
264 : /**
265 : * We offer the option to check convergence against the pre-SMO residual. This method handles the
266 : * logic as to whether we should perform such residual evaluation.
267 : *
268 : * @return A boolean indicating whether we should evaluate the pre-SMO residual.
269 : */
270 : bool shouldEvaluatePreSMOResidual() const;
271 :
272 : /**
273 : * Set whether to evaluate the pre-SMO residual and use it in the subsequent relative convergence
274 : * checks.
275 : *
276 : * If set to true, an _additional_ residual evaluation is performed before any
277 : * solution-modifying object is executed, and before the initial (0-th nonlinear iteration)
278 : * residual evaluation. Such residual is referred to as the pre-SMO residual. If the pre-SMO
279 : * residual is evaluated, it is used in the subsequent relative convergence checks.
280 : *
281 : * If set to false, no residual evaluation takes place before the initial residual evaluation, and
282 : * the initial residual is used in the subsequent relative convergence checks. This mode is
283 : * recommended for performance-critical code as it avoids the additional pre-SMO residual
284 : * evaluation.
285 : */
286 62595 : void setPreSMOResidual(bool use) { _use_pre_smo_residual = use; }
287 :
288 : /// Whether we are using pre-SMO residual in relative convergence checks
289 449663 : const bool & usePreSMOResidual() const { return _use_pre_smo_residual; }
290 :
291 : /// The reference residual used in relative convergence check.
292 : Real referenceResidual() const;
293 :
294 : /// The pre-SMO residual
295 : Real preSMOResidual() const;
296 :
297 : /// The initial residual
298 : Real initialResidual() const;
299 :
300 : /// Record the initial residual (for later relative convergence check)
301 : void setInitialResidual(Real r);
302 :
303 : void zeroVectorForResidual(const std::string & vector_name);
304 :
305 : void setInitialSolution();
306 :
307 : #ifdef MOOSE_KOKKOS_ENABLED
308 : void setKokkosInitialSolution();
309 : #endif
310 :
311 : /**
312 : * Sets the value of constrained variables in the solution vector.
313 : */
314 : void setConstraintSecondaryValues(NumericVector<Number> & solution, bool displaced);
315 :
316 : /**
317 : * Add residual contributions from Constraints
318 : *
319 : * @param residual - reference to the residual vector where constraint contributions will be
320 : * computed
321 : * @param displaced Controls whether to do the displaced Constraints or non-displaced
322 : */
323 : void constraintResiduals(NumericVector<Number> & residual, bool displaced);
324 :
325 : /**
326 : * Computes residual for a given tag
327 : * @param residual Residual is formed in here
328 : * @param the tag of kernels for which the residual is to be computed.
329 : */
330 : void computeResidualTag(NumericVector<Number> & residual, TagID tag_id);
331 :
332 : /**
333 : * Form multiple tag-associated residual vectors for all the given tags
334 : */
335 : void computeResidualTags(const std::set<TagID> & tags);
336 :
337 : /**
338 : * Form possibly multiple tag-associated vectors and matrices
339 : */
340 : void computeResidualAndJacobianTags(const std::set<TagID> & vector_tags,
341 : const std::set<TagID> & matrix_tags);
342 :
343 : /**
344 : * Compute residual and Jacobian from contributions not related to constraints, such as nodal
345 : * boundary conditions
346 : */
347 : void computeResidualAndJacobianInternal(const std::set<TagID> & vector_tags,
348 : const std::set<TagID> & matrix_tags);
349 :
350 : /**
351 : * Form a residual vector for a given tag
352 : */
353 : void computeResidual(NumericVector<Number> & residual, TagID tag_id);
354 :
355 : /**
356 : * Adds entries to the Jacobian in the correct positions for couplings coming from dofs being
357 : * coupled that
358 : * are related geometrically (i.e. near each other across a gap).
359 : */
360 : void addImplicitGeometricCouplingEntries(GeometricSearchData & geom_search_data);
361 :
362 : /**
363 : * Add jacobian contributions from Constraints
364 : *
365 : * @param jacobian reference to a read-only view of the Jacobian matrix
366 : * @param displaced Controls whether to do the displaced Constraints or non-displaced
367 : */
368 : void constraintJacobians(const SparseMatrix<Number> & jacobian_to_view, bool displaced);
369 :
370 : /**
371 : * Computes multiple (tag associated) Jacobian matricese
372 : */
373 : void computeJacobianTags(const std::set<TagID> & tags);
374 :
375 : /**
376 : * Method used to obtain scaling factors for variables
377 : * @returns whether this method ran without exceptions
378 : */
379 : bool computeScaling();
380 :
381 : /**
382 : * Associate jacobian to systemMatrixTag, and then form a matrix for all the tags
383 : */
384 : void computeJacobian(libMesh::SparseMatrix<Number> & jacobian, const std::set<TagID> & tags);
385 :
386 : /**
387 : * Take all tags in the system, and form a matrix for all tags in the system
388 : */
389 : void computeJacobian(libMesh::SparseMatrix<Number> & jacobian);
390 :
391 : /**
392 : * Computes several Jacobian blocks simultaneously, summing their contributions into smaller
393 : * preconditioning matrices.
394 : *
395 : * Used by Physics-based preconditioning
396 : *
397 : * @param blocks The blocks to fill in (JacobianBlock is defined in ComputeJacobianBlocksThread)
398 : */
399 : void computeJacobianBlocks(std::vector<JacobianBlock *> & blocks);
400 :
401 : void computeJacobianBlocks(std::vector<JacobianBlock *> & blocks, const std::set<TagID> & tags);
402 :
403 : /**
404 : * Compute damping
405 : * @param solution The trail solution vector
406 : * @param update The incremental update to the solution vector
407 : * @return returns The damping factor
408 : */
409 : Real computeDamping(const NumericVector<Number> & solution, const NumericVector<Number> & update);
410 :
411 : /**
412 : * Called at the beginning of the time step
413 : */
414 : void onTimestepBegin();
415 :
416 : using SystemBase::subdomainSetup;
417 : /**
418 : * Called from assembling when we hit a new subdomain
419 : * @param subdomain ID of the new subdomain
420 : * @param tid Thread ID
421 : */
422 : virtual void subdomainSetup(SubdomainID subdomain, THREAD_ID tid);
423 :
424 : /**
425 : * Called from explicit time stepping to overwrite boundary positions (explicit dynamics). This
426 : * will close/assemble the passed-in \p soln after overwrite
427 : */
428 : void overwriteNodeFace(NumericVector<Number> & soln);
429 :
430 : /**
431 : * Update active objects of Warehouses owned by NonlinearSystemBase
432 : */
433 : void updateActive(THREAD_ID tid);
434 :
435 : /**
436 : * Set transient term used by residual and Jacobian evaluation.
437 : * @param udot transient term
438 : * @note If the calling sequence for residual evaluation was changed, this could become an
439 : * explicit argument.
440 : */
441 : virtual void setSolutionUDot(const NumericVector<Number> & udot);
442 :
443 : /**
444 : * Set transient term used by residual and Jacobian evaluation.
445 : * @param udotdot transient term
446 : * @note If the calling sequence for residual evaluation was changed, this could become an
447 : * explicit argument.
448 : */
449 : virtual void setSolutionUDotDot(const NumericVector<Number> & udotdot);
450 :
451 : /**
452 : * Return a numeric vector that is associated with the time tag.
453 : */
454 : NumericVector<Number> & getResidualTimeVector();
455 :
456 : /**
457 : * Return a numeric vector that is associated with the nontime tag.
458 : */
459 : NumericVector<Number> & getResidualNonTimeVector();
460 :
461 : /**
462 : * Return a residual vector that is associated with the residual tag.
463 : */
464 : NumericVector<Number> & residualVector(TagID tag);
465 :
466 : virtual NumericVector<Number> & residualCopy() override;
467 : virtual NumericVector<Number> & residualGhosted() override;
468 :
469 : virtual NumericVector<Number> & RHS() = 0;
470 :
471 : virtual void augmentSparsity(libMesh::SparsityPattern::Graph & sparsity,
472 : std::vector<dof_id_type> & n_nz,
473 : std::vector<dof_id_type> & n_oz) override;
474 :
475 : /**
476 : * Sets a preconditioner
477 : * @param pc The preconditioner to be set
478 : */
479 : void setPreconditioner(std::shared_ptr<MoosePreconditioner> pc);
480 : MoosePreconditioner const * getPreconditioner() const;
481 :
482 : /**
483 : * If called with true this system will use a finite differenced form of
484 : * the Jacobian as the preconditioner
485 : */
486 120 : void useFiniteDifferencedPreconditioner(bool use = true)
487 : {
488 120 : _use_finite_differenced_preconditioner = use;
489 120 : }
490 :
491 : /**
492 : * If called with a non-null object true this system will use a field split preconditioner matrix.
493 : */
494 139 : void useFieldSplitPreconditioner(FieldSplitPreconditionerBase * fsp) { _fsp = fsp; }
495 :
496 : /**
497 : * @returns A field split preconditioner. This will error if there is no field split
498 : * preconditioner
499 : */
500 : FieldSplitPreconditionerBase & getFieldSplitPreconditioner();
501 :
502 : /**
503 : * If called with true this will add entries into the jacobian to link together degrees of freedom
504 : * that are found to
505 : * be related through the geometric search system.
506 : *
507 : * These entries are really only used by the Finite Difference Preconditioner and the constraint
508 : * system right now.
509 : */
510 1826 : void addImplicitGeometricCouplingEntriesToJacobian(bool add = true)
511 : {
512 1826 : _add_implicit_geometric_coupling_entries_to_jacobian = add;
513 1826 : }
514 :
515 : /**
516 : * Indicates whether to assemble residual and Jacobian after each constraint application.
517 : * When true, enables "transitive" constraint application: subsequent constraints can use prior
518 : * constraints' results.
519 : */
520 : void assembleConstraintsSeparately(bool separately = true)
521 : {
522 : _assemble_constraints_separately = separately;
523 : }
524 :
525 : /**
526 : * Attach a customized preconditioner that requires physics knowledge.
527 : * Generic preconditioners should be implemented in PETSc, instead.
528 : */
529 : virtual void attachPreconditioner(libMesh::Preconditioner<Number> * preconditioner) = 0;
530 :
531 : /**
532 : * Setup damping stuff (called before we actually start)
533 : */
534 : void setupDampers();
535 : /**
536 : * Compute the incremental change in variables at QPs for dampers. Called before we use damping
537 : * @param tid Thread ID
538 : * @param damped_vars Set of variables for which increment is to be computed
539 : */
540 : void reinitIncrementAtQpsForDampers(THREAD_ID tid, const std::set<MooseVariable *> & damped_vars);
541 :
542 : /**
543 : * Compute the incremental change in variables at nodes for dampers. Called before we use damping
544 : * @param tid Thread ID
545 : * @param damped_vars Set of variables for which increment is to be computed
546 : */
547 : void reinitIncrementAtNodeForDampers(THREAD_ID tid,
548 : const std::set<MooseVariable *> & damped_vars);
549 :
550 : ///@{
551 : /// System Integrity Checks
552 : void checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const;
553 : virtual bool containsTimeKernel() override;
554 : virtual std::vector<std::string> timeKernelVariableNames() override;
555 : ///@}
556 :
557 : /**
558 : * Return the number of non-linear iterations
559 : */
560 11651 : unsigned int nNonlinearIterations() const { return _n_iters; }
561 :
562 : /**
563 : * Return the number of linear iterations
564 : */
565 10748 : unsigned int nLinearIterations() const { return _n_linear_iters; }
566 :
567 : /**
568 : * Return the total number of residual evaluations done so far in this calculation
569 : */
570 316 : unsigned int nResidualEvaluations() const { return _n_residual_evaluations; }
571 :
572 : /**
573 : * Return the final nonlinear residual
574 : */
575 258 : Real finalNonlinearResidual() const { return _final_residual; }
576 :
577 : /**
578 : * Return the last nonlinear norm
579 : * @return A Real containing the last computed residual norm
580 : */
581 475770 : Real nonlinearNorm() const { return _last_nl_rnorm; }
582 :
583 : /**
584 : * Force the printing of all variable norms after each solve.
585 : * \todo{Remove after output update
586 : */
587 : void printAllVariableNorms(bool state) { _print_all_var_norms = state; }
588 :
589 40 : void debuggingResiduals(bool state) { _debugging_residuals = state; }
590 :
591 : unsigned int _num_residual_evaluations;
592 :
593 : void setPredictor(std::shared_ptr<Predictor> predictor);
594 66 : Predictor * getPredictor() { return _predictor.get(); }
595 :
596 : /**
597 : * Indicated whether this system needs material properties on boundaries.
598 : * @return Boolean if IntegratedBCs are active
599 : */
600 : bool needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const;
601 :
602 : /**
603 : * Indicated whether this system needs material properties on interfaces.
604 : * @return Boolean if IntegratedBCs are active
605 : */
606 : bool needInterfaceMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const;
607 :
608 : /**
609 : * Indicates whether this system needs material properties on internal sides.
610 : * @return Boolean if DGKernels are active
611 : */
612 : bool needInternalNeighborSideMaterial(SubdomainID subdomain_id, THREAD_ID tid) const;
613 :
614 : /**
615 : * Getter for _doing_dg
616 : */
617 : bool doingDG() const;
618 :
619 : //@{
620 : /**
621 : * Access functions to Warehouses from outside NonlinearSystemBase
622 : */
623 3934415 : MooseObjectTagWarehouse<KernelBase> & getKernelWarehouse() { return _kernels; }
624 20 : const MooseObjectTagWarehouse<KernelBase> & getKernelWarehouse() const { return _kernels; }
625 3928524 : MooseObjectTagWarehouse<DGKernelBase> & getDGKernelWarehouse() { return _dg_kernels; }
626 3928520 : MooseObjectTagWarehouse<InterfaceKernelBase> & getInterfaceKernelWarehouse()
627 : {
628 3928520 : return _interface_kernels;
629 : }
630 38401 : MooseObjectTagWarehouse<DiracKernelBase> & getDiracKernelWarehouse() { return _dirac_kernels; }
631 6052660 : MooseObjectTagWarehouse<IntegratedBCBase> & getIntegratedBCWarehouse() { return _integrated_bcs; }
632 56 : const MooseObjectTagWarehouse<ScalarKernelBase> & getScalarKernelWarehouse() const
633 : {
634 56 : return _scalar_kernels;
635 : }
636 4 : const MooseObjectTagWarehouse<NodalKernelBase> & getNodalKernelWarehouse() const
637 : {
638 4 : return _nodal_kernels;
639 : }
640 3867237 : MooseObjectTagWarehouse<HDGKernel> & getHDGKernelWarehouse() { return _hybridized_kernels; }
641 2883 : const MooseObjectWarehouse<ElementDamper> & getElementDamperWarehouse() const
642 : {
643 2883 : return _element_dampers;
644 : }
645 8701 : const MooseObjectWarehouse<NodalDamper> & getNodalDamperWarehouse() const
646 : {
647 8701 : return _nodal_dampers;
648 : }
649 : const ConstraintWarehouse & getConstraintWarehouse() const { return _constraints; }
650 :
651 : /**
652 : * Return the NodalBCBase warehouse
653 : */
654 3240368 : const MooseObjectTagWarehouse<NodalBCBase> & getNodalBCWarehouse() const { return _nodal_bcs; }
655 :
656 : /**
657 : * Return the IntegratedBCBase warehouse
658 : */
659 : const MooseObjectTagWarehouse<IntegratedBCBase> & getIntegratedBCWarehouse() const
660 : {
661 : return _integrated_bcs;
662 : }
663 :
664 : #ifdef MOOSE_KOKKOS_ENABLED
665 : ///@{
666 : /// Return the Kokkos residual object warehouses
667 : MooseObjectTagWarehouse<ResidualObject> & getKokkosKernelWarehouse() { return _kokkos_kernels; }
668 : MooseObjectTagWarehouse<ResidualObject> & getKokkosNodalKernelWarehouse()
669 : {
670 : return _kokkos_nodal_kernels;
671 : }
672 693 : MooseObjectTagWarehouse<ResidualObject> & getKokkosNodalBCWarehouse()
673 : {
674 693 : return _kokkos_nodal_bcs;
675 : }
676 : MooseObjectTagWarehouse<ResidualObject> & getKokkosIntegratedBCWarehouse()
677 : {
678 : return _kokkos_integrated_bcs;
679 : }
680 : ///@}
681 : #endif
682 :
683 : //@}
684 :
685 : /**
686 : * Weather or not the nonlinear system has save-ins
687 : */
688 3335767 : bool hasSaveIn() const { return _has_save_in || _has_nodalbc_save_in; }
689 :
690 : /**
691 : * Weather or not the nonlinear system has diagonal Jacobian save-ins
692 : */
693 526406 : bool hasDiagSaveIn() const { return _has_diag_save_in || _has_nodalbc_diag_save_in; }
694 :
695 562158256 : virtual libMesh::System & system() override { return _sys; }
696 1146774618 : virtual const libMesh::System & system() const override { return _sys; }
697 :
698 : virtual void setSolutionUDotOld(const NumericVector<Number> & u_dot_old);
699 :
700 : virtual void setSolutionUDotDotOld(const NumericVector<Number> & u_dotdot_old);
701 :
702 : virtual void setPreviousNewtonSolution(const NumericVector<Number> & soln);
703 :
704 84033 : TagID timeVectorTag() const override { return _Re_time_tag; }
705 55856 : TagID nonTimeVectorTag() const override { return _Re_non_time_tag; }
706 16601815 : TagID residualVectorTag() const override { return _Re_tag; }
707 1094855 : TagID systemMatrixTag() const override { return _Ke_system_tag; }
708 :
709 : /**
710 : * Call this method if you want the residual and Jacobian to be computed simultaneously
711 : */
712 : virtual void residualAndJacobianTogether() = 0;
713 :
714 : bool computeScalingOnce() const { return _compute_scaling_once; }
715 62587 : void computeScalingOnce(bool compute_scaling_once)
716 : {
717 62587 : _compute_scaling_once = compute_scaling_once;
718 62587 : }
719 :
720 : /**
721 : * Sets the param that indicates the weighting of the residual vs the Jacobian in determining
722 : * variable scaling parameters. A value of 1 indicates pure residual-based scaling. A value of 0
723 : * indicates pure Jacobian-based scaling
724 : */
725 62587 : void autoScalingParam(Real resid_vs_jac_scaling_param)
726 : {
727 62587 : _resid_vs_jac_scaling_param = resid_vs_jac_scaling_param;
728 62587 : }
729 :
730 17 : void scalingGroupVariables(const std::vector<std::vector<std::string>> & scaling_group_variables)
731 : {
732 17 : _scaling_group_variables = scaling_group_variables;
733 17 : }
734 :
735 : void
736 10 : ignoreVariablesForAutoscaling(const std::vector<std::string> & ignore_variables_for_autoscaling)
737 : {
738 10 : _ignore_variables_for_autoscaling = ignore_variables_for_autoscaling;
739 10 : }
740 :
741 80893 : bool offDiagonalsInAutoScaling() const { return _off_diagonals_in_auto_scaling; }
742 62587 : void offDiagonalsInAutoScaling(bool off_diagonals_in_auto_scaling)
743 : {
744 62587 : _off_diagonals_in_auto_scaling = off_diagonals_in_auto_scaling;
745 62587 : }
746 :
747 : libMesh::System & _sys;
748 : // FIXME: make these protected and create getters/setters
749 : Real _last_nl_rnorm;
750 : std::vector<unsigned int> _current_l_its;
751 : unsigned int _current_nl_its;
752 :
753 : /**
754 : * Setup the PETSc DM object (when appropriate)
755 : */
756 : void setupDM();
757 :
758 : using SystemBase::reinitNodeFace;
759 :
760 : /**
761 : * Create finite differencing contexts for assembly of the Jacobian and/or approximating the
762 : * action of the Jacobian on vectors (e.g. FD and/or MFFD respectively)
763 : */
764 0 : virtual void potentiallySetupFiniteDifferencing() {}
765 :
766 : /**
767 : * Destroy the coloring object if it exists
768 : */
769 : void destroyColoring();
770 :
771 : protected:
772 : /**
773 : * Compute the residual for a given tag
774 : * @param tags The tags of kernels for which the residual is to be computed.
775 : */
776 : void computeResidualInternal(const std::set<TagID> & tags);
777 :
778 : /**
779 : * Compute residual with Kokkos objects
780 : */
781 : #ifdef MOOSE_KOKKOS_ENABLED
782 : void computeKokkosResidual(const std::set<TagID> & tags);
783 : #endif
784 :
785 : /**
786 : * Enforces nodal boundary conditions. The boundary condition will be implemented
787 : * in the residual using all the tags in the system.
788 : */
789 : void computeNodalBCs(NumericVector<Number> & residual);
790 :
791 : /**
792 : * Form a residual for BCs that at least has one of the given tags.
793 : */
794 : void computeNodalBCs(NumericVector<Number> & residual, const std::set<TagID> & tags);
795 :
796 : /**
797 : * Form multiple tag-associated residual vectors for the given tags.
798 : */
799 : void computeNodalBCs(const std::set<TagID> & tags);
800 :
801 : /**
802 : * compute the residual and Jacobian for nodal boundary conditions
803 : */
804 : void computeNodalBCsResidualAndJacobian();
805 :
806 : /**
807 : * Form multiple matrices for all the tags. Users should not call this func directly.
808 : */
809 : void computeJacobianInternal(const std::set<TagID> & tags);
810 :
811 : /**
812 : * Compute Jacobian with Kokkos objects
813 : */
814 : #ifdef MOOSE_KOKKOS_ENABLED
815 : void computeKokkosJacobian(const std::set<TagID> & tags);
816 : #endif
817 :
818 : void computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian);
819 :
820 : void computeScalarKernelsJacobians(const std::set<TagID> & tags);
821 :
822 : /**
823 : * Enforce nodal constraints
824 : */
825 : void enforceNodalConstraintsResidual(NumericVector<Number> & residual);
826 : void enforceNodalConstraintsJacobian();
827 :
828 : /**
829 : * Do mortar constraint residual/jacobian computations
830 : */
831 : void mortarConstraints(Moose::ComputeType compute_type,
832 : const std::set<TagID> & vector_tags,
833 : const std::set<TagID> & matrix_tags);
834 :
835 : /**
836 : * Compute a "Jacobian" for automatic scaling purposes
837 : */
838 : virtual void computeScalingJacobian() = 0;
839 :
840 : /**
841 : * Compute a "residual" for automatic scaling purposes
842 : */
843 : virtual void computeScalingResidual() = 0;
844 :
845 : /**
846 : * Assemble the numeric vector of scaling factors such that it can be used during assembly of the
847 : * system matrix
848 : */
849 : void assembleScalingVector();
850 :
851 : /**
852 : * Called after any ResidualObject-derived objects are added
853 : * to the system.
854 : */
855 182028 : virtual void postAddResidualObject(ResidualObject &) {}
856 :
857 : /**
858 : * Reinitialize quantities such as variables, residuals, Jacobians, materials for node-face
859 : * constraints
860 : */
861 : void reinitNodeFace(const Node & secondary_node,
862 : const BoundaryID secondary_boundary,
863 : const PenetrationInfo & info,
864 : const bool displaced);
865 :
866 : /**
867 : * Perform some steps to get ready for the solver. These include
868 : * - zeroing iteration counters
869 : * - setting initial solutions
870 : * - possibly performing automatic scaling
871 : * - forming a scaling vector which, at least at some point, was required when AD objects were
872 : * used with non-unity scaling factors for nonlinear variables
873 : * @returns Whether any exceptions were raised while running this method
874 : */
875 : bool preSolve();
876 :
877 : /// ghosted form of the residual
878 : NumericVector<Number> * _residual_ghosted;
879 :
880 : /// Copy of the residual vector, or nullptr if a copy is not needed
881 : std::unique_ptr<NumericVector<Number>> _residual_copy;
882 :
883 : /// \f$ {du^dot}\over{du} \f$
884 : Number _du_dot_du;
885 : /// \f$ {du^dotdot}\over{du} \f$
886 : Number _du_dotdot_du;
887 :
888 : /// Tag for time contribution residual
889 : TagID _Re_time_tag;
890 :
891 : /// Vector tags to temporarily store all tags associated with the current system.
892 : std::set<TagID> _nl_vector_tags;
893 :
894 : /// Matrix tags to temporarily store all tags associated with the current system.
895 : std::set<TagID> _nl_matrix_tags;
896 :
897 : /// residual vector for time contributions
898 : NumericVector<Number> * _Re_time;
899 :
900 : /// Tag for non-time contribution residual
901 : TagID _Re_non_time_tag;
902 : /// residual vector for non-time contributions
903 : NumericVector<Number> * _Re_non_time;
904 :
905 : /// Used for the residual vector from PETSc
906 : TagID _Re_tag;
907 :
908 : /// Tag for non-time contribution Jacobian
909 : TagID _Ke_non_time_tag;
910 :
911 : /// Tag for system contribution Jacobian
912 : TagID _Ke_system_tag;
913 :
914 : ///@{
915 : /// Kernel Storage
916 : MooseObjectTagWarehouse<KernelBase> _kernels;
917 : MooseObjectTagWarehouse<HDGKernel> _hybridized_kernels;
918 : MooseObjectTagWarehouse<ScalarKernelBase> _scalar_kernels;
919 : MooseObjectTagWarehouse<DGKernelBase> _dg_kernels;
920 : MooseObjectTagWarehouse<InterfaceKernelBase> _interface_kernels;
921 : ///@}
922 :
923 : ///@{
924 : /// BoundaryCondition Warhouses
925 : MooseObjectTagWarehouse<IntegratedBCBase> _integrated_bcs;
926 : MooseObjectTagWarehouse<NodalBCBase> _nodal_bcs;
927 : MooseObjectWarehouse<DirichletBCBase> _preset_nodal_bcs;
928 : MooseObjectWarehouse<ADDirichletBCBase> _ad_preset_nodal_bcs;
929 : ///@}
930 :
931 : #ifdef MOOSE_KOKKOS_ENABLED
932 : ///@{
933 : /// Kokkos residual object warhouses
934 : MooseObjectTagWarehouse<ResidualObject> _kokkos_kernels;
935 : MooseObjectTagWarehouse<ResidualObject> _kokkos_integrated_bcs;
936 : MooseObjectTagWarehouse<ResidualObject> _kokkos_nodal_bcs;
937 : MooseObjectWarehouse<ResidualObject> _kokkos_preset_nodal_bcs;
938 : MooseObjectTagWarehouse<ResidualObject> _kokkos_nodal_kernels;
939 : ///@}
940 : #endif
941 :
942 : /// Dirac Kernel storage for each thread
943 : MooseObjectTagWarehouse<DiracKernelBase> _dirac_kernels;
944 :
945 : /// Element Dampers for each thread
946 : MooseObjectWarehouse<ElementDamper> _element_dampers;
947 :
948 : /// Nodal Dampers for each thread
949 : MooseObjectWarehouse<NodalDamper> _nodal_dampers;
950 :
951 : /// General Dampers
952 : MooseObjectWarehouse<GeneralDamper> _general_dampers;
953 :
954 : /// NodalKernels for each thread
955 : MooseObjectTagWarehouse<NodalKernelBase> _nodal_kernels;
956 :
957 : /// Decomposition splits
958 : MooseObjectWarehouseBase<Split> _splits; // use base b/c there are no setup methods
959 :
960 : /// Constraints storage object
961 : ConstraintWarehouse _constraints;
962 :
963 : /// increment vector
964 : NumericVector<Number> * _increment_vec;
965 : /// Preconditioner
966 : std::shared_ptr<MoosePreconditioner> _preconditioner;
967 :
968 : /// Whether or not to use a finite differenced preconditioner
969 : bool _use_finite_differenced_preconditioner;
970 :
971 : MatFDColoring _fdcoloring;
972 :
973 : /// The field split preconditioner if this sytem is using one
974 : FieldSplitPreconditionerBase * _fsp;
975 :
976 : /// Whether or not to add implicit geometric couplings to the Jacobian for FDP
977 : bool _add_implicit_geometric_coupling_entries_to_jacobian;
978 :
979 : /// Whether or not to assemble the residual and Jacobian after the application of each constraint.
980 : bool _assemble_constraints_separately;
981 :
982 : /// Whether or not a ghosted copy of the residual needs to be made
983 : bool _need_residual_ghosted;
984 : /// true if debugging residuals
985 : bool _debugging_residuals;
986 :
987 : /// true if DG is active (optimization reasons)
988 : bool _doing_dg;
989 :
990 : /// vectors that will be zeroed before a residual computation
991 : std::vector<std::string> _vecs_to_zero_for_residual;
992 :
993 : unsigned int _n_iters;
994 : unsigned int _n_linear_iters;
995 :
996 : /// Total number of residual evaluations that have been performed
997 : unsigned int _n_residual_evaluations;
998 :
999 : Real _final_residual;
1000 :
1001 : /// If predictor is active, this is non-NULL
1002 : std::shared_ptr<Predictor> _predictor;
1003 :
1004 : bool _computing_pre_smo_residual;
1005 :
1006 : /// The pre-SMO residual, see setPreSMOResidual for a detailed explanation
1007 : Real _pre_smo_residual;
1008 : /// The initial (i.e., 0th nonlinear iteration) residual, see setPreSMOResidual for a detailed explanation
1009 : Real _initial_residual;
1010 : /// Whether to use the pre-SMO initial residual in the relative convergence check
1011 : bool _use_pre_smo_residual;
1012 :
1013 : bool _print_all_var_norms;
1014 :
1015 : /// If there is any Kernel or IntegratedBC having save_in
1016 : bool _has_save_in;
1017 :
1018 : /// If there is any Kernel or IntegratedBC having diag_save_in
1019 : bool _has_diag_save_in;
1020 :
1021 : /// If there is a nodal BC having save_in
1022 : bool _has_nodalbc_save_in;
1023 :
1024 : /// If there is a nodal BC having diag_save_in
1025 : bool _has_nodalbc_diag_save_in;
1026 :
1027 : void getNodeDofs(dof_id_type node_id, std::vector<dof_id_type> & dofs);
1028 :
1029 : /// Flag used to indicate whether we have already computed the scaling Jacobian
1030 : bool _computed_scaling;
1031 :
1032 : /// Whether the scaling factors should only be computed once at the beginning of the simulation
1033 : /// through an extra Jacobian evaluation. If this is set to false, then the scaling factors will
1034 : /// be computed during an extra Jacobian evaluation at the beginning of every time step.
1035 : bool _compute_scaling_once;
1036 :
1037 : /// The param that indicates the weighting of the residual vs the Jacobian in determining
1038 : /// variable scaling parameters. A value of 1 indicates pure residual-based scaling. A value of 0
1039 : /// indicates pure Jacobian-based scaling
1040 : Real _resid_vs_jac_scaling_param;
1041 :
1042 : /// A container of variable groupings that can be used in scaling calculations. This can be useful
1043 : /// for simulations in which vector-like variables are split into invidual scalar-field components
1044 : /// like for solid/fluid mechanics
1045 : std::vector<std::vector<std::string>> _scaling_group_variables;
1046 :
1047 : /// Container to hold flag if variable is to participate in autoscaling
1048 : std::vector<bool> _variable_autoscaled;
1049 :
1050 : /// A container for variables that do not partipate in autoscaling
1051 : std::vector<std::string> _ignore_variables_for_autoscaling;
1052 :
1053 : /// Whether to include off diagonals when determining automatic scaling factors
1054 : bool _off_diagonals_in_auto_scaling;
1055 :
1056 : /// A diagonal matrix used for computing scaling
1057 : std::unique_ptr<libMesh::DiagonalMatrix<Number>> _scaling_matrix;
1058 :
1059 : private:
1060 : /**
1061 : * Finds the implicit sparsity graph between geometrically related dofs.
1062 : */
1063 : void findImplicitGeometricCouplingEntries(
1064 : GeometricSearchData & geom_search_data,
1065 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> & graph);
1066 :
1067 : /**
1068 : * Setup group scaling containers
1069 : */
1070 : void setupScalingData();
1071 :
1072 : /// Functors for computing undisplaced mortar constraints
1073 : std::unordered_map<std::pair<BoundaryID, BoundaryID>, ComputeMortarFunctor>
1074 : _undisplaced_mortar_functors;
1075 :
1076 : /// Functors for computing displaced mortar constraints
1077 : std::unordered_map<std::pair<BoundaryID, BoundaryID>, ComputeMortarFunctor>
1078 : _displaced_mortar_functors;
1079 :
1080 : /// The current states of the solution (0 = current, 1 = old, etc)
1081 : std::vector<NumericVector<Number> *> _solution_state;
1082 :
1083 : /// Whether we've initialized the automatic scaling data structures
1084 : bool _auto_scaling_initd;
1085 :
1086 : /// A map from variable index to group variable index and it's associated (inverse) scaling factor
1087 : std::unordered_map<unsigned int, unsigned int> _var_to_group_var;
1088 :
1089 : /// The number of scaling groups
1090 : std::size_t _num_scaling_groups;
1091 : };
|
