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