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