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 761344 : 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 53935 : 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 58927 : void setPreSMOResidual(bool use) { _use_pre_smo_residual = use; }
287 :
288 : /// Whether we are using pre-SMO residual in relative convergence checks
289 403615 : 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 : #ifdef MOOSE_KOKKOS_ENABLED
351 : void computeKokkosResidualAndJacobian(const std::set<TagID> & vector_tags,
352 : const std::set<TagID> & matrix_tags);
353 : #endif
354 :
355 : /**
356 : * Form a residual vector for a given tag
357 : */
358 : void computeResidual(NumericVector<Number> & residual, TagID tag_id);
359 :
360 : /**
361 : * Adds entries to the Jacobian in the correct positions for couplings coming from dofs being
362 : * coupled that
363 : * are related geometrically (i.e. near each other across a gap).
364 : */
365 : void addImplicitGeometricCouplingEntries(GeometricSearchData & geom_search_data);
366 :
367 : /**
368 : * Add jacobian contributions from Constraints
369 : *
370 : * @param jacobian reference to a read-only view of the Jacobian matrix
371 : * @param displaced Controls whether to do the displaced Constraints or non-displaced
372 : */
373 : void constraintJacobians(const SparseMatrix<Number> & jacobian_to_view, bool displaced);
374 :
375 : /**
376 : * Computes multiple (tag associated) Jacobian matricese
377 : */
378 : void computeJacobianTags(const std::set<TagID> & tags);
379 :
380 : /**
381 : * Method used to obtain scaling factors for variables
382 : * @returns whether this method ran without exceptions
383 : */
384 : bool computeScaling();
385 :
386 : /**
387 : * Associate jacobian to systemMatrixTag, and then form a matrix for all the tags
388 : */
389 : void computeJacobian(libMesh::SparseMatrix<Number> & jacobian, const std::set<TagID> & tags);
390 :
391 : /**
392 : * Take all tags in the system, and form a matrix for all tags in the system
393 : */
394 : void computeJacobian(libMesh::SparseMatrix<Number> & jacobian);
395 :
396 : /**
397 : * Computes several Jacobian blocks simultaneously, summing their contributions into smaller
398 : * preconditioning matrices.
399 : *
400 : * Used by Physics-based preconditioning
401 : *
402 : * @param blocks The blocks to fill in (JacobianBlock is defined in ComputeJacobianBlocksThread)
403 : */
404 : void computeJacobianBlocks(std::vector<JacobianBlock *> & blocks);
405 :
406 : void computeJacobianBlocks(std::vector<JacobianBlock *> & blocks, const std::set<TagID> & tags);
407 :
408 : /**
409 : * Compute damping
410 : * @param solution The trail solution vector
411 : * @param update The incremental update to the solution vector
412 : * @return returns The damping factor
413 : */
414 : Real computeDamping(const NumericVector<Number> & solution, const NumericVector<Number> & update);
415 :
416 : /**
417 : * Called at the beginning of the time step
418 : */
419 : void onTimestepBegin();
420 :
421 : using SystemBase::subdomainSetup;
422 : /**
423 : * Called from assembling when we hit a new subdomain
424 : * @param subdomain ID of the new subdomain
425 : * @param tid Thread ID
426 : */
427 : virtual void subdomainSetup(SubdomainID subdomain, THREAD_ID tid);
428 :
429 : /**
430 : * Called from explicit time stepping to overwrite boundary positions (explicit dynamics). This
431 : * will close/assemble the passed-in \p soln after overwrite
432 : */
433 : void overwriteNodeFace(NumericVector<Number> & soln);
434 :
435 : /**
436 : * Update active objects of Warehouses owned by NonlinearSystemBase
437 : */
438 : void updateActive(THREAD_ID tid);
439 :
440 : /**
441 : * Set transient term used by residual and Jacobian evaluation.
442 : * @param udot transient term
443 : * @note If the calling sequence for residual evaluation was changed, this could become an
444 : * explicit argument.
445 : */
446 : virtual void setSolutionUDot(const NumericVector<Number> & udot);
447 :
448 : /**
449 : * Set transient term used by residual and Jacobian evaluation.
450 : * @param udotdot transient term
451 : * @note If the calling sequence for residual evaluation was changed, this could become an
452 : * explicit argument.
453 : */
454 : virtual void setSolutionUDotDot(const NumericVector<Number> & udotdot);
455 :
456 : /**
457 : * Return a numeric vector that is associated with the time tag.
458 : */
459 : NumericVector<Number> & getResidualTimeVector();
460 :
461 : /**
462 : * Return a numeric vector that is associated with the nontime tag.
463 : */
464 : NumericVector<Number> & getResidualNonTimeVector();
465 :
466 : /**
467 : * Return a residual vector that is associated with the residual tag.
468 : */
469 : NumericVector<Number> & residualVector(TagID tag);
470 :
471 : virtual NumericVector<Number> & residualCopy() override;
472 : virtual NumericVector<Number> & residualGhosted() override;
473 :
474 : virtual NumericVector<Number> & RHS() = 0;
475 :
476 : virtual void augmentSparsity(libMesh::SparsityPattern::Graph & sparsity,
477 : std::vector<dof_id_type> & n_nz,
478 : std::vector<dof_id_type> & n_oz) override;
479 :
480 : /**
481 : * Sets a preconditioner
482 : * @param pc The preconditioner to be set
483 : */
484 : void setPreconditioner(std::shared_ptr<MoosePreconditioner> pc);
485 : MoosePreconditioner const * getPreconditioner() const;
486 :
487 : /**
488 : * If called with true this system will use a finite differenced form of
489 : * the Jacobian as the preconditioner
490 : */
491 86 : void useFiniteDifferencedPreconditioner(bool use = true)
492 : {
493 86 : _use_finite_differenced_preconditioner = use;
494 86 : }
495 :
496 : /**
497 : * If called with a non-null object true this system will use a field split preconditioner matrix.
498 : */
499 123 : void useFieldSplitPreconditioner(FieldSplitPreconditionerBase * fsp) { _fsp = fsp; }
500 :
501 : /**
502 : * @returns A field split preconditioner. This will error if there is no field split
503 : * preconditioner
504 : */
505 : FieldSplitPreconditionerBase & getFieldSplitPreconditioner();
506 :
507 : /**
508 : * If called with true this will add entries into the jacobian to link together degrees of freedom
509 : * that are found to
510 : * be related through the geometric search system.
511 : *
512 : * These entries are really only used by the Finite Difference Preconditioner and the constraint
513 : * system right now.
514 : */
515 1305 : void addImplicitGeometricCouplingEntriesToJacobian(bool add = true)
516 : {
517 1305 : _add_implicit_geometric_coupling_entries_to_jacobian = add;
518 1305 : }
519 :
520 : /**
521 : * Indicates whether to assemble residual and Jacobian after each constraint application.
522 : * When true, enables "transitive" constraint application: subsequent constraints can use prior
523 : * constraints' results.
524 : */
525 : void assembleConstraintsSeparately(bool separately = true)
526 : {
527 : _assemble_constraints_separately = separately;
528 : }
529 :
530 : /**
531 : * Attach a customized preconditioner that requires physics knowledge.
532 : * Generic preconditioners should be implemented in PETSc, instead.
533 : */
534 : virtual void attachPreconditioner(libMesh::Preconditioner<Number> * preconditioner) = 0;
535 :
536 : /**
537 : * Setup damping stuff (called before we actually start)
538 : */
539 : void setupDampers();
540 : /**
541 : * Compute the incremental change in variables at QPs for dampers. Called before we use damping
542 : * @param tid Thread ID
543 : * @param damped_vars Set of variables for which increment is to be computed
544 : */
545 : void reinitIncrementAtQpsForDampers(THREAD_ID tid, const std::set<MooseVariable *> & damped_vars);
546 :
547 : /**
548 : * Compute the incremental change in variables at nodes for dampers. Called before we use damping
549 : * @param tid Thread ID
550 : * @param damped_vars Set of variables for which increment is to be computed
551 : */
552 : void reinitIncrementAtNodeForDampers(THREAD_ID tid,
553 : const std::set<MooseVariable *> & damped_vars);
554 :
555 : ///@{
556 : /// System Integrity Checks
557 : void checkKernelCoverage(const std::set<SubdomainID> & mesh_subdomains) const;
558 : virtual bool containsTimeKernel() override;
559 : virtual std::vector<std::string> timeKernelVariableNames() override;
560 : ///@}
561 :
562 : /**
563 : * Return the number of non-linear iterations
564 : */
565 9797 : virtual unsigned int nNonlinearIterations() const { return _n_iters; }
566 :
567 : /**
568 : * Return the number of linear iterations
569 : */
570 8965 : virtual unsigned int nLinearIterations() const { return _n_linear_iters; }
571 :
572 : /**
573 : * Return the total number of residual evaluations done so far in this calculation
574 : */
575 152 : unsigned int nResidualEvaluations() const { return _n_residual_evaluations; }
576 :
577 : /**
578 : * Return the final nonlinear residual
579 : */
580 236 : virtual Real finalNonlinearResidual() const { return _final_residual; }
581 :
582 : /**
583 : * Return the last nonlinear norm
584 : * @return A Real containing the last computed residual norm
585 : */
586 435910 : Real nonlinearNorm() const { return _last_nl_rnorm; }
587 :
588 : /**
589 : * Force the printing of all variable norms after each solve.
590 : * \todo{Remove after output update
591 : */
592 : void printAllVariableNorms(bool state) { _print_all_var_norms = state; }
593 :
594 36 : void debuggingResiduals(bool state) { _debugging_residuals = state; }
595 :
596 : unsigned int _num_residual_evaluations;
597 :
598 : void setPredictor(std::shared_ptr<Predictor> predictor);
599 60 : Predictor * getPredictor() { return _predictor.get(); }
600 :
601 : /**
602 : * Indicated whether this system needs material properties on boundaries.
603 : * @return Boolean if IntegratedBCs are active
604 : */
605 : bool needBoundaryMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const;
606 :
607 : /**
608 : * Indicated whether this system needs material properties on interfaces.
609 : * @return Boolean if IntegratedBCs are active
610 : */
611 : bool needInterfaceMaterialOnSide(BoundaryID bnd_id, THREAD_ID tid) const;
612 :
613 : /**
614 : * Indicates whether this system needs material properties on internal sides.
615 : * @return Boolean if DGKernels are active
616 : */
617 : bool needInternalNeighborSideMaterial(SubdomainID subdomain_id, THREAD_ID tid) const;
618 :
619 : /**
620 : * Getter for _doing_dg
621 : */
622 : bool doingDG() const;
623 :
624 : //@{
625 : /**
626 : * Access functions to Warehouses from outside NonlinearSystemBase
627 : */
628 3584861 : MooseObjectTagWarehouse<KernelBase> & getKernelWarehouse() { return _kernels; }
629 18 : const MooseObjectTagWarehouse<KernelBase> & getKernelWarehouse() const { return _kernels; }
630 3578713 : MooseObjectTagWarehouse<DGKernelBase> & getDGKernelWarehouse() { return _dg_kernels; }
631 3578693 : MooseObjectTagWarehouse<InterfaceKernelBase> & getInterfaceKernelWarehouse()
632 : {
633 3578693 : return _interface_kernels;
634 : }
635 35512 : MooseObjectTagWarehouse<DiracKernelBase> & getDiracKernelWarehouse() { return _dirac_kernels; }
636 5486275 : MooseObjectTagWarehouse<IntegratedBCBase> & getIntegratedBCWarehouse() { return _integrated_bcs; }
637 51 : const MooseObjectTagWarehouse<ScalarKernelBase> & getScalarKernelWarehouse() const
638 : {
639 51 : return _scalar_kernels;
640 : }
641 20 : const MooseObjectTagWarehouse<NodalKernelBase> & getNodalKernelWarehouse() const
642 : {
643 20 : return _nodal_kernels;
644 : }
645 3519717 : MooseObjectTagWarehouse<HDGKernel> & getHDGKernelWarehouse() { return _hybridized_kernels; }
646 1170 : const MooseObjectWarehouse<ElementDamper> & getElementDamperWarehouse() const
647 : {
648 1170 : return _element_dampers;
649 : }
650 1801 : const MooseObjectWarehouse<NodalDamper> & getNodalDamperWarehouse() const
651 : {
652 1801 : return _nodal_dampers;
653 : }
654 : const ConstraintWarehouse & getConstraintWarehouse() const { return _constraints; }
655 :
656 : /**
657 : * Return the NodalBCBase warehouse
658 : */
659 3053091 : const MooseObjectTagWarehouse<NodalBCBase> & getNodalBCWarehouse() const { return _nodal_bcs; }
660 :
661 : /**
662 : * Return the IntegratedBCBase warehouse
663 : */
664 : const MooseObjectTagWarehouse<IntegratedBCBase> & getIntegratedBCWarehouse() const
665 : {
666 : return _integrated_bcs;
667 : }
668 :
669 : #ifdef MOOSE_KOKKOS_ENABLED
670 : ///@{
671 : /// Return the Kokkos residual object warehouses
672 19 : MooseObjectTagWarehouse<ResidualObject> & getKokkosKernelWarehouse() { return _kokkos_kernels; }
673 2 : MooseObjectTagWarehouse<ResidualObject> & getKokkosNodalKernelWarehouse()
674 : {
675 2 : return _kokkos_nodal_kernels;
676 : }
677 2265 : MooseObjectTagWarehouse<ResidualObject> & getKokkosNodalBCWarehouse()
678 : {
679 2265 : return _kokkos_nodal_bcs;
680 : }
681 : MooseObjectTagWarehouse<ResidualObject> & getKokkosIntegratedBCWarehouse()
682 : {
683 : return _kokkos_integrated_bcs;
684 : }
685 : ///@}
686 : #endif
687 :
688 : //@}
689 :
690 : /**
691 : * Weather or not the nonlinear system has save-ins
692 : */
693 3036277 : bool hasSaveIn() const { return _has_save_in || _has_nodalbc_save_in; }
694 :
695 : /**
696 : * Weather or not the nonlinear system has diagonal Jacobian save-ins
697 : */
698 472515 : bool hasDiagSaveIn() const { return _has_diag_save_in || _has_nodalbc_diag_save_in; }
699 :
700 483326616 : virtual libMesh::System & system() override { return _sys; }
701 726751697 : virtual const libMesh::System & system() const override { return _sys; }
702 :
703 : virtual void setSolutionUDotOld(const NumericVector<Number> & u_dot_old);
704 :
705 : virtual void setSolutionUDotDotOld(const NumericVector<Number> & u_dotdot_old);
706 :
707 : virtual void setPreviousNewtonSolution(const NumericVector<Number> & soln);
708 :
709 91042 : TagID timeVectorTag() const override { return _Re_time_tag; }
710 64407 : TagID nonTimeVectorTag() const override { return _Re_non_time_tag; }
711 15140990 : TagID residualVectorTag() const override { return _Re_tag; }
712 985880 : TagID systemMatrixTag() const override { return _Ke_system_tag; }
713 :
714 : /**
715 : * Call this method if you want the residual and Jacobian to be computed simultaneously
716 : */
717 : virtual void residualAndJacobianTogether() = 0;
718 :
719 : bool computeScalingOnce() const { return _compute_scaling_once; }
720 58921 : void computeScalingOnce(bool compute_scaling_once)
721 : {
722 58921 : _compute_scaling_once = compute_scaling_once;
723 58921 : }
724 :
725 : /**
726 : * Sets the param that indicates the weighting of the residual vs the Jacobian in determining
727 : * variable scaling parameters. A value of 1 indicates pure residual-based scaling. A value of 0
728 : * indicates pure Jacobian-based scaling
729 : */
730 58921 : void autoScalingParam(Real resid_vs_jac_scaling_param)
731 : {
732 58921 : _resid_vs_jac_scaling_param = resid_vs_jac_scaling_param;
733 58921 : }
734 :
735 15 : void scalingGroupVariables(const std::vector<std::vector<std::string>> & scaling_group_variables)
736 : {
737 15 : _scaling_group_variables = scaling_group_variables;
738 15 : }
739 :
740 : void
741 9 : ignoreVariablesForAutoscaling(const std::vector<std::string> & ignore_variables_for_autoscaling)
742 : {
743 9 : _ignore_variables_for_autoscaling = ignore_variables_for_autoscaling;
744 9 : }
745 :
746 69804 : bool offDiagonalsInAutoScaling() const { return _off_diagonals_in_auto_scaling; }
747 58921 : void offDiagonalsInAutoScaling(bool off_diagonals_in_auto_scaling)
748 : {
749 58921 : _off_diagonals_in_auto_scaling = off_diagonals_in_auto_scaling;
750 58921 : }
751 :
752 : libMesh::System & _sys;
753 : // FIXME: make these protected and create getters/setters
754 : Real _last_nl_rnorm;
755 : std::vector<unsigned int> _current_l_its;
756 : unsigned int _current_nl_its;
757 :
758 : /**
759 : * Setup the PETSc DM object (when appropriate)
760 : */
761 : void setupDM();
762 :
763 : using SystemBase::reinitNodeFace;
764 :
765 : /**
766 : * Create finite differencing contexts for assembly of the Jacobian and/or approximating the
767 : * action of the Jacobian on vectors (e.g. FD and/or MFFD respectively)
768 : */
769 0 : virtual void potentiallySetupFiniteDifferencing() {}
770 :
771 : /**
772 : * Destroy the coloring object if it exists
773 : */
774 : void destroyColoring();
775 :
776 : protected:
777 : /**
778 : * Compute the residual for a given tag
779 : * @param tags The tags of kernels for which the residual is to be computed.
780 : */
781 : void computeResidualInternal(const std::set<TagID> & tags);
782 :
783 : #ifdef MOOSE_KOKKOS_ENABLED
784 : /**
785 : * Compute residual with Kokkos objects
786 : */
787 : void computeKokkosResidual(const std::set<TagID> & tags);
788 : /**
789 : * Compute Kokkos nodal BCs
790 : */
791 : void computeKokkosNodalBCsResidual(const std::set<TagID> & tags);
792 : #endif
793 :
794 : /**
795 : * Enforces nodal boundary conditions. The boundary condition will be implemented
796 : * in the residual using all the tags in the system.
797 : */
798 : void computeNodalBCsResidual(NumericVector<Number> & residual);
799 :
800 : /**
801 : * Form a residual for BCs that at least has one of the given tags.
802 : */
803 : void computeNodalBCsResidual(NumericVector<Number> & residual, const std::set<TagID> & tags);
804 :
805 : /**
806 : * Form multiple tag-associated residual vectors for the given tags.
807 : */
808 : void computeNodalBCsResidual(const std::set<TagID> & tags);
809 :
810 : /**
811 : * Compute the Jacobian for nodal boundary conditions
812 : */
813 : void computeNodalBCsJacobian(const std::set<TagID> & tags);
814 :
815 : /**
816 : * Compute the residual and Jacobian together for nodal boundary conditions
817 : */
818 : void computeNodalBCsResidualAndJacobian(const std::set<TagID> & vector_tags,
819 : const std::set<TagID> & matrix_tags);
820 :
821 : /**
822 : * Form multiple matrices for all the tags. Users should not call this func directly.
823 : */
824 : void computeJacobianInternal(const std::set<TagID> & tags);
825 :
826 : /**
827 : * Compute Jacobian with Kokkos objects
828 : */
829 : #ifdef MOOSE_KOKKOS_ENABLED
830 : void computeKokkosJacobian(const std::set<TagID> & tags);
831 : #endif
832 :
833 : void computeDiracContributions(const std::set<TagID> & tags, bool is_jacobian);
834 :
835 : void computeScalarKernelsJacobians(const std::set<TagID> & tags);
836 :
837 : /**
838 : * Enforce nodal constraints
839 : */
840 : void enforceNodalConstraintsResidual(NumericVector<Number> & residual);
841 :
842 : /**
843 : * Enforce nodal constraints in the Jacobian
844 : * @param jacobian The Jacobian to read from while constructing the Jacobians corresponding to the
845 : * nodal constraints
846 : * @returns Whether there were active nodal constraints
847 : */
848 : bool enforceNodalConstraintsJacobian(const SparseMatrix<Number> & jacobian);
849 :
850 : /**
851 : * Do mortar constraint residual/jacobian computations
852 : */
853 : void mortarConstraints(Moose::ComputeType compute_type,
854 : const std::set<TagID> & vector_tags,
855 : const std::set<TagID> & matrix_tags);
856 :
857 : /**
858 : * Compute a "Jacobian" for automatic scaling purposes
859 : */
860 : virtual void computeScalingJacobian() = 0;
861 :
862 : /**
863 : * Compute a "residual" for automatic scaling purposes
864 : */
865 : virtual void computeScalingResidual() = 0;
866 :
867 : /**
868 : * Assemble the numeric vector of scaling factors such that it can be used during assembly of the
869 : * system matrix
870 : */
871 : void assembleScalingVector();
872 :
873 : /**
874 : * Called after any ResidualObject-derived objects are added
875 : * to the system.
876 : */
877 170191 : virtual void postAddResidualObject(ResidualObject &) {}
878 :
879 : /**
880 : * Reinitialize quantities such as variables, residuals, Jacobians, materials for node-face
881 : * constraints
882 : */
883 : void reinitNodeFace(const Node & secondary_node,
884 : const BoundaryID secondary_boundary,
885 : const PenetrationInfo & info,
886 : const bool displaced);
887 :
888 : /**
889 : * Perform some steps to get ready for the solver. These include
890 : * - zeroing iteration counters
891 : * - setting initial solutions
892 : * - possibly performing automatic scaling
893 : * - forming a scaling vector which, at least at some point, was required when AD objects were
894 : * used with non-unity scaling factors for nonlinear variables
895 : * @returns Whether any exceptions were raised while running this method
896 : */
897 : bool preSolve();
898 :
899 : /// ghosted form of the residual
900 : NumericVector<Number> * _residual_ghosted;
901 :
902 : /// Copy of the residual vector, or nullptr if a copy is not needed
903 : std::unique_ptr<NumericVector<Number>> _residual_copy;
904 :
905 : /// \f$ {du^dot}\over{du} \f$
906 : Number _du_dot_du;
907 : /// \f$ {du^dotdot}\over{du} \f$
908 : Number _du_dotdot_du;
909 :
910 : /// Tag for time contribution residual
911 : TagID _Re_time_tag;
912 :
913 : /// Vector tags to temporarily store all tags associated with the current system.
914 : std::set<TagID> _nl_vector_tags;
915 :
916 : /// Matrix tags to temporarily store all tags associated with the current system.
917 : std::set<TagID> _nl_matrix_tags;
918 :
919 : /// residual vector for time contributions
920 : NumericVector<Number> * _Re_time;
921 :
922 : /// Tag for non-time contribution residual
923 : TagID _Re_non_time_tag;
924 : /// residual vector for non-time contributions
925 : NumericVector<Number> * _Re_non_time;
926 :
927 : /// Used for the residual vector from PETSc
928 : TagID _Re_tag;
929 :
930 : /// Tag for non-time contribution Jacobian
931 : TagID _Ke_non_time_tag;
932 :
933 : /// Tag for system contribution Jacobian
934 : TagID _Ke_system_tag;
935 :
936 : ///@{
937 : /// Kernel Storage
938 : MooseObjectTagWarehouse<KernelBase> _kernels;
939 : MooseObjectTagWarehouse<HDGKernel> _hybridized_kernels;
940 : MooseObjectTagWarehouse<ScalarKernelBase> _scalar_kernels;
941 : MooseObjectTagWarehouse<DGKernelBase> _dg_kernels;
942 : MooseObjectTagWarehouse<InterfaceKernelBase> _interface_kernels;
943 : ///@}
944 :
945 : ///@{
946 : /// BoundaryCondition Warhouses
947 : MooseObjectTagWarehouse<IntegratedBCBase> _integrated_bcs;
948 : MooseObjectTagWarehouse<NodalBCBase> _nodal_bcs;
949 : MooseObjectWarehouse<DirichletBCBase> _preset_nodal_bcs;
950 : MooseObjectWarehouse<ADDirichletBCBase> _ad_preset_nodal_bcs;
951 : ///@}
952 :
953 : #ifdef MOOSE_KOKKOS_ENABLED
954 : ///@{
955 : /// Kokkos residual object warhouses
956 : MooseObjectTagWarehouse<ResidualObject> _kokkos_kernels;
957 : MooseObjectTagWarehouse<ResidualObject> _kokkos_integrated_bcs;
958 : MooseObjectTagWarehouse<ResidualObject> _kokkos_nodal_bcs;
959 : MooseObjectWarehouse<ResidualObject> _kokkos_preset_nodal_bcs;
960 : MooseObjectTagWarehouse<ResidualObject> _kokkos_nodal_kernels;
961 : ///@}
962 : #endif
963 :
964 : /// Dirac Kernel storage for each thread
965 : MooseObjectTagWarehouse<DiracKernelBase> _dirac_kernels;
966 :
967 : /// Element Dampers for each thread
968 : MooseObjectWarehouse<ElementDamper> _element_dampers;
969 :
970 : /// Nodal Dampers for each thread
971 : MooseObjectWarehouse<NodalDamper> _nodal_dampers;
972 :
973 : /// General Dampers
974 : MooseObjectWarehouse<GeneralDamper> _general_dampers;
975 :
976 : /// NodalKernels for each thread
977 : MooseObjectTagWarehouse<NodalKernelBase> _nodal_kernels;
978 :
979 : /// Decomposition splits
980 : MooseObjectWarehouseBase<Split> _splits; // use base b/c there are no setup methods
981 :
982 : /// Constraints storage object
983 : ConstraintWarehouse _constraints;
984 :
985 : /// increment vector
986 : NumericVector<Number> * _increment_vec;
987 : /// Preconditioner
988 : std::shared_ptr<MoosePreconditioner> _preconditioner;
989 :
990 : /// Whether or not to use a finite differenced preconditioner
991 : bool _use_finite_differenced_preconditioner;
992 :
993 : MatFDColoring _fdcoloring;
994 :
995 : /// The field split preconditioner if this sytem is using one
996 : FieldSplitPreconditionerBase * _fsp;
997 :
998 : /// Whether or not to add implicit geometric couplings to the Jacobian for FDP
999 : bool _add_implicit_geometric_coupling_entries_to_jacobian;
1000 :
1001 : /// Whether or not to assemble the residual and Jacobian after the application of each constraint.
1002 : bool _assemble_constraints_separately;
1003 :
1004 : /// Whether or not a ghosted copy of the residual needs to be made
1005 : bool _need_residual_ghosted;
1006 : /// true if debugging residuals
1007 : bool _debugging_residuals;
1008 :
1009 : /// true if DG is active (optimization reasons)
1010 : bool _doing_dg;
1011 :
1012 : /// vectors that will be zeroed before a residual computation
1013 : std::vector<std::string> _vecs_to_zero_for_residual;
1014 :
1015 : unsigned int _n_iters;
1016 : unsigned int _n_linear_iters;
1017 :
1018 : /// Total number of residual evaluations that have been performed
1019 : unsigned int _n_residual_evaluations;
1020 :
1021 : Real _final_residual;
1022 :
1023 : /// If predictor is active, this is non-NULL
1024 : std::shared_ptr<Predictor> _predictor;
1025 :
1026 : bool _computing_pre_smo_residual;
1027 :
1028 : /// The pre-SMO residual, see setPreSMOResidual for a detailed explanation
1029 : Real _pre_smo_residual;
1030 : /// The initial (i.e., 0th nonlinear iteration) residual, see setPreSMOResidual for a detailed explanation
1031 : Real _initial_residual;
1032 : /// Whether to use the pre-SMO initial residual in the relative convergence check
1033 : bool _use_pre_smo_residual;
1034 :
1035 : bool _print_all_var_norms;
1036 :
1037 : /// If there is any Kernel or IntegratedBC having save_in
1038 : bool _has_save_in;
1039 :
1040 : /// If there is any Kernel or IntegratedBC having diag_save_in
1041 : bool _has_diag_save_in;
1042 :
1043 : /// If there is a nodal BC having save_in
1044 : bool _has_nodalbc_save_in;
1045 :
1046 : /// If there is a nodal BC having diag_save_in
1047 : bool _has_nodalbc_diag_save_in;
1048 :
1049 : void getNodeDofs(dof_id_type node_id, std::vector<dof_id_type> & dofs);
1050 :
1051 : /// Flag used to indicate whether we have already computed the scaling Jacobian
1052 : bool _computed_scaling;
1053 :
1054 : /// Whether the scaling factors should only be computed once at the beginning of the simulation
1055 : /// through an extra Jacobian evaluation. If this is set to false, then the scaling factors will
1056 : /// be computed during an extra Jacobian evaluation at the beginning of every time step.
1057 : bool _compute_scaling_once;
1058 :
1059 : /// The param that indicates the weighting of the residual vs the Jacobian in determining
1060 : /// variable scaling parameters. A value of 1 indicates pure residual-based scaling. A value of 0
1061 : /// indicates pure Jacobian-based scaling
1062 : Real _resid_vs_jac_scaling_param;
1063 :
1064 : /// A container of variable groupings that can be used in scaling calculations. This can be useful
1065 : /// for simulations in which vector-like variables are split into invidual scalar-field components
1066 : /// like for solid/fluid mechanics
1067 : std::vector<std::vector<std::string>> _scaling_group_variables;
1068 :
1069 : /// Container to hold flag if variable is to participate in autoscaling
1070 : std::vector<bool> _variable_autoscaled;
1071 :
1072 : /// A container for variables that do not partipate in autoscaling
1073 : std::vector<std::string> _ignore_variables_for_autoscaling;
1074 :
1075 : /// Whether to include off diagonals when determining automatic scaling factors
1076 : bool _off_diagonals_in_auto_scaling;
1077 :
1078 : /// A diagonal matrix used for computing scaling
1079 : std::unique_ptr<libMesh::DiagonalMatrix<Number>> _scaling_matrix;
1080 :
1081 : private:
1082 : /**
1083 : * Finds the implicit sparsity graph between geometrically related dofs.
1084 : */
1085 : void findImplicitGeometricCouplingEntries(
1086 : GeometricSearchData & geom_search_data,
1087 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> & graph);
1088 :
1089 : /**
1090 : * Setup group scaling containers
1091 : */
1092 : void setupScalingData();
1093 :
1094 : /// Functors for computing undisplaced mortar constraints
1095 : std::unordered_map<std::pair<BoundaryID, BoundaryID>, ComputeMortarFunctor>
1096 : _undisplaced_mortar_functors;
1097 :
1098 : /// Functors for computing displaced mortar constraints
1099 : std::unordered_map<std::pair<BoundaryID, BoundaryID>, ComputeMortarFunctor>
1100 : _displaced_mortar_functors;
1101 :
1102 : /// The current states of the solution (0 = current, 1 = old, etc)
1103 : std::vector<NumericVector<Number> *> _solution_state;
1104 :
1105 : /// Whether we've initialized the automatic scaling data structures
1106 : bool _auto_scaling_initd;
1107 :
1108 : /// A map from variable index to group variable index and it's associated (inverse) scaling factor
1109 : std::unordered_map<unsigned int, unsigned int> _var_to_group_var;
1110 :
1111 : /// The number of scaling groups
1112 : std::size_t _num_scaling_groups;
1113 : };
|
